## 1 Introduction

Image registration along with image segmentation are two of the most important tasks in imaging sciences problems. Here we focus on image registration. Much research work has been done; for an extensive overview of registration techniques see [28, 18, 20]

. The methods can be classified into parametric and non-parametric image registration based on the geometric transformation. For the first category, the transformation is governed by a finite set of image features or by expanding a transformation in terms of basis functions. The second category (which is our main concern in this paper) of non-parametric image registration methods is not restricted to a certain parametrisable set. The problem is modelled as a functional minimisation via the calculus of variations. Given two images, the reference

and template , the functional consists of a distance measure and a regularisation term whereis the sought displacement vector at pixel

. The term removes the ill-posedness of the minimisation problem. We use the squared norm of the distance measure to quantify the differences between and as follows(1) |

The distance measure in equation (1) is the sum of the squared difference (SSD) which is commonly used and optimal for mono-modal image registration with Gaussian noise. For multi-modal image registration where cannot be compared directly, other distance measures must be used [8]. Generally, the regularisation terms are inspired by physical processes such as elasticity, diffusion and motion curvature. As such, elastic image registration was introduced in [1] which assumed that objects are deformed as a rubber band.

In previous works, higher order regularisation models [7, 3] were found to be the most robust while the diffeomorphic demon model [22] offers the most physical transform in terms of (nearly) bijective mapping. Diffusion and total variation regularisation models based on first order derivatives are less complicated to implement but are at a disadvantage compared to higher order regularisation models based on second order derivatives due to two reasons. First, the former methods penalise rigid displacement. They cannot properly deal with transformations with translation and rotation. Second, low order regularisation is less effective than high order one in producing smooth transformations which are important in several applications including medical imaging. The work of [5, 6, 7] proposed a high order regularisation model named as a linear curvature, which is an approximation of the surface (mean) curvature and the model is invariant to affine registration. This work was later refined in [11, 9, 10] where a related approximation to the sum of the squares of the principal curvatures was suggested and higher order boundary conditions were recommended. Without any approximation to the mean curvature, the works in [3, 2] developed useful numerical algorithms for models based on the nonlinear mean curvature regularisation and observed advantages over the linear curvature models for image registration; however the effect of mesh folding (bijective maps) was not considered. The diffeomorphic demon model [23] is widely used due to its property of bijective maps; however the bijection is not precisely imposed. Another useful idea of enforcing bijection, beyond the framework we consider, is via minimising the Beltrami coefficient which measures the distortion of the quasi-conformal map of registration transforms [12].

In this paper we propose a high order registration model based on Gaussian curvature and hope to achieve large and smooth deformation without mesh folding. Although the Gaussian curvature is closely related to the mean curvature, it turns out our new model based on the Gaussian curvature is much better. The motivation of the proposed model comes from two factors. Firstly, we are inspired by the work of Elsey and Esedoglu [4] in geometry processing where Gaussian curvature of the image surface is used in a variational formulation. The authors proposed the Gaussian curvature as a natural analogue of the total variation of Rudin, Osher and Fatemi (ROF) [19] model in geometry processing. Aiming to generalise the ROF model to surface fairing where the convex shapes in 3D have similar interpretation to the monotone functions in 1D problems for the ROF model, they showed that, based on the Gauss Bonnet theorem, the complete analogy of the total variation regularisation for surface fairing is the energy functional of the Gaussian curvature. A very important fact pointed out in [4] stated that the mean curvature of the surface is not a suitable choice for surface fairing because the model is not effective for preserving important features such as creases and corners on the surface (although the model is still effective for removing noise). Their claims are also supported by the work of [13] where the authors illustrated several advantages of Gaussian curvature over mean curvature and total variation in removing noise in 2D images. First, Gaussian curvature preserves important structures such as edges and corners. Second, only Gaussian curvature can preserve structures with low gradient. Third, the model is effective in removing noise on small scale features. Thus, we believe that Gaussian curvature is a more natural physical quantity of the surface than mean curvature. Here we investigate the potential of using Gaussian curvature to construct a high order regularisation model for non-parametric image registration of mono-modal images.

The outline of this paper is as follows. In §2 we review the existing models for non-parametric image registration with focus on the demon, linear and mean curvature models. In §3 we introduce the mathematical background of Gaussian curvature for surfaces. In §4 we introduce a Gaussian curvature model and a numerical solver to solve the resulting Euler-Lagrange equations. We show in §5 some numerical tests including comparisons. Finally, we discuss the parameters’ selection issue for our model in §6 and present our conclusions in §7.

## 2 Review of Non-parametric Image Registration

In image registration, given two images and (which are assumed to be compactly supported and bounded operators ), the task is to transform to match as closely as possible. Although we consider throughout this paper, with some extra modifications, this work can be extended to . In non-parametric image registration, the transformation is denoted by where is a vector valued function

where . To separate the overall mapping from the displacement field, we will define

where is the displacement field. Thus, finding is equivalent to finding . A non-parametric image registration model takes the form

(2) |

where the distance measure is given as in (1) and the choice of regulariser differentiates different models. Here is searched over a set of admissible functions that minimise in (2). Usually, the set is a linear subspace of a Hilbert space with Euclidean scalar product.

The force term , to be used by all models, is the gradient of (1) with respect to the displacement field

(3) |

which is non-linear. Different regularisation terms will lead to different non-parametric image registration models. Below we list three popular models selected for tests and comparisons.

Model LC. The first one is the linear curvature model by [6, 7, 15, 5], where

(4) |

This term is an approximation of the surface curvature through the mapping where

(5) |

when . The Euler Lagrange equation for (2) with as the regularisation term is given by a fourth order PDE

(6) |

with boundary conditions and the unit outward normal vector. The model consists of the second order derivative information of the displacement field which results in smoother deformations compared to those obtained using first order models based on elastic and diffusion energies. It is refined in [11, 9, 10]

with nonlinear boundary conditions. The affine linear transformation belongs to the kernel

which is not the case in elastic or diffusion registration.Model MC. Next is the mean curvature model [3, 2]

where and is as defined in (5). The Euler Lagrange equation for (2) with as the regularisation term is given by:

(7) |

with boundary condition . One can use the multigrid method to solve equation (7) as in [3]; refer also to [2] for multi-modality image registration work.

Model D. Finally Thirion [21] introduced the so-called demon registration method where every pixel in the image acts as the demons that force a pulling and pushing action in a similar way to what Maxwell did for solving the Gibbs paradox in thermodynamics. The original demon registration model is a special case of diffusion registration but it has been much studied and improved since 1998; see [17, 15, 25, 14]. The energy functional for the basic demon method is given by

(8) |

where is the current displacement field, and account for noise on the image intensity and the spatial uncertainty respectively. Equation (8) can be linearised using first order Taylor expansion,

(9) |

where is given by

for an efficient second order minimisation. The first order condition of (9) leads to the new update for

The additional use of for helps to achieve a nearly diffeormorphic transformation (mapping), where is the stationary velocity field of the displacement field ; see [24]. It should be remarked that the three main steps of the model cannot be combined into a single energy functional.

## 3 Mathematical Background of the Gaussian curvature

In differential geometry, the Gaussian curvature problem seeks to identify a hypersurface of as a graph over so that, at each point of the surface, the Gaussian curvature is prescribed. Let denote the Gaussian curvature which is a real valued function in . The problem is modelled by the following equation

(10) |

where is the first order derivative operator. Equation (10) is one of the Monge-Ampere equations. For , we have

(11) |

In [4], the authors define a regularisation term using the Gaussian curvature of a closed surface based on the Gauss-Bonnet theorem.

###### Theorem 3.1

Gauss-Bonnet Theorem. For a compact surface , we have

where is the length element to the surface and is the Euler characteristic of the surface.

Using this Theorem, it was shown in [4] that the complete analogy of the total variation regularisation for surface fairing is the energy functional of the Gaussian curvature. The analogous term , to the total variation of a function, that appears in the ROF model [19], is given by

where is the length element to the surface .

Gaussian curvature is one of the fundamental second order geometric properties of a surface. According to the Gauss’s Theorema Egregium, Gaussian curvature is intrinsic. For a local isometric mapping between two surfaces, Gaussian curvature remains invariant i.e. if and , then and the mapping is smooth and diffeomorphic.

We can also use a level set function to define the Gaussian curvature. Denote by the zero level set of the surface generated through the mapping . Then and where and . The Gaussian curvature of the level set is given by

(12) |

where , , is the Hessian matrix and is the adjoint matrix for . We have

Then,

This is why we set in equation (11). We shall use to measure the Gaussian curvature as in [4] for a monotonically increasing function (since the functional should be nonnegative).

## 4 Image Registration based on Gaussian Curvature

Before introducing our new image registration model, we first illustrate some facts to support the use of Gaussian curvature.

### 4.1 Advantages of a Gaussian curvature

The total variation and Gaussian curvature. We use the volume based analysis introduced in [13] to compare two denoising models, respectively based on Gaussian curvature and the total variation:

(13) |

(14) |

Consider, for each , the time change of the volume which is enclosed by the surface and the plane . Assume with either positive () or negative () at all points. Denote by the closed curve defined by the level set and accordingly by the 2D region enclosed by . The volume change in in time is given by

where is the area element. We now consider how changes from evolving (13) or (14).

If is the solution of equation (14), then from Gauss’ theorem

where is the length element and is the unit normal vector to the curve which is represented as . Then

where is the length of the curve . Furthermore, the volume variation in time is

where denotes the area of the region . We can see that the change in from to is proportional to the ratio . So, when this ratio is large (indicating possibly a noise presence), the total variation model reduces it and hence removes noise. However, important features of which have a large level set ratio are removed also and so not preserved by the total variation model (14).

Using similar calculations to before for the Gaussian curvature scheme (13), we have

(15) | ||||

From here, we observe that the quantity for the subdomain is dependent on the product of the variation and the Gaussian curvature on the level curve. The function in (15) controls and scales the speed of the volume change in contrast to the total variation scheme where depends only on the variation of the level curve. Consider a point where . Gaussian curvature based on two principal curvatures and where is the curvature of the level curve passing the point and is the curvature of the path which passes the point and is orthogonal to the level curve. If the Gaussian curvature on one level curve is zero then there is no change in regardless of variation on the level curves. In contrast, with the total variation, if there is a variation in the level curve, then there is a change in . Based on this observation, we believe that the Gaussian curvature model is better than the total variation model for preserving features on surfaces.

The mean curvature and Gaussian curvature. The mean curvature (MC) is also widely used. Next, we show that, though closely related, Gaussian curvature (GC) is better than mean curvature for surfaces in three ways.

First, Gaussian curvature is invariant under rigid and isometric transformations. In contrast, mean curvature is invariant under rigid transformations but not under isometric transformations. Rigid transformations preserve distance between two points while isometric transformations preserve length along surfaces and preserve angles between curves on surfaces. To illustrate invariance, consider a surface

whose Gaussian curvature and mean curvature are respectively

If we flip the surface upside down (isometric transformation) where , we will have the same value for the Gaussian curvature and a different value for the mean curvature. Thus, Gaussian curvature is invariant under isometric transformation.

Second, Gaussian curvature can be used to localise the tip of a surface better than mean curvature. Consider

When we compute the mean and Gaussian curvature for the surface, we see that the value given by the Gaussian curvature is sharper than that of the mean curvature. The highest point of the Gaussian curvature is better distinguished from its neighbourhood compared to the highest point of the mean curvature.

Third, Gaussian curvature can locate saddle points better than mean curvature. Take

as one example. We can again compute its mean and Gaussian curvatures analytically. The mean curvature for this surface appears complex where the largest value is not at the saddle point and the saddle point cannot be easily located. However, Gaussian curvature gets its highest value at the saddle point and is therefore able to accurately identify the saddle point within its neighbourhood.

In addition to these three examples and observations, a very important fact pointed out in [4] states that the mean curvature of the surface is not a suitable choice for surface fairing because the model is not effective for preserving important features such as creases and corners on the surface (although the model is effective for removing noise). This is true when we refer to surface fairing (surface denoising) but not necessarily true for 2D image denoising. From the recent works done in image denoising [4, 13], we observe several advantages of Gaussian curvature over total variation and mean curvature. Therefore, we conjecture that Gaussian curvature may outperform existing models in image registration. To our knowledge there exists no previous work on this topic.

### 4.2 The proposed registration model

Now we return to the problem of how to align or register two image functions . Let the desired and unknown displacement fields between and be the surface map where and with . We propose our Gaussian curvature based image registration model as

(16) |

where

The above model (16) leads to two Euler Lagrange equations:

(17) |

where

and boundary conditions , where is the normal vector at the boundary . Derivation of the resulting Euler-Lagrange equations for this model can be found in Appendix A. The equations are fourth order and nonlinear with anisotropic diffusion. Gradient descent type methods are not feasible and one way to solve them is by a geometric multigrid as in [3]. Here, instead, we present a fast and efficient way to solve the model (16) on a unilevel grid.

### 4.3 Augmented Lagrangian Method

The augmented Lagrangian method (ALM) is often used for solving constraint minimisation problems by replacing the original problem with an unconstrained problem. The method is similar to the penalty method where the constraints are incorporated in the objective functional and the problem is solved using alternating minimisation of each of the sub-problems. However, in ALM, there are additional terms in the objective functional known as Lagrange multiplier terms arising when incorporating the constraints. Similar works on the augmented Lagrangian method in image restoration can be found in [26, 27].

To proceed, we introduce two new dual variables and where and . Consequently we obtain a system of second order PDEs which are more amendable to effective solution.

We obtain the following refined model for Gaussian curvature image registration

and further reformulate to get the augmented Lagrangian functional

(18) | ||||

where are the Lagrange multipliers, the inner products are defined via the usual integration in and is a positive constant. We use an alternating minimisation procedure to find the optimal values of and where the process involves only two main steps.

Step 1. For the first step we need to update for any given . The objective functional is given by

This sub-problem can be solved using the following Euler Lagrange equations:

(19) |

where , and . We have a closed form solution for this step, if solving alternatingly, where

Similarly, we solve from

(20) |

where , and .

Step 2. For the second step we need to update for any given and with the following functional

Thus, we have the following Euler Lagrange equations:

(21) |

with Neumann boundary conditions . To solve equation (21), first, we linearise the force term using Taylor expansion

where

Second, we approximate and with

The discrete version of equation (21) is as follows

(22) |

where

is the discrete version of the Laplace operator and is the discrete version of

Third, we solve the system of equation (22) using a weighted pointwise Gauss Siedel method

where and we choose .

The iterative algorithm to solve (18) is now summarised as follows.

## 5 Numerical Results

We use two numerical experiments to examine the efficiency and robustness of the Algorithm 1 on a variety of deformations. To judge the quality of the alignment we calculate the relative reduction of the similarity measure

and the minimum value of the determinant of the Jacobian matrix of the transformation, denoted as

(23) |

We can observe that when , the deformed grid is free from folding and cracking.

All experiments were run on a single level. Experimentally, we found that works well for several types of images. As for the stopping criterion, we use for the residual of the Euler-Lagrange equations (19)-(21) and the maximum number of iterations is 30. Experiments were carried out using Matlab R2014b with Intel(R) core (TM) i7-2600 processor and 16G RAM.

### 5.1 Test 1: A Pair of Smooth X-ray Images

Images for Test 1 are taken from [16] where X-ray images of two hands of different individuals need to be aligned. The size of images are and the recovered transformation is expected to be smooth. The scaled version of the transformation and the transformed template image is given in Figures 1 (d) and (e) respectively. The transformation is smooth and the model is able to solve such a problem.

For comparison, the transformed template images for the diffeormorphic demon method, linear, mean and Gaussian curvatures are shown in Figures 2 (a), (b), (c) and (d) respectively. We can observe that there are some differences of these images inside the red boxes where only Gaussian curvature delivering the best result of the features inside the boxes.

Measure | Model D | Model LC | Model MC | GC | ||
---|---|---|---|---|---|---|

ML | SL | ML | SL | SL | SL | |

1.6 | 1.6 | 0.1 | 0.5 | 0.0001 | 0.0001 | |

Time (s) | 15.19 | 186.48 | 84.33 | 12.98 | 275.3 | 953.15 |

0.1389 | 0.1229 | 0.0720 | 0.3780 | 0.0964 | 0.0582 | |

0.0600 | 0.1082 | 0.3894 | 0.1973 | 0.6390 | 0.3264 |