# Fractional Regularization Term for Variational Image Registration

Authors: Rafael Verdú-Monedero; Jorge Larrey-Ruiz; Juan Morales-Sánchez; José Luis Sancho-Gómez
Journal: Mathematical Problems in Engineering
 (2009)

---

## 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.


---

## Body

## 1. Introduction

Image registration is the process of finding out the global and/or local correspondence between two images, templateT 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.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 discreted-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

Thed-dimensional reference and template datasets are defined R,T:Ψ⊂ℝd→ℝ, where Ψ:=]0,1[d, and d represents the spatial dimension of the datasets. The registration produces a displacement field u:ℝd→ℝd that makes the transformed template dataset similar to the reference dataset, T(x-u(x))≈R(x), where u(x)=(u1(x)⋯ud(x))⊤ and x is the spatial position x=(x1,…,xd)∈ℝd.The nonparametric registration can be approached in terms of the variational calculus, by defining the joint energy functional to be minimized:(2.1)𝒥[u]:=𝒟[R,T;u]+α𝒮[u].
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 α>0 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 u:
(2.2)𝒮hybr[u]:=12∑l=1d∫Ψ∥Δσ/2ul∥2dx,
with σ∈[1,2]⊂ℝ, where Δ≡∇2≡∇⊤∇ is the d-dimensional Laplace operator. Note that if σ=1, then Δσ/2=∇, 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
(2.3)𝒮hybr[u]:=12a[u,u]=12∫Ψ〈ℬ[u],ℬ[u]〉ℝddx,
where a[u,u]∈ℝ is a positive bilinear form, 〈·,·〉ℝd denotes the dot (or inner) product in ℝd, and ℬ is the following partial differential operator:
(2.4)ℬ[u]:=(∂σu1∂x1σ+⋯+∂σu1∂xdσ,…,∂σud∂x1σ+⋯+∂σud∂xdσ)⊤.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):
(2.5)𝒮hybr[u]≈12∑l=1d1N∑∀n|∑m=1dd-[nm]*⋯*d-[nm]︷σ-times*ul|2,
where N=N1N2…Nd is the cardinal of the datasets to be registered, n=(n1,…,nd)∈ℕd is the discrete spatial variable, ∑∀n(·)=∑n1=0N1-1⋯∑nd=0Nd-1(·), d-[nm] performs the discrete backward difference in the mth dimension (i.e., d-[nm]=δ[nm]-δ[nm-1], δ[nm] 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 ford-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, 𝒥[u]=𝒥̃[ũ], where
(2.6)𝒥̃[ũ]:=𝒟̃[R̃,T̃;ũ]+α𝒮̃[ũ],
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 𝒟̃ and the regularization term 𝒮̃ are now defined in the frequency domain.Taking into account that the Fourier transform of a Kronecker's delta shifted ton0,  δ[nm-n0], is e-jωmn0 [13], and applying Parseval's theorem to (2.5), the regularization term is obtained in the frequency domain:
(2.7)𝒮̃hybr[ũ]=ν2∑l=1d∫Ω|∑m=1d|1-e-jωm|2|σ|ũl|2dω,
where ν=((2π)dN)-1 is a positive real constant and Ω:=]-π,π[d is the d-dimensional frequency domain. Previous equation can be expressed in matrix form and as a bilinear form as follows:
(2.8)𝒮̃hybr[ũ]=ν2∫Ω∑l=1d∑m=1d𝒜̃lm(ω)ũmũl*dω=ν2∫Ω(ũ1⋯ũd)(𝒜̃11(ω)⋯𝒜̃d1(ω)⋮⋱⋮𝒜̃1d(ω)⋯𝒜̃dd(ω))(ũ1*⋮ũd*)dω=ν2∫Ωũ⊤𝒜̃⊤(ω)ũ*dω=ν2∫Ω〈𝒜̃(ω)ũ,ũ〉ℂddω=12ã[ũ,ũ],
where 〈r̃,s̃〉ℂd=r̃⊤s̃* is the complex inner product in ℂd,  ã[ũ,ũ]∈ℝ is a positive bilinear form in the frequency domain, and 𝒜̃ 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𝒜̃(ω) is therefore a diagonal matrix which produces a displacement field whose components are decoupled. Then, 𝒜̃(ω) can be written as
(2.9)𝒜̃(ω)=Id⊗𝒦̃(ω),
where Id is the d×d identity matrix, ⊗ denotes the Kronecker product, and the lth element of the diagonal of 𝒜̃(ω), the so-named 𝒦̃(ω), is given by
(2.10)𝒦̃(ω)=(2∑m=1d(1-cosωm))σ.The regularization term (2.7) finally results in
(2.11)𝒮̃hybr[ũ]=ν2∫Ω𝒦̃(ω)∥ũ∥2dω.
It should be noted that the regularization terms for diffusion and curvature registration are particular cases of (2.10), and therefore the smoother 𝒮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[.


### 2.2. Minimization in the Frequency Domain of the Energy Functional

According to the variational calculus, a necessary condition for a minimizerũ of the joint energy functional (2.6) is that the first variation of 𝒥̃[ũ] in any direction (also known as the Gâteaux derivative) vanishes for all suitable perturbations z̃, that is,
(2.12)δ𝒥̃[ũ;z̃]=δ𝒟̃[R̃,T̃;ũ;z̃]+αδ𝒮̃[ũ;z̃]=0,∀z̃∈ℂd.For the Gâteaux derivative of𝒟̃, we find
(2.13)δ𝒟̃[R̃,T̃;ũ;z̃]=limϵ→01ϵ(𝒟̃[R̃,T̃;ũ+ϵz̃]-𝒟̃[R̃,T̃;ũ])=ν∫Ω〈f̃(ω),z̃〉ℂddω,
where the so-called force field in the frequency domain, f̃(ω), depends on the particular choice of the distance measure [16], f̃(ω)=d-ℱ𝒯{∇𝒟[R,T;u]}∈ℂd.For the Gâteaux derivative of𝒮̃, the following expression is obtained (detailed in the appendix):
(2.14)δ𝒮̃[ũ;z̃]=limϵ→012ϵ(ã[ũ+ϵz̃,ũ+ϵz̃]-ã[ũ,ũ])=ã[ũ,z̃]=ν∫Ω〈𝒜̃(ω)ũ,z̃〉ℂddω,
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
(2.15)δ𝒥̃[ũ;z̃]=ν∫Ω〈f̃(ω)+α𝒜̃(ω)ũ,z̃〉ℂddω=0,∀z̃∈ℂd,
which leads to the Euler-Lagrange equation in the frequency domain:
(2.16)f̃(ω)+α𝒜̃(ω)ũ(ω)=0.
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:
(2.17)∂tũ(ω,t)+f̃(ω,t)+α𝒜̃(ω)ũ(ω,t)=0,
where ∂tũ(ω,t)=(∂tũ1(ω,t),…,∂tũd(ω,t))⊤ (in the steady-state ∂tũ(ω,t)=0 and (2.17) holds (2.16)). In order to solve (2.17), the time t is discretized, t:=ξτ,τ>0 being the time-step and ξ∈ℕ 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:
(2.18)ũ(ξ)(ω)=(I+τα𝒜̃(ω))-1(ũ(ξ-1)(ω)-τf̃(ξ-1)(ω)),
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
(2.19)ũl(ξ)(ω)=H(ω)(ũl(ξ-1)(ω)-η-1f̃l(ξ-1)(ω)),
where
(2.20)η:=1τ,K(ω):=α𝒦̃(ω),H(ω):=ηη+K(ω).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+η-1K(ω), 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).Spectra ofd-dimensional filters, d=2,  η=1,  α=1. (a) H(ω) for diffusion (σ=1), (b) H(ω) for the hybrid case (σ=1.5), (c) H(ω) for curvature (σ=2), (d) comparison of the frequency responses.
(a)(b)(c)(d)The practical implementation takes into account that datasets are discrete and then the spatial variablex gives rise to the discrete spatial index n=(n1,…,nd), where nl=0,…,Nl-1 with l=1,…,d. In the same way, the frequency domain is also discretized and index k=(k1,…,kd), with kq=0,…,Nq-1, and q=1,…,d, corresponds with the frequencies ((2π/N1)k1,…,(2π/Nd)kd) where the spectra are evaluated. Then, the expression that implements the iterative scheme is the following:
(2.21)ul(ξ)[n]=IFFT{H[k]·FFT{ul(ξ-1)[n]-η-1fl(ξ-1)[n]}},
where H[k]=(1+η-1α𝒦̃[k])-1, and 𝒦̃[k]=(2∑q=1d(1-cos((2π/Nq)kq)))σ. 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.


## 2.1. Translation into the Frequency Domain

This paper proposes a novel regularization term ford-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, 𝒥[u]=𝒥̃[ũ], where
(2.6)𝒥̃[ũ]:=𝒟̃[R̃,T̃;ũ]+α𝒮̃[ũ],
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 𝒟̃ and the regularization term 𝒮̃ are now defined in the frequency domain.Taking into account that the Fourier transform of a Kronecker's delta shifted ton0,  δ[nm-n0], is e-jωmn0 [13], and applying Parseval's theorem to (2.5), the regularization term is obtained in the frequency domain:
(2.7)𝒮̃hybr[ũ]=ν2∑l=1d∫Ω|∑m=1d|1-e-jωm|2|σ|ũl|2dω,
where ν=((2π)dN)-1 is a positive real constant and Ω:=]-π,π[d is the d-dimensional frequency domain. Previous equation can be expressed in matrix form and as a bilinear form as follows:
(2.8)𝒮̃hybr[ũ]=ν2∫Ω∑l=1d∑m=1d𝒜̃lm(ω)ũmũl*dω=ν2∫Ω(ũ1⋯ũd)(𝒜̃11(ω)⋯𝒜̃d1(ω)⋮⋱⋮𝒜̃1d(ω)⋯𝒜̃dd(ω))(ũ1*⋮ũd*)dω=ν2∫Ωũ⊤𝒜̃⊤(ω)ũ*dω=ν2∫Ω〈𝒜̃(ω)ũ,ũ〉ℂddω=12ã[ũ,ũ],
where 〈r̃,s̃〉ℂd=r̃⊤s̃* is the complex inner product in ℂd,  ã[ũ,ũ]∈ℝ is a positive bilinear form in the frequency domain, and 𝒜̃ 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𝒜̃(ω) is therefore a diagonal matrix which produces a displacement field whose components are decoupled. Then, 𝒜̃(ω) can be written as
(2.9)𝒜̃(ω)=Id⊗𝒦̃(ω),
where Id is the d×d identity matrix, ⊗ denotes the Kronecker product, and the lth element of the diagonal of 𝒜̃(ω), the so-named 𝒦̃(ω), is given by
(2.10)𝒦̃(ω)=(2∑m=1d(1-cosωm))σ.The regularization term (2.7) finally results in
(2.11)𝒮̃hybr[ũ]=ν2∫Ω𝒦̃(ω)∥ũ∥2dω.
It should be noted that the regularization terms for diffusion and curvature registration are particular cases of (2.10), and therefore the smoother 𝒮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[.


## 2.2. Minimization in the Frequency Domain of the Energy Functional

According to the variational calculus, a necessary condition for a minimizerũ of the joint energy functional (2.6) is that the first variation of 𝒥̃[ũ] in any direction (also known as the Gâteaux derivative) vanishes for all suitable perturbations z̃, that is,
(2.12)δ𝒥̃[ũ;z̃]=δ𝒟̃[R̃,T̃;ũ;z̃]+αδ𝒮̃[ũ;z̃]=0,∀z̃∈ℂd.For the Gâteaux derivative of𝒟̃, we find
(2.13)δ𝒟̃[R̃,T̃;ũ;z̃]=limϵ→01ϵ(𝒟̃[R̃,T̃;ũ+ϵz̃]-𝒟̃[R̃,T̃;ũ])=ν∫Ω〈f̃(ω),z̃〉ℂddω,
where the so-called force field in the frequency domain, f̃(ω), depends on the particular choice of the distance measure [16], f̃(ω)=d-ℱ𝒯{∇𝒟[R,T;u]}∈ℂd.For the Gâteaux derivative of𝒮̃, the following expression is obtained (detailed in the appendix):
(2.14)δ𝒮̃[ũ;z̃]=limϵ→012ϵ(ã[ũ+ϵz̃,ũ+ϵz̃]-ã[ũ,ũ])=ã[ũ,z̃]=ν∫Ω〈𝒜̃(ω)ũ,z̃〉ℂddω,
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
(2.15)δ𝒥̃[ũ;z̃]=ν∫Ω〈f̃(ω)+α𝒜̃(ω)ũ,z̃〉ℂddω=0,∀z̃∈ℂd,
which leads to the Euler-Lagrange equation in the frequency domain:
(2.16)f̃(ω)+α𝒜̃(ω)ũ(ω)=0.
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:
(2.17)∂tũ(ω,t)+f̃(ω,t)+α𝒜̃(ω)ũ(ω,t)=0,
where ∂tũ(ω,t)=(∂tũ1(ω,t),…,∂tũd(ω,t))⊤ (in the steady-state ∂tũ(ω,t)=0 and (2.17) holds (2.16)). In order to solve (2.17), the time t is discretized, t:=ξτ,τ>0 being the time-step and ξ∈ℕ 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:
(2.18)ũ(ξ)(ω)=(I+τα𝒜̃(ω))-1(ũ(ξ-1)(ω)-τf̃(ξ-1)(ω)),
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
(2.19)ũl(ξ)(ω)=H(ω)(ũl(ξ-1)(ω)-η-1f̃l(ξ-1)(ω)),
where
(2.20)η:=1τ,K(ω):=α𝒦̃(ω),H(ω):=ηη+K(ω).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+η-1K(ω), 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).Spectra ofd-dimensional filters, d=2,  η=1,  α=1. (a) H(ω) for diffusion (σ=1), (b) H(ω) for the hybrid case (σ=1.5), (c) H(ω) for curvature (σ=2), (d) comparison of the frequency responses.
(a)(b)(c)(d)The practical implementation takes into account that datasets are discrete and then the spatial variablex gives rise to the discrete spatial index n=(n1,…,nd), where nl=0,…,Nl-1 with l=1,…,d. In the same way, the frequency domain is also discretized and index k=(k1,…,kd), with kq=0,…,Nq-1, and q=1,…,d, corresponds with the frequencies ((2π/N1)k1,…,(2π/Nd)kd) where the spectra are evaluated. Then, the expression that implements the iterative scheme is the following:
(2.21)ul(ξ)[n]=IFFT{H[k]·FFT{ul(ξ-1)[n]-η-1fl(ξ-1)[n]}},
where H[k]=(1+η-1α𝒦̃[k])-1, and 𝒦̃[k]=(2∑q=1d(1-cos((2π/Nq)kq)))σ. 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 (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.3D view of two cylinders of an engine block: (a) reference dataset and (b) template dataset.
(a)(b)Figure3 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.Table 1
Parameters and numerical results for the registration of 3D industrial images. Each row corresponds, respectively, to registered datasets shown in Figure3.
RegularizerαoξoPSNRCRα𝒮̃Diffusion (σ=1)50199539.43 dB99.73%0.185Curvature (σ=2)65195041.55 dB99.85%0.162Fractional (σ=1.75)40133541.64 dB99.86%0.159Registration 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.
(a)(b)(c)Registration of industrial images using the proposed fractional regularization term. First column: views of the reference dataset. Second column: views of the template dataset. Third column: views of the registered template (σ=1.75). First row: front views. Second row: top views. Third row: side views.
(a)(b)(c)(d)(e)(f)(g)(h)(i)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).
(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 thed-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 thed-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.



---

*Source: 707026-2009-10-08.xml*