Iterative Image Reconstruction Algorithm with Parameter Estimation by Neural Network for Computed Tomography

Abstract

Recently, an extended family of power-divergence measures with two parameters was proposed together with an iterative reconstruction algorithm based on minimization of the divergence measure as an objective function of the reconstructed images for computed tomography. Numerical experiments on the reconstruction algorithm illustrated that it has advantages over conventional iterative methods from noisy measured projections by setting appropriate values of the parameters. In this paper, we present a novel neural network architecture for determining the most appropriate parameters depending on the noise level of the projections and the shape of the target image. Through experiments, we show that the algorithm of the architecture, which has an optimization sub-network with multiplicative connections rather than additive ones, works well.

1. Introduction

Computed tomography (CT) can be thought of in terms of a linear inverse problem wherein one estimates tomographic images from a measured projection and a known projection operator [1,2,3,4,5,6]. Although filtered back-projection, which is an analytical method, is widely used for image reconstruction in medical CT scanners, iterative image reconstruction [2,5,7,8,9] better meets modern requirements for practical use. The iterative image reconstruction method is based on an optimization algorithm and can obtain an image with relatively few artifacts from noisy projection data. The maximum-likelihood expectation-maximization (MLEM) algorithm [2], a typical iterative image reconstruction algorithm, is known to minimize the Kullback–Leibler (KL) divergence [10] of the measured and forward projections.

Since the objective function of a minimization problem for CT affects the quality of the reconstructed images, a rational design of the function is essential. Our research group previously presented a two-parameter extension of the power-divergence measure (PDM) [11,12,13,14], called EPDM, and proposed an iterative image reconstruction method [15] based on minimization of EPDM as an objective function, i.e., the PDEM algorithm. Note that, because the KL-divergence is a special case of the EPDM family, PDEM is a natural extension of MLEM. The PDEM algorithm can produce a high-quality reconstructed image by adjusting the values of the parameters depending on the degree of noise in the measured projection.

PDEM is derived from a nonlinear differential equation using the dynamical method [16,17,18,19,20,21] applied to tomographic inverse problems [22,23,24,25,26,27], which means that a discretization of the differential equation by using a first-order expansion of the vector field leads to an equivalent iterative formula. A theoretical result from the previous study [15] is that the global stability of an equilibrium corresponding to the desired solution observed in the continuous-time system is guaranteed by the Lyapunov theorem [28]. Furthermore, numerical experiments have shown that the appropriate parameters depend on the level of projection noise. As a natural conjecture, the appropriate parameters also depend on the shape of the image to be reconstructed, so one may consider whether it is possible to develop a method for estimating the parameters from the measured projection.

In recent years, many studies have incorporated deep-learning approaches for parameter estimation in iterative algorithms [29,30,31,32,33,34,35]. The well-known methods of ADMM-Net [32] and ISTA-Net [34] embed an iterative algorithm into a deep neural network by using the alternating direction method of multipliers [36] and the iterative shrinkage thresholding algorithm [37], respectively. Some of these studies have used recurrent neural networks (RNNs) in which convolutional sub-networks for estimating parameters are used in each iterative calculation and the iterative calculations are intertwined [29,33]. The complexity of the networks increases in proportion to the number of iterative calculations, but the risk of the vanishing gradient problem [38] occurring is low because the calculations are implemented as residual blocks [39]. The ideas behind this neural network architecture (referred to below as learned optimization; LO) have been organized by Andrychowicz et al. [40], and LO has been actively improved in terms of its design, training method, performance, and operation [41,42,43,44,45,46,47]. In particular, this technique has been applied to the field of medical imaging, especially the problem of compressed sensing in magnetic resonance imaging [48,49,50]. Zhang et al. [51] used it together with Nesterov’s acceleration [52] to estimate parameters associated with the regularization term for compensating for a lack of signals.

In this paper, we propose a novel neural-network architecture that estimates the values of parameters in the PDEM algorithm by using projection data. The proposed architecture is similar to LO because it has two subsystems for estimating parameters and iterative operations. Following the naming convention of Andrychowicz et al., we call these subsystems the optimizer and optimizee, respectively. Unlike conventional LO, the subsystems are completely separated in our system, and the network architecture is more straightforward than that of RNNs. The traditional LO estimates parameters by using the results of each iterative calculation, whereas our architecture uses only projection data. Nevertheless, the experimental results illustrate that our method works well. Moreover, the optimizee has a different structure than the traditional optimizee, with multiplicative rather than additive connections.

2. Problem Description

We assume that the measured projection $y\in R_{+}^I$ in computed tomography is acquired according to the linear model

$$y=Ax+\delta$$

(1)

where $A\in R_{+}^{I\times J}$, $x\in R_{+}^J$, and $\delta \in R^I$ are a projection operator, target image, and noise, respectively, with $R_{+}$ denoting the set of nonnegative real numbers. The purpose of iterative reconstruction is to obtain a high-quality image similar to the target image by solving an initial value problem consisting of a difference equation or an equivalent iterative algorithm. Recently, a novel iterative image reconstruction algorithm [15], which is based on minimization of the EPDM as an objective function of the measured and forward projections, has been developed:

$$\begin{array}{c}{z}_j(k+1)=z_j\left(k\right)f_j(z\left(k\right);y,\lambda )\end{array}$$

(2)

with

$$\begin{array}{c}{f}_j(z;y,\lambda ):=\frac{{\displaystyle \sum _{i=1}^I{A}_{ij}{\left(\frac{y_i}{{\left(A_{i}z\right)}^{\alpha }}\right)}^{\gamma }}}{{\displaystyle \sum _{i=1}^I{A}_{ij}{\left(\frac{A_{i}z}{{\left(A_{i}z\right)}^{\alpha }}\right)}^{\gamma }}}\end{array}$$

(3)

for $j=1,2,\dots ,J$ and $k=0,1,2,\dots ,K-1$ with $z\left(0\right)=x^0\in R_{+}^J$, where $A_i$ and $A_{ij}$ indicate the ith row vector and the element in the ith row and jth column of the matrix A, respectively, and $\gamma >0$ and $\alpha \ge 0$ are parameters that are collectively referred to as $\lambda :=\left(\gamma ,\alpha \right)$. The PDEM algorithm described in Equation (2) can be rewritten as

$$z(k+1)=z\left(k\right)\odot f\left(z\right(k);y,\lambda ),$$

(4)

with $f:={(f_1,f_2,\dots ,f_J)}^{\top }$, where the symbol ⊙ denotes the Hadamard product and the superscript ⊤ stands for the transpose.

The purpose of this paper is to construct a system that estimates the most appropriate values of $\lambda$ in the PDEM algorithm from the projection y by using a deep-learning approach.

3. Proposed Method

Our goal is to design a system for estimating an appropriate set of parameters $\lambda$ from the projection y by using supervised learning on a deep neural network, described by

$$\mathrm{\Lambda }:R_{+}^I\to R^2;y\mapsto \lambda .$$

(5)

Since it is complicated to prepare a large set of appropriate parameters as a teaching dataset, we chose to use the following network architecture,

$$\begin{array}{c}\mathrm{\Theta }\left(y\right):=\mathrm{\Omega }(y,\mathrm{\Lambda }(y\left)\right)\end{array}$$

(6)

where

$$\begin{array}{c}\mathrm{\Omega }(y,\lambda ):=\left(\underset{K}{\underbrace{p\circ p\circ \cdots \circ p}}\right)\left(x^0\right)\end{array}$$

with

$$\begin{array}{c}p\left(z\right):=z\odot f(z;y,\lambda )\end{array}$$

as illustrated in Figure 1, and convert the neural network’s output into a reconstructed image. This architecture allows us to use true images or high-quality reconstructed images as teaching data. The subsystem $\mathrm{\Omega }$ as an optimizee produces a reconstructed image $z\left(K\right)$ by iterating the formula in Equation (2) K times, as shown in Figure 2. Note that it does not include any training parameters. The whole system $\mathrm{\Theta }$ is implemented as a deep neural network and the training parameters of the subsystem $\mathrm{\Lambda }$ are automatically adjusted through the backpropagation procedure.

Although the depth of the computational graph of $\mathrm{\Theta }$ increases in proportion to the number of iterations K, the optimizee $\mathrm{\Omega }$ can be considered to be the residual blocks that can then be configured with multiplicative connections, as shown in Figure 2 and Figure 3, which is expected to reduce the risk of the vanishing gradient problem occurring. The well-known residual blocks can be described in terms of a difference equation

$$w_j(k+1)=w_j\left(k\right)+g_j(k,w\left(k\right)),$$

(7)

for $j=1,2,\cdots ,J$ and $k=0,1,\dots$, which is obtained by applying the (additive) Euler method with a step size of 1 to the differential equation

$$\frac{d{v}_j\left(t\right)}{dt}=g_j(t,v\left(t\right))$$

(8)

with $w_j\left(0\right)=v_j\left(0\right)$. On the other hand, applying the multiplicative Euler method [53] to Equation (8) yields the difference equation

$$u_j(k+1)=u_j\left(k\right)exp\left(\frac{g_j(k,u\left(k\right))}{u_j\left(k\right)}\right)$$

(9)

for $j=1,2,\cdots ,J$ and $k=0,1,\dots$ with $u_j\left(0\right)=v_j\left(0\right)$. Both residual blocks have shortcut edges that may help reduce the risk of the vanishing gradient problem. Indeed, the system in Equation (2) can be derived by discretizing the differential equation

$$\begin{array}{c}\frac{d{x}_j\left(t\right)}{dt}=x_j\left(t\right)log\left(f_j(x\left(t\right);y,\lambda )\right),\end{array}$$

(10)

for $j=1,2,\dots ,J$ and $t\ge 0$ with $x\left(0\right)=x^0$ by using the multiplicative Euler method. The differential equation in Equation (10) is guaranteed by the Lyapunov theorem to converge globally to the desired solution when the algebraic equation in Equation (1) with $\delta =0$ has a solution.

4. Results and Discussion

In this section, we present the results of an experiment on the proposed system using 5000 pairs of input and teaching data generated by numerical simulations. We generated teaching data in the form of $64\times 64$-pixel images deformed by changing the attributes of each oval, i.e., its intensity, size, position, and rotation angle, from the Shepp–Logan phantom [54] shown in Figure 4a. The attributes were varied in the ratio ranges shown in

Table 1. Deformation ratios from base phantom. Symbol * indicates that the deformation ratio for the attribute uses a common value for each oval in the same group. See Figure 4b for “oval numbers”.

Group	Oval Numbers	Intensity	Minor Axis Length	Major Axis Length	Horizontal Position	Vertical Position	Angle of Rotation
1	#1, #2	$\pm 10\%$	$\pm 20\%$ *	$\pm 15\%$ *	$\pm 0\%$	$\pm 0\%$	$\pm 0\%$
2	#3, #4, #5	$\pm 30\%$	$\pm 30\%$	$\pm 20\%$	$\pm 25\%$	$\pm 25\%$	$\pm 20\%$
3	#6, #7	$\pm 30\%$	$\pm 50\%$	$\pm 50\%$	$\pm 50\%$	$\pm 50\%$	$\pm 50\%$
4	#8, #9, #10	$\pm 30\%$ *	$\pm 25\%$ *	$\pm 25\%$ *	$\pm 20\%$ *	$\pm 20\%$ *	$\pm 50\%$

from that of the base phantom by using uniform random numbers. The MATLAB (MathWorks, Natick, MA, USA) function phantom [55] generated the images. To avoid generating extremely large or small ovals that should not be considered, some common deformation rates were chosen for one attribute in some of the sets of ovals in one phantom. The input data were projections calculated by Equation (1) by using a common projection operator A, projection noise $\delta$, and the target image x as the corresponding teaching data. The projection operator A, represented as a matrix, was simulated on the basis of the Radon transform assuming parallel projections from 90 directions with 95 detectors, so $I=90\times 95$ and $J=64\times 64$. The signal-to-noise ratio (SNR) for each projection was selected so as to be uniformly distributed in the range of 20 dB to 50 dB. Figure 5 shows the resulting experimental histogram of SNRs, while Figure 6 shows samples of learning data generated by the above procedure.

In the experiment, we designed the optimizer $\mathrm{\Lambda }$ to be a simple and typical network (Figure 7). The batch normalization layer [56] is supposed to stabilize learning, reduce the initial value dependence of the learning parameters, and speed up learning. The dropout layer [57] is used to suppress overfitting.

Since the first hidden layer is dense, it requires a massive number of training parameters, resulting in a vast computational cost for training. Therefore, designing a more efficient optimizer $\mathrm{\Lambda }$ is a future challenge.

The parameters $\gamma$ and $\alpha$ are power exponents, and giving them excessive values makes the numerical operations unstable. However, in the process of automatically updating the learning parameters by using the backpropagation procedure, $\gamma$ and $\alpha$ are given large values that temporarily exceed the tolerance of a stable calculation . We avoided this problem by setting an upper limit on the output values of the optimizer $\mathrm{\Lambda }$.

We define the activation function of the output layer of $\mathrm{\Lambda }$ to be

$$h(x_1,x_2):=1.8\times {(\mathrm{sgm}\left(x_1\right),\mathrm{sgm}\left(x_2\right))}^{\top },$$

(11)

where the function $\mathrm{sgm}$ is the standard sigmoid function. The coefficient of 1.8 for the sigmoid function, given as an upper bound for both parameters, is small enough not to cause numerical instability or affect the parameter estimation.

The learning conditions and hyperparameters were as follows. We used 4000 pairs of data for training and 1000 pairs different from them for validation. For the optimizee $\mathrm{\Omega }$ responsible for image reconstruction, the number of iterations was $K=200$ and an initial value was $x^0={(0.5,0.5,\dots ,0.5)}^{\top }$. Note that the value of the maximum iteration number K was chosen to ensure a sufficient reconstruction quality and to properly observe the difference in the convergence characteristics between the PDEM and MLEM algorithms, as shown in the previous study [15]. The number of epochs and the batch size were 5000 and 64, respectively. We used the squared L2 norm as the loss function. The Adam algorithm [58] was used to minimize the loss function; its learning rate was set to ${10}^{-6}$. The learning network was implemented using Python 3.9.7 and TensorFlow 2.8.0.

4.1. Results of Learning

Here, we present the experimental results of learning and show that the proposed network successfully learned the optimizer $\mathrm{\Lambda }$.

Figure 8 shows the values of the loss function at each epoch. At the final epoch (the 6000th), the loss function for the training data shows a decreasing trend, while the one for the validation data converged, indicating that the learning was sufficient.

To evaluate the performance of the learned optimizer $\mathrm{\Lambda }$, we validated it on a test dataset consisting of data not included in the training and validation datasets. The properties described below are common to all of the test datasets. We will discuss the results for five of the test data pairs shown in Figure 9. The first test data item ($\ell =1$) was a pair of a projection $y^1$ with an SNR of 20 dB and the expected reconstructed image $e^1$. In the second and third test data pairs ($\ell =2$ and 3), the expected reconstructed images $e^2$ and $e^3$ equaled $e^1$ but with different Sn Rs of the projection. On the other hand, in the fourth and fifth test data pairs ($\ell =4$ and 5), the expected reconstructed images $e^4$ and $e^5$ were different from $e^1$ and from each other, but the noise levels of the projections were 20 dB.

Figure 10 plots the evaluation function values

$$D(\ell ,z):=\parallel e^{\ell }{-z\parallel }_2$$

(12)

defined by the L2 norm of the difference between the ℓth true image $e^{\ell }$ and reconstructed image z when using the parameter ${\widehat{\lambda }}^{\ell }=({\widehat{\gamma }}^{\ell },{\widehat{\alpha }}^{\ell })$ estimated by the learned model $\mathrm{\Lambda }$ from the projection $y^{\ell }$ with $\ell =1,2,\cdots ,5$. The results for the MLEM algorithm, i.e., the PDEM algorithm with $\lambda =(1,1)$, are displayed for comparison.

For the first test data pair ($\ell =1$), due to the high noise level, the evaluation functions for both the MLEM and the proposed algorithm, which is the PDEM using the parameter values ${\widehat{\lambda }}^1=(0.40,1.05)$ estimated by the optimizer $\mathrm{\Lambda }$, do not decrease monotonically but rather start to increase at different iteration numbers, around $k=50$ and 140, respectively. The proposed method decreases the evaluation function significantly in the last iteration compared with that of the MLEM algorithm. This result indicates that the optimizer $\mathrm{\Lambda }$ can estimate parameter values that, while not optimal, significantly improve image quality. The results for the second test data pair ($\ell =2$), which has a different noise level than the first test data pair ($\ell =1$), shows that the proposed method also reduces the evaluation function compared with the MLEM algorithm, but the difference is slight. We discuss the reason in the following subsection. The proposed method can significantly reduce the evaluation function when using the third test data pair ($\ell =3$), which has weak noise. The reason for this is presumably that the power exponent parameter $\gamma$ is a value that affects the step size in the discretization of the continuous dynamical system and that, when the noise is weak, a relatively large parameter is chosen for fast convergence to a value close to the true one. In fact, the estimated parameter ${\widehat{\gamma }}^3$ is 1.64, which is larger than the others. The results for the fourth and fifth test data pairs ($\ell =4$ and 5), where the projection contains a lot of noise, show a slight decrease in the value of the evaluation function, contrary to the results for the first test data pair ($\ell =1$), where the SNRs of the projections are comparable. This result suggests that the optimal parameters also depend on the shape of the reconstruction target.

Figure 11 shows the contour maps of the evaluation function values of the reconstructed image $z\left(K\right)$ for changes in the parameters $\gamma$ and $\alpha$ according to the projection $y^{\ell }$. On all the test data, the optimizer $\mathrm{\Lambda }$ could estimate parameters that reduce the evaluation function value more than the MLEM algorithm. In the case of $\ell =1$ and $\ell =3$, where significant improvements were observed, the region of suitable parameters is located far from $(\gamma ,\alpha )=(1,1)$ corresponding to the parameter for MLEM, and the optimizer $\mathrm{\Lambda }$ succeeds in estimating the values around it. On the other hand, for $\ell =2,4$, and 5, where the differences in the evaluation functions are slight, the optimal parameters are relatively close to $(1,1)$. In other words, these cases are problem settings in which the MLEM algorithm can achieve high performance. Even in such cases, our optimizer can estimate the parameters better than it can. The previous study suggested the possibility of a rough parameter estimation by using the noise level. Indeed, the results for $\ell =1,4$, and 5 show similar trends in the contour plots for the same noise level, and the results for $\ell =1,2$ and 3 show a significant change in trend with the noise level. However, the optimal parameter values and contour maps depend not only on the noise level but also on the shape of the reconstruction target.

4.2. Application of Estimated Parameters to Reconstruction of Larger Images

Considering larger images and sinograms in practical CT applications, the significant increase in network size and training cost leads to problems of increased computer memory and computation time. To avoid these problems, we propose the following solution. First, a small reconstructed image and the corresponding projection are prepared as training data. These small-sized data are used to construct the optimizer $\mathrm{\Lambda }$ in a training manner. When estimating parameters from the original projection numbers, the projection is reduced to the same size as the training data and input to the optimizer. Then, the estimated parameters are used to reconstruct the image from the original projection. In this subsection, we verify the feasibility of this procedure by showing that the estimated parameters can be used to reconstruct a larger projection than that used for training.

We generated projections $y^6,y^7,\cdots ,y^{10}$, which were magnified versions of projections $y^1,y^2,\cdots ,y^5$. These projections were generated from phantoms $e^6,e^7,\cdots ,e^{10}$ shown in Figure 12, which were enlargements of $e^1,e^2,\cdots ,e^5$. The projection noise level to be applied was the same as in the previous subsection. The number of projection directions was increased to 180 so that the inverse problem would not be unduly difficult due to the large image size $128\times 128$. In the proposed use case, the parameters are estimated on the basis of a smaller projection reduced from the larger projection, but here the procedure is inverted. That is, the image reconstruction of projection $y^{\ell }$ for any $\ell =6,7,\dots ,10$ uses the parameter ${\lambda }^{\ell -5}$ estimated from its smaller version, projection $y^{\ell -5}$.

Figure 13 plots the values of the evaluation function during reconstruction. The same trend as in Figure 10 is obtained when the parameters estimated from the reduced projection are diverted. The contour maps of the parameters shown in Figure 14 are also consistent with the trend in Figure 11, suggesting that it is possible to use the parameters estimated from the miniaturized projection. Although developing a method to reduce the size of the large original projection appropriately is required for actual use, the results show the possibility of a practical application.

The reconstructed and subtraction images from the true value on the test data $\ell =6$ and $\ell =8$ are shown in Figure 15 and Figure 16. The images reconstructed by the MLEM algorithm and their subtraction images are also shown for comparison. Comparing the reconstructed images for $\ell =6$, it can be seen that the images of the proposed method significantly reduce the overall presence of granular artifacts. In the case of $\ell =8$, the reproducibility of the contour area is high, which is why the evaluation function value is smaller than that of the MLEM algorithm. Thus, the proposed method worked exceptionally well.

5. Conclusions

We developed a system for estimating appropriate parameters for the projection to be reconstructed by employing deep learning on the PDEM algorithm, which has two parameters and is a generalization of various iterative image reconstruction methods. As a learning network, we proposed a structure that connects the optimizer, i.e., the estimation model, and the PDEM algorithm. This architecture is less complex than traditional learned optimization architectures employing RNNs and achieves satisfactory estimations of the parameters of the PDEM algorithm. Although, the appropriate parameters depend on the noise level of the projection and the shape of the phantom, the learned estimation system can estimate them for any projection.

Furthermore, we showed that an optimizer trained on small projections can estimate appropriate parameters for image reconstruction of a projection of a practical size, indicating that the proposed algorithm can be applied to practical problems with larger images and sinograms. This fact enables us to avoid a considerable increase in training costs when reconstructing a clinical CT image.