Segmentation of artifacts and anatomy in CT metal artifact reduction

12 Pages·2012·1.63 MB·English

Segmentation of artifacts and anatomy in CT metal artifact reduction
SeemeenKarimia) andPamelaCosman
DepartmentofElectricalandComputerEngineering,UniversityofCalifornia,SanDiego,
LaJolla,California92093
ChristophWald
DepartmentofRadiology,LaheyClinic,Burlington,Massachusetts01805
HarryMartz
LawrenceLivermoreNationalLaboratories,Livermore,California94551
(Received25April2012;revised13August2012;acceptedforpublication14August2012;
published11September2012)
Purpose:Metalobjectspresentinx-raycomputedtomography(CT)scansareaccompaniedbyphysi-
calphenomenathatrenderCTprojectionsinconsistentwiththelinearassumptionmadeforanalytical
reconstruction.Theinconsistenciescreateartifactsinreconstructedimages.Metalartifactreduction
algorithmsreplacetheinconsistentprojectiondatapassingthroughmetalswithestimatesofthetrue
underlyingprojectiondata,butwhenthedataestimatesareinaccurate,secondaryartifactsaregener-
ated.Thesecondaryartifactsmaybeasunacceptableastheoriginalmetalartifacts;therefore,better
projectiondataestimationiscritical.Thisresearchusescomputervisiontechniquestocreatebetter
estimates of the underlying projection data using observations about the appearance and nature of
metalartifacts.
Methods:Theauthorsdevelopedamethodofestimatingunderlyingprojectiondatathroughtheuse
ofanintermediateimage,calledthepriorimage.Thismethodgeneratesthepriorimagebysegment-
ingregionsoftheoriginallyreconstructedimage,anddiscriminatingbetweenregionsthatarelikely
tobemetalartifactsandthosethatarelikelytorepresentanatomicalstructures.Regionsidentifiedas
metalartifactarereplacedwithaconstantsoft-tissuevalue,whilestructuressuchasboneorairpock-
etsarepreserved.Thispriorimageisreprojected(forwardprojected),andthereprojectionsguidethe
estimationoftheunderlyingprojectiondatausingpreviouslypublishedinterpolationtechniques.The
algorithm is tested on head CT test cases containing metal implants and compared against existing
methods.
Results:Usingthenewmethodofpriorimagegenerationontestimages,metalartifactswereelimi-
natedorreducedandfewersecondaryartifactswerepresentthanwithpreviousmethods.Theresults
apply even in the case of multiple metal objects, which is a challenging problem. The authors did
notobservesecondaryartifactsthatwerecomparabletoorworsethantheoriginalmetalartifacts,as
sometimesoccurredwiththeothermethods.Theaccuracyofthepriorwasfoundtobemorecritical
thantheparticularinterpolationmethod.
Conclusions: Metals produce predictable artifacts in CT images of the head. Using the
new method, metal artifacts can be discriminated from anatomy, and the discrimination can
be used to reduce metal artifacts. © 2012 American Association of Physicists in Medicine.
[http://dx.doi.org/10.1118/1.4749931]
Key words:metal artifact,computed tomography, image segmentation, sinogram in-painting, beam
hardening
I. INTRODUCTION centprogress,butcomparisonsacrossdifferentmethodshave
notbeenmadeotherthanwithlinearinterpolationacrossthe
When metal objects are present in x-ray computed tomog- metaltraces.1–13 Asaresult,thereisnorobustorwidelyac-
raphy (CT) scans, they are accompanied by bright and dark ceptedsolution,anditcontinuestobeachallengingresearch
shadowsandstreaks,collectivelycalledmetalartifacts.These problem.
artifactsareduetophysicalprocessesthatcausetheassump- The separation of artifact from real tissue is not a trivial
tion of linearity in the reconstruction to break down, so that task. The CT number (voxel intensity) ranges of metal arti-
the scanner projections cannot be accurately reconstructed facts and anatomical structures overlap, as do their gradient
using filtered backprojection or Radon transform inversion. ranges.Metalartifactreduction(MAR)algorithmsoperatein
Theartifactsobscureinformationaboutanatomicalstructures, projectionspace,sothatthescannerprojectionsamples,i.e.,
making it difficult for radiologists to correctly interpret the measuredray-sums,thatdonotpassthroughmetalcontribute
imagesorforcomputerprogramstoanalyzethem.Theprob- tothefinalimage,alongwithestimatedray-sumsthatreplace
lem has existed for many years,1,2 and there has been re- the measured ray-sums passing through the metal. In this
5857 Med.Phys.39(10),October2012 0094-2405/2012/39(10)/5857/12/$30.00 ©2012Am.Assoc.Phys.Med. 5857
5858 Karimietal.:SegmentationinCTMAR 5858
research,weimproveuponaclassofMARalgorithmsbased positionbasedmethodshaveshowngoodpreliminaryresults,
on sinogram in-painting.1–13 Sinogram in-painting methods sinogram in-painting is the most practical approach because
are vulnerable to the misrepresentation of anatomical edges itiscomputationallysimple,andneedsonlyoneenergyspec-
that overlap the metal traces. Various algorithms have tried trum.Itsmostseriouschallengeisaccuratedataestimation.
toaddressthisproblem.3–8,12,13Wesegmentregionscontami- In sinogram in-painting methods, the metal objects are
natedbymetalartifactsfromanatomicalstructuresintheorig- identified in the projection space. The regions in the projec-
inallyreconstructedimage,henceforthcalledtheoriginalim- tion space occupied by the metal are called metal traces. In
age.Voxelvaluesinregionsaffectedbyartifactsarereplaced somemethods,1–6metalobjectsarelocatedintheoriginalim-
withvaluesfromsurroundingregions.Thus,a“prior”image agebythresholding,andthetracesarelocatedbycalculation
isbuilt.Thepriorimagerepresentspriorknowledgethatcan or reprojection. We use the term scanner projections to de-
be used to guide ray-sum estimation within the metal trace. note log-attenuation projections measured from the scanner,
Asinothersinogramin-paintingmethods,themeasuredray- oralinearlyprocessedversionofthem,forexample,afterre-
sums outside the traces are left unmodified. The corrected binning to parallel projections. In other methods,7 the metal
projections are reconstructed to provide the artifact-reduced traces are located directly in the scanner projections. Early
image. work1,2interpolatedthescannerprojectiondataoneitherside
ofthetraces.Thismethod,whichwecalldirectinterpolation
(DI),missesedgesofrealstructuresasshowninFig.1.Inthe
II. BACKGROUND
picture,asofttissuebackgroundcontainsapieceofmetaland
abonystructure.Projectionsattwoanglesareshown.Inone
II.A. Causesofmetalartifacts
angle, the bone does not interfere with the metal projection,
Metalartifactsarecausedbybeamhardening (theprefer- andtheprojectioncanbeinterpolatedwithoutlossofedgein-
ential attenuation of low energy photons in a polyenergetic formation.Intheotherangle,interpolationofthedatawould
x-ray beam), photon scatter, partial volume effects, photon resultinanedgebeingblurredaway.Themissingedgeinfor-
starvation, and data sampling errors. Data sampling errors mationresultsinsecondaryartifactsinimages,whichmaybe
can be caused by inexact detector or view positions, cone assevereasthemetalartifacts.
beameffects,orpatientmotion.14 Beamhardeningandscat- It was proposed that edges be recorded from the original
tercauselow-frequencyartifacts,whichappearasbrightand image.3 The edges could be reprojected and combined with
dark shadows around the metal, obscuring the anatomy ac- the scanner projections to create a smoother projection on
tually present in the affected areas. The appearance of these which interpolation could be performed. However, a method
artifacts will be explained in more detail in Sec. III. The for recording edge information in reprojections was not de-
otherphenomenalistedabovecreatehighfrequencyartifacts, fined,butiscriticalintheperformanceofaMARalgorithm.
namely, streaks. Low frequency artifacts are difficult to re- It was proposed that the true values of image voxels be
move,buthighfrequencyartifactscanusuallyberemovedby estimated using k-means clustering and thresholding.5 The
nonlinearfiltering.15,16 reprojectionsoftheestimatedvoxelvalueswouldbeusedas
replacementdatainthemetaltrace.Thismethodwasreported
to fail when the artifacts were as bright as bone or as dark
II.B. Currentapproaches
asair,whichoccurredinalmostallofourtestimages.Other
Besides sinogram in-painting, other CT reconstruction existing methods built upon the notion of a prior knowledge
methodscanbeusedforMAR.Iterativereconstructionalgo- image that contains edge information.6–8 These methods
rithmshavebeenusedtoreconstructfromnoise-limiteddata. differinhowtheyutilizetheprior,usingvariousinterpolation
These can be adapted for MAR.17–19 Iterative reconstruction techniques.Onemethodusestheratioofscannerprojections
requires accurate modeling of the x-ray generation and at-
tenuation processes in the body and the detector, which is
difficult to accomplish. Another drawback of iterative meth-
odsisthattheyareslowbecausetheyrequirerepeatedrecon-
structions. Multiple-energy decomposition methods are used
todecomposematerialsintobasismaterials,andarelesssus-
ceptible to beam hardening artifacts.20,21 Therefore, energy
decomposition can be used to reduce metal artifacts.12,21,22
Energydecompositionmethodstakeintoaccounttheenergy-
dependent attenuation coefficients of different materials. As
aresult,theycancompensateforartifactsfrombeamharden-
ing.Twoormorex-rayspectraarerequiredforenergydecom-
position, but multienergy imaging is not standard in clinical
scanning protocols,nordoesitcompensate forscatter.Itera-
tive reconstruction of dual energy data21 has the potential to
provide excellent images if the dual spectra and models are
FIG.1. Aschematicofanobjectandtwoprojections.Directinterpolation
available. While iterative reconstruction and energy decom- resultsinmissededgesinoneprojection.
MedicalPhysics,Vol.39,No.10,October2012
5859 Karimietal.:SegmentationinCTMAR 5859
topriorreprojectionsinsteadofthedifference.6 Othermeth- 8000
odsrequirerepeatedreconstructions,eachofwhichimproves
e
tthhreesphroioldr.i7n,8gtHoopwroedvuecr,eathlletphreioserimmeatgheo,dasndreilnytenosnityintthernessihty- nit tim 6000
u
oldingleadstovoxelmisclassificationintheprior.Knowledge ea/
of the materials and accurate modeling helps retire some of nit ar 4000
this risk.12 Voxel misclassification leads to false edges or u
s/
missed edges and therefore to secondary artifacts. We call on 2000
the thresholded-prior methods TP-MAR.5–8 Some TP-MAR hot
P
methods suggest doing a first pass MAR with DI, and then
0
thresholdingthiscorrectedimagetocreateaprior.However, 30 60 90 120 150
thesecondaryartifactsfromDI-MARcanbeintenseenough Energy (keV)
thatthresholdingstillleadstovoxelmisclassification.
FIG.2. X-rayspectrumusedinbeamhardeningsimulations.
Weseethatwhilethenotionofthepriorimageisnotnew,
little attention has been paid to generating an accurate prior
image.Thesinogramin-paintingmethodswehavementioned photons(cid:2)atallenergiesE,
focus on the interpolation technique. Our experiments, de-
I = S(E)dE,
scribedinSecs.VandVI,showthattheaccuracyoftheprior 0
hasgreaterimpactonartifactreductionthantheparticularin-
andthepolyenergeticattenuationprocessisdescribedas
terpolation technique. There is much room for improvement (cid:2)
and for the use of image segmentation, in particular, in the I(s,θ)= S(E)e−∫s+lθ∈Lμ(s+lθ,E)dldE.
generation of a prior image. Strictly speaking, prior knowl-
edge consists of the observations we make about metal arti- In this case, the relationship between measured projec-
facts, but we use the term prior image for a combination of tionspandlog-attenuationnolongerapplies,leadingtobeam
theobservationsandmeasuredimage. hardeningartifacts.24 Compensationforbeamhardeningdue
toanatomicaltissuesisappliedduringscannercalibrationor
sometimesinpostprocessing.25
III. SIMULATIONSTUDYOFMETALARTIFACTS Wesimulatetwo-dimensional(2D)axialprojectionsofel-
CAUSEDBYBEAMHARDENING lipsoids in air of varying eccentricity, using a polyenergetic
spectrum.ThesimulatedspectrumisshowninFig.2.Itwas
Inordertosegmentareasofartifactandanatomy,wemust
obtained using the program XSPECW2.26 The tube voltage
understand the process of artifact generation. To help us un-
is 150 kVp, the beam filtration is 4 mm of aluminum, and
derstand the appearance of the metal artifacts, we simulate
thereisa10◦tungstenanodetarget.Thesimulatedmaterialis
beam hardening because it leads to low-frequency artifacts.
iron, whose energy-dependent attenuation cross sections we
Beam hardening artifacts are similar in appearance to scat-
obtainwithXCOM.27 Thereare1400projections,thedetec-
terartifacts,whichwehavenotsimulated.Subsequently,we
torspacingis0.5mmatisocenter,andthesource-to-isocenter
will make observations that we use to discriminate between
distanceis500mm.Theprojectionsarereconstructedwitha
regionsofartifactandanatomicalstructures.
Hanningfilterwithacutoffof10lp/cm.Ourimagevaluesare
The monoenergetic attenuation process can be described
clampedbetween−1024and15383.Amonoenergeticsimu-
usingtheBeer-Lambertlaw.23 IfI istheincidentnumberof
0 lationwasalsodoneforcomparison,withanominalCTvalue
photons with energy E , and μ(x,E ) is the linear attenua-
0 0 of30000HU(μcorrespondingto80keV).
tion coefficient at energy E for the material at some point
0 Whilesimulationsofmetalartifactshavebeenundertaken
representedbythevectorx,thenthenumberofphotonsatthe
before,14theshape-dependenceoftheartifactswasnotinves-
detectorisgivenby
tigated,norwastheaccuracyofreprojectionsthroughmetal.
I(s,θ)=I0e−∫s+lθ∈Lμ(s+lθ,Eo)dl. Our simulations, discussed in this section, demonstrate the
shape dependence of the beam hardening artifacts, explain
theirappearancewithrespecttoshape,andshowhowrepro-
In the above equation, the integration is done over the
jectionthroughthemetalcannotquantifythemetal.
scanned space L between source and detector, s is a vector
representingthesourceposition,andθ isaunitvectorinthe Figure 3 illustrates the shape dependence of the beam
hardening artifact. When the object cross section is circular,
directionfromthesourcepositiontoadetectorelement.
the beam hardening artifact does not exist outside the ob-
During reconstruction, we reconstruct the log attenuation
ject.However,beamhardeningisvisiblewithintheobjectin
projections:
(cid:2) (cid:3) (cid:4) the form of cupping. As the object’s eccentricity increases,
I the beam hardening artifact amplitude outside the object in-
p(s,θ)= μ(s+lθ,E )dl =ln 0 .
s+lθ∈L o I(s,θ) creases.Thedarkartifactisalongthelongaxisoftheellipse,
andthebrightartifactisalongtheshortaxis.
Now,ifwehaveapolyenergeticspectrumS,thentheinci- Figure 4 illustrates the formation of the bright and dark
dentnumberofphotonsisgivenbyintegratingthenumberof parts of the artifact. Ideal (monoenergetic) and hardened
MedicalPhysics,Vol.39,No.10,October2012
5860 Karimietal.:SegmentationinCTMAR 5860
ideal case, the undershoots of the filtered projection data in
geticgetic500/0500/0 aplrloajencgtleedsapcerrofesscttlhyeciommapgeen.sBautetfionrththeehmaredteanletrdacceasdea,ttahbeareckis-
erer==
onoenonoenW/WWW//LW/WL adrdoipsstolretsisonreolaftiavmeptolittuhdeeisd.eTalheinatmheplsihtuodrteeropfatthhe(c0e◦n)ttrhaalnrainy
mmWW thelongerpath(90◦).Theundershootsareevensmalleralong
thelongerpathsrelativetotheidealcase(duetomoreharden-
energeticenergeticWL=500/0WL = 500/0 iranerltgai)ftaitvcheta,lnyantlhdeestsshhenoelrgoteanrtgipveaerthpcsao.tmhTsphewisnisltaehtairvoeenlsa,ttlihveeaedlsyihnotgortoetormptuhacethhsbcrwoigmithh-t
olyolyW/WW//W/ pensation,leadingtothedarkartifact.
ppWW
Note that the metal attenuation cannot be quantified by
00
0000 simply reprojecting the metal voxels in a reconstructed im-
eticetic0/120/12 age.Thehardeningofthebeamcausesoverestimationofthe
gg00
olyenerolyenerWL = 60WL=60 rdeeprreosjteimctaiotinosniinntthheeddiirreeccttiioonnooffhloigwheessttaatttteennuuaattiioonn.,Faingduruen5-
ppW/W/WW// shows that the discrepancy between original hardened pro-
WW jectionsandreprojectionsincreaseswitheccentricity.Thisis
eelllliippssooiidd 11 eelllliippssooiidd 22 ellipsoid 3 because the metal object voxels are reconstructed using pro-
jectionswithdifferentamountsofhardeningineachview,but
FIG.3. Illustrationofthebeamhardeningartifactdependencewithobject thevoxelssummedineachreprojectedviewarethesamevox-
shape.Theellipseeccentricityincreasesfromlefttoright.Thetoprowshows els.Consideringthecentralrayalongthelongaxis,thevox-
anideal(monoenergetic)simulationandthebottomtworowsshowsimulated
elsalongitwerereconstructedusingrelativelyless-hardened
beamhardening.Forroundobjectsbeamhardeningisvisibleonlywithinthe
projectionsalongwiththemore-hardenedprojections.When
object.Thisisthewell-knowncuppingartifact.Asthemetalobjecteccen-
tricityincreases,theartifactsincreaseaswellandarevisibleoutsidetheob- compared with the most hardened projections, the reprojec-
ject.Thedarkartifactsarealongthelongaxis(maximumprojections).Top tionswillthereforebegreater.Asimilarreasoningexistsfor
tworows:WindowWidth(WW)/WindowLevel(WL)=500/0,bottomrow: every other ray through the metal. Due to this shape depen-
WW/WL=6000/12000HU.
dence, a simple beam hardening inversion cannot be per-
formedasisdoneforbonecorrectioninheadimagesinwhich
(polyenergetic) projection data that were filtered with the
theskullisroughlyring-shapedineachaxialslice.25
Ram-Lakkernelareshownnexttotheimage.Withagreater
We make three observations about metal artifacts in CT
ray length through the hardening material, there is a bigger
images, which are confirmed by the above simulations: the
discrepancy between hardened and ideal projections. In the
artifacts are adjacent to metal pieces, the amplitude of the
artifacts decreases as the distance from the metal increases,
90 degrees
1200 andthelocalmaximathroughmetalinprojectionspacecor-
respondtodarkartifactsintheimage.
1000
800
IV. METHOD
600
Ourmethodtobuildapriorimageoperatesontheoriginal
400
image. The goal is to discriminate artifacts from real struc-
200 tures in the original image so that we can replace artifact-
contaminated regions of the original image with tissue val-
-3 -2 -1 0 1 0 degrees ues, thus generating a prior. Our method for segmenting ar-
more drop 0.5 tifactsusesthethreeobservationsaboutmetalartifactsgiven
less drop
intheprevious paragraph. Wedescribethemainconcepts in
n
monochromatic atio 0 the paragraphs below, and then give details for each step in
nu subsections.
polychromatic e
att -0.5 Basedonthefirstobservation,localimageextremathatare
d
ere adjacenttothemetalvoxelsandsmallerthanagivenscaleare
filt -1 identified and interpreted as metal artifacts. The grouping of
voxels forartifactidentification (labeling)isdone viaregion
smaller undershoot
200 400 600 800 10001200 growing. The region growing uses a distance-based parame-
detector elements ter for voxel inclusion, which makes use of the second ob-
servation. The distance-based parameter limits the inclusion
FIG. 4. Illustration of the formation of the beam hardening artifact for a
of anatomy into the artifact labels. As will be discussed in
metalobjectwithellipticalcrosssection.Filteredprojectionsareshownat
detail below, these image regions are removed by replacing
twoanglesthatcorrespondtothelongestandshortestpathlengthsthrough
thecenteroftheobject. theirvoxelvalueswithsurroundingvoxelvalues.Alsobased
MedicalPhysics,Vol.39,No.10,October2012
5861 Karimietal.:SegmentationinCTMAR 5861
ellipsoid 1 ellipsoid 2 ellipsoid3
10
measured
reprojected
5
n
o
ati
u
n 0
e
att
g-
o 10
l
5
0
400 500 600 400 500 600 400 500 600
detectors
FIG. 5. ReprojectionsofthethresholdedironellipsoidsshowninFig.3(dottedlines)alongwiththeoriginalhardenedprojections(solidlines).Thetop
rowshowsreprojectionsalongtheshortaxisoftheellipsoids,andthebottomalongthelongaxis.Reprojectionsareunderestimatedalongtheshortaxisand
overestimatedalongthelongaxisoftheellipsoid.
onthesecondobservation,whenalocalmaximumlabelcon- lowedbyopening-by-reconstruction(OBR).28CBRisamor-
tains both artifact and bone, we use a discriminant curve to phological operation that performs grayscale dilation with
classify voxels within it as belonging either to artifact or to a structuring element followed by iterative erosion that is
bone.Werestorethebonevoxelstotheimage.Basedonthe constrained by the original image. Similarly, OBR first per-
third observation, when a local maximum through metal is formsgrayscaleerosionfollowedbyiterationsofdilationcon-
foundinprojectionspace,wecanexpect tofindlocalimage strained by the original image. CBR eliminates dark regions
minimacorrespondingtoit.Asaresult,wecaninterpretim- (localimageminima),andOBReliminatesbrightregions(lo-
ageminimaasartifactsevenwhentheyaresplitofffromthe calimagemaxima)thataresmallerthanthescaledetermined
metal. The dark artifacts can split off from the metal in the bythestructuringelement.Largerextremaorregionswithout
presence of high density CT objects, or multiple pieces of extrema are left alone. CBR and OBR, respectively, replace
metal. The approach, implemented in MATLAB version 7.10 voxelvalues intheclosedoropened regionswithvalues de-
(TheMathWorksInc.,Natick,MA),hasthefollowingsteps. rivedfromvoxelssurroundingtheseregions.Thestructuring
Figure6showsintermediateoutputs. elementshouldbeatleasttwiceaslargeasanymetalpiecein
the image to use replacement values outside the artifacts. In
order to prevent a too-large structuring element from flood-
IV.A. Segmentationofmetal
ing a local minimum with diffuse bright artifact values, we
Metalvoxelsaresegmentedbyregiongrowing.Seedsfor recommend clamping the image at a high soft tissue value
regiongrowingarevoxelsabove7000Hounsfieldunits(HU). (100HU)andrestoringtheimagevaluesafterCBR.
Neighboring voxels are successively added to the region if The OBR and CBR operations will also remove anatom-
they are above 3000 HU. Teeth, which have the highest CT ical structures that are smaller than the scale of the structur-
intensitiesforhumantissue,areusuallylessthan3000HU,so ing element. Figure 6(b) shows the result of OBR and CBR,
anatomyisnotincludedinmetallabels.Labelsaregenerated whereanatomicalstructureshavebeeneliminatedalongwith
foreachconnectedmetalregion. artifacts.Werestoretheanatomicalstructurestothepriorby
Wethencontourtheanatomytopreventdarkartifactsfrom using the following steps to discriminate between anatomy
blendingintothesurroundingair.Figure6(a)showsthecon- andartifacts.
tour by a dotted line. If dark artifacts blended into the sur-
rounding air (as shown in the rectangle), they would not be
IV.C. Recoveryofnonadjacentanatomicalstructures
interpretedaslocalminima.Thecontouringmethodweused
isdescribedinAppendixA. TheOBRandCBRprocessedimageissubtractedfromthe
originalimage.Inthisdifferenceimage[showninFig.6(c)],
smallintensitydifferences,attributedtonoiseorartifact,are
IV.B. Removaloflocalmaximaandminima
eliminatedbythresholding.Weusedathresholdofthreetimes
Metal artifacts create local maxima and minima around theimagenoise,whichwasabout20HUforourheadimages.
the metal. The removal of maxima and minima are, respec- Next, the positive and negative differences are consid-
tively,performedusingclosing-by-reconstruction(CBR)fol- eredseparately.Regiongrowingisperformedonthenegative
MedicalPhysics,Vol.39,No.10,October2012
5862 Karimietal.:SegmentationinCTMAR 5862
(a). original with contour (b). after OBR/CBR (c). a - b
(d). positive labels (e). Non-adj. restored (f). Prior
FIG.6. StagesinthegenerationofapriorimageofaCTheadscan.Thecontourisrepresentedbythedottedlinein(a),theimageafterCBRandOBRin
(b),thedifferenceimage(a–b)in(c),labelsfromthepositivedifferencesarein(d),therecoveryofnonadjacentlabelsisshownin(e),andthefinalpriorimage
isin(f).WW/WL=200/0HUin(a),(b),(e),and(f).WW/WL200/−1000HUin(c).
voxels of the difference image, using an intervoxel intensity are not neighboring the metals. The artifact labels may bor-
thresholdthatdependsondistancefromthemetal.Thejusti- dermetalpiecesorbeseparatedbyinterferencepatterns.We
ficationforadistance-dependentthresholdisoursecondob- defineaneighborhoodsizeof10mm.
servation, from which we infer that intervoxel variations in
artifacts decrease as distance from the metal increases. The IV.D. Recoveryofadjacentanatomicalstructures
distance-dependent threshold limits the grouping of artifacts
If a region of positive artifact grows into bone, voxels
withanatomy.Wegenerateadistancetransformoftheimage,
containing bone would be included in the labeled region
which is the smallest distance from each voxel to any metal
and incorrectly replaced with soft tissue values. We exploit
voxel.Ourdistance-dependentthresholdisdefinedasafunc-
the observation that artifact amplitude drops as a function
tionofthedistancetransformas
of distance from metal. In the positive labeled regions
T (x)=max(Te−aD(x),T ), that remain after the recovery of nonadjacent labels, voxel
δ min
values of the original image are plotted against the distance
whereD(x)isthevalueofthedistancetransformatlocationx, transform values. Figure 7(a) shows an example plot. The
andT,Tmin,andaareconstants.WechooseTtobe5000HU artifactsgenerateaclusterintheplot.Toseparatetheartifact
to accommodate the large variations in or near the metal, a cluster from nonartifact voxels, a set of exponential curves
=0.05,andTmin =100.Allvaluesweredeterminedempiri- is generated with different parameters. The equation for the
cally,butarenotcriticalasshownbyexperiments(varyinga familyofcurvesis
from0to0.2)andTminbetween50and200HU. I (D)=(I −I )e−cD +I ,
Similarly,regiongrowingisperformedonthepositivedif- c max min min
ferences.Weusethesameequationandparametersforregion wherecisthecurveparameter,I isthemaximumvaluein
max
growingofthepositivevoxels.Fromregiongrowing,weget the region, and D is the distance (distance-transform value).
labels of positive or negative polarity. For example, positive I istheminimumoftheregionvaluesandanoutlierbound
min
polaritylabelsareshowninFig.6(d). of200HU.Theoutlierboundvaluewaschosenbecausethe
If the labeled regions are not in the neighborhood of a minimum CT number within the artifact cluster was always
metallabel,theyareinterpretedasanatomicalstructures,and about 200 HU in our test cases. In any case, lower artifact
thevoxelvaluesoftheoriginalimagearerestoredinthosere- intensities are removed in a subsequent step. We used no
gions.Figure6(e)showstherecoveryoflabeledregionsthat outlierboundsfortheupperCTvalue.
MedicalPhysics,Vol.39,No.10,October2012
5863 Karimietal.:SegmentationinCTMAR 5863
s
w
e
Vi
Detector Samples
FIG. 8. Reprojectionsthroughmetalvoxels.Thelocalmaximaareshown
byredcircles.
negativelabels.Foreachnegativelabel,wereprojectitscen-
troidinthedirectionofitslargesteigenvector.Ifthecentroid
projects onto a point that is in the neighborhood of the local
maximum of the metal trace, then the label is considered an
artifact.Wehaveusedaneighborhoodsizeof10◦ and5mm.
Todeterminethelocalmaximainprojectionspace,were-
project only the metal. Figure 8 shows the metal traces. In
each projection view, we find the maxima in the sample di-
rection.Thereareoneormoremaximaineachview.Foreach
traceandforeachview,weextractthelocalmaximavalues.
We fit a sliding polynomial to the local maxima values ex-
tracted at each view. Then we locate the maxima of the fit-
ted curve. The sliding polynomial reduces spurious maxima
FIG.7. Therelationshipbetweentheimageintensityandthedistancetrans- whichmayappearduetonoiseorsamplingerrors.
formfortheimageinFig.6isshownin(a).Thelineshowstheexponen-
We do not use the correspondence of each negative label
tialcurvewiththebestseparationofartifactfromanatomy.Thenumberof
with a local maximum in projections as a requirement for
pointsundereachexponentialcurve,normalizedbytheareaunderthecurve
isshownin(b). classifyingthatlabelasanartifact,becausewithmultiplepro-
jectionmaxima,negativeregionsmaybecombinedbyregion
For each curve, the number of voxels under that curve is growing.However, thestepallowsustoaddthenegative la-
normalizedbytheareaunderthecurve.Thenormalizednum- bels that have broken off from the metal. We cannot apply
ber of voxels drops past the curve that includes the cluster thisruleforthebrightartifacts(imagelocalmaxima),because
[Fig.7(b)].Wechooseacurvethatispastthepeak,toensure they do not correspond to local maxima in projection space,
we have captured the cluster, even at the expense of some rather,theycorrespondtoregionsbetweenthelocalmaxima,
anatomy.Voxelsabovethiscurvearerecoveredbecausethey andaremorediffuseinappearancethanthedarkartifacts.
aremoreinfluencedbyanatomythanartifact. Theprocessedimagecreatedthusfaristhresholdedbythe
In order to determine whether artifact has indeed grown method in Appendix B. The segmented metal voxels are re-
into anatomy, we compute from the intensity-distance plot, stored,becauseeachmetalpieceisarealstructure,notanarti-
the highest CT value at each distance. We use distance bins fact.Thiscompletesthegenerationofthepriorimage,shown
of1mm.Thenwecomputetheskewnessoftheresultingdis- inFig.6(f).Thispriorisreprojected,andthereprojectionsare
tribution. If the skewness is less than 0.15, we assumed that usedininterpolation,asdescribedinAppendixC.
theselabelsmusthavegrownintobone.Theskewnessthresh-
oldwasselectedbecauseinourimages,theskewnessvalues
wereabout0.3inimageswheretheartifactdidnotgrowinto V. RESULTS
bone,andfrom−0.1to0.1whenitdid.
OurmethodwastestedonaxialheadCTscans.Figures9
and 10 each show four sets of images with metal artifacts
IV.E. Deletionofdarkartifacts
produced from metal coils in cerebral aneurysms, from a
The local maxima of the projection data are located and deep brain stimulator, and from dental fillings. The original
matchedwithnegativelabels.FromSec.III,wehaveseenthat images and images corrected by our MAR algorithm are
thelocalmaximainprojectionspacecorrespondtodarkarti- shown. For comparison, results are also shown for DI-
factsintheimage.Wecomputecentroidsandeigenvectorsof MAR (i. e., no prior) and an exemplified TP-MAR (i. e., a
MedicalPhysics,Vol.39,No.10,October2012
5864 Karimietal.:SegmentationinCTMAR 5864
FIG.9. Foursetsofimagesareshownincolumns.Eachsetcontainstheoriginal,theproposedmethod,DI-MAR,andTP-MARcorrectedimages.Inallcases,
theproposedpriorresultsinimprovementsovertheoriginalimage,overDI,andoverTP-MAR.WW/WL=200/40HUforcases1,2,and4,and=500/0HU
forcase3.Thearrowspointtoexampleswediscussinthetext.
MedicalPhysics,Vol.39,No.10,October2012
5865 Karimietal.:SegmentationinCTMAR 5865
FIG.10. Anotherfourexamplesareshown.WW/WL=200/40HUforcases5–7and500/40HUforcase8.
thresholded prior, using the same interpolation described in There are residual artifacts in the dental images, especially
AppendixC).TogeneratetheTP-MARprior,wehavethresh- in Fig. 10 (arrow 3), resulting from the imperfect separation
olded the DI-MAR corrected image. Although TP-MAR of the artifact from teeth. DI-MAR produces secondary arti-
algorithms vary, we used CT number values recommended factscomparabletotheoriginalmetalartifacts.TheTP-MAR
in Refs. 5 and 6 for threshold values. Then we used the imagesarebetterthanDI-MARbutnotasgoodasthosewith
difference interpolation method for data replacement in the ourprior.
metal traces. In the case of DI-MAR, the metal piece was
added back to the corrected image. Since the interpolation
VI. DISCUSSION
techniquewasthesameforalltheMARalgorithms,theprior
imagedeterminedtheimprovement. VI.A. Analysisofresults
Artifacts are removed by our algorithm even for multiple Ouralgorithmpreservesanatomicalstructuresintheprior
metalpieces,andfromlargemetalpieceswhichproducedark image, which is why secondary artifacts are reduced. DI-
artifactsbelow−1000HU(arrow1),andbrightartifactswith MARlosesedgeinformationinthemetaltraces,soestimated
the CT intensity of bone or cartilage (arrow 2). There were data are inconsistent with the rest of the sinogram and sec-
fewer secondary artifacts with our method than the others. ondaryartifactsaregenerated,asshownbyarrow4inFig.9.
MedicalPhysics,Vol.39,No.10,October2012
5866 Karimietal.:SegmentationinCTMAR 5866
TP-MAR also loses edge information as described later, but
toasmallerextentthanDI-MAR.
Our method operates at a region level while TP-MAR
methodsoperateatavoxellevel.Aregionismoreinformative
thanavoxelindistinguishingartifactsfromanatomybecause
theregioncanbeexaminedagainstmorecriteriathanthresh-
olds. We identify a region with constant polarity, and test if
theregionorapartofitfitsthecriteriaforartifacts.Consid-
eringthepositiveregions,weusetheintensity/distancerela-
tionship(observation2)toseparatebonefrombrightartifacts.
Negative regions can be matched with local maxima in the
Radontransform(observation3)toidentifythemasartifacts.
A TP-MAR prior is generated by thresholding an image
into air, soft tissue, or bone. With misclassification of arti-
fact voxels as bone or air, large errors may be produced in
the final corrected image, especially when the ratio interpo-
lation (discussed below) is used. Therefore, some TP-MAR
methods recommend a first-pass DI-MAR to provide an im-
age with smaller amplitude artifacts, so that thresholding of
this corrected image is more likely to create a good prior.
However, our experiments show that the secondary artifacts
withDI-MARmaybecomparabletoorworsethantheorigi-
nalmetalartifacts.UsingpoorqualityDI-MARimagesyields
poor priors. In the deep-brain stimulator case with multiple
metals (Fig. 9, case 2), the DI-MAR secondary artifacts are FIG.12. Imagesreconstructedwithratiointerpolation.Thepriorsforthetop
worsethanthemetalartifactsintheoriginalimage(arrow4). rowimageswerethresholdedDI-MARimages(Fig.9,row3),andthepri-
Thisleadstoapriorthatisworsethanthepriorfromthresh- orsforthebottomrowimageswerefromourproposedmethod.Alongwith
Fig. 9, these images suggest that the prior is more critical than the inter-
oldingtheoriginalimage.Inthedentalcases(cases4and8),
polationtechnique.However,formultiplepiecesofmetal,theinterpolation
the secondary artifacts are misclassified as anatomical struc-
techniqueitselfbeginstoplayaroleasseenbythedegradedimagequality
turesbythresholding,andpreservedorenhancedinthefinal oftherightcolumncomparedtoFig.9.
TP-MARcorrectedimage.
Figure11showsthedifferentpriorsresponsiblefortheim-
age quality in case 4. The original image is corrected with lossofsomeedgestructure,especiallyifthosestructuresare
DI-MAR.TheDI-MARimageisthenthresholdedtoproduce closetothemetalpieces.
aTP-MARprior.Notethepartiallossofboneandairpocket Inordertostudytheimpactoftheinterpolationtechnique
structuresinthethresholdedprior(arrow5),andbetterpreser- relative to the accuracy of the prior, we also use the ratio
vationintheproposedprior.Thepartiallossiswhytheimage interpolation with our prior and with the thresholded prior.6
qualityoftheTP-MARimageisinbetweenthatofDI-MAR This interpolation uses the ratio of scanner projections and
and our method. The DI-MAR algorithm itself results in the reprojections, so that the reprojections themselves need not
be substituted in the trace. This is a potential improvement
on methods that directly use the reprojections in the metal
trace,e.g.,5becausereprojectionsmaythemselvesbeinaccu-
rate.TheinterpolationtechniqueresultsareshowninFig.12.
ComparingwiththedifferenceinterpolationshowninFig.9,
we see that in the case of a single metal object, the image
quality of the difference interpolation and ratio interpolation
methods is nearly the same when our prior was used. The
image quality of the ratio interpolation image is nearly the
sameasthedifferenceinterpolationalsowhenthethresholded
prior was used. These two results indicate that the prior has
greaterimpactthantheparticularinterpolationmethod.Note
that the ratio interpolation however only works well when
themetalobjectisremovedfromtheprior.Thereasonisthe
FIG.11. Priorscorrespondingtocase4.TheleftimageshowstheTP-MAR sameonethatwediscussinSec.III,i.e.,thatbeamhardening
prior,createdbythresholdingtheDI-MARimage(Fig.9,row3).Thelossof coupled with thresholding creates overestimated or underes-
anatomicaledgesnearthemetalleadstoafinalTP-MARimagethatisnot
timated reprojections. These reprojections are naturally con-
muchbetterthanDI-MARinthisexample.ThisisbecausebothDI-MAR
sistentacrossviewangles,butnotoncetheyaremultipliedby
andTP-MARhavemissednearlythesamestructures.Therightimageisthe
proposedprior,andpreservesmoreoftheanatomicalstructure. the ratio of projections, which are not a constant unless the
MedicalPhysics,Vol.39,No.10,October2012