Automatic estimation of aortic and mitral valve displacements in dynamic CTA with 4D graph-cuts

The location of the mitral and aortic valves in dynamic cardiac imaging is useful for extracting functional derived parameters such as ejection fraction, valve excursions, and global longitudinal strain, and when performing anatomical structures tracking using slice following or valve intervention’s planning. Completely automatic segmentation methods are still challenging tasks because of their fast movements and the different positions that prevent good visibility of the leaflets along the full cardiac cycle. In this article, we propose a processing pipeline to track the displacement of the aortic and mitral valve annuli from high-resolution cardiac four-dimensional computed tomographic angiography (4D-CTA). The proposed method is based on the dynamic separation of left ventricle, left atrium and aorta using statistical shape modeling and an energy minimization algorithm based on graph-cuts and has been evaluated on a set of 15 electrocardiography-gated 4D-CTAs. We report a mean agreement distance between manual annotations and our proposed method of 2.52±1.06 mm for the mitral annulus and 2.00±0.69 mm for the aortic valve annulus based on valve locations detected from manual anatomical landmarks. In addition, we show the effect of detecting the valvular planes on derived functional parameters (ejection fraction, global longitudinal strain, and excursions of the mitral and aortic valves). Fifteen ECG-Gated 4D-CTA studies processed. the pixel size, field of and cardiac cycle duration of these studies. IDs #1 to #11 are healhy subjects, while IDs #12 to #15 are patients with low ventricular


Introduction
Detecting and tracking cardiac valve annuli and cardiac valvular planes in dynamic cardiac imaging is a necessary task to derive functional parameters, such as the ventricular ejection fraction, valve excursions, or the global longitudinal strain. In addition, it allows slice following to track dynamically the same anatomic slice through the cardiac cycle and allows the automatic selection of data around the valves, e.g., for surgery planning, or for threedimensional (3D) printing of a model of the valve (Giannopoulos et al., 2016;Vukicevic et al., 2017). Furthermore, the geometry, tracking, and movement analysis of the cardiac valves are crucial to assess the functional implications of several cardiovascular conditions in early stages, such as hypertension and dilated cardiomyopathy (Carlsson et al., 2007;Helein and Gibson, 1999;Romano et al., 2018).
To provide an accurate assessment, the cardiac valves must be anatomically and functionally characterized. This can be achieved by applying a heart model that integrates the anatomy of the valve located by the best-fitting plane as has been widely used in manual and semiautomatic segmentation and quantification protocols (Lang et al., 2015). A more elaborate solution would be a dynamic model that includes the complete anatomy of the valve, including the movements of leaflets along the entire cycle (Sun et al., 2014).
Regarding the mitral valve, an intermediate approach is the modeling of the mitral annulus, defined as the junctional zone providing attachment to the mitral valve, separating the left atrium (LA) from the left ventricle (LV) (McCarthy et al., 2010). The mitral annulus is traditionally modeled as a saddle-shaped ring (Levine et al., 1989), although alternatives such as the D-shaped model (Blanke et al., 2014) have also been proposed. As with the mitral valve, there are different approaches to modeling the aortic valve. It is commonly approximated by a ring delimited by different anatomical structures, such as the sinotubular junction, ventriculo-arterial junction, or virtual ring joining the hinge points (Piazza et al., 2008). Furthermore, the aortic valve annulus is also described as the semilunar crown-like structure demarcated by the hinges of the leaflets (Piazza et al., 2008). Regardless of the approach taken to model and characterize them, accurate detection and delineation of mitral and aortic valves in dynamic cardiac imaging is a complex task because they contain rapidly moving structures. Specifically, valve leaflets experience rapid movements through the cardiac phases that are not always visible.
The challenge of segmenting or tracking cardiac valves in 3D dynamic images (denoted as 3D+t or 4D) can be found as an individual objective or as a part of a more general proposal that also encompasses the automatic segmentation of heart chambers. Because mitral and aortic valves limit the upper part of the LV, a segmentation of the LV leads directly to the associated valvular model. However, generally, LV models do not follow the LV's complex shape in the area near the valves (Buckberg et al., 2008); thus, the LV outflow tract (LVOT) is left outside the model. Therefore, the state-of-the-art methods for shape modeling the LV do not track mitral and aortic valves independently because the upper part of the LV is modeled as a single opening with circular topology (Suinesiaputra et al., 2018). To have a complete model, it is necessary to include at least the LA and ascending aorta (AA) to be able to model the mitral and aortic valves as a secondary result of a multichamber segmentation approach. The multichamber cardiac segmentation has traditionally been performed by active shape models (ASMs) and deformable methods to extract the anatomy and contractile function of the heart. A complete survey can be found in (Tavakoli and Amini, 2013). More recently, convolutional neural networks (CNNs) have also been proposed for multichamber cardiac segmentation because of their outstanding performance in segmentation problems Payer et al., 2018).
In addition, 4D computed tomographic angiography (4D-CTA) is a very useful tool to complement echocardiography and magnetic resonance imaging (MRI) in cardiac imaging (Maurer et al., 2012;McVeigh et al., 2018;Ohnesorge et al., 2000;Storz et al., 2017), providing 3D+t high-resolution scans covering a whole heartbeat cycle, making it especially recommended in the task of assessing the heart valves (Chheda et al., 2010). However, although cardiac structures can be best visualized in 4D-CTA relative to echocardiography and cardiac MRI, providing a superior depiction of its anatomy and dynamics, segmentation, and tracking of cardiac chambers and valves in 4D-CTA usually requires user interactions, as anatomical landmarks need to be manually located to preregister the shape model (Von Berg and Lorenz, 2005). Moreover, high sensitivity to initialization and subjective landmarks annotations in training may reduce the performance of automatic methods (Ecabert et al., 2008;Peters et al., 2010;Zheng et al., 2008). More recently, a semiautomated method based on the registration of a template mesh has been proposed for canine studies (Vigneault et al., 2018). Because of the arduous task of manual segmentation of cardiac valves, especially when complex models are involved, several works have been proposed to achieve this goal automatically despite its inherent difficulty. A patient-specific constrained shape model of both mitral and aortic valves was proposed by Ionasec et al. (Ionasec et al., 2010) for 4D-CTAs transesophageal echocardiographic images that started from a landmark model refined by a full physiological surface valve model constrained by the previously defined landmarks and also used steerable features (Zheng et al., 2008). This work was later extended to the location and study of the motion of the four heart valves (Grbic et al., 2012). The mitral valve alone has been studied in (Voigt et al., 2011) from transesophageal echocardiographic images, where the valvular model is achieved by combining discriminative learning and leaflet biomechanics. The dynamic segmentation of aortic valve, including the root wall and leaflets, has also been proposed (Auricchio et al., 2011) and uses a combination of elliptical shapes to model the leaflet-sinus complexes (Rankin et al., 2013).
These aforementioned methods require basic training that must be manually determined across the cardiac phases, but given the complex nature and variability of cardiac data, methods based on ASM or active appearance models have problems with the accurate 3D+t representation of size variations of cardiac chambers, as well as the limited spatiotemporal representation of the 4D behavior. In addition, models have limitations in the representation of pathological cases (Hoogendoorn et al., 2013). A database correctly labeled by experts is also necessary to create the model and should include pathological cases if the latter were also to be segmented. New proposals based on CNNs and deep learning need even larger labeled data sets to train the network, but these large data sets are not always available.
To overcome these limitations, graph-cut methods (Boykov et al., 2001) have been used to solve different segmentation problems in cardiac imaging. Graph-cuts were developed to address a wide range of problems in computer vision using graph theory. The main idea is to minimize an energy cost function by finding the minimum cut of a graph. For an in-depth understanding of the graph-cut segmentation approach, see (Boykov and Funka-Lea, 2006).
As mentioned above, graph-cuts have been used to segment the endocardial wall of the LV in MRI and echocardiographic 3D images (Bernier et al., 2017). If we focus on 4D data, several methods have been proposed: a coarse 4D-graph-cut, later refined using competitive region growing, was used for segmentation of cardiac chambers in dynamic 4D MRI data (Lombaert et al., 2011). Other 4D methods use a spatiotemporal watershed cut for LV myocardium segmentation (Cousty et al., 2010), as well as LV segmentation in MRI with smoothness and interslice constraints (Albà et al., 2014). Apart from cardiac imaging, 4Dgraph-cuts have also been proposed for the segmentation of brain structures in longitudinal MRI data (Wolz et al., 2010).
The main contribution of this work is the proposal of a processing pipeline to locate and track the displacement of the aortic and mitral valve annuli and associated planes from cardiac electrocardiography (ECG)-gated 4D-CTAs based on a 4D-graph-cut algorithm. Concretely the pipeline targets to locate the posterior segment of the mitral valve (D-shape plane) (Blanke et al., 2014) and the virtual ring defined by the aortic hinges (Piazza et al., 2008). The acquisition protocol of 4D-CTA ensures that the injected contrast is limited to the left chambers and aorta; therefore, a mesh can be built by isolating regions with contrast. A 4D graph is then constructed by encompassing these meshes and the temporal information between time frames. We can detect and track the mitral and aortic valvular annuli after segmenting LV, LA, and aorta using a unique 4D-graph-cut algorithm for the whole cardiac cycle. Our method aims to be fully automatic and not dependent on a large training database, implying that no ASM should be used in the final step (the refinement stage). However, because of the need for robust initialization, we use an ASM for the initial location of the cardiac chambers.
We have tested our processing pipeline with 15 patients comprising a total of 278 3D scans. To validate the usefulness of our proposal, functional derived parameters that are influenced by the correct valve position during the cardiac cycle are computed and compared with the results obtained from manually detected valves. In addition, the error in the valve location is assessed using manually selected landmarks.

Methods
In this section, we present a new method to locate and automatically estimate the displacement of mitral and aortic valve annuli in ECG-gated 4D-CTAs. Left chambers are also segmented as an auxiliary result. The proposed workflow follows the main steps summarized in Fig. 1. The following sections detail the different steps in the process.

GMM and binary mask
The contrast agent highest concentration is found in the oxygenated arterial blood when the ECG-gated images are acquired. Therefore, left chambers, aorta, pulmonary veins, and coronary arteries present high values of Hounsfield units (HU). After removing background and lungs regions with Hounsfield values (also referred as CT numbers) closer to air values (HU = −1000), the remaining voxels are grouped in clusters values associated with soft tissues and blood without contrast, oxygenated blood with contrast, and bone. To differentiate between tissues and regions with the contrast agent, the image is modeled as a GMM.
GMMs have been extensively used for the characterization of tissues in ultrasound imaging (Vegas-Sanchez-Ferrero et al., 2014). Recently, the GMMs have also been studied for the characterization of CT numbers in CT imaging and were proven to be excellent models to describe the statistical properties of CT numbers for different doses and reconstruction kernels (Vegas-Sánchez-Ferrero et al., 2017). The characterization with GMMs allows us to calculate the posterior probability of belonging to each component of the mixture model. In our context, we take advantage of the GMM characterization to calculate the probability distribution associated with the contrast agent instead of the original CT numbers. For brevity, in the rest of the proposed workflow, we will refer to this posterior probability as the gamma map.
The proposed workflow uses a triangle mesh of the surface of the left chambers (LV and LA) and AA obtained from the binary mask. This mask is calculated with an automatic thresholding algorithm of the gamma map (Otsu, 1979) that was used to separate highintensity structures (chambers and arteries with contrast agent and bone) from the rest of the tissues.
Some postprocessing is required to obtain the binary mask comprising the LV, LA, and AA: first, a morphological closing is performed to soften the surface to eliminate irregularities caused by image noise and complex topology caused by the LV trabeculae and papillary muscles; second, the largest segment containing the left chambers and the AA is selected, discarding other subregions that include bones and other structures with high values in the gamma map; third, a morphological opening removes the coronary arteries. Lastly, all possible holes in the remaining mask are filled to get a closed and topologically simple isosurface prior to the extraction of a triangle mesh.

ASM of the LV
To obtain a first automatic location of the LV, the blood pool of the LV is detected and delineated using a statistically based shape model of the endocardium of the LV, and this model is derived from the model of the whole human heart (Hoogendoorn et al., 2013).
We first locate the LV model by using Horn's quaternion-based method (Horn, 1987) to find the rotation, translation, and scaling that optimally map the vertices of the shape model to candidate locations extracted from the gamma map. The fitting is performed by the weighted least squares method where the candidate points are selected in the normal directions of the surface that show the maximum gradient of the gamma map (i.e., the most likely border).
The weights are also derived from the gamma map. For a more detailed explanation, we refer the reader to (Haak et al., 2015). The fitting method is performed iteratively until the estimated error reaches a given threshold. In the next step, the modes are projected using a weighted regularized ASM (Zhao et al., 2004) to estimate the LV pose.

Triangle mesh
In the following step, we extract the surface of the binary mask described in Section 2.1 parametrized as an unordered triangular mesh (Qianqian and Boas, 2009). The meshing algorithm implementation is described in (Rineau and Yvinec, 2007) based on the previous work of ( Boissonnat and Oudot, 2005), choosing 1.5 mm as the maximum radius of the Delaunay sphere. The mesh is corrected to remove existing duplicated, isolated, and nonmanifold nodes to obtain a single watertight triangle mesh without singularities, selfintersections, and degenerate elements (Attene, 2010;Attene et al., 2013). A mesh of approximately 2 mm edges is later obtained using a remeshing algorithm that iteratively cuts small edges and subdivides larger ones.
The resultant mesh may also contain remaining structures associated with pulmonary veins and coronary arteries with contrast. These structures, although small in volume, can contain a high percentage of triangles; therefore, these structures are removed before the main graph-cut algorithm to ensure robustness and efficiency. These structures are segmented using the shape diameter function (SDF) introduced in (Shapira et al., 2008), which estimates the neighborhood diameter of the mesh as an approximation of the medial axis transform. On a smooth surface, the diameter can be defined by the distance to the antipodal surface point along the inward-normal direction. To reduce the variance because of noisy and low-resolution surfaces, the SDF (Shapira et al., 2008) takes the mean value of intersection distances over directions inside a cone centered in the inward-normal direction. A clustering algorithm then segments the surface regions with small shape diameter values, e.g., coronary arteries.
A coherent point drift (CPD) algorithm for nonrigid point set registration (Myronenko and Song, 2010) was used to register pairs of meshes in consecutive frames, being the barycenter of the faces the points used in the CPD algorithm.

Graph-cut formulation
To solve the separation of the different structures in the sequence, a 4D (also referred as 3d +t) graph-cut approach is proposed, establishing both spatial and temporal edges to impose both spatial and temporal coherence by combining the information of all the frames. An undirected graph G = ⟨V,E⟩ is composed of a set of nodes V and a set of edges E , where each edge connects two nodes {p,q}∈E : p,q∈V. Therefore, the topology of a triangle mesh can be formulated as an undirected graph if the triangle faces are identified as nodes and the neighboring ones are connected by edges. In addition, the triangle meshes of all the time frames can be joined in a single undirected graph to achieve the 4D structure in such a way that a set of nodes V is composed of all triangle faces from all N frames: V ≡ ⋃ n = 1 N V n , where V n is the set of nodes (i.e., triangle faces) obtained from the time frame n.
We define the set of edges E by considering two types of connections: edges connecting two neighboring faces of the same triangle mesh associated to a time frame (E s ), and edges connecting two triangle faces belonging to two consecutive time frames that were registered by using the previously described CPD algorithm (E T ). For clarity, we call them spatial and temporal edges, respectively, and we will use superscript S for spatial terms and T for temporal terms. Decomposing by time frames, the set of nodes and edges G = ⟨V,E⟩ can be expressed as: wℎere E n S ≡ {{p, q} ∈ E : p, q ∈ V n } and E n, n + 1 T ≡ {{p, q} ∈ E : p ∈ V n , q ∈ V n + 1 } . (1) In the preceding expression, E n s is the subset of spatial edges connecting nodes of the time frame n, and E n, n + 1 T is the subset of temporal edges connecting nodes of consecutive time frames n and n+1.
The complete graph G = ⟨V,E⟩ encodes a 4D structure because it includes spatial edges associated with the 3D triangular meshes and temporal edges encoding the temporal relation between 3D meshes extracted from the different frames. Fig. 2 schematizes this 4D graph.
The segmentation of the graph G is formulated as a graph-cut algorithm (Boykov et al., 2001) to find a labeling f of some finite set that assigns a label f p to each p∈V by minimizing an energy function E(f) with two cost terms: The boundary term, denoted as E bound , accounts for the edges connecting the nodes in which f is not piecewise smooth, while the regional term E reg describes the disagreement between f and the observed data. These terms can be calculated as a summation of costs over edges and nodes: where B {p,q} (f p , f q ) is the boundary cost of assigning labels {f p , f q } to the edge {p,q} and R p (f q ) is the regional cost of assigning the label f p to the node p. The detailed sum of boundary costs, distinguishing between temporal and spatial edges, is expressed as: where the boundary cost terms B {p, q} S are for spatial edges and the boundary cost terms B {p, q} T are for temporal edges. Using expressions (4) and (3), the components of (2) can be expressed as: Expression (5) is general for graph-cut algorithms with boundary and regional costs defined over 4D graphs (i.e., graphs with spatial and temporal edges In the following sections, the spatial components B {p, q} curv (shape curvature cost) and B {p, q} gamma (gamma value cost), as well as the regional costs R p , are described. The temporal boundary term is constant and only dependent on the CPD algorithm output.

Shape curvature cost function term-
The curvature of the triangle mesh can be used to identify the valve locations, as shown in Fig. 3. Note that the valve location is usually characterized by an indentation of the endocardial surface. Therefore, the cost associated with curvature should be minimal in concave areas of the surface. Each term B {p, q} curv is calculated as follows: first, we estimate the curvature tensor in the surface, as described in (Alliez et al., 2003); then, we regularize spatially the resulting tensor field though a local average of the curvature tensor (elementwise) in a local neighborhood defined around each vertex; lastly, the lowest valued eigenvalue of the curvature tensor is scaled to the range [0,1] in the entire triangle mesh, and the value in each edge is then calculated as the mean value in its two vertices. The minimum curvature was chosen as it presented a greater contrast at the valve indentation areas (see supplementary Fig. S3).

Gamma value cost function term-
The information around the surface can provide useful information to identify the location of the valve. Therefore, we consider the surrounding samples of the gamma map.
To retain the advantages of the fast segmentation method provided by the triangle mesh strategy, we perform a fast sampling strategy of the gamma map around each edge, searching for minimal values along directions randomly distributed around the normal of the edge. This search aims at finding valve leaflets in open or closed positions, which have small values in the gamma map. An example of this procedure is shown in Fig. 4, and the detailed procedure is as follows:

a.
Given the outward-pointing normal vectors n p and n q of the faces of the triangle mesh associated with the edge {p,q}, the vector n {p,q} = −(n p + n q ) is normal to the edge and points toward the interior of the mesh.

b.
The central point of the edge {p,q}, calculated from the coordinates of its two associated vertices in the mesh, and the vector n {p,q} define the vertex axis of a cone with an aperture angle δ between 30 degrees and 60 degrees.

c.
N vectors r n are randomly calculated inside the cone to avoid undesirable patterns derived from uniform sampling when N is small. Then, the gamma map is sampled across N segments of length d following directions of the vectors from the central point of the edge.

d.
The average of the samples is calculated per segment, and the cost is set to the minimal average across edges per node.
2.4.3 Regional cost function term-To calculate the regional terms R p in expression (5), a seed belonging to each chamber (i.e., AA, LA, and LV) located away from the valve planes that are going to be determined is placed automatically given the ASM fit. The triangle faces selected as seeds and their local neighborhood contribute to a negligible regional cost when the correct chamber label is used, and the maximum regional cost is when the wrong label is assigned corresponding to one of the classes (sink or source). The regional cost for the rest of the triangle mesh is calculated using the geodesic distance from each seed. Geodesic distances are computed using the fast marching method (Sethian, 1996) over a triangle mesh (Kimmel and Sethian, 1998), where the metric is induced by the embedding of the 3D surface in ℝ 3 .
We strive to provide a robust method to select the chamber seeds, so we use the ASM to define the most likely location of the chambers and valves. Then, we locate the seeds according to the location and pose of the ASM. The adjusted mesh of the LV using the ASM model has an open aperture composed of a loop of points, as represented in Fig. 5. A plane is fitted to those points. Because the model is expected to be reasonably adjusted, this plane is close to the mitral valve, while AA and LA are supposed to fall above this plane, whereas the ventricular apex is in the opposite direction. Accordingly, the point of the triangle mesh with maximum distance to this plane below the plane is chosen as the ventricular apex of the triangle mesh, which is the seed of the LV.
The rigid transformation obtained in the estimation of the ASM model of the LV is also applied to place a separation plane between LA and AA, and this transformation was previously calculated in the ASM model (Fig. 5). The parts of the triangle mesh that have a large distance to that plane in the right direction and simultaneously are very distant from the ventricular apex have a high probability of belonging to AA. In addition, this part of the triangle mesh, because of the tubular shape of the AA, is intersected by the plane defined by the Z-slices of the CT with an elliptical well-defined shape and is therefore adjusted with an ellipse fitting to ensure that the AA is correctly detected. The central point of the fitted elliptical shape is selected as the AA seed.
Once LV and AA seeds are detected, the upper LV plane and LA-AA planes of the ASM model are also used to choose the seed of the LA as the point that maximizes the sum of the squares of distances from both planes in the right direction, excluding points where the geodesic distance from the AA seed is smaller than a threshold to avoid selecting another point in the AA as the LA seed. (Boykov et al., 2001;Kolmogorov and Zabih, 2004) is used to segment the 4D graph to minimize the energy expression (5). The graph-cut algorithm optimizes multilabel energies via the α-β swap (Boykov and Kolmogorov, 2004).

Graph-cut algorithm-A graph-cut algorithm
To simplify the problem of multilabel segmentation further, the triangle meshes in LV, AA, and LA are divided by two independent bilabel segmentations: first, an AA-non-AA graphcut algorithm is performed in which the LV and LA are merged together as a unique label. In this way, the graph-cut detects the aortic valve as the separation between segments. This is followed by a second stage where the subgraph of the AA can be removed, and an LV-LA bilevel graph-cut segmentation is performed to detect the mitral valve.

Extraction of valvular annuli location
As a result of the 4D-graph-cut, the surfaces of the left chambers are segmented. The mitral and aortic valvular annuli locations are identified in a secondary step and extracted as the vertices of the mesh connecting two faces belonging to two different chambers. Because this contour has a 3D shape, the valvular plane is approximated by fitting this contour to a plane.
In the case of the mitral valve, its characteristic saddle shape causes that a fitted plane does not properly represent the annulus location. Instead, the posterior annulus segment is better approximated by a plane, according to the D-shape model (Blanke et al., 2014). We use this segment and the fitted plane to determine the displacement of the mitral valve along the cardiac cycle and to extract derived parameters.
The posterior annulus segment is calculated from the obtained mitral and aortic complete annuli estimated locations, extracting the approximated D-shape of the mitral annulus by selecting the principal axis of the two annuli, and the two most distant points in the mitral annulus to this axis. Fig. 6 illustrates this procedure.

Experiments and results
The evaluation of the proposed methodology was conducted in 4D-CTA scans acquired from 11 subjects with normal LV function and 4 patients with low ventricular ejection fraction. The acquisition was performed during an inspiratory breath-hold, full R-wave to R-wave, covering a complete cardiac cycle. To reduce the heart rate, beta blockers were administered to those subjects with no contraindications. The average intensity values of the ascending aorta given in HU were used in the bolus-tracking method to trigger the acquisition. For this process, 60-100 mL iodinated contrast (Iopamidol, Isovue 370, Bracco) was injected at a rate of 4.5-5.0 mL/s, followed by a 50 mL saline flush.
The studies where reconstructed using the vendor-provided methods. CT scanners used, acquisition data and reconstruction parameters using vendor-provided methods are shown in Table 1. Detailed pixel size, field of view (FOV), and cardiac cycle duration of the studies are listed in Table 2. Study identifiers (IDs) are maintained in all the figures and tables of this work to facilitate the identification of the reported results. IDs #1 to #11 are identifiers of healthy subjects, while IDs #12 to #15 identify pathological cases.
The CT data from subjects with normal cardiac function were conducted in the Department of Cardiology of the National Heart, Lung, and Blood Institute under IRB approved protocols. The CT data from subjects with reduced LV function were conducted at the UCSD Thornton Hospital (La Jolla, CA) under IRB approved protocols.
The number of distributions of the GMM used for the estimation of the gamma map was initially fixed to two for healhy subjects, and was tuned to three in the case of patients with contrast in the right cavities. Fig. 7 shows the segmented surfaces of the left chamber and their intersections with a slice of the LVOT. This cardiac view allows us simultaneously to visualize the location of the mitral and aortic valves as the points where the color changes. The rendering of segmented chambers is also superimposed onto Fig. 7. The result shown corresponds to study ID #1 in end-systolic and end-diastolic frames. The color codes used are red for the LV, green for the LA, and blue for the AA.
The complete result of the 4D-graph-cut segmentation, covering all the time frames of the cardiac cycle is depicted in Fig. 8 for study ID #5 and in Fig. 9 for all studies covering the end-systolic and end-diastolic frames. The surface mesh segmentation is labeled with red, green, and blue corresponding to the LV, LA, and AA, respectively. The animations of all the studies have been submitted as supplementary material.
The main parameters of the 4D-graph-cut algorithm were optimized using a small cohort of 3 initial healthy subjects and a simulated annealing algorithm. The table S1 of supplementary material specifies the parameters used in the experiments. Average time to process each case is less than 20 minutes, using a not fully optimized MATLAB ® implementation running on an Intel(R) Core(TM) i7-6700 CPU, 3.40GHz.
To test the performance of the preprocessing steps and the role of the gamma map, we compare the ASM with the gamma map to the ASM on the original CT numbers for all the cases. As a result, the ASM fitting using the gamma map achieved 100% success in all frames, while trying with the original CT numbers, only 39 frames out of 278 were adjusted acceptably.

Evaluation of the location estimation of the mitral and aortic valves annuli
The error in the location of the mitral and aortic valve annuli is defined as the agreement between the manual annotations and our automated method. To define the reference location of the mitral valvular annulus, all acquisitions were inspected manually, and three anatomical point landmarks along the mitral annulus were selected in all frames and studies. Two landmarks were selected near the anterior lateral and posterior medial commissures of the mitral valvular annulus, while the remaining landmark was selected approximately in the middle of the posterior leaflet. These three manually selected anatomical landmarks define the mitral valvular plane, as represented in Fig. 10.
The mean Euclidean distance from the extracted segment of the posterior mitral annulus and the plane defined by the manually selected landmarks is considered in the reported errors. Fig. 10 represents the method of error estimation, as described previously, illustrated with the details of mitral annular valve localization.
Because the contours are constituted by line segments connecting points, we calculate the Euclidean distance between each point of the contour with respect to the plane obtained from manually selected landmarks, reporting the mean value of this set of distances. We define the mean value of this measurement over all time frames as the mean absolute error (MabsE) in the localization of the mitral valvular annulus, which is reported in Table 3. The errors along the time frames are indicated in the box plots of Fig. 11.
In contrast, for the manual location of the aortic valvular annulus, the basal attachments of the three aortic cusps were manually selected as landmarks, defining a virtual ring that is slightly below the anatomic atrioventricular junction but clearly detectable in the image. The error of the location method, in this case, was calculated as the mean Euclidean distance between the manually selected landmarks and the corresponding closest points in the contour extracted from the segmented mesh. In this case, we did not use the distance between plane and contour, as reported for the mitral valvular annulus location because the aortic valvular annulus is far from being an ideal plane because of the anatomical shape of the cusps. Aortic valve localization MabsE is reported in Table 3, and agreement distances along the time frames are reported in Fig. 11. The total mean value of MabsE was 2.52±1.06 mm for the mitral annulus and 2.00±0.69 mm for the aortic valve annulus. Data are expressed as the mean values ± standard deviation.
In order to assess interobserver variability we asked a second observer to delineate the posterior mitral annulus in the frames #8 and #16 (approximately corresponding to endsystolic and end-diastolic moments) in studies IDs #1 to #11. The delineation was performed between the posteromedial commissure and the anterolateral commissure, covering the posterior leaflet the mitral valve. This part can be approximated by a plane using the Dshaped model as previously reported (Blanke et al., 2014). The agreement between plane defined by the three landmark points marked by the first observer and this delineation resulted in a final distance of 1.62±0.67 mm. For comparison, our reported agreement in the mitral valve is of 2.72±1.30 mm in the same 11 cases and two frames. The paired mean difference between the reported errors and the interobserver variability is 1.11 mm [95%CI: [0.57, 1.86] mm]. The P value of the two-sided permutation t-test is 0.0022. Therefore, these results assess a low (1.11mm) and significant difference with respect to the interobserver variability.
On the other hand, the overall quantitative effect of the temporal edges has been evaluated by calculating the 3D-graph-cuts segmentations independently for each frame using the same terms and parameters. The comparison respect to the 4D result is shown in Fig. S1 and Fig. S2 of supplementary material. It is emphasized that the improvement is mainly noted in the decrease of outliers thanks to the continuity obtained by the temporary edges. In the mitral annulus, the 10th percentile of worst frames have a mean agreement distance of 5.8 mm (worst frame, 7.2 mm) for the 4D-graph-cut method, while in the 3D case, the mean agreement distance of the 10th percentile of worst frames is 7.4 mm (worst frame, 13.7 mm). In the case of the aortic annulus, the comparison is more beneficial for the 4D method: the mean agreement distance for 10% percentile of worst frames is 4.4 mm (4D-graph-cut) and 11.0 mm (3D-graph-cut).

Excursion in mitral and aortic valves
The mitral annular plane systolic excursion (MAPSE) represents the total displacement of the mitral annular plane toward the apex between systole and diastole and assesses the size variation of the LV in the direction of the apex (Hu et al., 2013). Given the wide use of this descriptor, we performed an automatic measurement experiment to test the performance of our method. Apart from the displacement at systolic frames, the time evolution of the displacement between frames has also been quantified. We denote this dynamic mitral annular excursion as MAE. We have compared the measured MAE from the manually selected point landmarks with respect to the estimated MAE calculated from the mitral annulus extracted using the proposed 4D-graph-cut method.
Because of the complex movements of the mitral valve annulus and angular variations of MAE, we follow the approximation of measuring the linear excursion projected on the axis normal to the mean mitral valvular plane (considering the whole cardiac cycle) and consider the relative local coordinate along such axis.
The mean value of the projected excursion of the manual landmarks is considered as the measured MAE at a certain time frame, while the mean value of the projected estimated excursion of the anatomically equivalent points in the estimated mitral annulus is considered as the estimated MAE. Once the excursion estimations are obtained for all time frames, the time frame corresponding with the lowest mean value is taken as reference (diastolic frame), and the rest of the frames imply a negative displacement in mm. Fig. 12 shows this dynamic excursion, frame-by-frame, for all cases. The systolic excursion corresponds to the total mitral excursion (MAPSE) and is reported in Table 4.
A second estimate of MAE and its total displacement (MAPSE) has been computed using the ASM model used in the initialization, approximating the upper ring of the modeled LV as the mitral annulus. The results are also shown in Table 4 and Fig. 12.
The aortic annulus excursion (AAE) is calculated in the same manner as MAE for the mitral annulus, i.e., measured along the normal axis to the mean aortic plane. In Fig. 13, the frameby-frame AAE using manually selected landmarks is compared with AAE extracted from automatically estimated aortic annulus. The maximum displacement of AAE along the cardiac cycle is denoted as total AAE (TAAE) and is reported in Table 4.
Estimated total MAE differences between systolic and diastolic phases (i.e., MAPSE) ranged from −7.25 mm for study ID #3 to −15.1 mm for study ID #8. Comparing with the results from manually marked landmarks, the mean absolute error of the total mitral excursion is 1.58±1.63 mm, with only one case with the error above 5 mm. As a comparison, the error using the ASM fitting was 5.55±6.37 mm, with two cases of patients reporting errors around 20 mm. TAAE measurements are of lower range than MAPSE, as expected: varying from −5.3 mm for study ID #10 to −10.0 mm for study ID #3. The mean absolute error is 0.81±0.57 mm.
We analyze the differences between the MAE estimates for the different methods depicted in Fig. 12 by considering a mixed effects model described in the following equation: MAE i, j, k = Subject i + (Metℎod) k + (Metℎod ⋅ F rame) j(i), k + ε j(i), k where Subject accounts for the random effect given by the subject sampling (i=1,…, 15), Method is the fixed effect given by the different approaches (proposed method, manual landmarks, and ASM), and Frame is the temporal phase (where j(i) is the frame index for the i-th subject whose range is described in Table 2). The fixed effects test showed no statistical significance in the difference between methods (p-value = 0.1278). Similarly, we performed the same mixed effects model analysis for the aortic AAE estimates shown in Fig. 13. The fixed effect test showed no statistically significant differences between the proposed method and the measurements using manual landmarks (p-value = 0.6151). In

Global longitudinal strain of the LV
The estimated mitral valve annulus locations along the cardiac cycle have been used to report measurements of another relevant image-based descriptor: the global longitudinal strain of the LV (LVGLS) (Reisner et al., 2004), which is defined here as: where L ED and L ES represent the length of the LV over the myocardium in the end-diastolic frame and end-systolic frames, respectively. This definition is an approximation of the LVGLS typically defined as the mean value of the local longitudinal strain over the myocardium. L ED and L ES are approximated as the distance between the mitral valve and the ventricular apex.
In this experiment, we evaluated the validity of the proposed method to compute this descriptor. In our 4D-CTA, the ventricular apex is approximated as the most distant point of the ventricular cavity to the mitral annulus. The mitral valve location considered for estimating L ED and L ES is the center of mass of the estimated posterior segment of the mitral annulus.
The ventricular apex is calculated using the gamma map, highlighted by the administered contrast inside the ventricular cavity, which makes detecting the most distant point in the vicinity of the ventricle mesh a robust process. Therefore, manual marking was not used.
We denote the distance between the valve location and the ventricular apex as L MA . Therefore, in the end-systolic frame, it follows that L MA (ES) = L ES , while in the enddiastolic frame, it follows that L MA (ED) = L ED . We can extend the LVGLS to all time frames as the longitudinal strain (LS) as follows: In Fig. 14, the evolution of LS along the cardiac cycle is depicted, comparing the result obtained when the mitral location is estimated by the proposed segmentation using 4Dgraph-cuts with respect to the mitral valve location obtained using manually selected landmarks.
The mixed effects analysis performed for MAE and AAE was also employed for the LS measure of Fig. 14. As in the previous analyses, the fixed effect test showed no statistically significant differences between the proposed method and the manual landmarks (p-value = 0.8903).
Based on the maximum and minimum values of the evolution of L MA in Fig. 14, LVGLS is reported for both the automatic estimated method and for the manual measurements in Table  5. The mean absolute error is 1.57±1.24%.

Left ventricular ejection fraction
The left ventricular ejection fraction (LVEF) measured as LVEF = 1−V ES / V ED , where V ED corresponds to the maximal ventricular cavity volume, achieved at the expansion phase in the end-diastolic frame, and V ES is the minimal cavity volume reached at the contraction phase at end-systole. The LVEF is usually measured based on the LV cavity and a fixed selected plane at the base, but this value can be refined by including the total ventricular volume the mitral and aortic valvular annuli.

Discussion and conclusion
In this work, we propose an automatic method based on 4D-graph-cuts to determine the displacement of the mitral and aortic valvular annuli in 4D-CTA images. This method has been tested on a set of 15 dynamic studies, comprising a total of 278 3D images, without substantial errors. The proposed pipeline has been designed to be fully automatic and not dependent on a large training database. Initial alignment and scaling of the ASM model of the LV were employed in one reference study, while the remaining patients were correctly segmented using the same initialization regardless of substantial differences in the heart sizes and FOV of the images. The robustness in the initialization provided by the ASM shows that the method is sufficient without needing manual alignment.
The displacement is successfully completed throughout the time frames even when only partial LA and a small portion of the AA is found within the FOV, as shown in study ID #2 and in Fig. 9. The proposed procedure succeeded in both healthy subjects and pathological cases.
The robustness of the ASM fitting to each frame was ensured with the use of the gamma map instead of the original CT values present in the data because there is a clear difference between the probabilities of the gamma map with high values (close to one) in the blood pool with injected contrast, and the probabilities of the gamma map with low values (close to zero) in the surrounding tissues. Therefore, a clear gradient can be detected around the cavity of the LV filled with contrasted blood during the ASM fitting process. If the original CT numbers were used instead of the gamma map, there would also be clear additional gradients, e.g., between the epicardium and lungs, that would hinder the ASM fitting. The ASM model of the LV uses improved weighted optimization algorithms in both registration and projection steps. This way, points with negligible gradient values along the scanned line segments have close-to-zero weights and are not really involved in the ASM fitting process. In contrast, the most relevant points have higher weights and reliably guide the ASM fitting.
Our experiments on pathological cases demonstrated that, even though the ASM was built from a healthy population and its performance might be sub-optimal for these cases (e.g., Fig. 12, study IDs #13, #14 and #15), it was enough to detect the initial location of the ventricle and therefore the workflow performed properly in the unhealthy patients.
The role of the gamma map extracted from the GMM is not only confined to the fitting of the ASM model of the LV. It also allows us to create a binary mask containing the contrast regions (e.g., right chambers, aorta, and pulmonary veins). The unprocessed mask has an irregular surface and complex topology because of the presence of image noise and small structures, such as papillary muscles, LV trabeculae, or coronary arteries, and needs to be softened by morphological closing and opening operations to obtain a surface without singularities, degenerate elements, or topological artifacts. In contrast, excessively smoothed surface meshes are counterproductive because information on the surface curvature would have been lost and unable to adjust faithfully to the edges of the gamma map.
The proposed method could have benefited from a multiscale scheme in which the first coarser segmentation would be refined in an iterative process, obtaining a series of increasingly finer meshes, but the total processing time is increased. An alternative use of the increase in processing time is to use the multiscale approach only in a region around the valvular annuli segmented in coarser scales. We did not implement this multiscale approach because our dataset gives accurate estimates applying just one scale. Nevertheless, a multiscale implementation will be considered in future work.
The surface-based weighting factors of the cost function in the 4D-graph-cut, denoted as B {p, q} curv , were designed to include information about local curvature of the surface mesh because the valvular annuli can, in some cases, be identified by concave regions, as shown in the right picture of Fig. 3. This fact is particularly true for the aortic valve because of the presence of clearly visible semilunar cusps, although they are less prominent in some patients (see, for instance, study ID #8 in Fig. 9. In other cases, the entire perimeter of the mitral valve cannot be recognized by this characteristic concave zone. Moreover, other structures, especially the ventricular trabecula, produce concavities. Therefore, a graph-cut based solely on curvature factors would not be robust enough for our purpose. To overcome this problem, we consider the information related to the gray level of the image. Similar to the ASM fitting, instead of using the original CT numbers of the image, we use the gamma map and the sampling scheme in conic regions, as shown in Fig. 4. These gray-level-based weighting factors B {p, q} gamma significantly change during the cardiac cycle. In the systolic phase, when the mitral valve is closed and falls inside the sampling cone, The CPD registration between time frames, which adds new temporal dynamic edges to create a 4D graph, helps to localize the valvular annuli in the cardiac phases correctly where there are no low weights in the curvature components B {p, q} curv or in the gamma components B {p, q} gamma of the cost function. In those cases, the vertices that are difficult to segment are connected to anatomically equivalent points in other frames in which these characteristics appear to present temporal consistency. This triple contribution to the weights of the cost function (curvature, gamma map value, and temporal registration) increases the robustness of the 4D-graph-cut algorithm.
We have compared our results against manual measurements to report the error. However, this was a difficult task: 1668 points had to be manually selected on the 4D-CTAs (three markers per valve in 13 to 20 time frames for 15 studies). Completely delineating the valves to measure the localization error along its entire contour would have been even more arduous, although it would have made possible using a saddle-shaped mitral annular assessment and reduce the reported errors. Without this complete delineation, the error can only be measured at the manually marked points on the valve, as was done for the aortic valve in Table 3. However, for the case of the mitral valve, this measurement based on only three points could lead to apparently good estimates in cases in which the resultant mitral annulus segmentations differ substantially from the real mitral annulus (saddle-shaped) approximated by a plane if it was coincident in the three marked points but not in the rest of the valves. Therefore, the error was reported based on the mitral annulus defined by the intersection of the mesh and a plane derived from the three manually marked points with respect to the resultant annulus from the proposed pipeline. It is important to note that this procedure has some limitations because of the manual marking of the reference points and the mentioned saddle-shape of the mitral annulus, and these limitations both imply an overestimation of the error.
Regarding the parameters directly related with mitral excursions, MAE and AAE, we used for simplicity the typical approach that considers the movement of the valve as a single plane, although mitral annulus motion can be broken down into three components (Silbiger, 2012): translation toward and away from the LV apex, folding across the intercommissural axis, and circumferential contraction.  (Yingchoncharoen et al., 2013). The obtained LVGLS values in our small cohort of patients were substantially smaller: −9.7±2.9%, which is associated with higher risk of heart failure (Biering-Sørensen et al., 2017).
In conclusion, we have proposed a robust and fully automatic pipeline for the estimation of the location and displacement of the aortic and mitral valve annuli in high-resolution 4D-CTA. The method includes the calculation of a GMM to obtain a binary mask that is then processed to obtain a triangulated mesh for each time frame. Then, a robust initialization using an ASM model of the LV determines a priori locations, and a CPD registration algorithm connects the meshes of all time frames to yield a unique graph on which graphcuts are calculated to obtain a temporally consistent segmentation of the left heart cavities that identify the mitral and aortic annuli. The method has been evaluated on a set of 11 healthy subjects and 4 patients, and our results have been compared against the manual identification of anatomical landmarks at the valve annuli. The reported agreement between the manual annotations and our method (MabsE) was 2.52±1.06 mm for the mitral annulus and 2.00±0.69 mm for the aortic valve annulus. We have also shown the impact of the valve plane detections on derived functional parameters (LVVF, LVGLS, MAE, and AAE) that are evaluated throughout the cardiac cycle. Valve excursions and the LVGLS were compared with manually detected valves, reporting mean errors of 1.58±1.63 mm for MAPSE, 0.81±0.57 mm for TAAE, and 1.57±1.24% for the LVGLS. The MAE results were also compared with results obtained from the automatic ASM fitting: our method reported a mean error below 1.6 mm, while the ASM mean error was larger than 5 mm. These results suggest that we have achieved a robust method of localization and determination of the displacements of the mitral and aortic valve to evaluate derived functional parameters without the need for a large training database. Its use in both pathological and healthy data is therefore direct without the requirement for extra training.

Supplementary Material
Refer to Web version on PubMed Central for supplementary material. Basic steps of the proposed workflow: A) a bimodal tissue vs. blood gamma mixture model (GMM) is estimated from a 4D-CTA to differentiate tissues vs. contrasted blood in the cardiac region, followed by the segmentation, for every time frame, of a binary mask comprising the blood pool of the left chambers, with high values of intensity because of the contrast; B.1) an active shape model (ASM) is used to fit the left ventricle (LV) and atriumaorta separation plane to ensure a robust initialization of the graph-cut; B.2) triangular meshes of the binary mask (gray color) are calculated for all the time frames in the sequence (shown for one single frame) that are processed to remove structures not belonging to the left chambers (red color); C) a surface-based temporal registration of the meshes allows the construction of a 4D graph containing spatial edges from the meshes and temporal edges from the registration. Seeds for the graph-cut are shown as color points; D) a final 4D-graphcut segmentation is performed, separating ascending aorta, left atrium, and LV, to obtain the mitral and aortic valve annuli in all frames and associated planes (shown for one single frame). (a) Active shape model (ASM) of the left ventricle (LV) depicted as a red triangle mesh, left atrium (LA) as a green mesh, and ascending aorta (AA) as a blue mesh, also representing the fitted plane to the upper LV aperture (violet plane) and separation plane between LA and AA (brown plane); (b) when the ASM model of the LV is fitted to the study (gray mesh), both planes are useful when selecting the initial seeds of the LV (red point) LA (green point), and AA (blue point).  4D-graph-cut segmentation of surface mesh, covering end-diastolic and end-systolic time frames of the cardiac cycle for all studies. The left ventricle is represented in red, the left atrium is represented in green, and the ascending aorta is represented in blue. Example of posterior mitral annulus localization and method of error estimation. Three manually selected landmarks (blue dots) are used to define a mitral plane (translucent plane in the image on the right) which intersects with the triangular mesh of the left chambers forming a contour, drawn as a blue polyline in the image on the right. The contour obtained by the proposed method is depicted as a red polyline in the figure. The performance in each frame is measured by the mean Euclidean distance between this red contour and the blue plane. Left: agreement distances in the determination of the location of the mitral valve annulus, measured with respect to the valvular plane obtained by three manually-selected landmarks. The distance is calculated for each frame and averaged along the sequence. Right: agreement distances in the localization of the aortic valve annulus. The distance is calculated as the mean Euclidean distance between the manually-selected landmarks and the corresponding closest points in the contour extracted from the segmented mesh. Mitral annulus excursion (MAE) along the normal axis to the mean mitral plane, comparing estimated automatic results using the 4D-graph-cut (blue lines), the active shape model (ASM) depicted in red lines, and the manual landmark-based result (dotted gray lines) for all the cases under study. Aortic annulus excursion (AAE) along the normal axis to the mean aortic valvular plane, comparing estimated automatic results (blue lines) using the 4D-graph-cut and the manual landmark-based result (dotted gray lines) for all the cases under study.  Evolution of left ventricle volume fraction (LVVF) along the cardiac cycle, as a percentage of maximal volume with respect to the end-diastolic frame. The blue lines represent the calculated the ventricular volumes from the automatically derived mitral and aortic valve planes along the cardiac cycle, while the red lines consider the estimated upper ventricular plane of the active shape model (ASM).   Table 3 Mean absolute error (MabsE) and standard deviation (SD) in the localization of the mitral and aortic valve annuli. Mitral valve: for every frame, the estimated error is calculated as the mean Euclidean distance between all points in the estimated contour and the contour determined by the intersection between the plane defined by three manual markers and the mesh. Aortic valve: the estimated error by frame is calculated as the mean Euclidean distance between three manually selected landmarks and the corresponding closest points in the estimated contour. For both valves, the MabsE is then calculated as the mean error distance for all frames. The last row contains the overall values for all 15 studies.  Systolic excursion of mitral and aortic annuli: results from the manually selected landmarks, from 4D-graphcut (GC) method and from the ASM in the case of the mitral valve. The absolute error is the unsigned difference between the manual values and the automated values. The total excursion is obtained along the normal axis to the mean annular averaging the whole cardiac cycle.