SiSSR: Simultaneous subdivision surface registration for the quantification of cardiac function from computed tomography in canines

Highlights

• A pipeline for quantifying cardiac function from computed tomography is proposed.

• Patient-specific meshes modelling cardiac anatomy are automatically constructed.

• Cardiac mesh models are fitted to the endocardium in all frames simultaneously.

• Model fitting is formulated so as to encourage biologically plausible motion.

• Functional cardiac parameters are validated against the prior state-of-the-art.

Abstract

Recent improvements in cardiac computed tomography (CCT) allow for whole-heart functional studies to be acquired at low radiation dose (<2mSv) and high-temporal resolution (<100ms) in a single heart beat. Although the extraction of regional functional information from these images is of great clinical interest, there is a paucity of research into the quantification of regional function from CCT, contrasting with the large body of work in echocardiography and cardiac MR. Here we present the Simultaneous Subdivision Surface Registration (SiSSR) method: a fast, semi-automated image analysis pipeline for quantifying regional function from contrast-enhanced CCT. For each of thirteen adult male canines, we construct an anatomical reference mesh representing the left ventricular (LV) endocardium, obviating the need for a template mesh to be manually sculpted and initialized. We treat this generated mesh as a Loop subdivision surface, and adapt a technique previously described in the context of 3-D echocardiography to register these surfaces to the endocardium efficiently across all cardiac frames simultaneously. Although previous work performs the registration at a single resolution, we observe that subdivision surfaces naturally suggest a multiresolution approach, leading to faster convergence and avoiding local minima. We additionally make two notable changes to the cost function of the optimization, explicitly encouraging plausible biological motion and high mesh quality. Finally, we calculate an accepted functional metric for CCT from the registered surfaces, and compare our results to an alternate state-of-the-art CCT method.

Introduction

Regional cardiac dysfunction commonly precedes changes in global metrics such as ejection fraction (EF), and is useful in the characterization of subtle impairment and monitoring of treatment response (Tee et al., 2013). Techniques developed for cardiac magnetic resonance imaging (CMR), including harmonic phase (HARP) imaging (Osman et al., 1999), strain encoded (SENC) imaging (Osman et al., 2001), and displacement encoding with stimulated echos (DENSE) imaging (Aletras et al., 1999), have enabled the measurement of biomarkers for regional cardiac function, including strain, strain rate, and torsion. However, due to their complex and time-consuming acquisition and analysis, further complicated by contraindication of CMR in patients with pacemakers, these techniques have been confined to research studies in academic centers.

Echocardiographic techniques for measuring function are widely deployed in the clinic due to low cost and real-time, bedside acquisitions. However, in terms of cardiac function, echocardiography is subject to certain limitations. Tissue Doppler echocardiography (TDI) is used to measure cardiac motion, but is highly angle-dependent. Speckle-tracking, which has been developed both for 2-D and 3-D sequences, currently provides measurements which are less sensitive to probe angle than Doppler techniques (Tee et al., 2013). Transthoracic echocardiography (TTE) is non-invasive, but unfortunately is limited by the available acoustic window. Transesophageal echocardiography (TEE) avoids these issues, but is an uncomfortable, semi-invasive procedure. All echocardiographic techniques are limited by operator dependence, and 2-D techniques are additionally subject to through-plane motion errors.

Historically, cardiac computed tomography (CCT) has been clinically infeasible for functional studies due to high radiation dose and low temporal resolution (Suinesiaputra et al., 2016). However, recent advances in CCT have dramatically reduced the necessary radiation dose (<2 mSv) and increased the temporal resolution (<50 ms) for whole-heart functional studies (Tee et al., 2015), spurring interest in the use of CCT in assessing cardiac function. Moreover, CCT is an attractive modality for functional studies given its high spatial resolution (on the order of 0.25 mm), inherently volumetric data, fast scan times (single heart beat), ubiquity in clinic, and increasing inclusion in cardiac imaging protocols due to the depiction of coronary anatomy (Tee et al., 2013). Therefore, fast, automated, and robust algorithms for extracting regional functional information from CCT are of immediate interest.

Initial efforts to recover strain from CCT have largely focused on application of tools originally developed for 2-D ultrasound (Ogawa et al., 2006) to resliced CCT sequences (Helle-Valle, Yu, Fernandes, Rosen, Lima, 2010, Tee, Noble, Bluemke, 2013). This approach is limited, as it relies on myocardial features which are largely absent from CCT and discards useful volumetric information. An alternative metric and approach has been proposed by Pourmorteza et al. (2012), who employed the coherent point drift (CPD) method (Myronenko, Song, Carreira-Perpinan, 2007, Myronenko, Song, 2010) to register a template mesh representing the left ventricular (LV) endocardium at end diastole to each subsequent frame. In CPD, the “moving” point set is represented as a Gaussian Mixture Model (GMM). Using Expectation Maximization (EM), the GMM is deformed so as to maximize the probability of the “fixed” point set having been sampled from it. “Coherent” motion is enforced by constraining the motion of the GMM centroids relative to one another. CPD can be used for non-rigid registration, and does not prescribe any specific transformation model. Having performed the registration, Pourmorteza et al. (2012) calculated for each face the square root of the ratio between the face’s area and the area of the corresponding face at end diastole. This metric, which they termed Stretch QUantifier for Endocardial Engraved Zones (SQUEEZ), may be interpreted as a surrogate for cardiac contraction (<1) or expansion ( > 1) in the circumferential-longitudinal plane. SQUEEZ is well suited to volumetric CCT data, and was subsequently shown to correlate with strain obtained by tagged CMR (Pourmorteza et al., 2015); to have clinical utility in the planning of cardiac resynchronization therapy (Behar et al., 2017); and to be robust to noise in low-dose CCT (Pourmorteza et al., 2018). After detailing the technical contributions of our proposed method, we calculate SQUEEZ and compare our results to those obtained from the CPD method.

The objective of our method is conceptually similar to Pourmorteza et al. (2012) in that we derive SQUEEZ from a series of registered meshes. However, rather than registering meshes using CPD, we adapt and extend a method which has been successfully used for cardiac modeling and segmentation in the context of 3-D echocardiography (Stebbing, 2014, Stebbing, Namburete, Upton, Leeson, Noble, 2015). In Stebbing (2014), a subdivision surface is registered to an arbitrary point set obtained using an off-the-shelf edge detector. A notable advantage of this method is that the subdivision surfaces are registered across all frames simultaneously, contrasting with the conceptually simpler “frame-to-frame” and “reference-frame” formulations. In the frame-to-frame formulation, consecutive pairs of frames are registered, leading to propagation of error throughout the cardiac cycle. In the reference-frame formulation, each frame is registered to the end diastolic frame, resulting in independent motion estimates. While this may be desirable in some scenarios, it may also lead to discontinuities in the displacements measured between successive frames (as the transformation estimated by a registration algorithm is not in general a smooth function of its inputs). By formulating the problem as a single, global optimization, the method constructed by Stebbing (2014) avoids both drawbacks.

In this paper, we present the Simultaneous Subdivision Surface Registration (SiSSR) method, which extends the state-of-the-art method for Loop subdivision surface registration developed by Stebbing (2014) in order to register an anatomical reference mesh representing the endocardium of the LV to a 3-D+time cardiac image sequence (Fig. 1). The most notable differences between SiSSR and its predecessor are:

• An anatomical reference mesh is constructed directly from the (segmented) image data, obviating the need for a bespoke template mesh representing the structure of interest, and accounting for the wide anatomical variability between canines.

• Two notable changes are made to the cost function. First, the velocity, rather than acceleration, of the (coarse) mesh vertices are regularized, more plausibly modeling real cardiac motion. Second, (coarse) mesh faces with high aspect ratios are penalized, thereby explicitly discouraging solutions with low mesh quality.

• Subdivision surfaces naturally lend themselves to a multiscale approach, which may decrease computation time while directing the optimization away from local minima. Therefore, we perform an initial registration, explicitly subdivide the result, and use the subdivided result to initialize a second pass.

• Stebbing (2014) evaluates the method primarily as a segmentation algorithm, and stops short of calculating functional parameters. We calculate SQUEEZ, and compare our results with the CPD approach.

• The subdivision method was originally developed in the context of 3-D echocardiography. Here, we instead investigate 3-D+time CCT, which has received relatively less attention and which we feel will become increasingly common clinically in the near future.

The structure of this paper is as follows. In Section 2 we present our methodology in detail, specifically highlighting points of difference between our and prior methods. In Section 3, we provide technical details pertaining to the implementation. In Section 4, we describe the canine infarct model, imaging parameters, and validation against the CPD method. In Section 5, we present qualitative and quantitative results. In Section 6 we offer a discussion of our findings, its limitations, and future directions. In Section 7, we provide a summary of the paper.

Unless otherwise noted or obvious from context, images and meshes are denoted by uppercase script letters (e.g., $\mathcal{I}$); matrices by uppercase, boldface letters (e.g., X); tuples and vectors by lowercase boldface letters (e.g., t); real and integer scalars by lowercase, unweighted letters (e.g., r); and constants by uppercase, unweighted letters (e.g., F).