ISRN Biomathematics

Volume 2013 (2013), Article ID 534891, 6 pages

http://dx.doi.org/10.1155/2013/534891

Research Article

Total Variation Filtered Demons for Improved Registration of Sliding Organs

Damir Demirović,¹ Amira Šerifović-Trbalić,¹ Naser Prljača,¹ and Philippe C. Cattin²

¹Department of Electrical and Computer Engineering, University of Tuzla, 75000 Tuzla, Bosnia and Herzegovina

²Medical Image Analysis Center, University of Basel, 4031 Basel, Switzerland

Received 16 September 2013; Accepted 10 November 2013

Academic Editors: J. Chow and X. Wan

Copyright © 2013 Damir Demirović 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 present a new approach to regularize the displacement field of the accelerated Demons registration algorithm. The accelerated Demons algorithm uses Gaussian smoothing to penalize oscillatory motion in the displacement fields during registration. This regularization approach is often applied and ensures a smooth deformation field. However, when registering images with discontinuities in their motion field such as from organs sliding along the chest wall, the assumption of a smooth deformation field is invalid. In this work, we propose using total variation based smoothing that is known to better retain the discontinuities in the deformation field. The proposed approach is a first step towards automatically recovering breathing induced organ motion with good accuracy.

1. Introduction

Nonrigid image registration methods, Maintz and Viergever [1], play an important role in today’s modern medical diagnostics and Image-Guided Therapy (IGT) systems. They can be broadly split into two main categories, namely, the feature and intensity based approaches. Feature based methods extract corresponding image points as a part of the preprocessing, whereas the intensity based methods make direct use of the image voxels intensities. One of the well-known intensity based methods is optical flow by Horn and Schunk [2]. Because of the ill-posedness of image registration, additional information in a form of regularization rules is necessary. They usually come in form of heuristical rules with underlying physical reasoning.

Recent developments in medical imaging, such as 4DMRI developed by Von Siebenthal et al. [3], allow capturing 4-dimensional (3D+time) breathing organ motion. Today’s registration methods, however, generally are not suitable for registering these data sets, as the heuristical regularization rules produce invalid motion fields across sliding organ boundaries by enforcing smooth deformations fields. A need for properly handling discontinuities is thus a prerequisite when registering these data sets.

Lots of research has been devoted to discontinuity preserving registration in the last couple of years. For example, some studies preserve motion discontinuities by using joint registration/segmentation. Kiriyanthan et al. [4] utilize robust regularization energy formulations. In Ruan et al. [5], authors propose a regularization energy that encodes a discriminative treatment of different type of motion discontinuities, by using Helmholtz-Hodge decomposition. In another approach, Heinrich et al. [6] present a framework using variational formulation of the optical flow model for 3D medical image registration, using robust discontinuity preserving nonquadratic regularization. This allows globally smooth deformations with local discontinuities, such as can be seen in sliding motion of the abdominal organs during respiration.

Ding et al. [7] propose method to evaluate the sliding motion of the lobar surfaces during respiration using lobe-by-lobe mass-preserving nonrigid image registration. They previously segment the lobes in the three parts.

In Schmidt-Richberg et al. [8], diffusion-based model for incorporating physiological knowledge in image registration is presented. They estimate slipping motion at the organ borders by decoupling normal- and tangential-directed smoothing for estimation of respiratory lung motion.

Deformation field regularization term based on anisotropic diffusion that accommodates the deformation field discontinuities during sliding motion is presented in Pace et al. [9]. Forced continual field of movement which results from the today registration methods in general is not valid across borders of an organ (tissue). These discontinuities are formed, for example, in breathing cycle, when a liver moves across an abdomen. Taking into account discontinuities has equally role in registration of pre- and intra-operative images. Today methods of registration in general do not take into account discontinuities and are possibly adequately registered images when, for example, part of a tissue is removed by a surgeon. This constrain leads to poor or unaccepted registrations, and therefore obtained results can become useless. One can use alternative methods which are in general complex and need manual segmentation of tissues, which needs registration of previous segmented parts. Process of combining obtained deformation fields is complex task which takes into account manual intervention, and therefore these methods are semiautomatic.

In this paper, we propose an improvement to the well-known Demons algorithm that is better suited to handle such discontinuities in the deformation field.

2. Methods

The goal of registration is to obtain mapping from one coordinate system of one image to the coordinate system of the other image. In other words, it tries to find a displacement field that warps each pixel of the moving image onto the static reference image in an optimal sense. The optimality criteria are defined by the similarity measure and depend on the nature and modality of the input images. As image registration is an ill-posed problem, it needs further physically inspired constraints which ensure realistic deformations and penalize unwanted oscillations.

2.1. Thirion’s Demons

The Demons algorithm as proposed by Thirion [10] draws an analogy to thermodynamical principles found by Maxwell in the 19th century when describing a thermodynamical paradox. The aim of the Demons algorithm is to find a displacement field that maps all pixels of the moving image onto the corresponding pixels of the static image . The algorithm is iterative and new forces are calculated at each step by evaluating

where , are the respective pixel gray values, is the gradient of the static image, and is the displacement vector. This force is oriented according to the values of the image gradients and the differences between the images. represents the internal force and is the differential force between the static and the current moving image and hence is equivalent to an external force. When converging, the forces get smaller with each iteration of (1). As image registration is an ill-posed problem and the displacement field is calculated only from local information, regularization is essential. Thus the current displacement field requires a regularization after each iteration with a Gaussian filter.

2.2. Accelerated Demons Algorithm

Different extensions introducing additional forces have been proposed to speed up convergence of the original Demon registration algorithm. Wang et al. [11], for example, proposed an additional active force term that is added to the standard passive forces during calculation of the displacement field. To further improve convergence of the algorithm, Wang combined his active force approach with the normalization factor proposed earlier by Cachier et al. [12].

The factor can be interpreted as a step size that is adaptively changed after each iteration. They proposed to continuously reduce the step size until convergence. Equation (2) shows how the adaptive normalization factor is integrated into Wang’s Demon force calculation:

Wang’s integration of the moving image gradient and Cachier’s factor significantly improved convergence in comparison to Thirion’s Demon implementation. Another well-known strategy to reduce computational complexity is the multiresolution approach in which image resolution of the images is increased from being coarse to fine. The multiresolution approach is straightforward and all these enhancements form the basis for our further developments.

2.3. Total Variation Filtering of the Deformation Field

Filtering based on the total variation has been introduced as a method for edge preserving smoothing in image processing by Rudin et al. [13] and is widely known as the ROF model. In a recent paper, Wedel et al. [14], the Total Variation norm has also been applied when calculation the optical flow. Thanks to the TV norm, the proposed method preserves discontinuities in the flow field. A further advantage of using the -norm is the increased robustness against outliers.

The total variation is generally defined as a variational and can be described as follows:

where is the image domain, is the observed image function which is assumed to be corrupted by Gaussian noise, and is the sought solution. The balancing parameter is used to control the level of smoothing in . The aim of the ROF model is to minimize the total variation of with an image similar to the original. The first term of the variational in (3) is called regularizer and can be written as follows:

Here the main property is to allow image discontinuities and still being convex in . In contrast to the Tikhonov model [15], the squared term has been replaced by the norm. This change brings about some improvements with regard to robustness and discontinuity preservation but at the price of the norm not being a convex optimization problem. On the computational side, the TV norm is more challenging to optimize as it is not differentiable at zero. This problem opened a wide field of research in mathematics and many solutions have been proposed to handle these difficulties. To minimize the variational term, we employ the Euler-Lagrange equation with its associated differential equation. For the ROF model, this equation is given by

The main difficulty in this equation is that the term makes it highly nonlinear and the equation is furthermore not defined for . To address the problem at the location 0, we replace the term term by a regularized version . The value of has to be well balanced, as the ability to get discontinuities is limited for large values and conversely with lower , we have problems with stability and convergence.

The problem from (5) can be solved with the solvers such as Gauss-Seidel or Jacobi. In the original paper Rudin et al. [13] proposed to perform an explicit time marching until the steady state is reached for the following PDE:

In fact, explicit time marching is equal to the method of steepest descent which is given by the following scheme:

The ROF model (6) can also be solved by linearizing the non-linear term from the previous iteration and then we get the following iterative process:

If we compare the above equation with the Tikhonov model, we can see the similarities between them. The only difference lies in the fact that the quadratic norm is replaced by the TV norm. Although being a minor change, it leads to some strong improvements of the denoising performance. Unlike the Tikhonov model, this approach is able to remove noise and being robust against outliers while preserving sharp edges.

2.4. Overview of the Proposed Algorithm

Our proposed registration algorithm is based on the accelerated Demons algorithm. Instead of Gaussian smoothing, however, we regularize the deformation field with total variation filter to better recover the discontinuous deformation field. The proposed algorithm can be described in Algorithm 1.

Reset deformation field:

Set to a total of num iterations

Compute for each pixel

for to number of iterations

    for each pixel

        Compute intensity value in deformed image:

         where

        Compute additional displacement vector:

        Compute total displacement vector

    end

    Smooth v with a Total Variation filter

end

Algorithm 1

The Demons method, in general, converges fast and approximately 50 iterations were enough for our experiments. For smoothing the summed deformation fields at each iteration, the total variation method described in Chambolle [16] was used. This method approximates the ROF denoising model (3). The stopping criteria need to be specified, and they are usually enough to set it to 0.01.

3. Results

In first experiment, we evaluate the proposed method with synthetic images (Figures 1(a) and 1(b)) where the ground-truth deformation field is known. Here the floating image was deformed by shifting the right half of the image by five pixels down. Although the left half of the picture should present no motion field, the original Gaussian filtered Demons registration method produces a nonzero motion field, whereas the TV norm filtered Demons correctly provides a close to zero deformation field on the left side. In addition, the motion field along the discontinuity (Figure 1(d)) is much crisper with our proposed approach as compared to the original implementation using a Gaussian filtering stage.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 1: Registration results for synthetic image pair: (a) reference slice and (b) floating slice: second row; deformation field and magnitude after registration with Demons: third row; deformation field and magnitude after registration with TV-Demons.

Ground-truth data in this experiment is known, that is all deformation field vectors are normal to -axis. We can evaluate the performance of the proposed method by calculating the mean and the standard deviation of the deformation field with respect to the ground-truth motion. We use the following simulation parameters: Kernel of size , for Demon and and stopping criteria to be 0.01 for TV Demon. As a result, we get the following angles: and for Demon and TV Demon, respectively.

In the next experiment (Figure 2), we tested our approach on pairs of 2D liver images, taken from a dataset of 4DMRI Von Siebenthal et al. [3] obtained from their website. The temporal resolution of the volumes is 362 ms. Each volume consists of voxel of size  mm³. Eleven registration pairs were extracted from the dataset. These pairs represent the liver in different breathing phases. All image pairs undergo motion which show a discontinuous motion field along the chest wall. For evaluation purposes, we manually create eleven image masks which consist of the region of interest and than apply them to get input images (see example mask on Figure 2(c)). Then we registered each image pair with the original Demon and the proposed TV Demon approach.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 2: Comparison of the Demons and TV Demons registration. (a) reference slice, (b) floating slice, (c) masked out region of the image, (d) Demons deformation field, (e) close-up of Demons deformation field, (f) Magnitude of the calculated Demons deformation field, (g) TV Demons deformation field, (h) close-up of TV Demons deformation field, and (i) Magnitude of the calculated TV Demons deformation field.

We run the simulations for all eleven pairs for Demons and our proposed TV Demons. Both implementations run at 50 iterations with a value of . For the Demons, a Kernel of size and was used for all image pairs. We calculated SSD and CC before and after registration for both algorithms. For simulation purposes, we chose the and stoping criteria to be 0.001. Results are given in Table 1.

Table 1: SSD and CC improvements over Demons after registration.

Pair no.
Demons
TV Demons

SSD
CC
SSD
CC

1
0.0113
0.9748
0.0105
0.9782

2
0.0172
0.9543
0.0152
0.9645

3
0.0169
0.9589
0.0149
0.9683

4
0.0181
0.9568
0.0153
0.9692

5
0.0119
0.9825
0.0116
0.9836

6
0.0172
0.9577
0.0162
0.9625

7
0.0174
0.9616
0.0118
0.9824

8
0.0113
0.9822
0.0102
0.9853

9
0.0282
0.8950
0.0263
0.9086

10
0.0097
0.9775
0.0080
0.9848

11
0.0167
0.9523
0.0133
0.9697

From Table 1, we can see improvements for the TV Demons in both SSD and CC. It can be clearly seen in Figures 2(h) and 2(g) that the TV Demon produces more realistic deformation fields.

To verify that the improvement is significant, we statistically analyzed the data. All of these statistical tests were performed using the software package R (release 2.9.2 for Linux, http://www.r-project.org/). The normal distributed SSD and CC results were compared with a paired -test. Normality of the distributions was tested with the Kolmogorov-Smirnov test. A value was considered statistically significant for all tests. Statistical analysis of the obtained results showed that both registration approaches perform significantly better () than not registering the images at all. Taking the Demons algorithm with Gaussian filtering as the standard, it can be shown that both SSD and CC significantly improved () when using the proposed TV Demon approach. With regard to the computational complexity, the running times of the Demons and proposed TV Demon were 2.96 s and 64.21 s, respectively.

4. Conclusion

The Gaussian smoothing in Thirion’s Demon approach does not provide satisfying results when registering images showing discontinuous motion fields, that is, when registering images showing respiratory organ motion. In this paper, total variation regularization for the deformation field for Demos algorithm is proposed. The proposed algorithm, TV Demon, with total variation filtering as the regularizer significantly improves the results in the tested cases showing discontinuities. We showed its superior performance on synthetic as well as on clinical images.

Conflict of Interests

Authors declare no conflict of interest with the software package R (http://www.r-project.org/).

Acknowledgment

This work has been supported by the SCOPES Project of the Swiss National Science Foundation (http://www.snf.ch/).

References

A. J. B. Maintz and M. A. Viergever, “A survey of medical image registration,” Medical Image Analysis, vol. 2, no. 1, pp. 1–36, 1998.

B. K. P. Horn and B. G. Schunk, “Determining optical flow,” Artifical Inteligence, vol. 17, pp. 1185–3203, 1981.

M. Von Siebenthal, G. Székely, U. Gamper, P. Boesiger, A. Lomax, and P. Cattin, “4D MR imaging of respiratory organ motion and its variability,” Physics in Medicine and Biology, vol. 52, no. 6, article 1, pp. 1547–1564, 2007.

S. Kiriyanthan, K. Fundana, and P. C. Cattin, “Discontinuity preserving registration of abdominal MR images with apparent sliding organ motion,” in Abdominal Imaging. Computational and Clinical Applications, vol. 7029, pp. 231–239, 2012.

D. Ruan, J. A. Fessler, and S. Esedoglu, “Discontinuity preserving regularization for modeling sliding in medical image registration,” in Proceedings of the IEEE Nuclear Science Symposium Conference Record (NSS/MIC '08), pp. 5304–5308, October 2008.

M. P. Heinrich, M. Jenkinson, M. Brady, and J. A. Schnabel, “Discontinuity preserving regularization for variational optical- ow registration using the modified L-p norm,” in Proceedings of the MICCAI Workshop on Evaluation of Methods for Pulmonary Image Registration (EMPIRE '10), 2010.

K. Ding, Y. Yin, K. Cao et al., “Evaluation of lobar biomechanics during respiration using image registration,” in Proceedings of the 12th International Conference, Medical Image Computing and Computer-Assisted Intervention (MICCAI '09), pp. 739–746, London, UK, 2009.

A. Schmidt-Richberg, J. Ehrhardt, R. Werner, and H. Handels, “Slipping objects in image registration: improved motion field estimation with direction-dependent regularization,” in Proceedings of the Medical Image Computing and Computer-Assisted Intervention (MICCAI '09), vol. 5761, no. 1, pp. 755–762, 2009.

D. F. Pace, A. Enquobahrie, H. Yang, S. R. Aylward, and M. Niethammer, “Deformable image registration of sliding organs using anisotropic diffusive regularization,” in Proceedings of the 8th IEEE International Symposium on Biomedical Imaging: From Nano to Macro (ISBI '11), pp. 407–413, usa, April 2011.

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.

H. Wang, L. Dong, J. O'Daniel et al., “Validation of an accelerated “demons” algorithm for deformable image registration in radiation therapy,” Physics in Medicine and Biology, vol. 50, no. 12, pp. 2887–2905, 2005.

P. Cachier, X. Pennec, and N. Ayache, “Fast non-rigid matching by gradient descent: study and improvement of the “Demons” algorithm,” Tech. Rep. RR-3706.

L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Physica D, vol. 60, no. 1–4, pp. 259–268, 1992.

A. Wedel, T. Pock, C. Zach, H. Bischof, and D. Cremers, “An improved algorithm for TV-L1 optical flow,” in Statistical and Geometrical Approaches to Visual Motion Analysis, pp. 23–45, Springer, Berlin, Germany, 2009.

A. N. Tikhonov and V. Y. Arsenin, Solutions of Ill Posed Problems, Edited by F. John, V. H. Winston and Sons, New York, NY, USA, 1977.

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.