1 Introduction
Acording to the World Health Organization, cardiovascular diseases are the leading cause of death in the world. A key factor for early detection and prevention of these diseases is to analyze the cardiac motion to diagnose, for example, valve conditions and motion abnormality. Moreover, the cardiac mechanics can be studied and analyzed through the heart deformation [1]. In most of the medical laboratories, diagnosis is based on the visual inspection of the heart’s motion [2]. However, the results are conditioned to the experience of an expert, and in consequence, they are highly variable and subjective. Thereby, the need of having objective and understandable measures emerge, and uptodate cardiac motion estimation is a central topic in cardiac imaging [3].
Estimation of cardiac motion is a very active problem to be tackled. In the search of achieving a good motion estimation, different authors have used diverse biomedical imaging modalities including: magnetic resonance imaging (MRI), computed tomography (CT), single photon emission computed tomography (SPECT), and positron emission tomography (PET) (e.g. [4, 5, 6]). However, the lack of resolution in modalities such as PET and SPECT ( mm [7]), together with the exposure to radiation (i.e. CT, PET/SPECT), and the magnetic interference and cost of MRI make their use complicated.
An alternative modality to estimate cardiac motion is ultrasound imaging (US). US is very popular due to its low cost, high accessibility, realtime interaction, nonionization, portability, and rapid assessment [8]. It has become routinely used in multiple clinical scenarios including diagnostic and heart disease prevention. US allows capturing, for example, the heart’s size and shape, strain rate, ventricular deformation, and abnormal motions. Moreover, its feasibility for tissue tracking and estimated motion analysis has been demonstrated [9, 10].
Despite these benefits, ultrasound has disadvantages related to the presence of noise and occasional artifacts and as well a limited acquisition speed. The poor temporal resolution of conventional US hinders retrieving different mechanical events of the heart [11], [12]. Nonetheless, recent advances in ultrafast ultrasound (UUS) imaging have overcome previous issues, particularly temporal resolution, as they have a much higher frame rate (greater than 1000fps) [13], which is advantageous for cardiac motion estimation. An UUS system is capable of computing many lines in parallel, generating in this way a full image from one single transmitting event. Different applications of UUS emerge, such as tissue and blood motion estimation and imaging of micro bubbles or neurovascular coupling. Moreover, it facilitates advancements in disease prevention, diagnosis, and therapeutic monitoring [12]. In this work, we study the application of UUS to estimate the cardiac motion.
In the literature, different works for US cardiac motion estimation have been reported. Most of them use conventional ultrasound (e.g. [14, 15, 16, 18]), while few works have been reported using ultrafast ultrasound imaging (e.g. [19, 20]
). Despite the use of conventional or ultrafast ultrasound, the approaches proposed to retrieve the cardiac motion can be classified into those using the characteristics of the radio frequency (RF) signal, and those based on computer vision techniques applied to an image sequence.
Motion estimation techniques in the first category use the natural acoustic reflections of the radio frequency (RF) signal. When this option is used, the cardiac motion can be computed by either applying speckle tracking techniques, which use the amplitude of the signal, or radiofrequencybased correlation techniques, which use the phase information (e.g. [20, 21, 22]). This option is a promising approach to estimate the cardiac motion. However, the result greatly depends on the characteristics of the RF.
In the second category of motion estimation techniques, we can find solutions based on Optical Flow (OF), in which the relative motion of the heart is computed from the velocities of patterns’ brightness (examples can be found in [23, 24, 25]). Although different authors have proved the feasibility of OF for motion estimation, its main drawback is that it only works with relatively noncomplex deformations. Another common solution, usually adapted for more complex motions, is nonrigid registration, in which displacements of the tissue can be tracked by computing the spatial correspondences between frames [14, 26, 27, 28].
Another option for estimating the dynamics of the heart is by tracking the heart’s borders using deformable models (e.g. [29, 30, 32]). However, cardiac motion estimation can be inaccurate when the displacements are parallel to the edge or when there is a lack of welldefined borders, which is a common problem when using US [14].
In this work, we propose a new approach to estimate the cardiac motion using ultrafast ultrasound modality (see Figure. 1). Based on a variational formulation for nonrigid registration in , we include a maximum likelihood type estimator to increase the robustness of the solution in the sense of being able to deal with outliers. While this is an important part of the solution, the main contribution is: combining lowrank data representation with a topologypreserving approach. Particularly:

We promote lowrank data representation. As a standalone tool, it already offers several advantages, such as speeding up the global solution by reducing the computational time and decreasing the noise in the image sequence.

Another key point is topology preservation. A penalization term for the Jacobian determinant is used in order to guarantee a diffeomorphic transformation. We use a regularizer to rule out distortions while controlling the magnitude of expansion and compression.

The combination of the two previous tools turns out to be powerful as it allows computing an accurate displacement field that is mathematically wellmotivated and computationally efficient.
The remainder of this paper is organized as follows. Section 2 presents previous literature related to the topologypreserving problem. Our solution is tackled in Section 3. Particularly, in Subsection 3.1 we present the lowrank data representation strategy, while in Subsection 3.2 we describe our variational framework to recover the deformation over time. Moreover, in Subsection 3.3 we describe our penalization term for achieving topology preservation. In Section 4, we validate our proposal presenting experiments on both, a realistic and a simulated dataset. Section 5 provides a final conclusion and directions for future works.
2 Related Work
During cardiac motion estimation, unpreserved topologies result in violations of region convexity and are reflected in penetration of boundaries and overlapped or distorted mesh elements. Thereby, topology preservation is important in order to ensure connectivity between structures, maintain the relations between neighboring elements, and avoid distortion of existing structures. When the deformations are small, topology is preserved by the smoothness offered by the regularization term, but this is not the case when large deformations appear, as it happens when the complex dynamics of the heart is to be retrieved.
From the conservation principles of continuum mechanics, it is clear that the elastic displacement field in cardiac motion can be modeled as isomorphism resp. diffeomorphism. Necessary and sufficient conditions are that its deformation gradient tensor exists and is nonsingular at every point in the deformable body
[33]. In terms of the Jacobian , this means that exists and that the at every point in the deformable body. For a positive volume, it is required that throughout the deformable body [33].A wellknown approach to achieve topology preservation is controlling the Jacobian determinant. Dacorongna in [34] presented a detailed discussion related to the Jacobian determinant equation in which he demonstrated its ability to achieve topology preservation. Jacobian determinant has been used in problems involving deformable structures in order to achieve more realistic transformations (e.g. [35, 36, 37]). However, in this section we cover only those that have relation with the modeling of deformable objects. In a multidimensional elastic registration framework, authors in [31] explored a barrier function for penalizing locally noninvertible functions.
The problem of achieving topology preservation has been reported in different works addressing MRI and US data. The authors in [35] ensured topology preservation by defining a threshold of 0.5 for the Jacobian determinant. Then, for values lower than this threshold, they generate a new template, equal to the previous deformed template, to continue with the registration process. Similarly, Ashburner et al. in [36], [38] penalized singular values of the Jacobian having lognormal distribution. To enforce the Jacobian positivity, Musse and colleagues [39] proposed a 2D parametric approach based on the constraint of the continuous hierarchical modeling of the deformation field.
Later on, Noblet et al. [40] reported an extension of Musse’s work in the threedimensional space. They presented a hierarchical deformation field model in which the Jacobian determinant was conditioned when it had negative values. The main difference between the two works is that the nontrivial optimization problem in [40] was obtained in the 3D space. Topology preservation was treated as a hard constraint in [41], where they restricted the Jacobian determinant within a set of intervals in a grid region. If these conditions were not met, then they enforced topology preservation in terms of gradients. Explicit control of the deformation in terms of the determinant of the Jacobian was reported in [42]. In comparison with similar works, here, the authors proposed the use of pointwise inequality constraints (i.e. they achieved topology preservation by controlling voxel by voxel instead of using integral measures).
Another approach was presented in [43], in which topology was preserved by quantifying the magnitude of deformations and examining the statistical distributions of Jacobian maps in the logarithmic space. A twostep solution was proposed in [44]
. Authors first corrected the gradient vectors of the deformation and then they reconstructed the deformation based on a minimization problem on a convex subset of the underlying Hilbert space. As a result, they achieved a welldefined Jacobian in the image domain. An extension of that work was presented in
[37], in which authors proposed a solution based on independent problems of small dimension that allow parallel computation.Zhang et al. [15, 18] developed a temporally diffeomorphic motion estimation approach for conventional cardiac ultrasound sequences. In that work, the authors addressed the topologypreservation problem by using a smooth velocity field with a differential operator in a Sobolev space. The resulting transformation defines a group of diffeomorphisms. In [45], authors used the Beltrami coefficient (BC) to represent an orientationpreserving diffeomorphism. To deal with the computational cost of the BC method, they presented a splitting algorithm, one part solved the BC whereas the other involved the quasiconformal map. However, the BC was reduced to constrain the Jacobian. Moreover, a biophysically constrained framework for large deformation diffeomorphic image registration was proposed in [46]. They achieved topology preservation by controlling the Jacobian determinant and the amount of shear in the deformation map using a nonlinear Stroke regularization scheme.
3 Cardiac Motion Recovery
In this section, we present our solution to estimate the cardiac motion. The section contains three parts: First we describe the lowrank representation of the data. Then we give the variational framework to recover the deformation field. Finally, we present our regularizer for topology preservation to obtain a realistic diffeomorphic solution.
3.1 LowRank Data Representation
As a recent promising tool, lowrank data representation has been promoted in a variety of areas, including but not limited to: computer vision, machine learning for fitting problems, and computing the lowrank approximation of a matrix
[47, 48]. Lowrank representation has been useful for processing big data because in many datasets, the relevant information lies in a lowdimensional space [49]. In particular, lowrank representations have been used in motion recovery in optical flow, e.g. in [50], to recover the motion in videos of facial movements.Our motivation for using lowrank data is threefold: First, we aim to increase the computational speed of the solution convergence. Second, it constitutes a means of denoising the ultrasound data and avoiding artifacts in the recovered deformation field. Third, we aim to investigate the synergy of the lowrank representation with the preservation of the topology.
For a precise description of our lowrank data representation, consider an ultrafast ultrasound image sequence, , with frames of size . Then, a new structure of the data is given in the form of a single matrix, called Casorati matrix C, which columns are the frames in a vectorized way:
(1) 
where is the scalar value of the sequence in frame , at a given pixel location .
U  S  V  

13824x13824  13824x1814  1814x1814  
13824x100  100x100  1814x100 
Let us now turn to the theoretical prerequisites to produce a lowrank representation of the Casorati representation C, exploiting the high correlation in the columns of the matrix.
It is wellknown (see, for instance, [51]) that any matrix admits a decomposition in the way:
(2) 
where and are orthogonal matrices; the matrix S is diagonal with and the scalar quantities are called the singular values of the matrix A. Then, the triplet USV
is called a singular value decomposition (SVD).
A major issue in medical applications is the large amount of data to be processed, which is the case when the cardiac motion is estimated. An option to deal with this issue is to compute the set of dominant singular values instead of computing the singular value decomposition of a large and dense matrix. This allows keeping the most relevant information in a subspace which is smaller than the original one. Thus, the problem of building a lowrank representation can be given by finding the dominant singular values of . Mathematically, finding the rank approximation of the matrix can be described as:
(3) 
where is the squared Frobenius norm of a matrix. It is used in lowrank based problems since it is invariant to rotations and to the rank. The importance of the rankapproximation in (3) is given by the EckhartYoung theorem (see, for instance, [52]). Take a matrix with a SVD as in (2), and let . Let be the rank approximation in equation (3). Then, we have:
(4) 
For our data described in Section 4, we compute the Casorati matrix and different rankapproximations in equation (3). The resulting matrix decomposition is displayed in Table 1, in which the original matrix (of size 13824x1814) was reduced to rank 100. Both the error and the computational time are plotted against the rank in Figure 2. This figure shows that only an average of seconds was needed to obtain , and that the approximation error for is .
As we stated before, another main motivation to promote lowrank representation instead of using the full Casorati matrix is to reduce noise. This is accomplished by eliminating the subspace where the noise relies, which results in retrieving a subspace with only relevant information. Noise, which normally relies on another subspace due to its characteristics, is rejected from the solution (see Figure 3). This eliminates artifacts in the subsequent deformation computation.
In the last step, we use the invertibility in equation (1) to invert the lowrank representation back to the denoised video sequence . We will work with that sequence in the following section to extract the mechanical deformation field of the heart.
Let us now expose the variational framework to recover the deformation.
3.2 Deformation Recovery
Let be the image sequence over frames, where each image is a function over the bounded domain .
We will find a deformation vector field defined on the domain as a minimizer of an energy functional:
(5) 
where the three terms used have the following purposes: discrepancy measure, regularization term and topology preservation.
We will now go through our variational framework and explain each of the terms in the energy functional. We begin with the representation of the deformation field , then go on to the discrepancy term and the regularization term . In Subsection 3.3, we describe the topology preservation term necessary to achieve realistic deformations.
How to represent the deformation field ? The deformation model is an essential factor that defines how fast and accurate the approach is. In order to find a compromise between computational cost and accuracy, we will handle the changes over time using a lattice as in the following definition:
Definition 1
A mdimensional lattice is the linear span of a set of linearly independent vectors in .
We will then use a lattice in which its points are characterized by the tensor product of the bsplines [53]. These are widely used in medical applications. The advantage of this lattice deformation model is that it demands low running time, allows multiresolution, have optimal mathematical properties, and keeps affine invariance. Using bsplines has the additional advantage of being able to handle complex deformations.
Consider a given position in . Let be a basis of spline functions and let be control points. Then, we express the deformation vector at point through the model:
(6) 
In this work, we use cubic basis splines:
(7) 
The deformation model for in (6) is in and it shows that within our framework, we actually reconstruct the control points in order to get the deformation . – In our application, we will exploit this deformation model for dimension .
We now turn to the discrepancy term . Since the images are acquired by the same sensor, it is not expected that they have a big intensity variation between them. Therefore, an iconic method is a perfect match for this application. One possible option is to use the Sum of Squared Differences (SSD) method:
(8) 
SSD offers a low computational cost but it has the disadvantage of not dealing well with outliers.
For the actual expression for in (5), we use:
(9) 
Here, the function is a maximum likelihood type estimator motivated by robust statistics:
Definition 2
An Mestimator is a symmetric and positive definite function with a unique minimum at zero.
The Mestimator substitutes the minimization of , where is the residual error, with , in order to deal with the effect of outliers. This increases the robustness and accuracy of the result.
In this work, we use the Turkey estimator for :
(10) 
where is a tunning parameter. The discrepancy term with the Turkey estimator in (10) is known to offer a hard rejection of outliers [54].
For the regularization term , we use the Tikhonov method to impose stability to the energy functional. Let be the regularization parameter. Then we use for the term:
(11) 
as the regularizer in our variational framework (5).
3.3 Topology Preservation
A common drawback of most of the solutions working with complex deformations is that topologypreserving is not guaranteed, which leads to unrealistic deformations becoming a source of error. An example of such deformations is shown in Figure 4 in which topology is not preserved. In medical applications, this issue is of huge importance in order to maintain the anatomical structures. In order to preserve the topology, the deformation must be a diffeomorphism. One way to achieve topology preservation is using the Jacobian determinant. This option makes sense since it allows measuring the changes in the area/volume produced by the deformation at each patch.
Condition  Local type of deformation 

topology is destroyed, overlapping, distortion and penetration of boundaries may occur  
diffeomorphism, topology preservation  
contraction  
volume preservation  
expansion 
Let be the Jacobian determinant of the deformation in . The Jacobian determinant is described as:
(12) 
Table 2 shows the characteristics of the deformation encoded in the Jacobian and an illustration of these behaviors can be seen in Figure 4.
From our deformation model in (6), the partial derivatives of can be easily evaluated, as they come from the tensor product of independent functions. Using Equation (7), we have for the derivatives of the cubic basis splines:
(13) 
Then, the determinant in (12) can be straightforwardly evaluated as a function of the lattice points characterizing the deformation .
Once we have stated how to compute the Jacobian determinant, we turn to formulate our penalization term for the energy functional (5). We use a weak constraint, but do not penalize values lying near 1:
(14)  
Here offers a balance in our penalization, and is the margin of acceptance for values close to one.
The first term
heavily penalizes negative values of the deformation. Thus it prevents the field from having distortions or penetration of boundaries. The term
with the parameter controls the magnitude of the expansions and contractions.
Unlike most of the Jacobian determinant constrains, for example (see [36]) and (as in [31]), we do not only guarantee the positivity of the Jacobian determinant but also enforce its value to stay near one. Moreover, we penalize big expansions in order to achieve more realistic deformations.
There are different options for solving the regularized class. Traditional methods include Gradient Descent, Newton’s method, Nonlinear Conjugate Gradient, and Evolutionary Optimization Algorithms. However, they can get stuck at local minima, might need an infinite number of iterations to converge, or have a slow rate of convergence. A better alternative is the wellknown LevenbergMarquardt (LM) which offers better results with less computational time. For these reasons, in this work we use LM method to minimize our energy functional.
4 Experimental Results
This section describes in detail the experiments that we conducted to validate our proposal.
4.1 Subjects and acquisition
We used two ultrafast ultrasound datasets (see Figure 5) to evaluate our proposed approach.
The first is a realistic dataset from one patient. During the acquisition, the patient was placed in the supine position or left lateral decubitus. Then, the probe was placed on the left parasternal line at the fourth intercostal space with the marker pointing toward the right shoulder of the patient. The images were thus taken from the parasternal long axis view. This view is useful for global assessment of the motion of the heart’s wall and the function of different areas including the right and left ventricle, the mitral and aortic valves and the interventricular septum. This dataset is composed of frames with size of x, and a scaled version of size x. This data was acquired with an UUS device with a (fc) transducer with a bandwidth of 6MHz and 192 elements using coherent compounding of plane waves. The output is a 2D plane wave image sequence of the long axis view of a healthy heart.
The second dataset (see [19] for details) is a simulated ultrafast ultrasound sequence that has a realistic cardiac deformation field and describes the mechanics of a healthy left ventricle (see [55] for the description of the mechanics). This data was generated using (fc) transducer with a frame rate of 5000 Hz (single transmits) and effective frame rate of 500 Hz. The output of this simulated data is a set of 2D apical imaging planes. This view is helpful to study the left ventricle and the mitral inflow. The sequence is composed of frames with size of x. This simulated data is very helpful to assess the performance of our approach since it includes a ground truth of the displacements which can be compared against our estimation.
All the measurements and reconstructions in this section are taken from these two ultrafast ultrasound cardiac sequences. All test and comparison were run under the same condition on a CPUbased implementation. We used an Intel(R) Core i7 6700 CPU at 3.40GHz32GB RAM, and a Nvidia GeForce GT 610.
Preprocessing 

Minimum 

Minimum 

Minimum  

LowRank (SVD) 


Wavelet  
Gaussian Smoothing 
4.2 Validation scheme
We divided our validation scheme into two parts. The first part makes use of the realistic dataset and relies on the following measurements to evaluate the performance of our topologypreserving technique:

Comparison between our proposed topology regularizer and two from the literature: Table 3;

Assessment of using lowrank as a preprocessing step: Figure 6;

Inspection of the displacement field: Figure 7 (A)/(B);

Numerical results offered by the Jacobian determinant: Figure 7 (A.2)/(B.2); Table 4;

Careful comparison of the residual error for both the lowrank tool and the topology preservation tool and study of their synergy: Table 4;
In the second part, we use a simulated dataset with a provided ground truth in which we perform the following evaluations:

Numerical visualization and comparison of the displacement field: Figures 8 and 9;

Numerical visualization of mean accumulated displacement of the seven segments of the left ventricle Figure 9;

Nonparametric statistical analysis between the real and the estimated displacements;

Illustration of the computed strain as a reasonable clinical measure: Figure 10.
4.3 Results
In order to prove the benefits of using lowrank (SVD) as a preprocessing step, we compared it against two common preprocessing techniques: Gaussian smoothing (kernel and ) and Wavelets (Biorthogonal Spline Wavelet, 4 levels). We carried out the comparison of the three preprocessing techniques using our topology regularizer and two more from the state of the art [17, 16].
According to the results (Table 3 and Figure 6), we found that lowrank was able to find the best minima in our case study in a computationally efficient manner. The results showed that Wavelet was able to find an acceptable minima but, it needed an average of 30 iterations per frame to converge compared to the 16 needed by lowrank (see plots in Figure 6). Gaussian smoothing on the other hand performed the worst in terms of minima and average iterations per frame.
Overall, out of the three prepossessing techniques, lowrank offered a good tradeoff between accuracy and computational time since it requires less iterations per frame. This is further reflected in the overall CPU time as illustrated in the boxplot at right side of Figure 6 where lowrank only needed an average of 2.3217 seconds of computational time while Gaussian and Wavelet both required more than 11 seconds in average. Moreover, it performed the best across the different regularizers giving the best minima each time. In general, our topology regularizer approach was the best compared to the regularizers offered by Rohlfing [17] and Heyde [16]. It has the advantage of enforcing the value to be close to one and penalizing very strong expansions and contractions.
Exp.  Data 


Minimum 




1)  FullRank  12.0830  0.1069  0.0016  [2.8896 3.6884]  
2)  FullRank  (14) with  12.0672  [0.0617 1.0096]  
3)  FullRank  (14) with  11.3953  [0.9034 1.0023]  
4)  FullRank  (14) with  11.0304  [0.9546 1.0059]  
5)  LowRank  3.8771  7.9505e06  [1.0144 3.7375]  
6)  LowRank  (14) with  3.4895  [0.0213 1.0031]  
7)  LowRank  (14) with  3.0214  [0.9184 1.0028]  
8)  LowRank  (14) with  2.32170  [0.9950 1.0100] 
The performance of our proposed topologypreserving technique is illustrated in Figure 7 with and without topology preservation. First we show the performance without topology preservation ( in the energy functional (5)) where the resulting displacement fields are displayed in (A) and the corresponding Jacobian determinant are shown in (A.2). The example frames in the columns show critical details of the respective displacements fields and Jacobian. Row (A.1) shows zoomin parts where different violations of the topology occur in part (A), such as overlapping, boundaries penetration, and mesh elements distortion.
The plots of the Jacobian clearly show huge variations in the value of the determinant which results in an unstable representation of the anatomical structure. Not only that, but in some parts the Jacobian determinant presents big values (greater than 3), which indicates that some transformations produced very big expansions, while in others the Jacobian determinant gave negative values, which indicates that new structures were formed. These violations create new structures and result in an unrealistic representation of the complex deformation of the heart’s motion. This is not acceptable particularly in medical applications where preservation of the anatomical structure is remarkably needed.
We then ran the same tests after applying the proposed topology preserving approach ( in (14) with ) and show the results in part B. Looking at the zoom in parts in B.1, and comparing it with part A.1, we can verify that our approach allows controlling expansions and contractions, maintaining region convexity, and avoiding foldings. The corresponding plots of the Jacobian determinant are shown in part B.2. In comparison to the Jacobians without topology preservation, we can see stabilities on the values as they do not suffer large variations (mostly stay on 1) with guaranteed positivity. These are realistic values, as over two pairs of timeframes, the volume should be approximately preserved. From these illustrations one can therefore conclude that the approach was successful in avoiding topology violation even with complex deformations of the cardiac motion.
For a more detailed quantitative analysis, we evaluated the global performance of our approach by comparing its performance with lowrank and fullrank data using frames of the sequence (see Figure 5 dataset I). Unlike works where speckles were useful as indicators for tracking the heart’s motion, in the framework we proposed, results from Table 4 show that promoting lowrank offered a positive effect on the solution as it significantly reduced the computational time and allowed faster convergence of the energy functional (less iterations per frame).
Table 4, Exp. 1 and Exp. 5 show that without topology preservation in the functional (5), the residuum was about with and without the lowrank constraint. Applying lowrank, the Jacobian determinant had higher values but still suffered from heavy distortions. With topology preservation in Exps. 24, and 68, reasonable values of the Jacobian determinant were obtained, avoiding any penetration of boundaries. Thus, topology preservation is a necessary tool to get realistic deformation results for cardiac motion estimation.
The primary role of the lowrank representation seems to be the radical decrease in computational time, as seen from a comparison of Exps. 14 and Exps. 58, where the computational time was decreased by about 75 %.
The residual errors in Exps. 24 in the topologypreserved fullrank case were in the order of magnitude resp. . Contrary to that, topology preservation in the lowrank case in Exps. 68 yielded minima in the order of magnitude to . Moreover, the discrepancy error showed that lowrank achieved an order of magnitude in comparison with given by the fullrank case. In practice, topology preservation and lowrank constraint act together to get a more realistic deformation field in less time.
As a second part of our validation scheme and for an extended evaluation of our proposal, we illustrate the axial and lateral accumulated displacement of both the ground truth and our estimation of a single heart cycle (left side of Fig. 8). It is clear by visual inspection of the colored bar that our estimation is very close to the ground truth. In order to support this statement, we computed the RootMean Square error (RMSE) for the axial and lateral displacement and plotted the results in left side of Figure 9 where we can see RMSE values less than 1mm for the axial direction and less than 1.2mm for the lateral direction. The plots also show a concentration of values much lower than 1mm in both displacement directions.
To complement the analysis of the estimated displacement, in right side of Figure 9 we offer an analysis of the seven segments of the left ventricle during one cycle of the heart. Since the inferior and anterior parts are symmetric, we can evaluate their behavior for the axial and lateral displacements. As expected, axial displacements (blue circle) reported positive values acting in a similar way in both sides while lateral displacements (red squared) showed similar behavior with contrary signs since they go to opposite directions during heart motion.
We used the nonparametric Wilcoxon signed rank sum test to answer whether there is statistical significant difference between the real values (ground truth) and the estimated values. We found that the null hyphotesis was not rejected with significance level. This, lead us to conclude that we obtained a good estimation of the displacement.
To further support the results obtained in Figure 7 – (A.1 and A.2) and Table 4, the right side of Fig. 8 shows the corresponding Jacobian determinant of the illustrated frames where we can see that most of the values meet the desired criteria (stay at 1) for achieving topology preservation, which proofs the efficiency of our proposed term.
Finally, we provide strain images as these are often used clinically. It is wellknown that the strain can be calculated in terms of the components of the displacement vector field:
where is the Green deformation tensor, is the displacement gradient, and is the identity.
Strain is useful to evaluate the heart muscle and to identify subtle changes in heart function [56]. Moreover, it allows representing the percentage change in dimension from a resting state to a stressed state (after applying a force). Figure 10 shows the radial and longitudinal strains related to the left ventricle. The resulted plots can be evaluated according to the sign of the strain in which negative values indicate shortening and positive values denote stretching. According to the results at the upper part of the figure, radial strain reported a stretching behavior while longitudinal strain a shortening behavior. For illustration purposes, the radial strain of few frames from the sequence are displayed at the lower part of the same figure. The strain profiles and displacements of the left ventricle exhibit distinct features and clinically meaningful motion patterns.
5 Conclusion
In this paper, we presented a new approach to estimate the cardiac motion using ultrafast ultrasound data. In a variational framework, we combined a penalizer for topology preservation with a lowrank data representation. Together with the better temporal resolution of ultrafast ultrasound, our proposed approach overcame challenges of nonrigid registration, including noise and complex heart motion and inaccurate results exhibiting distortions. While keeping the computational time relatively low, a realistic and clinically meaningful displacement field was produced, with the diffeomorphic features and preserved structures.
In our variational framework, the displacement was represented by a lattice with splines, and a maximum likelihood estimator was used in the discrepancy term to provide robustness against outliers. The regularizer for the topology has two features: eliminating radically negative values and carefully controlling the volume expansion and compression. We validated the accuracy of our approach and showed that it offers a RMSE lower than 1 mm in comparison to the ground truth.
While this variational framework gives good results and is strong individually, the consuming CPU time and artifacts in the ultrafast ultrasound sequence motivated us to promote a lowrank data representation, as it has proved to be useful in other areas of imaging and computer vision. We represented the data in a single Casorati matrix and used the dominant singular values to compute the deformation. Apart from removing the noise in the ultrasound data, this technique greatly reduces the computational time and produces, together with the topology penalization term, significantly less discrepancy in the results.
While we wanted to show the potentials of combining ultrafast ultrasound with lowrank techniques and topology preservation from a technical point of view, the objective of this work is to have a first study for proof of concept and to open a new line of research for further clinical investigation. In this work, we evaluated the technique using data from one subject. Future work will include a more extensive evaluation with more subjects to examine the clinical potentials of the approach. Moreover, the technique promises to be useful for analyzing organs experiencing complex motion other than the heart, as for example the movement of the lungs in respiration.
Acknowledgment
We would like to thank Chris de Korte for discussion and coordination and Peter Bovendeerd for the initial code for the simulated data.
Support from the ERC Advanced Grant Project MULTIMOD 267184 is greatly acknowledged. Also, support from a FPU national scholarship with reference AP20121943.
References
 [1] N. Duchateau, B. Bijnens, J. D’hooge, and S. M., “Cardiac motion and deformation,” Chapter Book in 3D Echocardiogra phy, Second Edition, 2013.
 [2] B.D. Hoit, “Strain and strain rate echocardiography and coronary artery disease,” Circ Cardiovasc Imaging, no. 4, pp. 179–190, 2011.
 [3] A. LopezPerez, R. Sebastian, and J. Ferrero, “Threedimensional cardiac computational modelling: methods, features and applications,” BioMedical Engineering OnLine, pp. 1–31, 2015.
 [4] M. Prummer, J. Hornegger, G. Lauritsch, L. Wigstrom, E. GirardHughes, and R. Fahrig, “Cardiac carm ct: A unified framework for motion estimation and dynamic ct,” IEEE Transactions on Medical Imaging, pp. 1836–1849, 2009.
 [5] L. Wang, A. Basarab, P. Girard, P. Croisille, P. Clarysse, and P. Delachartre, “Analytic signal phasebased myocardial motion estimation in tagged MRI sequences by a bilinear model and motion compensation,” Medical Image Analysis, pp. 149–162, 2015.
 [6] J. Royueladel Val, L. CorderoGrande, F. SimmrossWattenberg, M. MartínFernández, and C. AlberolaLópez, “Nonrigid groupwise registration for motion estimation and compensation in compressed sensing reconstruction of breathhold cardiac cine MRI,” Magnetic Resonance in Medicine, 2015.
 [7] C. Nappi, W. Acampa, T. Pellegrino, M. Petretta, and A. Cuocolo, “Beyond ultrasound: advances in multimodality cardiac imaging,” Journal in Internal and Emergency Medicine, vol. 10, pp. 9–20, 2015.
 [8] K. Spencer, “Focused cardiac ultrasound,” Book Chaper in ASE’s Comprehensive Echocardiography, 2015.
 [9] V. De Luca, G. Székely, and G. Tanner, “Estimation of largescale organ motion in bmode ultrasound image sequences: A survey,” Journal in Ultrasound in Medicine and Biology, pp. 3044 – 3062, 2015.
 [10] K. Shung, Ultrasound: Imaging and Blood Flow Measurements, Second Edition. Taylor and Francis Group, 2015.
 [11] M. Cikes., L. Tong, G.R. Sutherland, and J. D’hooge, “Ultrafast cardiac ultrasound imaging: Technical principles, applications, and clinical benefits,” JACC: Cardiovascular Imaging, vol. 7, no. 8, pp. 812–823, 2014.
 [12] M. Tanter and M. Fink, “Ultrafast imaging in biomedical ultrasound,” Ultrasonics, Ferroelectrics, and Frequency Control, IEEE Transactions on, vol. 61, no. 1, pp. 102–119, 2014.
 [13] M. Tanter, “The advent of ultrafast imaging in biomedical ultrasound,” Med. Phys., no. 42, 2015.
 [14] M. J. LedesmaCarbayo, J. Kybic, M. Desco, A. Santos, M. Sühling, P. Hunziker, and M. Unser, “Spatiotemporal nonrigid registration for ultrasound cardiac motion estimation,” Medical Imaging, IEEE Transactions on, pp. 1113–1126, 2005.
 [15] Z. Zhang, D. Sahn, and X. Song, “Diffeomorphic cardiac motion estimation with anisotropic regularization along myofiber orientation,” WBIR 2012. LNCS, pp. 199–208, 2012.
 [16] B. Heyde, M. Alessandrini, J. Hermans, D. Barbosa, P. Claus and J. D’hooge, “Anatomical Image Registration Using Volume Conservation to Assess Cardiac Deformation From 3D Ultrasound Recordings,” in IEEE Transactions on Medical Imaging, pp. 501–511, 2016.
 [17] T. Rohlfing, C. Maurer, D. Bluemke and M. Jacobs, “Volumepreserving nonrigid registration of MR breast images using freeform deformation with an incompressibility constraint,” in IEEE Transactions on Medical Imaging, pp. 730–741, 2003.
 [18] Z. Zhang, M. Ashraf, D. J. Sahn, and X. Song, “Temporally diffeomorphic cardiac motion estimation from threedimensional echocardiography by minimization of intensity consistency error,” Medical Physics, pp. 1–16, 2014.
 [19] M. M. Nillesen, A. E. C. M. Saris, H. H. G. Hansen, S. Fekkes, F. J. van Slochteren, P. H. M. Bovendeerd, and C. L. De Korte, “Cardiac motion estimation using ultrafast ultrasound imaging tested in a finite element model of cardiac mechanics,” Functional Imaging and Modeling of the Heart, pp. 207–214, 2015.
 [20] S. Salles, A. Chee, D. Garcia, A. Yu, D. Vray, and H. Liebgott, “2d arterial wall motion imaging using ultrafast ultrasound and transverse oscillations,” Ultrasonics, Ferroelectrics, and Frequency Control, IEEE Transactions on, pp. 1047–1058, 2015.
 [21] I. Dydenko, D. Friboulet, J. Gorce, J. D’hooge, B. Bijnens, and I. Magnin, “Towards ultrasound cardiac image segmentation based on the radiofrequency signal,” Medical Image Analysis, vol. 3, pp. 353 – 367, 2003.
 [22] R. G. Lopata, M. M. Nillesen, J. M. Thijssen, L. Kapusta, and C. L. de Korte, “Threedimensional cardiac strain imaging in healthy children using rfdata,” Ultrasound in Medicine & Biology, vol. 9, pp. 1399 – 1408, 2011.
 [23] D. Tenbrinck, S. Schmid, X. Jiang., K. Schafers., and J. Stypmann, “Histogrambased optical flow for motion estimation in ultrasound imaging,” Journal Math Imaging Vis, pp. 138–150, 2013.
 [24] S. Vyas, J. Gammie, and P. Burlina, “Computing cardiac strain from variational optical flow in fourdimensional echocardiography,” 2014 IEEE 27th International Symposium on ComputerBased Medical Systems, pp. 149–152, 2014.
 [25] H. Gao, N. Bijnens, D. Coisne, M. Lugiez, M. Rutten., and J. D’hooge, “2d left ventricular flow estimation by combining speckle tracking with NavierStokesbased regularization: An in silico, in vitro and in vivo study,” Journal in Ultrasound in Medicine and Biology, pp. 99 – 113, 2015.
 [26] C. Metz, S. Klein, M. Schaap, T. Van Walsum, and W. Niessen, “Nonrigid registration of dynamic medical imaging data using nd + t bsplines and a groupwise optimization approach,” Medical Image Analysis, no. 15, pp. 238 – 249, 2011.
 [27] O. Oktay, A. Schuh, M. Rajchl, K. Keraudren, A. Gomez, M. Heinrich, G. Penney, and D. Rueckert, “Structured decision forests for multimodal ultrasound image registration,” Medical Image Computing and ComputerAssisted Intervention, pp. 363–371, 2015.
 [28] R. Morin, A. Basarab, S. Bidon, and D. Kouamé, “Motion estimationbased image enhancement in ultrasound imaging,” Elsevier Journal in Ultrasonics, no. 60, pp. 19–26, 2015.
 [29] C. Ahn and J. Seo, “Myocardial motion tracking method integrating localtoglobal deformation for echocardiography,” 2012 IEEE International Ultrasonics Symposium (IUS), pp. 1948–5719, 2012.
 [30] X. Huang, D. Dione, C. Compas, X. Papademetris, B.A Lin, A. Bregasi, A. Sinusas, L. Staib, and J. Duncan, “Contour tracking in echocardiographic sequences via sparse representation and dictionary learning,” Medical Image Analysis, pp. 253 – 271, 2014.
 [31] J. Kybic, M. Unser, and J. Marques, “Multidimensional elastic registration of images using splines,” Image Processing, 2000. Proceedings. 2000 International Conference on, vol. 2, pp. 455–458, 2000.
 [32] F. Milletari, A. S.A., D. Kroll, C. Hennersperger, F. F Tombari, A. Shah, A. Plate, K. Boetzel, and N. Navab, “Robust segmentation of various anatomies in 3d ultrasound using hough forests and learned data representations,” Medical Image Computing and ComputerAssisted Intervention, no. 15, pp. 111–118, 2015.
 [33] W. Slaughter, The Linearized Theory of Elasticity. Springer ScienceBusiness Media, LLC, 2002.

[34]
B. Dacorogna and J. Moser
On a partial differential equation involving the Jacobian determinant
. Ann. Inst. H. Poincaré Anal. Non Linéaire, 2002.  [35] G. Christensen, R. Rabbitt, and M. Miller, “Deformable templates using large deformation kinematics,” IEEE Transactions on Image Processing, 1996.
 [36] J. Ashburner, J. Andersson, and K. J. Friston, “Highdimensional image registration using symmetric priors,” 1999.
 [37] S. Ozere and C. Le Guyader, “Topology preservation for imageregistrationrelated deformation fields,” COMMUN. MATH. SCI., vol. 13, no. 5, pp. 1135–1161, 2015.
 [38] J. Ashburner and K. J. Friston, “Voxelbased morphometry–the methods,” NeuroImage, pp. 805 – 821, 2000.

[39]
O. Musse, F. Heitz, and J. Armspach, “Topology preserving deformable image matching using constrained hierarchical parametric models,”
Image Processing, IEEE Transactions on, pp. 1081–1093, 2001.  [40] V. Noblet, C. Heinrich, F. Heitz, and J.P. Armspach, “3d deformable image registration: a topology preservation scheme based on hierarchical deformation models and interval analysis optimization,” Image Processing, IEEE Transactions on, pp. 553–566, 2005.
 [41] B. Karacali and C. Davatzikos, “Estimating topology preserving and smooth displacement fields,” Medical Imaging, IEEE Transactions on, vol. 23, no. 7, pp. 868–880, 2004.
 [42] E. Haber and J. Modersitzki, “Image registration with guaranteed displacement regularity,” Int. J. Comput. Vision, vol. 71, no. 3, pp. 361–372, 2007.

[43]
I. Yanovsky, P. Thompson, S. Osher, and A. Leow, “Topology preserving
logunbiased nonlinear image registration: Theory and implementation,” in
Computer Vision and Pattern Recognition, 2007. CVPR ’07. IEEE Conference on
, 2007, pp. 1–8.  [44] C. Le Guyader, D. Apprato, and C. Gout, “On the construction of topologypreserving deformation fields,” Image Processing, IEEE Transactions on, vol. 21, no. 4, pp. 1587–1599, 2012.
 [45] L. Lui and T. Ng, “A splitting method for diffeomorphism optimization problem using beltrami coefficients,” Journal of Scientific Computing, vol. 63, no. 2, pp. 573–611, 2015.
 [46] A. Mang and G. Biros, “Constrained regularization schemes for diffeomorphic image registration,” arXiv:1503.00757 [math.OC], pp. 1–29, 2015.
 [47] E. Candès, C. SingLong, and J. D. Trzasko, “Unbiased risk estimates for singular value thresholding and spectral estimators,” IEEE Transactions on Signal Processing, pp. 4643–4657, 2013.
 [48] A. Newson, M. Tepper, and G. Sapiro, “Lowrank spatiotemporal video segmentation,” in Proceedings of the British Machine Vision Conference (BMVC), 2015.
 [49] B. Haeffele, E. Young, and R. Vidal, “Structured lowrank matrix factorization: Optimality, algorithm, and applications to image processing,” In Proceedings of the 31st International Conference on Machine Learning , pp. 2007–2015, 2014.
 [50] R. Garg, A. Roussos, and L. Agapito, “A variational approach to video registration with subspace constraints,” International Journal of Computer Vision, pp. 286–314, 2013.
 [51] G. Strang, “Linear algebra and its applications,” Thomson Brooks/Cole, 2005.
 [52] G. H. Golub and C. F. Van Loan, “Matrix computations,” TheJohns Hopkins University Press, Baltimore, MD, 1996.
 [53] M. Unser, “Splines: A perfect fit for signal and image processing,” IEEE Signal Processing Magazine, pp. 22–38, 1999.
 [54] C. V. Stewart, “Robust parameter estimation in computer vision,” SIAM Rev., vol. 41, no. 3, 1999.
 [55] P. H. M. Bovendeerd, W. Kroon and T. Delhaas,“Determinants of left ventricular shear strain” Am. J. Physiol. Heart Circ. Physiol. 297, pp. 1058–1068, 2009.
 [56] T. Abraham and R. Nishimura, “Myocardial strain: can we finally measure contractility?” Journal of the American College of Cardiology, pp. 731 – 734, 2001.
Comments
There are no comments yet.