*3.1. Definition*

The proposed methods for obtaining a solution to the tomographic inverse problem can be described as an iterative algorithm and dynamical system.

We present an iterative reconstruction method with a relaxation or scaling parameter *h* > 0:

$$z\_j(n+1) = z\_j(n) \left( f\_j(z(n)) \right)^h \tag{3}$$

with

$$f\_{\vec{f}}(w) := \frac{\sum\_{i=1}^{I} A\_{ij} \left(\frac{y\_i}{\left(A\_i w\right)^{\vec{\alpha}}}\right)^{\gamma}}{\sum\_{i=1}^{I} A\_{ij} \left(\frac{A\_i w}{\left(A\_i w\right)^{\vec{\alpha}}}\right)^{\gamma}} \tag{4}$$

for *j* = 1, 2, ... , *J* and *n* = 0, 1, 2, ... , *N* − 1, where *γ* > 0, *α* ≥ 0, and *z*(0) = *z*0 ∈ *<sup>R</sup>J*++, with *R*++ representing the set of positive real numbers. The accompanying system derived from a continuous analog based on the dynamical method is described by a dynamical system:

$$\frac{d\mathbf{x}\_{\circ}(t)}{dt} = \mathbf{x}\_{\circ}(t)\log\left(f\_{\circ}(\mathbf{x}(t))\right) \tag{5}$$

for *j* = 1, 2, ... , *J* at *t* ≥ 0, where the function *fj* is in Equation (4) and *x*(0) = *z*0. The system in Equation (5) can be equivalently written as

$$\frac{d\mathbf{x}(t)}{dt} = \begin{array}{c} \mathrm{X}\Big(\mathrm{Log}(\boldsymbol{A}^{\top}\operatorname{Exp}(\gamma(\mathrm{Log}(\boldsymbol{y}) - \boldsymbol{a}\operatorname{Log}(\mathrm{Ax}(t))))\Big) \\\\ -\mathrm{Log}(\boldsymbol{A}^{\top}\operatorname{Exp}(\gamma(1-\boldsymbol{a})\operatorname{Log}(\boldsymbol{A}\mathbf{x}(t)))) \Big) \end{array} \tag{6}$$

where *X* := diag(*x*). The iterative formula in Equation (3) is obtained by discretizing the differential equation of Equation (5) by using the multiplicative Euler method [41,42] with a step-size of *h*. Note that the iterative formula in Equation (3) with *h* = 1 is equivalent to the algorithm presented by Zeng [19] when *γ* = 1, to the algorithm in Reference [20] when *α* = 1, and to the MLEM algorithm when (*<sup>γ</sup>*, *α*)=(1, <sup>1</sup>).

We apply the proposed divergence in Equation (1) to the tomographic objective function consisting of measured and forward projections. Namely, we define

$$V(\mathbf{x}) := \Phi\_{\gamma, \mathbf{a}}(y, \mathbf{A}\mathbf{x}) = \sum\_{i=1}^{I} \int\_{y\_i}^{A\_i \mathbf{x}} \frac{s^{\gamma} - y\_i^{\gamma}}{s^{\gamma a}} ds \tag{7}$$

which can be written as

$$\begin{aligned} V(\mathbf{x}) &= \sum\_{i=1}^{I} \int\_{y\_i}^{A\_i \mathbf{x}} \frac{s^\gamma - y\_i^\gamma}{s} ds \\ &= \frac{1}{\gamma} \sum\_{i=1}^{I} y\_i^\gamma \left( \log \left( \frac{y\_i}{A\_i \mathbf{x}} \right)^\gamma + \left( \frac{A\_i \mathbf{x}}{y\_i} \right)^\gamma - 1 \right) \end{aligned}$$

if *γα* = 1;

$$\begin{aligned} V(\mathbf{x}) &= \sum\_{i=1}^{I} \int\_{y\_i}^{A\_i \mathbf{x}} \frac{\mathbf{s}^\gamma - y\_i^\gamma}{\mathbf{s}^{1+\gamma}} ds \\ &= \frac{1}{\gamma} \sum\_{i=1}^{I} \log \left( \frac{A\_i \mathbf{x}}{y\_i} \right)^\gamma + \left( \frac{y\_i}{A\_i \mathbf{x}} \right)^\gamma - 1 \end{aligned}$$

if *γα* = 1 + *γ*; and

$$\begin{aligned} V(\mathbf{x}) &= \sum\_{i=1}^{I} \int\_{y\_i}^{A\_i \mathbf{x}} \frac{\mathbf{s}^\gamma - y\_i^\gamma}{\mathbf{s}^{\gamma \alpha}} d\mathbf{s} \\ &= \sum\_{i=1}^{I} \frac{1}{1 - \gamma \alpha} y\_i^{1 + \gamma(1 - \alpha)} \left( 1 - \left( \frac{A\_i \mathbf{x}}{y\_i} \right)^{1 - \gamma \alpha} \right) \\ &+ \frac{1}{1 + \gamma(1 - \alpha)} y\_i^{1 + \gamma(1 - \alpha)} \left( \left( \frac{A\_i \mathbf{x}}{y\_i} \right)^{1 + \gamma(1 - \alpha)} - 1 \right) \end{aligned}$$

otherwise.

#### *3.2. Theoretical Results*

This section provides a theoretical result on the dynamical system defined in the preceding section. We show that any solution to the continuous analog converges to the desired solution of the system in Equation (2) with *δ* = 0 when the inverse problem is consistent.

**Theorem 1.** *Assume there exists e* ∈ *<sup>R</sup>J*++ *satisfying y* = *Ae. Then, e is an equilibrium observed in the continuous-time system in Equation (6) and is asymptotically stable.*

**Proof.** We see that *e* is an equilibrium of the system and the solutions to the system are in *<sup>R</sup>J*++ because the initial state value at *t* = 0 belongs to *<sup>R</sup>J*++ and the flow cannot pass through the invariant subspace *xj* = 0 for *j* = 1, 2, ... , *J* in the state space according to the uniqueness of solutions for the initial value problem. The nonnegative function *<sup>V</sup>*(*x*) of *xj* > 0 in Equation (7) is well-defined as a candidate of a Lyapunov function. Then, we have the derivative of *V* along the solutions to Equation (6):

$$\begin{split} \frac{dV}{dt}(\mathbf{x})\Big|\_{\left(\delta\right)} &= \quad -\sum\_{i=1}^{I} \Big( \Big( \frac{y\_i}{(A\_i \mathbf{x})^a} \Big)^{\gamma} - (A\_i \mathbf{x})^{\gamma(1-a)} \Big) A\_i \frac{d\mathbf{x}}{dt} \\ &= \quad - (\boldsymbol{\xi} - \boldsymbol{\xi})^{\top} X (\mathrm{Log}(\boldsymbol{\xi}) - \mathrm{Log}(\boldsymbol{\xi})) \\ &\leq \quad 0 \end{split} \tag{8}$$

where

$$\begin{array}{rcl} \xi & := & A \, ^\top \operatorname{Exp}(\gamma(\operatorname{Log}(y) - \mathfrak{a} \operatorname{Log}(A\mathfrak{x}))) \\ \zeta & := & A \, ^\top \operatorname{Exp}(\gamma(1-\mathfrak{a}) \operatorname{Log}(A\mathfrak{x})) .\end{array}$$

Therefore, *V* is a Lyapunov function and the equilibrium *e* is asymptotically stable.

This theorem guarantees that the proposed difference system in Equation (3) as a first-order approximation to the differential equation in Equation (6) has a stable fixed point *e* when the chosen step-size *h* is sufficiently small to ensure numerical stability.

#### **4. Experimental Results and Discussion**

We will illustrate the effectiveness of the extended MLEM algorithm based on the EPDM family in Equation (3) with the parameter set (*<sup>γ</sup>*, *α*) (in what follows, the iterative algorithm except for MLEM with (*<sup>γ</sup>*, *α*)=(1, 1) is referred to as PDEM) by using examples from numerical and physical CT experiments. The proposed systems were executed using a 6-core processor and computing tools provided by MATLAB (MathWorks, Natick, MA, USA).

We set *h* = 1 and a constant initial value *z*0 *j* for *j* = 1, 2, ... , *J*. Note that variation of *h* is out of the scope of this paper, although setting *h* > 1 would accelerate convergence. In addition, in the numerical simulation, the PDEM algorithm in Equation (3) with *h* = 1 as a simple forward Euler discretization qualitatively approximates the solutions to the differential equation in Equation (6), which were calculated by a standard MATLAB ODE solver ode113 implementing a variable step-size variable order method.

#### *4.1. Reconstruction Using Numerical Phantom*

We used a numerical phantom image consisting of *e* ∈ [0, 1] *J* with 128 × 128 pixels (*J* = 16,384), as shown in Figure 1. For our experiment, a Shepp–Logan phantom [43], which is a popular test image for developing reconstruction algorithms, was modified by changing the density values for ellipses so that the resulting image had better visual perception with high contrast. The noise-free and noisy projections *y* ∈ *R<sup>I</sup>* + derived from the phantom image were, respectively, simulated using Equation (2) without and with *δ* denoting white Gaussian noise such that the signal-to-noise ratio (SNR) was 30 dB and by setting the number of view angles and detector bins to 180 and 184 (*I* = 33,120) with 180-degree sampling.

**Figure 1.** Image of numerical phantom.

For directly evaluating the quality of the reconstructed images, we defined functions for comparing the reconstructed image compared against the true image, *e*, as

$$D\_j(z) := |\mathfrak{e}\_j - z\_j|\tag{9}$$

for *j* = 1, 2, . . . , *J* and

$$E(z) := ||\mathfrak{e} - \mathfrak{z}||\_2 = \left(\sum\_{j=1}^{J} \left(D\_j(z)\right)^2\right)^{\frac{1}{2}}.\tag{10}$$

First, we considered the case of a noise-free projection. Figure 2 shows the evaluation functions *<sup>E</sup>*(*z*(*n*)) of the iterative points *z*(*n*) for MLEM and PDEM with the sets of parameters (*<sup>γ</sup>*, *α*) being (0.3, 1.2), (0.5, 1.2), (0.8, 1.2), and (1.3, 1.2) for *n* = 0, 1, 2, ... , 200. All algorithms monotonically decreased the evaluation function, as supported by the theoretical result that the solutions of the continuous analog converge to the true value. Indeed, another experiment confirmed that the monotonic decrease continued as the number of iterations exceeded 200 iterations. We can see that PDEM with the parameter set (1.3, 1.2) reduces the evaluation function much more than MLEM does. To put it another way, the PDEM algorithm takes less computation time than MLEM for obtaining the same evaluation values.

**Figure 2.** Evaluation functions for MLEM and PDEM algorithms at each iteration in the experiment using numerical phantom with noise-free projection. Note that because the values of PDEM with (*<sup>γ</sup>*, *α*)=(0.8, 1.2) and MLEM are very similar, the plotted points for PDEM are almost invisible.

Figure 3 shows contour plots of the evaluation values on a logarithmic scale, log10(*E*(*z*(*N*))) for *N* = 50, 100, and 200, in the parameter plane (*<sup>γ</sup>*, *<sup>α</sup>*). The parameters *γ* and *α* were, respectively, sampled from 0.1 to 1.5 and 0 to 1.4 with a sampling interval of 0.1. We can see that, at least in the range examined, the evaluation function becomes smaller as the values of *γ* ≥ 1 and *α* ≥ 1 increase.

**Figure 3.** Contour plots of evaluation functions log10(*E*(*z*(*N*))) with *N* being (**a**) 50, (**b**) 100, and (**c**) 200 using numerical phantom with noise-free projections. The white dot indicates the position of MLEM.

Figure 4 illustrates images reconstructed by MLEM and PDEM with (*<sup>γ</sup>*, *α*)=(1.3, 1.2) at the 200th iteration and the corresponding subtraction images *Dj*(*z*(200)) (displayed in the range from 0 to 0.2) at every *j*th pixel, for *j* = 1, 2, ... , *J*. By comparing the values of the subtraction between MLEM and PDEM, e.g., the edges of the high-density objects in the image, we can see that PDEM produces high-quality reconstructions, as is quantitatively indicated by its small evaluation value between the reconstructed and phantom images.

**Figure 4.** Reconstructed images (**upper panel**) and images of the subtraction (**lower panel**) for MLEM and PDEM with (*<sup>γ</sup>*, *α*)=(1.3, 1.2) at 200th iteration using numerical phantom with noisefree projection.

Next, let us consider the effect of the measured noise. Figure 5 is a graph of the evaluation *<sup>E</sup>*(*z*(*n*)) as a function of iteration number *n* with *n* = 0, 1, 2, ... , 200. Given noisy projection data, the algorithm with each parameter set decreases the evaluation function in the early iterations. However, the time course does not show a monotonic decrease in further iterations. Similar characteristics are known to exist and have been considered for the alternative MLEM [19] that is described as Equation (3) with *γ* = 1. We can see that a set of parameters (*<sup>γ</sup>*, *α*) exists at which the PDEM algorithm reduces the evaluation function more than MLEM does for any iteration number. Additionally, the smallest value of the evaluation function among the iteration numbers for a fixed set of the parameters becomes smaller with decreasing *γ* in the set {0.3, 0.5, 0.8, 1, 1.3} considered for this example. The parameter dependence of the evaluation function is clearly visible in contour plots of Figure 6, showing the values of log10(*E*(*z*(*N*))) for *N* = 50, 100, and 200 in the parameter plane. When designing a parameterized PDEM algorithm, a relatively large value of *α* and a small value of *γ* compared with the reference values of (*<sup>γ</sup>*, *α*)=(1, 1) provide the best performance in early and sufficient iterations, respectively. The best choices of (*<sup>γ</sup>*, *α*) depending on the termination iteration number are approximately (0.8, 1.2) at the 50th iteration, (0.5, 1.2) at the 100th iteration, and (0.3, 1.2) at the 200th iteration. The evaluation values under these conditions are indicated in Table 1, showing that PDEM with each parameter set gives a smaller value than MLEM does. The reconstructed images

and subtraction images (displayed in the range from 0 to 0.3) are shown in Figure 7. The figure reveals lots of artifacts in the reconstructed image due to the presence of noise in the measured projection. In terms of a quantitative evaluation, the structural similarity index measure [44] (SSIM) between the reconstructed and the true images was calculated and summarized in Table 2. A higher value of SSIM, which is a perception-based quality metric, provides higher image quality. By comparing the images reconstructed by MLEM and PDEM at the 100th and 200th iterations (see Figure 7 and Table 2), we can see that the PDEM with a proper set of parameters is able to reconstruct high-quality images while reducing the effects of noise in the projections, which means that PDEM is more robust to noise than MLEM.

**Figure 5.** Evaluation functions for MLEM and PDEM algorithms at each iteration in the experiment using numerical phantom with noisy projection.


**Table 1.** Values of the evaluation function for MLEM and PDEM with (*<sup>γ</sup>*, *α*) equal to (0.8, 1.2) at 50th iteration, (0.5, 1.2) at 100th iteration, and (0.3, 1.2) at 200th iteration in the experiment using numerical phantom with noisy projection.

**Figure 6.** Contour plots of evaluation functions log10(*E*(*z*(*N*))) with *N* equal to (**a**) 50, (**b**) 100, and (**c**) 200 in the experiment using numerical phantom with noisy projection. The white dot indicates the position of MLEM.

**Figure 7.** *Cont.*

**Figure 7.** Reconstructed images (**upper panel**) and subtraction images (**lower panel**) for MLEM and PDEM with (*<sup>γ</sup>*, *α*) equal to (**a**) (0.8, 1.2) at 50th iteration, (**b**) (0.5, 1.2) at 100th iteration, and (**c**) (0.3, 1.2) at 200th iteration in the experiment using numerical phantom with noisy projection.

**Table 2.** SSIM for MLEM and PDEM with the same parameters, as shown in Table 1 at *N*th iteration in the experiment using numerical phantom with noisy projection.


#### *4.2. Reconstruction Using Physical Phantom*

A physical experiment was carried out to further validate the effectiveness of the proposed method, although the true image is not available for a quantitative evaluation. The projections were physically acquired from an X-ray CT scanner (Canon Medical Systems, Tochigi, Japan) with a body-simulated phantom [45] (Kyoto Kagaku, Kyoto, Japan) using 80 kVp tube voltage, 30 mA tube current, and an exposure time of 0.75 s per rotation. Figure 8 represents the sinogram, a two-dimensional array of data containing the projections *y* ∈ *<sup>R</sup>I*+, with *I* = 430,200 (956 acquisition bins and 450 projection directions in 180 degrees) and a reconstructed image created by FBP using a Shepp–Logan filter with *J* = 454,276 (674 × 674 pixels). Images reconstructed by MLEM and PDEM with (*<sup>γ</sup>*, *α*)=(0.5, 1.2) are shown in Figure 9. The parameter values were referred to as the results of the numerical phantom with noisy projection. Figure 10, which shows the density profiles along horizontal lines (indicated by white) in the reconstructed images of Figures 8b and 9, verifies that the PDEM has a lower density deviation on a flat distribution

of the X-ray absorption in the physical phantom than either MLEM or FBP. The parameter values of the power exponents in the PDEM algorithm make it more robust to noise in spite of the higher noise level due to the low-dose X-ray exposure condition. This fact implies that the proposed method contributes to reducing patient doses during X-ray CT examinations in clinical practice by adjusting the parameter values depending on the noise levels of the projection data.

 **Figure 8.** (**a**) Sinogram and (**b**) reconstructed image by FBP in the experiment using physical phantom.

**Figure 9.** Reconstructed images for (**a**) MLEM and (**b**) PDEM with (*<sup>γ</sup>*, *α*)=(0.5, 1.2) at 200th iteration in the experiment using physical phantom.

**Figure 10.** Density profiles for (**a**) FBP, (**b**) MLEM, and (**c**) PDEM of reconstructed images along horizontal line with *L* = 674 × 224 and - = 1, 2, . . . , 674.

## **5. Conclusions**

We presented an extension of the PDM family with two parameters, *γ* and *α*, and proposed a novel iterative algorithm based on minimization of the divergence measure as an objective function of the reconstructed images. The theoretical results show the convergence of solutions to the continuous analog of the iterative algorithm owing to the objective function decreasing as the time proceeds. Numerical experiments illustrated that the proposed algorithm, which is considered to be an extended MLEM with two power exponents *γ* and *α*, has advantages over MLEM, which is the most popular and suitable iterative method of image reconstruction from noisy measured projections. The algorithm is of practical importance because its image quality is superior to that of MLEM. Our results sugges<sup>t</sup> that a larger value of *α* accelerates convergence and a smaller value of *γ* improves its robustness to measured noise. An investigation of the direct relation between the parameter variation in the EPDM family and the quality of images reconstructed by the proposed algorithm based on minimization of the EPDM is a future work to be considered. Moreover, we will use techniques such as machine learning to determine the most appropriate parameter depending on the noise level of the projections, number of projections, number of pixels, etc.

**Author Contributions:** Conceptualization, T.Y.; data curation, Y.Y. and T.Y.; formal analysis, R.K., Y.Y., O.M.A.A.-O. and T.Y.; methodology, R.K., Y.Y., T.K., O.M.A.A.-O. and T.Y.; software, R.K., Y.Y., T.K. and T.Y.; supervision, T.Y.; validation, O.M.A.A.-O. and T.Y.; writing—original draft, R.K., Y.Y. and T.Y.; writing— review and editing, R.K., Y.Y., T.K., O.M.A.A.-O. and T.Y. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was partially supported by JSPS KAKENHI, gran<sup>t</sup> number 21K04080.

**Data Availability Statement:** All data used to support the findings of this study are included within the article.

**Conflicts of Interest:** The authors declare no conflict of interest regarding the publication of this paper.
