Intensity-based image registration using scatter search

13 Pages·2014·1.89 MB·English

by  Andrea Valsecchi

ArtificialIntelligenceinMedicine60(2014)151–163

ContentslistsavailableatScienceDirect

Artificial Intelligence in Medicine

journal homepage: www.elsevier.com/locate/aiim

Intensity-based image registration using scatter search

AndreaValsecchia,∗,SergioDamasa,JoséSantamaríab,LindaMarrakchi-Kacemc,d

aApplicationsofFuzzyLogicandEvolutionaryAlgorithmsResearchUnit,EuropeanCentreforSoftComputing,CalleGonzaloGutiérrezQuirósS/N,33600

Mieres,Spain
bDepar tmentofComputerScience,UniversityofJaén,EdificioTecnológicoydeIngenieríasA-3,ParajeLasLagunillasS/N,23071Jaén,Spain
cNeuroSpin,F ren chAlterna tiveEne rgiesandA to micE nergyC ommission,B â tim ent145,Ce ntre d’étude sde Saclay,91 191 Gif-sur -Yvet te,France
dCentredeR echerch edel’Instit utduCer veau etdela Moelle épinière,Hô pitalPitié Salpê trière, Boulevar dd el’Hôp ital47, 75013Paris,Fr ance

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

Articlehistory: Objective:Wepresentanovelintensity-basedalgorithmformedicalimageregistration(IR).
Receiv ed6November2012 Methodsa ndm aterials : TheIR problemisform ulatedas ac ontinuou sopti mizationtas k,andourwork
RAeccceepivteedd i2n 8 rJeavnisueadry fo2 r0m1 413 January 2014 focuses o n th e develop men t o f the optim i zation comp on en t. Our meth od is designe d ove r an adva nced
scattersearchtemplate,anditusesacombinationofrestartanddynamicboundarymechanismsinte-

gratedwithinamulti-resolutionstrategy.
Keywords: Results :Theex p erimentalvalidat ionisperformedovertwodatasetsofhumanbrainmagneticresonance
Globaloptimization
imaging.Thealgorithmisevaluatedinbothastand-aloneregistrationapplicationandanatlas-based
Heuristics
segmentationprocesstargetedtothedeepbrainstructures,consideringatotalof16and18scenarios,
Scattersearch
IAmtlaags-eb r aesgeidstsreagtimonentation rpeasrpiseocntivpeulryp. oFsivees, easntadbglirsohuendd I-Rtr tuetchhdnaiqtaueiss, ubsoetdh tfoeaqtuuaren-t iatantdiv ienltyenasssiteys-sbtahseedq,u aarleit cyoonfsitdheerreeds ufoltrs c.oOmu-r
Magneticre sonanceimaging approac hrankedfi rsti nbothstudies and it isable to outperforma llcomp etit orsin12 o f16 registra tion

scenariosandin14of18registration-basedsegmentationtasks.Astatisticalanalysisconfirmswithhigh

confidence(p<0.014)theaccuracyandapplicabilityofourmethod.

Conclusions:Withaproper,problem-specificdesign,scattersearchisabletoprovidearobust,global

optimization.Theaccuracyandreliabilityoftheregistrationprocessaresuperiortothoseofclassic

gradient-basedtechniques.

©2014ElsevierB.V.Allrightsreserved.

1. Introduction technology over the last decades determine a constant demand

forbetterimageprocessingandanalysistechniques.Dealingwith

In its most general formulation, image registration (IR) [1] is novel,morediverse,andincreasinglyaccuratesourcesofimaging

the task of aligning two or more images in order to establish a dataisthemainchallengeinIRanditexplainswhyitisstillavery

spatialcorrespondenceoftheircommoncontent.Suchimagesusu- activeresearchfield.

ally have the same or a similar subject but have been acquired The alignment between two images is specified as a spatial

underdifferentconditions,suchastimeandviewpoint,orbymul- transformation, mapping the content of one image to the corre-

tiple sensors. In medical image analysis, IR is a key technology spondingareaoftheother.ApopularstrategyamongIRmethods

that allows to “fuse” visual information from different sources istoperformthealignmentbyconsideringonlysalientanddistinc-

[2]. Applications include combining images of the same subject tivepartsoftheimage,suchaslines,cornersandcontours,called

fromdifferentmodalities,aligningtemporalsequencesofimages features. This strategy has the advantage of greatly reducing the

to compensate for motion between scans, image guidance dur- complexityoftheproblem,butreliesontheabilitytodetectthe

ing interventions and aligning images from multiple subjects in featurescorrectly.However,thisapproachislimitedtothecases

cohortstudies.Theremarkabledevelopmentsinmedicalimaging inwhichfeaturesaloneareabletocharacterizetheimagecontent.

IRmethodsfollowingthisapproacharecalledfeature-based[2,1],

whilethetermintensity-based(orvoxel-based)namesthemethods

inwhichthewholeimagedataisused.
∗ Correspondingauthor.Tel.:[phone]. Regar dles s of th is divi sion, th e core of every IR technique is

E-mailaddresses:[email protected], an optimization process that explores the space of geometrical
valsecchi.a ndrea@gm ail.com (A. Valsecchi), [email protected] tra nsformations .Twostr ategi esareava ilabl e.Inpa ram eters-based
(S.Damas),[email protected](J.Santamaría),[email protected]
approaches the search is directly performed in the space of the
(L. Marrakchi-Kacem).
http://dx.doi.org/10.1016/j.artmed.2014.01.006
0933-3657/©2014ElsevierB.V.Allrightsreserved.
152 A.Valsecchietal./ArtificialIntelligenceinMedicine60(2014)151–163

transformationparameters.Hence,asolutionisavectorofvalues

fortheparametersoftheregistrationtransformation.Inmatching-

based approaches, features are matched through a search in the

space of feature correspondences; once a suitable matching has

beenfound,thetransformationparametersarederivedaccordingly

bynumericalmethods.Inbothcasesthesearchisguidedbyasim-

ilaritymetric,afunctionthatmeasuresthedegreeofresemblance

betweentheinputimages.Thiscanbedoneeitherbycomparing

thewholeimagesorjusttheircorrespondingfeatures.Traditional

parameters-based methods use classic gradient-based optimiza-

tion algorithms, while matching-based methods use matching

algorithmslikeiterativeclosestpoint(ICP)[3].

ManyfeaturesofIRproblems,suchasnoise,discretizationand

differencesintheorderofmagnitudeofthetransformationparam-

eters still pose a challenge to traditional optimization methods.

A number of alternative approaches are based on metaheuristics

[4], which have proven their ability to deal with complex real-

world problems in a large number of fields, including computer

vision and image processing. In particular, metaheuristic-based

registration approaches have demonstrated to be a promising

solution to overcome the drawbacks of traditional optimization Fig.1. Theinteractionsamongthecomponentsofaregistrationtechnique.

algorithms[5,6].Scattersearch(SS)[7–9]isaprominentexample

ofsuchtechniques.Ithasalreadybeenappliedtoimageregistra- ThesecondcomponentofanyIRmethodisthesimilaritymet-
tio npro blemsasop ti mize roffea ture-b asedap pro aches. InIR,SS ric, a fu nction F(I , I ) tha t m eas ure s the de gr ee o f resembl ance
hasb eensucce ssf ulbothwh en usedtofindm atchingsam on gfe a- bet w eentwoi ma1ge s2.A sth efinalper form anceof an yIRmethod
ture s[10 ]aswellas insea rching fort he tran sformation parame ters depends onth eaccura tee stim ation ofthealignm en tof the images,
direct ly[1 1] . thisisac ruc ial IRcompo nent[13].T he qu alityofatr an sfor mation
Inth iswork,SSisusedasbaseforanovelintensity-basedIR fiso bt a inedby co mputingthe sim ilarit ymetri co v erthemodelI
meth od.T healg ori thm iss pec ifical lyd e signed totakeadvanta ge a nd thetran sfo rmedscene f(I ).Theactu alevalu ation m echanis mM
of the ch arac teristics o f t he IR proces s to impr ov e the optimiza- dep ends on the natu re of thSe r egis tration approach. In feature-
tio n.T oevaluateitsef fec tive nes s,ourm eth odiscom par edwithan basedme tho dst hesimil ari tym etricusually measurest he distance
heter og eneousg rou pofcompetit ors intwoe xp erimental stud ies betwe encorres pon dingfeatu res[14 ].Forin stance,th eal ignment
involvingsimu lateda nd realmedical im ages .Athorougha nalysis canbeev aluatedusingm eansqua reerr or(M SE)betw een thepoints
of the res ults is per form ed, and thei r signific an ce is asse ssed by ofth e modeland those ofth etrans forme dscen e.Ifthem od elhas
m eans ofdiffe re ntstatistical test s. rf eatu repoin ts,e achpo in tx i sassignedto theclo se stp ointin the
The pa perisstr ucturedas follows.InSection2,wereviewthe tr ansform edsce nec, andt hei M SEisgive n by:
imager egistra ti onproblem an dpresen ts everalte ch niq uestoso lve i
(cid:2)r
it. Section 3 introduces the basics of SS and the design of our IR 1
m ethod.In S ection4,w epr esentt he ex perim en talstud ies alo ng MSE = r (cid:3)xi− ci(cid:3)2
withthe an alysisof th eir results.F inal ly,conclusions arepro vided i=1

inSection5.
In intensity-based approaches, instead, the resemblance of the

intensityvaluesinthetwoimagesareconsidered.Sumofsquared

differences(SSD),normalizedcorrelation(NC)andmutualinforma-
2. Imageregistration
tion(MI)[15,16]aretypicallyused.Thesubjectoftheimagesand

theiracquisitiontechniquedeterminethekindoftherelationship
AtypicalIRprobleminvolvestwoimages,conventionallycalled
between the intensity distributions, which in turn decides what
model (I ) and scene (I ), with different roles in the registration
M S similaritymetricsareappropriate.Forinstance,whentwoimages
process. The model is the reference (or target) image, while the
havebeenacquiredusingdifferentsensors,ascenariocalledmulti-
sceneistheimagethatistransformedtoreachthegeometryofthe
modalregistration,therelationshipbetweentheintensityvaluesin
other.Theregistrationaimstofindageometrictransformationf
theimagescanbestronglynon-linear.WhileNCcanhandlealinear
thatalignsthescenetothemodel;inotherwords,fissuchthatthe
relationship,metricsbasedoninformationtheory,suchasMI,are
modelI andthetransformedscenef(I )areassimilaraspossible.
M S bettersuitedforthisscenario.MIisdefinedas
SeveralcomponentscharacterizeanIRmethod.Firstwehave (cid:2)
p(m,s,f)
ftohrem taratinosnfocrmanatbioenu mseoddetol, athliagtn dtehteerimmiangeess .wThhiicshc hkoinicde odfe ptreanndss- MI = p(m, s, f )log2pM(m) pS ( s,f)
s∈LS,m∈LM
entirelyontheconcreteapplication;verysimplemodelssuchas

translation transform can be enough in certain contexts such as whereL andL aresetsofregularlyspacedintensitybinscenters,p
M S

remote sensing [12]. At the other end of the spectrum there are isthediscretejointprobabilityandp ,p arethemarginaldiscrete
M S

non-rigid(alsocalledelastic)transformations,suchasB-splineand probabilitiesofthemodelandsceneimages.

thin-platesplinestransformations,abletorepresentlocaldeforma- ThethirdmaincomponentofanIRmethodistheoptimizer.Itis

tions(warpings)usinghundredsoreventhousandsofparameters. responsibleforfindingthebesttransformation,intermsofsimilar-

Othercommonchoicesincluderigidtransform,whichallowstrans- itymetric,amongthetransformationsinourtransformationmodel.

lationandrotation,similaritytransform,whichalsoadmitsscaling, Fig.1showstheflowchartofthewholeregistrationprocess.

andaffinetransformation,whichcanalsorepresentshearing.These Eachoptimizerhasadifferentsearchstrategy,whichdepends

areexamplesofglobaltransformationshavingrespectively6,7and also on the nature of the algorithm. One approach is to perform

12degreesoffreedomfor3Dimages. thesearchdirectlyinthespaceofthetransformationparameters.
A.Valsecchietal./ArtificialIntelligenceinMedicine60(2014)151–163 153

Table1 Adaptivestochasticgradientdescent. Adaptivestochasticgra-
The IR algorithms included in our experimental study. dient desce nt (ASGD) [31] is an optim ization me thod desig ned

Nature Strategy Optimizer for intensity-based IR. The algorithm is an extension of the

Robbins–Monrostochasticgradientdescentmethodinwhichsome
I-ICP Feature Matching Gradientdescent
Dyn- GA Feature Paramete rs Genetica lgorithms of the parameters are automatically computed, in particular the

SS* Feature Matching Scattersearch step size. In the comparison of intensity-based IR methods pre-
ASG D Intensit y Paramete rs Gradien tdescent sent edin [32 ],A SGDoutperf orm edtheotheral gor ithmsin both
GA+ Intensity Parameters Genetica lgorithms
SS+ Intensity Parameters Scatters earch affine and elastic registration problems. The algorithm uses a ran-
domimagesamplerandamulti-resolutionstrategyinwhichthe

imagesinthepyramidsareobtainedbyapplyingbothdownsam-

Thisturnstheregistrationinacontinuousoptimizationproblem, plingandGaussiansmoothing.

thereforeclassicnumericaloptimizationalgorithmscanbeused. Dyn-GA. Dyn-GA[22]isaparameter-basedIRtechniquebased

Gradientdescent,Newton’smethod,Powell’smethodanddiscrete onareal-codedGeneticAlgorithm(GA)[33].Thealgorithmuses

optimization[17]areamongthemostcommonchoicesalongwith fitness-proportionateselectionandanovelcrossoveroperatorthat

approaches based on evolutionary computation (EC) and other swapsanumberofgenesbetweenthetwoindividuals.Themuta-

metaheuristics[18–26].IRalgorithmsthatfollowthisapproachare tionreplacesthevalueofarandomlyselectedgene;thenewvalue

then called parameter-based. An alternative approach consists in isdrawnfromarangethatdependsonthefitnessoftheindividual:

searchingforamatchingbetweenfeaturesorareasoftheimage. thelargerthefitness,thelargertherange.Also,astheoptimization

Fromthematchonecanderivetheparametersofthecorresponding progresses,therangesoftransformationparametersarerestricted

transformationusingnumericalmethods.Thisclassofalgorithms aroundthoseoftheindividualsinthepopulation,sothatthesearch

iscalledmatching-based.Theiterativeclosestpointalgorithmisa isfocusednearthecurrentsolutions.
fa mous example followin g th e latter approac h [27 –29]. ICP u se s GA+. GA+ has been p roposed in [34] as an intensity-based

theclosest-pointassignmentrule(i.e.amodelpointisassignedto method whose optimizer is a GA. The registration is performed

theclosesttransformedscenepoint)toperformthematch.Once throughasearchinthespaceoftransformationparameters,there-

amatchisestablished,thecorrespondingtransformationiscom- fore the GA is real-coded and uses real-coded operators such as
p utedus in gleastsquar ese stimationorot hermorerobus tm odel BLX- ˛cr osso ve r.GA+canh and leanu mberoftra nsformatio nmo d-

fittingtechniques[30]. elsandsimilaritymetrics.Italsosupportmultipleresolutions.A

Finally, we mention two minor components that play a role restartmechanismisusedattheendofthefirstresolution.Ifthe

in the registration. In intensity-based methods, computing the fitnessofthebestindividualisgreaterthanafixedthreshold,the

similarity metric on the whole images is usually unfeasible and registrationmovestothesecondresolution,otherwisethefirstres-

unnecessary,thereforeasamplingstrategydetermineshowmany olutionisperformedagain.Experimentsshowedthismechanism

andwhichvoxelsareactuallyused.Thoseareusuallyselectedat helpsthealgorithmtoachieveabetterperformanceatalowextra

randomwithuniformprobabilityorsampledalongaregulargrid. computationalcost.

Second, it is common to perform the registration in multiple SS*. SS*[10]isamatching-basedapproachusingtheSSopti-

stages.Increasinglylargerandmoredetailedversionsoftheinput mization algorithm (see Section 3). The proposal is specifically

imagesareusedateachstage.Inthefirststage,theoptimizerfinds designedtoworkwithaclassoffeatures,crestlinespoints,and

asolutiontoacoarseversionoftheregistrationproblem.Ineach exploitstheknowledgeofthelocalcurvatureofthepointstoper-

ofthefurtherresolutions,partofthedetailsoftheinputimagesis formthematching.Inaddition,theauthorsproposedanadvanced

restoredandtheoptimizeraimstoadaptthesolutionoftheprevi- codingscheme,inwhichamatchingisrepresentedasapermuta-

ousphasetofitthenew,moredetaileddata.Themulti-resolution tionofpoints,andanoveldesignforsomeofthecomponentsof

strategydetermineswhichkindofprocessingisperformedonthe SS.

imagesineachdifferentstageoftheregistration;usuallythepro-

cedureincludesdown-samplingandsmoothing.Thesequenceof
3. Scattersearchdesignforintensity-basedIR
imagesusedduringtheregistrationiscalledpyramid.

In the rest of this section we present a selection of remark-
SSwasoriginallyproposedbyGlover[7]inthecontextofinteger
ableIRalgorithmsthatwillbelatercomparedinanexperimental
programming.ThemainideabehindSSistorecombinesystemat-
study.Thegroupisquiteheterogeneousintermsofnatureofthe
ically a reference set of high-quality solutions, possibly obtained
approaches,searchstrategiesandoptimizers(seeTable1),encom-
usingdifferenttechniques.MostSSimplementationstodatehave
passingmostofthecombinationsofapproachespresentedinthe
beenbasedontheSStemplatespresentedin[8,9];ourexposition
previous section. We begin with two modern versions of classic
followsthewell-knownfive-methodstemplateintroducedinthe
IRmethodsandthenmovetoalgorithmsbasedonECandother
latter.Thesemethodsare:
metaheuristics.
I-ICP. This extension of the original ICP proposal has been
introduc ed in L iu [28]. In com par ison wit h its predece ssor, I-ICP • adiversificationgenerationmethodthatgeneratesasetofdiverse
introduces the use ofco lli nearitytoass esst heq ualityofam atch- tr ialsolutions;
ing.Consid ertw op oi ntsx,yofth em odela ndt heircor re sp onding • an i mprovement method, used to enhance a solution, usually
poin tsx(cid:5),y(cid:5) in the scene .If t he ma tching isgo od,x (cid:5) andy(cid:5) should app lyinganheur isticmet hod;
bothbe cl oseto the straigh t line lpassingt hro ughx ,y.The distances • areferenc es etupdate methodtobuildandmaintainareference
d(x(cid:5), l)a ndd( y(cid:5), l)p rovidea ddit io nalinfo rmation o n the matching se tofsolut ion sselect edforth eir quali tyor diversity ;
that I-I CPe xploi ts tofindb ettersolut ions.Moreov er, the algorithm • asu bs etgenerat ionmeth od, byw hichse tso fsolutionsfromthe
has amor ecompl ex term inatio nmechan ism.Once the algorithm re ferenc esetarecr eated;
has co nverg edtoaso lution,atiny ,randomper turba tion isapplied • asolution com bin ationmethodthatcombinesasetofsolutions
toit ,thenthe alg o rithmcon t inues .Thisste pisrepeated u ntilthe in tooneo rmorenews olutions .

algorithm converges to the same solution found before the per-

turbation.ThishelpsI-ICPtotraverselocaloptima,increasingits An outline of the complete SS procedure [9] is shown in

robustness. Algorithm1.Atthebeginning,anumberofsolutionsarecreated,
154 A.Valsecchietal./ArtificialIntelligenceinMedicine60(2014)151–163

improved and stored in a temporary container P of size PSize. someofthetransformationparametersinphysicalunits,e.g.trans-

Thereferencesetisthenupdatedbyselectingthebestsolutions lations can be specified in millimeters rather than in number of

betweenPandthereferencesetitself.Next,theprocessentersa voxels.Analogously,differentsimilaritymetricsaresupportedand

loop(line8)thatiteratesthecoreSSprocedureuntilastoppingcon- itschoiceislefttotheuser.Inthealgorithm,theactualsimilarity

ditionismet,controllingthedurationoftheoptimization.Then,an metricvalueisdirectlyusedasfitness.

innerloopbegins.First,subsetsofsolutionsfromthereferenceset Our efforts have been devoted to the design of the opti-

arecreated.Eachsubsetiscombinedintoanewsolution,whichis mizer component. Our proposal is based on a SS algorithm

thenimprovedandstoredinanothercontainer,Pool.Next,theref- with 2-tier design. First, we present the five SS methods espe-

erencesetisupdatedwiththebestsolutionsamongthoseinPool cially designed to deal with IR as a continuous optimization

andthecurrentreferencesets.Ifnonewsolutionhasenteredthe problem.

referenceset,theinnerloopends.Finally,anewsetPiscreatedas Diversification generation method. We adopt an approach

inthebeginningofthealgorithm. based on frequency memory [35] to ensure the search space is

ApopularvariantofthecanonicalSSdesignisthe2-tierdesign exploredinauniformmanner.Therangeofeachtransformation

[9],inwhichthereferencesetisdividedintwotiers.Onetier,called parameter is divided in four sub-ranges of equal size and a fre-

thequalityreferenceset,storestheb mosthighqualitysolutions, quencycounterisassociatedtoeachofthem.Asolutionisbuilt
1

whiletheother,thediversityreferenceset,containstheb solutions intwosteps.First,foreachtransformationparameter,asub-range
2

havinghighdiversitywithrespecttothoseinthefirsttier.Each ischosenatrandomwithaprobabilityinverselyproportionalto

tierisorderedaccordingtoeitherqualityordiversity,andduring its frequency count. Then, the actual values of the parameters

thereferencesetupdateprocedurethenewsolutionshavingthe areselectedatrandomwithintheselectedsub-rangeswithuni-

highestvaluesofqualityordiversityareconsideredforinclusion formprobability.Finally,thecorrespondingfrequencycountersare

intothereferenceset.Theadvantageofthisdesignisthatitadds increased.

diversity to the search instead of focusing only around the best Subsetgenerationmethod. Ourmethodgeneratessubsetshav-

solutionsfound,avoidingprematureconvergence. ingtwoelements.First,thesolutionsfromthequalityreferenceset
are con sidered, yi eldin g b (b −1)/2 subs ets . Then, the divers ity
Algorithm 1 (A generic optimization procedure based on SS.). refe rencesetis used,so th1ato1the rb2 (b2−1)/2 subse tsa regener-

ated.Finally,wecombinesolutionsfromthequalityanddiversity
21 RPe←fe∅re;nceSet ←∅; refere nce set , ad ding b1b2 /2 sets to the re sult .
43 W xh←ileD |Piv| e<r PsSiifziec adtoionGeneration(); BLXS-o˛luctrioosns coovmerboinpaetrioatno mr[e3th6o],da. coTmom coomncbhinoeic seoilnutrieoanls-c wode eudseE vthoe-
5 x (cid:5)← ImprovementMethod(x); lutiona ryAlgorit hms.Two solu t ionsxand yare pr ovidedasin put.
6 ifx(cid:5)∈/ Pthenaddx(cid:5)toP For each position i of thei r encodin g, the o pera tor comp ut es the
78 eWnhdile ¬ stop c ond ition do vinatleu rev adli[= m |xiin−(x y,i|y a )n −d ˛ tdh,emn ai tx (rxan,ydo)m+˛l yd ]gwe nitehrautneisf oar mvaplureo bza ibni ltithye.
9 UpdatetheReferenceSetbyselectingthebestbsolutionsin i i i i
ReferenceS et∪ P; Thevalu ezisas sign e dto theith po si tion ofth eresulti ngsolution.
10 NewEleme nt s←true; The param e te r˛contro ls the wid thofthe ra nge sinwhich thenew
11 Pool ←∅; valu esaredraw n .
12 WhileNewElementsdo
13 Subs ets←SubsetG eneration(); Improvement method. Our improvement method is a local
14 NewEle me nts←false; searchbasedonacrossoveroperator.Toimproveasolutionx,we
15 WhileSubset s =/ ∅do perform cros sov er withara ndomlych os ensolutio n yfrome it her
16 S←p op(Sub sets ); therefer enceseto rP.T h isresultsi nanum berofn e wsol utions
11217908 enxxIdf (cid:5)← ←x(cid:5) S∈/ Iom Plopurotloi bovfne mtChoeenmntb MAindeadtth ixoo(cid:5)ndto((Sx P))o;;ol zba1ne,sd . t .ts h.o,e zlumwt,ih odonel pe aepmnrdoo nicneg gs xso, i nsy ,rt hez1pe , e a.a c.t t.eu, dzamlf ocrirs oa ts nhsouevnme cbro eonrpso eifdraie ttreoerrda ti nai os un ossue..Nt pTouhttee,

21 UpdatetheReferenceSetbyselectingthebestbsolutionsin thatifallthesolutionsgeneratedbythecrossoverareworsethan
ReferenceSet ∪ P ool; theo ri gin als olution,th eimprove me ntm ethodsim ply return sits
22 IfReferenceSethasnewelementsthen
2234 enNdewElement s ← true inpWut ex.use the “parent-centric” version of the BLX-˛ crossover
25 BuildanewsetPusingthediversificationgenerationand called PMX -˛[3 7].Incontrastw iththeB LX -˛,t herang esforthe
improvem e ntm eth od s; newso lutions param e tersare[ x −˛ d,x +˛d]an d[y −˛d,y + ˛d],
26 end i i i i
whichindeedresultsinsolutionsthatareclosertotheirparents.
27 end
ThisoperatorhasalreadybeenusedinIRandityieldedthebest

resultsinthecomparisonofmemeticapproachestorangeIRcarried

In what follows we present the design of the proposed IR outin[38].
meth od, ca lled SS+. The effectiv enes s of any SS im plementati on R ef erencesetupdatemethod. Theupdatepolicyaimstomain-

dependsonaproperdesignofthepreviousfivemethodsthatmust tainthehighestqualitysolutionsinthequalityreferencesetand

bespecifictotheoptimizationproblemathand.Aswereviewedin the most diverse solutions in the diversity one. The update is

Section2themaincomponentsofanIRapproacharethetransfor- performed in three steps. First, a new quality reference set is

mationmodel,thesimilaritymetricandtheoptimizer.Ingeneral, created using the b having the highest quality values in both
1

the choice of the transformation model and the similarity met- the current quality reference set and the pool. Then, the algo-

ricarehighlyproblem-dependent.Ourmethodhasbeendesigned rithm computes a diversity measure over each solution in the
to sup port m ultiple transformatio n m odels. S S+ i s able to han- pool a nd in the di versity re ference s et. Fin ally, the b so lut ions
2

dle all parametric transformation models, from the simple Euler having the highest distance values in both the pool and the

transformuptoB-andthin-platesplines.Regardlessofthecon- diversity reference set are placed in the new diversity reference

cretetransformationmodel,asolutionisstoredasavectorofreal set.

values, corresponding to its parameters. In medical imaging the The diversity value of a solution x measures the distance

imagesalsostorethephysicalsizeofthevoxels,sowecanencode betweenthesolutionandthequalityreferencesetS.Itisdefinedas
A.Valsecchietal./ArtificialIntelligenceinMedicine60(2014)151–163 155

theminimumofthesquaredifferencebetweenxandeachsolution resolution. Therefore, we can improve the intensification of the

yinthequalityreferenceset,asin search by limiting the range of the transformation parameters

inside this area. This approach is called dynamic boundary [22].
diversity(x,S)=min{distance(x,y)|y∈S} (1)
Letxbethebestsolutionfoundatthecurrentresolutionandl,u

(cid:2)n be,respectively,thevectorofthelowestandhighestvaluesinthe
1
distance(x, y) = n (xi− yi)2 (2) fcourrrtehnetn reaxntgreess. oTlhuetino, nthaer era[nl(cid:5),gue(cid:5)s] owf htheree transformation parameters
i=1
wherenisthenumberofcomponentsofasolution. li(cid:5) =max{xi−(ui−li)/(cid:2),li}
SS+ al so inc ludesan um berofspeci fic c omponentsthatexploit
ui(cid:5) =min{xi+(ui−li)/(cid:2),ui}
the features of the IR problem and are essential to improve the

performanceofSS. and (cid:2)>1 is the shrinking factor. Basically, dynamic boundary
Duplicatio nc ontrol. Acommonissueinreal-codedevolution- rest ric ts t he para meter ran ges aro und the be st solutio n, i.e. to a
aryalgorithmi sthatex trem elysimi larsol ut ionscanbe generated. (cid:2)-times sma llerinterval centere daroun dx.N otet hatevenf orm od -
Des pitebeing di ffere ntintheir represe ntation,s uch so lutionscan est valu es of (cid:2), this can drastica lly redu ce the size of th e s earch
besoclo setoe achother th atfor allpurposesthe yenc odethesa me spa ce,dep end in gon then umberoft ransform ati onpa ra met ers.For
ob jec t.This ca ncau sethe refe ren ce settocon taino nlyalm ost iden- instan ce,affinetra ns form ationha s1 2parameters,m eaningthat for
tical co pies of the sa me solution, wh ich both s trong ly focu s the (cid:2)=2,dyn amic boundaryreduc est he sizeofthes earchspa ceb ya
searc haroun d thec urrent bestsolu tionan dmak espointle ssthe use fa ct or of212=4 ,096.
ofmul tiplesol utio ns.Toa void thisbeh avio r,weu tilizead upl ica- The o verall structureofSS+,showingalsothespecificIRcompo-
tio ncontrol mechanis m. Weco nsid ertwosol utio nstob e different nents,i sgiven inAlgorit hm 2. Duplicati onco ntr olisimp le mented
only iftheir distance(Eq .(2) )isabove agi venthresh ol d. Thiscon- directl yi nside the SSproced ur e.Thealgorit hmbeg ins withthefirst
ditio n ische ckedbefo ret hein cl usion of asolu tioninase t(lin es6 resolutio nand sw itc hestothes econ dresoluti oninl ine8 .Du ring
and19 in Algorith m1). thechange ,dyn amicbou nd ary isappli ed(line10 ).I nthe fi rstres-
M ul ti-r esolution s trategy. Before the registration, the input olut ion, the restart m echanism i s used ( lines 2–7 ), a nd the best
imagesareprocess edapplyin gbothdo wn -samplingand Ga ussian solution sfou nddur ingthediffe ren trun sares aveda ndla ter used
smooth ing inorderto createtw oima gepyramids.Afe wo perations inthesec ondre solution .
intheoptim iz erallo w ittoha ndle thech angefrom on ere solutionto
th ene xt.Thetw orefe re nc esetsa rele ftunto uched ;as weencode d Algorithm2(PseudocodeofSS+,highlightingtheintegrationofthe

transformationsusingphysicalunitsratherthanvoxels,thesolu- generalSSoptimizationprocedurewithadvancedcomponentsspecific

tionsarenotaffectedbythechangeoftheimage’ssize.However, ofIR.).

thequalityofasolutionchangesbetweenoneresolutionandthe

next due to the new details that are now visible in the images, 1 Computetheimagepyramidsoftheinputimagesbyapplying

thereforeweneedtocomputethefitnessofthesolutionsagain. downsamplingandGaussiansmoothing;
Restart me chani sm . The m otiv ation be hi nd t he use of restart is 32 USasvee tdhSeo lfiurtsiot nlesv← el∅ o;f the imag e pyramids;
sfoimunpdlea. Avte trhyel eonwd- qouf athliety firtrsat nresfsoorlmutaiotino,n t.hRe eafilgnoinrigthsmuc mhitgrhatn hsfaovre- 45 FoRr ur n←t h1e TSoc Natu tmerbSeeraorfRchesotaprttism diozationprocedure;

mationisunlikelytoproduceagoodfinalsolution,thereforeitis 6 AddthecurrentsolutiontoSavedSolutions;
moreap p ropriatet op erforma g ainth esea rchfora suitablein iti al 7 end

8 Usethesecondleveloftheimagepyramids;
registrationbyrestartingthealgorithm.Tocheckwhetherthebest
9 UpdatefitnessvalueofSavedSolutions;
solutionobtainedattheendofthefirstresolutionisacceptableor 10 Applyd ynamic boun da ry;
not,one mightco nsi der tos et ath resho ldonitsfi tn essvalue(i .e. 11 Runth eScatter SearchplacingtheSavedSolutionsintotheinitial

itssimilaritymetric).However,thefitnessvalueofanappropriate referenceset;

solutiondependsontheactualcontentoftheinputimagesandis

hardtopredict.Weadoptedanalternativeapproach.Thefirstres- 4. Experimentalstudy
oluti on isperfor med afixed num beroftim es,indepen den tlyo fits

outcome.Attheendofthisprocess,thebestsolutionsfoundare Theaimoftheexperimentationistocarryoutanobjectivecom-
considere dfo rth ese con dre solution. parison ofo ur SSp roposalandtheo th er state- of-t he- artIRmet hods
Thecom put atio nalcost ofperformingarestartatthefirstreso- describe d inS ec tion2.To tha ten d,we designedtwoe xp eriments
lutioni squitelow,ast heim a gesinvolved a restill sm all. Note that involving sy nthetica n dr eal-w orld me dicalima ges. Tomakethe
ifpist he prob abilit yo fob taining alowqu ality sol utiond uring the compariso n as obje ctive as possibl e, in add ition to per formin g a
fir s t r eso lution, perfo rm ing n res ta rts reduces this pro bability to visualcomp ari sonofthe res ults,thee ffe ctiveness of eachmetho d
pn+1 ,i.e.anexpo nentialdecr ea se.This observat ionp rovidesaru le isasse ssedusinga qu anti tativeva lida tionmeasure sp ecific toeach
ofth umb to estimatean appropria ten umberofres tarts.Tob o und ex perimen t.Furth e rmore,asm ostofthea lgorithms involve da reof
to thepro ba bilityofo bta ininganina ppropria te solution toa value non-determi nisticnature, we carri ed out anumbero findepe nde nt
r,j ust set runsoneachscena rio.Ou ran alysisin ves ti gatesse ve ralaspectsof

n≥llnn((pr))− 1 othne e raecshu sltcse.n Fairrisot ,b wy ec ommepaustuirneg tmhee apne arfnodr mstaanndcea rodf dtehvei aatligoonr iotfh tmhes
validationmeasureandrankingthealgorithmsaccordingly.Next,

Dynamic boundary. This component is used to further take weassesstheoverallperformanceofthealgorithmsintwoways:

advantageoftheuseofmultipleresolutions.Theinitialrangesof bycomputingtheper-scenariomeanrankofeachalgorithmand

thetransformationparametersarehighlydependentontheappli- by counting the number of scenarios in which one outperforms

cation, therefore they should be provided by the user. However, another,calledwins.

aswechangefromoneresolutiontothenextone,wecanassume In the last part of the analysis, statistical tests are performed

thattheoptimalsolutionforthenextresolutionliesinthesame to determine which results are significantly different. We used

areaofthesearchspaceasthebestsolutionfoundatthecurrent thetestsandtheproceduresrecommendedin[39]forcomparing
156 A.Valsecchietal./ArtificialIntelligenceinMedicine60(2014)151–163

Table2 Table3

ThecharacteristicsofthefourbrainMRIimagesusedinthefirstexperimentalstudy Parametersofthesimilaritytransformationsweusedintheexperiments:rotation
(Fig . 2). The noise v al ue r epre sents the r atio bet ween th e s tand ard deviation of the angle ((cid:3)), ro ta tion axis (ax, a y, az), translation ve ctor ( tx, ty, t z) and uniform scaling
whiteGaussiannoiseandthesignalofthebrightesttissue.Thenumberofcrestline factors.

points,usedasfeatures,isalsoreported.
(cid:3) ax ay az tx ty tz s
Image Lesion Noise #offeatures
III123 NNYeoos N11%%o ne 533894 338 TTTT1234 112216375856 .9 −−−0000....683876073632 −−0000....282490538096 −0000....465237015168 −−1216662 −15555....5555 −−−−24444....6666 1110..28
I4 Yes 5% 248
It is important to remark this difference: while the input of

intensity-basedmethodsconsistsofthewholeimagesdata(incase,
algorithms over multiple problems. We used non-parametric tests twoimagesmad eof60× 181×2 17 =2, 356,6 20voxe lsha vin gan
toavoidmaking(ortesting)anyassumptionaboutthedistribution
of there sults.Th ep erforma nceo fSS+iscom pared wit hthatofthe 8-bit intensity value), that of feature-based approaches is a set of
justafewhundredpoints(Table2).
remainingalgorithms(i.e.,amultiplecomparisonagainstacontrol
SixteenIRprobleminstanceswerecreatedbychoosingpairsof
method),aprocedurethathasmorepowerthanapairwisecom-
differentimagesamongthefouravailableandapplyingoneofthe
parisonofallalgorithms.ThetestsweusedareNemenyi’stest[40]
foursimilaritytransformationsshowninTable3.Similaritytrans-
andsigntest.ThefirstisaposthocprocedureofFriedman’srank
formationsinvolverotation,translation,anduniformscaling.The
sumtest[41]andisbasedontheranksofthealgorithms.Sign-test,
parameters values of the transformations were chosen to obtain
instead,comparesthealgorithmsusingonlythenumberofwins
largechangesintheobjectlocation,orientationandscale.Changes
andlosses.Asmultiplecomparisonsareperformed,thep-valuesof
ofsuchmagnitudeareusuallychallengingforIRalgorithms.The
thetestshavebeenadjustedusingHolm’smethod[42]inorderto
con trolt hefam ily-w iseerror rate. scenarios we considered in the experiments are I1versus Ti(I2), I1
Fora llal gorithms,w eused theoriginalimplementationsbythe versus Ti(I3), I1versus Ti(I4) and I2versus Ti(I4), for i = 1, 2, 3, 4.
author s. I -ICP, Dyn-GA an d SS* are written in C, while SS+, GA + and comInp oonrdenert stion thriogdhulicgehdt itnhSee ccrtuiocnia3l oronleth pelapyeerfdo rbmy atnhcee sopfeScSifi+,cw IRe
ASGDareimplementedinC++.Theirrunningtimescanbefairly
first present the preliminary study in Fig. 3. We selected a sim-
compared, as no bias due to the programming language or the
pleandacomplexscenario,intermsofqualityoftheimagesand
executionenvironmentisintroduced.Allexperimentshavebeen
magnitudeofthetransformationsinvolved,andplottedtheresults
performedusinganIntelCorei5-2400CPUwith4GBofmemory.
SS+was integr ate din Elasti x[32],a toolb oxfo r inte ns ity-based obtained incrementally enabling all the components. In both sce-
narios,eachcomponentbringsanoticeableimprovementonthe
medicalimageregistration.Elastixisfree,open-sourceandithas
results,whichcanspanacrossafewordersofmagnitude.Inpartic-
been used in over one hundred publications in medical imaging
ular,restarthasamajoreffectonthevariabilityoftheresults,sois
[43].Thesoftwareisbuiltontopofthepopular“InsightSegmen-
beneficialintermsofrobustness,whiledynamicboundaryaffects
tationandRegistrationToolkit”(ITK)[44].
mostlytheaccuracy.Ananalogouspatternwasobservedthrough-
outthe stu dy.Note that ,excludin gthispr elim inarystu dy,SS+ is

4.1. Firstexperiment:registrationofsimulatedbrainmagnetic alwaysusedwithallthecomponentsenabled.
reso nance images Inor dert opro vid eau niformcomp arisonofSS+withrespectto
ourre cent res ults[34 ], wecons ideredthesa m ealgorit hms:GA +,

Thefirstexperimentissimilartotheonescarriedoutin[10,34]; ASGD,SS*,I-ICPandDyn-GA.Theparametersettingsweusedfor

thecurrentproposalextendsthestudyoffeature-andintensity- thesemethodsaretheonescorrespondingtothebestconfigura-

based methods performed in the two previous publications. For tionsofourpreviousstudy.Theyweremanuallyadjustedforthis

thisexperimentweusedfoursimulatedbrainmagneticresonance experiment, starting from the recommended/default values pro-
ima ges(MRIs)fr om apub licda tabase.At otalo f16regist rationsce- videdbyaut hors.As forSS +,th eparametersvalueswe readju sted

narioswereartificiallycreatedbyapplyingtotheimagesasetof through a preliminary study, that was performed using an addi-

fourlargetransformations.OnthoseIRinstances,weperformeda tionalregistrationscenarionotincludedintheexperimentation.
com pariso nconsideringal arg e,hete ro geneousg rou pofIRalgo - In SS+ theregis trationi spe rformedi n two resolutions;atthe

rithms. first resolution the images are smoothed (Gaussian smoothing,
The images used in this experiment were obtained from the (cid:4)=4 )anddown sam pledbya fact orof4inea chdimens ion.Thefirst

Brainweb database at McGill University [45]. Brainweb provides resolutionisrepeatedtwotimes(i.e.onerestart)independentlyof

simulatedbrainMRIalongwithground-truthdata,thereforeitcan theresults.Thealgorithmusedthesameconfigurationinbothres-

beeasilyusedtoevaluatetheperformanceofvariousimageanaly- olutions:fiveelementsinthequalityandinthediversityreference
sis metho ds.In d eed,Brain we bhasbeenfre qu entlyu sedby theIR set,PSize sett o32,blen df acto r˛=0.3 for th eso lutioncom bination
res earchcom munity [46].Toc reat escen arioswith differ en tdif fi- met hoda nd 50 ite rations ofPM X w ith ˛= 0.5 forthei mprovement

culties,weaddednoiseandmultiplesclerosislesionstosomeof method.Mutualinformationwasusedassimilaritymetric.

theimages,asdetailedinTable2.TheimagesareshowninFig.2; Forallalgorithms,thetransformationmodelissimilaritytrans-
eac himage ha ssize60 ×1 81×2 17 vox els. form,a nd thetransfo rma tionparameters range sa re[−30,3 0]for

This experiment compares both feature- and intensity-based thetranslationcomponentand[0.75,1.25]forthescalingfactor.No

algorithms, thus some features need to be extracted from the restrictionwasappliedtotherotationaxisortoitsmagnitude.In
imagestop rovide anin putforfe ature- bas ed algorithms .Asin the SS+,atrans form isrepre sen ted asarealv ecto rw ith se venelemen ts:

original comparison, we used crest line points, i.e. points where threespecifytheversoroftherotationv,threethetranslationtand
the surf ace normal h as a sha rp va riatio n, dete cted throu gh the oneth escalin gs. Validso lu tion srequir ev x,vy, vz∈ [−1,1]and s >0.

approachdescribedin[47].SS*alsousestheprincipalcurvatures Thestoppingcriterionisrunningtime.Itischallengingtodesign

ofthecrestlinepointsasheuristicinformationtoguidethefea- afaircomparisonbetweenalgorithmshavingdifferentinputs,in

turematching.Therefore,thisinformationhasbeencomputedand particularinputswithverydifferentsize.Indeed,recallthatwhile

providedtothealgorithm. intensity-basedmethodsusethewholeimagesdata(oratleasta
A.Valsecchietal./ArtificialIntelligenceinMedicine60(2014)151–163 157

Fig.2. ThefourMRIbrainimagesusedinthefirstexperiment.Fromlefttoright,axial,sagittalandcoronalviews,alongwiththecorrespondingcrestlinepointsusedas

features.
0
0
0
0
1
0
0
0
2
E
MS 00
5
0
0
1
0
4
none dupl. multi restart dynamic none dupl. multi restart dynamic
control resol. bound. contro l resol. bound.
FFrigo.m 3.l eTfhteto efrfiegchtt o,ft htheec sopmepciofince InRt scoamrepeonnaebnltes ionf aSnS+inoncr ietms peenrtfaolrfmasahnicoen. ,Teh.ge. rtehseulbtso xrepfleort tloa bthelee dfir“sdt yenxapmeriicmbeonutnald satruyd”ys,h socwensatrhioes rIe1svusl.t sT1o(fI2S)S (+leufts)i nagndd uI2pvlisc.a Tt4io(In4)c (ornigthrot)l.,

multipleresolutions,restartanddynamicboundary.Logarithmicscaleisusedtoshowtheimprovementbroughtbyeachcomponentdespitethedifferencesintheirorder

ofmagnitude.
158 A.Valsecchietal./ArtificialIntelligenceinMedicine60(2014)151–163

0
5
3
300 I-ICP 1500 I-ICP
Dyn-GA
50 SS* Dyn-GA
2 SS*
0
MSE200 100
0
5
1
0
00 50
1
0
5
0 10 20 30 40 50 60 0 10 20 30 40 50 60
Time (s) Time (s)

Fig. 4. The MSE scored by I-ICP, Dyn-GA and SS* against running time. Results are averaged over 15 runs and refer to scenarios I1vs. T1(I2) (left) and I2vs. T4(I4) (right).
largesubset),onlyasmallfractionoftheinputdataisactuallyused particularalgorithmswithdifferentsimilaritymetrics,wecannot

infeature-basedtechniques.Inthiscase,intensity-basedmethods simplycontrasttheirvalues.Instead,wehavetoagreeonacom-

use a random subset of 25,000 voxels to compute the similarity monmeasuretoevaluateallsolutions.WeusedtheMSEoverthe

metric,whiletheaveragenumberoffeaturesintheimagesis393. crestlinepoints.Forthefeature-basedalgorithmsinthecompar-

Thisgivesusaproportionofroughly60,whichisusedtosetthe ison, this is simply the similarity metric used by the algorithms.

twotimelimits:20sforthefeature-basedapproachesand20min Thesolutionsfoundbyintensity-basedalgorithmswereevaluated

fortheintensity-basedones.Inaddition,westudiedthebehaviorof inthesameway,i.e.byapplyingtheobtainedtransformationto

thealgorithmsovertime,toensurethesetimelimitsallowthemto the scene’s features and computing the MSE with respect to the

reachconvergenceanddeliverrepresentativeresults.Fig.4shows model’sfeatures.Weexpectthischoicetointroduceasmallbias

examplesofthisanalysisforthefeature-basedmethods,whilethe infavoroffeature-basedalgorithms.However,usingasimilarity

resultsforintensity-basedonesarereportedinFig.5a. metricbasedonintensitiesmightfavorintensity-basedmethods,

thereforeasweareproposinganalgorithmfromthelatterclass,it

4.1.1. Results seemsmoreappropriatetofavorthecompetitorsratherthanour

As in previous works with this dataset, for each registration approach.

scenario we performed 15 independent runs of each algorithm. Table4reportstheresultsofthefirstexperiment.Foreachsce-

Sincewearedealingwithalgorithmsofdifferentnatures,andin nario,wereportedmeanandstandarddeviationoftheMSEvalues

(a)
0 0
50 GA+ 50 GA+
SS+ SS+
0 0
0 0
4 4
0 0
SE 30 30
M
0 0
0 0
2 2
0 0
0 0
1 1
0 5 10 15 20 25 0 5 10 15 20 25

Time (m) Time (m)

(b)
9 9
0. ASGD 0. ASGD
GA+ GA+
SS+ SS+
8 8
0. 0.
p
a
verl 0.7 0.7
O
6 6
0. 0.
5 5
0. 0.
0 2 4 6 8 10 12 0 2 4 6 8 10 12

Time (m) Time (m)

Fig.5. TheresultsofASGD,GA+andSS+againstrunningtime.Resultsareaveragedover15runsandrefertothefirsttwoscenariosinthefirst(toprow)andsecond(bottom

row)studies.ThescoreofASGDinthetoprowisabove30,000andithasbeenomitted.
A.Valsecchietal./ArtificialIntelligenceinMedicine60(2014)151–163 159

Table4 Table4(Continued)

Detailedresultsofthefirstexperiment.Foreachscenario,thetablereportsthe
average M SE, sta nd ard devia tion and rank ing of the algorithm s in the c omparis on. Algorithm MMSeEan SD Rank

Algorithm MSE Rank
12 ASGD 885.8 356.7 6
Mean SD
Dyn-GA 414.8 258.2 4
1 ASGD 61,816.4 795.7 6 GA+ 47.3 0.4 2

Dyn-GA 194.9 50.5 4 I-ICP 416.6 5
GA+ 36.4 0.3 1 SS* 153.9 86.1 3
I-ICP 344.4 5 SS+ 47.2 0.0 1

SS* 37.0 1.5 3
SS+ 36.7 0.0 2 13 ASGD 56,932.0 568.6 6
Dyn-GA 179.8 59.5 3
2 ASGD 34,773.9 238.0 6 GA+ 35.0 0.2 2

Dyn-GA 107.5 52.1 4 I-ICP 237.6 5
GA+ 36.7 0.4 1 SS* 193.1 62.0 4
I-ICP 130.7 5 SS+ 35.0 0.1 1

SS* 43.4 3.6 3
SS+ 36.9 0.0 2 14 ASGD 31,521.1 6.6 6
Dyn-GA 105.7 50.8 4
3 ASGD 111,870.0 154.6 6 GA+ 30.7 0.3 1

Dyn-GA 211.0 137.3 3 I-ICP 341.3 5
GA+ 1736.5 6216.7 5 SS* 74.9 41.1 3
I-ICP 894.3 4 SS+ 31.0 0.1 2

SS* 63.2 2.9 2
SS+ 41.1 0.0 1 15 ASGD 112,134.4 1027.4 6
Dyn-GA 192.2 115.8 3
4 ASGD 1233.3 166.8 6 GA+ 1104.8 5128.7 5

Dyn-GA 302.0 121.4 4 I-ICP 608.8 4
GA+ 32.7 0.2 2 SS* 103.8 66.6 2
I-ICP 631.7 5 SS+ 40.4 0.0 1

SS* 53.9 2.6 3
SS+ 32.7 0.0 1 16 ASGD 512.8 233.2 5
Dyn-GA 298.1 144.8 4
5 ASGD 61,063.5 309.1 6 GA+ 29.5 0.3 2

Dyn-GA 299.3 144.1 4 I-ICP 1587.8 6
GA+ 51.4 0.2 2 SS* 150.2 78.3 3
I-ICP 517.7 5 SS+ 29.2 0.0 1

SS* 112.2 12.4 3
SS+ 51.4 0.0 1

Table5
6 ASGD 34,796.2 223.8 6 Firste xperiment:resultofNemenyi’stestcomparingSS+withtheremainingalgo-
Dyn-GA 154.0 114.2 4
GA+ 43.8 0.2 1 rithms. The table reports the average rankings of the algorithms and the adjusted
p-valueforeachcomparison.
I-ICP 330.3 5

SSSS*+ 5464..70 40..15 32 Algorithm Mean rank p-Value
SS+ 1.25
7 ASGD 110,131.2 1022.7 6 GA+ 2.38 0.0137
Dyn-G A 326.5 174.0 3 SS* 2.75 0.0020
GA+ 1091.8 4965.8 5 Dyn -GA 3.62 0.0000
I-ICP 437.8 4 I-ICP 4.81
SS* 63.8 46.2 2 ASGD 5.94
SS+ 56.2 0.0 1

8 ASGD 1017.4 252.7 6

Dyn-GA 354.3 146.9 4 obtained by the algorithms along with their ranks. The average
GA+ 44.5 0.3 2 ranks (Ta ble 5) and the cou nt of wins (Table 6) pro vide a basic
I-ICP 478.0 5
view of the results of the comparison. We also include a visual
SSSS*+ 14242..37 80..02 31 comp ari son ofthea ver age resultsinthe thir dsce nario(Fi g. 6).

FromthehighesttothelowesterrorrankingisASGD,I-ICP,Dyn-
9 ASGD 58,146.8 661.0 6 GA, SS*, GA+ and SS +. ASG D scor ed th e larges t M SE va lues in all
Dyn-GA 255.4 228.2 4
GA+ 52.9 0.3 2 but one of the sc enari os. Its perform an ce varie s gre atly de pen d-

I-ICP 704.3 5 ingonthescenario,butingeneralthemeanMSEisatleastone
SSSS*+ 15823..96 330..10 31 ord er o f m agnitude awa y fr om the best soluti ons. I-I CP delive red
abetter,moresteadyperformance,butstillwithverylargeMSE
10 ASGD 35,695.3 2465.3 6 v alues. C ompa ring th e mean values , Dy n-G A sco red b etter than
Dyn-G A 163.1 57.5 3
GA+ 476.4 3648.3 4

I-ICP 1493.2 5
Table6
SS* 89.2 40.8 2
SS+ 46.7 0.0 1 First experiment: the number of scenarios in which the algorithm on the row has a
bettermeanMSEvaluethanthatonthecolumn.

11 ASGD 111,384.4 574.2 6
ASGD Dyn-GA GA+ I-ICP SS* SS+
Dyn-GA 224.9 87.3 3
GA+ 2823.9 7863.5 5 ASGD – 0 0 1 0 0
I-ICP 951.3 4 Dyn-G A 16 – 5 16 1 0
SS* 82.2 45.1 2 GA+ 16 11 – 12 11 4
SS+ 57.7 0.0 1 I-ICP 15 0 4 – 0 0

SS* 16 15 5 16 – 0
SS+ 16 16 12 16 16 –
160 A.Valsecchietal./ArtificialIntelligenceinMedicine60(2014)151–163

4.2. Secondexperiment:atlas-basedsegmentationofreal-world

MRIs
In the second experiment we used real brain MRI images

withoutapplyinganytransformation.Theregistrationisusedto

perform atlas-based segmentation of deep brain structures [48].

Thequalityofthesegmentationobtainedinthisphaseisusedto

assesstheeffectivenessoftheregistrationmethods.

Atlas-basedsegmentationisaprocedurethataimstoautomat-

icallydelineatearegionofanimageusinganatlas(oran“average”

image)ofasimilarsubjectinwhichthedesiredregionhasbeen

alreadysegmented.Thefirststepistoregistertheatlas(thescene)
Fig.6. Firstexperiment.Avisualizationoftheoverlappingbetweenthescene(yel-
totheinputimage(themodel).Thetransformationresultingfrom
luoswed) atondc rtehaete mthoedefil g(ubrleuse)a. rTehteh ofisgeuhreasv irnegfetrh teoc tlohsee sthtiMrdS Escteonathrieo.m Tehaen sroelsuutilotnosf th isph aseis thenu sedt ooverla pthe segmentedreg ionofthe atlas

thecorrespondingalgorithm.InparenthesisistheMSEvalueofeachsolution.(For tothescene.Theregionofthescenethatoverlapsthesegmented
inte rpretationofth ereferenc es tocolorinth is figu rele gend, the read erisrefe rred reg ion ofthe atla sisthe re sult ofthe segm entation pro cess.Fig.7
tothewebver sio nof thearticle .)
illustratestheprocess.Often,atlas-basedsegmentationisusedas

preliminarystepinamorecomplexsegmentationapproach.
Table7 Thirteen T -w ei gh tedb rainMRI wereretrieved fromtheNMR
Firste xperiment:resultofsigntestcomparingSS+withtheotheralgorithms.The 1
database[49].Thedeepnucleistructuresineachimagehavebeen
tableliststhealgorithmsalongwiththeirnumberofscenariosinwhichtheyhave
been outp erfo rmedbySS+ (Table 6,bo ttom row)and th eassociat ed adjuste dp-v alue. manuallydelineatedbyanexpertinordertocreatetheground-
truthdata usedtoeva lua te theregi str ation. Fig .7show so neofthe
Algorithm Losses p-Value
imagesalongwiththecorrespondingdeepnuclei.Eighteenregis-

ASGD 16 0.0002 trationscenarioswerecreatedbyselectingpairsofdifferentimages
DGAyn+-G A 1162 00..00070628 at rand om.
Giventhenatureofthisexperiment,wecompareonlyintensity-
IS-SIC*P 1166 00..00000022 basedalg orit hms:SS +, GA+ andASGD.T he transform ation modelis

affinetransform,whichinvolvesrotation,translation,scalingand

shearing,anditcanberepresentedusing12realparameters.Affine

I-ICPinallscenarios,withlessvariabilitybetweendifferentsce- transformisapopularchoiceinregistrationofmedicalimages[50].

narios,butthegapwiththebestresultsislargenevertheless.Also, Itisflexibleenoughtorepresentawiderangeoftransformations

thestandarddeviationisoftencomparablewiththemean,mean- anditdoesnotproduceanatomicallyunrealisticresults,asitcould

ingthatthequalityofthesolutionsvariesaloteveninthesame happenwithdeformablemodels.Anaffinetransformisspecifiedby
sce nario .SS* andGA + deli veredsim ilarpe rfo rm ances in ter msof alinear trans formation( a3×3m at rixA= {a })and a translati on
ranking: 2.75 and 2.38, respecti vely. SS * scored const an tly qu ite v ectort .InSS+,bothele m en ts a recom bi ne dii,jna rea lv ector(a ,
1,1
closetot hebe stre sults, rankingthird or betteri nallbuton esce- a1,2,.. .,a 3,2 ,a3,3 ,tx,t y,tz)with 12e lements.Th e va lidr angeof each
nario .G A+ h ada lesssta blebeha vior: ite itherr ea che dv eryg ood matr i x e leme nta i sc om puted fr omthelim itsi mpos edove r rota-
i,j
solutions(11scenarios)orfailedtoconvergeandscoredareally tion,scalingandshearing.Astheconcreteparametervaluesofthe

highMSEvalue(theremaining5ones).OurSSapproachhasthe optimalregistrationtransformationsarenotknown,weestimated

lowestrank,1.25.Itfoundthebestsolutionsinmostofthescenarios parameters intervals considering a big enough range to include

(12over16)andcameextremelycloseintheothers,withdiffer- allregistrationsolutionsforthisapplication.Weallowedrotations
enc esbel ow 1.0.T heMS Evaluesva ryina sm all range,b etwee n29.2 bet ween−90◦a nd90◦,sc alin gin therange[0 .9,1 .1],shear inginthe
and57 .7. interval[ −0.1,0.1] and transla ti ons betwe en− 15cm and15 cm .
T able 5 reports the p-value of Nemenyi’s test comparing SS+ ForG A+ and SS+ wek eptthesame configur ation us edin th efirst
againstG A +,SS*an dD yn-GA.W e includedon lyt hebestrank ing experim ent.For ASGD wete sted seve ralconfigurati onsv ar ying the

algorithmstoavoidloweringthepowerofthetest.Inallthreecases numberofresolutions(2,3and4)anditerations(500,1000and
thetestcon fir msth eperform anc eofSS+ is sig nifica nt ly better than 2000).In w hatfollowsw e re port on lyth eresultsob taine dwith the
that oft hecompe tito rs,withtheh igh estp- valuebeing thatof GA+, bestco nfi gurat ion,wh ich uses4 resol utio nsand 1000iter ation s.

0.0137.Thesign-testcomparingthenumberofwins(Table7)has The stopping criterion is a time limit of 10min. Note that no
similar resu lts. All al gorithms b ut G A+ have a p-val ue of 0 .00 02, synthet ictransfo rmationi sa pp liedo nthe im age s;th eregi strat ion
whileG A+’sisc ons iderablyhi gher ,0.0768.T h etestco nfi rmsthe aimstoco mpensateforth ev ariabili tyin the poseof the patientdur-
differe nce b etw een SS+ and the ot hers alg orith ms is signific ant, ingth e acquisitiono fth eim ages.Ther efo re, with re spec ttothe first

althoughtoalowerdegree. experiment,the“magnitude”ofthetransformationstobefoundis

Fig.7. Theatlas-basedsegmentationofdeepbrainstructuresinbrainMRI.First,theatlasisregisteredtotheinputimageandtheresultingtransformationisappliedto

thelabeledregionoftheatlas(inblue).Theresultingregionisoverlappedontheinputimage(inyellow)todeterminetheoutputoftheprocess.(Forinterpretationofthe

referencestocolorinthisfigurelegend,thereaderisreferredtothewebversionofthearticle.)