Gibbs-Ringing Artifact Removal Based on Local Subvoxel-shifts

1 MB·English

by  Elias Kellner

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