Optimizing Parameters for Enhanced Iterative Image Reconstruction Using Extended Power Divergence

Abstract

In this paper, we propose a method for optimizing the parameter values in iterative reconstruction algorithms that include adjustable parameters in order to optimize the reconstruction performance. Specifically, we focus on the power divergence-based expectation-maximization algorithm, which includes two power indices as adjustable parameters. Through numerical and physical experiments, we demonstrate that optimizing the evaluation function based on the extended power-divergence and weighted extended power-divergence measures yields high-quality image reconstruction. Notably, the optimal parameter values derived from the proposed method produce reconstruction results comparable to those obtained using the true image, even when using distance functions based on differences between forward projection data and measured projection data, as verified by numerical experiments. These results suggest that the proposed method effectively improves reconstruction quality without the need for machine-learning techniques in parameter selection. Our findings also indicate that this approach is useful for enhancing the performance of iterative reconstruction algorithms, especially in medical imaging, where high-accuracy reconstruction under noisy conditions is required.

1. Introduction

X-ray computed tomography (CT) scanners and nuclear medicine CT scanners are technologies used to reconstruct tomographic images from projection data acquired through radiation transmission and emission, respectively. Image reconstruction methods can be broadly categorized into transform methods and iterative methods [1,2]. A transform method generates images analytically from projection data using Fourier transforms (or convolution integrals) and back-projection. When the CT inverse problem is consistent, an image equal to the true value of the inverse problem is obtained; however, when it is not consistent, artifacts may occur in the image. Meanwhile, iterative methods aim to reduce the differences between forward projections and measured projections through iterative calculations of pixel values. These methods are based on optimization principles and produce high-quality images even when the inverse problem cannot be analytically solved. Iterative methods include algebraic reconstruction technique [3], maximum-likelihood expectation-maximization (MLEM) [4], and multiplicative algebraic reconstruction technique [5], each with its unique characteristics. MLEM, in particular, is widely used in nuclear medicine diagnostic equipment due to its compatibility with noise models of nuclear disintegration and detection. In diagnostic X-ray CT, although iterative methods enable highly accurate image reconstruction, they are not mainstream due to the computational burden associated with the large number of pixels [6]. However, in recent years, iterative reconstruction has received much attention because it requires fewer radiation doses [6,7,8,9,10] while maintaining image quality. Additionally, as machine-learning techniques continue to evolve, there is a growing need for high-quality training data [11,12], which can be generated through advanced imaging methods. Therefore, the utility of high-performance iterative methods is expected to gain even more attention in the future.

In iterative reconstruction algorithms, convergence and image quality are influenced by the iterative calculations based on the objective function. Notably, the MLEM method is well-known for deriving its objective function from the Kullback–Leibler (KL) divergence [13]. Some of the authors previously proposed an iterative algorithm, the power divergence-based expectation-maximization (PDEM) method [14], which uses an extended power-divergence (EPD) measure as the objective function and demonstrated that it provides higher convergence performance and image quality than the traditional MLEM method. In the PDEM method, appropriately selecting the adjustable parameters in the iterative process is crucial for obtaining effective reconstruction results based on the projection conditions, which include factors such as the noise levels arising from the imaging conditions, as well as the concentration values of the object being scanned. Previous studies have proposed methods such as machine learning and sequential estimation for determining these parameters [15,16], but each has its pros and cons, and a complete solution has not yet been achieved. In this study, we propose a method to determine effective parameters before executing the reconstruction procedure by optimizing the evaluation function for projection data with common characteristics at a specified number of iterations.

In the numerical experiments, we examined projection noise and evaluated the results using the structural similarity index measure (SSIM) between the reconstructed images and the true images. The results revealed that using Euclidean distance or KL-divergence as distance functions did not yield satisfactory quality, but using EPD and weighted EPD (WEPD) resulted in high-quality reconstructed images. Although the Euclidean distance between forward and measured projections is commonly used for evaluating reconstruction performance, the results of this study suggest that this approach may need to be reconsidered. The effectiveness of the proposed method was also validated with projection data obtained from a physical X-ray CT scanner. The findings demonstrate the advantage of efficiently determining parameters in the preprocessing stage of reconstruction, allowing for high-quality images in iterative reconstruction using the PDEM method.

2. Image Reconstruction

Tomographic image reconstruction can be formulated as an inverse problem, where pixel values $x\in R_{+}^J$ are estimated from projections $y\in R_{+}^I$, following the model:

$$y=Ax+\sigma .$$

(1)

Here, $A\in R_{+}^{I\times J}$ is the projection operator, $\sigma \in R_{+}^I$ represents measurement noise, and $R_{+}$ denotes the set of non-negative real numbers. When all elements of the noise vector $\sigma$ are zero and Equation (1) is satisfied by x, the inverse problem is considered consistent. Hereafter, when the set

$$\mathcal{E}:=\left\{x\in R_{+}^J|y=Ax\right\}$$

(2)

is non-empty, we denote an element of this set by e.

A transform method based on the projection-slice theorem is commonly used when the tomographic inverse problem is consistent, i.e., if projections are measured over an angular range of 180 degrees, information on the 2D Fourier transform of the pixel values can be obtained, allowing for the reconstruction of pixel values from the projections through Fourier transform (or convolution integrals) and back-projection. However, in cases of ill-posed [17] inverse problems, such as when the measurements contain noise, certain adjustments are required, such as selecting a reconstruction filter (reconstruction function) that suppresses high frequencies instead of the standard ramp function used in convolution integrals.

Iterative methods are described by recursive algorithms, where desirable reconstructed images are obtained by solving an initial value problem of difference equations. Iterative algorithms have been proposed which vary depending on the objective function or the derivation principle of the tomographic inverse problem. In this paper, we define the MLEM and PDEM methods, which play an important role in our study.

2.1. MLEM Method

The MLEM algorithm, well-known for image reconstruction in nuclear medicine CT, is derived under the assumption that both the nuclear disintegration process and measurement noise follow a Poisson distribution. It has been proven that the generalized KL-divergence between the forward projection $Az\left(n\right)$ and the measured projection y,

$$\mathrm{KL}(y,Az\left(n\right)):=\sum _{i=1}^I{y}_{i}log{\displaystyle \frac{y_i}{A_{i}z\left(n\right)}}+A_{i}z\left(n\right)-y_i$$

(3)

decreases as the iteration n progresses [18,19]. Here, $z\left(n\right)$ represents the pixel values at iteration n, $A_i$ is the ith row vector of the projection operator A, and $y_i$ is the ith element of the projection vector y ($i=1,2,\dots ,I$). The iterative algorithm of the MLEM method can be described as the following recurrence relation for the elements $z_j$ of the pixel vector z ($j=1,2,\dots ,J$):

$$\begin{array}{ccc}& & z_j(n+1)=z_j\left(n\right){\displaystyle \frac{{\displaystyle \sum _{i=1}^I{A}_{ij}\frac{y_i}{A_{i}z\left(n\right)}}}{{\displaystyle \sum _{i=1}^I{A}_{ij}}}}\hfill \end{array}$$

(4)

for $n=0,1,2,\dots ,N-1$, where $A_{ij}$ denotes the element in the ith row and jth column of the matrix A ($i=1,2,\dots ,I$, $j=1,2,\dots ,J$), and the initial pixel values $z\left(0\right)$ are chosen as a constant vector $x^0$ with positive elements.

2.2. PDEM Method

Some of the authors extended the one-parameter family of the power-divergence measure [20,21,22,23] by introducing an additional parameter and applying it to the objective function in optimization, therefore deriving a new iterative reconstruction method [14]. Since the extended power-divergence (EPD) measure family includes KL-divergence, the PDEM method based on EPD is a natural extension of the MLEM method. It has been shown that by setting the two parameters in accordance with the level of measurement noise, the PDEM method can reconstruct high-quality images that are robust against noise, even under significant measurement noise. In other words, the proposed method outperforms conventional methods in terms of image quality and is particularly effective in reducing radiation exposure in X-ray CT diagnostics.

The PDEM method is an iterative algorithm that optimizes the following EPD measure family, which includes two parameters:

$${\Phi }_{\gamma ,\alpha }(y,Az):=\sum _{i=1}^I{\int }_{y_i}^{A_{i}z}{\displaystyle \frac{s^{\gamma }-y_i^{\gamma }}{s^{\gamma \alpha }}}ds$$

(5)

and is described by the following recurrence relation:

$$\begin{array}{ccc}& & z_j(n+1)=z_j\left(n\right){\displaystyle \frac{{\displaystyle \sum _{i=1}^I{A}_{ij}{\left(\frac{y_i}{{\left(A_{i}z\left(n\right)\right)}^{\alpha }}\right)}^{\gamma }}}{{\displaystyle \sum _{i=1}^I{A}_{ij}{\left(\frac{A_{i}z\left(n\right)}{{\left(A_{i}z\left(n\right)\right)}^{\alpha }}\right)}^{\gamma }}}}\hfill \end{array}$$

(6)

for $n=0,1,2,\dots ,N-1$, where $j=1,2,\dots ,J$, and the initial value is set to $z\left(0\right)=x^0$. The parameters $\gamma >0$ and $\alpha \ge 0$ are included in Equation (6). The PDEM algorithm in Equation (6) coincides with the MLEM algorithm in Equation (4) when $\gamma =1$ and $\alpha =1$. Furthermore, when $\alpha =1$, it becomes an iterative algorithm with the generalized Hellinger distance [24] as the objective function. Since Equation (5) encompasses many distance or divergence measures, including KL-divergence, Neyman’s ${\chi }^2$ distance, and generalized Hellinger distance, the parameter family derived through the optimization of Equation (5) is expected to exhibit diverse performance characteristics for image reconstruction. In fact, experiments have demonstrated that by selecting the appropriate parameters, the PDEM method can exhibit robustness against measurement noise and achieve higher-quality reconstructed images than the MLEM method [14]. It is important to note that the degree of improvement provided by PDEM over MLEM is likely to be dataset-dependent. Specifically, the advantage of PDEM becomes more pronounced when dealing with noisy data or under conditions where the number of projections is limited. In cases where the data quality is high, and noise is minimal, the primary advantage of increasing the parameter values in PDEM lies in achieving faster computation compared to MLEM, with no significant improvement in image quality.

3. Proposed Method

The proposed method focuses on iterative reconstruction algorithms with adjustable parameters and provides a way to explore parameter values that yield the most effective performance for a set of projection data sharing common characteristics, such as similar noise levels and internal structures.

First, given K projections $Y\in R_{+}^{I\times K}$ with the same number of elements I, we define an iterative algorithm with the target system f to determine the corresponding pixel values $Z\in R_{+}^{J\times K}$ as follows:

$$\begin{array}{ccc}& & Z_k(n+1)=f(Z_k\left(n\right);\lambda )\hfill \end{array}$$

(7)

for $n=0,1,2,\dots ,N-1$ and $k=1,2,\dots ,K$. Here, $\lambda \in R_{+}^M$ represents the parameters to be adjusted. For example, when the target f is the PDEM algorithm used in this paper, $M=2$ and

$$\lambda =\left(\begin{array}{c}\gamma \\ \alpha \end{array}\right).$$

(8)

Next, starting from the initial value $Z_k\left(0\right)$ ($k=1,2,\dots ,K$) according to Equation (7), after n iterations ($n=1,2,\dots ,N$), we define

$$Z_k\left(n\right)=:f^n(Z_k\left(0\right);\lambda )$$

(9)

for $k=1,2,\dots ,K$. We then seek the parameters that minimize the evaluation function

$$g(\lambda ;N):={\displaystyle \frac{1}{K}}\sum _{k=1}^{K}d\left(f^N(Z_k\left(0\right);\lambda )\right)$$

(10)

where d is the distance function, defined below, evaluated at the Nth iteration. Our goal is to find the elements of the set

$$\underset{\lambda }{\mathrm{argmin}}g(\lambda ;N).$$

(11)

Note that the function g depends on $\lambda$ through the computation process of $Z\left(N\right)$, which is obtained by iterating $f(·;\lambda )$, where f includes $\lambda$ as a parameter.

In the experiments conducted in this study, the distance function d is defined as follows. First, the Euclidean distance using the true value $e\in \mathcal{E}$ of the image is defined as:

$$d\left(x\right):=\Vert e-x\Vert .$$

(12)

Next, for the distance between the forward projection and the measured projection $Y_k$ ($k=1,2,\dots ,K$), we consider the following four types:

• Euclidean distance:

  $$d\left(x\right):=\Vert Y_k-Ax\Vert ,$$

  (13)
• KL-divergence:

  $$d\left(x\right):=\mathrm{KL}(Y_k,Ax),$$

  (14)
• EPD:

  $$d\left(x\right):={\Phi }_{0.1,0.1}(Y_k,Ax)$$

  (15)

  where

  $${\Phi }_{{\gamma }^0,{\alpha }^0}(p,q)=\sum _{i=1}^I{\int }_{p_i}^{q_i}{\displaystyle \frac{s^{{\gamma }^0}-p_i^{{\gamma }^0}}{s^{{\gamma }^0{\alpha }^0}}}ds,$$

  (16)
• WEPD:

  $$d\left(x\right):={\Psi }_{0.1,0.1}(Y_k,Ax)$$

  (17)

  where

  $${\Psi }_{{\gamma }^0,{\alpha }^0}(p,q)=\sum _{i=1}^I{\int }_{p_i}^{q_i}{\displaystyle \frac{s^{{\gamma }^0}-p_i^{{\gamma }^0}}{s^{{\gamma }^0{\alpha }^0}}}ds\sum _{j=1}^J{A}_{ij}.$$

  (18)

Equation (12) is used as a distance function to accurately evaluate the quality of the reconstructed image, while Equations (13)–(15), and (17) serve as candidate functions for applying the proposed method to actual projections. The parameters ${\gamma }^0$ and ${\alpha }^0$ in Equations (16) and (18) are both set to $0.1$ in Equations (15) and (17), respectively, based on the results from Ref. [14], which suggest that smaller parameter values are more desirable for achieving robustness against noise, and fixed values ensure fairness without dependency on other conditions. The initial values for optimization associated with the set in (11) are assigned to $\gamma \left(0\right)=1$ and $\alpha \left(0\right)=1$. Note that the iterative rule with these fixed parameters coincides with MLEM.

4. Experiments

Through numerical and physical experiments, we verify whether the proposed method operates as expected and discuss effective distance functions for evaluating the comparison between forward and measured projections. To demonstrate the advantages of iterative methods over transformation methods, we set conditions with fewer projection directions relative to the number of pixels and with a relatively high level of noise.

4.1. Experimental Method

We consider the adjustment of parameters $\lambda =(\gamma ,\alpha )$ for the PDEM iterative algorithm in Equation (6).

In the numerical experiments, we simulate a first-generation CT scan with 180-degree rotation and reconstruct images ($J=128\times 128$) from projections with 184 bins and 90 projection directions ($I=184\times 90$). We use the disc or modified Shepp–Logan phantom shown in Figure 1 as the true image $e\in R_{+}^J$. For each k of $K=8$ projections with the same noise level but different noise patterns, we create projections using

$$Y_k=Ae+{\sigma }_k\in R_{+}^I$$

(19)

for $k=1,2,\dots ,K$. Here, the normally distributed random noise ${\sigma }_k$ is set to have a variance such that the signal-to-noise ratio of the projection is 20 dB and is independently generated for each k. This choice is justified because the noise in X-ray CT is primarily quantum noise and electronic noise, which are well-represented by a normal distribution. The intention is to seek parameters that are common to the average projections, which share the same uncertainty in terms of signal-to-noise ratio, rather than focusing on a specific random variable. The number of iterations N for optimization is set to 30.

In the physical experiments, projection data were acquired using an X-ray CT scanner (Canon Medical Systems, Tochigi, Japan) and a phantom simulating the human body [25]. To investigate the noise suppression effect, which is a characteristic of the PDEM reconstruction method, we used an abdominal phantom with relatively simple density variations. The scanner was operated with a tube voltage of 80 kVp, tube current of 10 mA, and scan time of 0.75 s per rotation. This setup corresponds to imaging conditions under relatively high levels of measurement noise. Projections were obtained from 900 directions covering 360 degrees with 957 bins, and data were subsampled to achieve 90 projection directions ($I=\mathrm{86,130}$) to simulate sparse projections. Optimization was performed using projection data from four consecutive slices ($K=4$) and an evaluation function with $N=100$ iterations. The reconstructed image size was $674\times 674$ ($J=\mathrm{454,276}$). Figure 2 shows images created from one of the projection data slices using the filtered back-projection (FBP) method with a Shepp–Logan filter and the iterative reconstruction method using MLEM.

We solved the problem of finding the elements of the set in (11) using the solver fmincon provided in MATLAB R2023a (MathWorks, Natick, USA) for finding the minimum of constrained nonlinear multivariable functions. The initial parameter values were set to $\lambda \left(0\right)=(1,1)$. The lower bound for parameter variation was set to $0.01$ and the upper bound to $10.0$. Although a relatively large value is set as the exponent for the upper limit necessary for numerical accuracy, the internal state is also checked in the numerical experiments to ensure that no illegal values or divergence occurs.

For the numerical calculations, we used a computer equipped with a 12-core Apple M2 Max processor and 64 GB of memory, and the program utilized parallel computing capabilities.

4.2. Experimental Results

First, we present the results of numerical experiments when the distance function is chosen as the Euclidean distance between the reconstructed image and the true image (Equation (12)). Figure 3 shows the evolutions of the parameters and the evaluation function g obtained through solving our optimization problem associated with the set in (11) for two types of phantoms. The horizontal axis represents the iterations of the optimization algorithm $t=0,1,2,\dots$, while the green points represent the parameters, with their values referable on the left vertical axis. The brown points show the evaluation function $g\left(\lambda \right(t);N)$ along with its values on the right vertical axis. The algorithm stops when either the maximum number of iterations is reached, or the parameter changes fall below a predetermined tolerance, with the iteration count at which it stops denoted as T. Regardless of the phantom type, $\lambda \left(t\right)$ changes due to the iteration process, and $g\left(\lambda \right(t);N)$ decreases monotonically. However, the parameter set $\lambda \left(T\right)$ at iteration T differs by phantom type: $(0.17,3.45)$ for the disc and $(0.39,2.48)$ for the Shepp–Logan phantom. Figure 4 shows the iteration process of the evaluation function values when the initial parameter $\lambda \left(0\right)$ and the final parameter $\lambda \left(T\right)$ are given to the PDEM iterative system, represented by blue and red points, respectively. Since $\lambda \left(0\right)=(1,1)$ was chosen, the blue points provide results identical to the MLEM method. The function values $g\left(\lambda \right(0);N)$ and $g\left(\lambda \right(T);N)$ at $N=30$ match the values shown in Figure 3 at the initial iteration $t=0$ and final iteration $t=T$, respectively. Figure 5 shows the images reconstructed using MLEM and PDEM at iteration N and the corresponding subtraction images from the true values, where the parameter for PDEM is $\lambda \left(T\right)$. The density profiles along the column direction of the images are plotted in Figure 6.

Next, we present the results of numerical experiments when the distance function is chosen as the Euclidean distance (Equation (13)), KL-divergence (Equation (14)), EPD (Equation (15)), and WEPD (Equation (17)) between forward and measured projections. Figure 7 shows the evolution of the parameters obtained from the optimization associated with the set in (11) and the transition of the evaluation function through iterations, with the graph representation following the same format as in Figure 3. In each of the four distance functions, the evolution of the evaluation function $g\left(\lambda \right(t);N)$ decreases monotonically, demonstrating that the proposed algorithm operates as expected. Figure 8 illustrates the iterative process of the evaluation function value when the initial value $\lambda \left(0\right)$ and the final value $\lambda \left(T\right)$ are applied to the PDEM system. The reconstructed images at iteration N are shown in Figure 9.

Table 1. Performance evaluation of images reconstructed using MLEM and PDEM (with parameters $\lambda \left(T\right)$) for disc phantom. The distance functions are (a) Euclidean distance, (b) KL-divergence, (c) EPD, and (d) WEPD between forward and measured projections, and (e) Euclidean distance between the reconstructed and true images.

Method	Dist. Func.	$\mathit{\gamma }$	$\mathit{\alpha }$	SSIM	Std. Dev.
MLEM	—	1	1	0.816	0.111
	(a)	1.219	0.744	0.813	0.122
	(b)	0.957	1.39	0.798	0.118
PDEM	(c)	0.139	9.21	0.983	0.0719
	(d)	0.143	9.17	0.982	0.0727
	(e)	0.166	3.45	0.981	0.0517

and

Table 2. Performance evaluation of images reconstructed using MLEM and PDEM (with parameters $\lambda \left(T\right)$) for Shepp–Logan phantom. The distance functions are (a) Euclidean distance, (b) KL-divergence, (c) EPD, and (d) WEPD between forward and measured projections, and (e) Euclidean distance between the reconstructed and true images.

Method	Dist. Func.	$\mathit{\gamma }$	$\mathit{\alpha }$	SSIM
MLEM	—	1	1	0.900
	(a)	1.34	0.764	0.889
	(b)	0.959	1.40	0.889
PDEM	(c)	0.156	8.31	0.928
	(d)	0.169	7.86	0.932
	(e)	0.393	2.48	0.939

summarize the parameter values $\lambda \left(T\right)$, SSIM, and standard deviation (for the disc phantom) between the reconstructed images and the true images for each distance function when the optimization of the evaluation function reaches iteration T.

Finally, we present the results of the physical experiment. Figure 10 plots the optimization process and the transition of the evaluation function when the proposed method is applied to projections obtained from the X-ray CT scanner, using WEPD as the distance function between forward and measured projections as an example. The reconstructed images and density profiles, including those for other distance functions, are shown in Figure 11 and Figure 12, respectively.

5. Discussion

First, we discuss the results of the numerical experiments that enable quantitative evaluation using the true values. The algorithm applied to the PDEM system, as shown in Figure 3 and Figure 7, demonstrates that it operates appropriately as expected, regardless of whether the distance function compares images or projections or the type of function that was used. The parameters $\lambda \left(T\right)$ obtained using the true image significantly reduce the evaluation function compared to $\lambda \left(0\right)$, as seen in Figure 4. Additionally, visual assessment from Figure 5 and the density profiles in Figure 6 indicate that the reconstructed images exhibit less variability in density values, which should be flat, providing high-quality results. In practical applications, using distance functions based on forward and measured projections is necessary. It would be desirable if projections yield the same results as functions based on images. Although Euclidean distance and KL-divergence are commonly used as distance functions related to projections, effective parameters have not been obtained, as shown by the changes in distance functions (Figure 8a,b). No improvement in image quality was observed from MLEM (Figure 9). In contrast, distance functions using EPD and WEPD significantly reduce the evaluation function compared to MLEM (Figure 8c,d), and visual assessments (Figure 9) demonstrate that the reconstructed images are high quality.

Table 1 summarizes the results of the numerical experiments for the disc phantom. Focusing on the parameter values, previous experimental studies [14] have shown that to improve the image reconstruction quality with PDEM, it is preferable to set a large value for $\alpha$ to quickly reduce the evaluation function related to the true image with fewer iterations, and to set a small value for $\gamma$ to prevent the evaluation function from increasing in subsequent iterations. When the distance functions related to projections were selected as EPD (Equation (15)) and WEPD (Equation (17)), the values of the parameter sets $\lambda \left(T\right)$ were $(0.139,9.21)$ and $(0.143,9.17)$, respectively. These values, like the parameter set $(0.166,3.45)$ obtained when using the Euclidean distance for images (Equation (12)), resulted in a sufficiently small $\alpha$ and a sufficiently large $\gamma$, which clearly differs from the case when the Euclidean distance (Equation (13)) or KL-divergence (Equation (14)) was chosen for projections. Due to these parameter estimates, the quantitative evaluation of the reconstructed images also shows that, for projection-based EPD and WEPD distance functions, the results are of the same or even higher quality than those obtained with the image-based distance function. This is similarly observed in the results for the Shepp–Logan phantom in Table 2. However, in the experiment with the phantom, as shown in Figure 5b and Figure 6b, the edges of structures where the density values change sharply are blurred. The PDEM shows favorable results in the objective function and SSIM because the noise suppression effect outweighs the smoothing performance. Equation (16) indicates that the parameter sets $({\gamma }^0,{\alpha }^0)$ corresponding to $(1,0)$ and $(1,1)$ align with the square of the Euclidean distance and KL-divergence, respectively. In this study, the parameter set $(0.1,0.1)$ chosen for the distance functions EPD and WEPD is considered effective because it makes the reconstruction more robust to noise. The effectiveness of the distance function related to EPD may be due to the fact that PDEM is a reconstruction method based on the optimization of EPD. However, it is of significance that we demonstrated the use of a non-Euclidean distance function, rather than the commonly used Euclidean distance, for evaluating reconstruction methods, including MLEM, as a function of forward and measured projections.

As discussed, the quality of reconstruction varies significantly depending on the choice of distance function, with EPD and WEPD providing higher quality than the conventional Euclidean distance and KL-divergence. The experiments conducted in this study were designed under conditions of high projection noise and fewer projection directions than the number of pixels. This setting was motivated by the aim of proposing methods for reducing radiation exposure. It is known that with low noise levels, providing a high value of $\gamma$ to PDEM yields higher-quality images than MLEM [14]. Moreover, the advantages of using iterative methods diminish when noise is low. It should be noted that, under low noise conditions, the outcomes produced by the proposed method, depending on the distance function, may vary from the results presented in this paper.

6. Conclusions

In this study, we proposed a method of optimizing the evaluation function to derive the parameter values that yield the highest performance for systems involving iterative methods with adjustable parameters. Numerical experiments were conducted for PDEM, optimizing two power exponents included in the iterative rule as adjustable parameters. The results indicated that optimizing the evaluation function—specifically, using EPD and WEPD as the distance functions defined by the proposed algorithm—yields the most effective estimates. One of the key results of this study is that the proposed algorithm could yield similar performance to that of distance functions based on the reconstructed image and the true image despite using distance functions based on forward and measured projections. Physical experiments verified that results similar to those of the numerical experiments could be obtained, demonstrating that the proposed method contributes to improving image quality in iterative reconstruction using PDEM.