Fractional Regularization Term for Variational Image Registration

Image registration is a widely used task of image analysis with applications in many fields. Its classical formulation and current improvements are given in the spatial domain. In this paper a regularization term based on fractional order derivatives is formulated. This term is defined and implemented in the frequency domain by translating the energy functional into the frequency domain and obtaining the Euler-Lagrange equations which minimize it. The new regularization term leads to a simple formulation and design, being applicable to higher dimensions by using the corresponding multidimensional Fourier transform. The proposed regularization term allows for a real gradual transition from a diffusion registration to a curvature registration which is best suited to some applications and it is not possible in the spatial domain. Results with 3D actual images show the validity of this approach.


Introduction
Image registration is the process of finding out the global and/or local correspondence between two images, template T and reference R, of a scene in such a way that the transformed template matches 1 .This process is needed in various computer vision applications, such as stereo depth perception, motion analysis, remote sensing, change detection, object localization, object recognition, image fusion, and so forth.In these applications, a nonlinear transformation is necessary to correct the local differences between the images.This nonrigid registration is an ill-posed problem, therefore, it is necessary to constrain the estimated transformation as much as possible by using some prior knowledge on the deformation model.A regularization term is then used to preferentially obtain more likely solutions.

Variational Formulation
The d-dimensional reference and template datasets are defined R, T : Ψ ⊂ R d → R, where Ψ : 0, 1 d , and d represents the spatial dimension of the datasets.The registration produces a displacement field u : R d → R d that makes the transformed template dataset similar to the reference dataset, T x − u x ≈ R x , where u x u 1 x • • • u d x and x is the spatial position x x 1 , . . ., x d ∈ R d .The nonparametric registration can be approached in terms of the variational calculus, by defining the joint energy functional to be minimized: J u : D R, T ; u αS u .

2.1
The energy term D measures the distance between the deformed template dataset and the reference dataset; S is a penalty term which acts as a regularizer and determines the smoothness of the displacement field; and α > 0 weights the influence of the regularization.The distance measure D is chosen depending on the datasets to be registered.When the intensities of the two datasets are similar monomodal registration , the sum of squared differences of the datasets is commonly used 9 .When dealing with datasets from different sources or modalities multimodal registration , statistical based measures, for example, mutual information 10 or correlation ratio 11 , are more appropriate.
The regularization term S gives the smoothness characteristics to the displacement field 6 .In this paper, we propose the following term, which is given by the energy of σorder derivatives of u: ∇, and a diffusion registration is being performed; if σ 2, then Δ σ/2 Δ, and a curvature registration is in this case performed.A precise definition of the operator Δ σ/2 can be found, for example, in 12 , and references therein.The regularization term 2.2 can be also written as where a u, u ∈ R is a positive bilinear form, •, • R d denotes the dot or inner product in R d , and B is the following partial differential operator:

2.4
Because digital datasets are handled, which are typically encoded by uniformly distributed picture elements, the finite difference method is the natural approach to approximate 2.2 : where N N 1 N 2 . . .N d is the cardinal of the datasets to be registered, n n 1 , . . ., n d ∈ N d is the discrete spatial variable, ∀n • performs the discrete backward difference in the mth dimension i.e., d − n m δ n m − δ n m − 1 , δ n m being the Kronecker's delta for discrete signals 13 , and * denotes linear convolution of discrete signals.Note that the previous equation is meaningless in the spatial domain for noninteger values of σ.However, this problem can be approached from a frequency domain point of view, as it is detailed in following section.

Translation into the Frequency Domain
This paper proposes a novel regularization term for d-dimensional image registration, approached in terms of a variational formulation in the Fourier domain.As a starting point, the joint energy functional 2.1 has to be expressed in the frequency domain, following the procedure described in 14 .This translation is done by means of Parseval's theorem 15 , which relates the total energy of the signal in both spatial and frequency domain, that is, J u J u , where with u ω u 1 ω , . . ., u d ω being the frequency counterpart of the displacement field, ω ω 1 , . . ., ω d is the d-dimensional variable in the frequency domain, and where the distance measure D and the regularization term S are now defined in the frequency domain.
Taking into account that the Fourier transform of a Kronecker's delta shifted to n 0 , δ n m − n 0 , is e −jω m n 0 13 , and applying Parseval's theorem to 2.5 , the regularization term is obtained in the frequency domain: where frequency domain.Previous equation can be expressed in matrix form and as a bilinear form as follows: where r, s C d r s * is the complex inner product in C d , a u, u ∈ R is a positive bilinear form in the frequency domain, and A is a d × d matrix whose elements are scalar functions which implement the spatial derivatives in the frequency domain, allowing for their computation by means of products.The resulting frequency operator A ω is therefore a diagonal matrix which produces a displacement field whose components are decoupled.Then, A ω can be written as where I d is the d × d identity matrix, ⊗ denotes the Kronecker product, and the lth element of the diagonal of A ω , the so-named K ω , is given by

2.10
The regularization term 2.7 finally results in It should be noted that the regularization terms for diffusion and curvature registration are particular cases of 2.10 , and therefore the smoother S hybr can actually be seen as a generalized regularization term which allows for a registration technique which is between the diffusion σ 1 and curvature σ 2 cases, because it simultaneously includes partial features of both schemes if σ ∈ 1, 2 .

Minimization in the Frequency Domain of the Energy Functional
According to the variational calculus, a necessary condition for a minimizer u of the joint energy functional 2.6 is that the first variation of J u in any direction also known as the Gâteaux derivative vanishes for all suitable perturbations z, that is, For the Gâteaux derivative of D, we find where the so-called force field in the frequency domain, f ω , depends on the particular choice of the distance measure 16 , f ω d − FT{∇D R, T ; u } ∈ C d .For the Gâteaux derivative of S, the following expression is obtained detailed in the appendix : 2.14 where 2.8 has been used.Note that any energy-based smoother S can be incorporated into this framework, as long as it can be expressed in terms of 2.8 .Finally, we can write 2.12 as which leads to the Euler-Lagrange equation in the frequency domain:

2.16
Solving the previous equation in the frequency domain provides a stable implementation for the computation of a numerical solution for the displacement field, and in a more efficient way than existing approaches if the multidimensional fast Fourier transform is used.

Frequency Implementation of the E-L Equations
To solve the Euler-Lagrange equations 2.16 formulated in the frequency domain, a timemarching scheme can be employed, yielding the following iteration:  16 .In order to solve 2.17 , the time t is discretized, t : ξτ, τ > 0 being the timestep and ξ ∈ N being the iteration index, and then the time derivative of u ω, t is replaced by its discrete approximation first backward difference .Using the notation u ξ ω : u ω, ξτ , the iterative scheme is the following: where I denotes the identity on the domain Ω, and where u ξ ω is usually initialized to zero, u 0 ω : 0. As shown in 2.9 , the frequency components of the displacement field are independent i.e., they are not coupled .In this case, the matrix inversion in 2.18 does disappear due to the fact that the multiplication of a circulant block matrix and a column block vector becomes a Hadamard i.e., pointwise product of their respective spectra in the frequency domain 17 .Then, the iteration for the lth component is given by where

2.20
H ω is a d-dimensional low pass filter, it has its maximum at the frequencies of the DC component.The values of H ω are less or equal than one and are the reciprocal of 1 η −1 K ω , therefore the matrix inversion necessary for solving 2.18 has become a pointwise division.H ω is then the pseudoinverse filter of K ω , which corresponds to a d-dimensional high pass filter that contains the frequency representation of the spatial derivatives, and the constant η is related to the width of the transition band of filter H ω .The frequency point of view allows to understand the internal forces, with the restrictions imposed on the displacement field by the regularizer, as a low pass filtering.In 2.19 , each component of the displacement field as well as the driving external forces, weighted by the value η −1 , is low pass filtered.Figure 1 depicts the frequency spectra of filters H ω for the diffusion, curvature and hybrid cases.Filter H ω for curvature registration Figure 1 c shows a wider passband and a narrower transition band than filter H ω for diffusion registration Figure 1 a .As expected, filters obtained with the hybrid approach show an intermediate behavior between the diffusion and curvature cases, see Figures 1 b  and 1 d .
The practical implementation takes into account that datasets are discrete and then the spatial variable x gives rise to the discrete spatial index n n 1 , . . ., n d , where n l 0, . . ., N l − 1 with l 1, . . ., d.In the same way, the frequency domain is also discretized and index k k 1 , . . ., k d , with k q 0, . . ., N q − 1, and q 1, . . ., d, corresponds with the frequencies 2π/N 1 k 1 , . . ., 2π/N d k d where the spectra are evaluated.Then, the expression that implements the iterative scheme is the following: where H k 1 η −1 α K k −1 , and K k 2 d q 1 1 − cos 2π/N q k q σ .It should be noted that the proposed frequency domain implementation has the same complexity regardless the value of σ, that is, the same for diffusion, curvature, or intermediate registration scenarios.This implementation is, in terms of efficiency, two times faster than the fastest implementation available in the spatial domain 18 , which is the DCT-based algorithm included in the FLIRT toolbox 19 for the diffusion and curvature registration methods.

Results
In this section, the proposed regularization term is tested on a experiment involving threedimensional d 3 industrial images obtained from the CT scan of two mechanical pieces 20 .The reference and template volumes are 256 × 256 × 128-sized and are shown on the left and on the right, respectively, of Figure 2. To assess the validity of the proposed formulation, convenient quantitative measures are chosen: the peak signal-to-noise ratio PSNR and the correlation ratio CR .Before the registration process, the values for these measures are PSNR 15.82 dB, and CR 48.30%, respectively.The distance measure chosen for the registration is the sum of squared differences SSDs , the smoothing term is the proposed hybrid regularizer and the parameter η is set to a value of η 1.
Figure 3 shows the results obtained with the proposed method using σ 1, σ 2 and a fractional σ.Table 1 summarizes the simulation parameters α o and ξ o i.e., optimal regularization parameter and optimal number of iterations, resp., in terms of the  variational calculus, which are obtained as addressed in 21 , the PSNR, as well as CR, and the regularization energy for each registration scheme.In these simulations, the highest value of the similarity measure, the lowest value of the regularization energy, and the minimum iterations required correspond to the novel approach for the PSNR measure, an improvement of 0.5 dB becomes visible, and a value higher than 30 dB is usually considered a good match of the datasets, as shown in Figures 4 and 5; for the CR, a value of 100% implies a perfect match of the datasets .An advantage of using σ 1.75 is that in this way global rigid transformations are partially corrected because the regularizer is close to the curvature registration, moreover this value of σ avoids unlikely local curvatures because in this case the regularizer includes a slight component of diffusion registration.Note that the improvement in accuracy is not visually perceptible in most real cases as in this example , but it actually does exist, at least numerically see Table 1 .Additionally, the proposed registration scheme with noninteger values of σ can perform the optimal registration in a significantly lower number of iterations than border cases.
One point to take into account is about the boundary conditions considered, since spatial domain-based schemes impose Von Neumann boundary conditions, and the frequency domain-based scheme imposes periodic boundary conditions actually, due to the use of the d-dimensional discrete Fourier transform, periodic boundary conditions arise naturally when computing a numerical solution for the displacement field .Anyway, when dealing with images where the information is typically contained within a uniform domain by means of Parseval's theorem, and the minimization of the resulting variational equations is performed entirely in this domain, providing the Euler-Lagrange equation in the frequency domain with the internal forces for the hybrid regularizer.The proposed regularization term implemented in the frequency domain allows for a gradual transition between diffusion registration and curvature registration, thus providing better registration results in terms of both similarity of the images and smoothness of the transformation, and in a lower number of iterations of the algorithm.The use of the frequency domain especially if the d-FFT is taken into account reduces considerably the numerical

Figure 2 :
Figure 2: 3D view of two cylinders of an engine block: a reference dataset and b template dataset.

Figure 3 :
Figure 3: Registration results of two cylinders of an engine block: a diffusion σ 1 ; b curvature σ 2 ; c hybrid with σ 1.75.Table 1 summarizes the simulation parameters and quantifies numerically the results of registration.

Figure 5 :
Figure 5: Registration of industrial images using the proposed fractional regularization term.First column: slices of the reference dataset.Second column: slices of the template dataset.Third column: slices of the registered template σ 1.75 .First row: coronal slices slice no. 100 .Second row: axial slices slice no.64 .Third row: sagittal slices slice no.30 .

Table 1
summarizes the simulation parameters and quantifies numerically the results of registration.

Table 1 :
Parameters and numerical results for the registration of 3D industrial images.Each row corresponds, respectively, to registered datasets shown in Figure3.