Nonrigid Medical Image Registration by Finite-Element Deformable Sheet-Curve Models

Image-based change quantitation has been recognized as a promising tool for accurate assessment of tumor's early response to chemoprevention in cancer research. For example, various changes on breast density and vascularity in glandular tissue are the indicators of early response to treatment. Accurate extraction of glandular tissue from pre- and postcontrast magnetic resonance (MR) images requires a nonrigid registration of sequential MR images embedded with local deformations. This paper reports a newly developed registration method that aligns MR breast images using finite-element deformable sheet-curve models. Specifically, deformable curves are constructed to match the boundaries dynamically, while a deformable sheet of thin-plate splines is designed to model complex local deformations. The experimental results on both digital phantoms and real MR breast images using the new method have been compared to point-based thin-plate-spline (TPS) approach, and have demonstrated a significant and robust improvement in both boundary alignment and local deformation recovery.


INTRODUCTION
Chemoprevention is an established approach to cancer prevention based on 25 years of research [1][2][3]. The breast cancer prevention trial (BCPT) in the United States recently led the Food and Drug Administration (FDA) to approve tamoxifen as a chemopreventive drug to reduce breast cancer risk [4]. Although the overall benefits of tamoxifen outweigh the overall risks, the risks were found significant in the trials, including increased risks of endometrial cancer, problems related to blood clots, and cataracts [5]. Moreover, many clinical studies have also reported that not all women are protected against breast cancer by tamoxifen. Given the risks and the lack of full protection against breast cancer, some high-risk women are reluctant to use tamoxifen because they want to know whether tamoxifen is effective on their breasts, rather than relying on statistical evidence. Hence, it is clinically important to de-velop novel approaches to accurately assess early response to tamoxifen.
Previous studies have shown that estrogens increase breast density [6,7], but the effect on breast density of selective estrogen receptor modulators (SERMs), such as tamoxifen and raloxifene, is unknown. In a previous clinical trial [6], one author of this paper has evaluated the effects of two years of treatment with raloxifene, estrogen, or placebo on breast density using a digitized analysis of mammograms. The main findings of the study can be summarized as the following: (1) the mean breast density was statistically significantly greater in the conjugated estrogens group than it was in the other three groups, and (2) volunteers receiving raloxifene did not increase breast density after two years of treatment. In the study, breast density, as measured by mammography, is determined by the relative amounts of fat and fibroglandular tissue present. As discussed in [6], since estrogen action works through estrogen receptors that are present in glandular tissue, breast density measured on glandular tissue would more directly indicate the response to SERM chemoprevention. However, mammography dose not allow us to accurately measure the density of glandular tissue in practice.
In this study, we use gadolinium-enhanced MR imaging method to differentiate glandular tissue from fibroglandular tissue, and measure the breast density on the glandular tissue. In the light of the antiestrogenic properties of tamoxifen, that is, reduction of circulating growth factors and inhibition of angiogenesis [8], we propose to quantitatively measure changes of breast density and breast vascular density in glandular tissue to assess early response to tamoxifen. With such quantitative measurements, we can determine whether the drug is effective for a patient, early and individually. To compare the glandular tissue volumes between pre-and postcontrast MR images, nonrigid registration should be performed to align two MR images and recover local deformations.
Many nonrigid registration methods have been proposed in recent years [9][10][11], and can be categorized into intensitybased or feature-based approaches [12][13][14][15]. When using gadolinium-enhanced MR images, feature-based approaches offer advantages over intensity-based registration methods since geometric features-points, curves, and surfaces-are consistent between pre-and postcontrast MR images while intensity values are not. By establishing feature correspondences, deformations can be recovered by nonlinear interpolation methods such as the thin-plate-spline (TPS) method [13].
Registration accuracy often depends on the availability of a large number of features. Point-based registration methods typically use anatomical landmark points as the features to register images [13]. Although anatomical landmarks are the most robust features, the number of landmarks is usually limited in MR breast images. Moreover, if local deformations are confined in certain regions, the landmarks have to be well placed over those regions. It is therefore reasonable to use the boundaries of anatomical objects as features to recover local deformations [14,16].
Several approaches have been proposed that used boundaries or curves as matching features [14,[16][17][18]. Moshfeghi proposed an elastic matching method for registration of multimodality images [14]. The method is based on Burr's elastic contour matching algorithm and a Gaussian smoothed deformation model that estimates the deformations over the whole image [19,20]. To reduce the computational complexity, Davatzikos et al. proposed a twostep approach to recover deformations based on boundary mapping, that is, a homothetic mapping of boundaries followed by an elastic deformation transformation (EDT) [16].
However, elastic matching method often suffers from the following problems: (1) since no physical model (or constraint) is imposed on the curve deformation, the resulted curve matching is considered only an ad hoc solution; and (2) the Gaussian smoothed deformation model is an oversimplified model when facing the complex local deformations existed in MR breast images. For example, although Davatzikos' method is a very efficient algorithm resulting in a great reduction in computational complexity, the key assumption behind is homothetic mapping, that is, a mapping by a uniform scaling of its length and arbitrary lengthpreserving bending [16]. Such an assumption cannot hold in MR breast images in that the participated deformations are relatively large, therefore homothetic mapping could produce large curve matching error.
In this paper, we propose deformable sheet-curve models to overcome the problems associated with the existing methods. Deformable curves are used to obtain a reliable matching of the curves using physically constrained deformations, and a deformable sheet with the energy functional of thinplate splines is used to model complex local deformations between the images.

THEORY AND METHOD
In image registration a nonlinear transformation U(x) = (U x (x, y), U y (x, y)) is used to map the first image to the second image. We refer to the function U(x) as the deformation, which can be modeled as either thin-plate splines [13] or elastic sheets [16]. Under the assumption that we have K pairs of corresponding boundaries (or curves), (C 1i , C 2i ), i = 0, . . . , K − 1, we use deformable curves to model the deformation of curves in the two images (denoted as U c ), and a deformable sheet to model the deformation between two images (denoted as U s ). We propose a deformable sheetcurve model shown in Figure 1 to recover the deformation U s through matching the deformable curves and deforming the deformable sheet iteratively.
A deformable sheet-curve model can be presented in variational and finite-element forms. We denote the deformation of the sheet by U s (x, y, t) = (U x s (x, y, t), U y s (x, y, t)), where (x, y) is the bivariate material coordinate and t is the time index. Similarly, the deformation of the curve can be denoted by , where s is the univariate material coordinate of the curve. The strain energy E can be found to characterize the deformable material of the sheet or the curve as an instance of the spline function. At the center of our method, the continuum mechanical equation [21] μ governs the nonrigid motion of the sheet (curve) in response to the extrinsic force f(U), where μ is the mass density function of the deformable sheet (curve) and γ is the viscosity function of the ambient medium. The third term on the lefthand side of the equation is the variational derivative of the strain energy functional E, the internal elastic force of the sheet (curve).
Jianhua Xuan et al. 3 Deformable sheet Deformable curves Figure 1: Configuration of the deformable sheet-curve model.

Variational modeling
The deformable energy for the deformation U s (x, y, t) of the sheet is defined by where the weights w 01 , w 11 , and w 10 control the tensions of the sheet (stretching energy), and w 02 , w 22 , and w 20 control its "rigidities" by weighting the bending energy (see [22] for a detailed discussion on maximal rigidity). Accordingly, the deformable energy for the deformation U c (s, t) of the curve is given by The weight w 1 controls the tension along the curve (stretching energy), and w 2 controls its rigidity (bending energy). The dynamics of the deformable sheet-curve model is defined by the following Lagrangian motion equations: where f s and f i c are the external forces applied on the sheet and the curve, respectively. The registration process is an iterative procedure by repeating the following two steps.
(1) Boundary deformation recovery by deformable curve matching: deformable curves use the corresponding curves in the second image to define the external forces The final deformation U s is obtained when the energies of the deformable sheet-curve model reach their minima.
We define the external forces f i c to reflect the mismatch between the two corresponding curves (C 1i , C 2i ) : is the Gaussian weighted Euclidean distance of each point on the first curve to the nearest point on the second curve. After we recover the curve deformations U i c , we obtain the new curve positions as C new 1i = C 1i + U i c , i = 0, . . . , K − 1. Then we use Gaussian weighted Euclidean distance to define the external forces for the deformable sheet, that is,

Finite element modeling
We use a finite-element method to compute the numerical solutions of U s and U c . We first tessellate the continuous material domain, (x, y) for the sheet and s for the curve, into a mesh of m element subdomains D j . We then approximate U by a weighted sum of continuous basis functions N i (socalled shape functions):  (4) can then be discretized as where M is the mass matrix, C is the damping matrix, K is the stiffness matrix, and F is the forcing matrix. M, C, and F can be obtained as follows: To compute K, we have the following equation: where

Deformable sheet element
The deformable sheet consists of a set of connected triangular elements. Barycentric coordinates are the natural choice for defining shape functions over a triangular domain, since a unifying representation of different triangles can be achieved [23]. Barycentric coordinates (L 1 , L 2 , L 3 ) are defined by the following mapping with material coordinate (x, y): where (x 1 , y 1 ), (x 2 , y 2 ), and (x 3 , y 3 ) are the coordinates of three vertex locations of a triangle. We construct a 9-degree-of-freedom (DOF) triangular element with its position and first parametric partial derivatives at each triangle vertex. The shape functions for the 9-DOF triangle are [24] The triangle's symmetry in barycentric coordinates can be used to generate the shape function for the second and third nodes in terms of the first. To generate N 9 2 use the above equations but add a 1 to each index so that 1→2, 2→3, and 3→1. The N 9 3 functions are generated by adding another 1 to each index. By defining 9 DOFs of the triangular element as a T = [U 1 , U 1x , U 1y , U 2 , U 2x , U 2y , U 3 , U 3x , U 3y ], we can approximate the U s as U h s = N 9 a = N 9 1 N 9 2 N 9 3 a.

Deformable curve element
The finite element of the curve has 4 DOFs for the two nodes located atthe ends of the curve segment. The DOFs at each node are the position and tangent of the node. The deformation of the curve segment can be approximated as the weighted sum of a set of Hermite polynomials [23]: where U i c , i = 0, . . . , 3, are nodal variables and N i , i = 0, . . . , 3, are given as follows: where h is the parametric element length.

Numerical integration
The deformable sheet-curve model can be stabilized during the matching process if its motion is critically damped to minimize vibrations. Critical damping can be achieved by appropriately balancing the mass and damping distributions. A simple way of eliminating vibration while preserving useful dynamics is to set the mass density in (1) to zero, thus reducing (5) to This first-order dynamic system governs the model that has no inertia and comes to rest as soon as all the forces balance. We integrate (14) using an explicit first-order Euler method [24]. The method begins with a simple forward difference approximation. By considering extrapolation from time with the forward difference approximation (14) then becomes Jianhua Xuan et al.

5
Thus we obtain the updating formula as follows: It is well known that finite difference methods for initialvalue systems yield expressions very similar to the above results obtained by finite-element schemes [25]. A noteworthy distinction is that the coefficient matrix for finite differencing is diagonal in the usual difference approximation. This leads, in the forward difference approximation, to an efficient algorithm for solving the problem. In the finite-element method, due to the sparseness of C, it is difficult to compute C −1 since complicated singular decomposition methods have to be used. Fortunately, a physical solution exists in computational mechanics, namely, the "lumping" procedure, to overcome such difficulties. The idea can be physically interpreted as replacing the continuous material with the distributed mass, that is, concentrated material with "lumped" mass ("beads") at the nodes. In practice, there are several ways to perform such a "lumping" procedure-for example, using modified shape functions or different numerical integral methods [24]. Among those, the easiest way is to keep only the diagonal coefficients and is adopted here in solving the motion equations of the deformable sheet-curve model.

EXPERIMENTAL RESULTS
Deformable sheet-curve models were implemented using a finite-element method using 9-DOF triangular sheet elements and 4-DOF curve elements. The deformable curves are used to match boundaries, and the deformable sheet is modeled as a thin-plate spline to recover the local deformation. The registration process is an iterative procedure by repeating the following two steps: (1) the deformations of the curves (U i c ) are obtained by solving Lagrangian motion equations of the curves with the external forces (f i c ) defined by the corresponding curves, and (2) the deformation of the sheet (U s ) is obtained by solving Lagrangian motion equation of the deformable sheet with the external forces (f s ) defined by the new curve positions (derived from U i c ). The final deformation of the whole image is recovered when the energies of the deformable sheet-curve model reach their minima. In our experiments, we applied the deformable sheet-curve models to both simulated data and real breast MR images to demonstrate the performance of the method in recovering boundary deformation and image deformation.

Simulation result
Digital phantoms are generated and used to validate the proposed boundary-based registration approach in a controlled manner with known truth. As shown in Figure 2, the boundaries in the first image are depicted in solid line while those in the second image are in dashed line. In our approach, we use deformable curves to model the boundaries and match them by recovering the deformation iteratively. The recovered boundary deformation demonstrates an excellent matching between the initial and target boundaries with a mean square error (MSE) less than 2%, where the parame- ters (i.e., the weights in (3)) are experimentally chosen to obtain a good overall matching performance. The deformable sheet works cooperatively with deformable curves to infer the image deformation iteratively. As we can see from Figure 2, the recovered image deformation is constrained by the thinplate-spline functional. The digital phantom study therefore provides compelling evidence that deformable sheet-curve models can achieve accurate boundary matching and complex image deformation modeling.

Registration of breast MR images
We have further applied our deformable sheet-curve model to register real contrast-enhanced sequential MR breast images to assess early response to chemoprevention. The data set of this study consists of serial studies in 29 women at high risk of breast cancer who received either tamoxifen or no drug. MR images are acquired with contrast agent (gadolinium) at 0, 3, and 6 months along with biopsies for histopathology, RNA, and gene analysis at Lombardi Cancer Center (LCC) of the Georgetown University Medical Center. We use fibroglandular tissue as the reference object to register pre-and postcontrast MR breast images, see Figure 4 for an example. The regions of fibroglandular tissue are extracted using a stochastic model-based segmentation method based on the standard finite normal model (SFNM) and a contextual Bayesian relaxation labeling (CBRL) procedure [3]. The segmentation result is shown in Figure 5 with all the regions of fibroglandular tissue identified in both preand postcontrast images. Then all the corner points on the boundaries of fibroglandular tissue are extracted (see Figure 5). By establishing the corner point correspondence, we can partition the boundaries of fibroglandular tissue into corresponding curves in pre-and postcontrast images. These curves serve as the features in the registration of pre-and postcontrast MR breast images. The deformable curves with 4-DOF finite elements are used to model the extracted boundaries of fibroglandular tissue, and the deformable sheet with 9-DOF finite elements is then used to model image deformation. The registration process is a cooperative procedure by repeating boundary matching and image deformation inferring. The deformable curve matching gives us the boundary deformation, and the recovered boundary deformation is then used to deform the deformable sheet to derive image deformation. In our finite-element implementation, we use numerical integration method described in Section 2 to solve the Lagrangian motion equations. The final deformation of the image is recovered when the energies of the deformable sheet-curve model reach their minima. Figure 6 shows the extracted regions of glandular tissue by subtracting the registered precontrast image from the postcontrast image.
As a comparison, we show the extracted regions of glandular tissue using the registration result of the point-based TPS method in Figure 7. As we can see, many false regions of glandular tissue are resulted from the misalignment of the boundaries of fibroglandular tissue due to the limited number of the available feature points (25 feature points in this experiment). The boundary misalignment can be clearly seen in Figure 7(b), where two regions of fibroglandular tissue are registered by the point-based TPS method. Figure 7(a) shows the boundary registration result using our deformable sheetcurve model, which demonstrates a great improvement in boundary alignment compared with that of the point-based TPS method. The experimental results show that our deformable sheet-curve model can provide an expected improvement in registration accuracy regarding both boundary alignment and local deformation recovery.

CONCLUSION
We have developed deformable sheet-curve models for nonrigid medical image registration. The registration method supported by our deformable sheet-curve models is a feature-based approach in that (1) deformable curves are used to model curve features-the boundaries of objects, and (2) the thin-plate-spline sheet is used to model the local deformation over the image. The deformable sheetcurve model is also a dynamic system that is governed by Lagrangian motion equations. The model has been implemented in finite-element method with 9-DOF triangular sheet elements and 4-DOF curve elements. By solving Lagrangian motion equations iteratively, we can recover the deformation over the whole image. We have applied the deformable sheet-curve model to register MR breast images for the assessment of early response to chemoprevention. The experimental results have demonstrated a significant improvement in the registration accuracy with reliable boundary alignment and accurate local deformation recovery. We believe that our comparative studies provide useful information on the utility of the proposed method for nonrigid image registration. Given the difficulty of the task, while the optimality of the method may be data or modality dependent, we would expect it to be an important tool in change detection across temporal image sequences.

ACKNOWLEDGMENT
This work was supported in part by the Department of Defense under Grant DAMD17-03-0448.