University of Southern California Dissertations and Theses

Computational methods for fluorescence molecular tomography

(USC Thesis Other)

Computational methods for fluorescence molecular tomography

PDF

Download

Open document

Flip pages

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content COMPUTATIONALMETHODSFORFLUORESCENCEMOLECULARTOMOGRAPHY
by
JoyitaDutta
ADissertationPresentedto
FACULTYOFTHEUSCGRADUATESCHOOL
UNIVERSITYOFSOUTHERNCALIFORNIA
InPartialFulfillmentofthe
RequirementsfortheDegree
DOCTOROFPHILOSOPHY
(ELECTRICALENGINEERING)
May2011
Copyright 2011 JoyitaDutta
Dedication
toBaba,Ma,andShubhroz..
ii
Acknowledgments
Mygraduatecareerhasbenefitedfromtheexpertiseandadviceofmanyscientists,col-
leagues, friends, and family who have tangibly or intangibly provided inputs into this
doctoralthesis.
The Biomedical Imaging Research Lab at USC, of which I have been a part, has
provided a wonderful environment for transforming ideas into scientific publications.
Thecollaborativespirit,theinteractiveenvironment,andthekeeninterestinbiomedical
research prevailing in this lab have been immensely beneficial to my personal growth
as a researcher. I would like to thank my advisor, Dr. Richard Leahy, whose spirit
of research lends this lab its unique character. Dr. Leahy has been the perfect mentor
who offered me intellectual freedom in research, while keeping himself available any
time I needed his help or advice. This balance between independence and supervision
has helped me tremendously in maturing and evolving as a researcher. I have always
found Dr. Leahy’s outlook on research very inspiring. I am also thankful to him for
always encouraging and supporting my travel to collaborative research meetings and
conferences. This was crucial, given the interdisciplinary nature of my research. His
guidancethroughallphasesofmyPhDworkhasbeeninvaluable.
I would like to highlight the contribution of Dr. Sangtae Ahn, my former colleague
at the Biomedical Imaging Research Lab. I am thankful to him for his valuable inputs
onallaspectsofmyworkandforthemanyenrichingdiscussionsIhavehadwithhim.
iii
ItwasapleasurecollaboratingwiththelabofDr. SimonCherryattheUniversityof
California, Davis. On many occasions, Dr. Cherry’s insight has offered much-needed
direction to this project. I am thankful to Dr. Changqing Li and Dr. Gregory Mitchell
fortheirhardworkwithourprototypicimagingsystemsandforhelpingmeplanexper-
iments and acquire experimental data. I would also like to thank our collaborators at
Cambridge Research and Instrumentation Inc., in particular, Dr. Craig Gardner, Dr.
Dylan Bulseco, and Mr. Russell Gershman for their constant efforts toward coordinat-
ingtheinterdisciplinaryproject.
Iwouldliketo thankmycolleague andcollaboratorDr. AnandJoshi, whoseexper-
tisewithimageregistrationtechniqueshasbeenofgreathelp. IwouldliketothankDr.
Abhijit Chaudhari for helping me familiarize with optical imaging techniques during
my early days in the lab. I am immensely thankful to Dr. Sangeetha Somayajula, Dr.
Quanzheng Li, Dr. Dimitrios Pantazis, Dr. Sanghee Cho, Syed Ashrafulla, Yanguang
Lin, Wentao Zhu, Ran Ren, and all other present and past members of the Biomedical
ImagingResearchLabforbeingwonderfulfriendsandcolleaguesatalltimes.
Ifeelthankfulandfortunatetohavehadsomeverywonderfulfriendsandneighbors
at USC, who made living here a lot of fun and gave me memories to cherish. I would
like to thank all my friends from the Vidushak Improv Comedy Group and the Natural
PathMeditationGroupatUSC,whohelpedmelaughitawaythroughhardertimes.
I would like to thank all my family members, who made growing up fun and living
a pleasure. In particular, I would like to mention my grandmother Nilima Majumder, a
liberalwomanfaraheadofhertimewhosevisionmadeahugedifference. Iwouldlike
to thank my father, Dr. Binay K. Dutta, a chemical engineer and a professor, who is
the most influential teacher in my life. He always took time out of his busy schedule to
tutormethroughschoolandcontinuestooffermeacademicadviceeventoday. Iwould
like to thank my mother, Dr. Ratna Dutta, a professor and a Sanskrit scholar, who also
iv
tutoredmeduringmyschooldays,whoputmyeducationaheadofeverythingelse,and
whose pragmatism and resilience I still continue to depend on. Finally, I am thankful
to my husband, Shubhroz Gill, a biologist and an enthusiastic researcher, who offered
crucial inputs to all aspects of my work, provided me with much-needed support and
encouragementatalltimes,andhasalwaysremindedmetoaimhigh.
v
TableofContents
Dedication ii
Acknowledgments iii
ListofTables ix
ListofFigures x
Abstract xv
Chapter1: Introduction 1
1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 OrganizationofthisDissertation . . . . . . . . . . . . . . . . . . . . . 3
Chapter2: Overview 5
2.1 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.1.1 Surfaceprofiling . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.1.2 Fluorescencedataacquisition . . . . . . . . . . . . . . . . . . . 6
2.1.3 Atlaswarping . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.1.4 Datamapping . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.1.5 Forwardmodelcomputation . . . . . . . . . . . . . . . . . . . 7
2.1.6 Reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.2 ForwardProblem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.3 InverseProblem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
Chapter3: FastForwardSolver 14
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.2.1 Theory. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.2.2 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.3 ForwardModelComputation . . . . . . . . . . . . . . . . . . . . . . . 17
3.3.1 Correlationanalysis . . . . . . . . . . . . . . . . . . . . . . . . 17
3.3.2 Efficiencycomparison . . . . . . . . . . . . . . . . . . . . . . 17
vi
3.4 PhantomExperiment . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
Chapter4: IlluminationPatternOptimization 23
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
4.2 FormulationoftheOptimizationProblem . . . . . . . . . . . . . . . . 24
4.2.1 Fisherinformationmatrix . . . . . . . . . . . . . . . . . . . . . 25
4.2.2 Optimizationproblemformulation . . . . . . . . . . . . . . . . 26
4.2.3 Dimensionalityreduction . . . . . . . . . . . . . . . . . . . . . 28
4.2.4 Modifiedoptimizationproblem . . . . . . . . . . . . . . . . . . 29
4.2.5 Extensiontomultispectraldetection . . . . . . . . . . . . . . . 30
4.3 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
4.3.1 Computationoftheforwardmodel . . . . . . . . . . . . . . . . 31
4.3.2 Optimizationprocedure . . . . . . . . . . . . . . . . . . . . . . 32
4.3.3 Inversesolutionandperformancemetrics . . . . . . . . . . . . 34
4.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
4.4.1 Patternsonthecylinder . . . . . . . . . . . . . . . . . . . . . . 38
4.4.2 PatternsontheDigimouseatlas . . . . . . . . . . . . . . . . . 38
4.4.3 Pointsourcelocalization . . . . . . . . . . . . . . . . . . . . . 38
4.4.4 Performanceevaluation . . . . . . . . . . . . . . . . . . . . . . 40
4.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
4.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
Chapter5: ImagingSystems 44
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
5.2 PyramidalMirrorSetupwithPatternedIllumination . . . . . . . . . . . 44
5.3 ConicalMirrorSetup . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
5.3.1 Systemdescription . . . . . . . . . . . . . . . . . . . . . . . . 46
5.3.2 Mappingcalibration. . . . . . . . . . . . . . . . . . . . . . . . 47
Chapter6: PhantomStudies 53
6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
6.2 CubicPhantom . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
6.2.1 Phantomdescription . . . . . . . . . . . . . . . . . . . . . . . 53
6.2.2 Dataacquisition . . . . . . . . . . . . . . . . . . . . . . . . . . 54
6.2.3 Forwardmodelcomputationandimagereconstruction . . . . . 57
6.3 Mouse-ShapedPhantom . . . . . . . . . . . . . . . . . . . . . . . . . 58
6.3.1 Phantomdescription . . . . . . . . . . . . . . . . . . . . . . . 58
6.3.2 Dataacquisition . . . . . . . . . . . . . . . . . . . . . . . . . . 59
6.3.3 SurfacewarpingusingDigiWarp . . . . . . . . . . . . . . . . . 60
6.3.4 Forwardmodelcomputationandimagereconstruction . . . . . 61
vii
Chapter7: RegularizationUsingJointL
1
andTotalVariationNormPenal-
ties 64
7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
7.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
7.2.1 PreconditionedConjugateGradientalgorithm . . . . . . . . . . 68
7.2.2 SeparableParaboloidalSurrogatesalgorithm. . . . . . . . . . . 68
7.2.3 Compoundapproach . . . . . . . . . . . . . . . . . . . . . . . 74
7.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
7.3.1 Experimentalsetup . . . . . . . . . . . . . . . . . . . . . . . . 75
7.3.2 Simulationresults . . . . . . . . . . . . . . . . . . . . . . . . . 76
7.3.3 Experimentalresults . . . . . . . . . . . . . . . . . . . . . . . 78
7.4 SummaryandDiscussion . . . . . . . . . . . . . . . . . . . . . . . . . 79
Chapter8: In VivoStudies 84
8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
8.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
8.2.1 Experimentalsetup . . . . . . . . . . . . . . . . . . . . . . . . 84
8.2.2 Dataacquisition . . . . . . . . . . . . . . . . . . . . . . . . . . 85
8.2.3 Forwardmodelingandimagereconstruction . . . . . . . . . . . 85
8.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
Chapter9: SummaryandFutureWork 89
9.1 SummaryofContributions . . . . . . . . . . . . . . . . . . . . . . . . 89
9.1.1 Fastforwardsolver . . . . . . . . . . . . . . . . . . . . . . . . 89
9.1.2 Illuminationpatternoptimization . . . . . . . . . . . . . . . . . 90
9.1.3 RegularizationusingL
1
andTVnormpenalties . . . . . . . . . 90
9.1.4 Imagingsystemdevelopmentandcalibration . . . . . . . . . . 91
9.1.5 Phantomandinvivostudies . . . . . . . . . . . . . . . . . . . 91
9.2 FutureWork . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
9.2.1 Theoreticalgoals . . . . . . . . . . . . . . . . . . . . . . . . . 92
9.2.2 Experimentalgoals . . . . . . . . . . . . . . . . . . . . . . . . 93
Bibliography 94
viii
ListofTables
3.1 Comparisonofaccuracyandefficiencyofdifferentmonochromaticemis-
sionforwardmodels . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
4.1 FullwidthathalfmaximumofpointspreadfunctionsinFigure4.4 . . . 39
4.2 Conditionnumberofthesystemmatrixfordifferentilluminationschemes 40
7.1 Performance metrics for comparison of theL
2
,L
1
, TV,L
1
-TV penalty
functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
ix
ListofFigures
2.1 CompleteprocesschainforFMT. . . . . . . . . . . . . . . . . . . . . . 5
3.1 Comparison of correlation coefficients. The correlation coefficients of
the photon density at each detector due to all the source points are dis-
playedfor(a)ahomogeneoustissuemodelandFEM(b)ourperturbative
modelandFEM. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.2 Comparison of the increase in computational costs of our perturbative
forward model and the FEM-based forward model with increase in the
numberofdetectornodes. . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.3 (a)Schematicoftherectangularphantomshowingthefrontlayerthatis
4timesmorescatteringthantheremainingbulk. (b)The15pointsource
locationsonthefrontsurfaceusedforexcitation. (c)Apyramidalmirror
setupthatprovidestop,bottom,andsideviews. . . . . . . . . . . . . . 20
3.4 Horizontalsectionsofthe32mm×32mm×29mmphantomshowing
(a)theoriginalsourcepositionandthereconstructedsourcedistributions
obtainedusing(b)ourperturbativeforwardmodeland(c)anFEM-based
forwardmodel. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
4.1 (a) Front and back views of the curved surface of a tessellated cylin-
der (flattened for display purposes with the cylinder axis aligned hor-
izontally) and (b) top and bottom views of the Digimouse atlas show-
ing eigenfunctions corresponding to the 20 smallest eigenvalues of the
Laplace-Beltramioperator. Thebasisfunctionsareorthogonaland,there-
fore,havebothpositiveandnegativevaluesatvariouspoints. . . . . . . 33
4.2 A set of 3 optimal patterns on a tessellated cylinder. The individual
patterns are displayed in (a), (b), and (c). The sum of the 3 patterns
displayed in (d) shows rough overall uniformity of illumination. For
ease of display, the curved surface of the cylinder is laid out flat so that
the vertical axis represents the height of the cylinder and the horizontal
axisrepresentsthetransverseangle. . . . . . . . . . . . . . . . . . . . 36
x
4.3 ProjectionsofthetopandbottomsurfacesoftheDigimouseatlasshow-
ingoptimalpatternsfor(a)p = 1,(b)p = 3,(c)p = 6,and(d)p = 12. . 37
4.4 Point spread functions for (a) optimal patterns forp = 12, (b) optimal
patterns forp = 6, (c) optimal patterns forp = 3, (d) uniform patterns
forp = 3,and(e)optimalpatternforp = 1. . . . . . . . . . . . . . . . 39
4.5 Singularvaluedistributionsfordifferentilluminationschemes. . . . . . 40
4.6 Averageresolutionvs. variancecurvesfordifferentilluminationschemes. 41
5.1 (a)Acuboidalphantomplacedonthetransparentstageinsidethepyra-
midal mirror setup. (b) CCD camera image of a cubic phantom. The
top,bottom,left,andrightviewsareasindicated. . . . . . . . . . . . . 45
5.2 (a) A cubic phantom placed on the transparent stage inside the conical
mirror setup. (b) CCD camera image of a cubic phantom. (c) CCD
cameraimageofamouse. . . . . . . . . . . . . . . . . . . . . . . . . . 46
5.3 Radial section of the setup at an angleϕ = ϕ
p
showing an object point
P,(r
p
,ϕ
p
,z
p
)andthecorrespondingimagepointQ,(r
i
,ϕ
i
,z
i
). . . . . . 48
5.4 CalibrationparametersR,D,δθ,θ
xc
,θ
yc
,θ
xs
,andθ
zs
. . . . . . . . . . . 50
5.5 Corner-pointsinthecheckeredgridonthecubesurfaceusedforfeature
matching. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
5.6 (a) Photograph of the original checkered cube, (b) CCD camera image,
and(c)digitalrenderingofthecubewiththereconstructedcheckerboard
pattern.. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
6.1 (a) White-light image showing simultaneous views of four sides of the
cube phantom. (b) Log-scale CCD image of the fluorescence pattern
at a detection wavelength of 710 nm due to the embedded source. A
6.5 x 6.5 mm square patch of the bottom surface was illuminated at an
excitationwavelengthof635nm. . . . . . . . . . . . . . . . . . . . . . 54
6.2 Set of 17 marching square patterns illuminating the bottom side of the
cube. Thefrontandbackfaceshavebeenremovedfromthecubesurface
foreaseofviewing. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
6.3 Set of fluorescence patterns corresponding to the 17 patterns in figure
6.2 for the 635/710 nm excitation/emission pair. The front and back
faceshavebeenremovedfromthecubesurfaceforeaseofviewing. . . 56
xi
6.4 Reconstructed image displayed as (a) an isosurface using a 99th per-
centileintensitythresholdfortheimagevolumeand(b)transverseslices
throughtheaxisofthefluorescenttube. . . . . . . . . . . . . . . . . . 58
6.5 Mouse-shaped phantom. (a) Full-light photograph. (b) A coronal and a
transverse section from the CT image showing the twotubes filled with
fluorescentDiDdye. . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
6.6 The three steps of the DigiWarp procedure: posture correction, surface
warping by means of asymmetric L
2
distance minimization, and volu-
metricregistrationbasedonelasticenergyminimization. . . . . . . . . 60
6.7 Mouse-shaped phantom. (a) Transformation of the Digimouse surface
usingtheDigiWarpprocedure(b)Photographofthemouse-shapedphan-
tomshowingitssurfacegeometryforreference. . . . . . . . . . . . . . 61
6.8 Coronal sections of the mouse phantom showing the reconstruction of
two embedded line sources using the L
2
penalty with a regularization
parameter of 10
−4
. The slices are ordered from bottom to top starting
topleftintheimageandendingbottomright. . . . . . . . . . . . . . . 62
6.9 (a) Coronal and (b) transverse sections of the CT image of the mouse-
shapedphantomshowingthetwoembeddedfluorescentlinesources. (c)
Coronaland(d)transverseoverlayofCTandFMTimages. (e)Coronal
and (f) transverse sections of the FMT image showing the two fluores-
centlinesourcesreconstructedusingtheL
2
penaltywitharegularization
parameterof10
−4
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
7.1 Comparison of convergence curves of PCG, OSSPS, and OSSPS-PCG
with the transition from OSSPS to PCG occurring after iteration num-
bers 4, 6, 8, 10, and 12. These curves were obtained for the L
1
-TV
penaltyandforarandominitialization. . . . . . . . . . . . . . . . . . . 74
7.2 Mouse-shaped phantom. (a) Full-light photograph. (b) Coronal and
transversesectionsfromtheCTimageshowingthetwotubesfilledwith
fluorescentDiDdyeinsidethephantom. . . . . . . . . . . . . . . . . . 76
7.3 Coronalsectionsofthemousephantomshowingthetwosimulatedline
sources. The slices are ordered from bottom to top starting top left in
theimageandendingbottomright. . . . . . . . . . . . . . . . . . . . . 77
7.4 Coronalsectionsofthemousephantomshowingthetwosimulatedline
sourcesreconstructedusingtheL
2
penaltywitharegularizationparam-
eterof10
−3
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
xii
7.5 Coronalsectionsofthemousephantomshowingthetwosimulatedline
sourcesreconstructedusingtheL
1
penaltywitharegularizationparam-
eterof10
−4
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
7.6 Coronalsectionsofthemousephantomshowingthetwosimulatedline
sourcesreconstructedusingtheTVpenaltywitharegularizationparam-
eterof2×10
−6
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
7.7 Coronalsectionsofthemousephantomshowingthetwosimulatedline
sources reconstructed using bothL
1
and TV penalties with regulariza-
tionparametersof5×10
−5
and2×10
−6
respectively. . . . . . . . . . 79
7.8 Coronal sections of the mouse phantom showing the reconstruction of
two embedded line sources using the L
2
penalty with a regularization
parameterof10
−4
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
7.9 Coronal sections of the mouse phantom showing the reconstruction of
two embedded line sources theL
1
penalty with a regularization param-
eterof20. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
7.10 Coronal sections of the mouse phantom showing the reconstruction of
two embedded line sources using the TV penalty with a regularization
parameterof0.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
7.11 Coronal sections of the mouse phantom showing the reconstruction of
two embedded line sources using bothL
1
and TV penalties with regu-
larizationparametersof10and1respectively. . . . . . . . . . . . . . . 82
7.12 (a) Coronal and (b) transverse sections of the CT image of the mouse-
shapedphantomshowingthetwoembeddedfluorescentlinesources. (c)
Coronaland(d)transverseoverlayofCTandFMTimages. (e)Coronal
and (f) transverse sections of the FMT image showing the two fluores-
cent line sources reconstructed using both L
1
and TV penalties with
regularizationparametersof10and0.7respectively. . . . . . . . . . . . 83
8.1 Full-light photograph showing the anesthetized mouse placed on the
stageinsidetheconicalmirrorsystem. Thedetectornodesonthewarped
mouse atlas mapped onto the conical mirror surface are indicated in
green. The illumination locations mapped onto the mirror surface are
showninred. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
xiii
8.2 (a) Coronal and (b) transverse sections of a CT image of the mouse
showingthetumorbulgingontherightflank. (c)Coronaland(d)trans-
verse overlay of CT and FMT images. (e) Coronal and (f) transverse
sections of the FMT image showing the tumor on the right flank and
fluorescencefromthebladder. . . . . . . . . . . . . . . . . . . . . . . 88
xiv
Abstract
Fluorescence molecular tomography is an imaging modality which exploits the speci-
ficityoffluorescentbiomarkerstogeneratevolumetricimagesusingnear-infraredlight,
which is safe and non-ionizing, at a substantially lower cost than competing modalities
like positron emission tomography (PET) or single photon emission computed tomog-
raphy (SPECT). While these very attractive features make it well-suited for preclinical
research,the3Dreconstructionoffluorescentsourcesisconfoundedbythehighdegree
of absorption and scattering of photons propagating through tissue, which make the
inverseproblemill-posed. ThefocusofmyPh.D.researchhasbeenthedevelopmentof
computationaltechniquestotacklethischallengeinordertoreconstructhighresolution
3D images of molecular targets in vivo in small animals. To accurately and efficiently
predict the photon fluence on the animal surface, we have developed a fast approach
based on the Born approximation for computing the photon propagation model with
reasonable accuracy. We have also developed a systematic approach for designing effi-
cientspatial patternsof illumination forfluorescence tomography. Thekeyideabehind
our approach is to formulate a constrained optimization problem to improve the con-
ditioning of the Fisher information matrix for the fluorescence forward model. The
ill-conditioned nature of the inverse problem warrants the incorporation of some prior
information to improve the conditioning of the problem and make it solvable. We have
xv
developed a regularization approach that incorporates this information as a combina-
tion of sparsity- and smoothness-enforcing penalty functions. In collaboration with Dr.
Simon Cherry’s lab at UC Davis, and Cambridge Research and Instrumentation Inc.,
we have designed and built two fluorescence imaging systems capable of generating
a variety of illumination patterns and with optics and filters designed for 360
◦
surface
viewingandmultispectraldatacollection. Weconductedstudiesusingtissuephantoms
with embedded fluorescent sources. We also performed an in vivo study on an anes-
thetized mouse with a xenograft melanoma and validated it by overlaying it against the
CTimageofthesameanimal.
xvi
Chapter1
Introduction
1.1 Background
Opticalimagingtechniqueshavelongbeenpopularforgeneratingcontrastrepresenting
molecular processes in biological tissue. Over the past decade, their applications have
extended from microscopic to macroscopic imaging of deep-tissue optical sources by
exploitingthenear-infrared(NIR)windowwherewater,oxyhemoglobin,anddeoxyhe-
moglobin, the primary absorbers in tissue, have relatively low absorption coefficients
allowing light to penetrate several centimeters inside tissue [WN03]. The non-ionizing
nature of NIR light, the availability of a variety of highly specific fluorescent probes,
and,finally,lowcostofferthesemethodssomeleverageoverexistingradiologicaltech-
niques for molecular imaging [HAD97]. However, a penetration depth of only a few
centimeters is insufficient for most clinical deep-tissue imaging applications. Thus the
success of optical imaging methods in clinical diagnostics is restricted to only a few
applications. These include NIR spectroscopy and diffuse optical tomography used for
mappingbrainstructureandfunctioninneonates[KYM
+
03,VC97,HGY
+
02]andopti-
cal mammography for breast cancer screening [NC01, CCH
+
03]. In comparison, in
vivodeep-tissueopticalimagingofsmallanimalsrequiresphotonstotravelshorterdis-
tances leading to higher photon counts and better signal-to-noise ratio (SNR). Hence
optical tomographic methods show great promise in preclinical research, which is a
valuable translational tool between in vitro studies and clinical applications. Fluo-
rescence molecular tomography (FMT) and bioluminescence tomography (BLT) have
1
emerged as promising low-cost alternatives to positron emission tomography (PET)
and single photon emission computed tomography (SPECT) for functional imaging
in small animals, thus greatly impacting diagnostics, drug discovery, and therapeutics
[NRWW05,GHA05].
TheobjectiveofFMTandBLTistocomputethe3Ddistributionofaphotonsource
inside an animal volume from the photon density detected on the surface of the ani-
mal. Thephotonsourceistypicallyabioluminescentorafluorescentprobetaggingthe
molecule of interest. Unlike bioluminescence, which is triggered chemically, fluores-
cence involves absorption and emission of light with a downward shift in frequency
and thus requires an external light source for excitation. Although plagued by tis-
sue autofluorescence, which tends to reduce the signal-to-background ratio affecting
probe specificity [TJMSR04], FMT has many advantages over BLT. Compared to most
bioluminescent probes, commonly used fluorophores emit light at longer wavelengths
(wheretissuesarelessabsorbing)consequentlyofferinghigherdetectedsignalstrengths
[CB02,Hie05]. Additionallymultipleilluminationpatternscanbeusedtogeneratedif-
ferent mappings from the source space to the detector space making the FMT problem
lessill-posedthantheBLTproblem[CDB
+
05].
ThesuccessofFMTcanbeattributedtoanumberofbreakthroughs.
• The availability of a variety of new NIR fluorescent dyes, active or activatable
fluorescent biomarkers, and fluorescent proteins expressed by reporter genes has
enabledvisualizationofgeneexpressionandseveralcellularandsubcellularpro-
cessesinvivo[MG03,SRL
+
09].
• Advances in instrumentation have led to the development of a range of imaging
systems for time domain, frequency domain, and continuous wave (CW) FMT
[KRD
+
08,GSME05,ZVM
+
06]. Someofthenewersystemsfeaturenon-contact
tomographic detection using CCD cameras with free-space detection geometries
2
and/or innovative optics for full-surface visualization, thus eliminating the need
foropticalfibersandmatchingfluids[LMD
+
09b,GRWN03,PBAC05].
• Finally,anumberoftheoreticalandcomputationaladvanceshaveledtothedevel-
opmentofrealisticforwardmodelsandrobustinversemethods. Popularmethods
employed to model photon propagation through tissue for solving the forward
problemincludeanalyticalandnumericalsolutionstotheradiativetransportequa-
tion[KNBH02]andthediffusionequation[RCN01,DAL
+
08b,ASHD93]subject
to different boundary conditions [HST
+
94]. These forward modeling schemes
coupled with fast inversion techniques [RSM01, ACD
+
08, ZSA
+
09] have made
itfeasibletoreconstructFMTimagesaccuratelyandefficiently.
Despite their tremendous potential and increasing popularity, fluorescence tomo-
graphictechniquesareconfoundedbyhighdegreesofabsorptionandscatteringofpho-
tonspropagatingthroughtissue,makingtheFMTproblemill-posed. Inthisdissertation,
we will explore various computational methods aimed toward alleviating this problem
andimprovingsourcelocalization.
1.2 OrganizationofthisDissertation
Theunderlyingtheoreticalconceptsandexperimentalproceduresaredescribedinchap-
ter 2. Toward the end of this chapter, we describe some schemes for handling large
volumes of FMT data. In chapter 3, we present a fast, perturbative approach based
of the Born approximation to solve the forward problem of FMT for a realistic ani-
mal geometry and optically heterogeneous tissue types. Chapter 4 addresses one of
the most intriguing problems of FMT, which is to determine how to optimally illumi-
nate the animal surface so as to maximize the information content in the acquired data.
3
Thekeyideahereistoparameterizetheilluminationpatternsandtosolveanoptimiza-
tion problem that seeks to improve the conditioning of the Fisher information matrix.
In collaboration with Dr. Simon Cherry’s lab at University of California, Davis, and
Cambridge Research and Instrumentation Inc., we have built two fluorescence imaging
systemscapableofgeneratingavarietyofilluminationpatternsandwithopticsandfil-
ters designed for 360
◦
surface viewing and multispectral data collection. In chapter 5,
we present details of instrumentation and calibration. In chapter 6, we present recon-
struction results for preliminary studies on tissue phantoms with embedded fluorescent
tubesperformedusingthetwosystems. Owingtotheill-conditionednatureoftheFMT
inverse problem, image reconstruction approaches rely on the incorporation of prior
informationintheformofsuitableregularizers. Wederiveinspirationfromthefactthat
very often fluorescent probes and contrast agents are concentrated within specific areas
ofinterest(e.g.,insidelocalizedtumorsandexcretoryorgans)andpresentaregulariza-
tion approach in chapter 7 that uses a combination of L
1
and total variation penalties
in the FMT inverse problem to simultaneously enforce sparsity and smoothness in our
reconstructedimages. Inchapter8,wepresentreconstructionresultsfor in vivostudies
performed using one of the systems described in chapter 5 on an anesthetized mouse
with a xenograft melanoma. Chapter 9 presents a summary of all these contributions
andoutlinesresearchdirectionsforthefuture.
4
Chapter2
Overview
2.1 Methodology
TheoverallprocessflowofFMTissummarizedinfigure2.1. Steps1and2arethekey
experimental steps, while steps 3 to 6 are computational. These component tasks are
describedbelow:
2.1.1 Surfaceprofiling
Priorknowledgeoftheanimalsurfacegeometryandtissueopticalpropertiesisrequired
to solve the FMT forward problem. Anatomical imaging modalities such as CT or
MRI provide information about surface topography and internal anatomy. But, in
Figure2.1: CompleteprocesschainforFMT.
5
order to perform an anatomical scan, the anesthetized animal needs to be transferred
from the FMT system to a CT or MRI scanner. To avoid this cumbersome step, a
variety of all-optical surface topography setups have been integrated with FMT sys-
tems. Examples include photogrammetric techniques [RSN03], structured light based
schemes [RXK06], shadow photogrammetric techniques [MGZ
+
07], and 3D volume
carving [LSRN08]. However, estimation of the internal anatomy using all-optical tech-
niquesisnon-trivial. Thisnecessitatesstep3offigure2.1.
2.1.2 Fluorescencedataacquisition
In non-contact FMT, CCD cameras are used for data collection. The CCD camera by
itselfgeneratesa2Dplanarimagethatprovidesapartialviewofthemousesurface. The
FMTinverseproblemisextremelyill-conditionedforsuchpartialsurfaceviews. There-
fore, for tomographic reconstruction, a variety of optical elements are used to place the
complete animal surface within the field-of-view of the CCD camera. The illumination
anddetectionwavelengthsandspatialpatternsofilluminationarekeydesignparameters
forthisstepandarealsorequiredtogeneratetheforwardmodelinstep5offigure2.1.
2.1.3 Atlaswarping
For accurately solving the FMT forward problem, we require knowledge of the animal
geometryandtissueopticalproperties. Asmentionedbefore,thisinformationisdifficult
to extract using all-optical methods. An alternative approach is to warp an anatomical
atlas to the surface profile data and assign optical properties to the atlas nodes organ-
wisebasedonpublishedvalues.
6
2.1.4 Datamapping
Themappingfromthe3Dobjectspacetothe2DCCDimageplanedependsontheoptics
usedforfull-surfaceviewing. Thismappingestablishesthecorrespondencebetweenthe
detector nodes on the warped atlas and the pixel locations in the CCD camera images.
Using this mapping information, the fluorescence data can be mapped back from the
CCDcameraimagestotheatlasdetectornodes.
2.1.5 Forwardmodelcomputation
The FMT forward model maps a unit fluorescent source at any 3D location inside the
animal volume to a surface fluorescence pattern. The solution to the forward prob-
lem requires prior knowledge of the surface geometry and the tissue optical properties.
Hence,itmustbeprecededbystep3. Also,thesolutiondependsontheilluminationand
detection wavelengths and the spatial patterns of illumination. This step is discussed in
furtherdetailinsection2.2.
2.1.6 Reconstruction
ThegoaloftheFMTinverseproblemistoestimatethe3Ddistributionofafluorescent
source inside the animal volume. The FMT inverse problem is linear as described in
2.3. Oncetheforwardmodelmatrixandthedatavectorshavebeencomputed,avariety
oflinearinversionschemescanbeappliedtoreconstructthe3Dsourcedistribution.
2.2 ForwardProblem
In continuous wave fluorescence molecular tomography (CW FMT), the 3D biodistri-
bution of the fluorophore inside an animal volume is computed from the steady-state
7
surface fluorescence photon density measured using a CCD camera. With the diffu-
sion approximation, the CW FMT problem can be described by a set of coupled partial
differentialequations(PDEs)[HLSM95,MSO
+
04,PP06]:
[∇·(−κ(r,λ
ex
)∇)+µ
a
(r,λ
ex
)]Φ(r,λ
ex
) =w(r,λ
ex
), (2.1)
[∇·(−κ(r,λ
em
)∇)+µ
a
(r,λ
em
)]Φ(r,λ
em
) =η(λ
em
)q(r)Φ(r,λ
ex
). (2.2)
Hereκ(r,λ) = 1/[3(µ
a
(r,λ)+µ
′
s
(r,λ))]isthediffusioncoefficientandµ
a
(r,λ)and
µ
′
s
(r,λ)respectivelyaretheabsorptionandreducedscatteringcoefficientsataposition
randawavelengthλ. Forasurfaceilluminationpatternw(r,λ
ex
)atanexcitationwave-
lengthλ
ex
,thecorrespondingexcitationfield,Φ(r,λ
ex
),overananimalvolume,Ω,can
be computed by solving (2.1). For a known fluorophore distributionq(r) and emission
spectral coefficientη(λ
em
), the emission field, Φ(r,λ
em
), can be computed by solving
(2.2)forthepreviouslycomputedexcitationfield. InFMT,weassumepriorknowledge
oftheanimalsurfacegeometryandtissueopticalproperties,µ
a
(r,λ)andµ
′
s
(r,λ). The
diffusion approximation is valid under the condition µ
′
s
≫ µ
a
, an assumption that is
valid for most tissue types at near-infrared wavelengths. Appropriate boundary con-
ditions must be imposed while solving (2.1) and (2.2) [HST
+
94, Aro95]. The Robin
boundary condition, for example, requires that the total inwardly-directed photon cur-
rentattheboundarybezero:
Φ(r,λ
ex
)+2Gˆ ν·∇(κ(r,λ
ex
)Φ(r,λ
ex
)) = 0 (2.3)
Φ(r,λ
em
)+2Gˆ ν·∇(κ(r,λ
em
)Φ(r,λ
em
)) = 0. (2.4)
Here ˆ ν denotes the outward unit normal at position r on the boundary, ∂Ω, while G
dependsontherefractiveindexmismatchattheboundary[SAHD95].
8
ThePDEs(2.1)and(2.2)canbedecoupledandsolvedsubjecttoconditions(2.3)and
(2.4)respectivelybyreplacingthesourcetermsontherighthandsidesof(2.1)and(2.2)
bypointsources(deltafunctions)atdifferentlocations. Inadiscretizeddomainwithn
d
detectornodesontheanimalsurfaceandn
s
pointsourcelocationsdistributedinsidethe
volume, the excitationforwardmodel at an excitationwavelengthλ
ex
can be expressed
asamatrixA
ex
∈R
ns×n
d
,withthejthcolumnofthismatrixrepresentingtheexcitation
fieldcausedbyilluminatingthejthsurfacenode. Similarly,theemissionforwardmodel
at an emission wavelengthλ
em
can be expressed as a matrixA
em
∈ R
n
d
×ns
. Thejth
columnofthismatrixrepresentsthesurfacefluorescencepatternduetothejthinternal
pointsource. Forasetofpilluminationpatterns,wediscretizetheintensitydistribution
atthesurface,w(r,λ
ex
),asasetofvectorsw
k
∈R
n
d
,wherek = 1 :p. Thesolutionto
theforwardproblemisasystemmatrixdependentonA
ex
,A
em
,andallthew
k
vectors.
Theblockofthesystemmatrixcorrespondingtothekthilluminationpatterncanbe
obtainedbydiagonallyscalingtheemissionforwardmatrixbytheexcitationintensities
ateachinternalpoint:
A
k
=A
em
D
ex
k
. (2.5)
HereD
ex
k
∈R
ns×ns
isadiagonalmatrixrepresentingtheexcitationfieldduetow
k
and
isgivenby:
D
ex
k
= diag
i
(
d
i
k
)
, (2.6)
whered
i
k
istheithcomponentofvectord
k
computedfrom:
A
ex
w
k
=d
k
. (2.7)
Many FMT systems harness spectral variations of tissue optical properties to
improve depth resolution and source localization by using multispectral illumination
and/or detection [ZKS
+
05, CDB
+
05, CAL
+
09]. We will now look at how the block of
9
thesystemmatrixcorrespondingtothekthilluminationpatterngetsmodifiedwhenmul-
tispectraldetectionisdone. Weassumethat,foreachilluminationpattern,fluorescence
data is collected over s wavelength bins. We represent the monochromatic emission
forwardmodelatthelthwavelengthbinasA
em
l
. Thefluorophoreisassumedtohavean
emissionspectralcoefficientη
l
atthelthwavelengthbin. Then,theblockofthesystem
matrixcorrespondingtothekthilluminationpatternisgivenby:
A
k
=









η
1
A
em
1
D
ex
k
η
2
A
em
2
D
ex
k
.
.
.
η
s
A
em
s
D
ex
k









. (2.8)
The full system matrix,A ∈ R
psn
d
×ns
, is obtained by concatenating the individual
forward model matrices for a set of p different illumination patterns and s detection
wavelengths:
A =
[
A
′
1
A
′
2
... A
′
p
]
′
. (2.9)
2.3 InverseProblem
The goal of the inverse problem is to compute the unknown fluorophore distribution,
q∈R
ns
,bysolvingthefollowinglinearsystemofequations:
Aq =b. (2.10)
10
HereAisthesystemmatrixcomputedfrom(2.9)andb∈R
psn
d
isthemeasuredsteady-
statesurfacefluorescencedatacorrespondingtopdifferentilluminationpatternsstacked
upasasingledatavector,suchthat:
b =
[
b
′
1
b
′
2
... b
′
p
]
′
, (2.11)
whereb
k
∈R
sn
d
isthemultispectraldatavectorforthekthilluminationpattern.
The system matrixA in FMT is rank-deficient. There exist an infinite number of
solutions. Different constraints and penalties can be used to introduce additional infor-
mation and restrict ourselves to a particular solution-set. One common approach is to
useanL
2
-normpenaltythatgeneratestheminimumenergysolution. Theimagerecon-
struction problem then is an optimization problem that seeks to compute the following
regularizedleastsquaressolution:
^ q = argmin
q
(
1
2
||Aq−b||
2
+
α
2
||q||
2
)
. (2.12)
Here ^ q is the reconstructed 3D fluorescent source distribution and α is the Tikhonov
regularization parameter [TA77]. The first term in (2.12) is a data-fitting term, while
the second term is a regularization penalty term that smooths the reconstructed image.
Increasingα produces smoother and less noisy images but also tends to increase local-
izationerror.
11
Onewaytosolvetheproblemin(2.12)istousesingularvaluedecomposition. Given
the singular value decompositionA = UV
′
of the system matrix, this problem can
besolvedusingtheregularizedpseudoinverse:
^ q = (A
′
A+αI)
−1
A
′
b
= V (αI +
2
)
−1
U
′
b.
(2.13)
Alternatively, a variety of iterative methods can be used. One advantage of using itera-
tivemethodsisthattheyallowustoimposeanon-negativityconstraintonthesolution.
The system matrix,A, is typically very large in size, with many more rows (psn
d
)
than columns (n
s
). It is often infeasible to store this large matrix in computer memory.
Instead of solving the inverse problem in 2.10, we can solve the equivalent problem
basedonweightedleastsquares:
A
′
RAq =A
′
Rb. (2.14)
The matrixA
′
RA∈R
ns×ns
is of the same size irrespective of the number of patterns,
p, or the number of detection wavelengths, s. The matrixR is a diagonal weighting
matrix obtained by inverting the covariance matrix and can be written in terms of itsp
diagonalblockscorrespondingtotheppatterns.
R =









R
1
0 ··· 0
0 R
2
··· 0
.
.
.
.
.
.
.
.
.
.
.
.
0 0 ··· R
p









(2.15)
12
TheA
′
RAmatrixcanbecomputedasfollows:
A
′
RA =
∑
k
D
ex
k
(A
em
)
′
R
k
A
em
D
ex
k
. (2.16)
HereR
k
is thekth diagonal block ofR. The data vector is computed using equation
2.17below,where◦representsanentry-wisematrixproduct:
A
′
Rb =
∑
k
D
ex
k
(A
em
)
′
R
k
b
k
=
∑
k
d
k
◦((A
em
)
′
R
k
b
k
). (2.17)
Theapproachdescribedin(2.15),(2.16),and(2.17)willbeusedwhereverthedata-sets
aresolargethatitisimpossibletocomputethefullsystemmatrix,A.
13
Chapter3
FastForwardSolver
3.1 Introduction
As described in section 2.2, the FMT forward problem seeks to determine, for a given
3D source distribution, the photon density on the surface of an animal. Photon trans-
port through a turbid medium in which scattering is much stronger than absorption can
be modeled by applying the diffusion approximation to the radiative transport equation
[Arr99]. The key challenge is to accurately and efficiently solve the diffusion equation
for a realistic animal geometry and heterogeneous tissue types. Fast analytical solvers
are available that can be applied to arbitrary geometries but assume homogeneity of
tissue optical properties and hence have limited accuracy [RCN01, GRWN03, SN06].
These works assume simplified semi-infinite geometries and use appropriate boundary
conditions which prevent the violation of the diffusion approximation near the surface
[HST
+
94,Aro95,RNCNV01,RNVWN02]. Amoreaccurateapproachfortacklingthe
forward problem which allows reasonably accurate modeling of both animal geometry
and tissue heterogeneity is to use the finite element method (FEM) with volume tessel-
lation [ASHD93, RSM01, CDB
+
05]. But this approach is computationally intensive.
Thecomputationalchallengeisheightenedwhenoneisworkingwithmultispectraldata
to improve source localization and conditioning of the inverse problem. Thus, while
solving the forward problem, we are faced with a trade-off between accuracy and effi-
ciency. We present here a fast forward model based on the Born approximation that
falls in between these two approaches. Our model introduces tissue heterogeneity as
14
perturbations in diffusion and absorption coefficients at rectangular grid points inside a
mouseatlas.
The Born approximation is a well-known perturbation technique that can be used
to generate an analytical solution to the diffusion equation without ignoring tissue het-
erogeneity and has been used in DOT [SZ05] and BLT [CDWW06]. We present here
a perturbative forward modeling scheme based on the Born approximation that takes
into account both animal geometry and tissue heterogeneity and, at the same time, is
computationallyefficient.
3.2 Methods
3.2.1 Theory
The photon density at each surface detector node due to each unit point source can be
represented by a Green’s function, which, for a homogeneous tissue model and a real-
istic animal geometry, can be analytically computed using a semi-infinite slab approx-
imation and the extrapolated boundary condition [RCN01]. Extending this method to
account for tissue heterogeneity, we compute a perturbed model in which the diffusion
and absorption coefficients at each pointr inside the animal volume have base values
D
0
andµ
a0
andperturbationsδD(r)andδµ
a
(r)respectively. UsingthefirstorderBorn
approximation, the Green’s function, G(r
d
,r
s
), from a source point r
s
to a detector
pointr
d
isgivenby:
G(r
d
,r
s
) =G
h
(r
d
,r
s
)+
∫∫∫
V
G
h
(r
d
,r)[∇·(δD(r)∇G
h
(r,r
s
))−δµ
a
(r)G
h
(r,r
s
)]d
3
r,
(3.1)
where G
h
(r,r
s
) denotes the Green’s functions at a point r due to the source point
r
s
inside homogeneous animal volume, V. Whenr is an internal point, we compute
15
G
h
(r,r
s
) using the infinite, homogeneous medium assumption. For pointsr
d
on the
surface,wecomputeG
h
(r
d
,r
s
)usingasemi-infiniteslabapproximationalongwiththe
methodofimages. Theimagesourcepoint,r
si
,istheimageofr
s
acrossanextrapolated
boundaryabovetheanimalsurface. Thus,weconsiderG
h
(r
d
,r
s
)tobeoftheform:
G
h
(r
d
,r
s
) =
1
4πD
0
(
e
−
eff
|r
d
−rs|
|r
d
−r
s
|
−
e
−
eff
|r
d
−r
si
|
|r
d
−r
si
|
)
, (3.2)
whereµ
eff
isafunctionoftheabsorptionandreducedscatteringcoefficients. Basedon
theseassumptions,wehavesimplifiedequation3.1tothefollowingform:
G(r
d
,r
s
) =G
h
(r
d
,r
s
)+
∫∫∫
V
G
h
(r
d
,r)G
h
(r,r
s
)[µ
2
eff
δD(r)−δµ
a
(r)]d
3
r. (3.3)
Thisequationcanbeusedtocomputeeachentryoftheexcitationandemissionforward
modelmatrices,A
ex
andA
em
,describedinsection2.2.
3.2.2 Implementation
Instead of assigning optical properties to a large number of tetrahedral elements (as is
typically done for an FEM solver), we assign these properties to a smaller set of uni-
formly spaced grid points inside the mouse atlas. The optical properties were assigned
organ-wise based on published results [CPW90]. The set of internal source locations
for which the forward model was computed was the same as the set of grid points over
which the integral in equation 3.3 was evaluated. Using the formulation described in
section2.2,theA
ex
andA
em
matriceswerecomputedforFMTstudies.
16
3.3 ForwardModelComputation
Simulations were performed using the Digimouse atlas (http://neuroimage.
usc.edu/Digimouse.html) [DSCL07, SCS
+
02], a labeled atlas based on co-
registered CT and cryosection images of a 28g nude normal male mouse. First, we
examinetheaccuracyofourmodelbyperformingacorrelationanalysisbetweensurface
photon densities computed by our approach and an FEM-based forward model. Then,
we present efficiency comparisons of our approach with the more accurate FEM-based
approach.
3.3.1 Correlationanalysis
Our reference for the correlation study is an FEM-based forward model that uses vol-
ume tessellation of the Digimouse atlas. The atlas is tessellated into 306,773 tetrahe-
drons. Opticalpropertiesareassignedorgan-wisebasedonpublishedvalues[CPW90].
We generate monochromatic emission forward models (A
em
) for 9192 internal source
pointsarrangedasauniform,rectangulargridwitha1.2mmspacing,562surfacedetec-
tor nodes, and at a wavelength of 660 nm using our approach and FEM. Figure 3.1(b)
shows the correlation coefficients between the FEM model and our model. A similar
analysis is done to compare an analytical model based on homogeneous tissue proper-
ties and the extrapolated boundary condition with FEM. The result is shown in figure
3.1(a). This study reveals that our perturbative correction yields a better average corre-
lationcoefficientthanthehomogeneousmodelasshownintable3.1.
3.3.2 Efficiencycomparison
Table 3.1 provides the computation times on an AMD OPTERON 250 dual CPU com-
puter for generating monochromatic forward models based on the homogeneous tissue
17
Figure 3.1: Comparison of correlation coefficients. The correlation coefficients of the
photondensityateachdetectorduetoallthesourcepointsaredisplayedfor(a)ahomo-
geneoustissuemodelandFEM(b)ourperturbativemodelandFEM.
Table3.1: Comparisonofaccuracyandefficiencyofdifferentmonochromaticemission
forwardmodels
Criterion Homogeneous
tissuemodel
Perturbative
model
FEM-based
model
Average
correlation
coefficient
0.9293 0.9547 1
Run time
(seconds)
13 42 920
assumption based approach, our approach, and the FEM-based approach. MATLAB
R ⃝
(TheMathworksInc.,Natick,MA,USA)wasusedforallcomputations. Sinceincreas-
ingthenumberofdetectornodesimprovesresolution,theincreaseincomputationcosts
with the number of detector nodes was studied. Our studies indicate that the computa-
tioncostforourapproachincreaseslinearlyasthenumberofdetectornodesisincreased
whereas,forFEM,thecostvariationisnegligible. However,theruntimefortheformer
18
Figure 3.2: Comparison of the increase in computational costs of our perturbative for-
wardmodelandtheFEM-basedforwardmodelwithincreaseinthenumberofdetector
nodes.
continues to be significantly lower than that for FEM even when over 16,000 detector
nodesareused. Theseresultsareshowninfigure3.2.
3.4 PhantomExperiment
Finally, we test our forward modeling scheme by using it to reconstruct a fluorescent
source embedded inside an inhomogeneous phantom. As shown in figure 3.3(a), the
phantom is a 32 mm× 32 mm× 29 mm block with a 7.5 mm thick front layer that is
approximately 4 times more scattering than the remaining bulk and is designed to have
optical properties close to that of mouse tissue [CPT
+
97]. The back part of the block
is prepared from a solution of 1% intralipid and 1.5% agar in deionized water, while
the front layer is derived from the same contents, but with 4 times as much intralipid.
Theintralipidactsastheprimaryscatterer,withthereducedscatteringcoefficientvary-
ingwithfrequency[PP06], whereasthewatermoleculesthatremainaftercuringofthe
phantom are responsible for absorption of light. A 14 mm long capillary tube with an
19
Figure 3.3: (a) Schematic of the rectangular phantom showing the front layer that is 4
timesmorescatteringthantheremainingbulk. (b)The15pointsourcelocationsonthe
frontsurfaceusedforexcitation. (c)Apyramidalmirrorsetupthatprovidestop,bottom,
andsideviews.
insidediameterof1mmfilledwithfluorescentDiDdyewithameasuredemissionspec-
tral peak at 730 nm is embedded inside the phantom. We use an excitation wavelength
of 650 nm and 15 different excitation points on the front surface of the block as shown
infigure3.3(b). Thephantomisplacedonatransparentstage. Apyramidalmirrorsetup
that provides top, bottom, and side views of the block, as shown in figure 3.3(c) and a
CCDcameraplacedinfrontofitareusedtoacquirefluorescencedatafrom5surfacesof
thephantom. Multispectraldatawascollectedfor20nmbinswithwavelengthsranging
from710nmto730nm.
MultispectraldataatthreewavelengthsalongwithourforwardmodelandanFEM-
based model were used for reconstruction. Reconstruction was based on regularized
least squares using the gradient projection algorithm. Horizontal sections of the phan-
tom illustrating the original source location and the reconstructed source distributions
are shown in figure 3.4. The reconstruction results from our model and the FEM-based
model look similar, although the latter exhibits better resolution. The horizontal full
width at half maximum (FWHM) of the source distribution reconstructed using our
model was 5.5 mm. The results indicate that, with increasing depth, the increase in
20
Figure3.4: Horizontalsectionsofthe32mm×32mm×29mmphantomshowing(a)
theoriginalsourcepositionandthereconstructedsourcedistributionsobtainedusing(b)
ourperturbativeforwardmodeland(c)anFEM-basedforwardmodel.
localization error along the horizontal plane is more pronounced in our model than in
themoreaccuratetheFEMmodel.
3.5 Summary
Wehavedevelopedafast,perturbativeforwardmodelingschemeformultispectralFMT
based on the Born approximation in conjunction with an extrapolated boundary condi-
tion. Our model takes into account both animal geometry and heterogeneity of tissue
optical properties. We have presented simulation studies comparing the accuracy and
21
efficiency of our model with a homogeneous model and an FEM-based model. Corre-
lation analyses performed on the Digimouse atlas using our model and a homogeneous
model, both being with respect to an FEM-based model, indicate that our model has
better accuracy than the homogeneous model. Our studies show that the computational
savings with our model relative to FEM are significant. Finally, we have successfully
testedourmodelinamultispectralFMTphantomexperimentusinganinhomogeneous
phantom.
22
Chapter4
IlluminationPatternOptimization
4.1 Introduction
The high degree of absorption and scattering of photons propagating through tissue
makestheFMTproblemill-conditionedandposesasignificantchallengefor3Dimage
reconstruction. We seek to alleviate this problem and improve source localization by
exploitingthedegreeoffreedom offeredby externalillumination inFMT.In thischap-
ter, we present an approach to design an efficient set of spatial illumination patterns
that improve the conditioning of the forward model matrix [DAJL09, DAJL10]. FMT
setups typically acquire data-sets corresponding to different surface illumination pat-
terns. These patterns generate different excitation fields over the volume, which tune
the system matrix. FMT setups available today employ illumination schemes chiefly
guided by the availability and simplicity of the light source. Several of these use
laser sources with focusing or diffuser lenses to generate point or distributed patch
patterns [GRWN03, ZVM
+
06, LMD
+
09b]. Other approaches include raster scanning
[PBAC05, JBSM06] and structured light or spatially modulated illumination patterns
[LMS09]. Most of these approaches are ad hoc. Although various performance met-
rics could be used to theoretically compare these standard approaches, it is impossible
to make an exhaustive set of comparisons, since, for a given number of illumination
patterns being used for an experiment, infinitely many designs exist. Therefore, the
question we address in here is: given a fixed number of illumination patterns, how do
wedesignthesepatternssoastomaximizetheinformationintheacquireddata?
23
With the availability of Texas Instruments Digital Light Processor (DLP
R ⃝
) chips
[Hor96] which work in the near-infrared range and give us precise control over the
spatialintensitydistribution[DDS03],itisfeasibletogenerateanysetofspatialillumi-
nation patterns with grayscale intensity variation. Our focus here is to compute the set
ofilluminationpatternsforCWFMTthatmaximizetheinformationcontentinthedata
byimprovingtheconditioningoftheFisherinformationmatrix. Weformulateourprob-
lem as a constrained optimization problem that minimizes a cost function derived from
the Fisher information matrix and computes the parameterized set of optimal spatial
patterns.
4.2 FormulationoftheOptimizationProblem
According to the Cram´ er Rao inequality, the inverse of the Fisher information sets a
lower bound on the variance of an unbiased estimator [Kay93]. This establishes a reci-
procity between the variance and the information content of an estimator. For a scalar
estimator, the Fisher information is a scalar and estimator variance will be small when
Fisher information is large. However, when we are dealing with a vector estimator, the
Fisherinformationisamatrixofelements,andthereisnosinglescalaroptimalitycrite-
rion. In the field of optimal design of experiments, a variety of different functionals of
theeigenvaluesoftheFisherinformationmatrix(FIM)aretraditionallyusedasoptimal-
ity criteria. For example, A-optimality minimizes the trace of the inverse of the FIM,
D-optimality maximizes the determinant of the FIM, and E-optimality maximizes the
smallest eigenvalue of the FIM [Puk93]. The approach that we present here is inspired
byE-optimality. InFMT,simplybyfeedinginmoreilluminationpower,wesimultane-
ously increase all eigenvalues of the system, not just the smallest one. Maximizing the
smallesteigenvalueforafixedpowerlevelamountstominimizingtheratioofthelargest
24
to the smallest eigenvalues here. This amounts to improving the condition number of
theFIMand,hence,ofthesystemmatrix.
4.2.1 Fisherinformationmatrix
For a vector unknownq and observed datab with probability densityp(b), the FIM is
givenby[Kay93]:
F =−E
[
∂
2
∂q
2
lnp(b)
]
. (4.1)
We assume an additive white Gaussian noise model with noise vectorn∼N(0,σ
2
n
I)
anddatavectorb =Aq+n. Withoutlossofgenerality,weignoretheconstantmultiple,
σ
2
n
, in the covariance. Using (2.5) and (2.8), the Fisher information matrix for this
systemcanbeexpressedas:
F =A
′
A =
p
∑
k=1
D
ex
k
(A
em
)
′
A
em
D
ex
k
. (4.2)
Using(2.6),(2.7),and(4.2),theFIMcanbeexpressedasafunctionofthesetofillumi-
nation patterns,W = [w
1
w
2
...w
p
], whereW ∈R
n
d
×p
. We denote theith row vector
ofA
ex
by (a
ex
i
)
′
and theith column vector ofA
em
bya
em
i
. Then theijth term of the
FIMisgivenby:
F
ij
=
p
∑
k=1
(w
′
k
a
ex
i
)
[
(a
em
i
)
′
a
em
j
]
(w
′
k
a
ex
j
)
=
p
∑
k=1
F
em
ij
·
[
(a
ex
i
)
′
w
k
w
′
k
a
ex
j
]
= F
em
ij
·
[
(a
ex
i
)
′
p
∑
k=1
(w
k
w
′
k
)a
ex
j
]
= F
em
ij
·
[
(a
ex
i
)
′
WW
′
a
ex
j
]
, (4.3)
25
where we have used the notationsF
em
ij
= ((A
em
)
′
A
em
)
ij
= (a
em
i
)
′
a
em
j
. This implies
thattheFIMcanbeconvenientlyrepresentedasaHadamard(entrywise)matrixproduct
of two matrices, one dependent on the emission model and the other on the excitation
modelandtheilluminationpatterns:
F = [(A
em
)
′
A
em
]◦[A
ex
WW
′
(A
ex
)
′
]. (4.4)
Introducing the notationF
em
= [(A
em
)
′
A
em
] andG
ex
= [A
ex
WW
′
(A
ex
)
′
], we can
rewrite(4.4)as:
F =F
em
◦G
ex
. (4.5)
4.2.2 Optimizationproblemformulation
For the system matrix, A = UV
′
, to have a perfect condition number of 1, the
diagonal matrix of singular values should be σI where σ is a singular value. Note
that the singular value σ does not affect the condition number. Then the FIM for the
idealsystemis:
F
ideal
=A
′
A =V
2
V
′
=σ
2
VV
′
=σ
2
I. (4.6)
26
Thus,forperfectconditioning,F shouldapproachF
ideal
. Wefindasetofillumination
patterns,W,whichminimizestherelativedifferenceoftheresultantFIM,F(W),and
theidealone,F
ideal
:
Θ(W,σ) =
∥F(W)−F
ideal
∥
2
F
∥F
ideal
∥
2
F
=
∥F(W)−σ
2
I∥
2
F
∥σ
2
I∥
2
F
=
1
n
s

1
σ
2
F(W)−I

2
F
=
1
n
s

F
(
W
σ
)
−I

2
F
, (4.7)
whereF(W) is defined in (4.4) and∥·∥
F
represents the Frobenius norm. The mini-
mizer,(
ˆ
W,ˆ σ),ofthecostfunctionΘ(W,σ)withthenonnegativityconstraint,W ≥0,
isnotunique,sinceforanyγ, Θ(γ
ˆ
W,γˆ σ) = Θ(
ˆ
W,ˆ σ). Wethereforefirstsolveasim-
pleroptimizationproblem:
W
opt
= arg min
W≥0
||F(W)−I||
2
F
. (4.8)
Then we determine the scaling factor by imposing a power constraint,∥W∥
F
= β, on
theilluminationpatterns. Thefinalsolutionisgivenby:
(
ˆ
W,ˆ σ) = arg min
;W≥0; ∥W∥
F
=
Θ(W,σ)
=
(
β
∥W
opt
∥
F
W
opt
,
β
∥W
opt
∥
F
)
. (4.9)
In practice, the illumination is scaled to exploit the maximum available source power
and thus ensure high signal-to-noise ratio. Since the optimal solution
ˆ
W to minimize
27
thecostfunctionΘissimplyascaledversionofthesolutionW
opt
,wenowfocusonthe
simplerproblemgivenin(4.8).
4.2.3 Dimensionalityreduction
The formulated problem in (4.8) is fourth-order and non-convex with respect to the
patterns, W. W is of size n
d
×p, which is typically large, and hence solving this
problem is computationally intensive. Additionally, the solutions are not necessarily
spatially smooth and therefore not convenient for implementation. To allow smooth
solutionsandtoenablespeediercomputation,weparameterizetheilluminationpatterns
using a small number of spatially smooth basis functions. The transformation from the
basisfunctiondomaintothespatialpatterndomainisgivenby:
w
k
=Ly
k
,
W =LY. (4.10)
Herey
k
∈ R
m
is a vector of the linear coefficients of the basis functions for the kth
illumination pattern,L ∈ R
n
d
×m
is a matrix whose columns are the basis vectors, and
Y ∈ R
m×p
and Y = [y
1
y
2
... y
p
], m being the number of basis vectors used to
encodeeachspatialpattern. Then,theFIMcanberepresentedintermsofY asfollows:
F = ((A
em
)
′
A
em
)◦(A
ex
LYY
′
L
′
(A
ex
)
′
). (4.11)
28
4.2.4 Modifiedoptimizationproblem
With the inclusion of basis functions for dimensionality reduction, the modified opti-
mizationproblemisasfollows:
Y
opt
= arg min
LY≥0
||F −I||
2
F
,
W
opt
= LY
opt
. (4.12)
Agradient-basedapproachcanbeusedtosolvetheoptimizationproblemin(4.12). Our
costfunctionasafunctionoftheargumenty
k
,wherek = 1 :p,isasfollows:
Φ =||F −I||
2
F
=
ns
∑
i=1
ns
∑
j=1
(F
ij
−δ
ij
)
2
=
ns
∑
i=1
ns
∑
j=1
Ψ
2
ij
, (4.13)
whereδ
ij
istheKroneckerdeltafunctionandΨ
ij
isgivenby:
Ψ
ij
=F
em
ij
(a
ex
i
)
′
L
(
p
∑
k=1
y
k
y
′
k
)
L
′
a
ex
j
−δ
ij
. (4.14)
Weobservethatthesummationin(4.13)isoverthesetofdiscretepointsources. Thus,
if the biodistribution of a fluorophore is known to be within a specific local region of
interest,anoptimalsetspecifictothatcasecouldbeobtainedbysummingoverthepoint
sources within that region. The gradient,∇
k
Φ, of the cost function with respect to the
kthargumentvector,y
k
,isgivenby:
∇
k
Φ =
ns
∑
i=1
ns
∑
j=1
2Ψ
ij
∇
k
(F
ij
−δ
ij
)
=
ns
∑
i=1
ns
∑
j=1
2Ψ
ij
F
em
ij
(∇
k
G
ex
ij
). (4.15)
29
∇
k
G
ex
ij
in(4.15)canbecomputedasfollows:
∇
k
G
ex
ij
= ∇
k
(
(a
ex
i
)
′
L
(
p
∑
k=1
y
k
y
′
k
)
L
′
a
ex
j
)
= ∇
k
(
(a
ex
i
)
′
Ly
k
y
′
k
L
′
a
ex
j
)
= L
′
(a
ex
j
(a
ex
i
)
′
+a
ex
i
(a
ex
j
)
′
)Ly
k
. (4.16)
Substituting(4.16)in(4.15),wehave:
∇
k
Φ =
ns
∑
i=1
ns
∑
j=1
2Ψ
ij
F
em
ij
L
′
(a
ex
j
(a
ex
i
)
′
+a
ex
i
(a
ex
j
)
′
)Ly
k
. (4.17)
4.2.5 Extensiontomultispectraldetection
Ourformulation can beextendedfor applicationto multispectral detection. Weassume
thatforeachilluminationpattern,fluorescencedataiscollectedoverswavelengthbins.
The fluorophore is assumed to have an emission spectral coefficientη
l
at thelth wave-
length bin. Then, the block of the system matrix corresponding to thekth illumination
patternisgivenbyequation2.8. TheFisherinformationmatrixforthiscaseis:
F =
s
∑
l=1
p
∑
k=1
(
η
2
l
D
ex
k
(A
em
l
)
′
A
em
l
D
ex
k
)
=
(
s
∑
l=1
η
2
l
(A
em
l
)
′
A
em
l
)
◦(A
ex
WW
′
(A
ex
)
′
)
= F
em−ms
◦G
ex
, (4.18)
whereF
em−ms
=
(∑
s
l=1
η
2
l
(A
em
l
)
′
A
em
l
)
. TheFIMretainstheHadamardproductform
with the illumination-dependent matrix, G
ex
, unchanged. This implies that once the
30
multispectral emission model dependentF
em−ms
matrix is precomputed, the optimiza-
tion problem can be solved for the multispectral emission case without any additional
computationalcost.
4.3 Methods
The problem formulated in section 4.2 is first solved for an optically homogeneous,
symmetric solid shape – a cylinder – to gather geometric intuition. Next, we solve the
same problem for an optically inhomogeneous realistically shaped mouse atlas. We
perform source localization studies using the optimal patterns on the mouse atlas and
evaluatetheperformanceofthepatternsbasedonappropriatemetrics.
4.3.1 Computationoftheforwardmodel
The finite element method (FEM) was used to solve the coupled PDEs (2.1) and (2.2)
subjecttoboundaryconditions(2.3)and(2.4)respectivelytoprecomputetheexcitation
andemissionforwardmodelmatricesA
ex
andA
em
.
Tessellatedcylinder
For simulations on a cylinder, we use a tessellated cylinder of height 40 mm and diam-
eter 20 mm with 2373 tetrahedrons and 588 tessellation nodes. A set of 361 detector
nodesonthesurfaceareusedforilluminationanddetectionand1564pointsourceloca-
tionsonauniform,volumetricgridwitha1.5mmspacingareusedforsourcelocaliza-
tion. We use monochromatic excitation atλ
ex
= 650 nm and monochromatic detection
atλ
em
= 730 nm. Opticalpropertiesareassumedtobehomogeneouslydistributedand
areassumedtoresemblethoseformuscletissue.
31
Digimouseatlas
For more realistic simulations, we use the Digimouse atlas (http://neuroimage.
usc.edu/Digimouse.html) [DSCL07, SCS
+
02], a labeled atlas, based on co-
registeredCTandcryosectionimagesofa28gnudenormalmalemouse. Tissueoptical
properties are assumed to be spatially inhomogeneous and assigned organ-wise, for 17
different organs, according to published values [ARC05]. The tessellated atlas volume
consists of 306,773 tetrahedrons and 58,244 nodes. For generating the patterns, 810
surface nodes and 1129 internal grid points with a uniform 2.4 mm grid spacing are
used. We assume that all surface nodes except those lying on the limbs and the tip of
the snout are used for illumination. The AlexaFluor 700 dye, with an emission peak at
719 nm is used as the fluorescent probe. Accordingly, we use an excitation wavelength
ofλ
ex
= 650nmandperformmultispectraldetectionat710nm,730nm,and750nm.
4.3.2 Optimizationprocedure
The fourth-order cost function in (4.13) is minimized using fmincon, an inbuilt func-
tion for constrained optimization in MATLAB
R ⃝
(The Mathworks Inc., Natick, MA,
USA). This function chooses the active set method for optimization [GMW82]. The
non-convex nature of the cost function implies that several local minima could exist,
and our gradient-based optimization approach could get stuck in one of these. Hence,
weuseamulti-startapproachwithrandominitializations. Thisisaheuristicbutstraight-
forward approach to the global minimization problem in which the local minimization
algorithmisrunfrommanydifferentstartingpoints,andthenthebestsolutionispicked
[LM09].
Fordimensionalityreduction,wedefinebasisfunctionsonthecylinderandatlassur-
faces by treating them as smooth 2D manifolds. We use as our basis set eigenfunctions
32
Figure4.1: (a)Frontandbackviewsofthecurvedsurfaceofatessellatedcylinder(flat-
tened for display purposes with the cylinder axis aligned horizontally) and (b) top and
bottom views of the Digimouse atlas showing eigenfunctions corresponding to the 20
smallesteigenvaluesoftheLaplace-Beltramioperator. Thebasisfunctionsareorthogo-
naland,therefore,havebothpositiveandnegativevaluesatvariouspoints.
oftheLaplace-Beltramioperator[QBM06,Jos08,GM07]. TheLaplace-Beltramioper-
atorisageneralizationoftheLaplaceoperatoronRiemannianandpseudo-Riemannian
manifolds. Thesebasisfunctionsallowustogenerateaspatiallysmoothsolution. Their
orthonormality allows efficient representation of the solution. Since the eigenvalues of
33
a spatial differential operator are measures of the spatial rate of change of the corre-
sponding eigenfunctions, the smaller the eigenvalue the smoother is the corresponding
eigenfunction. Therefore, eigenfunctions corresponding to them = 20 smallest eigen-
values of the operator were used as basis functions. Figure 4.1 shows front and back
views of the tessellated cylinder and top and bottom views of the Digimouse atlas with
thebasisfunctionsdisplayedonthem.
4.3.3 Inversesolutionandperformancemetrics
The inverse problem is solved using the regularized pseudoinverse in (2.13). The illu-
minationpatternsontheDigimousesurfaceareinterpolatedontoadensersurfacetrian-
gulationwith3234nodes. Forsourcelocalizationstudies, auniformgridof9192point
sources with a 1.2 mm spacing is used. Our simulation setup uses surfacedata from all
nodes on the atlas surface except those lying on the limbs and the tip of the snout. We
use two approaches for comparing the performance of different illumination patterns.
Thefirstofthese,thesingularvaluedistribution,throwslightontheconditioningofthe
forwardmatrix. Thesecondapproachistolookataverageresolution-variancecurvesto
examinethebehavioroftheinversesolution.
Conditionnumber
The condition number of a matrix is given by the ratio of its largest to its smallest
singularvalue. Althoughthisisagoodfigure-of-meritfortherobustnessofthesystem,
the quality of the actual reconstruction results depends on the singular vectors as well.
Sothismetric,althoughinformative,isnotsufficient.
34
Resolution-varianceanalysis
Imaging systems exhibit an inherent trade-off between the spatial resolution and the
noisevarianceofthereconstructedimages. Fortheregularizedpseudoinverseapproach,
themeanvalueoftheestimatorforaregularizationparameterαisgivenintermsofthe
singularvaluedecomposition,A =UV
′
,ofthesystemmatrixby:
E[^ q] =V
2
(αI +
2
)
−1
V
′
q. (4.19)
For the jth unit point source, q = e
j
, a unit vector with 1 as the jth element. Then
(4.19)givesusthepointspreadfunction(PSF)forthepointsource,q =e
j
:
PSF
j
=V
2
(αI +
2
)
−1
V
′
e
j
. (4.20)
The spatial resolution is the full width at half maximum (FWHM) of the point spread
function. ForanadditivewhiteGaussiannoisemodelwithcovarianceσ
2
n
I,theestimator
varianceaveragedovernpointsourcelocationsis:
Var(^ q) =
pσ
2
n
n
n
∑
j=1
(
e
′
j
V
2
(αI +
2
)
−2
V
′
e
j
)
. (4.21)
Using these equations, the resolution and the variance can be computed for different
pointsourcelocationsinsidethevolumefordifferentvaluesoftheregularizationparam-
eter,α. For comparison of illumination setups with different numbers of patterns (p =
1, 3, 6, and 12), the effect of different scan times is eliminated by assuming that, when
p patterns are used, the acquisition time for data for each pattern is less by a factor ofp
relativetothep = 1case. Thisassumptionreflectsasascalingfactorpintheexpression
for the noise variance in (4.21). The curves of resolution vs. variance are obtained by
averagingeachofthesetwoquantitiesoverdifferentsourcelocationsandplottingthem
35
against each other for different values of the regularization parameter. Comparison of
these curves for different illumination patterns throws light on the relative magnitudes
ofmean-squarederrorforestimatorscorrespondingtodifferentilluminationpatterns.
Itmustbenotedthat(4.20)and(4.21)dependononlyV andandnotonU. Thus,
for the purpose of this analysis it is sufficient to compute the SVD ofA
′
A =V
2
V
′
.
Since, for a large number of illumination patterns and/or a large number of detection
wavelengths,thenumberofrowsinA,andhencethesizeofthematrixU,ishuge,the
computationalsavingsincurredusingthisapproacharesubstantial.
Figure4.2: Asetof3optimalpatternsonatessellatedcylinder. Theindividualpatterns
are displayed in (a), (b), and (c). The sum of the 3 patterns displayed in (d) shows
rough overall uniformity of illumination. For ease of display, the curved surface of the
cylinder is laid out flat so that the vertical axis represents the height of the cylinder and
thehorizontalaxisrepresentsthetransverseangle.
36
Figure 4.3: Projections of the top and bottom surfaces of the Digimouse atlas showing
optimalpatternsfor(a)p = 1,(b)p = 3,(c)p = 6,and(d)p = 12.
37
4.4 Results
4.4.1 Patternsonthecylinder
A set of three optimal patterns computed for the optically homogeneous, tessellated
cylinder is shown in figure 4.2. It is assumed that only the curved surface is used for
illumination. The patterns are displayed on the curved surface by laying it out flat. We
observethatthepattern(b)canbeshiftedby120
◦
(one-thirdofafullcircle)rightorleft
toobtainpatterns(a)and(c)respectively. Thesumofthethreepatternsin(d)looksquite
uniformoverallandadditionallyindicateroughbilateralsymmetryaboutthehorizontal
planecuttingthecylinderlengthwisemidway. Thus,theoptimalsetofpatternsexhibits
radialandbilateralsymmetries,characteristicofthecylindricalshape.
4.4.2 PatternsontheDigimouseatlas
OptimalpatternsontheDigimouseatlasforp = 1,3,6,and12areshowninfigure4.3.
Foreachset,theintensitiesarenormalizedwithrespecttothemaximumintensityvalue
for that set for display purposes. Unlike the cylinder, the mouse atlas is an optically
inhomogeneous, irregular shape. Accordingly, the patterns do not exhibit any clear
symmetry.
4.4.3 Pointsourcelocalization
We looked at the PSF of a sample deep source in the trunk of the mouse at a depth of
10.4 mm below the top surface and 7.6 mm from the bottom surface. The PSFs were
computedusing(4.20)fordifferentsetsofilluminationpatterns. Figure4.4showsslices
of the Digimouse atlas displaying a sample deep source and the point spread functions
for the different cases. To generate these PSFs, we picked values of the regularization
38
Patterntype FWHM(mm)
Optimal,p = 12 1.20
Optimal,p = 6 1.20
Optimal,p = 3 1.63
Uniform,p = 3 3.02
Optimal,p = 1 5.91
Table4.1: FullwidthathalfmaximumofpointspreadfunctionsinFigure4.4
parameter that equalize the estimator variances as computed using (4.21) for the dif-
ferent cases. We observe improvement in the measured FWHM from the casesp = 1
through12. Forreference,wehaveusedasetofp = 3patternseachofwhichuniformly
illuminates a different longitudinal section of the mouse. We observe that the recon-
structed point source for the optimal case for p = 3 appears to have better resolution
than thep = 3 uniform case for the same estimator variance. The FWHMs of the PSFs
forthedifferentcasesareprovidedintable4.1.
Figure 4.4: Point spread functions for (a) optimal patterns for p = 12, (b) optimal
patterns forp = 6, (c) optimal patterns forp = 3, (d) uniform patterns forp = 3, and
(e)optimalpatternforp = 1.
39
Patterntype Conditionnumber
Optimal,p = 12 2.5575×10
9
Optimal,p = 6 1.3466×10
10
Optimal,p = 3 7.1619×10
11
Uniform,p = 3 2.8376×10
12
Optimal,p = 1 3.3014×10
13
Table4.2: Conditionnumberofthesystemmatrixfordifferentilluminationschemes
4.4.4 Performanceevaluation
Figure4.5showsthesingularvaluedistributionsforthedifferentilluminationschemes.
Thesingularvaluesareplottedindescendingorderandnormalizedbythelargestsingu-
lar value. The condition numbers of the system matrices are provided in table 4.2. The
condition numbers clearly improve as p is increased from 1 through 12. The optimal
illumination pattern forp = 3 is observed to generate a better conditioned system than
thep = 3uniformilluminationscheme.
Figure4.5: Singularvaluedistributionsfordifferentilluminationschemes.
40
The plots of average resolution vs. variance for the different illumination schemes
areshowninfigure4.6. Theseareobtainedbyaveragingresolutionandvarianceforn =
250 randomly picked sources for different values of the regularization parameter. The
best possible FWHM is 1.2 mm, which is the grid spacing. The idea underlying these
curves is the trade-off between resolution and variance, which implies that, a decrease
in estimator variance is accompanied by an increase in FWHM and vice versa. The
intensities of the different sets were normalized to ensure that all sets of patterns have
thesameaverageintensity. Asexpected,aspgoesfrom1to12,theresolution-variance
properties improve. Also, the set of optimal patterns for p = 3 offers better average
resolutionthanthesetofuniformpatternsforp = 3forthesamevariance.
Figure4.6: Averageresolutionvs. variancecurvesfordifferentilluminationschemes.
4.5 Summary
We have developed an optimization framework for generating optimal spatial illumina-
tion patterns for CW FMT based on an approach that seeks to improve the condition
41
number of the Fisher information matrix. We formulated our problem as a constrained
optimization problem which, for a given number of patterns, can be solved to com-
putethesetofoptimalpatternswhichmaximizetheinformationcontentintheacquired
data. Our formulation assumes an i.i.d. Gaussian noise model. The fourth order, non-
convexcostfunctionwasminimizedinthepresenceofanon-negativityconstraintusing
a gradient-based approach. A set of random initializations was used to ensure global
convergence. The dimensionality of the optimization problem was reduced using a set
of geometrical basis functions defined on the 2D manifold representing the surface of
theobjectbeingimaged. EigenfunctionsoftheLaplace-Beltramioperatorwereusedas
basis functions owing to their spatial smoothness and orthogonality. To ensure that this
numerical approach is geometrically meaningful, we looked at optimal patterns gener-
atedforanopticallyhomogeneous,tessellatedcylinderandobservedradialandbilateral
symmetries that are intrinsic to the cylindrical shape. Finally, we used our method to
compute optimal patterns for the optically inhomogeneous, realistically shaped Digi-
mouseatlas. Theobtainedpatternslooksmoothandphysicallyrealizable. Forthesame
number of patterns,p = 3, the optimal set was shown to perform better than a uniform
illumination scheme on the basis of condition number and average resolution vs. vari-
ancecurves. SincetypicalFMTexperimentalsetupsusealargernumberofillumination
patterns,wehavegeneratedoptimalpatternsforuptop = 12andpresentedtheirsingu-
lar value spectra and average resolution-variance characteristics. For fair comparison,
thevariancewascomputedbasedontheassumptionofequalityoftotaldataacquisition
timesfordifferentvaluesofp.
42
4.6 Discussion
In our simulations we have assumed that the entire surface of the mouse except for the
limbs and the tip of the snout is being illuminated. However, in an experimental set-
ting, it may not be feasible to illuminate the entire mouse surface at one go. For such
applications, our formulation can be easily modified to introduce spatial constraints on
the illumination patterns. Another useful modification of our formulation would be its
extension to include multispectral excitation. Owing to the larger variability of optical
properties at the excitation wavelengths of commonly used fluorescent dyes and pro-
teins, for the same number of wavelength bins, multispectral excitation typically leads
toabetterconditionedsystemmatrixthanmultispectraldetection[CAL
+
09]. Wethere-
fore believe that optimal illumination coupled with multispectral excitation will allow
us to pack more information into the acquired data, for a given number of illumination
patternsandagivennumberofwavelengths.
The framework for generating optimal patterns assumes prior knowledge of mouse
surfacetopographyandtissueopticalpropertiesandrequirestheforwardproblemtobe
solved prior to the optimization procedure. It might not be feasible to repeat the opti-
mization procedure for each animal in between the surface profiling and fluorescence
data acquisition steps of an experiment. An alternative approach would be to use the
atlasasasurrogateandtowarpthesurfaceoftheatlastomatchthetargetmousebeing
imaged [ACL
+
09], and in doing so, we can warp the optimal patterns onto the surface
ofthetarget.
43
Chapter5
ImagingSystems
5.1 Introduction
A typical CW FMT system consists of a surface profiling setup, an excitation source, a
camera for detecting fluorescence light, optics for full-surface viewing, and excitation
andemissionfilters. Someexamplesofsurfaceprofilingschemeswereprovidedinsec-
tion2.1. Lasersourceswerebrieflydiscussedinsection4.1. Frequency-selectivefilters
for excitation and fluorescence light are required for both single-wavelength and multi-
spectraldataacquisition. Thechoiceofthesefiltersdependsonthedegreeofselectivity
required and the excitation and emission spectra of the dye used. CCD cameras, which
integrate incoming photons, are commonly used for contact-free acquisition of fluores-
cencedatainCWFMT.Toenablefull-surfacevisualizationwithoutthehassleofusing
multiplecameras,additionalopticsarerequiredtobringtheentiresurfaceoftheanimal
into the field-of-view of the CCD camera. In this chapter, we will look at two different
imagingsetupsandexploresomecalibrationissues.
5.2 Pyramidal Mirror Setup with Patterned Illumina-
tion
This system, developed in collaboration with the Cherry Lab at University of Califor-
nia,Davis,andCambridgeResearchandInstrumentationInc.,Woburn,MA,canproject
44
digitally controlled spatial illumination patterns [GDM
+
10]. These patterns are gener-
ated by a Sharp XR-10XR-10S projector that uses a Texas Instruments DLP
R ⃝
chip.
The projector light is passed through an excitation filter (Semrock, Rochester NY and
Chroma Technology Corp., Rockingham VT) contained in a computer-controlled filter
wheel and directed toward the sample stage using a steerable mirror. The sample sits
on a glass stage and can be illuminated from the top, bottom, left, and right sides by
orienting the 45
◦
steerable mirror under computer control. Various spatial patterns of
illuminationcanbeprojectedontotheselectedsideofthesamplebysendingappropriate
bitmap images to the projector from the control computer. The surface fluorescence is
simultaneouslyimaged from top, bottom, left, and rightusing a pyramidalmirrorsetup
as shown in figure 5.1. The intermediate image is passed through an emission filter
(Omega Optical, Brattleboro VT) on another controlled wheel. Additional optics are
used to increase the number of useable pixels on the CCD camera by 290% compared
to direct capture. The final image is captured by a 4 megapixel cooled CCD camera
(PhotometricsCoolSNAPK4,TusconAZ).
Figure5.1: (a)Acuboidalphantomplacedonthetransparentstageinsidethepyramidal
mirror setup. (b) CCD camera image of a cubic phantom. The top, bottom, left, and
rightviewsareasindicated.
45
5.3 ConicalMirrorSetup
5.3.1 Systemdescription
Thissystem,developedincollaborationwiththeCherryLabatUniversityofCalifornia,
Davis, uses a silver-coated aluminum conical mirror with a half angle of 45
◦
and an
outerdiameter of 200mm. This mirror setupwasdemonstrated tohavea betterphoton
collection efficiency than the more widely-used plane mirror-based imaging systems
[LMD
+
09b]. An EMCCD camera (C9100-3, Hamamatsu Corp., Bridgewater, NJ) was
used for image capture. Lasers of wavelengths 650 or 785 nm (BWF-OEM-650-200
and BWF-OEM-785-1.0-100-0.2, B&W Tech Inc.) with collimation lenses at the end
of their pigtail fibers were used to generate excitation light. Two motorized scanning
mirrors were used to direct the laser beam to different points on the animal surface. A
photograph of the conical mirror and CCD camera images of a cubic phantom and a
mouseareshowninfigure5.2.
Figure5.2: (a)Acubicphantomplacedonthetransparentstageinsidetheconicalmirror
setup. (b)CCDcameraimageofacubicphantom. (c)CCDcameraimageofamouse.
46
5.3.2 Mappingcalibration
When plane mirrors are used for data acquisition, the task of mapping the fluorescence
datafromthe2DCCDcameraimagesbacktothe3Dobjectspaceisrelativelystraight-
forward. This task is somewhat more involved for the conical mirror setup. We have
developed a calibration scheme for data mapping using the conical mirror setup. The
procedureisdescribedbelow:
Mappingfromthe3Dobjectspacetothe2Dimageplane
Ifweplaceanopaqueobjectofknownshapeataknownpositioninsideaconicalmirror,
we should be able to uniquely determine the location of the image point in the mirror
corresponding to an unoccluded point on the surface of the object. Consider the point
P with cylindrical coordinates (r
p
,ϕ
p
,z
p
) on the surface of an object and its image,Q,
givenby(r
i
,ϕ
i
,z
i
)intheconicalmirror. Fromradialsymmetry,weknow:
ϕ
i
=ϕ
p
. (5.1)
Figure 5.3 shows a radial section of the system at an angleϕ
p
. Distances are measured
with respect to the origin, O, which represents the vertex of the cone. The camera is
placedatpointC ontheaxisoftheconeatadistanceD fromtheorigin.
Weobservefrom∆PSQ,
tanα =
PS
QS
=
z
p
−z
i
r
i
−r
p
. (5.2)
From∆TQO,
tanθ =
TQ
TO
=
r
i
z
i
. (5.3)
47
Figure 5.3: Radial section of the setup at an angleϕ = ϕ
p
showing an object pointP,
(r
p
,ϕ
p
,z
p
)andthecorrespondingimagepointQ,(r
i
,ϕ
i
,z
i
).
The half-angle of the conical mirror is 45
◦
as per our design. In order to account for
any manufacturing imperfection, we assume the half-angle to beθ = 45
◦
+δθ where
δθ→ 0. Then:
1
tanθ
=
1−δt
1+δt
≈ 1−2δt, (5.4)
whereδt = tanδθ. Thisimplies:
z
i
=r
i
(1−2δt). (5.5)
) tanα =
PS
QS
=
z
p
−r
i
(1−2δt)
r
i
−r
p
. (5.6)
Basicgeometrytellsusthat\TQC = 2θ−α. Then,from∆TQC,
cot(2θ−α) =
TQ
TC
=
r
i
D−z
i
. (5.7)
48
Now,cot(2θ−α) = tan(π/2−(2θ−α)) = tan(α−2δθ) = (tanα−2δt)/(1+2δttanα).
SubstitutingthisintheLHSof(5.7),wecanshowthat:
tanα =
r
i
+2δt(D−r
i
)
D−r
i
. (5.8)
Solving(5.6)and(5.8)forr
i
,weget:
r
i
=
Dz
p
+2δtDr
p
D+z
p
−r
p
(1−2δt)
. (5.9)
Thus,using(5.1),(5.5),and(5.9),anygivensurfacepointP,(r
p
,ϕ
p
,z
p
),canbemapped
toanimagepointQ,(r
i
,ϕ
i
,z
i
). Alternatively,ifweknowthexy coordinatesofapoint
inspace(and,hence,r
p
andϕ
p
),wecanestimateitsz coordinatefromitsimageusing:
z
p
=
(D−r
p
(1−2δt))r
i
−2δtDr
p
D−r
i
. (5.10)
Estimationofmappingparameters
Inordertomapobjectpointstoimagepointsusing(5.5)and(5.9),weneedtocorrectly
estimateDandδθ. Additionally,weneedtoknowtheobjectspacetoimagespace(mm-
to-pixel) scaling factor, the camera orientation, and the orientation of the supporting
stage. We compute the scaling factor from the ratio of the pixel measurement of the
outer radius of the cone to the approximate true value,R, in mm. We, therefore, have
7 parameters to calibrate the system for: D,δθ, the outer radius of the cone,R, camera
angles,θ
xc
andθ
yc
,andstageangles,θ
xs
andθ
zs
,asshowninfigure5.4.
A 30 mm× 30 mm× 30 mm cube covered with a black-and-white checkerboard
patternisplacedonthestage. Thecorner-pointsonthecheckeredgridonthecubesur-
face, shown in figure 5.5, are used as feature points. The cube is positioned by roughly
49
Figure5.4: CalibrationparametersR,D,δθ,θ
xc
,θ
yc
,θ
xs
,andθ
zs
.
placing one of its faces along the transverse plane corresponding to the largest (outer-
most) circular cross-section of the cone with one of its bottom edges roughly aligned
with one of the edges of the stage (figure 5.4). The front vertex on this edge is used as
a reference point. Thexy coordinates, of this reference point, (x
ref
,y
ref
,z
ref
), (shown in
greeninfigures5.4and5.5)areestimatedfromthefrontviewofthecubeinthecamera
image, while itsz coordinate is computed from the corresponding image point on the
conicalmirrorfrom(5.10),allintermsofthepixel-to-distancescalingfactorsubjectto
calibration. The feature points were matched by performing a least squares fit on the 7
parametersfollowingthestepsdescribedbelow:
1. Chooseasetoffeaturepointsinobjectspace,(x
p;k
,y
p;k
,z
p;k
),wherek = 1 :N.
2. Transformthepoint-settoaccountforstagerotation:





x
ps;k
y
ps;k
z
ps;k





=R
x
(θ
xs
)R
z
(θ
zs
)





x
p;k
y
p;k
z
p;k





, (5.11)
50
Figure 5.5: Corner-points in the checkered grid on the cube surface used for feature
matching.
whereR
x
(θ
xs
)andR
z
(θ
zs
)arerotationmatrices.
3. Mapthetransformedpoint-settoimagespaceusing(5.1)and(5.9)togetasetof
imagepoints(x
is;k
,y
is;k
),wherek = 1 :N.
4. Identifythesetoffeaturepointsinimagespace,(x
i;k
,y
i;k
,z
i;k
),wherek = 1 :N.
5. Transformthepoint-settoaccountforcamerarotation:





x
ic;k
y
ic;k
z
ic;k





=R
x
(θ
xc
)R
y
(θ
yc
)





x
i;k
y
i;k
0





, (5.12)
whereR
x
(θ
xc
)andR
y
(θ
yc
)arerotationmatrices.
51
6. Use the Nelder-Mead simplex search method to compute the set of parameters
thatminimizethecostfunction:
Φ =
∑
k
(
(x
ic;k
−x
is;k
)
2
+(y
ic;k
−y
is;k
)
2
)
. (5.13)
Owingtotherefractionanddistortioncausedbythestage,welimitedoursetoffeature
pointstotheonesforwhichthestagedoesnotintercepttheraypaths. Figure5.6shows
the original cube, the CCD camera image, and the reconstructed checkerboard pattern
obtained using the optimal values of the 7 calibration parameters. Using these parame-
ters, we can now map points on an arbitrary object with a known height map onto the
imageplane.
Figure 5.6: (a) Photograph of the original checkeredcube, (b) CCD camera image, and
(c)digitalrenderingofthecubewiththereconstructedcheckerboardpattern.
52
Chapter6
PhantomStudies
6.1 Introduction
Asasteppingstonetoacquisitionandanalysisofinvivodata,wevalidatedthesystems
describedinchapter5byperformingsomepreliminaryphantomexperimentsusingboth
the pyramidal mirror setup with patterned illumination described in section 5.2 and the
conical mirror setup described in 5.3. The phantoms were designed to mimic tissue
opticalpropertiesandfluorescentsourceswereembeddedinsidethem. Theexperimen-
taldetailsandthereconstructionresultsarediscussedinthesubsequentsections.
6.2 CubicPhantom
Thefirstexperimentwasperformedusingthepyramidalmirrorsetup. Thephantomwas
cubic in shape and had an implanted fluorescent source. The details of this experiment
arediscussedbelow.
6.2.1 Phantomdescription
The cubic phantom (30 mm× 30mm× 30 mm) was composed of 1% Intralipid, 2%
agar, 20 µM bovine hemoglobin (H2625, Sigma-Aldrich Inc., St. Louis, MO), 1%
sodium azide (VW3465-2, VWR, West Chester, PA), and tri-buffered saline. Intralipid
was used to mimic tissue scattering, whereas hemoglobin was used to mimic tissue
53
Figure6.1: (a)White-lightimageshowingsimultaneousviewsoffoursidesofthecube
phantom. (b) Log-scale CCD image of the fluorescence pattern at a detection wave-
length of 710 nm due to the embedded source. A 6.5 x 6.5 mm square patch of the
bottomsurfacewasilluminatedatanexcitationwavelengthof635nm.
absorption. Additionally, hemoglobin is the dominant source of the wavelength depen-
dence of tissue optical properties. Sodium azide was added for phantom preservation
and prevention of hemoglobin oxygenation. This design ensures that the scattering and
absorption coefficients as well as spectral variation of these optical properties for the
phantomareclosetothosefortissue. A15mmlongcapillarytube(1mminnerdiame-
ter)filledwith6µMDiD(D307,InvitrogenCorporation)solutionwasembeddedinside
the phantom as the fluorescent target. The tube was placed 10 mm from and parallel to
each of two adjacent surfaces (the bottom and the right faces of the cube when placed
onthestage).
6.2.2 Dataacquisition
Multispectralfluorescencedatawasacquiredforatotalof68illuminationpatternsshone
on the top, bottom, left, and right faces of the cube. Each face was illuminated using
17 patterns. Of these, the first 16 patterns were 6.5 mm× 6.5 mm square patches of
54
Figure6.2: Setof17marchingsquarepatternsilluminatingthebottomsideofthecube.
Thefrontandbackfaceshavebeenremovedfromthecubesurfaceforeaseofviewing.
light stepped across the cube face, and the 17th pattern was full face illumination. Flu-
orescence data was collected for 10 excitation/emission filter pairs: 595/650, 595/670,
595/690, 595/710, 615/670, 615/690, 615/710, 635/690, 635/710 and 655/710 (all in
units of nm) for each of the 68 patterns. The entire dataset therefore consisted of 680
fluorescenceimages. Using4×4pixelbinning,thecameraexposuretimeswere5sand
55
Figure6.3: Setoffluorescencepatternscorrespondingtothe17patternsinfigure6.2for
the 635/710 nm excitation/emission pair. The front and back faces have been removed
fromthecubesurfaceforeaseofviewing.
1sforthesmallandlargesquarepatterns,respectively. Thesetof17patternsilluminat-
ing the bottom surface of the cube and the corresponding 17 fluorescence patterns for
the635/710nmexcitation/emissionpairareshowninfigures6.2and6.3respectively.
56
6.2.3 Forwardmodelcomputationandimagereconstruction
Theforwardproblemwassolvedusingthefiniteelementmethodwithvolumetessella-
tionforthe635/710nmexcitation/emissionpairandforthep = 68illuminationpatterns
described in section 6.2.2. A tessellated cube with 93,868 tessellation nodes was used
for computing the finite element model. The fluorescent source distribution was recon-
structedoverauniformgridwith15,625nodesandaspacingof1.2mm. Theexcitation
and emission forward model matrices, A
ex
andA
em
, were computed as described in
section 2.2 forn
d
= 3752 surface detector nodes andn
s
= 15,625 internal grid points.
Measuredvaluesoftheopticalpropertieswereusedforcomputingtheforwardmodels.
These values were: µ
a
(r,λ
ex
) = 0.0222,µ
′
s
(r,λ
ex
) = 1.2788,µ
a
(r,λ
em
) = 0.0.0084,
andµ
′
s
(r,λ
em
) = 1.1368,allinunitsofmm
−1
.
The size of the system matrix for this case is 255,136× 15,625. We used the for-
mulationin(2.14)toavoidcomputationofthefullsystemmatrix. Theregularizedleast
squares approach in (2.12) was used to solve the inverse problem. The gradient projec-
tion algorithm was used to minimize the cost function with a non-negativity constraint
[ACD
+
08]. TheinverseofthediagonalelementsoftheHessianmatrixofthecostfunc-
tion was used for preconditioning, a step that substantially accelerated the convergence
speed[ACD
+
08].
The reconstruction results are displayed in figure 6.4 both as transverse slices
through the axis of the fluorescent tube and as a three-dimensional isosurface. The
signal (99th percentile) to background (50th percentile) ratio in the reconstructed vol-
umewas82. Fromtheisosurface,thereconstructedfluorescentsourcewas17mmlong,
positioned11mmfromtheclosestcubefacesandhadanaveragediameterof4mm.
57
Figure 6.4: Reconstructed image displayed as (a) an isosurface using a 99th percentile
intensitythresholdfortheimagevolumeand(b)transverseslicesthroughtheaxisofthe
fluorescenttube.
6.3 Mouse-ShapedPhantom
The second experiment was performed using the conical mirror setup. The phantom
wasshapedlikeamouseandhadtwoimplantedfluorescentsources. Thedetailsofthis
experimentarediscussedbelow.
6.3.1 Phantomdescription
The mouse-shaped phantom (shown in figure 6.5 (a)) was made using a mold created
from a euthanized mouse. The phantom material consisted of 1% Intralipid, 2% agar,
and20µMhemoglobin. Twocapillarytubes, each20mmlong, 1mmindiameter, and
filled with 1µM DiD solution, were embedded inside the mouse-shaped phantom. The
embeddedlinesourcesarevisibleintheCTimage,showninfigure6.5(b).
58
Figure 6.5: Mouse-shaped phantom. (a) Full-light photograph. (b) A coronal and a
transversesectionfromtheCTimageshowingthetwotubesfilledwithfluorescentDiD
dye.
6.3.2 Dataacquisition
A 650 nm point laser source was used for excitation. Fluorescence data was collected
overtheentiresurfaceofthephantomatawavelengthof720nmusingtheconicalmirror
and the EMCCD camera. Data was acquired for 54 different illumination positions at
variouslocationsonthephantomsurface.
59
6.3.3 SurfacewarpingusingDigiWarp
Thesurfacegeometryofthemouse-shapedphantomwasextractedfromtheCTimage,
shown in figure 6.5 (b). In order to complete step 3 of figure 2.1, we used the Digi-
Warp technique [ACL
+
09, JCL
+
10]. This is a three-step procedure used to warp the
Digimouse atlas to fit the surface profile data. This method involves posture matching,
asymmetricL
2
distanceminimizationbetweenthetargetandtemplatesurfaces,andvol-
umetric registration using elastic energy minimization, as illustrated in figure 6.6. The
Figure 6.6: The three steps of the DigiWarp procedure: posture correction, surface
warpingbymeansofasymmetricL
2
distanceminimization,andvolumetricregistration
basedonelasticenergyminimization.
original Digimouse surface and the warped surface obtained by using DigiWarp on the
CTdataforthemouse-shapedphantomareshowninfigure6.7alongsideaphotoofthe
mouse-shapedphantom.
60
Figure6.7: Mouse-shapedphantom. (a)TransformationoftheDigimousesurfaceusing
the DigiWarp procedure (b) Photograph of the mouse-shaped phantom showing its sur-
facegeometryforreference.
6.3.4 Forwardmodelcomputationandimagereconstruction
The forward problem was solved using the finite element method with volume tessel-
lation for an excitation wavelength of 650 nm and an emission wavelength of 720 nm
for thep = 54 illumination patterns. The warped Digimouse atlas with 306,771 tetra-
hedrons and 58,244 tessellation nodes was used to generate the finite element model.
The fluorescent source distribution was reconstructed over a uniform grid with 11,036
nodes and a spacing of 1 mm. The excitation and emission forward model matrices,
A
ex
andA
em
, were computed forn
d
= 6762 surface detector nodes andn
s
= 11,036
61
internal grid points. The phantom was assumed to be optically homogeneous with
µ
a
(r,λ
ex
) = 0.0063,µ
′
s
(r,λ
ex
) = 0.88,µ
a
(r,λ
em
) = 0.0.009, andµ
′
s
(r,λ
em
) = 0.86,
allinunitsofmm
−1
.
We used the regularized least squares approach in (2.12) to solve the inverse prob-
lem. The gradient projection algorithm was used to minimize the cost function with a
non-negativity constraint. To speed up convergence, the inverse of the diagonal of the
Hessianmatrixwasusedforpreconditioning.
The reconstruction result is displayed in figure 6.8 as coronal slices of the phan-
tom volume. For validation, the FMT image was overlaid against the CT image and
the result is displayed in figure 6.9. The images were registered and displayed using
the AMIDE tool for viewing and analyzing medical imaging data (http://amide.
sourceforge.net).
Figure 6.8: Coronal sections of the mouse phantom showing the reconstruction of two
embedded line sources using the L
2
penalty with a regularization parameter of 10
−4
.
The slices are ordered from bottom to top starting top left in the image and ending
bottomright.
62
Figure 6.9: (a) Coronal and (b) transverse sections of the CT image of the mouse-
shaped phantom showing the two embedded fluorescent line sources. (c) Coronal and
(d) transverse overlay of CT and FMT images. (e) Coronal and (f) transverse sections
of the FMT image showing the two fluorescent line sources reconstructed using theL
2
penaltywitharegularizationparameterof10
−4
.
63
Chapter7
RegularizationUsingJointL
1
and
TotalVariationNormPenalties
7.1 Introduction
Inverse operators for ill-posed (or ill-conditioned) problems tend to be unbounded (or
have a very large norm). Regularization is the use of bounded (or low-norm) approx-
imations for such inverse operators that allow us to generate meaningful and numer-
ically stable solutions [DFM04]. This typically involves incorporating some prior
knowledge usually in the form of a penalty function that controls some desired prop-
erty of the unknown quantity or the prior probability distribution of an unknown ran-
dom variable. In the context of image reconstruction, for example, penalty func-
tions are often chosen to enforce smoothness [FH95], preserve edges [Lan90, SC03],
promote sparsity [MEHA07, CNJ07, GZ10a], or incorporate anatomical informa-
tion [LY91, FCR92, HMBN10]. The most widespread method for handling the ill-
conditioning of inverse problems is Tikhonov regularization [TA77], in which the cost
functioncontainsanL
2
normpenaltyterminadditiontothedata-fittingterm:
Φ(x) =
1
2
∥Ax−b∥
2
2
+
λ
2
∥Γx∥
2
2
, (7.1)
where Γ is an appropriately chosen Tikhonov matrix. In the simplest case, in which
Γ = I, this regularization technique tends to favor solutions with lower L
2
norms.
64
TheBayesianimplicationofTikhonovregularizationistheassumptionofamultivariate
Gaussianprobabilitydistributionontheunknownrandomvector.
While the simple Tikhonov regularizer generates satisfactory results when applied
to FMT reconstruction, we would like to design an improved regularizer that specifi-
cally takes into account prior knowledge about FMT images. A variety of regularizers
havebeenreportedinthefieldsofDOT,FMT,andBLT[HMBN10,PMP
+
99,ASAE07,
LBZ
+
05, GYIC05]. We focus on the fact that NIR probes used in FMT are designed
to preferentially accumulate in specific areas of interest, e.g., tumors or cancerous tis-
sue. Therefore reconstructed FMT images commonly exhibit fluorophore concentra-
tions localized within tumors and the major excretory organs [BLM
+
05, WTMB99].
Consequentlytheseimagestendtobeverysparsewithsomelocallysmoothhighinten-
sity regions. This inspires us to investigate a combination of regularizers that enforce
conditionsofsparsityandsmoothnessonreconstructedFMTimages.
The most fundamental sparsity metric is the L
0
norm, which is the total number
of non-zero elements in a vector. However, the problem with the L
0
norm penalty, at
least in the underdetermined case, is NP-hard [Nat95]. Instead, the L
1
norm penalty,
whichisaconvexrelaxationoftheL
0
normpenalty,isoftenusedtoenforcesparsityin
images and is particularly popular in the field of compressed sensing [Don06, CW08].
Approaches commonly used to solved the L
1
problem include iterative thresholding,
LASSO, and Basis Pursuit [DFM04, FNW07, BT09, Tib96, CDS98]. In this work, the
wechoosetheL
1
penaltyforenforcingsparsityinFMTimages. Avarietyofsmoothing
priorshavebeenappliedtotomographicreconstruction[HL89,Gre90]. Amongstthese,
priors involving quadratic penalties are particularly common, since they are easy to
handle. However,theytendtosmoothoutedgesinimages. We,therefore,choosetouse
the total variation (TV) penalty, which promotes smoothness while preserving edges
in images [ROF92]. In this work, we define the total variation as the L
1
norm of the
65
differencesbetweenneighboringpixels. ThisparticularformoftheTVpenaltyenforces
sparsity on pixel differences and consequently tends to generate images with piecewise
constant regions and sharp boundaries. A variety of methods have been employed for
handling the TV penalty, including dual-based and interior point approaches [Cha04,
VO98, HNW08, GZ10b]. In this work, we derive a method based on the Separable
ParaboloidalSurrogates(SPS)algorithmforminimizingtheTVpenalty.
Section 7.2 describes a compound approach that uses a combination of the SPS
method with the Preconditioned Conjugate Gradient (PCG) algorithm for handling the
joint L
1
and TV penalties. We use ordered subsets to accelerate the SPS approach.
In section 7.3 we describe an experimental setup using a mouse-shaped phantom for
testingthejointpenalties. Wevalidateourmethodbyapplyingittobothsimulatedand
experimentaldata. Adiscussionandanalysisofoursimulationandexperimentalresults
ispresentedinsection7.4.
7.2 Methods
Thecostfunctionweseektominimizecontainsthreeparts–adata-fittingterm,asparsi-
fyingpenaltyterm,andasmoothingpenaltyterm. Asmentionedbefore,thesparsifying
penaltyunderconsiderationistheL
1
normoftheunknownimagex. Forsmoothingthe
3Dimagex,wepenalizeitstotalvariation,definedasfollows:
TV(x) =
∑
k
|x
m(k)
−x
n(k)
|, (7.2)
66
wherem(k) andn(k) are pixel indices corresponding to thekth neighboring pixel pair
inthe3Dimagex. Theresultingoptimizationproblemisasfollows:
ˆ x = argmin
x≥0
1
2
∥Ax−b∥
2
2
+λ
L
1∥x∥
1
+λ
TV
∑
k
|x
m(k)
−x
n(k)
|. (7.3)
Here λ
L
1 and λ
TV
are regularization parameters. The non-differentiability of the
sparsity-enforcing L
1
norm penalty near zero makes it impossible to directly apply
gradient-basedmethodsforminimization. InFMT,however,wearehelpedbythenon-
negativity constraint on the unknown vector. We use a strictly positive initialization
for the reconstruction procedure and restrict the iterates to the non-negative quadrant.
This allows us to simplify the L
1
norm penalty to a linear term. To handle the non-
differentiabilityoftheTVpenalty,weapproximateitasfollows:
TV(x)≈ Φ
TV
(x),
∑
k
√
[Cx]
2
k
+δ
TV
, (7.4)
where the transformation v = Cx computes nearest neighbor differences and C ∈
R
nn×ns
, n
s
being the number of pixels inx and n
n
being the number of neighboring
pixel pairs. The quantityδ
TV
is an appropriately chosen thresholding parameter with a
smallvalue. Themodifiedcostfunctionis:
Φ(x) =
1
2
∥Ax−b∥
2
2
+λ
L
1x+λ
TV
∑
k
√
(x
m(k)
−x
n(k)
)
2
+δ
TV
. (7.5)
In the subsequent subsections, we describe two approaches for minimizing this cost
function.
67
7.2.1 PreconditionedConjugateGradientalgorithm
The Preconditioned Conjugate Gradient (PCG) method can used to minimize the cost
functionin(7.5). Thismethodusesthegradientgivenby:
∇Φ(x) =A
′
(Ax−b)+λ
TV
C
′
z +λ
L1
1. (7.6)
Herewedefinethevectorfunctionz as:
z(v) =
[
ξ(v
1
) ξ(v
2
) ... ξ(v
m
)
]
′
, (7.7)
wherethenotationξ(v)isusedtorepresentthescalarfunction:
ξ(v) =
1
√
v
2
+δ
TV
. (7.8)
To accelerate convergence, we use the inverse of the diagonal terms of the Hessian
A
′
A for preconditioning. The step size is determined using an Armijo line search. To
enforce non-negativity, every iterate is projected onto the non-negative quadrant. This
is achieved by computing a feasible direction within the non-negative orthant and an
optimalstep-sizealongthatdirection.
7.2.2 SeparableParaboloidalSurrogatesalgorithm
While the PCG algorithm is straightforward to implement, the computational time per
iterationisincreasedbytheexpensivelinesearch. Asanalternative,weexploreadiffer-
ent approach based on the optimization transfer principle, also known as Majorization-
Minimization(MM)[LHY00]. Ineveryiteration,wegenerateandminimizeasurrogate
function which satisfies a set of majorization conditions that guarantee monotonicity.
68
We design surrogates that are separable, using the SPS approach. Separable surrogate
functionsareeasiertominimizeandeliminatetheneedforanexpensivelinesearch.
Surrogatefunctionproperties
Optimizationtransferallowsustoreplaceacostfunction,Φ(x),thatisdifficulttomin-
imizebyaneasiertominimizesurrogatefunction,ϕ(x;x
n
)constructedateveryiterate,
x
n
. Theresultingupdateequationis:
x
n+1
= argmin
x∈D
ϕ(x;x
n
). (7.9)
This surrogate function must satisfy the following majorization conditions over the
domain,D,ofx:
ϕ(x
n
;x
n
) = Φ(x
n
) (7.10)
ϕ(x;x
n
) ≥ Φ(x), ∀x∈D (7.11)
∇
10
ϕ(x;x
n
)|
x=x
n = ∇Φ(x)|
x=x
n. (7.12)
Combining(7.11)and(7.12),weget:
Φ(x)−Φ(x
n
)≤ϕ(x;x
n
)−ϕ(x
n
;x
n
), ∀x∈D. (7.13)
Equation(7.13)ensuresthattheupdatein(7.9)allowsΦ(x)todecreasemonotonically.
69
SPSforthedata-fittingterm
We first construct a surrogate function that is convex, separable, and tangential to the
data-fitting cost function and lies above it using a method similar to [EF99]. We may
writethedata-fittingtermas:
Φ
DF
(x) =
1
2
(Ax−b)
′
(Ax−b)
=
1
2
∑
i
f
i
([Ax]
i
), (7.14)
where
f
i
([Ax]
i
), ([Ax]
i
−b
i
)
2
. (7.15)
Theterm[Ax]
i
canalsobewrittenas:
[Ax]
i
=
∑
j
a
ij
x
j
=
∑
j
a
ij
(x
j
−x
n
j
+x
n
j
)
=
∑
k
a
ik
x
n
k
n
∑
j=1
α
ij
+
∑
j
α
ij
a
ij
α
ij
(x
j
−x
n
j
)
=
∑
j
α
ij
(
a
ij
α
ij
(x
j
−x
n
j
)+[Ax
n
]
i
)
, (7.16)
where
∑
j
α
ij
= 1. (7.17)
Theα
ij
sareconstantsintroducedheretofacilitatethecomputationofthesurrogate. By
constraining these constants in (7.16) to satisfy (7.17) and by exploiting the convexity
off
i
definedin(7.15),wecancomputeafunctionthatliesabovef
i
:
f
i
([Ax]
i
)≤
∑
j
α
ij
f
i
(
a
ij
α
ij
(x
j
−x
n
j
)+[Ax
n
]
i
)
. (7.18)
70
Wepickα
ij
=a
ij
/
∑
k
a
ik
. Thisgivesusourseparablesurrogateforthedata-fitterm:
ϕ
DF
(x;x
n
) =
1
2
∑
i
(
∑
j
a
ij
∑
k
a
ik
f
i
(
a
ij
α
ij
(x
j
−x
n
j
)+[Ax
n
]
i
)
)
(7.19)
SPSforthetotalvariationpenaltyterm
The SPS function for the TV penalty is derived in two steps. We first approximate the
total variation penalty in (7.4) by a non-separable paraboloidal surrogate function as
follows:
ψ
n
(v) =
1
2
ξ(v
n
)v
2
. (7.20)
The function ξ(.) used here was defined in (7.8). We then use the convexity of the
functionψ
n
(v)toderiveaseparablesurrogateasfollows[LF95]:
ψ
n
(x
i(k)
−x
j(k)
) =ψ
n
(
1
2
[2x
i(k)
−x
n
i(k)
−x
j(k)
n]+
1
2
[2x
j(k)
−x
n
i(k)
−x
j(k)
n]
)
≤
1
2
ψ
n
(
2x
i(k)
−x
n
i(k)
−x
n
j(k)
)
+
1
2
ψ
n
(
2x
j(k)
−x
n
i(k)
−x
n
j(k)
)
.
(7.21)
TheresultantSPSfunctionfortheTVpenaltyis:
ϕ
TV
(x;x
n
) =
∑
k
1
2
ψ
n
(
2x
i(k)
−x
n
i(k)
−x
n
j(k)
)
+
∑
k
1
2
ψ
n
(
2x
j(k)
−x
n
i(k)
−x
n
j(k)
)
.
(7.22)
71
Updateequation
The update equation computes the minimum corresponding to the SPS cost function in
one step. Due to its separable nature, its Hessian is a diagonal matrix. It can be shown
thattheidealpreconditionerusestheinverseoftheHessiangivenby:
D
pre
(x) = diag
j
[
1
||A
′
Ae
j
||
1
+2λ
TV
(|C|
′
z(Cx))
j
]
. (7.23)
Owing to its paraboloidal (quadratic) nature, the minimum of the SPS function can be
computed in a single iteration. To ensure non-negativity, we project this result onto the
non-negativequadrant. Thus,thecompleteupdateequationis:
x
n+1
= [x
n
−D
pre
(x
n
)∇Φ(x
n
)]
+
. (7.24)
Thenotation[.]
+
usedin(7.24)representsprojectionontothenon-negativequadrant.
Orderedsubsetsimplementation
SincetheFMTinverseproblemtypicallyusesverylargedatasets,theSPSmethodcan
be accelerated using an ordered subsets (OS) (also known as the incremental gradient)
approach[Ber99,EF99]. TheFMTsystemmatrixistypicallyverytallandoftheform:
A =









A
em
D
ex
1
A
em
D
ex
2
.
.
.
A
em
D
ex
p









. (7.25)
72
ForOSimplementation,wedivideA
em
intosblocksorsubsets:
A
em
=









A
em
1
A
em
2
.
.
.
A
em
s









. (7.26)
Correspondingly,forpilluminationpatterns,theithblockofthesystemmatrixis:
A
i
=









A
em
i
D
ex
1
A
em
i
D
ex
2
.
.
.
A
em
i
D
ex
p









. (7.27)
Thgradientfortheithblockandforthenthiterationthenisgivenby:
g
i
(x) =∇r
i
(x)r
i
(x), (7.28)
where
r(x) =A
i
x−b
i
. (7.29)
Here, b
i
represents the block of the data vector corresponding to A
i
. The modified
updateequationfortheithblockandforthenthiterationthenis:
x
n;i+1
=
[
x
n;i
−α
n
D
pre
(x
n;i
)g
i
(x
n;i
)
]
+
, (7.30)
whereα
n
istherelaxationparameterorstepsize.
73
Figure7.1: ComparisonofconvergencecurvesofPCG,OSSPS,andOSSPS-PCGwith
thetransitionfromOSSPStoPCGoccurringafteriterationnumbers4,6,8,10,and12.
ThesecurveswereobtainedfortheL
1
-TVpenaltyandforarandominitialization.
7.2.3 Compoundapproach
We compared the convergence speeds of the SPS algorithm with ordered subsets
(OSSPS) with the PCG algorithm using a simulation setup that will be described in
section7.3.1. Theconvergencespeedishighlydependentontheinitialization. Thegen-
eral observationwasthat, for starting points veryfarfrom the true solution, the OSSPS
approachwasmuchfasterthan PCG.In thevicinity of thesolution, however, PCGwas
significantlyfaster. Thismotivatedustouseacompoundapproachwhereweuseafew
74
iterations of OSSPS to initialize the optimization problem and then use PCG to deter-
mine the final solution [LAL05]. Figure 7.1 demonstrates the speed-up achieved by
using the compound approach. We obtained convergence curves for OSSPS-PCG for
different transition points as shown in figure 7.1. Based on these curves, we chose to
switch between the methods after the 10th iteration. It should be noted that one com-
plete cycle through all the subsets is considered one iteration here. All convergence
plotswereobtainedforarandomlypickedpositivestartingpoint.
7.3 Results
We tested the L
1
-TV regularization method and compared it with the L
1
, TV, and L
2
penalties on simulated and experimental data. In this section, we describe the exper-
imental setup, present some simulations on a setup similar to the experimental setup,
evaluate the different penalties on the basis of some performance metrics, and, finally,
presentreconstructionresultsusingexperimentaldata.
7.3.1 Experimentalsetup
Weacquireddatausingtheexperimentalsetupdescribedinsection6.3. Theexperiment
descriptionisrepeatedhereforclarity. Weusedamouse-shapedphantomcreatedfrom
a mold (shown in figure 7.2 (a)). The phantom material consisted of 1% Intralipid, 2%
agar,and20µMhemoglobin. Twocapillarytubes,each20mmlong,1mmindiameter,
and filled with 1µM DiD solution, were embedded inside the mouse-shaped phantom,
asindicatedintheCTimageshowninfigure7.2(b).
A650nmpointlasersourcewasusedforexcitation. Fluorescencedatawascollected
over the entire surface of the phantom at a wavelength of 720 nm using the conical
mirrorsetupdescribedinsection5.3. Datawasacquiredforilluminationat54different
75
Figure 7.2: Mouse-shaped phantom. (a) Full-light photograph. (b) Coronal and trans-
versesectionsfromtheCTimageshowingthetwotubesfilledwithfluorescentDiDdye
insidethephantom.
points on the phantom surface. The surface geometry of the mouse-shaped phantom
was extracted from the CT image. The Digimouse atlas surface was warped to match
thephantomsurfaceusingtheDigiWarpprocedure.
7.3.2 Simulationresults
WesimulatedtwofluorescentlinesourcesinsidethewarpedDigimouseatlas. Theexci-
tation and emission wavelengths were assumed to be 650 nm and 720 nm respectively
– same as what was used in the experiment. The same 54 illumination locations from
the experiment were used for excitation. White Gaussian noise with an SNR of 5 was
addedtothefluorescencedatatomakethesimulationsrealistic. Coronalsectionsofthe
warpedmouseatlasindicatingthepositionsofthetwosimulatedlinesourcesareshown
infigure7.3.
76
Figure 7.3: Coronal sections of the mouse phantom showing the two simulated line
sources. The slices are ordered from bottom to top starting top left in the image and
endingbottomright.
Metric L
2
L
1
TV L
1
-TV
MSE 12.0879 11.1082 11.5229 11.7874
ROImean 0.0617 0.1770 0.1331 0.0919
ROIstandarddeviation 0.0148 0.0544 0.0380 0.0481
Backgroundmean 0.0120 0.0120 0.0130 0.0085
Backgroundstandarddeviation 0.0183 0.0346 0.0310 0.0232
Table7.1: PerformancemetricsforcomparisonoftheL
2
,L
1
,TV,L
1
-TVpenaltyfunc-
tions
The sources were reconstructed using theL
2
,L
1
, TV, andL
1
-TV penalties and the
results are shown in figures 7.4, 7.5, 7.6, and 7.7 respectively for indicated regulariza-
tion parameters. We computed the mean squared error (MSE), the mean and standard
deviation of the signal over the region of interest (ROI) defined by the two simulated
line sources, and finally the mean and standard deviation of the signal over the back-
ground pixels. Table 7.1 lists these of performance metrics for the different penalties.
TheL
1
,TV,andL
1
-TVpenaltieswereobservedtoyieldlowerMSEvaluesandstronger
meansignallevelsovertheROIthanL
2
. TheL
1
-TVapproachgeneratedtheleastmean
backgroundsignallevelandalowerbackgroundstandarddeviationthantheL
1
andTV
penaltiesindividually.
77
Figure 7.4: Coronal sections of the mouse phantom showing the two simulated line
sourcesreconstructedusingtheL
2
penaltywitharegularizationparameterof10
−3
.
Figure 7.5: Coronal sections of the mouse phantom showing the two simulated line
sourcesreconstructedusingtheL
1
penaltywitharegularizationparameterof10
−4
.
7.3.3 Experimentalresults
Wereconstructedtheexperimentaldatadescribedin7.3.1usingtheL
2
,L
1
,TV,andL
1
-
TV penalties and the results are shown in figures 7.8, 7.9, 7.10, and 7.11 respectively.
The result obtained using theL
1
-TV penalty is overlaid against the CT image and dis-
playedinfigure7.12. Thetwoimages wereregisteredanddisplayedusingtheAMIDE
tool for viewing and analyzing medical imaging data. We observe a reasonable degree
78
Figure 7.6: Coronal sections of the mouse phantom showing the two simulated line
sourcesreconstructedusingtheTVpenaltywitharegularizationparameterof2×10
−6
.
Figure 7.7: Coronal sections of the mouse phantom showing the two simulated line
sourcesreconstructedusingbothL
1
andTVpenaltieswithregularizationparametersof
5×10
−5
and2×10
−6
respectively.
of overlap between the reconstructed FMT image and the ground truth revealed by the
CTimage.
7.4 SummaryandDiscussion
We have used a combination of L
1
and TV penalties in the FMT inverse problem
to simultaneously enforce sparsity and smoothness in our reconstructed images. We
79
Figure 7.8: Coronal sections of the mouse phantom showing the reconstruction of two
embeddedlinesourcesusingtheL
2
penaltywitharegularizationparameterof10
−4
.
Figure 7.9: Coronal sections of the mouse phantom showing the reconstruction of two
embeddedlinesourcestheL
1
penaltywitharegularizationparameterof20.
derived inspiration from the fact that very often fluorescent probes and contrast agents
areconcentratedwithinspecificareasofinterest(e.g.,insidelocalizedtumorsandexcre-
toryorgans). WehavepresentedacompoundmethodthatusesacombinationofOSSPS
andPCGalgorithmstogenerateaconvergentsolution.
80
Figure7.10: Coronalsectionsofthemousephantomshowingthereconstructionoftwo
embeddedlinesourcesusingtheTVpenaltywitharegularizationparameterof0.5.
Wewereabletoapplyagradient-basedapproachinspiteofthenon-differentiability
of theL
1
norm penalty near zero by using a strictly positive initialization and restrict-
ing the iterates to the non-negative quadrant. This allowed us to simplify theL
1
norm
penalty to a linear term. SPS functions were derived for the data-fitting term and the
TV penalty. The MM algorithm guarantees monotonicity. However, the OS approach
that was incorporated to increase convergence speed does not guarantee convergence.
Typicallythisishandledbyusingadiminishingstepsize. Butourcompoundapproach,
wherethefinalsolutionisobtainedusingPCG,guaranteesconvergence. Thisapproach
was inspired by the observation that, far away from the solution, the SPS-OS method
hasafasterrateofconvergence,whilePCGisfasterinthevicinityofthesolution. The
advantage of the the OSSPS approach over PCG is that it eliminates the need for an
expensivelinesearchandrequiresfewernumbersofflopsperiteration. Ona2.93GHz
Intel
R ⃝
Core i7 CPU, it took about 40 s to complete one iteration of PCG while it took
about 20 s to cycle through all 5 subsets and complete one iteration of OSSPS. 5 data
blocks were used for the OS implementation, and the order of the blocks was chosen
randomlyatthestartoftheoptimizationprocedure.
81
Figure7.11: Coronalsectionsofthemousephantomshowingthereconstructionoftwo
embedded line sources using bothL
1
and TV penalties with regularization parameters
of10and1respectively.
Finally, we validated the derived procedure using simulated and experimental data
based on a mouse-shaped phantom with two embedded line sources. The simulation
results showed that theL
1
, TV, andL
1
-TV penalties generated lower MSE values and
stronger mean signal levels over the ROI than the L
2
regularizer. The L
1
-TV penalty
generated a smooth solution with a sparse background. TheL
1
-TV result had the least
mean background signal level and a lower background standard deviation than theL
1
andTVpenaltiesindividually. Finally,thereconstructedresultfortheexperimentaldata
obtainedusingtheL
1
-TVpenaltywasvalidatedbyoverlayingitagainstaCTimageof
thephantom.
82
Figure 7.12: (a) Coronal and (b) transverse sections of the CT image of the mouse-
shaped phantom showing the two embedded fluorescent line sources. (c) Coronal and
(d) transverse overlay of CT and FMT images. (e) Coronal and (f) transverse sections
oftheFMTimageshowingthetwofluorescentlinesourcesreconstructedusingbothL
1
andTVpenaltieswithregularizationparametersof10and0.7respectively.
83
Chapter8
In VivoStudies
8.1 Introduction
We performed an in vivo experiment using the conical mirror setup described in 5.3.1.
Fluorescence data was collected from a mouse with a xenograft tumor and was used to
reconstruct the fluorescent sources in 3D. The details of this experiment and the results
obtainedarepresentedthesubsequentsections.
8.2 Methods
Theexperimentalsetupandthedataacquisitionprocedurearediscussedbelow.
8.2.1 Experimentalsetup
A 30 g nude mouse with a xenograft melanoma (DX-6) dorsally located on the right
flankwasusedforthestudy. Thebulgeatthesiteofthetumormeasured 7mmacross.
The animal was anesthetized with 2% isoflurane and placed in the center of the field of
view of the conical mirror system. The mouse was injected intravenously with 10 nM
2-deoxyglucose (2-DG) conjugated with IRDye
R ⃝
800CW (LICOR, Lincoln, NE) 24 h
before imaging. The IRDye
R ⃝
800CW conjugated with 2-DG has an excitation peak
wavelengthof774nmandemissionpeakwavelengthof791nm.
84
8.2.2 Dataacquisition
A collimated laser beam with a spot size of 1 mm in diameter was used for illumi-
nating the mouse surface. The excitation wavelength was 785 nm. Data acquired by
illuminating the mouse at 48 locations spanning the entire mouse surface was used for
reconstruction. Thesepointsmappedontotheconicalmirrorsurfaceareindicatedinred
in figure 8.1. The full-light photograph obtained using the EMCCD camera also shows
the mouse on the stage inside the conical mirror system. For each illumination posi-
tion,afluorescenceimagewascollectedatanemissionwavelengthof820nm. Surface
topographic data for the mouse was acquired using the setup described in [LMD
+
09b].
After acquisition of fluorescence and surface profile data, the mouse was transferred to
amicro-CTscanner(whilestillundertheeffectofanethesia).
8.2.3 Forwardmodelingandimagereconstruction
Priortogeneratingtheforwardmodel,theDigiWarpmethodwasusedtowarptheDigi-
mouseatlaswith306,773tetrahedronsand58,244tessellationnodestothesurfacepro-
file data points following a procedure similar to that for the mouse-shaped phantom
describedin6.3.3.
The forward problem was solved using the finite element method with volume tes-
sellation for an excitation wavelength of 785 nm and an emission wavelength of 820
nmforthep = 48illuminationpatterns. Thefluorescentsourcedistributionwasrecon-
structed over a uniform grid with 11,232 nodes and a spacing of 1 mm. The excitation
andemissionforwardmodelmatrices,A
ex
andA
em
,werecomputedforn
d
= 7300sur-
face detector nodes andn
s
= 11,232 internal grid points. The detector node locations
85
Figure 8.1: Full-light photograph showing the anesthetized mouse placed on the stage
insidetheconicalmirrorsystem. Thedetectornodesonthewarpedmouseatlasmapped
onto the conical mirror surface are indicated in green. The illumination locations
mappedontothemirrorsurfaceareshowninred.
mapped to the conical mirror image are indicated in green in figure 8.1. The phan-
tom was assumed to be optically homogeneous withµ
a
(r,λ
ex
) =µ
a
(r,λ
em
) = 0.022,
µ
′
s
(r,λ
ex
) =µ
′
s
(r,λ
em
) = 1.41,bothinunitsofmm
−1
.
We used the regularized least squares approach in (2.12) to solve the inverse prob-
lem. The gradient projection algorithm was used to minimize the cost function with a
non-negativity constraint. To speed up convergence, the inverse of the diagonal of the
Hessianmatrixwasusedforpreconditioning.
86
8.3 Results
CoronalandtransversesectionsofthereconstructedFMTvolumeandtheCTimageof
the mouse are shown in figure 8.2. The coronal section from the CT shows the bulge
on the right flank indicating the tumor. The FMT image shows high intensity pixels at
this location as clearly revealed in the overlay image. The FMT image also shows high
intensitypixelswherethebladderislocatedowingtotheaccumulationofthedyeinthe
bladder24hafteritsinjection,whentheexperimentwasconducted.
87
Figure8.2: (a)Coronaland(b)transversesectionsofaCTimageofthemouseshowing
the tumor bulging on the right flank. (c) Coronal and (d) transverse overlay of CT and
FMT images. (e) Coronal and (f) transverse sections of the FMT image showing the
tumorontherightflankandfluorescencefromthebladder.
88
Chapter9
SummaryandFutureWork
9.1 SummaryofContributions
The research presented in this dissertation covers a number of aspects of FMT, some
theoretical and some practical. In course of this work, I put together a set of com-
putational tools for the component tasks described in section 2.1, I developed novel
theoretical approaches to address several challenges in the field of FMT, we collabora-
tivelydevelopedtwoprototypicFMTimagingsystems,and,finally,Isuccessfullyused
the computational toolbox to reconstruct phantom and in vivo data acquired using our
imagingsystems. Themaincontributionsofthisworkaresummarizedbelow:
9.1.1 Fastforwardsolver
We have developed a fast forward modeling scheme based on the Born approximation
[DAL
+
08b]. The performance of our method lies in between that of the slower but
moreaccuratefinite-elementmethodandfasterbutlessaccurateanalyticalmodels. Our
approach introduces tissue heterogeneity as perturbations in diffusion and absorption
coefficientsatrectangulargridpointsinsideamouseatlas. Thesereflectasacorrection
termaddedtothehomogeneousforwardmodel. Wehavetestedourmodelbyperform-
ing source localization studies first with a bioluminescence simulation setup and then
with an experimental setup using a fluorescent source embedded in an inhomogeneous
phantom that mimics tissue optical properties. This adds to our repertoire of forward
solversandthuscontributestostep5offigure2.1.
89
9.1.2 Illuminationpatternoptimization
We have developed a systematic approach for designing optimal illumination patterns
forFMT[DAL08a,DAJL09,DAJL10]. Weachievethisbyimprovingtheconditioning
of the Fisher information matrix. We parameterize the spatial illumination patterns and
formulateourproblemasaconstrainedoptimizationproblemthat,forafixednumberof
illumination patterns, yields the optimal set of patterns. For geometric insight, we used
our method to generate a set of 3 optimal patterns for an optically homogeneous, regu-
lar geometrical shape and observed expected symmetries in the result. We also gener-
atedoptimalilluminationpatternsfortheopticallyinhomogeneous,realisticallyshaped
Digimouseatlasfordifferentgivennumbersofpatterns. Theregularizedpseudoinverse
matrix generated using the singular value decomposition was employed to reconstruct
thepointspreadfunctionforeachsetofpatternsinthepresenceofasamplefluorescent
target deep inside the mouse atlas. We have further demonstrated the performance of
our method by examining the singular value spectra as well as plots of average spatial
resolutionvs. estimatorvariancecorrespondingtodifferentilluminationschemes.
9.1.3 RegularizationusingL
1
andTVnormpenalties
Owing to the ill-conditioned nature of the FMT inverse problem, image reconstruction
approachesrelyontheincorporationofpriorinformationintheformofsuitableregular-
izers. Veryoftenfluorescentprobesandcontrastagentsareconcentratedwithinspecific
areas of interest (e.g., inside localized tumors and excretory organs). Consequently,
the resultant reconstructed FMT images tend to be sparse yet locally smooth. We have
developed a regularization approach that takes this prior information into account and
applies a combination ofL
1
and TV penalties to the FMT inverse problem, the former
forthepurposeofenforcingsparsityandthelatterforthepurposeofgeneratingsmooth
90
andpiecewiseconstantreconstructedimages[DAL10]. TosolvetheL
1
-TVregularized
inverse problem, we have derived a compound approach that generates a convergent
solution by using a combination of the OSSPS and PCG algorithms. We validated our
techniquebyapplyingittosimulatedandexperimentalfluorescencedatafromamouse-
shapedphantomwithtwoembeddedfluorescentlinesources.
9.1.4 Imagingsystemdevelopmentandcalibration
In collaboration with the Cherry Lab at the University of California, Davis, and Cam-
bridge Research and Instrumentation Inc., we have designed and constructed two fluo-
rescence imaging systems capable of generating a variety of illumination patterns and
containing optics and filters designed for 360
◦
surface viewing and multispectral data
collection. One of these systems uses pyramidally arranged plane mirrors and is capa-
ble of generating computer-controlled spatial patterns for illumination using a DLP
R ⃝
-basedprojector[GDM
+
10,MDB
+
10,MGD
+
10]. Thesecondsystemfeaturesanovel
conical mirror geometry for full surface viewing [LMD
+
08, LMD
+
09b, LMD
+
09a].
These systems are capable for performing the tasks described in steps 1 and 2 of figure
2.1. Additionally, for the conical mirror-based imaging system, we have developed an
automated calibration scheme that allows us to map the fluorescence data from the 2D
CCDcameraimagesbacktotheanimalsurface. Thisisstep4offigure2.1.
9.1.5 Phantomand in vivostudies
Anumberofexperimentswereconductedusingtheimagingsystems. Thefirstphantom
study discussed here was performed using the pyramidal mirror setup [GDM
+
10]. The
studyinvolvedasimplecubicalphantommimickingtheopticalpropertiesofbiological
tissue with a capillary tube filled with fluorescent dye planted inside. Fluorescence
data was collected for several illumination patterns. The fluorescent line source was
91
reconstructedwithreasonableaccuracyusinganiterativeapproachundertheconstraint
ofnon-negativitywithoutexplicitcomputationofthesystemmatrix.
Thesecondstudyinvolvedamorerealisticmouse-shapedphantomwithtwoembed-
ded fluorescent tubes and was performed using the conical mirror setup [DAL10]. The
data was reconstructed using different regularizers: theL
2
norm, theL
1
norm, the TV
norm, and finally the jointL
1
-TV norm penalties. The results were validated by over-
layingandcomparingtheFMTreconstructionwithaCTimageofthephantom.
The final study was an in vivo study conducted using the conical mirror setup. The
studywasperformedonananesthetizedmousewithaxenograftmelanoma. Onceagain,
the data was reconstructed using an iterative approach and the result was validated by
overlayingagainstaCTimageofthediseasedmouse.
9.2 FutureWork
Inthecourseofthisproject,wehavedevelopedanumberofcomputationalmethodsfor
CWFMT.Wehavealsodesigned,built,andcalibratedtwoprototypicimagingsystems
and conducted phantom and in vivo studies using these setups. As future work, we
would like to incorporate some key innovations into our theoretical approaches as well
astheexperimentalprocedures.
9.2.1 Theoreticalgoals
FMT reconstruction results are highly dependent on the choice of the regularization
parameter, α. When α is increased, the reconstructed images look smoother and less
grainy,butwhenαistoolarge,significantlocalizationerrorsareintroduced. Thechoice
ofα, when working with experimental data, is often arbitrary. Also, because the sensi-
tivity varies significantly with source depth, a spatially-invariant regularization scheme
92
tends to generate reconstructed images with nonuniform spatial resolution properties.
Therefore, it is important to investigate spatially variant regularization approaches that
allow us to generate a desired point spread function and yield uniform spatial reso-
lution [SF00]. We would also like to investigate and document the properties of our
reconstruction methods through a combination of analytic and Monte Carlo simulation
studies.
9.2.2 Experimentalgoals
While developing the prototypic imaging systems and exploring their performance, it
becameclearthattherearenowidelyusedstandardsdefinedforFMTthatareappropri-
ateforcharacterizationofsystemperformanceintermsofresolution,quantitativeaccu-
racy, and sensitivity. In contrast, nuclear medicine systems use the NEMA (National
Electrical Manufacturers Association) standards [Nat07, Nat08, BNC
+
09] in order to
quantifyandcompareperformanceacrosssitesandmanufacturers. Itis,therefore,nec-
essary to develop a set of long-lived solid phantoms suitable for use in FMT and to
defineappropriateperformancemetricsforthem. Thefuturegoalforthisproject,onthe
experimental front, is to develop and use these phantoms to characterize commercially
availableCWFMTsystemsaswellasourownprototypicsystems. Wewillmakethese
phantomsanddetaileddesignspecificationsavailabletootherresearchersforcharacter-
izationoftheirimagingsystemsaswell.
93
Bibliography
[ACD
+
08] S.Ahn,A.J.Chaudhari,F.Darvas,C.A.Bouman,andR.MLeahy. Fast
iterative image reconstruction methods for fully 3D multispectral biolu-
minescencetomography. Phys.Med.Biol.,53(14):3921–3942,2008.
[ACL
+
09] A. Joshi A, A. J. Chaudhari, C. Li, D. Shattuck, J. Dutta, R. M. Leahy,
and A. W. Toga. Posture matching and elastic registration of a mouse
atlas to surface topography range data. In Biomedical Imaging: From
NanotoMacro,2009.ISBI’09.IEEEInternationalSymposiumon,pages
366–369,June28-July12009.
[ARC05] G. Alexandrakis, F. R. Rannou, and A. F. Chatziioannou. Tomographic
bioluminescence imaging by use of a combined optical-PET (OPET)
system: a computer simulation feasibility study. Phys. Med. Biol.,
50(17):4225–4241,2005.
[Aro95] R.Aronson. Boundaryconditionsfordiffusionoflight. J.Opt.Soc.Am.
A,12(11):2532–2539,1995.
[Arr99] S. Arridge. Optical tomography in medical imaging. Inverse Problems,
15:41–93,1999.
[ASAE07] J. Axelsson, J. Svensson, and S. Andersson-Engels. Spatially varying
regularization based on spectrally resolved fluorescence emission in flu-
orescence molecular tomography. Opt. Express, 15(21):13574–13584,
Oct.2007.
[ASHD93] S. R. Arridge, M. Schweiger, M. Hiroaka, and D. T. Delpy. A finite
element approach for modeling photon transport in tissue. Med. Phys.,
20:299–309,1993.
[Ber99] D.P.Bertsekas. NonlinearProgramming. AthenaScientific,1999.
94
[BLM
+
05] S. Bloch, F. Lesage, L. McIntosh, A. Gandjbakhche, K. Liang, and
S. Achilefu. Whole-body fluorescence lifetime imaging of a tumor-
targeted near-infrared molecular probe in mice. Journal of Biomedical
Optics,10(5):054003,2005.
[BNC
+
09] Q.Bao,D.Newport,M.Chen,D.B.Stout,andA.F.Chatziioannou. Per-
formance evaluation of the Inveon dedicated PET preclinical tomograph
basedonthenemanu-4standards. JNuclMed,50(3):401–408,2009.
[BT09] A. Beck and M. Teboulle. A fast iterative shrinkage-thresholding algo-
rithm for linear inverse problems. SIAM Journal on Imaging Sciences,
2(1):183–202,2009.
[CAL
+
09] A. J. Chaudhari, S. Ahn, R. Levenson, R. D. Badawi, S. R. Cherry, and
R. M. Leahy. Excitation spectroscopy in multispectral optical fluores-
cence tomography: methodology, feasibility and computer simulation
studies. Phys.Med.Biol.,54(15):4687–4704,2009.
[CB02] C.H.ContagandM.H.Bachmann.Advancesininvivobioluminescence
imaging of gene expression. Annu. Rev. Biomed. Eng., 4(1):235–260,
2002.
[CCH
+
03] J.P.Culver,R.Choe,M.J.Holboke,L.Zubkov,T.Durduran,A.Slemp,
V.Ntziachristos,B.Chance,andA.G.Yodh. Three-dimensionaldiffuse
optical tomography in the parallel plane transmission geometry: Evalu-
ation of a hybrid frequency domain/continuous wave clinical system for
breastimaging. Med.Phys.,30(2):235–247,2003.
[CDB
+
05] A. J. Chaudhari, F. Darvas, J. R. Bading, R. A. Moats, P. S. Conti, D. J.
Smith, S. R. Cherry, and R. M. Leahy. Hyperspectral and multispectral
bioluminescence optical tomography for small animal imaging. Phys.
Med.Biol.,50:5421–5441,2005.
[CDS98] S. S. Chen, D. L. Donoho, and M. A. Saunders. Atomic decomposition
by basis pursuit. SIAM Journal on Scientific Computing, 20(1):33–61,
1998.
[CDWW06] W. Cong, K. Durairaj, L. V. Wang, and G. Wang. A Born-type approxi-
mationmethodforbioluminescencetomography.Med.Phys.,33(3):679–
686,2006.
[Cha04] A. Chambolle. An algorithm for total variation minimization and appli-
cations. JournalofMathematicalImagingandVision,20:89–97,2004.
95
[CNJ07] N. Cao, A. Nehorai, and M. Jacobs. Image reconstruction for dif-
fuse optical tomography using sparsity regularization and expectation-
maximizationalgorithm. Opt.Express,15(21):13695–13708,2007.
[CPT
+
97] R.Cubeddu,A.Pifferi,P.Taroni,A.Torricelli,andG.Valentini. Asolid
tissuephantomforphotonmigrationstudies. Phys.Med.Biol.,42:1971–
1979,1997.
[CPW90] W. F. Cheong, S. A. Prahl, and A. J. Welch. A review of the optical
properties of biological tissues. IEEE J. Quantum Electron., 26:2166–
2185,1990.
[CW08] E.J.CandesandM.B.Wakin. Anintroductiontocompressivesampling.
SignalProcessingMagazine,IEEE,25(2):21–30,2008.
[DAJL09] J. Dutta, S. Ahn, A. A. Joshi, and R. M. Leahy. Optimal illumination
patterns for fluorescence tomography. In Biomedical Imaging: From
NanotoMacro,2009.ISBI’09.IEEEInternationalSymposiumon,pages
1275–1278,June2009.
[DAJL10] J. Dutta, S. Ahn, A. A. Joshi, and R. M. Leahy. Illumination pattern
optimizationforfluorescencetomography: theoryandsimulationstudies.
PhysicsinMedicineandBiology,55(10):2961,2010.
[DAL08a] J. Dutta, S. Ahn, and R. M. Leahy. Optimized illumination patterns for
fluorescencetomography.InWorldMolecularImagingConference,2008
(WMIC’08),Nice,France,September2008. PresentationNumber0345.
[DAL
+
08b] J. Dutta, S. Ahn, C. Li, A. J. Chaudhari, S. R. Cherry, and R. M. Leahy.
Computationally efficient perturbative forward modeling for 3D mul-
tispectral bioluminescence and fluorescence tomography. In Medical
Imaging2008: PhysicsofMedicalImaging,volume6913,page69130C.
SPIE,2008.
[DAL10] J. Dutta, S. Ahn, and R. M. Leahy. Fluorescence tomographic image
reconstruction using joint L1 and total variation norm penalties. In
World Molecular Imaging Conference, 2010, (WMIC’10), Kyoto,Japan,
September2010. PresentationNumber0384B.
[DDS03] D.Dudley,W.M.Duncan,andJ.Slaughter. Emergingdigitalmicromir-
ror device (DMD) applications. In MOEMS Display and Imaging Sys-
tems,volume4985,pages14–25.SPIE,2003.
[DFM04] I. Daubechies, M. De Friese, and C. De Mol. An iterative thresholding
algorithmforlinearinverseproblemswithasparsityconstraint. Commu-
nicationsinPureandAppliedMathematics,page14131457,2004.
96
[Don06] D. L. Donoho. Compressed sensing. Information Theory, IEEE Trans-
actionson,52(4),2006.
[DSCL07] B.Dogdas,D.Stout,A.Chatziioannou,andR.M.Leahy. Digimouse: A
3D whole body mouse atlas from CT and cryosection data. Phys. Med.
Biol.,52:577–587,2007.
[EF99] H. Erdogan and J. A. Fessler. Ordered subsets algorithms for transmis-
sion tomography. Physics in Medicine and Biology, 44(11):2835–2851,
1999.
[FCR92] J. A. Fessler, N.H. Clinthorne, and W.L. Rogers. Regularizedemission
imagereconstructionusingimperfectsideinformation. NuclearScience,
IEEETransactionson,39(5):1464–1471,Oct.1992.
[FH95] J. A. Fessler and III Hero, A. O. Penalized maximum-likelihood
imagereconstructionusingspace-alternatinggeneralizedEMalgorithms.
Image Processing, IEEE Transactions on, 4(10):1417–1429, October
1995.
[FNW07] M.A.T.Figueiredo,R.D.Nowak,andS.J.Wright. Gradientprojection
for sparse reconstruction: Application to compressed sensing and other
inverseproblems. SelectedTopicsinSignalProcessing,IEEEJournalof,
1(4):586–597,Dec.2007.
[GDM
+
10] C. Gardner, J. Dutta, G. S. Mitchell, S. Ahn, C.g Li, P. Harvey, R. Ger-
shman, S. Sheedy, J. Mansfield, S. R. Cherry, R. M. Leahy, and R. Lev-
enson. Improved in vivo fluorescence tomography and quantitation in
smallanimalsusinganovelmultiview,multispectralimagingsystem. In
BiomedicalOptics,OSATechnicalDigest(CD)(OpticalSocietyofAmer-
ica,2010).OSA,2010. PaperBTuF1.
[GHA05] A.P.Gibson,J.C.Hebden,andS.R.Arridge. Recentadvancesindiffuse
opticalimaging. Phys.Med.Biol.,50(4):R1–R43,2005.
[GM07] U.GrenanderandM.Miller.Patterntheory: fromrepresentationtoinfer-
ence. OxfordUniversityPress,USA,2007.
[GMW82] P.E.Gill,W.Murray,andM.H.Wright,editors. PracticalOptimization.
AcademicPress,1982.
[Gre90] P. J. Green. Bayesian reconstructions from emission tomography data
using a modified EM algorithm. Medical Imaging, IEEE Transactions
on,9(1):84–93,March1990.
97
[GRWN03] E. E. Graves, J. R., R. Weissleder, and V. Ntziachristos. A submillime-
ter resolution fluorescence molecular imaging system for small animal
imaging. Med.Phys.,30(5):901–911,2003.
[GSME05] A. Godavarty, E. M. Sevick-Muraca, and M. J. Eppstein. Three-
dimensional fluorescence lifetime tomography. Med. Phys., 32(4):992–
1000,2005.
[GYIC05] M. Guven, B. Yazici, X. Intes, and B. Chance. Diffuse optical tomog-
raphy with a priori anatomical information. Physics in Medicine and
Biology,50(12):2837,2005.
[GZ10a] H. Gao and H. Zhao. Multilevel bioluminescence tomography based
on radiative transfer equation part 1: L1 regularization. Opt. Express,
18(3):1854–1871,Feb2010.
[GZ10b] H. Gao and H. Zhao. Multilevel bioluminescence tomography based on
radiativetransferequationpart2: totalvariationandL1datafidelity. Opt.
Express,18(3):2894–2912,Feb2010.
[HAD97] J.CHebden,S.RArridge,andD.T.Delpy.Opticalimaginginmedicine:
I.experimentaltechniques. Phys.Med.Biol.,42(5):825–840,1997.
[HGY
+
02] J. C. Hebden, A. Gibson, R. M. Yusof, N. Everdell, E. M. C. Hillman,
D.T.Delpy,S.R.Arridge,T.Austin,J.H.Meek,andJ.S.Wyatt. Three-
dimensional optical tomography of the premature infant brain. Phys.
Med.Biol.,47(23):4155–4166,2002.
[Hie05] A. H. Hielscher. Optical tomographic imaging of small animals. Curr.
Opin.Biotechnol.,16(1):79–88,2005.
[HL89] T.HebertandR.M.Leahy. AgeneralizedEMalgorithmfor3-dbayesian
reconstruction from poisson data using gibbs priors. Medical Imaging,
IEEETransactionson,8(2):194–202,June1989.
[HLSM95] C. L. Hutchinson, J. R. Lakowicz, and E. M. Sevick-Muraca. Fluores-
cencelifetime-basedsensing intissues: a computational study. Biophys.
J.,68(4):1574–1582,1995.
[HMBN10] D. Hyde, E. L. Miller, D. H. Brooks, and V. Ntziachristos. Data specific
spatially varying regularization for multimodal fluorescence molecular
tomography. Medical Imaging, IEEE Transactions on, 29(2):365 –374,
Feb.2010.
98
[HNW08] Y. Huang, M. K. Ng, and Y.-W. Wen. A fast total variation minimiza-
tion method for image restoration. Multiscale Modeling & Simulation,
7(2):774–795,2008.
[Hor96] L.J.Hornbeck. DigitallightprocessingandMEMS:reflectingthedigital
display needs of the networked society. In Micro-Optical Technologies
forMeasurement,Sensors,andMicrosystems,volume2783,pages2–13.
SPIE,1996.
[HST
+
94] R. C. Haskell, L. O. Svaasand, T. T. Tsay, T. C. Feng, M. S. McAdams,
and B. J. Tromberg. Boundary conditions for the diffusion equation in
radiativetransfer. J.Opt.Soc.Am.A,11(10):2727–2740,1994.
[JBSM06] A. Joshi, W. Bangerth, and E. M. Sevick-Muraca. Non-contact fluo-
rescence optical tomography with scanning patterned illumination. Opt.
Express,14(14):6516–6534,2006.
[JCL
+
10] A. A. Joshi, A. J. Chaudhari, C. Li, J. Dutta, S. R. Cherry, D. W. Shat-
tuck,A.W.Toga,andR.M.Leahy. DigiWarp: amethodfordeformable
mouse atlas warping to surface topographic data. Physics in Medicine
andBiology,55(20):6197,2010.
[Jos08] A. A. Joshi. Geometric Methods for Image Registration and Analysis.
PhD thesis, University of Southern California, Los Angeles, California,
2008.
[Kay93] S. M. Kay. Fundamentals of Statistical Signal Processing: Estimation
Theory. PrenticeHall,EnglewoodCliffs,NJ,1993.
[KNBH02] A.D.Klose,U.Netz,J.Beuthan,andA.H.Hielscher. Opticaltomogra-
phyusingthetime-independentequationofradiativetransfer–part1: for-
wardmodel. J.Quant.Spectrosc.Radiat.Transf.,72(5):691–713,2002.
[KRD
+
08] A. T. N. Kumar, S. B. Raymond, A. K. Dunn, B. J. Bacskai, and D. A.
Boas. A time domain fluorescence tomography system for small animal
imaging. IEEETrans.Med.Imaging,27(8):1152–1163,Aug.2008.
[KYM
+
03] H. Koizumi, T. Yamamoto, A. Maki, Y. Yamashita, H. Sato,
H.Kawaguchi,andN.Ichikawa. Opticaltopography: Practicalproblems
andnewapplications. Appl.Opt.,42(16):3054–3062,2003.
[LAL05] Q. Li, S. Ahn, and R. M. Leahy. Fast hybrid algorithms for PET image
reconstruction. InNuclearScienceSymposiumConferenceRecord,2005
IEEE,volume4,page5pp.,2005.
99
[Lan90] K. Lange. Convergence of em image reconstruction algorithms with
Gibbs smoothing. Medical Imaging, IEEE Transactions on, 9(4):439–
446,Dec.1990.
[LBZ
+
05] A. Li, G. Boverman, Y. Zhang, D. Brooks, E. L. Miller, M. E. Kilmer,
Q. Zhang, E. M. C. Hillman, and D. A. Boas. Optimal linear inverse
solution with multiple priors in diffuse optical tomography. Appl. Opt.,
44(10):1948–1956,Apr.2005.
[LF95] K. Lange and J. A. Fessler. Globally convergent algorithms for max-
imum a posteriori transmission tomography. Image Processing, IEEE
Transactionson,4(10):1430–1438,October1995.
[LHY00] K. Lange, D. R. Hunter, and I. Yang. Optimization transfer using surro-
gateobjectivefunctions. JournalofComputationalandGraphicalStatis-
tics,9(1):1–20,2000.
[LM09] L. Liberti and N. Maculan, editors. Global Optimization: From Theory
toImplementation. SpringerUS,2009.
[LMD
+
08] C. Li, G. S. Mitchell, J. Dutta, S. Ahn, R. Leahy, and S. Cherry. A
conical mirror based multi-spectral three-dimensional fluorescence opti-
cal tomography system for small animal imaging. In Biomedical Engi-
neeringSociety(BMES)Annual FallMeeting,’08,St.Louis,MO,USA,
October2008.
[LMD
+
09a] C. Li, G. S. Mitchell, J. Dutta, S. Ahn, R. M. Leahy, and S. R. Cherry.
A high sensitivity multi-spectral three-dimensional fluorescence optical
tomographysystemforsmallanimalimaging.InSPIEBiomedicalOptics
Symposium (BiOS): Tomography and Spectroscopy of Tissue VIII, vol-
ume7174.SPIE,2009.
[LMD
+
09b] C.Li,G.S.Mitchell,J.Dutta,S.Ahn,R.M.Leahy,andS.R.Cherry. A
three-dimensional multispectral fluorescence optical tomography imag-
ing system for small animals based on a conical mirror design. Opt.
Express,17(9):7571–7585,2009.
[LMS09] V. Lukic, V. A. Markel, and J. C. Schotland. Optical tomography with
structuredillumination. Opt.Lett.,34(7):983–985,2009.
[LSRN08] T. Lasser, A. Soubret, J. Ripoll, and V. Ntziachristos. Surface recon-
structionforfree-space360
◦
fluorescencemoleculartomographyandthe
effects of animal motion. IEEE Trans. Med. Imaging, 27(2):188–194,
Feb.2008.
100
[LY91] R. M. Leahy and X. Yan. Incorporation of anatomical MR data for
improved functional imaging with PET. In Alan Colchester and David
Hawkes, editors, Information Processing in Medical Imaging, volume
511 of Lecture Notes in Computer Science, pages 105–120. Springer
Berlin/Heidelberg,1991.
[MDB
+
10] G.S.Mitchell,J.Dutta,D.Bulseco,C.Gardner,R.Gershman,P.Harvey,
S.Sheedy,J.R.Mansfield,R.M.Levenson,S.Ahn,R.M.Leahy,C.Li,
and S. R. Cherry. Improved in vivo fluorescence tomography and quan-
titation in small animals using a novel multiview, multispectral imag-
ing system. In Biophotonics Week, Optics Within Life Sciences, 2010,
(OWLS’10),Quebec,Canada,September2010.
[MEHA07] P. Mohajerani, A. A. Eftekhar, J. Huang, and A. Adibi. Optimal sparse
solutionforfluorescentdiffuseopticaltomography: theoryandphantom
experimentalresults. Appl.Opt.,46(10):1679–1685,2007.
[MG03] T. F. Massoud and S. S. Gambhir. Molecular imaging in living sub-
jects: seeing fundamental biological processes in a new light. Genes
Dev.,17(5):545–580,2003.
[MGD
+
10] J. R. Mansfield, C. Gardner, J. Dutta, G. Mitchell, S. Ahn, C. Li, P. Har-
vey,R.Gershman,S.Sheedy,S.R.Cherry,R.M.Leahy,andR.M.Lev-
enson. Improved in vivo fluorescence tomography and quantitation in
smallanimalsusinganovelmultiview,multispectralimagingsystem. In
World Molecular Imaging Conference, 2010, (WMIC’10), Kyoto,Japan,
September2010. PresentationNumber0366B.
[MGZ
+
07] H. Meyer, A. Garofalakis, G. Zacharakis, S. Psycharakis, C. Mamalaki,
D.Kioussis,E.N.Economou,V.Ntziachristos,andJ.Ripoll.Noncontact
opticalimaginginmicewithfullangularcoverageandautomaticsurface
extraction. Appl.Opt.,46(17):3617–3627,2007.
[MSO
+
04] A. B. Milstein, J. J. Stott, S. Oh, D. A. Boas, R. P. Millane, C. A.
Bouman, and K. J. Webb. Fluorescence optical diffusion tomography
using multiple-frequency data. J. Opt. Soc. Am. A, 21(6):1035–1049,
2004.
[Nat95] B. K. Natarajan. Sparse approximate solutions to linear systems. SIAM
JournalonComputing,24(2):227–234,1995.
[Nat07] National Electrical Manufacturers Association, Rosslyn, VA. NEMA
Standards NU 2-2007: Performance measurements of positron emission
tomographs(PETs),2007.
101
[Nat08] National Electrical Manufacturers Association, Rosslyn, VA. NEMA
Standards NU 4-2008: Performance measurements of small animal
positronemissiontomographs(PETs),2008.
[NC01] V.NtziachristosandB.Chance. Probingphysiologyandmolecularfunc-
tion using optical imaging: applications to breast cancer. Breast Cancer
Res.,3(1):41–46,2001.
[NRWW05] V. Ntziachristos, J. Ripoll, L. V. Wang, and R. Weissleder. Looking and
listening to light: the evolution of whole-body photonic imaging. Nat.
Biotechnol.,23(3):313–320,2005.
[PBAC05] S. Patwardhan, S. Bloch, S. Achilefu, and J. Culver. Time-dependent
whole-bodyfluorescencetomographyofprobebio-distributionsinmice.
Opt.Express,13(7):2564–2577,2005.
[PMP
+
99] B. W. Pogue, T. O. McBride, J. Prewitt, U. L.
¨
Osterberg, and K. D.
Paulsen. Spatiallyvariantregularizationimprovesdiffuseopticaltomog-
raphy. Appl.Opt.,38(13):2950–2961,May1999.
[PP06] B.W.PogueandM.S.Patterson. Reviewoftissuesimulatingphantoms
for optical spectroscopy, imaging and dosimetry. J. Biomed. Opt., 11,
2006.
[Puk93] F. Pukelsheim. Optimal Design of Experiments. Wiley series in proba-
bilityandmathematicalstatistics.JohnWiley&Sons,1993.
[QBM06] A. Qiu, D. Bitouk, and M. I. Miller. Smooth functional and structural
maps on the neocortex via orthonormal bases of the Laplace-Beltrami
operator. IEEETrans.Med.Imaging,25(10):1296–1306,2006.
[RCN01] B. W. Rice, M. D. Cable, and M. B. Nelson. In vivo imaging of light-
emittingprobes. J.Biomed.Opt.,6:432–440,2001.
[RNCNV01] J. Ripoll, V. Ntziachristos, R. Carminati, and M. Nieto-Vesperinas.
Kirchhoffapproximationfordiffusivewaves. Phys.Rev.E,64,2001.
[RNVWN02] J. Ripoll, M. Nieto-Vesperinas, R. Weissleder, and V. Ntziachristos.
Fast analytical approximation for arbitrary geometries in diffuse optical
tomography. Opt.Lett.,27(7):527–529,2002.
[ROF92] L.I.Rudin,S.Osher,andE.Fatemi.Nonlineartotalvariationbasednoise
removal algorithms. Physica D: Nonlinear Phenomena, 60(1-4):259–
268,1992.
102
[RSM01] R. Roy and E. M. Sevick-Muraca. Three-dimensional unconstrained
andconstrainedimage-reconstructiontechniquesappliedtofluorescence,
frequency-domainphotonmigration. Appl.Opt.,40:2206–2215,2001.
[RSN03] J. Ripoll, R. B. Schulz, and V. Ntziachristos. Free-space propagation of
diffuse light: Theory and experiments. Phys. Rev. Lett., 91(10):103901,
Sep2003.
[RXK06] B. W. Rice, H. Xu, and C. Kuo. Surface construction using combined
photographicandstructuredlightinformation. 2006. PatentApplication
0268153.
[SAHD95] M.Schweiger,S.R.Arridge,M.Hiraoka,andD.T.Delpy. Thefiniteele-
ment method for the propagation of light in scattering media: Boundary
andsourceconditions. Med.Phys.,22(11):1779–1792,1995.
[SC03] D. Strong and T. Chan. Edge-preserving and scale-dependent properties
oftotalvariationregularization. InverseProblems,19(6):S165,2003.
[SCS
+
02] D. Stout, P. Chow, R. Silverman, R. M. Leahy, X. Lewis, S. Gambhir,
and A. Chatziioannou. Creating a whole body digital mouse atlas with
PET,CTandcryosectionimages. Molec.Imag.Biol.,4(4)):S27,2002.
[SF00] J. W. Stayman and J. A. Fessler. Regularization for uniform spatial res-
olution properties in penalized-likelihood image reconstruction. IEEE
Trans.Med.Imaging,19(6):601–615,June2000.
[SN06] A. Soubret and V.Ntziachristos. Fluorescence molecular tomographyin
the presence of background fluorescence. Phys. Med. Biol., 51:3983–
4001,2006.
[SRL
+
09] X. Shu, A. Royant, M. Z. Lin, Todd A. Aguilera, V. Lev-Ram, P. A.
Steinbach, and R. Y. Tsien. Mammalian Expression of Infrared Flu-
orescent Proteins Engineered from a Bacterial Phytochrome. Science,
324(5928):804–807,2009.
[SZ05] E. Scherleitner and B. G. Zagar. Optical tomography imaging based on
higherorderBornapproximationofdiffusephotondensitywaves. IEEE
Trans.Instrum.Meas.,54(4):1607–1611,2005.
[TA77] A.N.TikhonovandV.Y.Arsenin. SolutionsofIll-PosedProblems. John
Wiley&Sons,Washington,D.C.,1977.
[Tib96] R.Tibshirani. RegressionshrinkageandselectionviatheLasso. Journal
of the Royal Statistical Society. Series B (Methodological), 58(1):267–
288,1996.
103
[TJMSR04] T. Troy, D. Jekic-McMullen, L. Sambucetti, and B. Rice. Quantitative
comparisonofthesensitivityofdetectionoffluorescentandbiolumines-
centreportersinanimalmodels. Mol.Imaging,3(1):9–23,2004.
[VC97] A. Villringer and B. Chance. Non-invasive optical spectroscopy and
imaging of human brain function. Trends Neurosci., 20(10):435 – 442,
1997.
[VO98] C. R. Vogel and M. E. Oman. Fast, robust total variation-based recon-
structionofnoisy,blurredimages. ImageProcessing,IEEETransactions
on,7(6):813–824,June1998.
[WN03] R. Weissleder and V. Ntziachristos. Shedding light onto live molecular
targets. NatureMed.,9:123–128,2003.
[WTMB99] R. Weissleder, C.-H. Tung, U. Mahmood, and A. Bogdanov. In
vivoimagingoftumorswithprotease-activatednear-infraredfluorescent
probes. NatureBiotechnology,17(4):375–378,1999.
[ZKS
+
05] G. Zacharakis, H. Kambara, H. Shih, J. Ripoll, J. Grimm, Y. Saeki,
R. Weissleder, and V. Ntziachristos. Volumetric tomography of fluores-
cent proteins through small animals in vivo. Proc. Natl. Acad. Sci. USA,
102(51):1825218257,2005.
[ZSA
+
09] A. D. Zacharopoulos, P. Svenmarker, J. Axelsson, M. Schweiger, S. R.
Arridge,andS.Andersson-Engels. Amatrix-freealgorithmformultiple
wavelength fluorescence tomography. Opt. Express, 17(5):3042–3051,
2009.
[ZVM
+
06] G.Zavattini,S.Vecchi,G.Mitchell,U.Weisser,R.M.Leahy,B.J.Pich-
ler, D. J. Smith, and S. R. Cherry. A hyperspectral fluorescence system
for3dinvivoopticalimaging. Phys.Med.Biol.,51(8):2029–2043,2006.
104

Abstract (if available)

Abstract Fluorescence molecular tomography is an imaging modality which exploits the specificity of fluorescent biomarkers to generate volumetric images using near-infrared light, which is safe and non-ionizing, at a substantially lower cost than competing modalities like positron emission tomography (PET) or single photon emission computed tomography (SPECT). While these very attractive features make it well-suited for preclinical research, the 3D reconstruction of fluorescent sources is confounded by the high degree of absorption and scattering of photons propagating through tissue, which make the inverse problem ill-posed. The focus of my Ph.D. research has been the development of computational techniques to tackle this challenge in order to reconstruct high resolution 3D images of molecular targets in vivo in small animals. To accurately and efficiently predict the photon fluence on the animal surface, we have developed a fast approach based on the Born approximation for computing the photon propagation model with reasonable accuracy. We have also developed a systematic approach for designing efficient spatial patterns of illumination for fluorescence tomography. The key idea behind our approach is to formulate a constrained optimization problem to improve the conditioning of the Fisher information matrix for the fluorescence forward model. The ill-conditioned nature of the inverse problem warrants the incorporation of some prior information to improve the conditioning of the problem and make it solvable. We have developed a regularization approach that incorporates this information as a combination of sparsity- and smoothness-enforcing penalty functions.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Hyperspectral and multispectral optical bioluminescence and fluorescence tomography in small animal imaging

PDF

Fast solvers in hyperspectral optical bioluminescence tomography for small animal imaging

PDF

Fast iterative image reconstruction for 3D PET and its extension to time-of-flight PET

PDF

Transmission tomography for high contrast media based on sparse data

PDF

New theory and methods for accelerated MRI reconstruction

PDF

The wavelet filter: a new method for nonlinear, nongaussian filtering in one dimension

PDF

Random access to compressed volumetric data

PDF

Statistical lesion detection in dynamic positron emission tomography

PDF

Signal processing methods for interaction analysis of functional brain imaging data

PDF

A treatise on cascaded computer generated holograms

PDF

Time-resolved laser induced fluorescence spectroscopy for detection of brain tumors

PDF

Advanced liquid simulation techniques for computer graphics applications

ZIP

Development of a multi-mode optical imaging system for preclinical applications in vivo [program]

PDF

Multi-softcore architectures and algorithms for a class of sparse computations

PDF

Shift-invariant autoregressive reconstruction for MRI

PDF

Model based PET image reconstruction and kinetic parameter estimation

PDF

3D inference and registration with application to retinal and facial image analysis

PDF

Correspondence-based analysis of 3D deformable shapes: methods and applications

PDF

Molecular imaging data grid (MIDG) for multi-site small animal imaging research based on OGSA and IHE XDS-i

PDF

Contributions to structural and functional retinal imaging via Fourier domain optical coherence tomography

Asset Metadata

Creator Dutta, Joyita (author)

Core Title Computational methods for fluorescence molecular tomography

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Electrical Engineering

Publication Date 04/26/2011

Defense Date 03/11/2011

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag cancer,fluorescence,imaging,near infrared,OAI-PMH Harvest,optical,preclinical,tomography

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Leahy, Richard M. (committee chair), Jadvar, Hossein (committee member), Jenkins, B. Keith (committee member)

Creator Email Dutta.Joyita@mgh.harvard.edu,joyita.dutta@gmail.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-m3773

Unique identifier UC1417997

Identifier etd-Dutta-4387 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-445926 (legacy record id),usctheses-m3773 (legacy record id)

Legacy Identifier etd-Dutta-4387.pdf

Dmrecord 445926

Document Type Dissertation

Rights Dutta, Joyita

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Repository Name Libraries, University of Southern California

Repository Location Los Angeles, California

Repository Email cisadmin@lib.usc.edu