FLIRT with Rigidity—Image Registration with a Local Non-rigidity Penalty

DOI 10.1007/s11263-007-0079-3

P O S I T I O N PA P E R

FLIRT with Rigidity—Image Registration with a Local Non-rigidity Penalty

Jan Modersitzki

Received: 24 November 2006 / Accepted: 23 July 2007 / Published online: 15 September 2007

© Springer Science+Business Media, LLC 2007

Abstract Registration is a technique nowadays commonly used in medical imaging. A drawback of most of the cur- rent registration schemes is that all tissue is being consid- ered as non-rigid (Staring et al., Proceedings of the SPIE 2006, vol. 6144, pp. 1–10,2006). Therefore, rigid objects in an image, such as bony structures or surgical instruments, may be transformed non-rigidly. In this paper, we integrate the concept of local rigidity to the FLexible Image Registra- tion Toolbox (FLIRT) (Haber and Modersitzki, in SIAM J.

Sci. Comput. 27(5):1594–1607, 2006; Modersitzki, Numer- ical Methods for Image Registration,2004). The idea is to add a penalty for local non-rigidity to the cost function and thus to penalize non-rigid transformations of rigid objects.

As our examples show, the new approach allows the main- tenance of local rigidity in the desired fashion. For example, the new scheme can keep bony structures rigid during regis- tration.

We show, how the concept of local rigidity can be in- tegrated in the FLIRT approach and present the variational backbone, a proper discretization, and a multilevel optimiza- tion scheme. We compare the FLIRT approach to the B- spline approach. As expected from the more general setting of the FLIRT approach, our examples demonstrate that the FLIRT results are superior: much smoother, smaller defor- mations, visually much more pleasing.

Keywords Image processing·Image registration· Warping·Fusion·Rigidity constraints·

Variational techniques·Constrained optimization

J. Modersitzki (

Institute of Mathematics, University of Lübeck, Wallstraße 40, 23560 Lübeck, Germany

e-mail: modersitzki@math.uni-luebeck.de

1 Introduction

1.1 Medical Image Registration

Registration is the determination of a reasonable, geometri- cal transformation that aligns points in one view of an ob- ject with corresponding points in another view of the same or a similar object. There exist many instances in imaging which demand for registration, particularly in medical imag- ing. Examples include the treatment verification of pre- and post-intervention images, the study of temporal series of car- diac images, and the monitoring of the time evolution of an agent injection subject to patient motion. Another important application is the combination of information from multiple images, acquired using different modalities, like for example computer tomography (CT) and magnetic resonance imag- ing (MRI), a technique also known as fusion. The problem of fusion and registration arises whenever images acquired from different subjects, at different times, or from different scanners need to be combined for analysis or visualization.

In the last two decades, computerized non-rigid image regis- tration has played an increasingly important role in medical imaging; see, e.g., Maintz and Viergever (1998), Glasbey (1998), Peckar et al. (1999), Fitzpatrick et al. (2000), Hajnal et al. (2001), Zitová and Flusser (2003), Yoo (2004), Mod- ersitzki (2004) and references therein.

Although registration is a generic image processing task like segmentation or enhancement, the tools to be used strongly depend on the particular application. The main dif- ficulty in registration arises from ill-posedness (Droske and Rumpf2004; Modersitzki2004): for every point in a view we give one information (grayvalue) but ask for two (or even three) unknowns (the transformation). Therefore, reg- ularization becomes a key issue. Figure1shows reference

Fig. 1 Registration results forR(c) andT (a);

transformed templateT(y)(first column), difference|T(y)−R| (second column), and

visualization ofy(third column)

and template images and visualizes three different trans- formations, obtained by weak regularization, strong regu- larization, and low regularization with additional penalty.

Although the transformations are quite different, the trans- formed images almost look the same. Without additional, external knowledge, it is impossible to rank these transfor- mations. More recent approaches to registration therefore

aim to incorporate additional pre-knowledge, as for example the correspondence of anatomical landmarks (Johnson and Christensen2002; Fischer and Modersitzki2003a,2003b), volume preservation of tissue (Rohlfing et al.2003; Haber and Modersitzki2004), limitation of volume changes (Haber and Modersitzki 2006b), or local rigidity (Staring et al.

2006; Modersitzki2007).

1.2 Local Rigidity

Most of the current registration schemes (except for Little et al.1997) treat all displayed tissue either as globally rigid or globally non-rigid (Staring et al.2006). Therefore, rigid ob- jects in an image, such as bony structures or surgical instru- ments, may be non-rigidly deformed. Other consequences are that tumor growth between follow-up images may be concealed, or that structures containing contrast material in only one of the images may be compressed by the registra- tion scheme. In order to improve the registration results, we thus aim to integrate a concept of local rigidity as illustrated in Fig.1. Reference and template images provide a natural segmentation in two gray squares on a white background.

The new objective is to find a transformation maintaining the inner structure of the segments, i.e., being locally rigid.

In a medical application, the concept of local rigidity ap- plies to perfection to the registration of bones but can also lead to more realistic results for the generic case of images displaying soft and hard tissues.

1.3 State of the Art

There are currently two main numerical approaches to im- age registration. The first one is based on an expansion of the wanted transformation in terms of B-splines (Rueckert et al.

1999) and the other one is based on the more general varia- tional framework (Modersitzki2004) (FLexible Image Reg- istration Toolbox (FLIRT)). Both approaches principally al- low for the integration of additional pre-knowledge in terms of a penalty, like, e.g. local rigidity. For a B-spline or an even more general radial basis function approach this has been implemented in (Little et al.1997; Loeckx et al.2004;

Staring et al.2006), while the objective of this paper is the integration of the local rigidity concept into FLIRT.

In contrast to (Keeling and Ring2005), where the regu- larizer is constructed to penalize departure from rigidity on the global average, our approach penalizes departure from rigidity on average over locally selected regions; see in par- ticular Fig. 9 in (Keeling and Ring2005). Here, we do not use a spatially variant regularization parameter (e.g., local elasticity parameters like Kabus et al. 2006) but explicitly penalize local non-rigidity.

1.4 FLIRT

The FLexible Image Registration Toolbox (FLIRT) (Fis- cher and Modersitzki2003c; Modersitzki2004, and in par- ticular Haber and Modersitzki 2004, 2006a, 2006b) pro- vides a general and efficient framework for image regis- tration. This approach treats registration as an optimization problem where the objective is to minimize a joint func- tionalJ with respect to the wanted transformationy. This

functional is assembled from different building blocks: a distance measure D quantifying image similarity, a regu- larizerS quantifying reasonability of the transformation, a penalty P penalizing unwanted transformations, and even additional hard constraints, if wanted. Our numerical ap- proach uses an efficient multilevel strategy which combines different discretizations of the continuous problem. On a fixed level, a Gauss-Newton type optimization scheme is used, where a multigrid scheme serves as a solver for the linear system.

1.5 Outline

In this paper, we integrate the concept of local rigidity to the FLexible Image Registration Toolbox (FLIRT) (Fis- cher and Modersitzki 2003c; Modersitzki2004). The idea is to add a penalty for local non-rigidity to the cost func- tion and thus to penalize non-rigid transformations of rigid objects. In Sect. 2 we briefly summarize the FLIRT ap- proach and present a variational setting for local rigidity, whereas Sect.3summarizes the discretization of the varia- tional FLIRT approach and presents our discretization of lo- cal rigidity. The optimization scheme is explained in Sect.4.

In Sect.5we present numerical results. We demonstrate that the new approach enables the maintenance of local rigid- ity. Moreover, we show that the FLIRT approach is superior to the B-spline approach (as expected from the additional smoothing properties of B-splines). We finally show the per- formance in a medical setting: in an X-ray application we keep a middle finger of a hand rigid and in a MRI applica- tion we keep the bones in a leg rigid.

2 Mathematical Framework

We briefly summarize the variational framework of the FLIRT approach (see Haber and Modersitzki2004,2006a, 2006bfor details). Moreover, we also present a functional setting for the local rigidity approach. The basic idea is to measure rigidity (by means of linearity, orthogonality and orientation preservation of the transformation) at every spa- tial position and to finally weight this measure by a moving weighting functionQ(y).

2.1 Summarizing FLIRT

The FLIRT approach is based on the minimization of the joint functional

J(y):=D(T(y),R)+αS(y−y^kern)+βC(y)

(1) subject to y∈M.

Here, y=(y1, . . . , yd):R^d→R^d is the wanted transfor- mation,R,T :→Rare the given reference and template images, andT(y)is the transformed image, i.e.,

T(y)(x):=T(y(x)) for allx∈, (2) where ⊂R^d is a given domain (typically =]0,1[^d) anddis the image dimensionality (typicallyd=2,3).

The functionalDmeasures the distance between the im- ages and could be the sum of squares difference (SSD), mu- tual information (Collignon et al.1995; Wells et al.1996), normal gradient field (Pluim et al.2000; Haber and Moder- sitzki2005), or any other distance measure at hand. In this paper, we focus on SSD,

D(T(y),R)=1

2T(y)−R²_L2(). (3)

The regularizerS is designed to yield a unique transforma- tion y. Common choices are summarized in (Modersitzki 2004). Here, we use the most common so-called elastic po- tential

S(y)=ω 2

j

∇y^j2_L2()+λ+ω

2 ∇ ·y²_L2(), (4) with Lamé-constantsλ, ω; see Modersitzki (2004) for de- tails. The regularization parameter α >0 compromises be- tween similarity and regularity. Note that the regularization in (1) is applied to the displacement y −y^kern which al- lows an elegant integration of any pre-registration. The con- straintsC andMare used to supply additional application conform pre-knowledge. Particular choices for local rigidity are discussed in the next subsection.

2.2 Local Rigidity

The goal of this subsection is to formalize and facilitate local rigidity. The naive approach would be based on a segmen- tation and may lead to localization errors in a multi-level framework. We therefore prefer a weighted approach, where the penalty is computed everywhere and then weighted by a differentiable function Q:→R; see Fig. 4j and Fig.5j for examples. The idea is to segment structures that are sup- posed to transform rigidly and to pickQas a smoothed char- acteristic function for the a priori given segments. In the ex- amples shown, a manual segmentation has been performed andQis based on a linear interpolation scheme. Note that the rigid transformations on different segments can be differ- ent as it is shown in the knee example. A second advantage of the weighted approach is that it easily allows for a correct location of the penalty in a multilevel framework. Since we like to maintain the rigidity of a displayed structure and not of a fixed subset of the domain, the weight has to be trans- formed in the same way as the structure. As a consequence, the appropriate weight isQ(y), whereQ(y)(x)=Q(y(x)).

We start with the Jacobian ofy,

∇y=

⎛

⎝

∂₁y₁ . . . ∂_dy₁ ... · · · ...

∂1yd . . . ∂dyd

⎞

⎠, (5)

where natural boundary conditions associated with∇y are used, i.e., that no additional artificial boundary condition constraints are imposed. A formal definition of rigidity reads as follows, see also Gurtin (1981), Staring et al. (2006).

Definition 1 A transformationy∈C²(,R^d)is rigid, if it is linear (∂i1,i2y_j=0 for alli₁, i₂, j=1, . . . , d), orthogonal (∇y∇y=I_d), and orientation preserving (det∇y=1) for allx∈.

Lemma 1 A transformationy∈C²(,R^d)is rigid, if and only ify(x)=Qx+b whereQ∈R^d,d is orthogonal and orientation preserving, andb∈R^d.

Proof From∂_i1,i2y_j=0 it follows that the transformation is linear, i.e.y(x)=Ax+band∇y(x)=A, from∇y∇y= Id it follows that Ais orthogonal and from∇y∇y =Id

orientation preservation follows. The opposite direction is

obvious.

With this definition, we haveM=C²(,R^d)and the penaltyCin (1) is given by

C(y)= d i1,i2=1,i2≤i1

1

2∂_i1,i2y_j²_Q(y) +

d i1,i2=1,i₂≤i1

1

2(∇y∇y−Id)i₁,i₂²_Q(y) +1

2det∇y−1²_Q(y), (6)

with the weightedL2-norm f²_Q(y)=

Q(y(x))²f (x)²dx. (7) The first term in (6) measures departure from linearity, the second departure from orthogonality, and the third departure from orientation preservation; the weighted norm puts a fo- cus on the moving structures to be kept rigid.

Particularly ford=2, and · = · _Q(y)we have 2C(y)= ∂_1,1y₁²+ ∂_1,2y₁²+ ∂_2,2y₁²

+ ∂1,1y2²+ ∂1,2y2²+ ∂2,2y2² + (∂1y1)²+(∂2y1)²−1²

+ (∂1y2)²+(∂2y2)²−1²

+ (∂1y1)(∂1y2)+(∂2y1)(∂2y2)²

+ (∂₁y₁)(∂₂y₂)−(∂₂y₁)(∂₁y₂)−1², (8) and ford=3 the formula is along the same lines but lengthy.

3 Discretization

Our numerical approach focusses on a discretize then opti- mize strategy, see Haber and Modersitzki (2004), Haber and Modersitzki (2006a) for details. The discretization is based on pixel/voxel (see Fig.2). An important issue is to guaran- tee anh-elliptic discretization for the regularizer. We there- fore use staggered grids for the discretization of the vector fieldy. We describe the discretization ford=2, the exten- sion tod =3 is straightforward; more details are given in Haber and Modersitzki (2004). In Sect.3.1we briefly sum- marize the FLIRT discretization (Haber and Modersitzki 2004,2006a) and in Sects.3.2and3.3we explain the dis- cretization of the rigidity penalty in detail.

3.1 Discretizing the Images, the Distance, and the Regularizer

We briefly summarize the discretization suggested in Haber and Modersitzki (2004) for dimension d = 2. For the discretization of y we assume that =]0, ω1[×]0, ω2[ and m=(m1, m2) are given. For any dimension j, we seth_j=ω_j/m_j and end up with the following grids (see also Fig.2):

x^h,•=(i1h1, i2h2)i_j=1/2:m_j−1/2, cell-centered, x^h,=(i1h1, i2h2)i1=0:m1,i2=1/2:m2−1/2, face staggered, x^h,=(i₁h₁, i₂h₂)_{i1=1/2:m1−1/2,i2=0:m2}, face staggered.

We discretizey=(y1(x1, x2), y2(x1, x2))by values on the staggered grids, i.e. y₁^h =y¹(x^h,) and y^h₂ =y¹(x^h,), where for convenience, we assemble the unknowns in lex- icographical order, i.e. y₁^h∈Rⁿ¹ andy₂^h ∈Rⁿ² withn₁= (m1+1)m2,n2=m1(m2+1), respectively.

The ultimate unknown is the vector y^h = ((y₁^h), (y₂^h))∈Rⁿ, withn=n1+n2. Explicit formulae for the

discretization of T, R, D and S are given in Haber and Modersitzki (2004) and are summarized as follows:

R=R(x^h,•)∈Rⁿ, T (y^h)=T(P y^h)∈Rⁿ, D(T (y^h), R)=1

2T (y^h)−R², S(y^h)=1

2By^h2, where

P=

P1 0 0 P2

, By^h2=λ+ω

2 ∇^h·y^h2+ω 2

d j=1

∇_j^hy_j^h2,

∇^h· =(∂₁^h,1, ∂₂^h,2), ∇_j^h=((∂₁^h,j), (∂₂^h,j)),

and the matricesP1,P2, and∂_k^h,j are specified in the next section.

3.2 Discretizing Derivatives

For the regularizer and the rigidity constraints, we need to discretize first and second order derivatives of y. The goal and the reason for staggered grids is to use short dif- ferences for the approximations of ∂_jy_j, which are then located at the cell-centered grid. Unfortunately, all other derivatives involve averaging operators that keep a formal description lengthy. Our description is based on the band matrices and Kronecker-products ⊗ (Brewer 1978). Let band(a₋k, . . . , a_k;k₁, k₂) denotes a k₁-by-k2 band matrix with diagonal bands a₋k, . . . , ak, wherea0 is on the main diagonal. For example,

band(0,−1,1;3,4)=

⎛

⎝−1 1 0 0

0 −1 1 0

0 0 −1 1

⎞

⎠.

Fig. 2 Staggered grid discretization ofy

With these matrices, we set Pj=band(0,1,1;mj, mj+1)/2,

D^short,j=band(0,−1,1;m_j, m_j+1)/ hj, D^long,j=BCband(−1,0,1;m_j, m_j+1)/(2hj), D^short,j,j=band(1,−2,1;m_j, m_j)/(h²_j),

D^long,j,j=BCband(0,1,−1,−1,1;m_j, m_j+1)/(2h²_j), where “=BC” means “up to boundary conditions”. Partic- ularly, the boundary conditions are chosen such that dis- cretizations of constants (for first order derivatives) and affine linear (second order derivatives) belong to the kernel.

The first and second order derivatives (linear parts in (6)) are approximated by

∂1y1(x^h,•)≈∂₁^h,1y₁^h, ∂₁^h,1=Im₂⊗D^short,1, (9a)

∂2y1(x^h,•)≈∂₂^h,1y₁^h, ∂₂^h,1=D^long,2⊗P1, (9b)

∂₁y₂(x^h,•)≈∂₁^h,2y₂^h, ∂₁^h,2=P₂⊗D^long,1, (9c)

∂₂y₂(x^h,•)≈∂₂^h,2y₂^h, ∂₂^h,2=D^short,2⊗I_m1, (9d)

∂1,1y1(x^h,•)≈C1(y^h)=∂_1,1^h,1y₁^h,

∂_1,1^h,1=Im₂⊗D^long,1,1, (10a)

∂1,2y1(x^h,•)≈C2(y^h)=∂_1,2^h,1y₁^h,

∂_1,2^h,1=D^long,2⊗D^short,1, (10b)

∂_2,2y₁(x^h,•)≈C₃(y^h)=∂_2,2^h,1y₁^h,

∂_2,2^h,1=D^short,2,2⊗P1, (10c)

∂1,1y2(x^h,•)≈C4(y^h)=∂_1,1^h,2y₂^h,

∂_1,1^h,2=P2⊗D^short,1,1, (10d)

∂1,2y2(x^h,•)≈C5(y^h)=∂_1,2^h,2y₂^h,

∂_1,2^h,2=D^short,2⊗D^long,1, (10e)

∂_2,2y₂(x^h,•)≈C₆(y^h)=∂_2,2^h,2y₂^h,

∂_2,2^h,2=D^long,1,1⊗I_m1, (10f)

and the non-linear part (∇y∇y −I_d) and (det∇y−1) in (6) are discretized by

C7(y^h)=(∂₁^h,1y₁^h)(∂₁^h,1y^h₁)

+(∂₂^h,1y₁^h)(∂₂^h,1y₁^h)−e, (10g) C8(y^h)=(∂₁^h,1y₁^h)(∂₁^h,2y^h₂)

+(∂₂^h,1y₁^h)(∂₂^h,2y₂^h), (10h) C₉(y^h)=(∂₁^h,2y₂^h)(∂₁^h,2y^h₂)

+(∂₂^h,2y₂^h)(∂₂^h,2y₂^h)−e, (10i)

C10(y^h)=(∂₁^h,1y₁^h)(∂₂^h,2y₂^h)

−(∂₂^h,1y₁^h)(∂₁^h,2y₂^h)−e, (10j) where e=(1, . . . ,1) and denotes the Hadamard- or pointwise product ((xy)i=xiyi).

3.3 Discretizing Local Rigidity

Using the discretizations of the previous subsection, we dis- cretize the penaltyCby anl2-norm of an appropriate chosen residualr:

C(y^h)=1

2r(y^h)r(y^h), r=(r₁, . . . , r₁₀),

andrj(y^h)=Q(P y^h)Cj(y^h). (11) For the derivative ofC, we thus get

dC(y^h)=r(y^h)dr(r^h), (12)

where the derivativedr_j is given by drj(y^h)=diag(Q(P y^h))dCj(y^h)

+diag(Cj(y^h))dQ(P y^h)P . (13) Note that for j =1, . . . ,6, C_j is linear anddC_j follows directly from (9). The parts C7, . . . , C10 are quadratic and their derivatives are given by

dC₇(y^h)=(2diag(∂₁^h,1y₁^h)∂₁^h,1+2diag(∂₂^h,1y₁^h)∂₂^h,1,0) (14a) dC8(y^h)=(diag(∂₁^h,2y^h₂)∂₁^h,1+diag(∂₂^h,2y₂^h)∂₂^h,1,

diag(∂₁^h,1y^h₁)∂₁^h,2+diag(∂₂^h,1y₁^h)∂₂^h,2), (14b) dC₉(y^h)=(0,2diag(∂₁^h,2y₂^h)∂₁^h,2+2diag(∂₂^h,2y₂^h)∂₂^h,2),

(14c) dC10(y^h)=(diag(∂₂^h,2y₂^h)∂₁^h,1−diag(∂₁^h,2y₂^h)∂₂^h,1,

diag(∂₁^h,1y₁^h)∂₂^h,2−diag(∂₂^h,1y₁^h)∂₁^h,2). (14d) 3.4 Image Registration Discretized

To summarize, our discrete version of (1) reads:

minimizeJ^h(y^h)=1

2T (y^h)−R² +α

2B(y^h−y^kern,h)²+β

2r(y^h)², (15) with analytic gradient

dJ^h(y^h)=(T (y^h)−R)dT (y^h)

+αBB(y^h−y^kern,h)+βr(y^h)dr(y^h).

(16)

Algorithm 1 (Multilevel Image Registration) 1. for level=0:Ldo

2. h=h^(level);

3. if(level=0), % coarsest level

4. perform rigid pre-registration and obtainy⁽⁰⁾=y^rigid;

5. else

6. sety⁽⁰⁾=I_h^2hy^2h,∗;

7. endIf

8. compute minimizery^h,∗ofJ^h(15) usingy⁽⁰⁾as starting guess 9. endFor

10. returny^h,∗.

4 Optimization

For the numerical optimization of (1) we use the techniques summarized in Haber and Modersitzki (2006a). The main ingredients are a multilevel strategy and a Gauss-Newton type optimizer for a fixed level. In the Gauss-Newton ap- proach the Hessian is approximated by a sum of the approx- imations tod²Dandd²C and the exact Hessiand²Sof the regularizer.

4.1 Multilevel

In order to compute a solution for the registration prob- lem (1), we use a hierarchy of different discretizations h⁽⁰⁾, . . . , h^(L), where h⁽⁺¹⁾ =2⁻h⁽⁰⁾, see Algorithm 1.

Starting with the coarsest grid=0, we initially perform a rigid pre-registration yieldingy^kernand compute a numer- ical minimizery^h,∗of J^h, whereh=h⁽⁰⁾ andy=y^kern,h serves as a starting guess. If we are not yet on the finest level ( < L), we switch to a higher level (←+1). Starting withy=I_h^2hy^2h,∗, whereI_h^2h is a standard multigrid inter- polation scheme (Trottenberg et al.2001), we compute a nu- merical minimizery^h,∗of J^h, whereh=h⁽⁾. These steps are repeated until we reach the finest resolutionh=h^(L). 4.2 Numerical Optimization ofJ^h

It remains to describe our numerical scheme for the mini- mizing ofJ=J^h. We use a Gauss-Newton type approach where the Hessian is approximated as follows

d²J (y)≈diag(dT (y)dT (y))

+αBB+βdr(y)dr(y), (17) see Algorithm2. Note that for the regularizerSthe second order derivativesBBis used, while for the distanceDand the penaltyC only first order derivatives are used. Though the dimensions ofH can be dramatically large (particularly for d =3), it has only a very small number of non-zero entries per row. As a consequence, multiplication by H is

Algorithm 2 (MinimizingJ=J^h(15)) 1. fork=1, . . .do

2. computeT (y),dT (y),D= T (y)−R²/2, dD=(T (y)−R)dT,

S= B(y−y^kern)²/2,dS=BB(y−y^kern), C= r(y)²/2,dC=r(y)dr(y),

computeJ=D+αS+βC, computedJ=dD+αdS+βdC;

3. check stopping criteria (see, e.g., Gill et al.1981) 4. solveH δ_y= −dJ, where

H=diag(dTdT )+αBB+βdrdr,

5. use Armijo line search to compute a “minimizer”t^∗of φ (t )=J (y+t δ_y),

6. updatey←y+t^∗δ_y; 7. endFor

8. returny.

anO(n)operation and efficient multigrid techniques might be adapted. However, in our current 2D MATLAB (Math- Works1992) implementation, we simply use the backslash operator.

5 Results

We tested our implementation on a variety of examples. Due to page limitations, we can only present a few intuitive and representative examples.

5.1 One Square

Our first example is a repetition of the experiment performed in (Staring et al. 2006), from which we also took the re- sults for the B-spline approaches, see Fig.3. From these re- sults, we see the effect of local rigidity constraints placed at non-zero locations in the moving template itself (Q=T).

The B-spline implementation uses a backward interpolation scheme, while the FLIRT implementation uses a forward scheme (in fact, for the FLIRT registration we interchanged reference and template in order to make the grids compara- ble).

As expected, both penalized approaches keep the square rigid. However, the B-spline approach introduces large de- formations in some parts of the image. This is partly due to the inappropriate vanishing Dirichlet boundary conditions (BC) but may also be related to the locations of the control points.

FLIRT takes advantage of a regularization kernel which handles global transformations. Moreover, FLIRT involves vanishing natural second order Neumann boundary condi- tions. As a consequence, the FLIRT result—with and with- out penalty—is a global rigid transformation leading to a

Fig. 3 Example from Staring et al. (2006): reference and template (first column), B-spline results taken from Staring et al.

(2006) without and with penalty (middle column), and FLIRT results without and with penalty (right column); all four transformations lead to a perfect match

perfect copy of the reference. This is probably also the most intuitive solution to this toy problem. The results are not shown here. In order to show the effects of the additional penalty term, we artificially removed the kernel information and these results are shown in Fig.3.

All results (B-spline/Flirt, without/with penalty) yield a perfect match, i.e. D(y) ≈0. However, we find the FLIRT transformations superior since they are visually much smoother, more local, and introduce less distortions.

5.2 Two Squares

Our next example shows registration results for images with two rigid objects which has already been introduce in the introduction; see Fig.1. In this exampleQ=T.

Our first experiment is done with small regularization and without penalty (α=1, β =0). As a result, the transfor- mation is heavily non-linear deformed and does not look appealing, see Fig. 1f. Therefore, the second experiment is done with large regularization but again without penalty (α=10,β=0). As expected, the resulting transformation is much smoother (Fig.1i) and looks more pleasing. However, especially the big square is transformed non-linearly and in addition, the transformed image is not even close to the reference. The third example use the local rigidity penalty (α=1,β=10⁻⁴). As a result, we get a very pleasing trans- formation (Fig.1l) with a perfect match (Fig.1j).

5.3 X-ray of Hands

A more realistic but still intuitive example is presented in Fig. 4 (see also Amit1994; Modersitzki 2004). Note that

the template image shows a global rotation of approxi- mately 25 degrees which outrules the B-spline approach with Dirichlet BC. In this example we make the middle finger of the hand to be rigid (see Fig. 4j), i.e. Q is a smoothed version of the bone in the middle finger. Figure4 shows FLIRT results without (β=0) as well as with penalty (β=0.01). For both variants, we pickedα=500. As it ap- parent, the penalized approach does keep the finger rigid while the unconstrained does not; see particularly the maps of det(∇y)(see Fig.4m, n).

5.4 Two Knees

This example shows the registration of computer tomogra- phy images from a human knee in a bent and straight po- sition (image courtesy of T. Netsch, Philips Research Ham- burg, Germany). Here, we pick two 2D subsections of the 3D reference and template images, respectively and perform the registration without (α=200,β=0) and with penalty (α=100, β=10⁻⁸; the small value ofβ is related to the height weight in this example: max{Q(x)} =255).

For the approach without penalty, a large valueα=200 is needed to keep the transformation one-to-one, particu- larly in the neighborhood of the joint. As it is apparent from Fig.5b, even with this strong regularization the femur is in- tolerably deformed. For the penalized approach, we find a much smoother overall transformation (even with a smaller regularization parameterα=100), that leads to a realisti- cally deformed template (Fig. 5c). The maps of det(∇y) (Fig.5m, n) clearly demonstrate the benefits of local rigidity for this application.

Fig. 4 Registration results forR, (g) andT, (a): data (first column), results forα=500 andβ=0 (second column) andα=500 andβ=0.01 (third column)

6 Conclusions

Registration is a technique nowadays commonly used in medical imaging. A drawback of most of the current reg-

istration schemes is that all tissue is being considered as non-rigid (Staring et al.2006). Therefore, rigid objects in an image, such as bony structures or surgical instruments, may be transformed non-rigidly. In this paper, we integrate

Fig. 5 Registration results forR, (g) andT, (a): data (first column), results forα=200 andβ=0 (second column) andα=100 andβ=10⁻⁸ (third column)

the concept of local rigidity to the FLexible Image Registra- tion Toolbox (FLIRT) (Modersitzki2004; Haber and Mod- ersitzki2006a). The idea is to add a penalty for local non- rigidity to the cost function and thus to penalize non-rigid

transformations of rigid objects. As our examples show, the new approach allows the maintenance of local rigidity in the desired fashion. For example, the new scheme can keep bony structures rigid during registration. We show, how the

concept of local rigidity can be integrated in the FLIRT ap- proach and present the variational backbone, a proper dis- cretization, and a multilevel optimization scheme. We com- pare the FLIRT approach to the B-spline approach. As ex- pected from the more general setting of the FLIRT approach, our examples demonstrate that the FLIRT results are supe- rior: much smoother, smaller deformations, visually much more pleasing.

Next steps are the automatization of the parameter choices for regularization and penalty and of course a 3D implementation. The FLIRT approach principally allows an integration hard constraints (like, e.g., landmarks Fischer and Modersitzki 2003bor volume preservation Haber and Modersitzki2004). However, the integration of local rigid- ity as a hard constrained requires further analysis.

References

Amit, Y. (1994). A nonlinear variational problem for image matching.

SIAM Journal on Scientific Computing, 15(1), 207–224.

Brewer, J. W. (1978). Kronecker products and matrix calculus in system theory. IEEE Transactions on Circuits and Systems, 25, 772–780.

Collignon, A., Vandermeulen, A., Suetens, P., & Marchal, G. (1995).

3d multi-modality medical image registration based on infor- mation theory. In Computational Imaging and Vision (Vol. 3, pp. 263–274).

Droske, M., & Rumpf, M. (2004). A variational approach to non-rigid morphological registration. SIAM Journal on Applied Mathemat- ics, 64(2), 668–687.

Fischer, B., & Modersitzki, J. (2003a). Combination of automatic non-rigid and landmark based registration: the best of both worlds. In M. Sonka & J. M. Fitzpatrick (Eds.), Proceedings of the SPIE: Vol. 5032. Medical imaging 2003: Image processing (pp. 1037–1048).

Fischer, B., & Modersitzki, J. (2003b). Combining landmark and in- tensity driven registrations. PAMM, 3, 32–35.

Fischer, B., & Modersitzki, J. (2003c). FLIRT: a flexible image regis- tration toolbox. In J. C. Gee, J. B. A. Maintz, & M. W. Vannier (Eds.), Lecture notes in computer science: Vol. 2717. 2nd interna- tional workshop on biomedical image registration 2003 (pp. 261–

270). Berlin: Springer.

Fitzpatrick, J. M., Hill, D. L. G., & Maurer Jr., C. R. (2000). Image registration. In M. Sonka, & J. M. Fitzpatrick (Eds.), Handbook of medical imaging: Vol. 2. Medical image processing and analysis (pp. 447–513). Bellingham: SPIE.

Gill, P. E., Murray, W., & Wright, M. H. (1981). Practical optimization.

London: Academic Press.

Glasbey, C. (1998). A review of image warping methods. Journal of Applied Statistics, 25, 155–171.

Gurtin, M. E. (1981). An introduction to continuum mechanics. Or- lando: Academic Press.

Haber, E., & Modersitzki, J. (2004). Numerical methods for volume preserving image registration. Inverse Problems, 20(5), 1621–

1638.

Haber, E., & Modersitzki, J. (2005). Beyond mutual information:

A simple and robust alternative. In H. P. Meinzer, H. Handels, A. Horsch, & T. Tolxdorff (Eds.), Bildverarbeitung für die Medi- zin 2005 (pp. 350–354). Berlin: Springer.

Haber, E., & Modersitzki, J. (2006a). A multilevel method for image registration. SIAM Journal on Scientific Computing, 27(5), 1594–

1607.

Haber, E., & Modersitzki, J. (2006b). Image registration with a guaran- teed displacement regularity. International Journal of Computer Vision, 1 DOI:10.1007/s11263-006-8984-4

Hajnal, J., Hawkes, D., & Hill, D. (2001). Medical image registration.

Boca Raton: CRC Press.

Johnson, H. J., & Christensen, G. E. (2002). Consistent landmark and intensity-based image registration. IEEE Transactions on Medical Imaging, 21(5), 450–461.

Kabus, S., Franz, A., & Fischer, B. (2006). Variational image registra- tion with local properties. In F. A. Gerritsen, J. P. W. Pluim, &

B. Likar (Eds.), Lecture notes in computer science. Third inter- national workshop, WBIR 2006: Biomedical image registration (pp. 92–100). Berlin: Springer.

Keeling, S. L., & Ring, W. (2005). Medical image registration and interpolation by optical flow with maximal rigidity. Journal of Mathematical Imaging and Vision, 23(1), 47–65.

Little, J. A., Hill, D. L. G., & Hawkes, D. J. (1997). Deformations in- corporating rigid structures. Computer Vision and Image Under- standing, 66(2), 223–232.

Loeckx, D., Maes, F., Vandermeulen, D., & Suetens, P. (2004). Non- rigid image registration using free-form deformations with a lo- cal rigidity constraint. In Lecture notes in computer science:

Vol. 3216. MICCAI (pp. 639–646).

Maintz, J. B. A., & Viergever, M. A. (1998). A survey of medical image registration. Medical Image Analysis, 2(1), 1–36.

MathWorks (1992). Matlab user’s guide. Natick.

Modersitzki, J. (2004). Numerical methods for image registration. Ox- ford: Oxford University Press.

Modersitzki, J. (2007). Image registration with local rigidity con- straints. In H. Handels, H. P. Meinzer, A. Horsch, T. M. Deserno,

& T. Tolxdoff (Eds.), Informatik Aktuell: Bildverarbeitung für die Medizin (pp. 444–448). Berlin: Springer.

Peckar, W., Schnörr, C., Rohr, K., & Stiehl, H. S. (1999). Parameter- free elastic deformation approach for 2d and 3d registration using prescribed displacements. Journal Mathematical Imaging and Vi- sion, 10(2), 143–162.

Pluim, J. P. W., Maintz, J. B. A., & Viergever, M. A. (2000). Image registration by maximization of combined mutual information and gradient information. IEEE TMI, 19(8), 809–814.

Rohlfing, T., Maurer Jr., C. R., Bluemke, D. A., & Jacobs, M. A.

(2003). Volume-preserving nonrigid registration of MR breast im- ages using free-form deformation with an incompressibility con- straint. IEEE Transactions on Medical Imaging, 22(6), 730–741.

Rueckert, D., Sonoda, L., Hayes, C., Hill, D., Leach, M., & Hawkes, D.

(1999). Non-rigid registration using free-form deformations.

IEEE Transactions on Medical Imaging, 18(1), 712–721.

Staring, M., Klein, S., & Pluim, J. P. W. (2006). Nonrigid registration using a rigidity constraint. In J. M. Reihnardt, & J. P. W. Pluim (Eds.), Proceedings of the SPIE 2006: Medical imaging, 2006 (Vol. 6144, pp. 1–10). Boca Raton: SPIE.

Trottenberg, U., Oosterlee, C., & Schüller, A. (2001). Multigrid. Lon- don: Academic Press.

Wells III, W. M., Viola, P., Atsumi, H., Nakajima, S., & Kikinis, R.

(1996). Multi-modal volume registration by maximization of mu- tual information. Medical Image Analysis, 1(1), 35–51.

Yoo, T. S. (2004). Insight into images: Principles and practice for seg- mentation, registration, and image analysis. AK Peters Ltd.

Zitová, B., & Flusser, J. (2003). Image registration methods: a survey.

Image and Vision Computing, 21(11), 977–1000.