A Novel 3D-to-3D Diffeomorphic Registration Algorithm With Applications to Left Ventricle Segmentation in MR and Ultrasound Sequences

Left ventricular segmentation is a difficult and time-consuming task performed by clinicians, requiring the use of manual contours. We propose a novel three-dimensional diffeomorphic registration algorithm for endocardial segmentation of the left ventricle in magnetic resonance images and ultrasound temporal sequences. The proposed diffeomorphic registration method computes a voxel-to-voxel correspondence in three-dimensional space and is parameterized by one radial and three curl components to emulate the cardiac deformations. In addition, the proposed method allows for enforcing constraints to control the amount of deformation. The method was evaluated on 521 temporal frames from 20 patients from the Automated Cardiac Diagnosis Challenge magnetic resonance imaging dataset and 213 frames from 10 patients undergoing ultrasound scans from the Mazankowski Alberta Heart Institute. The algorithm was compared against six other registration methods, the Symmetric Normalization algorithm from the Dipy package, two variants of the Demons algorithm from the Insight Toolkit software package, two versions of RealTiTracker, and Elastix. The proposed method yielded overall Dice scores of 98.10 (0.90)% for the MRI dataset and 92.90 (2.42)% for the ultrasound dataset. The robustness of the algorithm is demonstrated by the high performance on multiple imaging modalities and various patient abnormalities.


I. INTRODUCTION
Assesment of the left ventricle (LV) is frequently performed for a wide variety of disorders, including cardiomyopathy, hypertrophy, myocardial infarction, and congenital diseases [1]. Magnetic resonance imaging (MRI) and ultrasound (US) imaging are frequently used to evaluate LV function, where MRI is often considered the gold standard for diagnosis. A number of studies have shown that US can provide similar results to MRI [2] with advantages including low The associate editor coordinating the review of this manuscript and approving it for publication was Jinhua Sheng . cost, portability, and high spatial and temporal resolution. Analysis of the geometry and function of the LV over the cardiac cycle provides key indicators of many diseases. For instance, the volume of the ventricle can be used to calculate the ejection fraction, a measure of the efficiency of the heart at pumping blood. Other common diagnosis measures include dimensions of the chamber and global and longitudinal strain.

A. MRI REGISTRATION
Several approaches have been developed for the registration of 3D cardiac MRI volumes. These can be divided into those that are non-deep learning-based and those using deep learning techniques to compute the 3D displacement fields. One approach not based on deep learning describes a method to create spatially and temporally smooth displacement fields in 4D [3], [4]. To compute a set of displacement fields temporally, it is common to perform pairwise image registration, but in this manner, errors are propagated through the sequence. In order to register starting from the end-diastolic (ED) frame, authors create a second 4D sequence of only the ED frames, and registered this to the original temporal sequence. To perform the image matching, an attribute vector is formed for each voxel, which includes information concerning intensity, boundary as well as geometric moment invariants, calculated for various neighborhood sizes. A selection of the most distinct attribute vectors are matched together between the two sequences of volumes in a hierarchical fashion.
One set of authors in their work developed a deformable registration algorithm technique termed hyperelastic warping [5], [6], [7] which has the ability to measure strain in a set of medical images. The method performs image registration resulting in a diffeomorphic transform, and incorporates the use of fiber structure definitions and material properties. Authors have used the technique to calculate the strain during diastole [5] as well as during the systolic phase [8].
Others have developed methods to jointly perform segmentation and registration for delineation of the blood pools of the left and right ventricles and the myocardium [9]. The authors combine a level-set segmentation along with registration of two volumes, where multiple anatomical regions are able to be segmented. They demonstrate two uses of their method -temporal intra-patient segmentation and interpatient registration from an atlas.
A number of deep learning methods for the computation of displacement fields in a 3D-to-3D manner have been developed. The ability to create reference deformations needed to train a deep learning network is difficult. Authors [10] created a unique way to generate reference deformations from regions of interest that were previously segmented. Therefore the neural network they have developed learns to understand more global information instead of matching the intensity of voxels as in classic registration algorithms. Another supervised approach was developed by [11] where the authors developed a multi-resolution convolutional neural network that outputs a 3D displacement field taking as input two volumes to be registered. They introduce the idea of a smoothness constraint by use of the norm of the Jacobian matrix.
Other authors have developed unsupervised techniques for the formation of the displacement fields. One method developed by [12] is an unsupervised approach based on the use of a conditional variational autoencoder. Within this paradigm, constraints are applied to the transformation in order to enforce symmetry and to be diffeomorphic. Another set of authors developed an unsupervised technique where the network learns the transformation parameters that are used to form displacements. The displacements are then used for resampling of the moving image and this warped image is compared to the fixed image [13]. This method was extended [14] to use mutual information as a loss function instead of the typical mean squares and normalized cross-correlation metrics. Other methods that have been proposed calculate displacement fields over the entire cardiac cycle [15], where the authors use a multi-resolution patchbased approach based on the U-Net architecture [16].
As slice misalignment is an issue with 3D cine MR volumes, in one proposed method the authors [17] first employed a correction to align each of the 2D slices in a single frame based on [18]. The approach uses a regression-based method using convolutional neural networks to predict the center of the blood pool. After the slice correction, the authors used the VoxelMorph approach [19] to produce the deformation fields.
Other deep learning approaches that have been developed are based on free-form deformation. For instance, one set of authors [20] developed a group-wise registration approach where deformation is represented by a 1D convolution cascade. Authors perform registration and optimization of all frames simultaneously to a reference image, where the introduction of the 1D convolutional cascade greatly decreases the execution time.

B. ULTRASOUND REGISTRATION
There have been several approaches proposed in the literature for where 3D registration of the LV over a temporal sequence is used for computing the deformation fields [21] for US sequences. These can be divided into areas including intensity-based, regularization model-based approaches, and also deep learning methods. Authors [22] used spatiotemporal registration based on an elastic approach to measure the strain of the LV. A classic multi-resolution approach with a B-spline transformation model and mutual information as a similarity metric was used for the registration. Others have tried to use the envelope-detected image, instead of the processed ultrasound data, to perform registration in the spherical coordinate system [23]. One drawback of using B-spline transformation models as in [22] is that a rectangular grid is used in Cartesian space for the control points, which may not be appropriate for capturing the motion of the heart. In this formulation by [24], criteria for spatial smoothness are not enforced in the radial direction, which is necessary for cardiac images. Therefore basis functions are chosen that allow for the appropriate representation of the heart anatomy. This approach has been expanded by [25] in order to include the constraint of volume conservation in the anatomical free-form deformation model. Temporal diffeomorphic free-form deformation approaches have been developed [26] to perform an assessment of the motion and the strain of the LV. These have been extended by numerous authors by the introduction of new similarity metrics [27], [28]. The concept of sparsity has also been introduced [29] in order to represent deformations that are locally discontinuous. VOLUME 11, 2023 One deep learning approach for registration and segmentation of echocardiography sequences performs both tracking of the LV motion as well as segmentation of the chamber [30]. The authors proposed the use of an unsupervised 3D U-Net neural network architecture for obtaining the motion, similar to the method of VoxelMorph [31]. A 3D U-Net type architecture is again used for the segmentation network in a weakly supervised fashion [32]. The use of incompressibility constraints is incorporated within the combination of these networks in order to produce anatomically feasible displacement fields.
A deep learning technique for performing both motion and strain estimation has been developed by [21]. The authors first use a multi-layer perceptron in a supervised manner in order to learn the features between the input and ground truth deformation fields. The authors used three approaches to produce the starting estimates of the fields, a radio-frequency block matching method, flow network tracking, and a free form deformation approach for nonrigid registration. The displacement fields utilized were 4D, therefore including both spatial and temporal patterns. The method was further extended by the use of biomechanical constraints in order to produce deformation fields that were anatomically possible.

C. PUBLICLY AVAILABLE REGISTRATION ALGORITHMS
There are a number of registration algorithms that are publicly available and have been applied to a variety of image modalities and patient characteristics. Dipy's Symmetric Normalization (SyN) algorithm [33] is one widely used algorithm that ensures a diffeomorphic transformation is generated. The algorithm has been applied to a significant number of applications, most notably brain registration [34], [35]. Insight ToolKit (ITK) is a commonly used open-source toolbox for analysis of medical images, ranging from various filters to segmentation and registration algorithms [36]. One set of deformable registration algorithms in ITK includes an implementation of variants of the Demons algorithm, including the classic Demons algorithm and an extension to include fast symmetric forces. This set of algorithms has been used in many areas, such as head and neck CBCT registration [37] and registration of CT images of the prostate.
RealTiTracker [38], [39] is another registration toolbox that implements two types of optical flow based approaches. For the conservation of intensity term, the L 1 and L 2 norm have been implemented. It has been adapted for use in MR-based registration of the brain between T1 and T2 volumes [40], and for other applications such as for radiotherapy dosimetry estimation [41]. The software package Elastix [42], [43] consists of a multitude of transformation models, optimization methods and cost functions. It has been applied by a variety of researchers and used as a standard benchmark, such as for capturing the deformation between MR T1 brain volumes [34]. The wide use of the above algorithms and the public availability of the methods allows us to compare the proposed algorithm to the above approaches.

D. PROPOSED METHOD
We propose a semi-automated segmentation algorithm based on a moving mesh approach [44] to compute the spatial transformations across a cardiac cycle. The user provides manual segmentations of the LV for the end-diastolic and end-systolic frames, and the proposed algorithm produces a voxel-to-voxel correspondence for each of the frames in the sequence. After obtaining the segmentation of the LV for each frame, further analysis can include volumetric, regional, and strain computations, which are crucial for diagnosis of the patient [45].
The proposed algorithm has a number of advantages. First, the method performs a 3D-to-3D registration, as opposed to previous related methods which have only captured 2D motion [46], [47], [48]. This is a significant improvement as 3D-to-3D registration is crucial in capturing the true motion of the heart, as 2D is unable to detect out of plane motion. This is especially important for cardiac imaging as the motion of the heart consists of rotational, radial and torsional movement. The ability of the algorithm to represent the deformation field in terms of a single radial and three rotational components in 3D space makes it appropriate for analyzing cardiac data.
Similar to other deformable registration algorithms, the proposed method ensures that topology is preserved and that grid lines of the same family do not fold over each other. This is crucial as image registration should be able to capture realistic and accurate deformations for the cardiac chamber. The method also allows for the control of the amount of deformation, by setting explicit constraints on the determinant of the transformation Jacobian. In order to demonstrate these advantages, an experiment was performed using the diffeomorphic constraints, which demonstrates how modifying these values influence the presence of mesh folding. The proposed method also does not require the use of an explicit regularization term with its associated weight in the cost function.
The proposed method is evaluated on two datasets, one containing cardiac MRI sequences, and the other consisting of US sequences. The first dataset, the Automated Cardiac Diagnosis Challenge (ACDC) dataset [1], consists of MRI sequences from 100 patients who suffered from previous myocardial infarctions, cardiomyopathy, etc. The proposed algorithm was evaluated on a subset of 521 temporal frames acquired from 20 patients from the cohort of 100. In the second dataset, the algorithm was evaluated over 213 temporal 3D US volumes acquired from 10 patients. The proposed algorithm was then compared against six publicly available registration algorithms, the Symmetric Normalization diffeomorphic registration package from Dipy [33], two variants of the Demons algorithm from Insight Toolkit (ITK) [36], two variants of RealTiTracker [38], [39], and the Elastix package [42], [43].
The proposed method implements a 3D-to-3D diffeomorphic registration method, which requires a number of 3146 VOLUME 11, 2023 significant changes and improvements from previous work described by [46], [47], and [48]. 1) The divergence and curl operations were formulated for use in 3D 2) The Fourier transform based Poisson equation solver was developed for the 3D case 3). The div-curl system was significantly augmented by the development of the 3D formulation 4) The proposed method is evaluated on both MRI and US temporal sequences consisting of patients of varying cardiac conditions 5) The diffeomorphic registration method was compared to six other registration methods 6) Experiments were performed to examine the role of the diffeomorphic constraints.

II. METHODOLOGY
The proposed registration method is based on a moving mesh grid generation approach [44]. Instead of using traditional grid displacements to describe the deformation fields, they are parameterized using the transformation Jacobian determinant and the curl of end velocity field.

A. THEORETICAL OVERVIEW
The method computes the voxel-to-voxel correspondences for a series of 3D volumes in a temporal sequence from the nth image T n to the T n+1 image defined over ⊂ R 3 . The general registration problem is stated as an optimization of a similarity measure [49], where the similarity measure is defined as the squared L 2 norm: where φ : → represents the transformation function and ξ ∈ represents the voxel locations. This problem does not have a unique solution and hence requires more constraints, where the goal is to find a permissible deformation field. We define this deformation field using the monitor function µ and the curl of end velocity field γ . We define a continuous monitor function µ(ξ ): The purpose of registration is to obtain a transformation φ : → and ∂ → ∂ so that: where the Jacobian determinant of the transformation φ is given by J φ . In order to compute the transformation φ, the following steps are used.
Step 1: A vector field ρ(ξ ) is defined such that: Step 2: From the vector field ρ(ξ ), a velocity vector field is generated. t is artificially introduced time and is in the range [0,1]: To obtain φ the transformation, the ordinary differential equation below is solved, where t is again within the range [0,1] and ψ(ξ,t = 0) = ξ : If φ equal to ψ and is evaluated at t = 1 results in φ(ξ ) = ψ(ξ, t = 1).
As the problem above might have multiple solutions, a constraint is added to the curl of the vector field ρ(ξ ). The Dirichlet boundary condition is used, and an intermediate vector field ρ(ξ ) is solved from the div-curl system given below. By formulating it in such a manner, a unique solution can be formed.
The radial and rotational components of the transformation field are represented by µ and γ respectively. In order to have a unique solution, constraints need to be added to the vector field µ(ξ ). The problem can thus be reformulated as: (8) which is a constrained optimization problem using the parameterization from (1). τ lb and τ ub represent the lower and upper bounds for the allowable deformation, where τ lb > 0. This ensures that φ is a diffeomorphism. The first constraint enforces the fact that the domains before and after the transformation are equal while the second constraint ensures that the amount of compressibility is limited.
A step-then-correct optimization can be used to solve the above problem efficiently. Figure 1 displays a flowchart of the proposed registration method and Algorithm 1 displays the detailed steps. Given two 3D volumes to register, an iterative approach is taken to calculate the deformation field. First the gradients of both µ and γ are calculated. If the iteration number is not less than the maximum and the step size is not below a specified threshold, the optimization continues. The gradients are first updated and constraints are imposed for each voxel location to make sure the transformation is diffeomorphic. The vector field is then computed and the deformation field formed. The cost is then computed by using the resulting deformation field to warp the moving volume, which is then compared to the template volume by means of the similarity metric. If the current cost is less than the previous cost value, optimization continues by recomputing the gradients and the iteration number increases. If not, the step size is first decreased before continuing with the optimization, hence the name step-then-correct optimization.

B. NUMERICAL IMPLEMENTATION
The div-curl system given in Equation (7) can be written in the following manner, where Equation (9) provides the div-curl system in 3D. The divergence represents the radial motion, VOLUME 11, 2023 FIGURE 1. Flowchart displaying the components of the proposed 3D-to-3D diffeomorphic image registration algorithm. and the curl, comprised of three terms along each of the x, y, and z axes represent the rotational motion.
The derivatives for each of the f 1 , f 2 , f 3 , and f 4 terms are computed with respect to x, y, and z. Three Poisson equations of the form are given below, where the relevant derivative terms have been combined from (9), where f i k = df i dk , with i = 1,2,. . . ,4 and k = x,y,z: Under the null boundary condition, these Poisson equations can be efficiently solved by using the fast Fourier transform (FFT) based Poisson solver.

C. REGISTRATION FOR A SEQUENCE OF VOLUMES
The formulation previously detailed describes the registration between two volumes. In order to obtain the displacement fields over the entire sequence, registration is performed sequentially over the cardiac cycle. Since registration is performed sequentially over the temporal sequence in a pairwise manner, it may be prone to errors that can accumulate over

Algorithm 1
Step-Then-Correct Optimization Given two 3D volumes, comprised of the fixed volume T n and the moving volume T n+1 , the following steps are computed in order to calculate the deformation field φ: Step 1: Compute the gradients of µ and γ , which are given by ∇µ(T n ,T n+1 ,φ) and ∇γ (T n ,T n+1 ,φ) while δ > δ th and i < max_iter do Step 2: Update gradients: Step 3: Impose constraints from (8) for each voxel location ξ ∈ : Step 4: Compute a vector field ρ(ξ ) that satisfies (7) and compute the deformation field φ. Step Step 1 of recomputing the gradients else Decrease step size δ Start from Step 2 of updating the gradients end if end while time. A strategy was therefore devised to overcome these tracking errors, where registration is performed in both the forward and reverse directions and the displacement fields are combined using a weighting function: where DM f n,n+1 represents the forward displacement vectors and DM b n,n+1 the reverse displacement vectors between the n and n + 1th frames, and N represents the total number of frames in the sequence. This formulation assigns a higher weight to the displacement field computed using forward registration and less to the displacement field computed using reverse registration if the frame number is closer to the beginning of the sequence.

D. IMPLEMENTATION NOTES
The proposed algorithm was implemented using Python 3.6.2 with PyTorch version 1.7.1 using a NVIDIA Tesla P100 GPU. The other registration algorithms were run using Python 3.6.2 on a NVIDIA Tesla C2075 graphics card. Visualization Toolkit (VTK) version 8.1.2 was used for creation of the meshes and Paraview 5.7.0 was used to create the figures displaying the difference in the meshes between the ground truth and the registration methods. The average time to perform registration for a single frame was 0.59 seconds for the ACDC dataset, and 2.69 seconds for the Mazankowski dataset.

III. METHODOLOGY FOR EXPERIMENTAL RESULTS
The proposed registration method was evaluated on two patient groups, 20 patients from the Automated Cardiac Diagnosis Competition (ACDC) dataset consisting of MR scans [1], and ten patients who underwent US scans at the Mazankowski Alberta Heart Institute. Registration was performed temporally across the cardiac cycle in the forward and reverse directions. In order to compute the segmentation of the LV for each frame of the cardiac cycle, the ground truth ED and ES meshes representing the LV endocardium are used as input. Using the forward and reverse displacement fields, and the ground truth ED and ES meshes, a weighted average was computed and applied to form a single mesh for each of the remaining frames in the cardiac cycle.

A. PARAMETER SELECTION
To reduce the amount of time spent during the registration process, a bounding box was defined around the ED segmentation with a margin of 10 voxels. The maximum number of iterations was set to 20, and to allow for large deformations in the cardiac tissue, the values of the transformation Jacobian determinant τ lb and τ ub were set to 0.25 and 6.0 respectively. The other parameters include the step size h for the Runge-Kutta method was set to 0.05 (20 steps). For the step-thencorrect optimization method, the initial value of the step size δ was set to 0.5, and the factor to reduce δ was set to 0.66. The minimum δ threshold δ th was set to 0.0001. The values for the parameters were chosen based on the previous publication for the 2D-to-2D diffeomorphic registration algorithm development [47]. Slight modifications were made to the Jacobian determinant τ lb and τ ub values, which were modified from 0.1 to 0.25 and 4.0 to 6.0 respectively, where τ ub was increased in order to allow for larger tissue deformations.

B. ALGORITHMS USED FOR COMPARISON
The proposed registration method was compared to six standard registration methods, Dipy's Symmetric Normalization (SyN) algorithm [33], two variants of the Demons algorithm from ITK [36], two versions of RealTiTracker [38], [39], and Elastix [42], [43]. These algorithms were chosen for their ease of use and public availability. The Symmetric Normalization (SyN) algorithm from Dipy [33] ensures that the transformation is diffeomorphic, where the function and the inverse are both smooth. The algorithm extends the formulation described in [50], which used a Lagrangian diffeomorphic registration method. Symmetry between the two images is guaranteed in this method, where the invertibility constraints are used directly during the optimization process. The sum of squared differences (SSD) was used as the similarity metric. The method performs multi-resolution registration, using a Gaussian pyramid of three levels. The three levels of the Gaussian pyramid used 10, 10, and 5 iterations for each level respectively. For the Demons algorithm, two variants were employed from ITK to compare against the proposed method, first the classical Demons algorithm, and secondly a variant that incorporates fast symmetric forces. For both, histogram matching was first used to ensure that the 3D volumes to register were similar. 50 iterations were used for the first implementation, and 200 for the second approach. A Gaussian kernel was used to smooth the displacement field, with a value of 1.0 for the standard deviation for both approaches. Parameters were chosen placed on example codes provided in the documentation.
The RealTiTracker software package consists of a set of algorithms based on optical flow methods [51]. One term in optical flow ensures that the relative movement between the two frames is small, while another constraint enforces that the intensity difference between subsequent frames is conserved. For the conservation of intensity term, there are two forms that differ in how the spatial and temporal derivatives are employed. For the L 2 L 1 functional the L 1 norm is used, and for the L 2 L 2 functional the L 2 norm is used. One parameter, the weighting term alpha, determines how sensitive to gray level intensity variation the motion is. The value was set to 0.4 for the L 2 L 1 functional, and to 0.1 for L 2 L 2 . These were adapted from example codes and found to be reasonable values in terms of compromising for precision versus accuracy. Elastix [42], [43] is a software package that is comprised of a set of transformation models, cost functions and various optimization methods. The parameters used consisted of the following: B-spline transformation with a multi-resolution approach of five levels, advanced Mattes mutual information [52] for the similarity metric and adaptive stochastic gradient descent for the optimizer. The parameters were chosen based on the default values suggested for multiresolution B-spline registration.

C. METRICS
The proposed algorithm computes the displacements over the cardiac cycle using both forward and reverse registration VOLUME 11, 2023 and combined using a weighting function. Therefore, the same approach was used for each of the six other registration methods. For the validation of the methods, a set of standard distance and overlap measures were employed. The metrics were calculated for all frames except for the ED and ES frames, as these were used as input to the algorithm. The mean absolute distance (MAD) is the mean distance between each point in the ground truth mesh S and the closest point in the mesh from the proposed algorithm T , where the result is given in mm [53]. The Hausdorff distance (HD) measures the local maximum distance between the ground truth mesh S and the mesh from the proposed algorithm T , where the result is reported in mm [54]. The HD is defined in (13), where d(s, t) is the Euclidean distance between two points of each of the meshes.
The Dice metric (14) is a measure for overlap between two binary volumes V s and V t , with 1 indicating complete overlap, and 0 indicating no overlap between the regions [55].
The determinant of the Jacobian provides information about the local transformation at a point. Therefore in order to determine if mesh folding occurs, the percentage of voxels with a determinant of Jacobian less than zero is reported, given by J <0 %.

A. ACDC DATASET
The ACDC public dataset consists of 100 patients that underwent cine MR imaging at the University Hospital of Dijon, France [1]. A selection of 20 patients was chosen from the dataset for evaluation of the proposed algorithm. Scans were acquired using a 1.5T scanner (Siemens Arena, Siemens Medical Solutions, Germany) and a 3T scanner (Siemens Trio Tim, Siemens Medical Solutions, Germany). The volumes ranged from 184 × 216 × 8 to 256 × 256 × 11 voxels. The resolutions of the voxels ranged from 1.367 × 1.367 × 10.0 mm to 1.875 × 1.875 × 10.0 mm. The number of frames spanned from 13 to 35 frames within the cardiac cycle, with a total of 521 frames over the 20 patients. One clinical expert delineated the ground truth for the ED and ES frames. The set of patients consisted of normal subjects, those with a previous myocardial infarction, myopathy, and those with an abnormal right ventricle.
The ACDC dataset only includes the ground truth for the ED and ES frames, but the proposed method is designed to be evaluated over a sequence of frames in the cardiac cycle. In order to reduce the amount of time it would take an expert to manually annotate each short-axis slice for each frame and patient, an automated method was used. For each short-axis slice sequence, the method of [48] was used to obtain a set of ground truth contours using the provided ED and ES phases.
The contours were reviewed and edited when necessary by an expert radiologist. A subset of 20 out of the 100 patients were used for the validation of the proposed algorithm. This number was chosen as the contour assessment performed by the radiologist was a time consuming process. The first 20 patients were chosen from the ACDC dataset of 100.
In order to compare the ground truth contour points to the other registration methods, the 3D points were converted to a mesh. For each frame, spline interpolation was used to ensure each contour (for each short axis slice) had N = 241 points. Contours for each short axis slice were then aligned to ensure the start points were the same. Faces for the mesh were created by automatically creating triangular faces between adjacent pairs of short axis slice contours. The mesh for each frame was then used for computing metrics for the quantitative performance evaluation. In order to calculate the volume of the mesh for each frame, the method of [56] was implemented. Table 1 displays the evaluation of the proposed method and the six other registration algorithms using the mean absolute distance d m , Hausdorff distance d H , the Dice score Dice, and the percentage of voxels with a determinant of Jacobian less than zero, J <0 %. The robust performance can be seen for the proposed registration method, as it performs better than the others in terms of the d m and the d H metrics. Elastix performs slightly better than the proposed method in terms of only the Dice score, with values of 98.21 versus 98.10 respectively. The J <0 % value for the proposed algorithm is equal to zero, indicating that no mesh folding has occurred. It can be observed that for Elastix, the J <0 % is above zero, therefore displaying some mesh folding. Kruskal-Wallis H significance tests were performed for the d m , d H and Dice metrics (for each frame) in order to determine if there is a significant difference among the proposed method and the other registration methods. It can be seen from the table that the p values calculated are less than 0.0001, indicating that there exists a significant difference among the registration methods. The average time to perform the registration for a single frame was 0.59 seconds for the ACDC dataset.

b: VOLUME ANALYSIS
Examining the volume of the left ventricle over the cardiac cycle is a crucial part of diagnosing the patient. One method of analyzing the difference in the volumes between the ground truth and the proposed method, along with the other registration methods is to use a Bland-Altman plot. The Bland-Altman plot provides information about the bias, or the mean of the differences between the ground truth and the registration method. Limits of agreement, set to 1.96 times the standard deviation of the differences provides information about the span of the differences. Figure 2 provides  the Bland-Altman plot for the proposed registration method compared to the ground truth volumes. The red line indicates the bias, while the two yellow lines indicate the limits of agreement. The black line is the zero reference line, indicating that the volumes from the two methods match. It can be seen that the points are well clustered around the reference line, with few outliers beyond the two limits of agreement.
The bias and limits of agreement were also calculated for the six other registration methods in Table 2. It can be observed that the bias for Demons, 0.35 mL, and Elastix, 0.92 mL are both slightly smaller than the bias for the proposed method at 1.77 mL, indicating they are more accurate in terms of only volume estimation. The values for the bias for all methods is positive, indicating a slight underestimation of the volume.
The volume curves can also be analyzed, as they provide the clinician with information about the systolic and diastolic function of the heart. Figure 3 displays an example of volume curves for one patient from the ACDC MRI dataset. The ground truth is shown in red-filled circles while the proposed method is displayed in neon green squares. For both cases, it can be seen that the proposed method closely follows the  ground truth over the entire cardiac cycle. It can be observed that the Elastix method also closely follows the ground truth, but that towards the end of the diastolic phase, the proposed method has a higher accuracy. Figure 4 gives a visual representation of the difference between the ground truth mesh and the proposed method and VOLUME 11, 2023 the other registration methods for a frame in the systolic phase in the first row, and a frame in the diastolic phase in the second row. The areas highlighted in red indicate a larger difference between the ground truth (in mm), and the areas in blue indicate an almost complete overlap with the ground truth. The ground truth mesh is shown in gray. For the proposed method, it can be seen that the mesh closely follows the true ground truth mesh.

b: BINARY MASK COMPARISON
An alternate method of comparing the ground truth to the proposed segmentation is to view the overlap of the binary masks. Figure 5 displays the underlying MRI image overlaid with the ground truth segmentation in green and the mask from the proposed segmentation in red. Results are displayed for slices at the apex, mid-cavity and the base of the left ventricle for a frame in the systolic phase and a frame in the diastolic phase. It can be seen that overall there is a large overlap between the proposed segmentation and the ground truth.

B. ULTRASOUND DATASET
Subjects were approved to be scanned by the human research ethics committee at the University of Alberta. Using an X5-1 transducer, subjects were scanned at the Mazankowski Alberta Heart Institute (Edmonton, Alberta, Canada) by an expert sonographer on a Philips iE33 US machine (Philips Healthcare, Best, Netherlands). To achieve a volume rate of greater than 20 volumes per second, a 3D sector angle of 70 × 80 degrees was used. The set of subjects underwent US scans to assess and diagnose LV function. The dataset consists of subjects with the number of frames ranging from 16 to 26 for a total of 213 frames, with volumes ranging from 224 × 176 × 208 to 256 × 176 × 208 voxels. The resolutions of the voxels ranged from 0.617 × 0.787 × 0.533 mm to 0.810 × 1.005 × 0.681 mm. The TomTec Arena software (TomTec Imaging Systems, Unterschleissheim, Germany) was used by an expert cardiologist to provide the ground truth segmentation meshes and volume information for the left ventricle. In order to calculate the volume of the mesh for each frame, the method of [56] was implemented. Table 3 reports the evaluation using the d m , d H , Dice, and the J <0 % metric for the proposed method and the other registration algorithms. The proposed method performs better than the other registration methods in terms of the d m and the Dice metric, where for the d H metric, Elastix performs slightly better. Kruskal-Wallis H significance tests were performed for the d m , d H and Dice metrics (for each frame) in order to determine if there is a significant difference among the proposed method and the other registration methods. It can be seen from the table that the p values calculated are less than 0.0001, indicating that there exists a significant difference among the registration methods. The average time to perform registration for a single frame was 2.69 seconds for the Mazankowski dataset.

b: VOLUME ANALYSIS
A Bland-Atlman analysis was performed for the US dataset from the Mazankowski as seen in Figure 6. The red line indicates the bias, while the two yellow lines indicate the limits of agreement. The black line is the 0 reference line, indicating that the volumes from the two methods match. From the plot it can be seen that many of the points are clustered around the reference line, but that some outliers do exist beyond the two limits of agreement.
The Bland-Altman analysis was also performed for the other registration methods, with results concerning the bias and limits of agreement given in Table 4. It can be observed that the proposed method yields the smallest bias of 0.73 mL, indicating robust performance in terms of volume compared to the ground truth. The positive values indicate the underestimation of the volume, while negative values indicate the overestimation. Figure 7 displays an example of volume curves for one patient from the US dataset. The ground truth is shown in red-filled circles while the proposed method is displayed in neon green squares. It can be seen that the proposed method closely follows the ground truth over the entire cardiac cycle. In the volume curve for the patient from the Mazankowski, all of the methods perform well during systole, but it can be seen that for diastole there is high agreement between the proposed method and the ground truth. Figure 8 gives a visual representation of the difference between the ground truth mesh and the proposed method and the other software packages for registration for the systolic phase in the first row, and for the diastolic phase in the second row. The areas shown in red indicate a larger difference between the ground truth (in mm), and the areas highlighted in blue indicate an almost complete overlap. The ground truth mesh is shown in gray. For the proposed method, it can be seen that the mesh closely follow the true ground truth mesh.

b: BINARY MASK COMPARISON
An alternate method of comparing the ground truth to the proposed segmentation is to view the overlap of the binary masks. Figure 9 displays the underlying MRI image overlaid with the ground truth segmentation in green and the mask from the proposed segmentation in red. Results are displayed for slices at the apex, mid-cavity and the base of the left ventricle for a frame in the systolic phase and a frame in the diastolic phase. It can be seen that overall there is a large overlap between the proposed segmentation and the ground truth.

C. DIFFEOMORPHIC CONSTRAINTS EXPERIMENT
One of the main advantages of the proposed registration algorithm is that it forces the transformation to be diffeomorphic. This in turn enables the generation of anatomically plausible deformation within the cardiac tissue. An experiment was therefore performed using the ACDC MRI and Mazankowski US datasets in order to test the effect of the diffeomorphic constraints on the presence of mesh folding. Table 5 displays the values of the percentage of voxels with a determinant of Jacobian less than 0. It can be seen that including the diffeomorphic constraints results in no voxels with a determinant less than 0, indicating no mesh folding. When the constraints are removed, a small percentage of voxels exhibit mesh folding. Figure 10 displays the determinant of the Jacobian with and without the diffeomorphic constraints for the two datasets. The first rows displays figures from the MRI dataset, while the second row displays figures from the US dataset. The  TABLE 3. Quantitative evaluation results: The mean absolute difference d m in mm, Hausdorff distance d H in mm, Dice score Dice, and the percentage of voxels with a determinant of Jacobian less than zero, J <0 %, are provided for quantitative evaluation of the accuracy between the ground truth segmentation and the proposed segmentation method, Dipy SyN, ITK Demons (classic and fast symmetric forces), RealTiTracker (L 2 L 1 and L 2 L 2 ) and Elastix over 10 subjects using US data from the Mazankowski Alberta Heart Institute. The lower the values of d m , d H , and the higher the Dice, the more accurate the segmentation. Values above zero for J <0 % indicate the presence of mesh folding. Kruskal-Wallis H significance tests were performed for the d m , d H and the Dice metrics to determine if there exists a significant difference among the proposed method and the other registration methods. The average of the time it takes to perform registration of one frame is provided in seconds. Standard deviation values for all metrics are provided in parentheses. Values highlighted in bold indicate the metric with the best performance.  first column indicates the warped image, the middle column displays the determinant of Jacobian with the diffeomorphic constraints, and the last column displays the determinant of the Jacobian with no diffeomorphic constraints. Areas of interest are annotated with an arrow pointing to a circle. In these areas, it can be seen for the cases without the constraints, some discontinuities exist. This is evident by the black pixels (pixels with a determinant of Jacobian less than or equal to zero).

V. DISCUSSION
In this study, we have demonstrated the development of a 3D-to-3D registration algorithm using a moving mesh correspondence method [44]. Previous methods were developed only for 2D-to-2D image registration [46], [47], [48]. Applications of the algorithm include the delineation of the endocardial borders for the left ventricle in MRI and US datasets, where the algorithm was validated on a subset of publicly available data from the ACDC Challenge [1], and a set of US sequences from patients scanned at the Mazankowski Alberta Heart Institute. The algorithm was compared to six other available registration algorithms, Symmetric Normalization diffeomorphic registration package from Dipy [33], two variants of the Demons algorithm from ITK [36], two versions of RealTiTracker [38], [39], and the Elastix software package [42], [43].
The proposed registration algorithm is diffeomorphic, which does not allow for folding of the mesh. This is especially crucial when describing the actual tissue movement of the heart. Comparing the percentage of voxels with  a Jacobian determinant less than zero, it can be seen that all of the methods produced a value greater than zero for either the MRI or US datasets. Therefore, these methods may not be able to capture the true cardiac deformation of the tissue. The diffeomorphic algorithm is also unique in the fact that the optimization of the deformation field is performed on the radial and rotational representation, making it suitable for cardiac analysis.
There are several additional advantages of the proposed method compared to other algorithms for LV segmentation. One advantage is that the method does not require the use of a manually created training set, which could be difficult VOLUME 11, 2023 TABLE 5. Diffeomorphic constraint experiment results: For the proposed registration method, the diffeomorphic constraints were removed in order to see the effect on the determinant of the Jacobian. The percentage of voxels with a determinant of the Jacobian less than zero are reported J <0 %, indicating mesh folding. The mean absolute difference d m in mm, Hausdorff distance d H in mm, Dice score Dice are provided.

FIGURE 10.
Examples from a patient from the MRI dataset (first row) and a patient from the US data (second row). The first column shows the warped image, the second column displays the determinant of Jacobian map from using the diffeomorphic constraints, and the third column displays the determinant of Jacobian map without using the constraints. Areas of interest are annotated by an arrow pointing towards a black circle, where discontinuities exist in the case where no diffeomorphic constraints are used. Pixels with a determinant of Jacobian less than zero are black.
to obtain depending on the type of patients analyzed. A second advantage is that the proposed algorithm does not make geometric assumptions about the shape of the chamber; this is crucial if the algorithm is to capture markedly different anatomical differences compared to a normal patient, for instance in the case of patients that have cardiomyopathy, myocardial infarction or congenital heart disease. The algorithm has also been tested on MRI as well as US sequences, making it robust to the imaging modality.
There are a number of potential issues concerning the proposed algorithm. The voxel-to-voxel correspondence mapping may be affected by the presence of large abnormalities, artifacts, and other sources of noise. A second drawback is the method of obtaining the gold standard for both the MRI and US datasets. For the US dataset, a single expert cardiologist provided the ground truth annotations using the TomTec software. For the MRI dataset from the ACDC challenge, the ground truth was only provided at the end diastolic and end systolic frames. To generate the gold standard for the other frames in the cardiac cycle, the method of [48] was used for each 2D slice from apex to base, to propagate the contours temporally across the cardiac cycle. An expert radiologist modified the sets of contours to account for any possible issues. Another drawback of the proposed approach is the use of manual segmentations for both the end-diastolic and endsystolic frames.

VI. CONCLUSION
A diffeomorphic 3D-to-3D registration has been developed by the authors for use in MRI and US temporal sequences. The algorithm models the deformation field by use of the radial and rotational components, making it suitable for modeling plausible cardiac motion. We have demonstrated the use of the registration algorithm for segmentation of the endocardium in the left ventricle. As input the method requires two meshes delineating the left ventricle at enddiastole and end-systole. The image registration algorithm is applied in a sequential manner for a temporal sequence, where both forward and reverse deformation fields are computed and weighted to compute the final deformation field. Applying the deformation fields to the meshes results in a segmentation for each frame of the cardiac cycle.
The algorithm has been validated on two modalities, the publicly available ACDC MRI dataset and a set of US patients acquired from the Mazankowski Alberta Heart Institute. The proposed method was compared to a set of standard medical image registration software packages, the proposed algorithm yielded Dice scores of 98.10 (0.90)% and 92.90 (2.42)% on the MRI and US datasets respectively. These results demonstrate the robustness of the algorithm on a variety of datasets and patients.
A number of additions could be included to further increase the viability of the proposed algorithm. Preprocessing could be performed to reduce the influence of speckle noise and improve tracking. The ground truth segmentations could also be improved, as they could be created by manual contouring instead. Additionally, it would be beneficial in the future to obtain annotations of the endocardium by more than one clinical expert in order to perform inter-observer and intraobserver studies. In the future, the algorithm will be improved to only use the end-diastolic frame. In order to assess strain of the myocardium, some clinicians may prefer to annotate the myocardium instead of the endocardium and the epicardium. It would be useful to test this approach for obtaining the myocardial segmentation.
The sum of square differences as the distance metric was chosen as the registration is performed over a sequence of images obtained from the same modality. In the future, other similarity metrics can be used, such as mutual information if multi-modal registration is required. A variable incompressibility constraint could also be incorporated into the registration framework. This would avoid both the mesh folding issue and would allow for anatomically plausible myocardium deformations. In particular, the values of J lb and J ub can be set to [1-eps,1 + eps], respectively within the myocardium, and a wider range outside of the myocardium.
Future applications of the registration algorithm include the segmentation of other chambers of the heart, regional and strain analysis, as well as computed tomography or multi-modality registration. In addition to the endocardium, it would be useful to validate the algorithm on the epicardial and myocardial borders.