The improved differential demon algorithm

Abstract

BACKGROUND:

Differential demon is a fast and efficient registration algorithm. It drives the floating image to deform using the force based on the gradient between the reference and floating image. But it will cause abnormal deformation when the driving force approaches zero,which limits its practical applications.

OBJECTIVE:

This paper proposed an improved differential demon algorithm, which aimed to enhance the registration performance of the existing demon algorithm.

METHODS:

Firstly, we review the original differential demon algorithm. Then, we propose an improved differential demon algorithm and the process of mathematical deduction. Finally, we use experiment to prove that the improved differential demon algorithm is effective and it can improve the accuracy of registration.

RESULTS:

We tested our method on data sets provided by Xuanwu Hospital Capital Medical University. The registration performance proved to be better than the original demon algorithm in terms of mutual information, normalized correlation coefficient, mean square error and iteration number.

CONCLUSIONS:

Experiment results demonstrate the superiority of method proposed in this paper to the original demon algorithm.

1. Introduction

Medical image registration has become one of the most popular topics in the field of medical image analysis. All anatomical points, or at least the points of diagnostic significance and surgical interest of two images can be successfully matched [1].

In the past decade years, various studies of medical imaging registration had been conducted. The previous work on medical image registration can be divided into two categories: one is based on feature points and the other is based on physical model. The former method needs to build the transformation model based on the corresponding feature points extracted on both reference and floating images. Then this model is used to drive the floating image transform toward reference image [2, 3, 4, 5]. However, it’s really hard to extract corresponding feature points accurately, so the registration is not always correct. Normally, algorithm based on physical model needs to build a similarity criteria of statistics between the reference image and the floating image, and then complete the registration with one kind of physical model [6, 7, 8]. Thirion proposed the demon algorithm in 1998 [9]. Demon algorithm is a non-parametric, non-rigid registration method with high accuracy and fast convergence. However, it may cause unreasonable physical deformation or registration error when it deals with larger deformation [10]. In order to solve the problem, Vercauteren et al. proposed the differential demon algorithm [11] and Fritscher et al. used the “anatomically correct” [13] displacements vector fields to replace the tradition displacement field in demon [12].

2. The improved differential demon algorithm

Differential demon algorithm is proposed by Vercauteren et al. [11]. It takes each pixel of the floating image as a genie, and then these genies will drive the floating image to deform to the reference image. Under the assumption that the gray value of the image pixels is constant during the process of the transformation, the driving force comes from the theory of optical flow equation, namely:

$F(x,y,t)=F(x+\mathrm{\Delta }x,y+\mathrm{\Delta }y,t+\mathrm{\Delta }t)$(1)

$p$ is the pixel of the floating image and $M$ is the gray value of $p$ . $F$ is the gray value of the $p$ in reference image. From Eq. (1), we can get Eq. (2)

$f={\displaystyle \frac{(M(p)-F(p))\cdot (\nabla F(p)+\nabla M(p))}{{\parallel (\nabla F(p)+\nabla M(p))\parallel }^2+\frac{{\sigma }_i^2}{{\sigma }_x^2}I}}$(2)

In the Eq. (2),

$$f=(\frac{dx}{dt},\frac{dy}{dt}),\nabla F(p)={(\frac{\partial F}{\partial x},\frac{\partial F}{\partial y})}^T$$

and

$$\nabla M(p)={(\frac{\partial M}{\partial x},\frac{\partial M}{\partial y})}^T,$$

${\sigma }_i$ represents the noise at $p$ and ${\sigma }_x$ represents the uncertainty between the two geometric transformation. When $\nabla F(p)+\nabla M(p)$ is zero, the registration will be erroneous. Therefore, an differential demon algorithm to combine the normalized correlation coefficient is improved in this paper.

In this paper, energy function is proposed as Eq. (3):

$E_t^{\text{𝑐𝑜𝑟𝑟}}(u)={\parallel F-M\circ t\circ \exp (u)\parallel }^2+{\displaystyle \frac{{\sigma }_i^2}{{\sigma }_x^2}}\text{𝑑𝑖𝑠𝑡}{(t,t\circ \exp (u))}^2+{\gamma }^2(1-{\rho }_{F,M\circ t\circ \exp (u))})$(3)

$F$ stands for the reference image, $M$ stands for the floating image, $t$ is the prior deformation field. $u$ is the current deformation field, $\text{𝑑𝑖𝑠𝑡}{(t,t\circ \exp (u))}^2$ means the two adjacent displacement difference , $\gamma$ is the regular factor, ${\rho }_{F,M\circ t\circ \exp (u))}$ represents the normalized correlation between the image $F$ and the image $M\circ t\circ \exp (u)$ . Equation (3) is the differential transformation in the updated displacement field, which can maintain the invariance of topology during the deformation. We can get Eqs (4) and (5) by Taylor expansion:

$f_{1t}^p(u)=F(p)-M\circ t\circ \exp (u)(p)\approx f_{1t}^p(0)+J_1^p\cdot u(p)=F(p)-M\circ t(p)+J_1^p\cdot u(p)$(4) $f_{2t}^p(u)=\gamma \cdot (1-{\rho }_{F(p),M\circ t\circ \exp (u)(p))})\approx f_{2t}^p(0)+J_2^p\cdot u(p)=\gamma \cdot (1-{\rho }_{F(p),M\circ t(p)})+J_2^p\cdot u(p)$(5) $E_t^{corr}(u)={\displaystyle \frac{1}{2|\mathrm{\Omega }|}}{\displaystyle \sum _{p\in \mathrm{\Omega }}}{||\left(\begin{array}{c}\hfill F(p)-M\circ t(p)\hfill \\ \hfill 0\hfill \\ \hfill \gamma (1-{\rho }_{F(p),M\circ t(p)})I\hfill \end{array}\right)+\left(\begin{array}{c}\hfill J{}_{1p}\hfill \\ \hfill {\displaystyle \frac{{\sigma }_i(p)}{{\sigma }_x}}I\hfill \\ \hfill \gamma \cdot J_{2p}\hfill \end{array}\right)u(p)||}^2$(6)

From Eqs (4) and (5), we can get Eq. (6). Minimizing Eq. (6), we can get Eq. (7):

$u(p)={\displaystyle \frac{(F(p)-M\circ t(p))J_{1p}^T+r^{\wedge }2\cdot (1-{\rho }_{F,M\circ t})\cdot J_{2p}^T}{{\parallel J_{1p}^T\parallel }^2+\frac{{\sigma }_i^2}{{\sigma }_x^2}I+r^{\wedge }2\cdot {\parallel J_{2p}^T\parallel }^2}}$(7)

In the following, ESM is used to optimize the algorithm. In Eq. (7),

$$J_{1p}(u)={\frac{\partial }{\partial v^T}f(t\circ \exp (v))|}_{v=u}$$

and Taylor expansion is as Eq. (8), we join Eq. (8) to get Eq. (9).

$J_{1p}(u)=J_{1p}(0)+u^T\cdot H_{1p}+O({\parallel u\parallel }^2)$(8) $f(t\circ \exp (u))=f(t)+{\displaystyle \frac{1}{2}}(J_{1p}(u)+J_{1p})\cdot u+O({\parallel u\parallel }^3)$(9)

However, it is difficult to minimize the closed form of $J_{1p}(u)$ . Because of the particularity of the registration task, we can assume that $u$ is already optimized. The updating floating image is Eq. (10):

$M\circ t^{opt}=M\circ t\circ \exp (u^{opt})$(10) $J_{1p}\triangleq -{\displaystyle \frac{1}{2}}(\nabla F+\nabla M\circ t)$(11)

Figure 1.

(a)–(h) are 1–8 sets of data respectively. The first row is the reference image. The second row is floating image. The third row is the result image using this improved differential Demons algorithm. The fourth row is the result image using original algorithm of differential Demons. The red dashed represents the difference of the same parts.

Similarly,

$J_{2p}=\nabla F+I\ast {\rho }_{\nabla F,\nabla M\circ t}$(12)

Based on Eqs (11) and (12), we can get $u(p)$ .

$u(p)={\displaystyle \frac{(F(p)-M\circ t(p))J_p^T+r^{\wedge }2\ast (1-{\rho }_{F,M\circ t})\ast (\nabla F+I\ast {\rho }_{\nabla F,\nabla M\circ t})}{{\parallel J_p\parallel }^2+\frac{{\sigma }_i^2}{{\sigma }_x^2}I+r^{\wedge }2\ast {\parallel \nabla F+I\ast {\rho }_{\nabla F,\nabla M\circ t}\parallel }^2}}$(13)

Figure 2.

(a) is the reference image’s coronal plane. (b) is the floating image’s coronal plane. (c) is the result’s coronal plane using this improved differential Demons algorithm. (d) is the result’s coronal plane using original differential Demons algorithm. (e) is the result’s three dimensional figure using this improved differential Demons algorithm. (f) is the result’s three dimensional figure using original differential Demons algorithm. The red dashed represents the difference of the same parts.

The steps of the improved differential demon algorithm are listed as follow:

• Step1.

Choose the initial geometric transformation $t$ . The dense displacement $t$ is identical to zero vector fields;

Step2.

The current transform $t$ is given to obtain $u$ by minimizing the updated dense velocity field

$$E_t^{corr}(u),E_t^{corr}(u)={\parallel F-M\circ t\circ \exp (u)\parallel }^2+\frac{{\sigma }_i^2}{{\sigma }_x^2}dist{(t,t\circ \exp (u))}^2+{\gamma }^2{(1-{\rho }_{F,M\circ t\circ \exp (u))})}^2;$$

Step3.

For a fluid-like regularization, let $u$ $\leftarrow$ Fluid $\star u$ . The convolution kernel will be a typical Gaussian kernel.

Step4.

Update the transform $t$ by using the index mapping and compound operation $t\stackrel{}{⟵}t\circ \exp (u)$ ;

Step5.

Check the convergence. If not, return Step2; else, output the optimal transformation $t^{opt}$ .

Figure 3.

(a) is the comparison of mutual information. (b) is the comparison of normalized correlation coefficient. (c) is the comparison of mean square error. (d) is the comparison of iteration number.

3. Experiments

The experimental data sets are provided by Xuanwu Hospital Capital Medical University. The experimental machine has the CPU core i5-3570, 4 GB memory and works on the 64bits Win 7 operating system. In order to demonstrate that the proposed method is practicable, we choose eight groups of MRI brain data which size is 512*512*192.

In this paper, three evaluation indexes – mutual information (MI), normalized correlation coefficient (CC) and mean square error (MSE) are chosen to evaluate the image registration result. Firstly, Yang’s method is adopted to remove skull [14], and the multi-resolution strategy is used to improve the speed of convergence. The comparison results of our proposed method and the original demon algorithm on 8 different data sets are shown in Fig. 1, it is obviously observed that our method can effectively keep the texture and intensity details of the brain tissue. The corresponding 3D registration results for the proposed method and original demon algorithm are shown in Fig. 2. In order to evaluate the effectiveness and robustness of the proposed algorithm, we compute three different indexs (MI, CC and MSE) for traditional demon method and our proposed method, the comparison results are shown in Figs 3a–c, Fig. 3d shows the iteration comparison result. It can be observed that the greater mutual information and normalized correlation coefficient of two images are, and the smaller the mean square error of two images is, the more similar two images are. From Fig. 3, we can conclude that the precision of the improved diffeomorphic demon algorithm is greater than the original algorithm, and from Fig. 3d, the speed of the improved algorithm is faster than the original algorithm.

4. Conclusion

Differential demon algorithm is a fast and efficient registration algorithm. By using the index mapping and compound operations, it can keep the topology of the image. Original differential demon algorithm adopts the gradient of reference image and floating image as a driving force to drive floating image to deform. However, when driving force is close to zero, there may be an error (for example, gray value of the corresponding pixels in the reference image and floating image are different, but the driving force is zero). To solve the problem, we present an improved differential demon method in conjunction with normalized relations. Experiments demonstrate that the proposed method can achieve more accurate registration result than the original method, especially for the large deformation.

Conflict of interest

None to report.