Gibbs-Ringing Artifact Removal Based on Local Subvoxel-shifts EliasKellner1,BibekDhital1,MarcoReisert1, 5 Correspondingauthor: EliasKellner 1 0 2 1DepartmentofRadiology,MedicalPhysics,UniversityMedicalCenterFreiburg,Germany, n a J 0 Abstract 3 where I(x) denotes the image value at voxel index x andc(k)aretheN expansioncoefficients. Inpractice, ] h Gibbs-ringing is a well known artifact which mani- only a finite number of expansion coefficients can be p acquired. This truncation of Fourier space introduces fests itself as spurious oscillations in the vicinity of - artifactsiftheexpansioncoefficientsdonotdecayfast d sharp image transients, e.g. at tissue boundaries. enoughwithincreasingk(1). Thisisthecaseforsharp e The origin can be seen in the truncation of k-space m image transitions, where all higher frequency compo- during MRI data-acquisition. Consequently, cor- nentsarerequiredtoproperlyreconstructtheedge. The . rection techniques like Gegenbauer reconstruction s artifact can nicely be illustrated when reconstructing c or extrapolation methods aim at recovering these si missingdata. Here,wepresentasimpleandrobust the image on a fine grid using zero-filling by setting the missing high frequency components for k > N to y methodwhichexploitsadifferentviewontheGibbs- h zero. Thisoperationcorrespondstoamultiplicationof phenomena. Thetruncationink-spacecanbeinter- p the‘true’k-spacewitharectangularwindow, whichis pretedasaconvolutionwithasinc-functioninimage [ equivalent to a convolution with a sinc function in the space.Hence,theseverityoftheartifactsdependson 1 howthesinc-functionissampled. Weproposetore- imagedomain.Thesidelobesofthesincresultinoscil- v lations (‘ringing’) in the neighborhood of sharp edges. interpolatetheimagebasedonlocal,subvoxelshifts 8 TheissueisillustratedinFig.1fortheimageofasingle to sample the ringing pattern at the zero-crossings 5 edge. 7 of the oscillating sinc-function. With this, the arti- 7 fact can effectively and robustly be removed with a 0 minimalamountofsmoothing. Several approaches have been proposed to reduce . ringingartifacts. Moststraightforwardly,imagefilters, 1 Key words: Gibbs-Ringing — Ringing-Artifact — e.g.theLanczosσ-approximation(2,3)oramedianfil- 0 Truncation-Artifact 5 tercanbeusedtosmooththeoscillations. However,as 1 filtering implies global blurring, the spatial resolution : is effectively reduced, leading to a loss of fine image v i Introduction details. More advanced methods have been developed X basedonpiecewisere-reconstructionofsmoothregions r using Gegenbauer-Polynomials (2, 4). A drawback of a In MRI, images are not gained directly, but with a re- these methods is the requirement of an edge detection constructionfromacquisitionsofcorrespondingFourier and potential instabilities in some applications (5) due expansioncoefficientsink-space, totheinvolvedchoiceofparameters. Anothermodern approach consists of combined de-noising and ringing N−1 1 (cid:88) −i2πkx removalusingusingtotalvariationconstraineddataex- I(x)= c(k)·e N , (1) N ploration(6). k=0 1 absolutedifferencesintheneighborhood. Itseemsad- A B visabletomeasuretheringingtobothsidesofthepixel individually,andselectthesmallervalue. Thisway,in theneighborhoodofanedge,theringingisalwaysmea- sured in the direction opposite to the edge (see Fig.2. Thus, the step itself does not contribute, and interfer- ences from closely located edges are minimized. We Figure1: Animagewithasinglediscontinuity(blackedge)isrecon- structedfromtruncatedk-spacedata. Theresultingimage(bluedots) decided to use a simple absolute differences approach isdiscreteandexhibitsringingartifacts. Theamplitudeoftheringing tomeasuretheoscillations. dependsonwhethertheunderlyingsincpatternarisingfromthewin- dowingink-space(redcurve)issampledatitsextrema(A),oratthe zero-crossings(B). (cid:88)k2 TV+(x) = |I (x+n)−I (x+(n−1))| s s s n=k1 The method we propose in this work is based on a (cid:88)k2 different view on the effect: In fact, due to the finite TV−(x) = |I (x−n)−I (x−(n−1))| s s s number of expansion coefficients, the image is not re- n=k1 constructed on a continuous domain, but on a discrete where K = [k ,k ] defines the window size which is grid. Obviously, the strength of the ringing in the re- 1 2 used for measuring the amount of oscillations. There constructed image depends on the precise location of are reasons to exclude the central point itself, using a theedgerelativetothesamplinggrid,i.e. howthesinc- window of the form K = [1,k ]. This leads to more functionissampled. Ifthesidelobesofthesinc-pattern 2 stabilityforpointsdirectlyontheedge,andminimizes are sampled at its extrema, the ringing amplitude be- increaseinnoisecorrelation,astheoscillationmeasure- comes maximal, whereas it disappears when sampled mentitselfisnotcorrelatedwiththeactualpointweare atthezerocrossings(seeFig.1). Thisfeaturehasbeen interestedin. ThechoiceofthewindowK constitutes demonstratedinexperimentaldatainthecontextofthe theonlyparametersofthemethod. Theresultswillin- darkrimartifactincardiovascularimaging(7). Finding dicate that the proposed algorithm is relatively robust the optimal subvoxel-shift for pixels in the neighbor- againstthischoice. hood of sharp edges in the image can therefore min- Now,foreachpointxtheoptimalshiftisdetermined imize the oscillations. However, as there are multi- byfindingtheshiftssuchthatthetotalvariationTV is ple edges present in an image, the correction cannot minimal. First,independentlyfortheright(+)andleft be achieved by a single, global shift, but must be per- (-)side: formedonalocalbasis.Inthenextsection,wedescribe anon-iterativemethodwhichconductsthistask. r+(x) = argminTV+(x) (3) s s r−(x) = argminTV−(x) (4) Methods s s Finally, we decide whether the overall minimum One-dimensionalCase LetI (x)betheoriginal,dis- 0 min(TV+ (x),TV− (x))comesfromleftorright crete image, and c0(k) its Fourier expansion coeffi- r+(x) r−(x) cients. From this, a set of 2M images I (x), where andthisshiftisdeterminedtobetheoptimalshiftr(x). s s=−M...M −1,withsubsequentsubvoxel-shiftsis Now, we know the optimal shift at each image loca- tion, which minimizes TV and hence also the ringing createdbymultiplicationwithphase-rampsinFourier- artifact. But, finally we have to go back to original space: grid, i.e. evaluating I (x(cid:48)) at the non-integer posi- r(x) N−1 tion x(cid:48) = x−r(x)/(2M). In order to do so we can Is(x)= N1 (cid:88) c0(k)·e−Ni2πk(x+2Ms ) (2) useanykindofalternativeinterpolationwhichisnota k=0 sinc-interpolation. Formally,wehave From this dataset, for each pixel x, the optimal shift I (x):=I (x−r(x)/(2M)) unring r(x) whichminimizespotentialoscillationsintheneighbor- hoodisfound. Thecorrespondingmeasurecanbecal- We decided to do this final ’back interpolation’ by a culated with any oscillation-sensitive kernel, e.g. the simple linear interpolation. In fact, one could also use 2 alternativehigh-orderinterpolationschemes(likepoly- canbeachievedwiththeLanczossigma-approximation nomial interpolation), but they usually lead to ringing (2,3),wherethefilteredimageisgivenby artifactsagain. 1 N(cid:88)−1(cid:18) k (cid:19)p −i2πkx I(x)= sinc ·c(k)·e N , (7) Two-dimensionalCase Theextensiontothetwodi- N N k=0 mensionalcaseisnotstraightforwardasdiagonaledges produce checkerboard-like ringing patterns, as can be where the parameter p controls the filter strength. In observed in the phantom image in Fig.4). Hence, it this study, we use two different settings for p. On one is not possible to find the optimal shift in both dimen- hand, weapplythestandardchoiceofp = 1. Thisin- sionssimultaneously. Asasolution,wecorrecttheim- duces, however, a rather strong smoothing of the im- ageinbothdimensionsseparatelyresultingintwoone- age. Hence, for a sound comparison with the pro- dimensionally corrected images J and J . These are posed method, we further set p such that the increase x y thencombinedinFourierspacetothefinalimageJ via in the correlation of the noise are the same for both, the proposed method and the filtering approach. The (cid:110) (cid:111) J =FT−1 FT{J }·G +FT{J }·G (5) corresponding reference noise correlation was mea- x x y y sured in a pure-noise region on the mean-free image where FT{·} denotes the Fourier transform. We pro- I˜(x)=I(x)−I¯via posetousethe‘weightingfunctions’G andG witha x y (cid:80) I˜(x)I˜(x−1) saddle-likestructureinFourierSpaceoftheform r = x (8) (cid:80) I˜2(x) x G = 1+cosky x (1+cosky)+(1+coskx) Another popular, non-linear filtering approach which (6) preservesedgesisgivenbythemedianfilter. Forcom- Gy = (1+cos1k+y)c+os(1k+xcoskx) parison, we also applied a median filter, where the ‘width’ofthefilterwasfixedtoa2x2neighborhood. Thisway,thehighfrequencycomponentsalongthedi- rection of the correction are enhanced, while they are dampened along the non-corrected direction. Due to Numerical Phantoms The method was applied to the normalization G +G = 1, artifact-free images, twonumericalphantoms(Fig.4). Forthefirstphantom, x y whereJ = J , areleftunchangedbyEq.6. Thisen- a polygonal shape with some stripes and small struc- x y suresthatminimalsmoothingisintroduced. tures was simulated. Starting with a high-resolution image,theartifactwassimulatedbyreconstructingthe image from a truncated k-space with a reduction fac- tor of 20. Second, for a more realistic brain phantom, data from a T -weighted post-contrast MRI measure- 1 ment of the brain was used. The artifact was artifi- cially enhanced by re-reconstructing the image from a smaller k-space with a reduction factor of 4. In both cases,a‘groundtruth’imagewithoutartifact,butwith the same decreased spatial resolution was generated byconvolvingthehigh-resolutionimagewithaboxcar function,andsamplingtheresultonthecorresponding Figure3: WeightingfunctionsGx(left)andGy(right)formultiplication low-resolution grid. Gaussian noise with a signal-to- withthe1D-correctedimagesinFourierSpace. Highfrequenciesare noiseratioof100wasadded. enhancedalongthedirectionofthecorrectionandsuppressedalong theotherdirection. MRI Measurements The method was applied to diffusion-weighted-images (DWI) with 70 directions, Reference methods We compare the proposed b = 1000s/mm2, usingagradientechoEPIsequence method against other, non-iterative filtering methods. with TE = 107ms, matrix size 104x104, resolution Inordertominimizeunwantedsmoothing,low-passfil- 2mm3, performed on a 3T scanner (Siemens TIM terwithahighcutoff-frequencyseempreferable. This TRIO, Siemens, Erlangen, Germany). No distortion 3 A) ∆ = 0 pixels B) ∆ = 0.25 pixels C) ∆ = 0.5 pixels D) ∆ = 0.75 pixels Figure2: 1-Dimagewithtwoedges. Asetofimageswithincreasingsubvoxel-shiftsiscreated. Theoptimalshiftisdifferentforthetwoedges. Fortherightedge,imageA)isoptimalwhereasitisC)fortheleftedge. Consequently,theoptimizationmustbeperformedonalocalbasis, i.e. foreachpixelindividually. Theoptimizationcriteriaisgivenbytheamountofringing,whichcanbecharacterizedbycalculatingtheabsolute differencesinacertainrangeforeachpixel.InC),thisisexemplaryillustratedfortheredpixel.Thetotalvariationiscalculatedseparatelytoboth, theleftandrightsideofthecenterpixel. Thisensuresthattheringingismeasuredaway fromtheedge,andtheedgeitselfisnotcontributing. Further,thecenterpixelitselfisexcludedinordertominimizetheintroducedcorrelationsbetweenneighboringpixels. Finally,thepixelvalueis recalculatedbyinterpolatingtheshiftedimageD)attheoriginalpixelposition,whichinthisexampleis0.5pixelstotheleftside(seereddotat bottom).Mathematicaldetailsaregiveninthetext. correctionwasapplied,astheinvolvedcorrectionmeth- of slightly reduced artifact removal. A kernel size of ods already lead to significant filtering of the artifact, K = [1,3]seemstobeaappropriatecompromisebe- especiallyinphasedirection.Duetothedifferentimage tween artifact removal and noise correlation. This set- contrast for different b-values, the artifact might even tingisusedfortheapplicationtotheMRIimages. be amplified during post-processing of the diffusion These findings are basically the same for the phan- parameters. Therefore, we also calculated diffusion tomconstructedfromtheT -weightedimageinFig.5. 1 maps, without and with artifact correction. We further Alsohere,theartifactcanmosteffectivelyberemoved applied the method to a T2-weighted image acquired usingtheproposedmethod,whilepreservingfineimage with a turbo-spin-echo sequence, TE = 109s, resolu- details. tion1x1x5mm2ona1.5Tscanner(SiemensSONATA, Siemens, Erlangen, Germany). Both dataset were ac- MRI images The results for DWI measurements are quiredinthecontextofclinicalroutine,writtenconsent shownforonesliceinFig.6.Apparently,theb -images wasobtainedtousethedataforscientificuse. 0 exhibitstrongringingartifacts,whichisevenmoreem- phasizedafterdiffusioncalculation. Theartifactcanbe Results reducedwithboth,themedianfilterandtheLanczosap- proximationwithp = 1,however,atthecostofstrong Numerical Phantoms The results of the phantom smoothing. With the proposed method on the other simulations are shown in Fig.4. We show results us- hand, the artifact can virtually completely be removed ing the median filter, the Lanczos approximation with with minimal filtering. Results from the T2-weighted p=1,resultsobtainedwiththeproposedmethodusing image given in Fig.7. The findings are basically the differentparameters,andresultsusingtheLanczosap- sameasfortheDWImeasurement. proximationwithfilterparameterspadaptedtoyieldan equalnoisecorrelationastheproposedmethod. Discussion Obviously,themedianfilterpreservestheedgesbet- ter than the Lanczos-approximation with p = 1, at a smallerincreaseintheaveragesmoothing,indicatedby Even though the Gibbs-ringing artifact is omnipresent the smaller increase in noise correlation. However, it in MRI, most vendors do not include removal tech- showsastrongerresidualoftheartifact,andfineimage niquesinthestandardimagereconstructingpipelineto details like small, peak-like structures are destroyed. date.Onereasonforthismightbethatstandardfiltering With the proposed method, on the other hand, the ar- approaches inevitably reduce the effective image reso- tifactcaneffectivelyberemovedwithminimalsmooth- lution,andmoreadvancedmethodslikeGegenbauerre- ing of edges. The method is rather robust against the reconstruction(4,2,8,9)ordataextrapolationmethods choice of the kernel parameter. A larger neighbor- (10, 11) are practically difficult to handle due to their hood results in less smoothing, but comes at the price complexity and, the requirement of an edge detection, 4 Figure4: Phantomimagewithmultipleedgesandnoise. Insertscorrespondtoascaledareawithpurenoise,correlation. Theamountofnoise correlation(givenbyr)reflectsthestrengthofthesmoothingintroducedbythemethodsandtheirparameters. Both,medianfilterandLanczos approximationwithp = 1leadtoratherstrongfiltering. Themedianfilteradditionallyintroducesartifactsonpoint-likestructures(seee.g. the blackdotsintheimages). Resultsusingtheproposedmethodaregiveninthesecondrow,theLanczosapproximationwithcorrespondingequal noisecorrelationinthebottomrow.Withtheproposedmethod,theartifactcanvirtuallycompletelyberemovedwithpreservationoftheedgesand withoutdestroyingfineimagedetails,whereastheLanczosapproximationwithequalnoisecorrelationstillsuffersfromsignificantartifacts. andpotentialinstabilitiesinducedbythedependencyon methodcanthereforeappliedtoanyimagewithauni- thechoiceoftheparameters. versal value for K. We found that K = [1,3] consti- Another approach consists of optimizations based tutesagoodcompromisebetweenartifactremovaland on total variation (6). These methods have proven smoothing. to effectively remove the artifact, however, they treat noise and artifact equally, in contrast to the proposed method, which explicitly aims at separating both con- Conclusions tributions. Further, the outcome of total variation ap- plicationsstronglydependsonthestrengthofthefilter, Inthiswork,wepresentedanon-iterativemethodforre- whichmustbeadaptedtotherespectiveapplication. movalofringingartifactsbasedonre-samplingtheim- The proposed method is rather robust to the choice agesuchthatthesourceoftheringingpattern,thesinc- of its parameter, the kernel width K. Further, this pa- function, is sampled at its zero crossings. We demon- rameterisindependentoftheimagesize,asoscillation stratedthatthemethodeffectivelyremovestheartifact pattern of the ringing occurs always in the distance of whileintroducingminimalsmoothing. Themethodhas one voxel, and hence scales with the matrix size. The a low computational cost, consists of a rather simple 5 Figure5: T1weightedpost-contrastMRIimageswithartificiallyenhancedartifact. Here,thesamecharacteristicsasinFig.4canbeobserved. Both,medianfilterandLanczosapproximationwithp = 1leadtoblurringofimagedetails(comparee.g. areasatredarrowsingroundtruth). Again,intheLanczosapproximationwithadaptedparameters,thereisstillresidualringingvisible(seee.g.bluearrows),eventhoughthedifference totheproposedmethodisherelesspronouncedcomparedtothephantominFig.4. References mathematicalframeworkandisverystableagainstthe choiceofitsfewparameters. Thisrobustnesssuggests itasasuitablecandidateforarobustimplementationin 1. Gottlieb D, Shu CW, Solomonoff A, Vandeven H. the standard image processing pipeline in clinical rou- OntheGibbsphenomenonI:recoveringexponen- tine. EventhoughdesignedinthecontextofMRI,the tialaccuracyfromtheFourierpartialsumofanon- methodmightalsoproveitsapplicabilityinotherareas periodic analytic function. Journal of Computa- suchasFourier-baseddatacompressionalgorithms. tional and Applied Mathematics. 1992;43(1):81– 98. 2. Gottlieb D, Shu CW. On the Gibbs phenomenon anditsresolution. SIAMreview.1997;39(4):644– Acknowledgements 668. 3. Jerri AJ. Lanczos-like σ-factors for reducing the We are grateful to Valerij G. Kiselev for very helpful Gibbs phenomenon in general orthogonal expan- discussions. sionsandotherrepresentations. JournalofCompu- 6 b=0images Radialdiffusivity Figure6: DWImeasurement. Theoriginalb = 0imageexhibitsstrongringingartifacts,whichareevenmoreamplifiedinthecalculatedradial diffusivitymaps(rightmaps).Theartifactcanmosteffectivelyberemovedwiththeproposedmethodataminimalimagesmoothing.TheLanczos approximationshowsresidualartifacts. tationalAnalysisandApplications.2000;2(2):111– 8. Shizgal BD, Jung JH. Towards the resolution of 127. the Gibbs phenomena. Journal of Computational andAppliedMathematics.2003;161(1):41–65. 4. Archibald R, Gelb A. A method to reduce the 9. FENG Qj, HUANG X, FENG Yq. A Fast Algo- GibbsringingartifactinMRIscanswhilekeeping rithm to Reduce Gibbs Ringing Artifact Based on tissueboundaryintegrity. MedicalImaging, IEEE theChebyshevPolynomials. JournalofImageand Transactionson.2002;21(4):305–319. Graphics.2006;8:013. 10. AmarturS,HaackeEM. Modifiediterativemodel 5. Boyd JP. Trouble with Gegenbauer reconstruc- based on data extrapolation method to reduce tionfordefeatingGibbsphenomenon: Rungephe- Gibbs ringing. Journal of Magnetic Resonance nomenoninthediagonallimitofGegenbauerpoly- Imaging.1991;1(3):307–317. nomialapproximations. JournalofComputational Physics.2005;204(1):253–264. 11. ConstableR,HenkelmanR. Dataextrapolationfor truncationartifactremoval. Magneticresonancein 6. BlockKT,UeckerM,FrahmJ.SuppressionofMRI medicine.1991;17(1):108–118. truncationartifactsusingtotalvariationconstrained dataextrapolation.Internationaljournalofbiomed- icalimaging.2008;2008. 7. Ferreira P, Gatehouse P, Kellman P, Bucciarelli- DucciC,FirminD. Variabilityofmyocardialper- fusion dark rim Gibbs artifacts due to sub-pixel shifts. Journal of Cardiovascular Magnetic Reso- nance.2009;11(1):1–10. 7 Figure7: T2-weightedimagewithobviousringingartifacts. Withtheproposedmethod,theartifactcaneffectivelyberemoved. TheLanczos approximationshowsresidualartifacts(seee.g.bluearrow). 8