Structure-adaptive CBCT reconstruction using weighted total variation and Hessian penalties

Qi Shi; Nanbo Sun; Tao Sun; Jing Wang; Shan Tan

doi:10.1364/BOE.7.003299

1. Introduction

On-board cone-beam computed tomography (CBCT) is extensively used in clinical image-guided radiation therapy [1]. Repeated exposure to CBCT scanning during radiation therapy may increase the risk of cancer [2] and genetic defects [3].

The radiation dose may be reduced by lowering the mAs levels during the CBCT acquisition process. However, this strategy also decreases the number of x-ray photons detected by the scanner. As a result, the conventional FDK method will produce reconstructed CBCT images with increased noise from the projections acquired with lower mAs levels [4].

Several algorithms have been proposed to improve the quality of CBCT imaging from lower mAs level projections [5]. Statistical iterative reconstruction algorithms have shown to perform better than the FDK method. Algorithms based on compressed sensing have also been applied to CBCT reconstruction [6–9]. An adaptive steepest-descent projections on convex sets (ASD-POCS) iterative algorithm with the total variation (TV) penalty has been developed with encouraging reconstruction results [10]. The penalized weighted least-squares (PWLS) iterative algorithm with TV (PWLS-TV) has demonstrated its advantage in suppressing noises and preserving edges, showing state-of-the-art performance in CBCT reconstruction [11, 12].

Although the TV penalty is successful from many points of view, one disadvantage is the so-called staircase effect [13, 14]. Characterized by the first-order differential operator, TV penalizes the image gradient regardless of the image structures, and potentially produces reconstruction results with piecewise-constant artifacts that negatively affect the correct interpretation of images in medical imaging.

Recently, there has been growing interest in using penalties with higher-order differential operators, aiming to suppress the staircase effect for various inverse problems in image processing [15–20]. Higher-order penalties can potentially prevent the oversharpening of regions with smooth intensity transitions because of their piecewise-vanishing property, which produces piecewise-linear solutions that better fit smooth intensity changes. Most high-order penalties involve second-order differential operators. Bredies et al. proposed the Total Generalized Variation (TGV) with symmetric tensors, involving and balancing image higher-order derivatives [15]. TGV is characterized by convexity and weak semi-continuity, and shows appealing properties for natural image denoising particularly in reducing the staircase effect. You and Kaveh proposed a fourth-order partial differential equation (PDE) to optimize the trade-off between noise removal and edge preservation for image denoising problems [16]. Hu et al. developed both isotropic and anisotropic higher degree TV penalties for compressed sensing and image denoising [17]. Lefkimmiatis et al. proposed regularizers involving the Schatten norms of the Hessian matrix computed at every image pixel [18]. Higher-order penalties can effectively suppress the staircase effect, but may also lead to the oversmoothing of object edges [20].

In a previous study [21], we adopted the Hessian penalty to reduce the staircase effect in PWLS iterative CBCT reconstruction. Unlike TV, Hessian is a second-order penalty [19], which allows for linear smooth intensity transition between neighboring pixels. The Hessian penalty has a number of favorable properties including convexity, homogeneity, and rotation and translation invariance. It can also suppress the staircase effect of the TV penalty for PWLS CBCT reconstruction [21]. However, like most higher-order penalties, Hessian tends to slightly blur the edges of the reconstructed image.

A combination of Hessian and TV penalties may overall benefit CBCT reconstruction. Many approaches have been proposed to combine penalties for various image processing applications [22, 23]. Chambolle and Lions proposed the combination of first- and second-order convex regularizations through inf-convolution for image denoising, where the noisy image was decomposed into three parts [24]. Lysaker and Tai proposed the combination of TV with a fourth-order PDE to recover gray-scale images from noisy observations [25]. They separately solved two objective energy functions, one with a first-order penalty (TV) and one with a higher-order penalty, and then used the two solutions to generate a new result through a convex combination. Do et al. proposed a variational approach to combine the first- and higher-order TV penalties by using an equivalent constant weighting strategy for low-dose CT image reconstruction [26]. In their method, the weights were manually set to balance the first- and second-order terms based on visual inspection. Dutaa et al. presented a combination of L¹ and TV penalties for fluorescence molecular tomography; L¹ was used to suppress spurious background signals and enforce sparsity, while TV preserved local smoothness and piecewise constancy in the reconstructed images [27]. They manually selected the best constant weights for both penalties in an objective energy function through a trial-and-error approach. Bioucas-Dias et al. considered the regularizers through linear combinations of “simpler” regularizers [28]. In their method, the weighting parameters for individual regularizers were manually tuned for optimal performance. Chen et al. proposed the combination between the tight framelet and TV to restore images corrupted by blur and impulsive noise [29]. The resulting objective function included two penalties: the TV penalty on the framelet coefficients and the L¹ norm on the restored image. In the same method, the weighting parameters were also manually tuned for different images.

In this study, we proposed a new TV-H penalty combining TV and Hessian in a structure-adaptive way for CBCT reconstruction without introducing extra parameters. The proposed TV-H penalty can automatically adjust the weight parameters between TV and Hessian according to local image structures. We designed an exponential term to automatically recognize edges, gradual transition and uniform regions based on the image gradient. For uniform regions, the exponential term was about 1; for edges, the term approached 0; for gradual transition regions, values varied between 0 and 1. Using the majorization-minimization (MM) approach [30] and the Gauss-Seidel update strategy, we developed an effective algorithm to minimize the objective function with the proposed TV-H penalty. We conducted a systematic comparison between TV, Hessian and TV-H for CBCT reconstruction using two simulated phantoms and two physical phantoms. Our results indicated that the TV-H penalty not only outperformed the TV penalty in suppressing the staircase effect, but also overcame the drawbacks of the Hessian penalty, without blurring object edges.

2. Methods

2.1 PWLS reconstruction

Noise in the CBCT projection data after logarithm transformation follows an approximate Gaussian distribution; noise variance can be determined by the following exponential formula [31]

σ_{i}^{2} = \exp ({\bar{p}}_{i}) / N {}_{i 0}

where

{\bar{p}}_{i}

is the mean of the projection

p_{i}

, and

N_{i 0}

is the incident photon number, both at the detector bin i. Following the noise model in Eq. (1), the cost function of the PWLS algorithm can be written as [32]

Φ (μ) = \frac{1}{2} {(p - A μ)}^{T} Σ^{- 1} (p - A μ) + β R (μ)

The first term in Eq. (2) corresponds to the weighted least-squares (WLS) criterion, and p is the vector of the projection measurements. A represents the projection matrix,

μ

is the attenuation map and

Σ

is a diagonal matrix with its ith element

σ_{i}^{2}

. The symbol T denotes the transpose operator. The second term is a penalty term, where

β

is a parameter that balances the WLS and the penalty terms.

2.2. Adaptive structure reconstruction

The TV penalty consists of the first-order derivative of the reconstructed images. Although the TV penalty in CBCT reconstruction presents advantages in noise suppression and edge preservation [31], it may also lead to the staircase effect. We extended the first-order derivative to a second-order derivative and defined the Hessian penalty in our previous study for CBCT reconstruction [21].

The Hessian penalty successfully suppresses the staircase effect in CBCT reconstruction, but also reduces image resolution and blurs object edges [21]. In contrast, the TV penalty successfully preserves those edges. A combination between the two penalties may eliminate the TV staircase effect while preserving edges. In this study, we proposed the TV-H penalty that combines both penalties as follows:

R_{T V - H} (μ) = \int_{Ω} (1 - α_{X}) | \nabla μ | | d x + \int_{Ω} α_{X} | | H_{μ} (x) | |_{F} d x .

The first term in Eq. (3) is the weighted TV, where

\nabla μ = {(\frac{\partial μ}{\partial x}, \frac{\partial μ}{\partial y}, \frac{\partial μ}{\partial z})}^{T} = {(μ_{x}^{}, μ_{y}^{}, μ_{z}^{})}^{T}

,

Ω \subset ℝ^{3}

is the support of the reconstruction image, and X is the voxel coordinate

(x, y, z)

. The second term in Eq. (3) is a modified Hessian penalty, where

H_{μ} = (\begin{matrix} \frac{\partial^{2} μ}{\partial x^{2}} & \frac{\partial^{2} μ}{\partial x \partial y} & \frac{\partial^{2} μ}{\partial x \partial z} \\ \frac{\partial^{2} μ}{\partial y \partial x} & \frac{\partial^{2} μ}{\partial y^{2}} & \frac{\partial^{2} μ}{\partial y \partial z} \\ \frac{\partial^{2} μ}{\partial z \partial x} & \frac{\partial^{2} μ}{\partial z \partial y} & \frac{\partial^{2} μ}{\partial z^{2}} \end{matrix}) = (\begin{matrix} μ_{x x} & μ_{x y} & μ_{x z} \\ μ_{y x} & μ_{y y} & μ_{y z} \\ μ_{z x} & μ_{z y} & μ_{z z} \end{matrix}),

and

| | \cdot | |_{F}

denotes the Frobenius norm of a matrix.

α_{X}

is an adaptive variable. To take advantage of both penalties,

α_{X}

was designed to automatically balance TV and Hessian according to the local image structure, and defined as

α_{X}^{} = \exp (- \frac{{| \nabla μ_{X}^{} |}_{}^{2}}{η^{2}}),

where

η

determines the strength of each penalty. For edges within an image,

\nabla μ_{X}^{}

is large and leads to a small

α_{X}

, and the TV penalty plays a major role in preserving edges. In contrast, for uniform regions or gradual transition regions, the difference between the two pixels is small. In this situation,

\nabla μ_{X}^{}

becomes small,

α_{X}

increases, and the Hessian penalty plays a major role in suppressing the staircase effect of the TV penalty.

In discrete cases, the TV-H penalty in Eq. (3) can be re-written as

\begin{array}{l} R_{T V - H} = \sum_{x, y, z} {\begin{cases} (1 - α_{(x, y, z)}) \sqrt{μ_{x}^{2} (x, y, z) + μ_{y}^{2} (x, y, z) + μ_{z}^{2} (x, y, z)} \\ + α_{(x, y, z)} \sqrt{μ_{x x}^{2} (x, y, z) + μ_{y y}^{2} (x, y, z) + μ_{z z}^{2} (x, y, z) + 2 μ_{x y}^{2} (x, y, z) + 2 μ_{x z}^{2} (x, y, z) + 2 μ_{y z}^{2} (x, y, z)} \end{cases}} \\ = \sum_{x, y, z} {(1 - α_{(x, y, z)}) \sqrt{μ_{x}^{2} (x, y, z) + μ_{y}^{2} (x, y, z) + μ_{z}^{2} (x, y, z)} + α_{(x, y, z)} \sqrt{\sum_{i = 1}^{6} {[L_{i} (x, y, z) * μ (x, y, z)]}^{2}}}, \end{array}

where

{L_{i} (x, y, z); i = 1, \dots, 6}

denotes the second-order derivative filters; their definitions can be found in Appendix A.

Let $α = {α_{X}}_{\forall X \in Ω}$ . After substituting Eq. (6) into Eq. (2), we obtained the cost function for the structure-adaptive CBCT reconstruction as follows:

\begin{array}{l} Φ (μ, α) = \frac{1}{2} {(p - A μ)}^{T} Σ^{- 1} (p - A μ) + β R_{T V - H} (μ) \\ = \frac{1}{2} {(p - A μ)}^{T} Σ^{- 1} (p - A μ) + β \sum_{x, y, z} {\begin{cases} (1 - α_{(x, y, z)}) \sqrt{μ_{x}^{2} (x, y, z) + μ_{y}^{2} (x, y, z) + μ_{z}^{2} (x, y, z)} \\ + α_{(x, y, z)} \sqrt{\sum_{i = 1}^{6} {[L_{i} (x, y, z) * μ (x, y, z)]}^{2}} \end{cases}} . \end{array}

2.3. Optimization method

Both TV and Hessian penalties are convex. However, the cost function $Φ (μ, α)$ is no longer convex over the two variables $(μ, α)$ . To minimize $Φ (μ, α)$ , we used a strategy similar to the alternating minimization (AM) method. In the standard AM method, the cost function is iteratively minimized with two alternating minimization steps. Each step minimizes the problem over one variable while keeping the second variable fixed. In our strategy, the minimization of $Φ (μ, α)$ was conducted iteratively with a minimization step and an update step. The minimization step minimizes $Φ (μ | α)$ over $μ$ with a fixed $α$ , while the update step updates $α$ using the result from the minimization step. We noted that when $α$ is fixed, $Φ (μ | α)$ is convex over $μ$ .

2.3.1 Minimization step

This step minimizes $Φ (μ | α)$ . We observed that the proposed TV-H penalty in Eq. (6) is not continuously differentiable. Gradient-based optimization methods such as the gradient descent and the conjugate gradient [33] cannot be used for minimizing $Φ (μ | α)$ . In this study, we adopted the MM approach [30]. The MM approach takes a series of quadratic forms of the upper-bound of the original objective function $Φ (μ | α)$ . Solutions that minimize these quadratic forms converge into the solution of the original minimization problem. The minimization of a quadratic form objective function can be easily implemented using methods such as the Gauss-Seidel [5] update strategy.

In the MM approach, $μ^{(t)}$ represents a fixed $μ$ value at the tth iteration and $S (μ)$ is the objective function. $Q (μ | μ^{(t)})$ is a majorizer of $S (μ)$ at $μ^{(t)}$ so that for $\forall μ$ :

Q (μ | μ^{(t)}) \geq S (μ), Q (μ^{(t)} | μ^{(t)}) = S (μ^{(t)}) .

Ordinarily,

μ^{(t)}

represents the current iteration searching for the surface

S (μ)

. In the MM approach, we minimized the majorization

Q (μ | μ^{(t)})

rather than the actual function

S (μ)

.

μ^{(t + 1)}

is the minimizer of

Q (μ | μ^{(t)})

. The MM approach forces

S (μ)

downhill [30]:

S (μ^{(t + 1)}) \leq Q (μ^{(t + 1)} | μ^{(t)}) \leq Q (μ^{(t)} | μ^{(t)}) = S (μ^{(t)}) .

In our study, the following inequality was used to find the quadratic majorizer of the TV-H penalty:

\begin{array}{l} \sqrt{g (x)} \leq \frac{\sqrt{g (y)}}{2} + \frac{g (x)}{2 \sqrt{g (y)}}, \\ \forall y : g (y) \neq 0. \end{array}

If

g (x) = g (y)

, the equality holds. From Eq. (10), we calculated that

\begin{array}{l} \sqrt{μ_{x}^{2} (x, y, z) + μ_{y}^{2} (x, y, z) + μ_{z}^{2} (x, y, z)} \leq {\begin{matrix} \frac{1}{2} \sqrt{{[μ_{x}^{(t)} (x, y, z)]}^{2} + {[μ_{y}^{(t)} (x, y, z)]}^{2} + {[μ_{z}^{(t)} (x, y, z)]}^{2}} \\ + \frac{μ_{x}^{2} (x, y, z) + μ_{y}^{2} (x, y, z) + μ_{z}^{2} (x, y, z)}{2 \sqrt{{[μ_{x}^{(t)} (x, y, z)]}^{2} + {[μ_{y}^{(t)} (x, y, z)]}^{2} + {[μ_{z}^{(t)} (x, y, z)]}^{2}}} \end{matrix}, \end{array}

and

\sqrt{\sum_{i = 1}^{6} {[L_{i} (x, y, z) * μ (x, y, z)]}^{2}} \leq \frac{1}{2} \sqrt{\sum_{i = 1}^{6} {[L_{i} * μ^{(t)} (x, y, z)]}^{2}} + \frac{\sum_{i = 1}^{6} {[L_{i} * μ (x, y, z)]}^{2}}{2 \sqrt{\sum_{i = 1}^{6} {[L_{i} * μ^{(t)} (x, y, z)]}^{2}}} .

Equations (11) and (12) can construct majorizers for TV and Hessian, respectively. The sum of the TV and Hessian majorizers is obviously equal to the TV-H majorizers, satisfying the condition defined in Eq. (10). By replacing Eq. (6) with Eq. (11) and Eq. (12), we have the upper bound of $R_{T V - H}$ :

\begin{array}{l} Q_{R_{T V - H}}^{} (μ | μ^{(t)}) = \sum_{x, y, z} {\begin{matrix} (1 - α_{(x, y, z)}) [\begin{array}{l} \frac{1}{2} \sqrt{{[μ_{x}^{(t)} (x, y, z)]}^{2} + {[μ_{y}^{(t)} (x, y, z)]}^{2} + {[μ_{z}^{(t)} (x, y, z)]}^{2}} \\ + \frac{μ_{x}^{2} (x, y, z) + μ_{y}^{2} (x, y, z) + μ_{z}^{2} (x, y, z)}{2 \sqrt{{[μ_{x}^{(t)} (x, y, z)]}^{2} + {[μ_{y}^{(t)} (x, y, z)]}^{2} + {[μ_{z}^{(t)} (x, y, z)]}^{2}}} \end{array}] \\ + α_{(x, y, z)} [\frac{1}{2} \sqrt{\sum_{i = 1}^{6} {[L_{i} * μ^{(t)} (x, y, z)]}^{2}} + \frac{\sum_{i = 1}^{6} {[L_{i} * μ (x, y, z)]}^{2}}{2 \sqrt{\sum_{i = 1}^{6} {[L_{i} * μ^{(t)} (x, y, z)]}^{2}}}] \end{matrix}}, \\ = \frac{1}{2} μ^{T} (R_{T V}^{(t)} + R_{H e s s i a n}^{(t)}) μ + c o n s t \end{array}

where

R_{T V}^{(t)}

and

R_{H e s s i a n}^{(t)}

are symmetric square matrixes. The majorizers in Eq. (13) can be equivalently expressed as:

\begin{array}{l} Q_{R_{T V - H}} (μ | μ^{(t)}) = \frac{1}{2} \sum_{j} \sum_{i \in N_{j}} [{(w_{j i}^{T V})}^{(t)} + {(w_{j i}^{H e s s i a n})}^{(t)}] {(μ_{j} - μ_{i})}^{2} + c o n s t \\ = \frac{1}{2} \sum_{j} \sum_{i \in N_{j}} w_{j i}^{(t)} {(μ_{j} - μ_{i})}^{2} + c o n s t, \end{array}

where

N_{j}

is the set of neighbor voxels of voxel j. Both

{(w_{j i}^{T V})}^{(t)}

and

{(w_{j i}^{H e s s i a n})}^{(t)}

are the weighting parameters between neighboring voxels determined by

μ^{(t)}

, with

{(w_{j i}^{T V})}^{(t)} = - {(R_{T V}^{(t)})}_{i, j} / 2

and

{(w_{j i}^{H e s s i a n})}^{(t)} = - {(R_{H e s s i a n}^{(t)})}_{i, j} / 2

. In Eq. (14) for a fixed voxel j, the cross terms with index ji arise from the derivative operations, and only exist when i is no more than one voxel and two pixels away from j for TV and Hessian, respectively. More details about

w_{j i}^{(t)}

and its update during the reconstruction process are shown in Appendix B.

By substituting Eq. (14) into Eq. (7), we obtain the following quadratic objective function:

\begin{array}{l} Φ (μ | α) = \frac{1}{2} {(p - A μ)}^{T} Σ^{- 1} (p - A μ) + β (\frac{1}{2} \sum_{j} \sum_{i \in N_{j}} w_{j i}^{(t)} {(μ_{j} - μ_{i})}^{2}) + c o n s t \\ = \frac{1}{2} {(p - A μ)}^{T} Σ^{- 1} (p - A μ) + β (\frac{1}{2} μ^{T} R_{T V - H}^{(t)} μ) + c o n s t . \end{array}

Minimizing the objective function leads to

\frac{\partial Φ (μ | α)}{\partial μ} = 0.

The minimization of the resulting majorizers is equivalent to solving the following series of linear equations

(A^{T} Σ^{- 1} A + β R_{T V - H}^{^{(t)}}) μ = A^{T} Σ^{- 1} p .

Equation (17) can be efficiently solved using the Gauss-Seidel update strategy.

2.3.2 Update step

This step updates $α$ after the minimization step. In this step, each component in $α$ was updated as ${(α)}_{j} = α_{j}^{} = \exp (- {| \nabla μ_{j}^{} |}_{}^{2} / η^{2})$ , according to Eq. (5). Our experiments showed that $η$ can be reasonably set as 40% of the mean of the FDK reconstructed image gradient magnitude. In other words, $η$ was fixed during the iterative reconstruction. More details about our minimization strategy can be found in Appendix C.

2.4. Details in calculating projection matrix A

Calculating the 3D CBCT reconstruction is time-consuming, hence it is important to find an effective method to determine projection matrix $A$ in Eq. (17). To implement the forward projection, we adopted separable footprint (SF) projectors [34]. We specifically used the SF projector with a trapezoid/rectangle function (SF-TR) [34]. SF-TR approximates the voxel footprint functions as 2D separable functions, and uses the trapezoid function in the transaxial direction and the rectangular function in the axial direction.

3. Experiments

We tested our reconstruction algorithm on a Compressed Sensing (CS) phantom [35], a modified 3D Shepp-Logan phantom [36], a commercial calibration phantom CatPhan 600 (The Phantom Laboratory, Inc., Salem, NY) and an anthropomorphic head phantom. The first two are simulated digital phantoms, while the others are physical phantoms.

The projections of the two digital phantoms were simulated using the SF algorithm. There were 360 projections over 360°; the dimension of each projection was 620.8 × 155.2 mm², including 800 × 200 pixels. The source-to-axial distance was 100 cm and the source-to-detector distance was 150 cm. The noise of different levels was added to the projections according to Eq. (1). The incident photon number in Eq. (1) was set to 5 × 10³, 1 × 10⁴ and 5 × 10⁴ to simulate low-dose, medium-dose and high-dose protocols, respectively. The reconstructed image size was 350 × 350 × 16 while the voxel size was 0.776 × 0.776 × 0.776 mm³.

In both physical phantoms, the CBCT projection data were acquired by a Varian Acuity system (Varian Medical System, Palo Alto, CA) with 125 kVp voltage, 10 mA/10 ms current for the low dose and 125 kVp voltage 80 mA/12 ms for the high dose. There were 634 projections in a circle with individual dimensions of 397 × 298 mm², including 1024 × 768 pixels. The projection data were downsampled by a factor of 2 and only 16 out of 768 projection data were chosen for reconstruction along the axial direction. The source-to-axial distance was 100 cm and the source-to-detector distance was 150 cm. For CatPhan 600, the image size was 350 × 350 × 16 and the voxel size was 0.776 × 0.776 × 0.776 mm³. For the anthropomorphic head phantom, the image size was 550 × 550 × 32 mm³ and the voxel size was 0.388 × 0.388 × 0.388 mm³.

3.1. Evaluation criteria

To compare the performance of the TV, Hessian and TV-H penalties we adopted the following criteria: peak signal-to-noise ratio (PSNR), increase in the signal-to-noise ratio (ISNR), structural similarity (SSIM) [37], full width at half maximum (FWHM), edge-spread function (ESF) [5] and contrast-to-noise ratio (CNR).

The PSNR is defined as

P S N R = 10 \log_{10} (\frac{μ_{\max}^{2}}{M S E}),

where

μ_{\max}

is the maximum possible pixel value of the image and

M S E

is the mean-squared error between the reconstructed image and the original image. A higher PSNR indicates a reduced error between the reconstructed image and the original image.

The ISNR is defined as

I S N R = 10 \log (\frac{M S E_{i n}}{M S E_{o u t}}),

where MSE_in is the mean-squared error between the FDK reconstructed image and the original image. MSE_out is the mean-squared error between the reconstructed image and the original image. A higher ISNR indicates better quality of the reconstructed image as compared with the FDK reconstructed image.

The SSIM is defined as

S S I M (a, b) = \frac{(2 μ_{a} μ_{b} + C_{1}) (2 λ_{a b} + C_{2})}{(μ_{a}^{2} + μ_{b}^{2} + C_{1}) (λ_{a}^{2} + λ_{b}^{2} + C_{2})},

where a and b are two windows of 11 × 11 pixels in the same position of two images.

μ_{a}

and

μ_{b}

are the average voxel intensities,

λ_{a}^{2}

and

λ_{b}^{2}

are the voxel intensity variances, all in windows a and b.

λ_{a b}

is the intensity covariance between the two windows.

C_{1}

and

C_{2}

were chosen as

C_{1} = {(0.01 μ_{\max})}^{2}

and

C_{2} = {(0.03 μ_{\max})}^{2}

, respectively. A higher SSIM indicates that the regions of these two images are more similar to each other.

To calculate the FWHM, the gradient was first fitted along a selected profile perpendicular to an object’s edge using the following Gaussian function

f (x) = \frac{1}{ζ \sqrt{2 π}} \exp [- \frac{{(x - x_{0})}^{2}}{2 ζ^{2}}],

where

x

is the voxel location,

ζ

is the standard deviation and

x_{0}

is the mean value of the fitting Gaussian function. The FWHM is then defined as,

F W H M = 2 \sqrt{2 \ln 2} ζ .

A lower FWHM indicates that the object’s edge is sharper.

The resolution of an image can also be described by ESF along profiles through a region of interest (ROI). The ESF is defined as:

y = r + H \cdot erf ((x - \bar{x}) / κ),

where

x

is the voxel location and

y

represents the voxel intensity, erf is the error function, r is the shift of erf and H is the magnitude of erf, and

\bar{x}

is the center of an edge. The parameter

κ

represents the resolution of the image, and a higher

κ

indicates a blurrier edge.

CNR is defined as

C N R = \frac{| μ_{s} - μ_{b} |}{\sqrt{ρ_{s}^{2} + ρ_{b}^{2}}},

where

μ_{s}

and

μ_{b}

are the mean intensity of the ROI and background region, and

ρ_{s}^{}

and

ρ_{b}^{}

are the standard deviation of the ROI and the background region, respectively. A higher CNR suggests better preservation of the signal-to-noise contrast.

These criteria evaluated the performance of different penalties from various standpoints, though some of them presented potential limitations. For example, MSE-based criteria (ISNR and PSNR) do not always accurately model perceptual quality, CNR can become meaningless when applied to the ROI and background with constant intensity because it achieves infinity, while SSIM is not stable enough in regions of low intensity variance [38].

3.2. CS phantom

The CS phantom was used to measure the resolution properties of different penalties and the staircase effect of the TV penalty. This phantom has a uniform background with an attenuation coefficient of 0.0125 mm⁻¹. An octahedron and a ball with linearly gradual transitioning intensity were used to evaluate the potential staircase effect of different penalties. A set of line objects and circular cylinders were used to measure the resolution and CNR of the image, respectively. A representative slice of the CS phantom is shown in Fig. 1. Five green boxes were used to calculate the noise of the phantom, while the red box indicates the ROI.

Fig. 1 One representative slice of the CS phantom

Download Full Size | PDF

3.2.1 Convergence analysis

To examine the convergence of the proposed algorithm, we used the difference in the attenuation coefficient between the successive iterations to indicate whether the algorithm was convergent. Figure (2) shows the convergence curves of the CS phantoms reconstructed by the PWLS algorithm with TV, Hessian and TV-H penalties, where the residual was defined as ${{‖ μ^{(k + 1)} - μ^{(k)} ‖}_{2} / ‖ μ^{(k)} ‖}_{2}$ with k as the iteration number. In this study 30 iterations were used for the reconstruction in PWLS with different penalties. As shown, the three penalties had similar convergence speeds for all noise levels (Fig. 2).

Fig. 2 Convergence curves of the iterative reconstruction for the CS phantom by the three penalties with incident photon numbers: (a) 5 × 10³, (b) 1 × 10⁴, and (c) 5 × 10⁴, respectively.

Download Full Size | PDF

3.2.2 Evaluation

To compare the performance of the three penalties, we adjusted the parameter $β$ in Eq. (2) by trial-and-error so that the reconstruction results had the same noise level. Most evaluation criteria used in this work were related to the noise level. Particularly, image resolution and noise level is always a trade-off. By matching the noise level in the reconstructed image, we can directly compare the relative performance of different penalties [39].

The evaluation index values of the whole image are listed in Table 1. The proposed TV-H and Hessian penalties had similar performances, and were significantly better than the TV penalty and the FDK reconstruction for all incident photon numbers (Table 1). For example, when the incident photon number was 5 × 10³, the PSNR (44.09 dB), ISNR (17.52 dB) and SSIM (0.9812) of the TV-H penalty were higher than those of the TV penalty.

Table 1. A summary of the evalutation indexes of the CS phantom at different dose levels ^a

View Table | View all tables in this article

Figure 3 shows the SSIM maps between the reconstructed images and the original images for the three penalties. The regions near the octahedral and the ball in the reconstructed images using TV-H and Hessian had higher values than those reconstructed using TV, indicating a better preservation of the intensity gradual transition regions. In contrast, the annular and strip-type regions reconstructed using TV-H and TV exhibited higher values than those of the Hessian penalty, suggesting an improved preservation of object edges. In addition, PWLS reconstruction algorithms with the three penalties significantly outperformed the FDK algorithm. As expected, these results indicate that the TV-H penalty successfully enhanced the individual benefits of TV and Hessian.

Fig. 3 SSIM maps of the CS phantom reconstructed by PWLS and the original image at different dose levels. FDK reconstruction with (a) high-, (b) medium- and (c) low-dose projections; PWLS reconstruction with high-dose projection using (d) TV, (g) Hessian and (j) TV-H, respectively; PWLS reconstruction with medium-dose projection using (e) TV, (h) Hessian and (k) TV-H, respectively; PWLS reconstruction with low-dose projection using (f) TV, (i) Hessian and (l) TV-H, respectively.

Download Full Size | PDF

3.3. Modified 3D Shepp-Logan phantom

We modified the original 3D Shepp-Logan phantom by substituting the constant ellipsoid with a linear transitioning ellipsoid, in order to compare the reconstruction performance of different penalties for preserving both edges and transitioning regions in the image.

A representative slice of the modified 3D Shepp-Logan phantom obtained by different reconstruction algorithms is illustrated in Fig. 4, where the original (Fig. 4(a)), the analytical FDK reconstructed image (Fig. 4(b)), the reconstructed image using TV (Fig. 4(c)), the reconstructed image using Hessian (Fig. 4(d)) and the reconstructed image using TV-H (Fig. 4(e)) are shown. For improved illustration, a ROI (indicated by the green square in Fig. 4(a)) was enlarged at the ellipsoid for each image. The display window was set to 0-0.05mm⁻¹ for better illustration. The reconstructed images using TV-H and Hessian successfully preserved the smooth region while the reconstructed image using TV penalty showed an obvious staircase effect.

Fig. 4 A representative slice of the modified 3D Shepp-Logan phantom. (a) Original image; (b) FDK reconstructed image; (c) PWLS reconstructed image using TV; (d) PWLS reconstructed image using Hessian; (e) PWLS reconstructed image using TV-H.

Download Full Size | PDF

Figure 5 reports the profiles along the long red line in Fig. 4 (a) through the gradual transition regions in the modified ellipsoid of the reconstructed images. These plots clearly show that the profiles from TV-H (Fig. 5(e)) and Hessian (Fig. 5(d)) were smooth, while those from TV (Fig. 5(c)) had several obvious steps. This finding is consistent with our visual observation in Fig. 4.

Fig. 5 Profiles along the long red line through the ellipsoid object in Fig. 4. (a) Original image; (b) FDK reconstructed image; (c) PWLS reconstructed image using TV; (d) PWLS reconstructed image using Hessian; (e) PWLS reconstructed image using TV-H.

Download Full Size | PDF

To measure the blurring degree of the object edges, we calculated the derivative of the reconstructed image along the short blue line in Fig. 4(a). We then fitted the derivative with a Gaussian function to obtain the FWHM (Fig. 6). The short blue line goes through a sharp edge. The FWHMs of TV-H, TV and Hessian, were 1.75, 1.61 and 3.16, respectively. These results indicate that both TV-H and TV can preserve the edges better than Hessian.

Fig. 6 Image gradient along the short blue line in Fig. 4. (a) TV reconstruction (FWHM = 1.61); (b) Hessian reconstruction (FWHM = 3.16); (c) TV-H reconstruction (FWHM = 1.75).

Download Full Size | PDF

3.4. CatPhan 600 phantom

A typical slice of the reconstructed results from the CatPhan 600 phantom is reported in Fig. 7. The following information is shown: FDK reconstruction from projection data acquired with the low-dose protocol (10 mA/10 ms) (Fig. 7(a)); FDK reconstruction from projection data acquired with the high-dose protocol (80 mA/10 ms) (Fig. 7(b)); low-dose CBCT image reconstructed using TV (Fig. 7(c)); low-dose CBCT image reconstructed using Hessian (Fig. 7(d)); low-dose CBCT image reconstructed using TV-H (Fig. 7(e)). The three iterative reconstruction results were compared at the same noise level measured with the standard derivation of the reference region, indicated by a red square in Fig. 7(a). Among the three penalties, the noises of the original FDK reconstruction from the low-dose protocol in Fig. 7(a) were greatly reduced.

Fig. 7 A representative slice of the CatPhan 600 phantom: (a) FDK reconstruction with the low-dose (10 mA/10 ms), and (b) high-dose (80 mA/10 ms) protocols; (c) PWLS reconstruction with the low-dose protocol using TV; (d) PWLS reconstruction with the low-dose protocol using Hessian; (e) PWLS reconstruction with the low-dose protocol using TV-H

Download Full Size | PDF

To compare the CNR of the three penalties, we selected four low-contrast ROIs with a size of 17 × 17 pixels. Their exact locations are indicated by arrows in Fig. 7(a). The CNRs of the four regions from the three PWLS algorithms at different noise levels are listed in Table 2. Most of the CNR values from the TV-H reconstruction were similar to those of the TV and Hessian reconstruction, and were higher than those of the FDK (80mA) reconstruction (Table 2). Our findings indicate that these three penalties can preserve image contrast.

Table 2. CNRs of different ROI, as displayed in Fig. 7(a)

View Table | View all tables in this article

3.5. Anthropomorphic head phantom

A representative slice of the results reconstructed from the anthropomorphic head phantom is shown in Fig. 8. Reconstruction and images are shown as follows: FDK reconstruction from the projection data acquired with the low-dose protocol (10 mA/10 ms) (Fig. 8(a)); FDK reconstruction from the projection data acquired with the high-dose protocol (80 mA/10 ms) (Fig. 8(b)); low-dose reconstructions using TV, Hessian and TV-H, respectively (Fig. 8(c)-8(e)). The three iterative reconstruction results were approximately at the same noise level as the standard derivation of the region indicated by a blue square in Fig. 8(a), equal to 4.40 × 10⁻⁴, 4.29 × 10⁻⁴ and 4.49 × 10⁻⁴ for Fig. 8(c), Fig. 8(d) and Fig. 8(e), respectively. We focused on a selected ROI indicated by the red square in Fig. 8(a). We adopted the FDK result from the high-dose projections as a reference. The mean SSIMs of the ROI between the high-dose FDK result and the TV, Hessian and TV-H results were 0.9881, 0.9949 and 0.9935, respectively. The SSIM maps of the ROI for TV, Hessian, and TV-H are illustrated in Fig. 9(a), 9(b) and 9(c), respectively. Both Hessian and TV-H outperformed TV possibly because the head phantom was dominated by gradual intensity transition regions.

Fig. 8 A representative slice of the anthropomorphic head phantom. FDK reconstruction with (a) the low-dose (10 mA/10 ms) and (b) high-dose (80 mA/10 ms) protocols; PWLS reconstruction with the low-dose protocol using (c) TV, (d) Hessian and (e) TV-H penalties.

Download Full Size | PDF

Fig. 9 SSIM maps of the ROI for (a) TV, (b) Hessian, and (c) TV-H result

Download Full Size | PDF

A yellow line was drawn to fit the edge-spread function (ESF) as expressed in Eq. (23) (Fig. 8(a)). The parameter $κ$ was used to characterize the image resolution. The ESFs of the PWLS images reconstructed with the low-dose protocol using TV, Hessian, and TV-H penalty are shown in Fig. 10; $κ$ values were 0.9680, 4.4601 and 1.2143, respectively. This indicates that the resolution of the TV-H reconstruction was similar to that of the TV reconstruction, and higher than that of the Hessian reconstruction.

Fig. 10 Edge-spread function (ESF) along the yellow line in the anthropomorphic head phantom. PWLS reconstruction with the low-dose protocol using (a) TV penalty, (b) Hessian and (c) TV-H.

Download Full Size | PDF

We examined the adaptive parameter $α$ defined in Eq. (5) along the yellow line in the anthropomorphic head phantom for TV-H in the PWLS reconstruction process (Fig. 11). The parameter $α$ balanced the weight between TV and Hessian, and was automatically updated in each iteration during reconstruction. The values of $α$ gradually stabilized, with increases in the iteration number (Fig. 11).

Fig. 11 Update of the parameter $α$ along the yellow line in the anthropomorphic head phantom in the iteration process; (a) 1st iteration, (b) 10th iteration, (c) 20th iteration and (d) 30th iteration.

Download Full Size | PDF

The $α$ of a representative slice of the anthropomorphic head phantom is shown in Fig. 12. The parameter $α$ played the role of an edge detector, and gradually converged into the real boundary during the iteration process.

Fig. 12 Update of the parameter $α$ on a representative slice of the anthropomorphic head phantom in the iteration process; (a) 1st iteration, (b) 10th iteration, (c) 20th iteration, and (d) 30th iteration.

Download Full Size | PDF

4. Discussion

The state-of-the-art TV penalty successfully suppresses noise and preserves edges in low dose CBCT reconstruction. However, one of its drawbacks is the staircase effect, which modifies gradual intensity transition regions into piecewise constant regions in the reconstruction images.

We described a new TV-H penalty that balances the properties of TV and Hessian in the reconstruction process. As shown in the two digital simulation experiments, piecewise constant artifacts were observed over the regions with smooth intensity change in the TV reconstruction results, e.g., in the ellipsoid of the modified 3D Shepp-Logan phantom (Fig. 4(c)). The SSIM values also demonstrated the failure of the TV penalty in recovering the regions near the octahedral and the ball in the CS phantom (Fig. 3(d), 3(e), 3(f))). These TV reconstruction images were degraded by the piecewise constant regions and lost their original structures. The Hessian penalty suppresses the staircase effect by introducing the second-order derivative into the CBCT reconstruction process [21]. However, the Hessian penalty cannot properly preserve the edge of the line-style and the annular objects of the CS phantom (Fig. 3(g), 3(h), 3(i)), and the edges of the Shepp-Logan phantom (Fig. 6(b)). The proposed TV-H penalty can overcome the staircase effect of the TV penalty and also preserve object edges, as shown in Fig. 3(j), 3(k), 3(l)), Fig. 4(e) and Fig. 5(e).

The adaptive parameter $α$ was used to automatically detect edges, gradual transition and uniform regions in the reconstruction image, so that the proposed TV-H penalty could adaptively weigh different penalties in different local regions of the reconstruction images. A similar idea to explore local structure was previously published [39].Wang et al. constructed an anisotropic penalty by examining the intensity difference between a voxel and its direct neighbors. The anisotropic penalty in Wang’s method automatically adjusted the weight of the same penalty (quadratic penalty) over different local image structures.

In our method, $α$ was automatically updated during the iteration. We observed a significant fluctuation in $α$ at the beginning of the reconstruction process, as shown in Fig. 11(a)-11(c) and Fig. 12(a)-12(c). Noise is primarily responsible for this fluctuation. At the first iteration, $α$ was calculated from the FDK result which typically had a high noise level. With increasing numbers of iterations, the noise level in the reconstruction decreased, leading to $α$ with lower fluctuation (Fig. 11(d) and Fig. 12(d)). As the number of iterations increases, $α$ gradually stabilizes. In our experiment, we found that a maximum iteration number of 30 was enough to obtain satisfactory reconstruction.

Most current approaches based on mixed penalties introduce parameters that need extra manual adjustments [25–29]. We only observed one parameter, $η$ , in the TV-H penalty (See Eq. (5)) and found that its variation did not significantly change the reconstruction performance. We set $η$ as 40% of the mean of the FDK reconstructed image gradient magnitude in all experiments. Except for this parameter, the proposed reconstruction algorithm with TV-H did not introduce any other extra parameters. This is a potentially important advantage of our penalty. The additional parameter we needed to manually adjust in our study was the weighting parameter β. This parameter is required in a typical inverse problem, and is used to balance the fidelity and the penalty terms in the whole objective function, as in Eq. (2) or Eq. (7). For fair comparison, we manually adjusted this parameter to keep the images reconstructed with TV, Hessian and TV-H at the same noise level.

Bredies et. al previously extended TV to TGV [15]. TGV can incorporate the first and higher order derivatives, similarly to our proposed TV-H penalty. However, TGV and TV-H differ in many ways. We mentioned the second-order TGV as an example. The second-order TGV of image $μ$ can be written as $T G V_{α}^{2} (μ) = \min_{v \in Ω} α_{0} \int_{Ω} | \nabla μ - v | d x + α_{1} \int_{Ω} | ε (v) | d x,$ where $ε (v) = \frac{1}{2} (\nabla v + \nabla v^{T})$ is the matrix-valued weak symmetrized derivative, and $(α_{0}, α_{1})$ are two positive parameters [40]. First, the TGV does not explicitly incorporate the second derivative since the actual minimum $v$ may be located anywhere between 0 and $\nabla μ$ [40]. In contrast, our proposed TV-H combines the real TV and Hessian in a structure-adaptive way. Second, the first and second order terms in TGV are balanced via weights $(α_{0}, α_{1})$ . TV-H only includes one parameter $η$ (see Eq. (5)). We noted that the reconstruction was quite robust to the variation of $η$ , which was kept unchanged for all experiments in our study.

In this study, the PWLS algorithm (Eq. (2)) was adopted for CBCT reconstruction at a low-dose level, and the noise of the projection data was assumed to follow a Gaussian distribution with zero mean and signal dependent variance. Wang et al. studied the noise properties of the calibrated or preprocessed sinogram data in Radon space, and validated the nonlinear mean-variance relationship (i.e., Eq. (1)) of the projected data at different mAs levels [31]. The mean–variance relationship precludes the use of Poisson and Gamma distribution for image reconstruction. The authors argued that a PWLS cost function could be a suitable choice, fully utilizing the relationship between the first and the second statistical moments [31]. Nevertheless, the Gaussian distribution is only a rough approximation of the projection noise in CBCT [31, 41]. The reconstruction might benefit from a more sophisticated noise model on the raw projection data (i.e., before logarithm transform) such as Poisson-Gaussian [42, 43] and compound Poisson [44]. In addition, algorithms designed for signal-dependent noise may further improve the reconstruction performance [45]. In this study, the noise covariance matrix $Σ$ had a diagonal form, and no correlation was considered among the nearest bins of the CBCT projection. Considering the noise correlation can improve the reconstruction performance [46].

An important problem of statistical iterative reconstruction methods is the long computational time as compared to the FDK method [4]. In our experiments, 2 minutes were needed to complete one iteration to reconstruct CBCT images of 350 × 350 × 16 using a PC with 3.3GHz CPU. Thus, one image with 30 iterations will take at least 1 hour, which can be sped by various methods. As for the converging speed, many effective fast converging algorithms have recently been proposed including the FISTA [47] and the NESTA [6, 48]. In terms of hardware, graphics processing unit (GPU) can be used to reduce the computation time in CBCT reconstruction [12, 49].

5. Conclusion

In this study, we proposed the TV-H penalty for CBCT reconstruction, and developed an effective algorithm to minimize the objective function using the majorization-minimization strategy. The TV-H penalty balanced the weight between TV and Hessian in the objective function in a structure-adaptive way. Comparison studies among TV, Hessian and TV-H revealed that the proposed penalty effectively preserves both edges and other regions with smooth intensity transition, while restraining the staircase effect of the TV penalty.

Appendix A derivative filters

For a three-dimensional image, there are six commonly used second-order derivatives. The voxel index is represented by the triplet $(x, y, z)$ . For a discrete image $μ$ , the filtering operations can be defined as

\begin{array}{l} L_{1} (x, y, z) * μ (x, y, z) = μ (x + 1, y, z) - 2 μ (x, y, z) + μ (x - 1, y, z); \\ L_{2} (x, y, z) * μ (x, y, z) = μ (x, y + 1, z) - 2 μ (x, y, z) + μ (x, y - 1, z); \\ L_{3} (x, y, z) * μ (x, y, z) = μ (x, y, z + 1) - 2 μ (x, y, z) + μ (x, y, z - 1); \\ L_{4} (x, y, z) * μ (x, y, z) = \sqrt{2} [μ (x, y, z) - μ (x - 1, y, z) - μ (x, y - 1, z) + μ (x - 1, y - 1, z)]; \\ L_{5} (x, y, z) * μ (x, y, z) = \sqrt{2} [μ (x, y, z) - μ (x - 1, y, z) - μ (x, y, z - 1) + μ (x - 1, y, z - 1)]; \\ L_{6} (x, y, z) * μ (x, y, z) = \sqrt{2} [μ (x, y, z) - μ (x, y - 1, z) - μ (x, y, z - 1) + μ (x, y - 1, z - 1)] . \end{array}

Appendix B $w_{j i}^{(t)}$ and its update

For voxel $j$ , $N_{j}$ is the set of neighbor voxels which are at most 2-voxel away from it. From Eqs. (13) and (14), we can easily show that $w_{j i}^{(t)} : = - {(R_{T V - H}^{(t)})}_{j i} / 2, i \in N_{j}, i \neq j$ . $w_{j i}^{(t)}$ is updated in each iteration. For the proposed TV-H penalty, we defined $W_{(x, y, z)}^{(t)}$ as

W_{(x, y, z)}^{(t)} : = α_{(x, y, z)} / \sqrt{\sum_{i = 1}^{6} {[L_{i} (x, y, z) * μ^{(t)} (x, y, z)]}^{2}} .

Let

s

be a small constant number to avoid singularity of the derivative of TV. For each neighbor

i \in N_{j}

, the update formula of

w_{j i}^{(t)}

can be expressed as:

\begin{array}{l} if i = (x + 1, y, z), w_{j i}^{(t)} = {\begin{cases} \frac{1 - α_{(x + 1, y, z)}^{}}{\sqrt{{(μ_{x + 1, y, z} - μ_{x, y, z})}^{2} + {(μ_{x + 1, y, z} - μ_{x + 1, y - 1, z})}^{2} + {(μ_{x + 1, y, z} - μ_{x + 1, y, z - 1})}^{2} + s^{2}}} + \\ (W_{(x, y, z)}^{(t)} + 3 W_{(x + 1, y, z)}^{(t)} + W_{(x + 1, y + 1, z)}^{(t)} + W_{(x + 1, y, z + 1)}^{(t)}) \end{cases}}; \\ if i = (x, y + 1, z), w_{j i}^{(t)} = {\begin{cases} \frac{1 - α_{(x, y + 1, z)}^{}}{\sqrt{{(μ_{x, y + 1, z} - μ_{x - 1, y + 1, z})}^{2} + {(μ_{x, y + 1, z} - μ_{x, y, z})}^{2} + {(μ_{x, y + 1, z} - μ_{x, y + 1, z - 1})}^{2} + s^{2}}} + \\ (W_{(x, y, z)}^{(t)} + 3 W_{(x, y + 1, z)}^{(t)} + W_{(x + 1, y + 1, z)}^{(t)} + W_{(x, y +1, z +1)}^{(t)}) \end{cases}}; \\ if i = (x, y, z + 1), w_{j i}^{(t)} = {\begin{cases} \frac{1 - α_{(x, y, z + 1)}^{}}{\sqrt{{(μ_{x, y, z + 1} - μ_{x - 1, y, z + 1})}^{2} + {(μ_{x, y + 1, z + 1} - μ_{x, y, z + 1})}^{2} + {(μ_{x, y, z + 1} - μ_{x, y, z})}^{2} + s^{2}}} + \\ (W_{(x, y, z)}^{(t)} + 3 W_{(x, y, z +1)}^{(t)} + W_{(x +1, y, z + 1)}^{(t)} + W_{(x, y + 1, z + 1)}^{(t)}) \end{cases}}; \\ if i = (x - 1, y, z), w_{j i}^{(t)} = {\begin{cases} \frac{1 - α_{(x, y, z)}^{}}{\sqrt{{(μ_{x, y, z} - μ_{x - 1, y, z})}^{2} + {(μ_{x, y, z} - μ_{x, y - 1, z})}^{2} + {(μ_{x, y, z} - μ_{x, y, z - 1})}^{2} + s^{2}}} + \\ (W_{(x - 1, y, z)}^{(t)} + 3 W_{(x, y, z)}^{(t)} + W_{(x, y + 1, z)}^{(t)} + W_{(x, y, z + 1)}^{(t)}) \end{cases}}; \\ if i = (x, y - 1, z), w_{j i}^{(t)} = {\begin{cases} \frac{1 - α_{(x, y, z)}^{}}{\sqrt{{(μ_{x, y, z} - μ_{x - 1, y, z})}^{2} + {(μ_{x, y, z} - μ_{x, y - 1, z})}^{2} + {(μ_{x, y, z} - μ_{x, y, z - 1})}^{2} + s^{2}}} + \\ (W_{(x, y - 1, z)}^{(t)} + 3 W_{(x, y, z)}^{(t)} + W_{(x, y, z + 1)}^{(t)} + W_{(x + 1, y, z)}^{(t)}) \end{cases}}; \\ if i = (x, y, z - 1), w_{j i}^{(t)} = {\begin{cases} \frac{1 - α_{(x, y, z)}^{}}{\sqrt{{(μ_{x, y, z} - μ_{x - 1, y, z})}^{2} + {(μ_{x, y, z} - μ_{x, y - 1, z})}^{2} + {(μ_{x, y, z} - μ_{x, y, z - 1})}^{2} + s^{2}}} + \\ (W_{(x, y, z - 1)}^{(t)} + 3 W_{(x, y, z)}^{(t)} + W_{(x + 1, y, z)}^{(t)} + W_{(x, y + 1, z)}^{(t)}) \end{cases}}; \\ if i = (x, y -1, z -1) or (x - 1, y, z -1) or (x -1, y -1, z), w_{j i}^{(t)} = - W_{(x, y, z)}^{(t)}; \\ if i = (x, y +1, z +1), w_{j i}^{(t)} = - W_{(x, y + 1, z + 1)}^{(t)}; if i = (x, y +1, z -1), w_{j i}^{(t)} = - W_{(x, y + 1, z)}^{(t)}; if i = (x, y -1, z +1), w_{j i}^{(t)} = - W_{(x, y, z + 1)}^{(t)}; \\ if i = (x +1, y, z +1), w_{j i}^{(t)} = - W_{(x + 1, y, z + 1)}^{(t)}; if i = (x + 1, y, z -1), w_{j i}^{(t)} = - W_{(x + 1, y, z)}^{(t)}; if i = (x - 1, y, z +1), w_{j i}^{(t)} = - W_{(x, y, z + 1)}^{(t)}; \\ if i = (x - 1, y +1, z), w_{j i}^{(t)} = - W_{(x, y + 1, z)}^{(t)}; if i = (x +1, y +1, z), w_{j i}^{(t)} = - W_{(x + 1, y + 1, z)}^{(t)}; if i = (x +1, y -1, z), w_{j i}^{(t)} = - W_{(x + 1, y, z)}^{(t)}; \\ if i = (x +2, y, z), w_{j i}^{(t)} = - W_{(x + 1, y, z)}^{(t)} /2; if i = (x, y +2, z), w_{j i}^{(t)} = - W_{(x, y + 1, z)}^{(t)} /2; if i = (x, y, z +2), w_{j i}^{(t)} = - W_{(x, y, z + 1)}^{(t)} /2; \\ if i = (x -2, y, z), w_{j i}^{(t)} = - W_{(x - 1, y, z)}^{(t)} /2; if i = (x, y -2, z), w_{j i}^{(t)} = - W_{(x, y - 1, z)}^{(t)} /2; if i = (x, y, z -2), w_{j i}^{(t)} = - W_{(x, y, z - 1)}^{(t)} /2 . \end{array}

Appendix C structure adaptive CBCT reconstruction using TV-H penalty

Our task was to estimate the attenuation coefficient distribution map $μ$ from the projection data $p$ by minimizing the cost function $Φ (μ, α)$ defined in Eq. (7). To minimize $Φ (μ, α)$ , we used a strategy similar to the alternating minimization algorithm. We noted that $Φ (μ, α)$ is non-convex, and as a result the reconstruction depends on initialization. In our experiments, $μ$ was initialized using the FDK reconstruction $F D K {p}$ and $α$ was updated using Eq. (5). Let K and M be the maximum iteration number for the update and MM approach, respectively, N be the image size (the voxel numbers that need to be updated in the Gauss-Seidel method), and $T_{1}$ and $T_{2}$ be two small real number.

Initialization

\begin{array}{l} {\tilde{μ}}^{(1)} = F D K {p} \\ r = p - A {\tilde{μ}}^{(1)} \\ s_{j} = A_{j}^{T} Σ^{- 1} A_{j}, \forall j \end{array}

Iteration

For $k = 1 : K$ (Update step)

\begin{array}{l} {(α^{(k)})}_{j} : = α_{j}^{(k)} = \exp (- {| \nabla {\tilde{μ}}_{j}^{(k)} |}^{2} / η^{2}), \forall j (S e e E q . (5)) \\ μ : = {\tilde{μ}}^{(k)}, μ^{1} : = {\tilde{μ}}^{(k)}, α : = α^{(k)} \end{array}

For $t : = 1 : M$ (Minimization step)

For j=1: N (Gauss-Seidel)

$u p d a t e {w_{j i}^{(t)}} (S e e A p p e n d i x B)$

\begin{array}{l} μ_{j}^{o l d} : = μ_{j}^{} \\ μ_{j}^{n e w} : = \frac{(A_{j}^{T} Σ^{- 1} r + s {}_{j}μ_{j}^{o l d} + β \sum_{i \in N_{j}} w_{j i}^{(t)} μ_{i}^{})}{s_{j} + β \sum_{i \in N_{j}} w_{j i}^{'}} \\ μ_{j}^{} : = \max {μ_{j}^{n e w}, 0} \\ r : = r + A_{j} (μ_{j}^{} - μ_{j}^{o l d}) \end{array}

End

μ_{}^{t + 1} : = μ

If ${{‖ μ^{(t + 1)} - μ^{(t)} ‖}_{2} / ‖ μ^{(t)} ‖}_{2} \leq T_{1}$ , Break, End

End

{\tilde{μ}}^{(k + 1)} : = μ

If ${{‖ {\tilde{μ}}^{(k + 1)} - {\tilde{μ}}^{(k)} ‖}_{2} / ‖ {\tilde{μ}}^{(k)} ‖}_{2} \leq T_{2}$ , Break, End

End

In our experiment, we set K=30. We found that the parameter $α$ should be updated in each MM iteration (i.e., M=1) since in this setting we always obtained slightly better reconstruction results with lower objective costs and more satisfactory visual effects.

Acknowledgment

This work was supported in part by the National Natural Science Foundation of China (NNSFC), under Grant Nos.60971112 and 61375018, and Fundamental Research Funds for the Central Universities, under Grant No. 2012QN086. J. Wang was supported in part by grants from the Cancer Prevention and Research Institute of Texas (RP130109 and RP110562-P2), the National Institute of Biomedical Imaging and Bioengineering (R01 EB020366) and a grant from the American Cancer Society (RSG-13-326-01-CCE). We would like to thank Dr. Damiana Chiavolini for editing the paper.

References and links

1. D. A. Jaffray, J. H. Siewerdsen, J. W. Wong, and A. A. Martinez, “Flat-panel cone-beam computed tomography for image-guided radiation therapy,” Int. J. Radiat. Oncol. Biol. Phys. 53(5), 1337–1349 (2002). [CrossRef] [PubMed]

2. A. Berrington de González, M. Mahesh, K.-P. Kim, M. Bhargavan, R. Lewis, F. Mettler, and C. Land, “Projected cancer risks from computed tomographic scans performed in the United States in 2007,” Arch. Intern. Med. 169(22), 2071–2077 (2009). [CrossRef] [PubMed]

3. N. Wen, H. Guan, R. Hammoud, D. Pradhan, T. Nurushev, S. Li, and B. Movsas, “Dose delivered from Varian’s CBCT to patients receiving IMRT for prostate cancer,” Phys. Med. Biol. 52(8), 2267–2276 (2007). [CrossRef] [PubMed]

4. L. Feldkamp, L. Davis, and J. Kress, “Practical cone-beam algorithm,” J. Opt. Soc. Am. A 1(6), 612–619 (1984). [CrossRef]

5. L. Ouyang, T. Solberg, and J. Wang, “Effects of the penalty on the penalized weighted least-squares image reconstruction for low-dose CBCT,” Phys. Med. Biol. 56(17), 5535–5552 (2011). [CrossRef] [PubMed]

6. K. Choi, J. Wang, L. Zhu, T.-S. Suh, S. Boyd, and L. Xing, “Compressed sensing based cone-beam computed tomography reconstruction with a first-order method,” Med. Phys. 37(9), 5113–5125 (2010). [CrossRef] [PubMed]

7. H. Lee, L. Xing, R. Davidi, R. Li, J. Qian, and R. Lee, “Improved compressed sensing-based cone-beam CT reconstruction using adaptive prior image constraints,” Phys. Med. Biol. 57(8), 2287–2307 (2012). [CrossRef] [PubMed]

8. S. Lefkimmiatis, A. Bourquard, and M. Unser, “Hessian-based regularization for 3-D microscopy image restoration,” in Biomedical Imaging (ISBI),20129th IEEE International Symposium on, (IEEE, 2012), 1731–1734. [CrossRef]

9. J. C. Park, B. Song, J. S. Kim, S. H. Park, H. K. Kim, Z. Liu, T. S. Suh, and W. Y. Song, “Fast compressed sensing-based CBCT reconstruction using Barzilai-Borwein formulation for application to on-line IGRT,” Med. Phys. 39(3), 1207–1217 (2012). [CrossRef] [PubMed]

10. E. Y. Sidky and X. Pan, “Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization,” Phys. Med. Biol. 53(17), 4777–4807 (2008). [CrossRef] [PubMed]

11. E. Y. Sidky, C.-M. Kao, and X. Pan, “Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT,” J. XRay Sci. Technol. 14, 119–139 (2006).

12. X. Jia, Y. Lou, R. Li, W. Y. Song, and S. B. Jiang, “GPU-based fast cone beam CT reconstruction from undersampled and noisy projection data via total variation,” Med. Phys. 37(4), 1757–1760 (2010). [CrossRef] [PubMed]

13. T. Chan, A. Marquina, and P. Mulet, “High-order total variation-based image restoration,” SIAM J. Sci. Comput. 22(2), 503–516 (2000). [CrossRef]

14. J. Yuan, C. Schnörr, and G. Steidl, “Total-variation based piecewise affine regularization,” in Scale Space and Variational Methods in Computer Vision (Springer, 2009), pp. 552–564.

15. K. Bredies, K. Kunisch, and T. Pock, “Total generalized variation,” SIAM J. Imaging Sci. 3(3), 492–526 (2010). [CrossRef]

16. Y.-L. You and M. Kaveh, “Fourth-order partial differential equations for noise removal,” IEEE Trans. Image Process. 9(10), 1723–1730 (2000). [CrossRef] [PubMed]

17. Y. Hu and M. Jacob, “Higher degree total variation (HDTV) regularization for image recovery,” IEEE Trans. Image Process. 21(5), 2559–2571 (2012). [CrossRef] [PubMed]

18. S. Lefkimmiatis, J. P. Ward, and M. Unser, “Hessian Schatten-norm regularization for linear inverse problems,” IEEE Trans. Image Process. 22(5), 1873–1888 (2013). [CrossRef] [PubMed]

19. S. Lefkimmiatis, A. Bourquard, and M. Unser, “Hessian-based norm regularization for image restoration with biomedical applications,” IEEE Trans. Image Process. 21(3), 983–995 (2012). [CrossRef] [PubMed]

20. J. Liu, T.-Z. Huang, I. W. Selesnick, X.-G. Lv, and P.-Y. Chen, “Image restoration using total variation with overlapping group sparsity,” Inf. Sci. 295, 232–246 (2015). [CrossRef]

21. T. Sun, N. Sun, J. Wang, and S. Tan, “Iterative CBCT reconstruction using Hessian penalty,” Phys. Med. Biol. 60(5), 1965–1987 (2015). [CrossRef] [PubMed]

22. K. Papafitsoros and C.-B. Schönlieb, “A combined first and second order variational approach for image reconstruction,” J. Math. Imaging Vis. 48(2), 308–338 (2014). [CrossRef]

23. G. Steidl, “Combined first and second order variational approaches for image processing,” Jahresber. Dtsch. Math. Ver. 117(2), 133–160 (2015). [CrossRef]

24. A. Chambolle and P.-L. Lions, “Image recovery via total variation minimization and related problems,” Numer. Math. 76(2), 167–188 (1997). [CrossRef]

25. M. Lysaker and X.-C. Tai, “Iterative image restoration combining total variation minimization and a second-order functional,” Int. J. Comput. Vis. 66(1), 5–18 (2006). [CrossRef]

26. S. Do, W. C. Karl, M. K. Kalra, T. J. Brady, and H. Pien, “A variational approach for reconstructing low dose images in clinical helical CT,” in Biomedical Imaging: From Nano to Macro, 2010 IEEE International Symposium on, (IEEE, 2010), 784–787. [CrossRef]

27. J. Dutta, S. Ahn, C. Li, S. R. Cherry, and R. M. Leahy, “Joint L1 and total variation regularization for fluorescence molecular tomography,” Phys. Med. Biol. 57(6), 1459–1476 (2012). [CrossRef] [PubMed]

28. J. M. Bioucas-Dias and M. A. Figueiredo, “An iterative algorithm for linear inverse problems with compound regularizers,” in Image Processing,2008. 15th IEEE International Conference on, (IEEE, 2008), 685–688. [CrossRef]

29. F. Chen, Y. Jiao, G. Ma, and Q. Qin, “Hybrid regularization image deblurring in the presence of impulsive noise,” J. Vis. Commun. Image Represent. 24(8), 1349–1359 (2013). [CrossRef]

30. D. R. Hunter and K. Lange, “A tutorial on MM algorithms,” Am. Stat. 58(1), 30–37 (2004). [CrossRef]

31. J. Wang, H. Lu, Z. Liang, D. Eremina, G. Zhang, S. Wang, J. Chen, and J. Manzione, “An experimental study on the noise properties of x-ray CT sinogram data in Radon space,” Phys. Med. Biol. 53(12), 3327–3341 (2008). [CrossRef] [PubMed]

32. J. Wang, T. Li, Z. Liang, and L. Xing, “Dose reduction for kilovotage cone-beam computed tomography in radiation therapy,” Phys. Med. Biol. 53(11), 2897–2909 (2008). [CrossRef] [PubMed]

33. Z. Tian, X. Jia, K. Yuan, T. Pan, and S. B. Jiang, “Low-dose CT reconstruction via edge-preserving total variation regularization,” Phys. Med. Biol. 56(18), 5949–5967 (2011). [CrossRef] [PubMed]

34. Y. Long, J. A. Fessler, and J. M. Balter, “3D forward and back-projection for X-ray CT using separable footprints,” IEEE Trans. Med. Imaging 29(11), 1839–1850 (2010). [CrossRef] [PubMed]

35. W. Dong, G. Shi, and X. Li, “Nonlocal image restoration with bilateral variance estimation: a low-rank approach,” IEEE Trans. Image Process. 22(2), 700–711 (2013). [CrossRef] [PubMed]

36. D. Stsepankou, A. Arns, S. K. Ng, P. Zygmanski, and J. Hesser, “Evaluation of robustness of maximum likelihood cone-beam CT reconstruction with total variation regularization,” Phys. Med. Biol. 57(19), 5955–5970 (2012). [CrossRef] [PubMed]

37. Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment: from error visibility to structural similarity,” IEEE Trans. Image Process. 13(4), 600–612 (2004). [CrossRef] [PubMed]

38. J. F. Pambrun and R. Noumeir, “Limitations of the SSIM quality metric in the context of diagnostic imaging,” in Image Processing (ICIP),2015IEEE International Conference on, 2015, 2960–2963. [CrossRef]

39. J. Wang, T. Li, and L. Xing, “Iterative image reconstruction for CBCT using edge-preserving prior,” Med. Phys. 36(1), 252–260 (2009). [CrossRef] [PubMed]

40. K. Bredies and T. Valkonen, “Inverse problems with second-order total generalized variation constraints,” Proceedings of Sampta (2011).

41. B. R. Whiting, P. Massoumzadeh, O. A. Earl, J. A. O’Sullivan, D. L. Snyder, and J. F. Williamson, “Properties of preprocessed sinogram data in x-ray computed tomography,” Med. Phys. 33(9), 3290–3303 (2006). [CrossRef] [PubMed]

42. E. Chouzenoux, A. Jezierska, J. Pesquet, and H. Talbot, “A convex approach for image restoration with exact Poisson-Gaussian likelihood,” SIAM J. Imaging Sci. 8(4), 2662–2682 (2015). [CrossRef]

43. J. Li, Z. Shen, R. Yin, and X. Zhang, “A reweighted L² method for image restoration with Poisson and mixed Poisson-Gaussian noise,” Inverse Probl. Imaging (Springfield) 9(3), 875–894 (2015). [CrossRef]

44. B. R. Whiting, P. Massoumzadeh, O. A. Earl, J. A. O’Sullivan, D. L. Snyder, and J. F. Williamson, “Properties of preprocessed sinogram data in x-ray computed tomography,” Med. Phys. 33(9), 3290–3303 (2006). [CrossRef] [PubMed]

45. A. Repetti, E. Chouzenoux, and J. C. Pesquet, “A penalized weighted least squares approach for restoring data corrupted with signal-dependent noise,” in Signal Processing Conference (EUSIPCO),2012Proceedings of the 20th European, 2012, 1553–1557.

46. H. Zhang, L. Ouyang, J. Ma, J. Huang, W. Chen, and J. Wang, “Noise correlation in CBCT projection data and its application for noise reduction in low-dose CBCT,” Med. Phys. 41(3), 031906 (2014). [CrossRef] [PubMed]

47. A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” SIAM J. Imaging Sci. 2(1), 183–202 (2009). [CrossRef]

48. S. Becker, J. Bobin, and E. J. Candès, “NESTA: a fast and accurate first-order method for sparse recovery,” SIAM J. Imaging Sci. 4(1), 1–39 (2011). [CrossRef]

49. P. B. Noël, A. M. Walczak, J. Xu, J. J. Corso, K. R. Hoffmann, and S. Schafer, “GPU-based cone beam computed tomography,” Comput. Methods Programs Biomed. 98(3), 271–277 (2010). [CrossRef] [PubMed]

Incident photon number	Penalty	Noise level( × 10⁻⁴)	PSNR(dB)	ISNR(dB)	SSIM
5 × 10³	FDK	12.0	26.57	0.00	0.4057
	TV	1.76	41.06	14.48	0.9529
	Hessian	1.66	44.06	17.49	0.9811
	TV-H	1.65	44.09	17.52	0.9812
1 × 10⁴	FDK	9.01	29.75	0.00	0.5765
	TV	1.29	42.58	12.83	0.9661
	Hessian	1.34	45.99	16.24	0.9871
	TV-H	1.32	46.03	16.27	0.9873
5 × 10⁴	FDK	5.27	35.61	0.00	0.8350
	TV	0.55	45.67	10.05	0.9816
	Hessian	0.53	53.21	17.60	0.9974
	TV-H	0.52	53.25	17.63	0.9976

ROI	Noise Level ( × 10⁻⁴)	1.56	2.56	3.64	4.64	5.61
ROI1	FDK (80mA)	3.58
	FDK(10mA)	1.47
	TV	17.68	11.72	8.61	6.85	5.64
	Hessian	21.41	13.36	9.35	7.38	6.06
	TV-H	21.59	13.50	9.26	7.20	5.86
ROI2	FDK (80mA)	2.58
	FDK (10mA)	0.92
	TV	13.52	8.42	5.72	4.43	3.66
	Hessian	14.40	8.80	6.04	4.50	3.85
	TV-H	14.55	8.94	6.19	4.81	3.82
ROI3	FDK (80mA)	2.02
	FDK (10mA)	0.78
	TV	11.17	7.17	5.14	4.01	3.28
	Hessian	11.48	6.86	4.80	3.84	3.21
	TV-H	11.57	6.98	4.94	3.99	3.30
ROI4	FDK (80mA)	2.51
	FDK (10mA)	0.90
	TV	13.95	8.96	6.29	4.90	4.02
	Hessian	13.53	8.79	6.08	4.78	3.94
	TV-H	13.65	8.92	6.22	4.93	3.99

Incident photon number	Penalty	Noise level( × 10⁻⁴)	PSNR(dB)	ISNR(dB)	SSIM
5 × 10³	FDK	12.0	26.57	0.00	0.4057
	TV	1.76	41.06	14.48	0.9529
	Hessian	1.66	44.06	17.49	0.9811
	TV-H	1.65	44.09	17.52	0.9812
1 × 10⁴	FDK	9.01	29.75	0.00	0.5765
	TV	1.29	42.58	12.83	0.9661
	Hessian	1.34	45.99	16.24	0.9871
	TV-H	1.32	46.03	16.27	0.9873
5 × 10⁴	FDK	5.27	35.61	0.00	0.8350
	TV	0.55	45.67	10.05	0.9816
	Hessian	0.53	53.21	17.60	0.9974
	TV-H	0.52	53.25	17.63	0.9976

ROI	Noise Level ( × 10⁻⁴)	1.56	2.56	3.64	4.64	5.61
ROI1	FDK (80mA)	3.58
	FDK(10mA)	1.47
	TV	17.68	11.72	8.61	6.85	5.64
	Hessian	21.41	13.36	9.35	7.38	6.06
	TV-H	21.59	13.50	9.26	7.20	5.86
ROI2	FDK (80mA)	2.58
	FDK (10mA)	0.92
	TV	13.52	8.42	5.72	4.43	3.66
	Hessian	14.40	8.80	6.04	4.50	3.85
	TV-H	14.55	8.94	6.19	4.81	3.82
ROI3	FDK (80mA)	2.02
	FDK (10mA)	0.78
	TV	11.17	7.17	5.14	4.01	3.28
	Hessian	11.48	6.86	4.80	3.84	3.21
	TV-H	11.57	6.98	4.94	3.99	3.30
ROI4	FDK (80mA)	2.51
	FDK (10mA)	0.90
	TV	13.95	8.96	6.29	4.90	4.02
	Hessian	13.53	8.79	6.08	4.78	3.94
	TV-H	13.65	8.92	6.22	4.93	3.99

Structure-adaptive CBCT reconstruction using weighted total variation and Hessian penalties

Abstract

1. Introduction

2. Methods

2.1 PWLS reconstruction

2.2. Adaptive structure reconstruction

2.3. Optimization method

2.3.1 Minimization step

2.3.2 Update step

2.4. Details in calculating projection matrix A

3. Experiments

3.1. Evaluation criteria

3.2. CS phantom

3.2.1 Convergence analysis

3.2.2 Evaluation

3.3. Modified 3D Shepp-Logan phantom

3.4. CatPhan 600 phantom

3.5. Anthropomorphic head phantom

4. Discussion

5. Conclusion

Appendix A derivative filters

Appendix B $w_{j i}^{(t)}$ and its update

Appendix C structure adaptive CBCT reconstruction using TV-H penalty

Initialization

Iteration

Acknowledgment

References and links

Cited By

Figures (12)

Tables (2)

Equations (32)

Biomedical Optics Express

Abstract

1. Introduction

2. Methods

2.1 PWLS reconstruction

2.2. Adaptive structure reconstruction

2.3. Optimization method

2.3.1 Minimization step

2.3.2 Update step

2.4. Details in calculating projection matrix A

3. Experiments

3.1. Evaluation criteria

3.2. CS phantom

3.2.1 Convergence analysis

3.2.2 Evaluation

3.3. Modified 3D Shepp-Logan phantom

3.4. CatPhan 600 phantom

3.5. Anthropomorphic head phantom

4. Discussion

5. Conclusion

Appendix A derivative filters

Appendix B wji(t)and its update

Appendix C structure adaptive CBCT reconstruction using TV-H penalty

Initialization

Iteration

Acknowledgment

References and links

Cited By

Figures (12)

Tables (2)

Equations (32)

Biomedical Optics Express

Appendix B $w_{j i}^{(t)}$ and its update