Mathematical Problems in Engineering
Volume 2015 (2015), Article ID 235314, 16 pages
http://dx.doi.org/10.1155/2015/235314

Research Article

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

Xiaomei Yang,¹ Chaomin Shen,² Fang Li,¹ and Chunli Shen¹

¹Department of Mathematics, East China Normal University, Shanghai 200241, China
²Joint Laboratory for Imaging Science Technology and Department of Computer Science, East China Normal University, Shanghai 200241, China

Received 9 December 2014; Revised 31 March 2015; Accepted 31 March 2015

Academic Editor: Giuseppe Vairo

Copyright © 2015 Xiaomei Yang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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.

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

2. Preliminaries

A general framework for nonrigid image registration can be formulated in the following. We assume that is a bounded domain with Lipschitz boundary and satisfying the cone condition. Let us denote the reference image by and the template image by . In nonrigid image registration, the deformation field is usually described by displacement field denoted by

The problem is seeking for an deformation field such that the transformed template image is optimally correlated with the reference image ; 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 is the similarity measure between the transformed image and the reference image, is the regularization term, and is a weight parameter to balance the influence of and .

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 regularization term [17, 18]

and TV regularization term [19]

In image restoration problems, the following fourth-order regularization term is introduced:

3. 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 , where is the piecewise constant part and 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 is defined as

where and are the spatial and range (intensity) deviations and is a normalizing parameter to ensure that

in the neighborhood of . 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 and , 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 and consider the following minimization problem:

where the penalizing parameter is very small such that . To solve the optimization problem (16), we propose an iterative algorithm by alternating minimization method.

-Subproblem. Fixing , and denoting , the subproblem for can be written as the following minimization problem:

Define . Then, we can rewrite

Since is approximate to , we can utilize the first-order Taylor approximation to linearize near ; that is,

A point-wise minimization is employed for problem (17). Taking the first-order derivative of with respect to and setting it to zero, we get

Using the definition of , by direct computation, we get that (20) is equivalent to the following linear system:

where , , and . Since the determinant of the left-hand side matrix in (21) is always nonzero, this linear system can be solved efficiently.

-Subproblem. Fixing and , the subproblem for 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 in the following:

 Initialize: , here ; While

 end

where

And is the parameter introduced by split Bregman method. In our numerical experiments, the fast Fourier transform (FFT) is used to solve efficiently under the periodic boundary condition.

-Subproblem. Fixing and , the subproblem for 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 , space step , ,

 For , do Step 1. Find from

 where . Step 2. Find from

 Step 3. Test convergence. If it is convergent, then go to Step 4, otherwise go to Step 1. Step 4. Find from

 End

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.

Table 1: Discretization used in the implementation.

Derivative
Finite difference scheme

Finally, a coarse-to-fine multiresolution approach [30] is used to improve the performance for large deformation in our algorithm. The deformation field , , and 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.

(1) Set up pyramids , , and their derivatives at each layer.

(2) (a) If at the coarsest layer, initialize: .

      (a′) If not (a), upsample and from the previous layer to the current layer.

      Warp using bilinear interpolation.

      For

                        (see (17))

                    (see (22))

                  (see (25))

      end, go to the finer layer.

      end of the finest layer.

Algorithm 1: The proposed algorithm.

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 images and . The definitions of RMSE, CC, and MI are as follows.

(i)RMSE: , where is the size of .(ii)CC: , where is the covariance and and are the standard deviation of and , respectively.(iii)MI: , where and are the marginal entropies and is the joint entropy of and .

There are several regularization parameters , , , and in our model (15). The regularization parameters and in the -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³ in this paper. The parameter in solving 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.

(a) SSD

(b) Corr

Figure 1: The effect of different parameters and in our method: (a) the sum of squared difference SSD, (b) the correlation coefficient corr.

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 regularization (6), regularization (5) and fourth-order regularization (7). The three methods are abbreviated as , and fourth-order, respectively, where the -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 Figures 2(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.

(a) Template image

(b) Reference image

(c) ()

(d) ()

(e) (fourth-order)

(f) (our model)

Figure 2: Test 1: Registered image for a pair of synthetic images shown in (a) and (b) of size 473 × 449. (c) -diffusion model, (d) -diffusion model, (e) Fourth-order-diffusion model, and (f) our model. Parameters: , , , , , and .

(a)

(b)

(c) Fourth-order

(d) Our model

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

In this experiment, the correlation coefficient between the template image and the reference image before registered is corr.

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 image and the deformed images 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 and shown in Figures 5(a)–5(d), we can see that the -diffusion model and our model are more accurate for discontinuities, while the -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 [].

Image pairs
Method
RMSE ()
Mean ()
IM

Test 1: sliding rectangular
TV
380.0
72
1.1615

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

607
105
0.7701

Fourth-order
794
148
0.6815

Proposed
440
83
0.8164

(a) Corr

(b) SSD

Figure 4: Comparisons of the correlation coefficient corr and the sum of squared difference SSD for different models in Test 1.

(a)

(b)

(c) Fourth-order

(d) Our model

Figure 5: The residues after deformation using different methods in Test 1.

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

(a) Template image

(b) Reference image

(c) ()

(d) ()

(e) (fourth-order)

(f) (our model)

Figure 6: Test 2: Registered images for “circle to ellipse” (the smooth registration problem) shown in (a) and (b) of size 100 × 100. (c) -diffusion model, (d) -diffusion model, (e) Fourth-order-diffusion model, (f) Our model. Parameters: , , , , , and .

(a)

(b)

(c) Fourth-order

(d) Our model

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

(a)

(b)

(c) Fourth-order

(d) Our model

Figure 8: The residues after deformation using different methods in Test 2.

4.3. Test 3

In the third experiment, we apply our model on affine-linear transformation problem Figures 9(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).

(a) Template image

(b) Reference image

(c) ()

(d) ()

(e) (fourth-order)

(f) (our model)

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: , , , , , and .

(a)

(b)

(c) Fourth-order

(d) Our model

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

(a)

(b)

(c) Fourth-order

(d) Our model

Figure 11: The residues after deformation using different methods for the affine-linear transformation problem in Test 3.

4.4. Test 4

Finally, We consider a pair of 2D clinic images of dimension 320 × 317 that focuses on lung area, shown in Figures 12(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 and the reference image , 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.

(a) Template image

(b) Reference image

(c) ()

(d) ()

(e) (fourth-order)

(f) (our model)

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: , , , , , and .

(a)

(b)

(c) Fourth-order

(d) Our model

Figure 13: The residues after deformation using different methods with the clinic images.

In Figures 14(a)–14(d), we display the enlarged portions of the deformed images 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).

(a)

(b)

(c) Fourth-order

(d) Our model

Figure 14: The enlarged portions of the deformed images by different methods with the clinic images.

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.

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.

Conflict of Interests

The authors declared that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (11461064), the 973 Program (2011CB707104), and the College Funds of Xinjiang University (no. 201204011273).

References

R. Bajcsy and S. Kovačič, “Multiresolution elastic matching,” Computer Vision, Graphics, and Image Processing, vol. 46, no. 1, pp. 1–21, 1989.

C. Broit, Optimal registration of deformed images [Ph.D. dissertation], Department of Computer and Information Science, Philadelphia University, Philadelphia, Pa, USA, 1981.

J.-P. Thirion, “Image matching as a diffusion process: an analogy with Maxwell's demons,” Medical Image Analysis, vol. 2, no. 3, pp. 243–260, 1998.

P. Cachier, X. Pennec, and N. Ayache, “Fast non rigid matching by gradient descent: study and improvements of the ‘demons’ algorithm,” Tech. Rep. RR-3706, Institut National de Recherche en Informatique et en Automatique, 1999.

P. Cachier and D. Rey, “Symmetrization of the non-rigid registration problem using inversion-invariant energies: application to multiple sclerosis,” in Medical Image Computing and Computer-Assisted Intervention—MICCAI 2000, vol. 1935 of Lecture Notes in Computer Science, pp. 472–481, Springer, Berlin, Germany, 2000.

W. Lu, M.-L. Chen, H. G. Olivera, K. J. Ruchala, and T. R. Mackie, “Fast free-form deformable registration via calculus of variations,” Physics in Medicine and Biology, vol. 49, no. 14, pp. 3067–3087, 2004.

W. Lu, G. H. Olivera, Q. Chen, M.-L. Chen, and K. J. Ruchala, “Automatic re-contouring in 4D radiotherapy,” Physics in Medicine and Biology, vol. 51, no. 5, pp. 1077–1099, 2006.

L. G. Brown, “A survey of image registration techniques,” ACM Computing Surveys (CSUR), vol. 24, no. 4, pp. 325–376, 1992.

J. Modersitzki, Numerical Methods for Image Registration, Oxford University Press, New York, NY, USA, 2003.

D. L. Collins and A. C. Evans, “Animal: validation and application of nonlinear registration-based segmentation,” International Journal of Pattern Recognition and Artificial Intelligence, vol. 11, no. 8, pp. 1271–1294, 1997.

G. Hermosillo, C. Chefd, and O. Faugeras, “Variational methods for multimodal image matching,” International Journal of Computer Vision, vol. 50, no. 3, pp. 329–343, 2002.

B. Jian, B. C. Vemuri, and J. L. Marroquin, “Robust nonrigid multimodal image registration using local frequency maps,” in Information Processing in Medical Imaging, pp. 504–515, Springer, 2005.

J. P. W. Pluim, J. B. A. A. Maintz, and M. A. Viergever, “Mutual-information-based registration of medical images: a survey,” IEEE Transactions on Medical Imaging, vol. 22, no. 8, pp. 986–1004, 2003.

Y. He, A. Ben Hamza, and H. Krim, “A generalized divergence measure for robust image registration,” IEEE Transactions on Signal Processing, vol. 51, no. 5, pp. 1211–1220, 2003.

J. P. W. Pluim, J. B. A. Maintz, and M. A. Viergever, “f-information measures in medical image registration,” IEEE Transactions on Medical Imaging, vol. 23, no. 12, pp. 1508–1516, 2004.

M. Drulea and S. Nedevschi, “Total variation regularization of local-global optical flow,” in Proceedings of the 14th IEEE International Intelligent Transportation Systems Conference (ITSC '11), pp. 318–323, IEEE, Washington, DC , USA, October 2011.

Q. Zeng and Y. Chen, “Accurate inverse consistent non-rigid image registration and its application on automatic re-contouring,” in Bioinformatics Research and Applications, vol. 4983 of Lecture Notes in Computer Science, pp. 293–304, Springer, Berlin, Germany, 2008.

B. Fischer and J. Modersitzki, “Fast diffusion registration,” Contemporary Mathematics, vol. 313, pp. 117–128, 2002.

T. Pock, M. Urschler, C. Zach, R. Beichel, and H. Bischof, “A duality based algorithm for tv-l1 optical flow image registration,” in Medical Image Computing and Computer-Assisted Intervention—MICCAI 2007, pp. 511–518, Springer, 2007.

A. Chambolle and P.-L. Lions, “Image recovery via total variation minimization and related problems,” Numerische Mathematik, vol. 76, no. 2, pp. 167–188, 1997.

M. Lysaker, A. Lundervold, and X.-C. Tai, “Noise removal using fourth-order partial differential equation with applications to medical magnetic resonance images in space and time,” IEEE Transactions on Image Processing, vol. 12, no. 12, pp. 1579–1590, 2003.

H.-Z. Chen, J.-P. Song, and X.-C. Tai, “A dual algorithm for minimization of the LLT model,” Advances in Computational Mathematics, vol. 31, no. 1–3, pp. 115–130, 2009.

F. Li, C. Shen, J. Fan, and C. Shen, “Image restoration combining a total variational filter and a fourth-order filter,” Journal of Visual Communication and Image Representation, vol. 18, no. 4, pp. 322–330, 2007.

M. Lysaker and X.-C. Tai, “Iterative image restoration combining total variation minimization and a second-order functional,” International Journal of Computer Vision, vol. 66, no. 1, pp. 5–18, 2006.

A. N. Tikhonov and V. Y. Arsenin, Solutions of Ill-Posed Problems, Wiley, 1977.

F. Durand and J. Dorsey, “Fast bilateral filtering for the display of high dynamicrange images,” ACM Transactions on Graphics, vol. 21, no. 3, pp. 257–266, 2002.

T. Vercauteren, X. Pennec, A. Perchant, and N. Ayache, “Non-parametric diffeomorphic image registration with the demons algorithm,” in Medical Image Computing and Computer-Assisted Intervention—MICCAI 2007, vol. 4792 of Lecture Notes in Computer Science, pp. 319–326, Springer, Berlin, Germany, 2007.

T. Goldstein and S. Osher, “The split bregman method for l1-regularized problems,” SIAM Journal on Imaging Sciences, vol. 2, no. 2, pp. 323–343, 2009.

A. Chambolle, “An algorithm for total variation minimization and applications,” Journal of Mathematical Imaging and Vision, vol. 20, no. 1-2, pp. 89–97, 2004.

P. J. Kostelec, J. B. Weaver, and D. M. Healy Jr., “Multiresolution elastic image registration,” Medical Physics, vol. 25, no. 9, pp. 1593–1604, 1998.

N. Chumchob, K. K. Chen, and C. Brito-Loeza, “A fourth-order variational image registration model and its fast multigrid algorithm,” Multiscale Modeling & Simulation, vol. 9, no. 1, pp. 89–128, 2011.