Efficient low-dose CT artifact mitigation using an artifact-matched prior scan WeiXua) andKlausMuellerb) VisualAnalyticsandImagingLab,CenterofVisualComputing,ComputerScienceDepartment,StonyBrook University,StonyBrook,NewYork11794-4400 (Received10November2011;revised5June2012;acceptedforpublication12June2012; published19July2012) Purpose: Low-dose CT has attracted increasing attention due to growing concerns about radiation exposure in medical scans. However, the frugal use of x-ray radiation inevitably reduces the qual- ity of the CT images, introducing artifacts such as noise and streaks which make the reconstructed images difficult to read in clinical routine. For follow-up CT exams a prior scan is often available. It typically contains the same anatomical structures, just somewhat deformed and not aligned. This workdescribesatwo-steptechniquethatutilizesthispriorscantoachievehigh-qualitylow-doseCT imaging,overcomingdifficultiesarisingfromnoiseartifactsandmisalignment.Wespecificallyfocus onreducingthedosebyloweringthenumberofprojections.Thisgivesrisetoseverestreakartifacts whichpossiblylowerthereadabilityofCTimagestoalargerextentthanthefine-grainednoisethat resultsfromloweringthemAorkVsettings. Methods:Acommonapproachistoapplyimagefilteringtoreducethenoiseartifacts.Thesetech- niquestypicallyutilizepixelneighborhoodsinthedegradedimagetoestimatethetruevalueofapixel atthecenterofthisneighborhood.However,thiscanleadtopoorresultswhentheimageisseverely contaminatedunderverylowlow-dosesituations.Weproposeamethodthatutilizesthenondegraded, clean prior to determine higher quality pixel statistics to form the pixel estimates, supported by the matching scheme of the non-local means filter. To make this matching reliable, a good registration of prior and low-dose image is required. For this, we employ a state-of-the-art registration method, called SIFT-flow, which can tolerate the high amount of streak noise. But even for properly regis- teredimages,usinganartifactfreepriorforthematchingyieldsinferiorresults.Wehencedescribe a scheme that first constructs a tandem-prior with streak artifacts resembling those in the low-dose image,andthenemploysthisimageforthematching,butusesthecorrespondinghigh-qualityprior todeterminethepixelestimates. Results:Twoexperimentalstudiesarecarriedout,usingaheadphantomandahumanlungwithpro- jectionsgatheredviasimulation.Weassessthequalityoftheprocessedreconstructionwithvarious metrics:mathematicalandperceptual.Wefindthatthequalitythatcanbeobtainedwiththeartifact- matched prior-based scheme significantly exceeds that of all competing schemes. Even though the generalprior-basedapproachisabletoeliminatethestreakartifacts,onlytheartifact-matchedscheme canrestoresmalldetailandfeaturesharpness. Conclusions: The reduced-projection low-dose image reconstruction algorithm we present outper- forms traditional image restoration algorithms when a prior scan is available. Our method is quite efficientandassuchitiswellsuitedforfast-pacedclinicalapplicationssuchasimage-assistedinter- ventions,orthopedicalignmentscans,andfollow-ups.©2012AmericanAssociationofPhysicistsin Medicine.[http://dx.doi.org/10.1118/1.4736528] Keywords:artifact-matchedprior,artifactmitigation,low-doseCT,nonlocalmeans,reference-based I. INTRODUCTION adverseradiationeffectsofCTimaginghavedampenedthese developments.1,2 Due to these studies, the harmful effects Computed x-ray tomography has revolutionized modern of x-ray radiation in CT scans have become publicly heard, medicine and thanks to the rapid growth in scanner technol- threateningthefutureofthismodality.Tocounterthesecon- ogy the gamut of its applications has risen at an enormous cerns,campaignssuchasImageGently(Ref.29)andImage- rate. In this process, buoyed by the excitement of possibil- Wisely (Ref. 30) have been initiated that promote the opti- ities little attention was paid to the radiation dose adminis- mizationoftheradiationdoseusedinbothpediatricandadult tered to the patient. Scans with ever-improving spatial and medicalimaging. temporal resolutions were conducted on a routine basis and To reduce the radiation dose subjected to the patient one the associated CT reconstruction algorithms had the luxury can: (1) lower the number of scans, (2) lower the number of of an abundance of data collected at each exam. It was only x-ray projections per scan, and (3) lower the energy settings recently that the sobering results of long-term studies on the of the x-ray tube (kV, mA) per projection image. The first 4748 Med.Phys.39(8),August2012 0094-2405/2012/39(8)/4748/13/$30.00 ©2012Am.Assoc.Phys.Med. 4748 4749 W.XuandK.Mueller:Efficientlow-doseCTartifactmitigation 4749 measure, i.e., reducing the number of scans, is often left at at two different thicknesses, using the same acquired pro- thediscretionofthetreatingphysician.Thelattertwooptions jection data. They reconstruct the thicker slices from bin- are highly detrimental to image quality, resulting in images averagedprojectionswhichincreasesSNR,whilethethinner with significant noise artifacts. They greatly challenge the (andnoisier)slicesarereconstructedfromtheoriginalprojec- conventional CT reconstruction algorithms based on analyt- tiondata.Sinceregistrationisimplicit,itisrelativelystraight- icalformulationsrootedintheinverseradontransform,3 and forward to use the thicker slices for neighborhood-based de- theseshortcomingshaverecentlyinvigoratedresearchefforts noisingofthethinnerslices.Incontrast,ourmethodappliesto towards methods that seek alternatives to these conventional settingsinwhichthereferenceimagesarenotnecessarilyac- schemes. quiredsimultaneously.Yuetal.13presentthemethodprevious A popular approach to this end has been to enforce data scan-regularized reconstruction (PSRR). It replaces regions fidelityandimagequalityasajointoptimizationproblemand thatareunchangedinalow-doseCTreconstructionwiththeir solvethesetwopartsinaniterativeround-robinfashion.Data direct embodiments in a normal-dose CT reconstruction and fidelitycanbeassuredbywaysofanyCTreconstructional- uses a nonlinear diffusion approach for denoising in the re- gorithm,iterativeoranalytical,butmostusetheformer.The mainingregions.Thisapproachrequiresaneffectivestrategy reduction of noise artifacts, on the other hand, can be posed forfeaturerecognition,whichtheauthorsaccomplishviareg- as an image denoising problem. Many approaches use the istration. method of total variation minimization (TVM) (Ref. 4) for The approach most similar to ours is that of Ma et al.14 this task since it is often part of general compressive sens- They also use a registration algorithm—a combination of ingformulations5thatwereoriginallyprescribedtodealwith rigid principal component analysis (PCA) and nonrigid mu- sparsedata.ForCTreconstruction,anumberofsophisticated tual information optimization15—for rough alignment. They schemes have been developed that adapt the various param- then use a neighborhood-mechanism to locate suitable re- eters used in the process, such as ASD-POCS (Ref. 6) and placementcandidatesinapriorscanofthepatient.Ourwork soft-thresholdfiltering.7 differs from theirs in the following important ways. First, Inthiswork,wehaveattemptedtodeviseaframeworkthat whileMaetal.restoreimagesreconstructedfromprojections executesthedatafidelitystepandtheimagequalitystepeach generated at reduced mA settings, we treat artifacts caused exactly once. It is hence of lower computational complex- bythereductionofprojections.Theresultingstreakartifacts itythanthepresentschemeswhichperformthesestepsitera- are much more severe and irregular than the random noise- tively.Weachievethisbymakingcreativeuseofanartifact- artifacts caused by low-mA imaging. Second, instead of us- freeprior—constitutedbyanexistingregular-dosescanofthe ingacleanpriorformatchingweuseapriorwithsimulated patient. Such a clean prior scan is frequently available. For artifacts. We find that this affords much better accuracy, as example, it may be a regular-dose first scan acquired before wewillproofanddemonstrate.Apreliminaryversionofthe a low-dose follow-up scan or it may be a regular-dose diag- framework we describe here has been presented in Ref. 16, nostic scan preceding a low-dose setup scan for a surgical whichpredatestheworkbyMaetal.slightly. intervention such as orthopedic spine fixation, among other Our paper is organized as follows.Section II presents the scenarios. algorithms we have studied, Sec. III presents results, and Inthepresentwork,wefocusonthesecondformoflow- Sec.IVendswithconclusionsandpointstofuturework. dose CT, i.e., reducing the number of x-ray projections per scan.Thefirststepofourframeworkusesfilteredbackprojec- II. METHODSANDMATERIALS tion(FBP)toreconstructanimagewithsignificantstreakar- Tolocategoodpixelmatchesinthepriorweusethesim- tifactswhichresultfromthislownumberofprojections.The ilarity measures also employed by nonlocal means (NLM) useofFBPtoprovideaquickfirstestimateisacommonstrat- filtering.17 NLM filtering can be seen as a generalization of egy.Unlikethefirststepofaniterativescheme,suchasART Gaussian smoothing. It looks for structurally similar pixel (Ref.8)anditsderivatives,9,10 FBPtypicallyreconstructsall neighborhoodsinthesmoothingsite’sproximityandincludes imagefeaturesatgoodfidelitybutthehighimagenoisemakes them into the Gaussian filter statistics. This leads to a more themdifficulttoread.CommonapproachesthenfollowFBP robustestimateofthetruepixelvalueandconsequentlytoim- byaniterativepipelinefordenoising.Oursecondstep,onthe provedimagerestoration/denoisingresults.Wehaverecently other hand, uses a single prior-based image restoration that informally compared NLM filtering with TVM for low-dose eliminates the noise and so provides the desired viewing ex- imagingtasksandourresultshavebeenquiteencouraging.18 perience.Inthisstep,wefirstregister/alignourpriorwiththe In this current work, we do not employ NLM-filtering in a FBP-estimateusinganestablishedmultiscalefeatureregistra- conventional way as a nonlocal extension of Gaussian filter- tion algorithm, i.e., SIFT flow.11 Following, we simulate the ing.Rather,weonlyretainNLM’smechanismforsimilarity- low-dosestreakartifactsoftheFBP-estimateinthisregistered based pixel neighborhood matching in the artifact-matched cleanprior.Finally,foreachpixelinthetargetimageweusea prior. neighborhoodsimilaritymetrictodeterminethebestmatches inthecontaminatedpriorandthenreplaceitusingthecorre- II.A. StandardNLMfiltering spondingpixelsinthecleanprior. Using existing scans to support regularization is not new. The NLM algorithm was proposed by Buades et al.17 for Kelmetal.12 describeanapproachthatreconstructsvolumes image denoising. It takes advantage of the high degree of MedicalPhysics,Vol.39,No.8,August2012 4750 W.XuandK.Mueller:Efficientlow-doseCTartifactmitigation 4750 Target Target neighborhood pixel Candidate 1 neighborhood Candidate pixel 1 Candidate 2 neighborhood Candidate Search pixel 2 window for target pixel FIG.1. IllustrationofconventionalNLMfordenoising. (a) NLM-filtered (b) TVM-filtered redundancy that typically exists in an image. Given a tar- FIG.2. FilteringresultsobtainedwithNLMandTVM. get pixel subject to denoising, it defines a small Gaussian- weighted region around it, called a patch. It then searches be high. Nevertheless, our use of the NLM-based matching the entire image for similar patches and accumulated them mechanism relaxes the need for tight and laborious registra- weightedbytheirdegreeofsimilarity.Inpractice,onlyalo- tionofthepriorimageasitperformsthefineregistrationon cal neighborhood around the target pixel is searched, called the fly via its search mechanism. As such it is less sensitive searchwindow.Thishelpsperformancebutitalsobettertol- to spatial distortions than the PSRR approach. For registra- erates nonstationary noise processes. Further, since we first tion we have made use of the SIFT-flow algorithm, recently registerthepriortothetargetimage,wedonotrequirelarge published by Liu et al.11 This algorithm originates from the searchwindowsinanycase.Moreformally,theupdatedvalue optical-flow algorithm which produces dense, pixel-to-pixel p(cid:2)xofatargetp(cid:3)ixeliscomputedas (cid:4) correspondences between two images. It extends the match- (cid:2) (cid:2) ing from raw pixels to SIFT feature descriptors.20 A scale- exp − Ga(t)|px+t −py+t|2/h2 ·py invariantfeaturetransform(SIFT)featuredescriptorcaptures p(cid:2) = y∈Wx (cid:3)t∈P (cid:4) . (1) thehistogramofgradientorientationsinalocalneighborhood x (cid:2) (cid:2) exp − Ga(t)|px+t −py+t|2/h2 at a given scale. It is well suited to characterize salient local y∈Wx t∈P andtransform-invariantimagestructuresandatthesametime Here,xisthelocationofthetargetpixelandtheyaretheloca- encode contextual information. SIFT-flow has been specifi- tionsofthecandidatepixels,withvaluesp .W isthesearch cally designed for scene matching, where objects share sim- y x window around x, and P is the patch size of each pixel. The ilarscenecharacteristicsbutmayhavedifferentappearances patchsimilarityismeasuredbytheGaussianweightedL dis- andlocateatdifferentplaces.Thisisthecaseintheregistra- 2 tancebetweentwopatchvectorswithtrepresentingtheindex tion of the prior to the current scan. They are certainly from withinapatchandG beingaGaussiankernelwithstandard thesameperson,i.e.,theysharethescenecharacteristics,but a deviationa.Theexponentialfunctionconvertsthesedistances theywilllikelyhavedifferentSNRandundergonedistortions to weights, determined by a parameter h which controls the ofthefeatures. overallsmoothnessofthefiltering.Largervaluesofhwillre- The implementation of SIFT-flow has two parts: (i) gen- sultinmoresmoothing. erate dense SIFT features where each pixel has a 128- Figure 1 shows an illustration of this process when NLM dimensionalSIFTvector,and(ii)findthecorrespondenceof isusedinaconventionalwaytodenoiseareconstructionwith theseSIFTfeaturesviadiscreteoptimizationontheimagelat- severestreakartifacts.InthisparticularcaseweusedFBPto ticetoobtainthedisplacementfieldforalignment.Forthefirst reconstructaGEheadphantomfrom45fan-beamprojections part, Fig. 3 shows a typical SIFT feature descriptor summa- acquiredover360◦.Theillustrationshowsthesearchwindow, rizingthegradientorientationsina162 pixelarea(asplotted the target pixel in the center and two candidate pixels with insidetheredsquare).Thegradients(shownasbluearrows) similar neighborhoods as the target pixel. Figure 2(a) shows areGaussian-smoothedaccordingtotheirdistancetothearea theNLM-filteredresult,whileFig.2(b)showstheresultsob- center.Thisareaispartitionedinto42 blocks,eachofsize42 tained with TVM. We can observe that both of these filters pixels(shownasgreensquares).Thegradientorientationsare provide some amount of improvement over the original im- thenaccumulatedineachblockto8orientationbinsandare ageshowninFig.5(a).Itislikely(see,forexample,Ref.19) weightedbytheirgradientmagnitudes.Thereareatotalof16 that repeating the fidelity and denoising steps several times 8-bin orientation histograms (with each red arrow represent- would do substantially better, but since we have aimed for a ing one bin). Thus the dimension of a SIFT vector is 4 × 4 nonrepeating approach—one in which a prior is available to ×8=128overa16×16=256area.Forthesecondpart,to aidinthedenoising—wefocusonasinglestepscheme. estimateflow,theenergyfunctionforSIFTflowisdefinedas below (cid:5) II.B. RegistrationusingtheSIFT-flowalgorithm E(f(x,y)) = min((cid:5)s (x,y) 1 Acrucialelementinourprior-assistedframeworkisproper (x,y) registrationsincewithoutitthepossibilityformismatchescan −s (x+f (x,y),y+f (x,y))(cid:5) ,a) 2 x y 1 MedicalPhysics,Vol.39,No.8,August2012 4751 W.XuandK.Mueller:Efficientlow-doseCTartifactmitigation 4751 Target Target neighborhood pixel Candidate 1 neighborhood Candidate pixel 1 Candidate 2 neighborhood Candidate Search pixel 2 window for target pixel FIG.4. IllustrationofR-NLMwiththeregisteredcleanpriorusedforboth matchingandretrieval. more distortion-tolerant the algorithm becomes, but the po- tentialforlostdetailandoversmoothingalsorises.Inexperi- mentswefoundasizeof7×7pixelsforbothsearchwindow FIG.3. IllustrationofSIFTdescriptorsummarizingedgeorientationsover andpatchestorepresentagoodcompromise. 16×16pixelarea. Our R-NLM algorithm first uses SIFT-flow to align the prior with the current scan, call it target scan. Following, it (cid:5) visitseverypixelinthetargetscan,placesthesearchwindow + b(|f (x,y)|+|f (x,y)|) x y in the same location in the prior scan, and uses the NLM- (x,y) (cid:5) algorithm to determine the update. Figure 5 presents some + (min(c|f (x,y)−f (x(cid:2),y(cid:2))|,d) resultswehaveobtainedwithourR-NLMalgorithm,tomo- x x (x(cid:2),y(cid:2))∈W(x,y) tivate a further extension discussed in Sec. II.D. Figure 5(f) shows a regular-dose reconstruction obtained with 360 pro- + min(c|f (x,y)−f (x(cid:2),y(cid:2))|,d)), (2) y y jections,whileFig.5(a)showsalow-doseFBPreconstruction where (x, y) and (x(cid:2), y(cid:2)) are pixel locations, s and s are obtained from the same data but only 45 projections—about 1 2 SIFT descriptors, f is the displacement function with f in 1/8ofthedose.Figure5(b)showstheprior.Sincethiswasa x x-directionandf iny-direction,Wisthepixelneighborhood headphantomthatcouldnotbewarpedmechanically,weper- y anda,b,c,darefourthresholds(withdefaultsettings).This formed a digital warp—in this case a twirl distortion around function is designed according to three constraints: (1) the the center of the image. This deformation field is shown in matched pixel should have similar SIFT descriptors; (2) the Fig.5(h).Figure5(c)showstheregisteredpriorandFig.5(d) displacementshouldbeassmallaspossible;and(3)adjacent shows the result obtained with R-NLM filtering. It is clearly pixels should have similar displacements to maintain flow better than the NLM and TVM filtered results (see Fig. 2). smoothness. This discrete displacement function can be This of course is a comparison that is only partially fair be- estimated by optimizing the energy function with a belief causeR-NLMhadaccesstoacleanprior,whiletheothersdid propagation algorithm.11 Its time complexity is O(h2logh) not.ThisprocedurechangesEq.(1)intothefollowing: (cid:3) (cid:4) wherehisthewidth(height)oftheimage.Beforeregistration (cid:2) (cid:2) (cid:6) (cid:6) wesmooththeimageswitha7×7Gaussianfilter(standard exp − Ga(t)(cid:6)px+t−pycr+pt(cid:6)2/h2 ·pycrp wdMeeAvbTisaLittAieoB2n1 =itimto3po)l.kemTlehesinsstatythiiaoenlnd1s ommbiotnariefnoesrdtaobnlfeerorrmeegsuisltttrhsa.etiUonsaiunotgpheotrrha’es- px(cid:2) = y∈Wx(cid:2) exp(cid:3)t∈−P(cid:2)Ga(t)(cid:6)(cid:6)px+t−pycr+pt(cid:6)(cid:6)2/h2(cid:4) . (3) tionoftwo2562CTscansonaquad-coreDellXPS2.66GHz y∈Wx t∈P PCwith8GBofmemory.Intheirpaper,Liuetal.11alsopoint Here, the superscript crp indicates that the pixels originate out that a GPU implementation of their belief propagation fromthecleanregisteredpriorandnotfromthetarget.How- algorithmcouldyieldafurther(upto)50-timespeedupwhich ever,whencomparingthisresultwiththatobtainedatregular would bring the time required for the registration down to dose[Fig.5(f)],westillobservesomeamountofblurringin seconds. the image. Edges in general appear less defined, and small featuresarealsoweakenedorcompletelysuppressed.Forthe latter,compareforexampletheintricatedetailinthecenterof II.C. Reference-basedNLM(R-NLM)filtering theimage,totheleftofthepincushion-shapeddarkstructure, With the NLM-mechanism still being employed for whichisbarelyvisibleintheR-NLMresult.InSec.II.D,we the matching, we call our approach reference-based NLM describeanadvancedschemethatovercomestheseproblems. (R-NLM)filteringsinceitusesthepriorimageasareference Finally,intheeventthatnoreliableupdatecanbefoundin toguidethefiltering.Figure4providesanillustrationofthis thepriorforagiventargetpixel,wefallbacktoconventional process, now using a noticeable smaller search window than NLM-filteringusingthetargetimage.Weidentifythissitua- in the regular NLM-case. The larger the search window the tionbyalowsumofweightsinEq.(3).Inourexperiments, MedicalPhysics,Vol.39,No.8,August2012 4752 W.XuandK.Mueller:Efficientlow-doseCTartifactmitigation 4752 (a) low-dose FBP reconstruction (b) prior (c) prior registered to (a) (d) R-NLM filtering of (a) (e) MR-NLM filtering of (a) (f) regular-dose reconstruction for comparison (g) registered prior (c) with (h) deformation field, acting on (f), simulated artifacts of (a) to determine prior (b) FIG.5. Results-supportedillustrationofreference-basedNLM(R-NLM)andmatchedreference-basedNLM(MR-NLM). wehaveusedathresholdof0.001.Weusethiscriterionboth parameterhforscaling,itisaglobalparameterthatscalesall forR-NLMandfortheadvancedschemedescribednext. patchesatthesameweight.Thedifficultiesweencounterwith R-NLM cannot be solved just by adjusting the factor h, as wewilldemonstrateinSec.II.E.Justasthebestregistration II.D. Matchedreference-basedNLM(MR-NLM)filtering is achieved when the two scenes are similar in appearance, Since the features are generally weakened at a scale less thebestNLM-matchisobtainedwhentargetandpriorhavea thanthesizeoftheNLMsearchwindow,wecannotblamethe similarappearance.Thisisnotthecasewhenpairingaclean registrationalgorithmfortheseshortcomings.Rather,itisthe prior and a low-dose reconstruction with severe streak arti- qualityoftheNLM-matchingthatisattheheartoftheprob- facts. Hence, we require a method that transforms the clean lem.ConsidertheNLM-distancefunctionofEqs.(1)and(3) registered prior image into an image that bears similar arti- usedtodeterminethequalityofamatchforaspecificcandi- facts as the target image. We can achieve this by first sim- dateneighborhood(orpatch)P, ulating projections from the registered prior and then recon- (cid:5) structingitunderthesameconditionsasthetargetimage,i.e., Ga(t)|px+t −py+t|2. (4) with a lower number of projections and at the same viewing t∈P geometryasthetarget.ThisgivesrisetoFig.5(g)—whichas Here, |px+t – py+t| is the Euclidian distance of a correspond- we observe looks fairly close to the target image subject to ing pair of pixels parameterized by patch index t. The sum denoising[seeFig.5(a)]. of these distances determines the weight that the patch P The MR-NLM reconstruction procedure is illustrated in plays in determining the value of the target pixel, and thus formofpseudocodeinFig.6.Afterregisteringthepriorwith it is the patch’s structural similarity that is decisive for the thelow-dosetarget,adegradedregisteredprioriscreatedby scalingofitscontribution.WhileEqs.(1)and(3)alsohavea simulatingthelow-doseartifactalsopresentinthetarget.The MedicalPhysics,Vol.39,No.8,August2012 4753 W.XuandK.Mueller:Efficientlow-doseCTartifactmitigation 4753 Input: sharper,smallfeaturesarebettervisible,andwealsoseethat Low-dose scan L, normal-dose prior scan N, low-dose degradation D; theintricatedetailinthecenteroftheimage,totheleftofthe Preprocessing: pincushion-shapeddarkstructureisalsoclearlyrestored. 1. Register N to L using SIFT-flow and obtain N: R N ← SIFT_Flow_Registration (N, L); R 2. Generate projections P of N with the same low-dose degradation D as the input: R P← Forward_Projection(N, D); II.E. ComparingNLM,R-NLM,andMR-NLM R 3. Generate the degraded version of N – N – using FBP: N ← FBP(P, D); R DR Figure 7 compares the three different schemes, using a DR Filtering: case study at “microscopic” detail. In Fig. 7(a), we show a Apply MR-NLM filtering to L with <N, N > and return the denoised filtered result L: clean registered prior—the reference. For the matter of this R DR F LF← MR-NLM (NR, NDR, L); discussion, we shall focus on the 7 × 7 image cutout— equivalent to a NLM search window—within the black box FIG.6. Pseudocodeofthematchedreference-basedNLM(MR-NLM). inthelowerrighthalfofthisimage.Thiscutoutshowsapor- tion of a bony structure. In Fig. 7(b), we illustrate the data procedurethenusesthisdegradedregisteredpriorforNLM- flow and operations of all three schemes using the cutout as matching, but copies the corresponding candidate pixels in an example. In this schematic, the cutout labeled “clean ref- thecleanregisteredpriortotheweightedsum.Theresulting erence” is a copy of the black-boxed region, while the “tar- equationis,modifyingEq.(3), (cid:3) (cid:4) get” is the corresponding cutout in the low-dose scan which (cid:2) (cid:2) (cid:6) (cid:6) exp − Ga(t)(cid:6)px+t−pyd+rpt(cid:6)2/h2 ·pycrp issimmuulcahtedsetghreadleodw.-Adossdeisacrutsifsaecdt,sthineMthRe-cNleLaMn rperfoecreednucreepfirros-t px(cid:2) = y∈Wx(cid:2) (cid:3)t∈P(cid:2) (cid:6) (cid:6) (cid:4) . (5) ducing the “degraded reference.” It then uses this image for exp − Ga(t)(cid:6)px+t−pyd+rpt(cid:6)2/h2 NLM-matching, but retains the corresponding pixels in the y∈Wx t∈P clean reference to update the target, yielding the cutout la- Here the subscript crp denotes the clean registered prior, as beled “MR-NLM.” On the other hand, the R-NLM proce- before, while the subscript drp denotes the degraded regis- dure(dottedlines)usesthecleanreferenceformatchingand teredprior. updates the target directly, giving rise to the cutout labeled Figure5(e)showstheresultweobtainedforourtestexam- “R-NLM.”Finally,theconventionalNLMprocedureusesthe ple. We observe that the edges are now overall significantly targetforbothmatchandupdate,producingthecutoutlabeled FIG. 7. ComparingtheNLM,R-NLM,andMR-NLMfilteringschemesintermsoftheireffectona7×7imageregionequivalenttothesizeofasearch window.(a)Registeredcleanpriorwithblackboxinlowerrightimageregionindicatingtheregionstudied.(b)Thethreepipelinesillustratedbyexample. (c)Themapofpixelcontributionsforthissearchwindow.Eachpixelisassociatedwitha7×7patchcenteredonit. MedicalPhysics,Vol.39,No.8,August2012 4754 W.XuandK.Mueller:Efficientlow-doseCTartifactmitigation 4754 “NLM.”Therowofresultcutoutsdemonstratesanincreasing imentstheprioriscreatedbynonlinearlydistortingtheorigi- growth in quality from left to right. While the NLM cutout nalscan,thelow-doseimageiscreatedinalignmentwiththe is quite similar to the low-dose target subject to denoising, original scan, and the registration brings the prior back into theR-NLMcutouthassomewhatsharperdetail,inparticular approximatealignmentwiththeoriginalscan.Thus,themost inthecenter.Finally,theMR-NLMcutouthasthemostpro- appropriategoldstandardistheoriginallyobtainedscan).The nouncedsharpness—notquiteasstrongasthecleanreference μ andμ aretheaveragesofthelowandregulardoseimages l r butfairlyclose. landr,respectively,andNisthetotalnumberofpixels. Further insight can be obtained from visualizing the dis- Theadvantageofthesemetricsisthattheyareeasytocom- tance map for the NLM search window coinciding with the puteandhaveclearphysicalmeanings.However,theyreveal studiedimagecutout.Notethatthissearchwindowwillonly only little about the perceptual impact certain image differ- resolvethevalueforthepixelinthecenterofthecutout.This ences may have. The RMS metric computes the pointwise distance map is used for the matching—see Eq. (4). Plotted errors and pools them across the entire image—this ignores in Fig. 7(c) are the corresponding maps for the clean refer- anyspatialcoherenceandsocannotgaugethedifferencesin ence used in R-NLM and for the degraded reference used structure and contrast that may exist in local pixel neighbor- in MR-NLM, respectively. We can easily see that the dis- hoods.Ontheotherhand,whileCCdoesprovideastatistical tances in the latter map are much closer than those in the measureofimagedifferences,itcomputesitataglobalscale former which confirms the better correspondence. The third and it also considers only pixel intensities which are far less map, labeled “degraded/clean” shows the ratio of the two perceptuallysalientthanlocalcontrastsandedges. maps. We clearly see that this ratio is not constant across Asanattempttobetteraccountforhumanperceptionwhen the patch and thus a simple boosting of the h-parameter determining image quality, we have employed metrics that in the NLM equations would not be able to rectify this specificallygaugethepreservationofperceptuallysalientin- situation. formation, which we define as image content to which the As for the efficiency of the three algorithms, the compu- humanvisualsystemismostsensitiveto. tational complexity is O(NWP) where N, W, and P repre- The first such metric is E-CC, as defined in our earlier sent image size, search window size and patch size, respec- work.24 ItisidenticaltoCCbutoperatesontheedge-filtered tively. However, in a GPU implementation, due to the high images which we obtain using a Sobel mask. E-CC is still a pixel independence and therefore potential parallelism, the globaloperatorbutitconsidersmoreperceptuallysalientlow- speed could be greatly increased. As determined in Ref. 22 levelimagefeatures,i.e.,edgesthatdefinetheboundariesof forthe2Dcase,itonlytakes18mstodenoisea5122 image thereconstructedobjects. witha112 searchwindowanda72 patchsize.MR-NLM(R- Anothermetricweemployisthestructuralsimilarityindex NLM) incurs a small additional overhead for reading from (SSIM) devised by Wang et al.25 The SSIM is an enhance- two (one) other image(s) and for checking if falling back mentoftheuniversalimagequalityindex(UQI)(Ref.26)also to conventional NLM-filtering is more appropriate. Over- recently used by Bian et al.19 Both UQI and SSIM combine all, the entire algorithm, including FBP reconstruction (see thedifferencesinmeanintensity,contrast,andstructureinto Ref. 23), SIFT-flow registration (see Sec. II.B), and MR- asinglequalityfigure.TheSSIMiscomputedforeachimage NLM or R-NLM filtering, would most likely take on the pixel at position x over a sliding small image window—we j order of seconds when accelerated on a high-performance usean11×11mask—andthencombinedintoapooledindex GPU. SSIM byaveragingtheindividualSSIMmeasurements pooled (cid:3) (cid:4)(cid:3) (cid:4)(cid:3) (cid:4) 2μμ +c 2σσ +c σ +c SSIM= l r 1 l r 2 lr 3 II.F. Assessingimagequality μ2l +μ2r +c1 σl2+σr2+c2 σlσr +c3 To evaluate thequality performance of thevarious recon- 1 (cid:5)N struction schemes we have employed two groups of metrics. SSIMpooled = SSIM(xj). (7) N ThefirstgroupencompassesthetraditionalRMS(rootmean j=1 square)andCC(correlationcoefficient) measuresdefinedas Here,thesubscriptslandrdenotethelow-doseandregular- follows: dose images, respectively, and the μ and μ are the means (cid:7) l r (cid:2) of the pixels within these corresponding windows, while the N (p −p )2 RMS= i=1 l,i r,i σl and σr are their standard deviations and the σlr is their N covariance. The constants c , c , c are typically small (see (cid:2) 1 2 3 N (p −μ)(p −μ ) Ref. 25) and prevent numerical instabilities when a denom- CC= (cid:8)(cid:2) i=1 l,i (cid:2)l r,i r . (6) inator is close to zero—the UQI does not have these con- Ni=1(pl,i −μl)2 Ni=1(pr,i −μr)2 stantswhichcanleadtowrongestimateswhentheseadverse conditionsaremet.Finally,toavoidblockingartifactsWang In these metrics, the p are the pixels in the low-dose re- etal.recommendaGaussian-weightingofthesamplesunder l,i construction and the p are the pixels in the corresponding aSSIMwindow.SinceSSIMisthegenerallyacceptedname r,i regular-dose reconstruction, in our case constituted by the ofthemetric,wewilluseitthroughoutthepaperbutitisun- originally obtained scan image. (We note that in our exper- derstoodthatweuseitspooledversion. MedicalPhysics,Vol.39,No.8,August2012 4755 W.XuandK.Mueller:Efficientlow-doseCTartifactmitigation 4755 The first term in Eq. (7) is quite consistent with the just- TABLEI. Numericalcomparisonoftheresultsobtainedfortheheadphan- noticeable intensity difference (JND) metric often used in tomviavariousmetrics.Thepercentagefigureforametricmeasurestheim- provementwithrespecttothemethodtoitsimmediateleft.Totheleftof perceptualqualitystudies.Thesecondtermcomparesthelo- the%cell,abovethescores,welisttheoptimalparametersettingforeach calcontraststhatexistintheslidingwindow.Finally,thethird algorithmwhichweobtainedbymanualtuning. termevaluatesthestructuralsimilarityafterthedifferencesin means and contrasts have been accounted for. The SSIM is W/O TVM NLM R-NLM MR-NLM quitepowerful—studiesthataskhumanobserverstorankim- N/A λ=30 (%) h=220 (%) h=200 (%) h=120 (%) ages with identical scenes, but corrupted with different arti- facts,intermsofqualityshowthattheserankingcorrelateex- RMS 165 123 25.4 126 2.4 57 54.8 45 21.0 ceedinglywellwiththeSSIMoutcome.Furthermore,itisalso CC 0.96 0.98 2.1 0.98 0 0.99 1.0 0.99 0 interestingthatinthesestudiesallimageshadthesameRMS E-CC 0.62 0.73 17.7 0.74 1.4 0.94 27.0 0.96 2.1 error. Finally, large experiments27 have shown that SSIM is SSIM 0.43 0.53 23.2 0.53 0 0.95 79.2 0.97 2.1 particularlywellsuitedtodetectdistortionscausedbynoise. Italsodetectsspatiallycorrelatednoise,whichinCTimages couldmimicfalsefeatures. tunedforoptimalperformance.Thefirstobservationwemake isthatallmetricsshowsimilartrends(butwealsoobservethat CCismuchlesssensitivetothechangesinimagequality).In III. RESULTS general, for the CC, E-CC, and SSIM the maximum possi- We have run the algorithms described above on two ble value is 1.0, while for the RMS error the optimal value datasets:aheadphantomandahumanlung.Theheadphan- is0.Inthisparticularexperiment,allmetricsreachtheirbest tom is part of a body phantom scanned with a GE Light- valuesfortheMR-NLMalgorithm—around0.97fortheper- Speed scanner. The human lung scan was obtained from the ceptual metrics—and their worst values when no filtering is “Give a Scan” dataset collection.31 Specifically we used the applied.ItalsoappearsthatNLMandTVMreachquitesim- first dataset series 2 of patient p0015 obtained with a GE ilarscoresforeachmetric,withaslightadvantageforNLM. LightSpeed16 scanner. For all examples, we used the orig- This can be verified by comparing the images (see Fig. 2) inal floating-point reconstructions for three purposes. First, which look fairly similar. The improvement of R-NLM over theyservedasthebasisforahigh-qualityprojectionsimula- NLM is significant for both RMS (54%) and the perceptual tioninfan-beamgeometry(fanangle=20◦).Wethenpicked metrics (27% for E-CC and 79% for SSIM). The improve- a subset of these projections and reconstructed the reduced- ment for MR-NLM over R-NLM is another 21% for RMS, projections low-dose imagery studied in this paper. Second, and2%forthetwoperceptualmetrics. wealsousedthemtogeneratethepriors.Forthis,weapplied The absolute difference images presented in Fig. 9 show variousnonlineardistortionsonthemandsubsequentlyregis- similar trends. Marked improvements can be observed for teredthemtothelow-dosereconstructions.Third,theyrepre- R-NLM over the prior-less schemes, and more moderately sentedthegoldstandardforallnumericalqualityassessment for MR-NLM over R-NLM. The improvement achieved by viathevariousmetricsdescribedinSec.II.F. R-NLMismainlyinstreakremoval.MR-NLM,ontheother We begin with the head phantom already examined in hand, adds sharpness and detail definition which can be ap- Sec. II to illustrate the outcomes of the various algorithms. preciated by the overall much smaller errors, in particular at TableIandFig.8comparetheresultsobtainedforMR-NLM, edgesandsharpcorners. R-NLM,NLM,TVM,andnofiltering,asgaugedbytheRMS, In order to better explore the performance of the various CC,E-CC,andSSIMerrormetrics.TableIalsogivestheset- algorithmsatthelocaldetaillevelwehaveconductedaROI- tingsforthevariousalgorithmparameterswhichwemanually basedanalysis(seeFig.10).Figure10(a)depictsthelocations of four ROI regions and Fig. 10(b) plots the corresponding SSIMscores.Weobservea7%–8%improvementforR-NLM overNLMinallfourROIsandanother7%–8%inROI2,3, and 4 for MR-NLM over R-NLM. ROI 1 contains a fairly structured and high-contrast feature which is less in need of TABLEII. Numericalcomparisonoftheresultsobtainedforthehumanlung viavariousmetrics. W/O TVM NLM R-NLM MR-NLM N/A λ=20 (%) h=260 (%) h=180 (%) h=120 (%) RMS 222 154 30.6 153 0.6 128 16.3 120 6.2 CC 0.92 0.96 4.3 0.96 0 0.97 1.0 0.98 1.0 E-CC 0.6 0.69 15.0 0.71 2.9 0.81 14.1 0.85 4.9 FIG.8. Graphicalcomparisonoftheresultsobtainedfortheheadphantom SSIM 0.49 0.61 24.5 0.61 0 0.66 8.2 0.69 4.6 viavariouserrormetrics. MedicalPhysics,Vol.39,No.8,August2012 4756 W.XuandK.Mueller:Efficientlow-doseCTartifactmitigation 4756 (a) (b) (c) (d) (e) FIG.9. Absolutedifferenceimagesfortheheadphantomreconstructions.(a)MR-NLM;(b)R-NLM;(c)NLM;(d)TVM;and(e)withoutfiltering. the added fidelity of the MR-NLM scheme. Lastly, Fig. 11 plots overall evaluations with the various metrics. We make presents the ROI study visually in form of cutout details, similar observations as for the head phantoms but note that which best demonstrate the considerable benefits the MR- in this test case the quantitative improvements for R-NLM NLMrestorationmethodprovides.WeobservethatallROIs andMR-NLMaremorebalancedandtheMR-NLM/R-NLM showsignificantlymoredetailforMR-NLM,asopposedtoR- gainisaboutdoublethanthatfortheheadphantom.Finally, NLM.Infact,thereconstructionsarequiteclosetotheideal Fig. 15 depicts ROI-definitions and the SSIM-scores for the image. On the other hand, the differences of R-NLM vs the studied restoration schemes. Figure 16 shows the cutout de- prior-less methods are also significant, but not as marked as tails.Again,weseethatMR-NLMsignificantlyimprovesthe forMR-NLMvsR-NLM. fidelity of small detail and in fact it is even able to restore Next,Fig.12showthesamesequenceofresultsobtained someoftheoriginalCTimagenoisetexturethatwaspartof withthehumanlung,reconstructedfrom90projectionsover theprior. 360◦. The distortion applied was a fisheye warp. Figure 13 Finally, Fig. 17 explores the effect of different search presents difference images, and Table II lists and Fig. 14 window sizes, for the human lung. The size needed mainly 1 0.95 ROI1 ROI4 0.9 0.85 M ROI3 SI 0.8 S ROI2 0.75 0.7 ROI 1 ROI 2 ROI 3 ROI 4 0.65 0.6 W/O TVM NLM R-NLM MR-NLM Filters (a) (b) FIG.10. ROI-basedanalysisforonesliceoftheheadphantom.(a)ROIlocationsand(b)SSIMevaluation. MedicalPhysics,Vol.39,No.8,August2012 4757 W.XuandK.Mueller:Efficientlow-doseCTartifactmitigation 4757 1 OI R 2 OI R 3 OI R 4 OI R W/O TVM NLM R-NLM MR-NLM Ideal FIG.11. ComparingthefourROIsdefinedinFig.9visually. (a) prior (b) registered prior (c) degraded registered prior (d) low-dose reconstruction (e) TVM-filtered (f) NLM-filtered (g) R-NLM-filtered (h) MR-NLM-filtered (i) ideal phantom FIG.12. Resultsobtainedforthehumanlung. MedicalPhysics,Vol.39,No.8,August2012