Registration is a process of finding the optimal one-to-one correspondence between different data, such as images or surfaces. Applications can be found in various fields, including computer graphics, computer visions and medical imaging. For example, in medical imaging, finding accurate 1-1 correspondence between medical data is crucial for statistical shape analysis of the anatomical structures. While in computer graphics, surface registration is needed for texture mapping.
Different registration approaches have been developed. Existing algorithms can mainly be divided into three categories, namely, 1. landmark based registration, 2. intensity based registration and 3. hybrid registration using both landmark and intensity information. Landmark based registration computes a smooth 1-1 correspondence between corresponding data that matches important features. This kind of registration, with good feature alignment, is particularly crucial in medical imaging and computer graphics. For example, in computer graphics, landmark based registration is used to obtain the constrained texture mapping. The main advantage of the landmark based method is that larger deformations can be dealt with and intuitive user-interaction can be incorporated. Intensity based registration aims to match corresponding data without feature landmarks. Registration is usually obtained by matching intensity functions, such as image intensity for image registration or surface curvature for surface geometric registration. The main advantage of the intensity based registration is that more image information is taken into account and the delineation of feature landmarks is not required. However, it usually cannot cope with large geometric deformations. Recently, hybrid registration that combines landmark based and intensity based methods have gained increased attention. Hybrid approaches use both the landmark and intensity information to guide the registration. This type of approaches can usually obtain more accurate registration result, since the advantages of landmark based and intensity based registration can be combined. In this work, we will mainly focus on the landmark based registration and the hybrid registration.
Most existing algorithms can compute registration accurately and efficiently when the deformation is small. However, the registration problem becomes challenging when large deformations occur. Bijectivity can be easily lost and overlaps can usually be observed in the obtained registration. This causes inaccuracies in the registration. It is therefore necessary to develop an algorithm to obtain diffeomorphic registration with large deformations.
In this paper, we introduce a novel method to obtain diffeomorphic image or surface registrations via quasi-conformal maps, which can deal with large deformations. The key idea is to minimize an energy functional involving a Beltrami coefficient term, which measures the distortion of the quasi-conformal map. The Beltrami coefficient effectively controls the bijectivity and smoothness of the registration, even with very large deformations. By minimizing the energy functional, we obtain an optimal Beltrami coefficient associated to the desired registration, which is guaranteed to be bijective. Using the proposed algorithm, landmark-based registration between images or surfaces can be effectively computed. The obtained registration is guaranteed to be diffeomorphic (1-1 and onto), even with a large deformation or a large number of landmark constraints. The proposed algorithm can also be combined with matching intensity (such as image intensity or surface curvature) to improve the accuracy of the registration. Numerical results show that the combination of landmark constraints with intensity matching can significantly improve the accuracy of the registration. To test the effectiveness of the proposed algorithm, experiments have been carried out on both synthetic and real data. Results show that the proposed algorithm can compute diffeomorphic registration between images or surfaces effectively and efficiently.
In summary, the contributions of this paper are three-folded. Firstly, we propose a variational method to search for an optimized Beltrami coefficient associated to a diffeomorphic quasi-conformal map with large deformations, which minimizes the local geometric distortion. Secondly, we apply the model to compute the landmark based registration, which can deal with very large deformations and large amount of landmark constraints. Thirdly, we extend the landmark based registration model to a hybrid registration model, which combine both landmark and intensity information to obtain more accurate registration.
The rest of the paper is organized as follows. In section 2, we review some previous works closely related to this paper. In section 3, we describe some basic mathematical background related to our proposed model. In section 4, our proposed model for diffeomorphic registration with large deformations is explained in details. We describe the numerical implementation of the proposed algorithm in section 5. Experimental results are reported in section 6. Finally, we conclude our paper in section 7.
2 Previous works
In this section, we will review some related works closely related to this paper.
Intensity-based image registration has been widely studied. A comprehensive survey on the existing intensity-based image registration can be found in . One of the commonly used method is based on the variational approaches to minimize the intensity mismatching error. For example, Vercauteren et al.  proposed the diffeomorphic demons registration algorithm, which is a non-parametric diffeomorphic image registration algorithm based on Thirion’s demons algorithm. The basic idea is to adapt the optimization procedure underlying the demons algorithm to a space of diffeomorphic transformations. The obtained registration is smooth and bijective. Several algorithms for surface registration that matches geometric quantities, such as curvatures, have also been propsoed . For example, Lyttelton et al.  proposed an algorithm for surface parameterizations based on matching surface curvatures. Yeo et al.  proposed the spherical demons method, which adopted the diffeomorphic demons algorithm , to drive surfaces into correspondence based on the mean curvature and average convexity. Conformal surface registration, which minimizes angular distortions, has also been widely used to obtain a smooth 1-1 correspondence between surfaces [13, 5, 7, 6, 14, 15, 16, 17]. An advantage of conformal registrations is that they preserve local geometry well. Quasi-conformal surface registrations, which allows bounded amount of conformality distortion, have also been studied [18, 19, 20, 21]. For example, Lui et al.  proposed to compute quasi-conformal registration between hippocampal surfaces based on the holomorphic Beltrami flow method, which matches geometric quantities (such as curvatures) and minimizes the conformality distortion .
Landmark-based registration has also been widely studied and different algorithms have been proposed. Bookstein et al.  proposed to use a thin-plate spline regularization (or biharmonic regularization) to obtain a registration that matches landmarks as much as possible. Tosun et al.  proposed to combine iterative closest point registration, parametric relaxation and inverse stereographic projection to align cortical sulci across brain surfaces. These diffeomorphisms obtained can better match landmark features, although not perfectly. Wang et al. [27, 30, 28, 29] proposed to compute the optimized harmonic registrations of brain cortical surfaces by minimizing a compounded energy involving the landmark-mistmatching term [27, 30]. The obtained registration obtains an optimized harmonic map that better aligns the landmarks. However, landmarks cannot be perfectly matched, and bijectivity cannot be guaranteed under large number of landmark constraints. Later, Lin et al. 
propose a unified variational approach for registration of gene expression data to neuroanatomical mouse atlas in two dimensions that matches feature landmarks. Again, landmarks cannot be exactly matched. Inexact landmark-matching registrations are sometimes advantegous. In the case when landmark points/curves cannot be accurately delineated, this method is more tolerant of errors in labeling landmarks and gives better parameterization. In the situation when exact landmark matching is required, smooth vector field has been applied to obtain surface registration. Lui et al.[28, 29] proposed the use of vector fields to represent surface maps and reconstruct them through integral flow equations. They obtained shape-based landmark matching harmonic maps by looking for the best vector fields minimizing a shape energy. The use of vector fields to compute the registration makes the optimization easier, although it cannot describe all surface maps. An advantage of this method is that exact landmark matching can be guaranteed. Time dependent vector fields can also be used [8, 9, 10, 11, 12]. For example, Glaunés et al.  proposed to generate large deformation diffeomorphisms of a sphere, with given displacements of a finite set of template landmarks. The time dependent vector fields facilitate the optimization procedure, although it may not be a good representation of surface maps since it requires more memory. The computational cost of the algorithm is also expensive. Quasi-conformal mapping that matches landmarks consistently has also been proposed. Wei et al.  also proposed to compute quasi-conformal mappings for feature matching face registration. The Beltrami coefficient associated to a landmark points matching parameterization is approximated. However, either exact landmark matching or the bijectivity of the mapping cannot be guaranteed, especially when very large deformations occur.
Algorithms for hybrid registration, which combines both the landmark and intensity information to guide the registration, has also been proposed. For example, Christensen et al.  propsoed an algorithm for hybrid registration that uses both landmark and intensity information to guide the registration. The method utilizes the unidirectional landmark thin-plate spline (UL-TPS) registration technique together with a minimization scheme for the intensity difference to obtain good correspondence between images. Paquin et al.  proposed a registration method using a hybrid combination of coarse-scale landmark and B-splines deformable registration techniques. Chanwimaluang et al.  proposed a hybrid retinal image registration approach that combines both area-based and feature-based methods. Existing hybrid registration techniques can drive data into good correspondence when deformations are not too large. In this work, we propose a hybrid quasi-conformal registration method, which can deal with very large deformations.
3 Mathematical background
In this work, we apply quasi-conformal maps to obtain diffeomorphic registrations with large deformations. In this section, we describe some basic theories related to quasi-conformal geometry. For details, we refer to readers to .
A surface with a conformal structure is called a Riemann surface. Given two Riemann surfaces and , a map is conformal if it preserves the surface metric up to a multiplicative factor called the conformal factor. An immediate consequence is that every conformal map preserves angles. With the angle-preserving property, a conformal map effectively preserves the local geometry of the surface structure.
A generalization of conformal maps is the quasi-conformal maps, which are orientation preserving homeomorphisms between Riemann surfaces with bounded conformality distortion, in the sense that their first order approximations takes small circles to small ellipses of bounded eccentricity . Mathematically, is quasi-conformal provided that it satisfies the Beltrami equation:
for some complex-valued function satisfying . is called the Beltrami coefficient, which is a measure of non-conformality. It measures how far the map at each point is deviated from a conformal map. In particular, the map is conformal around a small neighborhood of when . Infinitesimally, around a point , may be expressed with respect to its local parameter as follows:
Obviously, is not conformal if and only if . Inside the local parameter domain, may be considered as a map composed of a translation to together with a stretch map , which is postcomposed by a multiplication of which is conformal. All the conformal distortion of is caused by . is the map that causes to map a small circle to a small ellipse. From , we can determine the angles of the directions of maximal magnification and shrinking and the amount of them as well. Specifically, the angle of maximal magnification is with magnifying factor ; The angle of maximal shrinking is the orthogonal angle with shrinking factor . Thus, the Beltrami coefficient gives us all the information about the properties of the map (See Figure 1).
The maximal dilation of is given by:
Quasiconformal mapping between two Riemann surfaces and can also be defined. Instead of the Beltrami coefficient, the Beltrami differential is used. A Beltrami differential on a Riemann surface is an assignment to each chart of an complex-valued function , defined on local parameter such that
on the domain which is also covered by another chart . Here, and .
An orientation preserving diffeomorphism is called quasi-conformal associated with if for any chart on and any chart on , the mapping is quasi-conformal associated with .
4 Proposed algorithm
In this section, we explain our proposed model for diffeomorphic registration with large deformation in details. The basic idea is to look for a quasi-conformal map to register two corresponding data, which can either be images or surfaces. The quasi-conformal map is obtained by minimizing an energy functional involving a Beltrami coefficient term, which measures the distortion of the quasi-conformal map. The Beltrami coefficient effectively controls the bijectivity and smoothness of the registration, even with very large deformations.
4.1 Proposed model
Let and be two corresponding images or surfaces. Our goal is to find a smooth 1-1 mapping between and satisfying certain prescribed criteria. For landmark-based registration, we look for a registration that matches corresponding feature landmarks. Let and be the sets of corresponding feature landmarks defined on and respectively. We search for a diffeomorphism subject to the landmark constraints that for all .
We propose a variational approach to obtain an optimized quasi-conformal map , which minimizes an energy functional involving the Beltrami coefficient terms. More specifically, we propose to solve the following minimization problem:
subject to the constraints that:
C(i) for (landmark constraint);
The first term of ensures the smoothness of . The second term of aims to minimize the conformality distortion of . The constraint C(i) is the landmark constraint, which enforces to match corresponding landmarks consistently.
If satisfies the constraint C(ii), then is a diffeomorphism. Suppose under some local coordinates. The Beltrami coefficient is given by:
Now, the Jacobian of , , is given by:
Since , . Also, . Hence, everywhere.
Since the Jacobian is postive everywhere, by the inverse function theorem, the mapping is locally invertible everywhere. In other words, is a diffeomorphism.
[Landmark-matching registration] Let:
Then: has a minimizer in . In fact, is compact.
Note that . In particular, the unique Teichmüller map is indeed in . We first prove that is complete. Let be a Cauchy sequence in under the norm . Then, is also a Cauchy sequence with respect to the norm. Since , which is complete. Thus, uniformly for some . Since is continuous, is also continuous. Besides, is convergent. Hence, and uniformly for some .
In addition, for all implies . implies . Since uniformly, locally uniformly. This implies for . Therefore, . Obviously, is totally bounded. We conclude that is compact.
Since is compact and is continuous in , has a minimizer in .
In order to improve the accuracy of the registration, one can combine the landmark-matching registration model with the intensity matching model. The intensities are functions defined on and . Usually, they are image intensities for image registration and surface curvatures for surface registration. Ideally, we want to obtain a landmark-matching diffeomorphism that matches the intensities as much as possible. We denote the intensities on and by and respectively. Our registration model can be modified as solving the following minimization problem:
subject to the constraints C(i) and C(ii).
[Landmark and intensity matching registration] has a minimizer in if and are continuous. The exisitence of minimizer depends on the continuity of and . If and are continuous, is continuous. Since is compact, has a minimizer in .
4.2 Energy minimization
In this subsection, we describe an algorithm to approximate the solutions of the above minimization problems.
4.2.1 Landmark based registration model
Given two corresponding sets of landmarks and on and respectively, our goal is to look for a diffeomorphism that satisfies () while minimizing the local geometric distortion. Our proposed model is to solve the variational problem (5) as described in the last subsection.
More specifically, our goal is to look for an optimal Beltrami coefficient , which is the Beltrami coefficient of some diffeomorphism , minimizing the following energy functional :
subject to the constraints that , for and , where is the Beltrami coefficient of .
We apply a splitting method to solve the constrained optimization problem. In particular, we consider to minimize:
subject to the constraints that and for .
We iteratively minimize subject to the constraints. Set . Suppose is obtained at the iteration. Fixing , we minimize over , subject to the constraint that (), to obtain . Once is obtained, fixing , we minimize over to obtain .
To minimize over fixing , it is equivalent to finding a landmark matching diffeomorphism , whose Beltrami coefficient closely resembles to and satisfies the landmark constraints . To obtain such , we propose to use the Linear Beltrami Solver (LBS) to find the descent direction for the Beltrami coefficients such that it approaches to and the corresponding satisfies the landmark constraints.
Let . From the Beltrami equation (1),
Let . We can write and as linear combinations of and ,
where ; ; .
Since , we obtain
In the discrete case, the elliptic PDEs (16) can be discretized into sparse positive definite linear systems. Given and the landmark constraints, one can solve the linear systems with the landmark constraints in the least square sense. A landmark matching quasi-conformal map , whose Beltrami coefficient closely resembles to , can then be obtained.
Once is obtained, we minimize over while fixing . In other words, we look for minimizing:
By considering the Euler-Lagrange equation, it is equivalent to solving:
In discrete case, equation (18) can be discretized into a sparse linear system and can be solved efficiently. However, by simply minimizing equation (18) does not ensure that the corresponding diffeomorphism of the resultant Beltrami coefficient satisfies the landmark constraints. Therefore, we use the Linear Beltrami Solver with input together with the landmark constraints to obtain a mapping with the corresponding Beltrami coefficients resemble to and satisfies the exact landmark matching requirement. Using as a descent direction, we update from the solution of the Equation (18) by for some small . This guarantees the resultant is smooth and the landmark mismatch decreases.
We keep the iteration going to obtain a sequence of pair . The iteration stops when for some small threshold . Theoretically, the conventional penalty method requires that increases in each iterations. In practice, we set to be a large enough constant and the algorithm gives satisfactory results.
In summary, the proposed landmark based registration model can be described as follows:
(Landmark based registration)
Images or surfaces: and ; cooresponding landmark sets and .
Optimal Beltrami coefficient and the landmark matching registration
Set . Use LBS to reconstruct from satisfying the landmark constraints.;
Given . Fixing , obtain by LBS satisfying the landmark constraints. Fixing , obtain by solving:
Use LBS to obtain from satisfying the landmark constraints. Obtain and update for some small .
If , continue. Otherwise, stop the iteration.
4.2.2 Hybrid registration model
The propsoed landmark based registration model can also be combined with matching intensity (such as image intensity for image registration or surface curvature for surface registration) to improve the accuracy of the registration result. More specifically, our goal is to look for an optimal Beltrami coefficient , which is the Beltrami coefficient of some diffeomorphism , minimizing the following energy functional :
subject to the constraints that and for . Here, and are the intensity functions defined on and respectively.
We again apply a splitting method to solve the above constrained optimization problem. We consider to minimize:
subject to the constraints that and is the quasi-conformal map with Beltrami coefficient satisfying for .
To solve the above optimization problem, we iteratively minimize subject to the constraints. Set . Suppose and is obtained at the iteration. Fixing , we minimize over , subject to the constraint that (), to obtain . Once is obtained, fixing , we minimize over to obtain .
We first discuss the minimization over fixing , subject to the constraint that (). This problem is equivalent to solving:
Using the gradient descent method, we compute the descent direction , which minimizes . is given by:
As is perturbed, its associated Beltrami coefficient is also adjusted by . The adjustment can be explicitly computed. Note that :
Thus, the adjustment can be obtained by:
Similiarly, we can obtain the descent direction , minimizes . is given by:
Therefore, the descent direction to solve the optimization problem (21) is given by:
Using the above formula for the descent direction, we obtain an updated Beltrami coefficient:
We then compute a quasi-conformal map , whose Beltrami coefficient closely resembles to , using LBS with the landmark constraints enforced. This step ensures a landmark matching registration can be obtained. We then update by: .
Once is obtained, fixing , we minimize over to obtain . In other words, we look for minimizing:
In the case when , by considering the Euler-Lagrange equation, it is equivalent to solving:
In discrete case, equation (30) can be discretized into a sparse linear system and can be solved efficiently. Similar to section 4.2.1, we use the Linear Beltrami Solver with input together with the landmark constraints to obtain a descent direction to update by for some small . This guarantees the resultant is smooth and the landmark mismatch decreases.
We keep the iteration going to obtain a sequence of pair . The iteration stops when for some small threshold . Again, the conventional penalty method requires that increases in each iterations. However, in practice, we set to be a large enough constant and the algorithm gives satisfactory results.
In summary, the proposed hybrid registration model can be described as follows:
Images or surfaces: and ; cooresponding landmark sets and ; intensity functions and defined on and respectively.
Optimal Beltrami coefficient and the landmark matching registration
Set . Use LBS to reconstruct from satisfying the landmark constraints. Set ;
Given and . Fixing , obtain by solving:
Using LBS, compute whose Beltrami coefficient closely resembles to with landmark constraints enforced. Let .
Fixing , obtain by solving:
Use LBS to obtain from satisfying the landmark constraints. Obtain and update for some small .
If , continue. Otherwise, stop the iteration.
5 Numerical implementation
The proposed models for landmark based and hybrid registration rely on the Linear Beltrami Solver(LBS) and solving the Euler-Lagrange(E-L) equations. In this section, we will describe the numerical implementation of the LBS and also the discretization of equations.
Practically speaking, 2D domains or surfaces in are usually represented discretely by triangular meshes. Suppose and are two surface meshes with the same topology representing and . We define the set of vertices on and by and respectively. Similarly, we define the set of triangular faces on and by and .
5.1 Numerical details of LBS
Suppose is an orientation preserving piecewise linear homeomorphism between and . We can assume and are both embedded in . In case and are surface meshes in , we first parameterize them conformally by and . The composition of with the conformal parameterizations, , is then an orientation preserving piecewise linear homeomorphism between and embedded in .
To compute the quasi-conformal mapping, the key idea is to discretize Equation 16 with two linear systems.
Given a map , we can easily compute its associated Beltrami coefficient , which is a complex-valued function defined on each triangular faces of . To compute , we simply need to approximate the partial derivatives on every face . We denote them by and respectively. Note that is piecewise linear. The restriction of on each triangular face can be written as:
Hence, , , and . Now, the gradient:
on each face can be computed by solving the linear system:
where: ; and