Harris Manchester College Oxford University Oxford, England

17 Pages·2016·1.03 MB·English

by  Antonietta Iadisernia Corsano

Galerkin-based meshless methods for
photon transport in the biological tissue
ChenghuQin,JieTian∗,XinYang,KaiLiu,GuoruiYan,Jinchao
Feng,YujieLv,MinXu
MedicalImageProcessingGroup
KeyLaboratoryofComplexSystemsandIntelligenceScience
InstituteofAutomationChineseAcademyofSciences
P.O.Box2728,Beijing,100190,China
[email protected]
Abstract: As an important small animal imaging technique, optical
imaging has attracted increasing attention in recent years. However, the
photonpropagationprocessis extremelycomplicatedfor highly scattering
propertyofthebiologicaltissue.Furthermore,thelighttransportsimulation
intissuehasasignificantinfluenceoninversesourcereconstruction.Inthis
contribution, we present two Galerkin-based meshless methods (GBMM)
to determine the light exitance on the surface of the diffusive tissue. The
twomethodsarebothbasedonmovingleastsquares(MLS)approximation
which requires only a series of nodes in the region of interest, so compli-
catedmeshingtaskcanbeavoidedcomparedwiththefiniteelementmethod
(FEM).Moreover,MLSshapefunctionsarefurthermodifiedto satisfythe
deltafunctionpropertyinonemethod,whichcansimplifytheprocessingof
boundaryconditionsincomparisonwiththeother.Finally,theperformance
of the proposed methods is demonstrated with numerical and physical
phantomexperiments.
© 2008 OpticalSocietyofAmerica
OCIScodes:(170.3660)Lightpropagationintissues;(170.3880)Medicalandbiologicalimag-
ing;(170.5280)Photonmigration
Referencesandlinks
1. V.Ntziachristos,J.Ripoll,L.V.Wang,andR.Weissleder,“Lookingandlisteningtolight:theevolutionofwhole
bodyphotonicimaging,”Nat.Biotechnol.23,313-320(2005).
2. G.Wang,W.Cong,H.Shen,X.Qian,M.Henry,andY.Wang,“Overviewofbioluminescencetomography–a
newmolecularimagingmodality,”Front.Biosci.13,1281-1293(2008).
3. S.BhaumikandS.S.Gambhir,“OpticalimagingofRenillaluciferasereportergeneexpressioninlivingmice,”
Proc.Natl.Acad.Sci.USA99,377-382(2002).
4. T.F.MassoudandS.S.Gambhir,“Molecularimaginginlivingsubjects:seeingfundamentalbiologicalprocesses
inanewlight,”GenesDev.17,545-580(2003).
5. W.Rice,M.D.Cable,andM.B.Nelson,“Invivoimagingoflight-emittingprobes,”J.Biomed.Opt.6,432-440
(2001).
6. E.E.Graves,J.Ripoll,R.Weissleder,andV.Ntziachristos,“Asubmillimeterresolutionfluorescencemolecular
imagingsystemforsmallanimalimaging,”Med.Phys.30,901-911(2003).
7. A. Cong and G. Wang, “A finite-element-based reconstruction method for 3D fluorescence tomography,”
Opt. Express 13, 9847-9857 (2005), http://www.opticsinfobase.org/abstract.cfm?URI=
oe-13-24-9847.
8. C. Contag and M. H. Bachmann, “Advances in bioluminescence imaging of gene expression,” Annu. Rev.
Biomed.Eng.4,235-260(2002).
#102140 - $15.00 USD Received 29 Sep 2008; revised 3 Nov 2008; accepted 11 Nov 2008; published 24 Nov 2008
(C) 2008 OSA 8 December 2008 / Vol. 16, No. 25 / OPTICS EXPRESS 20317
9. G. Wang, H. Shen, W. Cong, S. Zhao, and G. W. Wei, “Temperature-modulated bioluminescence tomog-
raphy,” Opt. Express 14, 7852-7871 (2006), http://www.opticsinfobase.org/abstract.cfm?
URI=oe-14-17-7852.
10. V.Y.Soloviev,“Tomographicbioluminescenceimagingwithvaryingboundaryconditions,”Appl.Opt.46,2778-
2784(2006),http://www.opticsinfobase.org/abstract.cfm?URI=ao-46-14-2778.
11. Y.Lv,J.Tian,W.Cong,G.Wang,W.Yang,C.Qin,andM.Xu,“Spectrallyresolvedbioluminescencetomogra-
phywithadaptivefiniteelementanalysis:methodologyandsimulation,”Phys.Med.Biol.52,4497-4512(2007).
12. A.P.Gibson,J.C.Hebden,andS.R.Arridge,“Recentadvancesindiffuseopticalimaging,”Phys.Med.Biol.
50,R1-R43(2005).
13. W. Cong, A. Cong, H. Shen, Y. Liu, and G. Wang, “Flux vector formulation for photon propagation in the
biological tissue,”Opt.Lett.32,2837-2839(2007),http://www.opticsinfobase.org/abstract.
cfm?URI=ol-32-19-2837.
14. Y.Lv,J.Tian,W.Cong,G.Wang,J.Luo,W.Yang,andH.Li,“Amultileveladaptivefiniteelementalgorithm
forbioluminescencetomography,” Opt.Express14,8211-8223(2006),http://www.opticsinfobase.
org/abstract.cfm?URI=oe-14-18-8211.
15. A. Joshi, W. Bangerth, and E. Sevick-Muraca, “Adaptive finite element based tomography for fluorescence
optical imaging intissue,” Opt. Express 12, 5402-5417 (2004), http://www.opticsinfobase.org/
abstract.cfm?URI=oe-12-22-5402.
16. A.Joshi,W.Bangerth,A.B.Thompson,andE.M.Sevick-Muraca,“Adaptivefiniteelementmethodsforfluo-
rescenceenhancedfrequencydomainopticaltomography:forwardimagingproblem,”IEEEInternationalSym-
posiumonBiomedicalImaging(ISBI2004)2,1103-1106(2004).
17. W.Cong,D.Kumar,Y.Liu,A.Cong,andG.Wang,“Apracticalmethodtodeterminethelightsourcedistribution
inbioluminescentimaging,”Proc.SPIE5535,679-686(2004).
18. L.H.Wang,S.L.Jacques,andL.Q.Zheng,“MCML-MonteCarlomodelingofphotontransportinmulti-layered
tissues,”Comput.Meth.Prog.Biomed.47,131-146(1995).
19. D.Boas,J.Culver,J.Stott,andA.Dunn,“ThreedimensionalMonteCarlocodeforphotonmigrationthrough
complex heterogeneous media including the adult human head,” Opt. Express 10, 159-169 (2002), http:
//www.opticsinfobase.org/abstract.cfm?URI=oe-10-3-159.
20. H.Li,J.Tian,F.Zhu,W.Cong,L.V.Wang,E.A.Hoffman,andG.Wang,“Amouseopticalsimulationenviron-
ment(MOSE)toinvestigatebioluminescentphenomenainthelivingmousewithMonteCarlomethod,”Acad.
Radiol.11,1029-1038(2004).
21. W.Cong, G.Wang,D.Kumar,Y.Liu,M.Jiang, L.V.Wang,E.A.Hoffman, G.McLennan, P.B. McCray,
J.Zabner,andA.Cong,“Practicalreconstructionmethodforbioluminescencetomography,” Opt.Express13,
6756-6771(2005),http://www.opticsinfobase.org/abstract.cfm?URI=oe-13-18-6756.
22. Y.Lv,J.Tian,H.Li,J.Luo,W.Cong,G.Wang,andD.Kumar,“Modelingtheforwardproblembasedonthe
adaptiveFEMsframeworkinbioluminescencetomography,”Proc.SPIE6318,63180I(2006).
23. S.R.Arridge,H.Dehghani,M.Schweiger,andE.Okada,“Thefiniteelementmodelforthepropagationoflight
inscatteringmedia:Adirectmethodfordomainswithnonscatteringregions,”Med.Phys.27,252-264(2000).
24. R.H.Bayford,A.Gibson,A.TizzardA,T.Tidswell,andD.S.Holder,“Solvingtheforwardprobleminelectrical
impedancetomographyforthehumanheadusingIDEAS(integrateddesignengineeringanalysissoftware),a
finiteelementmodellingtool,”Physiol.Meas.22,55-64(2001).
25. S. J. Koopman, A. C. Harvey, J.A. Doornick, and N. Shephard, Stamp 5.0:structural time series analyser,
modellerandpredictor,(TheManual.Chapman&Hall,London,1995).
26. I.V.Singh,K.Sandeep,andR.Prakash,“TheelementfreeGalerkinmethodinthreedimensionalsteadystate
heatconduction,”Int.J.Comput.Eng.Sci.3,291-303(2002).
27. I.V.Singh,“ParallelimplementationoftheEFGmethodforheattransferandfluidflowproblems,”Adv.Eng.
Software34,453-463(2004).
28. T.Belytschko,L.Gu,andY.Y.Lu,“Fractureandcrackgrowthbyelement-freeGalerkinmethods,”Modelling
Simul.Mater.Sci.Eng.2,519-534(1994).
29. G.Alexandrakis,F.R.Rannou,andA.F.Chatziioannou, “Tomographicbioluminescenceimagingbyuseofa
combinedoptical-PET(OPET)system:acomputersimulationfeasibilitystudy,”Phys.Med.Biol.50,4225-4241
(2005).
30. M.Schweiger,S.R.Arridge,M.Hiraoka,andD.T.Delpy,“Thefiniteelementmethodforthepropagationof
lightinscatteringmedia:Boundaryandsourceconditions,”Med.Phys.22,1779-1792(1995).
31. W.G.EganandT.W.Hilgeman,opticalpropertiesofinhomogeneousmaterials,(Academic,NewYork,1979).
32. T.Belytschko,Y.Y.Lu,andL.Gu,“Element-freeGalerkinmethod,”Int.J.Numer.MethodsEng.37,229-256
(1994).
33. J.DolbowandT.Belytschko, “Anintroduction toprogrammingthemeshlesselementfreeGalerkinmethod,”
Arch.Comput.MethodsEng.5,207-241(1998).
34. X.ZhangandY.Liu,Meshlessmethods,(TsinghuaUniversityPress,Beijing,2004).
35. J.S.ChenandH.P.Wang,“Newboundaryconditiontreatmentsinmeshfreecomputationofcontactproblems,”
Comput.MethodsAppl.Mech.Eng.187,441-468(2000).
#102140 - $15.00 USD Received 29 Sep 2008; revised 3 Nov 2008; accepted 11 Nov 2008; published 24 Nov 2008
(C) 2008 OSA 8 December 2008 / Vol. 16, No. 25 / OPTICS EXPRESS 20318
36. S.Li,W.Hao,andW.K.Liu,“Numericalsimulationsoflargedeformationofthinshellstructuresusingmeshfree
methods,”Comput.Mech.25,102-116(2000).
37. S.R.Arridge, M.Schweiger, M.Hiraoka, andD.T.Delpy, “Afinite elementapproach formodelingphoton
transportintissue,”Med.Phys.20,299-309(1993).
38. J.Scho¨berl,“Netgenanadvancingfront2D/3D-meshgeneratorbasedonabstractrules,”Comput.Visual.Sci.1,
41-52(1997).
39. D.Qin,H.Zhao,Y.Tanikawa,andF.Gao,“Experimentaldeterminationofopticalpropertiesinturbidmedium
byTCSPCtechnique,”Proc.SPIE6434,64342E(2007).
1. Introduction
Molecular imaging is a very promising and rapidly developing biomedical research field in
whichthemoderntoolsandmethodsarebeingmarriedtorepresentinvivocellularandmolec-
ularprocessesdirectly,sensitively,non-invasivelyandspecifically,suchasmonitoringprotein-
proteininteractions,geneexpression,celltraffickingandengraftment[1,2].Amongmolecular
imagingmodalities,opticalimaging,especiallyfluorescenceandbioluminescenceimaging,has
becomearesearchfocusoverthepastyearsforitsexcellentperformance,non-radiativityand
highcost-effectivenesscomparedwithtraditionalimagingtechniqueslikeX-raycomputedto-
mography(CT),magneticresonanceimaging(MRI),positronemissiontomography(PET)and
singlephotonemissioncomputedtomography(SPECT)[3,4,5].Influorescencetechnology,
thefluorophoreprobeinsidethetissueabsorbstheincidentexcitationphotonsproducedbythe
externallightsource,andthendecaystoitsgroundstate,emittingthephotonssynchronously
[6,7].Whereasbioluminescenceimagingemploysluciferaseenzymes,whichcancatalyzethe
biochemical reactions of substrate luciferin to generate bioluminescent photons in the pres-
ence of oxygen,ATP and Mg2+ [8, 9]. Although the photon generationschemes are various
indifferentopticalimagingmodality,theemissionspectrumintheopticalengineeringfieldis
generallyintheso-callednear-infraredlightwindowofthebiologicaltissue([650−900]nm),
which can travel several centimeters in tissue due to the low photon absorption in the above
spectralwindow[7,10,11].
The propagationofthe emission photonsin tissue can beaccuratelyrepresentedby thera-
diative transfer equation (RTE) approximated from Maxwell’s equations, but it is extremely
computationallyexpensiveforitsintegro-differentialnature[12,13].Therefore,thecommonly
usedmathematicalmodelinopticalimagingfield isthediffusionequationderivedfromRTE
in viewof highlyscatteringpropertyof thebiologicaltissue [12, 14, 15].Furthermore,many
algorithmshavebeenpresentedto simulate lighttransportin the turbidtissue andpredictthe
diffuselightfluxonthesurfaceofthesmallanimalbasedondiffusionmodel[16,17].Thecor-
respondingsimulationresultscanbeemployednotonlytoverifythecorrectnessofthephysical
model,butalsotogenerateasensitivitymatrixwhichrelatesthesurfacemeasurementstothe
internalopticalpropertiesandwillbeemployedintheinverselightsourcereconstruction[12].
Itiswellknownthatanalytical,statisticalandnumericaltechniquesarethreekindsofmeth-
odstosolvetheaforementioneddiffusionapproximationmodel[12,18,19,20].Amongthese
methods,numericaltechniquesarestudiedwidelybecauseofitshighefficiencyandgoodap-
plicability, and finite element method (FEM) is one of the most typical and successful algo-
rithms. For example, Cong etal. [17, 21] employed FEM to obtain the photon flux density
on the boundary of the homogeneous and heterogeneous phantoms. Lv etal. [22] used the
adaptiveFEMtocomputethephotonenergydistributiononthephantomsurface.Inaddition,
an improvedFEM was presented to handle light propagationin nonscatteringregionswithin
diffusing domains by Arridge etal. [23]. However, finite element mesh generation and data
pre-processingaredifficultandtime-consuming,especiallyforthree-dimensionalirregularob-
jectswithcomplexinternalstructureliketheheterogeneousbiologicaltissue.Forinstance,the
humanheadmodelwasdiscretizedinto155915elementsintheliterature[24].Whatismore,
#102140 - $15.00 USD Received 29 Sep 2008; revised 3 Nov 2008; accepted 11 Nov 2008; published 24 Nov 2008
(C) 2008 OSA 8 December 2008 / Vol. 16, No. 25 / OPTICS EXPRESS 20319
Shephardetal.[25]reportedameshwithtwomilliontetrahedralelements!Althoughremark-
able progresshasbeen made in generatingthe three-dimensionalstructuredmeshesfor FEM
analysis of solids and structures, it is generally recognized that the development of fast and
robustmeshingtechniquesof three-dimensionalobjectswith complexgeometricalshape and
internalstructureisstillachallengeinpractice.Furthermore,forthelinearFEManalysis,mesh
generationanddatapreparationbeforecalculationoftenneedmuchmoretimethantheassem-
bly andsolutionof the FEM equations.Therefore,meshlessmethodsare exploredasa novel
numericalanalysisapproachwhichcanavoidorgreatlysimplifymeshingtask,andtheyhave
been successfully applied to solve problemsof solid mechanics,heat transfer,fluid flow, etc.
[26,27,28].
In thiscontribution,two Galerkin-basedmeshlessmethods(GBMM) are proposedand de-
velopedforsimulatingthephotonpropagationprocessinthebiologicaltissue.Comparedwith
FEM, GBMM uses onlya set of discretized pointsand doesnotrequireany nodeconnectiv-
ity or element information, which helps not only to avoid the burdensome meshing but also
todescribecomplexinhomogeneousdomainsmoreaccurately.InGBMMalgorithms,moving
leastsquares(MLS)approximationplaysanimportantrole,andtwokindsofGBMMmethods
areprovidedinthispaperaccordingtowhethertheMLSshapefunctionsaremodified.Then,a
linearmatrixformlinkingthesourcedistributionandthephotonfluxdensitycanbeestablished
basedonGalerkinapproachandGausstheorywiththediffusionequationandRobinboundary
condition. Moreover,our two proposed algorithms should incorporate apriori knowledge of
thetissueopticalparameters,whichcanbeassignedonthebasisofavailabledatainliterature
[29] or invivo diffuse optical tomography (DOT) measurements. The paper is organized as
follows.ThenextsectionpresentsthedetailsoftheproposedGBMM, diffusionmodelbased
methodsforphotonpropagationinthebiologicaltissue. Inthethirdsection,theperformance
ofthetwomethodsisvalidatedusingnumericalandphysicalphantomsandcomparedwiththe
simulationdatabasedonFEMorMCandthemeasuredphotonfluxdensitybyacooledCCD
camera.Inthelastsection,relevantissuesarediscussedandconclusionsareprovided.
2. Methodology
2.1. Diffusionapproximationandboundarycondition
Under the assumption that light scattering dominatesover absorption, the propagationof the
emittingphotonsinthebiologicaltissuecanberepresentedbysteady-statediffusionequation
whenacontinuouswaveexternalexcitationlightsourceisusedinfluorescenceimagingexper-
imentorabioluminescentlightsourceisemployedinbioluminescencetechnology[7,14]:
−(cid:209) · D(x)F(cid:209) (x) +m (x)F (x)=S(x) (x∈W ) (1)
a
(cid:0) (cid:1)
where W is the region of interest; F (x) represents the photon flux density at location x
[Watts/mm2]; S(x) denotes the internal source density [Watts/mm3]; m (x) is the absorption
a
coefficient [mm−1]; D(x)=(3(m (x)+(1−g)m (x)))−1 is the optical diffusion coefficient,
a s
m (x)thescatteringcoefficient[mm−1]andgtheanisotropyparameter.
s
Inordertoeliminatethediffusionapproximationerrornearthesurfacewherelightdoesnot
propagate diffusively,an appropriateboundarycondition should be specified [12, 30]. When
the optical imaging experimentis carried out in a totally dark environment,Robin boundary
conditioncanbeemployed[21,30]:
F (x)+2A(x;n,n′)D(x) v(x)·F(cid:209) (x) =0 (x∈¶ W ) (2)
(cid:0) (cid:1)
where¶ W isthecorrespondingboundaryofW ;v(x)referstotheunitnormalvectoroutwardto
the boundary¶ W ; A(x;n,n′) isa functionto incorporatethemismatch betweenthe refractive
#102140 - $15.00 USD Received 29 Sep 2008; revised 3 Nov 2008; accepted 11 Nov 2008; published 24 Nov 2008
(C) 2008 OSA 8 December 2008 / Vol. 16, No. 25 / OPTICS EXPRESS 20320
indicesnwithinW andn′ inthesurroundingmedium.Furthermore,A(x;n,n′)canbeapprox-
imatelyexpressedasthefollowingformulawhentheimagingexperimentisperformedinair,
forwhichn′≈1.0:
A(x;n,n′)≈ 1+R(x) / 1−R(x) (3)
whereR(x)isaparameterrelativetothei(cid:0)nternalre(cid:1)fle(cid:0)ctionat¶ (cid:1)W .Accordingtothereference
[31], R(x) can be approximatedby R(x)≈−1.4399n−2+0.7099n−1+0.6681+0.0636n.In
ourstudy,themeasuredoutgoingfluxdensityQ(x)on¶ W is:
Q(x)=−D(x) v·F(cid:209) (x) =F (x)/ 2A(x;n,n′) (x∈¶ W ) (4)
(cid:0) (cid:1) (cid:0) (cid:1)
2.2. Galerkin-basedmeshlessmethods
2.2.1. Movingleastsquaresapproximation
The MLS approximation is the basis of the proposed GBMM algorithms. According to the
literatures[32]and[33],thephotonfluxdensityF (x)atnodexcanbeapproximatedbyF h(x)
intheregionofinterestW :
m
F (x)≈F h(x)= (cid:229) p (x)a (x)=pT(x)a(x) (5)
j j
j=1
where p (x)isthemonomialbasisfunctionofthespatialcoordinates,andmisthenumberof
j
thebasisfunctions;a (x)isthenon-constantcoefficientwhichcanbedeterminedbyminimiz-
j
ingthefollowingweighteddiscreteL normJ(x):
2
J(x)=(cid:229)Nn w(x−x)[pT(x)a(x)−F ]2 (6)
i i i
i=1
whereN isthenumberofthenodesintheregionofinterestW ;w(x−x)denotestheweight
n i
functionrelatedtothenodex,andx isapointinthesupportdomainofxforwhichw(x−x)6=
i i
0;F istheapproximationtothevalueF (x)atthenodex,whichiscalledasthegeneralized
i i
photonfluxdensity.
AftertheminimizationofJ(x),thefollowinglinearequationcanbeobtained:
A(x)a(x)=B(x)F (7)
where
A(x)=(cid:229)Nn w(x−x)p(x)pT(x) (8)
i i i
i=1
B(x)=[w(x−x )p(x ),w(x−x )p(x ),···,w(x−x )p(x )] (9)
1 1 2 2 Nn Nn
F =(F ,F ,···,F )T (10)
1 2 Nn
Solving a(x) from Eq. (7) and inserting it into Eq. (5), we can obtain the following MLS
approximationform:
F h(x)=(cid:229)Nn N(x)F (11)
i i
i=1
wheretheshapefunctionN(x)isdefinedby
i
N(x)=pT(x)A−1(x)B(x) (12)
i i
whereB(x)standsfortheithcolumnofthematrixB(x).
i
#102140 - $15.00 USD Received 29 Sep 2008; revised 3 Nov 2008; accepted 11 Nov 2008; published 24 Nov 2008
(C) 2008 OSA 8 December 2008 / Vol. 16, No. 25 / OPTICS EXPRESS 20321
ThepartialderivativesofN(x)aregivenby
i
N (x)=pTA−1B +pT[A−1(B −A A−1B)] (13)
i,s ,s i i,s ,s i
wheresrepresentsthespace variablex, y orz, andthecommaindicatesthepartialderivative
withregardtothespatialcoordinatethatfollows.
FromEq.(12),we cansee thattheperformanceoftheMLSapproximationisgovernedby
thebasisfunctionandtheweightfunction.Inthispaper,thequadraticbasisfunction
pT(x)=[1,x,y,z,x2,xy,y2,yz,z2,zx], m=10 (14)
andthefollowingquarticsplineweightfunctionareusedinthethree-dimensionalcaseaccord-
ingtothereferences[33]and[34].
1−6r2+8r3−3r4 0≤r≤1
w(r)= (15)
(cid:26) 0 r>1
where r = kx−xk/d is the ratio of the Euclidean distance between the evaluation point
i
x and the sampling node x to the radius of the support domain d. Removing all the vari-
i
ablescorrelatedwithz-componentinthree-dimensionalGBMMprocedure,wecanobtaintwo-
dimensionalGBMMprogrameasily.
2.2.2. ModifiedMLSshapefunction
Substitutingx=x backintoEq.(11),wehave
k
F h(x )=(cid:229)Nn N(x )F =NTF (16)
k i k i k
i=1
where F designates the generalized photon flux density at all nodes in W , and N =
k
[N (x ),N (x ),···,N (x )]T.Andthen,Eq.(16)canbefurtherwrittenasthefollowingmatrix
1 k 2 k Nn k
equation:
F =LF (17)
whereF denotesthenodalphotonfluxdensbity,andL isreferredtoasthetransformationmatrix
[35,36].Theycanbeexpressedas:
b
F =[F h(x ),F h(x ),···,F h(x )]T (18)
1 2 Nn
b N (x ) N (x ) ··· N (x )
1 1 2 1 Nn 1
L = N1(...x2) N2(...x2) ·.·..· NNn...(x2)  (19)
 
 N1(xNn) N2(xNn) ··· NNn(xNn) 
Therefore,thegeneralizedphotonfluxdensitycanbeobtainedfromEq.(17):
F =L −1F (20)
and b
F = (cid:229)Nn N(x)−1F (21)
i l i l
l=1
b
whereL −1istheinversematrixofL .IncorporatingEq.(21)withEq.(11),wehave
F h(x)=(cid:229)Nn N(x)(cid:229)Nn N(x)−1F = (cid:229)Nn M(x)F (22)
i l i l l l
i=1 l=1 l=1
b b
#102140 - $15.00 USD Received 29 Sep 2008; revised 3 Nov 2008; accepted 11 Nov 2008; published 24 Nov 2008
(C) 2008 OSA 8 December 2008 / Vol. 16, No. 25 / OPTICS EXPRESS 20322
whereM(x)=(cid:229) Nn N(x)N(x)−1 iscalledmodifiedMLSshapefunction,anditsatisfiesthe
l i=1 i l i
Kroneckerdeltafunctionproperty:
M(x )=(cid:229)Nn N(x )N(x)−1=d (23)
l k i k l i lk
i=1
2.2.3. Numericalimplementation
According to whether modifying the MLS shape function, two meshless methods based on
Galerkinapproacharedevelopedinthisarticle.
UsingGalerkinmethodandGausstheory,Eqs.(1)and(2)canbetransformedtothefollow-
ingweakform[30,37]:
D(x) F(cid:209) (x) · Y(cid:209) (x) +m (x)F (x)Y (x) dx
a
ZW (cid:16) (cid:0) (cid:1) (cid:0) (cid:1) (cid:17)
1
+ F (x)Y (x)dx= S(x)Y (x)dx (24)
Z¶ W 2A(x;n,n′) ZW
whereY (x)isatestfunctionfromSobolevspace.TheaforementionedEqs.(11)and(22)can
berewritteninthefollowingunifiedform:
F h(x)= (cid:229)Nn ¡ (x)G (25)
l l
l=1
where ¡ (x) represents N(x) in Eq. (11) or M(x) in Eq. (22); G denotes the generalized or
l l l
nodalphotonfluxdensity.
Substituting Eq. (25) into Eq. (24) and using ¡ (x) as the test function, and then we can
l
obtainthematrixequationasfollows:
(K+C+F)G =GG =S (26)
wherethecomponentsofthematricesK,C,FandthevectorSaregivenby
Kkl = W D(x) ¡(cid:209) k(x) · ¡(cid:209) l(x) dW
 Ckl =RW m a(x)(cid:0)¡ k(x)¡ l((cid:1)x)(cid:0)dW (cid:1) (27)
 Fkl = R¶ W ¡ k(x)¡ l(x)/ 2A(x;n,n′) d¶ W
Sk= RW S(x)¡ k(x)dW (cid:0) (cid:1)
SincethematrixGissymmetriRcandpositivedefinite[34],G canbeuniquelydeterminedfrom
G =G−1S (28)
When the MLS shape functions are employed, G that solved from Eq. (28) are only the
generalizedphotonfluxdensityF becausetheMLSapproximationdoesnotpassthroughthe
nodal function values according to the preceding subsection 2.2.1. In order to get the actual
photonfluxdensity,thatisthenodalphotonfluxdensityF ,atanypointinthegivendomainW ,
weneedtousetheMLSapproximationagainwiththeformula:
b
F =[N (x ),N (x ),···,N (x )][F ,F ,···,F ]T (29)
k 1 k 2 k Nn k 1 2 n
whereF isthenodalbphotonfluxdensityatpointx ;N (x ),N (x ),···,N (x )aretheMLS
k k 1 k 2 k Nn k
shapefunctionswithoutmodification.However,thenodalphotonfluxdensitycanbeobtained
b
directlythroughsolvingEq.(28)whenweusethemodifiedMLSapproximation.
Tosumup,theflowchartoftheabovetwoGBMMalgorithmsisshowninFig.1.
#102140 - $15.00 USD Received 29 Sep 2008; revised 3 Nov 2008; accepted 11 Nov 2008; published 24 Nov 2008
(C) 2008 OSA 8 December 2008 / Vol. 16, No. 25 / OPTICS EXPRESS 20323
Fig.1.(a)TheflowchartoftheproposedGBMMalgorithmwithoutmodification;(b)The
flowchartofthepresentedGBMMalgorithmwithmodification.
3. Experimentsandresults
To evaluate the developed GBMM algorithms in this paper, numerical and physical phan-
tom experiments were performed respectively. Furthermore, the computational results based
onGBMMwerecomparedwiththesimulationdatabyFEMorMCandthemeasuredphoton
fluxdensityusingaCCDcamera.Forthesakeofconvenience,GBMM1representstheGBMM
algorithmwithoutmodification,andGBMM2indicatestheothermethodinthissection.
#102140 - $15.00 USD Received 29 Sep 2008; revised 3 Nov 2008; accepted 11 Nov 2008; published 24 Nov 2008
(C) 2008 OSA 8 December 2008 / Vol. 16, No. 25 / OPTICS EXPRESS 20324
3.1. Numericalphantomexperiments
3.1.1. Homogeneousnumericalphantomexperiments
Firstly, a homogeneous tissue-like phantom was used in this experiment, and a light source
with a total power of 0.125nano−Watts was placed in the phantom with its center at
(6.5358,7.1537,15.0729).InGBMMstudy,2178randompointsweredistributedintheabove
phantom, as showed in Fig. 2(a). In addition, apriori optical parameters were specified as
m =0.035,m =6.0,g=0.9andn=1.37torepresentthediffusivebiologicaltissue,which
a s
couldbeobtainedfromtheliterature[29].Finally,thelightexitancemaponthephantomsur-
facewassolvedusingGBMM1andGBMM2proceduresrespectively,aspresentedinFig.2(b)
and2(c).Furthermore,thecomputationalresultofGBMM1algorithmiscompletelyidentical
tothatofGBMM2.
Fig.2.Homogeneousnumericalphantom.(a)Ahomogeneoustissue-likephantomwitha
seriesofnodesandalightsource;(b)and(c)ThesurfacelightpowersimulatedbyGBMM1
andGBMM2.
In order to demonstratethe accuracy of GBMM, FEM was also employedto calculate the
surfacephotonfluxdensity.InFEMframework,thephotonfluxdensitiesat24654discretized
nodes were simulated, and then we used the interpolation method to determine the photon
fluenceatthearrangedpointsinFig.2(a).Figure3(a)and3(b)givethevolumetricmeshused
in FEM and the correspondingsimulation result respectively.The computationalphotonflux
densitybasedonGBMM wasingoodagreementwiththenumericalresultbyFEM,withthe
averagerelativeerrorbeingabout0.9%,asshowedinFig.3(c).
Asweallknow,MCmethod,regardedasthe“goldstandard”,isrigorous,flexibleandpow-
erfultostudyphotontransportphenomenainthebiologicaltissue,withwhichothernumerical
techniques are often compared [12, 18, 19, 20, 30]. Therefore, MC approach was also used
to verify the performance of the above two GBMM algorithms in this paper. In the simula-
tion experiment,a cubic lightsource of 2.0mm side lengthand 1.0nano−Watts/mm3power
densitywasplacedin(11.0,11.0,11.0),andfourcubephantomscenteredat(10.0,10.0,10.0)
withdifferentsidelengthfrom5mmto15mmwereemployedtoobtaintheirrespectivesurface
photonfluxes.Table 1 lists the correspondingcomparativedata betweenGBMM1, GBMM2,
FEMandMCmethod.RE1,RE2andRE3inTable1aretherelativeerrors(RE)betweenthe
computationalresultsbyGBMM1,GBMM2,FEMandthesimulationdatausingMCmethod
respectively.ComparingtheGBMM,FEMsolutionwiththeMCsimulationresults,ithasfound
that theyhave the same tendency,while the GBMM resultsare in better agreementwith MC
simulationdatathantheFEMsolution.Furthermore,two GBMMmethodsinthispaperhave
thesamecalculationprecision.
#102140 - $15.00 USD Received 29 Sep 2008; revised 3 Nov 2008; accepted 11 Nov 2008; published 24 Nov 2008
(C) 2008 OSA 8 December 2008 / Vol. 16, No. 25 / OPTICS EXPRESS 20325
Fig.3.FEMsimulationforhomogeneousphantom.(a)ThevolumetricmeshusedinFEM
simulation; (b) The photon flux density on the phantom surface calculated by FEM. (c)
ComparisonofthecomputationalresultsbyGBMMandFEM.
Table1.Photonflux(nano−Watts)simulationforhomogeneousphantom.
Sidelength MC GBMM1 RE1 GBMM2 RE2 FEM RE3
5mm 45.26 45.86 1.32% 45.86 1.32% 45.89 1.39%
8mm 42.81 43.26 1.05% 43.26 1.05% 43.29 1.12%
12mm 39.13 39.35 0.56% 39.35 0.56% 39.39 0.66%
15mm 36.21 36.25 0.11% 36.25 0.11% 36.30 0.25%
3.1.2. Heterogeneousnumericalphantomexperiments
A cubical heterogeneous phantom with 8000mm3 volume was utilized to further test
the performance of the GBMM algorithms, and two cube light sources of 2.5mm side
length and 1.0nano−Watts/mm3 power density were located at (6.25,6.25,11.25) and
(13.75,13.75,11.25)intheabovephantomrespectively.Theopticalparameterswereassigned
toeachofthetwocomponents,aslistedinTable2.InFig.4and5,weillustratetheabilityof
our developedGBMM algorithmsto simulate photon transportin tissue in the multiple light
sourcescase.We cansee thatthesimulationdatabyGBMM1andGBMM2isidentical.Fur-
thermore, the average relative error of the results obtained using GBMM and FEM is about
2.40%. In order to better demonstrate the performance of GBMM1 and GBMM2, we com-
paredGBMMphotonfluxdensitywithFEMcalculationalresultalongthedetectionsquareon
thephantomsurfaceatheights5mm,10mmand15mmfromthebottomofthegeometricmodel,
asshowedinFig.5(c)-5(e)respectively.
Table2.Opticalparametersoftheheterogeneousphantom.
m (mm−1) m (mm−1) g n
a s
Region1 0.035 6.0 0.9 1.37
Region2 0.01 4.0 0.9 1.37
Similarly,fourheterogeneouscubicphantomswithdifferentsidelengthfrom11mmto20mm
weresetupfornumericalsimulationusingMCmethod.Acubelightsourcewithatotalvolume
of8.0mm3andapowerdensityof1.0nano−Watts/mm3wasembeddedintothephantomswith
itscenterat(11.0,11.0,11.0).Finally,thephotonfluxesontheboundaryofthephantomswere
solved by GBMM1, GBMM2, FEM andMC respectively,see Table 3. There are no any dif-
ferencebetweenthecomputationalresultsofGBMM1andthoseofGBMM2,andtheGBMM
#102140 - $15.00 USD Received 29 Sep 2008; revised 3 Nov 2008; accepted 11 Nov 2008; published 24 Nov 2008
(C) 2008 OSA 8 December 2008 / Vol. 16, No. 25 / OPTICS EXPRESS 20326