Cardiac motion tracking using a deformable 2D-mesh modeling

We present a framework for cardiac motion recovery using the adjustment of an electromechanical model of the heart to cine Magnetic Resonance Images (MRI). This approach is based on a constrained minimisation of an energy coupling the model and the data. Our method can be seen as a data assimilation of a dynamic system that allows us to weight appropriately the confidence in the model and the confidence in the data. After a short overview of the electromechanical model of the ventricles, we describe the processing of cine MR images and the methodology for motion recovery. Then, we compare this method to the methodology used in data assimilation. Presented results on motion recovery from given cine-MRI are very promising. In particular, we show that our coupling approach allows us to recover some tangential component of the ventricles motion which cannot be obtained from classical geometrical tracking approaches due to the aperture problem.


Introduction
The modelling of the heart's electromechanical activity is an active research area [5,9,1,13,4]. The simulation of the heart has received growing attention due to the importance of cardiovascular diseases in industrialised nations and to the high complexity of the cardiac function.
In order to help the clinical practice of cardiologists, it is important however that those models not only describe with some degree of realism the cardiac function but also be patient specific. Creating such personalised cardiac models implies that the anatomy of the patient is taken into account but also that the model parameters are tailored such that the simulated cardiac motion matches well with the observed cardiac motion. This represents a great challenge due to the intrinsic physiological complexity of the underlying phenomena which combine tissue mechanics, fluid dynamics, electrophysiology, energetic metabolism and cardiovascular regulation. Also only partial information can be derived from clinical data for a specific patient making the parameter estimation an ill-posed problem.
The objective of this paper is to propose a methodology that aims at creating personalised electromechanical model of the heart from cine MR images. Previous work [6,10,8,12] on the adjustment of a geometrical model of the heart on time series of medical images are mainly based on the concept of deformable models. In such a framework, a surface or volumetric mesh is fitted to the apparent boundaries of the heart by minimising the sum of two energies: an image term and a regularising or internal term. In such approaches, the model can be considered as a static system evolving under the minimization of an energy.
Conversely, electromechanical models of the heart are dynamic systems that evolve even in the absence of any image term. Adjusting such dynamic systems to time series of data (a method also known as "data assimilation") is fundamentally different from adjusting a static system since the parameters of the dynamic system are additional degrees of freedom that should be estimated. In the medical imaging community, P.C Shi and his group introduced data assimilation techniques by integrating cardiac models and Kalman filters for state and parameter estimation, see for instance [16] and [18]. However, such techniques, such as extended or unscented Kalman filtering, are often limited by the curse of dimensionality since they involve full covariance matrices whose size are equal to the square of the number of state variables augmented by the number of parameters to estimate. In the case of clinical applications, as cardiac electromechanical models are already complex dynamic systems with changing boundary conditions (cardiac phases), having a computationnally efficient estimation method is crucial.
In this paper, we propose an efficient method to estimate the state (i.e. the position and velocity) of an electromechanical model from cine MR images which is inspired from the deformable model framework used in medical image analysis. The goal of this paper is to show the formal equivalence between this approach and a filtering method introduced by Moireau et al. [7] used in data assimilation, which is different from Kalman-like filters such as the one used in [18]. The filtering approach proposed in [7] does not involve any matrix inversion (except the mass matrix which is a diagonal constant matrix), so that it allows much faster computations: the motion of a whole cardiac cycle on a mesh with 50 000 tetrahedral elements is estimated in about 10 minutes on a regular PC. This increases largely its potential future clinical application. The theoretical efficiency of this filter for mechanical systems has been demonstrated in [7]. The theoretical equivalence between the deformable model approach proposed here and this filtering approach leads to a better understanding of the trade-off between the electromechanical model and the image data.
We assume in this paper that model parameters are well known, in order to focus only on state estimation. Some preliminary results on parameter estimation are presented in conclusion, but this is not the goal of this paper. The proposed approach is first validated on synthetic time series of images and then applied to clinical cine MR images of a human heart.

Electromechanical model
We consider in this paper a fairly reduced electromechanical model since we want the complexity of the model to match the relatively sparse measures available from imaging data. Furthermore, this coarse level of modeling allows us to simulate a whole cardiac cycle on a mesh with 50 000 tetrahedral elements in about 5 minutes on a regular PC. Of course, the heart is a nonlinear material undergoing large strain. Thus, the assumptions of our simplified model are not realistic, but the global behavior of the heart is well represented. Furthermore, the limited computational time makes the estimation of the mechanical state and parameters tractable and allows us to test the behaviour of the model on series of heart beats.

Anatomy Description
The two ventricles are represented as a tetrahedral volumetric mesh including some anatomical information such as the myocardium geometry, the definition of some clinical anatomical regions (the American Heart Association regions), and the local orientation of fibres. We can build such a mesh from MR images, as explained below in section 3.1. The local fibre orientation can be either created from basic anatomical assumptions (elevation angle across the wall) or extracted from Diffusion Tensor MRI (DT-MRI) [11].
Several electrophysiological models have been proposed in the literature. Due to its efficiency, we use an Eikonal approach for the electrophysiology propagation, with a volumetric implementation of the algorithm described in [15]. The depolarisation time t d of the electrical wave for a given vertex of the volumetric mesh is computed by solving the anisotropic Eikonal equation v 2 (∇t T d D∇t d ) = 1, where v is the local conduction velocity parameter and D is the tensor defining the conduction anisotropy. In the fibre coordinates, D = diag (1, ρ, ρ), where ρ is the conduction anisotropy ratio between longitudinal and transverse directions. An anisotropic multi-front fast marching algorithm was developed in order to solve this model very efficiently.

Simulation of the myocardium contraction
The biomechanical model presented here is derived from a multi-scale modelling of the myocardium detailed in [2]. The mechanical model is composed of two elements, as shown on Fig. 1.a. The former is a parallel element which represents the passive properties of the tissue. This parallel element is anisotropic linear visco-elastic. The second element is an active contractile element controlled by the electrophysiology. More precisely, when the action potential is higher than a given threshold (i.e. when we reach the depolarisation time t d ), some calcium stored in the sarcoplasmic reticulum inside the cardiac cells is used for the ATP hydrolysis which provides energy to the molecular motors in the sarcomeres, generating the contraction of the fibre. The duration of this depolarisation is the action potential duration (APD). The electrical command u is then set to a constant k AT P which represents the rate of the hydrolysis of the ATP. After contraction, during the repolarisation, calcium moves back into the sarcoplasmic reticulum and this calcium decrease allows the relaxation of the muscle. The electrical command u is then set to another constant −k RS which represents the activity of the sarcoplasmic reticulum.
Thus, the contractile element is controlled by its corresponding command u through the differential equation:σ C + |u|σ C = |u| + σ 0 where σ C is the strength of the contraction, and σ 0 the maximum contraction. Then, with its associated command u described above, the strength of the contraction for each tetrahedron element is : where t r = t d + APD is the repolarisation time and HP the heart period. The command u and the intensity of the resulting contraction are represented on Fig. 1.b. Then, the active contractile element creates a stress tensor σ C f ⊗ f where f is the 3D fibre orientation and ⊗ the dyadic product. For each vertex of each element, this results in a 3D force vector f C = 1 4 R S (σ C f ⊗ f ) ndS with n the surface normal and S the element surface.
Finally, we represent the simplified dynamic law by a stiffness matrix K for the transverse anisotropic elastic part (parallel element), a diagonal mass matrix M, and a damping matrix C for the internal viscosity part, which is the Rayleigh damping matrix C = αM + βK, the contraction force vector F C created by the contractile elements, a force vector F P corresponding to the pressure forces in the ventricles and a force vector F B corresponding to other boundary conditions. The resulting law of motion is:  Let X = (Y,Ẏ ) T . Then, X is the state vector of the following dynamical system: where X 0 is the initial state vector, θ is the set of parameters of the model such as maximum contractility for example and where A (which depends of some parameters too) and R are defined by: We simulate the four cardiac phases (filling, isovolumetric contraction, ejection and isovolumetric relaxation) as described in [14]. The arterial pressures were computed using a 3-element Windkessel model described in [17].

Mesh Creation
4D (3D + t) cine MRI provides time series of high resolution images of the heart that describe in part or in total one (averaged) cardiac cycle. A cine-MRI typically consists in a sequence of 15 to 20 3D images for one cycle. The high intensity contrast between myocardium and ventricular blood pool allows a rough segmentation of the blood pools based on the combination of thresholding and connected component extraction. This segmentation is only used to demonstrate the possibilities of the method, a discussion on the various segmentation methods is out of the scope of this article. Fig. 2.c presents these two connected components for one image of the cardiac cycle. We need to build a computational mesh of the myocardium,

Model Initialisation
adjusted to the MRI image corresponding to the beginning of our simulation cycle. The first instant of our simulation cycle is the mid-diastole which corresponds to an instant when the ventricles are almost filled, just before the atrial contraction (P wave). We select for this the mid-diastole image, using the volume curves, detailed in the next paragraph. Then, the epicardium and left and right ventricles endocardia were delineated on this image using an interactive tool. These delineations generate three binary masks of the epicardium and the endocardia which are combined to obtain the binary mask of the myocardium used to create the mesh. This is done with isosurface extraction and tetrahedral mesh generation, using the INRIA software GHS3D (http://www.simulog.fr/mesh/gener2.htm).
We also need the local fibre orientation for this mesh. We generate synthetic fibre by linearly interpolating the elevation angle between the fibre and the short axis plane, from 80 o on the endocardium to −80 o on the epicardium. Fig. 2.b represents the obtained anatomical mesh with its synthetic fibre directions.  Electrical Model: As cardiac MRI is ECG-gated, we know the heart rate (here the heart period is 0.8 s) and the acquisition times of the 3D images related to the R-wave instant. This allows a first synchronisation between the image sequence and the simulation cycle. As the electrical information is not fully available, we need to extract additional information from the images. Due to the limited field of view, we only see part of the right ventricle in the MR images. Futhermore, the right ventricle blood pool has a grey level which varies along the cardiac cycle in cine MR images, thus thresholding is not reliable. Finally the trabeculae make the right ventricle segmentation difficult. For all these reasons, we have an important difference in volume between the two ventricles, as shown in Fig. 3. A more advanced segmentation method could overcome most of these difficulties, but this is out of the scope of this article. As our action potential propagation model only needs as inputs the time of the initialisation of the electrical wave and the action potential duration for each element, we extract average values from the volume curves. On these, one can observe the times of the beginning of the atrial contraction (P wave), of the ventricular contraction (R wave), and of the ventricular relaxation (T wave) independently for each ventricle (see Fig. 3). These times were set respectively to 0.0827 s, 0.125 s, 0.425 s. Then, we set the average value of the APD to the difference between the times of the beginning of the ventricular contraction and relaxation. Thus, for each vertex, APD is equal to 300 ms.

Model Initialisation
Mechanical Model: The passive mechanical parameters used are taken from the literature [16]. For the active component, we can use the volume curves to compute the ejection fraction, which is closely related to these parameters, in order to initialise it. However, due to the possible error on the right ventricle volume, we use only the left ventricle volume curve to calibrate the global contractility (the maximum contractility σ 0 constant for all the volumetric mesh) in order to obtain the same ejection fraction as the one computed from the left ventricle volume curve. For our data, σ 0 was set to 0.073 MPa/mm 2 . The rest position of the mechanical model is defined as the mid-diastole mesh created.

Coupling Model and Data: Methodology
In this section, we describe a method for coupling a dynamic system, the electromechanical model of the heart, and motion information from cine MRI. We start by discussing the choice of a metric to compare the simulated and observed motion and then describe formally the problem at hand: having a dynamic system that matches the available observations. Finally we show that motion tracking following a deformable model approach is equivalent to a data assimilation formulation where the error is minimised. This data assimilation formulation is directly inspired from the methodology of [7].

Metrics for comparing simulated and observed cardiac motion
Our objective is to minimize the discrepancy between the simulated cardiac motion and the actual one. One of the major difficulty is that in cine-MRI (which is the main dynamic modality in clinical routine MRI), only the apparent motion is visible. We see how the boundary moves, but we loose information on the tangential motion, which is important in the heart. We need to provide a metric to compare the model and the data taking this into account.
Since at each image instant the binary segmentation of the right and left blood pools are available it is reasonable to define the metric as the distance of the model endocardial surfaces to the blood pool surfaces as they should ideally match. Thus, for each point Y i of one endocardium surface of the mesh, we find the nearest point Y c i on the corresponding surface extracted in the MR image. Ideally, we want the distanced i between Y i and Y c i to be zero. This approach is illustrated in Fig. 4 (in green) in which n c is the normal to the blood pool surface at the point Y c i . Howewer, the distance maps must be either precomputed (storage costs) or computed during the estimation (computationnal costs). Thus, in this paper, we propose to use the reverse metrics: the distance of the blood pool voxels to the mesh vertices as shown in Fig. 4  Data interpolation: Due to limited temporal resolution, only a few MR images are available for a cardiac cycle. The time step used in the estimation is far smaller than the period between two MR images and we need information at each time step. Rather than interpolating the MR images, which would blur the contours, we prefer to interpolate the image forces described in the previous paragraph and computed at the previous and next images at each time step (see [14] for details).

Deformable Model Approach
This approach is based on segmentation by deformable models in which we minimise the sum of the energy of the dynamic system representing the heart and the energy corresponding to images forces, which are computed from contour images with distance maps for example. The introduction of the model in the minimised energy allows us to recover some movement which cannot be obtained from classical geometrical tracking approaches. Of course, the image forces have no physiological meaning, but if we couple the model and the data and if we estimate the model parameters (which is the next step of this work), the motion generated by the model should converge to the one observed in the images. Thus, the intensity of image forces should decrease along the estimation and the estimated motion should be more and more physiological.
The definition of image forces are consistent with the metrics chosen in the previous section. Namely, for each mesh point Y i , we seek the closest point Y img i along the normal direction N i of the mesh at Y i . Since the blood pool surfaces are roughly segmented as binary images, we compute Y img i as the intersection of the normal line at Y i with the isosurface I(x, y, z) = 127.5 for binary masks set to I = 255. This intersection can be computed fairly efficiently and with a subvoxel accuracy. More complex image forces involving intensity profiles, image blocks or textures could be used instead as shown in [3]. Here, we minimise the following energy:Ẽ where N i is the normal of the endocardium surface at the point Y i , m is the number of points of the endocardium surfaces (the points Y i are indexed from 1 to m for more simplicity) and γ i is the confidence in the measure Y img i . When we differentiate this energy with respect to Y , we obtain: Finally, this approach consists in adding the image forces 2γ i d(Z,Y i ) N i to the vertex Y i belonging to endocardium surfaces. This is similar to the pro-active deformable model described in [14].

Data assimilation approach
We will show in the following that this minimisation of energy can be related to a data assimilation approach. The methodology of this data assimilation is directly inspired from [7]. In this approach, two parts are taken into account: the electromechanical model described by Equation 3 with inputs consisting in the electrical command and different external loads, and the available observations. We assume that the parameters of the model are known, unlike the initial position condition X 0 on which we make an error of ξ X (X (0) = X 0 + ξ X ).
A new dynamical system called state observer takes as inputs the electrical command and the image data and returns the estimated state, written asX which should converge to the true state X . In classical data assimilation approach, the observation Z (measures) can be directly computed from the true state X , thanks to an observation operator H such that Z = H X . Then, the observations computed from the estimate state (Ẑ = HX) are compared to the measured observations (Z) and the difference (Ẑ − Z) called innovation is taken into account in the sate observer dynamics.
In our case, if we note Z the blood pool surfaces, we no longer have Z = H X since with cine MRI, we cannot track any material points during a cardiac cycle. Instead, we can compare the two surfaces X and Z through a distance map which can be formalized as H(X , Z) = 0. The observation operator is taken as the gradient of the square distance between the two surfaces H(X , Z) = ∇d 2 (Z, X ) = 0 The estimated stateX does not match perfectly with the observation, and therefore the error between the estimated state and the true state can be quantified with ∇d 2 (Z,X) = 2d(Z,X )∇d(Z,X). Note that ∇d 2 (Z,X) is a vector of the same size as X and its velocity components and its components which correspond to points that are not on endocardium surfaces are 0. For points on the endocardium, Furthermore, by definition of a distance map, ∇d i = N i where N i is the normal of the heart mesh at pointŶ i . Then, the built state observer is: with K d the gain associated with the data. We can see that with a high gain, the estimated state will rely more on image data information than in the electromechanical model. Conversely, with no gain, the observer do not take into account the data and is equivalent to the electromechanical model. Thus, the choice of the gain K d depends on the relative confidence in the model and the data.
It is of high interest to analyse the error between the estimated stateX and the true state X in order to choose the gain. With a proper choice of the gain, the error should converge towards zero. We write the error dynamics by subtracting the model (equation 3) from the observer (equation 7): After linearising the data and assuming that the estimated stateX is close to X : where H d (X ) a matrix n × n where n is the size of the state vector X . Its components corresponding to points on endocardium surfaces are the 3 ×3 Hessian matrix of the squared distance d i and are null otherwise. Since the real state X is supposed to coincide with the position and the movement of the apparent boundaries in the image Z, then ∇d 2 (Z, X ) = 0. The error dynamics is: A result of the control theory shows that this error converges to 0 if all eigenvalues of (A + K d H d ) matrix have negative real parts. This provides a criterion for selecting the gain matrix K d .
In practice, we choose the gain K d as in [7] : Indeed, if we decompose the error dynamics, we have: Therefore with this choice of K d , the stiffness of the error dynamics is increased. It implies an increase of the frequency and the damping of the eigenmodes, and therefore a better convergence toward zero. Here we see the difference between this filtering method and Kalman filtering methods such as the one proposed in [18]. The gain K d is not the Kalman gain, so that the result of the filter is not the optimal result in a stochastic way, but K d is chosen in order to ensure the convergence of the errorX toward zero. Although we do not ensure an optimal result, we avoid to compute the inverse of a combination of covariance matrices, thus leading to a much faster filter than the Kalman approach.
We use the Houbolt implicit scheme to integrate equation 7. Since the image term is also made implicit, the generalised stiffness matrix that is involved in the linear system of equations should change at each time step since the matrix H d depends on the position of endocardium vertices Y i . However, modifying the generalised stiffness matrix at each time step implies that a Cholesky decomposition or a preconditioning must be performed at each iteration which is computationally very expensive. Since the stiffness matrix K is constant, we chose to estimate the term γH T d H dŶt+dt numerically, by first computing the positionŶ t+dt as if there were no image forces and then multiplying it by γH T d H d . This proved to be a fairly efficient approach since the preconditioning of the generalised stiffness matrix is only done once. This also gives better results than a semi-implicit scheme where image forces are estimated explicitly.
Finally, one should note that N i is an eigenvector of the Hessian matrix of the distance map d i with eigenvalue 1. Therefore, when using the gain matrix as K d = γM −1 H T d , the dynamic law of the state observer is given by : This corresponds exactly to the formulation we obtained with the deformable model approach.

Validation with synthetic data
In order to validate our state estimation method in a quantitative way, we generated synthetic cine-MR images using the electromechanical model with standard values. We took 29 instants of the second simulated cycle and we generated the corresponding segmented 3D images, using rasterisation of the tetrahedra. As we assume here that the model is known, all parameters of the model used in the state estimation are the same than the ones used to generate the synthetic data. Thus the only error is on the initial position. We can then quantify the evolution of the mean position error in this ideal case.
State error analysis: We observed, as expected, that the root mean squared error (RMSE) decreases with time, under the action of the state estimation filter. Here, the gain γ was set to 0.8. Fig. 5.a shows the evolution of the position error along three cardiac cycles. Fig. 5.b shows the intensity of the contraction forces and the intensity of the image forces for one endocardial vertex and along three cardiac cycles. We can see that the image forces decrease rapidly in the first times of the first cycles and that the images forces remain small compared to the intensity of physical forces such as the contraction forces. We can see also that the image forces do not vanish exactly to zero. The decreasing of this RMSE depends on the spatial resolution of the images.

Results with clinical data
Effect of the spatial resolution of the MR images: The voxel sizes used in the synthetic images are respectively 1 mm, 0.75 mm and 0.5 mm in all three directions. The RMSE decreases if we increase the spatial resolution of the images and seems to converge to values which are smaller than the spatial resolution of the images and which should correspond to numerical approximation errors (see Fig. 5.a).
Effect of the temporal resolution of the MR images: For this we used real images (see details in next section). The first one was a complete cine-MRI sequence (30 images), the second and the third ones were subsamples of the cine-MRI sequence (respectively 15 and 5 images). Fig. 5.c shows that the left ventricle volume is better approximated in the case of sequences with 30 or 15 images than in the case of the sequence of 5 images. Nevertheless, as the contractility of the left ventricle was well calibrated, the knowledge of the model allows us to obtain good information on the left ventricle volume curve, and to compute good approximations of the ejection fraction. The left ventricle ejection fractions obtained respectively from the complete segmented sequence, from the estimations with complete MRI sequence, and with 15 and 5 images sequences are respectively: 59.20%, 59.34%, 57.56% and 56.84%.
Cardiac Function Estimation: Finally, in Fig. 6, the physiological curves obtained from the state estimation are compared with the ones given by the reference simulation. These physiological curves correspond to the right and left ventricular pressures ( Fig. 6.a), volumes ( Fig. 6.b) and flows (Fig. 6.c). In the isovolumic phases, pressures are computed to counterbalance external forces such as contraction forces and image forces in the case of the estimation in order to keep the volume constant. We can see that in these phases, and in the ejection phases in which the pressures depend on flows through the Windkessel model, the pressures are well recovered. We can see also that after a small period due to the initial position error, the volumes and the global evolution of the flow are well recovered. As flows are the derivative of volumes, errors on volumes due to the oscillation of image forces are magnified.

Results with clinical data
Several estimations were made with different values of the gain γ in order to see the effect of the gain on the state estimation. Fig. 7 shows the MRI segmentation at a time t i of the cardiac cycle. The superimposed lines represents the endocardium and epicardium surfaces of two heart meshes obtained with different values of γ.
The higher value of the gain gives more confidence in the data than in the model, then the image forces are larger in this case as we see in Figs. 7.b and 7.c. We can see that the left ventricle is well tracked in the two cases, while the right ventricle is better tracked in the case of the higher gain. It shows that the contractility parameter in the right ventricle does not equal the one in the left ventricle, which we calibrated with the left volume curve obtained from the cine-MRI. Thus, it allows us to detect differences in parameters, which can lead to parameter estimation.
In order to qualitatively evaluate the estimated motion, we used tagged MRI on the same subject to extract the projection of the 3D real cardiac motion in a number of short axis view ( Fig. 8.a). The qualitative comparison with the projection of the 3D estimated motion (Fig. 8.b) is promising, as we observe similar motion patterns. The estimated motion is much smoother due to the influence of the model. We are working on a more quantitative comparison with the estimated motion.

Conclusion
Coupling electromechanical models of the heart with clinical data in order to help diagnosis and therapy planning is still very challenging. This article presents the link between deformable models and data assimilation in order to estimate cardiac motion from cine-MRI. The proposed method allows to keep the low computational cost of deformable models while using a rigourous mathematical framework. Motion recovery is demonstrated on synthetic and real data. These promising preliminary results will be extended in order to perform parameter estimation, which is the ultimate goal of the approach.