Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC3004236

Formats

Article sections

- Abstract
- I. Introduction
- II. Background
- III. Local invertibility condition
- IV. Simple regularizer based on local invertibility condition
- V. Discussion
- References

Authors

Related links

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.2011116PMCID: PMC3004236

NIHMSID: NIHMS254895

Department of Electrical Engineering and Computer Science, University of Michigan, Ann Arbor, MI, 48109 USA

See other articles in PMC that cite the published article.

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.

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 *n*th-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.

A 3D nonrigid transformation *T*: **R**^{3} → **R**^{3} can be written

$$\underset{\_}{T}(\underset{\_}{r})=\underset{\_}{r}+\underset{\_}{d}(\underset{\_}{r}),$$

(1)

where *r* = (*x, y, z*) and *d*(*r*) is the deformation. We model the 3D deformation *d* = (*d ^{x}, d^{y}, d^{z}*) using a tensor product of

$${d}^{q}(\underset{\_}{r})=\sum _{\mathit{ijk}}{c}_{\mathit{ijk}}^{q}\beta (\frac{x}{{m}_{x}}-i)\beta (\frac{y}{{m}_{y}}-j)\beta (\frac{z}{{m}_{z}}-k),$$

(2)

where *q* {*x, y, z*}, *m _{q}* is knot spacing in

The goal in image registration is to estimate the deformation coefficients $\underset{\_}{c}=\{{c}_{i,j,k}^{q}\}$ by maximizing a similarity metric Ψ:

$$\underset{\_}{\widehat{c}}=arg\underset{\underset{\_}{c}}{max}\mathrm{\Psi}[g(\xb7),f(\underset{\_}{T}(\xb7;\underset{\_}{c}))]$$

(3)

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 *l*_{2} similarity metric for registering thorax CT images at different inhalations for the purpose of radiation therapy planning and monitoring.

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 *T*, is invertible everywhere, then near every point there exists a unique continuously differentiable local inverse. The determinant of the Jacobian for *T*, denoted |*T*|, must be non-zero for diffeomorphic nonrigid image registration. Also for *T* to be orientation preserving, we want |*T*| > 0.

Unfortunately, the condition |*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.

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

$$\underset{\_}{c}\in {C}_{0}\triangleq \{\underset{\_}{c}:\mid \nabla \underset{\_}{T}(\underset{\_}{r};\underset{\_}{c})\mid >0,\phantom{\rule{0.16667em}{0ex}}\forall \underset{\_}{r}\in {\mathbf{R}}^{3}\}.$$

(4)

In general this is an impractical constraint except when using linear deformation models [18], [19], [20] because *r* **R**^{3} 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]:

$${C}_{1}\triangleq \{\underset{\_}{c}:\mid \nabla \underset{\_}{T}(\underset{\_}{r};\underset{\_}{c})\mid >0,\phantom{\rule{0.16667em}{0ex}}\underset{\_}{r}\in \text{grid}\phantom{\rule{0.16667em}{0ex}}\text{points}\}.$$

(5)

However, because *C*_{0} *C*_{1}, this does not guarantee local invertibility *between* grid points. Nevertheless the smoothness of B-spline bases helps regularize *C*_{1} so using the constraint *C*_{1} often provides fairly good results [16]. However, computing |*T* (*r; c*)| at all the grid points is computationally expensive.

Simplifying the condition |*T* (*r; c*)| > 0 over **R**^{3} always involves smaller sets than *C*_{0}. Choi *et al.* [21] found box constraints for cubic B-spline deformation coefficients that ensure invertibility:

$${C}_{2}\triangleq \{\underset{\_}{c}:\mid {c}_{i,j,k}^{q}\mid <{m}_{q}/K,\phantom{\rule{0.16667em}{0ex}}\forall i,j,k\},$$

(6)

where *K* ≈ 2.05 in 2D and *K* ≈ 2.48 in 3D. The set *C*_{2} provides a sufficient condition for local invertibility because *C*_{2} *C*_{0}. However, *C*_{2} 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:

$${C}_{3}\triangleq \underset{q\in \{x,y,z\}}{\cap}\{\underset{\_}{c}:\mid {c}_{i+1,j,k}^{q}-{c}_{i,j,k}^{q}\mid <{m}_{q}{k}_{q},\mid {c}_{i,j+1,k}^{q}-{c}_{i,j,k}^{q}\mid <{m}_{q}{k}_{q},\mid {c}_{i,j,k+1}^{q}-{c}_{i,j,k}^{q}\mid <{m}_{q}{k}_{q},\phantom{\rule{0.38889em}{0ex}}\forall i,j,k\},$$

(7)

where *k _{x}* +

We first extend Kim’s sufficient conditions for local invertibility to overcome two limitations [27]. Firstly, a *n*th-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.

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

$$\mathbf{J}=\mathbf{I}+\left(\begin{array}{ccc}{x}_{1}& {x}_{2}& {x}_{3}\\ {x}_{4}& {x}_{5}& {x}_{6}\\ {x}_{7}& {x}_{8}& {x}_{9}\end{array}\right).$$

Then the corresponding determinant is given by

$$\mid \mathbf{J}\mid =(1+{x}_{1})(1+{x}_{5})(1+{x}_{9})+{x}_{2}{x}_{6}{x}_{7}+{x}_{3}{x}_{4}{x}_{8}-(1+{x}_{1}){x}_{6}{x}_{8}-(1+{x}_{5}){x}_{3}{x}_{7}-(1+{x}_{9}){x}_{2}{x}_{4}.$$

(8)

Suppose that the elements of the 3D Jacobian determinant satisfy x_{i} I_{i}, i = 1, …, 9 where I_{i} **R** are compact intervals. Then |**J**| achieves its global maximum and minimum values over I = I_{1}×···×I_{9} and those maximum and minimum values are achieved for a point
${x}_{i}^{\ast}$ for which
${x}_{i}^{\ast}\in \{max{I}_{i},min{I}_{i}\}$ 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 2^{9} 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 *x _{i}* value using Karush-Kuhn-Tucker conditions [25]. We suggest next a generalization using Lemma 1.

Suppose that
$\mid {x}_{i}\mid \le {k}_{{q}_{i}}<{\scriptstyle \frac{1}{2}}$ where q_{i} = x for i = 2, 3, q_{i} = y for i = 4, 6 and q_{i} = z for i = 7, 8. Also suppose that −k_{pi} ≤x_{i} ≤K_{pi}where p_{i} = x for i = 1, p_{i} = y for i = 5 and p_{i} = z for i = 9. Then min |**J**| = 1 − (k_{x} + k_{y} + k_{z}) and max |**J**| = (1 + K_{x})(1 + K_{y})(1 + K_{z}) + (1 + K_{x})k_{y}k_{z} + k_{x}(1 + K_{y})k_{z} + k_{x}k_{y}(1 + K_{z}). In other words,

$$1-({k}_{x}+{k}_{y}+{k}_{z})\le \mid \mathbf{J}\mid \le (1+{K}_{x})(1+{K}_{y})\xb7(1+{K}_{z})+(1+{K}_{x}){k}_{y}{k}_{z}+{k}_{x}(1+{K}_{y}){k}_{z}+{k}_{x}{k}_{y}(1+{K}_{z}).$$

(9)

Kim’s proposition was restricted to the case where *K _{x}* =

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 *n*th-order B-spline bases (*n* ≥ 1). We also generalize the bounds used by Kim *et al.* with Lemma 2 [27].

If ${b}_{m}\le {c}_{i+1,j,k}^{q}-{c}_{i,j,k}^{q}\le {b}_{M}$ for ∀i, j, k, then ${\scriptstyle \frac{{b}_{m}}{{m}_{x}}}\le {\scriptstyle \frac{\partial}{\partial x}}{d}^{q}(\underset{\_}{r})\le {\scriptstyle \frac{{b}_{M}}{{m}_{x}}}$ for ∀r where q {x, y, z} Similarly, if ${b}_{m}\le {c}_{i,j+1,k}^{q}-{c}_{i,j,k}^{q}\le {b}_{M}$ for ∀i, j, k,, then ${\scriptstyle \frac{{b}_{m}}{{m}_{y}}}\le {\scriptstyle \frac{\partial}{\partial y}}{d}^{q}(\underset{\_}{r})\le {\scriptstyle \frac{{b}_{M}}{{m}_{y}}}$ and if ${b}_{m}\le {c}_{i,j,k+1}^{q}-{c}_{i,j,k}^{q}\le {b}_{M}$ for ∀i, j, k, then ${\scriptstyle \frac{{b}_{m}}{{m}_{z}}}\le {\scriptstyle \frac{\partial}{\partial z}}{d}^{q}(\underset{\_}{r})\le {\scriptstyle \frac{{b}_{M}}{{m}_{z}}}$ for ∀r respectively.

This Lemma limits the range of values of the first derivative of *d*(*r*) over **R**^{3} 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.

Suppose $0\le {k}_{q}<{\scriptstyle \frac{1}{2}}$ for q {x, y, z}. Define:

$${C}_{4}\triangleq \{\underset{\_}{c}:-{m}_{x}{k}_{x}\le {c}_{i+1,j,k}^{x}-{c}_{i,j,k}^{x}\le {m}_{x}{K}_{x},-{m}_{y}{k}_{y}\le {c}_{i,j+1,k}^{y}-{c}_{i,j,k}^{y}\le {m}_{y}{K}_{y},-{m}_{z}{k}_{z}\le {c}_{i,j,k+1}^{z}-{c}_{i,j,k}^{z}\le {m}_{z}{K}_{z},\mid {c}_{i+1,j,k}^{q}-{c}_{i,j,k}^{q}\mid \le {m}_{q}{k}_{q}\phantom{\rule{0.16667em}{0ex}}\mathit{for}\phantom{\rule{0.16667em}{0ex}}q=y,z,\mid {c}_{i,j+1,k}^{q}-{c}_{i,j,k}^{q}\mid \le {m}_{q}{k}_{q}\phantom{\rule{0.16667em}{0ex}}\mathit{for}\phantom{\rule{0.16667em}{0ex}}q=x,z,\mid {c}_{i,j,k+1}^{q}-{c}_{i,j,k}^{q}\mid \le {m}_{q}{k}_{q}\phantom{\rule{0.16667em}{0ex}}\mathit{for}\phantom{\rule{0.16667em}{0ex}}q=x,y,\forall i,j,k\}.$$

In (2), if c C_{4} then |**J**| satisfies the bounds in (9) ∀r **R**^{3}. Moreover, if k_{x} + k_{y} + k_{z} < 1, then the transformation (2) is locally invertible everywhere.

This theorem applies to deformations based on any *n*th-order B-spline basis. We set the lower and upper bounds for |**J**| by setting appropriate *k _{q}* and

Theorem 1 establishes that *c* *C*_{4} is a simple sufficient condition for local invertibility. However, *C*_{4} does not allow all possible locally invertible deformations, *i.e.*, *C*_{4} *C*_{0}. Then one can ask how restrictive this sufficient condition is.

Although *C*_{4} 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 *C*_{0} because
$-1<{\scriptstyle \frac{\partial d(x)}{\partial x}}<\infty $. However if we impose the sufficient condition
$-0.33<{\scriptstyle \frac{\partial d(x)}{\partial 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 *k _{x}*,

Illustration of limitation of *C*_{4}. 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 *C*_{4} in terms of Lemma 2. Lemma 2 is trivial for a 2D Jacobian determinant |**J**| = (1+*a*)(1+*d*)−*bc* where
$\mathbf{J}=\left(\begin{array}{cc}a& b\\ c& d\end{array}\right)$. 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*).

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* = *k _{x}k_{y}* such that

Since *C*_{4} 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 *C*_{4}, *i.e.*, let *T* (*r*) = *T _{N}* (··· (

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

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.*, *m _{x}* =

Fig. 5 shows the unconstrained registration and results of using *C*_{1}, *C*_{2} and *C*_{4}. 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 *C*_{1}. This shows a more regular warp than Fig. 5(a). However, *C*_{1} allows a larger solution space than the ideal solution space *C*_{0}.

Fig. 5(c) and (d) show the limitation of using a single warp based on *C*_{2} and *C*_{4} respectively. The sufficient conditions *C*_{2} and *C*_{4} 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 *C*_{4}, as shown in Fig. 6(b). In contrast, to achieve a satisfactory warp by composing transformations that lie in the box constraint *C*_{2} [21], [22] required about 30 concatenations, as shown in Fig. 6(a). For larger and more complicated deformations, our proposed constraint *C*_{4} can be used as a simple elemental transformation to provide diffeomorphic composite transformations.

If we want to strictly ensure local invertibility, then we maximize a similarity metric subject to the linear constraints *c* *C*_{4}. 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:

$$p(t;{\zeta}_{1},{\zeta}_{2})=\{\begin{array}{ll}{\scriptstyle \frac{1}{2}}{(t-{\zeta}_{1})}^{2},\hfill & t<{\zeta}_{1}\hfill \\ 0,\hfill & {\zeta}_{1}\le t\le {\zeta}_{2}\hfill \\ {\scriptstyle \frac{1}{2}}{(t-{\zeta}_{2})}^{2},\hfill & {\zeta}_{2}<t,\hfill \end{array}$$

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

A variant of quadratic penalty function (solid) and real constraints (dashed) used with constraint set *C*_{4}.

$$R(\underset{\_}{c})=\sum _{q\in \{x,y,z\}}\sum _{i,j,k}[p({c}_{i+1,j,k}^{q}-{c}_{i,j,k}^{q};{\zeta}_{1}^{q,x},{\zeta}_{2}^{q,x})+p({c}_{i,j+1,k}^{q}-{c}_{i,j,k}^{q};{\zeta}_{1}^{q,y},{\zeta}_{2}^{q,y})+p({c}_{i,j,k+1}^{q}-{c}_{i,j,k}^{q};{\zeta}_{1}^{q,z},{\zeta}_{2}^{q,z})],$$

(10)

where
${\zeta}_{1}^{q,r}=-{m}_{q}{k}_{q}$ for ∀*r* {*x, y, z*},
${\zeta}_{2}^{q,r}={m}_{q}{k}_{q}$ for *r* ≠ *q* and
${\zeta}_{2}^{q,r}={m}_{q}{K}_{q}$ 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 *C*_{4}, 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 *C*_{4} 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 *C*_{0} or *C*_{1} is much more expensive for one transformation.

For diffeomorphic transformations using Theorem 1, the usual choice would be *k _{x}* =

We applied nonrigid image registration to the 256 × 256 images in Fig. 8 using no constraint, a Jacobian penalty based on *C*_{1}, a quadratic roughness penalty [7], a regularizer based on Kim’s constraint *C*_{3}, and our proposed penalty method based on *C*_{4}. Fig. 8 has an expanding circle and a shrinking ellipse to illustrate the difference between *C*_{3} and *C*_{4}. Since we have *a priori* knowledge about vertical motion, we investigated two sets of parameters in *C*_{4}: a symmetric way with *k _{x}* =

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 *l*_{0} 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.

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 *C*_{1} contains the original constraint *C*_{0} 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: *K _{x}* =

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.

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.

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 *k _{x}* =

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 *C*_{4} which is a smaller set than *C*_{1}. Our proposed method has smoothness property implicitly because it restricts the range of the differences between adjacent B-spline coefficients.

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.

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.

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* *C*_{4}. 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.

**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.

**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.

The global maximum and minimum values exist since |**J**| is continuous on the compact set *I*_{1}×···×*I*_{9}. Suppose that (*x*_{1}, ···, *x*_{9}) achieves the global minimum value of |**J**| and min *I _{k}* <

By Lemma 1, we need to evaluate |**J**| only on *x*_{1} {−*k _{x}, K_{x}*},

For *d*(*x*) = Σ* _{i} c_{i}β^{n}*(

$$\begin{array}{l}\frac{\partial}{\partial x}d(x)=\sum _{i}{c}_{i}\frac{\partial}{\partial x}{\beta}^{n}(x/{m}_{x}-i)\\ =\sum _{i}({c}_{i}-{c}_{i-1}){\beta}^{n-1}(x/{m}_{x}-i+1/2)/{m}_{x}.\end{array}$$

Using the constraints
${b}_{m}\le {c}_{i+1,j,k}^{q}-{c}_{i,j,k}^{q}\le {b}_{M}$ and the property Σ* _{i} β^{n}*(

$$\frac{\partial}{\partial x}{d}^{q}(\underset{\_}{r})=\sum _{i,j,k}({c}_{i,j,k}^{q}-{c}_{i-1,j,k}^{q}){\beta}^{n-1}(x/{m}_{x}-i+1/2){\beta}^{n}(y/{m}_{y}-j){\beta}^{n}(z/{m}_{z}-k)/{m}_{x}\le {b}_{M}/{m}_{x}\sum _{i}{\beta}^{n-1}(x/{m}_{x}-i+1/2)\xb7\sum _{j}{\beta}^{n}(y/{m}_{y}-j)\sum _{k}{\beta}^{n}(z/{m}_{z}-k)\le {b}_{M}/{m}_{x}.$$

Similarly,
${\scriptstyle \frac{\partial}{\partial x}}{d}^{q}(\underset{\_}{r})\ge {b}_{m}/{m}_{x}$. 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 *L*_{2}-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.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |