<%BANNER%>

Mathematical Models and Fast Numerical Algorithms in Medical Image Analysis

Permanent Link: http://ufdc.ufl.edu/UFE0042919/00001

Material Information

Title: Mathematical Models and Fast Numerical Algorithms in Medical Image Analysis
Physical Description: 1 online resource (121 p.)
Language: english
Creator: YE,XIAOJING
Publisher: University of Florida
Place of Publication: Gainesville, Fla.
Publication Date: 2011

Subjects

Subjects / Keywords: CONSISTENT -- CONVERGENCE -- MRI -- OPTIMIZATION -- RECONSTRUCTION -- REGISRATION -- STATISTICS
Mathematics -- Dissertations, Academic -- UF
Genre: Mathematics thesis, Ph.D.
bibliography   ( marcgt )
theses   ( marcgt )
government publication (state, provincial, terriorial, dependent)   ( marcgt )
born-digital   ( sobekcm )
Electronic Thesis or Dissertation

Notes

Abstract: A series of projects are presented to provide advanced mathematical models and numerical algorithms that improve the accuracy, robustness and efficiency of the compressed sensing reconstruction technique in magnetic resonance imaging and inverse consistent image regisration. Chapter one introduces a novel variational model that enforces the sparsity of the underlying image in terms of its spatial finite differences and representation with respect to a dictionary. The dictionary is trained using prior information to improve accuracy in reconstruction. In the meantime the proposed model enforces the consistency of the underlying image with acquired data by using the maximum likelihood estimator of the reconstruction error in partial k-space to improve the robustness to parameter selection. Moreover, a simple and fast numerical scheme is provided to solve this model. In chapters two and three we develop fast numerical algorithms for solving total variation and L1 (TVL1) based image reconstruction with application in partially parallel MR imaging. Our algorithms use variable splitting method to reduce computational cost. Moreover, the Barzilai-Borwein step size selection method is adopted in our algorithms for much faster convergence. Theoretical and experimental results on clinical partially parallel imaging data demonstrate that the proposed algorithm requires much fewer iterations and/or less computational cost than recently developed operator splitting and Bregman operator splitting methods, which can deal with a general sensing matrix in reconstruction framework, to get similar or even better quality of reconstructed images. Chapter four introduces a novel variational model for inverse consistent deformable image registration. The proposed model deforms both source and target images simultaneously, and aligns the deformed images in the way that the forward and backward transformations are inverse consistent. To avoid the direct computation of the inverse transformation fields, our model estimates two more vector fields by minimizing their invertibility error using the deformation fields. Moreover, to improve the robustness of the model to the choice of parameters, the dissimilarity measure in the energy functional is derived using the likelihood estimation. The experimental results on clinical data indicate the efficiency of the proposed method with improved robustness, accuracy and inverse consistency. These methods are aimed to benefit the practical usage of medical imaging techniques.
General Note: In the series University of Florida Digital Collections.
General Note: Includes vita.
Bibliography: Includes bibliographical references.
Source of Description: Description based on online resource; title from PDF title page.
Source of Description: This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
Statement of Responsibility: by XIAOJING YE.
Thesis: Thesis (Ph.D.)--University of Florida, 2011.
Local: Adviser: Chen, Yunmei.

Record Information

Source Institution: UFRGP
Rights Management: Applicable rights reserved.
Classification: lcc - LD1780 2011
System ID: UFE0042919:00001

Permanent Link: http://ufdc.ufl.edu/UFE0042919/00001

Material Information

Title: Mathematical Models and Fast Numerical Algorithms in Medical Image Analysis
Physical Description: 1 online resource (121 p.)
Language: english
Creator: YE,XIAOJING
Publisher: University of Florida
Place of Publication: Gainesville, Fla.
Publication Date: 2011

Subjects

Subjects / Keywords: CONSISTENT -- CONVERGENCE -- MRI -- OPTIMIZATION -- RECONSTRUCTION -- REGISRATION -- STATISTICS
Mathematics -- Dissertations, Academic -- UF
Genre: Mathematics thesis, Ph.D.
bibliography   ( marcgt )
theses   ( marcgt )
government publication (state, provincial, terriorial, dependent)   ( marcgt )
born-digital   ( sobekcm )
Electronic Thesis or Dissertation

Notes

Abstract: A series of projects are presented to provide advanced mathematical models and numerical algorithms that improve the accuracy, robustness and efficiency of the compressed sensing reconstruction technique in magnetic resonance imaging and inverse consistent image regisration. Chapter one introduces a novel variational model that enforces the sparsity of the underlying image in terms of its spatial finite differences and representation with respect to a dictionary. The dictionary is trained using prior information to improve accuracy in reconstruction. In the meantime the proposed model enforces the consistency of the underlying image with acquired data by using the maximum likelihood estimator of the reconstruction error in partial k-space to improve the robustness to parameter selection. Moreover, a simple and fast numerical scheme is provided to solve this model. In chapters two and three we develop fast numerical algorithms for solving total variation and L1 (TVL1) based image reconstruction with application in partially parallel MR imaging. Our algorithms use variable splitting method to reduce computational cost. Moreover, the Barzilai-Borwein step size selection method is adopted in our algorithms for much faster convergence. Theoretical and experimental results on clinical partially parallel imaging data demonstrate that the proposed algorithm requires much fewer iterations and/or less computational cost than recently developed operator splitting and Bregman operator splitting methods, which can deal with a general sensing matrix in reconstruction framework, to get similar or even better quality of reconstructed images. Chapter four introduces a novel variational model for inverse consistent deformable image registration. The proposed model deforms both source and target images simultaneously, and aligns the deformed images in the way that the forward and backward transformations are inverse consistent. To avoid the direct computation of the inverse transformation fields, our model estimates two more vector fields by minimizing their invertibility error using the deformation fields. Moreover, to improve the robustness of the model to the choice of parameters, the dissimilarity measure in the energy functional is derived using the likelihood estimation. The experimental results on clinical data indicate the efficiency of the proposed method with improved robustness, accuracy and inverse consistency. These methods are aimed to benefit the practical usage of medical imaging techniques.
General Note: In the series University of Florida Digital Collections.
General Note: Includes vita.
Bibliography: Includes bibliographical references.
Source of Description: Description based on online resource; title from PDF title page.
Source of Description: This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
Statement of Responsibility: by XIAOJING YE.
Thesis: Thesis (Ph.D.)--University of Florida, 2011.
Local: Adviser: Chen, Yunmei.

Record Information

Source Institution: UFRGP
Rights Management: Applicable rights reserved.
Classification: lcc - LD1780 2011
System ID: UFE0042919:00001


This item has the following downloads:


Full Text

PAGE 1

MATHEMATICALMODELINGANDFASTNUMERICALALGORITHMSINMEDICALIMAGINGAPPLICATIONSByXIAOJINGYEADISSERTATIONPRESENTEDTOTHEGRADUATESCHOOLOFTHEUNIVERSITYOFFLORIDAINPARTIALFULFILLMENTOFTHEREQUIREMENTSFORTHEDEGREEOFDOCTOROFPHILOSOPHYUNIVERSITYOFFLORIDA2011

PAGE 2

c2011XiaojingYe 2

PAGE 3

Tomyparents 3

PAGE 4

ACKNOWLEDGMENTS Foremost,mygreatestthankfulnessgoestomyadvisor,YunmeiChen,forherwiseguidanceandsupportthroughoutmytimeasagraduateresearcher.Iamverygratefulforherencouragementandeffortstomyacademiccareer.ManythankstoFengHuangforthefruitfuldiscussionontheinterestingproblemsinmagneticresonanceimaging.Severalprojectsrelatedtopartiallyparallelimagingreconstructionwereoriginallymotivatedbytheseclinicalapplications.IwouldliketothankLeiZhangforvaluablediscussionaboutmyacademiccareeraswellashisconsistentsupporttomystudyandlifetheseyears.ThankstoWilliamHagerandhisstudentsDungPhanandHongchaoZhangforthehelpfuldiscussiononoptimizationtechniqueswhichIhavebeenusingforseveralofmyprojects.ManythankstomyothercommitteemembersJayGopalakrishnan,ScottMcCullough,JonathanLiandMuraliRao,whohavebeenverysupportivetomygraduateresearch.IwouldliketothankHaniDoss,LiqingYan,andMarkYangforsupervisingmymasterthesisattheDepartmentofStatistics.Althoughthisissomehowanareaquitedifferentfrommyown,Ilearnedalotfromthem. 4

PAGE 5

TABLEOFCONTENTS page ACKNOWLEDGMENTS .................................. 4 LISTOFTABLES ...................................... 7 LISTOFFIGURES ..................................... 8 ABSTRACT ......................................... 10 CHAPTER 1ANovelMethodandFastAlgorithmforMRImageReconstructionviaLearntDictionaries ...................................... 12 1.1BackgroundsinCompressiveMagneticResonanceImaging ........ 12 1.1.1TrainedDictionariesasSparsifyingTransforms ........... 14 1.1.2LikelihoodEstimateasDataFidelityMeasure ............ 16 1.1.3FastNumericalAlgorithmsforSolvingCS-MRIModels ....... 18 1.1.4Organization .............................. 19 1.2ProposedModel ................................ 19 1.2.1SparseRepresentationUsingTrainedDictionary .......... 19 1.2.2LikelihoodEstimatefortheDataConsistency ............ 21 1.2.3VariationalModelforMRImageReconstructionfromUndersampledData ................................... 22 1.3Algorithm .................................... 23 1.3.1AFastAlgorithmforSolvingtheProposedModel .......... 23 1.3.2NumericalAlgorithmandConvergenceAnalysis ........... 26 1.4ExperimentalResults ............................. 27 1.4.1ImprovementonAccuracybyUsingDictionaries .......... 27 1.4.1.1ResultsofPhantomReconstruction ............ 28 1.4.1.2ResultsofBrainImageReconstruction .......... 29 1.4.1.3ResultsofChestImageReconstruction .......... 31 1.4.2ImprovementonRobustness(toParameterSelection)andEfciency 33 1.4.3RobustnessofDictionaryTrainingandReconstruction ....... 33 2ComputationalAccelerationforMRImageReconstructioninPartiallyParallelImaging ........................................ 37 2.1BackgroundsinTotalVariationBasedImageReconstruction ....... 37 2.1.1PartiallyParallelMRImaging ..................... 38 2.1.2PreviousWork ............................. 40 2.1.3Organization .............................. 44 2.2ProposedAlgorithm .............................. 44 2.3Method ..................................... 51 2.3.1DataAcquisition ............................ 51 5

PAGE 6

2.3.2TestEnvironment ............................ 53 2.4ComparisonAlgorithmsandResults ..................... 54 2.4.1ComparisonwithBOS ......................... 54 2.4.2ComparisonwithOS .......................... 56 2.5ConcludingRemarks .............................. 58 3AFastAlgorithmforTVImageReconstructionwithApplicationtoPartiallyParallelMRImaging ................................. 63 3.1BackgroundsinOptimizationMethodsforNonsmoothBasedImageReconstruction 63 3.2RelatedWork .................................. 65 3.3ConvergenceAnalysis ............................. 69 3.4AlgorithmsfortheTVandLSSubproblems ................. 76 3.5PartiallyParallelImaging(PPI) ........................ 78 3.6NumericalExperiments ............................ 81 3.6.1DataAcquisitionandExperimentalSetup .............. 81 3.6.2ComparisonAlgorithm ......................... 82 3.6.3ExperimentalResults .......................... 83 3.7ConcludingRemarks .............................. 85 4InverseConsistentDeformableImageRegistration ................ 88 4.1BackgoundofConsistentImageRegistration ................ 88 4.2ProposedMethod ............................... 93 4.2.1MotivationandIdeasofProposedMethod .............. 93 4.2.2AlternativeFormulationof( 4 )UsingDeformationFields ..... 95 4.2.3MLEbasedderivationfordis(S,T~) .............. 95 4.2.4Proposedmodel ............................ 97 4.3ExistenceofSolutions ............................. 99 4.4NumericalScheme ............................... 102 4.5ExperimentalResults ............................. 104 REFERENCES ....................................... 115 BIOGRAPHICALSKETCH ................................ 121 6

PAGE 7

LISTOFTABLES Table page 1-1Comparisonofresultsofphantomreconstructionsusingnonlinearconjugategradient(CG)formodel( 1 )withHaarwavelet,algorithmrecMRIwithHaarwavelet,andrecMRIformodel( 1 )withovercompleteDCTdictionary. .... 34 1-2Experimentalresultsof10runsofdictionarytrainingandbrainimagereconstructionasin 1.4.1.2 ..................................... 34 2-1Testsnumber,datainformationandparameters. ................. 53 2-2ResultsofBOSandTVL1recondata1. ...................... 55 3-1TestedDataandModelParameters ........................ 83 4-1NumberofiterationsusedforconvergenceandthenalCCobtainedbyproposedmodelwithupdated/xed. ............................. 108 7

PAGE 8

LISTOFFIGURES Figure page 1-1DictionarytrainedbyK-SVDalgorithm. ...................... 16 1-2Comparetheaccuracyofsparserepresentationsbywaveletandtraineddictionary. ............................................. 17 1-3Samplingmasksusedwithsamplingratios(SR). ................. 28 1-4Reconstructedphantomimagefromsimulatedpartialk-spacedata. ...... 30 1-5ReconstructionofbrainMRImage. ......................... 35 1-6ReconstructionsofchestMRimageusingmodel( 1 )and( 1 ). ....... 36 2-1Thek-spacemasksusedforthethreedatasetstestedintheexperiments. .. 52 2-2ComparisonofBOSandTVL1recondata1. .................... 59 2-3TestingdatausedforthecomparisonofOSandTVL1recfordata2anddata3. 60 2-4Reconstructionsofdata2anddata3byOSandTVL1rec. ............ 61 2-5ComparisonsofOSandTVL1recondata2(bluesolidlines)anddata3(blackdashedlines). ..................................... 62 3-1PPImask ....................................... 79 3-2Sensitivitymapsforan8-channelcoil ....................... 80 3-3AnillustrationofaheadPPIsystemwitheightreceivercoils. .......... 82 3-4ReconstructedimagesbyOSandAMforSB512. ................ 85 3-5ReconstructionsbyOSandAMforAB512. .................... 86 3-6ComparisonofOSandAMondatasetsAB256andSB512. ........... 87 4-1Inverseconsistentregistrationresultbyproposedmodel( 4 ). ........ 109 4-2Residueimageobtainedbyproposedmodel( 4 ). ............... 110 4-3Deformationeldsobtainedbyproposedmodel( 4 )inthezoomed-inareaappliedonregulargridwithhalfoforiginalresolutionofimages. ......... 110 4-4Autore-contouringresultusingthedeformationelduobtainedbyproposedmodel( 4 ). ..................................... 110 4-5Registrationresultofproposedmodel( 4 )appliedto3DchestCTimage.Thisgureshowsthez=33sliceataxialdirection. ............... 111 8

PAGE 9

4-6Registrationresultofproposedmodel( 4 )appliedto3DchestCTimage.Thisgureshowsthex=25sliceatsagittaldirection. ............. 112 4-7Registrationresultofproposedmodel( 4 )appliedto3DchestCTimage.Thisgureshowsthey=48sliceatcoronarydirection. ............. 113 4-8CCineachiterationobtainedbyfull-waymodel( 4 )andproposedmodel( 4 ). ........................................ 114 4-9Meanofinverseconsistenterrors(ICE)ofthenaldeformationeldsobtainedbyusingfull-waymodel( 4 )andproposedmodel( 4 ). ........... 114 9

PAGE 10

AbstractofDissertationPresentedtotheGraduateSchooloftheUniversityofFloridainPartialFulllmentoftheRequirementsfortheDegreeofDoctorofPhilosophyMATHEMATICALMODELINGANDFASTNUMERICALALGORITHMSINMEDICALIMAGINGAPPLICATIONSByXiaojingYeMay2011Chair:YunmeiChenMajor:MathematicsAseriesofprojectsarepresentedtoprovideadvancedmathematicalmodelsandnumericalalgorithmsthatimprovetheaccuracy,robustnessandefciencyofthecompressedsensingreconstructiontechniqueinmagneticresonanceimagingandinverseconsistentimageregisration.Chapteroneintroducesanovelvariationalmodelthatenforcesthesparsityoftheunderlyingimageintermsofitsspatialnitedifferencesandrepresentationwithrespecttoadictionary.Thedictionaryistrainedusingpriorinformationtoimproveaccuracyinreconstruction.Inthemeantimetheproposedmodelenforcestheconsistencyoftheunderlyingimagewithacquireddatabyusingthemaximumlikelihoodestimatorofthereconstructionerrorinpartialk-spacetoimprovetherobustnesstoparameterselection.Moreover,asimpleandfastnumericalschemeisprovidedtosolvethismodel.Inchapterstwoandthreewedevelopfastnumericalalgorithmsforsolvingtotalvariationand`1(TVL1)basedimagereconstructionwithapplicationinpartiallyparallelMRimaging.Ouralgorithmsusevariablesplittingmethodtoreducecomputationalcost.Moreover,theBarzilai-Borweinstepsizeselectionmethodisadoptedinouralgorithmsformuchfasterconvergence.Theoreticalandexperimentalresultsonclinicalpartiallyparallelimagingdatademonstratethattheproposedalgorithmrequiresmuchfeweriterationsand/orlesscomputationalcostthanrecentlydevelopedoperatorsplittingand 10

PAGE 11

Bregmanoperatorsplittingmethods,whichcandealwithageneralsensingmatrixinreconstructionframework,togetsimilarorevenbetterqualityofreconstructedimages.Chapterfourintroducesanovelvariationalmodelforinverseconsistentdeformableimageregistration.Theproposedmodeldeformsbothsourceandtargetimagessimultaneously,andalignsthedeformedimagesinthewaythattheforwardandbackwardtransformationsareinverseconsistent.Toavoidthedirectcomputationoftheinversetransformationelds,ourmodelestimatestwomorevectoreldsbyminimizingtheirinvertibilityerrorusingthedeformationelds.Moreover,toimprovetherobustnessofthemodeltothechoiceofparameters,thedissimilaritymeasureintheenergyfunctionalisderivedusingthelikelihoodestimation.Theexperimentalresultsonclinicaldataindicatetheefciencyoftheproposedmethodwithimprovedrobustness,accuracyandinverseconsistency.Thesemethodsareaimedtobenetthepracticalusageofmedicalimagingtechniques. 11

PAGE 12

CHAPTER1ANOVELMETHODANDFASTALGORITHMFORMRIMAGERECONSTRUCTIONVIALEARNTDICTIONARIESOutline Theaimofthisworkistoimprovetheaccuracy,robustnessandefciencyofthecompressedsensingreconstructiontechniqueinmagneticresonanceimaging.Weproposeanovelvariationalmodelthatenforcesthesparsityoftheunderlyingimageintermsofitsspatialnitedifferencesandrepresentationwithrespecttoadictionary.Thedictionaryistrainedusingpriorinformationtoimproveaccuracyinreconstruction.Inthemeantimetheproposedmodelenforcestheconsistencyoftheunderlyingimagewithacquireddatabyusingthemaximumlikelihoodestimatorofthereconstructionerrorinpartialk-spacetoimprovetherobustnesstoparameterselection.Moreover,asimpleandfastnumericalschemeisprovidedtosolvethismodel.Theexperimentalresultsonbothsyntheticandinvivodataindicatetheimprovementoftheproposedmodelinpreservationofnestructures,exibilityofparameterdecision,andreductionofcomputationalcost. 1.1BackgroundsinCompressiveMagneticResonanceImagingMagneticresonance(MR)imagingisatechniquethatallowsvisualizationofstructuresandfunctionsofabodybynon-invasiveandnon-ionizingmeans.Itprovidesbettercontrastbetweenthedifferentsofttissuesthanmostothermodalities.However,MRimagingtakesmuchlongeracquisitiontimethanmanyotherimagingmodalities,whichlimitsitsapplications.Toreduceacquisitiontime,themostcommonandfeasibleapproachisbyacquiringonlypartialk-spacedata,followedbyadequatereconstructiontechniquestoobtainimageswithwell-preservedquality.Theideaofreconstructingimagesfrompartialdatacoincideswiththecompressedsensing(CS),atechniqueusedinsignal/imageprocessing.CScanaccuratelyrecoverasignal/imageusingdatawithsignicantlylessmeasurementsthanregular,providedthesparsityoftheunderlyingsignal/imageandasophisticatedreconstructionprocedure. 12

PAGE 13

Recently,theapplicationofthistechniqueinmedicalimaginghasbecomeahotresearchtopic,andshownpromisingresults[ 16 18 21 22 27 42 46 50 68 ].Inparticular,theredundancyoftheMRdataacquiredinthefrequencydomain,i.e.thek-space,andimplicitsparsityinMRimageshavemotivatedmanyresearcherstostudytheapplicationofCStofastMRimaging(CS-MRI).CS-MRIhastheadvantageofproducinghighqualityreconstructionofMRimagesfrompartialFourierdata.RecentstudyhasshownthatthekeytothesuccessofCS-MRIisacombinationofthesparsityoftheunderlyingimageunderanappropriatedomain,thek-spacetrajectorythatprovidesincoherentundersamplingartifacts,andanadequatenonlinearreconstructionmethodthatenforcesboththesparsityanddataconsistencyoftheunderlyingimage[ 26 38 50 ].AgreatprogressofresearchesonCS-MRIhasbeenmade.However,forclinicalapplications,radiologistsoftendemandimprovementsonaccuracy,robustness,andefciencyofthecurrentCS-MRIalgorithms.Thedesiredimprovementsincludetheabilityofremovingartifactswhilepreservingimportantdiagnosticinformation(inparticular,sharpedgesandnestructures),therobustnesstothechoiceofparameters,andthespeedofreconstructions.Inthispaper,weproposeanovelvariationalmodelandafastnumericalalgorithmforMRimagereconstructionwithhighlyundersampleddata,whichtacklesthethreeproblemsmentionedaboveasfollows. AccuracyTheproposedmodelenforcesthesparsityoftheunderlyingimageintermsofitsspatialnitedifferencesandrepresentationbyadictionarytrainedusingpriorinformation.Thus,improvementonaccuracyofreconstructioncanbeachieved. RobustnessTheproposedmodelenhancesthedataconsistencybytheapproachofmaximumlikelihoodestimationforthediscrepancybetweenthereconstructionandacquireddataink-space.Thisleadstoanautomaticallyoptimizedweightingparameterwhichmakestheparameterselectionmoreexible. Efciency 13

PAGE 14

Tomaketheproposedmodelclinicallyapplicable,wealsoprovideasimpleandfastnumericalalgorithmtosolvethemodel.Themaincomputationsinvolveonlyshrinkage,matrix-vectormultiplicationandfastFouriertransform(FFT).Thebackgroundandbriefintroductionofourcontributionstotheseissuesareprovidedinthefollowingthreesubsections. 1.1.1TrainedDictionariesasSparsifyingTransformsSincesparsityisthekeytothesuccessofCSandconsequentreconstructions,manyresearchesexploitedthetransformsunderwhichimageshavetheirsparserepresentations[ 50 ].ThetheoryofCSindicatesthatoncesuchtransformswerefound,animagecanbeaccuratelyrecoveredusingasetofrandommeasurementswithcardinalitymuchlessthantheoriginalresolutionoftheimage[ 18 26 ].Inrecentyears,nitedifferenceoperatorandwavelettransformshavebeenwidelyusedassuchsparsifyingtransformsforMRimages[ 34 50 ].In[ 34 ],theauthorsproposedatotalvariation(TV)basedmodeltoreconstructMRimagesfrompartiallyacquiredk-spacedata.Theirmodelworkswellforpiecewiseconstantorveryhomogeneousimages[ 59 ].Forimageswithinhomogeneousintensity,TVbasedmodelsmaynotworkwellwhentheundersamplingrateishigh.In[ 50 ],Lustigetal.proposedamodelthatminimizestheBesovtogetherwithTVnormsoftheunderlyingimage,subjectedtoadataconsistencyconstraint: minuTV(u)+k>uk1,s.t.kFpu)]TJ /F4 11.955 Tf 11.96 0 Td[(fpk2<,(1)wherekk1isthe`1norm,TV(u),kDuk1isthe(anisotropic)TVsemi-normofu,istheHaarwavelettransform,thesuperscript>denotes(conjugate)transposeofmatrices.Intheconstraint,FpdenotestheundersampledFourieroperatorcorrespondingtothecustomizedk-spacesamplingpattern,fpisthepartiallyacquiredk-spacedata,andisanestimateofacquisitionerror.Asprovedin[ 18 26 ],minimizing`1normsubjectedtodataconsistencyyieldssparsesolutions.Therefore,model( 1 )infactleadstoareconstructedimagethathassparsegradientandwavelettransform 14

PAGE 15

coefcients.Simulationsin[ 50 ]showedverypromisingresultsusingmodel( 1 ).However,imagequalitydegradingandthelossofdiagnosticinformationmayhappeninreconstructionsusingTVandwavelettransformsastheymayeliminatesomenestructuresand/orusefullocalinformationintherecoveredimages.Asanalternate,weproposetousethedictionariestrainedusingpriorinformationassparsifyingtransformstotacklethisproblem.Arecentboostofthestudyondictionarydesignshowsthegreatpotentialofusingdictionariesinsignal/imageprocessing.Dictionaryisusuallyformedasasetofovercompletebasesanditselements/atomshavemuchsmallersizesthantheimagesize.Onthecontrary,wavelettransformhasasetofcompletebaseswithelementsofthesamesizeasimageitself,andthereforecanbetreatedasaspecialcaseofdictionary.Furthermore,adictionarycanbeproperlytrainedsuchthatitsprototypesignal-atomsaremoreadequatetosparselyrepresentobjectivesignalsthanwavelet.Anumberofresearcheshaveshownthebenetsofusingdictionariesforsparserepresentation;see,e.g.[ 46 62 ].Inthiswork,wetrainadictionarybyapplyingK-SVDalgorithmtoadatabaseconsistingofpatchesextractedfromimagesacquiredfromthesamesequencebutperhapsdifferentsubjects.ThenthetraineddictionaryA,asshowninFigure 1-1 ,isusedtoreconstructotherMRimagesunderthesameacquisitionsequencewithsimilarstructures.Thedatabaseusedfortrainingconsistsof4096patchesextractedfromfourMRbrainimages(butexcludestheimagetobereconstructed).Eachblockrepresentsanatomofsize88.Atomsaresortedbyascendingthevariancesofintensities.Detailsonthetrainingandreconstructionprocessesareprovidedinthefollowingsections.ComparisonoftheaccuracyofsparserepresentationsusingwavelettransformandthetraineddictionaryAisshowninFigure 1-2 .Resultsbyusingdictionaryhashasbetterpreservededgesandnestructuresbecausedictionaryabsorbedpriorknowledgebylearningfeaturesofthistypeofimagesduringthetrainingprocess[ 2 ]. 15

PAGE 16

Moreover,thedictionarytrainingandrepresentationprocessesarestableandrobustasshowninourexperimentalresultsinsection 1.4 Figure1-1. DictionarytrainedbyK-SVDalgorithm. InbothsparserepresentationsbywaveletandtraineddictionaryinFigure 1-2 ,imagesarerepresentedbypickingupthelargest12.5%transformcoefcients.Bottomimagesarecorrespondingzoomed-insquareareashownonthetopleftimage.Leftcolumn:reference.Middlecolumn:representationbywavelet,withRMSE5.74%andSNR22.3.Rightcolumn:representationbytraineddictionaryAshowninFigure 1-1 ,withRMSE4.36%andSNR24.7. 1.1.2LikelihoodEstimateasDataFidelityMeasureToimprovetherobustnesstotheparameterselection,weusethelikelihoodestimationofthereconstructionerrorasthedatadelitymeasure.ThereconstructionerroristhedifferencebetweenthepartiallyacquireddataandtheFouriertransformofthereconstructionatsampledk-spacelocations,i.e.Fpu)]TJ /F4 11.955 Tf 12.65 0 Td[(fpin( 1 ).InpreviouslyproposedCS-MRIalgorithms,leastsquares,i.e.thesumofsquareddifference(SSD)kFpu)]TJ /F4 11.955 Tf 13.06 0 Td[(fpk22,isacommonlyuseddataconsistencymeasure.Forinstance, 16

PAGE 17

A B C D E F Figure1-2. Comparetheaccuracyofsparserepresentationsbywaveletandtraineddictionary. theunconstrainedversionofmodel( 1 )solvesforthereconstructionby minuTV(u)+k>uk1+ 2kFpu)]TJ /F4 11.955 Tf 11.96 0 Td[(fpk22,(1)wheretheparameteriscrucialtothereconstructionresults:animproperlylargeweightforthedatadelitytermresultsinseriousresidualartifacts,whereasanimproperlysmallweightresultsindamagededgesand/ornestructures.Inthiswork,wetacklethisproblembytreatingthereconstructionerrorsatallpixelsassamplesindependentlydrawnfromaGaussiandistributionwithmeanzeroandvariance2tobeoptimized.Bymaximumlikelihoodestimate(MLE)approach,theweightonkFpu)]TJ /F4 11.955 Tf 12.69 0 Td[(fpk22=2becomes=2ratherthanaprescribed.Sinceistobeoptimized,itisupdatedduringiterations(infact,itisthestandarddeviationofthereconstruction 17

PAGE 18

error,see( 1 )below).Whenthereconstructionerrorreduces,theweight=2onkFpu)]TJ /F4 11.955 Tf 10.91 0 Td[(fpk22increases,andhencetheaccuracyisimproved.Thisautomaticallyoptimizedweightingfeaturemakesthechoiceofmuchmoreexible. 1.1.3FastNumericalAlgorithmsforSolvingCS-MRIModelsDespitethatdictionariesaremoreadequateinsignal/imagereconstructions,thecomputationalcostishigherthanthatusingwavelettransformduetotheredundancyofdictionariesandoverlappingofpatchestoberepresented.Also,thenon-differentiabilityofTVand`1normsbringsdifcultiestofastsolutionsofCS-MRImodels.TherehavebeenmanynumericalalgorithmsdevelopedtosolveTVand`1regularizedminimizationproblems,morerecentdevelopmentscanbefoundin[ 19 34 64 67 78 79 ]and[ 11 14 36 41 52 66 ]andreferencestherein.Ourapproachinthisworkiscloselyrelatedtothealgorithmdevelopedin[ 67 ],inwhichYangetal.introducedasimpleandfastmethodtosolvemodel( 1 )withisotropicTVnormofudenedby TV(u),NXi=1kDiuk2.(1)In( 1 )u2RNisthevectorformedbystackingallcolumnsoftheimagevertically,Nisthetotalnumberofpixelsintheimage,andDi2R2Nrepresentsthegradientoperatoratthei-thelementinthevectorofu.Toovercomethenon-differentiabilityofTVand`1norms,theyintroducedauxiliaryvariablesandusedaclassicalquadraticpenaltymethodwhichyieldsanalternatingminimizationscheme.BydiagonalizingthegradientoperatorusingFouriertransform,theymadethemaincomputationofthealgorithminvolvingonlysoftshrinkageandfastFouriertransform.However,theiralgorithmcannotbedirectlyappliedtothemodelsusingdictionaryAsinceitrequirestheorthogonalityofin( 1 ).Thereforethedevelopmentofefcientalgorithmsinvolvingtheuseofadictionaryisstillaremainingproblem.Inthispaper,weshowthatasimpletrickonselectingpatchesfromtheunderlyingimagecanbeusedtoovercomethisdifculty. 18

PAGE 19

Basedonthemethodin[ 67 ]andthistrick,weprovideasimpleandfastnumericalalgorithmthatcanbeappliedtoreconstructionmodelsinvolvingdictionaries. 1.1.4OrganizationTherestofthispaperisorganizedasfollows.Adetaileddescriptionoftheproposedmodelisgivenin 4.2 .Insection 4.4 ,afastalgorithmtosolvetheproposedmodelanditsderivationareprovided.Experimentalresultsonsyntheticandinvivodataarepresentedinsection 1.4 .Thelastsectionconcludesthepaper. 1.2ProposedModelBeforegoingintodetailsoftheproposedmodel,weaddressthenotationsusedthroughouttherestofthepaper.Firstofall,allvectorsinthispaperarecolumnvectors.Letu2RNbetheunderlyingreconstructionasin( 1 ),andFbethediscreteFouriertransform,whichcanbetreatedasanNNunitarymatrix.LetP2RpNdenotethebinarymatrixthatselectscertainrowsofFaccordingtothek-spacesamplingpattern.ThenFp,PFistheundersampledFouriertransform.Letfp2Cpbethepartiallyacquiredk-spacedataandusekktodenoteEuclideannormkk2ofvectorshenceforth.Thenotation(;)representsamatrixformedbystackingitsargumentsvertically.Inthispaper,allpatcheshavesizep np nandareoftentobetreatedasn-vectorsunlessotherwisenoted(inourexperimentsn=64).ThedictionaryA2RnKconsistsofKn-vectorsasatoms.BinarymatrixRj2RnNextractsthej-thpatchofu,andformsthepatchRjuasann-vector.AllpatchesfRjugJj=1covertheentireimageuandmayoverlap. 1.2.1SparseRepresentationUsingTrainedDictionaryToimprovetheaccuracyofreconstruction,especiallytheabilityofpreservingdiagnosticinformation,wehereproposetouseatraineddictionaryinsteadofwaveletasthesparsifyingtransforminMRimagereconstructions.Wechosetherecentlyproposeddictionarydesignmethod,termedasK-SVDalgorithm,toperformthedictionary 19

PAGE 20

training.K-SVDisaniterativemethodthatalternatesbetweensparsecodingoftheexamplesbasedonthecurrentdictionaryandaprocessofupdatingthedictionaryatomstobettertthegivendatabase.Theoutputisatraineddictionarythatcanrepresentallsignalsinthedatabaseunderstrictsparsityconstraintsanderrortolerance.Interestedreadersarereferredto[ 2 30 ]fordetails.OurprocedureofformingadatabaseandapplyingK-SVDalgorithmtotrainanadaptivedictionaryforMRimagereconstructionisdepictedasfollows. 1. CollectanumberofMRimagesacquiredusingthesamesequenceasthatfortheimagetobereconstructed,butfromdifferentsubjects.Thetrainingimagesandtheimagetobereconstructedarepreferredtobethesamebodypartstogetabettersparserepresentation.Usingthesameacquisitionsequenceensuresthattheyhavesimilarcontrast. 2. Decomposethetrainingimagestop np npatches,anddiscardthosepatcheswithconstantintensity.Thenrandomlychoose8Kpatchesfromtheremainingpatches,whereKisthenumberofatomsinthedictionarytobetrained. 3. TrainadictionaryAbyapplyingK-SVDalgorithmtothat8KpatcheswiththeovercompleteDCTmatrixastheinitial.TheresultingtraineddictionaryhasKelements,i.e.A2RnK.Inourexperiments,wesetnto64andKto512=256forbrain/chestMRdata.ThedictionarytrainedforbrainimagereconstructionisillustratedinFigure 1-1 .ThedictionaryAweobtainfromthistrainingprocedurecanadequatelyrepresentanypatchesofbrainMRimages(e.g.fromdifferentsubjectsorthesamesubjectindifferentperiods)acquiredunderthesamesequence.Inparticular,eachpatchRjuofucanbesparselyrepresentedbyA.Namely,thereexistrepresentationcoefcientsj2RKsuchthatkjk0
PAGE 21

reconstruction.Thatis,weenforcethefollowingintoourmodel minJXj=1kjk1+ 2kAj)]TJ /F4 11.955 Tf 11.96 0 Td[(Rjuk2,(1)where=(1;;J)2RKJ.Thisisinfacttherelaxedformofsparse-landproblemwith`0-normsubstitutedby`1-norm.Thereasonwhyweuse`1insteadof`0isthatminimizingthenon-convex`0isgenerallyaNP-hardproblemandhenceisnottractableinpractice.Moreover,ithasbeenprovedthatminimizationproblemswith`1and`0sharethesamesolutionundercertainconditions[ 18 26 ].NotethatifJ=1,R1=I(theidentitymatrix)andA=(thewavelettransform),then( 1 )reducestok>uk1asin( 1 )whenthedifferenceinthequadraticisexactly0.Namely,waveletisaspecialcaseofdictionary. 1.2.2LikelihoodEstimatefortheDataConsistencyOnedifcultyofapplyingtheunconstrainedenergyminimizationproblem( 1 )forMRimagereconstructionisindeterminingtheweightingparameterthatbalancesthedataconsistencyandimagesparsity.Thereconstructionresultsaresensitivetothechoiceofthisparameter.Totacklethisproblem,wederivethedataconsistencymeasure,theso-calleddatadelity,frommaximumlikelihoodestimate(MLE)approach.Let!=(!1,,!p)>2Cpbethereconstructionerrorink-space,whichisthedifferencebetweentheFouriertransformofthereconstructionuandtheacquireddatafpatthesampledk-spacelocations:fp=Fpu+!.Considerf!lgpl=1asindependentrandomsamplesdrawnfromanormaldistributionofmeanzeroandvariance2tobedetermined.Therefore,thejointprobabilitydensityfunction(pdf)off!lgpl=1,whichisalsothelikelihoodof,becomesL(j!)=pYl=11 p 2e)]TJ /F15 7.97 Tf 6.59 0 Td[(!2l=22=)]TJ /F3 11.955 Tf 5.48 -9.69 Td[(22)]TJ /F9 7.97 Tf 6.58 0 Td[(p=2ek!k2=22. 21

PAGE 22

Thus,thenegativelog-likelihoodis )]TJ /F3 11.955 Tf 11.96 0 Td[(logL(j!)=k!k2=22+plogp 2.(1)Substituting!byFpu)]TJ /F4 11.955 Tf 12.15 0 Td[(fp,andomittingtheconstantplogp 2,weobtainaMLEbasedconsistencyestimationwiththepartiallyacquireddata: F(u,,fp)=kFpu)]TJ /F4 11.955 Tf 11.96 0 Td[(fpk2=22+plog.(1)Thisisageneralizationoftheleastsquareestimation,whichisjustthecasewhere1.Wewilluse( 1 )asdatadelityterminourenergyfunctional. 1.2.3VariationalModelforMRImageReconstructionfromUndersampledDataNowwearereadytopresentourmodel.WeproposetouseTVandsparserepresentationbytraineddictionaryasregularizationandMLE( 1 )asdataconsistencymeasure.Ourmodelisformulatedasanunconstrainedminimizationproblem minu,,TV(u)+JXj=1kjk1+ 2kAj)]TJ /F4 11.955 Tf 11.95 0 Td[(Rjuk2+F(u,,fp), (1) whereTV(u)istheTVnormofudenedasin( 1 ),thesummationoveristheregularizationofuusingthesparsityunderrepresentationbydictionaryA,andF(u,,fp)isMLEdataconsistencymeasure( 1 ).ByusingMLEbasedapproach,isalsooptimizedalongwithu.In( 1 )theweightonkFpu)]TJ /F4 11.955 Tf 12.28 0 Td[(fpk2versusthesparsityoftheunderlyingimageis=2ratherthanonly.IntheEuler-Lagrange(EL)equationsassociatedwiththeproposedenergyfunctionbelow,onecanseethatisthestandarddeviationofthereconstructionerror!.Hence,whentheconstructionerror!decreases,theweight=2onminimizing`2normof!increasesautomatically,whichmakesthechoiceoftheinitialweightingparametermoreexible.ThisexibilitydramaticallyreducesthedifcultyofparameterdecisionandimprovestheapplicabilityofCS.Moreover,ourexperimentalresultsbelowshowthatthisautomaticallyupdatedweightingmakesfasterconvergence,andbetteraccuracyinreconstruction. 22

PAGE 23

1.3AlgorithmThereminimizationproblem( 1 )iscloselyrelatedtothewell-knownTVand`1basedsignal/imagereconstructionproblems.Sincethenon-differentiabilityofTVand`1termsbringadifcultyincomputations,therehavebeenanumberofnumericalalgorithmsdevelopedtoefcientlysolvethistypeofproblems.Thealgorithmprovidedinthissectionisinspiredbytheworkin[ 64 67 ],whichusesthevariablesplittingandclassicalquadraticpenaltytechniqueinoptimizationtomakethecomputationfastandstable. 1.3.1AFastAlgorithmforSolvingtheProposedModelWerstintroducetwoauxiliaryvariablesw=(w>1;w>2;;w>N)2RN2and=(1;2;;J)2RKJwherewi2R2andj2RKforalli=1,,Nandj=1,,J.Thenweconsidertheminimizationproblemequivalentto( 1 ): minu,w,,,NXi=1kwik+JXj=1kjk1+ 2kAj)]TJ /F4 11.955 Tf 11.96 0 Td[(Rjuk2+F(u,,fp)s.t.wi=Diu,j=j,8i=1,,N,j=1,,J. (1) Relaxingtheequalityconstraintandpenalizingtheirviolationsbyquadraticfunctions,weobtainanunconstrainedversionof( 1 ): minu,w,,,NXi=1(wi,Diu)+ (,)+JXj=1 2kAj)]TJ /F4 11.955 Tf 11.95 0 Td[(Rjuk2+F(u,,fp)(1)wherefunctionsand aredenedas(s,t)=ksk+ 2ks)]TJ /F4 11.955 Tf 11.95 0 Td[(tk2,s,t2R2and (s,t)=ksk1+ 2ks)]TJ /F4 11.955 Tf 11.95 0 Td[(tk2,s,t2RKJforgiven,>0.With,graduallyincreasing,solving( 1 )leadtoapproximationstothesolutionof( 1 ). 23

PAGE 24

Theminimization( 1 )canbecarriedoutinamuchfasterandmorestablemannerthan( 1 ):rst,forxeduand,theminimizationwithrespecttowandcanbecarriedoutinparallel: wi=S2(Diu),i=1,,N,(1)whereS2(t)isthetwo-dimensional(2D)shrinkagethatminimizes(s,t)forxedt:S2(t),maxktk)]TJ /F3 11.955 Tf 21.79 8.09 Td[(1 ,0t ktk,t2R2.Moreover,wehave =Sc()(1)whereSc(t)isthecomponentwiseshrinkagethatminimizes (s,t)forxedt=(t1,,tKJ)>2RKJ:Sc(t)=(S(t1),,S(tKJ))>andS(x)=maxfx)]TJ /F3 11.955 Tf 11.95 0 Td[(1=,0gsign(x),x2R,withassumption0(0=0)=0.BothcomputationalcostsforwandarelinearinN.Secondly,forxeduand,wecanhave=(1;;J)bysolvingthefollowingminimizationproblem: minJXj=1)]TJ /F5 11.955 Tf 5.48 -9.68 Td[(kj)]TJ /F5 11.955 Tf 11.95 0 Td[(jk2+kAj)]TJ /F4 11.955 Tf 11.96 0 Td[(Rjuk2.(1)Thesolutioncanbeobtainedbysettingjas j=V(I+))]TJ /F6 7.97 Tf 6.59 0 Td[(1V>)]TJ /F5 11.955 Tf 5.48 -9.68 Td[(j+A>Rju(1)wherethediagonalmatrixandorthogonalmatrixVcomefromtheeigendecompositionA>A=VV>.Thisdecompositiondoesnotdragthecomputationsincethedictionary 24

PAGE 25

Aispreparedbeforeanyexperiments,andhence,Vandcanbepre-computed.AlsothelargestdimensionKofAisusuallymuchlessthanN,anditcanbeseenthatthenumberofnonzeroeigenvaluescanneverexceedn.Asaresult,thecomputationsofj'scanbecarriedoutinparallel,andeachofthemonlyinvolvesmatrix-vectormultiplication.Thirdly,forxedw,and,theminimizationofuis minukwx)]TJ /F4 11.955 Tf 11.96 0 Td[(Dxuk2+kwy)]TJ /F4 11.955 Tf 11.95 0 Td[(Dyuk2+JXj=1kAj)]TJ /F4 11.955 Tf 11.95 0 Td[(Rjuk2+kFpu)]TJ /F4 11.955 Tf 11.96 0 Td[(fpk2,(1)HereDx,DyareN-squarematricesformedbythetopandbottomrowsofDi2R2N,i=1,,N,andhenceDxu,Dyurepresentthegradientofualongthexandydirections,respectively.wxandwyaretherstandsecondcolumnofw,respectively,and==,=()==2.ThustheEuler-Lagrangeequationof( 1 )yields Lu=r,(1)whereL=D>xDx+D>yDy+JXj=1R>jRj+F>pFpandr=D>xwx+D>ywy+JXj=1R>jAj+F>pfp.Undertheperiodicboundaryconditionforu,thenitedifferenceoperatorsDxandDyareblockcirculantmatriceswithcirculantblocksandhencecanbediagonalizedbyFouriermatrixF.Thus,^Dx=FDxF>and^Dy=FDyF>arediagonal.Also,periodicboundaryconditionenablesustoextractpatchesthatcovereachpixelmtimes,wherem=n=d2anddistheslidingdistancebetweenallconcatenatedpatches.Usually,wexthepatchsizeas88,i.e.n=64,andchoose(d,m)=(8,1)or(4,4).SincePjR>jRjisadiagonalmatrixwiththei-thdiagonalentrycountingthenumberoftimesthei-th 25

PAGE 26

pixelwascoveredbypatches,wehavePjR>jRj=mI,whereIistheidentitymatrix.SomultiplyingFonbothsidesof( 1 )gives ^LF(u)=^r,(1)where^L=^D>x^Dx+^D>y^Dy+mI+P>PisadiagonalmatrixsinceP>Pisdiagonal,and^r=^D>xF(wx)+^D>yF(wy)+F(u)+P>fpwhereu=PjR>jAjisanimageassembledusingpatchesthatarerepresentedbydictionaryAand.Finally,thecomputationofrstvariationofF(u,,fp)givesanupdateofineachiteration: =p kFpu)]TJ /F4 11.955 Tf 11.96 0 Td[(fpk2=p.(1)Therefore,similartothealgorithmdevelopedin[ 67 ],themaincomputationofouralgorithmfortheCSmodelusingadictionaryalsoinvolvesonlyshrinkage,matrix-vectormultiplicationandfastFouriertransform. 1.3.2NumericalAlgorithmandConvergenceAnalysisBasedonthediscussionabove,wearereadytoproposethefastalgorithmusedforsolvingmodel( 1 ).Notethatthesolutionto( 1 )canbeapproximatedbysolving( 1 )withcontinuationonthepenaltyparametersand[ 36 ].Forstoppingcriterion,weletresbethemaximumabsolute/normvalueofincrementsofw,,,u,andterminateeachinnerlooponceres
PAGE 27

of( 1 ).Basedonderivationsabove,wesummarizethealgorithmforourmodelasAlgorithm1. Algorithm1MRImageReconstructionviaSparseRepresentation(recMRI) InputP,fp,and,,,>0.Initializeu=F>pfp,==26and=0. while,<212do repeat Givenuand,computewandusing( 1 )and( 1 ). forj=1toJdo Givenuand,computejusing( 1 ). endfor Givenwand,computeubysolving( 1 )andupdateby( 1 ). untilres< returnu, u u,,(,) (2,2) endwhile Theproofoftheconvergenceoftheproposedalgorithm 5 issimilartotheonegivenin[ 64 ]withslightmodications,andthusisomittedhere. 1.4ExperimentalResultsInthissection,wepresenttheexperimentalresultsoftheproposedmodel( 1 )usingAlgorithm 5 andthecomparisonswiththatresultingfromusingwavelettransformonbothsyntheticandinvivoMRdata.AllimplementationsinvolvedintheexperimentswerecodedinMATLABrv7.3(R2006b),excepttheshrinkageandwavelettransformoperators,whichwerecodedinC++.ComputationswereperformedonaGNU/Linux(version2.6.16)workstationwithIntelrCore2CPUat1.86GHzand2GBmemory. 1.4.1ImprovementonAccuracybyUsingDictionariesToshowtheimprovementontheaccuracyofreconstructionsbyusingdictionariesassparsifyingtransforms,weappliedmodel( 1 ),whichuseswaveletasthesparsifyingtransform,andtheproposedmodel( 1 )tothreedatasets.Thethreedatasetsare:thedefaultShepp-LoganphantomimageprovidedbyMATLABR,a2DaxialbrainMRimageanda2DchestMRimage.ThesamplingmasksusedforthesethreedatasetsaredepictedinFigure 3-1 ,wherewhitepixelsindicatesampledlocations 27

PAGE 28

ink-space.InFigure 3-1 ,themasksarefor(a)phantom(b)brainimageand(c)chestimage,respectively.Weusedpseudoradialmaskforphantom,andCartesianmaskundersamplingphaseencoding(PE)linesforinvivodatatosimulaterandomacquisition.Allofthek-spacedatainthesimulatedpseudo-radialtrajectoryislocatedonCartesiangrid.Inpractice,CS-MRIalgorithmprefersrandomacquisitiontrajectorythatcanleadtoincoherentartifactsaliasing.However,thetrajectoryoftheacquireddataineachechotime(TE)islimitedbytheMRsystem,andhencetruerandomacquisitionisnotpossible.Inrecentyears,acquisitionschemesthatarefeasibleandcanproduceincoherentaliasingartifactsaredeveloped,e.g.radialandspiraltrajectories.Inourexperiments,forsimplicity,weusedpseudosamplingmaskswhichcansimulaterandomnessinacquisitionfordemonstrationpurpose. ASR=8.4% BSR=34.0% CSR=25.0% Figure1-3. Samplingmasksusedwithsamplingratios(SR). 1.4.1.1ResultsofPhantomReconstructionThedefaultShepp-Loganphantomofsize256256isshowninFigure 1-4A .Thenafullk-spacedatawassimulatedbythe2DFastFouriertransform(fft2inMATLAB)ofthephantom.Weusedthepseudoradialmask(showninFigure 1-3A )tothefullk-spacedataandaddedcomplexvaluedGaussianrandomnoisewithstandarddeviation(ofmagnitude)3.2tosimulatethepartiallyacquireddatafp.Directusing 28

PAGE 29

FFTofzerollingunscannedk-spacelocationsresultsinnotoriousartifactaliasing,asshowninFigure 1-4B .Thenweappliedmodel( 1 )withHaarwaveletandmodel( 1 )withanovercompletediscretecosinetransform(DCT)consistingof256atomsofsize88[ 2 ]asthedictionarytothepartialdatafp.Theparametersweusedforbothmodelswere(,,)=(1,103,10)]TJ /F6 7.97 Tf 6.58 0 Td[(3)andtheparameterwassetto1.Figure 1-4 showsthefollowing:(a)Referenceimage.(b)Zero-lling.(c)Reconstructionbymodel( 1 )withHaarwavelet(d)ReconstructionusingovercompleteDCTdictionary.TheresultsbyusingwaveletanddictionaryhavecorrespondingRMSEsare2.47%and2.18%,andSNRare32.9and34.5,respectively. 1.4.1.2ResultsofBrainImageReconstructionThesecondtestisonanaxialbrainMRimage.The2-dimentionalmulti-slicedatasetwascollectedona3TGEsystem(GEHealthcare,Waukesha,Wisconsin,USA)usingtheT1FLAIRsequence(FOV220mm,matrixsize512512,TR3060ms,TE126ms,ipangle90,slicethickness5mm,numberofaverages1)withan8-channelheadcoil(InvivoCorporation,Gainesville,FL,USA).Phaseencodingdirectionwasanterior-posterior.TheimagingprocessoutputahighresolutionandSNRimageofsize512512,whichwereusedasthereferenceimageinourexperiment,asshowninFigure 1-5A .Thesimulatedfullk-spacedatainourexperimentwasobtainedbytheFouriertransformofthisreferenceimage,andthenwasarticiallyundersampledusingtheCartesianmaskshowninFigure 1-3B ,whichledtothepartialk-spacedatafp.Thezoomed-inofthesquareareainFigure 1-5A isshowninFigure 1-5B .Figure 1-5 showsthefollowing:(a)Reference.(b)Zoomed-inofsquareareainthereferenceimage.(c)Zoomed-inofreconstructionbyzero-llingunscannedk-spacelocations.(d)Zoomed-inoflowresolution(LR)imagereconstructedby34.0%centralPElines,withRMSE10.32%andSNR17.2.(e)Zoomed-inofreconstructionobtained 29

PAGE 30

AReference BZero-Filling CWaveletRec. DDictionaryRec. Figure1-4. Reconstructedphantomimagefromsimulatedpartialk-spacedata. usingwaveletassparsifyingtransform,RMSEis8.52%andSNRis20.7.(f)Zoomed-inofreconstructionobtainedusingtraineddictionaryshowninFigure 1-1 ,RMSEis7.74%andSNRis22.0.Withonly34.0%dataforreconstruction,strongaliasingartifactscanbeobservedintheimagereconstructedbyzero-llingunscannedk-spacelocations,andthezoomed-inisshowninFigure 1-5C 30

PAGE 31

Inthisexperimentthedatabaseusedfortrainingadictionaryconsistsof409688patchesextractedrandomlyfromfour2DbrainMRimagesofdifferentnormalsubjects(excludingtheonetobereconstructedone)withthesameacquisitionsequence.ThetraineddictionaryA2R64512,asshowninFigure 1-1 ,consistsof512atomsofsize88.Weappliedmodel( 1 )withHaarwaveletandmodel( 1 )withthistraineddictionarytothepartialFourierdatafpforbrainMRimagereconstruction.Theparameters,andweresetto2e+3,1and5e-4inmodel( 1 )and( 1 ),respectively.Theparameterinmodel( 1 )wassetto106.Thezoomed-inareaofthereconstructedimagesbymodel( 1 )andproposedmodel( 1 )areshowninFigure 1-5E and 1-5F ,respectively.TheRMSEsofreconstructionsare8.52%formodel( 1 )and7.74%forproposedmodel( 1 ),andSNRsare20.7and22.0,respectively.Itcanbeseenthattheimagereconstructedbymodel( 1 )hasoil-paintingeffect.Onthecontrary,theimagereconstructedbytheproposedmodelhasbetterpreservednestructures.Thisfurtherconrmsthehigheraccuracyobtainedbytheproposedmethod.Wealsosimulatedalowresolution(LR)imagebyusingthe34.0%centralPElines(i.e.allwhiteverticallylinesinthemiddleinFigure 1-3B ),whichhasRMSE10.32%andSNR17.2. 1.4.1.3ResultsofChestImageReconstructionWealsovalidatetheproposedmethodonchestMRimages.Inthisexperimentthedictionarywastrainedbyslicesextractedfromathree-dimensional(3D)MRchestdataset,thatconsistsof19adjacent2Dimageslicesofsize256256nearthethoraxregion.Ourprocedureoftrainingadictionaryisasfollows:werandomlychosefourslicesanddecomposedthemintonon-overlappingpatchesofsize88,anddiscardedthosepatcheswithhomogeneousintensities,andthenuseK-SVDontheremainingpatchestotrainadictionaryA2R82256with256atoms.Itisworthnotingthatifadditionaldata(e.g.chestMRdatascannedfromdifferentsubjectsbutthesamesequence,which 31

PAGE 32

areusuallyavailableinclinicalapplications)wereavailable,onecanreadilyconstructadictionarythatiscomparabletotheonewetrainedusingadjacentslices,andobtainsimilarreconstructionresultsasweshowedbelow.Todemonstratetheimprovedaccuracybyusingdictionary,werandomlychosea2DsliceshowninFigure 1-6A ,whichisdifferentfromthefourusedastrainingslices.Thenwearticiallydownsampleditsk-spacedatausingaCartesianmaskwith25%samplingratio,asshowninFigure 1-3C .Zero-llingtheunscannedk-spacelocationsresultsinsevereartifactsasshowninFigure 1-6B withRMSE15.59%.Weagainsimulatedalowresolution(LR)image,asshowninFigure 1-6D ,byusingthe25.0%centralPElines,whichhasRMSE14.44%andSNR15.4.FromthecorrespondingerrormapFigure 1-6G ,i.e.theabsolutedifferencetothereferenceimage,wecanseethepotentiallossofedgesanddiagnosticinformation.Thereconstructionsperformedbyusingmodel( 1 )withHaarwaveletandmodel( 1 )withthetraineddictionaryareshowninFigure 1-6E and 1-6F ,andthecorrespondingRMSEsare12.04%and8.48%,andSNRsare17.3and20.1.Inbothcases,andweresetto1e+4,2.5and1e-5,respectively.Parameterinmodel( 1 )wassetto105.TheerrormapsofthesetworeconstructionsareshowninFigure 1-6H and 1-6I ,respectively.Itcanbeseenthattheimagereconstructedbytheproposedmodel( 1 )haslowerartifactslevelandbetterpreservededges.Thisexperimentdemonstratesagaintheadvantagesofusingpriorinformationtodenethesparsifyingtransformation,whichresultsinhigheraccuracyofreconstructions.InFigure 1-6 :(a)Referenceimage.(b)Zero-llingunscannedlocations,RMSEis15.59%andSNRis14.2.(d)LRimageobtainedbyusingcentralPElines,RMSEis14.44%andSNRis15.4.(e)Reconstructionbymodel( 1 )withHaarwavelet,RMSEis12.09%andSNRis17.1.(f)Reconstructionbyproposedmodel( 1 )withtraineddictionary,RMSEis8.48%andSNRis20.1.Figures(g),(h)and(i)arecorrespondingerrormapsof(d),(e)and(f)tothereferenceimage(a). 32

PAGE 33

1.4.2ImprovementonRobustness(toParameterSelection)andEfciencyTodemonstratetheimprovementontherobustnessoftheproposedmodel( 1 )withrespecttothechoiceofparameter,wetestedthereconstructionoftheShepp-Loganphantomusingmodel( 1 )withSSDandtheproposedmodel( 1 )usingMLEasasdataconsistencymeasuresonvariouschoicesof.TheresultingRMSEsareshowninTable 1-1 .FromthechangesofRMSEs,wecanseethemodelwithMLEasdataconsistencymeasuregeneratedsimilarlygoodresultswhereasthemodelwithSSDfailedwhenwentimproperlylarge.ThisresultshowsthattheproposedmodelwithMLEdataconsistencymeasureismuchlesssensitivetothechoiceofandhencemakesthereconstructionmorerobust.Table 1-1 alsoshowstheCPUtime(inseconds)forphantomreconstructionusingthreedifferentalgorithms:nonlinearconjugategradient(CG)algorithmformodel( 1 ),algorithm 5 recMRIformodel( 1 )withtheterminvolvingdictionaryJXj=1kjk1+ 2kAj)]TJ /F4 11.955 Tf 11.96 0 Td[(Rjuk2replacedbyHaarwavelettermk>uk1,andalgorithm 5 recMRIformodel( 1 )withovercompleteDCTasdictionary.HerethewavelettransformsweregeneratedusingoptimizedDWTpackageforMATLABRItcanbeseenthattheproposednumericalmethodwasover2.6timesfasterthanconjugategradientbasedmethod.ThedictionarybasedsparserepresentationconsistentlyproducedimageswithlowerRMSEthanwaveletbasedmethod,butittakeslongerreconstructiontimeduetotheredundancyofdictionaries.Moreover,usingoptimizeddiscretewavelettransform(DWT)packageinMATLABRmakesthecomputationforwaveletbasedmodelevenfaster.Therefore,weexpectanimprovementonspeedbycodeoptimizationwhenusingdictionaries. 1.4.3RobustnessofDictionaryTrainingandReconstructionIntheexperimentonbrainMRimagereconstructioninsection 1.4.1.2 ,thepatchesusedbyK-SVDalgorithmwererandomlyextractedfromthefourtrainingimages. 33

PAGE 34

Table1-1. Comparisonofresultsofphantomreconstructionsusingnonlinearconjugategradient(CG)formodel( 1 )withHaarwavelet,algorithmrecMRIwithHaarwavelet,andrecMRIformodel( 1 )withovercompleteDCTdictionary. MethodCG(Wavelet*)recMRI(Wavelet*)recMRI(Dictionary) RMSEObjCPURMSEObjCPURMSEObjCPU 1e+25.93%12.1386.37.93%13.9828.25.21%11.792111e+32.47%11.9271.62.52%12.1127.72.18%10.921991e+45.05%3.14771.44.98%3.02526.93.47%2.5401981e+525.9%2.27187.15.93%1.37527.03.67%1.1162011e+637.0%2.16581.26.16%1.12928.75.52%1.091212 DifferentpatchesmayleadtodifferenttraineddictionaryusingK-SVDalgorithm,andhencemayimpacttheconsequentreconstructions.ThereforeitisimportanttoverifythatthedictionarytrainingandreconstructionprocessarerobusttocertainlevelofvariabilityonthetrainingdatabaseusedinK-SVDalgorithm.Inthisexperiment,werepeated10timesoftheentireprocessfromformingadatasetoftrainingimagetousingtraineddictionaryinbrainMRimagereconstructionasdescribedinsection 1.4.1.2 .Thedifferenceisthat,ineachrun,the4096trainingimagesarerandomlychosenfromapoolofpatches(around50,000)extractedfromimagesacquiredunderthesamesequenceasthatusedfortheimagetobereconstructed.Therefore,thetrainingpatchesusedinonerunaredifferentfromthoseinanother.TheRMSEsandSNRsofreconstructionresultsareshowninTable 1-2 .Meanwhile,whenwedirectly Table1-2. Experimentalresultsof10runsofdictionarytrainingandbrainimagereconstructionasin 1.4.1.2 Runs12345678910 RMSE(%)7.747.777.657.817.837.787.797.797.727.68SNR22.021.822.621.621.621.821.721.722.122.5 usedtheovercompleteDCTinsteadofthetraineddictionaryinthereconstruction( 1 ),thereconstructionhadRMSE8.27%andSNR21.1.Table 1-2 indicatesthattheK-SVDalgorithmandtheconsequentreconstructionsusingtraineddictionariescanconsistentlygenerategoodresultsdespitethatthetrainingpatchesmayvary.Therefore,theproposedschemeusingdictionariestrainedbyK-SVDalgorithminMRimagereconstructionisstableandrobust,andhencehasgreatpracticalpotential. 34

PAGE 35

AReference BZoomed-InReference CZoomed-InZero-Filling DZoomed-InLR EZoomed-InWavelet FZoomed-InDictionary Figure1-5. ReconstructionofbrainMRImage. 35

PAGE 36

AReference BZero-Filling DLRRec. EWaveletRec. FDictionaryRec. GLRerrormap HWaveleterrormap IDictionaryerrormap Figure1-6. ReconstructionsofchestMRimageusingmodel( 1 )and( 1 ). 36

PAGE 37

CHAPTER2COMPUTATIONALACCELERATIONFORMRIMAGERECONSTRUCTIONINPARTIALLYPARALLELIMAGINGOutline Inthispaper,wepresentafastnumericalalgorithmforsolvingtotalvariationand`1(TVL1)basedimagereconstructionwithapplicationinpartiallyparallelMRimaging.Ouralgorithmusesvariablesplittingmethodtoreducecomputationalcost.Moreover,theBarzilai-Borweinstepsizeselectionmethodisadoptedinouralgorithmformuchfasterconvergence.Experimentalresultsonclinicalpartiallyparallelimagingdatademonstratethattheproposedalgorithmrequiresmuchfeweriterationsand/orlesscomputationalcostthanrecentlydevelopedoperatorsplittingandBregmanoperatorsplittingmethods,whichcandealwithageneralsensingmatrixinreconstructionframework,togetsimilarorevenbetterqualityofreconstructedimages. 2.1BackgroundsinTotalVariationBasedImageReconstructionInthispaperwedevelopanovelalgorithmtoacceleratethecomputationoftotalvariation(TV)and/or`1basedimagereconstruction.Thegeneralformofsuchproblemsis minukukTV+k>uk1+1 2kAu)]TJ /F4 11.955 Tf 11.96 0 Td[(fk2,(2)wherekkTVisthetotalvariation,kk1andkkkk2arethe`1and`2norms(Euclideannorms),respectively.Fornotationsimplicityweonlyconsidertwodimensional(2D)imagesinthispaper,whereasthemethodcanbeeasilyextendedtohigherdimensionalcases.Followingthestandardtreatmentwewillvectorizean(2D)imageuintoone-dimensionalcolumnvector,i.e.u2CNwhereNisthetotalnumberofpixelsinu.Then,the(isotropic)TVnormisdenedby kukTV=ZjDuj=NXi=1kDiuk(2) 37

PAGE 38

whereforeachi=1,,N,Di2R2Nhastwononzeroentriesineachrowcorrespondingtonitedifferenceapproximationstopartialderivativesofuatthei-thpixelalongthecoordinateaxes.In( 2 ),,0(+>0)areparameterscorrespondingtotherelativeweightsofthedatadelitytermkAu)]TJ /F4 11.955 Tf 12.39 0 Td[(fk2andthetermskukTVandkuk1.Model( 2 )hasbeenwidelyappliedtoimagereconstructionproblems.Solving( 2 )yieldsarestoredcleanimageufromanobservednoisyorblurredimagefwhenA=Iorablurringmatrix,respectively.Incompressivesensing(CS)applications,Aisusuallyalargeandill-conditionedmatrixdependingonimagingdevicesordataacquisitionpatterns,andfrepresentstheunder-sampleddata.InCS=[ 1,, N]2CNNisusuallyaproperorthogonalmatrix(e.g.wavelet)thatsparsiesunderlyingimageu. 2.1.1PartiallyParallelMRImagingTheCSreconstructionviaTVL1minimization( 2 )hasbeensuccessfullyappliedtoanemergingMRimagingapplicationknownaspartiallyparallelimaging(PPI).PPIusesmultipleRFcoilarrayswithseparatereceiverchannelforeachRFcoil.Asetofmulti-channelk-spacedatafromeachradiofrequency(RF)coilarrayisacquiredsimultaneously.Theimagingisacceleratedbyacquiringareducednumberofk-spacesamples.Partialdataacquisitionincreasesthespacingbetweenregularsubsequentread-outlines,therebyreducingscantime.However,thisreductioninthenumberofrecordedFouriercomponentsleadstoaliasingartifactsinimages.Therearetwogeneralapproachesforremovingthealiasingartifactsandreconstructinghighqualityimages:imagedomain-basedmethodsandk-spacebasedmethods.Variousmodelsintheframeworkof( 2 )havebeenemployedasimagedomain-basedreconstructionmethodsinPPI[ 15 29 54 57 60 70 73 ].Sensitivityencoding(SENSE)[ 54 55 ]isoneofthemostcommonlyusedmethodsofsuchkind.SENSEutilizesknowledgeofthecoilsensitivitiestoseparatealiasedpixelsresultedfromundersampledk-space. 38

PAGE 39

ThefundamentalequationforSENSEisasfollows:InaPPIsystemconsistingofJcoilarrays,theunder-sampledk-spacedatafjfromthej-thchannelrelatestotheunderlyingimageuby PF(Sju)=fj,j=1,,J,(2)whereFistheFouriertransform,Pisabinarymatrixrepresentingtheunder-samplingpattern(mask),andSj2CNisthesensitivitymapofthej-thchannelinthevectorformasu.ThesymbolistheHadamard(orcomponentwise)productbetweentwovectors.InearlyworksonSENSE,thereconstructionwasobtainedbysolvingaleastsquaresproblem minu2CNJXj=1kFp(Sju))]TJ /F4 11.955 Tf 11.95 0 Td[(fjk2,(2)whereFpistheundersampledFouriertransformdenedbyFp,PF.Denote A=0BBBBBBB@FpS1FpS2...FpSJ1CCCCCCCAandf=0BBBBBBB@f1f2...fJ1CCCCCCCA,(2)whereSj,diag(Sj)2CNNisthediagonalmatrixwithSj2CNonthediagonal,j=1,,J.Thenproblem( 3 )canbeexpressedas minu2CNkAu)]TJ /F4 11.955 Tf 11.96 0 Td[(fk2,(2)andthensolvedbyconjugategradient(CG)algorithm.However,duetotheill-conditioningoftheencodingmatrixA,ithasbeenshownin[ 57 ]thattheCGiterationsequenceoftenexhibitsasemi-convergencebehavior,whichcanbecharacterizedasinitiallyconvergingtowardtheexactsolutionandlaterdiverging.Moreover,theconvergencespeedislow,whentheaccelerationfactorishigh.Recently,totalvariation(TV)basedregularizationhasbeenincorporatedintoSENSEtoimprovereconstructedimagequalityandconvergencespeedoverthe 39

PAGE 40

un-regularizedCGmethod([ 15 73 ]).TVbasedregularizationcanbealsoconsideredasforcingthereconstructedimagetobesparsewithrespecttospatialnitedifferences.ThissparsityalongwiththesparsityofMRsignalsunderwavelettransformshavebeenexploitedin[ 44 ],wheretheframework( 2 )hasbeenemployedtoreconstructMRimagesfromunder-sampledk-spacedata.Therehavebeenseveralfastnumericalalgorithmsforsolving( 2 )thatwillbebrieyreviewedinthenextsection.However,computationalaccelerationisstillanimportantissueforcertainmedicalapplications,suchasbreath-holdingcardiacimaging.FortheapplicationinPPIthecomputationalchallengingisnotonlyfromthelackofdifferentiabilityoftheTVand`1terms,butalsotheinversionmatrixAin( 2 )whichhaslargesizeandisseverelyill-conditioned.Themaincontributionofthispaperistodevelopafastnumericalalgorithmforsolving( 2 )withgeneralA.TheproposedalgorithmincorporatestheBarzilai-Borwein(BB)methodintoavariablesplittingframeworkforoptimalstepsizeselection.Thenumericalresultsonpartiallyparallelimaging(PPI)problemsdemonstratemuchimprovedperformanceonreconstructionspeedforsimilarimagequality. 2.1.2PreviousWorkInreviewingthepriorworkonTVL1-basedimagereconstruction,wesimplify( 2 )bytaking=0.ItisworthpointingoutherethatTVhasmuchstrongerpracticalperformancethan`1inimagereconstructions,yethardertosolvebecausethegradientoperatorsinvolvedarenotinvertibleasinthe`1term.With=0,theimagereconstructionisequivalenttosolvingtheproblem minu2CN(NXi=1kDiuk+1 2kAu)]TJ /F4 11.955 Tf 11.95 0 Td[(fk2).(2)Earlyworkonalgorithmsfor( 3 )usedgradientdescentmethodswithexplicit[ 59 ]orsemi-implicitschemes[ 43 63 ]inwhichtheTVnormwasreplacedbyasmooth 40

PAGE 41

approximation kukTV,"=NXi=1p kDiuk2+".(2)However,thechoiceof">0iscrucialtothereconstructionresultsandconvergencespeed.Alarge"encouragesfastconvergencerate,butfailstopreservehighqualitydetailssuchasedgesintherestoredimage;asmall"betterpreservesnestructureinthereconstructionattheexpenseofslowconvergence.In[ 64 67 ],amethodisdevelopedbasedonthefollowingreformulationof( 3 ): minu,w(NXi=1kwik+1 2kAu)]TJ /F4 11.955 Tf 11.96 0 Td[(fk2:wi=Diu,8i)(2)Thenthelinearconstraintwastreatedwithaquadraticpenalty minu,w(NXi=1kwik+kDu)]TJ /F4 11.955 Tf 11.96 0 Td[(wk2+1 2kAu)]TJ /F4 11.955 Tf 11.95 0 Td[(fk2),(2)wherew2C2Nisformedbystackingthetwocolumnsof(w1,,wN)>,andD=(Dx;Dy)2C2NN.DxandDyarethehorizontalandverticalglobalnitedifferencematrices(N-by-N),i.e.theyconsistoftherstandsecondrowsofallDi's,respectively.Foranyxed,( 3 )canbesolvedbyalternatingminimizations.IfbothD>DandA>AcanbediagonalizedbytheFouriermatrix,astheywouldifAiseithertheidentitymatrixorablurringmatrixwithperiodicboundaryconditions,theneachminimizationinvolvesshrinkageandtwofastFouriertransforms(FFTs).Acontinuationmethodisusedtodealwiththeslowconvergencerateassociatedwithalargevaluefor.Themethod,however,isnotapplicabletomoregeneralA.In[ 34 ]GoldsteinandOsherdevelopedasplitBregmanmethodfor( 3 ).Theresultingalgorithmhassimilarcomputationalcomplexitytothealgorithmin[ 64 ];theconvergenceisfastandtheconstraintsareexactlysatised.LaterthesplitBregmanmethodwasshowntobeequivalenttothealternatingdirectionmethodofmultipliers 41

PAGE 42

(ADMM)[ 13 28 32 33 ]appliedtotheaugmentedLagrangianL(w,u;p)denedby NXi=1kwik+1 2kAu)]TJ /F4 11.955 Tf 11.96 0 Td[(fk2+hp,Du)]TJ /F4 11.955 Tf 11.96 0 Td[(wi+ 2kDu)]TJ /F4 11.955 Tf 11.96 0 Td[(wk2,(2)wherep2C2NistheLagrangianmultiplier.Nonetheless,thealgorithmsin[ 34 64 67 ]benetfromthespecialstructureofA,andtheyloseefciencyifA>Acannotbediagonalizedbyfasttransforms.TotreatamoregeneralA,theBregmanoperatorsplitting(BOS)method[ 75 ]replaceskAu)]TJ /F4 11.955 Tf 12.74 0 Td[(fk2byaproximal-liketermku)]TJ /F3 11.955 Tf 12.74 0 Td[((uk)]TJ /F5 11.955 Tf -431.03 -23.91 Td[()]TJ /F6 7.97 Tf 6.59 0 Td[(1A>(Auk)]TJ /F4 11.955 Tf 12.7 0 Td[(f))k2forsome>0.BOSisaninexactUzawamethodthatdependsonthechoiceof.TheadvantageofBOSisthatitcandealwithgeneralAanddoesnotrequiretheinversionofA>Aduringthecomputation.However,BOSisrelativelylessefcientthanthemethodpresentedinthispaper,evenifischosenoptimally.ThecomparisonofourmethodwiththeBOSalgorithmwillbepresentedinSection 2.4 .Therearealsoseveralmethodsdevelopedtosolvetheassociateddualorprimal-dualproblemsof( 3 )basedonthedualformulationoftheTVnorm: kukTV=maxp2Xhp,Dui,(2)whereX=fp2C2N:pi2C2,kpik1,8igandpiextractsthei-thand(i+N)-thentriesofp.Consequently,( 3 )canbewrittenasaminimaxproblem minu2CNmaxp2Xhp,Dui+1 2kAu)]TJ /F4 11.955 Tf 11.96 0 Td[(fk2.(2)In[ 20 ],Chanetal.proposedtosolvetheprimal-dualEuler-LagrangeequationsusingNewton'smethod.Thisleadstoaquadraticconvergencerateandhighlyaccuratesolutions;however,thecostperiterationismuchhighersincethemethodexplicitlyusessecond-orderinformationandtheinversionofaHessianmatrixisrequired.In[ 19 ],ChambolleusedthedualformulationoftheTVdenoisingproblem( 3 )withA=I,andprovidedanefcientsemi-implicitgradientdescentalgorithmforthedual.However,themethoddoesnotnaturallyextendtothecasewithmoregeneralA.Recently,Zhuand 42

PAGE 43

Chan[ 78 ]proposedaprimal-dualhybridgradient(PDHG)method.PDHGalternatelyupdatestheprimalanddualvariablesuandp.NumericalresultsshowthatPDHGoutperformsmethodsin[ 19 34 ]fordenoisinganddeblurringproblems,butitsefciencyagainreliesonthefactthatA>Acanbediagonalizedbyfasttransforms.Later,severalvariationsofPDHG,referredtoasprojectedgradientdescentalgorithms,wereappliedtothedualformulationofimagedenoisingproblemin[ 79 ]tomakethemethodmoreefcient.Furtherenhancementsinvolvedifferentstep-lengthrulesandline-searchstrategies,includingtechniquesbasedontheBarzilai-Borweinmethod[ 9 ].Anotherapproachthatcanbeappliedto( 3 )intheimagingcontext( 2 )withageneralAistheoperatorsplitting(OS)method.In[ 51 ]theOSideaof[ 47 ]isappliedtoimagereconstructionincompressedmagneticresonanceimaging.TheOSschemerewrites( 3 )as minuXih(Diu)+1 2kAu)]TJ /F4 11.955 Tf 11.95 0 Td[(fk2(2)whereh(),kk.Thentheoptimalconditionsfor( 2 )are wi2@h(Diu),1D>iwi+1A>(Au)]TJ /F4 11.955 Tf 11.95 0 Td[(f)=0,(2)where@h(z)isthesubdifferentialofhatsomepointzdenedbyaset@h(z),fd2CN:h(y))]TJ /F4 11.955 Tf 11.95 0 Td[(h(z)hd,y)]TJ /F4 11.955 Tf 11.96 0 Td[(zi,8yg.Thetheoryofconjugatedualitygivestheequivalencyy2@h(z),z2@h(y),8y,z,whereh(y),supvfhy,vi)]TJ /F4 11.955 Tf 19.26 0 Td[(h(v)g.Hencetherstconditionin( 2 )canbewrittenas 022h(wi)+(wi)]TJ /F4 11.955 Tf 11.96 0 Td[(ti),ti=2Diu+wi(2)andthentherstoneleadsto wi2@h((ti)]TJ /F4 11.955 Tf 11.95 0 Td[(wi)=2)=@h(ti)]TJ /F4 11.955 Tf 11.95 0 Td[(wi),(2) 43

PAGE 44

wheretheequalityisduetoh()=kk.( 2 )isequivalentto wi=argminwih(ti)]TJ /F4 11.955 Tf 11.96 0 Td[(wi)+1 2kwik2(2)thatprojectstiontotheunitballinR2.Then,combining( 2 )andthelastequalitiesin( 2 )and( 2 ),theOSschemeiteratesthefollowingforaxedpoint(whichisalsoasolutionto( 3 )):8>>>>>>><>>>>>>>:tk+1i=wki+2Diuk,8iwk+1i=argminwiktk+1i)]TJ /F4 11.955 Tf 11.95 0 Td[(wik+1 2kwik2,8iuk+1=1XiD>iwk+1i+1A>(Auk)]TJ /F4 11.955 Tf 11.95 0 Td[(f)+ukOSisefcientforsolving( 2 )withgeneralAwhenalltheparametersarecarefullychosen.Howeveritisstillnotasefcientasourmethodevenunderitsoptimalsettings.ThecomparisonofourmethodwiththeOSalgorithm[ 51 ]willbegiveninSection 2.4 2.1.3OrganizationTherestofthispaperisorganizedasfollows.InSection 2.2 wepresenttheproposedalgorithmwithdetailedderivations.Section 2.3 describesourexperimentdesignandtheclinicaldatausedinthispaper.Section 2.4 comparesouralgorithmtoBOS[ 75 ]andoperatorsplitting[ 51 ]onPPIdata.Finally,weconcludethepaperinSection 2.5 2.2ProposedAlgorithmInthispaper,wedevelopafastandsimplealgorithmtonumericallysolveproblem( 2 ).Notethatthecomputationalchallengeof( 2 )comesfromthecombinationoftwoissues:oneispossiblyhugesizeandoftheinversionmatrixA,andtheotheroneisthenon-differentiabilityoftheTVand`1terms.Asdiscussedearlier,despitethatthereweresomefastalgorithmsproposedrecentlytosolveimagerestorationproblemssimilarto( 2 ),theirefciencyreliesonaveryspecialstructureofAsuchthatA>Acanbediagonalizedbyfasttransforms, 44

PAGE 45

whichisnotthecaseinmostmedicalimagingproblems,suchasthatin( 2 )inPPIapplication.Totacklethecomputationalproblemof( 2 ),werstintroduceauxiliaryvariableswiandzitotransformDiuand >iuoutofthenon-differentiablenorms: minw,z,u(Xikwik+Xijzij+1 2kAu)]TJ /F4 11.955 Tf 11.95 0 Td[(fk2),wi=Diu,zi= >iu,8i=1,,N,(2)whichisclearlyequivalenttotheoriginalproblem( 2 )astheysharethesamesolutionsu.Todealwiththeconstraintsin( 2 )broughtbyvariablesplitting,weformtheaugmentedLagrangiandenedby L(w,z,u;b,c)=Xikwik)]TJ /F5 11.955 Tf 20.59 0 Td[(hbi,wi)]TJ /F4 11.955 Tf 11.95 0 Td[(Diui+ 2kwi)]TJ /F4 11.955 Tf 11.96 0 Td[(Diuk2+Xijzij)]TJ /F5 11.955 Tf 17.93 0 Td[(ci(zi)]TJ /F5 11.955 Tf 11.96 0 Td[( >iu)+ 2jzi)]TJ /F5 11.955 Tf 11.95 0 Td[( >iuj2+1 2kAu)]TJ /F4 11.955 Tf 11.96 0 Td[(fk2,(2)whereb2C2Nandc=(c1,,cN)>2CNareLagrangianmultipliers.Herebi2C2extractsthei-thand(i+N)-thentriesofb.Fornotationsimplicityweusedthesameparameter>0forallconstraintsin( 2 ).ThemethodofmultipliersiteratestheminimizationsofLagrangianLin( 2 )withrespectto(w,z,u)andtheupdatesofthemultipliersbandc: 8>>>>><>>>>>:(wk+1,zk+1,uk+1)=argminw,z,uL(w,z,u;bk,ck)bk+1i=bki)]TJ /F3 11.955 Tf 11.96 0 Td[((wk+1i)]TJ /F4 11.955 Tf 11.95 0 Td[(Diuk+1),8ick+1i=cki)]TJ /F3 11.955 Tf 11.96 0 Td[((zk+1i)]TJ /F5 11.955 Tf 11.95 0 Td[( >iuk+1),8i(2)Itisprovedthatthesequencef(wk,zk,uk)gkgeneratedby( 2 )convergestothesolutionof( 2 )withany>0. 45

PAGE 46

Sincetheupdatesofbkandckaremerelysimplecalculations,wenowfocusontheminimizationofL(w,z,u;bk,ck)in( 2 ).Firstweintroducefunctions 1(s,t)=jsj+(=2)js)]TJ /F4 11.955 Tf 11.95 0 Td[(tj2,s,t2C(2)and 2(s,t)=ksk+(=2)ks)]TJ /F23 11.955 Tf 11.96 0 Td[(tk2,s,t2C2.(2)Bycompletingthesquaresin( 2 ),wendtheequivalency argminw,z,uL(w,z,u;bk,ck)argminw,z,u(Xi2(wi,Diu+bki)+Xi1(zi, >iu+cki)+1 2kAu)]TJ /F4 11.955 Tf 11.96 0 Td[(fk2(2)becausetheobjectivefunctionsinthesetwominimizationsareequaluptoaconstantindependentof(w,z,u).Tosolve( 2 )werstrewritetheobjectivefunctioninasimplerway.Letx=(w;z;u)andB=(0,0,A),anddenefunctionsJk(x)Jk(w,z,u)by Jk(x),Xi2(wi,Diu+bki)+Xi1(zi, >iu+cki),(2)anddatadelityH(x)by H(x)=(1=2)kBx)]TJ /F4 11.955 Tf 11.96 0 Td[(fk2.(2)Thenproblem( 2 )(orequivalently,theminimizationsubproblemin( 2 ))canbeexpressedas xk+1=argminxfJk(x)+H(x)g,(2)WefurtherintroduceQ(x,y)denedby Q(x,y),H(y)+hrH(y),x)]TJ /F4 11.955 Tf 11.95 0 Td[(yi+ 2kx)]TJ /F4 11.955 Tf 11.96 0 Td[(yk2,(2) 46

PAGE 47

whichisalinearizationofH(x)atpointyplusaproximitytermkx)]TJ /F4 11.955 Tf 12.62 0 Td[(yk2=2penalizedbyparameter>0.Ithasbeenshownin[ 66 ]thattheiterativesequencefxk+1,lglgeneratedby xk+1,l+1=argminxJk(x)+Qk+1,l(x,xk+1,l)(2)convergestothesolutionxk+1of( 2 )withanyinitialxk+1,0andproperchoiceofk+1,lforl=0,1,1.Interestingly,wefoundthatinpracticetheoptimalperformancecanbeconsistentlyachievedifonlyiterating( 2 )oncetoapproximatethesolutionxk+1in( 2 ).Therefore,wesubstitutetherstsubproblemin( 2 )by xk+1=argminxJk(x)+Qk(x,xk),(2)wherekischosenbasedontheBarzilai-Borwein(BB)methodassuggestedin[ 66 ].BBmethodhandlesill-conditioningmuchbetterthangradientmethodswithaCauchystep[ 3 ].IntheBBimplementation,theHessianoftheobjectivefunctionisapproximatedbyamultipleoftheidentitymatrix.Weemploytheapproximation k=argmin)]TJ /F2 11.955 Tf 5.48 -9.69 Td[(rH(xk))-222(rH(xk)]TJ /F6 7.97 Tf 6.59 0 Td[(1))]TJ /F5 11.955 Tf 11.95 0 Td[((xk)]TJ /F4 11.955 Tf 11.95 0 Td[(xk)]TJ /F6 7.97 Tf 6.59 0 Td[(1)2,(2)andget k=hrH(xk))-222(rH(xk)]TJ /F6 7.97 Tf 6.58 0 Td[(1),xk)]TJ /F4 11.955 Tf 11.96 0 Td[(xk)]TJ /F6 7.97 Tf 6.59 0 Td[(1i=kxk)]TJ /F4 11.955 Tf 11.96 0 Td[(xk)]TJ /F6 7.97 Tf 6.59 0 Td[(1k2.(2)Thismakestheiteration( 2 )exhibitacertainlevelofquasi-Newtonconvergencebehavior. 1e.g.forxedk,anylimitpointoffxk+1,lglisasolutionof( 2 )whenk+1,lwaschosensuchthattheobjectivefunctionJk(xk+1,l)+H(xk+1,l)monotonicallydecreasesasl!1[ 66 ]. 47

PAGE 48

8>>>>>>>>>>>>>><>>>>>>>>>>>>>>:wk+1i=argminwikwik+ 2kwi)]TJ /F4 11.955 Tf 11.96 0 Td[(Diuk)]TJ /F4 11.955 Tf 11.96 0 Td[(bkik2+k 2kw)]TJ /F4 11.955 Tf 11.96 0 Td[(wkk2,8i;zk+1i=argminzijzij+ 2jzi)]TJ /F3 11.955 Tf 11.96 0 Td[(>iuk)]TJ /F4 11.955 Tf 11.96 0 Td[(ckij2+k 2jzi)]TJ /F4 11.955 Tf 11.95 0 Td[(zkij2,8i;uk+1=argminunkDu)]TJ /F4 11.955 Tf 11.96 0 Td[(wk+1k2+k>u)]TJ /F4 11.955 Tf 11.96 0 Td[(zk+1k2+ku)]TJ /F13 11.955 Tf 11.95 9.69 Td[()]TJ /F4 11.955 Tf 5.48 -9.69 Td[(uk)]TJ /F5 11.955 Tf 11.96 0 Td[()]TJ /F6 7.97 Tf 6.59 0 Td[(1kA>(Auk)]TJ /F4 11.955 Tf 11.96 0 Td[(f)2o;bk+1i=bki)]TJ /F3 11.955 Tf 11.96 0 Td[((wk+1i)]TJ /F4 11.955 Tf 11.95 0 Td[(Diuk+1),8i;ck+1i=cki)]TJ /F3 11.955 Tf 11.96 0 Td[((zk+1i)]TJ /F5 11.955 Tf 11.96 0 Td[( iuk+1),8i;k+1=kA(uk+1)]TJ /F4 11.955 Tf 11.96 0 Td[(uk)k2=)]TJ /F2 11.955 Tf 5.48 -9.68 Td[(kwk+1)]TJ /F4 11.955 Tf 11.95 0 Td[(wkk2+kzk+1)]TJ /F4 11.955 Tf 11.96 0 Td[(zkk2+kuk+1)]TJ /F4 11.955 Tf 11.95 0 Td[(ukk2.(2) FromthedenitionofJkandQk,( 2 )isequivalentto (wk+1,zk+1,uk+1)=argminw,z,uk(w,z,u)(2)wheretheobjectivek(w,z,u)isdenedby k(w,z,u),Xi2(wi,Diu+bki)+Xi1(zi, >iu+cki)+k 2kw)]TJ /F25 9.963 Tf 9.96 0 Td[(wkk2+kz)]TJ /F25 9.963 Tf 9.96 0 Td[(zkk2+u)]TJ /F25 9.963 Tf 9.96 0 Td[(uk+)]TJ /F28 6.974 Tf 6.22 0 Td[(1kA>(Auk)]TJ /F25 9.963 Tf 9.96 0 Td[(f)2(2)Theoretically,aniterativeschemecanbeappliedtoobtainthesolution(wk+1,zk+1,uk+1)of( 2 ).However,hereweproposeonlytodooneiterationfollowedbytheupdatesofbk,ckandkin( 2 ).ThisisananaloguetothesplitBregmanmethodandADMMappliedtotheaugmentedLagrangians,andleadstotheoptimalperformanceof( 2 ).Insummary,weproposeaschemeasin( 2 )forsolvingtheminimizationproblem( 2 ).Theupdatesofbk+1andck+1in( 2 )aremerelysimplecalculations.In( 2 ),k+1isderivedfrom( 2 )withHdenedin( 2 ),andalsohasanexplicitformthat 48

PAGE 49

canbequicklycomputed2.Next,weshowthatwk+1iandzk+1icanbeobtainedbysoftshrinkagesbythefollowingtheorem. Theorem2.1. Forgivend-vectorst1,t22Rdandpositivenumbersa1,a2>0,thesolutiontominimizationproblem mins2Rdnksk+a1 2ks)]TJ /F23 11.955 Tf 11.95 0 Td[(t1k2+a2 2ks)]TJ /F23 11.955 Tf 11.96 0 Td[(t2k2o(2)isgivenbytheshrinkageofaweightedsumoft1andt2: Sd(t1,t2;a1,a2),shrinkda1t1+a2t2 a1+a2,1 a1+a2(2)whereshrinkdisthed-dimensionalsoftshrinkageoperatordenedby shrinkd(t,),maxfktk)]TJ /F5 11.955 Tf 20.59 0 Td[(,0gt ktk.(2)withconvention0(0=k0k)=0. Proof. Bycompletingthesquares,theminimizationproblem( 2 )isequivalentto mins2Rd(ksk+a1+a2 2s)]TJ /F4 11.955 Tf 13.15 8.09 Td[(a1t1+a2t2 a1+a22),(2)becausetheobjectivefunctionsarethesameuptoaconstantindependentofs.Minimizationsofform( 2 )haveawellknownexplicitsolvershrinkdandhencetheconclusionfollows. AccordingtoTheorem 2.1 ,wk+1iandzk+1iin( 2 )canbeobtainedby wk+1i=S2)]TJ /F4 11.955 Tf 5.48 -9.68 Td[(Diuk+bki,wki;,k(2) 2Themaincomputationsforupdatingkarenormevaluations(noAoperationneededsinceAukhasbeencomputedintheu-stepandcanbesavedforusein-stepin( 2 )). 49

PAGE 50

and zk+1i=S1)]TJ /F5 11.955 Tf 5.48 -9.68 Td[( >iuk+cki,zki;,k(2)wherek=k=andk=k=.Thereforethecomputationalcostsfor( 2 )and( 2 )arelinearintermsofN.Theu-subproblemin( 2 )isaleastsquaresproblem.Theoptimalconditionofthisproblemreads Lku=rk(2)where Lk=D>D+I+kI(2)and rk=D>wk+1+>zk+1+kuk)]TJ /F4 11.955 Tf 11.96 0 Td[(A>(Auk)]TJ /F4 11.955 Tf 11.95 0 Td[(f).(2)Underperiodicboundarycondition,thematrixD>DisblockcirculantandhencecanbediagonalizedbyFouriermatrixF.Let=F>D>DFwhichisadiagonalmatrix,thenapplyFonbothsidesof( 2 )toobtain ^LkFu=^rk(2)where ^Lk=+I+kIand^rk=Frk.(2)Notethat^Lkcanbetriviallyinvertedbecauseitisdiagonalandpositivedenite.Therefore,uk+1canbeeasilyobtainedby uk+1=F>(^L)]TJ /F6 7.97 Tf 6.59 0 Td[(1kFrk).(2)Asallvariablesin( 2 )canbequicklysolved,weproposeAlgorithm 2 ,calledTVL1rec,tosolveproblem( 2 ).Asdiscussedabove,wandzcanbeupdatedusingsoftshrinkagesandhencethecomputationalcostsarelinearintermsofN.TheupdateofuinvolvestwofastFouriertransforms(FFTs)whichhascomputationalcomplexity 50

PAGE 51

Algorithm2TVL1ReconstructionAlgorithm(TVL1rec) Input,,and.Setu0=c=0,b=0,0=1,k=0. repeat Givenuk,computewk+1andzk+1using( 2 )and( 2 ); Givenwk+1andzk+1,computeuk+1using( 2 ); Updatebk,ckandkasin( 2 ); k k+1 untilkuk)]TJ /F4 11.955 Tf 11.95 0 Td[(uk)]TJ /F6 7.97 Tf 6.59 0 Td[(1k=kukk<. returnuk NlogNandtwooperationsofA(oneisA>).Therefore,unlikemostrecentlydevelopedalgorithms,ouralgorithmcandealwitharbitrarymatrixAandevenmoregeneralHwithnonlinearconstraint(aslongasHisconvexandrHiscomputable).Also,theperiterationcomputationoftheproposedalgorithmisverycheap,andthankstotheBBstepsizek,theconvergencespeedissignicantlyimprovedcomparedtoothertwomodernmethodsBOSandOS,asshowninSection 2.4 2.3MethodExperimentsweredesignedtotesttheeffectivenessoftheproposedalgorithmTVL1reconPPIreconstructions.Todemonstratethepotentialinclinicapplications,thethreedatasetsusedintheexperimentswereacquiredbycommerciallyavailable8-elementheadcoils.Forcomparison,wealsoimplementedtheBregmanoperatorsplittingalgorithm(BOS)[ 75 ]andacompressiveMRimagereconstructionalgorithmbasedonoperatorsplitting(OS)[ 51 ]forsolving( 2 ). 2.3.1DataAcquisitionTherstdataset(topleftinFigure 2-2 ),termedbydata1,isasetofsagittalCartesianbrainimagesacquiredona3TGEsystem(GEHealthcare,Waukesha,Wisconsin,USA).ThedataacquisitionparameterswereFOV220mm2,size5125128,TR3060ms,TE126ms,slicethickness5mm,ipangle90,andphaseencodingdirectionwasanterior-posterior.Theseconddataset(leftinFigure 2-4 )isaCartesianbraindatasetacquiredona3.0TPhilipsscanner(Philips,Best,Netherlands)usingT2-weightedturbospinecho(T2 51

PAGE 52

A B C Figure2-1. Thek-spacemasksusedforthethreedatasetstestedintheexperiments. TSE)sequence.TheacquisitionparameterswereFOV205mm2,matrix5125008,TR3000ms,TE85ms,andtheechotrainlengthwas20.Toavoidsimilarcomparisonplotduetothesamedatasize,wereducetheimageto2562508andobtainfullk-spacedataofthissamesize,termedbydata2.Thelastone(rightofFigure 2-4 ),denotedbydata3,isaradialbraindatasetacquiredona1.5TSiemensSymphonysystem(SiemensMedicalSolutions,Erlangen,Germany).TheacquisitionparameterswereFOV220mm2,matrix2565128(256radiallines),slicethickness5mm,TR53.5ms,TE3.4ms,andipangle75.Allthreedatasetswerefullyacquired,andthenarticiallydown-sampledusingthemasksinFigure 3-1 forreconstruction.Figure 3-1 shows,fromlefttoright,theCartesianmaskwithnetreductionfactor3,pseudorandommaskwithreductionfactor4,andradialmaskwith43(outof256)projections,i.e.reductionfactor6.Thereferenceimagesinourexperimentswereobtainedbyfullyacquiredk-space.AsummaryofthedatainformationisinTable 2-1 .InTable 2-1 ,Cart.Sag.meansCartesiansagittalbrainimage,andRad.Axi.standsforradialaxialbrainimage.ThecolumnPinTable 2-1 presentthemasknumber(refertoFigure 3-1 ). 52

PAGE 53

Table2-1. Testsnumber,datainformationandparameters. No.ImageAbbrev.Size(8)P(,) 1Cart.Sag.data15125121(1e-51e-2,0)2Cart.Sag.data22562502(1e-4,5e-5)3Rad.Axi.data32565123(1e-4,5e-5) 2.3.2TestEnvironmentAllalgorithmswereimplementedintheMATLABrprogrammingenvironment(VersionR2009a,MathWorksInc.,Natick,MA,USA).ThesparsifyingoperatorissettoHaarwavelettransformusingRicewavelettoolboxwithdefaultsettings.TheexperimentswereperformedonaDellOptiplexdesktopwithIntelRDualCore2.53GHzprocessors(only1corewasusedincomputation),3GBofmemoryandWindowsTMoperatingsystem.TheoreticallythechoiceofdoesnoteffecttheconvergenceofTVL1rec.Thisisalsodemonstratedbyourexperimentssincetheresultsarenotsensitivetoforalargerange.Thereforeinallexperimentswesettoamoderatevalue10.Algorithm 2 isautomaticallyterminatediftherelativechangeofukislessthanaprescribedtolerance.Inallofourexperiments,weset=10)]TJ /F6 7.97 Tf 6.59 0 Td[(3.Notethatsmallerleadstoslightlybetteraccuracyatthecostofmoreiterationsandlongercomputationaltime.Otherparametersettingsareshowninthenextsection.Forallalgorithmstestedinthispaper,thesensitivitymapsSj'swereestimatedfromthecentral3232k-spacedata(whichwasasubsetoftheacquiredpartialdata)andthenxedduringthereconstructions,andtheinitialu0wassetto0.Thereconstructionresultswereevaluatedqualitativelybyzoomed-inregionsofthereconstructedimages,andquantitativelybyrelativeerror(tothereferenceimage)andCPUtimes.Referenceandreconstructedimagescorrespondingtodata1anddata3werebrightenedby3times,andthosecorrespondingtodata2werebrightenedby2times,tohelpvisualjustications. 53

PAGE 54

2.4ComparisonAlgorithmsandResults 2.4.1ComparisonwithBOSIntherstexperiment,weusedata1withaCartesiansamplingpattern(leftinFigure 3-1 )toundersamplek-spacedata.WecompareTVL1recwithBOSwhichalsosolves( 2 )viaavariablesplittingframework( 2 ).Tosimplifycomparison,wehereset=0in( 2 )andfocusonthecomputationalefciencyoftwoalgorithmsinsolving( 2 ).TheBOSalgorithmsolves( 2 )byiterating 8>>>>>>>>><>>>>>>>>>:sk+1=uk)]TJ /F5 11.955 Tf 11.96 0 Td[()]TJ /F6 7.97 Tf 6.59 0 Td[(1A>(Auk)]TJ /F4 11.955 Tf 11.95 0 Td[(f)wk+1i=argminwinkwik+ 2kwi)]TJ /F4 11.955 Tf 11.96 0 Td[(Diuk)]TJ /F4 11.955 Tf 11.95 0 Td[(bkik2o,8iuk+1=argminufkDu)]TJ /F4 11.955 Tf 11.96 0 Td[(wk+1+bkk2+u)]TJ /F4 11.955 Tf 11.96 0 Td[(sk+12gbk+1i=bki)]TJ /F3 11.955 Tf 11.96 0 Td[((wk+1i)]TJ /F4 11.955 Tf 11.96 0 Td[(Diuk+1),8i(2)andconvergesifkA>Ak2,i.e.thelargesteigenvalueofA>A.InSENSEapplications,themagnitudesofsensitivitymapsareusuallynormalizedinto[0,1].ThereforefromthedenitionofAin( 2 ),wehavekA>Ak2=1andhenceset=1foroptimalperformanceofBOS.With=0,TVL1reconlyupdatesw,u,bandin( 2 ).Ascanbeseen,theperiterationcomputationalcostsforBOSandTVL1recarealmostidentical:themaincomputationsconsistofoneshrinkage,A,A>andtwoFFTs(includingoneinverseFFT).ThereforethecomputationcostforacompletereconstructionisnearlyproportionaltothenumberofiterationsrequiredbyBOSandTVL1rec.Inthispaper,wesetthestoppingcriterionofBOSthesameasTVL1rec,namelythecomputationwillbeautomaticallyterminatedwhentherelativechangeoftheiterateukislessthan=10)]TJ /F6 7.97 Tf 6.59 0 Td[(3.Table 2-2 showstheperformanceresultsofTVL1recandBOSondata1fordifferentvaluesofTVregularizationparameter.InTable 2-2 ,welistthefollowingquantities:therelativeerrorofthereconstructedimagestothereferenceimage(Err),thenalobjective 54

PAGE 55

functionvalues(Obj),thenumberofiterations(Iter),andtheCPUtimeinseconds(CPU). Table2-2. ResultsofBOSandTVL1recondata1. BOS TVL1rec ErrObjIterCPUErrObjIterCPU 1e-58.1%.2813375.17.2%.252718.61e-47.4%1.011738.97.1%.8601126.71e-37.4%6.003988.27.3%5.98716.01e-211.5%41.063142.110.6%40.7715.9 FromTable 2-2 ,wecanseethatbothBOSandTVL1areabletostablyrecovertheimagefrom34%k-spacedata.ThisisfurtherdemonstratedbyFigure 2-2 ,wherebothmethodgeneratedimagesveryclosetothereferenceimage.DespitethattherearefewobservablealiasingartifactsduetoCartesianundersampling,thedetailssuchasedgesandnestructureswerewellpreservedinbothreconstructions,ascanbeseeninthezoomed-insintherightcolumnofFigure 2-2 .Intermsofaccuracy,TVL1recgivesslightlybetterreconstructionqualityinthesenseoflowerrelativeerrorandobjectivevalues.Intermsofefciency,wefoundthatTVL1recsignicantlyoutperformsBOSbyrequiringmuchfeweriterations(andhencelessCPUtime)toobtainthesimilarorevenbetterimagequality,asshowninTable 2-2 .ComparedtoBOS,TVL1recisupto9timesfasterandhencehasmuchhigherefciency.Althoughtwoalgorithmshavealmostthesamecomputationalcostsperiteration,TVL1recbenetsfromtheadaptivechoiceofstepsizesandreadilyoutperformsBOSwhichusesxedstepsize=kA>Ak2throughoutthecomputations.TheadaptivestepsizeselectionmakesTVL1recexhibitsaquasi-NewtonconvergencebehaviorinsomesensebecausekIimplicitlyusespartialHessian(secondorder)information.TheadaptivestepsizeselectionnotonlyleadstohigherefciencybutalsobetterstablenessofTVL1rec.AsshowninTable 2-2 ,foralargerangeofin[10)]TJ /F6 7.97 Tf 6.58 0 Td[(5,10)]TJ /F6 7.97 Tf 6.58 0 Td[(2],TVL1recalwaysrequires11orfeweriterationstorecoverhighqualityimages.In 55

PAGE 56

comparison,BOSappearstobequitesensitivetothechoiceof:thisisexempliedbythelastrow(=10)]TJ /F6 7.97 Tf 6.59 0 Td[(2)ofTable 2-2 ,whereBOSrequiredmuchmoreiterationsthanusual;meanwhile,TVL1recbenetsfromtheoptimalstepsizeineachiterationandreadilyapproximatesthesolutioninonlyfewiterations.InFigure 2-2 ,therightcolumnshowsthezoomed-in(squareindata1)ofimagesintheleftcolumn.FromtoptobottomofFigure 2-2 ,theyarethereferenceimage,reconstructedimageusingBOS(Err=7.4%),andreconstructedimageusingTVL1rec(Err=7.1%),respectively.ThebetterperformanceofTVL1recoverBOSreliesontwophases:oneisthatTVL1recimposesproximitytermsnotonlyforubutalsoforwandzin( 2 ),whichleadtobetterchoicesoftheupdateswk+1andzk+1;theotheroneistheadoptionofBBmethodforoptimalpenaltyparameterskselection,whichaffectstheupdatesofallvariablesasin( 2 )andleadsmuchimprovedconvergencespeed. 2.4.2ComparisonwithOSFordata2anddata3,wecompareTVL1withOS[ 51 ]forsolving( 2 )withbothTVand`1terms(=10)]TJ /F6 7.97 Tf 6.59 0 Td[(4,==2).TheOSschemeof[ 51 ],withaminorcorrection,isasfollows: 8>>>>>>>><>>>>>>>>:sk+1=)]TJ /F4 11.955 Tf 5.48 -9.68 Td[(uk)]TJ /F3 11.955 Tf 11.96 0 Td[((d1=))]TJ /F4 11.955 Tf 5.48 -9.68 Td[(D>wk+A>(Auk)]TJ /F4 11.955 Tf 11.96 0 Td[(f),tk+1i=wki+d2Diuk,8i,uk+1=)]TJ /F3 11.955 Tf 5.48 -9.69 Td[(sign(sk+1)maxfjsk+1j)]TJ /F4 11.955 Tf 17.93 0 Td[(d1=,0g,wk+1i=min(1,ktk+1ik)tk+1i=ktk+1ik2,8i.(2)wheresk2CN,tkiandwki2C2,i=1,,N,wk2C2Nisformedbystackingthetwocolumnsofmatrix(wk1,,wkN)>,andthemaxandsignoperationsinthecomputationofuk+1arecomponentwiseoperationscorrespondingtoshrinkage.ThemaincomputationalcostperiterationintheOSschemecorrespondstothefollowingoperations:a2Dshrinkageduringthecomputationofwk+1,aprojectionduringthecomputationofuk+1,twowavelettransformsduringthecomputationofsk+1anduk+1,A 56

PAGE 57

andA>duringthecomputationofsk+1.In[ 51 ]itisshownthatford1,d2>0incertainranges,theOSschemeconvergestoaxedpointwhichisalsoasolutionof( 2 ).Theiterationswerestoppedwheneitherthefollowingconditionsweresatised: kuk+1)]TJ /F4 11.955 Tf 11.96 0 Td[(ukk2=maxf1,kukk2g<1(fk)]TJ /F4 11.955 Tf 11.95 0 Td[(fk+1)=maxf1,fkg<2p c=t,(2)wherefkistheobjectivevalueof( 2 )atuk,candtarethecurrentandtargetvaluesofrespectivelyand1and2areprescribedstoppingtolerances.SinceOShasmultipletuningparametersthataffecttheconvergencespeedandimagequality:largerdi'sandi'sleadtofasterconvergencebutresultinlargerrelativeerror,whereassmallerdi'sandi'syieldmonotonicdecreasesinobjectivevaluesandbetterimagequalityatthecostofmuchlongercomputation.Basedontheselectionbytheauthorsandseveraltries,wechosemoderatevaluesd1=d2=1,1=10)]TJ /F6 7.97 Tf 6.58 0 Td[(4and2=10)]TJ /F6 7.97 Tf 6.58 0 Td[(3whichappeartogiveabestcompromisebetweenconvergencespeedandimagequalityoftheOSscheme.Theresultsondata2anddata3areshowninFigure 2-4 ,andthecomparisononrelativeerrorsandobjectivevalues(bothinlogarithmic)areplottedinlogarithmicscaleinFigure 2-5 .InFigure 2-4 ,theleftcolumnandrightcolumncorrespondtotheresultsfordata2anddata3,respectively.Intheleftcolumn(resultsofdata2),theyarereference,reconstructedimagesbyOS(Err=7.6%)andTVL1rec(Err=4.6%).Intherightcolumn(resultsofdata3)fromtoptobottom,theyarereference,reconstructedimagesbyOS(Err=6.7%)andTVL1rec(Err=6.1%).ThehorizontallabelischosenasCPUtimebecausetheperiterationcomputationalcostsforOSandTVL1recareslightlydifferent.FromFigures 2-4 and 2-5 wecanseethatTVL1recconvergesmuchfasterthanOS,andachievedlowerrelativeerrorsandobjectivevaluesthanOSoverall.Therefore,itisevidentthatTVL1reccanoutperformOSschemeinefciencyastheformerrequiresmuchlesscomputationaltimetoreachthesimilarorevenbetterimagequality.Itisalso 57

PAGE 58

worthpointingoutthatbothalgorithmscanfurtherreducetherelativeerrorslightlybysettingatighterstoppingcriterionatthecostofmoreiterationnumbers.Nevertheless,theTVL1recstillcanmaintainlowerrelativeerrorandobjectivevaluethanOSduringthereconstructionprocess. 2.5ConcludingRemarksThispaperpresentsafastnumericalalgorithm,calledTVL1rec,forTVL1minimizationproblem( 2 )arisingfromCSreconstructionproblems.TheproposedalgorithmincorporatestheBarzilai-Borwein(BB)methodintoavariablesplittingframeworktooptimizetheselectionofstepsizes.TheoptimalstepsizesexploitpartialHessianinformationandhenceleadtoaquasi-NewtonconvergencebehaviorofTVL1rec.ExperimentalresultsdemonstratetheoutstandingefciencyoftheproposedalgorithminCS-PPIreconstruction.WecomparedTVL1rectoanothertworecentlydevelopedalgorithmsBOS[ 75 ]andOS[ 51 ]whichalsosolvetheminimizationproblem( 2 ).ThecommonpropertyofthesealgorithmsisthattheycandealwithgeneralsensingmatrixA,andevennonlineardatadelitytermH(u)otherthankAu)]TJ /F4 11.955 Tf 12.37 0 Td[(fk2aslongasHisconvexandrHiscomputable.Meanwhile,TVL1recsignicantlyoutperformstheothertwoalgorithmsbytakingadvantagesoftheoptimalstepsizeselectionbasedonBBmethod.WehopeTVL1reccanbebenecialtoPPIandothermedicalimagingapplications. 58

PAGE 59

A B C D E F Figure2-2. ComparisonofBOSandTVL1recondata1. 59

PAGE 60

A B Figure2-3. TestingdatausedforthecomparisonofOSandTVL1recfordata2anddata3. 60

PAGE 61

A B C D E F Figure2-4. Reconstructionsofdata2anddata3byOSandTVL1rec. 61

PAGE 62

A B Figure2-5. ComparisonsofOSandTVL1recondata2(bluesolidlines)anddata3(blackdashedlines). 62

PAGE 63

CHAPTER3AFASTALGORITHMFORTVIMAGERECONSTRUCTIONWITHAPPLICATIONTOPARTIALLYPARALLELMRIMAGINGOutline Thispaperpresentsafastalgorithmfortotalvariation-basedimagereconstruction.Theproposedmethodcombinesvariablesplittingandtheclassicpenaltytechnique.Thisreducestheimagereconstructionproblemtoanunconstrainedminimizationproblem,whichissolvedbyanalternatingproximalminimizationalgorithm.Onephaseofthealgorithmsolvesatotalvariation(TV)denoisingproblem,andsecondphasesolvesanill-conditionedlinearsystem.Linearandsublinearconvergenceresultsaregiven,andanimplementationbasedonasplitBregmanschemefortheTVproblemandaBarzilai-Borweinmethodforthelinearsystemisproposed.Thealgorithmisappliedtoimagereconstructionproblemsthatariseinpartiallyparallelmagneticresonanceimaging(PPI).Performanceiscomparedtothatofanoperatorsplittingscheme. 3.1BackgroundsinOptimizationMethodsforNonsmoothBasedImageReconstructionInthispaperwedevelopanewalgorithmfortotalvariation(TV)basedimagereconstruction.Thegeneralformofsuchproblemsis minu2CNJ(u)+H(u),(3)whereJisaconvexandpossiblynondifferentiablefunction,andHisconvexandcontinuouslydifferentiable.InTV-basedimagereconstructionproblems,JandHoftenhavetheform J(u)=kukTV+kuk1andH(u)=kAu)]TJ /F4 11.955 Tf 11.95 0 Td[(fk22,(3)wherekkTVisthetotalvariation,kk1isthe1-norm,kk2isthe2-norm(Euclideannorm),Aisapossiblylargeandill-conditionedmatrixdescribingtheimagingdeviceorthedataacquisitionpattern,fisthemeasureddata,and2CNNisanorthogonal 63

PAGE 64

sparsifyingmatrix.and0areparameterscorrespondingtotherelativeweightsofthedatadelitytermkAu)]TJ /F4 11.955 Tf 12.63 0 Td[(fk22andthetermskukTVandkuk1whichcontrolthesolutionsparsity.projectsualongasetofbasisfunctionsdescribedby,forexample,waveletsoradictionary[ 50 69 ].Itisexpectedthatmanycomponentsofuwillvanishatasolutionof( 4 )( 3 )formagneticresonanceimages.TV-basedregularizationwasoriginallyintroducedbyRudin,OsherandFetamiforimagedenoisingintheirpioneeringwork[ 59 ].AsignicantadvantageofTVregularizationisthatthesolutionyieldsarestoredcleanimagewithwell-preservededges.TheTVand1-normtermsin( 3 )leadtoanunderlyingsparsesolutionofAu=f.ThelackofsmoothnessinboththeTVand1-normtermsmakesthesolutionof( 4 )difcult.TocopewiththelackofsmoothnessinJ,weintroduceanauxiliaryvariablevtoobtaintheequivalentconstrainedproblem minJ(v)+H(u)subjecttou=v,u,v2CN.(3)Theequalityconstrainedproblemisconvertedtoanunconstrainedproblemusingaquadraticpenalty: minu,v2CNJ(v)+H(u)+kv)]TJ /F4 11.955 Tf 11.96 0 Td[(uk22,(3)where>0isaparameter.TheadditionalvariablevallowsustotreatthesmoothtermHandthenondifferentiabletermJsomewhatindependently.Startingfromaninitialguessu0,wesolvethepenalizedproblembyrstminimizingovervwithuxed,andthenminimizingoveruwithvxed: vk+1=T(uk),T(u),argminv2CNJ(v)+kv)]TJ /F4 11.955 Tf 11.96 0 Td[(uk22(TV)uk+1=L(vk+1),L(v),argminu2CNH(u)+kv)]TJ /F4 11.955 Tf 11.95 0 Td[(uk22(LS)9>=>;(3)SinceJisconvex,the(TV)subproblemisstronglyconvexinv.Likewise,sinceHisconvex,the(LS)subproblemisstronglyconvexinu.Hence,foranystartingguessu0, 64

PAGE 65

theiterationsequence(vk,uk),k1,existsandisunique.Intheimagingcontext,therstsubproblem,denoted(TV),hasthesameformasTV-waveletbasedimagedenoisingwhichhasbeenextensivelystudiedintheliterature,andsecondsubproblem(LS)isaleastsquaresproblem.Bothsubproblemscanbesolvedquickly.Intheliterature,algorithmsoftheform( 3 )arecalledalternatingproximalminimizationalgorithms.Referencesinclude[ 1 7 10 ].Theiteratesconvergetoasolutionof( 3 ),ifasolutionexists,accordingto[ 10 ,Cor.4.5],forexample.Ingeneral,oneneedstolettendtoinnitytoobtainthesolutionof( 4 ).However,ournumericalexperienceinpartiallyparallelimagereconstructionindicatesthatinthisapplication,asuitableapproximationtothesolutionof( 4 )isgeneratedusingaxed,notverylarge.Ourpaperisorganizedasfollows.InSection 3.2 wegiveanoverviewofTV-basedimagereconstructiontechniques.Section 3.3 studiestheconvergencerateof( 3 ).InSection 3.4 wepresentthealgorithmsthatwehaveusedtosolveeachoftheminimizationproblemsin( 3 )whenJandHaregivenby( 3 ).Section 3.5 givesanoverviewofanemergingmagneticresonanceimagingtechnologyknownaspartiallyparallelimaging(PPI).Finally,Section 3.6 comparesouralgorithmtooperatorsplitting[ 51 ]usingPPIgeneratedimages.Notation.Foradifferentiablefunction,rfdenotesthegradientoff,arowvector.Moregenerally,@J(x)denotesthesubdifferentialsetatx,asetofrowvectors.ForanymatrixM,N(M)isthenullspaceofM.xTdenotestheconjugatetransposeofthevectorxandhx,yi=xTyistheEuclideaninnerproduct.kkpisthep-norm,andkkTVisthediscretetotalvariationsemi-norm. 3.2RelatedWorkInreviewingthepriorworkonTV-basedimagereconstruction,wesimplify( 3 )bytaking=0.Withthissimplication,theimagereconstructionproblemisequivalentto 65

PAGE 66

solvingtheproblem minu2CNkukTV+kAu)]TJ /F4 11.955 Tf 11.95 0 Td[(fk22,(3)wherekkTVisthediscrete(isotropic)TVsemi-normdenedby kukTV,NXi=1kDiuk2,(3)whereDi2R2Nhastwononzeroentriesineachrowcorrespondingtonitedifferenceapproximationstopartialderivativesalongthecoordinateaxes,andNisthenumberofpixelsintheimage.Theearlyworkonalgorithmsfor( 3 )usedgradientdescentmethodswithexplicit[ 59 ]orsemi-implicitschemes[ 43 63 ]inwhichtheTVnormwasreplacedbyasmoothapproximation kukTV,=NXi=1q kDiuk22+.(3)Thechoiceof>0wascrucialtothereconstructionresultsandconvergencespeed.Alargeencouragesfastconvergencerate,butfailstopreservehighqualitydetailssuchasedgesintherestoredimage;asmallbetterpreservesnestructureinthereconstructionattheexpenseofslowconvergence.In[ 64 67 ],amethodisdevelopedbasedonthefollowingreformulationof( 3 ): minu,wNXi=1kwik2+kAu)]TJ /F4 11.955 Tf 11.96 0 Td[(fk22,subjecttowi=Diu,i=1,,N.(3)Thelinearconstraintwastreatedwithaquadraticpenalty minu,wNXi=1kwik2+kDu)]TJ /F4 11.955 Tf 11.96 0 Td[(wk22+kAu)]TJ /F4 11.955 Tf 11.95 0 Td[(fk22,(3)wherew2C2NandDisobtainedbystackingtheDimatrices.Foranyxed,( 3 )canbesolvedbyalternatingminimizations.IfbothD>DandA>AcanbediagonalizedbytheFouriermatrix,astheywouldifAiseithertheidentitymatrixorablurringmatrixwithperiodicboundaryconditions,theneachminimizationinvolvesshrinkageand 66

PAGE 67

afastFouriertransform(FFT).Acontinuationmethodisusedtodealwiththeslowconvergencerateassociatedwithalargevaluefor.Themethod,however,isnotapplicabletomoregeneralA.In[ 34 ]GoldsteinandOsherdevelopasplitBregmanmethodfor( 3 ).Theresultingalgorithmhassimilarcomputationalcomplexitytothealgorithmin[ 64 ];theconvergenceisfastandtheconstraintsareexactlysatised.LaterthesplitBregmanmethodwasshowntobeequivalenttothealternatingdirectionmethodofmultipliers(ADMM)[ 13 28 32 33 ]appliedtotheaugmentedLagrangian L(w,u,p),NXi=1kwik2+kAu)]TJ /F4 11.955 Tf 11.96 0 Td[(fk22+hp,Du)]TJ /F4 11.955 Tf 11.95 0 Td[(wi+kDu)]TJ /F4 11.955 Tf 11.96 0 Td[(wk22.(3)Nonetheless,thealgorithmsin[ 34 64 67 ]benetfromthespecialstructureofA,andtheyloseefciencyifATAcannotbediagonalizedbyfasttransforms.TotreatamoregeneralA,theBregmanoperatorsplitting(BOS)method[ 75 ]replaceskAu)]TJ /F4 11.955 Tf 12.4 0 Td[(fk22byaproximal-liketermku)]TJ /F3 11.955 Tf 12.6 0 Td[((uk)]TJ /F5 11.955 Tf 12.61 0 Td[(A>(Auk)]TJ /F4 11.955 Tf 12.61 0 Td[(f))k2=forsome>0.BOSisaninexactUzawamethodthatdependsonthechoiceof.ItisgenerallylessefcientthansplitBregman.Therearealsoseveralmethodsdevelopedtosolvetheassociateddualorprimal-dualformulationsof( 3 )basedonthedualformulationoftheTVnorm: kukTV=maxp2Xhp,Dui,whereX=fp=(p1,...,pN):pi2C2,1iNg(3)Consequently,( 3 )canbewrittenasaminimaxproblem minu2CNmaxp2Xhp,Dui+kAu)]TJ /F4 11.955 Tf 11.96 0 Td[(fk22.(3)In[ 20 ],Chanetal.proposedtosolvetheprimal-dualEuler-LagrangeequationsusingNewton'smethod.Thisleadstoaquadraticconvergencerateandhighlyaccuratesolutions;however,thecostperiterationismuchhighersincethemethodexplicitlyusessecond-orderinformationandtheinversionofaHessianmatrixisrequired.In[ 19 ], 67

PAGE 68

ChambolleusedthedualformulationoftheTVdenoisingproblem( 3 )withA=I,andprovidedanefcientsemi-implicitgradientdescentalgorithmforthedual.However,themethoddoesnotnaturallyextendtothecasewithmoregeneralA.Recently,ZhuandChan[ 78 ]proposedaprimal-dualhybridgradient(PDHG)method.PDHGalternatelyupdatestheprimalanddualvariablesuandp.NumericalresultsshowthatPDHGoutperformsmethodsin[ 19 34 ]fordenoisinganddeblurringproblems,butitsefciencyagainreliesonthefactthatATAcanbediagonalizedbyfasttransforms.Later,severalvariationsofPDHG,referredtoasprojectedgradientdescentalgorithms,wereappliedtothedualformulationofimagedenoisingproblemin[ 79 ]tomakethemethodmoreefcient.Furtherenhancementsinvolvedifferentstep-lengthrulesandline-searchstrategies,includingtechniquesbasedontheBarzilai-Borweinmethod[ 9 ].Anotherapproachthatcanbeappliedto( 4 )intheimagingcontext( 3 )withageneralAistheoperatorsplitting(OS)method.In[ 51 ]theOSideaof[ 47 ]isappliedtoimagereconstructionincompressedmagneticresonanceimaging.Theschemeisbasedontherst-orderoptimalityconditionatalocalminimizeru:02@J(u)+AT(Au)]TJ /F4 11.955 Tf 11.96 0 Td[(f).Thisisrewrittenintheform02@J(u)+1 (u)]TJ /F4 11.955 Tf 11.96 0 Td[(s),s=u)]TJ /F5 11.955 Tf 11.95 0 Td[(AT(Au)]TJ /F4 11.955 Tf 11.96 0 Td[(f).Theiterativeschemeissk=uk)]TJ /F5 11.955 Tf 11.95 0 Td[(AT(Auk)]TJ /F4 11.955 Tf 11.96 0 Td[(f),uk+1=argminuJ(u)+1 2jju)]TJ /F4 11.955 Tf 11.95 0 Td[(skjj22.Thecomputationofuk+1,givensk,isaTV-denoisingproblem.IfthisproblemissolvedusingasplitBregmanmethod[ 34 ],thenthisisequivalenttotheBregmanoperatorsplittingmethod[ 75 ].In[ 51 ]axedpointiterationisappliedtosolvetheusubproblem; 68

PAGE 69

theiterativeschemecombinestheOSideawithconjugateduality.Thecomparisonofourmethodwiththealgorithmof[ 51 ]isgiveninSection 3.6 3.3ConvergenceAnalysisInthissection,weexaminetheconvergencerateofthealternatingproximalminimizationscheme( 3 ).SinceHisconvex,thereexistsaconstant0suchthatthefollowingmonotonicityconditionholdsforalluandv2Cn: (rH(u))-222(rH(v))(u)]TJ /F4 11.955 Tf 11.96 0 Td[(v)ku)]TJ /F4 11.955 Tf 11.96 0 Td[(vk22(3)Here,rHdenotesthegradient,whichisarowvectorthroughoutthispaper.If>0,thenHisstronglyconvex.AsshownbelowinCorollary 1 ,strongconvexityofHandconvexityofJimplythattheobjectivefunctioninthepenalizedproblem( 3 )isstronglyconvex,whichensurestheexistenceofauniqueminimizer. Theorem3.1. If( 3 )hasminimizersvandu,thenforeachkwehave kvk+1)]TJ /F4 11.955 Tf 11.96 0 Td[(vk22 2+kvk)]TJ /F4 11.955 Tf 11.95 0 Td[(vk2andkuk+1)]TJ /F4 11.955 Tf 11.96 0 Td[(uk22 2+kuk)]TJ /F4 11.955 Tf 11.95 0 Td[(uk2.(3) Proof. Itiswell-knownthattheoperatorsTandLin( 3 )arenonexpansiverelativetotheEuclideannorm.Thatis,foralluandv,wehavekT(v))-222(T(u)k2kv)]TJ /F4 11.955 Tf 11.95 0 Td[(uk2andkL(v))-222(L(u)k2kv)]TJ /F4 11.955 Tf 11.95 0 Td[(uk2.Thisfollowsfromtherst-orderoptimalityconditionscharacterizingtheminimizersin( 3 ).Forexample,ifvi=T(ui)fori=1,2,then2(ui)]TJ /F4 11.955 Tf 11.99 0 Td[(vi)T2@J(vi),where@denotesthesubdifferential.BytheconvexityofJ,itfollowsthat J(v2)J(v1)+2(u1)]TJ /F4 11.955 Tf 11.96 0 Td[(v1)T(v2)]TJ /F4 11.955 Tf 11.96 0 Td[(v1).(3)Likewise,interchangingv1andv2gives J(v1)J(v2)+2(u2)]TJ /F4 11.955 Tf 11.96 0 Td[(v2)T(v1)]TJ /F4 11.955 Tf 11.96 0 Td[(v2).(3) 69

PAGE 70

Weadd( 3 )and( 3 )toobtain kv2)]TJ /F4 11.955 Tf 11.96 0 Td[(v1k2(u2)]TJ /F4 11.955 Tf 11.96 0 Td[(u1)T(v2)]TJ /F4 11.955 Tf 11.96 0 Td[(v1)ku2)]TJ /F4 11.955 Tf 11.96 0 Td[(u1k2kv2)]TJ /F4 11.955 Tf 11.96 0 Td[(v1k2.(3)Hence,kv2)]TJ /F4 11.955 Tf 12.53 0 Td[(v1k2=kT(u2))-270(T(u1)k2ku2)]TJ /F4 11.955 Tf 12.53 0 Td[(u1k2,whichyieldsthenonexpansiveproperty.Sincevanduachievetheminimumin( 3 ),wehavev=T(u).Subtractingthisidentityfromtheequationvk+1=T(uk)andutilizingthenonexpansivepropertygives kvk+1)]TJ /F4 11.955 Tf 11.96 0 Td[(vk2kT(uk))-222(T(u)k2kuk)]TJ /F4 11.955 Tf 11.96 0 Td[(uk2.(3)Therst-orderoptimalityconditionsforukanduarerH(uk))]TJ /F3 11.955 Tf 11.96 0 Td[(2(vk)]TJ /F4 11.955 Tf 11.96 0 Td[(uk)T=0,rH(u))]TJ /F3 11.955 Tf 11.95 0 Td[(2(v)]TJ /F4 11.955 Tf 11.95 0 Td[(u)T=0.Wesubtractthesecondequationfromtherstandmultiplyby(uk)]TJ /F4 11.955 Tf 11.96 0 Td[(u)toobtain (rH(uk))-221(rH(u))(uk)]TJ /F4 11.955 Tf 11.96 0 Td[(u)+2kuk)]TJ /F4 11.955 Tf 11.95 0 Td[(uk22=2(vk)]TJ /F4 11.955 Tf 11.96 0 Td[(v)T(uk)]TJ /F4 11.955 Tf 11.96 0 Td[(u) (3) 2kvk)]TJ /F4 11.955 Tf 11.96 0 Td[(vk2kuk)]TJ /F4 11.955 Tf 11.96 0 Td[(uk2.Utilizingthemonotonicitycondition( 3 )ontheleftsideof( 3 )gives(+2)kuk)]TJ /F4 11.955 Tf 11.96 0 Td[(uk222kvk)]TJ /F4 11.955 Tf 11.96 0 Td[(vk2kuk)]TJ /F4 11.955 Tf 11.95 0 Td[(uk2,whichyields kuk)]TJ /F4 11.955 Tf 11.95 0 Td[(uk22 +2kvk)]TJ /F4 11.955 Tf 11.95 0 Td[(vk2.(3)Combiningthiswith( 3 )giveskvk+1)]TJ /F4 11.955 Tf 11.95 0 Td[(vk22 +2kvk)]TJ /F4 11.955 Tf 11.95 0 Td[(vk2,therstinequalityin( 3 ).Combining( 3 ),withkreplacedbyk+1,andthenonexpansiveproperty( 3 )givesthesecondinequalityin( 3 ). 70

PAGE 71

Corollary1. If>0,thentheiteratesgeneratedby( 3 )convergelinearlytotheuniqueminimizerof( 3 ). Proof. Werstobservethatwhen>0,theobjectivefunctionin( 3 )stronglyconvex.LetF(u,v)=H(u)+kv)]TJ /F4 11.955 Tf 12.44 0 Td[(uk22bethepartoftheobjectivewhichexcludesJ.Bytheconvexityinequality( 3 ),wehave (rF(u1,v1))-222(rF(u2,v2))264uv375=(rH(u1))-222(rH(u2))(u1)]TJ /F4 11.955 Tf 11.95 0 Td[(u2)+2ku)]TJ /F5 11.955 Tf 11.96 0 Td[(vk22kuk22+2ku)]TJ /F5 11.955 Tf 11.95 0 Td[(vk22, (3) whereu=u1)]TJ /F4 11.955 Tf 11.27 0 Td[(u2andv=v1)]TJ /F4 11.955 Tf 11.27 0 Td[(v2.Thematrixcorrespondingtothequadraticin( 3 )is2264+=2)]TJ /F5 11.955 Tf 9.3 0 Td[()]TJ /F5 11.955 Tf 9.3 0 Td[(375.Sincetheeigenvaluesofthismatrixarestrictlypositive,Fisstronglyconvex.Theobjectivefunctionin( 3 )isthesumJ+FofaconvexfunctionJandastronglyconvexfunctionF.Hence,itisstronglyconvexandthereexistsauniqueminimizer(u,v).ByTheorem 3.1 ,theiteratesgeneratedby( 3 )convergeto(u,v)linearly. Inthecase=0,Theorem 3.1 onlyyields kvk+1)]TJ /F4 11.955 Tf 11.95 0 Td[(vk2kvk)]TJ /F4 11.955 Tf 11.96 0 Td[(vk2andkuk+1)]TJ /F4 11.955 Tf 11.95 0 Td[(uk2kuk)]TJ /F4 11.955 Tf 11.95 0 Td[(uk2,(3)whichdoesnotimplyconvergence.Ontheotherhand,bythetheoryforthealternatingproximalminimizationalgorithm,weknowthattheiteratesdoconverge.Wenowobservethattheinequalitiesin( 3 )arestrictexceptwhenconvergenceisachievedinanitenumberofsteps.Thisresultisbasedonthefollowingproperty. Lemma1. IfP:Cn!Cnsatises kP(u))-222(P(v)k22hP(u))-222(P(v),u)]TJ /F4 11.955 Tf 11.96 0 Td[(vi(3) 71

PAGE 72

foralluandv2Cn,then kP(u))-222(P(v)k2ku)]TJ /F4 11.955 Tf 11.95 0 Td[(vk2(3)foralluandv2CnwithequalityonlyifP(u))-222(P(v)=u)]TJ /F4 11.955 Tf 11.96 0 Td[(v.Operatorssatisfying( 3 )arecalledrmlynonexpansive.ThefactthattheproximalmapsTorLarermlynonexpansiveisimpliedby( 3 ). Proof. Theinequality( 3 )isaconsequenceoftheSchwarzinequalityappliedto( 3 ).Moreover,by( 3 )wehave k(u)]TJ /F4 11.955 Tf 11.95 0 Td[(v))]TJ /F3 11.955 Tf 11.95 0 Td[((P(u))-222(P(v))k22=ku)]TJ /F4 11.955 Tf 11.95 0 Td[(vk22)]TJ /F3 11.955 Tf 11.96 0 Td[(2hP(u))-222(P(v),u)]TJ /F4 11.955 Tf 11.96 0 Td[(vi+kP(u))-222(P(v)k22ku)]TJ /F4 11.955 Tf 11.95 0 Td[(vk22)-222(kP(u))-222(P(v)k22. (3) If( 3 )isanequality,thentherightsideof( 3 )vanishes,whichimpliesthattheleftsidevanishes:(u)]TJ /F4 11.955 Tf 11.96 0 Td[(v))]TJ /F3 11.955 Tf 11.95 0 Td[((P(u))-221(P(v))=0. Theorem3.2. Supposethatuandvareoptimalin( 3 ).Ifforsomek,theiteratesofthealternatingproximalminimizationalgorithm( 3 )satisfykuk+1)]TJ /F4 11.955 Tf 12.16 0 Td[(uk2=kuk)]TJ /F4 11.955 Tf 12.16 0 Td[(uk2,thenuj=ukandvj+1=vk+1forallj>k.Ifkvk+1)]TJ /F4 11.955 Tf 12.11 0 Td[(vk2=kvk)]TJ /F4 11.955 Tf 12.11 0 Td[(vk2forsomek,thenvj=vkanduj=ukforallj>k. Proof. Supposethatkuk+1)]TJ /F4 11.955 Tf 11.97 0 Td[(uk2=kuk)]TJ /F4 11.955 Tf 11.98 0 Td[(uk2.Sincevanduareoptimalin( 3 ),wehave (LT)(u)=L(T(u))=L(v)=u.(3) 72

PAGE 73

By( 3 ),itfollowsthatuk+1=(LT)(uk).Hence,theequalitykuk+1)]TJ /F4 11.955 Tf 12.07 0 Td[(uk2=kuk)]TJ /F4 11.955 Tf 12.07 0 Td[(uk2coupledwiththenonexpansivepropertiesofLandTyield kuk)]TJ /F4 11.955 Tf 11.96 0 Td[(uk2=k(LT)(uk))]TJ /F3 11.955 Tf 11.96 0 Td[((LT)(u)k2=kL(T(uk)))-222(L(T(u))k2kT(uk))-222(T(u)k2kuk)]TJ /F4 11.955 Tf 11.95 0 Td[(uk2. (3) Sincetherightandleftsidesof( 3 )areequal,alltheinequalitiesin( 3 )areequalities.TheequalitykT(uk))-222(T(u)k2=kuk)]TJ /F4 11.955 Tf 11.96 0 Td[(uk2andLemma 1 implythat T(uk))-221(T(u)=uk)]TJ /F4 11.955 Tf 11.96 0 Td[(u.(3)TheequalitykL(T(uk)))-221(L(T(u))k2=kT(uk))-221(T(u)k2andLemma 1 implythat (LT)(uk))]TJ /F3 11.955 Tf 11.95 0 Td[((LT)(u)=L(T(uk)))-222(L(T(u))=T(uk))-222(T(u).(3)Together,( 3 )and( 3 )yield (LT)(uk))]TJ /F3 11.955 Tf 11.95 0 Td[((LT)(u)=uk)]TJ /F4 11.955 Tf 11.96 0 Td[(u.(3)Wecombine( 3 )and( 3 )toobtainuk+1=(LT)(uk)=uk.Hence,ukisaxedpointof(LT)anduj=ukforallj>k.Sincevj+1=T(uj),weconcludethatvj+1=vk+1forallj>k.Theequalitykvk+1)]TJ /F4 11.955 Tf 12.05 0 Td[(vk2=kvk)]TJ /F4 11.955 Tf 12.05 0 Td[(vk2istreatedinthesamewayexceptthatLandTareinterchanged. Bytheconvergencetheoryforthealternatingproximalminimizationalgorithm,weknowthattheiteratesconvergetoasolution(u,v)of( 3 )providedasolutionexists.Theorem 3.2 impliesthatkuk+1)]TJ /F4 11.955 Tf 11.96 0 Td[(uk2=kuk)]TJ /F4 11.955 Tf 11.95 0 Td[(uk2<1 73

PAGE 74

exceptwhenuk=u.Likewisekvk+1)]TJ /F4 11.955 Tf 11.96 0 Td[(vk2=kvk)]TJ /F4 11.955 Tf 11.95 0 Td[(vk2<1exceptwhenvk=v.Thisimpliesatleastsublinearconvergenceofthealternatingproximalminimizationalgorithm( 3 ).Foranyxed,thesolutionof( 3 )generatesanapproximationtoasolutionof( 4 ).Letk,k0,denoteanincreasingsequenceofvaluesforthepenaltyparametertendingtoinnity,andlet(Uk,Vk)denoteassociatedsolutionsof( 3 ),assumingtheyexist.Bythetheorydescribingtheconvergenceofthepenaltyscheme(see[ 53 ,Thm.17.1]),convergentsubsequencesoftheiteratesapproachasolutionof( 4 ).Wenowshowinthecontext( 3 )ofimagereconstructionthattheiterates(Uk,Vk)arebounded. Theorem3.3. SupposethatJandHaregivenby( 3 ).If0,>0,andN(D)\N(A)=0,whereNdenotesnullspace,thenforeach0>0,thereexistsacompactsetKwhichcontainsthesolutionsof( 3 )forall0.Moreover,astendstoinnity,anyconvergentsubsequenceoftheiteratesapproachesasolutionofeither( 4 )ortheequivalentproblem( 3 ). Proof. Inthespecialcase( 3 ),J(0)=0andH(0)=kfk22.Let=kfk22bethevalueoftheobjectivefunctionvaluein( 3 )correspondingtou=v=0.Foranychoiceof,theoptimalobjectivefunctionvaluein( 3 )isboundedby.Hence,foranychoiceof,whenminimizingtheobjectivefunctionin( 3 ),weshouldrestrictourattentiontothoseuandvsatisfying J(v)+H(u)+kv)]TJ /F4 11.955 Tf 11.96 0 Td[(uk22.(3) 74

PAGE 75

SinceJ(v)=kvkTV+kvk10andH(u)=kAu)]TJ /F4 11.955 Tf 11.95 0 Td[(fk220,itfollowsfrom( 3 )that kv)]TJ /F4 11.955 Tf 11.95 0 Td[(uk2p =, (3) kvkTV, (3) kAu)]TJ /F4 11.955 Tf 11.95 0 Td[(fk2p =. (3) Decomposeu=un+upwhereun2N(A)andupisorthogonaltoN(A).By( 3 ),( 3 ),and( 3 ),wehave kvkTV=NXi=1kDivk2kDvk2kDuk2)-222(kD(v)]TJ /F4 11.955 Tf 11.95 0 Td[(u)k2kDunk2)-222(kDupk2)-222(kDk2kv)]TJ /F4 11.955 Tf 11.95 0 Td[(uk2kDunk2)-222(kDupk2)-222(kDk2p =. (3) SinceN(D)\N(A)=0,thereexists1>0suchthatkDuk21kuk2forallu2N(A).Hence,by( 3 ) kunk2+kDupk2+kDk2p ==1.(3)Similarly,thereexists2>0suchthatkAupk22kupk2.Hence,by( 3 ),wehave 2kupk2kAuk2kfk2+kAu)]TJ /F4 11.955 Tf 11.95 0 Td[(fk2kfk2+p =.(3)Combine( 3 )and( 3 )todeducethatu=un+upliesinacompactset.By( 3 ),wehavekvk2kuk2+p =, 75

PAGE 76

whichyieldsaboundforkvk2.Asincreases,thelevelsetof( 3 )correspondingtotheobjectivefunctionvaluecanonlyshrink.Hence,thislevelsetisboundedforany0.Letk,k=0,1,...,denoteanincreasingsequenceofvaluesforthepenaltytendingtoinnity,andlet(Uk,Vk)denoteassociatedsolutionsof( 3 ).By[ 53 ,Thm.17.1],everyconvergentsubsequenceoftheminimizers(Uk,Vk)approachesasolutionof( 3 ). Remark.IfHisstronglyconvex,then( 4 )hasauniquesolution;hence,anysequenceofsolutionsto( 3 )approachestheuniquesolutionof( 4 )astendstoinnity. 3.4AlgorithmsfortheTVandLSSubproblemsWenowprovideimplementationsforthe(TV)and(LS)subproblemsofthealternatingproximalminimizationalgorithm( 3 )intheimagingcontext( 3 ).Oneofthereasonsthatthesplitting( 3 )workedwellwasthateachofthesubproblemscouldbesolvedquickly.Forxedu,subproblem(TV)isaTV-waveletimagedenoisingproblem.Asdiscussedearlier,therearemanyfastalgorithmsforthisproblemthattakeadvantageofthesimplicityofthekv)]TJ /F4 11.955 Tf 11.01 0 Td[(uk22term.Recentworkincludesthedualapproachin[ 19 ],variablesplittingandcontinuation[ 64 67 ],splitBregman[ 34 ],primal-dualhybridgradient[ 78 79 ].InthenumericalexperimentsofSection 3.6 ,weusedasplitBregmanmethodwhichisamongthefastestmethodsforTV-waveletimagedenoising.WenowexplainindetailthesplitBregmanschemethatweusefortheTVsubproblemin( 3 ).Basedontherepresentation( 3 )fortheTVnorm,theTVsubproblemhastheformminvkvk1+kv)]TJ /F4 11.955 Tf 11.95 0 Td[(uk22+NXi=1kDivk2.Introducingadditionalvariableswandz,wherewi=Divandv=z,thisisrewrittenminv,w,zkzk1+kv)]TJ /F4 11.955 Tf 11.96 0 Td[(uk22+NXi=1kwik2subjecttoDv=w,v=z. 76

PAGE 77

WeapplythesplitBregmanschemegivenasAlgorithm(A1)in[ 76 ],butwithouttheproximalterms.LetbkandckdenoteapproximationstomultipliersfortheconstraintsDv=wandv=z,andlet1and2denotecorrespondingconstraintpenalties.Theiterationhasthefollowingform: Algorithm3SplitBregman[ 76 ]forTVSubproblem vk+1=argminv2CNkv)]TJ /F4 11.955 Tf 11.96 0 Td[(uk22+(Dv)]TJ /F4 11.955 Tf 11.95 0 Td[(wk)Tbk+(v)]TJ /F4 11.955 Tf 11.95 0 Td[(zk)Tck+1kDv)]TJ /F4 11.955 Tf 11.95 0 Td[(wkk22+2kv)]TJ /F4 11.955 Tf 11.96 0 Td[(zkk22 (3) wk+1i=argminwi2C2kwik2+(Divk+1)]TJ /F4 11.955 Tf 11.95 0 Td[(wi)Tbki+1kDivk+1)]TJ /F4 11.955 Tf 11.95 0 Td[(wik22,1iN (3) zk+1=argminz2CNkzk1+(vk+1)]TJ /F4 11.955 Tf 11.96 0 Td[(z)Tck+2kvk+1)]TJ /F4 11.955 Tf 11.95 0 Td[(zk22 (3) bk+1=bk+21(Dvk+1)]TJ /F4 11.955 Tf 11.95 0 Td[(wk+1) (3) ck+1=ck+22(vk+1)]TJ /F4 11.955 Tf 11.96 0 Td[(zk+1) (3) InAlgorithm1,( 3 )isawell-conditionedleastsquaresproblem;itcanbesolvedbyaniterativemethodsuchasGauss-SeidelorbyanFFTiftheassociatedimagesatisesperiodicboundaryconditions[ 64 ].Forthesteps( 3 )and( 3 ),thereareclosedformsolutions[ 34 64 67 ],the2Dandcomponentwiseshrinkageoperators.Thenalsteps( 3 )and( 3 )aretherst-ordermultiplierupdates.Basedontheresultsgivenin[ 34 ],Algorithm1isexpectedtobeveryefcient.TheLSsubproblemin( 3 )isaleast-squaresprobleminu.Thiscouldbesolvedbyaconjugategradientmethod,however,wehaveobtainedcomparableorbetterperformanceusingtheBarzilai-Borwein[ 9 ]method(BB),whichhandlesill-conditioningmuchbetterthangradientmethodswithaCauchystep[ 4 ].TheLSsubproblemhastheform minukAu)]TJ /F4 11.955 Tf 11.95 0 Td[(fk22+kv)]TJ /F4 11.955 Tf 11.95 0 Td[(uk2.(3)InthestandardimplementationoftheBBmethod,theHessianoftheobjectivefunctionisapproximatedbyamultipleoftheidentitymatrix.FortheLSproblem,however,theHessianofkv)]TJ /F4 11.955 Tf 12.32 0 Td[(uk2withrespecttoualreadyamultipleoftheidentity.Hence,weonly 77

PAGE 78

approximatetheHessianofkAu)]TJ /F4 11.955 Tf 11.68 0 Td[(fk22byamultipleoftheidentity.Moreprecisely,ifukisthecurrentBBiterate,thenweemploytheapproximation kAu)]TJ /F4 11.955 Tf 11.96 0 Td[(fk22kAuk)]TJ /F4 11.955 Tf 11.95 0 Td[(fk22+2(Auk)]TJ /F4 11.955 Tf 11.95 0 Td[(f)TA(u)]TJ /F4 11.955 Tf 11.96 0 Td[(uk)+kku)]TJ /F4 11.955 Tf 11.95 0 Td[(ukk22,(3)wherek=kA(uk)]TJ /F4 11.955 Tf 11.96 0 Td[(uk)]TJ /F6 7.97 Tf 6.58 0 Td[(1)k22=kuk)]TJ /F4 11.955 Tf 11.96 0 Td[(uk)]TJ /F6 7.97 Tf 6.58 0 Td[(1k22.SincethekAuk)]TJ /F4 11.955 Tf 12.46 0 Td[(fk22termin( 3 )doesnotdependonu,theBBmethodfortheLSsubproblemisasfollows: Algorithm4BBmethod[ 9 ]forLSSubproblem uk+1=argminu2CN)]TJ /F3 11.955 Tf 5.48 -9.68 Td[(2(Auk)]TJ /F4 11.955 Tf 11.96 0 Td[(f)TA(u)]TJ /F4 11.955 Tf 11.96 0 Td[(uk)+kku)]TJ /F4 11.955 Tf 11.96 0 Td[(ukk22+kv)]TJ /F4 11.955 Tf 11.95 0 Td[(uk2.(3) Theiteration( 3 )convergeslinearlytoasolutionof( 3 )by[ 24 25 31 ].EachiterationinvolvesmultiplicationbyAandAT.Intheapplicationtopartiallyparallelimagingdevelopedinthenextsection,AisrepresentedasaproductMFSjwherethematrixMistheidentitymatrixwithsomerowsremoved,FisaFouriertransform,andSjisadiagonalmatrix.ThetimetomultiplybyMorSjisboundedbyaconstanttimesN,whiletheFouriertransformcanbeperformedintimeproportionaltoNlog(N).Hence,eachiterationofAlgorithm 4 canbeperformedquicklyinourtargetapplication. 3.5PartiallyParallelImaging(PPI)Magneticresonance(MR)imagingisamedicalimagingtechniquecommonlyusedinradiologytovisualizetheinternalstructureandfunctionofthebodybynon-invasiveandnon-ionizingmeans.Itprovidesbettercontrastbetweenthedifferentsofttissuesthanmostothermodalities.MRimagesareobtainedthroughaninversionofFourierdataacquiredbythereceivercoil(s).Magneticresonanceimagesareobtainedbyplacinganobjectinastrongmagneticeldandthenturningonandoffaradiofrequencyelectromagneticeld.Differentbodypartsproducedifferentsignalswhicharedetectedbyreceivers.Theresultingdatais 78

PAGE 79

Figure3-1. PPImask theninvertedtoobtainanimageofthescannedobject.Toimprovethequalityoftheimage,anumberofreceiversarepointedatthescannedobjectfromdifferentdirections,anddataiscollectedinparallel.Toacceleratetheimagingprocess,onlypartoftheFouriercomponentsarerecordedbythereceiver.ThistechnologybasedoncollectinginparallelpartialFourierdatafromdifferentcoilarraysiscalledpartiallyparallelimagingorPPI.TheundersamplingpatternsoftheFouriercoefcientsareoftendescribedbyamask.Figure 3-1 showsasimulatedradialmaskfora2Dimage.ThewhitepixelscorrespondtotheFouriercomponentwhicharemeasured.ThewhiteregioninthecenterofthemaskindicatesthatthelowfrequencyFouriercomponentsareallmeasured,whilethewhiteraysinthesurroundingdarkerregionshowsthespacingbetweenthehigherfrequencyFouriercomponentsthataremeasured.Partialdataacquisitionincreasesthespacingbetweenregularsubsequentread-outlines,therebyreducingscantime,however,thisreductioninthenumberofrecordedFouriercomponentsleadstoaliasingartifactsinimages.Therearetwogeneralapproachesforremovingthealiasingartifactsandreconstructinghighqualityimages,imagedomain-basedmethodsandk-spacebasedmethods.Thek-spacebased 79

PAGE 80

Figure3-2. Sensitivitymapsforan8-channelcoil methodsusecoilsensitivityvariationstoreconstructthemissingk-spacedata,andthenapplytheFouriertransformtotheoriginalandreconstructeddatatoobtaintheunaliasedimage[ 6 35 48 ].Inthispaper,weemployimagedomainmethodsandcoilsensitivitymapstoreconstructtheunderlyingimage[ 15 29 54 57 60 70 73 ].SensitivityEncoding(SENSE)isthemostcommonimagedomain-basedparallelimagingmethod.Itisbasedonthefollowingequationwhichrelatesthepartialk-spacedatafj,acquiredbythej-thchannel,tothesensitivitymapSjandthemaskM: MF(Sju)=fj(3)HereistheHadamard(orcomponentwise)productbetweentwovectors,fjisthevectorofmeasuredFouriercoefcientsatreceiverj,MisthemaskwhichisobtainedbyextractingfromtheidentitythoserowscorrespondingtothemeasuredFouriercomponents,FistheFouriertransform,Sj2CNisthesensitivitymapforreceiverj,andu2CNistheunderlyingimagegottenbystackingallcolumnsoftheimagetoformaonedimensionalvector.ThesensitivitymapisanestimateoftheimpactofapixelintheimageonthemeasuredFouriercoefcients.Pixelsclosesttoareceivermayhavemoreimpactonthesignalthanpixelsfarawayfromthereceiver.Anexampleofthesensitivitymapforan8-channelcoilappearsinFigure 3-2 80

PAGE 81

Basedon( 3 ),thereconstructionoftheimageucouldbeaccomplishedbysolvingtheleastsquaresproblem minu2CNJXj=1kMF(Sju))]TJ /F4 11.955 Tf 11.95 0 Td[(fjk22,(3)whereJisthenumberofchannels.However,( 3 )oftendoesnothaveauniquesolutionandtheminimizationproblemcanbeill-conditioned.Toalleviatetheeffectoftheill-conditioning,theSENSEmodel( 3 )hasbeenimprovedrecentlybyincorporatingregularizationtermsintotheenergyfunctionaltotakeadvantageoftheunderlyingsparsityofMRimagesinthenitedifferencedomainandwavelettransformdomain[ 50 ].Theimagesarerecoveredbysolvinganoptimizationproblemoftheform minu2CNkukTV+kuk1+JXj=1kMF(Sju))]TJ /F4 11.955 Tf 11.96 0 Td[(fjk22.(3)Thersttwotermsin( 3 )correspondtoJin( 3 )whilethelasttermcorrespondstoH(u)=kAu)]TJ /F4 11.955 Tf 11.96 0 Td[(fk22. 3.6NumericalExperiments 3.6.1DataAcquisitionandExperimentalSetupInthissectionwegiveresultsforthreePPIreconstructionsbasedonthealgorithm( 3 ).Wecompareperformancetothatoftheoperatorsplitting(OS)in[ 51 ].Allk-spacedatawerefullyacquiredwith8-channelheadcoilasillustratedinFigure 3-3 .Byfullacquisitionwemeanthateachreceivercoilobtainsthecompletek-spacedataandhenceahighresolutionimage.Oneofthedatasets,denotedSB512,wasasetofsagittalCartesianbrainimagesacquiredona3TGEsystem(GEHealthcare,Waukesha,Wisconsin,USA).ThedataacquisitionparameterswereFOV220mm2,size5125128,TR3060ms,TE126ms,slicethickness5mm,ipangle90,andphaseencodingdirectionwasanterior-posterior.Theseconddatasetwasaradialbraindatasetacquiredona1.5TSiemensSymphonysystem(SiemensMedicalSolutions,Erlangen,Germany).TheacquisitionparameterswereFOV220mm2,size2565128, 81

PAGE 82

slicethickness5mm,TR53.5ms,TE3.4ms,andipangle75.ACartesiandatasetwithfullk-space,denotedAB512,withsize5125128,wasgeneratedbyGRAPPAoperatorgriding[ 61 ]whichcanshiftnon-CartesianradialdatatoCartesiangrids.Thethirddatasetofsize2562568,denotedAB256,hasfullk-spacesimulatedbythecompletecentralk-spaceofAB512.WesimulatedaPPIscanbyundersamplingtheactualdatausingaradialmasksimilartothatshowninFigure 3-1 .Weused88radiallinescorrespondingtosamplingratiosof33.5%forAB256andof16.8%forSB512andAB512. Figure3-3. AnillustrationofaheadPPIsystemwitheightreceivercoils. Thesensitivitymapswereestimatedusingcenterk-spacedataofsize3232,whichisasubsetofthepartiallysampleddata.TheestimatedsensitivitymapsforSB512areshowninFigure 3-2 .Inallexperiments,thesensitivitymapswereobtainedinthesameway;themapswerexedduringthereconstructionprocess.AlgorithmswereimplementedinMATLABr,VersionR2009b.AlltheexperimentswereperformedonaLenovolaptopwithIntelrDualCore2Duo2.53GHzprocessorsandaWindowsTMoperatingsystem.Only1corewasusedinthecomputations. 3.6.2ComparisonAlgorithmManyofthealgorithmsinSection 3.2 arenotveryeffectiveforPPIimagingduetothecomplicatedstructureofA.Forcomparison,wechosetheoperatorsplitting(OS) 82

PAGE 83

Table3-1. TestedDataandModelParameters DataImageAbbrev.SizeSampleRatio(,) 1AxialBrainAB256256256833.5%(.25,500)2SagittalBrainSB512512512816.8%(.25,500)3AxialBrainAB512512512816.8%(.25,500) schemefrom[ 51 ],whichisrelativelyefcientandonlyrequiresthecomputationsofAandA>ineachiteration.TheOSschemeof[ 51 ],withaminorcorrection,isasfollows:sk+1=uk)]TJ /F3 11.955 Tf 11.96 0 Td[((1=))]TJ /F4 11.955 Tf 5.48 -9.69 Td[(D>wk+A>(Auk)]TJ /F4 11.955 Tf 11.96 0 Td[(f),tk+1i=wki+2Diuk,1iN,uk+1=)]TJ /F3 11.955 Tf 5.48 -9.68 Td[(sign(sk+1)maxfjsk+1j)]TJ /F5 11.955 Tf 17.94 0 Td[(1=,0g,wk+1i=min(1,ktk+1ik)tk+1i=ktk+1ik21iN,wheresk2CN,tkiandwki2C2,i=1,,N,wk2C2Nisformedbystackingthetwocolumnsofmatrix(wk1,,wkN)>,andthemaxandsignoperationsinthecomputationofuk+1arecomponentwiseoperationscorrespondingtoshrinkage.ThemaincomputationalcostperiterationintheOSschemecorrespondstothefollowingoperations:a2Dshrinkageduringthecomputationofwk+1,aprojectionduringthecomputationofuk+1,twowavelettransformsduringthecomputationofsk+1anduk+1,andtwo(inverse)Fouriertransformsduringthecomputationofsk+1.In[ 51 ]itisshownthatfor1,2>0incertainranges,theOSschemeconvergestoaxedpointwhichisalsoasolutionof( 4 ).Theiterationswerestoppedwheneitherthefollowingconditionsweresatised: kuk+1)]TJ /F4 11.955 Tf 11.95 0 Td[(ukk2=maxf1,kukk2g<1or(fk)]TJ /F4 11.955 Tf 11.95 0 Td[(fk+1)=maxf1,fkg<2p c=t,(3)wherefkistheobjectivevalueof( 4 )atuk,candtarethecurrentandtargetvaluesofrespectivelyand1and2areprescribedstoppingtolerances. 3.6.3ExperimentalResultsTheparametersvaluesandtestproblemsaresummarizedinTable 3-1 .Inall 83

PAGE 84

experiments,weset1=2=1,1=10)]TJ /F6 7.97 Tf 6.58 0 Td[(4and2=10)]TJ /F6 7.97 Tf 6.58 0 Td[(3fortheOSscheme.Forthealternatingminimization(AM)scheme( 3 ),weset=25andweusedthestoppingcriterion(fk)]TJ /F4 11.955 Tf 12.17 0 Td[(fk)]TJ /F6 7.97 Tf 6.59 0 Td[(1)=fk10)]TJ /F6 7.97 Tf 6.58 0 Td[(3.Forthisstoppingcriterion,theimagequalityfortheOSandtheAMschemesarecomparable.ThereconstructedimagesforSB512areshowninFigures 3-4B and 3-4C .Therootmeansquarederror(RMSE)oftheimageu,givenbyku)]TJ /F3 11.955 Tf 12.29 0 Td[(uk2=kuk2whereuisthereferenceimagereconstructedfromfullyacquireddata,was13.1%fortheOSschemeand12.5%fortheAMscheme.InFigures 3-4E and 3-4F ,wezoomintothesquareshowninFigure 3-4A .ItisseenthatbothmethodsadequatelyrecoveredtheimagewhileAMhasslightlybetterpreservededges.ThereconstructionsforAB512areshowninFigures 3-5B and 3-5C .TheRMSEwas14.7%forOSand13.9%forAM.SinceAB512doesnotcontainobviousnestructureswhencomparedtoSB512(seeFigure 3-4D ),wecomparethedifferencesofthesetworeconstructionstothereferenceimage(Figure 3-5A )underthesamecontrastlevel.Figure 3-5A showsthefollowings:referenceimage,reconstructionbyOSwithRMSE=14.7%,reconstructionbyAMwithRMSE=13.9%,andthedifferencesofreconstructionsbyOS(left)andAM(right)tothereferenceimagedisplayedunderthesamecontrastlevel.FromFigure 3-5D ,wecanseethereconstructionbyAM(right)hasmuchlesssystematicerrorthanthatofOS(left),andhenceismorelikelytopreventlossofimportantdiagnosticinformationsuchasedges.ToexaminetheefciencyoftheproposedalgorithmAMcomparedtoOS,wetrackedtheirobjectivefunctionvaluesandreconstructionerrorsduringthecomputationprocessesfordatasetsAB256andSB512.AstheperiterationcostsfortwoalgorithmsOSandAMarequitedifferent,wecomparedthenormalizedRMSEandobjectivefunctionvaluesversusCPUtime,whichareplottedinFigures 3-6A and 3-6B ,respectively.FromFigures 3-6A and 3-6B ,wecanseebothAMandOSconvergedfasterforthesmallerimageAB256thanforSB512.ForbothofAB256and 84

PAGE 85

AReference BOS CAM DReferenceBox EOSBox FAMBox Figure3-4. ReconstructedimagesbyOSandAMforSB512. SB512,AMconsistentlyreachedandmaintainedthesameorlowerRMSEandobjectivefunctionsthanOSinmuchlessCPUtime. 3.7ConcludingRemarksAfastnumericalalgorithmfortotalvariation-basedimagereconstructionwasdevelopedandanalyzed.Theproposedmethodemploysvariablesplitting,aquadraticpenalty,andanalternatingproximalminimizationalgorithm(AM).Linearconvergencewasestablishedwhenthesmoothpartoftheobjectivefunctionwasstronglyconvex,whiletheconvergencewassublinearunderaweakerconvexityassumption.Onephaseofthealternatingproximalminimizationalgorithmrepresentsatotalvariationdenoisingproblem,andtheotherphaseisanill-conditionedlinearsystem.AnimplementationbasedonasplitBregmanschemefortheTVproblemandaBarzilai-Borweinmethod 85

PAGE 86

AReference BOS CAM DDifferencestoReferenceimage Figure3-5. ReconstructionsbyOSandAMforAB512. forthelinearsystemwasproposed.Numericalperformancewasevaluatedusingimagereconstructionproblemsthatarosefromclinicalapplicationsofpartiallyparallelmagneticresonanceimaging(PPI).Theperformanceofthealternatingproximalminimizationalgorithmwascomparedtothatofanoperatorsplittingalgorithm.Thenumericalresultsshowexcellentperformancefortheproposedalgorithmintermsofefciencyandaccuracyinreconstruction,whichsuggestsitsgreatpotentialforpracticaluse. 86

PAGE 87

ARelativeerrorsversusCPUtime. BObjectivefunctionsversusCPUtime. Figure3-6. ComparisonofOSandAMondatasetsAB256andSB512. 87

PAGE 88

CHAPTER4INVERSECONSISTENTDEFORMABLEIMAGEREGISTRATIONOutline Thispaperpresentsanovelvariationalmodelforinverseconsistentdeformableimageregistration.Theproposedmodeldeformsbothsourceandtargetimagessimultaneously,andalignsthedeformedimagesinthewaythattheforwardandbackwardtransformationsareinverseconsistent.Toavoidthedirectcomputationoftheinversetransformationelds,ourmodelestimatestwomorevectoreldsbyminimizingtheirinvertibilityerrorusingthedeformationelds.Moreover,toimprovetherobustnessofthemodeltothechoiceofparameters,thedissimilaritymeasureintheenergyfunctionalisderivedusingthelikelihoodestimation.Theexperimentalresultsonclinicaldataindicatetheefciencyoftheproposedmethodwithimprovedrobustness,accuracyandinverseconsistency. 4.1BackgoundofConsistentImageRegistrationImageregistrationisaveryimportantsubjectthathasbeenwidelyappliedinmedicalresearchandclinicalapplications.Thetaskofimageregistrationistondatransformationeldthatrelatespointsinthesourceimagetotheircorrespondingpointsinthetargetimage.Deformableimageregistrationallowslocalizedtransformations,andisabletoaccountforinternalorgandeformations.Therefore,ithasbeenincreasinglyusedinhealthcaretoassistdiagnosisandtreatments.Inparticular,deformableimageregistrationhasbecomeacriticaltechniqueforimageguidedradiationtherapy.Itallowsmoreprecisetumortargetingandnormaltissuepreservation.Acomprehensivereviewofimageregistrationinradiationtherapycanbefoundin[ 40 ].Adeformableimageregistrationiscalledinverseconsistent,ifthecorrespondencebetweentwoimagesisinvarianttotheorderofthechoiceofsourceandtarget.Moreprecisely,letSandTbethesourceandtargetimages,andhandgbetheforwardand 88

PAGE 89

backwardtransformations,respectively,i.e.Sh=TandTg=S,thenaninverseconsistentregistrationsatiseshg=idandgh=id,whereidistheidentitymap.Thiscanbeillustratedbythefollowingdiagramwithconstraintsg=h)]TJ /F6 7.97 Tf 6.59 0 Td[(1,h=g)]TJ /F6 7.97 Tf 6.59 0 Td[(1: S g// T hoo ,(4)whereeachofthetwosquaresin( 4 )representsthedomainonwhichthelabeledimageisdened.Byapplyinganinverseconsistentregistration,measurementsorsegmentationsononeimagecanbepreciselytransferredtotheother.Inimagingguidedradiationtherapy,theinverseconsistentdeformableregistrationtechniqueprovidesthevoxel-to-voxelmappingbetweenthereferencephaseandthetestphaseinfour-dimensional(4D)radiotherapy[ 49 ].Thistechniqueisreferredtoautomaticre-contouring.Inverseconsistentdeformableimageregistrationhasbeenanactivesubjectofstudyintheliterature.Therehasbeenagroupofworkdevelopedinthecontextoflargedeformationbydiffeomorphicmetricmapping,e.g.[ 8 12 37 39 ].Themainideaofthismethodismodelingtheforwardandbackwardtransformationsasaone-parameterdiffeomorphismgroup.Then,ageodesicpathconnectingtwoimagesisobtainedbyminimizinganenergyfunctionalsymmetrictotheforwardandbackwardtransformations.Thistypeofmodelsproduceaverygoodregistrationresults.However,ittakelongtimetocompute,sincestrongregularizationofthemappingsarerequired.Variationalmethodisoneofthepopularapproachesforinverseconsistentdeformableimageregistration.Thismethodminimizesanenergyfunctional(s)symmetrictotheforwardandbackwardtransformations,andingeneral,consistsofthreeparts:regularizationofdeformationelds,dissimilaritymeasureofthetargetanddeformedsourceimages,andpenaltyofinverseinconsistency[ 5 23 58 77 ].In[ 23 ], 89

PAGE 90

ChristensenandJohnsonproposedtominimizethefollowingcoupledenergyfunctionalswithrespecttohandgalternately: 8><>:E(h)=Es(Sh,T)+Er(u)+kh)]TJ /F4 11.955 Tf 11.95 0 Td[(g)]TJ /F6 7.97 Tf 6.59 0 Td[(1k2L2()E(g)=Es(Tg,S)+Er(v)+kg)]TJ /F4 11.955 Tf 11.95 0 Td[(h)]TJ /F6 7.97 Tf 6.59 0 Td[(1k2L2(),(4)whereuandvareforwardandbackwarddeformationeldscorrespondingtohandg,respectively,i.e.h(x)=x+u(x)andg(x)=x+v(x).ThedissimilaritymeasureEsandtheregularizationofthedeformationeldEraredenedbyEs(Sh,T)=kSh)]TJ /F4 11.955 Tf 11.96 0 Td[(Tk2L2(),Er(u)=kau+br(divu))]TJ /F4 11.955 Tf 11.96 0 Td[(cuk2L2()withpositiveconstantsa,b,c>0.Thelastterminbothenergyfunctionalsenforcestheinverseconsistencyofhandg.Thesolution(u,v)to( 4 )isobtainedbyiterativelysolvingasystemoftwoevolutionequationsassociatedwiththeirEuler-Lagrange(EL)equations.Thismodelgivesconsiderablygoodresultswithparameterschosencarefully.However,itneedstocomputetheinversemappingsg)]TJ /F6 7.97 Tf 6.59 0 Td[(1andh)]TJ /F6 7.97 Tf 6.59 0 Td[(1explicitlyineachiteration,whichiscomputationallyintensivecancausecumulatednumericalerrorsintheestimationofinversemappings.Thevariationalmodelsdevelopedin[ 5 ]and[ 77 ]havethesameframeworkasin[ 23 ],butwithdifferentrepresentationsofEs,Er,andinverseconsistentconstraints.In[ 5 ]and[ 77 ]thetermskhg(x))]TJ /F4 11.955 Tf 12.78 0 Td[(xk2L2()andkgh(x))]TJ /F4 11.955 Tf 12.79 0 Td[(xk2L2()areusedintheenergyfunctionaltoenforcetheinverseconsistency.Byusingthesetermstheexplicitcomputationoftheinversetransformsofhandgcanbeavoidedduringtheprocessofndingoptimalforwardandbackwardtransformations.Thesimilaritymeasurein[ 77 ]ismutualinformationformulti-modalimageregistration.TheEs(Sh,T)in[ 5 ]iskSh)]TJ /F4 11.955 Tf 12 0 Td[(Tk2L2()=maxjDTj.TheregularizationtermEr(u)in[ 77 ]isafunctionofDu,andthatin[ 5 ]isatensorbasedsmoothingwhichisdesignedtopreventthetransformationeldsfrombeingsmoothedacrosstheboundariesoffeatures.In[ 71 72 ]theproposed 90

PAGE 91

modelsincorporatedstochasticerrorsintheinverseconsistentconstraintsforbothforwardandbackwardtransformations.In[ 45 ],Leowetal.proposedanon-variationalapproachthatupdatestheforwardandbackwardtransformationssimultaneouslybyaforcethatreducesthersttwotermsinE(h)andE(g)in( 4 )andpreservestheinverseconsistency.However,inordertosimplifythecomputationthisalgorithmonlytakeslinearordertermsintheTaylorexpressiontoapproximatetheinverseconsistentconditionsforupdatedtransformationelds.Asaconsequence,thetruncatingerrorscanbeaccumulatedandexaggeratedduringiterations.Thiscanleadtolargeinverseconsistenterror,despitethatitcanproduceagoodmatchingquickly[ 74 ].Inthispaperweproposeanovelvariationalmodeltoimprovetheaccuracy,robustnessandefciencyofinverseconsistentdeformableregistration.AsanalternatetothecurrentframeworkofvariationalmethodswhichndstheforwardandbackwardtransformationsthatdeformasourceimageStomatchatargetimageTandviceversa,weproposetodeformSandTsimultaneously,andlettheregistrationalignthedeformedsourceanddeformedtargetimages.ItisclearthatthedisparitybetweendeformedSanddeformedTissmallerthanthatbetweendeformedSandxedTordeformedTandxedS.Therefore,thedeformationbythebidirectionalsimultaneousdeformationsisingeneralsmallerthanthedeformationbyunidirectionaldeformationthatdeformsSfullwaytoTorTfullwaytoS.Asshowninsection 4.5 ,deformingSandTsimultaneouslyleadstoafasterandbetteralignmentthandeformingStothexedTorviceversa.Letand~representthetransformationeldssuchthatSmatchesT~.Itisnotdifculttoverifythatifand~areinvertible,thentheregistrationsfromStoT,andTtoSareinverseconsistent.Toavoidthedirectcomputationoftheinversetransformationsofand~,ourmodelseeksfortwoadditionaldeformationelds ,~ suchthatand areinversetoeachother,andthesamefor~and~ .Moreover,theregistrationprocessenforcescertainregularizationof 91

PAGE 92

thesefourdeformationelds,andalignsthedeformedSanddeformedT.Then,theoptimalinverseconsistenttransformationsfromStoT,andTtoScanbeobtainedsimplybyappropriatecompositionsofthesefourtransformations.TheideaofdeformingSandTsimultaneouslyhasbeenadoptedinthemodelswheretheforwardorbackwardtransformationismodeledasaone-parameterdiffeomorphismgroup[ 8 ].However,ourmodelndsregularizedinvertibledeformationeldsbyminimizingtheL2normsofthedeformationeldsandinverseconsistenterrorsratherthanaone-parameterdiffeomorphismgroup,whosecomputationalcostisveryexpensiveandhencehindersitsapplicationinclinicaluse.Moreover,ourmodelallowsparallelcomputationsforallthedeformationeldstosignicantlyreducethecomputationaltime.Furthermore,toimprovetherobustnessofthemodeltonoisesandthechoiceoftheparameterthatbalancesthegoodnessofmatchingandsmoothnessofthedeformationelds(seetheinE(h)andE(g)of( 4 )),weadoptthemaximumlikelihoodestimate(MLE)thatisabletoaccommodatecertaindegreeofvariabilityinmatchingtoimprovetherobustnessandaccuracyoftheregistration.ByusingMLE,theratioofweightingparametersonthesumofsquareddistance(SSD)oftheresidueimageS)]TJ /F4 11.955 Tf 12.46 0 Td[(T~andtheregularizationtermisnotaxed,but=2(see( 4 )below).Thisresultsinaself-adjustableweightingfactorthatmakesthechoiceofmoreexible,andalsospeedsuptheconvergencetotheoptimaldeformationeld.Therestofthepaperisorganizedasfollows.Insection 4.2 ,wepresentadetaileddescriptionoftheproposedmodel.Theexistenceofsolutionstotheproposedmodelisshowninsection 4.3 .Thecalculusofvariationandanoutlineofafastalgorithmforsolvingtheproposedmodelnumericallyareprovidedinsection 4.4 .Insection 4.5 ,wepresenttheexperimentalresultsonclinicaldata,andtheapplicationinautore-contouring.Thelastsectionconcludesthepaper. 92

PAGE 93

4.2ProposedMethodLetSandTbethesourceandtargetimagesdenedonSandTinRd,respectively.Notethat,inrealapplications,SandTareusuallyfullyoverlapped.ForsimplicityweassumethatimagesSandTarereal-valuedfunctionswithcontinuousderivatives.Letjjdenotetheabsolutevalue(length)ofascaler(vector)inEuclideanspaces,andkkdenotekkL2()henceforth.Wealsoextendthisnotationtovector-valuedfunctionswhosecomponentsareinL2orH1:u=(u1,,ud)>witheachcomponentuj2H1(),j=1,,d,thereiskukH1(),)]TJ /F2 11.955 Tf 5.48 -9.69 Td[(kuk2+kDuk21=2andkuk, dXj=1kujk2!1=2,kDuk, dXj=1kDujk2!1=2,wherekujk=Zjuj(x)j2dx1=2andkDujk=ZjDuj(x)j2dx1=2,forj=1,,d. 4.2.1MotivationandIdeasofProposedMethodInthispaper,weproposeanovelvariationalmodelforinverseconsistentdeformableregistrationtoimproveitsefciencyandrobustness.OurideadiffersfromthecurrentframeworkwhichdeformssourceimageStotargetimageT,orviceversa:asanalternate,weproposetodeformSandTsimultaneously,andmatchbothdeformedimages.Thismeansthatideallywepursuitforapairofhalf-waytransforms:S!Mand~:T!MsuchthatS=T~,whereMistheregionwhereSandT~haveoverlap.ToensurethetransformationsfromStoTandTtoSareinverseconsistent,thetransformsand~arerequiredtobeinvertible(butnotnecessarilytobeinversetoeachother).Hence,ourpurposeistondthetransformationsand~such 93

PAGE 94

that S=T~,,~invertible.(4)Toavoiddirectcomputationofinversesofand~duringiterations,weenforcetheinvertibilityofand~byndinganothertwotransformations :M!Sand~ :M!Tsuchthat =id, =id, (4) ~ ~=id,~~ =id.Onceweobtainedsuch and~ ,wecanconstructtheobjectivefull-waytransformationshandgasfollows,h=~ ,g=~ .Itiseasytoseethathandgsatisfytheinverseconsistentconstraintshg=gh=id.Thisideaisillustratedbythefollowingdiagram,whereMisanintermediateimage. M ~ ?????????? S ?? g// T ~__ hoo (4)SincebydeformingSandTsimultaneouslythedifferencebetweendeformedSanddeformedTateachiteration,ingeneral,issmallerthanthatbetweendeformedSandxedT,ordeformedTandxedS,thecomputationalcostofdeformingbothSandTismuchlessthantheconventionalonethatdeformSallthewaytoTandTtoS.Inparticular,iftheunderlyingdeformationsofhandgarelarge,deformingbothSandTcanmaketheeachdeformationofand~intheproposedmodelalmosthalfsmallerthanthatofhandg,andachieveafasterconvergenceforthecomputationofand~.Also,seeking and~ alongwithand~avoidsdirectcomputationofinverse 94

PAGE 95

transformationsineachiterationasthatin( 4 ),whichusuallycausescumulatederrorsduringiterationsifusingapproximationsoftheinverses.Moreover,regularizingthedeformationeldsisveryimportanttoobtainphysicallymeaningfulandaccurateregistrations.Also,iftheenergyfunctionalconsistsofonlydissimilaritymeasuresandinvertibleconstraints,itisill-posedingeneral.Therefore,weproposethefollowingframeworkfordeformableinverseconsistentregistration: min,~, ,~ R(,~, ,~ )+dis(S,T~),s.t.condition( 4 )holds(4)whereRisaregularizationoperatorofitsarguments,dis(S,T~)measuresthedissimilaritybetweenSandT~. 4.2.2AlternativeFormulationof( 4 )UsingDeformationFieldsLetthefunctionsu,~u,vand~vrepresentthecorrespondingdeformationeldsofthetransformations,~, and~ ,respectively.Thatis, (x)=x+u(x),~(x)=x+~u(x), (4) (x)=x+v(x),~ (x)=x+~v(x).Then,theconstraintsin( 4 )canberewrittenas u+v(x+u)=v+u(x+v)=0, (4) ~u+~v(x+~u)=~v+~u(x+~v)=0. 4.2.3MLEbasedderivationfordis(S,T~)Toimprovetherobustnessofthealgorithmfordeformableimageregistration,weusethenegativelog-likelihoodoftheresidueimageasameasureofmismatching.ConsidervoxelintensitiesoftheresidueimagedenedbyW(x),S(x))]TJ /F4 11.955 Tf 11.96 0 Td[(T~(x),x2M, 95

PAGE 96

asindependentsamplesdrawnfromaGaussiandistributionofmeanzeroandvariance2tobeoptimized(seeremarkbelowforthereasonofthisassumption),whoseprobabilitydensityfunction(pdf)isdenotedbyP(j).ThenthelikelihoodoftheresidualimageW(x)canbecomputedas L(jfW(x),x2g)=Yx2P(W(x)j)=Yx21 p 2ejS)]TJ /F9 7.97 Tf 6.59 0 Td[(T~j2=22.(4)Then,bywritingthesummationoverallx2asanintegraloverthenegativelog-likelihoodfunctionisgivenasfollows:kS)]TJ /F4 11.955 Tf 11.96 0 Td[(T~k2=22+jjlogp 2.Omittingtheconstantlogp 2,wedenethedissimilaritytermas dis(S,T~),kS)]TJ /F4 11.955 Tf 11.95 0 Td[(T~k2=22+jjlog.(4)whichcanberewrittenasourMLEttingtermFbyusingcorrespondingdeformationeldsuand~u: F(u,~u,),dis(S(x+u),T(x+~u))=kS(x+u))]TJ /F4 11.955 Tf 11.95 0 Td[(T(x+~u)k2=22+jjlog.(4)remarkLet^PbetheestimationofthepdffortherandomvariableX,W(x),x2.Weshowbelowwhyitisreasonabletoassume^PtobeaGaussiandistributionofzeromeanandvariance2.Infact,^PisafunctioninC0(R),thespaceofallthecontinuousfunctionsonreallinevanishingatinnitywiththesupremenorm.LetH0(R)betheHilbertspaceconsistingofalllinearcombinationsof(xl,x)fornitemanyofxl2R,where (xl,x)=)]TJ /F3 11.955 Tf 5.48 -9.69 Td[(22)]TJ /F6 7.97 Tf 6.59 0 Td[(1=2e)]TJ /F6 7.97 Tf 6.58 0 Td[((xl)]TJ /F9 7.97 Tf 6.59 0 Td[(x)2=22,8x2R.(4) 96

PAGE 97

DeneaninnerproductonH0(R)by*mXi=1ai(xi,),nXj=1bj(yj,)+H0(R)=mXi=1nXj=1aibj(xi,yj).Weclaimthat H0(R)isdenseinC0(R).(4)Infact,iftheclaim( 4 )isnottrue,byHahn-BanachtheoremthereexistsaboundedsignedmeasureminthedualspaceofC0(R),suchthat ZR^Pdm6=0,(4)butRRfdm=0,forallf2H0(R).Inparticular,foranyx2R,ZR(x,y)dmy=0,where(,)isasin( 4 ),andhence,ZRR(x,y)dmxdmy=0.Thisimpliesm=0,whichcontradicts( 4 ).Therefore,theclaimholds.Bythisclaimitiseasytoseethat ^P(z)kXl=1l(xl,z)=)]TJ /F3 11.955 Tf 5.48 -9.68 Td[(22)]TJ /F6 7.97 Tf 6.59 0 Td[(1=2kXl=1le)]TJ /F6 7.97 Tf 6.58 0 Td[((xl)]TJ /F9 7.97 Tf 6.59 0 Td[(z)2=22(4)forsomefxl;lgkl=1.SinceagoodregistrationrequiresthetheintensitiesoftheresidueimageW(x)closetozero.Hence,in( 4 )theonlydominateterminthesumshouldbetheonecorrespondingtoxl=0,andothertermsarenegligible.Thismeansthat^PisapproximatelyN(0,2),theGaussiandistributionwithmean0andvariance2. 4.2.4ProposedmodelBaseonthediscussionabove,wearereadytopresenttheproposedmodel.WedenetheregularizationtermR(,~, ,~ )in( 4 )usingtheircorresponding 97

PAGE 98

deformationeldsas R(,~, ,~ )=R(u,~u,v,~v),kDuk2+kD~uk2+kDvk2+kD~vk2.(4)Byplugging( 4 )and( 4 )into( 4 ),andreplacingtheconstraintin( 4 )by( 4 ),theproposedmodelcanbewrittenas: minu,~u,v,~v,R(u,~u,v,~v)+F(u,~u,),s.t.condition( 4 )holds,(4)whereR(u,~u,v,~v)andF(u,~u,)aredenedin( 4 )and( 4 ),respectively.Tosolveproblem( 4 ),werelaxtheequalityconstraintsofinverseconsistency,andpenalizetheirviolationusingquadraticfunctions,thenwriteitasanunconstrainedenergyminimizationproblem minu,~u,v,~v,R(u,~u,v,~v)+F(u,~u,)+(I(u,v)+I(~u,~v)),(4)whereandI(u,v)isthecostofinverseinconsistencyofuandv: I(u,v)=Iv(u)+Iu(v),(4)with Iv(u)=ku+v(x+u)k2andIu(v)=kv+u(x+v)k2.(4)Similarly,wehaveI(~u,~v).Withsufcientlylarge,solving( 4 )givesanapproximationtothesolutionof( 4 ).ThetermF(u,~u,)isfromthenegativelog-likelihoodoftheresidualimage( 4 ).Minimizingthistermforcesthemeanoftheresidueimagetobezero,butallowsittohaveavariancetoaccommodatecertainvariability.Thismakesthemodelmorerobusttonoiseandartifacts,andlesssensitivetothechoiceoftheparameterthanthemodelusingtheSSD,i.e.thesquaredL2-norm,oftheresidueimageasadissimilaritymeasureasin( 4 ).Theparameterbalancesthesmoothnessofdeformationeldsandgoodnessofalignments,andaffectstheregistrationresultsignicantly.Inthe 98

PAGE 99

proposedmodel,theratiooftheSSDoftheresidueimageoverthesmoothingtermsis=2ratherthanaprescribed.Sinceistobeoptimized,andfromitsELequationisthestandarddeviationoftheresidueimage.Therefore,intheproposedmodeltheweightonthematchingtermupdatesduringiterations.Whenthealignmentgetsbetter,thestandarddeviationoftheresidueasshownin( 4 )decreases,andhencetheweightonthematchingtermautomaticallyincreases.Thisself-adjustablefeatureoftheweightnotonlyenhancestheaccuracyofalignment,butalsomakesthechoiceofexible,andresultsinafastconvergence.Asshownearlier,thenalforwardandbackwardtransformshandgcanbeobtainedbyh=~ =x+~v+u(x+~v)andg=~ =x+~u+v(x+~u).Thus,thecorrespondingnalfull-wayforwardandbackwarddeformationeldsuandvaregivenas u=~v+u(x+~v)andv=~u+v(x+~u),(4)respectively.Thentheinverseconsistentconstraints( 4 )canberepresentedusingu,vasfollows: u+v(x+u)=v+u(x+v)=0.(4) 4.3ExistenceofSolutionsInthissectionweprovetheexistenceofsolutions(u,~u,v,~v,)totheproposedmodel( 4 ).Forsimplicity,weassumethatbothSandTdenedonthesamedomain,whichissimplyconnected,closedandboundedinRdwithLipschitzboundary@.AlsoS,T2C1().Asinreality,deformationeldcannotbeunbounded,werestrictu,~u,v,~vtobeinaclosedsubsetofL1():B,fu2L1():kukL1()B,B2R+onlydependsong 99

PAGE 100

Then,weseeksolutions(u,~u,v,~v,)totheproblem( 4 )inthespacesu,~u,v,~v2H1()\Band2R+.Forshortnotations,weletwdenotethequaternion(u,~u,v,~v).Then,weshowtheexistenceofsolutionstothefollowingminimizationproblem: min(w,)2(H1\B)R+E(w,)(4)whereE(w,)=kDwk2+F(w,)+I(w)andFandIaredenedcorrespondinglyin( 4 )usingthesimpliednotationofw,i.e.kDwk2=kDuk2+kD~uk2+kDvk2+kD~vk2,F(w,)=kS(x+u))]TJ /F4 11.955 Tf 11.96 0 Td[(T(x+~u)k2=2+jjlog,I(w)=Iv(u)+Iu(v)+I~v(~u)+I~u(~v).andthetermsontherightsideofI(w)aredenedasin( 4 ).Theandareprescribedpositiveconstants. Theorem4.1. Theminimizationproblem( 4 )admitssolutions(w,)2)]TJ /F4 11.955 Tf 5.48 -9.68 Td[(H1\BR+. Proof. For(w,)2)]TJ /F4 11.955 Tf 5.48 -9.69 Td[(H1\BR+,E(w,)isboundedbelow.Hence,thereexistsaminimizingsequencef(wk,k)g1k=1)]TJ /F4 11.955 Tf 5.48 -9.69 Td[(H1\BR+suchthatlimk!1E(wk,k)=inf(H1\B)R+E(w,).ThereforefkDwkkg1k=1areuniformlybounded.Alongwithwk2Bweknowthatfwkg1k=1isaboundedsequenceinH1.BytheweakcompactnessofH1andthefactthatH1isprecompactinL2,thereexistsaconvergentsubsequence,whichisstilldenotedbyfwkg1k=1,andafunction^w2H1,suchthat wk*^wweaklyinH1,(4) wk!^wstronglyinL2,anda.e.in.(4) 100

PAGE 101

Moreover,sinceE(wk,k)!1ifk!0or1,thereisaconstantC>0suchthatfkg1k=1areboundedbelowandaboveby1=CandCrespectively.Hence,thereisasubsequenceoffkg1k=1andascaler^2R+,withoutchangingthenotationforthesubsequencewehave k!^2R+.(4)Fromtheweaklowersemi-continuityofnormsand( 4 ),weknow kD^wk2limk!1kDwkk2.(4)Also,asI(w)8Bforanyw2H1\Bandwk!^wa.e.in,weget,bydominantconvergencetheorem,that limk!1I(wk)=I(^w).(4)BythesameargumentwiththesmoothnessofSandT,theconvergenceoffkg1k=1,andthefactthatwk!^wa.e.in,wecanalsohave limk!1F(wk,k)=F(^w,^)(4)Combining( 4 ),( 4 )with( 4 ),weobtainthatE(^w,^)limk!1E(wk,k)=inf(H1\B)R+E(w,).Furthermore,sincefwkg1k=1BL1(),weknowwk*w^wweakly*inL1andhence^w2B.Therefore,(^w,^)2)]TJ /F4 11.955 Tf 5.48 -9.68 Td[(H1\BR+.HenceE(^w,^)=inf(H1\B)R+E(w,).whichimpliesthat(^w,^)isasolutiontotheminimizationproblem( 4 ). 101

PAGE 102

4.4NumericalSchemeInthissection,weprovidethenumericalschemeforsolving( 4 ).SincethecompositionsintheinverseconsistencyconstraintsIu(v)andIv(u)bringadifcultyingettinganexplicitformoftheELequationsforthedeformationeldsandtheirinverses,inourcomputation,insteadofdirectlysolving( 4 ),wesolvethefollowingtwocoupledminimizationproblemsalternately: 8>><>>:minu,~uEv,~v,(u,~u)minv,~vEu,~u(v,~v)(4)where Ev,~v,(u,~u,)=kDuk2+kD~uk2+F(u,~u,)+(Iv(u)+I~v(~u))(4)and Eu,~u(v,~v)=kDvk2+kD~vk2+(Iu(v)+I~u(~v)).(4)Bytakingrstvariationwithrespecttou,~u,v,~v,wegettheELequations: 8>>>>>>><>>>>>>>:)]TJ /F3 11.955 Tf 9.3 0 Td[(u+ 2Wu,~uDS(x+u)+hI+Dv(x+u),u+v(x+u)i=0)]TJ /F3 11.955 Tf 9.3 0 Td[(v+hI+Du(x+v),v+u(x+v)i=0)]TJ /F3 11.955 Tf 9.3 0 Td[(~u)]TJ /F5 11.955 Tf 15.74 8.09 Td[( 2Wu,~uDT(x+~u)+hI+D~v(x+~u),~u+~v(x+~u)i=0)]TJ /F3 11.955 Tf 9.3 0 Td[(~v+hI+D~u(x+~v),~v+~u(x+~v)i=0, (4) in,withfreeNeumannboundaryconditionsforeachofthemon@: hDu,ni=hD~u,ni=hDv,ni=hD~v,ni=0,on@,(4)whereWu,~u,S(x+u))]TJ /F4 11.955 Tf 12.3 0 Td[(T(x+~v),Iistheidentitymatrixofsized,andnistheouternormalof@.Also,therstvariationofgives =kS(x+u))]TJ /F4 11.955 Tf 11.96 0 Td[(T(x+~u)k=jj1=2.(4) 102

PAGE 103

ThesolutiontotheELequations( 4 )canbeobtainedbyndingthestationarysolutiontotheevolutionequationsassociatedwiththeELequations.Innumericalimplementation,weusesemi-implicitdiscreteformoftheevolutionequations.Theadditiveoperatorsplitting(AOS)schemewasappliedtosolvetheproblemfaster[ 65 ].AnalternativewayofAOStosolvethesemi-implicitdiscreteevolutionequationinthiscasecanbeobtainedbyapplyingdiscretecosinetransforms(DCT)todiagonalizetheLaplaceoperatorwiththeassumptionthatthedeformationeldshavesymmetricboundarycondition,whichiscompatiblewith( 4 ).Intwo-dimensional(2D)case,thesemi-implicitdiscreteformof( 4 )withxedstepsizesu,vfortheevolutionequationsofu(k+1)as u(k+1)i,j)]TJ /F4 11.955 Tf 11.95 0 Td[(u(k)i,j u=i,ju(k+1))]TJ /F4 11.955 Tf 11.95 0 Td[(Di,j)]TJ /F5 11.955 Tf 5.48 -9.69 Td[(F)]TJ /F4 11.955 Tf 5.48 -9.69 Td[(u(k),~u(k),(k)+Iv(k))]TJ /F4 11.955 Tf 5.48 -9.69 Td[(u(k),(4)andv(k+1)as v(k+1)i,j)]TJ /F4 11.955 Tf 11.96 0 Td[(v(k)i,j v=i,jv(k+1))]TJ /F5 11.955 Tf 11.96 0 Td[(Di,jIu(k))]TJ /F4 11.955 Tf 5.48 -9.68 Td[(v(k),(4)wherei,jandDi,jrepresentthediscreteLaplacianandgradientoperatorsatthepixelindexedby(i,j),respectively.The3Dcaseisasimpleanaloguewithonemoresubscriptinindices.Similarly,wehavethediscreteevolutionequationfor~uand~vwiththetwocomponentswithineachofthethreepairs(u,~u),(v,~v)and(S,T)switchedin( 4 )and( 4 ).WithAOSschemebeingapplied,thecomputationforeachupdateofuinvolvesofsolvingdtridiagonalsystemswhosecomputationalcostsarelinearinN,whereNisthetotalnumberofpixelsinS(orT).Also,ineachiterationofupdatinguandv,thereneeds2(d+1)interpolationswithsizeN.Itisimportanttopointoutthat,ineachiteration,thecomputationsofu,~u,v,~vcanbecarriedoutinparallel.WesummarizeicDIRinAlgorithm 5 ,wherethemaximuminverseconsistencyerror(ICE)cisdenedby c=maxxfju+v(x+u)j,jv+u(x+v)jg,(4) 103

PAGE 104

Algorithm5InverseConsistentDeformableImageRegistration(icDIR) InputS,T,andu,v,,>0,=.5,c=1.Initialize(u(0),~u(0),v(0),~v(0))=0,k=0. whilecdo repeat fAlltermsin(u(k+1),~u(k+1),v(k+1),~v(k+1))canbecalculatedinparallelg Calculate(u(k+1),v(k+1))using( 4 )and( 4 ). Calculate(~u(k+1),~v(k+1))using( 4 )and( 4 )with(u(k+1),v(k+1))replacedby(~u(k+1),~v(k+1)). update(k+1)by( 4 ). k k+1 untilconvergence return(u,~u,v,~v) (u,~u,v,~v) (u,~u,v,~v), 2. Computeuandvusing( 4 )andthencusing( 4 ). endwhile anduandvarethenalfull-waydeformationeldsshownin( 4 ).Thatis,itmeasuresthemaximumICEofdeformationsobtainedbyquaternion(u,~u,v,~v).Theparameterin( 4 )mayincreaseduringiterationstoensuresmallerICE.Ineachinnerloopwithxed,thecomputationisterminatedwhenthemeanofCC(S(x+u),T)andCC(T(x+v),S)converges.Wesetastoppingtolerance=.5andterminatethewholecomputationoncecislowerthan,inwhichcasethemaximumICEislessthanhalfofthegridsizebetweentwoconcatenatepixels/voxelsandhencetheinverseconsistencyisexactlysatisedwithrespecttotheoriginalresolutionoftheimages. 4.5ExperimentalResultsInthissection,wepresenttheexperimentalresultsofproposedmodelusingalgorithm1(icDIR).AllimplementationsinvolvedintheexperimentswerecodedinMATLABRv7.3(R2006b),excepttheThomastridiagonalsolver,whichwascodedinC++.Weusedbuild-infunctionsinterp2/interp3ofMatlabwithdefaultsettingsforinterpolations.AllComputationswereperformedonaGNU/Linux(version2.6.16)workstationwithIntelRCore2CPUat1.86GHzand2GBmemory.Wersttesttheaccuracyofregistrationandautore-contouringoftheproposedalgorithmonaclinicaldatasetof1002D-prostateMRimages.Eachimage,calleda 104

PAGE 105

phase,isa2Dimageofdimension288192thatfocusesontheprostatearea.TherstphaseisusedasasourceimageS,asshowninFigure 4-1A .Theboundariesoftheregionsofinterests(ROI)inSweredelineatedbycontoursandsuperimposedbymedicalexperts,asenlargedandshowninFigure 4-4A .Therest99phaseswereconsideredastargets.Inthisexperimentweappliedtheproposedmodel( 4 )withparameters(,,)settobe(.05,.2,.05)toSandTs.Fordemonstration,weonlyshowedtheresultusingthe21stphaseasT,asdepictedinFigure 4-1B .ThedeformedTanddeformedS,i.e.T(x+v)andS(x+u),areshownintheFigure 4-1C and 4-1D respectively,whereuandvaredenedin( 4 )usingtheoptimal(u,~u,v,~v)obtainedbymodel( 4 ).Theerrorsofthealignments,jT(x+v))]TJ /F4 11.955 Tf 12.74 0 Td[(SjandjS(x+u))]TJ /F4 11.955 Tf 12.74 0 Td[(Tj,onthesquaredarea(showninFigure 4-1A )aredisplayedinFigure 4-2A and 4-2C ,respectively.WithcomparisontotheoriginalerrorjS)]TJ /F4 11.955 Tf 12.02 0 Td[(TjshowninFigure 4-2B ,wecanseetheerrorsofalignmentsaresignicantlyreduced.Thisindicatesthattheproposedregistrationmodel( 4 )hashighaccuracyinmatchingtwoimages.Thenaloptimalforwardandbackwarddeformationeldsuandvaredisplayedbyapplyingthemtoadomainofregulargrids,showninFigure 4-3A and 4-3C ,respectively.Furthermore,tovalidatetheaccurateinverseconsistencyobtainedbyourmodel( 4 ),weappliedu+v(x+u)onadomainwithregulargrids,andplottedtheresultinggridsinFigure 4-3B .Theresultinggridsbyv+u(x+v)hadthesamepatternsoweomittedithere.FromFigure 4-3B ,wecanseethattheresultinggridsarethesameastheoriginalregulargrids.Thisindicatesthattheinverseconsistentconstraintsu+v(x+u)=v+u(x+v)=0arewellpreserved.WealsocomputedthemaximumICEcusingu,vand( 4 )andtheresultwas.46.ThemeanICE(ku+v(x+u)k+kv+u(x+v)k)=2jjversusthenumberofiterationsisplottedinFigure 4.5 ,whichshowstheinverseconsistencyispreservedduringtheregistration.Theseimplythattheproposedalgorithmprovidesanaccurateinverseconsistentregistration. 105

PAGE 106

Anaccurateinverseconsistentregistrationcantransformsegmentationsfromoneimagetoanotheraccurately.Oneoftheapplicationsisautore-contouring,thatdeformstheexpert'scontoursfromaplanningimagetonewimagesduringthecourseofradiationtherapy.Inthisexperiment,wehadexpert'scontourssuperimposedonthesourceimageSasshowninFigure 4-4A .ThenbyapplyingthedeformationelduonthiscontourswegetthedeformedcontoursonthetargetimageTasshowninFigure 4-4B .Theaccuracyinautore-contouringisevident.Figure 4-3 shows,fromlefttoright,thefollowings:u,u+v(x+u),whichdemonstratestheinverseconsistencyiswellpreserved,andv,respectively.Thesecondexperimentwasaimedtotesttheefciencyoftheproposedmodel( 4 )inregistering3Dimages.Weapplied( 4 )toapairof3DchestCTimagesofdimension648348takenfromthesamesubjectbutatdifferentperiods.Theparameters(,,)weresettobe(.05,.1,.004).Theregistrationwasperformedin3D,butfordemonstration,weonlyshowthecorrespondingaxial(xyplanewithz=33),sagittal(yzplanewithx=25)andcoronal(zxplanewithy=48)slices.TheregistrationresultsareplottedinFigures 4-5 4-6 and 4-7 ,respectively.Ineachgure,theimagesintheupperrowareSandT,respectively,andtheimagesinthemiddlerowaredeformedTandS,i.e.T(x+v)andS(x+u),respectively.ThebottomrowshowstheresidualimagesjS(x+u))]TJ /F4 11.955 Tf 12.14 0 Td[(Tj,jS)]TJ /F4 11.955 Tf 12.14 0 Td[(TjandjT(x+v))]TJ /F4 11.955 Tf 12.14 0 Td[(Sj.ThemeanofCC(S(x+u),T)andCC(T(x+v),S)reached.998after50iterations,andthemeanofinverseconsistencyerrorswas.015.Theresultsshowsthehighaccuracyofproposedmodel( 4 )andthewellpreservedinverseconsistency.Thethirdexperimentwasaimedtocomparetheeffectivenessofmodel( 4 )withthefollowingconventionalfull-wayinverseconsistentdeformableregistrationmodel: minu,v,u,vkDuk2+kDvk2+J(u,v,u,v)+(Iv(u)+Iu(v))(4) 106

PAGE 107

whereuandvareforwardandbackwarddeformationelds,respectively,andthetermJisdenedbyJ(u,v,u,v)=kS(x+u))]TJ /F4 11.955 Tf 11.95 0 Td[(Tk2=22u+kT(x+v))]TJ /F4 11.955 Tf 11.95 0 Td[(Sk2=22v+jjlogvv.Thecomparisonismadeontheefciencyandaccuracyofmatching,aswellasthepreservationofinverseconsistency.Theaccuracyofmatchingismeasuredbycorrelationcoefcients(CC)betweenthetargetimageanddeformedsourceimagewiththeoptimalforwardandbackwarddeformationsobtainedbymodel( 4 )andproposedmodel( 4 ),respectively.RecallthatforanytwoimagesSandTbothwithNpixels,theCCofSandTisdenedbyCC(S,T)=PNi=1(Si)]TJ /F3 11.955 Tf 12.84 2.66 Td[(S)(Ti)]TJ /F3 11.955 Tf 13.94 2.66 Td[(T) q PNi=1(Si)]TJ /F3 11.955 Tf 12.84 2.66 Td[(S)2PNi=1(Ti)]TJ /F3 11.955 Tf 13.95 2.66 Td[(T)2,whereSiandTiaretheintensitiesattheithpixelsofSandT,respectively,SandTarethemeanintensitiesofSandT,respectively.ThemaximumvalueofCCis1,inwhichcaseSandTare(positively)linearlyrelated.Inthisexperimentweappliedmodels( 4 )and( 4 )totheimagesintherstexperimentshowninFigure 4-1 withthesameparameters(,,)tobe(.05,.2,.05).InFigure 4.5 ,weplottedtheCCobtainedbymodel( 4 )andproposedmodel( 4 )ateachiteration.OnecanobservedthattheCCobtainedbymodel( 4 )ishigherandincreasesfasterthanmodel( 4 ).Thisdemonstratesthatproposedmodel( 4 )ismoreefcientthantheconventionalfull-waymodel.ThereasonisthatthedisparitybetweendeformedSanddeformedTissmallerthanthatbetweendeformedSandxedTordeformedTandxedS.WhenSandTaredeformedsimultaneously,thetwodirectionaldeformationeldsarenotnecessarilytobelargeeveniftheunderlyingdeformationeldislarge,whichusuallymakesitdifcultforthefull-waybasedregistrationmodeltoreachasatisfactoryalignmentinshorttime. 107

PAGE 108

Table4-1. NumberofiterationsusedforconvergenceandthenalCCobtainedbyproposedmodelwithupdated/xed. Update Fix CCIterCCIter 1e2.96248.955891e1.96297.9464201e0.960356.9331762 ThelastexperimentisaimtotesttherobustnessofthemodeltonoisesandthechoiceoftheparameterwiththeuseofMLEbasedapproach( 4 )formeasuringthegoodnessofmatching.TheimagesSandTinFigure 4-1 withadditiveGaussiannoises(standarddeviationis3%oflargestintensityvalueofS)wereusedinthisexperiment.TheCCbetweenSandTbeforeregistrationisCC(S,T)=.901.Weappliedmodel( 4 )withtobeupdated/optimizedbyitsELequation( 4 ),andtobeset=1,thatisthesameasusingSSDassimilaritymeasure,respectively,tothenoisedatamentionedabove.Weproceededtheregistrationwithvariousvaluesof,butkeptotherparametersxed.Thenthenumbersofiterations(Iter)forconvergenceandthenalCCwererecordedandshowninTable 4-1 .Onecanseethatwhiledecreases,theaccuracyofmodel( 4 )usingxedreducesasthenalCCbecomemuchsmaller,anditalsotakesmuchlongertimeforthealgorithmtoconverge.Ontheotherhand,withbeingupdated(whosecomputationalcostisextremelycheap)model( 4 )canobtaingoodmatchinginmuchlessiterationsforalargerangeof.ThisshowsthatmodelwithMLEttingismuchlesssensitivetonoiseandthechoiceof,andcanachievefastandaccurateresultscomparedwiththemodelusingSSDtomeasuremismatching. 108

PAGE 109

ASourceimageS BTargetimageT CDeformedT DDeformedS Figure4-1. Inverseconsistentregistrationresultbyproposedmodel( 4 ). 109

PAGE 110

AjS(x+u))]TJ /F25 9.963 Tf 9.97 0 Td[(Tj BjS)]TJ /F25 9.963 Tf 9.96 0 Td[(Tj CjT(x+v))]TJ /F25 9.963 Tf 9.97 0 Td[(Sj Figure4-2. Residueimageobtainedbyproposedmodel( 4 ). Au Bu+v(x+u) Cv Figure4-3. Deformationeldsobtainedbyproposedmodel( 4 )inthezoomed-inareaappliedonregulargridwithhalfoforiginalresolutionofimages. AOriginalContouronS BRe-contouringonT Figure4-4. Autore-contouringresultusingthedeformationelduobtainedbyproposedmodel( 4 ). 110

PAGE 111

Figure4-5. Registrationresultofproposedmodel( 4 )appliedto3DchestCTimage.Thisgureshowsthez=33sliceataxialdirection. 111

PAGE 112

Figure4-6. Registrationresultofproposedmodel( 4 )appliedto3DchestCTimage.Thisgureshowsthex=25sliceatsagittaldirection. 112

PAGE 113

Figure4-7. Registrationresultofproposedmodel( 4 )appliedto3DchestCTimage.Thisgureshowsthey=48sliceatcoronarydirection. 113

PAGE 114

Figure4-8. CCineachiterationobtainedbyfull-waymodel( 4 )andproposedmodel( 4 ). Figure4-9. Meanofinverseconsistenterrors(ICE)ofthenaldeformationeldsobtainedbyusingfull-waymodel( 4 )andproposedmodel( 4 ). 114

PAGE 115

REFERENCES [1] Acker,F.andPrestel,M.-A.Convergenced'unschemademinimisationalternee.Ann.Fac.Sci.ToulouseMath.(5)2(1980):1. [2] Aharon,Michal,Elad,Michael,andBruckstein,Alfred.TheK-SVD:AnalgorithmforDesigningofOvercompleteDictionariesforSparseRepresentation.IEEETrans.SignalProcess.54(2006).11:4311. [3] Akaike,H.Onasuccessivetransformationofprobabilitydistributionanditsapplicationtotheanalysisoftheoptimumgradientmethod.Ann.Inst.Statist.Math.Tokyo11(1959):1. [4] .Onasuccessivetransformationofprobabilitydistributionanditsapplicationtotheanalysisoftheoptimumgradientmethod.Ann.Inst.Statist.Math.Tokyo11(1959):1. [5] Alvarez,L.,Deriche,R.,Papadopoulo,T.,andSanchez,J.Symmetricaldenseopticalowestimationwithocclusionsdetection.ProceedingsofEuropeanConferenceonComputerVision(2002):721. [6] Arunachalam,A.,Samsonov,A.,andBlock,W.F.Self-calibratedGRAPPAmethodfor2Dand3Dradialdata.Magn.Reson.Med.57(2007).5:931. [7] Attouch,H.,Redont,P.,andSoubeyran,A.Anewclassofalternatingproximalminimizationalgorithmswithcost-to-move.SIAMJ.Optim.18(2007):1061. [8] Avants,B.B.,Grossman,M.,andGee,J.C.Symmetricdiffeomorphicimageregistration:Evaluatingautomatedlabelingofelderlyandneurodegenerativecortexandfrontallobe.ProceedingsofBiomedicalImageRegistration4057(2006):50. [9] Barzilai,J.andBorwein,J.M.Twopointstepsizegradientmethods.IMAJ.Numer.Anal.8(1988):141. [10] Bauschke,H.H.,Combettes,P.L.,andNoll,D.JointminimizationwithalternatingBregmanproximityoperators.Pac.J.Optim.2(2006):401. [11] Beck,AmirandTeboulle,Marc.AFastIterativeShrinkage-ThresholdingAlgorithmforLinearInverseProblems.SIAMJ.Imag.Sci.2(2009).1:183. [12] Beg,MirzaFaisalandKhan,Ali.SymmetricDataAttachmentTermsforLargeDeformationImageRegistration.IEEETransactionsonMedicalImaging26(2007).9:1179. [13] Bertsekas,D.ParallelandDistributedComputation.PrenticeHall,1989. 115

PAGE 116

[14] Bioucas-Dias,J.andFigueiredo,M.AnewTwIST:Two-stepiterativeshrinkage/thresholdingalgorithmsforimagerestoration.IEEETrans.ImagePro-cess.16(2007).12:2992. [15] Block,K.,Uecker,M.,andFrahm,J.UndersampledRadialMRIwithMultipleCoils:IterativeImageReconstructionUsingaTotalVariationConstraint.Magn.Re-son.Med.57(2007):1086. [16] Block,K.T.,Uecker,M.,andFrahm,J.UndersampledRadialMRIwithMultipleCoils.IterativeImageReconstructionUsingaTotalVariationConsraint.Magn.Re-son.Med.57(2007):1086. [17] Candes,EmmanuelJ.andRomberg,JustinK.SignalRecoveryfromRandomProjections.ProceedingsofSPIEComputationalImagingIII5674(2005):76. [18] Candes,EmmanuelJ.,Romberg,JustinK.,andTao,Terence.RobustUncertaintyPrinciples:ExactSignalReconstructionfromHighlyIncompleteFrequencyInformation.IEEETrans.Inform.Theory.52(2006).2:489. [19] Chambolle,A.Analgorithmfortotalvariationminimizationandapplications.J.Math.ImagingVis.20(2004):89. [20] Chan,T.F.,Golub,G.H.,andMulet,P.Anonlinearprimal-dualmethodfortotalvariationbasedimagerestoration.SIAMJ.Optim.20(1999):1964. [21] Chang,T.C.,He,L.,andFang,T.MRImageReconstructionfromSparseRadialSamplesUsingBregmanIteration.Proc.Intl.Soc.Mag.Reson.Med.(2006):696. [22] Chen,S.,Donoho,D.,andSaunders,M.AtomicDecompositionbyBasisPursuit.SIAMJ.Sci.Comput.20(1999):33. [23] Christensen,G.E.andJohnson,H.J.ConsistentImageRegistration.IEEETransactionsonMedicalImaging20(2001).7:721. [24] Dai,Y.H.,Hager,W.W.,Schittkowski,K.,andZhang,H.ThecyclicBarzilai-Borweinmethodforunconstrainedoptimization.IMAJ.Numer.Anal.26(2006):604. [25] Dai,Y.H.andYuan,Y.Alternateminimizationgradientmethod.IMAJ.Numer.Anal.23(2003):377. [26] Donoho,DavidL.ForMostLargeUnderdeterminedSystemsofLinearEquations,theMinimall1-normSolutionisAlsotheSparsestSolution.Com-mun.PureAppl.Math.59(2006).7:907. [27] Donoho,DavidL.,Elad,M.,andTemlyakov,V.StableRecoveryofSparseOvercompleteRepresenationinthePresenceofNoise.IEEETrans.Inf.Theory52(2006):6. 116

PAGE 117

[28] Eckstein,J.andBertsekas,D.OntheDouglas-Rachfordsplittingmethodandtheproximalpointalgorithmformaximalmonotoneoperators.MathematicalProgramming55(1992).1-3:293. [29] Eggers,H.andBoesiger,P.Gridding-andConvolution-BasedIterativeReconstructionwithVariablek-spaceResolutionforSensitivity-EncodedNon-CartesianImaging.Proc.Intl.Soc.Mag.Reson.Med..2003,2346. [30] Elad,MichaelandAharon,Michal.ImageDenoisingViaSparseandRedundantRepresentationsOverLearnedDictionaries.IEEETrans.ImageProcess.15(2006).12:3736. [31] Friedlander,A.,Martnez,J.M.,Molina,B.,andRaydan,M.Gradientmethodwithretardsandgeneralizations.SIAMJ.Numer.Anal.36(1999):275. [32] Gabay,D.andMercier,B.Adualalgorithmforthesolutionofnonlinearvariationalproblemsvianite-elementapproximations.Comput.Math.Appl.2(1976):17. [33] Glowinski,R.andMarrocco,A.Surl'approximationparelementsnisd'ordreun,etlaresolutionparpenalisation-dualited'uneclassedeproblemesdeDirichletnonlineaires.RAIROAnalyseNumerique9(1975):41. [34] Goldstein,T.andOsher,S.TheSplitBregmanMethodforL1RegularizedProblems.SIAMJ.Imag.Sci.2(2009):323. [35] Griswold,M.A.,Jakob,P.M.,Heidemann,R.M.,Mathias,N.,Jellus,V.,Wang,J.,Kiefer,B.,andHaase,A.GeneralizedAutocalibratingPartiallyParallelAcquisitions(GRAPPA).Magn.Reson.Med.47(2002):1202. [36] Hale,E.T.,Yin,W.,andZhang,Y.Fixed-PointContinuationMethodfor`1-MinimizationwithApplicationstoCompressedSensing.SIAMJ.Optim.19(2008).3:1107. [37] He,J.C.andChristensen,G.E.Largedeformationinverseconsistentelasticimageregistration.ProceedingsofInformationProcessinginMedicalImaging2732(2003):438. [38] He,Lin,Chang,Ti-Chiun,Osher,Stanley,Fang,Tong,andSpeier,Peter.MRImageReconstructionbyUsingtheIterativeRenementMethodandNonlinearInverseScaleSpaceMethods.Tech.Rep.06-35,CAMUCLA,2006. [39] Joshiand,S.,Davis,B.,Jomier,M.,andGerig,G.Unbiaseddiffeomorphicatlasconstructionforcomputationalanatomy.Neuroimage,Supplement23(2004).1:151. [40] Kessler,M.L.Imageregistrationanddatafusioninradiationtherapy.BritishJournalonRadiology79(2006):99. 117

PAGE 118

[41] Kim,Seung-Jean,Koh,K.,Lustig,M.,Boyd,Stephen,andGorinevsky,Dimitry.AnInterior-PointMethodforLarge-Scale`1-RegularizedLeastSquares.SIAMJ.Optim.1(2007).4:606. [42] Kim,S.J.,Koh,K.,Lustig,M.,andBoyd,S.AnEfcientMethodforCompressedSensing.ProceedingsofIEEEIntl.Conf.ImageProcess.(ICIP)3(2007):117. [43] Krishnan,D.,Lin,P.,andTai,X.Anefcientoperatorsplittingmethodfornoiseremovalinimages.Commun.Comput.Phys.1(2006):847. [44] Larson,A.C.,White,R.D.,Laub,G.,McVeigh,E.R.,Li,D.,andSimonetti,O.P.Self-gatedcardiaccineMRI.Magn.Reson.Med.51(2004):93. [45] Leow,A.D.,Huang,S.C.,Geng,A.,Becker,J.,Davis,S.,Toga,A.,andThompson,P.Inverseconsistentmappingin3Ddeformableimageregistration:itsconstructionandstatisticalproperties.ProceedingsofInformationProcessinginMedicalImaging(2005):493. [46] Liao,H.Y.andSapiro,Guillermo.SparseRepresentationforLimitedDataTomography.ProceedingsofIEEEIntl.Symp.Biomed.Imag.(ISBI)(2008):1375. [47] Lions,P.L.andMercier,B.Splittingalgorithmsforthesumoftwononlinearoperators.SIAMJ.Numer.Anal.16(1979):964. [48] Liu,C.,Bammer,R.,andMoseley,M.ParallelImagingReconstructionforArbitraryTrajectoriesUsingk-SpaceSparseMatrices(kSPA).Magn.Reson.Med.58(2007):1071. [49] Lu,W.,Olivera,G.H.,Chen,Q.,Chen,M.,andRuchala,K.Automaticre-contouringin4Dradiotherapy.PhysicsinMedicineandBiology51(2006):1077. [50] Lustig,Michael,Donoho,David,andPauly,JohnM.SparseMRI:TheApplicationofCompressedSensingforRapidMRImaging.Magn.Reson.Med.58(2007).6:1182. [51] Ma,Shiqian,Yin,Wotao,Zhang,Yin,andChakraborty,A.AnefcientalgorithmforcompressedMRimagingusingtotalvariationandwavelets.IEEEProc.Conf.onComp.Vis.Patt.Recog.(2008):1. [52] Nesterov,Yu.Gradientmethodsforminimizingcompositeobjectivefunction.COREDiscussionPapers2007/76,UniversitcatholiquedeLouvain,CenterforOperationsResearchandEconometrics(CORE),2007.URL http://ideas.repec.org/p/cor/louvco/2007076.html [53] Nocedal,J.andWright,S.J.NumericalOptimization.NewYork:Springer,1999. 118

PAGE 119

[54] Pruessmann,K.P.,Weiger,M.,Bornert,P.,andBoesiger,P.AdvancesinSensitivityEncodingWithArbitraryk-SpaceTrajectories.Magn.Reson.Med.46(2001):638. [55] Pruessmann,K.P.,Weiger,M.,Scheidegger,M.B.,andBoesiger,P.SENSE:SensitivityencodingforfastMRI.Magn.Reson.Med.42(1999):952. [56] Qian,Y.,Zhang,Z.,Stenger,V.A.,andWang,Y.Self-CalibratedSpiralSENSE.Magn.Reson.Med.52(2004):688. [57] Qu,P.,Zhong,K.,Zhang,B.,Wang,J.,andShen,G.-X.ConvergenceBehaviorofIterativeSENSEReconstructionwithNon-CartesianTrajectories.Magn.Re-son.Med.54(2005):1040. [58] Rogelj,P.andKovacic,S.Symmetricimageregistration.MedicalImageAnalysis10(2006).3:484. [59] Rudin,L.,Osher,S.,andFatemi,E.Non-linearTotalVariationNoiseRemovalAlgorithm.PhysicsD.60(1992):259. [60] Samsonov,A.andJohnson,C.Non-CartesianPOCSENSE.Proc.Intl.Soc.Mag.Reson.Med..2004,2648. [61] Seiberlich,N.,Breuer,F.A.,Blaimer,M.,Barkauskas,K.,andGriswold,P.M.JakobM.A.Non-CartesianDataReconstructionUsingGRAPPAOperatorGridding(GROG).Magn.Reson.Med.58(2007):1257. [62] Starck,J.,Elad,M.,andDonoho,D.ImageDecompositionviatheCombinationofSparseRepresentationsandaVariationalApproach.IEEETrans.ImageProcess.14(2005):1570. [63] Vogel,C.R.andOman,M.E.Iterativemethodsfortotalvariationdenoising.SIAMJ.Sci.Comput.17(1996):227. [64] Wang,Y.,Yang,J.,Yin,W.,andZhang,Y.ANewAlternatingMinimizationAlgorithmforTotalVariationImageReconstruction.SIAMJ.Imag.Sci.1(2008).3:248. [65] Weickert,Joachim,terHaarRomeny,BartM.,andViergever,MaxA.EfcientandReliableSchemesforNonlinearDiffusionFiltering.IEEETransactionsonImageProcessing7(1998).3:398. [66] Wright,StephenJ.,Nowak,RobertD.,andFigueiredo,M.SparseReconstructionbySeparableApproximation.IEEETrans.SignalProcess.57(2009).7:2479. [67] Yang,Junfeng,Zhang,Yin,andYin,Wotao.AFastTVL1-L2MinimizationAlgorithmforSignalReconstructionfromPartialFourierData.Tech.Rep.08-29,CAAMRiceUniv.,2008. 119

PAGE 120

[68] Ye,J.C.,Tak,S.,Han,Y.,andPark,H.W.ProjectionReconstructionMRImagingUsingFOCUSS.Magn.Reson.Med.57(2007):764. [69] Ye,X.,Chen,Y.,andHuang,F.MRImageReconstructionviaSparseRepresentation:ModelingandAlgorithm.InternationalConferenceonImageProcessing,ComputerVision,andPatternRecognition(IPCV)(2009):10. [70] Yeh,E.N.,Stuber,M.,McKenzie,C.A.,RM,R.M.Botnar,Leiner,T.,Ohliger,M.A.,Grant,A.K.,Willig-Onwuachi,J.D.,andSodickson,D.K.InherentlySelf-CalibratingNon-CartesianParallelImaging.Magn.Reson.Med.54(2005):1. [71] Yeung,Sai-Kit,Tang,Chi-Keung,Shi,Pengcheng,Pluim,JosienP.W.,Viergever,MaxA.,Chung,AlbertC.S.,andShen,HelenC.Enforcingstochasticinverseconsistencyinnon-rigidimageregistrationandmatching.ProceedingsofIEEEConferenceonComputerVisionandPatternRecognition(2008):1. [72] Yeung,S.K.andShi,P.C.Stochasticinverseconsistencyinmedicalimageregistration.InternationalConferenceonMedicalImageComputingandComputer-AssistedIntervention8(2005).2:188. [73] Ying,L.,Liu,B.,Steckner,M.C.,Wu,G.,andM.Wu,S.-J.Li.AStatisticalApproachtoSENSERegularizationWithArbitraryk-SpaceTrajectories.Magn.Re-son.Med.60(2008):414. [74] Zeng,Q.andChen,Y.AccurateInverseConsistentNon-rigidImageRegistrationandItsApplicationonAutomaticRe-contouring.ProceedingsofInternationalSymposiumonBioinformaticsResearchandApplications(2008):293. [75] Zhang,Xiaoqun,Burger,Martin,Bresson,Xavier,andOsher,Stanley.BregmanizedNonlocalRegularizationforDeconvolutionandSparseReconstruction.SIAMJ.Imag.Sci.3(2010).09-03:253. [76] Zhang,Xiaoqun,Burger,Martin,andOsher,Stanley.Auniedprimal-dualalgorithmframeworkbasedonBregmaniteration.Tech.Rep.09-99,CAMUCLA,2009. [77] Zhang,Z.,Jiang,Y.,andTsui,H.Consistentmulti-modalnon-rigidregistrationbasedonavariationalapproach.PatternRecognitionLetters(2006):715. [78] Zhu,M.andChan,T.F.AnEfcientPrimal-DualHybridGradientAlgorithmforTotalVariationImageRestoration.Tech.Rep.08-34,CAMUCLA,2008. [79] Zhu,M.,Wright,S.J.,andChan,T.F.Duality-BasedAlgorithmsforTotal-Variation-RegularizedImageRestoration.Comput.Optim.Appl.47(2010):377. 120

PAGE 121

BIOGRAPHICALSKETCH XiaojingYewasborninBeijing,China,onAugust7th,1983.Hereceivedhisbachelor'sdegreeinpureandappliedmathematicsfromtheSchoolofMathematicalSciencesatPekingUniversity,Beijing,China,inJuly2005.HereceivedtheMasterofStatisticsandDoctorofPhilosophyinMathematicsinDecember2009andMay2011,respectively. 121