Feature Tracking Cardiac Magnetic Resonance via Deep Learning and Spline Optimization

TL;DR

This paper introduces a fully automated deep learning and spline optimization method for feature tracking in cardiac MRI, enabling accurate segmentation and strain analysis for cardiac function assessment.

Contribution

It combines CNN-based segmentation with spline fitting for automated, reproducible cardiac feature tracking in MRI, improving over manual methods.

Findings

Segmentation accuracy ~97% pixel accuracy

IoU ~87% across classes

Significant difference in strain between control and HCM patients

Abstract

Feature tracking Cardiac Magnetic Resonance (CMR) has recently emerged as an area of interest for quantification of regional cardiac function from balanced, steady state free precession (SSFP) cine sequences. However, currently available techniques lack full automation, limiting reproducibility. We propose a fully automated technique whereby a CMR image sequence is first segmented with a deep, fully convolutional neural network (CNN) architecture, and quadratic basis splines are fitted simultaneously across all cardiac frames using least squares optimization. Experiments are performed using data from 42 patients with hypertrophic cardiomyopathy (HCM) and 21 healthy control subjects. In terms of segmentation, we compared state-of-the-art CNN frameworks, U-Net and dilated convolution architectures, with and without temporal context, using cross validation with three folds. Performance relative to expert manual segmentation was similar across all networks: pixel accuracy was ~97%, intersection-over-union (IoU) across all classes was ~87%, and IoU across foreground classes only was ~85%. Endocardial left ventricular circumferential strain calculated from the proposed pipeline was significantly different in control and disease subjects (-25.3% vs -29.1%, p = 0.006), in agreement with the current clinical literature.

Figures

Tables

Table 1. Table 1: CNN architecture variants considered. Note: ED = End Diastole; DC = Dilated Convolution.

Name	Variant	Input Size	Temporal Context
Network A	U-Net	$128\times 128\times 1$	None
Network B	DC	$128\times 128\times 1$	None
Network C	U-Net	$128\times 128\times 2$	ED Frame
Network D	DC	$128\times 128\times 2$	ED Frame

Table 2. Table 2: Network Performance Compared with Expert Manual Segmentations.

Name	Description	Metric	Pixel Accuracy	IoU (All)	IoU (Foreground)
Network A	U-Net, No Context	Mean	0.977	0.874	0.838
		Median	0.981	0.885	0.851
Network B	DC, No Context	Mean	0.977	0.876	0.840
		Median	0.981	0.886	0.853
Network C	U-Net, Context	Mean	0.976	0.874	0.838
		Median	0.981	0.887	0.855
Network D	DC, Context	Mean	0.976	0.873	0.837
		Median	0.981	0.888	0.855

Table 3. Table 3: Circumferential Strain Results. ( † : Significant at the p < 0.05 𝑝 0.05 p<0.05 level.)

Plane	Structure	Control ($\%$)	Overt ($\%$)	p
Base	LV Endocardium	$-26.3$	$-29.0$	$0.098$
	LV Epicardium	$-11.8$	$-11.7$	$0.952$
	RV Endocardium	$-11.9$	$-14.3$	$0.529$
Midslice	LV Endocardium	$-25.3$	$-29.1$	$0.006{}^{†}$
	LV Epicardium	$-8.8$	$-8.5$	$0.661$
	RV Endocardium	$-8.0$	$-12.7$	$0.152$
Apex	LV Endocardium	$-25.5$	$-28.6$	$0.110$
	LV Epicardium	$-10.7$	$-9.9$	$0.470$
	RV Endocardium	$-13.8$	$-12.4$	$0.566$

Equations

$$B_d(r)=\sum _{i=0}^d{b}_{i,d}(r)x_i,$$

$$b_{i,d}(r)=\left(\genfrac{}{}{0ex}{}{n}{i}\right){\left(1-r\right)}^{n-i}{r}^i,$$

$$\left(\genfrac{}{}{0ex}{}{n}{i}\right)=\frac{n!}{i!(n-i)!}.$$

$$B_2(r)=(1-r{)}^2{x}_0+2r(1-r)x_1+r^2{x}_2,$$

$$\begin{gathered} B_2(r)=\left[\begin{array}{ccc}{x}_0& x_1& x_2\end{array}\right]\left[\begin{array}{c}(1-r{)}^2\\ 2r(1-r)\\ r^2\end{array}\right] \\ =\left[\begin{array}{ccc}{x}_0& x_1& x_2\end{array}\right]\left[\begin{array}{ccc}1& -2& 1\\ 0& 2& -2\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}1\\ r\\ r^2\end{array}\right]. \end{gathered}$$

$$E_{cf}=\sum _{f,j,r}\Vert {\mathbf{u}}_{f,j,r}-\varphi \left({\mathbf{u}}_{f,j,r}\right){\Vert }^2.$$

$$E_{ac}=\sum _{f,c}{\Vert \left[\begin{array}{ccc}1& -2& 1\end{array}\right]\left[\begin{array}{c}{\mathbf{x}}_{(f+2)\mathrm{mod}K,c}^{\top }\\ {\mathbf{x}}_{(f+1)\mathrm{mod}K,c}^{\top }\\ {\mathbf{x}}_{f,c}^{\top }\end{array}\right]\Vert }^2$$

$$E_{cv}=\sum _{f,c}{\Vert \frac{d^2{B}_2(r)}{d{r}^2}\Vert }^2$$

$$E=\min \left({\rho }_{cf}{E}_{cf}+{\rho }_{ac}{E}_{ac}+{\rho }_{cv}{E}_{cv}\right)$$

Full text

1 Introduction

Quantification of regional cardiac function is of utmost importance in the characterization of subtle abnormalities which may precede changes in global metrics [1]. Harmonic phase (HARP) analysis [2] of tagged cardiac magnetic resonance (CMR) images is the gold-standard for regional function, but has not been adopted clinically due to lengthy acquisition and analysis. Moreover, HARP analysis is difficult to apply to chambers other than the left ventricle (LV) due to the thinness of the myocardial wall in the atria and right ventricle (RV). Recently, feature tracking (FT) has emerged as a promising alternative to tagging [3]. Because FT-CMR can be applied to balanced steady state free precession (SSFP) images, acquisition of specialized image sequences is avoided. Moreover, because FT-CMR primarily tracks myocardial borders and trabeculation, the thinness of the atrial and right ventricular myocardium does not hinder tracking.

Despite recent interest in FT-CMR, two important challenges remain. First, all current commercially available implementations (MTT, TomTec, CMR42) require manual contouring of one or more cine frames, preventing full automation and reducing reproducibility. Second, FT generally has been implemented by repeatedly applying methods designed to determine a displacement field between a single pair of images (e.g., optical flow, block matching, deformable registration), rather than an image sequence; either matching successive pairs, or matching each frame to a single reference, typically end diastole (ED, Fig 1). Each of these approaches has well-known potential drawbacks [4], which may be overcome by empirically optimizing over all frames simultaneously [5].

Here, we propose a method for FT-CMR analysis which overcomes the first of these challenges by using a deep learning approach in place of human contouring. Deep convolutional neural networks have been used to great effect in image classification [6, 7], and semantic segmentation [8, 9]. Recently, CNNs have also shown state-of-the-art performance in biomedical image analysis [10, 11]. CNN segmentation of short axis CMR has been applied to the LV blood-pool [12], the RV blood-pool [13], and both simultaneously [14]. Here, we perform segmentation of the LV myocardium, LV blood-pool, and RV blood-pool. Moreover, we apply the segmentation to patients with HCM, which increases the complexity of the problem due to the highly variable appearance of the LV in these patients.

The second of these challenges we address by fitting quadratic basis splines to the segmentation data jointly, rather than frame by frame, adapting the technique presented in [15, 5] in the context of cardiac ultrasound. In our pipeline, extraction is performed using a CNN, tracking is performed with simultaneous spline optimization, and cardiac strain is estimated from the registered splines.

2 Methods

Broadly, the automated analysis pipeline involves three steps: feature extraction (segmentation), feature tracking (spline fitting), and calculation of functional parameters (strain estimation). These steps are discussed in detail in the following sub-sections.

2.1 Segmentation

Following the work of [8, 10, 11], we designed our segmentation architecture as a fully convolutional network. In order to obtain a segmentation map with the same spatial resolution as the input image, up-sampling operators are used to replace the pooling operators in traditional classification networks. This strategy enables our network to segment arbitrarily large images.

As shown in Fig. 2, the architecture consists of a down-sampling path (left) followed by an up-sampling path (right). During the first several layers, the structure resembles the canonical classification CNN [6, 7], as a $3\times 3$ convolution, rectified linear unit (ReLU), and $2\times 2$ max pooling are repeatedly applied to the input image and feature maps. In the second half of the architecture, we “undo” the reduction in spatial resolution by performing $2\times 2$ up-sampling, ReLU activation, and $3\times 3$ convolution, eventually mapping the intermediate feature representation back to the original resolution. To provide accurate boundary localization, low-level feature representations from the down-sampling path are concatenated with the feature maps from the up-sampling path. For all layers, we apply 128 trainable kernels. We performed batch normalization [16], which has been shown to increase generalizability, between each pair of convolution and ReLU activation layers.

In addition to this basic architecture (Network A), we varied the amount of temporal context by either inputting the input image alone, or the input image and ED image together. We based this on the intuition that the papillary muscles, which frequently interfere with LV segmentation, are least compacted at ED and may guide the segmentation of the input frame. Additionally, in Networks B and D, the final up-sampling/down-sampling pass was replaced by a dilated convolution (DC). All architectures have $\sim 3.1$ million trainable parameters. The architectures tested are summarized in Table 1.

2.2 Quadratic Bézier Curve Registration

Following semantic segmentation, contours defining the boundaries of the LV endocardium, LV epicardium, and RV endocardium were extracted through standard morphological operations. In this work, the pixels belonging to these contours are known as “boundary candidates.” Unfortunately, these contours cannot be used directly to quantify cardiac function, because they lack anatomical correspondence between frames. It is the aim of this section to describe an optimization procedure for jointly registering a sequence of closed, quadratic Bézier curves to these boundary candidates.

A segment of a closed Bézier curve $B$ of degree $d$ parameterized by $r\in[0,1]$ is a linear combination of $d+1$ control points $x_{i}:0\leq i\leq d$,

[TABLE]

where $b_{i,d}(r)$ is the $i$th Bernstein polynomial of degree $d$,

[TABLE]

and $\binom{n}{i}$, often read aloud as “$n$ choose $i$”, is the binomial coefficient,

[TABLE]

For $d=2$, $B_{2}(r)$ expands to

[TABLE]

which may be more conveniently expressed in terms of the monomial basis as

[TABLE]

Importantly, the first and second derivatives of $B_{2}(r)$ with respect to $r$ are trivial to compute.

2.3 Formulating the Optimization

Levenberg-Marquardt least squares optimization [17] is used to register a set of closed, quadratic Bézier curves to the boundary candidates. The parameters $\Delta X\in\mathbb{R}^{2\times\left(C\times K\right)}$ (where $C$ is the number of control points in a single curve and $K$ is the number of cardiac phases) of the optimization are Cartesian displacements to the control points of all template curves across all frames. A fixed number of points $\mathbf{u}_{f,j,r}$ were sampled across each curve. At each step in the optimization, for each of these points, the nearest boundary candidate $\phi\left(\mathbf{u}_{f,j,r}\right)$ was calculated, where $\phi:\mathbb{R}^{2}\rightarrow\mathbb{R}^{2}$. This was computed efficiently by representing the boundary candidate point set at each frame as a $K_{d}$ tree. The Cartesian components of the distance between the points sampled from the curve and nearest boundary candidate were the residuals of $E_{cf}$, the first term of the cost function,

[TABLE]

Additionally, two regularizers were included to enforce physical constraints of anatomical deformation: control point acceleration and spline curvature.

In our cost function, the control point acceleration regularizer allows information to be shared between frames. At a minimum, regularizing against control point velocity as in [5] is necessary to maintain anatomical consistency (the assumption that, for fixed $j$ and $r$, $\mathbf{u}_{f,j,r}$ corresponds to the same material point $\forall f$). By regularizing against acceleration rather than velocity, our method additionally encourages smooth, biologically plausible motion. The control point acceleration regularizer $E_{ac}$ was defined as the Cartesian components of the second differences between corresponding vertices $\mathbf{x}_{f,c}$ in three adjacent frames, where $\mathbf{x}_{f,c}$ is the $(f\times C)+c$th column of $X$.

[TABLE]

Curvature for segment $j$ in frame $f$ is the second derivative of $B_{2}(r)$ with respect to $r$.

[TABLE]

The overall optimization problem may then be written in terms of Eqns. 1, 2, and 3 and corresponding scaling factors. Scaling factors $\rho_{cf}=10.0$, $\rho_{ac}=1.0$, and $\rho_{cv}=0.1$ were set empirically to prevent any single term from dominating the optimization.

[TABLE]

Two points relating to computational efficiency are worth noting. First, the Jacobians of Eqns. 1, 2, and 3 can all be calculated analytically. By providing explicit Jacobians, we avoid the need for numeric derivatives, which would slow computation precipitously. Moreover, for a given set of correspondences between surface positions and boundary candidates, the Jacobians of all residuals are linear with respect to the Cartesian displacements of the control points and therefore trivial to calculate. Second, each individual residual depends upon a very small number of parameters. Specifically, each individual residual depends on exactly three control points (six parameters). This sparsity is exploited during the optimization to limit the number of components of the Jacobian which must be evaluated, further reducing computational cost.

Following the initial fit, the spline is subdivided and used to initialize a second optimization, and this process is repeated one further time. This multiresolution approach has benefit over registering the highly subdivided spline directly, which can be sensitive to initialization.

3 Experiments

3.1 Segmentation

The LV myocardium, LV blood-pool, and RV blood-pool were manually segmented in 189 short axis 2D+time volumes (basal, equatorial, and apical cine series from each of 63 subjects). The papillary muscles of the LV were excluded from the myocardium. The subjects were partitioned into three folds of approximately equal size such that the images from any one subject were present in one fold only. The volumes were cropped to $128\times 128$ pixels in the spatial dimensions, and varied from $25{}$ to $50{}$ pixels in the time dimension, totaling $2706{}$, $3000{}$, and $2775{}$ images in the three folds, respectively. For each of the four architectures (U-Net and DC with and without temporal information), three models were trained on two folds and tested on the remaining fold. The images were histogram equalized and normalized to zero mean and unit standard deviation before being input into the CNN. The network weights were initialized with orthogonal weights [18], and were trained with standard stochastic gradient descent (SGD) with momentum (0.9) by optimizing categorical cross entropy. Learning rate was initialized to $0.01$ and decayed by $0.1$ every $32$ epochs. To avoid over-fitting, we used considerable data augmentation (horizontal and vertical flipping, random translations and rotation) and a weight decay of $10^{-4}{}$. Accuracy was measured as pixel accuracy between the prediction and manual segmentations. The model was implemented in the Python programming language using the Keras interface to Tensorflow [19], and trained on one NVIDIA Titan X graphics processing units (GPU) with 12 GB of memory. For all network architectures, it took roughly 200 seconds to iterate over the entire training set (1 epoch). At test time, the network predicted segmentations at roughly 75 frames per second (real-time).

3.2 Tracking

Short axis (SA) scans from 42 subjects with overt hypertrophic cardiomyopathy (HCM) and 21 control subjects were segmented as described above. For each scan, the LV endocardium, LV epicardium, and RV endocardium were tracked using the spline optimization method. The tracking algorithm was implemented in the C++ programming language using the Insight Toolkit (ITK) for reading, writing, and manipulating images and point sets, and using the Ceres Solver for least squares optimization. The three registration passes took $\sim 2.1s$ per cine sequence.

3.3 Regional Function

Following tracking, registered splines from the 63 subjects were used to calculate global strain. For each structure in each SA plane, global strain was compared between HCM and control subjects using the Student’s t test.

4 Results

In terms of segmentation, performance relative to expert manual segmentation was similar across all networks: mean pixel accuracy was $\sim 97\%$, intersection-over-union (IoU) across all classes was $\sim 87\%$, and IoU across foreground classes only was $\sim 84\%$. However, inspection of the images revealed that a single subject with severe, nonuniform illumination was incorrectly segmented by all networks, with a disproportionate effect on mean performance metrics. For this reason, median values are also reported. Broadly, performance improved with the addition of temporal context over the target frame alone, and with the dilated convolution (DC) networks compared with the U-Net networks (Table 2). Therefore, Network D was selected to provide segmentations for the tracking data.

Tracking was performed in three passes (Fig. 3), where the output of one pass was subdivided and passed to the next pass. In each pass, the contours tighten towards the segmentation, allowing for acute structures such as the RV insertion points to be better described. Compared with registering a highly subdivided spline directly, this technique avoids local minima and converges faster.

The relationship between strain values measured in the control and overt groups was consistent with other studies [20]. In particular, circumferential strain measured in the equatorial LV endocardium was higher (more negative) in overt subjects relative to control subjects (-29.1 vs -25.3, $p=$ 0.006). Detailed circumferential strain results are given in Table 3.

Representative segmentation and tracking results are shown for control and overt subjects (Fig. 4). The model learned to avoid the papillary muscles of the LV myocardium and performed well even in subjects with severe hypertrophy. Tracking visually followed the contours of the segmentation closely.

5 Conclusions

Measuring cardiac function in a fully automated way from SSFP CMR has the potential to simplify quantification of regional cardiac function, and expedite clinical adoption. We have presented a fully automated pipeline for cardiac segmentation, tracking, and estimation of cardiac strain. We obtained segmentation results with and without temporal context in U-Net and DC architectures, and found improvements with temporal context, as well as in DC architectures. The best-performing architecture (Network D) had a median pixel accuracy of $0.981{}$, all-class IoU of $0.887{}$, and foreground IoU of $0.855{}$. We then presented a feature tracking algorithm taking these segmentations as input and jointly optimizing a set of quadratic splines over all frames simultaneously. We applied this segmentation and tracking to the LV endocardium, LV epicardium, and RV endocardium of healthy and disease subjects, and found statistically significant differences between control and overt HCM subjects consistent with previous studies.

Our algorithm is novel in three principal ways. First, in terms of application, the wide anatomical variability observed in subjects with overt HCM make segmentation a particularly difficult problem; this work is the first to demonstrate that CNN-based segmentation is effective in these subjects. In addition, we have directly compared dilated convolution and U-Net architectures to select an appropriate state-of-the-art architecture solution to this problem. Second, a persistent problem in FT-CMR is the interference of the papillary muscles in cardiac segmentation; we have demonstrated that deep learning neatly solves this problem, and are unaware of deep learning segmentation being used for FT-CMR before. Third, the spline optimization method presented avoids the errors inherent to the various pairwise sequential and reference frame formulations ubiquitous in the feature tracking literature to date.

Notably, all networks failed to segment a single case with severe nonuniform illumination. Augmentation during training to counteract this effect will be the subject of future work. Moreover, because only edge features are tracked, our method suffers from the so-called “aperture-problem,” such that anatomical correspondence may not be reliable. In future work, we will incorporate features from our pre-trained CNN into the spline optimization to mitigate this effect. However, this may be a fundamental limitation of FT-CMR where trabeculation is minimal, such as when measuring LV endocardial strain in the basal slice.

In conclusion, we have presented a fully automated pipeline which addresses a number of longstanding challenges to the adoption of FT-CMR, and tested this pipeline successfully in the context of a difficult clinical problem.

6 Acknowledgements

D. Vigneault is supported by the NIH-Oxford Scholars Program and the NIH Intramural Research Program. W. Xie is supported by the Google DeepMind Scholarship, and the EPSRC Programme Grant Seebibyte EP/M013774/1.