3.1.1Current literature

Although the registration problem has been studied for almost two decades, there continue to be some fundamental limitations in the popular solutions that make them suboptimal, difficult to evaluate, and limited in scope. To explain these limitations let ℱ be a certain set of ℝn-valued functions on a domain D, made precise later. A pairwise registration between any two images f1, f2 ∈ ℱ is defined as finding a mapping γ, typically a diffeomorphism from D to itself, such that pixels f1(s) and f2(γ(s)) are optimallymatched to each other for all s ∈ D. To develop an algorithm for registration one needs: (1) An objective function for formalizing the notion of optimality, and (2) a numerical procedure for finding the optimal γ. Although the numerical techniques for optimization, that is item (2), have become quite mature over the last ten years, the commonly used objective functions themselves have several fundamental shortcomings. It is the choice of objective function that is the focus of this paper. The registration problems are commonly posed as variational problems, with the most common form of an objective function being

where |⋅| is the Euclidean norm, R is a regularization penalty on γ, typically involving its first and/ or second derivatives, λ is a positive constant, and Γ denotes the space of relevant deformations. Several variations of this functional are also used. We highlight the shortcomings of these methods using a broader discussion about desired properties of an objective function in the next section.

3.1.2Our approach

We propose a novel way to approach the image registration problem in this paper. Focusing on the deformations of image domain using diffeomorphisms of unimodal images, we introduce the concept of phase and amplitude of an image. This is motivated by progress in shape analysis of objects [1, 9, 11, 12, 14, 18, 22, 24] and phase-amplitude modes in functional data [16, 23, 25, 31]. We start by defining the amplitude of an image as the property that is unchanged under arbitrary diffeomorphic deformations. A deformation is said to change only the phase of an image. This description forms the basis for a formal development of a notion of phase and amplitude components of an image, and a proper computational framework for separating these components in the given image data. This separation effectively yields registration of the images involved. These developments rely on novel mathematical representations of images, termed q-maps, and the L2-norms between them. Alternatively, one can view this construction as based on a novel metric imposed on the image space, this metric is defined by pulling back the L2-metric under the q-map.

The main advantage of this viewpoint is a logical and interpretable separation between sources of variability in images. In cases where the images under study are primarily considered as deformations of each other, the main variability lies in their phase components and one can exploit this restriction to improve performance. In another situation where images are well registered already, the primary source of variability is amplitude and phase representation can be left out of the analysis. Most practical situations lie between these two extreme cases and, similar to analysis of variance in statistics, one is interested in decomposing the observed variance in these two components. Furthermore, one can impose separate metric structures on these components and study/ model them in an individual manner. In fact, we can interpret many of the past methods as imposing such metric structures without making them explicit.

The rest of this paper is organized as follows. In Section 3.2, we introduce the concepts of phase and amplitude of images, and describe the accompanying mathematical framework. In Section 3.3, we present a discussion on the theoretical properties of this framework, and its advantages relative to the current literature. We use Section 3.4 to outline a computational approach for separating phase-amplitude in a given dataset, and we present some experimental results in Section 3.5. We close the paper with a short summary and conclusion in Section 3.6.

3.2.1q-Map and amplitude distance

We start with this question: Suppose we want to compare the amplitudes of any two images, f1 and f2, using a proper metric, whatmetric can we use? If we try to use the L2-norm ‖f1 − f2‖ to quantify the difference in their amplitudes, we face a problem because ‖f1 − f2‖ ≠ ‖(f1, γ1) − (f2, γ2)‖ in general. In fact, it is not true even if γ1 = γ2. This is problematic because the action of Γ does not change the amplitudes of f1 and f2 and yet the norm is changing. Another way to express this behavior is that the action of Γ is not by isometries under the L2-metric. That is, for any two f1, f2 ∈ ℱ, and a general γ ∈ Γ, we have ‖f1 − f2‖ ≠ ‖f1 ∘ γ − f2 ∘ γ‖.

What are our other options? One potential metric can be derived as follows. This derivation requires a mathematical representation of images defined by a mapping called the q-map, that has been prompted by recent work in shape analysis of surfaces [13]. Here,we adopt it for analyzing images as follows. Let (x1, . . . , xm) : D → ℝm be coordinates on (a chart of) D and Jf(s) be the Jacobian matrix of f at s with the (j, i)-th element as ∂f j/∂xi(s). Define the generalized area multiplication factor of f at s for arbitrary n ≥ m as a(s) = ‖Jf(s)‖V where ‖Jf(s)‖V = |∂f/∂x1(s) ∧ ∂f/∂x2(s) ∧⋅ ⋅⋅ ∧ ∂f/ ∂xm(s)|. Here ∧ denotes the wedge product. The two special cases are: if m = n = 2, then a(s) = det Jf(s); if m = 2 and n = 3, then a(s) = ∂f/∂x1(s) × ∂f/∂x2(s). Given the area multiplication factor, we can define the q-map.

Definition 2.3. For an f ∈ ℱ, define a mapping Q: ℱ → L2 such that for any s ∈ D,

For any f ∈ ℱ, we will refer to q = Q(f) as its q-map; note that q : D → ℝm. Also, we remark that this is a general version of q-map for arbitrary ℝn that extends the work of [13]. If the original image is smooth, then the resulting q-map is square integrable. Intuitively, the q-map leaves homogeneous (constant) regions in an image as zeros while preserving edge information in such a way that it is compatible with change of variables, that is stronger edges get higher values. If we deform an image f to f ∘ γ, how does its q-map change? This is given by the corresponding action of Γ on L2 as follows.

Lemma 2.4. The right action of Γ on the space of q-maps, corresponding to the deformation of the original image, is given by the mapping L2 ×Γ → L2 as where Jγ(s) denotes the Jacobian matrix of γ at s.

Note that the mapping Q is equivariant, that is it can be shown that for an image f and any γ, Q((f, γ)) = (Q(f), γ), where the action on the right side is given by Lemma 2.4 and the action on the left is just the composition. Why is it useful to introduce a new mathematical representation of images. This question is answered by the following important property.

Proposition 2.5. The deformation group Γ acts on the L2-space by isometries under the L2-norm, that is ∀q1, q2 ∈ L2, ∀γ ∈ Γ, ‖(q1, γ) − (q2, γ)‖ = ‖q1 − q2‖.

Proof. Using a simple change of variables:

Setting q2 ≡ 0, ∀q1 ∈ L2 and γ ∈ Γ, we have ‖q1‖ = ‖(q1, γ)‖. Thus, warping of images is actually a unitary operator under this representation.

We return to the problem of comparing amplitudes of any two given images and use q-maps to define a new metric.

Definition 2.6 (Amplitude Distance). Define the distance between the amplitude components in any two images f1 and f2 as the minimum L2-norm between the two orbits in terms the q-maps:

In other words, find the two elements in their respective orbits such that their q-maps are as close to each other as possible with respect to the L2-distance. This distance is denoted by the letter B in Figure 1. Some properties of this distance are as follows:

Not only is the amplitude distance well defined, it agrees with the desired goals of this definition. It compares amplitude components of images and is invariant to arbitrary deformations of those images. Beyond comparing the amplitude components, the quantity da is of central importance in registration of images and in generalizing the notion of relative phase to arbitrary pairs of images. This is described in the next section.

3.2.2Relative phase and image registration

The introduction of q-map sets up the general concept of relative phases between images. As mentioned earlier, the notion of phases relates to the amount of deformation in an image. However, the specification of deformation requires two images – one before and one after deformation. In fact, as shown below, one can talk about the amount of deformation between any two images. We start with the registration problem.

Definition 2.7 (Registration of Images). The registration between any two images f1 and f2, represented by their q-maps, q1 and q2, is defined as the solution to the optimization problem:

This makes γ12 the optimal mapping from the domain of f1 to the domain of f2. With this definition, the pixel f1(s) is said to be registered to the pixel f2(γ12(s)) for any s ∈ D. This also leads to the definition of relative phase.

Definition 2.8 (Relative Phase between Images). The relative phase of an image f1, with respect to another image f2, is given by γ12 defined using the equation 2.4.

It is interesting to note that while the notion of amplitude is independent of any metric structure on the image space, the notion of relative phase is closely tied to the choice of the amplitude distance. That is da([f1], [f2]) = ‖q1 − (q2, γ12)‖. The actual value of relative phase may change substantially if another metric structure is used instead.

There are some interesting properties associated with this definition of relative phase:

for any c > 0.

3.4.1Basis on

In order to compute the gradient ∇E and to update γo we need to specify an orthonormal basis for . Since Γ = Diff+([0, 1]2), the set of all boundary preserving diffeomorphisms on [0, 1]2, the tangent space is defined as the set of smooth vector fields on [0, 1]2:

We start from constructing an orthonormal basis on L2([0, 1], ℝ) and then extend it to the 2D case. It is known that

forms an orthonormal basis for L2([0, 1], ℝ) under the L2-metric. We seek an orthonormal basis for L2([0, 1], ℝ) under the Palais metric due to some nice properties of this Riemannian metric [21], to get a basis set that vanishes at the boundaries. The 1D Palais metric is defined as for f, g ∈ L2([0, 1], ℝ). Under this metric, the set

provides an orthonormal basis of functions under L2 that vanish at t ∈ {0, 1}.

We will use the Cartesian product of b (with elements b) and b ̃ (with elements b̃) to construct an orthonormal basis for functions on [0, 1]2 under the 2D Palais metric given by the inner product:

By enumerating the 1D basis as and the 2D basis elements are constructed as The basis elements can be orthonormalized using the Gram–Schmidt process under the metric defined in equation 4.3.

3.4.2Mean image and group-wise registration

An important problem in image analysis, especially medical image analysis, is to compute a “typical” or an “average” of several images from the same class and use it as a template. Then, the individual images can be registered to the sample mean in a pairwise manner, resulting in a group registration. By registering member images to the group mean, one can analyze their variations from the typical template image. Suppose there is a set of n images, Their Karcher mean is defined as the image that minimizes the sum of squares of the distances to the given images, that is Since this quantity is an amplitude, it results in a set of images that have the same amplitude. In order to pick a central amplitude from this set, we use the notion of center of an orbit, as defined in [25] for real-valued functions on [0, 1]. It lays out a procedure to find an instance of mean image fμ, such that fμ ∈ [μ]. With this image, we can register groups of images to the mean image rather than to an arbitrarily chosen template, as is often done in current methods.

3.5.1Pairwise image registration

We first present some results on synthetic images to demonstrate the use of the registration framework introduced in equation 2.4. When given a pair of synthetic images, f1 and f2, we register them twice. First take f1 as the template image and estimate γ12 that optimally deforms f2, that is γ12 = arg minγ∈Γ ℒ(f1, (f2, γ)). Then, the roles are reversed and f2 is used as the template to obtain γ21, that is γ21 = arg minγ∈Γ ℒ (f2, (f1, γ)).

Figure 3 shows a synthetic example to illustrate the property of deformations from our algorithm. From left to right, we show images f1 and f2, followed by optimal deformations for forward and backward registration.

In Figure 4, we use a different set of images and show the original images f1 and f2 before and after warping. The diffeomorphisms γ12 and γ21, used to register the images, are also presented. By composing them in different orders we expect the resulting diffeomorphisms to be the identity map. We show the two converged objective functions, ‖(q1, γ21) − q2‖ and ‖q1 − (q2, γ12)‖, associated with the optimal γ12 and γ21 to verify symmetry. The cumulative diffeomorphisms γ21 ∘γ12 and γ12 ∘γ21 are also used to demonstrate the inverse consistency of the proposed metric. As mentioned above, the theory shows that γ12 and γ21 are expected to be inverses of each other. In order to better visualize that the composed diffeomorphisms are close to identity, we apply them to checkerboard images and we observe that the composed diffeomorphisms γ21 ∘ γ12 and γ12 ∘ γ21 are close to the identity map. Figure 5 shows similar results on registering a pair of binary images.

Next, we add simulated Gaussian noise to images. First, the original images are registered, as described previously, and then the corresponding noisy image pairs are registered again, as shown in Figure 6. When the registration results and the diffeomorphisms are compared, we find that the results are quite similar to each other. This result underlines the fact that although the image representation using q-maps involves taking a gradient, the resulting registration is relatively immune to the presence of high-frequency noise.

In Figure 7, we present registration results on 2D brain MR images. In order to illustrate our method, in each of the two experiments, we show (1) images, f1 and f2, (2) f1 and f2 ∘ γ21, and (3) f2 and f1 ∘ γ12. In each panel, the two concerned images are shown overlapped on a common canvas such that red and green denote positive and negative image differences, respectively. We also take brain MRIs from two modalities collected simultaneously (in this case n = 2 and our method directly applies) and register both modalities jointly using a common diffeomorphism. The results are shown in Figure 8. In both cases, the deformations are relatively small and the identity map provides a good initial condition for the gradient method.

Landmark-aided registration

Our framework can be extended to incorporate landmark information during registration and all of the nice mathematical properties of the objective function are preserved. Assume that there are a fixed number, say K, of distinct landmark points, ? = {p1, p2, . . . , p K}, in the image domain D. They are typically chosen according to the application but fixed within the analysis. The landmark-guided registration is

achieved by defining a subgroup of Γ, denoted by ΓP, which preserves the landmark points:

Given two images f1, f2 with landmark information ?, the images can be registered in two steps in an iterative way:

1. Register the landmark points ? and apply an initial deformation to f1 to form such that the landmarks are at the same locations in and f2.
2. Register to f2 using equation 2.4 restricting to the subgroup ΓP.

A similar technique of forming a landmark-constrained basis on ?2 has been used on closed surfaces as described in [14]. In the second step, searching over ΓP ensures correspondences of landmarks established in step 1 are preserved. The registration is refined without moving the landmarks. This search is based on a basis BP in Tγid (ΓP) constructed such that its elements, the vector fields on D, vanish at the landmark locations {pi}. We remark that ΓP forms a subgroup of Γ and, as a result, the desired properties discussed earlier are still satisfied. This approach may be termed landmarkguided registration where the landmarks are treated as hard constraints.

We show results on MRIs with landmarks in Figure 9. In the first example, presented in the first row of Figure 9, the optimally deformed f1 searched over the full group Γ is displayed at the end as (f1, γ). The deformation in the skull is so large that our original method fails to reach a global minimizer of the objective function. By adopting the landmark-guided registration, we at first get a deformed image with nicely matched landmarks and the skull deformed correspondingly. Then, is further deformed to register the intensity details without moving the landmarks achieved by optimizing over the restricted subgroup ΓP. The final result matches f2 with no artifacts near the skull. Two more examples are shown in the bottom two rows of Figure 9. Generally, the registration with landmarks outperforms the identity map as the initial condition of our procedure, especially when the deformations are large.

3.5.2Registering multiple images

We use part of the MNIST database of handwritten digits to illustrate our method to compute the mean image and multiple image registration. Images containing a certain digit all have the same pattern, or the same amplitude. Therefore, within the same group, the phase variation can be left out of analysis by groupwise registration.

In Figure 10 we present the mean image of each digit computed without and with registration, respectively. Two examples of registering images within digit groups are shown in Figures 11 and 12.

We also show an example using brain MRIs in Figure 13. Four brain images without alignment are shown in Figure 13a with the corresponding mean image. This mean without registration appears blurred due to misalignment. In Figure 13b, the images are aligned to the Karcher mean as described in Section 3.4.2. We can see that with multiple registration the mean image improves the blurriness issue.

3.5.3Image classification

The amplitude distance introduced in equation 2.3 is useful in many applications relating to pattern analysis of images such as clustering or classification. For the task of handwritten digit recognition, amplitude variation is of main interest: images containing the same digit are from the same orbit thus are supposed to have zero amplitude distance (or very small amplitude distance in practice); conversely, images containing distinct digits have different amplitudes so that their amplitude distance should be large. Distance-based classifiers can then be applied for recognition after the pairwise distances between images are computed. To illustrate this idea, we use a subset of the MNIST database of images of handwritten digits from 0 to 9. It contains ten images of handwriting for each digit. In addition to the baseline L2-distance, or sum of squared distances (SSD) (without any warping), we compare our method to three other methods – diffeomorphic demons [35], FAIR [20], and NiftyReg [19].

For computing distance matrices between all pairs of images, the digits are registered in a pairwise manner using each of the three methods and then the SSD is computed for the registered images as a measure of distance. In the case of our method the amplitude distance defined in equation 2.3 is used. Using the leave-one-out nearestneighbor (LOO-NN) classifier, the rates of correct classification are listed in Table 2. Provided as a baseline, the L2-distance without any registration provides a rate of classification of 76% with our method performing the best with 94%. Our method has misclassified three pairs of images. One misclassified pair in all methods is presented in Figure 14.

Table 2: Classification of MNIST Digits.