Robust Cardiac Motion Estimation using Ultrafast Ultrasound Data: A Low-Rank-Topology-Preserving Approach

10/26/2016 ∙ by Angelica I. Aviles, et al. ∙ Universitat Politècnica de Catalunya 0

Cardiac motion estimation is an important diagnostic tool to detect heart diseases and it has been explored with modalities such as MRI and conventional ultrasound (US) sequences. US cardiac motion estimation still presents challenges because of the complex motion patterns and the presence of noise. In this work, we propose a novel approach to estimate the cardiac motion using ultrafast ultrasound data. -- Our solution is based on a variational formulation characterized by the L2-regularized class. The displacement is represented by a lattice of b-splines and we ensure robustness by applying a maximum likelihood type estimator. While this is an important part of our solution, the main highlight of this paper is to combine a low-rank data representation with topology preservation. Low-rank data representation (achieved by finding the k-dominant singular values of a Casorati Matrix arranged from the data sequence) speeds up the global solution and achieves noise reduction. On the other hand, topology preservation (achieved by monitoring the Jacobian determinant) allows to radically rule out distortions while carefully controlling the size of allowed expansions and contractions. Our variational approach is carried out on a realistic dataset as well as on a simulated one. We demonstrate how our proposed variational solution deals with complex deformations through careful numerical experiments. While maintaining the accuracy of the solution, the low-rank preprocessing is shown to speed up the convergence of the variational problem. Beyond cardiac motion estimation, our approach is promising for the analysis of other organs that experience motion.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 5

page 9

page 11

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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 up-to-date 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, real-time interaction, non-ionization, 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 radio-frequency-based 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 non-complex deformations. Another common solution, usually adapted for more complex motions, is non-rigid 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 well-defined borders, which is a common problem when using US [14].

Figure 1: Overview of our proposed approach. (From left to right) an ultrafast ultrasound cardiac sequence is first acquired, this data is then represented in low-rank in order to speed up the solution and reduce noise. Later, cardiac motion is computed enforcing topology preservation which allows keeping the anatomical structure of the heart. Finally, analysis of the results is offered.

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 non-rigid 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 low-rank data representation with a topology-preserving approach. Particularly:

  • We promote low-rank data representation. As a stand-alone 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 well-motivated and computationally efficient.

The remainder of this paper is organized as follows. Section 2 presents previous literature related to the topology-preserving problem. Our solution is tackled in Section 3. Particularly, in Subsection 3.1 we present the low-rank 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 well-known 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 non-invertible 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 three-dimensional 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 point-wise 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 two-step 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 well-defined 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 topology-preservation 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 orientation-preserving 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 quasi-conformal 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 low-rank 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 Low-Rank Data Representation

As a recent promising tool, low-rank 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 low-rank approximation of a matrix

[47, 48]. Low-rank representation has been useful for processing big data because in many datasets, the relevant information lies in a low-dimensional space [49]. In particular, low-rank 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 low-rank 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 low-rank representation with the preservation of the topology.

For a precise description of our low-rank 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
Table 1: Decomposition of the Casorati matrix and the rank 100-approximation
Figure 2: Left: Singular values . Right: CPU time to compute the rank- approximation .

Let us now turn to the theoretical prerequisites to produce a low-rank representation of the Casorati representation C, exploiting the high correlation in the columns of the matrix.

It is well-known (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 low-rank 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 low-rank based problems since it is invariant to rotations and to the rank. The importance of the rank--approximation in (3) is given by the Eckhart-Young 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 rank--approximations 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 .

Figure 3: Top row shows an original and denoised frame after applying low-rank process along with the rejected space. Bottom row, A and B, show zoom in views of the same frames in which we can see that both noise and some artifacts were removed.

As we stated before, another main motivation to promote low-rank 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 low-rank 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 m-dimensional 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 b-splines [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 b-splines 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 M-estimator is a symmetric and positive definite function with a unique minimum at zero.

The M-estimator 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 topology-preserving 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.

Figure 4: When topology-preserving is not enforced, unrealistic transformations result. A way to ensure topology preservation is by checking the Jacobian determinant . When it is equal to 1 then the volume is preserved. Small positive or large positive numbers of result in contractions or expansions. But having can result in distortions, overlapping, and creation of new structures.
Condition Local type of deformation
topology is destroyed, overlapping, distortion and penetration of boundaries may occur
diffeomorphism, topology preservation
contraction
volume preservation
expansion
Table 2: Jacobian determinant conditions for topology preservation

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 well-known Levenberg-Marquardt (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 CPU-based implementation. We used an Intel(R) Core i7- 6700 CPU at 3.40GHz-32GB RAM, and a Nvidia GeForce GT 610.

Figure 5: Sample frames of the raw data extracted from the two datasets used for evaluating our approach.
Preprocessing
Topology Regularizer
Our Approach
Minimum
Topology Regularizer
Rohlfing et al. [17]
Minimum
Topology Regularizer
Heyde et al. [16]
Minimum
Low-Rank (SVD)
with as in Eq. [14]
Wavelet
Gaussian Smoothing
Table 3: Performance comparison between our proposed approach and other state of the art approaches

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 topology-preserving technique:

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

  • Assessment of using low-rank 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 low-rank 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 low-rank (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 low-rank 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 low-rank (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, low-rank 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 box-plot at right side of Figure 6 where low-rank 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.

Figure 6: The first row shows the convergence history of the complete sequence (1000 frames) using three different pre-processing techniques, while the second row shows the number of iterations it took to find the minimum for few of those frames. Box-plot at the right side shows the CPU time comparison of the different techniques.
Exp. Data
Energy functional in (5):
specifics of topology preservation
CPU time (seconds)
per frame
Minimum
Discrepancy
(mean) (8)
Average [min max]
of
1) Full-Rank 12.0830 0.1069 0.0016 [-2.8896 3.6884]
2) Full-Rank (14) with 12.0672 [0.0617 1.0096]
3) Full-Rank (14) with 11.3953 [0.9034 1.0023]
4) Full-Rank (14) with 11.0304 [0.9546 1.0059]
5) Low-Rank 3.8771 7.9505e-06 [-1.0144 3.7375]
6) Low-Rank (14) with 3.4895 [0.0213 1.0031]
7) Low-Rank (14) with 3.0214 [0.9184 1.0028]
8) Low-Rank (14) with 2.32170 [0.9950 1.0100]
Table 4: Performance analysis: low-rank vs. full-rank data and their reaction to different degrees of topology preservation (see text for discussion).

The performance of our proposed topology-preserving 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 zoom-in 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.

Figure 7: (A) Resulted transformations, during complex deformations, without applying topology preservation. Highlighted areas denote structure violations that are more clearly displayed in the zoom-in views (A.1). The resulted Jacobian determinant are shown in (A.2). Resulted transformations, after applying topology preservation, are shown in (B) and can be compared with (A) in which (B.1) and (B.2) show that they keep the mesh structures with most of the Jacobian determinant staying at 1.
Figure 8: (From right to left) Accumulated displacement for the apical view of the left ventricle. Few samples of the approximated axial and lateral displacements (top and bottom) are compared against the ground truth. Left side shows the Jacobian determinant of the same sample frames which reflects preservation of the anatomy.

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 time-frames, 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 low-rank and full-rank 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 low-rank 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).

Figure 9: Mean accumulated displacement of the seven segments of the left ventricle. Blue circles make reference to the axial displacement while red squared refer to the lateral displacement.
Figure 10: (Top) Radial and longitudinal strain profiles of the left ventricle. These profiles are evaluated by their sign where negative values reflect shortening and positive ones reflect stretching. (Bottom) Few frames of the cardiac cycle showing the radial strain.

Table 4, Exp. 1 and Exp. 5 show that without topology preservation in the functional (5), the residuum was about with and without the low-rank constraint. Applying low-rank, the Jacobian determinant had higher values but still suffered from heavy distortions. With topology preservation in Exps. 2-4, and 6-8, 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 low-rank representation seems to be the radical decrease in computational time, as seen from a comparison of Exps. 1-4 and Exps. 5-8, where the computational time was decreased by about 75 %.

The residual errors in Exps. 2-4 in the topology-preserved full-rank case were in the order of magnitude resp. . Contrary to that, topology preservation in the low-rank case in Exps. 6-8 yielded minima in the order of magnitude to . Moreover, the discrepancy error showed that low-rank achieved an order of magnitude in comparison with given by the full-rank case. In practice, topology preservation and low-rank 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 Root-Mean 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 well-known 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 low-rank data representation. Together with the better temporal resolution of ultrafast ultrasound, our proposed approach overcame challenges of non-rigid 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 low-rank 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 low-rank 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 AP2012-1943.

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. Lopez-Perez, R. Sebastian, and J. Ferrero, “Three-dimensional cardiac computational modelling: methods, features and applications,” BioMedical Engineering OnLine, pp. 1–31, 2015.
  • [4] M. Prummer, J. Hornegger, G. Lauritsch, L. Wigstrom, E. Girard-Hughes, and R. Fahrig, “Cardiac c-arm 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 phase-based myocardial motion estimation in tagged MRI sequences by a bilinear model and motion compensation,” Medical Image Analysis, pp. 149–162, 2015.
  • [6] J. Royuela-del Val, L. Cordero-Grande, F. Simmross-Wattenberg, M. Martín-Fernández, and C. Alberola-López, “Nonrigid groupwise registration for motion estimation and compensation in compressed sensing reconstruction of breath-hold 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 large-scale organ motion in b-mode 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. Ledesma-Carbayo, J. Kybic, M. Desco, A. Santos, M. Sühling, P. Hunziker, and M. Unser, “Spatio-temporal 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, “Volume-preserving nonrigid registration of MR breast images using free-form 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 three-dimensional 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, “2-d 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, “Three-dimensional cardiac strain imaging in healthy children using rf-data,” Ultrasound in Medicine & Biology, vol. 9, pp. 1399 – 1408, 2011.
  • [23] D. Tenbrinck, S. Schmid, X. Jiang., K. Schafers., and J. Stypmann, “Histogram-based 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 four-dimensional echocardiography,” 2014 IEEE 27th International Symposium on Computer-Based Medical Systems, pp. 149–152, 2014.
  • [25] H. Gao, N. Bijnens, D. Coisne, M. Lugiez, M. Rutten., and J. D’hooge, “2-d left ventricular flow estimation by combining speckle tracking with Navier-Stokes-based 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 b-splines 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 multi-modal ultrasound image registration,” Medical Image Computing and Computer-Assisted Intervention, pp. 363–371, 2015.
  • [28] R. Morin, A. Basarab, S. Bidon, and D. Kouamé, “Motion estimation-based 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 local-to-global 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 Computer-Assisted 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, “High-dimensional image registration using symmetric priors,” 1999.
  • [37] S. Ozere and C. Le Guyader, “Topology preservation for image-registration-related deformation fields,” COMMUN. MATH. SCI., vol. 13, no. 5, pp. 1135–1161, 2015.
  • [38] J. Ashburner and K. J. Friston, “Voxel-based 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, “3-d 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 log-unbiased 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 topology-preserving 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. Sing-Long, 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, “Low-rank spatio-temporal video segmentation,” in Proceedings of the British Machine Vision Conference (BMVC), 2015.
  • [49] B. Haeffele, E. Young, and R. Vidal, “Structured low-rank 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.