Nonrigid Medical Image Registration Using an Information Theoretic Measure Based on Arimoto Entropy with Gradient Distributions

Abstract

This paper introduces a new nonrigid registration approach for medical images applying an information theoretic measure based on Arimoto entropy with gradient distributions. A normalized dissimilarity measure based on Arimoto entropy is presented, which is employed to measure the independence between two images. In addition, a regularization term is integrated into the cost function to obtain the smooth elastic deformation. To take the spatial information between voxels into account, the distance of gradient distributions is constructed. The goal of nonrigid alignment is to find the optimal solution of a cost function including a dissimilarity measure, a regularization term, and a distance term between the gradient distributions of two images to be registered, which would achieve a minimum value when two misaligned images are perfectly registered using limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) optimization scheme. To evaluate the test results of our presented algorithm in non-rigid medical image registration, experiments on simulated three-dimension (3D) brain magnetic resonance imaging (MR) images, real 3D thoracic computed tomography (CT) volumes and 3D cardiac CT volumes were carried out on elastix package. Comparison studies including mutual information (MI) and the approach without considering spatial information were conducted. These results demonstrate a slight improvement in accuracy of non-rigid registration.

1. Introduction

Volume registration is an essential task of image processing, especially in medical field, such as aiding diagnosis, surgical applications, and image-guided radiation therapy [1,2]. The images to be registered are generally obtained at different times and from different imaging sensors, namely, multi-modality imaging. Different medical imaging modalities could provide various and complementary information. For instance, CT and MRI display the anatomic structures of an organ, while positron emission tomography (PET) and single photon emission computed tomography (SPECT) provide the functional and metabolic information. Therefore, in clinic, these multi-modal volumes are often registered and fused together, in this way, much complementary information derived from different modalities are supplied with the physicians to improve the diagnosis accuracy and assessment efficiency of lesion progression.

Image registration is related by the process to find the optimal mapping function between two images to be aligned [3]. In recent years, the registration algorithms using the similarity based on information theory have been attracted more and more attention in medical image registration, among which maximization of MI was early reported for registration of medical images from different modalities by Collignon et al. [4], Maes et al. [5], Wells et al. [6]. Studholme et al. [7] studied a normalized similarity measure called normalized MI (NMI) to tackle the problem of changing field of view (FOV). MI and NMI are both estimated by the probability distributions of images to be registered. Besides, the concept of cumulative probability distributions was introduced into image registration, and cumulative residual entropy (CRE) was investigated in [8]. Additionally, the relations between CRE and Shannon entropy were ulteriorly researched in [9], and a comparison with MI was reported in [10]. In this new measure, cumulative density functions (CDF) instead of probability density functions (PDF) were adopted to calculate values of the similarity measure, which illustrates a good robustness to noise.

The aforementioned measures based on information theory—such as MI, NMI, and CRE—are constructed by Shannon entropy. Nonetheless, the additivity of Shannon entropy signifies that it is an extensive entropy. However, Antolin et al. [11] pointed out that the extensive entropy does not consider the correlation of two variables. Consequently, they presented a similarity measure exploiting Tsallis entropy. Subsequently, this new divergence was employed to construct a non-rigid registration model for medical images [12,13].

Illuminated by the reference [11], a divergence measure based on Arimoto entropy was presented called the Jensen–Arimoto divergence (JAD) [14]. In [15], some properties, such as the concavity of Arimoto entropy and the boundedness of JAD, have been further investigated. However, the registration method based on JAD does not take the spatial information between voxels into account, which is of significance to medical image registration. This paper aims to present a novel nonrigid registration method of medical images adopting a normalized measure based on Arimoto entropy and gradient distributions. Firstly, the properties of Arimoto entropy and JAD are analyzed, and a distance of gradient distributions is constructed. Secondly, a nonrigid deformation model is chosen and the registration process is formulated by an optimization procedure. In the sequel, the continuous probability distributions are estimated using Parzen window method applying B-splines and the analytical gradient of objective function can be obtained. Finally, the L-BFGS optimization [16] is adopted to obtain the optimal deformation parameters. To assess the performance of our registration framework for medical images, several groups of non-rigid experiments on simulated volumes and real 3D data are implemented.

Our contributions are twofold. Firstly, the related measures based on Shannon entropy do not consider the correlation. Therefore, we present a normalized measure based on Arimoto entropy, a non-extensive entropy, as the dissimilarity measure. Secondly, in the existing measures, such as MI, NMI and JAD, the intensity values are directly exploited to calculate the similarity measure, while the spatial information has been not considered. To take the spatial information between voxels into account, a distance term of gradient distributions is constructed and incorporated into the objective function.

The rest of this work is arranged as follows. In Section 2, the knowledge of information theory was firstly reviewed, and then introduce the Arimoto entropy and JAD measure, constructing gradient distribution distance. We formulate the model of nonrigid registration and give the detailed description of the registration method in Section 3 adopting a normalized measure based on Arimoto entropy with gradient distribution. Section 4 demonstrates the nonrigid test results on 3D MR volumes and real 3D clinical datasets, with the compared results to other registration algorithms illustrated. Finally, we provide the conclusions and perspectives in Section 5.

2. Preliminaries

For this section, we briefly review the theoretical concept of information theory, and then introduce the Arimoto entropy and JAD, along with studying their properties. In addition, the gradient distributions of reference image and float image are constructed and a distance of them is derived.

2.1. Shannon Entropy and Mutual Information

For an arbitrary random variable X(x₁, x₂,…, x_N), with its probability distributions p(x₁, x₂,…, x_N), the Shannon entropy of X is used to measure the amount of average uncertainty included in this random variable,

$$H(X)=-{\displaystyle \sum _{i=1}^{N}p(x_i)\log p(x_i)}$$

(1)

which is also employed to measure the account of information provided by this random variable. Then, considering another random variable Y, H(X|Y) is remarked as the conditional entropy of X when Y is known. The reduction in uncertainty due to Y is called the MI. MI of X and Y is defined by

$$I(X,Y)=H\left(X\right)-H\left(X|Y\right)={\displaystyle \sum _x{\displaystyle \sum _{y}p(x,y)\log \frac{p(x,y)}{p(x)\cdot p(y)}}}$$

(2)

where p(x) and p(y) denote marginal probability of two random variables, as well as p(x, y) being the joint distribution of them. MI is applied as a measure of the dependence between X and Y. MI is symmetric in X and Y and always nonnegative [17].

2.2. Arimoto Entropy

Arimoto [18] introduced a generalized form of Shannon entropy, Arimoto entropy is defined by

$$A_{\alpha }(X)=\frac{\alpha }{\alpha -1}\left[1-{({\displaystyle \sum _{i=1}^N{p}_i^{\alpha }})}^{\frac{1}{\alpha }}\right]\alpha >0,\alpha \ne 1$$

(3)

Boekee et al. [19] investigated some significant properties of Arimoto entropy, here, we only exhibit several useful properties as follows.

Non-negativity:

$$A_{\alpha }(X)\ge 0\alpha >0,\alpha \ne 1$$

(4)

Pseudo-additivity:

$$A_{\alpha }(X,Y)=A_{\alpha }(X)+A_{\alpha }(Y)-\frac{\alpha -1}{\alpha }{A}_{\alpha }(X)A_{\alpha }(Y)\alpha >0,\alpha \ne 1$$

(5)

Concavity:

$$A_{\alpha }(t{X}_1+(1-t)X_2)\ge t{A}_{\alpha }(X_1)+(1-t)A_{\alpha }(X_2)t\in (0,1),\alpha >0,\alpha \ne 1$$

(6)

Symmetry:

$$A_{\alpha }\left(\cdots ,p_i,\cdots p_j\cdots \right)=A_{\alpha }\left(\cdots ,p_j,\cdots p_i\cdots \right)$$

(7)

Upper bound:

$$A_{\alpha }\left(p_1,p_2,\cdots ,p_N\right)\le A_{\alpha }\left(\frac{1}{N},\frac{1}{N},\cdots ,\frac{1}{N}\right)=\frac{\alpha }{\alpha -1}\left[1-N^{\frac{1-\alpha }{\alpha }}\right]$$

(8)

See [19] for the detailed proof of these properties. Property one ensures that Arimoto entropy is non-negative. Its pseudo-addivity illustrates that Arimoto entropy accounts for a non-extensive entropy. The parameter α in (3) accounts for the degree of non-extensivity.

2.3. Jensen Arimoto Divergence

For a random variable X, and $P({\mathit{p}}_1,{\mathit{p}}_2,\dots ,{\mathit{p}}_N)$ is probability distributions on X. JAD is defined to be [15]

$$J{A}_{\alpha }({\mathit{p}}_1,{\mathit{p}}_2,\dots ,{\mathit{p}}_N)=A_{\alpha }\left({\displaystyle \sum _{i=1}^N{\omega }_i{\mathit{p}}_i}\right)-{\displaystyle \sum _{i=1}^N{\omega }_i{A}_{\alpha }({\mathit{p}}_i)}\alpha >0,\alpha \ne 1$$

(9)

where $A_{\alpha }(\cdot )$ denotes Arimoto entropy and ${\omega }_i$ is a weighted vector to constrain ${\omega }_i\ge 0$ and ${\displaystyle {\sum }_{i=1}^N{\omega }_i}=1$. In the following, we review the properties of JAD.

Proposition 1.

JAD has these properties of non-negativity, symmetry. Also, JAD is identical to 0 when and only when all of the probability distributions are the same as each other. The proof has been reported in [15].

Li et al. pointed that a distance must fulfill four requirements [20]. JAD does not meet the triangle inequality, so JAD is not a true distance metric. Nonetheless, JAD can still be adopted to measure the disparity among the probability distributions of two random variables.

Proposition 2.

The JA divergence has the maximum when p₁, p₂,…, p_N are degenerate distributions, where ${\mathit{p}}_i={\delta }_{ij}=1$ when i = j and 0 otherwise. The proof had been provided in [21], and the maximum equals to $A_{\alpha }(\omega )$.

According to Proposition 1 and 2, it is obviously observed that JAD is bounded, $0\le J{A}_{\alpha }({\mathit{p}}_1,{\mathit{p}}_2,\dots ,{\mathit{p}}_N)\le A_{\alpha }(\omega )$.

2.4. Gradient Distributions Distance

Given the reference image R and float image F, ∇F, and ∇R represent the gradients of F and R, along with q(∇F) and p(∇R) representing the gradient distributions of F and R. Kullback–Leibler divergence (KLD) can be used to calculate the distance of q(∇F) and p(∇R) as

$$KLD\left(q\left|\right|p\right)={\displaystyle \sum _{\mathrm{x}}q\left(\nabla F\left(\mathrm{x}\right)\right)\log \frac{q\left(\nabla F\left(\mathrm{x}\right)\right)}{p\left(\nabla R\left(\mathrm{x}\right)\right)}}$$

(10)

where x denotes any point of image gradient, x = [x, y, z]^T. KLD is also known as relative entropy, measuring the diversity of two probability distributions. In (10), we use the convention that 0·log(0/0) = 0. In other word, when R and F are completely registered, KLD between gradient distributions defined in (10) achieves the minimum value. Considering the spatial transformation T_μ, μ denoted by the parameters of transformation model, the gradient distribution distance of F and the transformed F is given by

$$KLD\left(\nabla F\left(T_{\mu }\left(\mathrm{x}\right)\right)\left|\right|\nabla R\left(\mathrm{x}\right)\right)={\displaystyle \sum _d{\displaystyle \sum _{\mathrm{x}}q\left(\nabla F_d\left(T_{\mu }\left(\mathrm{x}\right)\right)\right)\log \frac{q\left(\nabla F_d\left(T_{\mu }\left(\mathrm{x}\right)\right)\right)}{p\left(\nabla R_d\left(\mathrm{x}\right)\right)}}}$$

(11)

where ∇F_d(T_μ(x)) and ∇R_d(x) represent the gradients of transformed float image and reference image, and d is the dimension of the images to be registered. Subsequently, the Parzen method is employed to estimate the gradient distributions.

Denote β⁽⁰⁾ and β⁽³⁾ by zero-order and three-order B-spline, respectively. In the process of image registration, the transformation parameters μ does not affect the reference image, the gradient of reference image is also constant in registration process. Consequently, the gradient of reference image can be calculated before registration to improve the implementation efficiency. In addition, a zero-order B-spline is exploited to estimate the gradient distributions ∇R_d(x) of reference image R(x), with the probability density function of ∇R_d(x) defined by

$$p\left(\nabla R_d\left(\mathrm{x}\right)\right)={\tilde{p}}_d(r_i)=\frac{1}{V}{\displaystyle \sum _{\mathrm{x}\in \mathsf{\Omega }}{\beta }^{(0)}\left(r_i-\frac{\nabla R_d(\mathrm{x})-\nabla R_d^0}{\Delta b_R}\right)}$$

(12)

where Ω is volume domain to estimate probabilities. V denotes the number of voxels of Ω domain, as well as d being the image dimension. r_i represents intensity levels of ∇R_d(x), and ‘∇’ is the gradient operator. The three-order B-spline is adopted to compute the gradient distributions ∇F_d(T_μ(x)) of transformed float image, with the probability density function of ∇F_d(T_μ(x)) shown as

$$q\left(\nabla F_d\left(T_{\mu }\left(\mathrm{x}\right)\right)\right)={\tilde{q}}_d(f_j)=\frac{1}{V}{\displaystyle \sum _{x\in \mathsf{\Omega }}{\beta }^{(3)}\left(f_j-\frac{\nabla F_d(T_{\mu }(\mathrm{x}))-\nabla F_d^0}{\Delta b_F}\right)}$$

(13)

In Equations (12) and (13), Δb_R and Δb_F are the widths of bins. $\nabla R_d^0$ and $\nabla F_d^0$ is the minimum values in two images, and f_j represents intensity levels of ∇F_d(T_μ(x)). Substituting (12) and (13) into (11), we can obtain the following formula.

$$\begin{array}{ll}KLD\left(\nabla F\left(T_{\mu }\left(\mathrm{x}\right)\right)\left|\right|\nabla R\left(\mathrm{x}\right)\right)=& {\displaystyle \sum _d{\displaystyle \sum _{\mathrm{x}}{\tilde{q}}_d(f_j)\log \frac{{\tilde{q}}_d(f_j)}{{\tilde{p}}_d(r_i)}}}\\ & ={\displaystyle \sum _d\left[H\left({\tilde{q}}_d\right)-{\displaystyle \sum _{\mathrm{x}}{\tilde{q}}_d(f_j)\log {\tilde{p}}_d(r_i)}\right]}\end{array}$$

(14)

where H is the Shannon entropy. In this paper, KLD of gradient distributions will be applied as a distance term in medical image registration, regularizing that the gradient distribution of float image is similar to the gradient distribution of reference image.

3. Description of Proposed Nonrigid Registration Method

The details of the registration method using a normalized information theoretic measure based on Arimoto entropy with gradient distributions is described. Firstly, an appropriate transformation model needs to be selected, where the registration criteria (objective function or cost function) and optimization algorithm is applied to optimize the criteria. Finally, the images to be registered are aligned using an optimal solution obtained by the optimization scheme. Figure 1 displayed the block diagram of our nonrigid registration framework.

3.1. Formulation

Figure 2a,b depict the corresponding planes in two 3D MR images, which account for T1-weighted MR and T2-weighted MR. Non-rigid registration is formulated as the process of searching for the optimal spatial deformation function of reference image R and float image F as

$$\mathrm{y}=g\left(\mathrm{x};\mu \right)$$

(15)

where g(x; μ) is the deformation function, μ denotes the vector of deformation parameters, with x and y being the coordinates of arbitrary point in R and F, respectively. We can formulate the non-rigid registration of F to R as

$$T^{*}=\underset{\mu }{\arg \min }D\left(F\left(\mathrm{x}\right)\circ g\left(\mathrm{x};\mu \right),R\left(\mathrm{x}\right)\right)=\underset{\mu }{\arg \min }D\left(F\left(g\left(\mathrm{x};\mu \right)\right),R\left(\mathrm{x}\right)\right)$$

(16)

where D represents a dissimilarity measure that can achieves its minimum in registration of R(x) and F(g(x; μ)).

However, the process of image registration is an ill-posed issue, and a penalty term need to be incorporated to obtain a smooth transformation. Considering the regularization term, the cost function E is expressed by

$$E=D\left(F\left(g\left(\mathrm{x};\mu \right)\right),R\left(\mathrm{x}\right)\right)+\lambda S\left(g\left(\mathrm{x};\mu \right)\right)$$

(17)

3.2. Transformation Model

Clinically, there exist some large deformations between medical images. Therefore, a nonrigid transformation is generally employed to deal with organ or tissue deformation. Figure 2c,d display the deformation fields and vectors of the nonrigid transformation between Figure 2a,b. The free-form deformations (FFD) model can preferably described local deformations between medical images. Therefore, cubic B-splines are employed to constructed FFD model [22] to simulate this elastic deformation.

In a 3D image, Φ is a mesh with the size of n_x × n_y × n_z, as well as the control points [ω_i, ω_j, ω_k]^T, and δ is the size of spacing. Then, 3D deformation function of x = [x, y, z]^T could be defined as

$$g\left(\mathrm{x};\mu \right)={\displaystyle \sum _{ijk}{\mu }_{ijk}{\beta }^{(3)}\left(\frac{x-{\omega }_i}{\delta }\right){\beta }^{(3)}\left(\frac{y-{\omega }_j}{\delta }\right){\beta }^{(3)}\left(\frac{z-{\omega }_k}{\delta }\right)}$$

(18)

where μ_ijk is the vector of deformation coefficients, and β⁽³⁾ is the third-order B spine function.

3.3. Registration Criteria

The dissimilarity measure D is the most significant part of objective functions, which is used to measure the difference between two images to be registered. When D achieve the minimum value, the similarity of two volumes is maximized, resulting in the two volumes completely registered. We have employed JAD as similarity measure to register medical images in the presence of non-rigid transformation, in which a negative sign is assigned to the JAD to construct the dissimilarity measure [15]. In this paper, according to Proposition 2 in Section 2.3, a normalized dissimilarity measure based on JAD is introduced. Its definition is given as

$$D\left(F(g(\mathrm{x};\mu )),R(\mathrm{x})\right)=1-\frac{J{A}_{\alpha }\left(F(g(\mathrm{x};\mu )),R(\mathrm{x})\right)}{A_{\alpha }\left(\omega \right)}$$

(19)

In [15], JAD is expressed by

$$\begin{array}{ll}J{A}_{\alpha }\left(F(g(\mathrm{x};\mu )),R(\mathrm{x})\right)& =\frac{\alpha }{1-\alpha }\{{\left[{{\displaystyle \sum _{j=1}^M\left[{\displaystyle \sum _{i=1}^{M}p(r_i)p(f_j|r_i)}\right]}}^{\alpha }\right]}^{\frac{1}{\alpha }}-{\displaystyle \sum _{i=1}^{M}p(r_i)}{\left[{\displaystyle \sum _{j=1}^{M}p{(f_j|r_i)}^{\alpha }}\right]}^{\frac{1}{\alpha }}\}\\ & =\frac{\alpha }{1-\alpha }\left\{{\left[{\displaystyle \sum _{j=1}^{M}p{\left(f_j\right)}^{\alpha }}\right]}^{\frac{1}{\alpha }}-{\displaystyle \sum _{i=1}^M{\left[{\displaystyle \sum _{j=1}^{M}p{(r_i,f_j)}^{\alpha }}\right]}^{\frac{1}{\alpha }}}\right\}\end{array}$$

(20)

where f = (f₁, f₂, …, f_M) and r = (r₁, r₂, …, r_M) are the intensity values in F(g(x; μ)) and R(x). Also, M is bins number. Also, p(f_j|r_i) is the conditional probability. JAD and A_α(ω) are substituted into (19). In consequence, the dissimilarity measure D in (19) is rewritten as

$$D\left(F(g(\mathrm{x};\mu )),R(\mathrm{x})\right)=1-\left\{{\displaystyle \sum _{i=1}^M{\left[{\displaystyle \sum _{j=1}^{M}p{(r_i,f_j)}^{\alpha }}\right]}^{\frac{1}{\alpha }}-{\left[{\displaystyle \sum _{j=1}^{M}p{\left(f_j\right)}^{\alpha }}\right]}^{\frac{1}{\alpha }}}\right\}/\left(1-M^{\frac{1-\alpha }{\alpha }}\right)$$

(21)

where M is the number of bins. To address the problem of nonrigid registration, the smooth deformation needs to be acquired by regularizing the deformation model. Incorporated with the regularization term, the objective function E is rewritten by

$$E(F(g(\mathrm{x};\mu )),R(\mathrm{x}))=D\left(F(g(\mathrm{x};\mu )),R(\mathrm{x})\right)+\lambda S(g(\mathrm{x};\mu ))$$

(22)

where D is the dissimilarity measure defined in (21), and S is the regularization term, with its expression given as follows.

$$\begin{array}{ll}S(g\left(\mathrm{x};\mu \right))& =\frac{1}{V}{\displaystyle {\int }_0^X{\displaystyle {\int }_0^Y{\displaystyle {\int }_0^Z[{\left(\frac{{\partial }^{2}g\left(\mathrm{x};\mu \right)}{\partial x^2}\right)}^2+{\left(\frac{{\partial }^{2}g\left(\mathrm{x};\mu \right)}{\partial y^2}\right)}^2+{\left(\frac{{\partial }^{2}g\left(\mathrm{x};\mu \right)}{\partial z^2}\right)}^2}}}\\ & +2{\left(\frac{{\partial }^{2}g\left(\mathrm{x};\mu \right)}{\partial x\partial y}\right)}^2+2{\left(\frac{{\partial }^{2}g\left(\mathrm{x};\mu \right)}{\partial x\partial z}\right)}^2+2{\left(\frac{{\partial }^{2}g\left(\mathrm{x};\mu \right)}{\partial y\partial z}\right)}^2]dxdydz,\end{array}$$

(23)

However, objective function shown in (22) does not take into account the spatial information between voxels. To deal with the issue, the distance between two gradient distributions q(∇F(g(x; μ))) and p(∇R(x)) displayed in (14) is introduced to (22). As a result, the nonrigid registration process is expressed by

$$\begin{array}{ll}{\mu }^{*}& =\underset{\mu }{\arg \min }E\left(F(g(\mathrm{x};\mu )),R(\mathrm{x})\right)\\ & =\underset{\mu }{\arg \min }\{D\left(F(g(\mathrm{x};\mu )),R(\mathrm{x})\right)+{\lambda }_{1}S(g(\mathrm{x};\mu ))+{\lambda }_{2}KLD\left(\nabla F\left(g\left(\mathrm{x};\mu \right)\right)\left|\right|\nabla R\left(\mathrm{x}\right)\right)\}\end{array}$$

(24)

where KLD represents gradient distribution distance, as well as λ₁ and λ₂ being weight parameters, balancing the tradeoff among a dissimilarity measure D, a regularization S, and a distance term KLD.

3.4. Optimization

Newton–Raphson algorithm has been widely exploited, in which second-order derivatives can show better convergence [23] compared with these strategies based on first-order gradient. L-BFGS [24] does not calculate second-order information. Thus, a high computation efficiency can be achieved. A second-order Taylor approximation [16] of E with respect to μ is given as

$$E(\mu +\Delta \mu )\approx E(\mu )+\Delta {\mu }^T\cdot \nabla E(\mu )+\frac{1}{2}\Delta {\mu }^T\cdot {\nabla }^{2}E(\mu )\cdot \Delta \mu ,$$

(25)

where Δμ is the increment of μ, ∇ is gradient operation. The deformation parameter μ of the L-BFGS optimization algorithm is updated as

$${\mu }^{(k+1)}={\mu }^{(k)}-{(H^{(k)})}^{-1}\cdot \nabla E({\mu }^{(k)}),$$

(26)

In the sequel, the derivative of objective function E with respect to μ need to computed.

$$\frac{\partial E}{\partial \mu }=\left[\frac{\partial E}{\partial {\mu }_1},\frac{\partial E}{\partial {\mu }_2},\cdot \cdot \cdot ,\frac{\partial E}{\partial {\mu }_n}\right],$$

(27)

The pseudo code of our registration approach is displayed in Algorithm 1.

To solve the optimization process, we need calculate the analytical gradient of the objective function E. Traditionally, the probability distributions expressed in (21) was not continuous. Hence, the continuous probability density function (pdf) needs to be estimated by Parzen-window method. The continuous marginal and joint pdfs of two images to be registered have been calculated [15]. The continuous expression of gradient distribution distance has been also provided in Section 2.4. Equalization (18) is substituted into (23), the continuity of smoothness term is acquired. Consequently, the objective function E is continuous and its analytical derivative with respect to μ can be calculated.

Algorithm 1. Nonrigid medical image registration with gradient distributions

Input: Reference image R, floating image F Output: Optimal deformation parameters μ* Set λ1, λ2, NMAX, α, M, N, δ, ε Compute the gradient of R, denote as ∇R(x) and gradient distributions p(∇R(x)) Initialize deformation parameters μ(0), iteration k = 0, F(g(x; μ(0))) = F, E(μ(0)) = 0 While |E(μ(k + 1)) − E(μ(k))|> threshold ε or k < =NMAX  Obtain the deformed float image F(g(x; μ(k+1))) and the regularization S(g(x, μ(k + 1)))  Compute ∇F(g(x; μ(k + 1))) and gradient distributions q(∇F(g(x; μ(k + 1))))  Estimate the dissimilarity measure D and gradient distributions distance KLD  Calculate objective function E(μ(k + 1)) = D(R(x), F(g(x; μ(k + 1)))) + KLD(q(k + 1)||p) + S(g(x, μ(k + 1)))  μ(k + 1) =μ(k) − (H(k)) −1·∇E(μ(k))  k = k + 1 end

Derivative of the Objective Function

The objective function defined in (24) includes dissimilarity measure D, a regularization term S, and a distance term KLD. The derivative of D is deduced as

$$\frac{d\left[D\left(F(g(\mathrm{x};\mu )),R(\mathrm{x})\right)\right]}{d\mu }=-\frac{1}{A_{\alpha }\left(\omega \right)}\frac{d\left[J{A}_{\alpha }\left(F(g(\mathrm{x};\mu )),R(\mathrm{x})\right)\right]}{d\mu }$$

(28)

According to (20), we obtain

$$\frac{d\left[D\left(F(g(\mathrm{x};\mu )),R(\mathrm{x})\right)\right]}{d\mu }=\frac{1}{\left(1-M^{\frac{1-\alpha }{\alpha }}\right)}\cdot \left\{{\displaystyle \sum _i{\displaystyle \sum _{j}Y}}\frac{\partial \tilde{p}(r_i,f_j)}{\partial \mu }\right\}$$

(29)

$$Y={\left({\displaystyle \sum _j\tilde{p}{(f_j)}^{\alpha }}\right)}^{\frac{1}{\alpha }-1}\tilde{p}{(f_j)}^{\alpha -1}-{\left({\displaystyle \sum _j\tilde{p}{(f_j|r_i)}^{\alpha }}\right)}^{\frac{1}{\alpha }-1}\tilde{p}{(f_j|r_i)}^{\alpha -1}$$

(30)

where $\partial \tilde{p}(f_j,r_i)/\partial \mu$ represents derivative of estimated joint probability. The derivative of $\widetilde{p(f_j,r_i;\mu )}$ is calculated by

$$\begin{array}{ll}\frac{\partial \widetilde{p(f_j,r_i;\mu )}}{\partial \mu }=& -\frac{1}{N\cdot \Delta b_F}{\displaystyle \sum _{x\in \mathsf{\Omega }}{\beta }^{(0)}\left(r_i-\frac{R(x)-R^0}{\Delta b_R}\right)}\\ & \times {{\beta }^{\prime }}^{(3)}\left(f_j-\frac{F(g\left(x;\mu \right))-F^0}{\Delta b_F}\right)\times \left(\frac{\partial F(s)}{\partial s}{|}_{s=g\left(x;\mu \right)}\right)\times \frac{\partial (g\left(x;\mu \right))}{\partial \mu }\end{array}$$

(31)

with β⁽⁰⁾ and β′⁽³⁾ being the zero-order B-Splines and the derivative of the three-order B-Splines, respectively. R⁰ and F⁰ are the minimal intensities in R(x) and F(g(x; μ)), as well as $\partial F(t)/\partial t$ being the gradient of the deformed float image F(g(x; μ)), $\partial (g\left(x;\mu \right))/\partial \mu$ can be estimated by FFD model.

To obtain the derivative of the penalty term S, we rewrite (23) as

$$S\left(g(\mathrm{x};\mu )\right)=\frac{1}{V}{{\displaystyle \sum _{\mathrm{x}}{\displaystyle \sum _{i,j=1}^3\left(\frac{{\partial }^{2}g(\mathrm{x};\mu )}{\partial x_i\partial x_j}\right)}}}^2$$

(32)

where x represents the points in image region, and V denotes the number of pixels. The derivative of S has been provided by Staring and Klein [25],

$$\frac{\partial S\left(g(\mathrm{x};\mu )\right)}{\partial \mu }=\frac{1}{V}{\displaystyle \sum _{\mathrm{x}}{\displaystyle \sum _{i,j=1}^{3}2\left(\frac{{\partial }^{2}g(\mathrm{x};\mu )}{\partial x_i\partial x_j}\right)}}\frac{\partial }{\partial \mu }\frac{{\partial }^{2}g(\mathrm{x};\mu )}{\partial x_i\partial x_j}$$

(33)

where ${\partial }^{2}g(\mathrm{x};\mu )/\partial x_i\partial x_j$ denotes Hessian matrix of deformation function g(x; μ), $\partial \left({\partial }^{2}T/\partial x_i\partial x_j\right)/\partial \mu$ is the Jacobi of Hessian matrix.

Next, we calculate the derivative of KLD,

$$\frac{d\left[KLD\left(\nabla F\left(g\left(\mathrm{x};\mu \right)\right)\left|\right|\nabla R\left(\mathrm{x}\right)\right)\right]}{d\mu }={\displaystyle \sum _d{\displaystyle \sum _{\mathrm{x}}\left(1+\log \frac{{\tilde{q}}_d(f_j)}{{\tilde{p}}_d(r_i)}\right)\frac{\partial {\tilde{q}}_d(f_j)}{\partial \mu }}}$$

(34)

where ${\tilde{p}}_d(r_i)$ and ${\tilde{q}}_d(f_j)$ represent the gradient distributions of the transforms float image and reference image, respectively. The derivative of ${\tilde{q}}_d(f_j)$ is calculated as

$$\begin{array}{ll}\frac{\partial {\tilde{q}}_d(f_j)}{\partial \mu }=-\frac{1}{V\cdot \Delta b_F}& {\displaystyle \sum _{x\in \mathsf{\Omega }}\cdot {{\beta }^{\prime }}^{(3)}\left(f_j-\frac{\nabla F_d(g(\mathrm{x};\mu ))-\nabla F_d^0}{\Delta b_F}\right)}\\ & \cdot \left({\nabla }^2{F}_d(s){|}_{s=g(x;\mu )}\right)\cdot \frac{\partial (g(x;\mu ))}{\partial \mu }\end{array}$$

(35)

where β′⁽³⁾ is derivative of three-order B-spline, and ${\nabla }^2{F}_d(t)$ represents the second-order gradient of F(g(x; μ)), as well as $\partial (g(\mathrm{x};\mu ))/\partial \mu$ being the derivative of deformation function g(x; μ) with respect to the parameter μ. Substituting (12), (13), and (35) into (34), we can obtain the derivative of gradient distribution distance. In the terms of (29), (33), and (34), the derivative of E will be easily calculated.

4. Experiments and Results

To evaluate the registration method using the normalized JAD with gradient distribution (NJAD-GD), we designed several groups of tests and performed on simulated and real 3D data, respectively. In Section 4.1, the experimental data is depicted, including simulated and real medical images. The non-rigid registration of simulated MR volumes is performed in Section 4.2. The tests on real 3D thoracic CT images and 3D cardiac data are implemented, and the experimental results are shown in Section 4.3 and Section 4.4, respectively. Our nonrigid registration algorithm employing JAD and gradient distributions was implemented in the elastix package [26].

4.1. Experimental Data

In this paper, simulated brain MR volumes, thoracic CT volumes and real 3D cardiac CT images were exploited as experimental data. The detailed descriptions of brain MR and 3D thoracic CT images have been reported in [15]. Additionally, non-rigid tests were also performed on twelve 4D cardiac CT sequences acquired from twelve patients. Each of 4D CT sequence consists of 10 3D cardiac CT images, which were obtained from one whole cardiac cycle of one patient. These CT images have 256 × 256 pixels along axial direction. Figure 3 exhibits 10 3D cardiac CT volumes of one 4D CT sequence. It is obviously observed that some elastic deformations are existed between 10 images.

4.2. Nonrigid Registration of Simulated Brain Images

Simulated brain volumes were firstly used to design the elastic alignment experiments. Furthermore, we employ a multiresolution hierarchical strategy with three levels to carry out these non-rigid tests. Also, a comparison with JAD without gradient distribution and MI is also reported.

We selected 60 warping indexes (see parameters m of the warping function in [15]), which were yielded randomly from the interval [1,7]. Consequently, 60 float images were produced based on the 60 deformations for each pair of test volumes and 540 nonrigid trials of three pairs of brain MR volumes for NJAD-GD, JAD, and MI algorithms in total.

To assess quantitatively test results of these trails, we exploit registration error as the evaluation standard. Here, the registration error is defined as the difference of true values that can be calculated by warping indexes and the obtained values by optimization strategy. In the registration trails of brain images, the involved parameters are set as follows: the nonextensive parameter α = 1.5, bins M = 16, the number of random samples N = 2000, δ = 20 × 20 × 20. The weighting parameters λ₁ = 0.005 and λ₂ = 0.001 can provide a good tradeoff among three terms: D, S, and gradient distribution term KLD in the objective function E.

Figure 4 shows the test results of all 540 non-rigid registrations. From Figure 4, the NJAD-GD registration algorithm could result in the lower errors of three pairs of test volumes compared to other two approaches.

4.3. Experiments of 3D Thoracic CT Images

3D thoracic CT volumes were chosen as the test images to carry out non-rigid registrations. These volumes consist of four 4D sequences, and each of them includes 10 3D volumes. A three-level implementation scheme was still employed to decrease registration accuracy and improve computation efficiency.

We denote the 10 volumes from each 4D sequence by T00-T90, in which the maximal inhalation and maximum exhalation are included, with indicated by T10 and T60, respectively. Then, we designed the following experiments: in each 4D sequence, the T60 frame is applied as reference image and the residual nine frames are chosen as float image, leading to nine non-rigid tests. Hence, 36 trails of elastic alignments were yielded in total for four 4D CT sequences. We also compared the results adopting NJAD-GD and JAD without considering spatial information. Finally, 72 elastic registration tests were conducted for two methods. In order to quantify the test errors, target registration error (TRE) and Hausdorff distance meansure (HDM) were calculated. HDM is a widely-used measure to calculate the distance of two clouds of points. In the 3D thoracic CT registration, the manually marked landmarks can be applied to calculate HDM of two images.

Figure 5 and Figure 6 demonstrate the registration results of 72 tests, along with TREs before registration and after alignment. Figure 7 illustrates the box-and-whisker plots of HDM values of four 4D CT sequences. It is observed from these results that the registration errors applying NJAD-GD algorithm are less than these obtained by the method based on JAD without gradient distribution.

In implementation of experiments, the nonextensive parameter α = 1.5, bins M = 16, the number of random samples N = 8000, the spacing of mesh points δ = 20 × 20 × 20. The weighting parameters λ₁ = 0.005 and λ₂ = 0.001.

4.4. Registration of 3D Cardiac CT Image

The cardiac CT data consists of 12 groups of 4D image sequence, and each of which includes 10 3D images acquired from one whole cardiac cycle. One cardiac cycle consists of the phase of systole and phase of diastole. In each 4D CT sequence, two 3D images with the maximum deformation were employed as the test images, and 12 nonrigid registration experiments were carried out adopting NJAD-GD approach. Figure 8 illustrates the checkboard of 12 examples adopting our non-rigid framework (α = 1.50) for 12 4D CT sequences. As it can be seen, the registration algorithm based on gradient distribution demonstrates the accuracy results.

In these tests, a multiresolution scheme with three levels was also exploited to implement these nonrigid registrations. Due to the large deformation, the number of random samples N and the spacing of mesh points δ were set to 10,000 and 10 × 10 × 10, with other parameters being as follows: bins M = 16, the weighting parameters λ₁ = 0.005 and λ₂ = 0.001, the maximum number of iterations of the limited memory BFGS scheme NMAX = 200.

5. Conclusions

In this work, we review the definition and properties of Arimoto entropy, with an information measure based on Arimoto entropy, called JAD. The gradient distributions of reference image and float image are constructed and a distance between them is derived. Additionally, a normalized dissimilarity measure based on JAD was presented. A nonrigid registration method exploiting the normalized measure with gradient distributions is proposed.

Arimoto entropy is regarded as a generalized form of the classical Shannon entropy. In the aforementioned section, it is proofed that the JAD measure is equal to MI when α tends to 1. We adopted FFDs as the parameter space for non-rigid registration, along with objective function E including three elements: the normalized JAD as the dissimilarity measure, a regularization to acquire the smooth deformation and a distance term of the gradient distributions.