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 endwhere 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 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.
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.
| ||||||||||||||||||||
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 103 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
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 [email protected] 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)
(a)
(b)
(c) Fourth-order
(d) Our model
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.
(a) Corr
(b) SSD
(a)
(b)
(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 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)
(a)
(b)
(c) Fourth-order
(d) Our model
(a)
(b)
(c) Fourth-order
(d) Our model
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)
(a)
(b)
(c) Fourth-order
(d) Our model
(a)
(b)
(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 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)
(a)
(b)
(c) Fourth-order
(d) Our model
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
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).