Variational Reconstruction of Left Cardiac Structure from CMR Images

Cardiovascular Disease (CVD), accounting for 17% of overall deaths in the USA, is the leading cause of death over the world. Advances in medical imaging techniques make the quantitative assessment of both the anatomy and function of heart possible. The cardiac modeling is an invariable prerequisite for quantitative analysis. In this study, a novel method is proposed to reconstruct the left cardiac structure from multi-planed cardiac magnetic resonance (CMR) images and contours. Routine CMR examination was performed to acquire both long axis and short axis images. Trained technologists delineated the endocardial contours. Multiple sets of two dimensional contours were projected into the three dimensional patient-based coordinate system and registered to each other. The union of the registered point sets was applied a variational surface reconstruction algorithm based on Delaunay triangulation and graph-cuts. The resulting triangulated surfaces were further post-processed. Quantitative evaluation on our method was performed via computing the overlapping ratio between the reconstructed model and the manually delineated long axis contours, which validates our method. We envisage that this method could be used by radiographers and cardiologists to diagnose and assess cardiac function in patients with diverse heart diseases.


Introduction
Cardiovascular Disease (CVD), the leading cause of death all over the world, accounts for 17% of overall deaths in the USA. CVD claims more lives than the combination of the next seven leading causes of death. Modern diagnosis and treatment of CVD strongly depend on the aid of medical imaging tools. Numerous advances in medical imaging techniques such as computed tomography (CT) and cardiovascular magnetic resonance imaging (CMR or cardiac MRI) greatly help clinicians to acquire both morphological and functional information regarding the heart.
Visual assessment based on experience used to be clinicians' primary method to diagnosis medical conditions. High intra-/inter-observer variability of visual assessment made clinicians search quantitative methods with high repeatability and reproducibility. Any such quantitative method has an invariable prerequisite, i.e., an accurate computational cardiac model, the study of which has started since the introduction of the imaging modalities. [1] surveyed the literature of modeling techniques from a variety of imaging modalities.
Most studies of cardiac modeling were focused on the left ventricle (LV). This exclusive preference was due to the physiopathological significance of LV. Early approaches [2] used idealized ellipsoidal shapes and variations [3] to modeling the LV. Those compromise methods were practical in the past, however they are quite oversimplified given nowadays advanced imaging techniques and numerical methods. Similar recent approaches [4][5][6][7] used more sophisticated superquadric shapes with more control points to replace the simple ellipsoidal shapes. Images were fitted to determine the shape parameters. The accuracy of results based on the pre-assumption may be limited when handling the highly variable subjects. A few recent approaches [8][9][10][11] used image data to fit a population-based (statistical) model, which also suffer from the above limitations.
Besides, most reconstructed models are not complete LV models. The basal short axis imaging plane was usually considered as the top of LV models. Two significant parts in the anatomic definition of LV, i.e., the atrioventricular junction (LV inflow tract) and the aorto-ventricular junction (LV outflow tract), are not included in the reconstructed models. However, modern analysis tools such as finite element methods and hemodynamic simulation requires a complete LV model for computation [12][13][14][15][16][17][18]. Two junctions are especially significant for the assessment precision considering the influence from mitral valve and aorta valve to the hemodynamic simulation.
Only a few studies addressed the LV inflow or outflow tract as well as the other chambers and vessels in the last decades. Four-chamber model was reconstructed in [9,19,20] and a model containing more components such as coronary arteries was established in [10]. As comprehensive as they are, these models have two limitations. Firstly, some of them are population-based (statistical) model. To adapt them to patient-specific modeling methods is not a trivial task. Secondly, the imaging modalities in both studies were 3D CT or high resolution MR or Diffusion Tensor MR providing a fair image resolution and 3D isotropy. However, CMR is considered as the gold-standard for assessing cardiac anatomy and function due to its ability to capture multiframes images throughout the cardiac cycle and widely used in current clinical practice.
The major difficulties of the whole left heart modeling based on CMR include the large imaging slice spacing, relating three dimensional anisotropy, and the uncertainty of morphologies of two junction orifices. Researchers in numerical engineering proposed various methods to reconstruct models from contours, which generally use triangulation and meshing technique to determine connections in between neighboring contours [21][22][23][24][25], which however suffer from the inevitable heuristics due to the topology variations or contours sparseness. The demand of a robust and precise reconstruction method motivates this study.
In this study, a novel method is proposed to address the left cardiac modeling problem from a variational approach. Contours were delineated to indicate LV, LA, and AO in both short and long axis images. Contours from different images were registered in the patient-based coordinate system. After a series of preprocessing including contour matching and interpolation, a variational surface reconstruction method was applied. Delaunay based tetrahedral meshes were generated to discretize the underlying space. Graph-cuts were applied to solve the variational problem. The surface model containing the LV, LA, and AO was then extracted from the tetrahedral mesh according to the min-cut. The intersection of the reconstructed surface with the long axis imaging plane was validated against to the manually delineated contours using both curve-based and region-based criteria.
The remainder of this article is organized as follows. Section 2 describes the methodology. Section 3 provides the experimental results and the validation. Section 4 concludes this article.

Methods
In this study, we tested algorithm on ten healthy volunteers. This study was approved by the Singhealth Centralised Institutional Review Board (CIRB No: 2009/705/C) for human research. All enrolled participants gave written informed consent. The MR data are deposited in hospital and are available for research and education purpose. Cardiac relating measurements for each volunteer/patient were given in Table 1.

Image Acquisition and Contour Delineation
A 1.5T Siemens scanner with ECG gating was utilized to acquire the cine MR images. Both long axis and short axis images were acquired. A parallel stack of short axis images was acquired from the left atrium to left ventricle apex (8mm inter-slice thickness with no interslice gap). Three long axis images were acquired orthogonal to the short axis images, i.e., two chamber view, three chamber view, and four chamber view. The TR/TE/flip angle is typically 68ms/1ms/70°. The field of view was typically 320 mm with spatial resolution of less than 1.5 mm (typically 1.43 mm). Each slice was acquired in a single breath hold, with 22 temporal frames per cardiac cycle.
Both long axis and short-axis MR images were processed in the CMRtools suite (Cardiovascular Solution, UK). Endocardium was delineated by experts for the end-diastole (ED) and end-systole (ES) of each image slice.
For the two chamber view, the LV and LA were delineated; for the three chamber view, the LV, LA, and AO were delineated; for the four chamber view, the LV and LA were delineated. Examples of such long axis images and delineations are shown in Fig 1. Identification of the boundaries of LA inflow tract and AO outflow tract could be of less significance. Since the subsequent registration and surface reconstruction procedure would use little information from those.
For short axis images containing only left ventricle, left ventricle endocardium was delineated, see and 2(e). All papillary muscles were excluded from the myocardial region and were considered as the blood pool instead.

Point Cloud Registration
The contours delineated from all images were projected into one three dimensional space, i.e., the patient-based coordinate system, according to the imaging specification, e.g., the pixel spacing, the image position, and the image orientation. These image specification are contained in the DICOM file meta information, and the transformation from 2D planar contours to 3D point clouds is as follows.
• C 2ch : the set from the contour on the two chamber left view image.
• C 3ch : the set from the contour on the three chamber view image.
• C 4ch : the set from the contour on the four chamber view image.
Simply collecting all contour points for reconstruction is by no means a robust method, considering the breath hold motion and other movements. Certain registration measures must be done before incorporating all contours from different imaging planes.
Recently, the registration methods addressing the misalignment problem could be categorized into two groups. One group of approaches is image-based approaches [26]. The other is geometry-based approaches [27][28][29]. The image-based approaches are inherently inaccurate due to the large slice spacing and the complex nature of images such as the inhomogeneity and non-uniformity as well as the existence of papillary muscles. Hence, we choose the geometrybased registration method in this stage.
Iterative Closest Point (ICP) algorithm proposed by Besl and McKay [30] and its variations are widely used to register two sets of points, which partially overlap and are usually sampled from one surface.
Directly applying the point cloud registration method to a pair of point clouds may cause artifact registering configuration due to the different data characteristics. Most point cloud registration algortithms were designed for two clouds with relatively dense overlappings, which was not true in our situation. To tackle this issue, we modified the algorithm. Instead of using the whole point clouds, we compute a subset for each point cloud for registration.
The subset filtering criteria are based on the distance between two point clouds. Let C X and C Y be two point clouds, which have limited overlapping part.
The resulting C 0 X and C 0 Y have the following properties: From Property 2 and 3, the relationship between a proper and the Haursdorf distance between the originial two point clouds is implied. Larger generates larger subsets. Hence, in our experiment, we use = αd H (C X , C Y ),0.2 α 0.5. The empirical α is used for the cases when the motion is not significant. The adaptive setting of the α parameter taking into account of the initial position of two point clouds is still in study. The algorithm to compute C 0 X and C 0 Y is described in Table 2.
After computing of such subsets C 0 X and C 0 Y , the classic ICP was applied to them. Assume C 0 X is fixed as the reference, the registration configuration could be obtained then: translation v, and rotation R. The registered subset C 0 Y could be represented as T v;R ðC 0 Y Þ. For the four obtained point clouds, i.e., C sax , C 2ch , C 3ch , and C 4ch , the above registration was performed between three pairs among them, namely {C sax , C 2ch }, {C sax , C 3ch }, and {C sax , C 4ch }. C sax was used for reference in each registration procedure. Note that the selected subset C 0 sax for C sax could be different for different procedure.
Three registration configuration were obtained then: (v 2 , R 2 ) for {C sax , C 2ch }, (v 3 , R 3 ) for {C sax , C 3ch }, and (v 4 , R 4 ) for {C sax , C 4ch }. This registration parameters was defined in earlier study. Briefly, R is the rotation matrix and v is the translate vector. The registration configurations were applied on each long axis point set respectively, obtaining C 0 2ch , C 0 3ch , and C 0 4ch , see Fig  4, where three pairs of registrations are shown respectively as well as a zoomed view. Table 2. Subset computation for two point clouds.

1.
Two point clouds C X , C Y

Variational Surface Reconstruction
The registered point clouds from long/short axis images were then utilized to reconstruct an endocardial surface of the left heart. The reconstruction task consisted of three steps: (a) interpolation between parallel contour points; (b) tetrahedral mesh generation; (c) variational mesh segmentation and surface extraction. Interpolation between Parallel Contour Points. Short axis contours C sax were identified as LV, LA, and AO contours C LV , C LA , C AO . The interpolation between parallel contours, i.e., short axis contours, was conducted in each contour group, respectively. The region of the left ventricular inflow tract and outflow tract, i.e., the bifurcation portion, was not interpolated due to the uncertainty of the correspondence between contour points.
The interpolation of contours in each contour group includes (a) contour re-orientation; (b) intra-contour interpolation; (c) contour matching; (d) inter-contour interpolation.
(a). The re-orientation step was adopted to ensure that all contours to be interpolated were counter-clockwise. This could be accomplished by a number of methods. In our approach, we adjusted all contours to be counter-clockwise by making their signed areas to be positive.
(b). The intra-contour interpolation was conducted to make all contours contain the same number of contour points. The interpolation method was chosen as piecewise cubic Hermite interpolating polynomial (PCHIP). Considering the fact that the contours of LV, LA, or AO are nearly circular shapes, smoother interpolation results with C 2 from a more computationally expensive interpolation method such as cubic spline could not produce a significantly different result compared with PCHIP method.
(c). Contour matching was applied to establish the point-to-point correspondence between two neighboring contours. To obtain a reasonable and most likely correct correspondence between points of two contours, one contour was used as the reference contour and a circular shift was applied on the other one such that the mean of the point-wise distance between two contours was minimized. After the contour matching was applied to all neighboring contours sequentially, all contours in a group (LV, LA, or AO) were reformulated as the format in Table 3.
In Table 3, C 1 , C 2 , . . ., C L are referred as the L parallel contours in a contour group. Each contour contains the same number of contour points, i.e., {C i1 , C i2 , Á Á Á, C iN }. Point correspondence was established in each column, i.e., contour points with the same second subscript {C 1j , C 2j , Á Á Á, C Lj }.
(d). Using the "vertical" correspondence shown in Table 3, inter-contour interpolation was then conducted in each vertical corresponded point set, i.e., {C 1j , C 2j , Á Á Á, C Lj }. The interpolation method was also chosen as PCHIP method. After this inter-contour interpolation, the contour point matrix in Table 3 was enlarged from L × N to M × N, where M > N.
Note that the intra-/inter-contour interpolation size, i.e., M and N was determined by the size/density of the tetrahedral mesh generated in the next step. Too large M and N could produce an almost same result as that of proper M and N, while introducing extra computation load.
The whole interpolation procedure was conducted on each short axis contour group C LV , C LA , and C AO , obtaining three interpolated and re-formulated point sets C 0 LV , C 0 LA , and C 0 AO . This interpolation step is shown in Fig 5. Tetrahedral Mesh Generation. The registered long axis contour points and interpolated short axis contour points form a new point set profiling the whole left cardiac model with a higher point density:  Table 3. Re-formulated contour points. (b) illustrates this preparation for mesh generation, in which the point cloud C new is annotated in red, while the auxiliary point is in bright yellow. The selection of auxiliary grid points was described in our previous work [31], where the usage of a Delaunay-based mesh was also justified.  Variational Mesh Segmentation and Surface Extraction. To reconstruct a triangular mesh surface from the tetrahedral mesh is equivalent to segmenting the tetrahedral mesh into two partitions, interior and exterior. Such task could be addressed as a variational problem, i.e., the weighted minimal surface energy [32].

Parallel contours
where d(x, C new ) = min y 2 C new d(x, y), d(x, y) is the Euclidean distance between x and y. The surface S minimizing this energy functional is the reconstructed surface.
After discretizing this energy functional on the underlying mesh space, it was noted that this minimization problem could be solved by graph-cuts technique [33], i.e., max-flow/mincut algorithm. In other words, this energy is graph-representable. Solving this problem via graph-cuts has some more technical details to concern, such as determining the solution space and establishing a proper boundary condition for the solver. The operation procedure was also described in our previous study. In this section, we illustrate these steps in Fig 7. Applying the graph-cuts on the problem, a min-cut was then obtained efficiently. A triangular surface mesh was then extracted from the tetrahedral mesh according to the min-cut. After certain post-processing such as smoothing [34] and remeshing [35], a processed left cardiac surface was obtained and an example is shown in Fig 8.

Results
The present method was applied on data from ten volunteers. Cardiac relating measurements for each volunteer were give in Table 1. The average time to reconstruct the left cardiac model for one frame is around 6 seconds on a 2.5GHz CPU Desktop. ED and ES frames were reconstructed for each case. The reconstructed surface model of Volunteer 1, ED frame was shown in Fig 8. All ten reconstructed models were shown in Fig 9. The gold standard to reconstruct the left cardiac model from MRI is not easily accessible. Commercial softwares such as CMRTools usually reconstruct the sole LV model, which is relatively trivial task. Alternatively we use the accurate manually delineated contours to validate the reconstructed model. Since contours were extracted from two-dimensional imaging planes, we projected the reconstructed model to the corresponding planes to produce the sectional contours. Both short and long axis contours should have been used for validation. Since the point cloud contains much more short axis contours than the long axis contours, and interpolation was performed in between short axis contours, the validation result between short axis contours and reconstructed models are extremely good, almost coinciding with each other. Hence, we choose the long axis contours as the validation references, which could evaluate the method more critically.
The overlapping ratio between the reconstructed model and long axis contours was used to evaluate the reliability and accuracy of the reconstruction method. The overlapping ratio is widely used in evaluating the performance of segmentation methods. In our study, the reconstruction task could be thought of as segmenting the left cardiac endocardial cavity in a highly anisotropic 3D image (slice spacing up to 8mm). The long axis contours could serve as the ground truth, against which the reconstruction results could be compared with.
In the experiments, the intersection of the reconstructed surface model and the long axis imaging planes were computed and validated against the contours drawn by experts in the beginning of the experiments, see   Three criteria were utilized in the evaluation, i.e., Hausdorff distance, Dice similarity Coefficient, Jaccard similarity coefficient. Hausdorff distance is a curve-based coefficient to measure the furthest displacement from the reconstructed model to the ground truth contour.
Meanwhile, the Dice and Jaccard similarity coefficients are region-based measurements of the overlapping ratio between the reconstructed model and the ground truth contour. The Dice (D) and Jaccard (J) coefficients are defined as follows.
where Re and Tr are areas bounded by the reconstructed model and the manual delineated contour, respectively. A value of 0.7 and above is considered an adequate overlap. An example of the validation against two chamber view long axis contour is shown in Fig 11. The statistics are given in Table 4 for the validation results of all ten volunteers. The average of all cases are 8.17 for Hausdorff distance, 0.96 for Dice coefficients, and 0.92 for Jaccard coefficients. These results validate our method as a reliable and accurate reconstruction tool.
This validation was not performed for the reconstructed model from un-registered point clouds (skipping the step in Section 2.2) to evaluate the impact of the registration on the accuracy. The registration step was designed to correct the motion relating to the breath hold. For the ten healthy volunteers in this study, the motions were relatively subtle. The impact of the registration step is not obvious in this study. In our earlier study, we have encountered data with large distortion from breath hold motions, which critically undermined the accuracy of function assessment. We proposed the method to recover the distorted shape in [29]. The registration step in this article could eliminate this distortion before shape modeling. We will validate this in future works.
Implication: Comprehensive reconstruction of left cardiac structure from CMR images holds several potential implications in heart modelling and function assessment. First, the  comprehensive reconstruction of left cardiac structure including left ventricle, left atrium and aorta will straightforward facilitate the downstream heart flow simulation, as reported in our previous publication [12][13][14][15][16][17][18]. Second, the reconstruction of left cardiac structure will facilitate rapid automated numerical characterization of point heart surface curvedness and thereby point function spatial-temporal fluctuations of curvedness reflect local heart muscle contraction, allowing comprehensive annotation of regional and global cardiac left ventricular structure and function [36][37][38][39] and aorto-ventricular matching before and after heart surgery [40]. Last, this method will contribute a sizeable contemporary CMR imaging atlas of left heart morphology and function in normal subjects and diverse diseased hearts in the near future [41]. These models from CMR images will also be used to validate LV three dimensional echocardiography measurements (i.e., LV volumes and ejection fraction) (3D-echo).

Conclusions
In this study, we proposed a novel method to semi-automatically reconstruct the whole left cardiac surface model from contours in short and long axis CMR images. Contours from multiple images were registered to each other in the patient coordinate system. Intra/Inter-contour interpolations were performed on the parallel short axis contours then. The resulting contour points were utilized in a variational surface reconstruction via Delaunay triangulation and graph-cuts. The method was validated via evaluating the overlapping ratio between the reconstructed model and multiple long axis contours on the corresponding image planes. High overlapping ratio indicates a reliable reconstruction of the left cardiac structure. This validation is an alternative method to evaluate our reconstruction method when the gold standard of the left cardiac model was not available in our study. More comprehensive validation will be studied in our future works, including validation against the model reconstructed from CT scanning of the same subjects (which has a higher spatial resolution), and against various medical measurements such as ejection fraction (EF) from CT or cardiac catheterization, or nuclear imaging, the prerequisite of which is an accurate localization and modelling of the mitral valve and aortic valve.