Goossens, Bart

Pizurica, Aleksandra

Staelens, Steven

Van Holen, Roel

Vandeghinste, Bert

Vandenberghe, Stefaan

Vanhove, Christian

English

Ghent University Academic Bibliography

Combined shearlet and TV regularization insparse-view CT reconstructionBert Vandeghinste, Bart Goossens, Roel Van Holen, Christian Vanhove, Aleksandra Pizˇurica,Stefaan Vandenberghe and Steven StaelensAbstract—Preclinical in vivo micro computerized tomographysuffers from high image noise, due to limitations on total scanningtime and the small pixel sizes. A lot of different noise mini-mization algorithms have already been proposed to reconstructimages acquired in low dose settings. Sparse-view reconstructionamongst others can reduce acquisition dose significantly, byacquiring only a small subset of projection views. Total Variationminimization has been used extensively to solve these problems.However, the performance of TV is suboptimal for compleximages, compared to simple images with little texture. This ismainly due to the underlying piecewise constant image modelimposed by TV.A recent efficient solver was developed for convex problems,able to incorporate regularization terms different from TV.The work presented here is a proof-of-concept study combiningboth TV as well as shearlets as regularization terms into onegeneral CT reconstruction algorithm. Shearlets, closely relatedto wavelets, take edges into account in a multitude of directionsat different scales, and have good compaction properties. Thismakes shearlets a better candidate than TV for compressed sens-ing problems. The resulting reconstructions were compared toTV minimization and to shearlet minimization. The combinationof both shows benefits for sparse-view CT imaging, and leadsto edge-preserved image denoising. Difference images show avery small loss in resolution, which may be caused by difficultparameter selection.Index Terms—Computed Tomography, Iterative Algorithms,Noise, Reconstruction AlgorithmsI. INTRODUCTIONPreclinical in vivo micro computerized tomography (µCT)suffers from high image noise, as a result of the small detectorpixel sizes, a high scatter-to-primary noise ratio [1], and thelimited time animals can be safely kept under anesthetics.This results in limited soft tissue contrast. Reducing thedose without sacrificing image quality could offer significantbenefits for longitudinal preclinical research, where the smallanimals receive a large dose within a timeframe of severaldays to weeks.This work was supported in part by a PhD grant to Bert Vandeghinste of theInstitute for the Promotion of Innovation through Science and Technology inFlanders (IWT-Vlaanderen). B. Goossens and R. Van Holen are postdoctoralresearch fellows with FWO, Flanders. C. Vanhove is supported by theGROUP-ID consortium. Asterisk indicates corresponding author.*B. Vandeghinste, R. Van Holen, C. Vanhove, S. Vandenberghe and S.Staelens are with the Medical Image and Signal Processing (MEDISIP)research group, Ghent University–IBBT, 9000 Gent, Belgium. (e-mail:bert.vandeghinste@ugent.be).B. Goossens and A. Pizˇurica are with the Image Processing and Interpre-tation (IPI) research group, Ghent University–IBBT, Sint-Pietersnieuwstraat41, 9000 Gent, Belgium.S. Staelens is also with the Molecular Imaging Centre Antwerp, Universityof Antwerp, 2650 Edegem, Belgium.Total variation (TV) minimization has been extensivelyinvestigated in the last decade for image denoising in generaland sparse-view reconstruction in particular [2]–[5]. Thesemethods have been shown to have superior denoising perfor-mance in simple classes of images. However, TV minimizationproduces cartoon-like approximations due to its underlyingimage model. This model biases the results towards distortedimages, which may be less suitable for medical images usedin diagnostics [3].In the search for objective functions different from TV, a re-cent efficient solver was developed, based on a split-Bregmanapproach. With this solver, other regularization terms caneasily be tested. One possibility is the shearlet [6], [7], closelyrelated to wavelets, which has better directional sensitivity,better `1-norm sparsity and, because of a different underlyingimage model, does not lead to the staircasing effect. Wehave previously shown that shearlet-regularized reconstructedimages show no staircasing and exhibit small aliasing artifacts.However, these reconstructions did not outperform TV-basedregularization for all datasets [8]. In the research presentedhere, we investigated if combining TV and shearlets cancombine their benefits and reduce the small artifacts inducedby using shearlets in CT reconstruction. This combination hasalready been shown effective in Magnetic Resonance Imagingreconstruction [9]–[11].The remainder of this paper is organized as follows. In Sec.II, we introduce the problem formulation and quickly reiteratethe mathematical background for using the split-Bregman ap-proach in CT. Section III describes the set-up for the evaluationon simulated and measured preclinical in vivo µCT data. InSec. IV we compare the combination of shearlets and TVto conventional algebraic reconstruction and to reconstructionwith only one of these terms, for sparse-view data. Theseresults are then discussed in Sec. V, where suggestions arealso made for further research. Our conclusions are in Sec.VI.II. PRELIMINARIESPreviously, we have developed the split-Bregman frameworkfor regularized CT reconstruction [8], [12], [13]. We denotean `1-norm by |.|1 and and `2-norm by ||.||2. The followingminimization problem is solved:xˆ = argminxE(x) + λ||C−1/2(y −Wx)||22, (1)with x the unknown reconstructed image, E(x) the penaltyterm, λ the Lagrangian multiplier, a constant which determinesthe amount of data-fitting and amount of regularization, C aprewhitening term, y the measured data and W the systemmatrix.The penalty term E(x) can include different regularizers.Previously, the `1-norm of the discrete gradient operator hasalready been used in anisotropic fasion [13], as well as a casewhere shearlets were used [8]. In this study, we will use apenalty term combining both:E(x) = γTV (|∇xx|1 + |∇yx|1) + γSH |Sx|1, (2)with ∇ the discrete gradient operator, S the shearlet transformand constants γTV and γSH weighting the influence of the twocomponents.This cost function can be minimized by using Bregmaniterations and variable splitting [12], effectively splitting the`1- and `2-norm into subproblems which are easier to solve[13]. This results in equations (3a) - (3g).III. MATERIALS AND METHODSA. Simulated dataFan-beam data was simulated using MC-GPU v1.2 [14]1, aGPU-accelerated x-ray transport simulator. A high resolutionphantom was built based on the work of the FORBILD group2(Fig. 1a). The resolution rods were simulated as containing air,surrounded by soft tissue. One 3.6864× 3.6864× 0.1395 cmthick slice was generated, containing 8192 × 8192 voxels, toget sufficient subsampling in the holes with smallest diameter.The detector was simulated as a perfect detector with 100%efficiency, and consists of 296 elements with a pixel pitchof 0.14 mm, acquiring 360 uniformly spaced projection viewsover 2pi. All data was generated using a 60 keV monoenergeticx-ray source with 108 photons per ray. The sinogram includesscattered photons.B. Measured dataThe X-O CT system (Gamma Medica Ideas, Northridge,California, USA) was used to obtain preclinical data of onein vivo contrast-enhanced mouse study. This flat-panel cone-beam system consists of a 1280×1120 detector with a 100µm pixel pitch. The tube current is determined automaticallyduring calibration to ensure that the dynamic range of thedetector is optimally used. Fan-beam data were generated byretaining only the central detector row. 2048 projection viewswere obtained over 2pi. A new dataset was generated fromthis projection data, by removing all but every 16th projection(128 views).C. Data reconstruction and analysisAll datasets were reconstructed using 4 methods: SIRT,split-Bregman using anisotropic TV (SpBR-TV), split-Bregman using shearlets (SpBR-SH) and split-Bregman usingboth shearlets and anisotropic TV (SpBR-SHTV). In SpBR,Eq. 3a was solved using 30 iterations of conjugate gradient1Freely available from http://code.google.com/p/mcgpu/2http://www.imp.uni-erlangen.de/phantoms/highcontrast/highcontrast.htmlresolution phantom preclinical referencea bFig. 1. Reference images for the simulated data (resolution phantomconverted to 60 keV attenuation values) and the preclinical data (convergedSIRT reconstruction of 2048 projection views).on the normal equations. Matrix W was implemented as the2D Distance Driven projector [15]. The shearlet transform Swas implemented as previously implemented by Goossens etal. [16], based on the Meyer wavelet. Equations 3b-3d wereimplemented using soft-shrinkage [12]. Equations 3e-3g aretrivial to solve. All regularized reconstructions converged andwere stopped at iteration 20. All SIRT reconstructions werestopped when ||x(i+1) − x(i)||22/||x(i+1)||22 < 10−4.The simulated data was reconstructed to a 2562-grid withvoxel pitch 0.16 mm, the preclinical data was reconstructed toa 2562-grid with 0.13 mm voxel pitch. The diagonal elementsof C (see Eq. 3a) were set to cii = e−yi [17] with yi themeasured counts, serving as an estimator for the mean numberof counts.For each regularized reconstruction, parameter λ was em-pirically chosen, by reconstructing with different λ valuesand gradually making the search interval smaller, fine tuningthe amount of denoising. This generally results in 5 to 10reconstructions needed to determine a good λ value for thecase of only one regularizer. When SH and TV were com-bined, γSH and γTV also had to be empirically determined.The parameter µ was always set to 0.585 × λ, which wasempirically determined.The peak signal-to-noise ratio (PSNR) was determined toobjectively evaluate the preclinical images. The SIRT recon-struction of 2048 projection views was used as the referencehigh dose image (Fig. 1b).IV. RESULTSFigure 2 shows the resolution phantom reconstructed withthe different methods, zoomed in to the low resolution partof the resolution phantom. some For SpBR-TV, noise patchesbecome apparent in the heavily denoised image (λ = 600),compared to less denoising (λ = 4000). However, a slightincrease in resolution can be noted when only a small amountof denoising is applied. Furthermore, the rods are shapedirregularly and are not perfectly round.In the SH reconstruction, denoising with a low λ factordoes not eliminate all noisy patches in between of the rods.The resolution increases when less denoising is used, alsoincreasing the noise in the background. However, when SHand TV are combined (γSH = 1.0 and γTV = 0.45), roundx(i+1) = argminxλ2||C−1/2(y −Wx)||22 +µ2||d(i)TV,x −∇x(x)− b(i)TV,x||22 (3a)+µ2||d(i)TV,y −∇y(x)− b(i)TV,y||22 +µ2||d(i)SH − Sx− b(i)SH ||22d(i+1)TV,x = argmind|d|1 + µ2 ||d−∇xx(i+1) − b(i)TV,x||22 (3b)d(i+1)TV,y = argmind|d|1 + µ2 ||d−∇yx(i+1) − b(i)TV,y||22 (3c)d(i+1)SH = argmind|d|1 + µ2 ||d− Sx(i+1) − b(i)SH ||22 (3d)b(i+1)TV,x = b(i)TV,x + (∇xx(i+1) − d(i+1)TV,x ) (3e)b(i+1)TV,y = b(i)TV,y + (∇yx(i+1) − d(i+1)TV,y ) (3f)b(i+1)SH = b(i)SH + (Sx(i+1) − d(i+1)SH ). (3g)SpBR-SH! 1300SpBR-SHTV! 4000! 1600SpBR-TVSIRT! 4000SpBR-TVSpBR-SH! 3000Fig. 2. Comparison of SIRT of noisy data, combined SH and TV(γSH = 1.0, γTV = 0.45), TV with low and high λ, SH with lowand high λ. All images were normalized to the same window.rods are achieved, whilst still perfectly minimizing the noisein the background.Figure 3 compares SpBR-TV, SpBR-SH and SpBR-SHTVto the reference high-dose SIRT reconstruction in the caseof measured preclinical data. Parameter λ was set to 2000for SpBR-SHTV, with γSH = 0.25 and γTV = 1.0. Plottedon the right are the absolute difference images, each timebetween the reconstruction and the reference image. SIRTshows streaking artifacts when only 128 views are used. Allregularized reconstructions lead to higher PSNR compared toSIRT.The difference image for SpBR-TV shows some resolutionloss at the edges of the animal bed and at the body contour, pri-marily at edges which are not mainly horizontally or verticallyoriented. For SH, these edges are not visible. However, thereare some streaking artifacts left, which could not be minimizedwith a different λ choice without sacrificing resolution. WhenSH and TV are combined, good denoising properties areobtained with a small amount of resolution loss. This isprimarily visible at the sternum, where the spongious bone ismore difficult to distinguish on SpBR-SHTV reconstructionsthan on SIRT or with TV alone.V. DISCUSSIONCombining TV and SH regularization into one algorithmshows benefits for sparse-view CT imaging. Previous researchhas shown that shearlets do not lead to any form of piecewise-constant behavior, but on the other hand do not tend toapproximate uniform regions as well as TV [8]. The prelim-inary results presented here show that a combination of bothregularization terms combines the benefits of both SH as wellas TV.One tricky and very sensitive area is parameter selection initerative algorithms. Next to determining λ empirically, thereis now also the difficulty of weighing the contribution of SHagainst the contribution of TV with γSH and γTV . Special carehas to be taken to not let TV overpower SH, as this will resultin piecewise constant behavior. In theory, we would like touse shearlets to minimize the cost function in general, and usea little bit of TV to minimize the introduced artifacts, such asGibbs phenomena next to jump discontinuities, or remainingstreaking artifacts. We have shown that this works sufficientlyin the case of preclinical data, although with a small loss ofresolution. This may be due to parameter selection.On the other hand, a larger γTV than γSH is needed whenphantom data is reconstructed. Shearlets can not reconstructSIRTSpBR-TVSpBR-SHSpBR-SHTV!1800!600!200034.73 dB40.15 dB40.34 dB43.08 dBref.Fig. 3. Preclinical sparsified data reconstructed from 128 projectionviews with SIRT, SpBR-TV, SpBR-SH and SpBR-SHTV (γSH = 0.25,γTV = 1.0). All images were normalized to the same window. Differenceis the absolute difference between reconstruction and reference.the uniform areas in phantoms accurately, as was demonstratedwith the resolution phantom in Fig. 2. However, these casesare not realistic when doing (pre)clinical measurements. Op-timal selection of these parameters will be subject to futureresearch, as well as task-based observer studies to determineif SpBR-SHTV has better diagnostic value than simple TVminimization.VI. CONCLUSIONWe have combined TV and shearlet minimization into onereconstruction algorithm, and have shown its benefits forsparse-view CT imaging in a proof-of-concept study. A smallloss of resolution is apparent, probably due to suboptimalparameter selection.REFERENCES[1] D. Gutierrez and H. Zaidi, “Assessment of scatter for the micro-CTsubsystem of the trimodality FLEX TriumphTM preclinical scanner,”Medical Physics, vol. 38, no. 7, p. 4154, 2011.[2] L. I. Rudin, O. Stanley, and E. Fatemi, “Nonlinear total variation basednoise removal algorithms,” Physica D: Nonlinear Phenomena, pp. 259–268, 1992.[3] G. T. Herman and R. Davidi, “Image reconstruction from a small numberof projections,” Inverse Problems, vol. 24, no. 4, 2008.[4] M. Lustig, D. Donoho, and J. M. Pauly, “Sparse MRI: The applicationof compressed sensing for rapid MR imaging,” Magnetic Resonance inMedicine, vol. 58, no. 6, pp. 1182–1195, 2007.[5] M. Defrise, C. Vanhove, and X. Liu, “An algorithm for total variationregularization in high-dimensional linear problems,” Inverse Problems,vol. 27, no. 6, p. 065002, May 2011.[6] G. Kutyniok and D. Labate, “Resolution of the Wavefront Set UsingContinuous Shearlets,” Transactions of the American MathematicalSociety, vol. 361, no. 5, pp. 2719–2754, 2009.[7] K. Guo, D. Labate, and W.-Q. Lim, “Edge analysis and identificationusing the continuous shearlet transform,” Applied and ComputationalHarmonic Analysis, vol. 27, no. 1, pp. 24–46, 2009.[8] B. Vandeghinste, B. Goossens, R. Van Holen, C. Vanhove, A. Pizˇurica,S. Vandenberghe, and S. Staelens, “Iterative CT reconstruction usingshearlet-based regularization,” Proc. of SPIE, pp. 1–7, 2012.[9] A. Pizˇurica, J. Aelterman, F. Bai, S. Vanloocke, H. Quang Luong,B. Goossens, and W. Philips, “On structured sparsity and selected ap-plications in tomographic imaging,” Proc. of SPIE Conference Waveletsand Sparsity XIV, p. 81381D, 2011.[10] L. He, T.-C. Chang, S. Osher, T. Fang, and P. Speier, “MR imagereconstruction from undersampled data by using the iterative refinementprocedure,” PAMM, vol. 7, no. 1, pp. 1 011 207–1 011 208, Dec. 2007.[11] J. Huang, S. Zhang, H. Li, and D. Metaxas, “Composite splittingalgorithms for convex optimization,” COMPUTER VISION AND IMAGEUNDERSTANDING, pp. 1–13, Sep. 2011.[12] T. Goldstein and S. Osher, “The split Bregman method for L1 regularizedproblems,” SIAM Journal on Imaging Sciences, vol. 2, no. 2, pp. 323–343, 2009.[13] B. Vandeghinste, B. Goossens, J. De Beenhouwer, A. Pizˇurica,W. Philips, S. Vandenberghe, and S. Staelens, “Split-Bregman-basedsparse-view CT reconstruction,” Proceedings of the 11th InternationalMeeting on Fully 3D Image Reconstruction, pp. 431–434, 2011.[14] A. Badal and A. Badano, “Accelerating Monte Carlo simulations ofphoton transport in a voxelized geometry using a massively parallelgraphics processing unit,” Medical Physics, vol. 36, no. 11, p. 4878,2009.[15] B. de Man and S. Basu, “Distance-driven projection and backprojectionin three dimensions.” Physics in Medicine and Biology, vol. 49, no. 11,pp. 2463–2475, 2004.[16] B. Goossens, J. Aelterman, H. Luong, A. Pizˇurica, and W. Philips,“Efficient design of a low redundant discrete shearlet transform,” Proc.of the 2009 Int. Workshop on Local and Non-Local approximation inImage Processing, 2009.[17] S. Ramani and J. Fessler, “A Splitting-Based Iterative Algorithm forAccelerated Statistical X-Ray CT Reconstruction.” IEEE Transactionson Medical Imaging, Nov. 2011.

Combined shearlet and TV regularization in sparse-view CT reconstruction

https://biblio.ugent.be/publication/3072110/file/4381380.pdf

Combined shearlet and TV regularization in sparse-view CT reconstruction

Abstract

Similar works

Full text

Available Versions

Ghent University Academic Bibliography