Noise correlation in CBCT projection data and its application for noise reduction in lowdose CBCT Hua Zhang, Luo Ouyang, Jianhua Ma, Jing Huang, Wufan Chen, and Jing Wang Citation: Medical Physics 41, 031906 (2014); doi: 10.1118/1.4865782 View online: http://dx.doi.org/10.1118/1.4865782 View Table of Contents: http://scitation.aip.org/content/aapm/journal/medphys/41/3?ver=pdfcov Published by the American Association of Physicists in Medicine Articles you may be interested in Iterative image-domain decomposition for dual-energy CT Med. Phys. 41, 041901 (2014); 10.1118/1.4866386 Database-assisted low-dose CT image restoration Med. Phys. 40, 031109 (2013); 10.1118/1.4790693 Noise suppression in reconstruction of low-Z target megavoltage cone-beam CT images Med. Phys. 39, 5111 (2012); 10.1118/1.4737116 Inverse determination of the penalty parameter in penalized weighted least-squares algorithm for noise reduction of low-dose CBCT Med. Phys. 38, 4066 (2011); 10.1118/1.3600696 Projection space denoising with bilateral filtering and CT noise modeling for dose reduction in CT Med. Phys. 36, 4911 (2009); 10.1118/1.3232004

Noise correlation in CBCT projection data and its application for noise reduction in low-dose CBCT Hua Zhang Department of Biomedical Engineering, Southern Medical University, Guangzhou, Guangdong 510515, China and Department of Radiation Oncology, University of Texas Southwestern Medical Center, Dallas, Texas 75390

Luo Ouyang Department of Radiation Oncology, University of Texas Southwestern Medical Center, Dallas, Texas 75390

Jianhua Ma,a) Jing Huang, and Wufan Chen Department of Biomedical Engineering, Southern Medical University, Guangzhou, Guangdong 510515, China

Jing Wanga) Department of Radiation Oncology, University of Texas Southwestern Medical Center, Dallas, Texas 75390

(Received 22 October 2013; revised 15 January 2014; accepted for publication 2 February 2014; published 19 February 2014) Purpose: To study the noise correlation properties of cone-beam CT (CBCT) projection data and to incorporate the noise correlation information to a statistics-based projection restoration algorithm for noise reduction in low-dose CBCT. Methods: In this study, the authors systematically investigated the noise correlation properties among detector bins of CBCT projection data by analyzing repeated projection measurements. The measurements were performed on a TrueBeam onboard CBCT imaging system with a 4030CB flat panel detector. An anthropomorphic male pelvis phantom was used to acquire 500 repeated projection data at six different dose levels from 0.1 to 1.6 mAs per projection at three fixed angles. To minimize the influence of the lag effect, lag correction was performed on the consecutively acquired projection data. The noise correlation coefficient between detector bin pairs was calculated from the corrected projection data. The noise correlation among CBCT projection data was then incorporated into the covariance matrix of the penalized weighted least-squares (PWLS) criterion for noise reduction of low-dose CBCT. Results: The analyses of the repeated measurements show that noise correlation coefficients are nonzero between the nearest neighboring bins of CBCT projection data. The average noise correlation coefficients for the first- and second-order neighbors are 0.20 and 0.06, respectively. The noise correlation coefficients are independent of the dose level. Reconstruction of the pelvis phantom shows that the PWLS criterion with consideration of noise correlation (PWLS-Cor) results in a lower noise level as compared to the PWLS criterion without considering the noise correlation (PWLS-Dia) at the matched resolution. At the 2.0 mm resolution level in the axial-plane noise resolution tradeoff analysis, the noise level of the PWLS-Cor reconstruction is 6.3% lower than that of the PWLS-Dia reconstruction. Conclusions: Noise is correlated among nearest neighboring detector bins of CBCT projection data. An accurate noise model of CBCT projection data can improve the performance of the statistics-based projection restoration algorithm for low-dose CBCT. © 2014 American Association of Physicists in Medicine. [http://dx.doi.org/10.1118/1.4865782] Key words: low-dose CBCT, noise correlation, PWLS

1. INTRODUCTION Onboard x-ray cone-beam computed tomography (CBCT) provides volumetric information of a patient in the treatment position using cone-beam acquisition geometry and a flat panel detector (FPD). There is increasing interest in using CBCT for patient treatment position setup and dose reconstruction for adaptive radiation therapy.1–3 Due to repeated CBCT acquisitions over a radiotherapy treatment course, there has also been concern about the extra CBCT imaging dose delivered to patients,4–6 especially for pediatric patients. One cost-effective way to reduce CBCT imaging dose is to 031906-1

Med. Phys. 41 (3), March 2014

acquire the projection data with a lower mAs protocol. However, the image quality of low-mAs CBCT degrades dramatically due to the excessive x-ray quantum noise. Much effort has been devoted to improving the image quality of CBCT acquired with low mAs protocol.7–9 Strategies based on statistical projection smoothing followed by analytical Feldkamp-Davis-Kress (FDK) (Ref. 10) reconstruction have shown promising results for low-dose CBCT noise reduction.11 Statistics-based projection smoothing algorithms try to estimate the true line integrals of CBCT through minimizing an objective function that is constructed from the noise statistical properties of CBCT projection data. The

0094-2405/2014/41(3)/031906/10/$30.00

© 2014 Am. Assoc. Phys. Med.

031906-1

Zhang et al.: Noise correlation in CBCT projection

031906-2

performance of the statistics-based algorithms strongly depends on the accurate noise modeling of measured projection data. One commonly used objective function is based on the penalized weighted least-squares (PWLS) criterion.12, 13 In the PWLS criterion, the second-order statistics of the projection data is incorporated into the weighting term, where the covariance matrix of the noise in projection data is used to weight the squared difference between the measured noisy data and ideal signal to be estimated. If the noise is not correlated among different detector bins, then the covariance matrix is reduced to a diagonal form, whose element is the inverse of the noise variance of each projection measurement. Previously developed statistics-based image restoration and reconstruction algorithms for CBCT (Refs. 11, 14, and 15) simply assume that there is no noise correlation among detector bins. Such an assumption has been validated in fan-beam multidetector CT (MDCT), where analyses on repeated projection data of a chest phantom using a Siemens MDCT scanner have shown that the noise of sinogram data in MDCT is not correlated among either bin- or view-direction.16 However, such an assumption is not valid for CBCT projection noise, where charge sharing may lead to noise correlation among nearest bins of the CBCT projection data acquired with a FPD. The purpose of this work is to systematically investigate the noise correlation among detector bins of CBCT projection data and to quantify the gain by incorporating noise correlation information into the PWLS algorithm for noise reduction of low-dose CBCT. We analyzed the noise correlation coefficient of CBCT projection data from repeated measurements of a pelvis phantom at different dose levels and varying projection views. Lag correction was performed on the transmitted raw data to avoid the influence of image lag on the noise modeling. With the obtained noise correlation coefficients, a nondiagonal noise covariance matrix is constructed and incorporated into the PWLS projection smoothing algorithm.12 The noise covariance matrix provides additional information of the measured noisy data, and the performance of a conventional PWLS algorithm can be improved with a more accurate noise model. An anthropomorphic pelvis phantom study is performed to quantitatively evaluate the gain of PWLS by using the nondiagonal covariance matrix as the weighting term.

resents the Poisson statistics of the photon counting with mean measurement λ¯ i (Ek ) collected by the flat panel detector (FPD). The second term in Eq. (1) represents the electronic noise, which is adequately modeled as a Gaussian distributed 2 . The offset random variable with mean di and variance σe,i mean di can be estimated prior to each scan and subtracted from the measurement.20 As shown in previous studies,17, 18 the linear relationship between the signal mean and variance still holds for the compound Poisson noise model. Additionally, the exact and closed-form solution is not available to compound Poisson likelihood function. Thus, in practice, the polychromatic nature of x rays is often not explicitly considered in present statistics-based image reconstruction for CT.21 Without considering the polychromatic nature of x-ray generation, the following simplified Poisson noise model is often used:   2 I˜i ∼ Poisson{λi } + Normal di , σe,i , (2)

031906-2

2. METHODS AND MATERIALS 2.A. CBCT measurement model

For polychromatic x-ray generation, the noise in CBCT transmission data acquired by an energy integrating detector can be modeled as the compound Poisson noise of photon counting plus the electronic Gaussian noise.17, 18 Mathematically, the transmitted CBCT measurement I˜i at a detector bin i can be described by the following model:17–19    2 I˜i ∼ , (1) G(Ek )Poisson{λ¯ i (Ek )} + Normal di , σe,i k

where k represents the energy-bin index, and G(·) is the mean energy dependent gain factor. The notation Poisson{·} repMedical Physics, Vol. 41, No. 3, March 2014

where λi represents the mean measure of the radiation intensity collected by the detector. By taking logarithm transformation of the transmitted data, the line integral of the linear attenuation coefficients along each ray path at detector bin i can be calculated by pi = ln

I0i = ln I0i − ln I˜i , I˜i

(3)

where the value of pi refers to the log-transformed projection datum at the detector bin i. I0i is the average number of photons emitted from the x-ray source and can be measured in the absence of the object in the FOV, i.e., by an air scan. I˜i represents the number of transmitted photons at the detector bin i with the mean value λi , thus, the calculated line integral pi is also a random variable. In this study, the variance σi2 of the line integral pi is estimated by the following mean-variance relationship proposed by Ma et al.:19 var(pi ) =

1 σ 2 − 1.25 + e 2 , λi λi

(4)

where σe2 is the background electronic noise variance. 2.B. Noise correlation analysis in CBCT projection data

2.B.1. Lag correction

Ideally, under the constant irradiation, the response of a FPD is expected to be constant over time. However, exponential rise in the rising step-response function (RSRF) (Ref. 22) can be seen in the repeated projection scans, as shown in Fig. 1. This is caused by the image lag. Lag affects the calculation of signal mean and variance of the projection data.23 In addition, in CBCT projections, lag causes temporal noise blurring, and leads to a range of reconstructed image artifacts, such as streaking, skin-line artifacts, or shading artifacts.22 Therefore, in this study, lag correction is first implemented on the transmitted raw data before the noise properties analysis and CBCT reconstruction.

Zhang et al.: Noise correlation in CBCT projection

031906-3

bins from the repeated projection measurements described in Sec. 2.D. The noise correlation coefficient is defined as cov(pi , pj ) cov(pi , pj )  , (10) =√ ρij = σi σj var(pi ) var(pj )

1.01

Normalized RSRF

1

0.99

where cov(pi , pj ) denotes the covariance of repeated logtransformed projection data p between the ith and the jth detector bins, while σ i and σ j are the standard deviations of the log-transformed projection data of the ith and the jth detector bins, respectively. The correlation coefficient is bounded by −1 ≤ ρ ≤ 1. ρ = 0 when there is no correlation and ρ = ±1 when pi and pj has perfectly positive correlation or negative correlation.

Uncorrected RSRF Correction

0.98

0.97

0.96

0.95

031906-3

0

50

100

150

200

250

300

350

400

450

500

Frame# F IG . 1. Normalized rising step responses of the measured data at 0.6 mAs dose level and its corresponding RSRF based corrected data.

In this study, an empirical method for lag correction is used.22, 24 The model of lag decay is described by an impulse response function (IRF) h(k) = b0 δ(k) +

N 

bn e−an k ,

(5)

n=1

where k is the discrete-time variable for x-ray frames and N = 2 is the number of exponential terms used in this study. δ(k) is the impulse function, and b0 δ(k) represents the portion of the input signal that is not affected by lag. bn represents the lag coefficients and an refers to the lag rates. Given h(k) defined in Eq. (5), the actual output L(k) of a FPD is L(k) = I (k)∗ h(k),

(6)

where I(k) denotes the ideal output signal without lag. In this work, the IRF is determined from the RSRF. The measured RSRF(k) is related to the IRF in Eq. (5) by   N  bn −an k RSRF(k) = 1 − e u(k). (7) 1 − e−an n=1 By fitting the normalized measurements to the model in Eq. (7), we get the lag coefficients bn and the lag rates an . Then the IRF is generated with the obtained parameters. The correction is a temporal iterative deconvolution of the measurements with the modeled IRF22, 24  −an L(k) − N n=1 bn Sn,k e I (k) = Ik = , (8) N n=0 bn Sn,k = Ik + Sn,k e−an .

(9)

Here, Sn, k is a state variable that contains the contribution from previous inputs. 2.B.2. Noise correlation analysis

To analyze the noise correlation of the projection data, we calculated the correlation coefficients of noise among detector Medical Physics, Vol. 41, No. 3, March 2014

2.C. Penalized weighted least-squares criterion incorporating noise correlation

2.C.1. Cost function of PWLS

Noise correlation in measured projection data can be incorporated into the PWLS criterion.12 In the projection domain, the PWLS objective function can be described as −1 1 β ˆT ˆ + R(p). (11) (p) = (yˆ − p) (yˆ − p) 2 2 The first term in Eq. (11) is a weighted least-squares criterion, where yˆ is the vector of the log-transformed projection data and pˆ is the vector of the ideal projection image to be estimated. Symbol T denotes the transpose operator. The ma trix is the noise covariance matrix. The covariance matrix offers knowledge of the second-order statistics. When the noise  correlation among projection data is ignored, the matrix is diagonal and its element is simply the variance of each measured data, while the noise variance determines the reliability and contribution of each measurement to the PWLS objective function. When the correlation information  of the noise in projection data is considered, the matrix becomes nondiagonal  var(p) = ⎤ ⎡ cov(p1 , p2 ) · · · cov(p1 , pN ) σ12 ⎥ ⎢ σ22 · · · cov(p2 , pN )⎥ ⎢ cov(p1 , p2 ) ⎥ ⎢ =⎢ ⎥. .. .. .. .. ⎥ ⎢ . . . . ⎦ ⎣ cov(pN , p1 ) cov(p2 , pN ) · · · σN2 (12) From the correlation coefficient ρ ij , the element cov(pi , pj ) of the covariance matrix can be estimated by cov(pi , pj ) = σij2 = ρij σi σj ,

(13)

where the standard deviations σ i of each √ measurement at the detector bin i can be estimated with σi = var(pi ). The noise correlation provides additional information of the measurements among neighboring detectors. Consideration of the noise correlation can potentially improve the performance of the PWLS criterion for noise reduction in low-dose CBCT. For simplification, we refer to the PWLS

Zhang et al.: Noise correlation in CBCT projection

031906-4

criterion with consideration of noise correlation as PWLSCor and the PWLS criterion ignoring noise correlation (i.e., with diagonal covariance matrix) as PWLS-Dia. The second term in Eq. (11) is a smoothness penalty or a priori constraint, where β is a parameter that controls the strength of the regularization. A quadratic form penalty has been used widely for iterative image reconstruction and sinogram restoration,12, 13   2 R(p) = R(pj ) = w(k, j )(pj − pk ) , (14)

form (p) is convex with linearity property. This implies that the minimizer to (p) can be found by searching for the stationary point of the corresponding gradient ∇(p). The gradient ∇(p) can be calculated as follows. The gradient of the quadratic penalty term can be rewritten by ∇R(p) = ∇ p 2 = Ap where A is a quint-diagonal matrix and satisfies ⎧ i=j ⎪ ⎨ 4, Ai,j = −1, i − j = ±1 or i − j = ±B . (15) ⎪ ⎩ 0, otherwise

031906-4

j

j

k∈Sj

where index j runs over all image elements. Sj represents the nearest neighbors around pixel j. The weight w(k, j ) is positive and symmetric, i.e., w(k, j ) ≥ 0 and w(k, j ) = w(j, k), which are usually designed to be an inverse proportion of the distance between pixels k and j in Sj . In this study, equal weight was used for the four nearest neighbors. 2.C.2. Optimization of the PWLS-Cor criterion

The nondiagonal covariance matrix used in Eq. (11) can be expressed as a real symmetric nine-diagonal matrix with  are nonzero when bandwidth B + 1, i.e., the entries σ ij of i = j, i − j = ±1, i − j = ±B, i − j = ±(B ± 1), and zero otherwise, as shown in the results of the noise correlation analysis in Sec. 3.C. B depends on the way in which the data taken from a two-dimensional detector array are rearranged to form a one-dimensional projection data vector. B is equal to the column dimension when columnwise rearrangement is used and equal to the row dimension when rowwise rearrangement is used. The resulting covariance matrix of the one-dimensional projection data vector from repeated measurement is very sparse and its structure is illustrated in Fig. 2. Recall that we are to minimize the cost function in Eq. (11) with nondiagonal covariance matrix where the penalty term R(p) is a convex function. Thus, the functional

Then the gradient of (p) is given by   −1 β 1 T (y − p) + ∇R(p) ∇(p) = ∇ (y − p) 2 2 −1 (p − y) + βAp. (16) = Solving the gradient equation is equivalent to solving the following linear system: −1 (p − y) + βAp = 0, (17) the solution of which is the minimizer of Eq. (11).  To fully exploit the structure of matrices and A, we can solve the linear system iteratively by utilizing the following algorithm with the available sparse solvers, for example, and A are large and sparse enough. The lsqr lsqr,25 as the method is based on the bidiagonalization procedure25 and it is analytically equivalent to the standard method of conjugate gradients, but possesses more favorable numerical properties. The present algorithm is implemented in MATLAB (Math Works, Inc., Natick, MA) on a 64-bit Windows PC R Core 2 Duo CPU of 3.2 GHz. For each projecwith Intel tion view, the computation time is about 10 s for 20 iterations with projection matrix size of 1024 × 768. A LGORITHM : Solving

−1

(p − y) + βAp = 0 iteratively.

Initialize with guess p(0) and set Iteration = 20; For n = 1: Iteration  Solve b = y − p(n) for b Solve βAp(n+1) = b for p(n+1) End

2.D. Data acquisition

F IG . 2. Nine-diagonal structure of covariance matrix of projection data vector. Medical Physics, Vol. 41, No. 3, March 2014

A TrueBeamTM onboard imaging system (Varian Medical Systems, Palo Alto, CA) with PaxScan 4030CB FPD (Varian Medical Systems, Palo Alto, CA) was used for the data acquisition. The FPD is an amorphous silicon kV image detector using a cesium iodide (CsI) scintillator. The active rectangular imaging area is 397 × 298 mm with a pixel pitch of 0.388 × 0.388 mm. The kV image detector is optimized for kV beam energies and provides an integrated 14 bit analog digital conversion and an optical interface. In this study, the half-fan mode with a half bow tie filter was used for the data acquisition. CBCT projection data were acquired on an

Zhang et al.: Noise correlation in CBCT projection

031906-5

anthropomorphic male pelvis phantom (CIRS Model 801-PF). The phantom is designed for use in diagnostic radiology and radiation therapy, and the anatomical dimensions of the phantom are based on the Visible Human Project (VHP) data sets. Data acquisition was performed with a tube voltage of 125 kVp. Six different dose levels were chosen for imaging the phantom from 0.1 to 1.6 mAs per projection, including 0.1, 0.3, 0.6, 0.9, 1.2, and 1.6 mAs. At each dose level, 500 repeated projections were acquired in the dynamic gain fluoro mode at three fixed gantry angles: 0◦ , 45◦ , 90◦ . CBCT scans of the same phantom were also performed at the same dose levels. Each scan included about 680 projection views evenly distributed over 360◦ . In the CBCT acquisition, the source-toaxis distance (SAD) was 1000 mm and the source-to-imager distance (SID) was 1500 mm.

acquisitions. After lag correction using the RSRF-calibrated IRF, a very flat RSRF was obtained, which indicated that excellent lag correction was achieved.

031906-5

2.E. Performance evaluation

The influence of different weight matrices on the PWLS projection restoration algorithm was quantified by a noise-resolution tradeoff measurement.26, 27 Axial-plane and coronal-plane image resolution were analyzed using an edge spread function (ESF) along profiles through the selected regions in the phantom. The ESF is a measure of the broadening of a step edge, and can be described by an error function (erf) parameterized by t, as ¯ y = r + Herf((x − x)/t),

From the repeated measurements of transmitted projection data at each mAs level, we calculated the sample mean and variance at each detector bin. Then we plotted the mean value against the variance value of each data point of a randomly selected region of interest (ROI) in the projection data image. These data points were fitted to a linear function, using a minimum least-squares fitting algorithm. The results are shown in Fig. 3. Figures 3(a)–3(d) show the results at different dose levels (0.3, 0.6, 0.9, and 1.2 mAs) overlaid with each one’s respective linear fitting. It can be observed that the mean and variance of measured raw data are linear, which reflects the compound Poisson noise nature of x-ray photons. It also can be seen that the intercept is nonzero with a positive value, which is caused by the background electronic noise. The noise floor is approximately constant, independent of the dose level. Thus, the variance of background electronic noise σe2 is set to 19 ADU for projection variance estimation in Eq. (3).

(18)

where x is the location of the pixel, and y represents the intensity of the pixel. r is the shift of the erf from baseline of zero, H is the magnitude of the half of the step, and x¯ is the center of the edge. By fitting the profiles to the error function, we can obtain parameter t which was used as the measure of the image resolution. The noise level of each reconstructed image was characterized by the standard deviation of a uniform region of 20 × 20 pixels. The sharp edges and the uniform regions were marked in Figs. 8(g) and 8(h). The smoothing parameter β controls the strength of the penalty term. A larger β value produces a more “smoothed” reconstructed image with less noise. A set of β values were chosen for each PWLS algorithm to get the noise-resolution tradeoff curves of the reconstructed pelvis phantom images. The range of β values was selected in such a way that the noise level in the PWLS-processed low-dose CBCT was close to the noise level in the high-dose CBCT (1.6 mAs) reconstructed by FDK. 3. RESULTS 3.A. Lag correction

In Fig. 1, the normalized rising step responses of the uncorrected and RSRF-corrected image intensity were both plotted as a function of the frame number. Before lag correction, the average intensity of projection data increased over time although the acquisition protocols were kept the same. This indicates that lag occurs during the consecutive projection data Medical Physics, Vol. 41, No. 3, March 2014

3.B. Linear relationship between the signal mean and variance of transmitted projection data

3.C. Noise correlation coefficients of projection data

A region of interest of size 21 × 21 pixels was chosen in the CBCT projection data to calculate the noise correlation coefficients. Figure 4(a) shows the noise correlation coefficients of CBCT projection from repeated scans at the 0.6 mAs dose level. For a better illustration, Fig. 4(b) shows a zoomed-in image of Fig. 4(a). It can be observed that there are nonzero values in addition to the principal diagonal values. The nonzero values are mainly along the nearest neighbors in the two-dimension projection data. This observation is different from that of the conventional fan-beam CT where noise is not correlated among neighboring detector bins.12 To quantitatively measure the noise correlation coefficients, a profile was obtained by averaging the correlation coefficients matrix along its diagonal direction at each dose level. The profiles of the lowest (0.1 mAs) and the highest (1.6 mAs) dose levels are shown in Figs. 5(a) and 5(b) as examples. Nonzero noise correlation coefficient values can be found between the eight nearest neighboring pixels. For the first order neighbors, the horizontal noise correlation coefficients are slightly larger than the vertical ones. The average noise correlation coefficient for the first order neighbors is 0.20. The average correlation coefficient for the second-order neighbors is 0.06, which is much smaller than that of the firstorder neighbors. Figure 6 shows amplitudes of the eight neighboring noise correlation coefficients as a function of dose level. It can be seen that the noise correlation coefficients are almost constant at different dose levels. Figure 7 shows amplitudes of the eight neighboring noise correlation coefficients as a function

Zhang et al.: Noise correlation in CBCT projection

031906-6

031906-6

3000

6000

0.3 mAs Linear fit of measured data Intercept=19.8 R-squared=0.9870

2500

4000

Variance

Variance

2000

1500

3000

1000

2000

500

1000

0

0.6 mAs Linear fit of measured data Intercept=19.9 R-squared=0.9878

5000

0

500

1000

1500

2000

0

2500

0

500

1000

1500

2000

Mean (a) 10000

3000

3500

4000

4500

8000

9000 10000

12000

0.9 mAs Linear fit of measured data Intercept=18.9 R-squared=0.9881

9000 8000 7000

1.2 mAs Linear fit of measured data Intercept=18.8 R-squared=0.9866

10000

8000

6000

Variance

Variance

2500

Mean (b)

5000 4000

6000

4000

3000 2000

2000

1000 0

0

1000

2000

3000

4000

5000

6000

7000

8000

Mean (c)

0

0

1000

2000

3000

4000

5000

6000

7000

Mean (d)

F IG . 3. Linear relationship between mean and variance of the transmitted data. (a)–(d) correspond to the dose level dose at 0.3, 0.6, 0.9, and 1.2 mAs.

of view angle for three view angles at 0◦ , 45◦ , 90◦ . Constant correlation coefficient values are observed among different view angles, where the same pixel in the projection image at different view angles has different intensity levels. These results suggest that the noise correlation coefficients are independent of the dose level and image intensity. Thus in this study, the noise correlation coefficient is assumed to be constant across the detector for the same neighboring detector bin. 3.D. Reconstruction analysis

Figures 8(a) and 8(b) show FDK reconstructed images from projection data acquired with 1.6 mAs protocol. Figures 8(c)–8(h) show FDK reconstructed CBCT images from projection data acquired with 0.6 mAs protocol. The left column shows the axial view and the right column shows the coronal view. Figures 8(c) and 8(d) are reconstructed from the noisy projection data without noise treatment, where excessive noise makes it difficult to identify soft tissue structures from the background. Figures 8(e) and 8(f) are from the projection data processed by the PWLS-Dia criterion and Figs. 8(g) and 8(h) are from the projection data processed by the PWLS-Cor criterion. Noise is substantially reduced Medical Physics, Vol. 41, No. 3, March 2014

in CBCT from projection data processed by both PWLSCor and PWLS-Dia criteria. When the same smoothing parameter β was used, the CBCT noise level from the PWLSCor processed projection is lower than that from the PWLSDia processed projection. Zoomed ROIs are also shown for a better visual inspection. The performance difference between the PWLS-Cor and PWLS-Dia criteria is further quantified by the noiseresolution tradeoff measure. Figure 9 shows both the axial-plane [Fig. 9(a)] and coronal-plane [Fig. 9(b)] noise resolution curves by analyzing the axial-view image and coronalview image reconstructed from the PWLS processed projection data, either with the PWLS-Cor criterion or with the PWLS-Dia criterion. The noise-resolution data point of the high dose image (1.6 mAs) reconstructed by the FDK algorithm was used as a reference for the choice of β values and performance comparison. At matched resolution, PWLS-Cor has a lower noise level than PWLS-Dia in both curves. At the 2.0 mm resolution level in the axial-plane noise-resolution tradeoff analysis, the noise level of PWLS-Cor reconstruction is 6.3% lower than that of PWLS-Dia reconstruction. These results show that the PWLS-Cor criterion outperforms PWLS-Dia criterion by incorporating noise correlation information among neighboring detector bins.

Zhang et al.: Noise correlation in CBCT projection

031906-7

Correlation Coefficient

031906-7

0.25 0.2 0.15 0.1 0.05 0 0

0.5

1

1.5

Dose(mAs)

Correlation Coefficient

F IG . 6. Noise correlation coefficients between a pixel and its eight nearest neighboring pixels at different dose levels.

0.25 0.2 0.15 0.1 0.05

F IG . 4. (a) Correlation coefficients matrix of noise among detector bins from repeated measurements at 0.6 mAs dose level. (b) Zoomed in image of (a) for a better illustration.

0 0

45

90

View Angle F IG . 7. Noise correlation coefficients between a pixel and its eight nearest neighboring pixels at different view angles.

1.2 Dose level: 0.1 mAs

1 Correlation Coefficient

Correlation Coefficient

1 0.8 0.6 0.4 0.2 0 -0.2 350

Dose level: 1.6 mAs

0.8 0.6 0.4 0.2 0

400

450

(a)

500

550

-0.2 350

400

450

500

550

(b)

F IG . 5. Profiles of averaged correlation coefficients matrix along its diagonal direction. The lowest [(a) 0.1 mAs] and the highest [(b) 1.6 mAs] dose levels were shown for illustration. Medical Physics, Vol. 41, No. 3, March 2014

031906-8

Zhang et al.: Noise correlation in CBCT projection

031906-8

(a)

(b)

(c)

((d)

((e)) (e

((f)

(g)

(h)

F IG . 8. (a) and (b) show FDK reconstructed CBCT images from projection data acquired with 1.6 mAs protocol. (c)–(h) show the results from projection data acquired with 0.6 mAs protocol: (c) and (d) are the reconstructions from the noisy projection data; (e) and (f) are the reconstructions from the projection data processed by the PWLS-Dia criterion (β = 0.5 × 103 ); (g) and (h) are the reconstructions from the projection data processed by PWLS-Cor criterion (β = 0.5 × 103 ). The left column shows axial view and the right column shows the coronal view. The display window is [−620 680] HU.

2.5

2.6 PWLS-Cor PWLS-Dia FDK-1.6 mAs

2.4

2.4

Resolution (mm)

Resolution (mm)

2.3 2.2 2.1 2 1.9

2.3 2.2 2.1 2 1.9

1.8 1.7 1.5

PWLS-Cor PWLS-Dia FDK-1.6 mAs

2.5

1.8 1.7 2

2.5

3

3.5

4

Standard Deviation(attn.coeff(1/mm))

(a)

4.5

5 -4

x 10

3

3.5

4

4.5

5

Standard Deviation(attn.coeff(1/mm))

5.5 -4

x 10

(b)

F IG . 9. Axial (a) and coronal (b) noise-resolution tradeoff curves of the reconstructed images for the two PWLS based projection restoration algorithms at 0.6 mAs dose level. Medical Physics, Vol. 41, No. 3, March 2014

031906-9

Zhang et al.: Noise correlation in CBCT projection

(a)

031906-9

(b)

(c)

F IG . 10. FDK reconstructed CBCT images before and after lag correction from projection data acquired with 0.6 mAs protocol. (a) Image reconstructed from the uncorrected projection data; (b) image reconstructed from the corrected projection data; (c) different image between images (a) and (b). Display window for (a) and (b) is [−620 680] HU, and for (c) is [−45 50] HU.

4. DISCUSSION AND CONCLUSION Statistics-based image reconstruction and restoration algorithms for CBCT (Refs. 11, 14, and 15) usually model the noise in CBCT projection data at each detector bin as conditionally independent variable. However, the indirect flat plane detector in kilovoltage (kV) CBCT uses the active readout mechanisms28 and cross talk (i.e., charge sharing) that occurs among neighboring detector bins. Thus, the signal of one detector bin may alter the noise properties of its neighbors, resulting in noise correlation among neighboring detector bins. The noise correlation contains additional information of the measurements and such information can be used to improve accuracy from the noisy measurements. The primary objective of this work is to systematically analyze noise correlation in CBCT projection data. In this study, we performed lag correction on the acquired projection data before any further data analysis. The main purpose of lag correction is to make the calculation of the signal mean, variance, and the correlation coefficient unbiased and more accurate from the repeated scans. Without lag correction, the increasing trend of the signal during repeated scans will introduce a bias in the mean and variance calculation. Additionally, lag correction can reduce shading artifacts in the final reconstructed images. An example in Fig. 10 illustrates the benefit of lag correction. Figure 10(a) is the image reconstructed from the uncorrected projection data while Fig. 10(b) is the image reconstructed from the projection data after lag correction. Figure 10(c) shows the difference between images in Figs. 10(a) and 10(b), where strong artifacts are observed. For the PWLS criterion, the choice of the penalty function is also very important. In this study, to focus on the influence of the weighting matrix on the PWLS criterion and minimize the influence of the penalty term, a simple isotropic penalty function was used, i.e., simple differences of nearest neighbors were considered in the penalty term. Other edgerpreserving regulation terms, such as Huber penalty,29 total variation,30 and anisotropic quadratic penalty,14 can be coupled with the WLS with a covariance weighting matrix and the performance of PWLS is expected to be further improved. In addition to the penalty term, the performance of the PWLS criterion also strongly depends on the choice of the regularization parameter β, which controls the tradeoff between the data fidelity term and smoothness penalty term. In this study, β was empirically selected by matching the noise level Medical Physics, Vol. 41, No. 3, March 2014

of processed low-dose CBCT to the corresponding high dose images. Other automatic methods can be used for determining optimal β value, such as the generalized crossvalidation (GCV) (Ref. 31) and the L-curve criterion.32 In summary, we have performed experimental studies to systematically investigate the noise correlation properties of CBCT projection data. The noise correlation coefficients among detector bins were calculated with repeated scanned measurements. The analyses showed that noise correlation coefficients are nonzero between the nearest neighboring bins of CBCT projection data, and that they are independent of dose level. The obtained noise correlation coefficients were then applied to the PWLS criterion. The reconstructed results showed that when the off-diagonal correlation terms were included in the covariance matrix, the PWLS criterion outperformed the one with only diagonal weighted terms.

ACKNOWLEDGMENTS This work was supported in part by grants from the National Natural Science Foundation of China (Nos. 81371544, 81000613, and 81101046), the National Key Technology Research and Development Program of the Ministry of Science and Technology of China (No. 2011BAI12B03), and the 973 Program of China under Grant No. 2010CB732503. J. Wang was supported in part by grants from Cancer Prevention and Research Institute of Texas (RP110562-P2 and RP130109) and the American Cancer Society (RSG-13-326-01-CCE). The authors would like to thank Ms. Dee Hill for editing the paper. The authors would also like to thank the anonymous reviewers for their constructive comments and suggestions that greatly improved the quality of the paper.

a) Authors

to whom correspondence should be addressed. Electronic addresses: [email protected] and [email protected] 1 Y. Yang, E. Schreibmann, T. Li, C. Wang, and L. Xing, “Evaluation of on-board kV cone beam CT (CBCT)-based dose calculation,” Phys. Med. Biol. 52, 685–705 (2007). 2 L. Xing, B. Thorndyke, E. Schreibmann, Y. Yang, T. F. Li, G. Y. Kim, G. Luxton, and A. Koong, “Overview of image-guided radiation therapy,” Med. Dosim. 31, 91–112 (2006). 3 L. Lee, Q. T. Le, and L. Xing, “Retrospective IMRT dose reconstruction based on cone-beam CT and MLC log-file,” Int. J. Radiat. Oncol., Biol., Phys. 70, 634–644 (2008).

031906-10

Zhang et al.: Noise correlation in CBCT projection

4 D.

J. Brenner and E. J. Hall, “Current concepts – Computed tomography – An increasing source of radiation exposure,” N. Engl. J. Med. 357, 2277– 2284 (2007). 5 N. Wen, H. Q. Guan, R. Hammoud, D. Pradhan, T. Nurushev, S. D. Li, and B. Movsas, “Dose delivered from Varian’s CBCT to patients receiving IMRT for prostate cancer,” Phys. Med. Biol. 52, 2267–2276 (2007). 6 M. K. Islam, T. G. Purdie, B. D. Norrlinger, H. Alasti, D. J. Moseley, M. B. Sharpe, J. H. Siewerdsen, and D. A. Jaffray, “Patient dose from kilovoltage cone beam computed tomography imaging in radiation therapy,” Med. Phys. 33, 1573–1582 (2006). 7 L. Ouyang, T. Solberg, and J. Wang, “Noise reduction in low-dose cone beam CT by incorporating prior volumetric image information,” Med. Phys. 39, 2569–2577 (2012). 8 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, 5955–5970 (2012). 9 H. Zhang, Y. Liu, H. Han, J. Wang, J. H. Ma, L. H. Li, and Z. R. Liang, “A comparison study of sinogram- and image-domain penalized re-weighted least-squares approaches to noise reduction for low-dose cone-beam CT,” Proc. SPIE 8668, 86683E (2013). 10 L. A. Feldkamp, L. C. Davis, and J. W. Kress, “Practical cone-beam algorithm,” J. Opt. Soc. Am. A 1, 612–619 (1984). 11 J. Wang, T. F. Li, Z. R. Liang, and L. Xing, “Dose reduction for kilovoltage cone-beam computed tomography in radiation therapy,” Phys. Med. Biol. 53, 2897–2909 (2008). 12 J. Wang, T. Li, H. B. Lu, and Z. R. Liang, “Penalized weighted least-squares approach to sinogram noise reduction and image reconstruction for lowdose X-ray computed tomography,” IEEE Trans. Med. Imaging 25, 1272– 1283 (2006). 13 J. A. Fessler, “Penalized weighted least-squares image-reconstruction for positron emission tomography,” IEEE Trans. Med. Imaging 13, 290–300 (1994). 14 J. Wang, T. Li, and L. Xing, “Iterative image reconstruction for CBCT using edge-preserving prior,” Med. Phys. 36, 252–260 (2009). 15 J. Tang, B. E. Nett, and G. H. Chen, “Performance comparison between total variation (TV)-based compressed sensing and statistical iterative reconstruction algorithms,” Phys. Med. Biol. 54, 5781–5804 (2009). 16 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, 3327–3341 (2008). 17 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, 3290–3303 (2006).

Medical Physics, Vol. 41, No. 3, March 2014

031906-10 18 S.

Y. Huang, K. Yang, C. K. Abbey, and J. M. Boone, “A semiempirical linear model of indirect, flat-panel x-ray detectors,” Med. Phys. 39, 2108– 2118 (2012). 19 J. H. Ma, Z. R. Liang, Y. Fan, Y. Liu, J. Huang, W. F. Chen, and H. B. Lu, “Variance analysis of x-ray CT sinograms in the presence of electronic noise background,” Med. Phys. 39, 4051–4065 (2012). 20 J. Hsieh, Computer Tomography, Principles, Design, Artifacts, and Recent Advances (SPIE, Bellingham, WA, 2003). 21 J. B. Thibault, K. D. Sauer, C. A. Bouman, and J. Hsieh, “A threedimensional statistical approach to improved image quality for multislice helical CT,” Med. Phys. 34, 4526–4544 (2007). 22 J. Starman, J. Star-Lack, G. Virshup, E. Shapiro, and R. Fahrig, “Investigation into the optimal linear time-invariant lag correction for radar artifact removal,” Med. Phys. 38, 2398–2411 (2011). 23 K. Yang, S. Y. Huang, N. J. Packard, and J. M. Boone, “Noise variance analysis using a flat panel x-ray detector: A method for additive noise assessment with application to breast CT applications,” Med. Phys. 37, 3527– 3537 (2010). 24 J. Hsieh, O. E. Gurmen, and K. F. King, “A recursive correction algorithm for detector decay characteristics in CT,” Proc. SPIE 3977, 298–305 (2000). 25 C. C. Paige and M. A. Saunders, “LSQR: An algorithm for sparse linearequations and sparse least-squares,” ACM Trans. Math. Softw. 8, 43–71 (1982). 26 J. D. Evans, D. G. Politte, B. R. Whiting, J. A. O’Sullivan, and J. F. Williamson, “Noise-resolution tradeoffs in x-ray CT imaging: A comparison of penalized alternating minimization and filtered backprojection algorithms,” Med. Phys. 38, 1444–1458 (2011). 27 J. Y. Xu and B. M. W. Tsui, “Electronic noise modeling in statistical iterative reconstruction,” IEEE Trans. Image Process. 18, 1228–1238 (2009). 28 W. Zhao and J. A. Rowlands, “Digital radiology using active matrix readout of amorphous selenium: Theoretical analysis of detective quantum efficiency,” Med. Phys. 24, 1819–1833 (1997). 29 W. Chlewicki, F. Hermansen, and S. B. Hansen, “Noise reduction and convergence of Bayesian algorithms with blobs based on the Huber function and median root prior,” Phys. Med. Biol. 49, 4717–4730 (2004). 30 C. R. Vogel and M. E. Oman, “Iterative methods for total variation denoising,” SIAM J. Sci. Comput. 17, 227–238 (1996). 31 G. H. Golub, M. Heath, and G. Wahba, “Generalized cross-validation as a method for choosing a good ridge parameter,” Technometrics 21, 215–223 (1979). 32 P. C. Hansen and D. P. Oleary, “The use of the L-curve in the regularization of discrete ill-posed problems,” SIAM J. Sci. Comput. 14, 1487–1503 (1993).