Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
IEEE J Sel Top Signal Process. Author manuscript; available in PMC 2010 December 20.
Published in final edited form as:
IEEE J Sel Top Signal Process. 2009 February 1; 3(1): 159–169.
doi:  10.1109/JSTSP.2008.2011116
PMCID: PMC3004236

A simple regularizer for B-spline nonrigid image registration that encourages local invertibility

Se Young Chun, Student Member, IEEE and Jeffrey A. Fessler, Fellow, IEEE


Nonrigid image registration is an important task for many medical imaging applications. In particular, for radiation oncology it is desirable to track respiratory motion for thoracic cancer treatment. B-splines are convenient for modeling nonrigid deformations, but ensuring invertibility can be a challenge. This paper describes sufficient conditions for local invertibility of deformations based on B-spline bases. These sufficient conditions can be used with constrained optimization to enforce local invertibility. We also incorporate these conditions into nonrigid image registration methods based on a simple penalty approach that encourages diffeomorphic deformations. Traditional Jacobian penalty methods penalize negative Jacobian determinant values only at grid points. In contrast, our new method enforces a sufficient condition for invertibility directly on the deformation coefficients to encourage invertibility globally over a 3D continuous domain. The proposed penalty approach requires substantially less compute time than Jacobian penalties per iteration.

Index Terms: B-splines, nonrigid image registration, penalty method, local invertibility sufficient conditions, thorax CT images

I. Introduction

Image registration is a core tool in many medical imaging applications, including fusion of structural and functional images. Several image reconstruction schemes for MR, PET and CT incorporate motion correction or joint estimation of motion into the reconstruction process to improve image quality [1], [2], [3], [4], [5], [6]. Radiation treatments may be able to target cancer cells more accurately through motion correction [7]. Rigid or affine transformations can provide fast image registration. However, most of the human body does not conform to rigid or affine approximations [8]. Lamare et al. [3] used affine image registration for respiratory motion correction, but reported that it was sufficient only for a single organ and associated lesions. Effective motion correction usually requires nonrigid image registration, which enables more flexible matching of local details between two images than rigid registration.

There are many methods for nonrigid image registration [8], [9]. B-spline bases are used frequently for nonrigid image registration [10] because locally supported basis function expansions are convenient computationally and B-splines have the properties of smoothness, compact support, fast interpolation schemes and hierarchical structure for multi-resolution [9], [11], [12]. However deformations with high degrees of freedom can lead to unrealistic transformation results such as folding in the absence of appropriate constraints [8].

There have been some efforts to regularize nonrigid image registration based on B-splines by making certain reasonable assumptions. Rueckert et al. [13] penalized the bending energy of the deformation directly, assuming that the local deformation of tissues should be smooth. Sorzano et al. [14] proposed a regularizer based on the gradients of the divergence and the curl of the displacement field. Rohlfing et al. [15] used an incompressibility constraint: the Jacobian determinant of a transformation should be unity, assuming that local deformations are volume preserving. They applied this method after making an initial affine transformation. See [9] for other methods for constraining the transformation.

Another reasonable constraint is to impose local invertibility of the nonrigid transformation to ensure that image registration is topology-preserving or diffeomorphic.

One way to ensure local invertibility is to require the Jacobian determinant of the transformation to be positive everywhere, either as a hard constraint or by a penalty method [16]. However most such approaches constrain the Jacobian determinant of a transformation only at each discrete voxel grid point, so local invertibility is not strictly guaranteed on the whole continuous domain. Recently, Sdika [17] described a condition involving the gradient of the Jacobian determinant that encourages the local invertibility to be achieved everywhere even though that condition is invoked only at each discrete grid point. However, compared to unregularized image registration, calculating the Jacobian determinant or its gradient significantly increases computation time due to additional B-spline interpolations of the partial derivatives of a deformation.

Ensuring invertibility is somewhat easier when using 1st-order B-spline bases for deformations. Musse et al. [18] derived elegant linear constraints that provide necessary and sufficient conditions to ensure that the Jacobian determinant values of such transformations are positive everywhere. However, that 2D approach was restricted to 1st-order B-spline deformations. Karacali et al. [19] proposed a method to regularize 2D and 3D deformations to ensure that 1st-order B-splines are topology-preserving. Noblet et al. [20] generalized [18] for 3D B-spline deformations and illustrated their method with 1st-order B-splines, but enforcing the constraints requires much higher computation than regularization based on bending energy.

Lastly, one can ensure local invertibility by imposing sufficient conditions that are simpler than the necessary conditions. Choi et al. [21] suggested box constraints for cubic B-spline deformation coefficients that ensure invertibility, but those sufficient conditions preclude large deformations. Rueckert et al. [22] concatenated many transformations based on those box constraints to achieve large deformations. Rohde et al. [23] suggested a sufficient condition for local invertibility, derived using Neuman series for a transformation model that uses a sum of deformations. Motivated by [23], Kim et al. [24], [25], [26] suggested similar sufficient conditions for 3D transformations based on cubic B-splines and implemented a constrained minimization algorithm using Dykstra’s cyclic projection method. We recently extended Kim’s sufficient conditions for local invertibility of deformations so that we can use nth-order B-spline bases and so we can also assign the upper bound on the Jacobian determinant value independently from the lower bound choice. We implemented it with a simple and fast quadratic-like penalty function [27].

This paper elaborates on the method in [27] and compares it empirically with methods based on other sufficient conditions as well as with the traditional Jacobian penalty method that uses a discrete grid [16], [25]. This paper is organized as follows. Section II reviews some related work. Section III proposes a new simple sufficient condition for the local invertibility of transformations based on B-splines and compares it with the box constraint [21], [22] empirically. Section IV proposes a new simple regularizer based on the local invertibility sufficient condition and presents 2D and 3D results.

II. Background

A. Mathematical model for nonrigid transformation

A 3D nonrigid transformation T: R3R3 can be written


where r = (x, y, z) and d(r) is the deformation. We model the 3D deformation d = (dx, dy, dz) using a tensor product of nth-order B-splines as follows:


where q [set membership] {x, y, z}, mq is knot spacing in q direction and β is a nth-order B-spline basis.

The goal in image registration is to estimate the deformation coefficients c_={ci,j,kq} by maximizing a similarity metric Ψ:


where g(r) and f(r) denote two 3D images.

To help stabilize the estimation, and to have physically plausible deformations, often we would like to ensure that the estimated coefficients ĉ correspond to a diffeomorphic transformation T. The methods in this paper are applicable to any similarity metric; for a survey of such metrics, see [28]. Section IV focuses on the l2 similarity metric for registering thorax CT images at different inhalations for the purpose of radiation therapy planning and monitoring.

B. Invertibility and diffeomorphic transformations

Invertibility of a nonrigid transformation T is a necessary condition for it to be diffeomorphic. T is diffeomorphic if both T and T−1 are continuously differentiable. If we use a B-spline basis with n ≥ 2 in (1), then T is continuously differentiable. (Musse et al. [18] addressed the case where n = 1. ) By the implicit function theorem, if the Jacobian matrix of T, denoted [nabla]T, is invertible everywhere, then near every point there exists a unique continuously differentiable local inverse. The determinant of the Jacobian for T, denoted |[nabla]T|, must be non-zero for diffeomorphic nonrigid image registration. Also for T to be orientation preserving, we want |[nabla]T| > 0.

Unfortunately, the condition |[nabla]T| > 0 everywhere does not by itself ensure that T is globally one-to-one. One way to ensure that T is invertible globally is to ensure that transformation maps the boundary of the domain onto itself [18], [20]. However, we do not enforce such boundary conditions in this paper because the field of view for thorax inhale and exhale CT images does not contain the whole body and there is usually missing anatomy in the superior-inferior directions.

C. Related work

The goal of diffeomorphic nonrigid image registration with the parametric representation of deformation (2) is to maximize the similarity metric (3) subject to the constraint


In general this is an impractical constraint except when using linear deformation models [18], [19], [20] because r [set membership] R3 so there are uncountably many conditions. One way to simplify (4) is to replace the “∀r” requirement with a set of voxel grid points [16], [17]:


However, because C0 [subset or is implied by] C1, this does not guarantee local invertibility between grid points. Nevertheless the smoothness of B-spline bases helps regularize C1 so using the constraint C1 often provides fairly good results [16]. However, computing |[nabla]T (r; c)| at all the grid points is computationally expensive.

Simplifying the condition |[nabla]T (r; c)| > 0 over R3 always involves smaller sets than C0. Choi et al. [21] found box constraints for cubic B-spline deformation coefficients that ensure invertibility:


where K ≈ 2.05 in 2D and K ≈ 2.48 in 3D. The set C2 provides a sufficient condition for local invertibility because C2 [subset or is implied by] C0. However, C2 is a very restrictive constraint set that allows only very small deformations. To achieve large deformations, Rueckert et al. [22] composed several transformations that each satisfied this condition.

Kim et al. [24], [25], [26] suggested a sufficient condition for ensuring invertibility of cubic B-spline deformations that allows a larger family of deformations. Instead of restricting the absolute values of the coefficients as in (6), this condition limits the differences of adjacent B-spline coefficients:


where kx + ky + kz < 1. Although C3 [subset or is implied by] C0, this sufficient condition only allows large deformations with fairly small Jacobian determinant values. In particular, one can show that 1 −(kx + ky + kz) ≤ |[nabla]T (r; c)| ≤ (1 + kx)(1 + ky)(1 + kz) + (1+kx)kykz +kx(1+ky)kz +kxky(1+kz) ∀c [set membership] C3 [24], [25], [26]. This means that C3 does not allow acute volume changes locally. This is because the upper bound on the Jacobian determinant is determined by the lower bound design. For example, if we choose kq = 1/3 so that the lower bound for the Jacobian determinant |J| is 0, then the upper bound for the Jacobian determinant value would be automatically determined to 76/27 ≈ 2.8148 which is fairly small [25]. The next section provides new broader sets of sufficient conditions.

III. Local invertibility condition

A. Lemmas

We first extend Kim’s sufficient conditions for local invertibility to overcome two limitations [27]. Firstly, a nth-order B-spline basis (n ≥ 1) can be used instead of cubic B-spline basis for deformation modeling. Secondly, the upper bound of Jacobian determinant can be designed independently from the lower bound of Jacobian determinant.

Lemma 1

For concise notation, denote the Jacobian J = [nabla]T of a 3D transformation as


Then the corresponding determinant is given by


Suppose that the elements of the 3D Jacobian determinant satisfy xi [set membership] Ii, i = 1, …, 9 where Ii [subset or is implied by] R are compact intervals. Then |J| achieves its global maximum and minimum values over I = I1×···×I9 and those maximum and minimum values are achieved for a point xi for which xi{maxIi,minIi} for ∀i = 1, …, 9.

The Appendices have the proofs of these Lemmas. This Lemma implies that we can determine the global minimum and maximum of |J| over the compact set I “simply” by calculating the 29 possible values of |J| at the vertices of I. (It is trivial to apply this Lemma to 2D cases.)

Kim et al. provided a specific formula for the ‘possible’ maximum and minimum of |J| for given ranges of each xi value using Karush-Kuhn-Tucker conditions [25]. We suggest next a generalization using Lemma 1.

Lemma 2

Suppose that xikqi<12 where qi = x for i = 2, 3, qi = y for i = 4, 6 and qi = z for i = 7, 8. Also suppose that −kpi ≤xi ≤Kpiwhere pi = x for i = 1, pi = y for i = 5 and pi = z for i = 9. Then min |J| = 1 − (kx + ky + kz) and max |J| = (1 + Kx)(1 + Ky)(1 + Kz) + (1 + Kx)kykz + kx(1 + Ky)kz + kxky(1 + Kz). In other words,


Kim’s proposition was restricted to the case where Kx = kx, Ky = ky, and Kz = kz. To ensure local invertibility, kx + ky + kz should be less than 1, where each kq is positive, so that the lower bound in (9) is positive.

Kim et al. showed a second proposition about the relationship between the first partial derivative of deformation and adjacent deformation coefficients for the cubic B-spline basis case [25]. We show next that this relation is also valid for general nth-order B-spline bases (n ≥ 1). We also generalize the bounds used by Kim et al. with Lemma 2 [27].

Lemma 3

If bmci+1,j,kqci,j,kqbM for ∀i, j, k, then bmmxxdq(r_)bMmx for ∀r where q [set membership]{x, y, z} Similarly, if bmci,j+1,kqci,j,kqbM for ∀i, j, k,, then bmmyydq(r_)bMmy and if bmci,j,k+1qci,j,kqbM for ∀i, j, k, then bmmzzdq(r_)bMmz for ∀r respectively.

This Lemma limits the range of values of the first derivative of d(r) over R3 by restricting the differences of adjacent deformation coefficients. Combined, Lemmas 2 and 3 show that one can obtain a transformation T that is everywhere locally invertible by maximizing a similarity metric subject to constraints on the differences between adjacent deformation coefficients, as summarized in the following Theorem.

Theorem 1

Suppose 0kq<12 for q [set membership] {x, y, z}. Define:


In (2), if c [set membership] C4 then |J| satisfies the bounds in (9) ∀r [set membership] R3. Moreover, if kx + ky + kz < 1, then the transformation (2) is locally invertible everywhere.

This theorem applies to deformations based on any nth-order B-spline basis. We set the lower and upper bounds for |J| by setting appropriate kq and Kq values for q [set membership] {x, y, z}.

B. Restrictions

Theorem 1 establishes that c [set membership] C4 is a simple sufficient condition for local invertibility. However, C4 does not allow all possible locally invertible deformations, i.e., C4 [subset or is implied by] C0. Then one can ask how restrictive this sufficient condition is.

Although C4 allows for acute volume expansion, it precludes acute volume shrinkage. Fig. 1 illustrates this limitation for a 1D transformation. The desired transformation maps [0.0 0.6] to [0.3 0.6], i.e., T (x) = x + d(x) where d(x) = 0.3 − x/2 (acute volume shrinkage). This deformation belongs to C0 because 1<d(x)x<. However if we impose the sufficient condition 0.33<d(x)x, then Fig. 1 shows that acute volume shrinkage is precluded because the minimum derivative of the transformation is 0.67. The constrained transformation maps [0 0.6] to [0.3 0.7] instead of [0.3 0.6]. More generally, when we choose kx, ky and kz subject to kx + ky + kz < 1 to ensure invertibility, C4 imposes restrictions for acute volume changes in each direction.

Fig. 1
Illustration of limitation of C4. The constrained transformation maps [0 0.6] to [0.3 0.7] instead of [0.3 0.6].

The 2D case illustrates the solution space of C4 in terms of Lemma 2. Lemma 2 is trivial for a 2D Jacobian determinant |J| = (1+a)(1+d)−bc where J=(abcd). A deformation c d having a positive Jacobian determinant must satisfy (1+a)(1+ d) > bc. We can introduce a free parameter k such that |J| is always positive if (1 + a)(1 + d) > k and bc < k for any k. Fig. 2 visualizes the solution space for 2D invertible deformations in terms of a, b, c, d, and k. For fixed k, any values of (a, d) that lie above the upper line or below the lower line yield a positive Jacobian determinant if (b, c) lies between these lines. Lines vary as k varies. To allow acute volume shrinkage, we need k to be close to 0 as observed in Fig. 1. However smaller k values imply more restrictive sets for (b, c).

Fig. 2
Solution space for 2D positive Jacobian determinant. Smaller k values admit smaller a, d values but preclude more values of b, c.

Lemma 2 corresponds to fixing k = kxky such that kx + ky < 1 and kx ≥ 0, ky ≥ 0. This yields the rectangular areas for a, b, c, d shown in Fig. 3 (for kx = ky = 1/2 and k = 1/4). Thus Theorem 1 not only uses a fixed value for k, but also imposes restrictive box constraints on the deformation derivatives. However it still has a larger solution space than traditional box constraints on the B-spline coefficients such as [21]. Because k is fairly small, relaxing this sufficient condition may allow larger volume shrinkage [27].

Fig. 3
Local invertibility sufficient condition space in 2D, for k = 1/4, kx = ky = 1/2 and Kx = Ky = ∞. C4 corresponds to using a fixed k value. (a) a > −1/2 and d > −1/2. (b) |b| < 1/2 and |c| < 1/2.

C. Concatenating transformations

Since C4 is a restrictive sufficient condition, it may not contain all real deformations of interest. To allow larger deformations, we can concatenate multiple elemental transformations that belong to C4, i.e., let T (r) = TN (··· (T2(T1(r)))) where each Tk satisfies C4. Since each Tk is diffeomorphic, T is also diffeomorphic.

Rueckert et al. [22] used a box constraint C2 [21] to guarantee that each elemental transformation is diffeomorphic. We use C4 for our elemental transformations. This should require fewer elemental transformations because C4 allows a larger solution space than C2, as illustrated in the next section.

D. 2D simulation: warping a disk to a “C” shape

We applied several constrained nonrigid image registration methods to the challenging registration problem shown in Fig. 4. We placed deformation knot points every 4th pixel, i.e., mx = my = 4. The data fit term used sum of squared differences. For optimization we used augmented Lagrangian multipliers [17] with the conjugate gradient method. Line search step size was determined by one step of Newton’s method. We used fast B-spline interpolation [29], [30], [31] with a 4-level multiresolution scheme [32].

Fig. 4
Images for illustrating 2D nonrigid registration.

Fig. 5 shows the unconstrained registration and results of using C1, C2 and C4. The unconstrained result in Fig. 5(a) shows some unrealistic warping such as folding. Fig. 5(b) shows the regularized deformed images with a Jacobian penalty based on C1. This shows a more regular warp than Fig. 5(a). However, C1 allows a larger solution space than the ideal solution space C0.

Fig. 5
Deformed images (left) and their warped grids (right)

Fig. 5(c) and (d) show the limitation of using a single warp based on C2 and C4 respectively. The sufficient conditions C2 and C4 do not contain the complicated diffeomorphic transformation needed to map the source image to the target image in Fig. 4. However, this warp can be achieved satisfactorily by composing just 3 warps that each belong to C4, as shown in Fig. 6(b). In contrast, to achieve a satisfactory warp by composing transformations that lie in the box constraint C2 [21], [22] required about 30 concatenations, as shown in Fig. 6(a). For larger and more complicated deformations, our proposed constraint C4 can be used as a simple elemental transformation to provide diffeomorphic composite transformations.

Fig. 6
Proposed constraint requires much less transformations to achieve a satisfiable deformation than the box constraint.

IV. Simple regularizer based on local invertibility condition

A. Proposed simple regularizer

If we want to strictly ensure local invertibility, then we maximize a similarity metric subject to the linear constraints c [set membership] C4. However, to simplify the computation, we can relax the invertibility condition by using a penalty method [16], [27]. In a penalty method we maximize an objective function that is the similarity metric minus a penalty function that encourages the invertibility condition, but does not enforce it strictly.

We propose to construct a penalty function based on the following piecewise quadratic function:


which is illustrated in Fig 7. The argument t denotes a difference between two adjacent deformation coefficients. This function does not strictly constrain such differences, but its first and second derivatives are simple and convenient for use in optimization algorithms such as conjugate gradient. The final new penalty function is

Fig. 7
A variant of quadratic penalty function (solid) and real constraints (dashed) used with constraint set C4.


where ζ1q,r=mqkq for ∀r [set membership] {x, y, z}, ζ2q,r=mqkq for rq and ζ2q,r=mqKq for r = q. Note that choosing ζ1 = ζ2 = 0 would correspond to a quadratic roughness penalty over B-spline coefficients, which is akin to encouraging the volume preserving condition |J| = 1, ∀r.

Being based on the somewhat restrictive solution space C4, the new penalty method can encourage the local invertibility on the whole continuous domain with a fast and memory efficient implementation. This implementation is possible because C4 does not require additional B-spline interpolations beyond the interpolations needed for the data fitting term. It also encourages the smoothness of deformations inherently because it constrains the differences between adjacent deformation coefficients. In contrast, using C0 or C1 is much more expensive for one transformation.

B. Incorporating a priori knowledge of motions

For diffeomorphic transformations using Theorem 1, the usual choice would be kx = ky = kz = 1/3 − ε for some small ε. However, if we have a priori knowledge about the deformation, then we can assign each kq accordingly. For instance, for registering thorax inhale and exhale images, we can assign kx = ky = 1/4 − ε and kz = 1/2 − ε because the deformation in the z direction is larger due to diaphragm motion, whereas the deformations in the x and y directions are smaller. With this design, C4 allows 50% local shrinkage along z and 75% local shrinkage along x and y instead of allowing 67% shrinkage in each direction. We can use this sufficient condition for the proposed simple regularizer to encourage local invertibility of the deformation.

C. 2D simulation: expansion and shrinkage

We applied nonrigid image registration to the 256 × 256 images in Fig. 8 using no constraint, a Jacobian penalty based on C1, a quadratic roughness penalty [7], a regularizer based on Kim’s constraint C3, and our proposed penalty method based on C4. Fig. 8 has an expanding circle and a shrinking ellipse to illustrate the difference between C3 and C4. Since we have a priori knowledge about vertical motion, we investigated two sets of parameters in C4: a symmetric way with kx = ky = 1/2 − 0.01 × 1/2 as well as an asymmetric way with kx = 0.35 − 0.01 × 0.35 < ky = 0.65 − 0.65 × 0.01.

Fig. 8
Images for illustrating expansion and shrinkage.

We placed deformation knot points every 4th pixel. The data fit term used sum of squared differences. For optimization we used the conjugate gradient method. Line search step size was determined by one step of Newton’s method. We used fast B-spline interpolation and the 4-level multiresolution scheme as in III-D. We ran 200 iterations for each level or ran until the l0 norm of the gradient is less than the machine accuracy. We checked the local invertibility by computing Jacobian determinant values on a grid 10 times finer than the image resolution.

Fig. 9 quantifies the tradeoff between image similarity and local invertibility for the 5 different registration methods for a range of regularization parameters. The horizontal axis is the root mean square (RMS) difference between the deformed image and the target image (log scale) and the vertical axis is the number of the finer (10 times) voxel grid points having a non-positive Jacobian determinant (log scale). We took a log after adding 1 for the number of non-positive Jacobian determinant since the lowest number of it is 0.

Fig. 9
RMS error and negative Jacobian determinant trade-off for different regularization parameters. (log scale)

For the unconstrained case, the RMS difference was 0.0109 and the number of negative Jacobian determinants was 497644. As the regularization parameters decrease, the RMS differences and the number of negative Jacobian determinants of all other methods approached closely to these values. This is expected because the unconstrained case is the same as any other penalty method with regularization parameter 0.

As we increased the regularization parameters, the number of negative Jacobian determinants “generally” decreases, eventually towards zero, although not always monotonically. The RMS differences also “generally” increase as the regularization parameters increase for most methods except Jacobian penalty. This is because the Jacobian constraint C1 contains the original constraint C0 and it does not restrict the deformation so that it can achieve low RMS difference for strong penalty parameters. For properly chosen regularization parameters, symmetric/asymmetric proposed simple penalties show fairly good performance compared to Kim’s or quadratic penalty methods based on more restrictive sufficient conditions (Kim’s: Kx = kx, Ky = ky and quadratic: Kx = kx = 0, Ky = ky = 0).

Table I shows the best RMS difference of image for each method with zero non-positive Jabocian determinant values over the 10 times finer grid. The proposed simple penalty with a priori motion information performed well compared to Jacobian penalty. However, as the regularization method depends on more restrictive condition, the RMS difference is larger. It clearly shows that a quadratic penalty oversmoothes the deformations.

The best RMS error for each method with zero negative Jabocian determinant in 2D simulation.

Our proposed asymmetric penalty performed a little better than Jacobian penalty in these experiments. However, our proposed penalty may not always perform better. It depends on the convergence, regularization parameter, image structure and so on. Since the data fitting term is non-convex, local minima may affect the result, too. However, for simpler cases like in Fig. 8 our proposed regularization method may perform close to the Jacobian penalty method. In the next section, we apply both methods to the 3D real CT images of a patient.

D. 3D real CT images

We applied our proposed regularization method (10) to the problem of registering 3D breath-hold X-ray CT images of a real oncology patient scanned at inhale and at exhale. These images are useful for radiation treatment planning. The image size was 396 × 256 × 128 as shown in Fig 10. We chose kx = ky = 1/4 − 0.01 × 1/4 and kz = 1/2 − 0.01 × 1/2 because we expect the deformation in the z direction to be larger than the deformations in the x and y directions due to diaphragm motion.

Fig. 10
3D source (exhale) and target (inhale) X-ray CT images.

We used the same methods as in Section IV-C except for the multiresolution scheme. For the first 3 levels of multiresolution, the knot spacing was every 8 pixels for downsampled images, and for the last level of multiresolution the knot spacing was every 4 pixels. We ran 120 iterations at each level to see the convergence properties. The regularization parameter that multiplies (10) was chosen experimentally to achieve the minimum value of data fitting term such that all Jacobian values on the image grid were positive.

Fig. 11 shows the difference images between the target image and the deformed images. As expected, the difference image for unconstrained registration in Fig. 11(a) has smaller values than the constrained difference images in Fig. 11(b) and (c). The RMS difference for unconstrained registration was the smallest, which was 19.9 HU. The RMS errors of the Jacobian penalty (25.9 HU) and the proposed penalty (29.2 HU) were somewhat higher. However, Fig. 12(a) shows that unconstrained registration yields an unrealistic warped grid. The number of negative Jacobian determinant voxels was 316914 out of 12582912 voxels. Fig. 12(c) shows a smoother warp than Fig. 12(b) because our proposed penalty method is based on C4 which is a smaller set than C1. Our proposed method has smoothness property implicitly because it restricts the range of the differences between adjacent B-spline coefficients.

Fig. 11
Differences between 3D target and deformed images.
Fig. 12
Warped grids for 3D inhale-exhale registration.

The proposed penalty method was much faster and more memory efficient than the traditional Jacobian penalty method per iteration. If one uses the sum of squared error as the data fitting term and penalizes negative Jacobian determinant values on each image grid point in 3D with cubic B-splines, then the interpolations needed to compute the gradients of the direct Jacobian penalty function require about 1.8 times more operations than the interpolations needed for the gradient of the data fitting term. Table II shows the computational cost for one iteration at the last (finest) level of the multiresolution procedure. Our proposed method requires only slightly more time per iteration than unconstrained registration, and much less time than using a Jacobian penalty. Furthermore, in this simulation, the traditional Jacobian penalty method required about twice as much memory as our proposed method because it must store the interpolation results for the Jacobian gradient. Fig. 13 shows the convergence of each method.

Fig. 13
Convergence of each method.
Computational cost at the finest level

We could compose a coarse resolution warp based on (10) with one full sequence of coarse-to-fine warp based on (10) to reduce RMS differences further with only slight increase in computation.

V. Discussion

We proposed a new condition Theorem 1 that is sufficient to ensure the local invertibility of transformations based on B-splines. Its limitation can be overcome by using composite transformations. This proposed sufficient condition can be used with constrained optimization such as augmented Lagrangian multiplier method [17] or Dykstra’s cyclic projection method [25].

We showed that the proposed sufficient condition is more general than other simple sufficient conditions that ensure local invertibility everywhere such as box constraint [21]. When used in composite transformations, it requires many fewer transformations to achieve comparable deformations [22].

We also relaxed our local invertibility condition by a simple quadratic-like penalty. This approach achieves more flexible image matching compared to other penalty methods based on more restrictive local invertibility conditions. For practical use in a thorax image registration, we used a single transformation with a simple quadratic-like penalty that encourages c [set membership] C4. This gave a fairly good deformation with no negative Jacobian determinant values on image voxel grid points. This approach is much simpler and faster than the traditional Jacobian determinant penalty and is more memory efficient.

Some application areas require not only local invertibility, but also require computing the inverse transformation. One approach is to estimate both forward and backward image registration parameters with consistency regularizer [33]. Using both consistency regularizer and our proposed regularizer can be interesting future work.

In the 3D thorax registration example, we observed some bone warping in each of the deformed images; a natural direction for further research would be to use a rigidity penalty term [34], [35].


The authors thank Dr. Marc Kessler for the 3D CT data and Dr. James Balter for discussions about image registration. We also thank the reviewers for insightful questions that improved the paper.

This work was supported in part by NIH/NCI grant 1P01 CA87634.


An external file that holds a picture, illustration, etc.
Object name is nihms254895b1.gif

Se Young Chun Se Young Chun received his BSE from the School of Electrical Engineering, Seoul National University, Seoul, South Korea in 1999. He started his graduate study at the University of Michigan, Ann Arbor in 2003 and received his MSE in Electrical Engineering: Systems and MS in Mathematics in April and December 2005, respectively. He is currently pursuing his Ph.D. in Electrical Engineering: Systems under the supervision of Prof. Jeffrey A. Fessler. He was a recipient of a 2-year scholarship from the Ministry of Information and Communications, South Korea during his graduate study. His research interests are diffeomorphic nonrigid image registration, motion-compensated image reconstruction in medical imaging as well as statistical signal processing in brain signals.

An external file that holds a picture, illustration, etc.
Object name is nihms254895b2.gif

Jeffrey A. Fessler Jeff Fessler received the BSEE degree from Purdue University in 1985, the MSEE degree from Stanford University in 1986, and the M.S. degree in Statistics from Stanford University in 1989. From 1985 to 1988 he was a National Science Foundation Graduate Fellow at Stanford, where he earned a Ph.D. in electrical engineering in 1990. He has worked at the University of Michigan since then. From 1991 to 1992 he was a Department of Energy Alexander Hollaender Post-Doctoral Fellow in the Division of Nuclear Medicine. From 1993 to 1995 he was an Assistant Professor in Nuclear Medicine and the Bioengineering Program. He is now a Professor in the Departments of Electrical Engineering and Computer Science, Radiology, and Biomedical Engineering. He became a Fellow of the IEEE in 2006, for contributions to the theory and practice of image reconstruction. He received the Francois Erbsmann award for his IPMI93 presentation. He serves as an associate editor for IEEE Transactions on Medical Imaging and is a past associate editor for the IEEE Transactions on Image Processing and the IEEE Signal Processing Letters. He was co-chair of the 1997 SPIE conference on Image Reconstruction and Restoration, technical program co-chair of the 2002 IEEE International Symposium on Biomedical Imaging (ISBI), and general chair of ISBI 2007. His research interests are in statistical aspects of imaging problems, and he has supervised doctoral research in PET, SPECT, X-ray CT, MRI, and optical imaging problems.

Appendix A

Proof of Lemma 1


The global maximum and minimum values exist since |J| is continuous on the compact set I1×···×I9. Suppose that (x1, ···, x9) achieves the global minimum value of |J| and min Ik < xk < max Ik for some k. Fix all xi except xk, |J| is an affine function with respect to xk so |J| can achieve equal or better global minimum value on either xk = min Ik or xk = max Ik. The same argument can be applied for all xi such that min Ii < xi < max Ii and thus it generates a contradiction. The same argument can be applied to the global maximum case.

Appendix B

Proof of Lemma 2


By Lemma 1, we need to evaluate |J| only on x1 [set membership] {−kx, Kx}, x4 [set membership] {−ky, Ky}, x9 [set membership] {−kz, Kz} and xi [set membership] {−kqi, kqi} where qi = x for i = 2, 3, qi = y for i = 4, 6 and qi = z for i = 7, 8. For fixed xi except x1, |J(x1)| = (1 + x1){(1 + x5)(1 + x9) − x6x8} + c where c is a constant for x1 and (1 + x5)(1 + x9) − x6x8 is always positive under given conditions. So x1 = Kx for max |J| and x1 = −kx for min |J|. Similarly we determine x5 and x9. For fixed xi except x2, |J(x2)| = x2{x6x7 −(1+x9)x4}+c where c is a constant for x2. For min |J|, x2 = −kx if x4 = −ky and x2 = kx if x4 = ky. In other words, x2x4 = kxky for min |J|. Similarly, x2x4 = −kxky for max |J|. In this fashion, x6x8 and x3x7 will be determined for max |J| and min |J|. From these results, one can induce that x2x6x7 + x3x4x8 = 0 for max |J| and x2x6x7 = x3x4x8 = −kxkykz for min |J|.

Appendix C

Proof of Lemma 3


For d(x) = Σi ciβn(x/mxi), by using xβn(x)=βn1(x+1/2)βn1(x1/2) in [30]


Using the constraints bmci+1,j,kqci,j,kqbM and the property Σi βn(x/mxi) = 1, we have the bounds


Similarly, xdq(r_)bm/mx. The other directions y, z can be proved similarly.


1. Reyes M, Malandain G, Koulibaly PM, Gonzalez-Ballester MA, Darcourt J. Model-based respiratory motion compensation for emission tomography image reconstruction. Physics in Medicine and Biology. 2007 June;52(12):3579–3600. [PubMed]
2. Qiao F, Pan T, Jr, JWC, Mawlawi OR. A motion-incorporated reconstruction method for gated PET studies. Physics in Medicine and Biology. 2006 August;51(15):3769–3783. [PubMed]
3. Lamare F, Cresson T, Savean J, Rest CCL, Reader AJ, Visvikis D. Respiratory motion correction for PET oncology applications using affine transformation of list mode data. Physics in Medicine and Biology. 2007 January;52(1):121–140. [PubMed]
4. Lamare F, Carbayo MJL, Cresson T, Kontaxakis G, Santos A, Rest CCL, Reader AJ, Visvikis D. List-mode-based reconstruction for respiratory motion correction in PET using non-rigid body transformations. Physics in Medicine and Biology. 2007 September;52(17):5187–5204. [PubMed]
5. Mair BA, Gilland DR, Sun J. Estimation of images and nonrigid deformations in gated emission CT. IEEE Transactions on Medical Imaging. 2006;25(9):1130–1144. [PubMed]
6. Jacobson MW, Fessler JA. Joint estimation of respiratory motion and activity in 4D PET using CT side information. Proc IEEE Intl Symp Biomed Imag. 2006:275–8.
7. Zeng R, Fessler JA, Balter JM. Estimating 3-D respiratory motion from orbiting views by tomographic image registration. IEEE Transactions on Medical Imaging. 2007;26(2):153–163. [PMC free article] [PubMed]
8. Crum WR, Hartkrns T, Hill DLG. Non-rigid image registration: theory and practice. The British Journal of Radiology. 2004;77:S140–S153. [PubMed]
9. Holden M. A review of geometric transformations for nonrigid body registration. IEEE Transactions on Medical Imaging. 2008;27:111–128. [PubMed]
10. Szeliski R, Coughlan J. Spline-based image registration. Int J Comput Vision. 1997;22(3):199–218.
11. Unser M. Splines: A perfect fit for signal and image processing. IEEE Signal Processing Magazine. 1999 November;16(6):22–38.
12. Kybic J, Unser M. Fast parametric elastic image registration. IEEE Transactions on Image Processing. 2003 November;12(11):1427–1442. [PubMed]
13. Rueckert D, Sonoda LI, Hayes C, Hill DL, Leach MO, Hawkes DJ. Nonrigid registration using free-form deformations: application to breast MR images. IEEE Transactions on Medical Imaging. 1999;18(8):712–21. [PubMed]
14. Sánchez Sorzano C, Thévenaz P, Unser M. Elastic registration of biological images using vector-spline regularization. IEEE Transactions on Biomedical Engineering. 2005 April;52(4):652–663. [PubMed]
15. Rohlfing T, CRM, Bluemke DA, Jacobs MA. Volume-preserving non-rigid registration of MR breast images using free-form deformation with an incompressibility constraint. IEEE Transactions on Medical Imaging. 2003;22(6):730–741. [PubMed]
16. Kybic J, Thévenaz P, Nirkko A, Unser M. Unwarping of unidirectionally distorted EPI images. IEEE Transactions on Medical Imaging. 2000 February;19(2):80–93. [PubMed]
17. Sdika M. A fast nonrigid image registration with constraints on the Jacobian using large scale constrained optimization. IEEE Transactions on Medical Imaging. 2008 February;27(2):271–281. [PubMed]
18. Musse O, Heitz F, Armspach J. Topology preserving deformable image matching using constrained hierarchical parametric models. IEEE Transactions on Image Processing. 2001;10(7):1081–1093. [PubMed]
19. Karacali B, Davatzikos C. Estimating topology preserving and smooth displacement fields. IEEE Transactions on Medical Imaging. 2004;23(7):868–880. [PubMed]
20. Noblet V, Heinrich C, Heitz F, Armspach J. 3-D deformable image registration: A topology preservation scheme based on hierarchical deformation models and interval analysis optimization. IEEE Transactions on Image Processing. 2005;14(5):553–566. [PubMed]
21. Choi Y, Lee S. Injectivity conditions of 2D and 3D uniform cubic B-spline functions. Graph Models. 2000;62(6):411–427.
22. Rueckert D, Aljabar P, Heckemann RA, Hajnal JV, Hammers A. Diffeomorphic registration using b-splines. MICCAI; 2006. pp. 702–709. [PubMed]
23. Rohde GK, Aldroubi A, Dawant BM. The adaptive bases algorithm for intensity-based nonrigid image registration. IEEE Transactions on Medical Imaging. 2003;22(11):1470–1479. [PubMed]
24. Kim J, Fessler JA. Nonrigid image registration using constrained optimization. SIAM Conference on Imaging Science; 2004.
25. Kim J. PhD dissertation. University of Michigan; 2004. Intensity based image registration using robust similarity measure and constrained optimization: applications for radiation therapy.
26. Kim J. Non-rigid image registration using constrained optimization (korean) Journal of Korean Information and Communications. 2004;29(10C):1402–1413.
27. Chun SY, Fessler JA. Regularized methods for topology-preserving smooth nonrigid image registration using B-spline basis. Proc IEEE Intl Symp Biomed Imag. 2008:1099–1102.
28. Hill DLG, Batchelor PG, Holden M, Hawkes DJ. Medical image registration. Physics in Medicine and Biology. 2001;46(3):R1–R45. [PubMed]
29. Unser M, Aldroubi A, Eden M. Fast B-Spline transforms for continuous image representation and interpolation. IEEE Transactions on Pattern Analysis and Machine Intelligence. 1991 March;13(3):277–285.
30. Unser M, Aldroubi A, Eden M. B-Spline signal processing: Part I—Theory. IEEE Transactions on Signal Processing. 1993 February;41(2):821–833.
31. Unser M, Aldroubi A, Eden M. B-Spline signal processing: Part II—Efficient design and applications. IEEE Transactions on Signal Processing. 1993 February;41(2):834–848.
32. Unser M, Aldroubi A, Eden M. The L2-polynomial spline pyramid. IEEE Transactions on Pattern Analysis and Machine Intelligence. 1993 April;15(4):364–379.
33. Christensen GE, Johnson HJ. Consistent image registration. IEEE Trans Med Imag. 2001 Jul;20(7):568–82. [PubMed]
34. Ruan D, Fessler JA, Roberson M, Balter J, Kessler M. Nonrigid registration using regularization that accomodates local tissue rigidity. Proc of SPIE. 2006;6144:346–354.
35. Loeckx D, Maes F, Vandermeulen D, Suetens P. Nonrigid image registration using free-form deformations with a local rigidity constraint. Medical Image Computing and Computer-Assisted Intervention. 2004;LNCS3216:639–646.