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