Efficient Low Dose X-ray CT Reconstruction through Sparsity-Based MAP Modeling

1 Ef ficient Lo w Dose X-ray CT Reconstruction through Sparsity-Based MAP Modeling SayedMas oud Hashemi, Student Me mber , IEEE , Soosan Beheshti, Se nior Member , IEEE , Patrick R. Gill, Narinder S. Paul, Richard S.C. Cobbold, Life Membe r , IEEE Abstract —Ultra low radiation dose in X-ray Comput ed T o- mography (CT) is an important cli nical objective in order to minimize the risk of carcinogenesis. Compr essed Sensing (CS ) enables significant reductions in radiation dose to be achiev ed by producing diagnostic images from a limited number of CT projec- tions. Howev er , the excessiv e com putation t ime that conv entional CS-based CT reconstruction typically requires has limited clinical implementation. In this paper , we first demonstrate that a thorough analysis of CT reconstruction through a Maximum a Posteriori objective fun ction results in a weighted compressiv e sensing p roblem. This analysis enables us to formulate a low dose fan beam and helical cone beam CT reconstruction. Subsequently , we provide an efficient solution to the for mulated CS pr oblem based on a Fast Composite Splitting Algorithm-Latent Expected Maximization (FCS A-LEM) algorithm. In th e proposed method we use pseudo polar Fourier transf orm as the measurement matrix in order to decrease the computational complexity; and rebinning of the projections to p arallel rays in order to extend its application to fan beam and helical cone beam sca ns. The weight inv olved in the proposed weighted CS model, denoted by Error Adaptation W eight (EA W), is calculated based on the statistical characteristics of CT reconstruction and is a function of P oisson measurement noise and r ebinnin g interpolation error . Simulation results show that low computational complexity of the proposed method made t he fast re cov ery of the CT images possible and using EA W reduces the reconstruction error by one order of m agnitude. Recov ery of a high q uality 51 2 × 512 imag e was achiev ed in less than 2 0 sec on a desktop computer without numerical optimizations. Index T erms —Computed T omog raphy , Direct Fourier Recon- struction, Pseud o-Polar F ourier T ransform , Compre ssed S ensing, Statistical Iterative CT reconstruction I . I N T R O D U C T I O N The clinical use of Computed T omog raphy (CT) has dra- matically increased over the last two decades. This is primar ily due to its u nsurpassed speed and the fine d etails that can be obtained in cro ss-sectional views o f soft tissues and organs. Compared to co n ventional radiog raphy , CT results in a rela- ti vely large rad iation dose to patien ts. Studies o ver the pa st decade ha ve sho wn that the h igher r adiation d ose is of s erious long-ter m concern in its potential for inc reasing the risk of developing cancer [1], [2]. As a result, low dose CT ima ging Masoud Hashemi (Email : sayedmasoud.ha shemiamroaba di@mail.utoronto.ca), and Ric hard S.C Co bbold (co bbold@ec f.utoronto.ca ) are with the Institut e of Biomate rials and Biomedical Engineeri ng (IBBME), Unive rsity of T oronto (UofT), T oronto, Canada. Soosan Be heshti (Email: soosan@ee.ryerson.ca ) is with the Depa rtment of Elect rical and Compute r Engine ering, Ry erson Uni versi ty , T oronto, Can ada. Pa trick. R. Gill is with Ramb us Inc ., Sunn yvale, CA, USA. Narin der S. Paul (Narinder .Pau l@uhn.ca) i s with Uni versity Health Netwo rk, joint Department of Medical Imaging and the Institute of Biomate rials an d Biomedica l Engi neering, Univ ersity of T oronto, T oronto, Canada . that maintain s the reso lution and achie ves goo d contra st to noise ratio has been the goal of many CT developments over the p ast d ecade. Low dose CT images reconstructed with con ventional Filtered Back Pro jection (FBP), which directly calculates the image in a sin gle rec onstruction step, suffer from low contrast to noise ratios. A reduced radiation do se decreases either the number of emitted photo ns or their energy . This incr eases the amount of photon noise in CT images and degrades the image quality . Several methods have been prop osed for lowering th e relative a mount of noise in a lo w d ose CT scan. These m ethods can be cate gorized into the follo wing three different approaches: 1) improving scan proto cols [3 ], 2) add ing denoising algor ithms [4], [5], and 3) investigating new reconstructio n methods [ 6], [7], [8], [9], [ 10], [1 1]. Th e first appro ach is h ardware based. Ne w iterative reconstruction approa ches are proposed by com bining the go als of the seco nd and third appro aches. Th ese methods aim to improve the reconstruc tion quality and to dec rease image artifacts. Ite rati ve reconstruc tion m ethods can be categorized into two group s: Algebraic Reconstruction T echn iques (AR T) [ 7], [12], [13] and Statistical Iter ati ve Reconstruction (SIR) [6]. Wh ile SIR methods ar e m ore succe ssful in noisy (low do se) reconstruc- tions, AR T b ased method s ha ve advantages in dealing with incomplete data. Howe ver , compared to FBP bo th method s are computatio nally expensi ve enoug h to hinder their wid espread clinical a doption. Iterative reconstruction methods have pro gressed with the introdu ction of Compr essed Sensing (CS). CS is a relatively recent innov ation in signal processing that allows recovery of imag es from f ewer p rojections than that required by the Nyquist sampling theore m [14], [15]. The overall X-ray ex- posure in CT scanners is the pro duct of th e X-ra y exposure at each projection vie w and the number of projection vie ws. While con ventio nal iterative CT reconstru ction meth ods focus on red ucing th e amoun t of X-ray expo sure in the pr ojections, CS permits reconstructions from fe wer X-ray projection vie ws. Such metho ds are cap able of reconstruc ting high quality images from app roximately on e ten th of th e numb er of views needed in FBP [8], p ermitting a much lo we r dose scanning protoco l than than that n eeded in co n ventional reconstru ction methods. Ho wev er , CS-based recon struction/tomo graphy al- gorithms suffer fro m two drawback s: they are pr ohibitively computatio nally intensive for clinical use [10], [ 11], and they have not in corpor ated CT statistics and geo metries in pr oblem formu lation [16], [17], [18], [19]. As a result, it seems un likely that these m ethods could be used d irectly for the clinical CT 2 systems. CS p rescribes solving optimization pro blems su ch as those giv en by : ˆ x = argmin x 1 2 k y − Ax k 2 2 + λ R TV ( f ) (1) ˆ x = argmin x 1 2 k y − A x k 2 2 + λ R k W T x k 1 (2) where λ R acts as a regular ization parameter specifying a trad e- off between the image prior model and th e fid elity to o bser- vations, A is the measuremen t matrix, x is the colum n vector representatio n of the desired image ( f ), y is the measur ed data, W T is a sparsifying transform, k x k 2 2 = ∑ i | x i | 2 , k x k 1 = ∑ i | x i | , and TV ( f ) = ∑ i q ( ∇ x f ) 2 i + ( ∇ y f ) 2 i where ∇ x and ∇ y are th e first deriv ati ves in direction x and y acco rdingly . The main challenge in so lving this o ptimization problem within a reasonable amou nt of time is due to the size of the measuremen t m atrix A . Curren tly , in most available CS-based reconstruc tion m ethods used f or modern CT geo metries the measuremen t matrix A is a Radon sampling matrix which models the ra ys going thr ough the patient. For example, to reconstruc t a 512 × 512 pixel image fr om 900 sensors and 1200 projection ang les, A would be a 10 80000 × 26214 4 matr ix. As typical iterations each usually require two multiplica tions by A and A T , it takes several hours of co mputation on ty pical desktop c omputers to recon struct a 5 12 × 512 image with such methods [ 10], [11]. T o r educe the co mputation al comp lexity of the image re- construction , Fourie r based recon struction metho ds have been propo sed [20], [2 1], [22]. The Central Slice Th eorem (CST) or Direc t Fourier Reco nstruction (DFR) relates th e 1D Fourier transform of the projections to the 2D F ourier transform of the image. DFR reconstructions c omprise: 1) interp olation of polar data to a Cartesian grid and 2 ) ca lculation o f the in verse FFT on th at grid to reconstruct the CT im age. Mor eover , to achieve an accepta ble re construction quality , the interpolatio n step needs o versampling, which requires add itional radiation exposure. T o add ress the interpola tion pro blem in DFR b ased methods, Equally Sloped T omog raphy (EST) has be en p roposed [23], [24], [2 5]. EST is an iterative metho d u sing the Pseud o Polar Fourier transform (PPFT) [ 26]. The PPFT has three impor tant proper ties which makes it a good alternative to conventional DFR m ethods: 1) it is c loser to polar ( equiangu lar line ) grids compare d to Cartesian grids, 2) it c an be computed with a fast algorithm [26], and 3 ) unlike interpolating the polar data on a Cartesian grid in regular DFR methods it has an analytical conjuga te fun ction. Note th at all the above mentioned meth ods such as CST , DFR, and EST assume pa rallel beam geometr y and do not take account of the fan and cone beam geometries used in mo st current CT systems. Consequently , in order to use the m the projected rays sho uld first be transformed to parallel bea ms. This step , called rebin ning, includ es in terpolation tha t induces additional err or to the reconstructed image. This p roblem has received slight attention , altho ugh it has been ad dressed in the following two paper s. A n EST based method w as propo sed to rec onstruct fan beam an d helical con e beam images in [27]. In this method, to overcome the rebinn ing i nterpo la- tion pro blem, at each iteratio n a no n-local total variation minimization smoothing step is used. An ℓ 2 -TV optimiza tion scheme was used to reco nstruct the CT imag es fro m fan beam projection s in [ 28]. T o compensate the interpo lation error, a confidenc e matr ix is add ed to the CS scheme, which co ntrols the p ropagatio n o f the erro r in the iteratio ns. It should be n oted that wh ile the ge ometry can b e incorpo rated into CS-based recon struction by rebinnin g, the statistics of th e noise that typica lly occur s in CT data has not been u tilized in formu lating the prob lem, such as tho se given by (1) and ( 2). A modified CS formulatio n, called reweighted ℓ 1 minimization , was proposed in [ 29] where it w as discussed that using approp riate weights the q uality o f the recovered signal can be improved. Using the weig hts intr oduced in this modified formu lation, some statistical priors can be added to the model. For example, two weighting m ethods wer e propo sed in [30] and [31] for recovery of the signals with a partial known support and with a priori info rmation about the pro bability of each entry of the signal being n on-zero . In this pap er , we rigo urously explore the statistical ch arac- teristics of CT image r econstruction to model it based o n Maximum a Posteriori (MAP) formula tion. It is sh own that the MAP formulation is transformed into a weighted CS problem . The weights are direct co nsequence s of the geometr y and statistics o f CT itself. The r esulted weigh t, denoted by Error Adap tation W e ight (E A W), is a fu nction of Poisson measuremen t n oise and the interpolation error caused by rebinnin g. T he first pa rt of this paper l eads to a pro posed weighted CS problem , which is solved by the method proposed in th e secon d part. T o provide an efficient solu tion, we first break the optimization problem into two simpler ℓ 2 − ℓ 1 and ℓ 2 -TV problem s by using Fast Comp osite Splitting Alg orithm (FCSA) [32]. Next we so lve each optim ization pro blem with a latent Expectation-Ma ximization ( LEM) method. The overall solution, denoted by FCSA-LEM, is ab le to recon struct high quality images from fe wer projections and consequently lower - dose CT scans while using su bstantially less co mputation load than conv entional methods. The p aper is organized as follows . In section I I th e CT reconstruc tion pro blem is form ulated and a MAP model is introdu ced for CT im ages. In section III we discuss the pro ce- dure o f how to transfor m the regularize d CT inv erse p roblem into one that c an b e solved quickly and with few inter polation artifacts. The p roposed image r econstruction algorithm ba sed on the EM estimator is provid ed in sectio n IV, and section V contains the simulation results. I I . P RO B L E M F O R M U L A T I O N In this section we describe a CT reco nstruction method based on th e op timization of the maxim um a posterior i of the projection data an d g i ven p riors. W e use two classes of prior: sparsity of the wavelet coefficients an d piecewise linearity o f the image s, and will sho w that the proposed MAP model is similar to weighted CS model. The k ey innov ation of ou r method is the introduction of weights applied to th e data, which depen d on the ma gnitude o f the noise sou rces from both m easurement a nd in terpolation. 3 A. Maximum a P osteriori Model (MAP) of CT X-ray p rojections of the parallel beam CT can be expressed as the Radon transfo rm of the object. The Rado n tran sform is defined a s [33]: g ( l , ϕ ) = R ( f ) = Z ∞ − ∞ Z ∞ − ∞ f ( u , v ) δ ( u cos ϕ + v sin ϕ − l ) d u d v (3) which is th e integral alon g a ray at angle ϕ and at the distance l from th e o rigin, δ ( u , v ) is Dirac delta function, an d f ( u , v ) is th e image attenuation at ( u , v ) . Howe ver, this is not wha t the scanners d irectly m easure. Detectors of the scanner measure the n umber of photons which hit the d etector in different angles, λ ( l , ϕ ) , whic h is usually modeled by Poisson distribution with expected v alue of λ ( l , ϕ ) [34]. The relation between the p rojections, g ( l , ϕ ) , and λ ( l , ϕ ) is g i ven b y: λ ( l , ϕ ) = λ T exp ( − g ( l , ϕ )) (4) in w hich λ T is the nu mber of radiated photo ns fro m the X-ray source. T his leads to: g ( l , ϕ ) = − log ( λ ( l , ϕ ) λ T ) (5) λ ( l , ϕ ) is usually corrup ted with two k inds o f noise: electrical noise of the detectors (with variance of σ 2 n ) an d th e photon counting noise (observed coun ts are drawn from a Poisson distribution of mean ¯ λ ). If we consider the discrete fo rmulation in which y deno tes the vectorized g ( l , ϕ ) , x den otes the vector - ized f ( x , y ) , and A is the projectio n matrix, using the secon d order T aylor series expansion of the Poisson distribution and log likelihood of the m easurements, w e have [35], [3 6]: log p ( y | x ) ≈ − 1 2 ( y − Ax ) T D ( y − Ax ) + O ( y 3 ) (6) in which O ( y 3 ) is a fun ction wh ich depends upon measured data only a nd D is a diagona l matrix. The i t h diagona l element of D is den oted by d i . Ign oring O ( y 3 ) , (6) describes a simplified CT model: y = Ax + n (7) in which n is a Gaussian distributed n oise with a cov ariance matrix D − 1 and d i is pr oportion al to detector co unts wh ich are the maximu m likelihood of the in verse of the variance of th e projection measurem ents, 1 / σ 2 y . From (5) the relation for the i t h measured projection y i is: y i = log ( λ T λ i ) = lo g ( λ T λ i ) + log ( λ i λ i ) ≈ y i + ( 1 − λ i λ i ) (8) in which y i is n oiseless an d λ i follows the Po isson distribution with σ 2 λ = λ i . As a result the v ariance of proje ction data can be e stimated from: σ 2 y i ≈ ( σ 2 λ + σ 2 n )( 1 λ i ) 2 (9) Using λ i as an unbiased estimation of λ i the diago nal elements of D can b e expressed as: d i = 1 σ 2 y i = λ 2 i σ 2 n + λ i (10) Therefo re, the MAP estimator ca n be u sed to reco nstruct the image from the p rojections, which uses the follo wing equation: ˆ x = argmax x log p ( y | x ) + log p ( x ) (11) Here h ( x ) = log p ( x ) acts as a pen alty function , which will be used later in the p aper . I t has b een shown in many studies [37], [38] that the wa velet transfo rm of the natural and medical images, θ = W T x , can be modeled by Generalized Gaussian Distribution (GGD): p ( θ ) = p ( W T x ) = K ( s , q ) · exp ( −| θ s | q ) (12) where W T is a sparsifyin g transfor m such as the w av elet transform and W is its in verse, s , q are the p arameters of the GGD and K ( s , q ) is the normalizatio n parameter . When q = 1, the GGD is eq uiv alen t to Lapla cian distribution and when q = 2 it describe s a Gaussian d istrib ution. Using (1 2), (6) and (11) we hav e the following MAP model for CT images: ˆ x = argmin x 1 2 ( y − Ax ) T D ( y − Ax ) + λ R k W T x k q (13) T ypically , q is ch osen to be 0 < q ≤ 1, θ is a spar se rep resen- tation of the image x = W θ , an d k x k q = ∑ i | x i | q . Another prior o n p ( x ) is the piecewise linearity of the images. A p -variation distrib ution is p roposed to describe piecewise constant function s [39]. If x n ( t ) = ∑ n j = 1 x n j ψ n j ( t ) is a piece wise func tion spa nned by ψ n j ( t ) , the roof- top basis, the following class of probab ility distribution can b e u sed to describe it: p ( x n 1 , ..., x n n ) = c q , n exp ( − a n ( n + 1 ) 1 − q n + 1 ∑ j = 1 | x n j − x n j − 1 | q ) (14) where a n > 0, x n 0 = x n n + 1 = 0, c q , n is no rmalizing factor, and [ x n 1 , ..., x n n ] T is a R n -valued rando m vector . When q = 1, this yields the total variation norm . Using ( 14) with q = 1, (6) an d (11) become th e following MAP model for CT images: ˆ x = argmin x 1 2 ( y − Ax ) T D ( y − Ax ) + λ R TV ( f ) (15) As can be seen, (15) and (13) are g eneralized f orms of the CS mod els gi ven by (1) and (2 ), respectively . Consequ ently , applying CS for CT reconstru ction is equiv alent to a MAP estimation of the CT images. I I I . G E N E R A L I Z E D C S M O D E L F O R F A N B E A M A N D H E L I C A L C O N E B E A M G E O M E T R I E S T o r educe the comp utational c omplexity of the CT reco n- struction, the pseudo polar Fourier transform is used as the measuremen t matrix, A . Since pseudo polar grids ar e p laced on eq ually sloped ra dial lines, the projection rays sho uld be 4 measured or in terpolated on the equally slop ed r adial lines, as shown in Figure 1: ϕ BH = tan − 1 2 m / N , − N / 2 ≤ m < N / 2 ϕ B V = tan − 1 2 m / N + π / 2 , − N / 2 ≤ m < N / 2 ϕ = ϕ BH [ ϕ B V (16) (A) (B) Fig. 1: (A) Pseudo-Polar Grids: red lines ar e Basically Hori- zontal (BH) and black lines are Basically V ertical (BV). (B) Polar Grids (red dots) on Pseudo-Po lar Grids (gray dots). Although CT scanners typically collect data along equally spaced angles, they ha ve the flexibility to collect data along the angles of a p seudo polar grid in stead. Then , the equally d istant measured data should be interpolated to the p seudo polar girds, as shown in Figure 1-B. The resulting interpolatio n error can be limited by oversampling the Fourier data by zero-p adding the pro jections on the equally sloped radial lines. Fan be am and he lical geo metries need extra in terpolation s to estimate the measur ed p rojections o n the parallel e qually sloped radia l lines first, a pr ocess called re binning. At e ach interp olation step, the in terpolation error is tracked to be in cluded in th e EA W . Th e calcu lated weig hts an d the prepa red data are fed into the FCSA-LEM solver . Consequently , the pr oposed CT rec onstruction me thod, shown in Figure 2, can be summarized by the fo llowing two major stages: • Data p reparation and rebinn ing: fan or h elical projectio ns are mapped to p arallel eq ually sloped rad ial lines u sed in the p seudo polar Four ier tran sform. The o utput y o f this stage is the 1D Fourier transform of the ca lculated parallel rays. • Image Reconstruction : the CT im age is reconstru cted us- ing the pr oposed FCSA-LE M method . The measurem ent matrix is the fast pseudo p olar Four ier transform fun ction and the in put data is y fro m the first stage. A. Comple xity Reductio n thr ough th e use of the Pseudo P olar F ourier T ransform (PP FT) Pseudo polar grids contain two types of samples: basically horizon tal (BH) and b asically vertical (BV), as seen in Figu re Fig. 2 : Flowchart o f the r econstruction m ethod. 1. BV an d BH lines are describ ed by: B V = { ω y = π l N for − N ≤ l < N , ω x = ω y . 2 m N for − N / 2 ≤ m < N / 2 } BH = { ω x = π l N for − N ≤ l < N , ω y = ω x . 2 m N for − N / 2 < m ≤ N / 2 } (17) The Fourier transfo rm on the BV grid s can be fo und fr om: F ( ω x , ω y ) = F [ m , l ] = N − 1 ∑ i 1 = 0 ˆ f 1 [ i 1 , l ] e xp ( − i 2 π i 1 m N . l N ) (18) where ˆ f 1 [ i 1 , l ] = ∑ 2 N − 1 i 2 = 0 f Z [ i 1 , i 2 ]( − 1 ) i 2 exp ( − i 2 π i 2 l 2 N ) is the 1D Fourier tr ansform of the zero-pad ded column s o f the image ( f Z ). I n fact, the s ame equation can be written fo r BH by applying the same equation on rows rather than colum ns. As a r esult, (18) can b e interpr eted as the fractiona l F ourier transform of the 1D Fourier transform of the zero-padded columns o f the image weighted b y ( − 1 ) i 2 . T o reconstruct an N × N image from its PFFT coefficients, 4 N 2 samples are needed. A fast alg orithm is proposed in [26] to calculate the PPFT and its conjugate ; it is used in o ur pro posed algorithm as A and A T , respectively . 5 B. Rebinning Pr o cess T o be able to use the central slice theorem and direct Fourier reconstruction in fan beam and helical geometries, the projec tions should be rearrang ed to parallel rays. This redistribution of the ray s is called rebinning [40], [41]. Sinc e we u se PPFT as our measurement matrix A , all the parallel rays sh ould be p laced on eq ually sloped rad ial lines, ϕ in (1 6). 1- F an beam to Equally Sloped P arallel Beams: T w o in terpolation steps ar e ne eded fo r fan beam geome try . In the first step, pro jections are interpolated on equally sloped radial lin es, as sho wn in Figure 1-A. This step makes use o f the following relationships between fan and parallel beams: R ( γ , β ) = g ( R sin γ , β + γ ) l = R sin γ ϕ = β + γ (19) where γ , R , ϕ and β are geometry p arameters d efined in Figure 3 . R ( γ , β ) is the fan beam pro jected data and g ( R sin γ , β + γ ) is th e correspo nding r ebinned parallel ray . These radial lines are then zero-padded and the 1D F ourier transform s of the zero-padded radial lines are calculated. This is eq uiv alen t to oversampling in the Fourier do main. I n the second interpolatio n step, the oversampled radial Fourier coefficients are in terpolated on p seudo po lar grid s, sh own in Figure 1-B. Sinc e the rad ial coefficients are ov ersampled, the interp olation error in this st ep is man ageably small. The output of th is step is th e measur ed data y . (A) (B) Fig. 3 : ( A) Parallel b eam g eometry and (B) 3 rd generation Fan beam geometr y . 2- Helical Geometry to Equally S loped Beams: T o rec onstruct th e helica lly scanned objects, the scanned cone beam d ata are first con verted to fan-be am data and then the fan beam data are conv erted to parallel beams. This rebinn ing process is based on the m ethod introduced in [41], called Cone Beam Sin gle Slice ReBinning (CB-SSRB). CB-SSRB consists of the following two steps: 1) F or each sou rce p osition in the h elical trajectory , ψ , fix the z -sampling distances. 2) F or each z-slice, calc ulate the comp lete f an-beam set, from which the imag e can b e estimated. This step uses the fo llowing equation to interpo late the cone b eam scanned d ata on the fan beam po ints of in terest: p z ( ϕ , u ) ≃ √ u 2 + D 2 √ u 2 + v 2 + D 2 g ψ ( u , v ) , v = u 2 + D 2 RD ∆ z (20) where p z ( ϕ , u ) is estimated fan beam pro jection at source angle ϕ and axial p osition z , g ψ ( u , v ) is the co ne beam projections at helical position, D is the distance between the so urce and the o rigin of the d etector , and u , v , a nd R are g eometry p arameters de fined in Figure 4. Each interpolated fan b eam p rojection is weig hted by : w ( φ ss , u ) =      sin 2 ( π φ ss 2 ( 2 γ T + 2 γ ) ) if φ ss ∈ [ 0 , 2 γ T + 2 γ ] 1 if φ ss ∈ [ 2 γ T + 2 γ , π + 2 γ ] sin 2 ( π ( π + 2 γ T − φ ss ) 2 ( 2 γ T − 2 γ ) ) if φ ss ∈ [ π + 2 γ , π + 2 γ T ] where φ ss = ( π 2 + γ T )( 1 − ∆ z d ) , d = 0 . 5 P ( π / 2 + γ T ) / ( 2 π ) , P is the pitch o f the helical trajecto ry , and 2 γ T is th e maximu m fan a ngle. The parallel beams g ( l , ϕ ) are estimated from the weighted fan b eams p z ( ϕ , u ) usin g (1 9), fr om which y will be calculated by com puting the 1D Fourier transfor m of g ( l , ϕ ) s. (A) ο ο ߶ ௦௦ (B) Fig. 4: (A) Helical trajectory an d (B) th e fan beam s in parallel z-slices. C. Gener alizing the CS Model to Ad apt to Nonun iform Mea- sur e ment Noise an d I nterpolation Err or In section II-A we sho wed the MAP estimator o f CT is a form of weighted CS (1 3 and 15), in which th e weight is a function of n oise v ariance and is denoted as D in (1 0). W e assert that the effects o f n oise variance an d in terpolation error can be lump ed to gether into the form o f an Err or Ada ptation W eight (AEW) , deno ted by c : d i = 1 σ y i − → c i = 1 σ y i + e i (21) in wh ich e i is the interpo lation error . The gr eater the interpo - lation erro r , the g reater the uncertainty ab out the value of ¯ y i , so the effect of inter polation error is similar to the effect of the m easurement noise σ y i . AEW can be r e written as: c i = 1 σ y i + ε i σ y i = 1 σ y i × 1 1 + ε i (22) 6 Using th is d efinition, the generalized CS is as fo llows : ˆ x = argmin x 1 2 k c • ( y − Ax ) k 2 2 + λ R 1 k W T x k q + λ R 2 TV ( f ) , c =    c 1 . . . c ( kn ϕ ) × n 2    ∝ vec ( D ) • 1 1 + ε (23) where • denotes the element-b y-element multiplication, vec ( . ) converts th e matrix into a co lumn vector , and ε = [ ε 1 , ..., ε kn ϕ × n 2 ] , ε i ∈ [ 0 , ∞ ) repre sents the effect of interpolatio n er ror . The method used for calculating ε is illustrated in Figure 5. The value o f ε i for a lin e exactly b etween two p olar lin es is ∞ , since its distance fr om the true measured values ar e maximal and ther efore the er ror is maximal. Usin g (21) this can b e thoug ht as e i → ∞ an d consequen tly ε i → ∞ o r c i → 0. The closer the eq ually sloped lines are to the rays on which the measurem ents ar e do ne, the interpolatio n erro r g ets smaller and ε i ’ s on that line get clo ser to zero. Finally , if the desired eq ually sloped ra ys are exactly on the p olar lines, th e interpolatio n er ror e i is z ero wh ich is equiv alen t to ε i = 0. -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 ω dž ω LJ The line with    ∞ The line with    0 The desired Equally Slop ed lines with   ∈ 0, ∞ Fig. 5: Calculating the ef f ect of the in terpolation error for inclusion in the Er ror Ada ptation W eight (EA W). I V . E M S O L U T I O N F O R T H E G E N E R A L I Z E D C S M O D E L It has p reviously bee n shown that the qu ality of the recon - structed imag e can be improved by co mbining the sp arsity and to tal variation p enalty term s [1 6] both of wh ich are incorpo rated in (23). It is thus the equa tion used in th is p aper to r ecover CT images from und ersampled d ata. In solving this optimization pr oblem, we use th e Fast Co mposite Splitting Algorithm (FCSA) [32 ] to de compose (23) into two simpler sub-pro blems similar to (15) and (1 3), which are given b y: ˆ x 1 = argmin x 1 2 k c • ( y − Ax ) k 2 2 + λ R 1 k W T x k q ˆ x 2 = argmin x 1 2 k c • ( y − Ax ) k 2 2 + λ R 2 TV ( f ) (24) Calculating ˆ x 1 and ˆ x 2 , FCSA pro poses th at the solution to the pro blem can be obtained by a linear comb ination of th e solutions of the two sub- problem s, i.e. , ˆ x = δ ˆ x 1 + ( 1 − δ ) ˆ x 2 (25) in which δ is a weigh t that is defined as a fun ction of the value of the objective f unctions o f th e two subprob lems, denoted as f 1 and f 2 and is gi ven by δ = f 2 / ( f 1 + f 2 ) . Therefore, an FCSA based EM method is proposed to recover the CT images from the X-r ay p rojections and is called the FCSA- LEM algorithm. A. Latent variable and EM meth od for ℓ 2 − ℓ 1 and ℓ 2 -TV Optimization Here an ef ficien t method is p roposed to solve ℓ 2 − ℓ 1 and ℓ 2 -TV sub problems in (24). An Expectation- Maximization algorithm is u sed to solve these two op timizations prob lems [42], [43] for a CT modeled by (7). A latent variable z is defined su ch tha t th e problem in ( 7) can b e written as 1 : y = Az + n 1 (26) in wh ich z is chosen to be: z = x + α n 2 = W θ + α n 2 (27) The no ise is split in to two parts: n = α A n 1 + n 2 and p ( n 1 ) = N ( n 1 ; 0 , I ) , p ( n 2 ) = N ( n 2 ; 0 , D − 1 − α 2 AA T ) . Using this nota- tion, the E M algorithm is a s f ollows: E-step: Compute the conditional expec tation o f the log like- lihood, giv en the observed data and the current estimate θ ( t ) = W T x ( t ) : Q ( θ , θ ( t ) ) = E [ log p ( y , z | θ ) | y , θ ( t ) ] (28) M-Step: Up date th e estimate: Q ( t + 1 ) = arg max θ { Q ( θ , θ ( t ) ) − h ( θ ) } (29) In th e E-step, z ( t ) = E [ z | y , θ ( t ) ] sh ould be calculated a nd plugged in (28), in which we need to c alculate the like- lihood p ( y , z | θ ) = p ( y | z , θ ) p ( z | θ ) = p ( y | z ) p ( z | θ ) . Sin ce p ( y | z ) ∼ N ( y ; Az , D − 1 − α 2 AA T ) and p ( z | θ ) ∼ N ( z ; 0 , α 2 I ) , log p ( y , z | θ ) = − k x − z k 2 2 α 2 + K = − x T x − 2 x T z 2 α 2 + K ′ and therefore we have th e fo llowing equ ation fo r Q ( θ , θ ( t ) ) : Q ( θ , θ ( t ) ) = − k x − z ( t ) k 2 2 α 2 + K = − x T x − 2 x T z ( t ) 2 α 2 + K ′ (30) in which K and K ′ do no t depend on x and as a resu lt θ . Since p ( z | y , ˆ θ ( t ) ) ∝ p ( y | z ) p ( z | ˆ θ ( t ) ) , in which b oth p ( y | z ) and p ( z | ˆ θ ( t ) ) ar e Gaussian, p ( z | y , ˆ θ ( t ) ) is Gau ssian with mean value of z ( t ) = E [ p ( z | y , ˆ θ ( t ) )] [45], [46]: z ( t ) = x ( t ) + C z A T ( A C z A T + C n 2 ) − 1 ( y − Ax ( t ) ) (31) in wh ich C n 2 = D − 1 − α 2 AA T and C z = α 2 I . Therefor e, th e E-step can b e summarized by the calculation o f: z ( t ) = x ( t ) + α 2 A T D ( y − A x ( t ) ) = W θ ( t ) + α 2 A T D ( y − A W θ ( t ) ) (32) 1 The definition and the eff ect of this latent varia ble in the final solution is similar to the varia ble splitting strategy used in alternatin g direction methods (ADM) [44], [16]. 7 By inclusion of the rebinnin g e rror, i.e. using the EA W introdu ced in (23), this step will be: z ( t ) = x ( t ) + α 2 A T ( c • y − c • A x ( t ) ) = W θ ( t ) + α 2 A T ( c • y − c • A W θ ( t ) ) (33) that will be followed by an M-step: θ ( t + 1 ) = argmax θ {− k W θ − z ( t ) k 2 2 α 2 − h ( θ ) } = argmin θ { k θ − W T z ( t ) k 2 2 α 2 + h ( θ ) } x ( t + 1 ) = argmin x { k x − z ( t ) k 2 2 α 2 + h ( x ) } (34) Equation (3 4) has a clo sed form solution for some special cases, e.g. soft thresholdin g if h ( θ ) = λ R k θ k 1 [47] and th e TV denoising prob lem if h ( θ ) = λ R TV ( W θ ) [48]. B. Solving the Gen eralized CS Mo del The pseu docode shown in Alg orithm 1 ou tlines the FCSA- LEM metho d used to solve the gen eralized CS problem derived in section III-C, using the meth od intro duced in ( 32)- (34). In this algorithm prox L { g ( x ) , z } = argmin x g ( x ) + L 2 k x − z k 2 2 . Th e optimization p roblem in step 2 of Algorith m 1: ˆ x 1 = W ( prox 1 / α 2 { λ R 1 k θ k 1 , W T z } has a closed form solution given by : ˆ x 1 = W ( sign ( W T z ) max { 0 , | W T z | − λ R 1 α 2 } ) (35) T o calculate ˆ x 2 in step 3, a total v ar iation minimization s cheme is u sed and solved by a split Bregman based method [48] 2 . Alg. 1 Pseudocode of the FCSA-LEM algorithm used t o solve the optimization problem. Initialize: α , λ R 1 , λ R 2 , c , r 1 = 0, t 1 = 1, maxiter, tol while k ˆ x k − ˆ x k − 1 k 2 k ˆ x k k > tol o r k < maxiter do 1 z = r k + α 2 A T ( c • y − c • A r k ) 2 ˆ x 1 = W ( prox 1 / α 2 { λ R 1 k W T x k 1 , W T z } ) 3 ˆ x 2 = prox 1 / α 2 { λ R 2 TV ( x ) , z } 4 ˆ x k = δ ˆ x 1 + ( 1 − δ ) ˆ x 2 5 t k + 1 = 1 + q 1 + 4 t 2 k 2 6 r k + 1 = ˆ x k + ( t k − 1 t k + 1 )( ˆ x k − ˆ x k − 1 ) 7 k ← k + 1 endwhile 2 The final steps are very similar to iterati ve soft t hresholdi ng based methods [19]. V . R E S U LT S In this section we present results obtained u sing the pro- posed algorithm with fan and helical cone beam geometries using a She pp-Log an software p hantom av ailab le in MA TLAB (MathW orks, Massachusetts, USA), a custom made p hantom which mimics d ifferent cardiac plaqu es, and a che st scan from a hospital patient. MA TLAB w as u sed to calculate the equiangular fan beam projection s thr ough the Shep p-Logan phan tom. Th e X-ray projection s of the p laque phanto m an d the p atient were taken using a T oshiba Aq uilion ONE c  scanner (T oronto Gener al Hospital, Can ada), wh ich ha s 89 6 detecto rs and 320 r ows of detectors. The scanner gathers data from 1200 pr ojections in each 360 ◦ rotation. When images were reconstructed from fewer than the 1200 projections the projection vie ws were selected equiangu larly . For all the scan proto cols used th e X-ray tube cu rrent-exposu re time pro duct was 50mAs and the p eak voltage was 120k V . Data from the central row o f a volumetric scan o n one sing le rotation served as th e fan beam data. T o s implify the EA W calcu lation, α was chosen to be a diagona l matrix wh ose elemen ts were α i ∝ 1 / d i so that the e lements o f EA W will be c i = 1 / ( 1 + ε i ) . Figure 6 com pares the Shepp-L ogan phantom recon structed from 1 28 pro jections using 1) the inverse pseudo po lar Fourier transform ( using the least squ ares metho d), 2) an iterative soft threshold- based method (T wIST) [1 9], and 3) the propo sed FCSA-LEM metho d. Based on th e same phantom Figu re 7 compare s the accuracy of the reconstru ction error f or all three methods as the nu mber of p rojections is varied from 5 0 to 1000. Both of th ese figures show that using EA W improves the recovery accuracy . In particular Figure 7 shows that the use of m ore than 256 projec tions for a 512 × 5 12 imag e does not significantly affect the reconstru ction accuracy . A major improvement in the proposed method , com pared to the other CS-based reconstructio n method s, is its mu ch reduced computational burden. Figure 8 compares the recovery time of 1) filtered back projection (FBP), 2) th e proposed method (FCSA-LEM), and 3) an AR T -TV based metho d [10], [11], which is basically an algebr aic r econstruction followed by a TV smoo thing at eac h step. The compu ter used for the simulation is a desktop i5 comp uter with 16GB of RAM. Using this computer we could n ot use AR T -TV m ethods with a resolution higher than 128 × 128 pixels due to memory constraints in MA TL AB. It can be seen that the recovery time for the p roposed method g ets closer to the time of FBP as the image approa ches 1 024 × 10 24 pixels. Figure 9 com pares the plaqu e phan tom reconstruc ted with FBP from 1200 projections an d the phanto m reconstru cted with the prop osed m ethod from 256 projections. It can b e seen that the image quality is almost the same with an e rror less th an 1 %. Reconstruction s of a ch est CT scan from a hospital patient using FBP from 1200 pro jections an d the propo sed meth od from 256 projectio ns is sho wn in Figure 10. I t is evident that the im age recon structed with the prop osed method has almost th e same qua lity as FBP , which has abo ut 5 time s more 8 (A) (B) (C) (D) Fig. 6: (A) Original Shepp-Log an p hantom Image. Recon- structions u sing 1 28 pr ojections with (B) inverse pseud o polar Fourier transform using the least squ ares m ethod (nor malized error ≈ 0 . 9 ), (C) an iterative soft thresholdin g based method (T wIST) (normalized error ≈ 10 − 1 ), and (D) the prop osed FCSA-LEM metho d ( normalized erro r ≈ 10 − 2 ). The reb inned parallel rays are used in all three methods to recover the image. 100 200 300 400 500 600 700 800 900 1000 10 −2 10 −1 Num ber of Pro jections Normal ized Error = k x − ˆ x k 2 k x k 2 Tw IST Least Square i nv ers e PPFT F C SA-LEM Fig. 7: Normalized recon struction err or fo r the Shep p-Logan phantom reconstructed with the in verse pseudo polar Fourier transform (PPFT) using least squ ares (LS) method, using an iterativ e soft thresho lding based meth od ( T wI ST), an d the propo sed FCSA-LEM method. Rebin ned parallel ray s were used in all three method s to recover the image. projection s, i.e. 5 times the radiation dose. While th e image from th e propo sed me thod is a little blur ry comp ared to FBP reconstruc ted ima ge, but all the details are preserved. Figure 11 sho ws a simple simulated phantom reconstructed by the prop osed helical re construction method . The h elix so urce position is defin ed as ψ = [ R cos ( ϕ ) , R sin ( ϕ ) , P ϕ 2 π ] an d in this 10 2 10 3 10 −3 10 −2 10 −1 10 0 10 1 10 2 Image Siz e (Pix els) Reconstr ucti on Time (s ec) Fi ltered B ack Pro jection F C SA-LEM AR T -TV based met ho d Fig. 8: Reconstru ction time comparison using a standard desktop com puter, for: ( 1) fan be am filtered b ack pro jection (FBP) rec onstruction, ( 2) the propo sed metho d ( FCSA-LEM), and (3) a fan b eam AR T -TV based m ethod. (A) (B) Fig. 9 : Compar ison of FBP an d the p roposed meth od for a car diac p laque ph antom made in T oronto Gener al Hospi- tal. The scan p rotocol was 50mAs and 1 20kVp. (A) Image reconstruc ted from 1200 pr ojections using FBP . (B) Image reconstruc ted from 2 56 p rojections with the propo sed method . (A) (B) Fig. 10: Com parison of FBP an d the proposed metho d f or the chest CT scan data from a patient. The scan proto col was 50 mAs and 120 kVp. (A) I mage reconstru cted f rom 1200 projection s with FBP . (B) Im age reconstructed fro m 256 projection s with the p roposed meth od. test pitch factor P = 0 . 5 . As can be seen, aside from the start or end o f the scan the reconstru ction is almost pe rfect. Howev er , when the image is close to o ne of the end points the error of 9 rebinnin g increases an d a s a result the imag e reco nstruction error increases. V I . C O N C L U S I O N It has been sho wn that CT recon struction can be statistically mode led as a weighted compressed sensing optimization p roblem. The propo sed mo del was derived from a MAP mode l o f CT imaging with sparsity and piec e wise linearity co nstraints. T o solve the pro posed model a fast CS recovery meth od was prop osed in which p seudo polar Fourier transform was used as the measureme nt fu nction to reduce the computatio nal co mplexity . Mor eover , to b e able to rec onstruct CT imag es fro m fan beam and helical cone beam p rojections, rebinnin g to parallel be ams was used. T o adapt the proposed CS recovery metho d to the interp olation er ror , a weigh ting approa ch ( EA W) was proposed, in which the weights accounted for the me asurement noise and inte rpolation errors. This enabled CT images to be recon structed f rom a reduced number of fan or helical con e be am X-r ay p rojections. It was shown that using EA W imp roves the reco nstruction q uality substantially . For instance , a 5 12 × 51 2 Sh epp-Log an pha ntom reconstruc ted with 128 projections using a conv entional CS method h ad ∼ 1 0% error . Howe ver, using the same data with our new method the recon struction error was as low as ∼ 1%. The pr oposed weighted CS-CT reconstruc tion model was solved with a pr oposed FCSA-EM based metho d, called FCSA-LEM. Th e low compu tational com plexity o f our FCSA-LEM method m ade f ast recov ery of the CT images possible. For e xample, we were able to recover a 512 × 51 2 image in less than 2 0 s ec on a desktop computer witho ut numerical optimiza tions, thus our prop osed meth od may be among the fir st CS-CT methods whose com putational complexity is within the rea lm of what could b e clinically relev ant today . Ackno wl edgements — W e thank the Canadian Mitacs-Accelerate program and T oshiba Canada for partial financial support. RSCC wishes to thank the N atural S ciences an d Engineering Council of Canada for support unde r grant #32 47-2012, and SMH is grateful for the award of a L oo Geok Graduate Scholarship. R E F E R E N C E S [1] D.J. Brenner and E.J. Hall, “Computed tomography - an increasi ng source of radiati on exposure, ” The New England Jou rnal of Medicine , vol. 57, pp. 2277–2284, Nov . 2007. [2] M.S. Pearce, J.A. Salott, M.P . L ittle , K. McHugh, C. Lee, K.P . Kim, N.L. Howe, C.M. Ronckers, P . Rajarama n, L. Parker A.W . Craft, and A.B. de Gonz ´ alezl, “Radiati on e xposure from CT scans in childhood and s ubsequen t risk of leuka emia and brai n tumours: a ret rospecti ve cohort study , ” The Lance t , v ol. 380, no. 9840, pp. 499– 505, Aug. 2012. [3] B. Desjardins and E.A . Kazerooni, “Revi e w: ECG-gated cardiac CT, ” American Journal of Roentgeno logy , vol. 128, no. 4, pp. 993–1010, Apr . 2004. [4] F . Zhu, T . Carpente r , D. Rodriguez Gonzale z, M. Atkinson, and J. W ard law , “Computed tomog raphy per fusion imaging denoising using gaussian process regre ssion, ” Physics in Medici ne and Biolo gy , v ol. 57, no. 12, pp. 183–198, Jun. 2012. [5] S.A. Hyder and R. Sukane sh, “ An effici ent algorit hm for denoising MR and CT images using digi tal c urvel et transform, ” Advances in Experiment al Medicine and Biology , pp. 471–480, Sept. 2011. [6] J. Bro wne and A.B. de Pierro, “ A row-a ction altern ati ve to the E M algorit hm for maximizing likel ihood in emission tomography , ” IEE E T ransaction on Medical Imaging , vol. 15, no. 5, pp. 687–699, Oct. 1996. [7] J. Ming and W . Ge, “Con ver gence of the simul taneous al gebraic recon- structio n techni que (SAR T), ” IEEE Tr ansacti on on Image Proce ssing , vol. 12, no. 8, pp. 957–961, Aug. 2003. [8] E.Y . Sidky , C.M. Kao, and X. Pan, “ Accura te image reconstruction from fe w-vie ws and limit ed-angl e data in di verg ent-bea m CT, ” J ournal of X- Ray Science and T echn ology , vol. 14, no. 2, pp. 119–139, Jan. 2006. [9] E.Y . Sidk y , X. P an, I.S. Reiser , R.M. Nishika wa, R.H. Moore, and D.B. Ko pans, “ Accurate image reconstru ction from few-vie ws and limited- angle data in di ver gent-be am CT, ” Medical Physics , vol . 36 , no. 11, pp. 4920–4932, 2009. [10] G.H. Chen, J. T ang, and S. Leng , “Prior image con strained compressed sensing (piccs): A met hod to accurately rec onstruct dyna mic CT imag es from highly undersampled projection data sets, ” Medic al Physics , vol. 35, no. 2, pp. 660–663, Feb . 2008. [11] H. Lee, L. Xing, R. Davi di, R. Li, J. Qian, and R. Lee, “Improved compressed sensing -based cone-beam CT reconstruct ion using adapti ve prior image constrai nts, ” Physic s in Medicine and Biology , vol. 57, no. 8, pp. 2287–2307, Feb . 2012. [12] T . Strohmer and R. V ershynin, “ A randomized kac zmarz a lgorithm w ith expo nential conv ergence, ” Journal of F ourier Analysis and Applicatio ns , vol. 15, no. 2, pp. 262–278, Apr . 2009. [13] G.T . Herman and L.B. Meye r , “ Algebraic reconstructi on techni ques can be made computatio nally effic ient, ” IEEE T ransaction on Medical Imagin g , vol. 12, no. 3, pp. 600–609, Sept. 1993. [14] E.J. Cand ` es, J. Romberg, and T . T ao, “Robust uncertain ty princ iples: exa ct signal reconstructio n from highly incomplete frequency informa- tion, ” IEEE T rans. Info. Theory , vol . 52, no. 2, pp. 489– 509, Feb . 20 06. [15] D.L. Donoho, “Compressed sensing, ” IE EE Tr ansacti on on Information Theory , vol. 52, no. 4, pp. 1289–1306, Apr . 2006. [16] Y . Junfeng, Z. Y in, and Y . W otao, “ A fast alternat ing direction method for TVL1-L2 signal reconstructi on from parti al fouri er data, ” IEEE J ournal of Sele cted T opics in Signal Pr ocessing , vol. 4, no. 2, pp. 288– 297, Apr . 2010. [17] B.P . Fa himian, Y . Mao, P . Cloete ns, and J. Miao, “Low- dose x-ray phase- contrast and absorption ct using equall y sloped tomography , ” Physics in Medicine and Biology , vol. 55, no. 18, pp. 5383–5400, Aug. 2010. [18] M. Jiang, J. Jin, F . Liu, Y . Y u, L. Xia, Y . W ang, and S. Cro zier , “Sparsit y- constrai ned SENSE reconstruc tion: an effic ient implement ation using a fast composite splitting algorithm, ” M agnet ic Resonance Imaging , vol. 31, no. 7, pp. 1218–1227, Sept. 2013. [19] J.M. Bioucas-Dias and M.A.T . Figueir edo, “ A ne w T wIST : T wo-ste p iterat i ve shrinkage /threshol ding algorithms for image restoration, ” IE EE T ransaction on Image Pr ocessing , vol. 16, no. 12, pp. 2992–3004, Dec. 2007. [20] D. Got tleib, B. Gustafsson, and P . Forssen, “On the direct fourier method for c omputer tomograp hy , ” IEEE T ransact ion on Medical Imaging , v ol. 19, no. 3, pp. 223–232, Mar . 2000. [21] H. Stark, J . W oods, I. Paul, and R. Hingorani, “Direct fourier recon- structio n in comput er tomography , ” IEEE T ransactions on Acoustics, Speec h and Signal P r ocessing , vol. 29, no. 2, pp. 237–245, Apr . 1981. [22] H. Schomber g and J. T immer , “The gri dding method for image reconstru ction by fourier tra nsformation , ” IEEE T ransac tions on Me dical Imagin g , vol. 14, no. 3, pp. 596–607, Sept. 1995. [23] Y . Mao, B.P . Fahimian , S.J. Osher , and J. Miao, “De velo pment and opti- mizatio n of regu larize d tomographic reconstructi on algorithms utilizi ng equall y-sloped tomogra phy , ” IEEE T ransaction on Image Proc essing , vol. 19, no. 5, pp. 1259–1268, May 2010. [24] J. Miao, F . F ¨ orster , and O. Levi , “Equally sloped tomograp hy with ov ersampling reconstruct ion, ” Physical Revie w B , vol. 72, no. 5, pp. 052103, Aug. 2005. [25] E. Lee, B .P . Fahi mian, C.V . Ianc u, C. Sulow ay , G.E. Mur phy , E.R. Wright, D. Castao-Dez, G.J. Jensen, and J. Miao, “Radia tion dose re- duction and image enhancement in biolog ical imaging through equall y- sloped tomography , ” J ournal of Structura l Biology , vol. 164, no. 2, pp. 221–227, 2008. [26] A. A verbuc h, R.R. Coifmanb, D. L. Donoho, M. Elad, and M. Israeli, “Fast and accurate pol ar fourier tra nsform, ” Applied and Co mputation al Harmonic Analysis , vol. 21, pp. 145–167, 2006. [27] B.P . Fahimia n, Y . Zhao, Z. Huang, R. Fung, Y . Mao, C. Zhu, M. Kha- tonabad i, J.J. DeMarc o, S. J. Osher , M. F . McNitt-Gr ay , and J. Miao, “Radia tion dose re duction in medical x-ray ct vi a four ier-b ased iterat i ve reconstru ction, ” Medical Physics , vol. 40, no. 3, pp. 031914–1–031914– 10, Mar . 2013. [28] M. Hashemi, S. Beheshti, P .R. Gill, N.S. Pa ul, and R.S.C. Cobbold, “Fast fan/p aralle l bea m CS-based lo w-dose CT reconstructi on, ” in Internati onal Confer ence on Acoustics, Speec h, and Signal Proc essing (ICASSP) . IEEE, May 2013. [29] E. J. Cand ` es, M. W ak in, and S. Boyd, “Enhancing sparsit y by re weighted ℓ 1 minimizat ion, ” Journal of F ourier Analysis and Applications , vol. 14, pp. 877–905, 2008. [30] M.A. Khajehneja d, X. W eiyu, A.S. A vestimehr , and B. Hassibi, “ Ana - lyzing weighted ℓ 1 minimizat ion for sparse recov ery with nonuniform 10 (A) (B) (C) Fig. 11: Helical scan tested on a simple simulated p hantom. Pitch factor is 0. 5 in th is p hantom data. ( A) Th e o riginal ph antom. (B) Image reconstructe d with the prop osed m ethod. ( C) Difference between th e true image and the r econstructed image. sparse models, ” IEEE T ransactions on Signal P r ocessing , vol. 59, no. 5, pp. 1985–2001, 2011. [31] N. V asw ani a nd L. W ei, “Modified-CS: Modifying compressiv e s ensing for problems with par tially kno wn support, ” in IEEE Int ernation al Symposium on Information Theory , 2009. ISIT 2009. , 2009, pp. 488– 492. [32] A. Be ck and M. T eboull e, “ A f ast iterati ve shrinkage-t hresholdin g algorit hm for linear inv erse problems, ” SIAM Journal on Imaging Scienc es , vol. 2, no. 1, pp. 183–202, Mar . 2009. [33] S.R. Deans, The Radon T ransform and Some of Its Applications , John W ile y & Sons, 1983. [34] J. Hsieh, Compu ted T omograph y: Princ iples, Design, Artifacts, an d Recent Advances , vol. PM188, SPIE Press Book, Bellingham, W A, 2 editi on, 2009. [35] J.B. Thiba ult, D.K. Sauer , C.A. Bouman, and J. Hsieh, “ A three- dimensiona l statisti cal approach to improv ed image quality for multisl ice helic al CT, ” Medical Physic s , vol. 34, no. 11, pp. 4526–4544, Aug. 2007. [36] S. Maeda, W . Fukuda, K. Atsunori, and S. Ishii , “Maximum a posteriori X-ray compute d tomography using gra ph cuts, ” 2010, v ol. 233, pp. 4526–4544. [37] M. Hashemi and S. Behesht i, “ Adapti ve bayesian denoising for ge neral gaussian distribute d (ggd) signals i n wav elet domai n, ” , Jul. 2012. [38] M. Basse ville , A. Ben veniste, K.C. Chou, S .A. Golden, R. Nikouk hah, and A.S . Wil lsky , “Modeling and e stimation of multiresolutio n stochas- tic processes, ” IEEE T ransction on Information Theory , vol. 38, no. 2, pp. 766–784, Mar . 1992. [39] M. Lassas and S. Silt anen, “Can one use total v ariat ion prior for edge- preservin g bayesia n in version?, ” In verse Probl ems , vol. 20, no. 5, pp. 1537–1563, Aug. 2004. [40] G. Besson, “CT image reconstru ction from fan-pa rallel data, ” Medical Physics , vol. 26, pp. 415–426, 1999. [41] F . Noo, M. Defrise, and R. Cl ackdo yle, “Single-sli ce rebinning method for helical cone-beam CT, ” Physic s in Medicin e and Biology , vol. 44, no. 2, pp. 561–570, Feb . 1999. [42] Q. Kun and A. Dogan dzic, “Sparse signal reconstruct ion via ECME hard thresholdi ng, ” IEE E Tr ansacti ons on Signal Pr ocessing , vol. 60, no. 9, pp. 4551–4569, Sept. 2012. [43] M.A. Figuei redo and R.D. Nowak , “ An EM algorithm for wa velet-base d image re storatio n, ” IEEE T ransaction on Image Proce ssing , v ol. 12, no. 8, pp. 906–916, Aug. 2003. [44] J.F . Y ang and Y . Zhang, “ Alternati ng direction algorit hms for ℓ 1 - problems in compressi ve sensing, ” SIAM Journal on Scientific Com- puting , vol. 33, no. 1, pp. 250–278, 2011. [45] L. Scharf and C. D emeure, Statistical si gnal proc essing: detection, estimati on, and time series analysis , Addison-W e sley Pub . Co., Reading, MA, 1991. [46] M.A. Figuei redo and R.D. Nowak , “ An EM algorithm for wa velet-base d image re storatio n, ” IEEE T ransaction on Image Proce ssing , v ol. 12, no. 8, pp. 906–916, Aug. 2003. [47] M.A.T . Figueiredo and R.D. Now ak, “W avele t-based image estimati on: an empirical bayes approach using J ef fre y’ s noni nformati ve pri or , ” IEEE T ransaction on Image Pro cessing , vol. 10, no. 9, pp. 1322–1331, Sept. 2001. [48] T . Goldst ein and S. Osher, “The split bregman method for ℓ 1 -regu larize d problems, ” SIAM Journa l on Imaging Sciences , vol. 2, no. 2, pp. 323– 343, Apr . 2009.