Physical Constraint Finite Element Model for Medical Image Registration

Due to being derived from linear assumption, most elastic body based non-rigid image registration algorithms are facing challenges for soft tissues with complex nonlinear behavior and with large deformations. To take into account the geometric nonlinearity of soft tissues, we propose a registration algorithm on the basis of Newtonian differential equation. The material behavior of soft tissues is modeled as St. Venant-Kirchhoff elasticity, and the nonlinearity of the continuum represents the quadratic term of the deformation gradient under the Green- St.Venant strain. In our algorithm, the elastic force is formulated as the derivative of the deformation energy with respect to the nodal displacement vectors of the finite element; the external force is determined by the registration similarity gradient flow which drives the floating image deforming to the equilibrium condition. We compared our approach to three other models: 1) the conventional linear elastic finite element model (FEM); 2) the dynamic elastic FEM; 3) the robust block matching (RBM) method. The registration accuracy was measured using three similarities: MSD (Mean Square Difference), NC (Normalized Correlation) and NMI (Normalized Mutual Information), and was also measured using the mean and max distance between the ground seeds and corresponding ones after registration. We validated our method on 60 image pairs including 30 medical image pairs with artificial deformation and 30 clinical image pairs for both the chest chemotherapy treatment in different periods and brain MRI normalization. Our method achieved a distance error of 0.320±0.138 mm in x direction and 0.326±0.111 mm in y direction, MSD of 41.96±13.74, NC of 0.9958±0.0019, NMI of 1.2962±0.0114 for images with large artificial deformations; and average NC of 0.9622±0.008 and NMI of 1.2764±0.0089 for the real clinical cases. Student’s t-test demonstrated that our model statistically outperformed the other methods in comparison (p-values <0.05).


Introduction
Widely available biomedical images are essential for more accurate diagnosis and management of patients with a variety of diseases. Registration of these medical images into a spatial correspondence and alignment is important for fusion and maximization of the underlying information in the imaging datasets. Image registration is achieved by estimating an optimal transformation or displacement field [1]. It is one of the most fundamental research areas for various clinical applications, including planning and delivery of radiotherapy treatment [2][3], image-guided surgery [4][5], cross modality image fusion [6], morphometric study [7] and treatment monitoring [8].
After more than three decades of intensive research, various registration methods have been proposed [9][10][11]. Among these methods, the intensity based registration has the advantage of direct usage of the image intensity information without necessary preprocessing of segmentation and extraction of features such as salient points, curves and surfaces as in the feature based registration. However, compared with feature based counterpart, these intensity based methods normally face the challenges of reducing high computational cost and avoiding local optimization.
Basically, two types of geometric transformations are used in the intensity based methods which are derived from interpolation theory or derived from physical models. Free-form deformation (FFD) method [12] and local affine transformation (LAT) method [13] are two most common types of interpolation strategies. In FFD method, a rectangular grid of control points is defined in order to determine the deformation. Displacements between control points are propagated by interpolation. This method requires few degrees of freedom to describe local deformations and is able to efficiently provide smooth results. However, the dependence on a regular grid of control points restricts their adaptability and it is difficult to change control points topology. Many extensions of FFD methods have been proposed [14][15]. LAT methods parameterize the transformation by locally linear deformation. The images are hierarchically partitioned into contiguous blocks and an affine transformation is recovered for each one of them. However, the underlying image content has not been considered in the splitting process, and the regularization providing for the global smoothness in LAT method may lead to the ambiguous matching [16].
The diffusion model and the elastic body model in the framework of finite element method are the currently recognized physical models for medical image registration. The diffusion model upon the optical flow (OF) constraint describes deformations by assuming a constant brightness constraint of floating voxels [17][18]. However, it lacked a sound theoretical justification [9]. In finite element elastic model, the deformation of the image was modeled as an elastic body that is described by the Navier-Cauchy partial differential equation [19]. The gradient of the similarity measure is used as an external force field which tries to deform the floating image to fit the reference image configuration, while the internal elastic forces of the solid oppose the deformation. Thus, the floating image is deformed until the internal and external forces reach an equilibrium state. On the basis of a linear assumption [20], the governing equation can be simplified and the nodal displacement matrix can be obtained under the nodal external force matrix in the framework of finite element analysis method. However, the linear assumption may be not valid for soft tissues with complex nonlinear behavior and cannot accurately simulate large organ deformations [21].
To address the challenging issue in the elastic registration, in this paper, we propose a new registration method based on the St. Venant-Kirchhoff (VK) model that is the simplest hyperelastic material model for representing the material behavior of soft tissues. Furthermore, we extend this model in a dynamic fashion by using the Newtonian differential equation [22][23].
The geometrical nonlinearity of the continuum is taken into account by considering the quadratic term of the deformation gradient in the strain tensor when deriving the elastic force. For a more efficient implementation of our proposed realistic model, we adopt a hierarchical (also referred to as a "global-to-local") strategy in our registration algorithm.

Methods and Materials
In the proposed deformable registration algorithm, the displacement field is obtained subject to external constraints imposed by the similarity metric of the image pair for registration. The equilibrium of our model is achieved by the minimization of the global energy functional that consists of elastic deformation energy and the energy enforced by the external image force.

Governing Equation for the Nonlinear Elastic Model
In the framework of the finite element analysis method, the displacement of any image pixel can be interpolated with the nodal displacement and the shape function. On the basis of an isotropic linear assumption, a matrix equation describing the motion of the elastic model can be expressed as: where the stiffness matrix K is obtained from the mesh and elastic parameters.R is the nodal external force matrix obtained by optimizing the spatially encoded mutual information (SEMI), and u is the displacement vector of element nodes. The linear elastic finite element model (LFEM) has an important limitation when coping with large deformations. To account for large deformations, we propose a nonlinear elastic model governed by the second order differential equation as below: where M and C are the matrices of the mass and damping, respectively. These two matrices are interrelated to each other with C = αM where α is a scale factor. The mass matrix M is chosen to be a diagonal matrix with diagonal elements as M el ii ¼ 1 3 rV el , where V el is the volume of element el, and ρ is the mass density.
The elastic force matrixF constrains the solution space and together with the external force, this elastic force guides optimization process to converge. Under the linear elastic condition,F ¼ KU according to Eq (1), thus Eq (2) is deduced to the linear elastic dynamic finite element model (DFEM) as below: In our proposed nonlinear finite element model (NFEM), the elastic force matrixF is derived from the elastic deformation energy and innovatively models the geometrical nonlinearity based on hyper-elastic material properties. Fig 1 illustrates the components governing equations for LFEM, DFEM and NFEM, and elaborates the differences of these three FEM registration methods.
The nonlinear governing Eq (2) is solved with an explicit integration scheme [24], where the displacement vector at time t+1 is computed from the elastic force and external force estimated at time t as follows: where 4t is the integration time step.

Modeling of Elasticity via Green-St.Venant Strain Measure
The hyperelastic material model has been proven appropriate to represent the material behavior of soft tissues and has been successfully applied in surgery simulation [25][26]. We used the simplest isotropic hyperelastic material model as proposed by Wittek et al. [27]. In this paper, we derive the elastic forceF to model tissue deformation on the basis of St.Venant-Kirchhoff (VK) material model [28]. The VK material is the simplest hyper-elastic material and holds nonlinear property between the strain and the displacement gradient, and in the meantime obeys Hooke's law in constitutive relationships.
The deformation of the image pair for registration is assumed to be resulted from a large displacement such as the flexible plates or shell deformed with external forces. The strain relative to the initial (base) configure which is usually called the Green-St.Venant Strain and is defined as where J r CG is the right Cauchy-Green strain tensor, I is an identity tensor and J X is the deformation gradient. Physically, the right Cauchy-Green strain tensor leads to the square of local change in distances due to deformation. By substituting J r CG in Eq (6) to Eq (5), we have and with the assumption that |ru|<<1, the quadratic term ru T ru can be neglected, which is encountered in the linear elastic finite element model analysis.
The components of e can be rearranged as a 3-component strain vector ε for the finite element analysis of our elastic model as follows: where where h m is a 4×1 vector, H m is a 4×4 sparse symmetric matrix, and g is the displacement gradient vector, i.e., where u x and u y are the displacement in x and y directions. Assume an elastic continuum without initial stresses of strains, the elastic deformation energy W e of an elastic body can be expressed as: where O is the image domain, σ is the stress tensor which relates to the Green-St.Venant Strain in the following manner: Eq (14) can be rewritten in a more compact form as σ = Dε. For the material with isotropic characteristics, matrix D can be determined with the Young's modulus E and Poisson's ratio v as follows: By substituting Eq (8) and Eq (14) to Eq (13), the elastic energy can be rewritten as: where W mn is usually called the strain energy density and can be expressed as follows: where L is a differential operator given in Eq (12) and D mn is the element of matrix D.

Elastic Force
According to the Delaunay criteria [29], uniform meshes that are composed of triangular elements O el are generated for both the reference image I R and the floating image I T . If the displacement vectors u i el (i = 1,2,3) of the triangular nodes are obtained, the displacement vectors u(x) of other points in the triangular element can be obtained approximately by the following interpolations: where N i el (x)(i = 1,2,3) is the shape function of an element. Using linear interpolation, the shape functions are defined as the so-called natural coordinates L i of the element. i.e., With natural coordinates, the area A el and the coefficients a i el , b i el , c i el are all defined as following: where (x i , y i ) i = 1,2,3 is the global coordinates of element nodes, and the other coefficients are found by cyclic interchange of the indexes. Given the expression of the elastic energy, the elastic forceF el i acting on node i of element el can be derived as: and by substituting Eq (18) to Eq (17), @W mn @u el i can be written as: where B el i ¼ LN el i , and L is a differential operator given in Eq (12).F el as determined in Eq (21) can be assembled to the globalF in the image domain O.

External Force
In the framework of FEM, the nodal external force matrixR can be assembled by the external forceR el of all element el, which is defined as follow: where R(x) is the external force field subjected to the gradient flows of the matching metric, N el is the shape function vector assembled by N i el (x)(i = 1,2,3).
Our aim of registration is to find a deformation field so that the registration similarity, whch is represented as spatially encoded mutual information (SEMI) [30], can be maximized. Mathematically, the matching metric increases faster along its gradient direction. Therefore, a plausible choice of the external force would be the gradient flows of the matching metric which would then drive the tissue to deform in a trend increases the matching similarity.
The gradient flows of the similarity metric between reference image r(x) and floating image m(x) is defined as follow: where R is the intensity value domain and G s (x, x 0 ) is the weighting function for spatial information encoding, i.e., x 0 )dx is the normalization factor and ϕ(Á,Á) is the Parzen intensity kernel, i.e.,  where where β is the smoothing parameter for the Parzen estimates, p s (r ' ) and p s (m ' ) are the region marginal intensity probability density of the reference image and floating image respectively.

Hierarchical Registration Strategy
For a more efficient implementation of our proposed registration method for images with large deformation, a hierarchical strategy (often referred to as the "global-to-local" strategy) was adopted in our registration procedure. A schematic illustration of the block diagram of our hierarchical registration method is shown in Fig 2. The global registration is to find the transformation matrix maximizing the similarity metric that is mutual information in our paper [31]. After the global registration, a subsequent local non-rigid registration was performed based on our proposed NFEM. We constructed the coarse-to-fine registration pyramid by firstly performing the filtering of the images with a Gaussian low pass filter for noise reduction and then setting uniform mesh nodes using the predefined element size to generate the mesh pyramid.

Parameters Setting
We set three types of parameters in our method. The first one is the elastic material related parameters. For simplicity, we considered our model as a homogeneous isotropic hyperelastic model and we fixed the elastic modulus E = 100kPa and the Poisson's ratio v = 0.45 as in [32]. The second type is the parameters related to the optimization procedure including the number of hierarchical levels, the element sizes and the time step for the explicit integration scheme. In our experiments, the spatial resolution of the images is 256×256. Upon pyramid generation, excessively larger number of nodes would increase computational requirement while excessively smaller number of nodes would lead to unacceptable loss of information. In our experiments, meshes with 64×64, and 32×32 nodes were used to construct a two-level multiresolution pyramid for better balance of computational efficiency and registration accuracy.
Since the explicit integration scheme for the nonlinear elastic model is only conditionally stable, a thorough analysis about it has been presented by determining a suitable time step. Similar to the procedure as in [33], we chose 4t = 0.004, the value of the damping coefficient and the mass density were deduced according to the predefined spectral radius ρ Ã (ρ Ã = 0.99) [33]. These parameters are summarized in Table 1.
The third type of parameters is for SEMI metric including the size of the local region denoted as the diameter d, the standard deviation γ and the smoothing parameter β for the Parzen estimates. As analyzed in [30], there was a tradeoff between the global robustness and the local accuracy. The large local region contributes more information for estimating the joint intensity probability density, while the small local region contributes the better accuracy in the local region. For obtaining optimized values for the region diameter d, analysis of registration accuracy over 20 abdominal CT image pairs for different d values was conducted.

7 Materials
Evaluations of our nonlinear finite element model (NFEM) were performed on 30 clinical CT image pairs with artificial deformations. These abdominal and chest CT images were acquired from the First Affiliated Hospital of Soochow University and used as the reference images. 10 B-Spline based small synthetic deformation fields and 10 large synthetic deformation fields were used to warp the reference image to generate floating images for the nonrigid registration. In order to compare our method as well as the LFEM and the DFEM methods with the RBM method, 10 different non-uniform artificial fields were used to generate ten floating images, instead of the B-Spline based predefined fields used above. The deformation exceeding 15 pixels in local regions was considered as "large" deformation and was included to validate our proposed NFEM algorithm.
We also evaluated the performance of our proposed algorithm for clinical images including 10 chest CT image pairs and 20 brain MRI image pairs. The chest CT image pairs were acquired from two different treatment periods of 10 patients in the First Affiliated Hospital of Soochow University. The brain MRI images from different healthy volunteers were collected and made available by the CASILab at The University of North Carolina at Chapel Hill and distributed by the MIDAS Data Server at Kitware, Inc [34]. We randomly selected 21 T1-Flash Images and chose one T1-Flash image as the reference and the rest 20 T1-Flash images as floating images.

Ethics Statement
The study about CT images was carried out according to the Helsinki Declaration and approved by the ethical committee of the First Affiliated Hospital of Soochow University. The need for informed consent was waived, because the data used in this study had already been collected for clinical purposes. Furthermore, the present study did not interfere with the treatment of patients and the database was organized in a way that makes the identification of an individual patient impossible. The study about MR images was carried out according to the Helsinki Declaration and approved by the ethical committee of the University of North Carolina at Chapel Hill. All subjects provided signed consent allowing images to be made publicly available on the website. The database was organized in a way that makes the identification of an individual patient impossible.

Registration Accuracy Assessment
To evaluate the proposed NFEM algorithm for non-rigid registration, we compared it with 1) the conventional LFEM registration algorithm; 2) DFEM registration algorithm; and 3) the robust block matching (RBM) based registration software which was called "Nifty_Reg" and developed at University College London containing programs to perform rigid, affine and the RBM algorithms for nonlinear registration [35]. Evaluation was carried out by visual assessment in subtraction images between the reference image and the floating image after registration and also by three popular measures: the mean square difference (MSD), the normalized correlation (NC) and the normalized mutual information (NMI). Furthermore, for the registration of the images with artificial deformations, 40 random foreground points were automatically selected from the reference images. Then the mean and max distance of the corresponding points in the image pairs before and after registration was calculated to measure the registration accuracy. Fig 3 show the impact of parameter d on the average values of NC and NMI. We found d = 15 contributed the highest NC and NMI values. The standard deviation γ was set to a small value, typically γ = 0.25 Ã d.  Fig 4D and 4E show the color images of the distance of a point in the reference to its counterpart in the floating image to display the deformation degree in local regions.

Registration Results of Medical Images with Uniform Deformations
The visual assessment and comparison of the registration results by observing the subregions in the squares were performed as shown in Fig 5. The more corresponding and similar the registration results and the reference images are, the more capable the algorithm of tackling large deformations. As demonstrated in Fig 5D, our registration result was more corresponding to the reference image. Table 2 shows the quantitative evaluations of different methods on the ten pairs of images with large uniform deformations in terms of the mean distance (mean_x and mean_y) and maximal distance (max_x and max_y) between true and recovered deformation in the x and y direction. The mean values of the mean and max distances of our method were the lowest when compared to the mean values from other three methods, which indicated that our method achieved the highest registration accuracy. Student's t-test (p-values <0.05) validated this improvement was statistically significant.
The comparisons of our proposed NFEM method with other three methods in terms of MSD, NC and NMI (shown in Table 3) demonstrated that the mean values of NC and NMI using our NFEM was greater than the corresponding values using other three methods, and consistently, the mean value of MSD using NFEM was the lowest. Student's t-test (pvalues <0.05) consistently validated our method statistically outperformed its counterpart comparison methods. The quantitative validation on ten pairs of images with small uniform deformations is displayed in Table 4. Table 4 summarized the mean distances and different metrics for 10 image pairs with small deformation and demonstrated that the RBM achieved the best registration results followed by our NFEM. In comparison, although there were no obvious matching errors in the registration results shown in Fig 6, LFEM and DFEM algorithms slightly fell behind regarding the registration accuracy because of the infinitesimal deformation assumption.    Table 5. In addition, the Student's t-test (p-values <0.05) shown in Table 5 indicated that our NFEM method significantly outperformed other three methods. The visual assessment displayed as the subtraction images in Fig 7 further consolidated the same conclusion. The statistical performance in terms of NC and NMI for different methods in both chest CT and brain MRI cases is shown in Fig 9 where these two matching metrics are illustrated by means of statistic box-plots. The t-tests were implemented and results are illustrated in Table 6 and Table 7, which indicated that, at the 0.05 level, the mean value for NC and NMI of all registration validation using method NFEM is significantly greater than the ones using Global, LFEM, DFEM, RBM method.

Registration Results of Real Medical Images
All evaluations were implemented with C on Windows 7 operating system and performed on a DELL desktop with Intel(R) Core(TM) i7-4770 @ 3.4GHz CPU. Using the aforementioned two-level pyramid, the total freedom degree for computing the displacement field was 10240. The average computation time using LFEM was 3 minutes and 35 seconds and the one using DFEM was 4 minutes and 3 seconds. The average computation time using our NFEM method was 5 minutes and 26 seconds.

Discussion
Complicated nonlinear characteristics of the soft tissues or the tremendous deformations introduced by the internal organ movements or treatment interventions pose significant challenges for the current non-rigid image registration algorithms. Underlying linear hypotheses of these algorithms are one of the foremost reasons to underpin the requirement for further improvement of these nonlinear methods. We took into account the geometric nonlinearity of elastic body model. In our algorithm, the elastic force was modeled as the derivative of the deformation energy with respect to the nodal displacement vectors of the finite element; the external force was derived from the registration similarity gradient flow. We validated our algorithm firstly on image pairs with artificial deformations and then on real CT and MRI images. We compared our NFEM method with three other methods: 1) the conventional LFEM registration algorithm; 2) DFEM registration algorithm; and 3) the robust block matching (RBM) algorithm.
To validate the performance of our algorithm when registering images with large deformations, we conducted experiments on ten pairs of images with large artificial deformations. The visual assessment and comparison of the registration results were performed by observing the sub-regions in the squares as shown in Fig 5. As demonstrated in Fig 5D, our registration result was more corresponding to the reference image, which indicated that our proposed NFEM algorithm outperformed other three algorithms. The better registration results from our algorithm were mainly due to the contribution of the internal elastic force that took into account the geometrical nonlinearity to compensate for the local large displacements.
The quantitative evaluations of different methods (as shown in Table 2 and Table 3) on the ten pairs of images with large uniform deformations consistently demonstrated that our method statistically outperformed its counterpart comparison methods. These better statistical results from our algorithm also indicated that the nonlinear internal elastic force played the positive role in the registration.
The quantitative validation on ten pairs of images with small uniform deformations as displayed in Table 4 demonstrated that the RBM achieved the best registration results followed by our NFEM. The better accuracy from RBM method benefited from the fact that RBM was a B-Spline based method and it was therefore inherently more suitable for handling the small  uniform deformations generated by B-Spline interpolation. In comparison, although there were no obvious matching errors in the registration results shown in Fig 6, LFEM and DFEM algorithms slightly fell behind regarding the registration accuracy because of the infinitesimal deformation assumption. The registration results in Table 2 and Table 4 demonstrated that our method was more robust for the cases with both large deformations and small deformations. Meanwhile, 10 corresponding floating images were generated by pre-defined non-uniform artificial fields to investigate the feasibility of our method for estimating large non-uniform deformations. The performance of different methods was compared in terms of MSD, NC and NMI. The Student's t-test (p-values <0.05) shown in Table 5 indicated that our NFEM method significantly outperformed other three methods. The visual assessment displayed as the subtraction images in Fig 7 further consolidated the same conclusion. The nonlinear elasticity modeling via Green-St.Venant strain measure contributed to the improved deformation estimation and hence led to the improved registration accuracy for the images with non-uniform deformations.
In addition to the evaluations of our proposed algorithm on artificially deformed images, we validated our algorithm on clinical CT and MRI images as shown in Fig 8. The statistic boxplots illustrated in Fig 9 demonstrated that our method was able to effectively compensate the deformations. Statistical analysis as illustrated in Table 6 and Table 7 indicated that, at the 0.05 level, the mean value for NC and NMI of all ten registration validation using method NFEM is significantly greater than the ones using Global, LFEM, DFEM, RBM method, which demonstrated that our proposed NFEM algorithm outperformed other methods for registering real clinical image pairs and would be potentially usable for assessing disease progression or patient response to treatment. While our method achieved improved registration accuracy compared to other linear elastic registration models, it has a tradeoff in regard to the computation time. Theoretically, explicit (b) NMI parameters for box-plots. In CT image registration cases, these methods are represented by appending a suffix 1 (e.g. Global1) to the name abbreviations, while in MRI image cases, the suffix 2 is added to their name abbreviations (e.g. Global2).
doi:10.1371/journal.pone.0140567.g009 integration in our NFEM method should contribute to improve the computational efficiency. However, updating the relatively complicated elastic force in the iterations was computationally costly. The convergence towards the steady-state solution would be increased if the integration time step could be adaptively updated during iterations.

Conclusion
In this paper, a novel physics-based nonlinear registration algorithm is proposed based on the FEM elastic body model in the framework of hierarchical strategy. The novelty of this algorithm is the inclusion of nonlinearity between the strain and the displacement gradient based on St.Venant-Kirchhoff model. This nonlinear model is represented as a second order differential equation to obtain equilibrium between the internal and the external forces. Experimental validation on images with artificial deformations and real medical images demonstrated that our nonlinear FEM elastic model outperformed the traditional LFEM model, DFEM method, and RBM method in terms of the quantitative and qualitative analysis.