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.
| ||||||||||||||||||||||||||
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)
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 mm3. 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)
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.
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/).