Image reconstruction based on nonlinear diffusion model for limited-angle computed tomography

For simplicity, we take two-dimensional fan-beam CT scanning mode as an example to illustrate limited-angle CT imaging problem, as shown in figure 1.

Zoom In Zoom Out Reset image size

The line connecting the ray source S and the detector point D passes through the origin O of the coordinate system and is perpendicular to the line AB where the detector is located, satisfying AD = DB. The limited angular range is $\beta\in [\theta, \pi-\theta]$ , where $0\lt \theta \lt \frac{\pi}{2}$.

Using traditional reconstruction algorithms such as simultaneous algebraic reconstruction technique (SART), incomplete data of limited-angle projection will cause artifacts in the reconstruction results. Figure 2 shows the reconstruction results after 5 iterations of full-angle SART reconstruction and limited-angle SART reconstruction for real data, where the limited angle range is $[\frac{\pi}{6}, \frac{5\pi}{6} ]$ . In figure 2(b), the upper and lower ends of the circle (y direction) are blurred and an oblique edge of the lower right triangle is seriously missing.

Zoom In Zoom Out Reset image size

In this paper, we hope to improve traditional reconstruction algorithms with the aid of some diffusion model.

The diffusion model has found widespread application in the field of image processing and has recently been employed in addressing limited-angle CT imaging problems [30, 31]. In AEDS [30] algorithm, besides SART iteration, two regularization terms were incorporated that play the role of edge-preserving diffusion and smoothing along the x-direction and y-direction respectively. However, we find that if the directions are fixed (e.g. along the coordinate axes), some staircase effect would appear in the reconstruction images, see figure 7(c). In order to overcome this problem, we hope to construct a model whose direction is dynamic, variable, and adaptive in some sense.

Different from AEDS, the gradient direction is used as the diffusion direction. It is well known that gradients, which are perpendicular to image edges, are the steepest ascent directions. It is expected that the edge information can be diffused to other non-edge area quickly, so we choose the gradient direction as the direction of information diffusion. Note that the gradients at different image points are generally different, so the directions of information diffusion at different points are also different, which achieves an adaptive diffusion process.

Let u be a 2D image and suppose its gradient $\nabla u = ( \dfrac{\partial u}{\partial x},\dfrac{\partial u}{\partial y})\neq\textbf{0}$. Denote the unit gradient vector as $l = \dfrac{\nabla u}{\parallel\nabla u\parallel_{2}}$. The first and second order directional derivatives of u along direction l are as follows:

$$\begin{equation*} \dfrac{\partial u}{\partial l} = \sqrt{\left( \dfrac{\partial u}{\partial x}\right) ^{2}+\left( \dfrac{\partial u}{\partial y}\right) ^{2}} = |\nabla u|, \end{equation*}$$

$$\begin{equation} \dfrac{\partial^{2}u}{\partial l^{2}} = \dfrac{1}{\left(\frac{\partial u}{\partial x}\right) ^{2}+\left( \frac{\partial u}{\partial y}\right) ^{2}}\left[ \dfrac{\partial^{2}u}{\partial x^{2}}\left( \frac{\partial u}{\partial x}\right) ^{2}+2\dfrac{\partial^{2}u}{\partial x\partial y}\frac{\partial u}{\partial x}\frac{\partial u}{\partial y}+\dfrac{\partial^{2}u}{\partial y^{2}}\left( \frac{\partial u}{\partial y}\right) ^{2}\right]\triangleq \mathcal{G}u. \end{equation} \tag{ 2.1 }$$

Define the following PDE of u:

$$\begin{equation} \lambda\left(x,y\right)\left(u\left(x,y\right)-\phi\left(x,y\right)\right)-\mathcal{G}u = 0, \end{equation} \tag{ 2.2 }$$

where $\phi(x,y)$ is the initial image and $\lambda(x,y)$ plays the role of edge indicators, with larger values for edge points and smaller values for non-edge points. The detailed discussion about $\lambda(x,y)$ will be presented in section 2.2.2.

Remark 2.1. Notice that, we assume $\parallel\nabla u\parallel_{2}\neq0$ before deriving the above model, which means $\dfrac{1}{\parallel\nabla u\parallel_{2}^{2}}$ is reasonable. We set a heterogeneity parameter δ and solve (2.2) when $\parallel\nabla u\parallel_{2}\gt\delta$. Otherwise, the value of u does not change if $\parallel\nabla u\parallel_{2}\unicode{x2A7D} \delta$. We take $\delta = 10^{-25}$ in subsequent experiments. We point out that the above approach is not unique. For instance, one can also replace the term $\dfrac{1}{\parallel\nabla u\parallel_{2}^{2}}$ with an inexact term $\dfrac{1}{\parallel\nabla u\parallel_{2}^{2}+\hat{\delta}}$ with another small parameter $\hat{\delta}$ so that the assumption $\parallel\nabla u\parallel_{2}\gt\delta$ can be removed, and $\mathcal{G}u$ is given by:

$$\begin{equation} \mathcal{G}u = \dfrac{1}{\left(\frac{\partial u}{\partial x}\right) ^{2}+\left( \frac{\partial u}{\partial y}\right) ^{2}+\hat{\delta}} \left[ \dfrac{\partial^{2}u}{\partial x^{2}}\left( \frac{\partial u}{\partial x}\right) ^{2}+2\dfrac{\partial^{2}u}{\partial x\partial y}\frac{\partial u}{\partial x}\frac{\partial u}{\partial y}+\dfrac{\partial^{2}u}{\partial y^{2}}\left( \frac{\partial u}{\partial y}\right) ^{2} \right]. \end{equation} \tag{ 2.3 }$$

Since the unit gradient vector l is a function of u, the aforementioned PDE (2.2) is nonlinear. More specifically, it is a two-dimensional nonlinear steady-state equation. One can solve the model by iteration methods (like Newton's iteration method). Herein, we introduce another approach with the time-stepping technique which seems more convenient for GPU computing. By introducing an artificial variable t, the following time-dependent equation is obtained.

$$\begin{equation} \begin{cases} u_{t}\left(t;x,y\right) & = -\lambda\left(x,y\right)\left(u\left(t;x,y\right)-\phi\left(x,y\right)\right)+\mathcal{G}u,\\ u\left(0;x,y\right) & = \phi\left(x,y\right),\\ u\left(t;x,y\right)|_{\partial\Omega} & = 0, \end{cases} \end{equation} \tag{ 2.4 }$$

where the $\partial\Omega$ is the boundary of the image.

2.2.1. Difference discretization.

By carefully applying central difference formulas, we construct the following difference scheme to approximate the time-dependent equation (2.4):

$$\begin{align} &\dfrac{\partial u}{\partial t}\Big|_{\left(t_n, x_{i}, y_{j}\right)}\approx \dfrac{1}{\tau}\left(u_{i,j}^{n+1}-u_{i,j}^{n}\right), \end{align} \tag{ 2.5 }$$

$$\begin{align} &\dfrac{\partial u}{\partial x}\Big|_{\left(t_n, x_{i}, y_{j}\right)}\approx \dfrac{1}{2h}\left(u_{i+1,j}^n-u_{i-1,j}^n\right), \end{align} \tag{ 2.6 }$$

$$\begin{align} &\dfrac{\partial u}{\partial y}\Big|_{\left(t_n, x_{i}, y_{j}\right)}\approx \dfrac{1}{2s}\left(u_{i,j+1}^n-u_{i,j-1}^n\right), \end{align} \tag{ 2.7 }$$

$$\begin{align} &\dfrac{\partial^{2} u}{\partial x^{2}}\Big|_{\left(t_n, x_{i}, y_{j}\right)}\approx\dfrac{1}{h^{2}}\left(u_{i+1,j}^n-2u_{i,j}^n+u_{i-1,j}^n\right), \end{align} \tag{ 2.8 }$$

$$\begin{align} &\dfrac{\partial^{2} u}{\partial y^{2}}\Big|_{\left(t_n, x_{i}, y_{j}\right)}\approx\dfrac{1}{s^{2}}\left(u_{i,j+1}^n-2u_{i,j}^n+u_{i,j-1}^n\right), \end{align} \tag{ 2.9 }$$

$$\begin{align} &\dfrac{\partial^{2} u}{\partial y\partial x}\Big|_{\left(t_n, x_{i}, y_{j}\right)}\approx\dfrac{1}{2s}\left( \dfrac{u_{i+1,j+1}^n-u_{i-1,j+1}^n}{2h}-\dfrac{u_{i+1,j-1}^n-u_{i-1,j-1}^n}{2h} \right), \end{align} \tag{ 2.10 }$$

where h and s represent the discrete steps in the x and y directions respectively, and τ represents the discrete steps in the time direction. Nodes in the grid are denoted as $(x_{i}, y_{j})$ with $x_{i} = ih, y_{j} = js$, and the corresponding grey values are denoted as $u_{i,j}^n = u(t_n, x_{i}, y_{j})$, $n = 0,1,\ldots,N_t$, $i = 0,1,\ldots,N_1$, $j = 0,1,\ldots,N_2$. Then

$$\begin{align} u_{i,j}^{n+1}& = u_{i,j}^{n}+\tau\left[-\lambda\left(x_i,y_i\right)\left(u_{ij}^n-\phi\left(x_i,y_j\right)\right)\right.\nonumber\\ & \quad+\dfrac{1}{\dfrac{1}{4h^2}\left(u_{i+1,j}^n-u_{i-1,j}^n\right)^2+\dfrac{1}{4s^2}\left(u_{i,j+1}^n-u_{i,j-1}^n\right)^2}\left[\dfrac{\left(u_{i+1,j}^n-2u_{i,j}^n+u_{i-1,j}^n\right)\left(u_{i+1,j}^n-u_{i-1,j}^n\right)^2}{4h^2}\right.\nonumber\\ & \quad+\dfrac{\left(\left(u_{i+1,j+1}^n-u_{i-1,j+1}^n\right)-\left(u_{i+1,j-1}^n-u_{i-1,j-1}^n\right) \right)\left(u_{i+1,j}^n-u_{i-1,j}^n\right)\left(u_{i,j+1}^n-u_{i,j-1}^n\right)}{8h^2s^2}\nonumber\\ & \quad\left.\left.+\dfrac{\left(u_{i,j+1}^n-2u_{i,j}^n+u_{i,j-1}^n\right)\left(u_{i,j+1}^n-u_{i,j-1}^n\right)^2}{4s^2}\right]\right] .\end{align} \tag{ 2.11 }$$

By Taylor's expansion, it is not difficult to obtain the following truncation error results.

Theorem 2.1. The truncation error of the difference scheme (2.5)–(2.10) for the time-dependent equation (2.4) is $O(\tau+h^2+s^2)$ with a constant dependent on δ or $\hat{\delta}$ in remark 2.1 but independent of τ, h and s.

Proof. Denote $Lu = u_{t}(t;x,y)+\lambda(x,y)(u(t;x,y)-\phi(x,y))-\mathcal{G}u$. Without loss of generality, we consider $\mathcal{G}u$ to be the form (2.3) (see remark 2.1). By Taylor's expansion, we have

$$\begin{align} &\dfrac{\partial u}{\partial t}\Big|_{\left(t_n, x_{i}, y_{j}\right)} = \dfrac{1}{\tau}\left(u_{i,j}^{n+1}-u_{i,j}^{n}\right)+O\left(\tau\right)\\ &\dfrac{\partial u}{\partial x}\Big|_{\left(t_n, x_{i}, y_{j}\right)} = \dfrac{1}{2h}\left(u_{i+1,j}^n-u_{i-1,j}^n\right)+O\left(h^2\right)\\ &\dfrac{\partial u}{\partial y}\Big|_{\left(t_n, x_{i}, y_{j}\right)} = \dfrac{1}{2s}\left(u_{i,j+1}^n-u_{i,j-1}^n\right)+O\left(s^2\right)\\ &\dfrac{\partial^{2} u}{\partial x^{2}}\Big|_{\left(t_n, x_{i}, y_{j}\right)} = \dfrac{1}{h^{2}}\left(u_{i+1,j}^n-2u_{i,j}^n+u_{i-1,j}^n\right)+O\left(h^2\right)\\ &\dfrac{\partial^{2} u}{\partial y^{2}}\Big|_{\left(t_n, x_{i}, y_{j}\right)} = \dfrac{1}{s^{2}}\left(u_{i,j+1}^n-2u_{i,j}^n+u_{i,j-1}^n\right)+O\left(s^2\right)\\ &\dfrac{\partial^{2} u}{\partial y\partial x}\Big|_{\left(t_n, x_{i}, y_{j}\right)} = \dfrac{1}{2s}\left( \dfrac{u_{i+1,j+1}^n-u_{i-1,j+1}^n}{2h}-\dfrac{u_{i+1,j-1}^n-u_{i-1,j-1}^n}{2h} \right)+O\left(h^2+s^2\right) . \end{align}$$

Denote that L_h u is the discrete expression of Lu under the difference scheme (2.5)–(2.10). Then we have

$$\begin{align*} L_hu\left(t_n, x_{i}, y_{j}\right)& = \dfrac{1}{\tau}\left(u_{i,j}^{n+1}-u_{i,j}^{n}\right)+\lambda\left(x_i,y_i\right)\left(u_{ij}^n-\phi\left(x_i,y_j\right)\right)\\ & \quad-\dfrac{1}{\dfrac{1}{4h^2}\left(u_{i+1,j}^n-u_{i-1,j}^n\right)^2+\dfrac{1}{4s^2}\left(u_{i,j+1}^n-u_{i,j-1}^n\right)^2+\hat{\delta}}\\ & \quad \times\left[\dfrac{\left(u_{i+1,j}^n-2u_{i,j}^n+u_{i-1,j}^n\right)\left(u_{i+1,j}^n-u_{i-1,j}^n\right)^2}{4h^2}\right.\\ & \quad+\dfrac{\left(\left(u_{i+1,j+1}^n-u_{i-1,j+1}^n\right)-\left(u_{i+1,j-1}^n-u_{i-1,j-1}^n\right) \right)\left(u_{i+1,j}^n-u_{i-1,j}^n\right)\left(u_{i,j+1}^n-u_{i,j-1}^n\right)}{8h^2s^2}\\ & \quad\left.+\dfrac{\left(u_{i,j+1}^n-2u_{i,j}^n+u_{i,j-1}^n\right)\left(u_{i,j+1}^n-u_{i,j-1}^n\right)^2}{4s^2}\right] .\end{align*}$$

Then,

$$\begin{align*} &\left|Lu\left(t_n, x_{i}, y_{j}\right)-L_hu\left(t_n, x_{i}, y_{j}\right)\right|\unicode{x2A7D} O\left(\tau+h^2+s^2\right) .\end{align*}$$

 □

2.2.2. Edge recognition.

In the proposed method, image edges are recognized according to the gradient at each point, where points with larger values are identified as an edge points, while points with smaller values are regarded as non-edge points. For an discrete image u, we set

$$\begin{align*} \begin{aligned} &u_{0,j} = u_{i,0} = 0,\\ &u_{N_1,j} = u_{i,N_2} = 0, \end{aligned} \end{align*}$$

$i = 0,1,\dots,N_1, j = 0,1,\dots,N_2$, where N₁ and N₂ are the pixel numbers of the image u in the x and y directions respectively. The discrete gradient is typically obtained by calculating the pixel difference between adjacent pixels in the four horizontal and vertical directions, i.e. the discrete gradient g_ij of the (i, j)th pixel is usually taken as

$$\begin{equation} g_{ij} = |u_{l}|+|u_{r}|+|u_{u}|+|u_{d}| \end{equation} \tag{ 2.12 }$$

where

$$\begin{equation*} \begin{aligned} u_{l} = u_{i-1,j}-u_{i,j},\\ u_{r} = u_{i+1,j}-u_{i,j},\\ u_{u} = u_{i,j+1}-u_{i,j},\\ u_{d} = u_{i,j-1}-u_{i,j}, \end{aligned} \end{equation*}$$

$i = 1,2,\dots,N_1-1, j = 1,2,\dots,N_2-1$. Introducing a threshold parameter κ, one can identify whether the current pixel is an object edge in the following way.

$$\begin{equation} \begin{cases} g_{ij}\,\,\text{is the edge point},&\text{If }g_{ij}>\kappa,\\ g_{ij}\,\,\text{is the non-edge point},&\text{If }g_{ij}\unicode{x2A7D}\kappa.\\ \end{cases} \end{equation} \tag{ 2.13 }$$

The function $\lambda(x, y)$ serves as an edge indicator, with the values at edge points expected to be large and values at non-edge points expected to be small. We can define $\lambda(x, y)$ as a function related to $\phi(x, y)$, for example, $\lambda(x,y) = c_1 e^{c_2\vert\nabla u\vert^{2}}$ , where $c_1,c_2$ are scalar parameters. For industrial nondestructive testing, the image contrast is usually relatively high and the edges are relatively easy to identify. $\lambda(x, y)$ could be defined as a binary function, with large value at edge points and small value at non-edge points. According to the empirical values of numerical experiments, we set the large value as 10⁶ and the small value as 10⁻²³.

We find that if (2.12) with the difference in four directions is used to identify the edge in solving the proposed model, it may lead to a phenomenon of 'Light Leakage'. For example, let $\Omega = [0,1]\times[0,1],~\Omega_s = \{(x,y)|0.25\unicode{x2A7D} x,y\unicode{x2A7D}0.75\}$, and

$$\begin{equation} \phi\left(x,y\right) = \left\{ \begin{aligned} &2x+2y+\textrm{sin}\left[4\pi\left(x+0.25\right)\right]\textrm{sin}\left[4\pi\left(y+0.25\right)\right],&\left(x,y\right)\in\Omega_s,\\ &0,&\left(x,y\right)\in\Omega\backslash\Omega_s, \end{aligned} \right. \end{equation} \tag{ 2.14 }$$

where a trigonometric function with inter-slice aliasing is taken in the above. Figure 3(a) represents the result by solving the time-dependent equation (2.4) in one time step of using (2.12) to identify edges from the above initial image $\phi(x,y)$ in (2.14), where κ = 0.9 and the size of the image is $512\times512$. Figure 3(b) represents the local amplification results in the lower right corner of figure 3(a). One can find that the inner information of the object in the image leak out of the object.

Zoom In Zoom Out Reset image size

To address the above problem, we propose an improved method for edge recognition. We define the new discrete gradient by using the pixel values in 8 directions of the pixel point shown as

$$\begin{equation} g_{ij} = |u_{l}|+|u_{r}|+|u_{u}|+|u_{d}|+|u_{lu}|+|u_{ld}|+|u_{ru}|+|u_{rd}|, \end{equation} \tag{ 2.15 }$$

where

$$\begin{equation*} \begin{aligned} u_{lu} = u_{i-1,j+1}-u_{i,j},\\ u_{ld} = u_{i-1,j-1}-u_{i,j},\\ u_{ru} = u_{i+1,j+1}-u_{i,j},\\ u_{rd} = u_{i+1,j-1}-u_{i,j}. \end{aligned} \end{equation*}$$

In other words, we use 8 discrete directional derivatives instead of (2.12). We did the same experiment of (2.14) by just replacing (2.12) with the new discrete gradient (2.15). The corresponding result is shown in figure 4 and one can see the 'Light Leakage' phenomenon disappears.

Zoom In Zoom Out Reset image size

This section demonstrates the effectiveness of the proposed model (2.4) with a concrete example. From the initial value $\phi(x,y)$ (2.14), the image is evolved according to the equation (2.4). The settings of parameters are the same as those used in section 2.3. Figure 5 shows the evolution process, where figures 5(a)–(d) are the results of 0, 5, 10 and 50 time steps respectively and figures 5(e)–(h) respectively show the profiles of the red dotted lines shown in figures 5(a)–(d).The red dotted lines in figures 5(f)–(h) respectively show the gray value of the red dotted lines in figures 5(b)–(d) , where each blue solid lines represents the initial signal as a contrast. It is evident that the oscillation caused by trigonometric function in (2.14) in the initial image has been gradually eliminated, and edges have been well preserved. It is also worth noting that the final evolution result (d) approximates a piecewise linear function determined by its edge values, indicating the effectiveness of the nonlinear diffusion model.

Zoom In Zoom Out Reset image size

Based on the previous discussion, we propose the Nonlinear Diffusion Reconstruction (NDR) algorithm for limited-angle CT image reconstruction. The algorithm comprises two main components: solving optimization problems and solving differential equations.

The well used objective functional minimization in CT reconstruction is as follows

$$\begin{equation*}u^* = \text{arg}\min_u\left\{\Vert u^{\left(n\right)}-u\Vert+\boldsymbol{R}\left(u\right)\right\}\end{equation*}$$

where $\Vert u^{(n)}-u\Vert$ is data fidelity term; $\boldsymbol{R}(u)$ represents some traditional regularization term. In this paper, $\boldsymbol{R}(u)$ can be taken as l₀-norm or l_p -norm ($p>$ = 1).

The differential equation is the nonlinear diffusion equation above. The algorithm is implemented as follows

Steps 2–8 are the main ingredient of the algorithm 1, where steps 4–6 are to solve the nonlinear equation. k₁ and k₂ are parameters in the loop of SART, which are usually taken as small integer numbers. k₁ = 5 are taken in the following experiments section 3. k₂ equals 3 in the sections 3.1 and 3.2, and equals 5 in the section 3.3.

Algorithm 1. NDR algorithm.
Input: the initial image $\hat{u}$, the projection p, maximum iteration number parameter $k_\textrm{max}$, integer parameters $k_{1},k_{2}$, and initialize n = 0.
Output: reconstructed image u
1: Execute SART k1 times:
$$\begin{equation*}u^{\left(0\right)} = \text{SART}\left( \hat{u},p,k_{1} \right);\end{equation*}$$
2: while $n\lt k_\textrm{max}$ do
3:   Solving optimization problem:
$$\begin{equation*}u^{\left(n+\frac{1}{2}\right)} = \text{arg}\min_u\left\{\Vert u^{\left(n\right)}-u\Vert+\boldsymbol{R}\left(u\right)\right\}.\end{equation*}$$
4:   Determine the edge information of $u^{(n + \frac{1}{2})}$ according to (2.13) and (2.15);
5:   Determine the edge indicator function according to the edge information;
6:   Solve equation (2.4) by iteration scheme (2.11) taking the initial $\phi(x,y) = u^{(n+\frac{1}{2})}$ to get $u^{(n+1)}$;
7:   $n: = n+1$;
8: end while
9: Execute SART k2 times:
$$\begin{equation*}u = \text{SART}\left( u^{\left(n\right)},p,k_{2} \right).\end{equation*}$$

In this subsection, we establish a reconstruction algorithm based on linear diffusion models. For the reconstruction of two-dimensional CT images, a linear diffusion model is considered, which follows bellow

$$\begin{align} \lambda\left(x,y\right)\left(u\left(x,y\right)-\phi\left(x,y\right)\right)-\left(\frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2}\right) = 0. \end{align} \tag{ 2.16 }$$

The limited-angle CT image reconstruction algorithm based on the linear model is called linear diffusion reconstruction (LDR) algorithm. The pseudocode of two-dimensional LDR is given as follows:

Algorithm 2. LDR algorithm.
Input: The initial image $\hat{u}$, the projection p, maximum iteration number parameter $k_\textrm{max}$, integer parameters $k_{1},k_{2}$, and initialize n = 0.
Output: Reconstructed image u.
1: Execute SART k1 times:
$$\begin{equation*}u^{\left(0\right)} = \text{SART}\left( \hat{u},p,k_{1} \right)\end{equation*}$$
2: while $n\lt k_\textrm{max}$ do
3:   Solving optimization problem:
$$\begin{equation*}u^{\left(n+\frac{1}{2}\right)} = \text{arg}\min_u\left\{\Vert u^{\left(n\right)}-u\Vert+\boldsymbol{R}\left(u\right)\right\}.\end{equation*}$$
4:   Determine the edge information of $u^{(n + \frac{1}{2})}$ according to (2.13) and (2.15);
5:   Determine the edge indicator function according to the edge information;
6:   Solve equation (2.16) by taking the initial $\phi(x,y) = u^{(n+\frac{1}{2})}$ to get $u^{(n+1)}$;
7:   $n: = n+1$;
8: end while
9: Execute SART k2 times:
$$\begin{equation*}u = \text{SART}\left( u^{\left(n\right)},p,k_{2} \right).\end{equation*}$$

Similarly, using the time-stepping technique and introducing the artificial variable t for the elliptic equation (2.16), the corresponding parabolic regularization is obtained for two-dimensional case

$$\begin{equation} \begin{cases} \dfrac{\partial u}{\partial t}& = \left(\dfrac{\partial^2u}{\partial x^2}+\dfrac{\partial^2u}{\partial y^2}\right)-\lambda\left(x,y\right)\left(u\left(x,y\right)-\phi\left(x,y\right)\right),\\ u\left(0;x,y\right)& = \phi\left(x,y\right),\\ u\left(t;x,y\right)|_{\partial\Omega} & = 0. \end{cases} \end{equation} \tag{ 2.17 }$$

We present three numerical difference schemes here.

Format 1: solve the elliptic equation (2.16) directly by the following scheme

$$\begin{align} \lambda_{i,j}\left(u_{i,j}-\phi_{i,j}\right)-\left[\frac{1}{h^2}\left(u_{i+1,j}-2u_{i,j}+u_{i-1,j}\right)+\frac{1}{s^2}\left(u_{i,j+1}-2u_{i,j}+u_{i,j-1}\right)\right] = 0. \end{align} \tag{ 2.18 }$$

Format 2: solve the parabolic equation (2.17) by the following full explicit scheme

$$\begin{align} u_{i,j}^{n+1} = u_{i,j}^{n}+\tau\left[\frac{1}{h^2}\left(u_{i+1,j}^{n}-2u_{i,j}^{n}+u_{i-1,j}^{n}\right)+\frac{1}{s^2}\left(u_{i,j+1}^{n}-2u_{i,j}^{n}+u_{i,j-1}^{n}\right)-\lambda_{i,j}\left(u_{i,j}^{n}-\phi_{i,j}\right)\right]. \end{align} \tag{ 2.19 }$$

Format 3: solve the parabolic equation (2.17) by the following semi-explicit scheme

$$\begin{align} u_{i,j}^{n+1} = u_{i,j}^{n}+\tau\left[\frac{1}{h^2}\left(u_{i+1,j}^{n}-2u_{i,j}^{n}+u_{i-1,j}^{n}\right)+\frac{1}{s^2}\left(u_{i,j+1}^{n}-2u_{i,j}^{n}+u_{i,j-1}^{n}\right)-\lambda_{i,j}\left(u_{i,j}^{n+1}-\phi_{i,j}\right)\right]. \end{align} \tag{ 2.20 }$$

Among the above, $u_{i,j}$ in equation (2.18) denotes the grey value at $(x_{i}, y_{j})$. The definition of h, s, τ, $(x_{i}, y_{j})$ and $u_{i,j}^n$ are the same as in section 2.2.1.

Notice that the approximate solution of the elliptic equation (2.16) can be obtained by means of solving parabolic equation (2.17) with format 2 or format 3. It is not difficult to get the stability condition of format 2 is:

$$\begin{equation*}4r\left(\textrm{sin}^2\frac{\sigma_1\Delta x}{2}+\textrm{sin}^2\frac{\sigma_2\Delta y}{2}\right)+\tau\lambda_{i,j}\unicode{x2A7D} 2,\end{equation*}$$

and the stability condition of format 3 is:

$$\begin{equation*}4r\left(\textrm{sin}^2\frac{\sigma_1\Delta x}{2}+\textrm{sin}^2\frac{\sigma_2\Delta y}{2}\right)-\tau\lambda_{i,j}\unicode{x2A7D} 2,\end{equation*}$$

where $r = \frac{\tau}{h^2}$.

The NDR algorithm and the LDR algorithm can be extended to the three-dimensional case by replacing the two-dimensional nonlinear and linear diffusion equations with the corresponding three-dimensional counterparts. In addition, 26 directions of the pixel will be used to calculate the new discrete gradient for the three-dimensional case.

Experimental CT image reconstruction often employs discrete phantoms. However, simulating projections from discrete phantoms, especially when the phantom and the reconstructed image share the same dimensions, can lead to the inverse crime, resulting in overly optimistic experimental outcomes. In this study, to avoid the occurrence of inverse crime during the reconstruction process, we acquire the projection data analytically, i.e. the forward projection is calculated analytically on an analytic phantom. The phantom is an inclined rectangle with an inclination angle of 15 degrees. The long side is twice as long as short side. Table 1 lists the scanning parameters of the experiment, and figure 6 shows the full-angle SART reconstruction (5 iterations) results.

Zoom In Zoom Out Reset image size

Table 1. Scanning parameters.

Parameter	Values
Scanning type	Fan-beam scanning
Scanning angular interval	1 degree
Number of detector units	1000
Reconstructed image size	$512\times 512$
Scanning angular range	$\left[ \frac{\pi}{6},\frac{5\pi}{6} \right]$

3.1.1. Experiments of the inclined rectangular phantom.

For NDR, the regularization term $\boldsymbol{R}(\cdot)$ can be selected in a variety of ways, and l₀-norm is selected in the experiment of this paper. The parameter κ of NDR is taken as 0.01 and values of other parameters have been introduced in section 2.2. The parameter λ of l₀GM is taken as 10⁻³ while the parameters λ₁ and λ₂ in AEDS are taken as 10⁻⁴and $5\times10^{-5}$. In our experiments with noisy simulation data, the noise we add to the simulation data is Poisson noise specified by the incident intensity which is denoted as I₀, i.e.

$$\begin{equation} p_\textrm{noisy} = -\textrm{log}\left(\textrm{poissrnd}\left(I_0e^{-p}\right)/I_0\right), \end{equation} \tag{ 3.1 }$$

where the symbols $\textrm{log}(\cdot)$ and $\textrm{poissrnd}(\cdot)$ denote two functions for logarithm transform and Poisson random numbers generation, respectively, while p and $p_\textrm{noisy}$ denote the noise-free and noisy projection data. We take $I_0 = 1 \times 10^5$ in our experiments. Experients with noise-free data and noise data are executed on graphic processing cards NVIDIA Tesla P100 and RTX 3090 respectively in this part.

In figure 7, the columns from left to right are the reconstruction results by SART(5 iterations), l₀GM [34], AEDS and NDR algorithms. The first row shows the reconstruction images by noise-free data and the second row shows the reconstruction images by noisy data. Figure 8 shows the corresponding residual images of figrue 7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As shown in figures 7 and 8, the limited-angle SART reconstruction (5 iterations) fail to restore the upper and lower edges. The results of l₀GM algorithm and AEDS algorithm are superior to the results of SART, but staircase effect appears on the upper and lower inclined edges in different extent. It is evident to see the quality of the NDR result is better than others and no staircase effect appears in upper and lower inclined edges, as direction considered in NDR is gradient direction, which is no longer restricted in x axis and y axis. The gray value profiles along column 256 of the reconstruction results are given in figure 9. The image profiles reconstructed by SART has obvious errors across the upper and lower edges of the oblique rectangle. Other three reconstruction results are closer to the reference full-angle reconstruction image. As shown in table 2, we use the peak signal to noise ratio (PSNR) [35] and structural similarity (SSIM) [36] to compare the results quantitatively. The result of indices shows that the NDR algorithm outperforms other algorithms in this experiment.

Zoom In Zoom Out Reset image size

Table 2. PSNR and SSIM of the reconstructed results for the inclined rectangle.

	Measure	SART	l0GM	AEDS	NDR
noise-free	PSNR	28.2094	35.2946	33.6236	37.012
	SSIM	0.9107	0.9816	0.9761	0.9813
noisy	PSNR	27.1777	32.4113	32.348 06	37.0347
	SSIM	0.9035	0.954 66	0.9599	0.979 41

3.1.2. Experiments of the robustness testing.

The experiment aims to evaluate the NDR algorithm's robustness in recovering object edges across various angular orientations. This presents a notably challenging problem due to the alterations in visible and invisible edges with changes in the object's angular positioning. The significant discrepancy between the length and width of the rectangular model amplifies the difficulty in reconstruction, particularly when the longer side serves as the invisible edge. The experiment is conducted on a rectangle model, with the object rotating counterclockwise in 30-degree increments, where the limited angular range for experiments is consistently set at 120 degrees. Experients in this part are executed on a graphic processing card NVIDIA RTX 3090.The reconstruction results are illustrated in figure 10.

Zoom In Zoom Out Reset image size

Observing figure 10, it is apparent that the NDR algorithm effectively reconstructs edges across different angular ranges, seemingly unaffected by the variation in scanning angles. However, combined with table 3, upon integrating residual image analyses and assessment metrics, it becomes evident that as the scanning range rotates, the invisible edges transition between the longer side, the shorter side, and situations where both sides are not entirely invisible edges. The optimal reconstruction occurs when both sides are not completely invisible edges, as evidenced in the experimental results of the central column. Additionally, when the shorter side serves as the invisible edge, the reconstruction outperforms scenarios where the longer side is the invisible edge.

Table 3. PSNR and SSIM of the reconstructed results for robustness testing.

Measure	30 degrees	60 degrees	90 degrees	120 degrees	150 degrees	180 degrees
PSNR	36.7584	49.5221	39.5244	39.4303	46.7386	36.5177
SSIM	0.9839	0.9970	0.9906	0.9904	0.9965	0.9830

The following experiment is conducted to further investigate the stability of the proposed NDR algorithm in this study. The scanning angular range is reduced to 90 degrees, where the long side of a rectangle experiences more severe loss, making the reconstruction more challenging. The proposed NDR algorithm is also compared with the SART, l₀GM and AEDS algorithms, and the results are shown in figure 11. In figure 11, the columns from left to right are the reconstruction results of SART (5 iterations), l₀GM, AEDS and NDR algorithms. The first row displays the reconstruction images, while the second row exhibits the corresponding residual images. It is evident that the NDR algorithm outperforms other algorithms, and as shown in table 4 the values of PSNR and SSIM are also higher.

Zoom In Zoom Out Reset image size

Table 4. PSNR and SSIM of the reconstructed results by 90 limited-angle.

Measure	SART	l0GM	AEDS	NDR
PSNR	23.8528	31.3506	31.4405	34.2548
SSIM	0.8399	0.9660	0.9659	0.9673

Experiments on real projection data, which has noise interference, are performed to further demonstrate the behavior of the proposed NDR algorithm. The experimental phantom is a combination of objects with different shapes, whose boundary types include curved and inclined edges as shown in figure 12(a). Table 5 lists the system and geometric parameters of real data. A full angular CT scanning is carried out from which we can reconstruct image as the reference (figure 12(b)). The purpose of this experiment is to evaluate how well the NDR algorithm can reconstruct the image from the limited-angle projection data with noise.

Zoom In Zoom Out Reset image size

Table 5. System and geometrical scanning parameters.

Parameter	Values
The ray source	YXLON-FXE-225.48
Type of the detector	VarianPS2520V
Scanning type	Fan-beam scanning
Experimental voltage	140kV
Experimental current	100uA
Number of detector units	1920
Detector unit width	0.127 mm
Distance of x-ray source to rotation center	239 mm
Distance of x-ray source to to detector	459 mm
Reconstructed image size	$512\times 512$
Scanning angular range	$\left[ \dfrac{\pi}{6},\dfrac{5\pi}{6} \right]$

3.2.1. Experiments of geometric phantom.

The reconstruction results of the phantom with SART, l₀GM, AEDS and NDR are shown in figure 13, where the first row show their reconstruction results and the second row shows the corresponding residual images. The parameter $\lambda_{l_{0}}$ of l₀GM is taken as 10⁻⁶ while the parameters λ₁ and λ₂ in AEDS are taken as $5\times10^{-7}$and 10⁻⁶. The parameter κ of NDR is taken as 10⁻⁵ and values of other parameters have been introduced in section 2.2. Experients are executed on a graphic processing card NVIDIA Tesla P100 in this part.

Zoom In Zoom Out Reset image size

One can see from the first column of figure 13 that the limited-angle SART reconstruction (5 iterations) has very obvious artifacts. The disk is blurred along the vertical direction and the bottom inclined edge of the triangle at the right is difficult to identify. The l₀ algorithm improves the image quality, but there are deformations in the top of the disk, bottom of the left triangle, and the bottommost oblique edge of the right triangle. The results obtained by AEDS has suppressed the blur and artifacts to a certain extent, but obvious staircase effect is produced. As can be seen from the results of NDR, the blurring and artifacts at the edges are effectively suppressed and there is no staircase effect. The comparison of PSNR and SSIM indicators for above reconstruction results is shown in table 6.

Table 6. PSNR and SSIM of the reconstructed results for real data.

Measure	SART	$l_{0}GM$	AEDS	NDR
PSNR	56.3721	61.4911	61.3668	64.7735
SSIM	0.9927	0.9984	0.9983	0.9988

In order to further observe the process of the NDR algorithm, the image error in the kth cycle of NDR is defined as

$$\begin{equation} err^{\left(k\right)} = \parallel u^{\left(k\right)}-u_{\text{ref}}\parallel_2, \end{equation} \tag{ 3.2 }$$

where u_ref represents the reference image (full-angle reconstruction image). Based on this definition, two residual curves of NDR of the inclined rectanglular phantom experiment with noise-free analytic projection data and the geometric phantom experiment with real data are shown in figure 14 respectively. The overall trend of two residual curves goes to be stable in the iterative process. From our experiments, NDR, to some extent, overcomes the limitations of the existing reconstruction methods, and can significantly reduce the blur and artifacts in limited-angle CT image reconstruction.

Zoom In Zoom Out Reset image size

3.2.2. Experiments of the robustness testing.

The reconstruction difficulty within limited angles primarily stems from the limited data acquired. Consequently, variations in the quantity of acquired data exert a certain influence on the experimental outcomes. Further experimentation was conducted using the data obtained from the previous phase to enhance the analysis. Experients in this part are executed on a graphic processing card NVIDIA RTX 3090. The reconstruction results are illustrated in figure 15. The computed quantitative measures, including PSNR and SSIM indices, are shown in table 7.

Zoom In Zoom Out Reset image size

Table 7. PSNR and SSIM of the reconstructed results for robustness testing.

Measure	90 degree	100 degree	110 degree	120 degree	130 degree	140 degree	150 degree
PSNR	61.4884	63.2017	63.6623	65.1030	65.4009	65.6645	66.7495
SSIM	0.9976	0.9982	0.9986	0.9989	0.9990	0.9990	0.9991

As depicted in (a) through (g), the NDR reconstruction algorithm within the angular range of 90 degrees to 150 degrees can effectively reconstruct the structural form. It is observed that as the angular range expands, there is an improvement in the quality of reconstruction.

3.2.3. Experiments of LDR.

We use LDR algorithm with three numerical formats, discussed in section 2.5, respectively to reconstruct images for this experiment. Experients are executed on a graphic processing card NVIDIA Tesla P100 in this part. The reconstruction image results are shown in figure 16. The PSNR and SSIM measurements corresponding to the reconstruction results are listed in table 8. Based on the definition of image error in (3.2), the error of the results and the reference image is shown in figure 17. Even though the error curve has some fluctuations, the overall trend of three residual curves goes to be stable in the iterative process. Comparing the index values of the results reconstructed by the NDR and LDR algorithms, it is evident that NDR yields better performance.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 8. PSNR and SSIM of the results reconstructed by LDR for real data.

Measure	Format 1	Format 2	Format 3
PSNR	63.8195	64.0942	64.1003
SSIM	0.9986	0.9987	0.9987

The experiment conducted in this section involves the use of the PCB phantom, presenting increased complexity for the proposed method. Notably, the PCB features phantoms with non-piecewise constant intensities and subtle structures. The reconstruction of the PCB phantom from a full-angle scan is illustrated in figure 18. The relevant system and geometric parameters are detailed in table 5. Experients are executed on a graphic processing card NVIDIA RTX 3090 in this subsection.

Zoom In Zoom Out Reset image size

The image obtained through limited-angle SART reconstruction (5 iterations) is presented in figures 19(a) and (d) displays the corresponding residual image. It can be observed that noticeable artifacts exist in many positions in figure 19(a). The reconstruction results of the phantoms using l₀GM and AEDS are illustrated in other subfigures of figure 19, with the first row showcasing their reconstruction outcomes and the second row depicting the corresponding residual images.

Zoom In Zoom Out Reset image size

For l₀ GM, the parameter $\lambda_{l_{0}}$ is set to $5\times10^{-5}$, and for AEDS, the parameters λ₁ and λ₂ are set to $2\times10^{-5}$ and $5\times10^{-5}$, respectively. Both the l₀GM and AEDS algorithms have shown significant improvements in eliminating image artifacts. However, some staircase effects still persist in certain oblique edges.

In comparison, figure 20 presents the reconstruction results of the proposed NDR algorithm in this paper. It is worth noting that this paper utilizes the in-slice magnitude of the image gradient to identify edge points. Unlike the previous experiment where a constant threshold κ was used for the $\lambda(x, y)$ function, a more dynamic strategy is further employed here. A relatively small threshold is initially used and gradually increased for subsequent iterations as the reconstruction images become clearer. In fact, there is no need to set a fixed threshold value for the in-slice gradient norm. Instead, the threshold could be specified as a fixed percentage of pixels with the largest in-slice gradient norm, or it could vary within a certain interval. Figure 20(a) displays the reconstruction result of the NDR algorithm with the threshold $\kappa = 5\times10^{-5}$, while figure 20(b) shows the reconstruction result of the NDR algorithm with the threshold varied from $1\times10^{-8}$ to $1\times10^{-4}$. From the images, it is evident that the NDR algorithm has significantly improved the artifacts near the oblique edges for both constant and varied thresholds. Moreover, the reconstruction results with the varied threshold surpass those with the constant threshold. Similar conclusions can be drawn from table 9. Comparing the PSNR and SSIM values obtained from different algorithms, it can be observed that these metrics align well with the reconstruction results.

Zoom In Zoom Out Reset image size

Table 9. PSNR and SSIM of the reconstructed results for PCB data.

Measure	SART	$l_{0}GM$	AEDS	NDR(constant κ)	NDR(varied κ)
PSNR	48.7460	51.5598	50.1098	54.6685	54.9050
SSIM	0.9845	0.9955	0.9931	0.9970	0.9970