**1. Introduction**

Image reconstruction in computed tomography (CT) is the process of estimating unknown density images from measured projections. When the system of a tomographic inverse problem is ill-posed, iterative reconstruction algorithms [1,2] based on the optimization strategy provide significant improvements over transform methods, including the filtered back-projection [3,4] (FBP) procedure. In recent years, iterative reconstruction has received much attention because of its ability to reduce radiation doses [5–9] in X-ray CT. Iterative algorithms implemented in, e.g., the algebraic reconstruction technique [1], maximum-likelihood expectation-maximization [10] (MLEM) method, and multiplicative algebraic reconstruction technique, have been used for reconstructing CT images. The MLEM algorithm, which is the most popular method used in emission CT and is derived for the maximum likelihood estimation of a Poisson distribution, reconstructs high-quality images even for noisy projection data, but it is slow to converge [11–14] under iteration. The ordered-subsets EM algorithm [11], in which the EM iteration is performed in each subset by dividing the projection into subsets or blocks, is known to be useful for accelerating MLEM [13,15,16]. However, divergence or oscillation of solutions may occur in the iterative process when the subset balance is not satisfied. Because of the high quality of

**Citation:** Kasai, R.; Yamaguchi, Y.; Kojima, T.; Abou Al-Ola, O.M.; Yoshinaga, T. Noise-Robust Image Reconstruction Based on Minimizing Extended Class of Power-DivergenceMeasures. *Entropy* **2021**, *23*, 1005. https://doi.org/10.3390/e23081005

Academic Editor: Jiayi Ma

Received: 1 June 2021 Accepted: 24 July 2021 Published: 31 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

image reconstruction afforded by the MLEM algorithm, improved MLEM methods have been presented for accelerating convergence. Some schemes accelerate the convergence rate by increasing a relaxation parameter or the step-size in iterative operations [14,17,18] or by introducing a parameter with a power exponent related to the projection for controlling the noise model [19,20]. However, no theory has explained the divergence and oscillation phenomena affecting solutions when the step-size parameter is large.

The convergence of iterative solutions and the quality of images are governed by the underlying objective function that has to be minimized. Hence, the base objective function is one of the most important decisions when designing an iterative algorithm. In this paper, we present an extended class of power-divergence measures [21–24] (PDMs) and derive a novel iterative algorithm based on the minimization of the extended PDM (EPDM) as an objective function for the optimization strategy. Let us define the parameterized function <sup>Φ</sup>*γ*,*<sup>α</sup>*(*p*, *q*) of vectors *p* and *q* with nonnegative elements *pi* and *qi*, respectively, as

$$\Phi\_{\gamma,a}(p,q) := \sum\_{i} \int\_{p\_i}^{q\_i} \frac{\mathbf{s}^{\gamma} - p\_i^{\gamma}}{\mathbf{s}^{\gamma a}} ds \tag{1}$$

where *γ* and *α* indicate positive and nonnegative parameters, respectively. The extension is performed by incorporating the parameter *α* in the conventional class of PDMs, which includes a large set of distance and relative entropy measures. By fixing the parameter *α* = 1, <sup>Φ</sup>*γ*,<sup>1</sup> gives the family of PDMs with a single parameter *γ*. Therefore, the measure coincides with the Kullback–Leibler (KL), or relative entropy, divergence [25] if *γ* = 1, Neyman's *χ*2 distance if *γ* = 2, and the generalized Hellinger distance otherwise. Moreover, it corresponds to the squared *L*<sup>2</sup> norm when *γ* = 1 and *α* = 0 and the reverse KL-divergence when *γ* = 1 and *α* = 2. Thus, the parameters *γ* and *α* provide a smooth connection among the forward and reverse KL-divergences, the Hellinger distance, the *χ*2 distance, and the *L*<sup>2</sup> distance and can control the trade-off between robustness and asymptotic efficiency of the estimators, in a similar way as in other families of distance measures [26–29].

By exploiting the vectors *p* and *q* in Equation (1) as the measured and forward projections, respectively, for the tomographic inverse problem, it is expected that we can create a high-performance iterative reconstruction algorithm thanks to the high degree of freedom. For constructing this novel iterative algorithm, we introduce a nonlinear differential equation whose numerical discretization is equivalent to the iterative system. The purpose of applying a dynamical method [30–35] to tomographic inverse problems [36–39] is as follows: it enables us to prove the stability of the equilibrium corresponding to the desired solution of the system of differential equations by using the Lyapunov stability theorem [40] if a proper Lyapunov function can be found; then, since the step-size used to discretize the set of differential equations corresponds to the relaxation or scaling parameter of the system of difference equations, a family of iterative algorithms incorporating a parameter for acceleration is naturally derived. Moreover, it provides a methodology for systematically designing a new iterative reconstruction algorithm based on optimization of an objective function depending on the features of the image to be reconstructed.

Since the EPDM family includes the KL-divergence, the resulting iterative algorithm with power exponents corresponding to the parameters *γ* and *α* is a natural extension of the MLEM method with (*<sup>γ</sup>*, *α*)=(1, <sup>1</sup>). The convergence of solutions to the continuous analog of the proposed iterative algorithm is theoretically shown using the EPDM as a Lyapunov function when the tomographic inverse problem is consistent.

We conducted image reconstruction experiments using numerical and physical phantoms with noisy projections and found that the proposed algorithm outperformed the conventional MLEM method with respect to reconstructing high-quality images that are robust to measured noise when selecting a set of proper parameter values.

#### **2. Definitions and Notations**

Image reconstruction is a problem of obtaining unknown pixel values *x* ∈ *RJ*+ satisfying

$$y = A\mathfrak{x} + \delta \,\tag{2}$$

where *y* ∈ *<sup>R</sup>I*+, *A* ∈ *RI*×*<sup>J</sup>* + , and *δ* ∈ *R<sup>I</sup>* denote the measured projection, projection operator, and noise, respectively, with *R*+ representing the set of nonnegative real numbers. When the system in Equation (2) without noise, i.e., *δ* = 0, has a solution *e* ∈ *<sup>R</sup>J*+, it is consistent; otherwise, it is inconsistent. The tomographic inverse problem can be reduced to one of finding *x*, which can be accomplished using an optimization approach such as an iterative method or a continuous-time system by minimizing an objective function.

Here, we introduce the notation that will be used below. The superscript stands for the transpose of a matrix or vector, *θk* indicates the *k*th element of the vector *θ*, Θ*i* and <sup>Θ</sup>*ij* indicate the *i*th row vector and the element in the *i*th row and *j*th column of the matrix Θ, respectively, Log(*θ*) and Exp(*θ*) are, respectively, the vector-valued functions Log(*θ*) := (log(*<sup>θ</sup>*1), log(*<sup>θ</sup>*2), ... , log(*<sup>θ</sup>i*)) and Exp(*θ*) := (exp(*<sup>θ</sup>*1), exp(*<sup>θ</sup>*2), ... , exp(*<sup>θ</sup>i*)) of each element in vector *θ* = (*<sup>θ</sup>*1, *θ*2, ... , *<sup>θ</sup>i*), and diag(*θ*) indicates the diagonal matrix in which the diagonal entries are the elements of the vector *θ*.

#### **3. Proposed System**
