Abstract
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.
1. Introduction
Image registration is the process of finding out the global and/or local correspondence between two images, template and reference , 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.
Recently, new approaches have arisen for regularization terms in nonparametric image registration, for example, minimal curvature [2], the membrane energy [3], or the consistent symmetric approaches [4], along with efficient and stable implementations like using, for example, singular value decomposition [5], a DCT-type factorization [6], additive operator splitting schemes [7], or approaches which exploit the multiscale nature of the registration problem [2] or which take advantage of multigrid techniques [8]. In some applications, it could be interesting a registration scheme based on a regularization term with gradual stiffness properties, apart from the diffusion [7] and curvature [6] models, which are based, respectively, on firs-t and second-order derivatives. This paper proposes a hybrid regularization term based on fractional order derivatives, which can be seen as a generalization of the diffusion and curvature smoothing terms.
This paper is structured as follows: we start out with the variational framework for the registration problem, proposing the hybrid regularization term. In the following section, the energy functional is adapted to handle discrete -dimensional signal, then it is translated into the frequency domain, where it is minimized. Next, results show the performance of this regularizer and finally the conclusions close the paper.
2. Variational Formulation
The -dimensional reference and template datasets are defined , where , and represents the spatial dimension of the datasets. The registration produces a displacement field that makes the transformed template dataset similar to the reference dataset, , where and is the spatial position .
The nonparametric registration can be approached in terms of the variational calculus, by defining the joint energy functional to be minimized: The energy term measures the distance between the deformed template dataset and the reference dataset; is a penalty term which acts as a regularizer and determines the smoothness of the displacement field; and weights the influence of the regularization.
The distance measure 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 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 : with , where is the -dimensional Laplace operator. Note that if , then , and a diffusion registration is being performed; if , then , and a curvature registration is in this case performed. A precise definition of the operator can be found, for example, in [12], and references therein. The regularization term (2.2) can be also written as where is a positive bilinear form, denotes the dot (or inner) product in and is the following partial differential operator:
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 is the cardinal of the datasets to be registered, is the discrete spatial variable, , performs the discrete backward difference in the th dimension (i.e., , 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.
2.1. Translation into the Frequency Domain
This paper proposes a novel regularization term for -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, , where with being the frequency counterpart of the displacement field, is the -dimensional variable in the frequency domain, and where the distance measure and the regularization term are now defined in the frequency domain.
Taking into account that the Fourier transform of a Kronecker's delta shifted to , , is [13], and applying Parseval's theorem to (2.5), the regularization term is obtained in the frequency domain: where is a positive real constant and is the -dimensional frequency domain. Previous equation can be expressed in matrix form and as a bilinear form as follows: where is the complex inner product in , is a positive bilinear form in the frequency domain, and is a 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 is therefore a diagonal matrix which produces a displacement field whose components are decoupled. Then, can be written as where is the identity matrix, denotes the Kronecker product, and the th element of the diagonal of , the so-named , is given by
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 can actually be seen as a generalized regularization term which allows for a registration technique which is between the diffusion () and curvature () cases, because it simultaneously includes partial features of both schemes if .
2.2. Minimization in the Frequency Domain of the Energy Functional
According to the variational calculus, a necessary condition for a minimizer of the joint energy functional (2.6) is that the first variation of in any direction (also known as the Gâteaux derivative) vanishes for all suitable perturbations , that is,
For the Gâteaux derivative of , we find where the so-called force field in the frequency domain, , depends on the particular choice of the distance measure [16], .
For the Gâteaux derivative of , the following expression is obtained (detailed in the appendix): where (2.8) has been used. Note that any energy-based smoother 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: 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.
2.3. Frequency Implementation of the E-L Equations
To solve the Euler-Lagrange equations (2.16) formulated in the frequency domain, a time-marching scheme can be employed, yielding the following iteration: where (in the steady-state and (2.17) holds (2.16)). In order to solve (2.17), the time is discretized, being the time-step and being the iteration index, and then the time derivative of is replaced by its discrete approximation (first backward difference). Using the notation , the iterative scheme is the following: where denotes the identity on the domain , and where is usually initialized to zero, .
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 th component is given by where is a -dimensional low pass filter, it has its maximum at the frequencies of the DC component. The values of are less or equal than one and are the reciprocal of , therefore the matrix inversion necessary for solving (2.18) has become a pointwise division. is then the pseudoinverse filter of , which corresponds to a -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 .
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 , is low pass filtered. Figure 1 depicts the frequency spectra of filters for the diffusion, curvature and hybrid cases. Filter for curvature registration (Figure 1(c)) shows a wider passband and a narrower transition band than filter 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).
(a)
(b)
(c)
(d)
The practical implementation takes into account that datasets are discrete and then the spatial variable gives rise to the discrete spatial index , where with . In the same way, the frequency domain is also discretized and index , with and , corresponds with the frequencies where the spectra are evaluated. Then, the expression that implements the iterative scheme is the following: where , and . 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.
3. Results
In this section, the proposed regularization term is tested on a experiment involving three-dimensional () industrial images (obtained from the CT scan of two mechanical pieces [20]). The reference and template volumes are -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 dB, and , 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 .
(a)
(b)
Figure 3 shows the results obtained with the proposed method using , and a fractional . Table 1 summarizes the simulation parameters and (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 implies a perfect match of the datasets). An advantage of using 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.
(a)
(b)
(c)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
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 -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 background, this difference is hardly noticeable, as has been previously stated by other authors, see, for example, [5].
4. Conclusions
In this paper, a fractional regularization term is proposed for approaching the variational image registration problem. The joint energy functional is translated into the frequency 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 -FFT is taken into account) reduces considerably the numerical complexity and memory requirements of the overall iterative schemes, thus becoming efficient and fast variational registration techniques.
Appendix
Minimization of the Hybrid Regularization Term Defined in the Frequency Domain
The Gâteaux derivative of the proposed energy term is given by where it is taken into account that is real and symmetric, therefore, .
Acknowledgment
This work is partially supported by Ministerio de Educación y Ciencia, under Grant TEC2006-13338/TCM.