- Email: [email protected]
Resolution and noise performance of sparse view X-ray CT reconstruction via Lp-norm regularization
Physica Medica 52 (2018) 72–80
Contents lists available at ScienceDirect
Physica Medica journal homepage: www.elsevier.com/locate/ejmp
Original paper
Resolution and noise performance of sparse view X-ray CT reconstruction via Lp-norm regularization
T
⁎
Lin Zhanga, Huijuan Zhaoa,b, Wenjuan Mac, Jingying Jiangd, , Limin Zhanga,b, Jiao Lia,b, ⁎ Feng Gaoa,b, Zhongxing Zhoua,b, a
School of Precision Instrument and Optoelectronics Engineering, Tianjin University, Tianjin, People’s Republic of China Tianjin Key Laboratory of Biomedical Detecting Techniques and Instruments, Tianjin, People’s Republic of China c Tianjin Medical University Cancer Institute and Hospital, Tianjin, People’s Republic of China d Paul C Lauterbur Research Center for Biomedical Imaging, Shenzhen Institutes of Advanced Technology, Chinese Academy of Sciences, Shenzhen, Guangdong 518055, People’s Republic of China b
A R T I C LE I N FO
A B S T R A C T
Keywords: Sparse view Lp-norm regularization Modulation transfer function (MTF) Noise power spectrum (NPS) Noise equivalent quanta (NEQ)
Objectives: Adaptive steepest descent projection onto convex sets (ASD-POCS) algorithms with Lp-norm (0 < p ≤ 1) regularization have shown great promise in sparse-view X-ray CT reconstruction. However, the difference in p value selection can lead to varying algorithm performance of noise and resolution. Therefore, it is imperative to develop a reliable method to evaluate the resolution and noise properties of the ASD-POCS algorithms under different Lp-norm priors. Methods: A comparative performance evaluation of ASD-POCS algorithms under different Lp-norm (0 < p ≤ 2) priors were performed in terms of modulation transfer function (MTF), noise power spectrum (NPS) and noise equivalent quanta (NEQ). Simulation data sets from the EGSnrc/BEAMnrc Monte Carlo system and an actual mouse data set were used for algorithms comparison. Results: A considerable MTF improvement can be achieved with the decrement of p. L1 regularization based algorithm obtains the best noise performance, and shows superiority in NEQ evaluation. The advantage of L1norm prior is also confirmed by the reconstructions from the actual mouse data set through contrast to noise ratio (CNR) comparison. Conclusion: Although the ASD-POCS algorithms using small Lp-norm (p ≤ 0.5) priors yield a higher MTF than do the high Lp-norms, the best noise-resolution performance is achieved when p is between 0.8 and 1. The results are expected to be a reference to the choice of p in Lp-norm (0 < p ≤ 2) regularization.
1. Introduction X-ray cone-beam computed tomography (CBCT) is an important tool for biomedical research and preclinical applications. In spite of the remarkable advantages of CBCT, the exposure risk remains a major concern in clinical practice. To decrease the radiation dose, techniques have been developed to realize the limited-angle CT [1–3] or limitedview CT [4–6]. However, these dose reduction approaches will unavoidably increase the data inconsistency, leading to a challenging task for image reconstruction. To deal with such situation, much effort has been directed to developing iterative image reconstruction algorithms for applications in X-ray CBCT with data inconsistency [7,8]. Iterative reconstruction methods produce good quality images when the projection data is not theoretically sufficient for exact image
reconstruction [9,10]. The CBCT reconstruction problem can be solved using iterative algorithms by formulating the data consistency constraint with an additional regularization term [10–13]. One representative algorithm is the adaptive steepest descent projection onto convex sets (ASD-POCS) algorithm [10], which uses the adaptive steepest descent algorithm for total variance minimization—L1-norm of the gradient magnitude images, and employs projection onto convex sets algorithm to keep data fidelity. In spite of the advantages of using a L1 regularization term to solve the sparsity constrained problems, enhancing the sparsity constraint for better imaging performance is still a main concern. CT image reconstructions from ASD-POCS algorithms using Lp-norm (0 < p < 1) prior are expected to have higher spatial resolution than those using L1-norm or L2-norm prior [14–16]. However, quantitative
⁎
Corresponding authors at: School of Precision Instrument and Optoelectronics Engineering, Tianjin University, Tianjin 300072, People’s Republic of China (Z. Zhou). E-mail addresses: [email protected] (L. Zhang), [email protected] (H. Zhao), [email protected] (J. Jiang), [email protected] (L. Zhang), [email protected] (J. Li), [email protected] (F. Gao), [email protected] (Z. Zhou). https://doi.org/10.1016/j.ejmp.2018.04.396 Received 12 January 2018; Received in revised form 28 March 2018; Accepted 25 April 2018 1120-1797/ © 2018 Published by Elsevier Ltd on behalf of Associazione Italiana di Fisica Medica.
Physica Medica 52 (2018) 72–80
L. Zhang et al.
energy (PCUT) was set to 521 keV and 10 keV, respectively. Rayleigh scattering and bound Compton scattering were switched on, along with atomic relaxations and electron impact ionization. The prediction of Xray spectrum was derived from the BEAM Data Processor (BEAMDP) program. Specifically speaking, an 80 keV mono-energetic parallel circular electron beam in the simulation was impinged to the target, and then the multiple spectrum containing energy from 0 to 80 KeV was generated. The simulated spectral distribution was presented in the left of Fig. 2. The spectra have been assigned 0.4 keV bin widths and normalized to unit area. CBCT scanning simulation was performed using egs_cbct, an EGSnrc-based application for CBCT-related calculations. Photons were simulated down to 1 keV and no electron was transported. The latter was achieved by setting ECUT higher than the maximum photon energy used in the MC simulation. CBCT scans of the phantom were conducted with an isotropic point X-ray source modeled by BEAMnrc. The derived spectrum of the X-ray source was presented in Fig. 2. An ideal detector made of 512 × 512 pixels with a size of 1.5 mm is adopted in the simulation. The X-ray source was located at 630 mm from the phantom’s center, and the detector was centered at 970 mm from the source along the line connecting the source and the phantom’s center. The sourcedetector system can be rotated 360° around the phantom. The phantom used in CBCT scanning simulation was constructed according to the third module of physical ACR phantom, which is a solid cylinder phantom constructed primarily from a uniform waterequivalent material. It is 8 cm in depth and 20 cm in diameter made of solid water, a material in the standard 521ICRU library of EGSnrc (see the right of Fig. 2). The digital ACR phantom was then reconstructed by the ASD-POCS algorithm using varying Lp-norm (0 < p ≤ 2) priors, based on TIGRE, a MATLAB-GPU toolbox for CBCT image reconstruction [33]. Three simulation data sets under three different signal to noise ratios (SNRs) were generated by running different particle histories in EGS simulation for all projection angles. Table 1 shows the simulation and reconstruction parameters. For all Lp-norm (0 < p ≤ 2) priors, iteration numbers were set to 300, which were enough to get the convergence.
frequency comparisons on resolution and noise performance of ASDPOCS algorithms using different Lp-norm (0 < p ≤ 2) priors are not sufficient at the present time [10,14,17]. Some papers were based on global evaluation metrics such as signal to noise ratio (SNR), contrast to noise ratio (CNR), and root of mean square error (RMSE) [18–20], and although there have been reports evaluating the spatial resolutions using modulation transfer function (MTF) analysis, they presented results within limited number of Lp-norm priors [21,22]. Moreover, relatively few works were conducted on quantitatively comparing both the resolution and noise performance by spatial frequency dependent metrics for the ASD-POCS algorithms with Lp-norm (0 < p ≤ 2) regularization. There, the aim of this study is to employ an overall and comprehensive performance comparison of Lp-norm (0 < p ≤ 2) priors. The modulation transfer function (MTF) and the noise power spectrum (NPS) have been widely recognized as the most relevant spatial frequency dependent metrics of resolution and noise evaluation in radiographic imaging. The MTF is a measure of image sharpness and the NPS describes the noise in the image. These two metrics have been introduced into CBCT imaging system evaluation with varying measurement methods [23,24]. MTF and NPS can be further combined to derive the noise equivalent quanta (NEQ), which is a decisive indicator of the signal and noise transfer properties of an imaging system [25]. In CT, methods for measuring the presampling MTF were developed by using an aluminum edge [26] or a tungsten wire [27]. These methods typically use highly dense materials to simulate an infinite impulse to maximize CNR and improve measurement precision. However, these materials are atypical in biomedical research and preclinical applications. Moreover, the use of high-density metallic test object for MTF measurement may introduce metal artifacts that can introduce additional uncertainties owing to their CT artifacts. According to this problem, it is important to use test objects with lower contrast to avoid metal artifacts and yield MTFs relevant to realistic one used in clinical images. A currently widely adopted low contrast phantom for MTF measurement in CT is the American College of Radiology (ACR) computed tomography accreditation program (further denoted as ACR phantom) [28]. The ACR phantom is a solid cylinder phantom constructed primarily from a water-equivalent material. Note that in 2011, low-contrast cylinder test objects similar to ACR phantom were proposed by the American National Standard for MTF measurement of CT [29]. In this standard, the circumference of the cylinder is used as an edge to measure the MTF, and the inner part of it is used to calculate the NPS. In this work, a Monte Carlo (MC) simulation platform was built employing the EGSnrc/BEAMnrc code system. The CBCT scans were simulated under the condition of different noise levels. Then the simulated and experiment data were employed to quantitatively investigate and compare the resolution and noise performance of Lp-norm (0 < p ≤ 2) regularization. Several figures of merit including the MTF, NPS and NEQ, which account for the spatial resolution and the noise performance, were introduced to objectively evaluate the reconstruction performance at different noise levels.
2.2. Evaluation metrics Frequency dependent metrics including MTF, NPS and NEQ were calculated for algorithm performance evaluation. MTF was assessed by a radial oversampling methodology [28] from the reconstructed images as illustrated in Fig. 3. The oversampled edge spread function (ESF) was determined on increments of the distance from the edge [28], and then it was fitted with the linear combination of two Gaussian functions and a Boltzmann function to avoid instability [34]. ESF was then differentiated to obtain the line spread function (LSF), from which the MTF was calculated by using the Fourier transform. To investigate the frequency variation of noise under Lp-norm (0 < p ≤ 2) regularization, the NPS was calculated based on ensemble statistics from the ACR phantom images. The NPS utilizes the Fourier transform of noise images to determine the variance of noise power present at each spatial frequency. Noise-only images were obtained by the subtraction of the averaged image from each of 20 central reconstructed slices. Total 720 region of interests (ROIs) in 20 noise-only images were averaged to calculate the 2D NPS distribution. ROIs were created, centered along a circle, as illustrated in Fig. 3 [28]. Finally, the desired NPS profile was generated by averaging 8 radial profiles from center to the edge to avoid fluctuation. The image quality of the reconstructed ACR phantom was further assessed by the measurement of the noise equivalent quanta (NEQ), which is a measure of the signal to noise ratio of an imaging system. Mathematically, the NEQ can be represented as [29]
2. Materials and methods 2.1. Virtual CT system The full MC simulation of a kV CBCT imaging system was carried out with the EGSnrc/BEAMnrc MC code system [30,31]. The component modules in the BEAMnrc simulation consisted of an 'XTUBE' (the X-ray target) and a 'SLABS' (the filter). The SLABS included a filter of 4.0 mm of Aluminum plus 0.1 mm of Beryllium. The tungsten anode target (2 mm thick) was contained in a copper holder, with 9 degree from z-axis (see Fig. 1). Cross sections for all materials in the simulation were generated from the XCOM data set using the PEGS4 code system [32]. The global electron cut-off energy (ECUT) and photon cut-off
NEQ (f ) = 73
[S·MTF (f )]2 NPS (f )
(1)
Physica Medica 52 (2018) 72–80
L. Zhang et al.
Fig. 1. X-ray tube geometry in BEAMnrc for x-z projection.
Fig. 2. Left: X-ray spectral distribution derived from BEAMDP, right: digital ACR phantom.
where f is spatial frequency. S is the signal level, measured as the average CT value of the object. To validate the conclusion obtained from the simulation, a mouse data set was then adopted for reconstruction. The mouse was fixed in an imaging cavity and was scanned with 42 views distributed evenly in 360 degrees. The study protocol was approved by the local institution review board and performed according to the guidelines of the committee on animal use and care of Chinese Academy of Sciences. The related parameters for image acquisition were shown in Table 2. The acquired projections were sheared and rebinned to 482 × 482 pixels with a pixel size of 0.0977 mm × 0.0977 mm. An example projection view of the mouse was shown in Fig. 4. The voxel size of the reconstruction was 0.123 mm × 0.123 mm × 0.123 mm for 241 × 241 × 241 voxels. The iteration numbers were also set to 300 for all Lpnorm (0 < p ≤ 2) priors, which were enough to get the convergence. To carry out a quantitative research, contrast to noise ratios (CNRs) of the reconstructions using different Lp-norm (0 < p ≤ 2) priors were calculated from the areas as indicated in Fig. 5. The red rectangle area denotes the target area and the green one represents the background area, they were selected to calculate the CNR as following [18]:
Fig. 3. Diagram of ROIs selecting in NPS calculation. Table 2 Acquisition parameters for real CT data set. Parameters
Voltage
Current
Exposure time
Pixel size
Pixel numbers
Value
45KV
0.3 mA
375 s
0.0488 mm
1024 × 1024
Table 1 Parameters setting in phantom simulation and image reconstruction. Parameters
Value
Projection parameters
Reconstruction parameters
Pixel size(mm)
Pixel numbers
Projection numbers
Voxel size(mm)
Voxel numbers
1.5 × 1.5
256 × 256
20
1.96 × 1.96 × 1.96
128 × 128
74
Physica Medica 52 (2018) 72–80
L. Zhang et al.
difference can be found within the range of p ∈ (0, 0.5]. When the value of p increases in the range of (1.5, 2], MTF changes slowly. The same trends can be seen from the MTFs obtained from the simulation data sets under different SNRs. The NPSs associated with different Lp-norm (0 < p ≤ 2) regularization at varying noise levels were given in Fig. 8. The left, middle and right column shows the NPS of reconstructed images from projections with SNR = 35 dB, 40 dB and 45 dB respectively. As shown in the figure, L1 regularization resulted in the lowest NPS among all Lp-norm (0 < p ≤ 2) regularizations under three different projection noise levels. The measured NPS was increased when the value of p deviated from 1. Fig. 9 compares the NEQ measurements on log scale at different noise levels from the reconstructions with various Lp-norm priors. Within the range of p ∈ [1, 2], L1-norm prior achieves the highest NEQ, and an obvious NEQ reduction can be found with the increment of the value of p. The NEQ curves for Lp-norm (1 < p ≤ 2) priors have a local maximum around the intermediate frequency within cut-off frequency (half the Nyquist sampling frequency). Within the range of p ∈ (0, 1), NEQ is increased when the value of p is approaching 1. No obvious difference among NEQ curves can be found within the range of p ∈ (0,0.5] in terms that the projection SNR is 35 dB or 40 dB. Generally speaking, ASD-POCS algorithm with L1-norm prior achieves the best NEQ among all results except for the most high frequency part approaching the cut-off frequency.
Fig. 4. An example projection view of the mouse.
3.2. Experiment Fig. 10 shows the mouse head cross section of a display window [0, 0.1] reconstructed with different Lp-norm (0 < p ≤ 2) priors. The oval structure in the middle is the skull of the mouse, and two circle structures in front of the skull are the front claws. The CNRs of the reconstructions using varying Lp-norm priors (p = 0.3, 0.5, 0.8, 1.0, 1.1, 1.2, 1.3, 1.5, 1.8, 2.0) were calculated as 22.13, 25.25, 57.64, 58.59, 49.94, 47.07, 46.27, 42.25, 36.90 and 32.64 respectively. The results demonstrated that the best CNR performance was obtained when p equals 1, which is in consistent with the NEQ results in the simulation studies (see Fig. 9). The CNR is a metric for the overall evaluation of contrast and noise in the image, while the NEQ provides detailed contrast and noise assessment related to the spatial frequency. In order to observe the difference of each algorithm, the image profiles along the horizontal and vertical lines, as shown in Fig. 5, are plotted in Fig. 11. The ROI is selected from 1D profiles and marked by the red rectangles in the left-hand column of Fig. 11. The profiles of these selected pixels are plotted in the right-hand column to observe the differences more clearly. It can obviously be seen that the profiles of the ASD-POCS algorithm with smaller Lp-norm priors show sharper edges than those using higher Lp-norm priors. However, the ASD-POCS algorithm using small Lp-norm (p ≤ 0.5) priors have an effect of introducing spikes at edges. These spikes behave as pepper-like impulse noise as shown in Fig. 10 (a) and (b), which may lead to difficulties in identifying small structures of the sample.
Fig. 5. Illustration to the calculation of CNR.
CNR =
It −Ib σb
(2)
where It represents for the average pixel value in the target area, Ib stands for the average pixel value in the background and σb is the square root of the background variance. 3. Results 3.1. Simulation Fig. 6 displays the central axial slices reconstructed by the ASDPOCS algorithm using varying Lp-norm priors. The simulation data set under the SNR of 40 dB was selected as a typical example. The display gray scale window of the figure was set to [0.0044, 0.0054] for better showing the edge and noise performance. As shown in the figure, ASDPOCS algorithms using small Lp-norm (0 < p < 1) priors preserved the high contrast edges of the image, while high Lp-norm (1 < p < 2) priors blurred the edges. When p ∈ [0.5, 1), the smaller the value of p, the sharper edger obtained by the Lp regularization. When p ∈ (0, 0.5], the performance of Lp-norm priors has no significant difference on edge enhancement. When p ∈ (1, 2], artifacts can be found in the reconstructed images. Obvious artifacts can be found in the reconstructed slices when the value of p deviated from 1. The subsequent calculated MTF based on the fitted ESFs are shown in Fig. 7. It can be observed that the MTF is improved when reducing the value of p. This demonstrates that better spatial resolution is achieved by the ASD-POCS algorithm using small Lp-norm (0 < p < 1) priors. Specifically speaking, the MTF improvement is significant when decreasing p from 1.5 to 0.5. No significant MTF
4. Discussion and conclusion Many researches have proved that an exact reconstruction of the sparse signal is possible with substantially fewer measurements by replacing the L1-norm with the Lp-norm (0 < p < 1) [14,35–37]. A few researchers perform the reconstruction using ASD-POCS algorithm with a L0-norm prior approximated by the Lp-norm, and compared the image reconstruction performance with those using L2-norm and L1-norm prior [14,38]. In this paper, the resolution and noise performances of ASD-POCS algorithms using different Lp-norm (0 < p ≤ 2) priors are assessed and compared by simulation and experiment. In the phantom study, the comparison via MTF showed that a 75
Physica Medica 52 (2018) 72–80
L. Zhang et al.
Fig. 6. Reconstructions of the ACR phantom from different p with SNR of 40 dB, central axial slice reconstructed using (a) p = 0.3 (b) p = 0.5 (c) p = 0.8 (d) p = 1.0 (e) p = 1.3 (f) p = 1.5 (g) p = 1.8 (h) p = 2.0. The display gray scale window is [0.0044, 0.0054].
speaking, a better NEQ in low spatial frequency is obtained when p = 1 while a higher NEQ appears at the high frequency part when p = 0.8. The resolution and noise performances of ASD-POCS algorithms was also tested along with application to an actual mouse data set. Experimental data set study shows that relative better spatial resolution and noise performance could be obtained when p is between 0.8 and 1, which is in consistent with the NEQ results in the simulation studies. What should be pointed out is that the main drawback of Lp (p > 1) regularization is not sparsity preserving [39], and as a consequence the corresponding reconstructions potentially suffer from noise and the overall performance is lower than that of p ≤ 1. This phenomenon is verified by the simulation and experiment in this study. Besides, the reconstruction time was not taken into consideration in this work. Although iterative reconstruction is usually characterized by a high computational cost, the steady improvement of the compute capabilities of GPUs can decrease the reconstruction time significantly, making them increasingly suitable for iterative reconstruction in CT [40,41].
smaller value of p resulted in a higher resolution, and when p was smaller than 0.5 or larger than 1.5, the difference in reconstructed resolution were not obvious. Thus, it could be concluded that the ASDPOCS algorithm with smaller Lp-norm prior has a higher capability to preserve edge details compared to those using bigger Lp-norm. However, the ASD-POCS algorithms using small Lp-norm (p ≤ 0.5) priors may result in inconsistent artifacts, while those with high Lpnorm (p > 1.2) priors will lead to blurred reconstructions. In our phantom study, the best noise performance was obtained around p = 1 for the simulation data sets under any SNR condition. It should be noted that as the SNR of projection views increased to 40 dB and 45 dB, no marked improvement was obtained on the MTF performance compared to the reconstruction results from projection views of 35 dB. To take both MTF and NPS into consideration, the noise-resolution tradeoff study is employed by calculating the NEQ of the ASD-POCS algorithms. Although the ASD-POCS algorithms using small Lp-norm (p ≤ 0.5) priors yield a higher MTF than do the high Lp-norms, the best NEQ performance is achieved when p is between 0.8 and 1. Specifically
Fig. 7. MTF based on fitted ESFs for the projections SNR of 35 dB (left), 40 dB (middle) and 45 dB (right).
76
Physica Medica 52 (2018) 72–80
L. Zhang et al.
Fig. 8. NPS reconstructed using different Lp-norm priors under projection SNRs of 35 dB (left), 40 dB (middle) and 45 dB (right).
Fig. 9. NEQ reconstructed with different p of different SNRs (left:35 dB, middle:40 dB, right:45 dB).
Fig. 10. Reconstructed images from real CT data set with different p: (a) p = 0.3 (b) p = 0.5 (c) p = 0.8 (d) p = 1.0 (e) p = 1.1 (f) p = 1.2 (g) p = 1.3 (h) p = 1.5. (i) p = 1.8. (j) p = 2.0. The display gray scale window is [0, 0.1].
77
Physica Medica 52 (2018) 72–80
L. Zhang et al.
Fig. 11. One-dimensional profiles along the horizontal line (the 180st row of the reconstructed images) and the vertical line (the 121st column of the reconstructed images.) (a) The image profiles along the horizontal line, (b) the partial profiles of the selected ROI, (c) the image profiles along the vertical line, (d) the partial profiles of the selected ROI.
Acknowledgments
system for early cancer detection. XCT imaging is an important part of this system. This study serves the goal of finding a high-performance algorithm for XCT reconstruction with limited data.
This work was supported by the National Natural Science Foundation of China (No. 81401453, 81371602, 61475115, 61475116, 61575140, 81571723, 81671728), and Tianjin Municipal Government of China (15JCZDJC31800, 16JCZDJC31200, 17JCZDJC32700).
Conflicts of interest
Role of the funding source
The authors have no potential conflicts of interest related with this article.
The role of the study sponsor is to develop multi-model imaging Appendix Iterative reconstruction algorithms based on Lp-norm (0 < p ≤ 2) regularization For simplicity, the discretized Beer-Lambert law is selected as the forward projection model for monochromatic x rays in a partially attenuating medium.
⎡ ⎤ Yi = di exp ⎢− ∑ hij μj ⎥ j ⎣ ⎦
(3)
where i and j are linear indices in the image vectors and projection data, respectively. di is the incident photon count. Yi is the detector photon count after attenuation by the object along the jth X-ray path. hij is a discrete length path through voxel j and μj denotes the voxel’s linear attenuation coefficient. From Eq. (3), a logarithmic operation is performed to obtain the so-called projection data along jth X-ray path
Pi = ln
di = Yi
∑
hij μj = Hi μ
(4)
j
78
Physica Medica 52 (2018) 72–80
L. Zhang et al.
A discrete linear system for modeling the CT-imaging process can be derived from Eq. (4) and written as
P = Hμ
(5)
where the vector P = {pi |i = 1,2,…,M } denote the projection data and the vector μ = {μ|i = 1,2,…,n} is the vectorized image to be reconstructed; H is the system matrix modeling the forward projection. With this discrete data model, the image reconstruction problem is equivalent to inverting the linear system Eq. (3), finding μ from a given projection data set P. However, even under ideal conditions, it is difficult to directly invert Eq. (5) when the system is under-determined, i.e., M < N, for image reconstruction from a reduced number of projections under study. Regularization is an effective way to constraint an under-determined problem. Similar to the previous studies [10], the CT image reconstruction problem can be formulated as a constrained optimization problem below
μ∗ = arg min‖μ‖pp
(6)
Subject to the inequality constraints: (a) data fidelity
|Hμ−P| ⩽ ε
(7)
and (b) non-negativity
μj ⩾ 0
(8)
where the user-defined parameter ε is interpreted as the variance of the difference between the real and the predicted of p to the Lp-norm of the gradient magnitude image, which can be written as
‖μ‖pp = (∂x μijk ) p + (∂y μijk ) p + (∂z μijk ) p
projections. ‖μ‖pp
is the power (9)
When p = 0, it’s exactly the L0-norm priors which equal to the number of nonzero elements. While minimization of the L0-norm is a combinatorial optimization problem and not practical [42]. Replacing p = 0 with p = 1 yields the L1-norm priors, it’s a convex problem and can be solved efficiently. L1 priors has been applied in image reconstruction in CT and there are plenty of impressive results [10,43]. Another choice is Lp-norm (0 < p < 1) priors, has proposed in image processing and image reconstruction since it promotes sparser results [44,45]. Reconstructions with L2norm (p = 2) priors has been proposed early, while a main drawback is the results with L2-norm priors are oversmoothing [15]. The above constrained optimization problem is extended from the ASD-POCS algorithm introduced by Sidky and Pan [9], which is modified by replacing the conventional L1-norm priors by an Lp-norm (0 < p ≤ 2) regularized version in Eq. (9). The algorithm alternatively minimizes the Lpnorm using adaptive steepest descent (ASD) and data fidelity using projection onto convex sets (POCS).
[19] Lu Y, Chan HP, Wei J, Hadjiiski LM. Selective-diffusion regularization for enhancement of microcalcifications in digital breast tomosynthesis reconstruction. Med Phys 2010;37(11):6003–14. [20] Huang J, Ma J, Liu N, Zhang H, Bian Z, Feng Y, et al. Sparse angular CT reconstruction using non-local means based iterative-correction POCS. Comput Biol Med 2011;41(4):195–205. [21] Hashemi SM, Song WY, Sahgal A, Lee Y, Huynh C, Grouza V, et al. Simultaneous deblurring and iterative reconstruction of CBCT for image guided brain radiosurgery. Phys Med Biol 2017;62(7):2521. [22] Gonzales B, Spronk D, Cheng Y, Tucker AW, Beckman M, Zhou O, et al. Rectangular fixedgantry CT prototype: combining CNT x-ray sources and accelerated compressed sensingbased reconstruction. IEEE Access 2014;2:971–81. [23] Han M, Park S, Baek J. Effect of anatomical noise on the detectability of cone beam ct images with different slice direction, slice thickness, and volume glandular fraction. Opt Express 2016;24(17):18843–59. [24] Gao Y, Bian Z, Huang J, Zhang Y, Niu S, Feng Q, et al. Low-dose X-ray computed tomography image reconstruction with a combined low-mAs and sparse-view protocol. Opt Express 2014;22(12):15190–210. [25] Tward DJ, Siewerdsen JH. Cascaded systems analysis of the 3D noise transfer characteristics of flat-panel cone-beam CT. Med Phys 2008;35(12):5510–29. [26] Boone JM. Determination of the presampled MTF in computed tomography. Med Phys 2001;28(3):356. [27] Richard S, Husarik DB, Yadava G, Merphy SN, Samei E. Towards task-based assessment of CT performance: system and object MTF across different reconstruction algorithms. Med Phys 2012;39(7Part1):4115–22. [28] Friedman SN, Fung GSK, Siewerdsen JH, Tsui BMW. A simple approach to measure computed tomography (CT) modulation transfer function (MTF) and noise-power spectrum (NPS) using the American College of Radiology (ACR) accreditation phantom. Med Phys 2013;40(5). [29] American national standard for evaluating the image quality of x-ray computed tomography (CT) security-screening systems, IEEE Standard ANSI N42.45, 2011. [30] Rogers DWO, Faddegon BA, Ding GX, Ma CM, We J, Mackie TR. BEAM: a Monte Carlo code to simulate radiotherapy treatment units. Med Phys 1995;22(5):503–24. [31] BEAMnrc users manual, National Research Council of Canada, Ottawa, ON, 2011. [32] The EGS4 Code System, Stanford Linear Accelerator Center, Stanford, CA, 1985. [33] Biguri A, Dosanjh M, Hancock S, Soleisani M. TIGRE: a MATLAB-GPU toolbox for CBCT image reconstruction. 2016;2(5):055010. [34] Brunner CC, Kyprianou IS. Material-specific transfer function model and SNR in CT. Phys Med Biol 2013;58(20):7447. [35] Chartrand R. Exact reconstruction of sparse signals via nonconvex minimization. IEEE Signal Process Lett 2007;14(10):707–10. [36] Chartrand R, Staneva V. Restricted isometry properties and nonconvex compressive sensing. Inverse Prob 2008;24(3):035020. [37] Sidky EY, Chartrand R, Boone JM, Pan X. Constrained TpV minimization for enhanced exploitation of gradient sparsity: application to CT image reconstruction. IEEE J Transl Eng Health Med 2014;2:1–18.
References [1] Qi H, Chen Z, Wu S, Xu Y, Zhou L. Iterative image reconstruction using modified non-local means filtering for limited-angle computed tomography. Phys Med 2016;32(9):1041–51. [2] Chen Z, Jin X, Li L, Wang G. A limited-angle CT reconstruction method based on anisotropic TV minimization. Phys Med Biol 2013;58(7):2119–41. [3] Protonotarios NE, Fokas AS, Gaitanis A, Kastis GA. Inversion of the attenuated radon transform using cubic splines. Phys Med 2016;32(9):242–3. [4] Yan B, Zhang W, Li L, Zhang H, Wang L. Quantitative study on exact reconstruction sampling condition by verifying solution uniqueness in limited-view CT. Phys Med 2016;32(1):1321–30. [5] Huang J, Ma J, Liu N, Zhang H, Bian Z, Feng Y, et al. Sparse angular CT reconstruction using non-local means based iterative-correction POCS. Comput Biol Med 2011;41(4):195–205. [6] Chen Z, Qi H, Wu S, Xu Y, Zhou L. Few-view CT reconstruction via a novel non-local means algorithm. Phys Med 2016;32(10):1276–83. [7] Slavine NV, Guild J, McColl RW, Anderson JA, Oz OK, Lenkinski RE. An iterative deconvolution algorithm for image recovery in clinical CT: a phantom study. Phys Med 2015;31(8):903–11. [8] Egiazarian K, Foi A, Katkovnik V. Compressed sensing image reconstruction via recursive spatially adaptive filtering. IEEE Int Conf Image Process 2007;1:I-549–52. [9] Mueller K, Yagel R, Wheller JJ. Anti-aliased three-dimensional cone-beam reconstruction of low-contrast objects with algebraic methods. IEEE Trans Med Imaging 2002;18(6):519–37. [10] Sidky EY, Pan X. Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization. Phys Med Biol 2008;53(17):4777. [11] Han X, Bian J, Eaker DR, Kline TL, Sidky EY, Ritman EL, et al. Algorithm-enabled lowdose micro-CT imaging. IEEE Trans Med Imaging 2011;30(3):606. [12] Niu T, Zhu L. Accelerated barrier optimization compressed sensing (ABOCS) reconstruction for cone-beam CT: phantom studies. IEEE; 2012. p. 2354–7. [13] Lee H, Xing L, Davidi R, Li R, Qian J, Lee R, et al. Improved compressed sensing-based cone-beam CT reconstruction using adaptive prior image constraints. Phys Med Biol 2012;57(8):2287–307. [14] Sidky EY, Pan X, Reiser IS, Reiser IS, Niskikawaet RM, et al. Enhanced imaging of microcalcifications in digital breast tomosynthesis through improved image-reconstruction algorithms. Med Phys 2009;36(11):4920. [15] Morigi S, Reichel L, Sgallari F. Fractional Tikhonov regularization with a nonlinear penalty term. J Comput Appl Math 2017;324:142–54. [16] Fu SJ, Zhang CM, Tai XC. Image denoising and deblurring: non-convex regularization, inverse diffusion and shock filter. Sci China Inf Sci 2011;54(6):1184–98. [17] Sidky EY, Duchin Y, Pan X, Ullberg C. A constrained, total-variation minimization algorithm for low-intensity x-ray CT. Med Phys 2011;38(S1). [18] Shieh CC, Kipritidis J, O’Brien RT, Kuncic Z, Keall PJ. Image quality in thoracic 4D conebeam CT: a sensitivity analysis of respiratory signal, binning method, reconstruction algorithm, and projection angular spacing. Med Phys 2014;41(4).
79
Physica Medica 52 (2018) 72–80
L. Zhang et al.
[42] Sidky EY, Chartrand R, Pan X. Image reconstruction from few views by non-convex optimization. 5. IEEE; 2007. p. 3526–30. [43] Chartrand R, Yin W. Iteratively reweighted algorithms for compressive sensing. IEEE; 2008. p. 3869–72. [44] Okawa S, Hoshi Y, Yamada Y. Improvement of image quality of time-domain diffuse optical tomography with lp sparsity regularization. Biomed Opt Express 2011;2(12):3334–48. [45] Chartrand R, Sidky EY, Pan X. Nonconvex compressive sensing for X-ray CT: an algorithm comparison. IEEE; 2013. p. 665–9.
[38] Hu Y, Xie L, Luo L, Nunes JC, Toumoulin C. L0 constrained sparse reconstruction for multi-slice helical CT reconstruction. Phys Med Biol 2011;56(4):1173. [39] Theodoridis S, Kopsinis Y, Slavakis K. Sparsity-aware learning and compressed sensing: an overview. arXiv preprint arXiv. 2012;1211;5231. [40] Jia X, Lou Y, Li R, Song WY, Jiang SB. GPU-based fast cone beam CT reconstruction from undersampled and noisy projection data via total variation. Med Phys 2010;37(4):1757–60. [41] Després P, Jia X. A review of GPU-based medical image reconstruction. Phys Med: Eur J Med Phys 2017;42:76–92.
80
Contents lists available at ScienceDirect
Physica Medica journal homepage: www.elsevier.com/locate/ejmp
Original paper
Resolution and noise performance of sparse view X-ray CT reconstruction via Lp-norm regularization
T
⁎
Lin Zhanga, Huijuan Zhaoa,b, Wenjuan Mac, Jingying Jiangd, , Limin Zhanga,b, Jiao Lia,b, ⁎ Feng Gaoa,b, Zhongxing Zhoua,b, a
School of Precision Instrument and Optoelectronics Engineering, Tianjin University, Tianjin, People’s Republic of China Tianjin Key Laboratory of Biomedical Detecting Techniques and Instruments, Tianjin, People’s Republic of China c Tianjin Medical University Cancer Institute and Hospital, Tianjin, People’s Republic of China d Paul C Lauterbur Research Center for Biomedical Imaging, Shenzhen Institutes of Advanced Technology, Chinese Academy of Sciences, Shenzhen, Guangdong 518055, People’s Republic of China b
A R T I C LE I N FO
A B S T R A C T
Keywords: Sparse view Lp-norm regularization Modulation transfer function (MTF) Noise power spectrum (NPS) Noise equivalent quanta (NEQ)
Objectives: Adaptive steepest descent projection onto convex sets (ASD-POCS) algorithms with Lp-norm (0 < p ≤ 1) regularization have shown great promise in sparse-view X-ray CT reconstruction. However, the difference in p value selection can lead to varying algorithm performance of noise and resolution. Therefore, it is imperative to develop a reliable method to evaluate the resolution and noise properties of the ASD-POCS algorithms under different Lp-norm priors. Methods: A comparative performance evaluation of ASD-POCS algorithms under different Lp-norm (0 < p ≤ 2) priors were performed in terms of modulation transfer function (MTF), noise power spectrum (NPS) and noise equivalent quanta (NEQ). Simulation data sets from the EGSnrc/BEAMnrc Monte Carlo system and an actual mouse data set were used for algorithms comparison. Results: A considerable MTF improvement can be achieved with the decrement of p. L1 regularization based algorithm obtains the best noise performance, and shows superiority in NEQ evaluation. The advantage of L1norm prior is also confirmed by the reconstructions from the actual mouse data set through contrast to noise ratio (CNR) comparison. Conclusion: Although the ASD-POCS algorithms using small Lp-norm (p ≤ 0.5) priors yield a higher MTF than do the high Lp-norms, the best noise-resolution performance is achieved when p is between 0.8 and 1. The results are expected to be a reference to the choice of p in Lp-norm (0 < p ≤ 2) regularization.
1. Introduction X-ray cone-beam computed tomography (CBCT) is an important tool for biomedical research and preclinical applications. In spite of the remarkable advantages of CBCT, the exposure risk remains a major concern in clinical practice. To decrease the radiation dose, techniques have been developed to realize the limited-angle CT [1–3] or limitedview CT [4–6]. However, these dose reduction approaches will unavoidably increase the data inconsistency, leading to a challenging task for image reconstruction. To deal with such situation, much effort has been directed to developing iterative image reconstruction algorithms for applications in X-ray CBCT with data inconsistency [7,8]. Iterative reconstruction methods produce good quality images when the projection data is not theoretically sufficient for exact image
reconstruction [9,10]. The CBCT reconstruction problem can be solved using iterative algorithms by formulating the data consistency constraint with an additional regularization term [10–13]. One representative algorithm is the adaptive steepest descent projection onto convex sets (ASD-POCS) algorithm [10], which uses the adaptive steepest descent algorithm for total variance minimization—L1-norm of the gradient magnitude images, and employs projection onto convex sets algorithm to keep data fidelity. In spite of the advantages of using a L1 regularization term to solve the sparsity constrained problems, enhancing the sparsity constraint for better imaging performance is still a main concern. CT image reconstructions from ASD-POCS algorithms using Lp-norm (0 < p < 1) prior are expected to have higher spatial resolution than those using L1-norm or L2-norm prior [14–16]. However, quantitative
⁎
Corresponding authors at: School of Precision Instrument and Optoelectronics Engineering, Tianjin University, Tianjin 300072, People’s Republic of China (Z. Zhou). E-mail addresses: [email protected] (L. Zhang), [email protected] (H. Zhao), [email protected] (J. Jiang), [email protected] (L. Zhang), [email protected] (J. Li), [email protected] (F. Gao), [email protected] (Z. Zhou). https://doi.org/10.1016/j.ejmp.2018.04.396 Received 12 January 2018; Received in revised form 28 March 2018; Accepted 25 April 2018 1120-1797/ © 2018 Published by Elsevier Ltd on behalf of Associazione Italiana di Fisica Medica.
Physica Medica 52 (2018) 72–80
L. Zhang et al.
energy (PCUT) was set to 521 keV and 10 keV, respectively. Rayleigh scattering and bound Compton scattering were switched on, along with atomic relaxations and electron impact ionization. The prediction of Xray spectrum was derived from the BEAM Data Processor (BEAMDP) program. Specifically speaking, an 80 keV mono-energetic parallel circular electron beam in the simulation was impinged to the target, and then the multiple spectrum containing energy from 0 to 80 KeV was generated. The simulated spectral distribution was presented in the left of Fig. 2. The spectra have been assigned 0.4 keV bin widths and normalized to unit area. CBCT scanning simulation was performed using egs_cbct, an EGSnrc-based application for CBCT-related calculations. Photons were simulated down to 1 keV and no electron was transported. The latter was achieved by setting ECUT higher than the maximum photon energy used in the MC simulation. CBCT scans of the phantom were conducted with an isotropic point X-ray source modeled by BEAMnrc. The derived spectrum of the X-ray source was presented in Fig. 2. An ideal detector made of 512 × 512 pixels with a size of 1.5 mm is adopted in the simulation. The X-ray source was located at 630 mm from the phantom’s center, and the detector was centered at 970 mm from the source along the line connecting the source and the phantom’s center. The sourcedetector system can be rotated 360° around the phantom. The phantom used in CBCT scanning simulation was constructed according to the third module of physical ACR phantom, which is a solid cylinder phantom constructed primarily from a uniform waterequivalent material. It is 8 cm in depth and 20 cm in diameter made of solid water, a material in the standard 521ICRU library of EGSnrc (see the right of Fig. 2). The digital ACR phantom was then reconstructed by the ASD-POCS algorithm using varying Lp-norm (0 < p ≤ 2) priors, based on TIGRE, a MATLAB-GPU toolbox for CBCT image reconstruction [33]. Three simulation data sets under three different signal to noise ratios (SNRs) were generated by running different particle histories in EGS simulation for all projection angles. Table 1 shows the simulation and reconstruction parameters. For all Lp-norm (0 < p ≤ 2) priors, iteration numbers were set to 300, which were enough to get the convergence.
frequency comparisons on resolution and noise performance of ASDPOCS algorithms using different Lp-norm (0 < p ≤ 2) priors are not sufficient at the present time [10,14,17]. Some papers were based on global evaluation metrics such as signal to noise ratio (SNR), contrast to noise ratio (CNR), and root of mean square error (RMSE) [18–20], and although there have been reports evaluating the spatial resolutions using modulation transfer function (MTF) analysis, they presented results within limited number of Lp-norm priors [21,22]. Moreover, relatively few works were conducted on quantitatively comparing both the resolution and noise performance by spatial frequency dependent metrics for the ASD-POCS algorithms with Lp-norm (0 < p ≤ 2) regularization. There, the aim of this study is to employ an overall and comprehensive performance comparison of Lp-norm (0 < p ≤ 2) priors. The modulation transfer function (MTF) and the noise power spectrum (NPS) have been widely recognized as the most relevant spatial frequency dependent metrics of resolution and noise evaluation in radiographic imaging. The MTF is a measure of image sharpness and the NPS describes the noise in the image. These two metrics have been introduced into CBCT imaging system evaluation with varying measurement methods [23,24]. MTF and NPS can be further combined to derive the noise equivalent quanta (NEQ), which is a decisive indicator of the signal and noise transfer properties of an imaging system [25]. In CT, methods for measuring the presampling MTF were developed by using an aluminum edge [26] or a tungsten wire [27]. These methods typically use highly dense materials to simulate an infinite impulse to maximize CNR and improve measurement precision. However, these materials are atypical in biomedical research and preclinical applications. Moreover, the use of high-density metallic test object for MTF measurement may introduce metal artifacts that can introduce additional uncertainties owing to their CT artifacts. According to this problem, it is important to use test objects with lower contrast to avoid metal artifacts and yield MTFs relevant to realistic one used in clinical images. A currently widely adopted low contrast phantom for MTF measurement in CT is the American College of Radiology (ACR) computed tomography accreditation program (further denoted as ACR phantom) [28]. The ACR phantom is a solid cylinder phantom constructed primarily from a water-equivalent material. Note that in 2011, low-contrast cylinder test objects similar to ACR phantom were proposed by the American National Standard for MTF measurement of CT [29]. In this standard, the circumference of the cylinder is used as an edge to measure the MTF, and the inner part of it is used to calculate the NPS. In this work, a Monte Carlo (MC) simulation platform was built employing the EGSnrc/BEAMnrc code system. The CBCT scans were simulated under the condition of different noise levels. Then the simulated and experiment data were employed to quantitatively investigate and compare the resolution and noise performance of Lp-norm (0 < p ≤ 2) regularization. Several figures of merit including the MTF, NPS and NEQ, which account for the spatial resolution and the noise performance, were introduced to objectively evaluate the reconstruction performance at different noise levels.
2.2. Evaluation metrics Frequency dependent metrics including MTF, NPS and NEQ were calculated for algorithm performance evaluation. MTF was assessed by a radial oversampling methodology [28] from the reconstructed images as illustrated in Fig. 3. The oversampled edge spread function (ESF) was determined on increments of the distance from the edge [28], and then it was fitted with the linear combination of two Gaussian functions and a Boltzmann function to avoid instability [34]. ESF was then differentiated to obtain the line spread function (LSF), from which the MTF was calculated by using the Fourier transform. To investigate the frequency variation of noise under Lp-norm (0 < p ≤ 2) regularization, the NPS was calculated based on ensemble statistics from the ACR phantom images. The NPS utilizes the Fourier transform of noise images to determine the variance of noise power present at each spatial frequency. Noise-only images were obtained by the subtraction of the averaged image from each of 20 central reconstructed slices. Total 720 region of interests (ROIs) in 20 noise-only images were averaged to calculate the 2D NPS distribution. ROIs were created, centered along a circle, as illustrated in Fig. 3 [28]. Finally, the desired NPS profile was generated by averaging 8 radial profiles from center to the edge to avoid fluctuation. The image quality of the reconstructed ACR phantom was further assessed by the measurement of the noise equivalent quanta (NEQ), which is a measure of the signal to noise ratio of an imaging system. Mathematically, the NEQ can be represented as [29]
2. Materials and methods 2.1. Virtual CT system The full MC simulation of a kV CBCT imaging system was carried out with the EGSnrc/BEAMnrc MC code system [30,31]. The component modules in the BEAMnrc simulation consisted of an 'XTUBE' (the X-ray target) and a 'SLABS' (the filter). The SLABS included a filter of 4.0 mm of Aluminum plus 0.1 mm of Beryllium. The tungsten anode target (2 mm thick) was contained in a copper holder, with 9 degree from z-axis (see Fig. 1). Cross sections for all materials in the simulation were generated from the XCOM data set using the PEGS4 code system [32]. The global electron cut-off energy (ECUT) and photon cut-off
NEQ (f ) = 73
[S·MTF (f )]2 NPS (f )
(1)
Physica Medica 52 (2018) 72–80
L. Zhang et al.
Fig. 1. X-ray tube geometry in BEAMnrc for x-z projection.
Fig. 2. Left: X-ray spectral distribution derived from BEAMDP, right: digital ACR phantom.
where f is spatial frequency. S is the signal level, measured as the average CT value of the object. To validate the conclusion obtained from the simulation, a mouse data set was then adopted for reconstruction. The mouse was fixed in an imaging cavity and was scanned with 42 views distributed evenly in 360 degrees. The study protocol was approved by the local institution review board and performed according to the guidelines of the committee on animal use and care of Chinese Academy of Sciences. The related parameters for image acquisition were shown in Table 2. The acquired projections were sheared and rebinned to 482 × 482 pixels with a pixel size of 0.0977 mm × 0.0977 mm. An example projection view of the mouse was shown in Fig. 4. The voxel size of the reconstruction was 0.123 mm × 0.123 mm × 0.123 mm for 241 × 241 × 241 voxels. The iteration numbers were also set to 300 for all Lpnorm (0 < p ≤ 2) priors, which were enough to get the convergence. To carry out a quantitative research, contrast to noise ratios (CNRs) of the reconstructions using different Lp-norm (0 < p ≤ 2) priors were calculated from the areas as indicated in Fig. 5. The red rectangle area denotes the target area and the green one represents the background area, they were selected to calculate the CNR as following [18]:
Fig. 3. Diagram of ROIs selecting in NPS calculation. Table 2 Acquisition parameters for real CT data set. Parameters
Voltage
Current
Exposure time
Pixel size
Pixel numbers
Value
45KV
0.3 mA
375 s
0.0488 mm
1024 × 1024
Table 1 Parameters setting in phantom simulation and image reconstruction. Parameters
Value
Projection parameters
Reconstruction parameters
Pixel size(mm)
Pixel numbers
Projection numbers
Voxel size(mm)
Voxel numbers
1.5 × 1.5
256 × 256
20
1.96 × 1.96 × 1.96
128 × 128
74
Physica Medica 52 (2018) 72–80
L. Zhang et al.
difference can be found within the range of p ∈ (0, 0.5]. When the value of p increases in the range of (1.5, 2], MTF changes slowly. The same trends can be seen from the MTFs obtained from the simulation data sets under different SNRs. The NPSs associated with different Lp-norm (0 < p ≤ 2) regularization at varying noise levels were given in Fig. 8. The left, middle and right column shows the NPS of reconstructed images from projections with SNR = 35 dB, 40 dB and 45 dB respectively. As shown in the figure, L1 regularization resulted in the lowest NPS among all Lp-norm (0 < p ≤ 2) regularizations under three different projection noise levels. The measured NPS was increased when the value of p deviated from 1. Fig. 9 compares the NEQ measurements on log scale at different noise levels from the reconstructions with various Lp-norm priors. Within the range of p ∈ [1, 2], L1-norm prior achieves the highest NEQ, and an obvious NEQ reduction can be found with the increment of the value of p. The NEQ curves for Lp-norm (1 < p ≤ 2) priors have a local maximum around the intermediate frequency within cut-off frequency (half the Nyquist sampling frequency). Within the range of p ∈ (0, 1), NEQ is increased when the value of p is approaching 1. No obvious difference among NEQ curves can be found within the range of p ∈ (0,0.5] in terms that the projection SNR is 35 dB or 40 dB. Generally speaking, ASD-POCS algorithm with L1-norm prior achieves the best NEQ among all results except for the most high frequency part approaching the cut-off frequency.
Fig. 4. An example projection view of the mouse.
3.2. Experiment Fig. 10 shows the mouse head cross section of a display window [0, 0.1] reconstructed with different Lp-norm (0 < p ≤ 2) priors. The oval structure in the middle is the skull of the mouse, and two circle structures in front of the skull are the front claws. The CNRs of the reconstructions using varying Lp-norm priors (p = 0.3, 0.5, 0.8, 1.0, 1.1, 1.2, 1.3, 1.5, 1.8, 2.0) were calculated as 22.13, 25.25, 57.64, 58.59, 49.94, 47.07, 46.27, 42.25, 36.90 and 32.64 respectively. The results demonstrated that the best CNR performance was obtained when p equals 1, which is in consistent with the NEQ results in the simulation studies (see Fig. 9). The CNR is a metric for the overall evaluation of contrast and noise in the image, while the NEQ provides detailed contrast and noise assessment related to the spatial frequency. In order to observe the difference of each algorithm, the image profiles along the horizontal and vertical lines, as shown in Fig. 5, are plotted in Fig. 11. The ROI is selected from 1D profiles and marked by the red rectangles in the left-hand column of Fig. 11. The profiles of these selected pixels are plotted in the right-hand column to observe the differences more clearly. It can obviously be seen that the profiles of the ASD-POCS algorithm with smaller Lp-norm priors show sharper edges than those using higher Lp-norm priors. However, the ASD-POCS algorithm using small Lp-norm (p ≤ 0.5) priors have an effect of introducing spikes at edges. These spikes behave as pepper-like impulse noise as shown in Fig. 10 (a) and (b), which may lead to difficulties in identifying small structures of the sample.
Fig. 5. Illustration to the calculation of CNR.
CNR =
It −Ib σb
(2)
where It represents for the average pixel value in the target area, Ib stands for the average pixel value in the background and σb is the square root of the background variance. 3. Results 3.1. Simulation Fig. 6 displays the central axial slices reconstructed by the ASDPOCS algorithm using varying Lp-norm priors. The simulation data set under the SNR of 40 dB was selected as a typical example. The display gray scale window of the figure was set to [0.0044, 0.0054] for better showing the edge and noise performance. As shown in the figure, ASDPOCS algorithms using small Lp-norm (0 < p < 1) priors preserved the high contrast edges of the image, while high Lp-norm (1 < p < 2) priors blurred the edges. When p ∈ [0.5, 1), the smaller the value of p, the sharper edger obtained by the Lp regularization. When p ∈ (0, 0.5], the performance of Lp-norm priors has no significant difference on edge enhancement. When p ∈ (1, 2], artifacts can be found in the reconstructed images. Obvious artifacts can be found in the reconstructed slices when the value of p deviated from 1. The subsequent calculated MTF based on the fitted ESFs are shown in Fig. 7. It can be observed that the MTF is improved when reducing the value of p. This demonstrates that better spatial resolution is achieved by the ASD-POCS algorithm using small Lp-norm (0 < p < 1) priors. Specifically speaking, the MTF improvement is significant when decreasing p from 1.5 to 0.5. No significant MTF
4. Discussion and conclusion Many researches have proved that an exact reconstruction of the sparse signal is possible with substantially fewer measurements by replacing the L1-norm with the Lp-norm (0 < p < 1) [14,35–37]. A few researchers perform the reconstruction using ASD-POCS algorithm with a L0-norm prior approximated by the Lp-norm, and compared the image reconstruction performance with those using L2-norm and L1-norm prior [14,38]. In this paper, the resolution and noise performances of ASD-POCS algorithms using different Lp-norm (0 < p ≤ 2) priors are assessed and compared by simulation and experiment. In the phantom study, the comparison via MTF showed that a 75
Physica Medica 52 (2018) 72–80
L. Zhang et al.
Fig. 6. Reconstructions of the ACR phantom from different p with SNR of 40 dB, central axial slice reconstructed using (a) p = 0.3 (b) p = 0.5 (c) p = 0.8 (d) p = 1.0 (e) p = 1.3 (f) p = 1.5 (g) p = 1.8 (h) p = 2.0. The display gray scale window is [0.0044, 0.0054].
speaking, a better NEQ in low spatial frequency is obtained when p = 1 while a higher NEQ appears at the high frequency part when p = 0.8. The resolution and noise performances of ASD-POCS algorithms was also tested along with application to an actual mouse data set. Experimental data set study shows that relative better spatial resolution and noise performance could be obtained when p is between 0.8 and 1, which is in consistent with the NEQ results in the simulation studies. What should be pointed out is that the main drawback of Lp (p > 1) regularization is not sparsity preserving [39], and as a consequence the corresponding reconstructions potentially suffer from noise and the overall performance is lower than that of p ≤ 1. This phenomenon is verified by the simulation and experiment in this study. Besides, the reconstruction time was not taken into consideration in this work. Although iterative reconstruction is usually characterized by a high computational cost, the steady improvement of the compute capabilities of GPUs can decrease the reconstruction time significantly, making them increasingly suitable for iterative reconstruction in CT [40,41].
smaller value of p resulted in a higher resolution, and when p was smaller than 0.5 or larger than 1.5, the difference in reconstructed resolution were not obvious. Thus, it could be concluded that the ASDPOCS algorithm with smaller Lp-norm prior has a higher capability to preserve edge details compared to those using bigger Lp-norm. However, the ASD-POCS algorithms using small Lp-norm (p ≤ 0.5) priors may result in inconsistent artifacts, while those with high Lpnorm (p > 1.2) priors will lead to blurred reconstructions. In our phantom study, the best noise performance was obtained around p = 1 for the simulation data sets under any SNR condition. It should be noted that as the SNR of projection views increased to 40 dB and 45 dB, no marked improvement was obtained on the MTF performance compared to the reconstruction results from projection views of 35 dB. To take both MTF and NPS into consideration, the noise-resolution tradeoff study is employed by calculating the NEQ of the ASD-POCS algorithms. Although the ASD-POCS algorithms using small Lp-norm (p ≤ 0.5) priors yield a higher MTF than do the high Lp-norms, the best NEQ performance is achieved when p is between 0.8 and 1. Specifically
Fig. 7. MTF based on fitted ESFs for the projections SNR of 35 dB (left), 40 dB (middle) and 45 dB (right).
76
Physica Medica 52 (2018) 72–80
L. Zhang et al.
Fig. 8. NPS reconstructed using different Lp-norm priors under projection SNRs of 35 dB (left), 40 dB (middle) and 45 dB (right).
Fig. 9. NEQ reconstructed with different p of different SNRs (left:35 dB, middle:40 dB, right:45 dB).
Fig. 10. Reconstructed images from real CT data set with different p: (a) p = 0.3 (b) p = 0.5 (c) p = 0.8 (d) p = 1.0 (e) p = 1.1 (f) p = 1.2 (g) p = 1.3 (h) p = 1.5. (i) p = 1.8. (j) p = 2.0. The display gray scale window is [0, 0.1].
77
Physica Medica 52 (2018) 72–80
L. Zhang et al.
Fig. 11. One-dimensional profiles along the horizontal line (the 180st row of the reconstructed images) and the vertical line (the 121st column of the reconstructed images.) (a) The image profiles along the horizontal line, (b) the partial profiles of the selected ROI, (c) the image profiles along the vertical line, (d) the partial profiles of the selected ROI.
Acknowledgments
system for early cancer detection. XCT imaging is an important part of this system. This study serves the goal of finding a high-performance algorithm for XCT reconstruction with limited data.
This work was supported by the National Natural Science Foundation of China (No. 81401453, 81371602, 61475115, 61475116, 61575140, 81571723, 81671728), and Tianjin Municipal Government of China (15JCZDJC31800, 16JCZDJC31200, 17JCZDJC32700).
Conflicts of interest
Role of the funding source
The authors have no potential conflicts of interest related with this article.
The role of the study sponsor is to develop multi-model imaging Appendix Iterative reconstruction algorithms based on Lp-norm (0 < p ≤ 2) regularization For simplicity, the discretized Beer-Lambert law is selected as the forward projection model for monochromatic x rays in a partially attenuating medium.
⎡ ⎤ Yi = di exp ⎢− ∑ hij μj ⎥ j ⎣ ⎦
(3)
where i and j are linear indices in the image vectors and projection data, respectively. di is the incident photon count. Yi is the detector photon count after attenuation by the object along the jth X-ray path. hij is a discrete length path through voxel j and μj denotes the voxel’s linear attenuation coefficient. From Eq. (3), a logarithmic operation is performed to obtain the so-called projection data along jth X-ray path
Pi = ln
di = Yi
∑
hij μj = Hi μ
(4)
j
78
Physica Medica 52 (2018) 72–80
L. Zhang et al.
A discrete linear system for modeling the CT-imaging process can be derived from Eq. (4) and written as
P = Hμ
(5)
where the vector P = {pi |i = 1,2,…,M } denote the projection data and the vector μ = {μ|i = 1,2,…,n} is the vectorized image to be reconstructed; H is the system matrix modeling the forward projection. With this discrete data model, the image reconstruction problem is equivalent to inverting the linear system Eq. (3), finding μ from a given projection data set P. However, even under ideal conditions, it is difficult to directly invert Eq. (5) when the system is under-determined, i.e., M < N, for image reconstruction from a reduced number of projections under study. Regularization is an effective way to constraint an under-determined problem. Similar to the previous studies [10], the CT image reconstruction problem can be formulated as a constrained optimization problem below
μ∗ = arg min‖μ‖pp
(6)
Subject to the inequality constraints: (a) data fidelity
|Hμ−P| ⩽ ε
(7)
and (b) non-negativity
μj ⩾ 0
(8)
where the user-defined parameter ε is interpreted as the variance of the difference between the real and the predicted of p to the Lp-norm of the gradient magnitude image, which can be written as
‖μ‖pp = (∂x μijk ) p + (∂y μijk ) p + (∂z μijk ) p
projections. ‖μ‖pp
is the power (9)
When p = 0, it’s exactly the L0-norm priors which equal to the number of nonzero elements. While minimization of the L0-norm is a combinatorial optimization problem and not practical [42]. Replacing p = 0 with p = 1 yields the L1-norm priors, it’s a convex problem and can be solved efficiently. L1 priors has been applied in image reconstruction in CT and there are plenty of impressive results [10,43]. Another choice is Lp-norm (0 < p < 1) priors, has proposed in image processing and image reconstruction since it promotes sparser results [44,45]. Reconstructions with L2norm (p = 2) priors has been proposed early, while a main drawback is the results with L2-norm priors are oversmoothing [15]. The above constrained optimization problem is extended from the ASD-POCS algorithm introduced by Sidky and Pan [9], which is modified by replacing the conventional L1-norm priors by an Lp-norm (0 < p ≤ 2) regularized version in Eq. (9). The algorithm alternatively minimizes the Lpnorm using adaptive steepest descent (ASD) and data fidelity using projection onto convex sets (POCS).
[19] Lu Y, Chan HP, Wei J, Hadjiiski LM. Selective-diffusion regularization for enhancement of microcalcifications in digital breast tomosynthesis reconstruction. Med Phys 2010;37(11):6003–14. [20] Huang J, Ma J, Liu N, Zhang H, Bian Z, Feng Y, et al. Sparse angular CT reconstruction using non-local means based iterative-correction POCS. Comput Biol Med 2011;41(4):195–205. [21] Hashemi SM, Song WY, Sahgal A, Lee Y, Huynh C, Grouza V, et al. Simultaneous deblurring and iterative reconstruction of CBCT for image guided brain radiosurgery. Phys Med Biol 2017;62(7):2521. [22] Gonzales B, Spronk D, Cheng Y, Tucker AW, Beckman M, Zhou O, et al. Rectangular fixedgantry CT prototype: combining CNT x-ray sources and accelerated compressed sensingbased reconstruction. IEEE Access 2014;2:971–81. [23] Han M, Park S, Baek J. Effect of anatomical noise on the detectability of cone beam ct images with different slice direction, slice thickness, and volume glandular fraction. Opt Express 2016;24(17):18843–59. [24] Gao Y, Bian Z, Huang J, Zhang Y, Niu S, Feng Q, et al. Low-dose X-ray computed tomography image reconstruction with a combined low-mAs and sparse-view protocol. Opt Express 2014;22(12):15190–210. [25] Tward DJ, Siewerdsen JH. Cascaded systems analysis of the 3D noise transfer characteristics of flat-panel cone-beam CT. Med Phys 2008;35(12):5510–29. [26] Boone JM. Determination of the presampled MTF in computed tomography. Med Phys 2001;28(3):356. [27] Richard S, Husarik DB, Yadava G, Merphy SN, Samei E. Towards task-based assessment of CT performance: system and object MTF across different reconstruction algorithms. Med Phys 2012;39(7Part1):4115–22. [28] Friedman SN, Fung GSK, Siewerdsen JH, Tsui BMW. A simple approach to measure computed tomography (CT) modulation transfer function (MTF) and noise-power spectrum (NPS) using the American College of Radiology (ACR) accreditation phantom. Med Phys 2013;40(5). [29] American national standard for evaluating the image quality of x-ray computed tomography (CT) security-screening systems, IEEE Standard ANSI N42.45, 2011. [30] Rogers DWO, Faddegon BA, Ding GX, Ma CM, We J, Mackie TR. BEAM: a Monte Carlo code to simulate radiotherapy treatment units. Med Phys 1995;22(5):503–24. [31] BEAMnrc users manual, National Research Council of Canada, Ottawa, ON, 2011. [32] The EGS4 Code System, Stanford Linear Accelerator Center, Stanford, CA, 1985. [33] Biguri A, Dosanjh M, Hancock S, Soleisani M. TIGRE: a MATLAB-GPU toolbox for CBCT image reconstruction. 2016;2(5):055010. [34] Brunner CC, Kyprianou IS. Material-specific transfer function model and SNR in CT. Phys Med Biol 2013;58(20):7447. [35] Chartrand R. Exact reconstruction of sparse signals via nonconvex minimization. IEEE Signal Process Lett 2007;14(10):707–10. [36] Chartrand R, Staneva V. Restricted isometry properties and nonconvex compressive sensing. Inverse Prob 2008;24(3):035020. [37] Sidky EY, Chartrand R, Boone JM, Pan X. Constrained TpV minimization for enhanced exploitation of gradient sparsity: application to CT image reconstruction. IEEE J Transl Eng Health Med 2014;2:1–18.
References [1] Qi H, Chen Z, Wu S, Xu Y, Zhou L. Iterative image reconstruction using modified non-local means filtering for limited-angle computed tomography. Phys Med 2016;32(9):1041–51. [2] Chen Z, Jin X, Li L, Wang G. A limited-angle CT reconstruction method based on anisotropic TV minimization. Phys Med Biol 2013;58(7):2119–41. [3] Protonotarios NE, Fokas AS, Gaitanis A, Kastis GA. Inversion of the attenuated radon transform using cubic splines. Phys Med 2016;32(9):242–3. [4] Yan B, Zhang W, Li L, Zhang H, Wang L. Quantitative study on exact reconstruction sampling condition by verifying solution uniqueness in limited-view CT. Phys Med 2016;32(1):1321–30. [5] Huang J, Ma J, Liu N, Zhang H, Bian Z, Feng Y, et al. Sparse angular CT reconstruction using non-local means based iterative-correction POCS. Comput Biol Med 2011;41(4):195–205. [6] Chen Z, Qi H, Wu S, Xu Y, Zhou L. Few-view CT reconstruction via a novel non-local means algorithm. Phys Med 2016;32(10):1276–83. [7] Slavine NV, Guild J, McColl RW, Anderson JA, Oz OK, Lenkinski RE. An iterative deconvolution algorithm for image recovery in clinical CT: a phantom study. Phys Med 2015;31(8):903–11. [8] Egiazarian K, Foi A, Katkovnik V. Compressed sensing image reconstruction via recursive spatially adaptive filtering. IEEE Int Conf Image Process 2007;1:I-549–52. [9] Mueller K, Yagel R, Wheller JJ. Anti-aliased three-dimensional cone-beam reconstruction of low-contrast objects with algebraic methods. IEEE Trans Med Imaging 2002;18(6):519–37. [10] Sidky EY, Pan X. Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization. Phys Med Biol 2008;53(17):4777. [11] Han X, Bian J, Eaker DR, Kline TL, Sidky EY, Ritman EL, et al. Algorithm-enabled lowdose micro-CT imaging. IEEE Trans Med Imaging 2011;30(3):606. [12] Niu T, Zhu L. Accelerated barrier optimization compressed sensing (ABOCS) reconstruction for cone-beam CT: phantom studies. IEEE; 2012. p. 2354–7. [13] Lee H, Xing L, Davidi R, Li R, Qian J, Lee R, et al. Improved compressed sensing-based cone-beam CT reconstruction using adaptive prior image constraints. Phys Med Biol 2012;57(8):2287–307. [14] Sidky EY, Pan X, Reiser IS, Reiser IS, Niskikawaet RM, et al. Enhanced imaging of microcalcifications in digital breast tomosynthesis through improved image-reconstruction algorithms. Med Phys 2009;36(11):4920. [15] Morigi S, Reichel L, Sgallari F. Fractional Tikhonov regularization with a nonlinear penalty term. J Comput Appl Math 2017;324:142–54. [16] Fu SJ, Zhang CM, Tai XC. Image denoising and deblurring: non-convex regularization, inverse diffusion and shock filter. Sci China Inf Sci 2011;54(6):1184–98. [17] Sidky EY, Duchin Y, Pan X, Ullberg C. A constrained, total-variation minimization algorithm for low-intensity x-ray CT. Med Phys 2011;38(S1). [18] Shieh CC, Kipritidis J, O’Brien RT, Kuncic Z, Keall PJ. Image quality in thoracic 4D conebeam CT: a sensitivity analysis of respiratory signal, binning method, reconstruction algorithm, and projection angular spacing. Med Phys 2014;41(4).
79
Physica Medica 52 (2018) 72–80
L. Zhang et al.
[42] Sidky EY, Chartrand R, Pan X. Image reconstruction from few views by non-convex optimization. 5. IEEE; 2007. p. 3526–30. [43] Chartrand R, Yin W. Iteratively reweighted algorithms for compressive sensing. IEEE; 2008. p. 3869–72. [44] Okawa S, Hoshi Y, Yamada Y. Improvement of image quality of time-domain diffuse optical tomography with lp sparsity regularization. Biomed Opt Express 2011;2(12):3334–48. [45] Chartrand R, Sidky EY, Pan X. Nonconvex compressive sensing for X-ray CT: an algorithm comparison. IEEE; 2013. p. 665–9.
[38] Hu Y, Xie L, Luo L, Nunes JC, Toumoulin C. L0 constrained sparse reconstruction for multi-slice helical CT reconstruction. Phys Med Biol 2011;56(4):1173. [39] Theodoridis S, Kopsinis Y, Slavakis K. Sparsity-aware learning and compressed sensing: an overview. arXiv preprint arXiv. 2012;1211;5231. [40] Jia X, Lou Y, Li R, Song WY, Jiang SB. GPU-based fast cone beam CT reconstruction from undersampled and noisy projection data via total variation. Med Phys 2010;37(4):1757–60. [41] Després P, Jia X. A review of GPU-based medical image reconstruction. Phys Med: Eur J Med Phys 2017;42:76–92.
80