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

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 KullbackLeibler 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.


Introduction
Various kinds of iterative reconstruction (IR) [1][2][3][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][8][9][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' -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][16][17][18][19][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 continuoustime system using Lyapunov's stability theorem [21].The theory is extended to a case where the Lyapunov function can be an -skew -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 2 Mathematical Problems in Engineering 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.

ML-EM and MART Algorithms
The problem of image reconstruction from measured projections in CT is formulated as solving unknown variable  ∈   + for pixel values such that where  ∈   ++ ,  ∈  × + , and  ∈   represent the measured projection, projection operator, and noise, respectively, with  + and  ++ , respectively, denoting the sets of nonnegative and positive real numbers.If the system  =  has a nonnegative solution, it is consistent; otherwise, it is inconsistent.
We introduce the following definitions and notations.The variable denoted by the symbol   for  = 1, 2, . . .,  is defined by where   is the element in the th row and th column of .The notations   and   indicate the th element and the th row vector of the vector  and the matrix , respectively.The function KL(, ) of two nonnegative vectors  and  represents the generalized KL-divergence [5]: Known methods of reconstructing tomographic images in clinical CT are filtered back projection (FBP) as a transform method and IR as an optimization process.KLdivergence 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 sequences {KL(, ())} ∞ =0 and {KL((), )} ∞ =0 decrease for the respective iterative sequences {()} ∞ =0 of the ML-EM and the (simultaneous) MART algorithms defined by and for  = 1, 2, . . .,  with (0) =  0 ∈   ++ .

Proposed System
The IR algorithm proposed in this paper can be described using the following difference equation.

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 value  0 ∈   ++ for the system in (9), the solution (,  0 ) remains in the subspace   ++ for all  ≥ 0.
Proof.The vector field of the th element of the system is multiplied by   ; therefore, the derivative   / ≡ 0 holds for any  on the subspace   ≡ 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  ++ behaves in the substate space.

Let us define
() fl (1 − ) KL (, ) + KL (, ) with  in the interval [0, 1] as an -skew -divergence.The function is essential for a premetric measure of difference between  and  in accordance with   () ≥ 0 and   () = 0 if and only if  = .When  = 0.5, the divergence   becomes a symmetrical form of KL-divergence by averaging KL(, ) and KL(, ), 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 function   monotonically decreases along solution trajectories.

Theorem 2. If a unique solution 𝑒 ∈ 𝑅 𝐽
++ exists to the system  = , the equilibrium  of the dynamical system in (9) with any  ∈ [0, 1] is asymptotically stable.
Proof.A Lyapunov candidate function,   , is defined for applying Lyapunov's stability theorem in (10): which is well defined via Proposition 1 when choosing an initial value  0 in   ++ .We then calculate its derivative with respect to the dynamical system in (9) with any  ∈ [0, 1] as The derivative equals zero at  =  ∈   ++ .Consequently, the system has a Lyapunov function, so the equilibrium  of the system is asymptotically stable.with sufficiently smooth functions () and () satisfying () = () = 0, where  > 0 denotes an equilibrium of (13).

Hybrid Euler
Then, we have the geometric multiplicative derivative  ()  = exp ( ( ())) exp ( ( ())) , (0) =  0 , ( as a counterpart of the additive derivative.Analogous to the ordinary Euler method, the multiplicative first-order expansion [26,27]  ((1 + ) ) =  () (  ()  ) (15) with small  > 0 leads to the iterative formula of the multiplicative Euler discretization with iteration numbers  = 0, 1, 2, . ... However, when  ≡ 0 in (13), its additive Euler discretization corresponds to the multiplicative one in (16) with  ≡ 0. The term (1 + ()) is also considered as a first-order term in the Taylor series expansion of the function exp(()) with  in the neighborhood of the steady state .By replacing a part of the multiplicative term exp(()) in ( 16) with the term (1+()) while preserving the multiplication of , we obtain the formula of a hybrid Euler discretization  ( + 1) =  () (1 +  ( ())) exp ( ( ())) , for  = 0, 1, 2, . .., with a combination based on the additive and multiplicative calculus for the functions  and , 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 () and ().
The hybrid Euler discretization of the differential equation in (9) gives the IR algorithm in (6).

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 of  ∈   + shown in Figure 1 Let us first consider which is the most robust discretization method for numerical approximation.Figure 2 shows the time course of the objective function   (()) with  = 0.5, which we denote as  0.5 (()), for the continuous-time system with  = 0.5 at  and the discrete-time systems at  = , 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.
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  2 -norm of difference between the phantom image  and iterative state () obtained using the formula in (6) at the th iteration starting from the common initial value   (0) = ∑  =1   / ∑  =1 ∑    =1    ,  = 1, 2, . . ., , 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  is relatively small.In particular, the measure of SKLD is less than those of ML-EM and MART at, e.g.,  = 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][29][30][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.
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,  = 7, 680 (128 acquisition bins and 60 projection directions in 180 degrees) and  = 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. 10 4 , which is less than 4.61 × 10 4 obtained using FBP with a Shepp-Logan filter.
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 number  and the parameter  ∈ (0, 1) (or especially  = 0.5) in which there is a minimum value of the objective function  0.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).

Concluding Remarks
We proposed the IR algorithm based on the minimization of -skew -divergence between the measured and forward projections.The objective functions to be minimized with iterative process are KL(, ), KL(, ), and -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 -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.

Figure 5 Figure 2 :
Figure 2: Time course of  0.5 (()) as a function of t for continuous-time (magenta) and discrete-time (blue) systems using additive, multiplicative, and hybrid Euler methods with  = 0.01, 0.1, and 1.

Figure 3 :Figure 4 :
Figure 3: Values of RMS distance  using proposed algorithm in (6) at Nth iteration while varying  with (a) noisy and (b) noise-free projections.

Figure 5 :
Figure 5: (a) Sinogram and (b) reconstructed image using SKLD algorithm at 30th iteration in X-ray CT example.

Figure 6 :
Figure 6: Values of  0.5 using proposed algorithm in (6) at Nth iteration while varying  in (a) SPECT and (b) X-ray CT examples.