Elektrotehniški vestnik 84(1-2): 61-65, 2017 Izvirni znanstveni članek Vpliv zmanjšanja simulacijskega volumna na napako in hitrost izračunov reflektance z metodo Monte Carlo Miran BUrmen1 1 Univerza v Ljubljani, Fakulteta za elektrotehniko, Tržaška 25, 1000 Ljubljana, Slovenija E-pošta: miran.buermen@fe.uni-lj.si Povzetek. Širjenje nepolarizirane svetlobe po bioloških tkivih pogosto formuliramo v okviru enačbe sevalnega prenosa energije. Ta za večino uporabnih geometrij ni analitično rešljiva, zato jo najpogosteje rešujemo numerično v okviru simulacij Monte Carlo. Tovrstno metodologijo je mogoče uspešno uporabiti v številnih eksperimentalnih postavitvah, ki vključujejo tudi prostorsko razločene meritve povratno sipane svetlobe (reflektance) z optično vlaknenimi sondami. V okviru te študije smo podrobneje proučili vpliv omejenega simulacijskega volumna na čas in natančnost izračunov prostorsko razločene reflektanče zajete z linearno razporeditvijo optičnih vlaken v sondi. Rezultati kazejo, da je napaka pri izračunih reflektance s simulacijami Monte Carlo pri uporabi omejenega simulacijskega volumna največja za vzorce z nizkim absorpcijskim in reduciranim sipalnim koeficientom. Nasprotno, rezultati kazejo, da je napaka le malo odvisna od faktorja anizotropije sipalne fazne funkcije. Z zmanjšanjem simulacijskega volumna je mogoče izračune reflektance pohitriti tudi do dvakrat, a le pri vzorcih z nizkim absorpcijskim (« 1cm-1) in do zmernim reduciranim sipalnim koeficientom (« 30 cm-1). Poleg številnih bioloških tkiv se s tovrstnimi optičnimi lastnostmi pogosto srečujemo pri optičnih fantomih za vrednotenje in umerjanje reflektance izračunane s simulacijami Monte Carlo. Ključne besede: reflektančna spektroskopija, absorpcija, sipanje, Monte Carlo simulacije, napaka, hitrost izračunov, optične sonde The impact of a reduced simulation volume on the accuracy and computational time of the Monte Carlo reflectance simulations Light propagation in biological tissues formulated within the framework of the radiative transport equation is frequently simulated by stochastic Monte Carlo simulations. Such methodology is used in numerous experimental settings including spatially resolved reflectance measurements conducted by optical fiber probes. In this paper, we analyse the impact of the simulation volume on the computational time and relative error of the Monte Carlo simulated reflectance for an optical fiber probe with a linear fiber layout. The results show, that the largest relative reflectance error for a reduced simulation volume is observed for samples with low absorption and reduced scattering coefficients and for larger source detector separations. A subsequent analysis revealed that the reflectance error does not significantly depend on the anisotropy factor of the scattering phase function. A twofold reduction in the computational time of the Monte Carlo simulations can be attained without significantly affecting the accuracy of the computed reflectance for the given experimental setting. However, the computational time is improved only for the samples that exhibit low absorption (« 1 cm-1) and up to medium reduced scattering (« 30 cm-1) coefficients. Such optical properties can be found in many biological tissues, and in particular in purely scattering optical phantoms that are frequently used for validation and calibration of the Monte Carlo reflectance simulations. Prejet 12. december, 2016 Odobren 20. januar, 2017 1 Uvod Povratno sipana svetloba (reflektanca) iz biološkega tkiva je bogata z informacijo o njegovi kemični sestavi in strukturi. Pri tem lahko absorpcijo svetlobe poveZemo s sestavo glavnih tkivnih kromofoijev (kri, voda, melanin, maščobe, rumeni pigmenti), kompleksni proces sipanja svetlobe pa z lokalnimi spremembami lomnega kolicšnika, ki so posledica morfologije in mikrostruk-ture tkiva [10], [13], [14], [11]. Prostorsko in valovno razlocšene meritve povratno sipane svetlobe (reflektance) lahko ucšinkovito izvedemo z opticšnimi sondami z vecš vlakni [5], [12], [8], [2]. Širjenje nepolarizirane svetlobe v sipajocih medijih zadovoljivo dobro opišemo v okviru enacbe sevalnega prenosa energije (RTE). Zajeto reflektanco v tem primeru opišemo z absorpcijskim koeficientom pa in sipalnim koeficientom ps, ki dolocata atenuacijo kolimiranega snopa svetlobe, ter sipalno fazno (verjetnostno) funkcijo p(s, s'), ki doloca verjetnost sipanja fotona iz vpadne smeri s v smer s'. Pri bioloških tkivih pogosto predpostavimo izotropnost sipalne fazne funkcije, ki se tako poenostavi v p(cos 6), kjer 6 pomeni kot med s in s'. V obstojeci literaturi se sipalne lastnosti tkiv pogosto poenostavljeno podajajo z reduciranim si-palnim koeficientom p' = (1~9)^s, ki povezuje sipalni koeficient in prvi Legnderov moment g (povprecni kot sipanja) sipalne fazne funkcije [4]. Zal RTE v splošnem ni analiticno rešljiva [7], [6] zato se pri njenem 62 MIRAN BÜRMEN reševanju najpogosteje poslužujemo stohasticnih simulacij Monte Carlo (MC) [15]. Pri tem simulacije dodatno pohitrimo ž uporabo strojne opreme, najpogosteje ž vecjedrnimi grafičnimi pospeševalniki [1]. Ko izčrpamo razpoložljive možnosti za strojno pohitritev simulacij, nam ostanejo še tiste, ki so vežane na poenostavitve eksperimentalne geometrije. Pri tem moramo pažiti, da s tovrstnimi poenostavitvami ne vnasšamo dodatne ali prevelike napake v numericne ižracune reflektance. V okviru te študije smo podrobneje proucili možnosti ža dodatne pohitritev Monte Carlo simulacij prostorsko ražlocšene reflektance ž uporabo omejenega simu-lacijskega volumna. Pri tovrstni simulacijah prenehamo slediti vsem fotonom, ki žapustijo ižbrano okolico iž-vornega opticšnega vlakna. Pristop je žanimiv predvsem takrat, ko je vrednost absorpcijskega koeficienta vžorca nižka, in se posledicšno fotoni lahko bistveno oddaljijo od ižvornega opticšnega vlakna. 2 Metodologija Za namene študije smo uporabili polneskoncno geometrijo biološkega tkiva in opticno vlakneno sondo s sšestimi linearno ražporejenimi opticšnimi vlakni premera 200 pm in raždaljo med središci vlaken 220 pm. Tako smo v vsaki simulaciji ižracšunali reflektanco pri petih raždaljah med ižvornim in sprejemnim vlaknom. Pri tem smo skrajno levo vlakno uporabili kot ižvorno. Vrednost numericne odprtine (ang. numerical aperture) opticnih vlaken smo postavili na NA = 0,22, lomni kolicnik pa na nfib = 1,45. Mejo med opticno sondo in tkivom smo pri ižracšunih obravnavali lateralno uniformno. V simula- Slika 1: Prikaž eksperimentalne postavitve s precnim prerežom opticne sondo in pripadajoco ražporeditev oddajnega (rumeno) in sprejemnih (modra) opticnih vlaken na naležni površini cijah Monte Carlo smo žacetno požicijo fotonov vžorcili enakomerno s precšnega preseka ižvornega opticšnega vlakna, žacšetno smer pa enakomerno iž prostorskega kota, ki ga omejujeta numericšna odprtina vlakna in lomni kolicšnik tkiva. Postopek nakljucšnega vžorcšenja žacetnih požicij in smeri širjenja fotonov obicajno udejanjimo s pomocjo generatorja enakomerno ražporejenih psevdonakljucnih števil iž intervala [0,1]. V ta namen je vsakokrat treba uporabiti nova štiri psevdonakljucna števila, in sicer ^ do £4. Zacetne koordinate (x,y,z) in smer širjenja fotona (sx,sy,sz) nato dolocimo po naslednjih enačbah: x = xo + r^v^l cos(2n&), (1) y = yo + r sin(2n^2), (2) z = 0, (3) dNA = 1 - ^4(1 - cosarcsin(NA/nt)), (4) sx = cos(2n^3) sin(^NA), (5) sy = sin(2n^3) sin(^NA), (6) sz = cos(6 n a) , (7) kjer rv pomeni radij x0 in y0 pa koordinati središca ižvornega opticnega vlakna, nt pa lomni kolicnik opažo-vanega vžorca. Podrobnejšo geometrijo eksperimentalne postavitve in opticne sonde prikažuje slika 1. Opticne lastnosti tkiva, ki smo jih obravnavali v okviru te študije, smo ižbrali tako, da pokrivajo celoten ražpon vrednosti absorpcijskega in reduciranega sipalnega koeficienta, ki jih navaja obstojeca literatura [4]. Pri tem smo vrednost absorpcijskega koeficienta ^a spreminjali v ražponu od 0cm-1 do 25 cm-1, vrednosti reduciranega sipalnega koeficienta pa v ražponu od 5 cm-1 do 60 cm-1, in sicer oba s 25 enakomerno ražporejenimi koraki. Lomni kolicnik biološkega tkiva smo postavili na nm = 1, 33. V simulacijah MC smo uporabili Henyey-Greenstein (HG) [3] model sipalne fažne funkcije, ki ga doloca naslednja enacšba: p(cos 0) = ^ 1 - g2--, (8) 2 1 + g2 — 2g cos 0 kjer 0 pomeni kot med smerjo širjenja fotona pred sipanjem in po njem, parameter g pa faktor anižotropije sipalne fažne funkcije (povprecna vrednost cos(0)). Pri simulacijah ž omejenim volumnom tkiva smo prenehali slediti vsem fotonom, ki so se od ižvornega vlakna oddaljili ža vec kot rt. Referencne ižracune reflektance pri vrednosti absorpcijskega koeficienta 0cm-1 smo iž-vedli pri veliki vrednosti radija simulacijskega volumna rt = 30 mm. Vsi ižracuni reflektance v tej študiji so bili ižvedeni s pomocšjo simulacij Monte Carlo podprtih ž OpenClTM vmesnikom. Pri tem smo v eni simulaciji sledili 108 fotonom. Za dolocitev reflektance, ižmeijene skoži posamežna opticšna vlakna, smo ižkoristili radialno simetrijo reflektance in s tem dodatno ižboljšali kakovost signala. 3 Rezultati V prvem eksperimentu smo poskušali identificirati tista podrocšja opticšnih lastnosti, pri katerih nastane kot posledica omejenega simulacijskega volumna najvecja VPLIVA SIMULACIJSKEGA VOLUMNA NA NAPAKO IN HITROST IZRAČUNOV 63 napaka pri izračunih reflektance. V ta namen smo uporabili relativno majhen radij simulacijskega volumna rt = 1.5 mm, ki komajda se zaobjame vsa optična vlakna. Ta lahko v dani geometriji popolnoma zaob-jamemo s hemisfero radija 1.2 mm. Slika 2 prikazuje relativno napako izračunane reflektance pri vrednosti radija simulacijskega volumna rt = 1.5 mm, in sicer za najmanjšo (220 |im) in največjo (1100 |im) razdaljo med izvornim in sprejemnim opticnim vlaknom. Vrednost faktorja anizotropije g sipalne fazne funkcije HG smo postavili na g =0.85 Vidimo, da najvecja relativna napaka izracunane reflektance nastopi pri nizkih vrednostih absorpcijskega in reduciranega sipalnega koeficienta, in narašca z razdaljo med izvornim in sprejemnim vlaknom. Posledicno smo v drugem eksperimentu podrobneje proucili vpliv vrednost radija simulacijskega volumna rt na relativno napako izracunane reflektance. V ta namen smo vrednost absorpcijskega koeficienta postavili na 0.01 cm-1, vrednost reducira-nega sipalnega koeficienta pa na 5 cm-1, kar je najnizja vrednost, ki jo obicajno še srecamo med biološkimi tkivi v vidnem in zacetnem bliznjem infrardecem delu spektra [4]. Vrednost faktorja anizotropije sipalne fazne funkcije smo postavili na 0,85. Slika 3 prikazuje, kako se relativna napaka izracšunane reflektance spreminja z radijem simulacijskega volumna rt in razdaljo med izvornim in sprejemnim opticnim vlaknom. Vidimo, da je relativna napaka izracšunane reflektance pri vseh sprejemnih opticnih vlaknih manjša od 2%, ko je radij simulacijskega volumna rt priblizno vecji ali enak 8 mm. V naslednjem koraku smo proucšili, kako se relativna napaka reflektance spreminja v odvisnosti od faktorja anizotropije sipalne fazne funkcije. Slika 4 prikazuje relativno napako reflektance v odvisnosti od razdalje med izvornim in sprejemnim opticšnim vlaknom za vrednosti faktorja anizotropije med 0,4 in 0,95. Pri izracunih smo uporabili radij simulacijskega volumna 8 mm, vrednost absorpcijskega koeficienta 0.01 cm-1 ter vrednost reduciranega sipalnega koeficienta 5 cm-1. Za konec smo ocenili še faktorje pohitritve simulacij Monte Carlo, ki jih dosezemo pri simulacijskem volumnu rt = 8 mm, ko relativna napaka reflektance za biolosško zanimiv nabor opticšnih lastnosti ne presezše 2 %. Slika 5 prikazuje izracšunane faktorje pohitritve pri vrednosti faktorja anizotropije 0,85. 4 Razprava Rezultati eksperimentov kažejo, da lahko z zmanjšanjem simulacijskega volumna bistveno, tudi do dvakratno, pohitrimo MC simulacije, ne da bi pri tem bistveno vplivali na natancnost izracunov reflektance. Prakticna uporaba tovrstnih pohitritev v izbrani geometriji opticšnih vlaken je omejena na vzorce z nizkimi absorpcijskimi koeficienti (do « 1cm-1) in vrednostmi reduciranega sipalnega koeficienta do priblizno 30 cm-1. Pri tem je relativna napaka izracunane reflektance le malo odvisna od faktorja anizotropije sipalne fazne funkcije (slika 4), hkrati pa mocšno odvisna od razdalje med izvornim in sprejemnim opticšnim vlaknom (slika 3). V dani eksperimentalni geometriji lahko pri relativno nizki vrednosti absorpcijskega in reduciranega sipalnega koeficienta dosezemo okoli dvakratno pohitritev izracunov, ne da bi pri tem relativna napaka reflektance presegla 2 %. Pri narašcajocih vrednostih reduciranega sipalnega koeficienta, pa se pohitritev zmanjša in pri 30 cm-1 znaša priblizšno 30 %. Z opticšnimi lastnostmi, ki omogocšajo uporabo tovrstnih pohitritev, se pogosto srecujemo pri opticšnih fantomih, ki so namenjeni za vrednotenje in umerjanje simulacij Monte Carlo [9]. Slednje pogosto pripravimo iz vodnih suspenzij polistirenskih mikros-fericšnih delcev, ki brez dodatkov izkazujejo zelo nizek absorpcijski koeficient. Pri taksšnih opticšnih lastnostih je skoraj nujno uporabiti koncšni simulacijski volumen. Pred prakticšno uporabo omejenega simulacijskega volumna je za izbrano eksperimentalno geometrijo treba podrobno proucšiti vpliv na napako izracšunane reflek-tance. Pri tej se pogosto zadovoljimo z najvec 2 odstotno relativno napako. Glede na prikazane rezultate vidimo, da je napaka pri izracunih reflektance najmanjša za visoke vrednosti absorpcijskega koeficienta. Rezultati so smiselni, saj z narašcajoco absorpcijo fotoni pospešeno ugašajo in se vecina tistih fotonov, ki bi jih pri pol-neskoncšni obravnavi medija zaznali skozi sprejemna opticšna vlakna, ne uspe dovolj oddaljiti od izvornega opticšnega vlakna, da bi povzrocšili prekinitev simulacije. To je tudi razlog, da pri tovrstnih opticšnih lastnostih vzorca zmanjšanje simulacijskega volumna ne vodi do bistvenih pohitritev simulacij Monte Carlo. Nasprotno se pri nizki vrednosti absorpcijskega in reduciranega sipalnega koeficienta, ko medij postaja prosojen, fotoni uspejo mocšno oddaljiti od izvornega vlakna. Posledicšno lahko z zmanjšanjem simulacijskega volumna bistveno pohitrimo izracšune, a hkrati tudi najhitreje pridelamo nesprejemljivo veliko napako izracšunane reflektance, saj je delez fotonov, ki zapustijo simulacijski volumen in bi jih pri polneskoncšni obravnavi medija zaznali skozi sprejemna opticna vlakna, velik. Narašcanje relativne napake reflektance z razdaljo med izvornim in sprejemnim vlaknom je pricakovano, saj lahko fotoni pri vecji razdalji med izvornim in sprejemnim vlaknom uberejo daljše poti, še zlasti pri nizkih vrednostih absorpcijskega koeficienta. 5 Sklep V sštudiji smo pokazali, da lahko s smiselno izbiro simulacijskega volumna pohitrimo izracune simulacij MC pri nizkih vrednostih absorpcijskega koeficienta, ne da bi pri tem bistveno vplivali na natancšnost izracšunov reflektance. Faktor pohitritve, ki ga na ta nacin lahko dosezemo, je odvisen predvsem od opticnih lastnosti vzorca in geometrije eksperimentalne postavitve. Na splosšno je najbolje, da pred uporabo natancšno ovre- 64 MIRAN BÜRMEN 5 10 E o ¡10 -t-» m a k a 10 o a n a n iv 15 ro el R 20 2 % Razdalja med izvornim in sprejemnim vlaknom: . — 220 pm — 440 pm — 660 pm ' — 880 pm 1100 pm . g = 0,85 HI ~8 mm iua = 0,01 cm-1 = 5 cm-1 5 10 15 Radij simulacijskega volumna rt (mm) 20 C ra 2 3 0 3 ra cp ra c ra c > 0 a: Razdalja med izvornim in sprejemnim vlaknom: — 220 |jm — 440 |jm — 660 |jm — 880 jim — 1100 jim = 5 cm-1 = 0,01 cm' rt = 8 mm 0,4 0,5 0,6 0,7 0,8 Faktor anizotropije g 0,9 Red. sipalni koeficient 10 20 30 40 (cm-1) 50 60 15 & 20 o s _Q < 25 Slika 3: Relativna napaka izracunane reflektance v odvisnosti od radija simulacijskega volumna rt pri vrednosti absorpcijskega koeficienta 0.01 cm-1, reduciranega sipalnega koeficienta 5 cm-1 ter faktorja anizotropije 0,85 Slika 4: Relativna napaka reflektance v odvisnosti od razdalje med izvornim in sprejemnim opticšnim vlaknom za vrednosti faktorja anizotropije med 0,4 in 0,95 pri absorpcijskem koeficientu 0.01 cm-1, reduciranem sipalnem koeficientu 5 cm-1 in radiju simulacijskega volumna rt =8 mm Slika 5: Faktor pohitritve simulacij Monte Carlo kot posledica zmanjšanja simulacijskega volumna pri rt = 8 mm. Vrednost faktorja anizotropije g je bila pri izracunih 0,85 dnotimo vpliv zmanjšanega simulacijskega volumna za izbrano eksperimentalno geometrijo in pricakovani razpon opticnih lastnosti proucevanih vzorcev. Na podlagi rezultatov lahko nato dolocšimo najmanjsši simulacijski volumen, ki nam daje še zadovoljivo tocne izracune reflektance. Zahvala Raziskavo sta omogočila Ministrstvo za visoko šolstvo, znanost in tehnologijo Republike Slovenije v okviru programa P2-0232 in Javna agencija za raziskovalno dejavnost Republike Slovenije v okviru projektov J2-7211, J7-6781, J2-7211 in J2-7118. Literatura [1] Alerstam, E., Svensson, T., Andersson-Engels, S.: Parallel computing with graphics processing units for high-speed Monte Carlo simulation of photon migration. Journal of Biomedical Optics 13(6), 060504 (Dec 2008). [2] Bravo, J. J., Paulsen, K. D., Roberts, D. W., Kanick, S. C.: Subdiffuse optical biomarkers characterize localized microstructure 5 0 5 5 10 0 5 ra 2 k 1 VPLIVA SIMULACIJSKEGA VOLUMNA NA NAPAKO IN HITROST IZRAČUNOV 65 and function of cortex and malignant tumor. Optics Letters 41(4), 781-784 (Feb 2016). [3] Henyey, L. G., Greenstein, J. L.: Diffuse radiation in the Galaxy. The Astrophysical Journal 93, 70-83 (Jan 1941). [4] Jacques, S. L.: Optical properties of biological tissues: a review. Physics in Medicine and Biology 58(11), R37 (Jun 2013). [5] van Leeuwen-van Zaane, F., Gamm, U. A., van Driel, P. B. A. A., Snoeks, T. J. A., de Bruijn, H. S., van der Ploeg-van den Heuvel, A., Mol, I. M., Lowik, C. W. G. M., Sterenborg, H. J. C. M., Amelink, A., Robinson, D. J.: In vivo quantification of the scattering properties of tissue using multi-diameter single fiber reflectance spectroscopy. Biomedical Optics Express 4(5), 696 (May 2013). [6] Liemert, A., Kienle, A.: Light transport in three-dimensional semi-infinite scattering media. Journal of the Optical Society of America A 29(7), 1475 (Jul 2012). [7] Liemert, A., Kienle, A.: Novel analytical solution for the radiance in an anisotropically scattering medium. Applied Optics 54(8), 1963 (Mar 2015). [8] Naglic, P., Pernus, F., Likar, B., Biirmen, M.: Limitations of the commonly used simplified laterally uniform optical fiber probe-tissue interface in Monte Carlo simulations of diffuse reflectance. Biomedical Optics Express 6(10), 3973-3988 (Oct 2015). [9] Pogue, B. W., Patterson, M. S.: Review of tissue simulating phantoms for optical spectroscopy, imaging and dosimetry. Journal of Biomedical Optics 11(4), 041102 (Aug 2006). [10] Rogers, J. D., Radosevich, A. J., Yi, J., Backman, V.: Modeling light scattering in tissue as continuous random media using a versatile refractive index correlation function. IEEE journal of selected topics in quantum electronics : a publication of the IEEE Lasers and Electro-optics Society 20(2), 7000514 (2013). [11] Schmitt, J. M., Kumar, G.: Turbulent nature of refractive-index variations in biological tissue. Optics Letters 21(16), 1310-1312 (Aug 1996). [12] Utzinger, U., Richards-Kortum, R. R.: Fiber optic probes for biomedical optical spectroscopy. Journal of Biomedical Optics 8(1), 121-147 (Jan 2003). [13] Wang, R. K.: Modelling optical properties of soft tissue by fractal distribution of scatterers. Journal of Modern Optics 47(1), 103120 (Jan 2000). [14] Xu, M., Alfano, R. R.: Fractal mechanisms of light scattering in biological tissue and cells. Optics Letters 30(22), 3051-3053 (Nov 2005). [15] Zhu, C., Liu, Q.: Review of Monte Carlo modeling of light transport in tissues. Journal of Biomedical Optics 18(5), 50902 (May 2013). Miran Biirmen je zaposlen na Fakulteti za elektrotehniko Univerze v Ljubljani, kjer se pedagoško in raziskovalno ukvarja s spektroskopijo in hiperspektralnim slikanjem, modeliranjem spektroskopskih in hiperspektralnih slikovnih sistemov, obdelavo in analizo spektrov in hiperspektralnih slik ter modeliranjem širjenja svetlobe po sipajocih medijih.