Geometrical methods for level set based abdominal aortic aneurysm

16 Pages·2012·2.49 MB·English

Author's personal copy
computer methods and programs in biomedicine 107 (2012) 202–217

journal homepage: www.intl.elsevierhealth.com/journals/cmpb

Geometrical methods for level set based abdominal aortic

aneurysm thrombus and outer wall 2D image segmentation

Christos Zohiosa, Georgios Kossiorisa,b, Yannis Papaharilaoub,∗

aDepartmentofMathematics,UniversityofCrete,KnossouAvenue,P.O.Box2208,Heraklion71409,Greece

bInstituteofAppliedandComputationalMathematics,FoundationforResearchandTechnology-Hellas,VassilikaVouton,P.O.Box

1385,Heraklion71110,Greece

a r t i c l e i n f o a b s t r a c t

Articlehistory: Abdominalaorticaneurysm(AAA)isalocalizeddilatationoftheaorticwall.Accuratemea-

Received18December2009 surementsofitsgeometriccharacteristicsarecriticalforareliableestimateofAAArupture

Receivedinrevisedform risk.However,currentimagingmodalitiesdonotprovidesufficientcontrasttodistinguish

9December2010 thrombusfromsurroundingtissuethusmakingthetaskofsegmentationquitechalleng-

Accepted28June2011 ing.Themainobjectiveofthispaperistoaddressthisproblemandaccuratelyextractthe

thrombusandouterwallboundariesfromcrosssectionsofa3DAAAimagedataset(CTA).

Keywords: Thisisachievedbynewgeometricalmethodsappliedtotheboundarycurvesobtainedbya

Abdominalaorticaneurysm LevelSetMethod(LSM).Suchmethodsaddresstheproblemofleakageofamovingfrontinto

LevelSetMethods sectorsofsimilarintensityandthatofthepresenceofcalcifications.Theversatilityofthe

Thrombusandouterwall methodsistestedbycreatingartificialimageswhichsimulatetherealcases.Segmentation

segmentation qualityisquantifiedbycomparingtheresultswithamanualsegmentationoftheslicesof

Stoppingcriterion tenpatientdatasets.Sensitivitytotheparametersettingsandreproducibilityareanalyzed.

Wallthickness Thisisthefirstworktoourknowledgethatutilizesthelevelsetframeworktoextractboth

thethrombusandexternalAAAwallboundaries.

©2011ElsevierIrelandLtd.Allrightsreserved.

stress, wall stiffness, intraluminal thrombus thickness and
1. Introduction
walltensionhavebeenproposed(seee.g.[1]).Tostudythese

indicatorsandobtainamorereliablepatientspecificestimate

Theaortaisthelargestarteryinthehumanbodyandthemain of AAA rupture risk, accurate geometric characterization of

bloodvesselleadingawayfromtheheart.Anabdominalaor- theaneurysmiscritical.

ticaneurysm(AAA)isapermanentandirreversiblelocalized 3D reconstruction of the complex anatomy of an AAA

dilatationoftheabdominalsectionofthisvessel.Inclinical from medical imaging data can be achieved through either

practice,diagnosticinformationofthe3DanatomyofanAAA aslice-by-slice2Dsegmentationofthestructuresofinterest

isextractednon-invasivelyinvivo,throughComputedTomog- and the subsequent application of a 3D surface reconstruc-

raphyAngiography(CTA).AAAcangrowprogressivelylarger tion method [2,3] on the extracted boundary points, or by

and may eventually rupture if not diagnosed and treated. utilization of a 3D segmentation approach that extracts the

Accurate estimation of AAA rupture risk remains an open surfacesofinterestinonestep[4–7].Amanual2Dsegmenta-

problem. Numerous risk indicators beyond the peak trans- tionisarathertime-consumingprocedurewiththeadditional

verse diameter, which is not always reliable, such as wall drawback of large intra- and inter-observer variability. An

∗ Correspondingauthor.Tel.:[phone].

E-mailaddress:[email protected](Y.Papaharilaou).

0169-2607/$–seefrontmatter©2011ElsevierIrelandLtd.Allrightsreserved.

doi:10.1016/j.cmpb.2011.06.009
Author's personal copy
computer methods and programs in biomedicine 107 (2012) 202–217 203

automaticsegmentationmethod,whichisbothaccurateand fromneighboringregionsisproposed.Twogeometricalmeth-

robust, would alleviate these problems. Level Set Methods ods are presented to reconstruct the thrombus and wall

(LSM’s) [8,9] have been used in this direction due to their boundaries exploiting the strong image features associated

main advantages compared to other methods, such as the withthepresenceofcalcifications.Theproposedsegmenta-

formulationofthelevelsetwithoutaparametrizationresult- tionalgorithmisappliedtotenpatientdatasets(450slices)

ing in a free transform and change of its topology, and the andtheextractedresultsarecomparedtoamanualsegmen-

relative ease in extending the method to 3D problems. AAA tationobtainedfromamedicalexpert.

segmentationisachallengingtaskforLSM’sdueto(a)poor

contrastbetweenthrombusandsurroundingtissuethatcur-
2. Levelsetframework
rent imaging modalities produce and (b) strong edges that

manyneighboringstructurespresent.OtherissuesthatLSM’s
2.1. Themainequations
faceare(a)theleakageoftheadvancingfrontintoregionsof

similar intensity and (b) their high computational complex-
The application of the level set framework consists of two
itywhichmakesthemslow.TheseissuesexplainwhyLSM’s
phases. In the first phase we use the FMM. In this method,
havebeenprimarilyusedsofaronlyforlumenandouterwall
thepositionoftheexpandingfrontischaracterizedbycom-
AAAsegmentation.ExamplesofLSMapplicationforAAAseg-
putingthearrivaltimeT(x,y)ofthefrontasitcrosseseach
mentation are the works by Loncaric et al. [10] and Subasic
point(x,y)[8].If(cid:2) istheinitiallocationoftheinterfacethen
etal.[4]whoimplementeda3DLSMforAAAandintroduceda
T=0on(cid:2).Attimetthemovingfrontisgivenby
curvestoppingcriterionandtheworkbyMageeetal.[5]which

utilizes a 3D deformable model for an initial segmentation
combin ed wi thanefficie ntlevel set im plemen tation.Zhuge (cid:2)(t)={(x,y)∈R2|T(x,y)=t}. (1)

et al. [11] performed segmentationof the outer AAA wall in

five steps: preprocessing, global region analysis, surface ini- Thefrontismovingwithaspeedfunctioninthenormaldirec-

tialization,localfeatureanalysisandlevelsetsegmentation tiongI>0,sotheequationforthearrivalfunctionT(x,y)is

andSonkaetal.[7]appliedanovel4Doptimalborderdetection

algorithm for automatic surface segmentation of AAA lumen. gI|∇T| = 1. (2)
IntheworkofdeBruijneetal.[12]aninteractivemethodfor

AAAsacsegmentationispresentedwhichreliesonthefitting
Here, the proposed edge-indicator function, which restrains
ofashapemodeltopointswithhighcorrelationwiththeref-
theevolvingfrontfromleakingoutofthedesiredregionis
erencecontour.Thesamegroup[13]alsoproposedamethod

dsbeaeqslieundee naoctnieos n.A Oocltfai vbthea erSr oihauagtpaeere taManloe.d[u1er5yls]s imn(At rbSooMduusn)c de[ad1r4ay] nifneow rm mauuletttiohpmoleda tMfeoRdr gI(x, y) = 1 + (cid:3)(cid:2)(cid:2)∇Gs1∗ I(x, y)(cid:2)(cid:2)2. (3)

deformable model-based segmentation of thrombus in AAA
based on a 3D discrete d eformable mod el (DM), Bod ur et al. I is the image intensity and for ˛ ∈R2, Gs(˛) is the Gaussian
withwidths,
[16] presented an automatic segmentation of the aortic bor-

derbyobtainingorthogonaltothecenterlineslicesandusing
a no ve l variation of the isop eri met ric algorith m wh ich incor- Gs(˛) = 41(cid:4)se−(|˛|2/4s), (4)
poratescircularconstraintsandDePutteretal.[6]createdan

initial3Dactiveobjectfromthelumencenterlinewhichwas
which is used to reduce possible noise effects in the image.
theniterativelydeformedviaatime-discretizedsecondorder
(cid:3) is a positive parameter. gI is close to unity away from the
Newtonian-evolution equation to track the boundary. None,
boundariesanddropstozeronearsharpchangesintheimage
however, of these AAA segmentation methods is applied to
gradient.Thesechangespresumablycorrespondtotheedges
trackthethrombus-wallboundary.Toourknowledge,theonly
of the desired shape. Efficient numerical schemes (upwind
work that presents a semi automatic method to extract the
differencing)canbeusedtosolveEq.(2)(see[8]).
AAAwallthicknessfromCTAimagedataisthatofShumetal.
Inthesecondphaseweembedtheinitialpositionofthe
[17]whichdoesnot,however,relyontheLSMbutratheruti-
frontasthezerolevelsetofahigher-dimensionalfunction(cid:5)
lizesintensityhistogramsforthethrombussegmentationand
andwelinktheevolutionofthisfunctiontothepropagation
aneuralnetworkfortheouterwallsegmentation.
ofthefrontitselfthroughatime-dependentinitialvalueprob-
The objective of the present work is to use the level set
lem.Atanytime,thefrontisgivenbythezerolevelsetofthe
framework and introduce geometrical methods to generate time -de pen dentl evel setfu nc tion(cid:5) .S o,at time t≥0 ,the fr ont
a semi-automatic segmentation algorithm that accurately
is
extractsthethrombusandouterwallboundariesfromcross

sectionsof3DAAAimagedata.Thisisachievedbyexploiting
(cid:2)(t)={(x,y)∈R2|(cid:5)(x,y,t)=0}, (5)
anyavailablesourceofcontrastatthethrombus-wallbound-

ary including the presence of calcifications which appear in
where(cid:5)(x,y,0)istheresultoftheFMM.Thelevelsetequation
AAACTAimagesasveryhighintensityregionsandarecom-
usedinthisphaseis
monly part of the wall. In cases where there is insufficient

informationataregionoftheboundarytobeexploitedbythe
LSM, a meth od whereb y wal l thickness c an be interp ola ted (cid:5)t=gI(x,y)((cid:6)−c)|∇(cid:5)|+∇(cid:5)·∇P. (6)
Author's personal copy
204 computer methods and programs in biomedicine 107 (2012) 202–217

gIisthesameasbefore,(cid:6)istwicethemeancurvature,c>0a initializationwetaketheresultoftheFMM.Asanouterini-

parameterusedforfasterconvergence[8]and tializationweconsiderthesmallestcirclecontainingtheinner

(cid:2) (cid:2) initialization,extendedbyapropernumberofpixelsp(forfur-
(cid:2) (cid:2)
P(x,y)=− ∇Gs∗I(x,y) . (7) ther details on p see Section 5). The final result is the inner

curveatthetimeofcontactofthetwocurves(seeFig.1(e,f)).

Theabovelevelsetequationconsistsofthreeterms:gI(cid:6)|∇(cid:5)| Acomparisonofthethreestoppingcriteria,basedonthe

which acts as a smoothing term for the front away from speedoftheirresponsetoaleakage,ispresentedinSection

regions where the gradient magnitude is large, gI(−c)|∇(cid:5)| 4.7.Theresultingcurveisthenusedforthereconstructionof

whichisadrivingexpansionforceand∇(cid:5)·∇Pwhichhasasta- theouterwallboundary(Section3.2).

bilizingeffect,asPconstantlyattractsthefronttotheedges

in the image [8,9,18]. See Sapiro [18] and Antiga [19] for the
3. Geometricalmethodsforthe
derivationof(6),and[9]foritsnumericalapproximation.Two
re constructionofth rombusa nd outerwall
othervariationalformulationsoftheLSMcanbefoundinLi
boundaries
etal.[20]andGelasetal.[21].

2.2. DetectionofLeakageRegionmethodandstopping 3.1. Thrombusboundaryreconstructioninthe
crite ria pres enceofcalcifica tions

The main problem of the LSM in segmenting the outer aor- The presence of calcified deposits which are located within
tic b ound ary is the l eaka ge o f th e moving f ront (whic h is theA AAwallc om binedwi thpoorc ontrast betw eenthr ombus
the thrombus in itiali zation)in to surr ounding tissue inregio ns and oute rwal lresultina “too thed ”thromb usbound ary,asthe
with image in tensity simila r to that of the aortic w all (see evol ving c urve enter s t h e regions between t he calcific ati ons
they ellowc urveinFig .1(d)).B ef oread dr essi ngthis issue ,we (see Fig. 2). Ass uming tha t region s in-betw een calcifications
des cribe a metho d whi ch ca n be u sed to dete ct th e regio ns are c ons ide red to be w all t hey sho uld be exclu ded from the
where th e leakage takes p lace (se e Fig . 1( a)) wh ich i s called bou ndary.Toac hie ve this, weco nsider the inversion mapp ing
Detect ion ofLeaka geReg ion(D LR)m eth od.W efirst sp litthe withrespe ct tothec urve (cid:2) o btained byt heLSM.T heinver-
interiorof th ethromb usinitia lizatio ncurvei nto secto rs,us ing sion mapping is the usualt ra nsformat ion T(z) =1/(z −z0) inthe
radials sta rtin g from its center of m ass, w hic h are eq ually com plexplan e, whe rez0= (x0,y0)isthec ente r ofm a ss of the
spaced byause r-defi ned sector an gle(cid:5).T heinit ializ ationis domainb ounde dby(cid:2) .We co nsid er the doma in bound ed by
evolvin gw it hthelevelse tequa tiona nd are initialization as T((cid:2)) in the image pl an e un der T, an d w e find t he points of
asigned dista nce funct ion takesplac eev er yteniterations .A itsc o nv exh ull.To obtain theth rom bus bou nda ryw etaket he
ca lculati onofthe areaofe achs ectori salso ma deeveryte n im ageund erT− 1 o ftheb oun daryofth econstruc ted con vex
iterations.W e then calc ula teth erelati ve chan geofi tsarea by hull an d we apply t he L SM for a s m all n umber of ite rations
comparing itw ithth esectora rea attheen dofthe pr evi oust en (app roxim at ely 20) to o btain a ph ysiolo gically m or e realistic
iterations.T h eevo luti onisin terru p ted whe na loc alminim um result.
of the rela tive area chan g e has appea red fo r all th e sectors. By using the inversion mapping the existence of leakage
W hen thecurv esto ps,thes ecto rsforwhi cht he relat ivearea does not aff ect the recon struction of t hrombus a nd , there-
chang e is of ord er r (s ee S ection 5 fo r a disc uss ion on r ) for fore, thefi nalre sult isnotsensitive to theinitializ ation ofthe
allitera tio ns define t hele akagere gi oni n theimage.A lth ou gh meth od.

therateofchangeoftheareaenclosedbytheevolvingbound-

aryhasbeenpreviouslyused[22]forleakagedetection,herein 3.2. Reconstructionoftheouterwallboundaryinthe

weadoptamorepreciselocalapproach,wherebysmalland presenceofcalcifications

localizedleakageregionscanbeidentified.

The DLR method provides a stopping criterion which we ThecircumferentialdistributionofcalcificationsintheAAA

willcallSectorCriterion(SC)(seeFig.1(b)).Thesecondstop- wall can be exploited to segment both the outer wall and

pingcriterion,calledLengthDifferenceRateCriterion(LDRC), thrombusboundaries.Thefirststepintheprocedureforthe

stopsthecurvewhenthefirstminimumoftheabsolutevalue reconstruction of the outer wall boundary is a thresholding

of the length difference rate of the curve appears. Here, as of the image to extract the calcification regions (Fig. 3(a)).

lengthdifferencerateweconsiderthedifferenceofthetotal Thetotalareaofthepixelsabovethethresholdiscomputed

lengthofthecurveeveryteniterations.Asuddenincreaseof anddisplayedandtheuserselectsanupperandlowerbound

thisrateisduetotheleakage,sothecurvehastostop(see forthethreshold.Thenaniterativeprocessseekswithinthis

Fig.1(c,d)). rangethethresholdvaluethatminimizesthederivativeofcal-

For the third stopping criterion, called Internal-External cificationareawithrespecttothresholdvalue.Iftheresultis

ContactCriterion(IECC),weconsideraninnerandanouter(to satisfactory (i.e. the SNR (Signal to Noise Ratio) of the final

theaorticwall)initializationwhichmovetowardeachother thresholdedimageishigh,wherecalcificationsaresignaland

simultaneously;thisisachievedbyapplyingthelevelsetEq. alltheotherpartsconstitutenoise)itisacceptedbytheuser.

(6),withparametercpositiveinthecaseoftheinnerinitializa- Otherwiseanewthresholdrangeisspecifiedandtheprocess

tionandnegativeinthecaseoftheouterinitialization(with isrepeateduntilasatisfactoryresultisobtained.

thesameabsolutevalue).Theevolutionisinterruptedassoon Followingtheapplicationofastoppingcriteriontotheevo-

as the two evolving curves come into contact. As an inner lutionofthethrombusinitializationwiththeLSM,weconsider
Author's personal copy
computer methods and programs in biomedicine 107 (2012) 202–217 205

Fig.1– ResultsoftheDLRmethodandthethreestoppingcriteriaforthesamecrosssection:(a)theleakageregions
dete c te dbythe DL Rm etho daresh own inb lack (r=10−2). Theang leo fth esect orsis (cid:2)/24.(b) Th ere sultofSC whichstops

thecurveafter80iterations.(c)TheinitialcurveusedforLDRCobtainedbytheFMM.(d)Theblackcurveistheresultof

LDRC(40iterations)andtheyellowcurvetheresultafter100iterations,showingtheleakageifLDRChadnotbeenused.(e)

TheouterandtheinnerinitializationusedforIECC.(f)Thefinalouterandinnercurves,wherethelatteristheIECCresult

(60iterations).(Forinterpretationofthereferencestocolourinthisfigurelegend,thereaderisreferredtothewebversionof

thearticle.)

the extracted curve and the smallest circle containing this contrastbetweenthrombusandouterwallexists,therecon-

curve.Atleastonepointofeachcalcificationwillbeinterior structedthrombusandouterwallboundariesobtainedusing

to this circle (see Fig. 3(b)). We then apply the region grow- the methods described in Sections 3.1 and 3.2 may coin-

ingmethodwitheachcalcificationpixelinteriortothecircle cide in certain regions. However, although in our analysis

acting as a seed and we use the threshold extracted in the suchcaseswererare,toaddressthisissuewedevelopedthe

previous step as a rule of similarity to locate all the calcifi- TWBDmethodwhich,basedontheassumptionthatthewall

cationpixelsintheimage.Wethencomputetheconvexhull thickness distribution should present a smooth continuity,

oftheunionofallmid-edgepointsofthecalcificationpixels smoothlyinterpolatesinformationfromneighboringregions

andthepointsofthecurvewhichwastheresultoftheLSMon to compute wall thickness in regions where thrombus and

thethrombusinitializationcurve.Theboundaryoftheconvex outerwallboundariesaretooclosetoeachother.

hull(thedashedblackcurveinFig.3(e))willcertainlycontain The first step in the TWBD method is to consider many

allthecalcifications.Theprocedureisconcludedbyapplying raysstartingfromthecentroidofthethrombus,sothattwo

thelevelsetequationtothisboundary,buttakinganinward successiveraysformthesameangle(cid:5),andfindthedistance

motionthistime,aimingatasmootherresultwhichismore (wallthickness)betweenthetwointersectionpointsofeach

biologicallyplausible(yellowcurveinFig.3(e)).Fig.3(f)depicts ray with the thrombus and outer wall curves. We find the

both the thrombus and outer wall boundaries extracted by thrombuspointswhichareconnectedwithdistancesoforder

themethodspresentedinthisandtheprevioussectionwhich d/10mm (see Section 5 for details on d). For each point, we

exploitthepresenceofcalcifieddeposits. find the nearest n thrombus points in the counterclockwise

directionwhosedistances(wallthicknesses)areoforderdmm

3.3. ThrombusandWallBoundariesDistance(TWBD) and calculate the mean of these distances. Then, we find

method whichofthesenpointshasdistanceclosertothismean.We

denotethisdistanceasd andthedistancefoundfollowingthe
l
In images where calcified deposits are not present or their sameprocedureintheclockwisedirectionasdr.Thenalinear

AAA wall circumferential distribution density is low and poor interpolation between dland dris made to obtain a new wall
Author's personal copy
206 computer methods and programs in biomedicine 107 (2012) 202–217

Fig.2– ThestepsforthereconstructionofthethrombusboundaryofthecrosssectionofFig.1:(a)theyellowcurveisthe

thrombusinitializationobtainedbytheFMMandtheblackcurvetheresultoftheLSM(200iterations).(b)Theimageofthe

blackcurveundertheinversionmappingwithrespecttoitscentroid(inorange)andtheboundaryofitsconvexhull(in

black)aredepicted.(c)Thepreimagesofthepointsoftheconvexhullundertheinversionmappingformadomain,whose

boundaryistheyellowcurve.(d)Theyellowcurveisthefinalresult,afterapplyingthelevelsetequationfor20iterations.

(Forinterpretationofthereferencestocolourinthisfigurelegend,thereaderisreferredtothewebversionofthearticle.)

thickness (and a new thrombus point) for each point in the givesapairofpoints,astartandanendpoint(apoint,here

regionwithdistancesoforderd/10mm.Finally,weapplythe andintherestofthepaper,isthecenterofapixel).Thefront

levelsetequationtothesoreconstructedthrombusboundary is initialized with the start point inside the desired region

for a small number of iterations (e.g. 10) taking an outward and expands until the solution has reached the vicinity of

motion,toobtainasmootherandphysiologicallymorerealis- the desired edge. It stops expanding when it meets the end

ticresult. point for the first time. These points are specified either by

Numericaltestshaveshownthatforn∈{4,5,...,10}the theircoordinatesorinteractivelybycursorpositioningonthe

results are very similar, because of the very small angle (cid:5). image.

Fig.4showstheresultofthemethodfortwocrosssections, Toachieveagoodinitialization,whichisimportantfora

onewhereapproximately30%(Fig.4(a))andtheotherwhere satisfactory final result, the initial boundary should: (a) be

only approximately 8% (Fig. 4(d)) of the boundary perimeter close to the contour we want to detect, (b) not have leaked

includescalcifications. outofit,and(c)nothavebeencapturedbythelumeninsome

region. In order to obtain an initialization that satisfies the

aboveconstraintsthestartpointoftheFMMshouldbeeither

4. Experiments nearaconceivablecenterofmassofthethrombusornextto

aregionofhighcontrastandtheendpointshouldbeatalow

4.1. ImplementationoftheFMMandLSM contrastregion,inordertoavoidleakage.

However,forimagesthatareratherdifficulttosegmentand

WeappliedtheFMMusingtheMatlab-basedFastMarching obtainaproperthrombusinitialization(duetoleakageprob-

Toolbox[23]createdbyGabrielPeyrétoobtainathrombusini- lemsorstrongresponsesfromneighboringregionse.g.spine,

tialization for the LSM. We took (cid:3) = 50 for gIand we used a 5 × 5 lumen or calcified deposits) more start and end points should
Gaussian filter with s=1. In this Toolbox, the user typically beselected.Forthispurpose,wemodifiedtheFMMToolboxto
Author's personal copy
computer methods and programs in biomedicine 107 (2012) 202–217 207

Fig.3– (a)–(d)ThestepsforthereconstructionoftheouterwallboundaryinthecrosssectionofFig.1:(a)appropriate

thresholdingtolocatethecalcifications.(b)TheLSMextractedcurveinblack(usingLDRC,40iterations)andthesmallest

circlecontainingthiscurvearedepicted.(c)Detectionofthecalcificationpixels(showninblack)usingthiscircle.(d)The

dashedcurveistheboundaryoftheconvexhulloftheunionofthemid-edgepointsofthesepixelsandoftheblackcurve

depictedin(b).Theyellowcurveisthefinalresult,afterapplyingthelevelsetequationtothisboundaryfor20iterations.(e)

depictstheresultofthemethods(thrombusandouterwallboundariesinblackandyellow,respectively)forthesamecross

section.(Forinterpretationofthereferencestocolourinthisfigurelegend,thereaderisreferredtothewebversionofthe

article.)
allowtheusertoselectuptofourpairsofpoints.Thisnumber beenobtainedafterexperimentationwithvariousvaluesand

ofpoints,whichwasdeterminedafterextensiveexperimen- itisparameter-independent.Incontrasttothat,thenumber

talwork,wassufficienttoproduceagoodinitializationforthe of iterations of the LSM (applied again to the thrombus ini-

mostdifficulttosegmentcases.Inthiscasetheinitialization tialization)fortheouterwallsegmentationisdeterminedby

istheunionoftheresultsextractedforeachpair.TheFMM astoppingcriterion.

methodiscodedinC++andisthereforeveryfast.

FortheimplementationoftheLSMweusedtheLevelSet 4.2. Anoverviewofthealgorithm

Toolbox[24]createdbyBarisSumengen.InthisToolboxthere

aretwobasicroutines,thefirstisforthereinitializationofthe The basic steps of the segmentation algorithm proposed in

levelsetfunction(cid:5)asasigneddistancefunction(usuallyevery thisworkareshowninFig.5.Initially,asuitablecroppingis

fiveorteniterations)sothat(cid:5)remainssmoothasitevolves performed on the selected slice of the AAA CTA data set to

andthesecondisfortheevolutionofthecurvewhereallthree extracttheregionofinterestinordertoreducesegmentation

typesofmotionareapplied(motioninthenormaldirection, time. The FMM is first applied and the user is interactively

motioninvolvingmeancurvatureandmotioninanexternal asked to enter up to four pairs of start and end points in
forcefie ld).Weto ok(cid:3)=0 .1,c=5000 and a5×5 Ga uss ianfilter order too btain an opt imize dinit iali zation tob eus edfort he
with s=1.M ore over, w e too k r =10− 2 as th e o rd erofthe rela- throm bu sandt heo uterwalls egmentations .Fo rt hethr om bus

tiveareachangeintheDLRmethod,p=1.5pixelsintheIECC boundarysegmentation,theLSMisappliedtothethrombus

andd=1mmastheorderusedintheTWBDmethod.Foran initializationtoobtainaboundarycurvetowhichthethrom-

explanationofthevaluesofalltheparametersusedandfora busboundaryreconstructionmethodisthenapplied.

sensitivityanalysis,seeSection5. In the outer wall boundary reconstruction, the leakage

To perform the thrombus boundary segmentation, the issuesoftheadvancingfrontareresolvedbyusingtheLSM

LSM is applied to the thrombus initialization for 200 itera- tothethrombusinitializationaccompaniedbyoneormore,if

tionsbeforetheboundaryisreconstructedwiththemethod required,oftheSC,LDRCandIECCstoppingcriteria.Thestop-

describedinSection3.1.Theabovenumberofiterationshas ping criterion which detects the leakage region earlier thus
Author's personal copy
208 computer methods and programs in biomedicine 107 (2012) 202–217

Fig.4– ApplicationoftheTWBDmethodfortwocrosssectionswithfew(a)andno(d)calcifications(d=1mm,n=4and

(cid:3)=(cid:2)/200).Thrombusandouterwallboundariescoincideinsomeregionsin(a),duetothelackofcalcifications.Theresultof

theTWBDmethodisshownin(b)andthefinalresult,afterapplyingtheLSM,in(c).Aproperdistancebetweenthetwo

boundarieshasbeensettled.Acrosssectionwithveryfewcalcificationsisdepictedin(d).(e)presentstheintensityprofile

alongthedashedlineof(d).Thisshowsthatthereisintensitydifferenceinsidethewall(regionsA–BandE–F)comparedto

thrombus(regionsB–CandD–E).Thetwoboundariesmaycoincideinsomeregionsbuttheexistenceofregionswith

distancebetweentheboundariesduetothisdifferenceisexploitedviatheTWBDmethod.Thefinalthrombusboundary

(dashedline),afterapplyingtheTWBDmethodandtheLSM,andtheouterwallboundary(continuousline)areshownin(f).

withtheleastnumberofiterationsisselected.Theshapeof basedonthemeanintensityofrespectiveregionsintheAAA

theleakageregionplaysanimportantroleinthisselection(see CTimagesprocessed.Thereisalsoaregioninwhich“throm-

Section4.7).Theleakageregionofaslicecanalsobedetected, bus” and “outer wall” have similar intensity (right side) to

if desired by the user, applying the DLR method. The outer test leakage detection. The results for thrombus initializa-

wallboundaryreconstructionmethodisthenappliedtothe tionthroughFMMandthetrackedthrombuswallboundaries

stopped boundary curve. As part of the method, the user is showninFig.6(b)and(c)respectivelydemonstratetheeffec-

interactivelyaskedtoapplyanappropriatethresholdingtothe tivenessoftheproposedmethodologyinaddressingtheissues

imageinordertoextracttheexistingcalcifications.Inregions beingtested.

where,duetolackofcalcificationsandpoorintensitygradient, Fig.6(d)showsadifferenttestimage,wherewehavetaken

the resulting thrombus and outer wall boundaries coincide, the(extreme)caseofexactlythesameintensityfor“throm-

theTWBDmethodisappliedtothethrombusboundarytoset bus”and“outerwall”.Therearealsotwoouterleakageregions

anappropriateinterpolatedwallthickness. andeight“calcifications”.Contrastratioandintensityvalues

arethesameasinFig.6(a).Theuniformintensityof“throm-

4.3. Resultsforartificialimages bus”,“outerwall”andleakageregionsis0.46.InFig.6(f)the

maximumerrorinthetwoleakageregionsisapproximately6

To test the versatility of our methods and assess their pixelsor1.6mm,whilethebordersoftheseregionsareplaced

effectiveness in: (a) tracking the thrombus and outer wall 15pixelsoutsidethe“outerwall”boundary.Thisresultshows

boundariesand(b)detectingleakageoutsidetheouterwall, theeffectivenessofLDRCinexploitingthepresenceof“calcifi-

we generated the artificial images shown in Fig. 6(a, d, g). cations”togenerateasatisfactorysegmentationresultdespite

Fig.6(a)includesfiveareasofuniformbutdifferentintensity theabsenceofcontrastatthe“thrombus”,“wall”boundary.

which on 0–1 intensity scale, background is 0 and “lumen”, Noise was not added to these images since it was quite

“thrombus”,“outerwall”and“calcification”haveintensityval- low in the data sets of our experiments and the Gaussian

uesof0.5,0.46,0.48and0.9,respectively.Thesevalueswere filterusedproperlydiminisheditseffects.However,toassess
Author's personal copy
computer methods and programs in biomedicine 107 (2012) 202–217 209

Fig.5– Diagramofthesegmentationalgorithm.Novelaspectsofourapproacharepresentedinboxeswithblack

background.Fewiterationsinthecontextofouralgorithmare,basedonextensiveexperimentalwork,approximately20.

the effectiveness of the proposed method in cases where methods in leakage detection and wall boundary tracking

image noise is significant we added white Gaussian noise in images with an asymmetric and sparse circumferential

ofzeromeanandavarianceof0.0001toFig.6(d).Toassess distribution of calcifications we removed 4 ‘calcification’

the effects of random noise on the ability of our intensity regionsfromFig.6(d)andgeneratedFig.6(h).Thefinalresult

gradient based method to segment the thrombus and wall depicted in Fig. 6(i) demonstrates the effectiveness of our

boundaries in an AAA image we consider the Contrast to methodsinprocessingsuchchallengingcases.

NoiseRatio(CNR)betweenneighboringregionswhichdefine

theseboundaries.HereCNRiscomputedasthedifferenceof 4.4. Manualandsemi-automaticsegmentation
them eanpixelin tensit yoft he correspon din gim ageregio ns

to the standard deviation of the pixel intensities in these A quantitative validation of our algorithm was performed

regions.Withtheadditionofnoisetheresultingimagehada using CTA image data sets acquired from ten patients at

CNRofapproximately2.8betweenthelumenandthrombus the University Hospital - Heraklion. The image reconstruc-

andofapproximately3.5betweenthrombusandbackground. tion matrix size was 512×512 pixels with a pixel width

For reference, note that corresponding CNR values obtained between0.64mmand0.9mm.Thenumberofslicesthatcon-

from typical CTA AAA images in our data set were 9 and 8 tainedaneurysmaltissueineachdatasetvariedfrom29to

respectively.Thefinalresultofthemethodasappliedtothe 56.

noisyversionofFig.6(d)isdepictedinFig.6(g)anddemon- A vascular surgeon first selected the slices defining the

strates the low sensitivity of the proposed segmentation proximal and distal boundaries of the aneurysm and then

approachtonoise.Finally,toassesstheeffectivenessofour selected a sufficient number of points along the outer wall
Author's personal copy
210 computer methods and programs in biomedicine 107 (2012) 202–217

Fig.6– Artificialimages((a),(d)and(h))usedtotesttheflexibilityofthemethods.(a)ischaracterizedbyaregionwith

similarintensitybetween“thrombus”and“outerwall”.(b)showsthethrombusinitializationfromtheFMMand(c)depicts

theresultofourmethodsforthetwoboundaries.Inartificialimage(d),thesameintensitybetween“thrombus”and“outer

wall”hasbeensetandtherearetwoouterleakageregions.TheDLRmethodresultisshownin(e),wheretheleakage

regionsareinyellow(angle(cid:2)/45).(f)and(g)showthefinalresult(LDRCisused,30iterations)intheabsenceandpresence

respectivelyofwhiteGaussiannoiseofzeromeanandvariance=0.0001.Image(h)wasgeneratedbyremovalof4

“calcifications”fromimage(d).Thefinalresultforthewallboundaryisshownin(i).(Forinterpretationofthereferencesto

colourinthisfigurelegend,thereaderisreferredtothewebversionofthearticle.)

boundaryforeachslicewithinthatrangesothatdetailsofits datasets,tracingthethrombusandouterwallboundariesfor

shape are captured. Boundary contours were then obtained eachslice.Asuitablecroppingofeachslicewasmadebefore-

bylinearinterpolationoftheselectedpoints.Thisprocessis hand,selectingtheregionofinterest,toreducesegmentation

repeatedforthemorechallengingthrombusboundarybutin time.Togeneratetheresultswehaveusedtheparameterval-

thiscasethemedicalexpertwasrequestedtoprovideforeach ues presented in Section 4.1 and the LDRC as the stopping

sliceonlythosepointsthathewascertainoftheirpositionin criterionofchoice.

anefforttoreducetheuncertaintyofthemanualsegmenta-

tionresult. 4.5. Validationfortheouterwallboundary

Themanualsegmentationwasperformedbyapointselec-

tiontooloftheimageprocessingsoftwareImageJ[25].Wethen The quantities used to evaluate segmentation error for the

appliedoursegmentationmethodtoalltheseslicesoftheten outer wall boundary for a slice were: Absolute and Relative
Author's personal copy
computer methods and programs in biomedicine 107 (2012) 202–217 211

Area Error (AAE and RAE), Mean Distance (MD), Hausdorff
Table1–Summarizedresultsfortheouterwall:mean,
Distance (HD) and Area Overlap (AO). We also extracted the
standarddeviation,minimumandmaximumvalues
AbsoluteandRelativeVolumeErrors(AVEandRVE)foreach overtenp atientdat asetsarep rovi dedforthe quantities

patient.AAEandRAEaredefinedas usedtoevaluatethesegmentationerror.

Mean±s.d. [min,max]
AAE
AAE = |A(M) − A(A)| and RAE = A(M), (8) NAbos. oolfu stleicAerse aError(AAE,cm2) 415.1 ±47 7± 0.355 [[209.6, 5546,]1.729]
RelativeA reaE rror( RAE,% ) 6.0± 3.3 [2.5,13 .1]
where A(M) and A(A) are the areas of the manual and auto- Absolute Volu meE rror(A VE ,cm3) 4.5 63 ±2 .083 [0.23 7,7.079]
matic contours, respectively, measured in cm2. RelativeV olumeE rror( RVE,% ) 4.0± 3.7 [0.3,12 .4]
Mean and Hausdorff distances are used to quantify the MeanDi stance(M D,m m) 1.3 20 ±0.318 [0.97 8,1.910]
similarity of t he two cu rves. They are defin ed for the p oint Hausd orffDista nce( HD,mm) 4.160 ± 1.096 [2.656, 5.999]
clustersA an dB as AreaOverl ap(AO,% ) 94.6± 1 .8 [91.8,9 7.0]

MD(A,B)=max{m(A,B),m(B,A)}, (9) Table1summarizestheaverageresultsfortheouterwall

overthetenpatientdatasets.Themeanofaquantityinthis

HD(A,B)=max{h(A,B),h(B,A)}, (10) tableisthemeanofthetenM’s,i=1,...,10,whereM isthe
i i
meanofthisquantityovertheslicesoftheithdataset.Mini-

where mum and maximum values are also referring to these ten Mi’s.
Fig.7presentsanalyticallytheresultsforeachofthedatasets

m(A, B) = meana ∈ A{mb ∈i nB{d(a, b)}}, (11) fwoar ltlhbroeuen odfa trhyes aebgmoveen qtuatainotniteierrso ur,sneadm foerl yevthaeluMatDin,tgh tehHe Douatnedr

theAreaOverlap.Twosliceswiththemanualandautomatic
h(A, B) = maxa ∈ A{mb∈inB{d(a, b)}}, (12) out erwa llcontou rssu perim posed are shown inFig .8(a)and

(b)respectively.

whered(a,b)isaL2normforthepointsaandb.Themainissue

regardingthesedistancesisthattheyonlytakeintoaccount 4.6. Validationforthethrombusboundaryby

thesetsofpointsonbothcurvesanddonotreflectthecourse comparingwallthi ckn ess

ofthecurve.Weovercomethisproblembyfinding400points

on each of the two contours. These points are taken as the
Themanualsegmentationforthethrombusboundaryfollow-
intersections of the curve with rays which start at the cen-
ing linear interpolation of the selected points includes only
troid of the contour and each two successive ones form an
segmentsandnotaclosedcurveforsomeslices.Thisisdue
angleequalto(cid:4)/200.Inthisway,weconsideralargenumber
to our decision to include in the construction of the man-
ofpointswhichareuniformlydistributedalongtheseconvex
ualsegmentationboundaryonlythosepointsthattheexpert
or, in the worst case, star-shaped curves. We then calculate
is certain of their position (see Section 4.4). We introduce
theirMDandHD(inmm),makingthesequantitiestogethera
a method for approximating the mean wall thickness of an
veryreliablemeasureoftheresemblanceofthetwocontours.
AAA in a slice-by-slice fashion, which is used as a throm-
MDisameasureoftheirproximityandHDidentifiesanyout-
bus boundary validation. Following proper alignment of the
liers, which may be hidden in the averaging process of the
thrombusandouterwallcurves,wedrawraysfromthecen-
MeanDistance,byprovidingtheworstpossibledisagreement
troidofthewall,sothateverytwosuccessiveraysformthe
betweenthecurves.
sameangle(cid:5).Sixsuccessiveraysdefineasector,sothereare
AreaOverlapisdefinedas
2(cid:4)/5(cid:5)sectors(theangleofeachsectoris5(cid:5)and(cid:5)shouldbe

chosensothat(2(cid:4)/5(cid:5))∈N).Foraraywhichintersectsboththe
A(A∩M)
AO=2 , (13) thrombus and outer wall curves, we define Ray Wall Thick-
A(A)+A(M)
ness(RWT)tobethedistancebetweenthesetwointersection

points.ForasectorwedefineSectorWallThickness(SWT)to
whereA(A∩M)istheareaoftheintersectionofthemanual
bethemeanoftheexistentRWTovertheraysofthesector.
andautomaticcontours.Fortheperpatientresults,Absolute
ThemeanSWToverallthesectorsofalltheslicesofanAAA
andRelativeVolumeErrorsaredefinedasfollows:
scanprovidesanevaluationofthemeanwallthicknessofthe

AVE AAA.Fig.8(c)showstheresultsobtainedfromthemanualand
AVE=|V(M)−V(A)| and RVE= , (14)
V(M) semi-automaticsegmentationofthethrombusandouterwall

boundaries.
whereV(M)andV(A)arethevolumesofthemanualandthe Oneoftheadvantagesofthismethodforestimatingwall

semi-automaticsegmentation,respectively.Volumeforadata thicknessisthatsincewecanincreasethenumberofsectors

setiscalculatedastheproductofthesumoftheareasofthe byreducingtheangle(cid:5),wecancontrolthespatialresolutionin

curvesforalltheslicesofthesetandtheslicedistance[26].Let ouranalysisofthewallthicknessalongthethrombusbound-

usnotethatwedonotpresentthesevolumeerrorsasanerror ary.Furthermore,theresultisindependentofthenumberof

indicationoftheactual3Dvolume.Giventhefactthatthein- pointsofthethrombusandouterwallcurves.

slicedistanceisthesameforallslices,thevolumeerrorscould InmanualsegmentationSWTwasnotcalculatedinthose

beconsideredasaweightederrorofallslicesforeachpatient. sectorswerenoraycouldbefoundtointersectboththrombus