Gamma regularization based reconstruction for low dose CT

The standard CT reconstruction problem can be considered as an inverse problem formulated by equation (1):

$$\begin{eqnarray}&&Gf=y\end{eqnarray} \tag{ 1 }$$

where $f$ denotes the target image of discrete attenuation coefficients, and $y$ represents the calibrated and log-transformed projection data as measurements. $G$ is the system matrix operator reflecting the specific geometry of a given CT imaging system, which can be calculated via techniques like 'pixel-driven' (Peters 1981, Joseph 1982), 'ray-driven' (Siddon 1985, Zhuang et al 1994) or 'distance-driven' (Man and Basu 2004).

As pointed in Lu et al (2001) and Li et al (2004), for low dose CT with low tube currents, the calibrated log-transformed projection data is degraded by additive Gaussian noise, and the relation between the detected projection data $\tilde{\mathop{y}}\,$ and the true projection data $y$ for each channel $i$ is:

$$\begin{eqnarray}&&{{\tilde{\mathop{y}}\,}_{i}}={{y}_{i}}+{{\mu}_{i}}={{(Gf)}_{i}}+\text{Gaussian}\left(0,\sigma _{i}^{2}\right)\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\sigma _{i}^{2}=h\times \text{exp}\left(\frac{{{y}_{i}}}{T}\right)\end{eqnarray} \tag{ 3 }$$

Here,$\text{Gaussian}\left(0,\sigma _{i}^{2}\right)$ denotes the Gaussian noise term ${{\mu}_{i}}$ for each channel $i$ with zero mean and the variance specified by equation (3). $\sigma _{i}^{2}$ is the variance of projection measurement ${{y}_{i}}$ at channel $i$. $T$ and $h$ are object-independent parameters completely determined by the system and can be determined by fitting a large amount of projection data from repeated scanning. Through maximizing the likelihood estimation, we can reconstruct image $f$ by minimizing the weighted least square (WLS) model:

$$\begin{eqnarray}&&f=\underset{f}{\mathop{\text{arg}\,\text{min}}}\,\left\{\Phi(\,f)=\frac{1}{2}W||Gf-\tilde{\mathop{y}}\,|{{|}^{2}}\right\}\end{eqnarray} \tag{ 4 }$$

Where $W$ is a diagonal matrix, with the $i\text{th}$ entry $1/\sigma _{i}^{2}$. Solving $f$ via equation (4) is in fact an ill-posed problem because the rank of the system matrix $G$ is often lower than the rank of $f$, which results in infinite possible solutions (Tikhonov and Arsenin 1997). Some prior knowledge on image $f$ can be used to overcome this ill-posed problem by adding a regularization term $\Psi(\,f)$ into the equation (4), which leads to the following regularization model equation (5):

$$\begin{eqnarray}&&f=\underset{f}{\mathop{\text{argmin}}}\,\left\{\Phi(\,f)=\frac{1}{2}W||Gf-\tilde{\mathop{y}}\,|{{|}^{2}}+\lambda \Psi(\,f)\right\}\end{eqnarray} \tag{ 5 }$$

Here, $\lambda$ is a positive regularization parameter modulating the trade-off between the fidelity term and the regularization term. The prior knowledge can be introduced in $\Psi(\,f)$ using different functions of gradients. We can for instance define the anisotropic regularization ${{\Psi}_{a}}(\,f)$ and isotropic regularization ${{\Psi}_{i}}(\,f)$ (Beck and Teboulle 2009), (Guo and Yin 2012):

$$\begin{eqnarray}&&\Psi(\,f)\,\text{=}\,{{\Psi}_{a}}(\,f)\,=\sum_{i=1}^{m}\sum_{j=1}^{n}\left(\psi \left(|\nabla f_{i,j}^{v}|\right)\,+\,\psi \left(|\nabla f_{i,j}^{h}|\right)\right)\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\Psi(\,f)\,\text{=}\,{{\Psi}_{i}}(\,f)\,\text{=}\sum_{i=1}^{m}\sum_{j=1}^{n}\psi \left(\sqrt{{{\left(\nabla f_{i,j}^{v}\right)}^{2}}+{{\left(\nabla f_{i,j}^{h}\right)}^{2}}}\right)\end{eqnarray} \tag{ 7 }$$

In equations (6) and (7), $\psi$ is the regularization function acting on the gradient term, and ${{f}_{i,j}}$ represents each 2D pixel in the reconstructed 2D image with index $i=1,\ldots m$ and $j=1,\ldots n$. Here, $\nabla f_{i,j}^{{}}$ takes the form $\left(\nabla f_{i,j}^{v},\nabla f_{i,j}^{h}\right)$ including the vertical and horizontal gradients $\nabla f_{i,j}^{v}$ and $\nabla f_{i,j}^{h}$:

$$\begin{eqnarray} &&\nabla f_{i,j}^v = \left\{ {\begin{array}{*{20}{l}} {{f_{i + 1,j}} - {f_{i,j}}i < m}\\ {0i = m} \end{array}, \nabla f_{i,j}^h = } \right.\left\{ {\begin{array}{*{20}{l}} {{f_{i,j + 1}} - {f_{i,j}}j < n}\\ {0j = n} \end{array}, i = 1, \ldots m, j = 1, \ldots ,n} \right. \end{eqnarray} \tag{ 8 }$$

Different regularization function $\psi$ can be used in building $\Psi(\,f)$. Using $\psi (x)={{x}^{2}}$ in the models expressed by equations (6) or (7), we obtain the l₂-norm regularization based WLS reconstruction model (l₂-WLS) (Gonzalez et al 2004):

$$\begin{eqnarray}&&f=\text{arg}\underset{f}{\mathop{\text{min}}}\,\left\{\Phi(\,f)=\frac{1}{2}W||Gf-\tilde{\mathop{y}}\,||_{2}^{2}+\lambda \sum_{i=1}^{m}\sum_{j=1}^{n}\left(|\nabla f_{i,j}^{v}{{|}^{2}}+|\nabla f_{i,j}^{h}{{|}^{2}}\right)\right\}\end{eqnarray} \tag{ 9 }$$

If setting $\psi (x)=x$ in equation (6) for the isotropic model, we obtain the l₁-norm anisotropic regularization WLS reconstruction model (l_1a-WLS) (Beck and Teboulle 2009, Guo and Yin 2012):

$$\begin{eqnarray}&&f=\text{arg}\underset{f}{\mathop{\text{min}}}\,\left\{\Phi(\,f)=\frac{1}{2}W||Gf-\tilde{\mathop{y}}\,||_{2}^{2}+\lambda \left(\sum_{i=1}^{m}\sum_{j=1}^{n}|\nabla f_{i,j}^{h}|+\sum_{i=1}^{m}\sum_{i=1}^{n}|\nabla f_{i,j}^{v}|\right)\right\}\end{eqnarray} \tag{ 10 }$$

Similarly, with $\psi (x)=x$ in the isotropic model in equation (7), the l₁-norm isotropic regularization WLS reconstruction model (l_1i-WLS) (Rudin et al 1992) can be derived:

$$\begin{eqnarray}&&f=\text{arg}\underset{f}{\mathop{\text{min}}}\,\left\{\Phi(\,f)=\frac{1}{2}W||Gf-\tilde{\mathop{y}}\,||_{2}^{2}+\lambda \sum_{i=1}^{m}\sum_{j=1}^{n}\sqrt{|\nabla f_{i,j}^{h}{{|}^{2}}+|\nabla f_{i,j}^{v}{{|}^{2}}}\right\}\end{eqnarray} \tag{ 11 }$$

Through calculus transform, we can easily rewrite the WLS reconstruction model with the l₂-norm regularization equation (9) into equation (12):

$$\begin{eqnarray}&&f=\text{arg}\underset{f}{\mathop{\text{min}}}\,\left\{\Phi(\,f)=\frac{1}{2}W||Gf-\tilde{\mathop{y}}\,||_{2}^{2}+\lambda \sum_{i=1}^{m}\sum_{j=1}^{n}\left(\int_{0}^{|\nabla f_{i,j}^{v}|}2t\text{d}t+\int_{0}^{|\nabla f_{i,j}^{h}|}2t\text{d}t\right)\right\}\end{eqnarray} \tag{ 12 }$$

Similarly, equations (10) and (11) can be respectively rewritten by equations (13) and (14):

$$\begin{eqnarray}&&f=\text{arg}\underset{f}{\mathop{\text{min}}}\,\left\{\Phi(\,f)=\frac{1}{2}W||Gf-\tilde{\mathop{y}}\,||_{2}^{2}+\lambda \left(\sum_{i=1}^{m}\sum_{j=1}^{n}\int_{0}^{|\nabla f_{i,j}^{h}|}\text{d}t+\sum_{i=1}^{m}\sum_{i=1}^{n}\int_{0}^{|\nabla f_{i,j}^{v}|}\text{d}t\right)\right\}\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&f=\text{arg}\underset{f}{\mathop{\text{min}}}\,\left\{\Phi(\,f)=\frac{1}{2}W||Gf-\tilde{\mathop{y}}\,||_{2}^{2}+\lambda \sum_{i=1}^{m}\sum_{j=1}^{n}\int_{0}^{\sqrt{|\nabla f_{i,j}^{h}{{|}^{2}}+|\nabla f_{i,j}^{v}{{|}^{2}}}}\text{d}t\right\}\end{eqnarray} \tag{ 14 }$$

Based on equations (12)–(14), a generalized reconstruction model equation (15) can be obtained by introducing a integrand function $\phi (t)$:

$$\begin{eqnarray}&&f=\underset{f}{\mathop{\text{argmin}}}\,\left\{\Phi(\,f)=\frac{1}{2}W||Gf-\tilde{\mathop{y}}\,||_{2}^{2}+\lambda \sum_{i=1}^{m}\sum_{j=1}^{n}\left(\int_{0}^{|\nabla f_{i,j}^{v}|}\phi (t)\text{d}t+\int_{0}^{|\nabla f_{i,j}^{h}|}\phi (t)\text{d}t\right)\right\}\end{eqnarray} \tag{ 15 }$$

where $\phi (t)$ denotes the general kernel function with positive constraint ($\phi (t)\geqslant 0$). Note equation (15) will become equations (9) or (10) when $\phi (t)\text{=2}t\,or~1$.

We propose a new kernel function ${{\phi}_{\Gamma}}\left(t;\alpha ,\beta \right)$ based on the Gamma probability density function (PDF) in the forms of equations (16) and (17) (Hogg and Craig 1978):

$$\begin{eqnarray}&&{{\phi}_{\Gamma}}\left(t;\alpha ,\beta \right)=\frac{{{t}^{\alpha -1}}{{\beta}^{\alpha}}{{\text{e}}^{-\beta t}}}{\Gamma(\alpha )}\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&\Gamma(\alpha )=\int_{0}^{\infty}\frac{{{t}^{\alpha -1}}}{{{\text{e}}^{t}}}\text{d}t\end{eqnarray} \tag{ 17 }$$

In equation (16), $\alpha$ and $\beta$ denote the shape parameter and the scale parameter, respectively. The Gamma-distributed random variable $t$ in equation (16) has expectation $E$ = $\alpha /\beta$ and variance $V$ = $\alpha /{{\beta}^{2}}$, respectively (Hogg and Craig 1978). Then, the Gamma regularization function ${{\psi}_{\Gamma}}\left(x;\alpha ,\beta \right)$ can be defined as:

$$\begin{eqnarray}&&{{\psi}_{\Gamma}}\left(x;\alpha ,\beta \right)=\int_{0}^{x}\frac{{{t}^{\alpha -1}}{{\beta}^{\alpha}}{{\text{e}}^{-\beta t}}}{\Gamma(\alpha )}\text{d}t\end{eqnarray} \tag{ 18 }$$

Here, ${{\psi}_{\Gamma}}\left(x;\alpha ,\beta \right)$ is in fact the cumulative distribution function (CDF) of the Gamma distributed random variable $t$ with shape parameter $\alpha$ and scale parameter $\beta$ (Hogg and Craig 1978). Figure 1 depicts the regularization functions for l₂-norm regularization, l₁-norm regularization, l₀-norm regularization, and the proposed Gamma regularization function ${{\psi}_{\Gamma}}\left(x;\alpha ,\beta \right)$ with different shape parameter $\alpha \in${1.0, 1.2} and scale parameter $\beta$$\in${2, 4, 6, 8}. We can see in figure 1 that the Gamma regularization functions can be modulated by selecting different values of $\alpha$ and $\beta$ to get a flexible trade-off between l₁-norm regularization and l₀-norm regularization.

Zoom In Zoom Out Reset image size

Consequently, the cost function for the reconstruction with anisotropic Gamma regularization (Γ_a-WLS) can be given by equation (19):

$$\begin{eqnarray}&&f=\underset{f}{\mathop{\text{argmin}}}\,\left\{\Phi(\,f)=\frac{1}{2}W||Gf-\tilde{\mathop{y}}\,||_{2}^{2}+\lambda {{\Psi}_{{{\Gamma}_{a}}}}(\,f)\right\}\end{eqnarray} \tag{ 19 }$$

where

$$\begin{eqnarray}&&{{\Psi}_{{{\Gamma}_{a}}}}(\,f)=\sum_{i=1}^{m}\sum_{j=1}^{n}\int_{0}^{|\nabla f_{i,j}^{v}|}\frac{{{t}^{\alpha -1}}{{\beta}^{\alpha}}{{\text{e}}^{-\beta t}}}{\Gamma(\alpha )}\text{d}t+\sum_{i=1}^{m}\sum_{j=1}^{n}\int_{0}^{|\nabla f_{i,j}^{h}|}\frac{{{t}^{\alpha -1}}{{\beta}^{\alpha}}{{\text{e}}^{-\beta t}}}{\Gamma(\alpha )}\text{d}t\end{eqnarray} \tag{ 20 }$$

Accordingly, the cost function for isotropic Gamma regularization (Γ_i-WLS) is defined as:

$$\begin{eqnarray}&&f=\underset{f}{\mathop{\text{argmin}}}\,\left\{\Phi(\,f)=\frac{1}{2}W||Gf-\tilde{\mathop{y}}\,||_{2}^{2}+\lambda {{\Psi}_{{{\Gamma}_{\text{i}}}}}(\,f)\right\}\end{eqnarray} \tag{ 21 }$$

where

$$\begin{eqnarray}&&{{\Psi}_{{{\Gamma}_{i}}}}(\,f)=\sum_{i=1}^{m}\sum_{j=1}^{n}\int_{0}^{\sqrt{|\nabla f_{i,j}^{h}{{|}^{2}}+|\nabla f_{i,j}^{v}{{|}^{2}}}}\frac{{{t}^{\alpha -1}}{{\beta}^{\alpha}}{{\text{e}}^{-\beta t}}}{\Gamma(\alpha )}\text{d}t\end{eqnarray} \tag{ 22 }$$

The conjugate gradient method (CG) (Boyd and Vandenberghe 2004) was chosen to solve the reconstruction problems in equations (19) and (21). The CG algorithm is an improved steepest descent algorithm, with the descent direction determined by the current descent direction as well as the previous search direction (Fletcher and Reeves 1964). Since the optimization of equation (19) is similar, we only derive the algorithm solving equation (21). We first define the partial derivative ${{\left[\nabla \Phi(\,f)\right]}_{i,j}}$ of the Gamma regularization term in equation (21) with respect to each image point ${{f}_{i,j}}$:

$$\begin{eqnarray}&&{{\left[\nabla \Phi(\,f)\right]}_{i,j}}={{\left({{G}^{T}}W\left(Gf-\tilde{\mathop{y}}\,\right)\right)}_{i,j}}+\lambda \frac{\partial}{\partial {{f}_{i,j}}}\left({{\Psi}_{{{\Gamma}_{i}}}}(\,f)\right)\end{eqnarray} \tag{ 23 }$$

The partial derivative of the regularization term with respect to pixel ${{f}_{i,j}}$ is:

$$\begin{eqnarray}&&\frac{\partial}{\partial {{f}_{i,j}}}\left({{\Psi}_{{{\Gamma}_{i}}}}(\,f)\right)=-\frac{{{\left({{\left(\,{{f}_{i+1,j}}-{{f}_{i,j}}\right)}^{2}}+{{\left(\,{{f}_{i,j+1}}-{{f}_{i,j}}\right)}^{2}}\right)}^{\frac{\alpha}{2}-1}}{{\beta}^{\alpha}}{{\text{e}}^{-\beta {{\left({{\left(\,{{f}_{i+1,j}}-{{f}_{i,j}}\right)}^{2}}+{{\left(\,{{f}_{i,j+1}}-{{f}_{i,j}}\right)}^{2}}\right)}^{\frac{1}{2}}}}}}{\Gamma(\alpha )}\left(\,{{f}_{i+1,j}}-{{f}_{i,j}}\right) \\ &&-\frac{{{\left({{\left(\,{{f}_{i+1,j}}-{{f}_{i,j}}\right)}^{2}}+{{\left(\,{{f}_{i,j+1}}-{{f}_{i,j}}\right)}^{2}}\right)}^{\frac{\alpha}{2}-1}}{{\beta}^{\alpha}}{{\text{e}}^{-\beta {{\left({{\left(\,{{f}_{i+1,j}}-{{f}_{i,j}}\right)}^{2}}+{{\left(\,{{f}_{i,j+1}}-{{f}_{i,j}}\right)}^{2}}\right)}^{\frac{1}{2}}}}}}{\Gamma(\alpha )}\left(\,{{f}_{i,j+1}}-{{f}_{i,j}}\right) \\ &&+\frac{{{\left({{\left(\,{{f}_{i,j}}-{{f}_{i-1,j}}\right)}^{2}}+{{\left(\,{{f}_{i-1,j+1}}-{{f}_{i-1,j}}\right)}^{2}}\right)}^{\frac{\alpha}{2}-1}}{{\beta}^{\alpha}}{{\text{e}}^{-\beta {{\left({{\left(\,{{f}_{i,j}}-{{f}_{i-1,j}}\right)}^{2}}+{{\left(\,{{f}_{i-1,j+1}}-{{f}_{i-1,j}}\right)}^{2}}\right)}^{\frac{1}{2}}}}}}{\Gamma(\alpha )}\left(\,{{f}_{i,j}}-{{f}_{i-1,j}}\right) \\ &&+\frac{{{\left({{\left(\,{{f}_{i+1,j-1}}-{{f}_{i,j-1}}\right)}^{2}}+{{\left(\,{{f}_{i,j}}-{{f}_{i,j-1}}\right)}^{2}}\right)}^{\frac{\alpha}{2}-1}}{{\beta}^{\alpha}}{{\text{e}}^{-\beta {{\left({{\left(\,{{f}_{i+1,j-1}}-{{f}_{i,j-1}}\right)}^{2}}+{{\left(\,{{f}_{i,j}}-{{f}_{i,j-1}}\right)}^{2}}\right)}^{\frac{1}{2}}}}}}{\Gamma(\alpha )}\left(\,{{f}_{i,j}}-{{f}_{i,j-1}}\right)\end{eqnarray} \tag{ 24 }$$

The Gamma regularization term equation (22) can be approximated by equation (25):

$$\begin{eqnarray}&&{{\Psi}_{{{\Gamma}_{i}}}}(\,f)\approx \sum_{i=1}^{m}\sum_{j=1}^{n}\int_{0}^{{{\left(|\nabla f_{i,j}^{h}{{|}^{2}}+|\nabla f_{i,j}^{v}{{|}^{2}}+\varepsilon \right)}^{\frac{1}{2}}}}\frac{{{t}^{\alpha -1}}{{\beta}^{\alpha}}{{\text{e}}^{-\beta t}}}{\Gamma(\alpha )}\text{d}t\end{eqnarray} \tag{ 25 }$$

Here, $\varepsilon$ is a small constant used to avoid the zero denominator in equation (24) (when $\alpha <2$). Then, the corresponding partial derivative of the regularization term equation (25) with respect to pixel ${{f}_{i,j}}$ is:

$$\begin{eqnarray}&&\frac{\partial}{\partial {{f}_{i,j}}}\left({{\Psi}_{{{\Gamma}_{i}}}}(\,f)\right)=-\frac{{{\left({{\left(\,{{f}_{i+1,j}}-{{f}_{i,j}}\right)}^{2}}\,+{{\left(\,{{f}_{i,j+1}}-{{f}_{i,j}}\right)}^{2}}\,+\varepsilon \right)}^{\frac{\alpha}{2}-1}}{{\beta}^{\alpha}}{{\text{e}}^{-\beta {{\left({{\left(\,{{f}_{i+1,j}}-{{f}_{i,j}}\right)}^{2}}+{{\left(\,{{f}_{i,j+1}}-{{f}_{i,j}}\right)}^{2}}+\varepsilon \right)}^{\frac{1}{2}}}}}}{\Gamma(\alpha )}\left(\,{{f}_{i+1,j}}-{{f}_{i,j}}\right) \\ &&-\frac{{{\left({{\left(\,{{f}_{i+1,j}}-{{f}_{i,j}}\right)}^{2}}\,+{{\left(\,{{f}_{i,j+1}}-{{f}_{i,j}}\right)}^{2}}\,+\varepsilon \right)}^{\frac{\alpha}{2}-1}}{{\beta}^{\alpha}}{{\text{e}}^{-\beta {{\left({{\left(\,{{f}_{i+1,j}}-{{f}_{i,j}}\right)}^{2}}+{{\left(\,{{f}_{i,j+1}}-{{f}_{i,j}}\right)}^{2}}+\varepsilon \right)}^{\frac{1}{2}}}}}}{\Gamma(\alpha )}\left(\,{{f}_{i,j+1}}-{{f}_{i,j}}\right) \\ &&+\frac{{{\left({{\left(\,{{f}_{i+1,j-1}}-{{f}_{i,j-1}}\right)}^{2}}\,+{{\left(\,{{f}_{i,j}}-{{f}_{i,j-1}}\right)}^{2}}\,+\varepsilon \right)}^{\frac{\alpha}{2}-1}}{{\beta}^{\alpha}}{{\text{e}}^{-\beta {{\left({{\left(\,{{f}_{i+1,j-1}}-{{f}_{i,j-1}}\right)}^{2}}+{{\left(\,{{f}_{i,j}}-{{f}_{i,j-1}}\right)}^{2}}+\varepsilon \right)}^{\frac{1}{2}}}}}}{\Gamma(\alpha )}\left(\,{{f}_{i,j}}-{{f}_{i,j-1}}\right) \\ &&+\frac{{{\left({{\left(\,{{f}_{i,j}}-{{f}_{i-1,j}}\right)}^{2}}\,+{{\left(\,{{f}_{i-1,j+1}}-{{f}_{i-1,j}}\right)}^{2}}\,+\varepsilon \right)}^{\frac{\alpha}{2}-1}}{{\beta}^{\alpha}}{{\text{e}}^{-\beta {{\left({{\left(\,{{f}_{i,j}}-{{f}_{i-1,j}}\right)}^{2}}+{{\left(\,{{f}_{i-1,j+1}}-{{f}_{i-1,j}}\right)}^{2}}+\varepsilon \right)}^{\frac{1}{2}}}}}}{\Gamma(\alpha )}\left(\,{{f}_{i,j}}-{{f}_{i-1,j}}\right)\end{eqnarray} \tag{ 26 }$$

The overall solution is outlined in algorithm 1 below. We applied the strategy of backtracking line search to improve its convergence. In this strategy, the iteration step $\tau$ is recursively selected to determine the maximum step along a given search direction by means of a backtracking line search with the Armijo-Goldstein condition (Frank and Wolfe 1956, Armijo 1966):

$$\begin{eqnarray}&&\Phi(\,f+\tau d)<\Phi(\,f)+\eta \tau {{\left(\frac{\partial}{\partial f}\Phi(\,f)\right)}^{T}}d\end{eqnarray} \tag{ 27 }$$

where $\eta \in \left[0,\,0.5\right]$, $d$ is the descent direction and the superscript 'T' denotes the matrix transpose operator. The basic idea of the backtracking line search is to recursively decrease the step $\tau$ from a large preset value until a step $\tau$ is found to satisfy the constraint of equation (27). In the algorithm 1, the iteration stopping is jointly controlled by the iteration update limit ${{\varepsilon}_{0}}$ and the maximum iteration number ${{N}_{0}}$.

Algorithm 1. WLS reconstruction with gamma regularization

Input: projection data $\tilde{\mathop{y}}\,$

Output: reconstructed image $f$

Set $\lambda$, $\alpha$, $\beta$, $\eta$, $\rho$, $\tau$, ${{N}_{0}}$, ${{\varepsilon}_{0}}$, $k=1$

Initialize ${{f}_{0}}$, ${{g}_{0}}$, ${{d}_{0}}$, ${{W}_{0}}$

While $||{{f}_{k+1}}-{{f}_{k}}||/{{f}_{k}}>{{\varepsilon}_{0}}$ or $k<{{N}_{0}}$

while $\Phi\left(\,f+\tau {{d}_{k}}\right)>\Phi(\,f)+\eta \tau {{\left(\frac{\partial}{\partial f}\Phi(\,f)\right)}^{T}}{{d}_{k}}$

$\tau =\rho \tau$

end

${{f}_{k+1}}={{f}_{k}}+\tau {{d}_{k}}$

${{y}_{k+1}}=G\times {{f}_{k+1}}$

${{W}_{k+1}}=h\times \text{exp}\left({{y}_{k+1}}/T\right)$

${{g}_{k+1}}=V\Phi\left({{f}_{k+1}}\right)$

$\vartheta =\frac{||{{g}_{k+1}}||_{2}^{2}}{||{{g}_{k}}||_{2}^{2}}$

${{d}_{k+1}}=-{{g}_{k+1}}+\vartheta {{d}_{k}}$

$k=k+1$

end

From the above algorithm outline, we can easily obtain the following relation:

$$\begin{eqnarray}&&\Phi\left(\,{{f}_{k+1}}\right)=\frac{1}{2}{{W}_{k+1}}||G{{f}_{k+1}}-\tilde{\mathop{y}}\,|{{|}^{2}}+\lambda \left({{\psi}_{{{\Gamma}_{i}}}}\left(\sqrt{|\nabla f_{i,j}^{\text{h}}{{|}^{2}}+|\nabla f_{i,j}^{v}{{|}^{2}}}\right)\right)\end{eqnarray}$$

$$\begin{eqnarray}&&\leqslant \frac{1}{2}{{W}_{k}}||G{{f}_{k}}-\tilde{\mathop{y}}\,|{{|}^{2}}+\lambda \left({{\psi}_{{{\Gamma}_{i}}}}\left(\sqrt{|\nabla f_{i,j}^{\text{h}}{{|}^{2}}+|\nabla f_{i,j}^{v}{{|}^{2}}}\right)\right)=\Phi\left({{f}_{k}}\right)\end{eqnarray} \tag{ 28 }$$

Where ${{f}_{k}}$ denotes the ${{k}^{\text{th}}}$ iterated image in reconstruction. This analysis shows that a monotonic decrement of the cost function can be ensured over iterations for the Algorithm 1, which means that a local minimum can be obtained for the proposed reconstruction algorithm.

The regularization function shown in equation (18) monotonically increases with respect to the image gradients and the approximation (to 1) goes stable as gradient increases. The inherent assumption is that the regularization effect should decrease as the gradient values increase because large gradients are often related to image edges, and small gradients often come from noise and artifacts. The shape parameter $\alpha$ and the scale parameter $\beta$ jointly determine the regularization function in equation (18), and are crucial to the performance of the proposed reconstruction model. Figure 2(a) depicts the regularization function values for different $\beta$ when $\alpha$ is fixed to 1.2, and figure 2(b) the values of the regularization function equation (18) for different $\alpha$ in the case of fixed $\beta$. From figures 2(a) and (b), we can see that the regularization function approximates more closely to the l₀ norm as $\beta$ increases or $\alpha$ decreases. This observation means that the parameters $\alpha$ and $\beta$ affect the approximation to l₀ norm in opposite ways. Figures 2(c) and (d) illustrate the regularization functions when $\alpha \in${1.0, 1.1, 1.2, 1.3, 1.4, 1.5} and the ratio $E$ of $\alpha$ and $\beta$ is fixed to 0.15 and 0.2, respectively. We can see in figures 2(c) and (d) that the regularization functions take similar shapes for a fixed $E$ value when $\alpha$ is varying. The plots in figure 2(e) also show that the regularization function becomes closer to a binary function when the parameter $\alpha$ increases to 150 and the $E$ is fixed (the $\beta$ can be deterministically calculated as $\alpha /E$). Considering that a binary regularization function cannot reflect the spatially varying gradients in CT images, the $\alpha$ value is selected within the range [1.0, 1.5] in this study.

Zoom In Zoom Out Reset image size

The plots in figures 2(c) and (d) also show that we can modulate the shapes of regularization functions (or its approximation to 1) by selecting a suitable ratio $E$. The equation (29) can be derived in order to depict the relation between the regularization function and the ratio $E$:

$$\begin{eqnarray}&&\forall \beta ,\,~{{\psi}_{\Gamma}}\left(x=5E;\alpha ,\beta \right)={{\psi}_{\Gamma}}\left(5E;\alpha ,\beta \right)=\int_{0}^{5E}\frac{{{t}^{\alpha -1}}{{\beta}^{\alpha}}{{\text{e}}^{-\beta t}}}{\Gamma(\alpha )}\text{d}t\geqslant 0.99\,~~~~s\text{.t}\text{.}\text{}\alpha \geqslant 1\end{eqnarray} \tag{ 29 }$$

The proof of equation (29) is provided below.

Proof: by simply transforming the regularization function in equation (29), we can get:

$$\begin{eqnarray}&&\forall \beta ,\,~~{{\psi}_{\Gamma}}\left(5E;\alpha ,\beta \right)=\int_{0}^{5E}\frac{{{t}^{\alpha -1}}{{\beta}^{\alpha}}{{e}^{-\beta t}}}{\Gamma(\alpha )}\text{d}t=\int_{0}^{5\beta E}\frac{{{t}^{\alpha -1}}{{\text{e}}^{-t}}}{\Gamma(\alpha )}\text{d}t=\int_{0}^{5\alpha}\frac{{{t}^{\alpha -1}}{{\text{e}}^{-t}}}{\Gamma(\alpha )}\text{d}t\,~~~~s\text{.t}\text{.}\text{}\alpha \geqslant 1\end{eqnarray} \tag{ 30 }$$

Taking the derivative with respect to the variable $\alpha$, we get:

$$\begin{eqnarray}&&\forall \beta ,\,~\frac{\partial}{\partial \alpha}\left({{\psi}_{\Gamma}}\left(5E;\alpha ,\beta \right)\right)=\frac{5{{t}^{\alpha -1}}{{\text{e}}^{-t}}}{\Gamma(\alpha )}>0\,~~s\text{.t}\text{.}\text{}\alpha \geqslant 1\end{eqnarray} \tag{ 31 }$$

Equation (31) implies that the regularization function in equation (30) is a monotonically increasing function with respect to $\alpha$ and we have:

$$\begin{eqnarray}&&\text{}\text{}\forall \beta ,\text{}\text{}\text{}{{\psi}_{\Gamma}}\left(x>5E;\alpha ,\beta \right)>{{\psi}_{\Gamma}}\left(5E;\alpha ,\beta \right)=\int_{0}^{5\alpha}\frac{{{t}^{\alpha -1}}{{\text{e}}^{-t}}}{\Gamma(\alpha )}\text{d}t\geqslant \int_{0}^{5}\frac{{{\text{e}}^{-t}}}{\Gamma(1)}\text{d}t={{\psi}_{\Gamma}}\left(5E;1,\beta \right)\geqslant 0.99\text{}\text{ s}\text{.t}\text{.}\text{}\alpha \geqslant 1\end{eqnarray} \tag{ 32 }$$

Equation (32) shows that the gradients larger than $5E$ ($x>5E$) will cause nearly constant value 1 with derivatives near to zero for the regularization function, and should be considered as image edges to be preserved in the reconstructed image. In the present study, we set the shape parameter $\alpha$ to 1.2, and the value of $5E$ to 25% quantile of the gradient values (Hyndman and Fan 1996). Then the scale parameter $\beta$ can thus be simply calculated as $\alpha /E$. This strategy of parameter setting was found robust in this study. In practical situation, the image to be reconstructed is unavailable, and we just use the FBP reconstructed image to calculate the 25% quantile of the gradient amplitudes.

The two simulated phantoms contain the Modified Shepp-Logan (MSL) phantom (figure 3(a)) and a simulated non-uniform rational B-splines (NURBS) based cardiac-torso (NCAT) phantom (figure 3(b)) (Segars and Tsui 2002). These phantoms are composed by 256   ×   256 pixels with 1 mm  ×  1 mm pixel size. A monochromatic CT parallel scanner with a total of 180 scanning views and 367 radial bins per view is modeled. We simulated two CT dose levels by adding noise1 and noise 2 according to equations (2) and (3) with parameters $T=10000$, $h=5$ and $T=10000$, $h=10$ (Lu et al 2001, Li et al 2004), respectively. For convenience, the corresponding low dose sinograms are named M-sin1(sinogram of MSL phantom with added noise1), M-sin2 (sinogram of MSL phantom with added noise2), N-sin1 (sinogram of NACT phantom with noise1) and N-sin2 (sinogram of NACT phantom with noise2), respectively. All the simulated sinograms are illustrated in figure 4. Quantitative assessments are given in terms of PSNR (peak signal to noise ratio), SNR (signal to noise ratio), and SSIM (structural Similarity Index measuring) (Wang et al 2004). These quantitative metrics are defined by equations (33)–(35):

$$\begin{eqnarray}&&\text{PSNR}(P,I)=10\text{lo}{{\text{g}}_{10}}\left(\frac{P_{\text{max}}^{^{2}}}{\frac{1}{mn}\sum_{i=1}^{m}\sum_{j=1}^{n}{{\left({{P}_{i,j}}-{{I}_{i,j}}\right)}^{2}}}\right)\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}&&\text{SNR}(P,I)=10\text{lo}{{\text{g}}_{10}}\left(\frac{\sum_{i=1}^{m}\sum_{j=1}^{n}{{\left({{P}_{i,j}}\right)}^{2}}}{\sum_{i=1}^{m}\sum_{j=1}^{n}{{\left({{P}_{i,j}}-{{I}_{i,j}}\right)}^{2}}}\right)\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&\text{SSIM}(P,I)=\frac{1}{mn}\frac{\left(2{{\mu}_{p}}{{\mu}_{i}}+{{C}_{1}}\right)\left(2{{\sigma}_{pi}}+{{C}_{2}}\right)}{\left(\,\mu _{p}^{2}+\mu _{i}^{2}+{{C}_{1}}\right)\left(\sigma _{p}^{2}+\sigma _{i}^{2}+{{C}_{2}}\right)}\end{eqnarray} \tag{ 35 }$$

where $I$ is the reconstructed image and $P$, the ground-truth image. ${{P}_{\text{max}}}$ denotes the maximum intensity value in $P$. ${{\mu}_{p}}$ and ${{\mu}_{i}}$ are the mean values of the 8   ×   8 square window in $I$ and $P$. ${{\sigma}_{p}}$ and ${{\sigma}_{i}}$ are the corresponding standard deviations, and ${{\sigma}_{pi}}$ is the corresponding covariance. ${{C}_{1}}$, ${{C}_{2}}$ and ${{C}_{3}}$ are three constant parameters, which are set according to (Wang et al 2004). ${{C}_{1}}={{(0.01\times L)}^{2}}$, ${{C}_{2}}={{(0.03\times L)}^{2}}$, ${{C}_{3}}={{C}_{2}}/2$, with $L$ denoting the grayscale range of the image to reconstruct.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The reconstruction results from the proposed approaches were compared with the results obtained from the l₂-WLS, l_1a-WLS and l_1i-WLS methods. In the experiments conducted on the MSL phantom, with the shape parameter $\alpha$ set to 1.2, the scale parameter $\beta$ was calculated to 6.0 and 4.2 for the Γ_a-WLS and Γ_i-WLS methods, respectively. In the experiments performed on the NACT phantom, with the shape parameter $\alpha$ set to 1.2, the scale parameter $\beta$ was calculated to 16.7 and 11.8 for the Γ_a-WLS and Γ_i-WLS methods, respectively. Figure 5 gathers the reconstructed images for all the simulations carried out with M-sin1, M-sin2, N-sin1 and N-sin2, figures 5(a1) and (d1), (a2) and (d2), (a3) and (d3), (a4) and (d4) and (a5) and (d5) using l₂-WSL, l_1a-WSL, l_1i-WSL Γ_a-WLS and Γ_i-WLS, respectively. We can see that the proposed Γ_a-WLS and Γ_i-WLS methods perform better in preserving edges and suppressing noise than other methods. Tables 2 and 3 list the PSNR, SNR and SSIM values of the reconstructed images with respect to the reference phantom images. The $\lambda$ values are also listed in tables 2 and 3. We can note that the Γ_a-WLS and Γ_i-WLS approaches lead to reconstructions with higher PSNR, SNR and SSIM values. However, the computation time required by the proposed algorithm is higher than for the other methods as shown in tables 2 and 3. Figure 6 depicts the intensity profiles (specified as the red lines in figure 7) for the reconstructions in figure 5. A better match with the ground truth phantom profiles can be observed for the Γ_a-WLS and Γ_i-WLS methods. We can note that, though giving higher PSNR values than the Γ_a-WLS algorithm, the Γ_i-WLS algorithm can reconstruct images with no significant visual difference with respect to the results from the Γ_a-WLS algorithm.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Experiments on the Catphan 600 phantom under low dose protocol were also performed to evaluate the algorithm performance on the preservation of high resolution and low contrast features. The tube voltage and the current were set to 100 KVp and 10 mA. The images to reconstruct are composed by 256   ×   256 pixels with 1 mm  ×  1 mm pixel size. The imaging geometry was specified as: the source to isocenter distance was 500 mm and the distance from the source to the detector was set to 1500 mm. The data were acquired from a monochromatic CT fan beam scanner with a total of 661 scanning views and 1024 radial bins per view. In the experiments worked out using the Catphan 528 phantom, with the shape parameter $\alpha$ set to 1.2, the scale parameter $\beta$ was calculated based on the above parameter setting strategy as 550.46 and 322.58 for Γ_a-WLS and Γ_i-WLS, respectively. For the Catphan 404, $\alpha$ being still set to 1.2, following the same procedure, the scale parameter $\beta$ values were respectively calculated to be 419.58 and 402.68. The global parameter $\lambda$ was set to give the best results in terms of edge preservation and noise suppression.

The reconstructions of the central 2D slices in phantoms Catphan 528 and Catphan 404 are displayed in figures 8 and 9, respectively. In figures 8 and 9, the reconstructed images (a)-(f) correspond to the algorithms FBP (with ramp filter), l₂-WSL, l_1a-WSL, l_1i-WSL, Γ_a-WLS and Γ_i-WLS, respectively. We can see that the images reconstructed by FBP suffer from severe noise, and the l₂-WSL algorithm leads to noise suppression at the cost of obvious detail blurring. The Γ_a-WLS and Γ_i-WLS methods perform better in noise suppression and edge preservation than the others. From figures 8(c)–(f), it is found that the proposed methods provide reconstructions with higher resolution than the l₁-norm regularization methods (see the structures pointed out by red arrows in figure 8). In Figure 9, we can also observe that the proposed methods can preserve low contrast regions with clearer boundaries than the l₁-norm regularization based methods in figures 8(c) and (d) (see the structures defined by red arrows in figure 9).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The convergence property of the Γ_a-WLS and Γ_i-WLS algorithms is analyzed based on the cost functions in equations (19) and (21). M-sin1 was used for this analysis. The evolutions of the cost function values and the PSNR values over iterations are depicted in figure 10. We can observe monotonic decrease of the log-transformed cost function values and the monotonic increase of the PSNR values (of reconstructed images) as the iteration proceeds.

Zoom In Zoom Out Reset image size

Only the shape parameter $\alpha$ needs to be set in the proposed parameter setting strategy. Figure 11 shows the reconstructed images with M-sin1 and N-sin1 using Γ_a-WLS and Γ_i-WLS when the shape parameter $\alpha$ is set to 1.2, 1.8, 2.0 and 5.0, and the parameter $\beta$ calculated as $\alpha /E$ following the above parameter setting strategy. Figure 11 confirms the proposed parameter setting strategy and that the choice of $\alpha$ equal to 1.2 leads to reconstructed images with higher quality than others. In particular, as $\alpha$ is increased more noise from small gradients can be observed in the reconstructed images.

Zoom In Zoom Out Reset image size

In this section, we compared our method with another l₀-approximate regularization with tunable parameter, which was proposed by Hu et al (Hu et al 2011). This regularization also uses a tunable parameter $p$ to realize an approximation to the l₀ norm regularization, with the regularization function given by equation (36):

$$\begin{eqnarray}&&{{\psi}_{\text{Hu}}}(x;p)=\text{log}\left(\frac{x}{p}+1\right)\end{eqnarray} \tag{ 36 }$$

Here, $x$ denotes image gradient. The regularization function ${{\psi}_{\text{Hu}}}$ modulates the function by changing the scale through the parameter $p$. We can see in figure 12(a) that the regularization function ${{\psi}_{\text{Hu}}}$ exhibits a worse approximation to l₀ norm regularization.

Zoom In Zoom Out Reset image size

We also provide the reconstructions of M-sin1 and N-sin1 employing our methods (Γ_a-WLS and Γ_i-WLS) and the regularization function ${{\psi}_{\text{Hu}}}$ with the same weighted least square models used in above experiments. The methods using the regularization function ${{\psi}_{\text{Hu}}}$ are referred as H_a-WLS and H_i-WLS for the anisotropic and isotropic models, respectively. The reconstruction results are illustrated in figure 13 with different combinations of parameters $p$ and $\lambda$ given below. We can see that the proposed Gamma regularization leads to reconstruction results with better performance in both edge preservation and noise suppression than the reconstructions with H_a-WLS and H_i-WLS.

Zoom In Zoom Out Reset image size