# Tomographic Image Reconstruction Based on Minimization of Symmetrized Kullback-Leibler Divergence

Authors: Ryosuke Kasai; Yusaku Yamaguchi; Takeshi Kojima; Tetsuya Yoshinaga
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

Iterative reconstruction (IR) algorithms based on the principle of optimization are known for producing better reconstructed images in computed tomography. In this paper, we present an IR algorithm based on minimizing a symmetrized Kullback-Leibler divergence (SKLD) that is called Jeffreys’J-divergence. The SKLD with iterative steps is guaranteed to decrease in convergence monotonically using a continuous dynamical method for consistent inverse problems. Specifically, we construct an autonomous differential equation for which the proposed iterative formula gives a first-order numerical discretization and demonstrate the stability of a desired solution using Lyapunov’s theorem. We describe a hybrid Euler method combined with additive and multiplicative calculus for constructing an effective and robust discretization method, thereby enabling us to obtain an approximate solution to the differential equation. We performed experiments and found that the IR algorithm derived from the hybrid discretization achieved high performance.


---

## Body

## 1. Introduction

Various kinds of iterative reconstruction (IR) [1–4] algorithms based on the principle of optimization are known for producing better reconstructed images in computed tomography (CT). In accordance with the base objective function, each IR method has intrinsic properties in the quality of images and in the convergence of iterative solutions. The maximum-likelihood expectation-maximization (ML-EM) [4] algorithm decreases the generalized Kullback-Leibler (KL) divergence [5, 6] between the measured projection and the forward projection along with the iteration. However, the objective function, in which the multiplicative algebraic reconstruction technique (MART) [7–10] forces a decrease in the iterative process, is an alternative KL-divergence done by exchanging the measured and forward projections. Due to the asymmetry of the KL-divergence, both objective functions have generally different values.In this paper, we present an IR algorithm based on the minimization of a symmetrized KL-divergence (SKLD), which is formulated using the mean of the two mutually alternative KL-divergences and is called Jeffreys’J-divergence [11, 12]. Because the base optimization function is a symmetric premetric measure and it gives an upper bound of the Jensen-Shannon divergence [13, 14], one can expect a better performance while preserving good properties of ML-EM and MART algorithms. The convergence to an exact solution and the monotonic decreasing of the SKLD with each iterative step for a consistent inverse problem are guaranteed using the approach of the continuous dynamical method [15–20]. Specifically, we construct an autonomous differential equation for which the proposed iterative formula gives a first-order numerical discretization with some step size and demonstrate the stability of an equilibrium in a continuous-time system using Lyapunov’s stability theorem [21]. The theory is extended to a case where the Lyapunov function can be an α-skew J-divergence with a parameter α∈[0,1], which is the SKLD when α=0.5.We also propose a novel method for constructing an effective and robust method in the discretization, thereby enabling us to obtain an approximate solution to the differential equation. That is, a hybrid Euler discretization combined with additive and multiplicative calculus works well; however, neither additive nor multiplicative Euler method does.We conducted experiments and found that the IR algorithm derived from the hybrid Euler discretization of the continuous-time dynamical system achieved high performance.


## 2. ML-EM and MART Algorithms

The problem of image reconstruction from measured projections in CT is formulated as solving unknown variablex∈R+J for pixel values such that (1)y=Ax+σ,where y∈R++I, A∈R+I×J, and σ∈RI represent the measured projection, projection operator, and noise, respectively, with R+ and R++, respectively, denoting the sets of nonnegative and positive real numbers. If the system y=Ax has a nonnegative solution, it is consistent; otherwise, it is inconsistent.We introduce the following definitions and notations. The variable denoted by the symbolλj for j=1,2,…,J is defined by (2)λj≔∑i=1IAij-1,where Aij is the element in the ith row and jth column of A. The notations zj and Ai indicate the jth element and the ith row vector of the vector z and the matrix A, respectively. The function KL(p,q) of two nonnegative vectors p and q represents the generalized KL-divergence [5]: (3)KLp,q=∑i=1Ipilog⁡piqi+qi-pi.Known methods of reconstructing tomographic images in clinical CT are filtered back projection (FBP) as a transform method and IR as an optimization process. KL-divergence is a common measure for deriving IR algorithms, and it plays an important role in accordance with the axiom that minimizing the KL-divergence is equivalent to maximizing the likelihood function that is considered suitable for reconstruction modeled with a probability distribution. The sequencesKL(y,Az(n))n=0∞ and KL(Az(n),y)n=0∞ decrease for the respective iterative sequences {z(n)}n=0∞ of the ML-EM and the (simultaneous) MART algorithms defined by (4)zjn+1=zjnλj∑i=1IAijyiAiznand (5)zjn+1=zjnexp⁡λj∑i=1IAijlog⁡yiAiznfor j=1,2,…,J with z(0)=x0∈R++J.


## 3. Proposed System

The IR algorithm proposed in this paper can be described using the following difference equation.(6)zjn+1=zjn1+δfjznexp⁡δgjzn,where δ>0 is a parameter and where the functions fj and gj are, respectively, defined by (7)fjzn=1-αλj∑i=1IAijyiAizn-1and (8)gjzn=αλj∑i=1IAijlog⁡yiAiznfor j=1,2,…,J, with α∈[0,1] and z(0)=x0∈R++J.For investigating the dynamics of the difference system in (6), we give its continuous analog described using the differential equation: (9)dxjtdt=λjxjt∑i=1IAij1-αyiAixt-1+αlog⁡yiAixtfor j=1,2,…,J and t≥0 with x(0)=x0∈R++J. The term “continuous analog” means that a numerical discretization of the differential equation is the same as the difference equation. The idea is based on the approach for solving a constrained tomographic inverse problem using nonlinear dynamical methods [22–25].


### 3.1. Theoretical Analysis

The theoretical results on the continuous-time system in (9) are given. First, we show that a solution is possible to ensure positive values remain.Proposition 1.
If we choose an initial valuex0∈R++J for the system in (9), the solution ϕ(t,x0) remains in the subspace R++J for all t≥0.Proof.
The vector field of thejth element of the system is multiplied by xj; therefore, the derivative dϕj/dt≡0 holds for any j on the subspace xj≡0, which means the subspace is invariant and which means any flow cannot pass through invariant subspace on the basis of the uniqueness property of solutions to the initial value problem. Hence, the trajectory emanating from the initial value in R++ behaves in the sub-state space.Let us define(10)Jαx≔1-αKLy,Ax+αKLAx,y=∑i=1I1-αyilog⁡yiAix+Aix-yi+αAixlog⁡Aixyi+yi-Aixwith α in the interval [0,1] as an α-skew J-divergence. The function is essential for a premetric measure of difference between y and Ax in accordance with Jα(x)≥0 and Jα(x)=0 if and only if y=Ax. When α=0.5, the divergence Jα becomes a symmetrical form of KL-divergence by averaging KL(y,Ax) and KL(Ax,y), which we abbreviate as SKLD.Then, under the assumption that the inverse problem is consistent, the following shows that the solution converges to a desired equilibrium and that the nonnegative functionJα monotonically decreases along solution trajectories.Theorem 2.
If a unique solutione∈R++J exists to the system y=Ax, the equilibrium e of the dynamical system in (9) with any α∈[0,1] is asymptotically stable.Proof.
A Lyapunov candidate function,Jα, is defined for applying Lyapunov’s stability theorem in (10): (11)Jαx=∑i=1I1-α∫yiAix1-yiwdw+α∫yiAixlog⁡wyi=∑i=1I1-αw-yilog⁡wyiAix+αwlog⁡wyi-wyiAix,which is well defined via Proposition 1 when choosing an initial value x0 in R++J. We then calculate its derivative with respect to the dynamical system in (9) with any α∈[0,1] as (12)dJαdtx9=∑i=1I1-α1-yiAix+αlog⁡Aixyi∑j=1JAijdxjdt=-∑j=1J∑i=1IAij1-αyiAix-1+αlog⁡yiAixλjxj∑k=1IAkj1-αykAkx-1+αlog⁡ykAkx=-∑j=1Jλjxj∑i=1IAij1-αyiAix-1+αlog⁡yiAix2≤0.The derivative equals zero at x=e∈R++J. Consequently, the system has a Lyapunov function, so the equilibrium e of the system is asymptotically stable.


### 3.2. Hybrid Euler Discretization

Consider an autonomous differential equation for the state variablex(t)>0, t≥0, written by (13)dxtdt=xtfxt+gxt,x0=x0,with sufficiently smooth functions f(x) and g(x) satisfying f(e)=g(e)=0, where e>0 denotes an equilibrium of (13). Then, we have the geometric multiplicative derivative (14)πxtπt=exp⁡fxtexp⁡gxt,x0=x0,as a counterpart of the additive derivative. Analogous to the ordinary Euler method, the multiplicative first-order expansion [26, 27](15)x1+δt=xtπxtπtδwith small δ>0 leads to the iterative formula of the multiplicative Euler discretization (16)zn+1=znexp⁡δfznexp⁡δgzn,z0=x0,with iteration numbers n=0,1,2,….However, wheng≡0 in (13), its additive Euler discretization (17)zn+1=zn1+δfzncorresponds to the multiplicative one in (16) with g≡0. The term (1+δf(z)) is also considered as a first-order term in the Taylor series expansion of the function exp(δf(z)) with z in the neighborhood of the steady state e. By replacing a part of the multiplicative term exp(δf(z)) in (16) with the term (1+δf(z)) while preserving the multiplication of z, we obtain the formula of a hybrid Euler discretization (18)zn+1=zn1+δfznexp⁡δgzn,z0=x0,for n=0,1,2,…, with a combination based on the additive and multiplicative calculus for the functions f and g, respectively. The hybrid Euler method is effective for a practical calculation of an initial value problem in (13) from the viewpoint of choosing a larger value of the step-size δ, when either an additive or multiplicative calculus is better for each term of the partial vector fields xf(x) and xg(x).The hybrid Euler discretization of the differential equation in (9) gives the IR algorithm in (6).


## 3.1. Theoretical Analysis

The theoretical results on the continuous-time system in (9) are given. First, we show that a solution is possible to ensure positive values remain.Proposition 1.
If we choose an initial valuex0∈R++J for the system in (9), the solution ϕ(t,x0) remains in the subspace R++J for all t≥0.Proof.
The vector field of thejth element of the system is multiplied by xj; therefore, the derivative dϕj/dt≡0 holds for any j on the subspace xj≡0, which means the subspace is invariant and which means any flow cannot pass through invariant subspace on the basis of the uniqueness property of solutions to the initial value problem. Hence, the trajectory emanating from the initial value in R++ behaves in the sub-state space.Let us define(10)Jαx≔1-αKLy,Ax+αKLAx,y=∑i=1I1-αyilog⁡yiAix+Aix-yi+αAixlog⁡Aixyi+yi-Aixwith α in the interval [0,1] as an α-skew J-divergence. The function is essential for a premetric measure of difference between y and Ax in accordance with Jα(x)≥0 and Jα(x)=0 if and only if y=Ax. When α=0.5, the divergence Jα becomes a symmetrical form of KL-divergence by averaging KL(y,Ax) and KL(Ax,y), which we abbreviate as SKLD.Then, under the assumption that the inverse problem is consistent, the following shows that the solution converges to a desired equilibrium and that the nonnegative functionJα monotonically decreases along solution trajectories.Theorem 2.
If a unique solutione∈R++J exists to the system y=Ax, the equilibrium e of the dynamical system in (9) with any α∈[0,1] is asymptotically stable.Proof.
A Lyapunov candidate function,Jα, is defined for applying Lyapunov’s stability theorem in (10): (11)Jαx=∑i=1I1-α∫yiAix1-yiwdw+α∫yiAixlog⁡wyi=∑i=1I1-αw-yilog⁡wyiAix+αwlog⁡wyi-wyiAix,which is well defined via Proposition 1 when choosing an initial value x0 in R++J. We then calculate its derivative with respect to the dynamical system in (9) with any α∈[0,1] as (12)dJαdtx9=∑i=1I1-α1-yiAix+αlog⁡Aixyi∑j=1JAijdxjdt=-∑j=1J∑i=1IAij1-αyiAix-1+αlog⁡yiAixλjxj∑k=1IAkj1-αykAkx-1+αlog⁡ykAkx=-∑j=1Jλjxj∑i=1IAij1-αyiAix-1+αlog⁡yiAix2≤0.The derivative equals zero at x=e∈R++J. Consequently, the system has a Lyapunov function, so the equilibrium e of the system is asymptotically stable.


## 3.2. Hybrid Euler Discretization

Consider an autonomous differential equation for the state variablex(t)>0, t≥0, written by (13)dxtdt=xtfxt+gxt,x0=x0,with sufficiently smooth functions f(x) and g(x) satisfying f(e)=g(e)=0, where e>0 denotes an equilibrium of (13). Then, we have the geometric multiplicative derivative (14)πxtπt=exp⁡fxtexp⁡gxt,x0=x0,as a counterpart of the additive derivative. Analogous to the ordinary Euler method, the multiplicative first-order expansion [26, 27](15)x1+δt=xtπxtπtδwith small δ>0 leads to the iterative formula of the multiplicative Euler discretization (16)zn+1=znexp⁡δfznexp⁡δgzn,z0=x0,with iteration numbers n=0,1,2,….However, wheng≡0 in (13), its additive Euler discretization (17)zn+1=zn1+δfzncorresponds to the multiplicative one in (16) with g≡0. The term (1+δf(z)) is also considered as a first-order term in the Taylor series expansion of the function exp(δf(z)) with z in the neighborhood of the steady state e. By replacing a part of the multiplicative term exp(δf(z)) in (16) with the term (1+δf(z)) while preserving the multiplication of z, we obtain the formula of a hybrid Euler discretization (18)zn+1=zn1+δfznexp⁡δgzn,z0=x0,for n=0,1,2,…, with a combination based on the additive and multiplicative calculus for the functions f and g, respectively. The hybrid Euler method is effective for a practical calculation of an initial value problem in (13) from the viewpoint of choosing a larger value of the step-size δ, when either an additive or multiplicative calculus is better for each term of the partial vector fields xf(x) and xg(x).The hybrid Euler discretization of the differential equation in (9) gives the IR algorithm in (6).


## 4. Experimental Results and Discussion

We conducted experiments to demonstrate the effectiveness of not only the proposed IR method in (6) with α=0.5 and δ=1, say SKLD algorithm, based on the minimization of SKLD but also the hybrid Euler discretization of the continuous analog in (9).We used a modified Shepp-Logan phantom image consisting ofe∈R+J shown in Figure 1(a), which is composed of 128×128 pixels (J=16,384). The projection y∈R+I was simulated by using the model y=Ae+σ with σ denoting the white Gaussian noise such that the signal-to-noise ratio (SNR) was 30 dB unless otherwise specified and by setting the number of view angles and detector bins to 180 and 184, respectively (I=33,120), as illustrated in Figure 1(b) formatted into a sinogram. A standard MATLAB (MathWorks, Natick, USA) ODE solver ode113 implementing a variable step-size variable order Adams-Bashforth-Moulton method was used for a numerical calculation of solving an initial value problem of the differential equation in (9).Figure 1
(a) Phantom image and (b) sinogram with noisy projection.
(a) (b)Let us first consider which is the most robust discretization method for numerical approximation. Figure2 shows the time course of the objective function Jα(x(t)) with α=0.5, which we denote as J0.5(x(t)), for the continuous-time system with α=0.5 at t and the discrete-time systems at t=δn, which are discretizations of the differential equation using the additive, multiplicative, and hybrid Euler methods with the fixed step sizes δ of 0.01, 0.1, and 1. We see that the hybrid Euler discretization gives a similar good convergence to that of a continuous-time system; however, in the case δ=1, the additive Euler method causes, as time goes on, an error with oscillation and solutions using the multiplicative method divergence.Figure 2
Time course ofJ0.5(x(t)) as a function oftfor continuous-time (magenta) and discrete-time (blue) systems using additive, multiplicative, and hybrid Euler methods with δ = 0.01, 0.1, and 1.Next, we consider the property of the SKLD algorithm compared with the algorithms using the ML-EM method and MART. Note that the iterative formula in (6) with δ=1 becomes the ML-EM, SKLD, and MART algorithms when the parameter α is equal to 0, 0.5, and 1, respectively. The values of the root mean squared (RMS) distance or L2-norm of difference between the phantom image e and iterative state z(N) obtained using the formula in (6) at the Nth iteration starting from the common initial value zj(0)=∑i=1Iyi/∑i=1I∑j′=1JAij′, j=1,2,…,J, and variations in the parameter α are plotted in Figure 3(a). We observed a value of α exists in the open interval (0,1) minimizing the quantitative measure when the number N is relatively small. In particular, the measure of SKLD is less than those of ML-EM and MART at, e.g., N=10 and 12. The presence of noise is required for this property, which is supported by the fact that an experiment from another simulation with noise-free projection data gave different results, as shown in Figure 3(b), such that we had the minimum value of the measure at α=1. This was due to characteristics of the ML-EM and MART algorithms being disadvantageous. As a result of these characteristics, ML-EM has a slow convergence [28–31], which means that a relatively large number of iterations is required to obtain an acceptable value of the objective function, and MART is more susceptible to noise [32, 33], which may result in an increase of the convergence rate with increasing iterative steps due to the noisy projection, as typically seen at α=0 and 1 in Figure 3, respectively.Figure 3
Values of RMS distanced using proposed algorithm in (6) atNth iteration while varying α with (a) noisy and (b) noise-free projections.
(a) (b)We carried out other experiments to examine whether there exists anα not equal to 0 or 1, in which the value of an objective function obtained by using the algorithm in (6) with δ=1 takes a minimum. A set of SPECT projection data from a publicly accessible database [34] was examined first. The sinogram of a brain scan and an image reconstructed by the SKLD algorithm are shown in Figure 4, where the number of projections and pixels of the image are, respectively, I=7,680 (128 acquisition bins and 60 projection directions in 180 degrees) and J=7,569 (87×87 pixels). The projection data for the second example were acquired from an X-ray CT scanner (Toshiba Medical Systems, Tochigi, Japan) with a body phantom [35] (Kyoto Kagaku, Kyoto, Japan) using a tube voltage of 120 kVp and a tube current of 300 mA. Figure 5 represents the sinogram with I=86,130 (957 acquisition bins and 90 projection directions in 180 degrees) and a reconstructed image with J=454,276 (674×674 pixels). The value of the objective function J0.5 for the image reconstructed by SKLD algorithm at 30th iteration is 2.86×104, which is less than 4.61×104 obtained using FBP with a Shepp-Logan filter.Figure 4
(a) Sinogram and (b) reconstructed image using SKLD algorithm at 30th iteration in SPECT example.
(a) (b)Figure 5
(a) Sinogram and (b) reconstructed image using SKLD algorithm at 30th iteration in X-ray CT example.
(a) (b)In the physical experiments, as well as a feature of the object to be reconstructed and also an overdetermined or underdetermined inverse problem, there exists a nonempty set of the iteration numberN and the parameter α∈(0,1) (or especially α=0.5) in which there is a minimum value of the objective function J0.5 as shown in Figure 6. From Figure 6(a), we observe that the algorithm with α≥0.1 sufficiently converges to a local minimum at a small iteration number, whereas the ML-EM algorithm has a slow convergence rate when the state variable is far away from the local minimum in the state space. Considering the noisy nature of the measured projection and large datasets required for reconstructing large sized images in clinical X-ray CT, the conditions suitable for practical use include a relatively small number of iterations exclusively used in IR systems [36] under which our IR method is effective, as seen in Figure 6(b).Figure 6
Values ofJ0.5 using proposed algorithm in (6) atNth iteration while varying α in (a) SPECT and (b) X-ray CT examples.
(a) (b)


## 5. Concluding Remarks

We proposed the IR algorithm based on the minimization ofα-skew J-divergence between the measured and forward projections. The objective functions to be minimized with iterative process are KL(y,Az), KL(Az,y), and J-divergence after setting the value of parameter α to 0, 1, and 0.5 in the IR algorithm, respectively. We used Lyapunov’s theorem to construct a differential equation for which the IR algorithm is equivalent to the hybrid Euler discretization and demonstrated the monotonically decreasing convergence of the J-divergence with iterative steps to a desired solution of the continuous-time system. The hybrid Euler was found to have more robust performance than the additive and multiplicative Euler methods. The numerical and physical experiments revealed that the SKLD method is advantageous with respect to the minimization of a distance measure when the projection data are noisy and when the maximum iteration number is relatively small. Further investigation is required to clarify the performance of the proposed IR algorithm from the viewpoint of image quality and noise evaluation.



---

*Source: 8973131-2018-07-17.xml*