# Constrained Linear Curvature Image Registration Model and Its Numerical Algorithm

Authors: Jin Zhang
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

We propose a constrained linear curvature image registration model to explicitly control the deformation according to the transformed Jacobian matrix determinant using point-by-point inequality constraints in this paper. In addition, an effective numerical method is proposed to solve the resulting inequality constrained optimization model. Finally, some numerical examples are given to prove the obvious advantages of the curvature image registration model with inequality constraints.


---

## Body

## 1. Introduction

In image processing, people are interested not only in analyzing an image but also comparing or combining information from images which take different time, different places, different viewpoints, or different modalities. Thus, image registration is one of the most useful and challenging problems in the field of image processing. Its main idea is to find a geometric transformation which aligns points in one view of one object with corresponding points in another view of the same or similar object. At present, there are a large number of application areas which require image registration, such as computer vision, biological imaging, remote sensing, and medical imaging. For comprehensive surveys of these applications, refer to [1–5].The basic framework of image registration can be described as follows: given two images of the same object, which are called reference imageR and template image T, respectively, and our purpose is to find a vector value transformation φ as defined below:(1)φu⋅:ℝd⟶ℝd,φux:x⟶x+ux,or equivalently find the unknown displacement field u:(2)u:ℝd⟶ℝd,u:x⟶ux=u1x,u2x,…,udx⊤,so that the transformed template image Tφx=Tx+ux≜Tu is as similar to the reference image R as possible. Here, d∈ℕ denotes spatial dimension of the given images.Without loss of generality, here we focus ond=2 throughout this paper, but it is easy to generalize to d=3 with some additional modifications. The variational model is an important tool for studying image registration and has been widely concerned by many researchers [5–9]. This variational model treats the image registration problem as a minimization problem of the joint energy functional in the following form:(3)minuJαu=Du+αSu,where(4)Du=12∫ΩTu−R2dΩ,where Du denotes distance measure which quantifies distance or similarity of the transformed template image Tu and reference R, for other choices on Du, refer to [5, 7], Su is the deformation regularizer which constrains u and ensures the well-posedness of the problem, and α>0 is a regularization parameter which balances similarity and regularity of displacement.We all know that different regularizers will produce displacement fields with different degrees of smoothness and the selection of a regularizer is critical to the solution of the problem and its properties; for more details, refer to [5]. Usually the choice of regularizer can be classified into two main categories: the first type is to limit the displacement field ux to the parametric model [10–14], for example, rigid or affine transformations (parameterized by rotation, scaling, and translation) or linear combinations of a set of basis functions (B-splines) [3, 5, 15–18]; the second type is based on the derivative of the displacement field. At present, there are regularizers based on first-order derivatives, such as elastic regularizer [19–21], diffusion regularizer [22], total variational regularizer [23, 24], modified total variational regularizer [25, 26], total fractional-order regularizer [27], and the ones based on higher-order derivatives, such as linear curvature [28, 29], mean curvature [8, 30], and Gaussian curvature [31]. It has been proved that, in many cases, the selection method of the first kind of a regularizer is too strict, and the required transformation cannot be guaranteed to be included in the parametric model. Therefore, the second kind of method is a common method to select a regularizer. For the second method, it is easy to implement for low-order regularizers, while they are less effective than high-order ones in producing smooth displacement fields which are important in some applications including medical imaging. Although the registration models based on a higher-order regularizer can produce more satisfactory registration results visually, they do not take into account mesh folding.In fact, the regularity of the displacement field is also an important measure in image registration [32]. In many variational models (1) that currently exist, although they can produce satisfactory registration results visually, they cannot ensure that the transformation φx=x+ux found is reversible. The irreversibility of the transformation means that the displacement field is not regular. In this case, there will be mesh folding during the registration process, which is not allowed in practical applications. Therefore, it is necessary to avoid mesh folding during the registration process. Currently, a direct idea to avoid mesh folding is to use a larger regularization parameter α. However, such a value will cause the similarity between the transformed template image and the reference image to become worse. In order to avoid mesh folding phenomenon, some scholars have proposed to add an additional regular term Cu on the transformed Jacobian matrix determinant in the objective function formula (3) [33–36], i.e.,(5)minuJαu=Du+αSu+βCu−12,where Cu=detId+∇u represents the determinant of the Jacobian matrix of the transformation. However, this method only penalizes the irregular displacement field as a whole, while the local displacement field cannot be guaranteed to be regular [32]. In addition, this method is only effective for the smaller regularizer parameter β, and increasing the value of β usually leads to ill-posed optimization problems [37]. To solve this problem, Haber and Modersitzki proposed a new registration model by adding additional explicit volume inequality constraints [32]; however, this constrained method usually leads to solving a large-scale highly nonlinear inequality constrained optimization problem. Other methods to ensure the regularity of the displacement field can be found in the literature [20, 38–43]. However, some of them require more computation time due to the complexity of the regularizer.There are two purposes for image registration. One is to enhance some similarities between two images by geometrically transforming one of the given two images. The other is to ensure that this transformation is reasonable. In fact, it is equivalent to find geometric transformationφ and the displacement field ux in the framework of variational model. If the displacement field is irregular, the transformation is considered unreasonable, and then the mesh folding phenomenon will appear which is not allowed in practical applications. In this paper, we propose a new image registration model by integrating the evaluation criteria to measure the registration results directly into the basic framework of the variational model (3).The rest of the paper is organized as follows: in Section2, we propose a new constrained linear curvature image registration model. Then in Section 3, we discuss the numerical method for solving the new model by using a combination of the multiplier method and Gauss–Newton scheme with the Armijos line search and further combine with a multilevel method to achieve fast convergence. Next, some experimental results from syntectic and real images are illustrated in Section 4. Finally, conclusions and future work are summarized in Section 5.


## 2. Constrained Linear Curvature Image Registration Model

Firstly, we briefly review the Fischer–Modersitzki’s linear curvature image registration model [25, 27]. Choosing Su in (3) based on an approximation to the curvature of the surface of the displacement field ul is given by the following form:(6)SLCu=12∑l=12∫ΩΔul2dΩ.There are two major advantages to the particular choice of the regularizer. Firstly, it can penalize oscillations; secondly, without requiring an additional affine linear preregistration step, it can produce visually more satisfactory registration results than a diffusion model and an elastic model for smooth displacement fields. However, a mesh folding phenomenon is not considered in this linear curvature model. In order to avoid this, the evaluation criteria to measure the registration results are directly integrated into the basic framework of the variational model (3), and we propose a constrained linear curvature image registration model in the following form:(7)minuJαu=12∫ΩTu−R2dΩ+α2∑l=12∫ΩΔul2dΩs.t.Cu>0,where(8)Cu=detJφx=1+u1xu1yu2x1+u2y=1+u1x1+u2y−u1yu2x.Compared with model (5), our new model can ensure that the displacement field is regular both globally. In addition, the new model prevents mesh folding even for very small regularization parameters α. Finally, visually pleasing registration results can be obtained by using our new model with low computing time for smooth registration problems. The numerical solution of the new model (7) is given below.


## 3. Numerical Solution of the New Model

In general, it is difficult to solve the optimization problem (7) by the analytic method. Thus it is necessary to adopt the numerical method and appropriate discretization. In this paper, we choose the discretize-optimize method which aims to take advantage of efficient optimization techniques. In this section, we first discuss briefly the discretization we use and then describe the details of numerical algorithms.


### 3.1. Finite Difference Discretization

Assume that given discrete images havem1×m2 pixels. For simplicity, the image region is further assumed to be Ω=0,1×0,1⊂ℝ2, and then each side of these m1×m2 cell-centered images has width hi=1/mi,i=1,2. Thus the discrete domain can be denoted by(9)Ωh=x∈Ωx=xi,yj⊤=i−0.5h1,j−0.5h2⊤,i=1,2,…,n1;j=1,2,…,n2.


#### 3.1.1. Discretization of Regularizer

The discrete form of the continuous displacement fieldu=u1,u2⊤ can be represented by uh=u1h,u2h, where u1h and u2h are the discrete grid functions defined on the discrete region Ωh. For convenience, let ulhij=ulhxi,yj,i=1,2,…,m1,j=1,2,…,m2, and l=1,2. Since the curvature regularizer is expressed based on the Laplacian operator Δ which can be regarded as the product of gradient operator ∇ and divergence operator ∇, we introduce the symbols ∇h and ∇h⋅ to represent their discrete forms, respectively. The discrete gradient operator ∇h can be defined at each pixel i,j by the following form:(10)∇huhi,j=∇hu1hi,j,∇hu2hi,j⊤,where(11)∇hulhi,j=∂xhulhi,j,∂yhulhi,j⊤,∂xhulhi,j=1h1ulhi+1,j−ulhi,j,ifi<m1;0,ifi=m1;∂yhulhi,j=1h2ulhi,j+1−ulhi,j,ifj<m2;0,ifj=m2.The displacement fieldu satisfies the homogeneous Neumann boundary conditions on the boundary ∂Ω of the image region Ω:(12)∂ul∂ν=0,l=1,2.Through the analysis of continuous setting, we know that the discrete divergence operator is the negative conjugate transposition of the gradient operator, namely,∇⋅=−∇∗. Thus, it can be defined by the following form:(13)∇⋅ωi,j=1h1ω1i,j−ω1i−1,j1h1ω1i,j−1h1ω1i−1,j+1h2ω2i,j−ω2i,j−1,if1<i<n1,1<j<n2;1h2ω2i,j,ifi=j=1;−1h2ω2i,j−1,ifi=n1,j=n2,where ω=ω1,ω2 is a vector. For the convenience of calculation, the grid functions u1h and u2h can be changed into column vectors u1h and u2h according to lexicographical ordering, respectively:(14)u1h=u1,1,1h,…,u1,n1,1h,u1,1,2h,…,u1,n1,2h,…,u1,n1,n2h⊤,u2h=u2,1,1h,…,u2,n1,1h,u2,1,2h,…,u2,n1,2h,…,u2,n1,n2h⊤.We can getu1h∈ℝN, u2h∈ℝN, and Uh=u1h,u2h⊤∈ℝ2N, where N=n1n2. Discrete gradient operator ∇hulhi,j can also be expressed as the product of matrix Ak⊤∈ℝ2×Nk=1,2,…,N and the vector ulhl=1,2 in the following form:(15)Ak⊤ulh=ulhk+1−ulhk;ulhk+n2−ulhk,ifkmodn1≠0andk+n2≤N;0;ulhk+n2−ulhk,ifkmodn1=0andk+n2≤N;ulhk+1−ulhk;0,ifkmodn1≠0andk+n2>N;0;0,ifkmodn1=0andk+n2>N.Let(16)A=A1,A2,…,AN=A1,1,A1,2,…,AN,1,AN,2∈ℝN×2N;Ax=A1,1,A2,1,…,AN,1∈ℝN×N;Ay=A1,2,A2,2,…,AN,2∈ℝN×N.By this notation, we can get(17)∇hu1h=Ax⊤Ay⊤u1h≜Bu1h,∇hu2h=Ax⊤Ay⊤u2h≜Bu2h.LetC=B⊤B, ℂ=C00C, and(18)ℬu=∑l=12Δul2=Δu12+Δu22.Then the discrete form of (18) is as follows:(19)BhUh=−B⊤Bu1h⊤−B⊤Bu1h+−B⊤Bu2h⊤−B⊤Bu2h=−Cu1h⊤−Cu1h+−Cu2h⊤−Cu2h=u1h⊤C⊤Cu1h+u2h⊤C⊤Cu2h=u1h⊤,u2h⊤C⊤C00C⊤Cu1hu2h=u1h⊤,u2h⊤C00CTC00Cu1hu2h=Uh⊤ℂ⊤ℂUh.According to the midpoint quadrature formula, the linear curvature regularizerSLCu=1/2∫ΩℬudΩ has the following discrete form:(20)ShUh=12hdUh⊤ℂ⊤ℂUh,where hd=h1h2.


#### 3.1.2. Discretization of TemplateT and Reference R

For a given discrete image, if we want to know the gray value at any spatial location other than the grid point, then image interpolation is needed. In order to take full advantage of the fast and effective optimization method, a smooth cubicB-spline is used for interpolation. Next, T and ℛ are used to represent the continuous smooth approximation of template image T and reference image R, respectively. Let(21)xc=x1,1,…,xn1,1,x1,2,…,xn1,2,…,x1,n2,…,xn1,n2⊤,yc=y1,1,…,yn1,1,y1,2,…,yn1,2,…,y1,n2,…,yn1,n2⊤,and Xch=xc;yc.Thus the discrete reference image and transformed template image can be represented by the following form, respectively:(22)R⟶=ℛXch,(23)T⟶Uh=TXch+Uh,and further we can get the Jacobian of T⟶:(24)T⟶Uh=∂T⟶∂UhUh=∂T∂UchUch,where Uch=Xch+Uh and the Jacobian of T⟶ is a block matrix with diagonal blocks.


#### 3.1.3. Discretization of Distance MeasureD

Although it is in a continuous setting, it is not possible to compute integrals analytically. Thus it is necessary to use numerical integration. In discrete simulation, the midpoint quadrature formula can be used to approximate the integral. According to (22) and (23), the discrete form of distance measurement D2 can be written directly as follows:(25)DhUh=12h1h2T⟶Uh−R⟶⊤T⟶Uh−R⟶.In addition, the derivative of the discrete functionalDhUh on Uh can also be calculated and has the following form:(26)dDhUh=h1h2T⟶Uh⊤T⟶Uh−R⟶.Furthermore, we can calculate the second derivatived2DhUh of the distance measurement DhUh:(27)d2DhUh=h1h2T⟶Uh⊤T⟶Uh+h1h2∑i=1NdiUh∇2diUh,where dUh=T⟶Uh−R⟶∈ℝN. On one hand, it is consuming and numerically unstable to compute higher-order derivatives (27) in registering two images for practical applications. On the other hand, the difference between T⟶Uh and R⟶ will become smaller if the template image is well registered. To have an efficient and stable numerical algorithm as proposed in work [5], d2DhUh can be approximated by the following form:(28)d2DhUh=h1h2T⟶Uh⊤T⟶Uh.


#### 3.1.4. Discretization of Inequality Constraint FunctionalCu

In model (7), the inequality constraint functional Cu is defined by(29)Cu=1+u1x1+u2y−u1yu2x.According to the previous analysis, the discrete form of the partial derivative of the continuous displacement field elementul can be expressed as follows:(30)ulhx=Ax⊤ulh≜ml;ulhy=Ay⊤ulh≜wl,l=1,2.Obviously,ml∈ℝN, and wl∈ℝN, where N=n1×n2. Let(31)e=1,1,…,1⊤∈ℝN,c=e+m1⊛e+w2−w1⊛m2,where symbol ⊛ denotes the multiplication of the corresponding elements of two vectors and c∈ℝN. For the convenience of calculation, let ci denote the i-th element of c, i=1,2,…,N. Therefore, the continuous inequality constraint function Cu has the following discrete form:(32)ChUh=c1,c2,…,cN⊤.Since the first-order variation of the continuous inequality constraint functionCu on continuous displacement field u is as follows:(33)dCu=u2xy−u2yx,u1yx−u1xy⊤.Thus we can get the discrete form of the first-order variationdCu:(34)dChUh=0Ay⊤Ax⊤−Ax⊤Ay⊤Ax⊤Ay⊤−Ay⊤Ax⊤0u1hu2h≜AUh.Obviously,dChUh∈ℝ2N, 0∈ℝN×N, and A∈ℝ2N×2N.


### 3.2. Solving the Discrete Optimization Problem

According to the above analysis, inequality constrained functional (7) has the following discrete form:(35)minUhJαUh=DhUh+αShUhs.t.ChUh>0.Below we use the multiplier method to numerically solve the inequality constrained optimization problem (35). The basic idea of this method is to transform the original problem into a series of unconstrained optimization problems to solve and simultaneously estimate the Lagrangian multiplier. For more details on multiplier scheme, see [37]. Before solving (35), let us briefly review the multiplier method of inequality constrained optimization.


#### 3.2.1. Multiplier Method for Inequality Constrained Problems

Consider the following inequality constrained optimization problem:(36)minfx,s.t.gix≥0,i=1,2,…,m.Letyi≥0, and the above inequality constraint can be transformed into the following equivalent equality constraint problem:(37)minfx,s.t.gix−yi2=0,i=1,2,…,m.In this case, the augmented Lagrange function can be expressed as(38)φ˜x,y,λ,σ=fx−∑i=1mλigix−yi2+σ2∑i=1mgix−yi22.In order to eliminate the auxiliary variablesy, the minimization of φ˜ with respect to variable y can be considered. According to the first-order necessary condition, let(39)∇yφ˜x,y,λ,σ=0.We can get(40)2yiλi−2σyigix−yi2=0,i=1,2,…,m.Namely,(41)yiσyi2−σgix−λi=0,i=1,2,…,m.Therefore, whenσgix−λi>0,(42)yi2=1σσgix−λiσgix−λi>00σgix−λi≤0,i=1,2,…,m;that is to say,(43)gix−yi2=λiσσgix−λi>0gixσgix−λi≤0,i=1,2,…,m.Thus whenσgix−λi≤0, we have(44)−λigix−yi2+σ2gix−yi22=−λigix+σ2gix2=12σσgix−λi2−λi2.And whenσgix−λi>0, we can obtain(45)−λigix−yi2+σ2gix−yi22=−1σλi2+12σλi2=−12σλi2.According to the above two cases,(46)−λigix−yi2+σ2gix−yi22=12σmin0,σgix−λi2−λi2.Substituting it into formula (38), we can get the corresponding augmented Lagrange function of (36):(47)ϕ˜x,λ,σ=minyψ˜x,y,λ,σ=fx+12σ∑i=1mmin0,σgix−λi2−λi2.Since the multiplier vector needs to be updated to solve the inequality constrained optimization problems (36) by using the multiplier method, next we derive the multiplier iterative formula. Firstly, fix the penalty parameter σ to some value σk>0 at its k-th iteration, and fix λ at the current estimate λk. Secondly, perform minimization with respect to x. Using xk to denote the approximate minimizer of φ˜x,y,λ,σ, then we can get by the optimality conditions for unconstrained minimization that(48)0≈∇xφ˜xk,yk,λk,σk=∇fxk−∑i=1mλik−σkgixk−yik2∇gixk−yik2.Letx∗,y∗,λ∗ satisfy the KKT conditions for (37), then we have(49)∇fx∗−∑i=1mλi∗∇gix∗−yi∗2=0.By comparing (48) with (49), we can deduce that(50)λi∗≈λik−σkgixk−yik2,i=1,2,…,m.According to (50), to improve the current estimate λk of the Lagrange multiplier vectors, the multiplier iteration formula can be given by the following form:(51)λik+1=λik−σkgixk−yik2.Then, taking (43) into the multiplier iteration formula (51), we have(52)λik+1=0,σkgixk−λik>0,λik−σkgixk,σkgixk−λik≤0.Furthermore, it can be written as(53)λik+1=max0,λik−σkgixk≥0,i=1,2,…,m.Similarly, take (43) into the termination criterion(54)∑i=1mgixk−yik221/2≤ε.We can get(55)∑i=1mmingixk,λikσk21/2≤ε.


#### 3.2.2. Multiplier Method for Solving Model

Next, we use the multiplier method to solve the model (35). Firstly, we construct the augmented Lagrange function for solving model (35):(56)ψUh,λ,σ=JαUh+12σ∑i=1Nmin0,σCihUh−λi2−λi2.The corresponding multiplier iteration formula has the following form:(57)λik+1=max0,λik−σkCihUhk.And the corresponding stopping criterion is(58)βk=∑i=1NminCihUhk,λikσk21/2≤ε.Although the augmented Lagrangian function (57) of the model (35) contains the min function, it is still continuously differentiable; for details, see [37, 44]. The detailed steps of the multiplier method for solving the model (35) can be summarized by Algorithm 1.Algorithm 1: Multiplier scheme.
Step 1: input the initial value:Uh0∈ℝ2N, the objective function JαUh and its gradient dJαUh, inequality constrained vector ChUh, and the transpose of its Jacobian matrix dChUh; let maxk:=10, σ1:=1, ε:=10−5, ϑ:=0.3, η:=0.2, k:=0, λ1:=10−3,10−3,…,10−3⊤∈ℝN, and βk:=10.Step 2: solving the subproblem. WithUhk−1 as the initial point, solve the minimum value Uhk of the unconstrained subproblem (51) by using the Gauss–Newton scheme with Armijo line search.Step 3: check the termination condition. Ifβk≤ε or k>maxk, where βk is defined by (57), the iteration is stopped, and Uhk is output as the approximate minimum of the original problem; otherwise, go to Step 4.Step 4: update penalty parameters. Ifβk≥ϑβk−1, let σk+1:=ησk; otherwise, set σk+1:=σk.Step 5: update multiplier vector. Calculateλik+1=max0,λik−σkCihUhk,i=1,2,…,N.Step 6: setk:=k+1, and go to Step 1.In Algorithm1, the Gauss–Newton method is used to solve the unconstrained subproblem (56). And its main idea is to use a quadratic function ψ^ instead of ψ near the iteration value Uhk of the previous step by the Taylor expansion given below:(59)ψUhk+δUh≈ψ^Uhk+δUh=ψUhk+dψUhkδUh+12δUh⊤HδUh,where dψUhk is the Jacobian matrix of ψ at Uhk and H is the approximation of its Hessian. Due to d2DhUhk, ℂ⊤ℂ and MUhk⊤MUhk are both positive semidefine, and it is easy to prove that H is also positive semidefinite. Thus, we know that ψ^ is convex. In this way, the nonconvex problem can be transformed into a convex problem to be solved. For further detailed description, see [37]. The detailed steps are described below.Given the initial valueUhk, the Jacobian matrix dψUhk and the Hessian matrix H are calculated by the following forms:(60)dψUhk=dDhUhk+αhdℂ⊤ℂUhk+MUhk⊛σChUhk−λ,(61)H=d2DhUhk+αhdℂ⊤ℂ+MUhkMUhk⊤,and in each outer iteration step, respectively, where MUhk=dChUh∈ℝ2N is defined by (34). In formula (59), the disturbance value δUh can be obtained by finding the stability point of the quadratic function ψ^, namely,(62)HδUh=−dψUhk.Usually,H is the positive definite. Thus the equation (62) can be solved by the preconjugate gradient method. To ensure that the objective function (59) is descending, the standard Armijo line search method can be used. The specific steps for the Armijos line search can be summarized by using Algorithm 2. And the Gauss–Newton method mentioned above is described by using Algorithm 3. In order to provide a good initial value, the Gauss–Newton method with Armijo line search and the multilevel method are combined to solve model (56), which can reduce the risk of getting trapped at an unwanted minimizer and save computing time. Firstly, we use the initial value Uh0 and the Gauss–Newton method with the Armijos line search to solve (56) on the coarsest level. Secondly, the solution on the coarsest level is interpolated to the next finer level; next it is used as the initial value, and the same method is used to solve the model (56) on the finer level, where bilinear interpolation operator IHh is used. Finally, this process is repeated until the loop terminates. We summarize the multilevel method using Algorithm 4.Algorithm 2: Armijo line search method.
Step 1: input initial value:Uhk∈ℝ2N; δUh∈ℝN, ψUhk, and dψUhk are calculated by (56) and (60), respectively. Set t:=1, maxIter:=10, η:=10−4, and k:=0.Step 2: setUthk:=Uhk+tδUh, and compute ψUthk.Step 3: check the termination condition. IfψUthk<ψUhk+tηdψUhk⊤δUh, then the iteration is stopped; otherwise, go to Step 4.Step 4: Sett:=t/2; k:=k+1; and go to Step 2.Algorithm 3: The Gauss–Newton method.
Step 1: setk:=0, maxIter:=10, and ε=2.2e−13; input initial value: Uhk∈ℝ2N.Step 2: computeψUhk, dψUhk, and H by using (56), (60), and (61), respectively.Step 3: check the termination condition. IfdψUhk≤ε or k≥maxIter, stop, and Uh∗:=Uhk is the output.Step 4: Solve the Quasi–Newton equation (62) by Preconjugate Gradient method. If the equation (61) has a solution δUh which meets dψUhkδUh<0, then go to Step 5; otherwise, set δUh=−dψUhk, and go to Step 5.Step 5: the step factort is solved by the Armijos line search technique.Step 6: setUthk:=Uhk+tδUh, k:=k+1; then go to Step 2.Algorithm 4: Multilevel method.
Step 1: given initial values:Maxlevel:=ceillog2minm1,m2 and Minlevel:=3; given multilevel representation of the reference image R and the template image T.Step 2: setUh0=0 on the coarsest level l=Minlevel, then solve (51) by using Algorithm 1; otherwise, go to Step 3.Step 3: the initial valueUh0 on the finer level is obtained by interpolation operator IHh.Step 4: setl:=l+1; then go to Step 2.


## 3.1. Finite Difference Discretization

Assume that given discrete images havem1×m2 pixels. For simplicity, the image region is further assumed to be Ω=0,1×0,1⊂ℝ2, and then each side of these m1×m2 cell-centered images has width hi=1/mi,i=1,2. Thus the discrete domain can be denoted by(9)Ωh=x∈Ωx=xi,yj⊤=i−0.5h1,j−0.5h2⊤,i=1,2,…,n1;j=1,2,…,n2.


### 3.1.1. Discretization of Regularizer

The discrete form of the continuous displacement fieldu=u1,u2⊤ can be represented by uh=u1h,u2h, where u1h and u2h are the discrete grid functions defined on the discrete region Ωh. For convenience, let ulhij=ulhxi,yj,i=1,2,…,m1,j=1,2,…,m2, and l=1,2. Since the curvature regularizer is expressed based on the Laplacian operator Δ which can be regarded as the product of gradient operator ∇ and divergence operator ∇, we introduce the symbols ∇h and ∇h⋅ to represent their discrete forms, respectively. The discrete gradient operator ∇h can be defined at each pixel i,j by the following form:(10)∇huhi,j=∇hu1hi,j,∇hu2hi,j⊤,where(11)∇hulhi,j=∂xhulhi,j,∂yhulhi,j⊤,∂xhulhi,j=1h1ulhi+1,j−ulhi,j,ifi<m1;0,ifi=m1;∂yhulhi,j=1h2ulhi,j+1−ulhi,j,ifj<m2;0,ifj=m2.The displacement fieldu satisfies the homogeneous Neumann boundary conditions on the boundary ∂Ω of the image region Ω:(12)∂ul∂ν=0,l=1,2.Through the analysis of continuous setting, we know that the discrete divergence operator is the negative conjugate transposition of the gradient operator, namely,∇⋅=−∇∗. Thus, it can be defined by the following form:(13)∇⋅ωi,j=1h1ω1i,j−ω1i−1,j1h1ω1i,j−1h1ω1i−1,j+1h2ω2i,j−ω2i,j−1,if1<i<n1,1<j<n2;1h2ω2i,j,ifi=j=1;−1h2ω2i,j−1,ifi=n1,j=n2,where ω=ω1,ω2 is a vector. For the convenience of calculation, the grid functions u1h and u2h can be changed into column vectors u1h and u2h according to lexicographical ordering, respectively:(14)u1h=u1,1,1h,…,u1,n1,1h,u1,1,2h,…,u1,n1,2h,…,u1,n1,n2h⊤,u2h=u2,1,1h,…,u2,n1,1h,u2,1,2h,…,u2,n1,2h,…,u2,n1,n2h⊤.We can getu1h∈ℝN, u2h∈ℝN, and Uh=u1h,u2h⊤∈ℝ2N, where N=n1n2. Discrete gradient operator ∇hulhi,j can also be expressed as the product of matrix Ak⊤∈ℝ2×Nk=1,2,…,N and the vector ulhl=1,2 in the following form:(15)Ak⊤ulh=ulhk+1−ulhk;ulhk+n2−ulhk,ifkmodn1≠0andk+n2≤N;0;ulhk+n2−ulhk,ifkmodn1=0andk+n2≤N;ulhk+1−ulhk;0,ifkmodn1≠0andk+n2>N;0;0,ifkmodn1=0andk+n2>N.Let(16)A=A1,A2,…,AN=A1,1,A1,2,…,AN,1,AN,2∈ℝN×2N;Ax=A1,1,A2,1,…,AN,1∈ℝN×N;Ay=A1,2,A2,2,…,AN,2∈ℝN×N.By this notation, we can get(17)∇hu1h=Ax⊤Ay⊤u1h≜Bu1h,∇hu2h=Ax⊤Ay⊤u2h≜Bu2h.LetC=B⊤B, ℂ=C00C, and(18)ℬu=∑l=12Δul2=Δu12+Δu22.Then the discrete form of (18) is as follows:(19)BhUh=−B⊤Bu1h⊤−B⊤Bu1h+−B⊤Bu2h⊤−B⊤Bu2h=−Cu1h⊤−Cu1h+−Cu2h⊤−Cu2h=u1h⊤C⊤Cu1h+u2h⊤C⊤Cu2h=u1h⊤,u2h⊤C⊤C00C⊤Cu1hu2h=u1h⊤,u2h⊤C00CTC00Cu1hu2h=Uh⊤ℂ⊤ℂUh.According to the midpoint quadrature formula, the linear curvature regularizerSLCu=1/2∫ΩℬudΩ has the following discrete form:(20)ShUh=12hdUh⊤ℂ⊤ℂUh,where hd=h1h2.


### 3.1.2. Discretization of TemplateT and Reference R

For a given discrete image, if we want to know the gray value at any spatial location other than the grid point, then image interpolation is needed. In order to take full advantage of the fast and effective optimization method, a smooth cubicB-spline is used for interpolation. Next, T and ℛ are used to represent the continuous smooth approximation of template image T and reference image R, respectively. Let(21)xc=x1,1,…,xn1,1,x1,2,…,xn1,2,…,x1,n2,…,xn1,n2⊤,yc=y1,1,…,yn1,1,y1,2,…,yn1,2,…,y1,n2,…,yn1,n2⊤,and Xch=xc;yc.Thus the discrete reference image and transformed template image can be represented by the following form, respectively:(22)R⟶=ℛXch,(23)T⟶Uh=TXch+Uh,and further we can get the Jacobian of T⟶:(24)T⟶Uh=∂T⟶∂UhUh=∂T∂UchUch,where Uch=Xch+Uh and the Jacobian of T⟶ is a block matrix with diagonal blocks.


### 3.1.3. Discretization of Distance MeasureD

Although it is in a continuous setting, it is not possible to compute integrals analytically. Thus it is necessary to use numerical integration. In discrete simulation, the midpoint quadrature formula can be used to approximate the integral. According to (22) and (23), the discrete form of distance measurement D2 can be written directly as follows:(25)DhUh=12h1h2T⟶Uh−R⟶⊤T⟶Uh−R⟶.In addition, the derivative of the discrete functionalDhUh on Uh can also be calculated and has the following form:(26)dDhUh=h1h2T⟶Uh⊤T⟶Uh−R⟶.Furthermore, we can calculate the second derivatived2DhUh of the distance measurement DhUh:(27)d2DhUh=h1h2T⟶Uh⊤T⟶Uh+h1h2∑i=1NdiUh∇2diUh,where dUh=T⟶Uh−R⟶∈ℝN. On one hand, it is consuming and numerically unstable to compute higher-order derivatives (27) in registering two images for practical applications. On the other hand, the difference between T⟶Uh and R⟶ will become smaller if the template image is well registered. To have an efficient and stable numerical algorithm as proposed in work [5], d2DhUh can be approximated by the following form:(28)d2DhUh=h1h2T⟶Uh⊤T⟶Uh.


### 3.1.4. Discretization of Inequality Constraint FunctionalCu

In model (7), the inequality constraint functional Cu is defined by(29)Cu=1+u1x1+u2y−u1yu2x.According to the previous analysis, the discrete form of the partial derivative of the continuous displacement field elementul can be expressed as follows:(30)ulhx=Ax⊤ulh≜ml;ulhy=Ay⊤ulh≜wl,l=1,2.Obviously,ml∈ℝN, and wl∈ℝN, where N=n1×n2. Let(31)e=1,1,…,1⊤∈ℝN,c=e+m1⊛e+w2−w1⊛m2,where symbol ⊛ denotes the multiplication of the corresponding elements of two vectors and c∈ℝN. For the convenience of calculation, let ci denote the i-th element of c, i=1,2,…,N. Therefore, the continuous inequality constraint function Cu has the following discrete form:(32)ChUh=c1,c2,…,cN⊤.Since the first-order variation of the continuous inequality constraint functionCu on continuous displacement field u is as follows:(33)dCu=u2xy−u2yx,u1yx−u1xy⊤.Thus we can get the discrete form of the first-order variationdCu:(34)dChUh=0Ay⊤Ax⊤−Ax⊤Ay⊤Ax⊤Ay⊤−Ay⊤Ax⊤0u1hu2h≜AUh.Obviously,dChUh∈ℝ2N, 0∈ℝN×N, and A∈ℝ2N×2N.


## 3.1.1. Discretization of Regularizer

The discrete form of the continuous displacement fieldu=u1,u2⊤ can be represented by uh=u1h,u2h, where u1h and u2h are the discrete grid functions defined on the discrete region Ωh. For convenience, let ulhij=ulhxi,yj,i=1,2,…,m1,j=1,2,…,m2, and l=1,2. Since the curvature regularizer is expressed based on the Laplacian operator Δ which can be regarded as the product of gradient operator ∇ and divergence operator ∇, we introduce the symbols ∇h and ∇h⋅ to represent their discrete forms, respectively. The discrete gradient operator ∇h can be defined at each pixel i,j by the following form:(10)∇huhi,j=∇hu1hi,j,∇hu2hi,j⊤,where(11)∇hulhi,j=∂xhulhi,j,∂yhulhi,j⊤,∂xhulhi,j=1h1ulhi+1,j−ulhi,j,ifi<m1;0,ifi=m1;∂yhulhi,j=1h2ulhi,j+1−ulhi,j,ifj<m2;0,ifj=m2.The displacement fieldu satisfies the homogeneous Neumann boundary conditions on the boundary ∂Ω of the image region Ω:(12)∂ul∂ν=0,l=1,2.Through the analysis of continuous setting, we know that the discrete divergence operator is the negative conjugate transposition of the gradient operator, namely,∇⋅=−∇∗. Thus, it can be defined by the following form:(13)∇⋅ωi,j=1h1ω1i,j−ω1i−1,j1h1ω1i,j−1h1ω1i−1,j+1h2ω2i,j−ω2i,j−1,if1<i<n1,1<j<n2;1h2ω2i,j,ifi=j=1;−1h2ω2i,j−1,ifi=n1,j=n2,where ω=ω1,ω2 is a vector. For the convenience of calculation, the grid functions u1h and u2h can be changed into column vectors u1h and u2h according to lexicographical ordering, respectively:(14)u1h=u1,1,1h,…,u1,n1,1h,u1,1,2h,…,u1,n1,2h,…,u1,n1,n2h⊤,u2h=u2,1,1h,…,u2,n1,1h,u2,1,2h,…,u2,n1,2h,…,u2,n1,n2h⊤.We can getu1h∈ℝN, u2h∈ℝN, and Uh=u1h,u2h⊤∈ℝ2N, where N=n1n2. Discrete gradient operator ∇hulhi,j can also be expressed as the product of matrix Ak⊤∈ℝ2×Nk=1,2,…,N and the vector ulhl=1,2 in the following form:(15)Ak⊤ulh=ulhk+1−ulhk;ulhk+n2−ulhk,ifkmodn1≠0andk+n2≤N;0;ulhk+n2−ulhk,ifkmodn1=0andk+n2≤N;ulhk+1−ulhk;0,ifkmodn1≠0andk+n2>N;0;0,ifkmodn1=0andk+n2>N.Let(16)A=A1,A2,…,AN=A1,1,A1,2,…,AN,1,AN,2∈ℝN×2N;Ax=A1,1,A2,1,…,AN,1∈ℝN×N;Ay=A1,2,A2,2,…,AN,2∈ℝN×N.By this notation, we can get(17)∇hu1h=Ax⊤Ay⊤u1h≜Bu1h,∇hu2h=Ax⊤Ay⊤u2h≜Bu2h.LetC=B⊤B, ℂ=C00C, and(18)ℬu=∑l=12Δul2=Δu12+Δu22.Then the discrete form of (18) is as follows:(19)BhUh=−B⊤Bu1h⊤−B⊤Bu1h+−B⊤Bu2h⊤−B⊤Bu2h=−Cu1h⊤−Cu1h+−Cu2h⊤−Cu2h=u1h⊤C⊤Cu1h+u2h⊤C⊤Cu2h=u1h⊤,u2h⊤C⊤C00C⊤Cu1hu2h=u1h⊤,u2h⊤C00CTC00Cu1hu2h=Uh⊤ℂ⊤ℂUh.According to the midpoint quadrature formula, the linear curvature regularizerSLCu=1/2∫ΩℬudΩ has the following discrete form:(20)ShUh=12hdUh⊤ℂ⊤ℂUh,where hd=h1h2.


## 3.1.2. Discretization of TemplateT and Reference R

For a given discrete image, if we want to know the gray value at any spatial location other than the grid point, then image interpolation is needed. In order to take full advantage of the fast and effective optimization method, a smooth cubicB-spline is used for interpolation. Next, T and ℛ are used to represent the continuous smooth approximation of template image T and reference image R, respectively. Let(21)xc=x1,1,…,xn1,1,x1,2,…,xn1,2,…,x1,n2,…,xn1,n2⊤,yc=y1,1,…,yn1,1,y1,2,…,yn1,2,…,y1,n2,…,yn1,n2⊤,and Xch=xc;yc.Thus the discrete reference image and transformed template image can be represented by the following form, respectively:(22)R⟶=ℛXch,(23)T⟶Uh=TXch+Uh,and further we can get the Jacobian of T⟶:(24)T⟶Uh=∂T⟶∂UhUh=∂T∂UchUch,where Uch=Xch+Uh and the Jacobian of T⟶ is a block matrix with diagonal blocks.


## 3.1.3. Discretization of Distance MeasureD

Although it is in a continuous setting, it is not possible to compute integrals analytically. Thus it is necessary to use numerical integration. In discrete simulation, the midpoint quadrature formula can be used to approximate the integral. According to (22) and (23), the discrete form of distance measurement D2 can be written directly as follows:(25)DhUh=12h1h2T⟶Uh−R⟶⊤T⟶Uh−R⟶.In addition, the derivative of the discrete functionalDhUh on Uh can also be calculated and has the following form:(26)dDhUh=h1h2T⟶Uh⊤T⟶Uh−R⟶.Furthermore, we can calculate the second derivatived2DhUh of the distance measurement DhUh:(27)d2DhUh=h1h2T⟶Uh⊤T⟶Uh+h1h2∑i=1NdiUh∇2diUh,where dUh=T⟶Uh−R⟶∈ℝN. On one hand, it is consuming and numerically unstable to compute higher-order derivatives (27) in registering two images for practical applications. On the other hand, the difference between T⟶Uh and R⟶ will become smaller if the template image is well registered. To have an efficient and stable numerical algorithm as proposed in work [5], d2DhUh can be approximated by the following form:(28)d2DhUh=h1h2T⟶Uh⊤T⟶Uh.


## 3.1.4. Discretization of Inequality Constraint FunctionalCu

In model (7), the inequality constraint functional Cu is defined by(29)Cu=1+u1x1+u2y−u1yu2x.According to the previous analysis, the discrete form of the partial derivative of the continuous displacement field elementul can be expressed as follows:(30)ulhx=Ax⊤ulh≜ml;ulhy=Ay⊤ulh≜wl,l=1,2.Obviously,ml∈ℝN, and wl∈ℝN, where N=n1×n2. Let(31)e=1,1,…,1⊤∈ℝN,c=e+m1⊛e+w2−w1⊛m2,where symbol ⊛ denotes the multiplication of the corresponding elements of two vectors and c∈ℝN. For the convenience of calculation, let ci denote the i-th element of c, i=1,2,…,N. Therefore, the continuous inequality constraint function Cu has the following discrete form:(32)ChUh=c1,c2,…,cN⊤.Since the first-order variation of the continuous inequality constraint functionCu on continuous displacement field u is as follows:(33)dCu=u2xy−u2yx,u1yx−u1xy⊤.Thus we can get the discrete form of the first-order variationdCu:(34)dChUh=0Ay⊤Ax⊤−Ax⊤Ay⊤Ax⊤Ay⊤−Ay⊤Ax⊤0u1hu2h≜AUh.Obviously,dChUh∈ℝ2N, 0∈ℝN×N, and A∈ℝ2N×2N.


## 3.2. Solving the Discrete Optimization Problem

According to the above analysis, inequality constrained functional (7) has the following discrete form:(35)minUhJαUh=DhUh+αShUhs.t.ChUh>0.Below we use the multiplier method to numerically solve the inequality constrained optimization problem (35). The basic idea of this method is to transform the original problem into a series of unconstrained optimization problems to solve and simultaneously estimate the Lagrangian multiplier. For more details on multiplier scheme, see [37]. Before solving (35), let us briefly review the multiplier method of inequality constrained optimization.


### 3.2.1. Multiplier Method for Inequality Constrained Problems

Consider the following inequality constrained optimization problem:(36)minfx,s.t.gix≥0,i=1,2,…,m.Letyi≥0, and the above inequality constraint can be transformed into the following equivalent equality constraint problem:(37)minfx,s.t.gix−yi2=0,i=1,2,…,m.In this case, the augmented Lagrange function can be expressed as(38)φ˜x,y,λ,σ=fx−∑i=1mλigix−yi2+σ2∑i=1mgix−yi22.In order to eliminate the auxiliary variablesy, the minimization of φ˜ with respect to variable y can be considered. According to the first-order necessary condition, let(39)∇yφ˜x,y,λ,σ=0.We can get(40)2yiλi−2σyigix−yi2=0,i=1,2,…,m.Namely,(41)yiσyi2−σgix−λi=0,i=1,2,…,m.Therefore, whenσgix−λi>0,(42)yi2=1σσgix−λiσgix−λi>00σgix−λi≤0,i=1,2,…,m;that is to say,(43)gix−yi2=λiσσgix−λi>0gixσgix−λi≤0,i=1,2,…,m.Thus whenσgix−λi≤0, we have(44)−λigix−yi2+σ2gix−yi22=−λigix+σ2gix2=12σσgix−λi2−λi2.And whenσgix−λi>0, we can obtain(45)−λigix−yi2+σ2gix−yi22=−1σλi2+12σλi2=−12σλi2.According to the above two cases,(46)−λigix−yi2+σ2gix−yi22=12σmin0,σgix−λi2−λi2.Substituting it into formula (38), we can get the corresponding augmented Lagrange function of (36):(47)ϕ˜x,λ,σ=minyψ˜x,y,λ,σ=fx+12σ∑i=1mmin0,σgix−λi2−λi2.Since the multiplier vector needs to be updated to solve the inequality constrained optimization problems (36) by using the multiplier method, next we derive the multiplier iterative formula. Firstly, fix the penalty parameter σ to some value σk>0 at its k-th iteration, and fix λ at the current estimate λk. Secondly, perform minimization with respect to x. Using xk to denote the approximate minimizer of φ˜x,y,λ,σ, then we can get by the optimality conditions for unconstrained minimization that(48)0≈∇xφ˜xk,yk,λk,σk=∇fxk−∑i=1mλik−σkgixk−yik2∇gixk−yik2.Letx∗,y∗,λ∗ satisfy the KKT conditions for (37), then we have(49)∇fx∗−∑i=1mλi∗∇gix∗−yi∗2=0.By comparing (48) with (49), we can deduce that(50)λi∗≈λik−σkgixk−yik2,i=1,2,…,m.According to (50), to improve the current estimate λk of the Lagrange multiplier vectors, the multiplier iteration formula can be given by the following form:(51)λik+1=λik−σkgixk−yik2.Then, taking (43) into the multiplier iteration formula (51), we have(52)λik+1=0,σkgixk−λik>0,λik−σkgixk,σkgixk−λik≤0.Furthermore, it can be written as(53)λik+1=max0,λik−σkgixk≥0,i=1,2,…,m.Similarly, take (43) into the termination criterion(54)∑i=1mgixk−yik221/2≤ε.We can get(55)∑i=1mmingixk,λikσk21/2≤ε.


### 3.2.2. Multiplier Method for Solving Model

Next, we use the multiplier method to solve the model (35). Firstly, we construct the augmented Lagrange function for solving model (35):(56)ψUh,λ,σ=JαUh+12σ∑i=1Nmin0,σCihUh−λi2−λi2.The corresponding multiplier iteration formula has the following form:(57)λik+1=max0,λik−σkCihUhk.And the corresponding stopping criterion is(58)βk=∑i=1NminCihUhk,λikσk21/2≤ε.Although the augmented Lagrangian function (57) of the model (35) contains the min function, it is still continuously differentiable; for details, see [37, 44]. The detailed steps of the multiplier method for solving the model (35) can be summarized by Algorithm 1.Algorithm 1: Multiplier scheme.
Step 1: input the initial value:Uh0∈ℝ2N, the objective function JαUh and its gradient dJαUh, inequality constrained vector ChUh, and the transpose of its Jacobian matrix dChUh; let maxk:=10, σ1:=1, ε:=10−5, ϑ:=0.3, η:=0.2, k:=0, λ1:=10−3,10−3,…,10−3⊤∈ℝN, and βk:=10.Step 2: solving the subproblem. WithUhk−1 as the initial point, solve the minimum value Uhk of the unconstrained subproblem (51) by using the Gauss–Newton scheme with Armijo line search.Step 3: check the termination condition. Ifβk≤ε or k>maxk, where βk is defined by (57), the iteration is stopped, and Uhk is output as the approximate minimum of the original problem; otherwise, go to Step 4.Step 4: update penalty parameters. Ifβk≥ϑβk−1, let σk+1:=ησk; otherwise, set σk+1:=σk.Step 5: update multiplier vector. Calculateλik+1=max0,λik−σkCihUhk,i=1,2,…,N.Step 6: setk:=k+1, and go to Step 1.In Algorithm1, the Gauss–Newton method is used to solve the unconstrained subproblem (56). And its main idea is to use a quadratic function ψ^ instead of ψ near the iteration value Uhk of the previous step by the Taylor expansion given below:(59)ψUhk+δUh≈ψ^Uhk+δUh=ψUhk+dψUhkδUh+12δUh⊤HδUh,where dψUhk is the Jacobian matrix of ψ at Uhk and H is the approximation of its Hessian. Due to d2DhUhk, ℂ⊤ℂ and MUhk⊤MUhk are both positive semidefine, and it is easy to prove that H is also positive semidefinite. Thus, we know that ψ^ is convex. In this way, the nonconvex problem can be transformed into a convex problem to be solved. For further detailed description, see [37]. The detailed steps are described below.Given the initial valueUhk, the Jacobian matrix dψUhk and the Hessian matrix H are calculated by the following forms:(60)dψUhk=dDhUhk+αhdℂ⊤ℂUhk+MUhk⊛σChUhk−λ,(61)H=d2DhUhk+αhdℂ⊤ℂ+MUhkMUhk⊤,and in each outer iteration step, respectively, where MUhk=dChUh∈ℝ2N is defined by (34). In formula (59), the disturbance value δUh can be obtained by finding the stability point of the quadratic function ψ^, namely,(62)HδUh=−dψUhk.Usually,H is the positive definite. Thus the equation (62) can be solved by the preconjugate gradient method. To ensure that the objective function (59) is descending, the standard Armijo line search method can be used. The specific steps for the Armijos line search can be summarized by using Algorithm 2. And the Gauss–Newton method mentioned above is described by using Algorithm 3. In order to provide a good initial value, the Gauss–Newton method with Armijo line search and the multilevel method are combined to solve model (56), which can reduce the risk of getting trapped at an unwanted minimizer and save computing time. Firstly, we use the initial value Uh0 and the Gauss–Newton method with the Armijos line search to solve (56) on the coarsest level. Secondly, the solution on the coarsest level is interpolated to the next finer level; next it is used as the initial value, and the same method is used to solve the model (56) on the finer level, where bilinear interpolation operator IHh is used. Finally, this process is repeated until the loop terminates. We summarize the multilevel method using Algorithm 4.Algorithm 2: Armijo line search method.
Step 1: input initial value:Uhk∈ℝ2N; δUh∈ℝN, ψUhk, and dψUhk are calculated by (56) and (60), respectively. Set t:=1, maxIter:=10, η:=10−4, and k:=0.Step 2: setUthk:=Uhk+tδUh, and compute ψUthk.Step 3: check the termination condition. IfψUthk<ψUhk+tηdψUhk⊤δUh, then the iteration is stopped; otherwise, go to Step 4.Step 4: Sett:=t/2; k:=k+1; and go to Step 2.Algorithm 3: The Gauss–Newton method.
Step 1: setk:=0, maxIter:=10, and ε=2.2e−13; input initial value: Uhk∈ℝ2N.Step 2: computeψUhk, dψUhk, and H by using (56), (60), and (61), respectively.Step 3: check the termination condition. IfdψUhk≤ε or k≥maxIter, stop, and Uh∗:=Uhk is the output.Step 4: Solve the Quasi–Newton equation (62) by Preconjugate Gradient method. If the equation (61) has a solution δUh which meets dψUhkδUh<0, then go to Step 5; otherwise, set δUh=−dψUhk, and go to Step 5.Step 5: the step factort is solved by the Armijos line search technique.Step 6: setUthk:=Uhk+tδUh, k:=k+1; then go to Step 2.Algorithm 4: Multilevel method.
Step 1: given initial values:Maxlevel:=ceillog2minm1,m2 and Minlevel:=3; given multilevel representation of the reference image R and the template image T.Step 2: setUh0=0 on the coarsest level l=Minlevel, then solve (51) by using Algorithm 1; otherwise, go to Step 3.Step 3: the initial valueUh0 on the finer level is obtained by interpolation operator IHh.Step 4: setl:=l+1; then go to Step 2.


## 3.2.1. Multiplier Method for Inequality Constrained Problems

Consider the following inequality constrained optimization problem:(36)minfx,s.t.gix≥0,i=1,2,…,m.Letyi≥0, and the above inequality constraint can be transformed into the following equivalent equality constraint problem:(37)minfx,s.t.gix−yi2=0,i=1,2,…,m.In this case, the augmented Lagrange function can be expressed as(38)φ˜x,y,λ,σ=fx−∑i=1mλigix−yi2+σ2∑i=1mgix−yi22.In order to eliminate the auxiliary variablesy, the minimization of φ˜ with respect to variable y can be considered. According to the first-order necessary condition, let(39)∇yφ˜x,y,λ,σ=0.We can get(40)2yiλi−2σyigix−yi2=0,i=1,2,…,m.Namely,(41)yiσyi2−σgix−λi=0,i=1,2,…,m.Therefore, whenσgix−λi>0,(42)yi2=1σσgix−λiσgix−λi>00σgix−λi≤0,i=1,2,…,m;that is to say,(43)gix−yi2=λiσσgix−λi>0gixσgix−λi≤0,i=1,2,…,m.Thus whenσgix−λi≤0, we have(44)−λigix−yi2+σ2gix−yi22=−λigix+σ2gix2=12σσgix−λi2−λi2.And whenσgix−λi>0, we can obtain(45)−λigix−yi2+σ2gix−yi22=−1σλi2+12σλi2=−12σλi2.According to the above two cases,(46)−λigix−yi2+σ2gix−yi22=12σmin0,σgix−λi2−λi2.Substituting it into formula (38), we can get the corresponding augmented Lagrange function of (36):(47)ϕ˜x,λ,σ=minyψ˜x,y,λ,σ=fx+12σ∑i=1mmin0,σgix−λi2−λi2.Since the multiplier vector needs to be updated to solve the inequality constrained optimization problems (36) by using the multiplier method, next we derive the multiplier iterative formula. Firstly, fix the penalty parameter σ to some value σk>0 at its k-th iteration, and fix λ at the current estimate λk. Secondly, perform minimization with respect to x. Using xk to denote the approximate minimizer of φ˜x,y,λ,σ, then we can get by the optimality conditions for unconstrained minimization that(48)0≈∇xφ˜xk,yk,λk,σk=∇fxk−∑i=1mλik−σkgixk−yik2∇gixk−yik2.Letx∗,y∗,λ∗ satisfy the KKT conditions for (37), then we have(49)∇fx∗−∑i=1mλi∗∇gix∗−yi∗2=0.By comparing (48) with (49), we can deduce that(50)λi∗≈λik−σkgixk−yik2,i=1,2,…,m.According to (50), to improve the current estimate λk of the Lagrange multiplier vectors, the multiplier iteration formula can be given by the following form:(51)λik+1=λik−σkgixk−yik2.Then, taking (43) into the multiplier iteration formula (51), we have(52)λik+1=0,σkgixk−λik>0,λik−σkgixk,σkgixk−λik≤0.Furthermore, it can be written as(53)λik+1=max0,λik−σkgixk≥0,i=1,2,…,m.Similarly, take (43) into the termination criterion(54)∑i=1mgixk−yik221/2≤ε.We can get(55)∑i=1mmingixk,λikσk21/2≤ε.


## 3.2.2. Multiplier Method for Solving Model

Next, we use the multiplier method to solve the model (35). Firstly, we construct the augmented Lagrange function for solving model (35):(56)ψUh,λ,σ=JαUh+12σ∑i=1Nmin0,σCihUh−λi2−λi2.The corresponding multiplier iteration formula has the following form:(57)λik+1=max0,λik−σkCihUhk.And the corresponding stopping criterion is(58)βk=∑i=1NminCihUhk,λikσk21/2≤ε.Although the augmented Lagrangian function (57) of the model (35) contains the min function, it is still continuously differentiable; for details, see [37, 44]. The detailed steps of the multiplier method for solving the model (35) can be summarized by Algorithm 1.Algorithm 1: Multiplier scheme.
Step 1: input the initial value:Uh0∈ℝ2N, the objective function JαUh and its gradient dJαUh, inequality constrained vector ChUh, and the transpose of its Jacobian matrix dChUh; let maxk:=10, σ1:=1, ε:=10−5, ϑ:=0.3, η:=0.2, k:=0, λ1:=10−3,10−3,…,10−3⊤∈ℝN, and βk:=10.Step 2: solving the subproblem. WithUhk−1 as the initial point, solve the minimum value Uhk of the unconstrained subproblem (51) by using the Gauss–Newton scheme with Armijo line search.Step 3: check the termination condition. Ifβk≤ε or k>maxk, where βk is defined by (57), the iteration is stopped, and Uhk is output as the approximate minimum of the original problem; otherwise, go to Step 4.Step 4: update penalty parameters. Ifβk≥ϑβk−1, let σk+1:=ησk; otherwise, set σk+1:=σk.Step 5: update multiplier vector. Calculateλik+1=max0,λik−σkCihUhk,i=1,2,…,N.Step 6: setk:=k+1, and go to Step 1.In Algorithm1, the Gauss–Newton method is used to solve the unconstrained subproblem (56). And its main idea is to use a quadratic function ψ^ instead of ψ near the iteration value Uhk of the previous step by the Taylor expansion given below:(59)ψUhk+δUh≈ψ^Uhk+δUh=ψUhk+dψUhkδUh+12δUh⊤HδUh,where dψUhk is the Jacobian matrix of ψ at Uhk and H is the approximation of its Hessian. Due to d2DhUhk, ℂ⊤ℂ and MUhk⊤MUhk are both positive semidefine, and it is easy to prove that H is also positive semidefinite. Thus, we know that ψ^ is convex. In this way, the nonconvex problem can be transformed into a convex problem to be solved. For further detailed description, see [37]. The detailed steps are described below.Given the initial valueUhk, the Jacobian matrix dψUhk and the Hessian matrix H are calculated by the following forms:(60)dψUhk=dDhUhk+αhdℂ⊤ℂUhk+MUhk⊛σChUhk−λ,(61)H=d2DhUhk+αhdℂ⊤ℂ+MUhkMUhk⊤,and in each outer iteration step, respectively, where MUhk=dChUh∈ℝ2N is defined by (34). In formula (59), the disturbance value δUh can be obtained by finding the stability point of the quadratic function ψ^, namely,(62)HδUh=−dψUhk.Usually,H is the positive definite. Thus the equation (62) can be solved by the preconjugate gradient method. To ensure that the objective function (59) is descending, the standard Armijo line search method can be used. The specific steps for the Armijos line search can be summarized by using Algorithm 2. And the Gauss–Newton method mentioned above is described by using Algorithm 3. In order to provide a good initial value, the Gauss–Newton method with Armijo line search and the multilevel method are combined to solve model (56), which can reduce the risk of getting trapped at an unwanted minimizer and save computing time. Firstly, we use the initial value Uh0 and the Gauss–Newton method with the Armijos line search to solve (56) on the coarsest level. Secondly, the solution on the coarsest level is interpolated to the next finer level; next it is used as the initial value, and the same method is used to solve the model (56) on the finer level, where bilinear interpolation operator IHh is used. Finally, this process is repeated until the loop terminates. We summarize the multilevel method using Algorithm 4.Algorithm 2: Armijo line search method.
Step 1: input initial value:Uhk∈ℝ2N; δUh∈ℝN, ψUhk, and dψUhk are calculated by (56) and (60), respectively. Set t:=1, maxIter:=10, η:=10−4, and k:=0.Step 2: setUthk:=Uhk+tδUh, and compute ψUthk.Step 3: check the termination condition. IfψUthk<ψUhk+tηdψUhk⊤δUh, then the iteration is stopped; otherwise, go to Step 4.Step 4: Sett:=t/2; k:=k+1; and go to Step 2.Algorithm 3: The Gauss–Newton method.
Step 1: setk:=0, maxIter:=10, and ε=2.2e−13; input initial value: Uhk∈ℝ2N.Step 2: computeψUhk, dψUhk, and H by using (56), (60), and (61), respectively.Step 3: check the termination condition. IfdψUhk≤ε or k≥maxIter, stop, and Uh∗:=Uhk is the output.Step 4: Solve the Quasi–Newton equation (62) by Preconjugate Gradient method. If the equation (61) has a solution δUh which meets dψUhkδUh<0, then go to Step 5; otherwise, set δUh=−dψUhk, and go to Step 5.Step 5: the step factort is solved by the Armijos line search technique.Step 6: setUthk:=Uhk+tδUh, k:=k+1; then go to Step 2.Algorithm 4: Multilevel method.
Step 1: given initial values:Maxlevel:=ceillog2minm1,m2 and Minlevel:=3; given multilevel representation of the reference image R and the template image T.Step 2: setUh0=0 on the coarsest level l=Minlevel, then solve (51) by using Algorithm 1; otherwise, go to Step 3.Step 3: the initial valueUh0 on the finer level is obtained by interpolation operator IHh.Step 4: setl:=l+1; then go to Step 2.


## 4. Numerical Experiments

In this part, we use three experiments to show that our new model (CLC) has good performance by comparing it with diffeomorphic demons (DDemons) [43], linear curvature model (LC) [28], mean curvature model (MC) [8], hyperelastic regularizer (Hyper) [20], and Zhang–Chen model (ZC) [42]. In order to illustrate the capabilities of our model, we select two pairs of artificial images. In addition, considering the important application of image registration in biomedical images, we also use a pair of medical lung images for experiments.In order to quantify the quality of the registered image, the relative reductionε of the dissimilarity, which is given by [7], is shown as follows:(63)ε=DuDstop×100%.And the minimum valueℱ of the determinant of the Jacobian matrix J of the transformation φ(64)J=1+u1xu1yu2x1+u2y,ℱ=mindetJ,are used, where Du is defined by equation (4), u is the current iteration, and Dstop is the value of Du at u=0.


### 4.1. Test 1: A Pair of Lena Images

In the first experiment, we used a pair of Lena images with a size of256×256. The test images and registered ones using our new model are shown in Figure 1. The transformed template image obtained using the other five models and the image difference after registration are represented by Figure 2. For Example 1 shown in Figures 1(a) and 1(b), the registration results using our new model and linear curvature (LC [28]) model and other four models are recorded in Tables 1 and 2, respectively. For a smaller regularizer parameter α=3.12e−3, from Figures 1 and 2, we can see that our new model and hyperelastic regularizer-, linear curvature-, and mean one-based image registration models can produce visually satisfactory registration results. However, we find that the latter two have mesh folding by Tables 1 and 2. Although the demons-based registration model and ZC model do not produce folding, the registration effects are relatively poor. From Tables 1 and 2, we can further see that for the same image with different sizes, our new model can give very satisfactory registration results without folding.Figure 1
Registration results of testing Lena images using our new model. (a) Reference image, (b) template image, (c) the transformed template image using our new model, (d) image difference before registrationε=100%, (e) the transformation x+ux, and (f) image difference after registration ε=0.36%.
(a)(b)(c)(d)(e)(f)Figure 2
Registration results for a pair of Lena images with a size of256×256 using other five models. The transformed template image Tx+ux from (a) the LC model [28], (b) Hyper model [20], and (c) DDemons model [43]. The image difference after registration from (d) the LC model ε=1.48% [28], (e) Hyper model [20] ε=0.52%, and (f) DDemons model [43] ε=2.12%. The transformed template image Tx+ux from (g) the ZC model [42] and (h) MC model [8]. The image difference after registration from (i) the ZC model [42] ε=3.48% and (j) MC model [8] ε=1.41%.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)Table 1
Quantitative measurements using our new model and LC model [28] for processing Example 1 shown in Figures 1(a) and 1(b).
LayerModelLCCLCα=3.12e−3α=3.12e−3ε%Nℱε%Nℱh=1/2561.48315−1.22280.3600.5791h=1/1281.1692−0.88560.2200.4997h=1/641.0420−0.71840.3900.5946h=1/320.664−0.37300.4000.0168h=1/161.231−0.56291.0400.0486N represents the number of mesh folding of the transformation x+ux. ℱ>0 means that the transformation does not include folding or cracking.Table 2
Quantitative measurements using the DDemons model [43], ZC model [42], Hyper model [20], and MC model [8] for processing Example 1 shown in Figures 1(a) and 1(b).
LayerModelDDemonsZCHyperMCσfluid=1, σdiff=0.756α=3.12e−3α=3.12e−3α=3.12e−3ε%Nℱε%Nℱε%Nℱε%Nℱh=1/2562.1200.04343.4800.06490.5200.30381.418−0.1055h=1/1284.6300.17643.2900.06601.0500.64971.141−0.0273h=1/6411.9900.54326.4300.11914.1200.71163.0400.1613h=1/3223.6900.689910.3300.044015.8400.71932.912−0.2810h=1/1648.1400.865028.4500.379285.5500.82472.1700.1656N represents the number of mesh folding of the transformation x+ux. ℱ>0 means that the transformation does not include folding or cracking.In order to analyze how much our model will be affected when changing the regularizer parameterα, we use Algorithm 4 to test Example 1, and the corresponding quantitative measurement values are recorded in Tables 3–5. From Tables 3 to 5, we find that the registration quality becomes worse and worse with the increase of the regularizer parameter α. However, when alpha is greater than or equal to 3.12×10−3, our new model does not produce folding for Example 1. For this example, our new model is able to produce better registration results, especially when the regularizer parameter α is taken in an appropriate range 3.12×10−3,3×10−1.Table 3
Dependent results of our new model on the regularizer parameterα by testing Example 1 (Figures 1(a) and 1(b)).
α=e−33.124581030ε (%)0.35740.37340.38760.41690.43070.4996ℱ0.57910.59680.61190.63870.64930.6826ℱ>0 means that the transformation does not include folding or cracking.Table 4
Dependent results of our new model on the regularizer parameterα by testing Example 1 (Figures 1(a) and 1(b)).
α=e−380100200300400500ε (%)0.58360.61180.74880.89031.03801.1914ℱ0.67360.67190.66780.66690.66790.6699ℱ>0 means that the transformation does not include folding or cracking.Table 5
Dependent results of our new model on the regularizer parameterα by testing Example 1 (Figures 1(a) and 1(b)).
α=e−31e32e33e34e35e36e3ε (%)2.02853.84245.67387.45499.20010.8363ℱ0.68950.73500.75340.76880.77980.7878ℱ>0 means that the transformation does not include folding or cracking.


### 4.2. Test 2: A Pair of Medical Lung Images

In this part, we select a pair of medical lung images with a size of256×256 for testing. For each model, we choose the optimal parameters. The test images and the corresponding registered ones are shown in Figures 3 and 4, respectively. The registration results on different layers are recorded in Tables 6 and 7, respectively. According to Figures 3 and 4, we can observe that the other five models can produce visually relatively satisfactory registration effects except the DDemons model. From the image difference after registration and the results recorded in Tables 6 and 7, we find that although our model requires slightly more computing time than the LC model and MC model when the resolution 256×256 is larger, it produces more satisfactory registration effects than them. In addition, from the above tables we can also see that when the resolution ≤128×128 of the image is relatively small, our model can produce more satisfactory registration effects than the other models in a relatively short time. Although the ZC model and Hyper model produce slightly better results than our new model for the image with the size of 256×256, our new model is more robust with respect to mesh parameter h. As can be seen from Tables 8–10, with the increase of the regularizer parameter α, the registration effects gradually become worse. But when α is taken in a certain range 0.02846,0.1, our model can still produce satisfactory registration effects. From the tables above, we also find that when the resolution of the image changes and the regularizer parameter α is taken in a certain range, the change of registration results generated by our new model is not particularly significant. This shows that our new model is more robust.Figure 3
Registration results of testing a pair of medical lung images using our new model. (a) Reference image, (b) template image, (c) the transformed template image using our new model, (d) image difference before registrationε=100%, (e) the transformation x+ux, and (f) image difference after registration ε=1.23%.
(a)(b)(c)(d)(e)(f)Figure 4
Registration results for a pair of lung images with size of256×256 using other five models. The transformed template image Tx+ux from (a) the LC model [28], (b) Hyper model [20], and (c) DDemons model [43]. The image difference after registration from (d) the LC model [28] ε=2.41%, (e) Hyper model [20] ε=1.01%, and (f) DDemons model [43] ε=16.74%.The transformed template image Tx+ux from (g) the ZC model [42] and (h) MC model [8]. The image difference after registration from (i) the ZC model [42] ε=1.22% and (j) MC model [8] ε=2.15%.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)Table 6
Quantitative measurements using our new model and LC model [25] for processing Example 2 shown in Figures 3(a) and 3(b).
LayerModelLCCLCα=4.84e−2α=2.846e−2ε%Tℱε%Tℱh=1/2562.4114.50.00061.2344.10.0186h=1/1281.967.60.09100.768.90.2540h=1/641.174.00.13630.432.60.4581h=1/321.072.60.25790.291.10.5437h=1/160.811.2−0.32570.130.50.6884T represents the total run time which includes the output of the image (seconds). ℱ>0 means that the transformation does not include folding or cracking.Table 7
Quantitative measurements using the DDemons model [43], ZC model [42], Hyper model [20], and MC model [8] for processing Example 2 shown in Figures 3(a) and 3(b).
LayerModelDDemonsZCHyperMCσfluid=2, σdiff=2.336α=2, β=8α=10, β=260α=8e−2ε%Tℱε%Tℱε%Tℱε%Tℱh=1/25616.7412.10.27911.2252.10.00131.0113.60.02152.1515.40.0080h=1/12816.176.70.07191.689.40.00291.143.60.42661.525.60.0844h=1/6422.723.80.49602.384.60.00911.461.70.48590.873.10.2770h=1/3220.294.30.82794.341.00.01483.181.20.59130.502.00.4201h=1/1635.673.00.946416.480.40.045926.010.30.65720.371.20.5814T represents the total run time which includes the output of the image(seconds). ℱ>0 means that the transformation does not include folding or cracking.Table 8
Dependent results of our new model on the regularizer parameterα by testing Example 2 (Figures 3(a) and 3(b)).
α0.028460.030.040.050.060.07ε (%)1.22821.23641.34631.42681.49261.5496ℱ0.01860.00420.08940.18360.25240.3152ℱ>0 means that the transformation does not include folding or cracking.Table 9
Dependent results of our new model on the regularizer parameterα by testing Example 2 (Figures 3(a) and 3(b)).
α0.10.20.51510ε (%)1.68501.96002.32502.62873.64704.3160ℱ0.40190.51920.69010.73520.84250.8998ℱ>0 means that the transformation does not include folding or cracking.Table 10
Dependent results of our new model on the regularizer parameterα by testing Example 2 (Figures 3(a) and 3(b)).
α501005001000500010000ε (%)6.98378.952315.318818.216621.337923.6637ℱ0.99030.98870.98320.98180.98200.9802ℱ>0 means that the transformation does not include folding or cracking.


### 4.3. Test 3: A Pair of Artificial Image

A pair of artificial image with a size of256×256 is used in the third experiment. For each model, we still choose its optimal parameters. Figure 5 represents the test images and registration results using our new model (CLC). The transformed template images and image differences after registration from all other models are shown in Figure 6. The registration results on different layers using our new model (CLC) and LC model [28] and other four models are summarized in Tables 11 and 12, respectively. In Figures 5 and 6, we can see that all five models except the DDemons model are fine in producing satisfactory registration results. According to Tables 11 and 12, we find that although our new model requires slightly more computing time when the resolution 256×256 is large, it has the best value of ε. We also find that LC model and MC model have mesh folding during the registration process. In addition, our new model is more robust than the ZC model and Hyper model with the change of grid parameter h. For this example, an accurate regularizer parameter α is also not needed. From Tables 13 to 15, we also find our proposed new model can produce acceptable registration results for any α between 0.1 and 1.Figure 5
Registration results of testing a pair of artificial images using our new model. (a) Reference image, (b) template image, and (c) the transformed template image using our new model (CLC). The image difference (d) before registrationε=100%, (e) the transformation, and (f) after registration using our new model (CLC) ε=0.07%, respectively.
(a)(b)(c)(d)(e)(f)Figure 6
Registration results for a pair of lung images with size of256×256 using other five models. The transformed template image Tx+ux from (a) the LC model [28], (b) Hyper model [20], and (c) DDemons model [43]. The image difference after registration from (d) the LC model [28] ε=0.73%, (e) Hyper model [20] ε=0.10%, and (f) DDemons model [43] ε=5.26%. The transformed template image Tx+ux from (g) ZC model [42] and (h) MC model [8]. The image difference after registration from (i) the ZC model [42] ε=0.12% and (j) MC model [8] ε=0.19%.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)Table 11
Quantitative measurements using our new model (CLC) and LC model [25] for processing Example 3.
LayerModelLCCLCα=4.2α=1e−1ε%Tℱε%Tℱh=1/2560.7311.00.46830.0744.00.4410h=1/1280.885.10.49030.0810.50.4554h=1/641.584.00.50740.182.30.2405h=1/329.402.6−0.51340.231.10.3042h=1/169.431.7−0.08341.380.40.2617T represents the total run time which includes the output of the image(seconds). ℱ>0 means that the transformation does not include folding or cracking.Table 12
Quantitative measurements using the DDemons model [43], ZC model [42], Hyper model [20], and MC model [8] for processing Example 3 shown in Figures 5(a)–5(c).
LayerModelDDemonsZCHyperMCσfluid=1, σdiff=0.4α=35, β=40α=240, β=300α=4e−1ε%Tℱε%Tℱε%Tℱε%Tℱh=1/2565.2629.60.05810.121.80.14360.103.30.42330.199.50.4524h=1/1281.9624.10.04990.381.20.13760.522.30.51670.195.60.4416h=1/641.394.70.30792.350.90.05411.621.80.52160.783.8−0.5506h=1/329.347.40.44547.240.60.06336.040.60.56400.672.5−0.5497h=1/165.446.30.697912.170.50.068236.820.50.59400.661.60.1492T represents the total run time which includes the output of the image (seconds). ℱ>0 means that the transformation does not include folding or cracking.Table 13
Dependent results of our new model on the regularizer parameterα by testing Example 3 (Figures 5(a)–5(c)).
α0.10.20.30.40.50.6ε (%)0.06730.11040.13980.16340.18390.2025ℱ0.44100.45080.45660.46100.46470.4679ℱ>0 means that the transformation does not include folding or cracking.Table 14
Dependent results of our new model on the regularizer parameterα by testing Example 3 (Figures 5(a)–5(c)).
α0.70.80.91.02.03.0ε (%)0.21980.23640.25230.26780.41180.5519ℱ0.47070.47320.47170.47050.46740.4690ℱ>0 means that the transformation does not include folding or cracking.Table 15
Dependent results of our new model on the regularizer parameterα by testing Example 3 (Figures 5(a)–5(c)).
α4.05.06.010.020.030.0ε (%)0.69380.83860.98591.59803.03175.0536ℱ0.47110.47170.47210.47710.47910.5214ℱ>0 means that the transformation does not include folding or cracking.


## 4.1. Test 1: A Pair of Lena Images

In the first experiment, we used a pair of Lena images with a size of256×256. The test images and registered ones using our new model are shown in Figure 1. The transformed template image obtained using the other five models and the image difference after registration are represented by Figure 2. For Example 1 shown in Figures 1(a) and 1(b), the registration results using our new model and linear curvature (LC [28]) model and other four models are recorded in Tables 1 and 2, respectively. For a smaller regularizer parameter α=3.12e−3, from Figures 1 and 2, we can see that our new model and hyperelastic regularizer-, linear curvature-, and mean one-based image registration models can produce visually satisfactory registration results. However, we find that the latter two have mesh folding by Tables 1 and 2. Although the demons-based registration model and ZC model do not produce folding, the registration effects are relatively poor. From Tables 1 and 2, we can further see that for the same image with different sizes, our new model can give very satisfactory registration results without folding.Figure 1
Registration results of testing Lena images using our new model. (a) Reference image, (b) template image, (c) the transformed template image using our new model, (d) image difference before registrationε=100%, (e) the transformation x+ux, and (f) image difference after registration ε=0.36%.
(a)(b)(c)(d)(e)(f)Figure 2
Registration results for a pair of Lena images with a size of256×256 using other five models. The transformed template image Tx+ux from (a) the LC model [28], (b) Hyper model [20], and (c) DDemons model [43]. The image difference after registration from (d) the LC model ε=1.48% [28], (e) Hyper model [20] ε=0.52%, and (f) DDemons model [43] ε=2.12%. The transformed template image Tx+ux from (g) the ZC model [42] and (h) MC model [8]. The image difference after registration from (i) the ZC model [42] ε=3.48% and (j) MC model [8] ε=1.41%.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)Table 1
Quantitative measurements using our new model and LC model [28] for processing Example 1 shown in Figures 1(a) and 1(b).
LayerModelLCCLCα=3.12e−3α=3.12e−3ε%Nℱε%Nℱh=1/2561.48315−1.22280.3600.5791h=1/1281.1692−0.88560.2200.4997h=1/641.0420−0.71840.3900.5946h=1/320.664−0.37300.4000.0168h=1/161.231−0.56291.0400.0486N represents the number of mesh folding of the transformation x+ux. ℱ>0 means that the transformation does not include folding or cracking.Table 2
Quantitative measurements using the DDemons model [43], ZC model [42], Hyper model [20], and MC model [8] for processing Example 1 shown in Figures 1(a) and 1(b).
LayerModelDDemonsZCHyperMCσfluid=1, σdiff=0.756α=3.12e−3α=3.12e−3α=3.12e−3ε%Nℱε%Nℱε%Nℱε%Nℱh=1/2562.1200.04343.4800.06490.5200.30381.418−0.1055h=1/1284.6300.17643.2900.06601.0500.64971.141−0.0273h=1/6411.9900.54326.4300.11914.1200.71163.0400.1613h=1/3223.6900.689910.3300.044015.8400.71932.912−0.2810h=1/1648.1400.865028.4500.379285.5500.82472.1700.1656N represents the number of mesh folding of the transformation x+ux. ℱ>0 means that the transformation does not include folding or cracking.In order to analyze how much our model will be affected when changing the regularizer parameterα, we use Algorithm 4 to test Example 1, and the corresponding quantitative measurement values are recorded in Tables 3–5. From Tables 3 to 5, we find that the registration quality becomes worse and worse with the increase of the regularizer parameter α. However, when alpha is greater than or equal to 3.12×10−3, our new model does not produce folding for Example 1. For this example, our new model is able to produce better registration results, especially when the regularizer parameter α is taken in an appropriate range 3.12×10−3,3×10−1.Table 3
Dependent results of our new model on the regularizer parameterα by testing Example 1 (Figures 1(a) and 1(b)).
α=e−33.124581030ε (%)0.35740.37340.38760.41690.43070.4996ℱ0.57910.59680.61190.63870.64930.6826ℱ>0 means that the transformation does not include folding or cracking.Table 4
Dependent results of our new model on the regularizer parameterα by testing Example 1 (Figures 1(a) and 1(b)).
α=e−380100200300400500ε (%)0.58360.61180.74880.89031.03801.1914ℱ0.67360.67190.66780.66690.66790.6699ℱ>0 means that the transformation does not include folding or cracking.Table 5
Dependent results of our new model on the regularizer parameterα by testing Example 1 (Figures 1(a) and 1(b)).
α=e−31e32e33e34e35e36e3ε (%)2.02853.84245.67387.45499.20010.8363ℱ0.68950.73500.75340.76880.77980.7878ℱ>0 means that the transformation does not include folding or cracking.


## 4.2. Test 2: A Pair of Medical Lung Images

In this part, we select a pair of medical lung images with a size of256×256 for testing. For each model, we choose the optimal parameters. The test images and the corresponding registered ones are shown in Figures 3 and 4, respectively. The registration results on different layers are recorded in Tables 6 and 7, respectively. According to Figures 3 and 4, we can observe that the other five models can produce visually relatively satisfactory registration effects except the DDemons model. From the image difference after registration and the results recorded in Tables 6 and 7, we find that although our model requires slightly more computing time than the LC model and MC model when the resolution 256×256 is larger, it produces more satisfactory registration effects than them. In addition, from the above tables we can also see that when the resolution ≤128×128 of the image is relatively small, our model can produce more satisfactory registration effects than the other models in a relatively short time. Although the ZC model and Hyper model produce slightly better results than our new model for the image with the size of 256×256, our new model is more robust with respect to mesh parameter h. As can be seen from Tables 8–10, with the increase of the regularizer parameter α, the registration effects gradually become worse. But when α is taken in a certain range 0.02846,0.1, our model can still produce satisfactory registration effects. From the tables above, we also find that when the resolution of the image changes and the regularizer parameter α is taken in a certain range, the change of registration results generated by our new model is not particularly significant. This shows that our new model is more robust.Figure 3
Registration results of testing a pair of medical lung images using our new model. (a) Reference image, (b) template image, (c) the transformed template image using our new model, (d) image difference before registrationε=100%, (e) the transformation x+ux, and (f) image difference after registration ε=1.23%.
(a)(b)(c)(d)(e)(f)Figure 4
Registration results for a pair of lung images with size of256×256 using other five models. The transformed template image Tx+ux from (a) the LC model [28], (b) Hyper model [20], and (c) DDemons model [43]. The image difference after registration from (d) the LC model [28] ε=2.41%, (e) Hyper model [20] ε=1.01%, and (f) DDemons model [43] ε=16.74%.The transformed template image Tx+ux from (g) the ZC model [42] and (h) MC model [8]. The image difference after registration from (i) the ZC model [42] ε=1.22% and (j) MC model [8] ε=2.15%.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)Table 6
Quantitative measurements using our new model and LC model [25] for processing Example 2 shown in Figures 3(a) and 3(b).
LayerModelLCCLCα=4.84e−2α=2.846e−2ε%Tℱε%Tℱh=1/2562.4114.50.00061.2344.10.0186h=1/1281.967.60.09100.768.90.2540h=1/641.174.00.13630.432.60.4581h=1/321.072.60.25790.291.10.5437h=1/160.811.2−0.32570.130.50.6884T represents the total run time which includes the output of the image (seconds). ℱ>0 means that the transformation does not include folding or cracking.Table 7
Quantitative measurements using the DDemons model [43], ZC model [42], Hyper model [20], and MC model [8] for processing Example 2 shown in Figures 3(a) and 3(b).
LayerModelDDemonsZCHyperMCσfluid=2, σdiff=2.336α=2, β=8α=10, β=260α=8e−2ε%Tℱε%Tℱε%Tℱε%Tℱh=1/25616.7412.10.27911.2252.10.00131.0113.60.02152.1515.40.0080h=1/12816.176.70.07191.689.40.00291.143.60.42661.525.60.0844h=1/6422.723.80.49602.384.60.00911.461.70.48590.873.10.2770h=1/3220.294.30.82794.341.00.01483.181.20.59130.502.00.4201h=1/1635.673.00.946416.480.40.045926.010.30.65720.371.20.5814T represents the total run time which includes the output of the image(seconds). ℱ>0 means that the transformation does not include folding or cracking.Table 8
Dependent results of our new model on the regularizer parameterα by testing Example 2 (Figures 3(a) and 3(b)).
α0.028460.030.040.050.060.07ε (%)1.22821.23641.34631.42681.49261.5496ℱ0.01860.00420.08940.18360.25240.3152ℱ>0 means that the transformation does not include folding or cracking.Table 9
Dependent results of our new model on the regularizer parameterα by testing Example 2 (Figures 3(a) and 3(b)).
α0.10.20.51510ε (%)1.68501.96002.32502.62873.64704.3160ℱ0.40190.51920.69010.73520.84250.8998ℱ>0 means that the transformation does not include folding or cracking.Table 10
Dependent results of our new model on the regularizer parameterα by testing Example 2 (Figures 3(a) and 3(b)).
α501005001000500010000ε (%)6.98378.952315.318818.216621.337923.6637ℱ0.99030.98870.98320.98180.98200.9802ℱ>0 means that the transformation does not include folding or cracking.


## 4.3. Test 3: A Pair of Artificial Image

A pair of artificial image with a size of256×256 is used in the third experiment. For each model, we still choose its optimal parameters. Figure 5 represents the test images and registration results using our new model (CLC). The transformed template images and image differences after registration from all other models are shown in Figure 6. The registration results on different layers using our new model (CLC) and LC model [28] and other four models are summarized in Tables 11 and 12, respectively. In Figures 5 and 6, we can see that all five models except the DDemons model are fine in producing satisfactory registration results. According to Tables 11 and 12, we find that although our new model requires slightly more computing time when the resolution 256×256 is large, it has the best value of ε. We also find that LC model and MC model have mesh folding during the registration process. In addition, our new model is more robust than the ZC model and Hyper model with the change of grid parameter h. For this example, an accurate regularizer parameter α is also not needed. From Tables 13 to 15, we also find our proposed new model can produce acceptable registration results for any α between 0.1 and 1.Figure 5
Registration results of testing a pair of artificial images using our new model. (a) Reference image, (b) template image, and (c) the transformed template image using our new model (CLC). The image difference (d) before registrationε=100%, (e) the transformation, and (f) after registration using our new model (CLC) ε=0.07%, respectively.
(a)(b)(c)(d)(e)(f)Figure 6
Registration results for a pair of lung images with size of256×256 using other five models. The transformed template image Tx+ux from (a) the LC model [28], (b) Hyper model [20], and (c) DDemons model [43]. The image difference after registration from (d) the LC model [28] ε=0.73%, (e) Hyper model [20] ε=0.10%, and (f) DDemons model [43] ε=5.26%. The transformed template image Tx+ux from (g) ZC model [42] and (h) MC model [8]. The image difference after registration from (i) the ZC model [42] ε=0.12% and (j) MC model [8] ε=0.19%.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)Table 11
Quantitative measurements using our new model (CLC) and LC model [25] for processing Example 3.
LayerModelLCCLCα=4.2α=1e−1ε%Tℱε%Tℱh=1/2560.7311.00.46830.0744.00.4410h=1/1280.885.10.49030.0810.50.4554h=1/641.584.00.50740.182.30.2405h=1/329.402.6−0.51340.231.10.3042h=1/169.431.7−0.08341.380.40.2617T represents the total run time which includes the output of the image(seconds). ℱ>0 means that the transformation does not include folding or cracking.Table 12
Quantitative measurements using the DDemons model [43], ZC model [42], Hyper model [20], and MC model [8] for processing Example 3 shown in Figures 5(a)–5(c).
LayerModelDDemonsZCHyperMCσfluid=1, σdiff=0.4α=35, β=40α=240, β=300α=4e−1ε%Tℱε%Tℱε%Tℱε%Tℱh=1/2565.2629.60.05810.121.80.14360.103.30.42330.199.50.4524h=1/1281.9624.10.04990.381.20.13760.522.30.51670.195.60.4416h=1/641.394.70.30792.350.90.05411.621.80.52160.783.8−0.5506h=1/329.347.40.44547.240.60.06336.040.60.56400.672.5−0.5497h=1/165.446.30.697912.170.50.068236.820.50.59400.661.60.1492T represents the total run time which includes the output of the image (seconds). ℱ>0 means that the transformation does not include folding or cracking.Table 13
Dependent results of our new model on the regularizer parameterα by testing Example 3 (Figures 5(a)–5(c)).
α0.10.20.30.40.50.6ε (%)0.06730.11040.13980.16340.18390.2025ℱ0.44100.45080.45660.46100.46470.4679ℱ>0 means that the transformation does not include folding or cracking.Table 14
Dependent results of our new model on the regularizer parameterα by testing Example 3 (Figures 5(a)–5(c)).
α0.70.80.91.02.03.0ε (%)0.21980.23640.25230.26780.41180.5519ℱ0.47070.47320.47170.47050.46740.4690ℱ>0 means that the transformation does not include folding or cracking.Table 15
Dependent results of our new model on the regularizer parameterα by testing Example 3 (Figures 5(a)–5(c)).
α4.05.06.010.020.030.0ε (%)0.69380.83860.98591.59803.03175.0536ℱ0.47110.47170.47210.47710.47910.5214ℱ>0 means that the transformation does not include folding or cracking.


## 5. Conclusions

Avoiding mesh folding is a key issue to ensure the invertibility of transformation in the image registration model. In this paper, we propose a constrained linear curvature image registration model by integrating the evaluation criteria to measure the registration results directly into the basic framework of the variational model. In order to solve the new model, we use a combination of the multiplier method and Gauss–Newton scheme with the Armijos line search and further combine with a multilevel method to achieve fast convergence. To illustrate the good performance of our new model, we compare it with five representative models based on the LC model [28], MC model [8], DDemons model [43], Hyper model [20], and ZC model [42] using three numerical examples. Numerical experiments show that our proposed model is more efficient and robust than the competing models. Future works will consider the use of homotopy method to solve the corresponding inequality constraint registration model.



---

*Source: 9468215-2020-10-13.xml*