Constrained Linear Curvature Image Registration Model and Its Numerical Algorithm

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.


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. us, 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][2][3][4][5].
e basic framework of image registration can be described as follows: given two images of the same object, which are called reference image R and template image T, respectively, and our purpose is to find a vector value transformation φ as defined below: or equivalently find the unknown displacement field u: so that the transformed template image T(φ(x)) � T(x + u(x)) ≜ T(u) is as similar to the reference image R as possible. Here, d ∈ N denotes spatial dimension of the given images. Without loss of generality, here we focus on d � 2 throughout this paper, but it is easy to generalize to d � 3 with some additional modifications. e variational model is an important tool for studying image registration and has been widely concerned by many researchers [5][6][7][8][9]. is variational model treats the image registration problem as a minimization problem of the joint energy functional in the following form: where where D(u) denotes distance measure which quantifies distance or similarity of the transformed template image T(u) and reference R, for other choices on D(u), refer to [5,7], S(u) 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 u(x) to the parametric model [10][11][12][13][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][16][17][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][20][21], diffusion regularizer [22], total variational regularizer [23,24], modified total variational regularizer [25,26], total fractionalorder 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. erefore, 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 highorder 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 + u(x) found is reversible. e 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. erefore, 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 C(u) on the transformed Jacobian matrix determinant in the objective function formula (3) [33][34][35][36], i.e., where C(u) � det(I d + ∇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][39][40][41][42][43]. However, some of them require more computation time due to the complexity of the regularizer. ere are two purposes for image registration. One is to enhance some similarities between two images by geometrically transforming one of the given two images. e other is to ensure that this transformation is reasonable. In fact, it is equivalent to find geometric transformation φ and the displacement field u(x) 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). e rest of the paper is organized as follows: in Section 2, we propose a new constrained linear curvature image registration model. en 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.

Constrained Linear Curvature Image Registration Model
Firstly, we briefly review the Fischer-Modersitzki's linear curvature image registration model [25,27]. Choosing S(u) in (3) based on an approximation to the curvature of the surface of the displacement field u l is given by the following form: ere 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: where C(u) � det(J(φ(x))) 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. e numerical solution of the new model (7) is given below.

Numerical Solution of the New Model
In general, it is difficult to solve the optimization problem (7) by the analytic method. us 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.1. Discretization of Regularizer. e discrete form of the continuous displacement field u � (u 1 , u 2 ) ⊤ can be represented by u h � (u h 1 , u h 2 ), where u h 1 and u h 2 are the discrete grid functions defined on the discrete region Ω h . For conve- . . , m 2 , 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. e discrete gradient operator ∇ h can be defined at each pixel (i, j) by the following form: where e displacement field u satisfies the homogeneous Neumann boundary conditions on the boundary zΩ of the image region Ω: rough the analysis of continuous setting, we know that the discrete divergence operator is the negative conjugate transposition of the gradient operator, namely, ∇· � −∇ * . us, it can be defined by the following form: where ω � (ω 1 , ω 2 ) is a vector. For the convenience of calculation, the grid functions u h 1 and u h 2 can be changed into column vectors u h 1 and u h 2 according to lexicographical ordering, respectively:

Mathematical Problems in Engineering
We can also be expressed as the product of matrix A ⊤ k ∈ R 2×N (k � 1, 2, . . . , N) and the vector u h l (l � 1, 2) in the following form: By this notation, we can get en the discrete form of (18) is as follows: According to the midpoint quadrature formula, the linear curvature regularizer S LC (u) � (1/2) Ω B[u]d Ω has the following discrete form: where

Discretization of Template T 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 cubic B-spline is used for interpolation. Next, T and R are used to represent the continuous smooth approximation of template image T and reference image R, respectively. Let and us the discrete reference image and transformed template image can be represented by the following form, respectively: and further we can get the Jacobian of T → : where U h c � X h c + U h and the Jacobian of T → is a block matrix with diagonal blocks.

Discretization of Distance Measure D.
Although it is in a continuous setting, it is not possible to compute integrals analytically. us 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 D(2) can be written directly as follows: In addition, the derivative of the discrete functional D h (U h ) on U h can also be calculated and has the following form: Furthermore, we can calculate the second derivative 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 → (U h ) 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], d 2 D h (U h ) can be approximated by the following form:

Discretization of Inequality Constraint Functional
C(u). In model (7), the inequality constraint functional C(u) is defined by According to the previous analysis, the discrete form of the partial derivative of the continuous displacement field element u l can be expressed as follows: Obviously, m l ∈ R N , and w l ∈ R N , where N � n 1 × n 2 . Let where symbol ⊛ denotes the multiplication of the corresponding elements of two vectors and c ∈ R N . For the convenience of calculation, let c i denote the i-th element of c, i � 1, 2, . . . , N.
erefore, the continuous inequality constraint function C(u) has the following discrete form: Since the first-order variation of the continuous inequality constraint function C(u) on continuous displacement field u is as follows: us we can get the discrete form of the first-order variation dC(u): Obviously, 3.2. Solving the Discrete Optimization Problem. According to the above analysis, inequality constrained functional (7) has the following discrete form: Below we use the multiplier method to numerically solve the inequality constrained optimization problem (35). e 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.

Multiplier Method for Inequality Constrained
Problems. Consider the following inequality constrained optimization problem: Let y i ≥ 0, and the above inequality constraint can be transformed into the following equivalent equality constraint problem: In this case, the augmented Lagrange function can be expressed as In order to eliminate the auxiliary variables y, the minimization of φ with respect to variable y can be considered. According to the first-order necessary condition, let We can get Namely, Mathematical Problems in Engineering that is to say, And when σg i (x) − λ i > 0, we can obtain According to the above two cases, Substituting it into formula (38), we can get the corresponding augmented Lagrange function of (36): 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 x k to denote the approximate minimizer of φ(x, y, λ, σ), then we can get by the optimality conditions for unconstrained minimization that Let (x * , y * , λ * ) satisfy the KKT conditions for (37), then we have By comparing (48) with (49), we can deduce that 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: en, taking (43) into the multiplier iteration formula (51), we have Furthermore, it can be written as Similarly, take (43) into the termination criterion We can get

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): e corresponding multiplier iteration formula has the following form: And the corresponding stopping criterion is Although the augmented Lagrangian function (57) of the model (35) contains the min function, it is still continuously differentiable; for details, see [37,44]. e detailed steps of the multiplier method for solving the model (35) can be summarized by Algorithm 1.
In Algorithm 1, 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 U h(k) of the previous step by the Taylor expansion given below: where dψ (U h(k) ) is the Jacobian matrix of ψ at U h(k) and H is the approximation of its Hessian. Due to ) are both positive semidefine, and it is easy to prove that H is also positive semidefinite. us, 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]. e detailed steps are described below.
Given the initial value U h(k) , the Jacobian matrix dψ (U h(k) ) and the Hessian matrix H are calculated by the following forms: and in each outer iteration step, respectively, where M(U h(k) ) � dC h (U h ) ∈ R 2N is defined by (34). In formula (59), the disturbance value δ U h can be obtained by finding the stability point of the quadratic function ψ, namely, Usually, H is the positive definite. us 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. e 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 U h(0) 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 I h H is used. Finally, this process is repeated until the loop terminates. We summarize the multilevel method using Algorithm 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: And the minimum value F of the determinant of the Jacobian matrix J of the transformation φ are used, where D(u) is defined by equation (4), u is the current iteration, and D stop is the value of D(u) at u � 0.  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. 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
Step 2: solving the subproblem. With U h(k− 1) as the initial point, solve the minimum value U h(k) 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 > max k, where β k is defined by (57), the iteration is stopped, and U h(k) is output as the approximate minimum of the original problem; otherwise, go to Step 4.

then the iteration is stopped; otherwise, go to
Step 4.
Step 5: the step factor t is solved by the Armijos line search technique.
Step 3: the initial value U h(0) on the finer level is obtained by interpolation operator I h H .

Test 2: A Pair of Medical Lung Images.
In this part, we select a pair of medical lung images with a size of 256 × 256 for testing. For each model, we choose the optimal parameters. e test images and the corresponding registered ones are shown in Figures 3 and 4, respectively. e 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. is shows that our new model is more robust.

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 T represents the total run time which includes the output of the image(seconds). F > 0 means that the transformation does not include folding or cracking.    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.

Data Availability
e image datasets used to support the findings of this study are included within the article.

Conflicts of Interest
e author declares no conflicts of interest.