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

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 similaritymeasure is used in ourmethod to improve the accuracy and robustness for imagematching. 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.


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 inversioninvariant 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][12][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 localglobal 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][22][23].It is reported that fourthorder diffusion damps oscillations much faster than secondorder 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 Section 2, 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.

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 and the template image by  : Ω → R. In nonrigid image registration, the deformation field is usually described by displacement field u : Ω → Ω denoted by The problem is seeking for an deformation field u such that the transformed template image is optimally correlated with the reference image (x); that is, 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: 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 Section 1, the frequently used similarity measure is SSD measure [4] Meawhile, the widely used regularization term is quadratic  2 regularization term [17,18] and TV regularization term [19] In image restoration problems, the following fourth-order regularization term is introduced:

The Proposed Model and Algorithm
We assume that the deformation field, which transform the template image  to the reference image , 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 Firstly, we present the local-global similarity measure based on the bilateral filter [26]: Here,  denotes the bilateral filter operator.Recall that the bilateral weight of function (p) is defined as where   and   are the spatial and range (intensity) deviations and  p is a normalizing parameter to ensure that in the neighborhood N(p) of p.The output of the bilateral filter is then given by 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 where ,  are positive balance parameters.To enhance numerical stability, we add L2-norm regularization of u 1 and u 2 , that is, where ,  are positive balance parameters.Combining ( 9), (13), and ( 14), the proposed model is formulated as 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 k = (V 1 , . . ., V 2 ) and consider the following minimization problem: where the penalizing parameter  is very small such that u 1 (x) + u 2 (x) ≈ k(x).To solve the optimization problem ( 16), we propose an iterative algorithm by alternating minimization method.
k-Subproblem.Fixing u 1 , u 2 and denoting u = u 1 + u 2 , the subproblem for k can be written as the following minimization problem: Define (k) := (x + k(x)) − (x).Then, we can rewrite Since u(x) is approximate to k(x), we can utilize the first-order Taylor approximation to linearize (k) near u(x); that is, A point-wise minimization is employed for problem (17).Taking the first-order derivative of E 11 with respect to k and setting it to zero, we get Using the definition of (k), by direct computation, we get that ( 20) is equivalent to the following linear system: where u = ( 1 ,  2 )  , k = (V 1 , V 2 )  , and  0 = (x + u) − (x) − ⟨u, ∇(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 k and u 2 , the subproblem for u 1 is 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  = 1, 2 in the following: where And  > 0 is the parameter introduced by split Bregman method.In our numerical experiments, the fast Fourier transform (FFT) is used to solve  1 efficiently under the periodic boundary condition.
u 2 -Subproblem.Fixing k and u 1 , the subproblem for u 2 is 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 size  ∈ (0, 1/64], space step ∇ = 1, ∇ = 1, For  = 0, 1, 2, . .., do Step 1. Find   , from where div End Table 1: Discretization used in the implementation.

Derivative
Finite difference scheme Here, 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.
Finally, 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 k 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.

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 images (x) and (x).The definitions of RMSE, CC, and MI are as follows.
where  is the size of (x).
end, go to the finer layer.end of the finest layer.There are several regularization parameters , , , and  in our model (15).The regularization parameters  and  in the  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 10 3 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.
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 includes TV regularization (6),  2 regularization (5) and fourth-order regularization (7).The three methods are abbreviated as TV,  2 and fourthorder, respectively, where the TV-diffusion model was solved efficiently using the Bregman split method and the fourthorder 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.40GHz and 8 GB memory.The programming language is MATLAB.

Test 1.
We first consider matching a pair of 2D synthetic images shown in Figures 2(a In this experiment, the correlation coefficient between the template image (x) and the reference image (x) before registered is corr((x), (x)) = 0.8.  2. From the measures, we find that our method achieves the best performance among all.2. Again we find that our method is the best.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  By 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.

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

The L 2 Figure 4 :
Figure 4: Comparisons of the correlation coefficient corr((u + x), (x)) and the sum of squared difference SSD((u + x), (x)) for different models in Test 1.

Figure 7 :
Figure 7: The transformation presentations for different methods in Test 2.

Figure 10 :
Figure 10: The transformation presentations for different methods in Test 3.

Algorithm 1 :Figure 11 :
Figure 11: The residues |(x + u) − (x)| after deformation using different methods for the affine-linear transformation problem in Test 3.
) 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 Figures2(c)-2(f) and 3(a)-3(d), respectively.

Figure 13 :
Figure 13: The residues |(x + u) − (x)| after deformation using different methods with the clinic images.

4. 2 .
Test 2. In Test 2, a smooth registration problem "circle to ellipse" is considered, where the source images are shown in Figures 6(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 TV 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

Figure 14 :
Figure 14: The enlarged portions of the deformed images (x + u) by different methods with the clinic images.

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 Figures12(c)-12(f) are almost identical.Comparing the dissimilarity between the deformed images (x + u) and the reference image (x), we find that our method is more robust than the other methods.It is obvious that in Figure13(a) much information is lost in smooth region.Meanwhile, Figure13(c) fails to preserve the edge of object.In Figures14(a)-14(d), we display the enlarged portions of the deformed images (x + u) from the four methods for detail comparison.Obviously, staircase effect occurs in Figure14(a).The deformed images by our method in Figure14(d) is the most close to the reference image in Figure 12(b).