# A Combination of the Total Variation Filter and a Fourth-Order Filter for Image Registration

Authors: Xiaomei Yang; Chaomin Shen; Fang Li; Chunli Shen
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

We introduce a novel method for nonrigid image registration which combines the total variation filter and a fourth-order filter. We decompose the deformation field into two components, that is, a piecewise constant component and a smooth component. The total variation filter is used for the first component and the fourth-order filter is used for the second one. Then, we present a new PDE-based image registration model suitable for both smooth and nonsmooth deformation problem. Meanwhile, the local-global similarity measure is used in our method to improve the accuracy and robustness for image matching. By applying the split Bregman algorithm and dual algorithm, we present a fast and stable numerical scheme. The numerical experiments and comparisons on both synthetic images and real images demonstrate the effectiveness of our method in nonrigid image registration.


---

## Body

## 1. Introduction

Image registration is playing an important role in image analysis, and having application in various fields (e.g., image fusion, atlas matching, and pathological diagnosis). The purpose of image registration is to find an optimal geometric transformation that aligns points in one view of an object with corresponding points in another view of the same object or a similar one. The transformation can be either rigid or nonrigid deformation. Particularly, nonrigid image registration method is a challenging subject in today’s modern medical diagnostics and image-guided therapy systems.Lots of research has been devoted to nonrigid image registration in the last couple of years. We focus on variational and PDE-based image registration methods in this paper, which have been proven to be very successful techniques in recent years. Bajcsy and Kovačič [1] utilized the Navier-Stoke equation to represent the local deformation and Broit [2] proposed the elastic model for nonrigid registration. Thirion proposed the Demon’s algorithm by considering image registration as a diffusion process [3]. Improvements of Demon’s algorithm for nonrigid registration was studied by Cachier et al. [4]. Cachier and Rey [5] introduced a method to symmetrize the registration problem using an inversion-invariant energy.In many variational methods, the energy includes two terms. One is similarity measure term which describes the similarity between the deformed source image and the target image. The other is regularization term which is related to the smoothness of deformation field. The choices of the similarity measures depend on the particular images to be registered. When the intensities of the images are similar, the sum of square distance (SSD) [6, 7] or the sum of absolute differences (SAD) [8, 9] of images is commonly used. When the images come from different sources or modalities, the cross correlation (CC) [10, 11], mutual information (MI) [11–13], or other information-theoretic measures [14, 15] are considered. For the mono-modality registration, one of the most common and simple similarity measure is SSD measure [4]. A local-global similarity measure is introduced in [16] to enhance the matching accuracy for very large displacements. The classical regularization term uses the quadratic term of gradient [17, 18] which leads to smooth estimation of deformation field and does not allow for discontinuities. The second-order total variation (TV) regularization is considered in [19] which is better for preserving discontinuities of the deformation fields. However, it has been shown that TV regularization leads to staircase effect in image restoration problems since it favors piecewise constant solutions which seriously decreases the visual image quality [20]. To reduce the staircase effect of the second-order TV regularization, fourth-order regularization methods are introduced in [21–23]. It is reported that fourth-order diffusion damps oscillations much faster than second-order diffusion and is more suitable for restoring smooth contents.In fact, for the nonrigid registration problems, the deformation field have both discontinuities and smooth structures as images [24]. So it is necessary to design a deformation model which is suitable for deformation field with both smooth and nonsmooth structures. In this paper, we propose a novel registration model, which combines the total variation regularization with the fourth-order regularization. The combined algorithm takes the advantage of both regularization methods and overcomes their demerits. Additionally, we draw support from the local-global similarity measure to enhance the matching accuracy for very large displacements. Finally, we adopt an iterative reweighted minimization scheme to solve our proposed model. All numerical experiments demonstrate the efficiency and stability of our proposed model.The remainder of the paper is organized as follows: in Section2, we give some preliminaries. In Section 3, we propose our model and algorithm for image registration. In Section 4, we show experimental results on both synthetic images and real images. We also compare our method with some closely related methods. Finally, we give concluding remarks in Section 5.


## 2. Preliminaries

A general framework for nonrigid image registration can be formulated in the following. We assume thatΩ
⊂
R
2 is a bounded domain with Lipschitz boundary and satisfying the cone condition. Let us denote the reference image by R
:
Ω
→
R and the template image by T
:
Ω
→
R. In nonrigid image registration, the deformation field is usually described by displacement field u
:
Ω
→
Ω denoted by (1)
u
:
x
⟶
u
x
=
u
1
x
,
u
2
x
T
.The problem is seeking for an deformation field u such that the transformed template image is optimally correlated with the reference image R
(
x
); that is, (2)
T
x
+
u
x
≈
R
x
,
x
∈
Ω
.Meanwhile, in order to avoid the ill conditioned problem [25] in the sense of Hadamard, it is necessary to impose an appropriate regularizer, which will penalize the unstable and nonsmooth solutions. In variational methods, the general model for nonrigid image registration is formulated as the following minimization problem: (3)
m
i
n
u
E
u
=
F
u
+
λ
R
u
,where F
(
u
) is the similarity measure between the transformed image and the reference image, R
(
u
) is the regularization term, and λ
>
0 is a weight parameter to balance the influence of F
(
u
) and R
(
u
).In the existing works related to image registration in Section1, the frequently used similarity measure is SSD measure [4] (4)
F
S
S
D
u
=
1
2
∫
Ω
T
x
+
u
x
-
R
x
2
d
x
.Meawhile, the widely used regularization term is quadratic L
2 regularization term [17, 18] (5)
R
u
=
∑
i
=
1
d
∫
Ω
∇
u
i
2
d
xand TV regularization term [19] (6)
R
u
=
∑
i
=
1
d
∫
Ω
∇
u
i
d
x
.In image restoration problems, the following fourth-order regularization term is introduced: (7)
R
u
=
∑
i
=
1
d
∫
Ω
∇
2
u
i
d
x
.


## 3. The Proposed Model and Algorithm

We assume that the deformation field, which transform the template imageT to the reference image R, is decomposed into two parts u
=
u
1
+
u
2, where u
1 is the piecewise constant part and u
2 is the smooth part. We denote(8)
u
i
x
=
u
i
1
x
,
u
i
2
x
T
,
i
=
1,2
.Firstly, we present the local-global similarity measure based on the bilateral filter [26]: (9)
F
ω
u
=
∫
Ω
ω
∘
T
x
+
u
1
+
u
2
-
R
x
2
d
x
.Here, ω denotes the bilateral filter operator. Recall that the bilateral weight of function T
(
p
) is defined as (10)
W
p
,
q
=
1
Z
p
exp
⁡
-
p
-
q
2
σ
s
2
exp
⁡
-
T
p
-
T
q
2
σ
r
2
,where σ
s and σ
r are the spatial and range (intensity) deviations and Z
p is a normalizing parameter to ensure that (11)
∑
q
∈
N
p
W
p
,
q
=
1in the neighborhood N
(
p
) of p. The output of the bilateral filter is then given by (12)
ω
∘
T
p
=
∑
q
∈
N
p
W
p
,
q
T
q
.Remark that the local-global similarity measure (9) is a generalization of SSD measure. Since the bilateral weight depends not only on spatial distance but also on the intensity difference, the local-global similarity measure can improve the matching accuracy and robustness for very large displacement [26].Secondly, we construct the regularization term. As aforementioned in the introduction, TV regularization is good for piecewise constant part and fourth-order regularization is good for smooth part. Hence, it is natural to propose the combined regularization term(13)
R
1
u
=
μ
∑
j
=
1
2
∫
Ω
∇
u
1
j
d
x
+
τ
∑
j
=
1
2
∫
Ω
∇
2
u
2
j
d
x
,where μ, τ are positive balance parameters. To enhance numerical stability, we add L2-norm regularization of u
1 and u
2, that is, (14)
R
2
u
=
α
2
∑
j
=
1
2
∫
Ω
u
1
j
2
d
x
+
β
2
∑
j
=
1
2
∫
Ω
u
2
j
2
d
x
,where α, β are positive balance parameters.Combining (9), (13), and (14), the proposed model is formulated as (15)
min
u
1
,
u
2
⁡
E
u
1
,
u
2
=
F
ω
u
+
R
1
u
+
R
2
u
.Generally, variational formulation would be used for this minimization problem directly. However, it is a challenging task for the nonlinear term in (15). Motivated by the technique in [27], we introduce an additional variable v
=
(
v
1
,
…
,
v
2
) and consider the following minimization problem: (16)
min
u
1
,
u
2
,
v
⁡
F
ω
v
+
R
1
u
+
R
2
u
+
1
2
θ
∫
Ω
u
1
x
+
u
2
x
-
v
x
2
d
x
,where the penalizing parameter θ is very small such that u
1
(
x
)
+
u
2
(
x
)
≈
v
(
x
). To solve the optimization problem (16), we propose an iterative algorithm by alternating minimization method.v-Subproblem. Fixing u
1, u
2 and denoting u
=
u
1
+
u
2, the subproblem for v can be written as the following minimization problem: (17)
min
v
⁡
E
11
v
=
F
ω
v
+
1
2
θ
∫
Ω
u
x
-
v
x
2
d
x
.Define r
(
v
)
:
=
T
(
x
+
v
(
x
)
)
-
R
(
x
). Then, we can rewrite (18)
F
ω
v
=
∫
Ω
ω
∘
r
v
2
d
x
.Since u
(
x
) is approximate to v
(
x
), we can utilize the first-order Taylor approximation to linearize r
(
v
) near u
(
x
); that is, (19)
r
v
=
r
u
+
v
-
u
,
∇
T
x
+
u
.A point-wise minimization is employed for problem (17). Taking the first-order derivative of E
11 with respect to v and setting it to zero, we get (20)
ω
∘
2
r
v
∇
T
x
+
u
x
+
1
θ
v
x
-
u
x
=
0
.Using the definition of r
(
v
), by direct computation, we get that (20) is equivalent to the following linear system:(21)
1
+
2
θ
ω
T
x
2
x
+
u
2
θ
ω
T
x
x
+
u
T
y
x
+
u
2
θ
ω
T
y
x
+
u
T
x
x
+
u
1
+
2
θ
ω
T
y
2
x
+
u
v
1
v
2
=
u
1
-
2
θ
ω
T
x
x
+
u
r
0
u
2
-
2
θ
ω
T
y
x
+
u
r
0
,where u
=
u
1
,
u
2
T, v
=
v
1
,
v
2
T, and r
0
=
T
(
x
+
u
)
-
R
(
x
)
-
〈
u
,
∇
T
(
x
+
u
)
〉. Since the determinant of the left-hand side matrix in (21) is always nonzero, this linear system can be solved efficiently.u
1-Subproblem. Fixing v and u
2, the subproblem for u
1 is (22)
min
u
1
⁡
E
12
u
1
=
μ
∑
j
=
1
2
∫
Ω
∇
u
1
j
d
x
+
α
2
∑
j
=
1
2
∫
Ω
u
1
j
2
d
x
+
1
2
θ
∫
Ω
u
1
x
+
u
2
x
-
v
x
2
d
x
.For the anisotropic total variation term, many efforts have devoted to obtain fast numerical schemes and overcome the nondifferentiability of the term. In this paper, we apply the efficient split Bregman algorithm [28] for each d
=
1,2 in the following: Initialize:u
1
d
0
=
v
d
-
u
2
d, here d
x
0
=
d
y
0
=
b
x
0
=
b
y
0
=
0;
Whileu
1
d
k
-
u
1
d
k
+
1
2
>
t
h
r
e
s
h
o
l
d
(23)
u
1
d
k
+
1
=
G
k
;
d
x
k
+
1
=
s
h
r
i
n
k
∇
x
u
1
d
k
+
1
+
b
x
k
,
μ
λ
;
d
y
k
+
1
=
s
h
r
i
n
k
∇
y
u
1
d
k
+
1
+
b
y
k
,
μ
λ
;
b
x
k
+
1
=
b
x
k
+
∇
x
u
1
d
k
+
1
-
d
x
k
+
1
;
b
y
k
+
1
=
b
y
k
+
∇
y
u
1
d
k
+
1
-
d
y
k
+
1
;
endwhere (24)
G
k
i
,
j
=
θ
λ
μ
1
+
α
θ
+
4
λ
μ
θ
u
1
d
,
i
+
1
,
j
k
+
u
1
d
,
i
-
1
,
j
k
+
u
1
d
,
i
,
j
+
1
k
+
u
1
d
,
i
,
j
-
1
k
+
d
1
,
i
-
1
,
j
k
-
d
1
,
i
,
j
k
+
d
2
,
i
,
j
-
1
k
-
d
2
,
i
,
j
k
-
b
1
,
i
-
1
,
j
k
+
b
1
,
i
,
j
k
-
b
2
,
i
,
j
-
1
k
+
b
2
,
i
,
j
k
+
1
1
+
θ
α
+
4
λ
v
d
,
i
,
j
-
u
2
d
,
i
,
j
.And λ
>
0 is the parameter introduced by split Bregman method. In our numerical experiments, the fast Fourier transform (FFT) is used to solve u
1
d efficiently under the periodic boundary condition.u
2-Subproblem. Fixing v and u
1, the subproblem for u
2 is (25)
min
u
2
⁡
E
13
u
2
=
τ
∑
j
=
1
2
∫
Ω
∇
2
u
2
j
d
x
+
β
2
∑
j
=
1
2
∫
Ω
u
2
j
2
d
x
+
1
2
θ
∫
Ω
u
1
x
+
u
2
x
-
v
x
2
d
x
.This minimization problem is similar to LLT model [21] in image denoising. For LLT model, Chambolle’s dual algorithm for TV denoising [29] is generalized to get an efficient algorithm [22] for the corresponding fourth-order PDE. It overcomes the numerical difficulties related to the nondifferentiability of the first regularization term. Using the same technique as in the derivation of dual algorithm of LLT model, we can easily get the dual algorithm for problem (25). The details of the numerical algorithm are given in the following: Initialize: time step sizet
∈
(
0
,
1
/
64
], space step ∇
x
=
1, ∇
y
=
1, (26)
p
0
=
p
11
p
12
p
21
p
22
=
0
.
Forn
=
0,1
,
2
,
…, do
Step 1. FindA
i
,
j
n from (27)
A
i
,
j
n
=
D
x
x
d
i
v
2
p
i
,
j
n
-
1
μ
τ
θ
v
d
,
i
,
j
-
u
1
d
,
i
,
j
D
x
y
+
d
i
v
2
p
i
,
j
n
-
1
μ
τ
θ
v
d
,
i
,
j
-
u
1
d
,
i
,
j
D
y
x
-
d
i
v
2
p
i
,
j
n
-
1
μ
τ
θ
v
d
,
i
,
j
-
u
1
d
,
i
,
j
D
y
y
d
i
v
2
p
i
,
j
n
-
1
μ
τ
θ
v
d
,
i
,
j
-
u
1
d
,
i
,
j
,
whered
i
v
2
p
i
,
j
n
=
D
x
x
p
i
,
j
11
+
D
x
y
+
p
i
,
j
21
+
D
y
x
-
p
i
,
j
12
+
D
y
y
p
i
,
j
22.
Step 2. Findp
i
,
j
n from(28)
p
i
,
j
n
+
1
=
p
i
,
j
n
-
t
A
i
,
j
n
1
+
t
A
i
,
j
n
.
Step 3. Test convergence. If it is convergent, then go to Step 4, otherwise go to Step 1.
Step 4. Findu
2
d from (29)
u
2
d
,
i
,
j
=
v
d
,
i
,
j
-
u
i
d
,
i
,
j
-
τ
μ
θ
d
i
v
2
p
i
,
j
n
+
1
.
EndHere, the first-order and second-order difference operators are defined in Table 1. Note that the computation of each component can be carried out parallelly in each iteration.Table 1
Discretization used in the implementation.
Derivative
Finite difference scheme
D
x
±
(
u
i
,
j
)
±
1
Δ
x
u
i
±
1
,
j
-
u
i
,
j
D
y
±
(
u
i
,
j
)
±
1
Δ
y
u
i
,
j
±
1
-
u
i
,
j
D
x
x
(
u
i
,
j
)
1
Δ
x
D
x
+
u
i
,
j
-
D
x
+
u
i
-
1
,
j
D
x
y
±
(
u
i
,
j
)
±
1
Δ
y
D
x
±
u
i
,
j
±
1
-
D
x
±
u
i
,
j
D
y
x
±
(
u
i
,
j
)
±
1
Δ
x
D
y
±
u
i
±
1
,
j
-
D
y
±
u
i
,
j
D
y
y
(
u
i
,
j
)
1
Δ
y
D
y
+
u
i
,
j
-
D
y
+
u
i
,
j
-
1Finally, a coarse-to-fine multiresolution approach [30] is used to improve the performance for large deformation in our algorithm. The deformation field u
1, u
2, and v are propagated from the low resolution image derived from the original image to the next higher resolution. Then, our algorithm is summarized in Algorithm 1.Algorithm 1:The proposed algorithm.
(1)  Set up pyramidsT
(
x
), R
(
x
), and their derivatives at each layer.
(2)  (a)  If at the coarsest layer, initialize:u
1
0
=
u
2
0
=
0.
(a′)  If not (a), upsampleu
1 and u
2 from the previous layer to the current layer.
WarpR
,
T
,
T
x
,
T
y using bilinear interpolation.
Fork
=
0,1
,
…
v
k
+
1
=
arg
⁡
min
⁡
v
ε
11
v
,
u
1
k
,
u
2
k   (see (17))
u
1
k
+
1
=
arg
⁡
min
⁡
u
1
ε
12
v
k
+
1
,
u
2
k
,
u
1  (see (22))
u
2
k
+
1
=
arg
⁡
min
⁡
u
2
ε
13
v
k
+
1
,
u
2
,
u
1
k
+
1   (see (25))
end, go to the finer layer.
end of the finest layer.


## 4. Numerical Experiments

In this section, we present experimental results of our method and compare them with some closely related variational methods.For quantitative comparison, we compute the pixel-wise root mean square error (RMSE), the correlation coefficient (CC), and the mutual information (MI) between two imagesT
(
x
) and S
(
x
). The definitions of RMSE, CC, and MI are as follows.(i)
RMSE: R
M
S
E
(
T
(
x
)
,
S
(
x
)
)
=
∑
‍
(
(
T
(
x
)
-
S
(
x
)
)
2
/
n
), where n is the size of T
(
x
).
(ii)
CC: c
o
r
r
(
T
(
x
)
,
S
(
x
)
)
=
c
o
v
(
T
,
S
)
/
σ
S
σ
T, where c
o
v is the covariance and σ
T and σ
S are the standard deviation of T
(
x
) and S
(
x
), respectively.
(iii)
MI: M
I
(
T
(
x
)
,
S
(
x
)
)
=
H
(
T
)
+
H
(
S
)
-
H
(
T
,
S
), where H
(
T
) and H
(
S
) are the marginal entropies and H
(
T
,
S
) is the joint entropy of T
(
x
) and S
(
x
).There are several regularization parametersα, β, μ, and τ in our model (15). The regularization parameters α and β in the L
2-norm regularization term (14) affect the stability of the proposed algorithm. Larger parameters α, β ensure more stability of the algorithm. However, larger parameters lead to slower convergence, as shown in Figures 1(a)-1(b) (which corresponds to Test 3 in Section 4.3). The regularization parameters μ, τ control the smoothness of the estimated deformation field. By increasing μ or τ, smoother deformation field can be obtained. Parameters μ and τ are ranging from 0 to 103 in this paper. The parameter λ in solving u
1 is tuned in each experiment to obtain optimal result. Theoretically the penalizing parameter θ should be selected to be as small as possible. However, as shown in our experiments, θ in the range [1/1000 1/10] is good for the registration problems.Figure 1
The effect of different parametersα and β in our method: (a) the sum of squared difference SSD(
T
(
u
+
x
)
,
R
(
x
)
), (b) the correlation coefficient corr(
T
(
u
+
x
)
,
R
(
x
)
).
(a)
SSDT
u
+
x
,
R
x
(b)
Corr(
T
(
u
+
x
)
,
R
(
x
)
)To verify the robustness and accuracy of our method, a series of experiments are presented. Moreover, we also implement the registration models with SSD measure fidelity term and three kinds of regularization for comparison. The regularization method includesT
V regularization (6), L
2 regularization (5) and fourth-order regularization (7). The three methods are abbreviated as T
V, L
2 and fourth-order, respectively, where the T
V-diffusion model was solved efficiently using the Bregman split method and the fourth-order diffusion model was solved using the Chambolle’s dual algorithm. Note that all the involved parameters of the above three methods are tuned for each experiment to get optimal result.All the experiments are performed under Windows 7 and MATLAB R2013a with Intel Core i3-4130 CPU@3.40 GHz and 8 GB memory. The programming language is MATLAB.


### 4.1. Test 1

We first consider matching a pair of 2D synthetic images shown in Figures2(a) and 2(b) for “sliding rectangular,” which are almost piecewise constant. The registration results and the transformation of coordinate grid by the four methods are displayed in Figures 2(c)–2(f) and 3(a)–3(d), respectively.Figure 2
Test 1: Registered image for a pair of synthetic images shown in (a)T and (b) R of size 473 × 449. (c) T
V-diffusion model, (d) L
2-diffusion model, (e) Fourth-order-diffusion model, and (f) our model. Parameters: α
=
500, β
=
500, μ
=
10, τ
=
1, λ
=
300, and θ
=
0.1.
(a)
Template imageT
(b)
Reference imageR
(c)
T
(
x
+
u
) (T
V)
(d)
T
(
x
+
u
) (L
2)
(e)
T
(
x
+
u
) (fourth-order)
(f)
T
(
x
+
u
) (our model)Figure 3
The transformation presentations for different methods in Test 1.
(a)
T
V
(b)
L
2
(c)
Fourth-order
(d)
Our modelIn this experiment, the correlation coefficient between the template imageT
(
x
) and the reference image R
(
x
) before registered is corr(
T
(
x
)
,
R
(
x
)
)
=
0.8.It is obvious that all the four methods are effective in this test. The correlation coefficients and the sum of squared differences between the reference imageR
(
x
) and the deformed images T
(
x
+
u
) from four models are all close to the best results after 200 iterations, which are shown in Figures 4(a)-4(b). Meanwhile, from the dissimilarity between R
(
x
) and T
(
x
+
u
) shown in Figures 5(a)–5(d), we can see that the T
V-diffusion model and our model are more accurate for discontinuities, while the L
2-diffusion model and the fourth-order model yield oversmooth on the edge of the object. Furthermore, in order to get a comprehensive evaluation for comparison, the RMSE values and MI values for all the registration results with the four registration methods are reported in Table 2. From the measures, we find that our method achieves the best performance among all.Table 2
Summary of RMSE, mean, and MI. The images are rescaled to [0,1].
Image pairs
Method
RMSE  (1
0
-
4)
Mean (1
0
-
4)
IM
Test 1: sliding rectangular
TV
380.0
72
1.1615
L
2
274.1
36
1.1610
Fourth-order
343.2
26
1.1358
Proposed
248
24
1.1620
Test 2: circle to ellipse
TV
708
111
0.7040
L
2
607
105
0.7701
Fourth-order
794
148
0.6815
Proposed
440
83
0.8164Figure 4
Comparisons of the correlation coefficient corr(
T
(
u
+
x
)
,
R
(
x
)
) and the sum of squared difference SSD(
T
(
u
+
x
)
,
R
(
x
)
) for different models in Test 1.
(a)
Corr(
T
(
u
+
x
)
,
R
(
x
)
)
(b)
SSD(
T
(
u
+
x
)
,
R
(
x
)
)Figure 5
The residues|
T
(
x
+
u
)
-
R
(
x
)
| after deformation using different methods in Test 1.
(a)
T
V
(b)
L
2
(c)
Fourth-order
(d)
Our model


### 4.2. Test 2

In Test 2, a smooth registration problem “circle to ellipse” is considered, where the source images are shown in Figures6(a)-6(b). The results and the deformation fields by the four models are displayed in Figures 6(c)–6(f) and 7(a)–7(d), respectively. Carefully comparing Figures 8(a)–8(d), we notice that the result of T
V model is poor for the smooth deformation field, and our model performs well on smooth deformations. Furthermore, the performance measures RMSE and MI are reported in Table 2. Again we find that our method is the best.Figure 6
Test 2: Registered images for “circle to ellipse” (the smooth registration problem) shown in (a)T and (b) R of size 100 × 100. (c) T
V-diffusion model, (d) L
2-diffusion model, (e) Fourth-order-diffusion model, (f) Our model. Parameters: α
=
100, β
=
10, μ
=
100, τ
=
10, λ
=
100, and θ
=
0.01.
(a)
Template imageT
(b)
Reference imageR
(c)
T
x
+
u (T
V)
(d)
T
(
x
+
u
) (L
2)
(e)
T
(
x
+
u
) (fourth-order)
(f)
T
(
x
+
u
) (our model)Figure 7
The transformation presentations for different methods in Test 2.
(a)
T
V
(b)
L
2
(c)
Fourth-order
(d)
Our modelFigure 8
The residues|
T
(
x
+
u
)
-
R
(
x
)
| after deformation using different methods in Test 2.
(a)
T
V
(b)
L
2
(c)
Fourth-order
(d)
Our model


### 4.3. Test 3

In the third experiment, we apply our model on affine-linear transformation problem Figures9(a)-9(b). Although the regularization method is widely used for nonrigid models, it is not usually used in affine-linear transformation problem. For the affine-line displacement, it is difficult to be penalized in the interior of the image domain [31]. Experiment results show that our method performs better than the other three methods, see Figures 9(c)–9(f), 10(a)–10(d), and 11(a)–11(d).Figure 9
Test 3: Test on affine-linear transformations problem. (a) and (b) are a pair of images to be registered, size 414 × 414. (c)–(f) are the registered results of different methods. Parameters:α
=
50, β
=
50, μ
=
1, τ
=
10, λ
=
1, and θ
=
0.1.
(a)
Template imageT
(b)
Reference imageR
(c)
T
(
x
+
u
) (T
V)
(d)
T
(
x
+
u
) (L
2)
(e)
T
(
x
+
u
) (fourth-order)
(f)
T
(
x
+
u
) (our model)Figure 10
The transformation presentations for different methods in Test 3.
(a)
T
V
(b)
L
2
(c)
Fourth-order
(d)
Our modelFigure 11
The residues|
T
(
x
+
u
)
-
R
(
x
)
| after deformation using different methods for the affine-linear transformation problem in Test 3.
(a)
T
V
(b)
L
2
(c)
Fourth-order
(d)
Our model


### 4.4. Test 4

Finally, We consider a pair of 2D clinic images of dimension 320 × 317 that focuses on lung area, shown in Figures12(a)-12(b). The registered images obtained from the different methods shown in Figures 12(c)–12(f) are almost identical. Comparing the dissimilarity between the deformed images T
(
x
+
u
) and the reference image R
(
x
), we find that our method is more robust than the other methods. It is obvious that in Figure 13(a) much information is lost in smooth region. Meanwhile, Figure 13(c) fails to preserve the edge of object.Figure 12
Test 4: test on the real clinic images. (a) and (b) are a pair of images to be registered, size 320 × 317. (c)–(f) are the registered results of different methods. Parameters:α
=
100, β
=
10, μ
=
100, τ
=
10, λ
=
100, and θ
=
0.1.
(a)
Template imageT
(b)
Reference imageR
(c)
T
(
x
+
u
) (T
V)
(d)
T
(
x
+
u
) (L
2)
(e)
T
(
x
+
u
) (fourth-order)
(f)
T
(
x
+
u
) (our model)Figure 13
The residues|
T
(
x
+
u
)
-
R
(
x
)
| after deformation using different methods with the clinic images.
(a)
T
V
(b)
L
2
(c)
Fourth-order
(d)
Our modelIn Figures14(a)–14(d), we display the enlarged portions of the deformed images T
(
x
+
u
) from the four methods for detail comparison. Obviously, staircase effect occurs in Figure 14(a). The deformed images by our method in Figure 14(d) is the most close to the reference image in Figure 12(b).Figure 14
The enlarged portions of the deformed imagesT
(
x
+
u
) by different methods with the clinic images.
(a)
T
V
(b)
L
2
(c)
Fourth-order
(d)
Our modelBy the above four experiments, it can be seen that our method is effective not only for smooth data but also nonsmooth data. It takes the advantage of both total variational regularization and fourth-order regularization method while overcomes their demerits.


## 4.1. Test 1

We first consider matching a pair of 2D synthetic images shown in Figures2(a) and 2(b) for “sliding rectangular,” which are almost piecewise constant. The registration results and the transformation of coordinate grid by the four methods are displayed in Figures 2(c)–2(f) and 3(a)–3(d), respectively.Figure 2
Test 1: Registered image for a pair of synthetic images shown in (a)T and (b) R of size 473 × 449. (c) T
V-diffusion model, (d) L
2-diffusion model, (e) Fourth-order-diffusion model, and (f) our model. Parameters: α
=
500, β
=
500, μ
=
10, τ
=
1, λ
=
300, and θ
=
0.1.
(a)
Template imageT
(b)
Reference imageR
(c)
T
(
x
+
u
) (T
V)
(d)
T
(
x
+
u
) (L
2)
(e)
T
(
x
+
u
) (fourth-order)
(f)
T
(
x
+
u
) (our model)Figure 3
The transformation presentations for different methods in Test 1.
(a)
T
V
(b)
L
2
(c)
Fourth-order
(d)
Our modelIn this experiment, the correlation coefficient between the template imageT
(
x
) and the reference image R
(
x
) before registered is corr(
T
(
x
)
,
R
(
x
)
)
=
0.8.It is obvious that all the four methods are effective in this test. The correlation coefficients and the sum of squared differences between the reference imageR
(
x
) and the deformed images T
(
x
+
u
) from four models are all close to the best results after 200 iterations, which are shown in Figures 4(a)-4(b). Meanwhile, from the dissimilarity between R
(
x
) and T
(
x
+
u
) shown in Figures 5(a)–5(d), we can see that the T
V-diffusion model and our model are more accurate for discontinuities, while the L
2-diffusion model and the fourth-order model yield oversmooth on the edge of the object. Furthermore, in order to get a comprehensive evaluation for comparison, the RMSE values and MI values for all the registration results with the four registration methods are reported in Table 2. From the measures, we find that our method achieves the best performance among all.Table 2
Summary of RMSE, mean, and MI. The images are rescaled to [0,1].
Image pairs
Method
RMSE  (1
0
-
4)
Mean (1
0
-
4)
IM
Test 1: sliding rectangular
TV
380.0
72
1.1615
L
2
274.1
36
1.1610
Fourth-order
343.2
26
1.1358
Proposed
248
24
1.1620
Test 2: circle to ellipse
TV
708
111
0.7040
L
2
607
105
0.7701
Fourth-order
794
148
0.6815
Proposed
440
83
0.8164Figure 4
Comparisons of the correlation coefficient corr(
T
(
u
+
x
)
,
R
(
x
)
) and the sum of squared difference SSD(
T
(
u
+
x
)
,
R
(
x
)
) for different models in Test 1.
(a)
Corr(
T
(
u
+
x
)
,
R
(
x
)
)
(b)
SSD(
T
(
u
+
x
)
,
R
(
x
)
)Figure 5
The residues|
T
(
x
+
u
)
-
R
(
x
)
| after deformation using different methods in Test 1.
(a)
T
V
(b)
L
2
(c)
Fourth-order
(d)
Our model


## 4.2. Test 2

In Test 2, a smooth registration problem “circle to ellipse” is considered, where the source images are shown in Figures6(a)-6(b). The results and the deformation fields by the four models are displayed in Figures 6(c)–6(f) and 7(a)–7(d), respectively. Carefully comparing Figures 8(a)–8(d), we notice that the result of T
V model is poor for the smooth deformation field, and our model performs well on smooth deformations. Furthermore, the performance measures RMSE and MI are reported in Table 2. Again we find that our method is the best.Figure 6
Test 2: Registered images for “circle to ellipse” (the smooth registration problem) shown in (a)T and (b) R of size 100 × 100. (c) T
V-diffusion model, (d) L
2-diffusion model, (e) Fourth-order-diffusion model, (f) Our model. Parameters: α
=
100, β
=
10, μ
=
100, τ
=
10, λ
=
100, and θ
=
0.01.
(a)
Template imageT
(b)
Reference imageR
(c)
T
x
+
u (T
V)
(d)
T
(
x
+
u
) (L
2)
(e)
T
(
x
+
u
) (fourth-order)
(f)
T
(
x
+
u
) (our model)Figure 7
The transformation presentations for different methods in Test 2.
(a)
T
V
(b)
L
2
(c)
Fourth-order
(d)
Our modelFigure 8
The residues|
T
(
x
+
u
)
-
R
(
x
)
| after deformation using different methods in Test 2.
(a)
T
V
(b)
L
2
(c)
Fourth-order
(d)
Our model


## 4.3. Test 3

In the third experiment, we apply our model on affine-linear transformation problem Figures9(a)-9(b). Although the regularization method is widely used for nonrigid models, it is not usually used in affine-linear transformation problem. For the affine-line displacement, it is difficult to be penalized in the interior of the image domain [31]. Experiment results show that our method performs better than the other three methods, see Figures 9(c)–9(f), 10(a)–10(d), and 11(a)–11(d).Figure 9
Test 3: Test on affine-linear transformations problem. (a) and (b) are a pair of images to be registered, size 414 × 414. (c)–(f) are the registered results of different methods. Parameters:α
=
50, β
=
50, μ
=
1, τ
=
10, λ
=
1, and θ
=
0.1.
(a)
Template imageT
(b)
Reference imageR
(c)
T
(
x
+
u
) (T
V)
(d)
T
(
x
+
u
) (L
2)
(e)
T
(
x
+
u
) (fourth-order)
(f)
T
(
x
+
u
) (our model)Figure 10
The transformation presentations for different methods in Test 3.
(a)
T
V
(b)
L
2
(c)
Fourth-order
(d)
Our modelFigure 11
The residues|
T
(
x
+
u
)
-
R
(
x
)
| after deformation using different methods for the affine-linear transformation problem in Test 3.
(a)
T
V
(b)
L
2
(c)
Fourth-order
(d)
Our model


## 4.4. Test 4

Finally, We consider a pair of 2D clinic images of dimension 320 × 317 that focuses on lung area, shown in Figures12(a)-12(b). The registered images obtained from the different methods shown in Figures 12(c)–12(f) are almost identical. Comparing the dissimilarity between the deformed images T
(
x
+
u
) and the reference image R
(
x
), we find that our method is more robust than the other methods. It is obvious that in Figure 13(a) much information is lost in smooth region. Meanwhile, Figure 13(c) fails to preserve the edge of object.Figure 12
Test 4: test on the real clinic images. (a) and (b) are a pair of images to be registered, size 320 × 317. (c)–(f) are the registered results of different methods. Parameters:α
=
100, β
=
10, μ
=
100, τ
=
10, λ
=
100, and θ
=
0.1.
(a)
Template imageT
(b)
Reference imageR
(c)
T
(
x
+
u
) (T
V)
(d)
T
(
x
+
u
) (L
2)
(e)
T
(
x
+
u
) (fourth-order)
(f)
T
(
x
+
u
) (our model)Figure 13
The residues|
T
(
x
+
u
)
-
R
(
x
)
| after deformation using different methods with the clinic images.
(a)
T
V
(b)
L
2
(c)
Fourth-order
(d)
Our modelIn Figures14(a)–14(d), we display the enlarged portions of the deformed images T
(
x
+
u
) from the four methods for detail comparison. Obviously, staircase effect occurs in Figure 14(a). The deformed images by our method in Figure 14(d) is the most close to the reference image in Figure 12(b).Figure 14
The enlarged portions of the deformed imagesT
(
x
+
u
) by different methods with the clinic images.
(a)
T
V
(b)
L
2
(c)
Fourth-order
(d)
Our modelBy the above four experiments, it can be seen that our method is effective not only for smooth data but also nonsmooth data. It takes the advantage of both total variational regularization and fourth-order regularization method while overcomes their demerits.


## 5. Conclusion

In this work, we have proposed a novel variational approach for nonrigid image registration by employing a local-global similarity measure and a combinational regularizer. We have assumed that the deformation field can be decomposed into a piecewise constant component and a smooth component. Then, by taking use the advantages of the total variational diffusion and the fourth-order diffusion, each component was described appropriately. A combination of LK method, the split Bregman algorithm and the dual algorithm was introduced to get an efficient algorithm. The proposed method has provided satisfactory results for matching both non-smooth and smooth structures in nonrigid image registration problem.In the future work, we will study the automatical choice of parameters in our model and extend the proposed method for multimodality image registration.



---

*Source: 235314-2015-05-31.xml*