Metoda Monte Carlo pro transport fotonů - Monte Carlo method for photon transport

Modelování šíření fotonů pomocí Metody Monte Carlo je flexibilní, ale přísný přístup k simulaci transportu fotonů. V této metodě jsou místní pravidla přenosu fotonů vyjádřena jako rozdělení pravděpodobnosti, která popisují velikost kroku pohybu fotonů mezi místy interakce foton-tkáň a úhly vychýlení v trajektorii fotonu, když dojde k události rozptylu. To je ekvivalentní analytickému modelování přenosu fotonů metodou radiační přenosová rovnice (RTE), který popisuje pohyb fotonů pomocí diferenciální rovnice. Uzavřená řešení RTE však často nejsou možná; pro některé geometrie je difúzní aproximace lze použít ke zjednodušení RTE, i když to zase přináší mnoho nepřesností, zejména v blízkosti zdrojů a hranic. Naproti tomu simulace Monte Carlo lze libovolně přesně zvýšit zvýšením počtu sledovaných fotonů. Podívejte se například na film, kde je simulace Monte Carlo a tužkový paprsek incident na a napůl nekonečný střední modely jak počáteční tok balistického fotonu, tak i pozdější difúzní šíření.

Metoda Monte Carlo je nutně statistická, a proto k dosažení přesnosti vyžaduje značnou dobu výpočtu. Kromě toho mohou simulace Monte Carlo sledovat více fyzických veličin současně s libovolným požadovaným prostorovým a časovým rozlišením. Díky této flexibilitě je modelování Monte Carlo mocným nástrojem. I když jsou výpočtově neefektivní, metody Monte Carlo jsou často považovány za standard pro simulovaná měření přenosu fotonů pro mnoho biomedicínských aplikací.

Monte Carlo simulace dopadu paprsku tužky na polo nekonečné rozptylové médium.

Biomedicínské aplikace metod Monte Carlo

Biomedicínské zobrazování

Optické vlastnosti biologické tkáně nabízejí vzrušující přístup k biomedicínskému zobrazování. Existuje mnoho zajímavých endogenních kontrastů, včetně absorpce z krve a melaninu a rozptylu z nervových buněk a jader rakovinných buněk. Kromě toho lze fluorescenční sondy cílit na mnoho různých tkání. Mikroskopické techniky (včetně konfokální, dva fotony, a optická koherentní tomografie ) mají schopnost zobrazovat tyto vlastnosti s vysokým prostorovým rozlišením, ale protože se spoléhají na balistické fotony, je jejich hloubková penetrace omezena na několik milimetrů. Zobrazování hlouběji do tkání, kde byly fotony mnohonásobně rozptýleny, vyžaduje hlubší pochopení statistického chování velkého počtu fotonů v takovém prostředí. Metody Monte Carlo poskytují flexibilní rámec, který byl použit různými technikami k rekonstrukci optických vlastností hluboko ve tkáni. Zde je uveden krátký úvod k několika z těchto technik.

• Fotoakustická tomografie V PAT je absorbováno difúzní laserové světlo, které generuje místní nárůst teploty. Tato lokální teplotní změna zase generuje ultrazvukové vlny prostřednictvím termoelastické expanze, které jsou detekovány pomocí ultrazvukového měniče. V praxi se mění nejrůznější parametry nastavení (tj. Vlnová délka světla, numerická apertura snímače) a ve výsledku je modelování Monte Carlo cenným nástrojem pro predikci tkáňové reakce před experimentálními metodami.
• Difúzní optická tomografie DOT je zobrazovací technika, která využívá k měření optických vlastností biologických tkání řadu blízkých infračervených světelných zdrojů a detektorů. Lze měřit různé kontrasty, včetně absorpce v důsledku oxy- a deoxy-hemoglobinu (pro funkční neuroimaging nebo detekci rakoviny) a koncentrace fluorescenčních sond. Aby bylo možné rekonstruovat obraz, je třeba znát způsob, jakým světlo putovalo z daného zdroje do daného detektoru a jak měření závisí na distribuci a změnách optických vlastností (známé jako dopředný model). Vzhledem k vysoce rozptylové povaze biologické tkáně jsou takové cesty komplikované a funkce citlivosti jsou rozptýlené. Dopředný model je často generován pomocí metod Monte Carlo.

Radiační terapie

Cílem radiační terapie je dodávat energii, obvykle ve formě ionizujícího záření, do rakovinné tkáně, zatímco šetří okolní normální tkáň. Modelování Monte Carlo se běžně používá v radiační terapii k určení periferní dávky, kterou pacient pocítí v důsledku rozptylu, a to jak z tkáně pacienta, tak z rozptylu z kolimace proti proudu v lineárním urychlovači.

Fotodynamická terapie

v Fotodynamická terapie (PDT) světlo se používá k aktivaci chemoterapeutických látek. Vzhledem k povaze PDT je ​​užitečné použít metody Monte Carlo pro modelování rozptylu a absorpce v tkáni, aby se zajistilo dodávání vhodných úrovní světla k aktivaci chemoterapeutických látek.

Implementace přenosu fotonů v rozptylovacím médiu

Zde je představen model fotonové metody Monte Carlo v homogenním nekonečném médiu. Tento model je však snadno rozšířitelný pro vícevrstvá média. U nehomogenního média musí být brány v úvahu hranice. Kromě semi-nekonečného média (ve kterém jsou fotony považovány za ztracené, pokud opustí horní hranici) je třeba věnovat zvláštní pozornost. Další informace naleznete na odkazech v dolní části stránky. Problém vyřešíme pomocí nekonečně malého bodového zdroje (analyticky znázorněného jako a Diracova delta funkce v prostoru a čase). Odpovědi na libovolné zdrojové geometrie lze sestavit pomocí metody Greenovy funkce (nebo konvoluce, pokud existuje dostatečná prostorová symetrie). Požadované parametry jsou absorpční koeficient, koeficient rozptylu a funkce fáze rozptylu. (Pokud jsou brány v úvahu hranice, musí být také uveden index lomu pro každé médium.) Odezvy s časovým rozlišením se naleznou sledováním celkové uplynulé doby letu fotonu pomocí délka optické dráhy. Odpovědi na zdroje s libovolnými časovými profily lze poté modelovat pomocí konvoluce v čase.

V našem zjednodušeném modelu používáme následující techniku ​​snižování odchylek ke snížení výpočetního času. Místo individuálního šíření fotonů vytváříme fotonový paket se specifickou hmotností (obvykle inicializovanou jako jednota). Jak foton interaguje v zakaleném médiu, bude ukládat váhu v důsledku absorpce a zbývající váha bude rozptylována do dalších částí média. Během cesty lze zaznamenat libovolný počet proměnných, v závislosti na zájmu konkrétní aplikace. Každý fotonový paket opakovaně projde následujícími očíslovanými kroky, dokud nebude buď ukončen, odražen nebo přenesen. Proces je schématicky znázorněn vpravo. Lze spustit a modelovat libovolný počet fotonových paketů, dokud výsledná simulovaná měření nemají požadovaný poměr signál-šum. Všimněte si, že jelikož modelování Monte Carlo je statistický proces zahrnující náhodná čísla, budeme proměnnou ξ v celém textu používat jako a pseudonáhodné číslo pro mnoho výpočtů.

Schéma modelování toku fotonů v nekonečném rozptylovacím a absorpčním médiu pomocí simulací Monte Carlo.

Krok 1: Spuštění fotonového paketu

V našem modelu ignorujeme počáteční zrcadlovou odrazivost spojenou se vstupem do média, které neodpovídá indexu lomu. S ohledem na to jednoduše musíme nastavit počáteční polohu fotonového paketu a počáteční směr. Je vhodné použít globální souřadnicový systém. Použijeme tři Kartézské souřadnice k určení polohy, spolu se třemi směrové kosiny k určení směru šíření. Počáteční počáteční podmínky se budou lišit v závislosti na aplikaci, ale pro tužkový paprsek inicializovaný v počátku můžeme počáteční kosinus polohy a směru nastavit následujícím způsobem (izotropní zdroje lze snadno modelovat náhodným výběrem počátečního směru každého paketu):

${displaystyle {egin {aligned} x & = 0 {ext {Position:}} y & = 0 z & = 0 mu _ {x} & = 0 {ext {Direction cosines:}} mu _ {y} & = 0 mu _ {z} & = 1end {aligned}}}$

Krok 2: Výběr velikosti kroku a pohyb fotonového paketu

Velikost kroku, s, je vzdálenost, kterou fotonový paket urazí mezi interakčními místy. Existuje řada metod pro výběr velikosti kroku. Níže je uvedena základní forma výběru velikosti kroku fotonu (odvozená pomocí metoda inverzní distribuce a Pivo – Lambertův zákon ) ze kterého používáme pro náš homogenní model:

${displaystyle s = - {frac {ln xi} {mu _ {t}}}}$

kde ${displaystyle xi}$ je náhodné číslo a ${displaystyle {mu _ {t}}}$ je celkový koeficient interakce (tj. součet absorpčních a rozptylových koeficientů).

Jakmile je vybrána velikost kroku, fotonový paket se šíří o vzdálenost s ve směru definovaném kosiny směru. Toho lze snadno dosáhnout aktualizací souřadnic následujícím způsobem:

${displaystyle {egin {aligned} x & leftarrow x + mu _ {x} s y & leftarrow y + mu _ {y} s z & leftarrow z + mu _ {z} odeslat {aligned}}}$

Krok 3: Absorpce a rozptyl

Část hmotnosti fotonu je absorbována v každém interakčním místě. Tento zlomek hmotnosti se stanoví takto:

${displaystyle Delta W = {frac {mu _ {a}} {mu _ {t}}} W}$

kde ${displaystyle {mu _ {a}}}$ je absorpční koeficient.

Hmotnostní podíl lze poté zaznamenat do pole, pokud je pro konkrétní studii zajímavé distribuční absorpce. Hmotnost fotonového paketu musí být poté aktualizována následovně:

${displaystyle Wleftarrow W-Delta W,}$

Po absorpci je fotonový paket rozptýlen. Vážený průměr kosinu úhlu rozptylu fotonu je známý jako rozptylová anizotropie (G), který má hodnotu mezi −1 a 1. Pokud je optická anizotropie 0, obecně to znamená, že rozptyl je izotropní. Li G se blíží hodnotě 1, což znamená, že rozptyl je primárně ve směru dopředu. Abychom určili nový směr fotonového paketu (a tedy kosiny směru fotonů), musíme znát fázovou rozptylovou funkci. Často se používá fázová funkce Henyey-Greenstein. Poté se pomocí následujícího vzorce určí úhel rozptylu, θ.

${displaystyle cos heta = {egin {cases} {frac {1} {2g}} left [1 + g ^ {2} -left ({frac {1-g ^ {2}} {1-g + 2gxi}} ight) ^ {2} ight] & {ext {if}} geq 0 1-2xi & {ext {if}} g = 0end {cases}}}$

A polární úhel φ obecně se předpokládá, že je rovnoměrně rozloženo mezi 0 a ${displaystyle 2pi}$. Na základě tohoto předpokladu můžeme nastavit:

${displaystyle varphi = 2pi xi {frac {} {}}}$

Na základě těchto úhlů a původních směrových kosinů můžeme najít novou sadu směrových kosinů. Nový směr šíření lze v globálním souřadnicovém systému reprezentovat takto:

${displaystyle {egin {aligned} mu '_ {x} & = {frac {sin heta (mu _ {x} mu _ {z} cos varphi -mu _ {y} sin varphi)} {sqrt {1-mu _ {z} ^ {2}}}} + mu _ {x} cos heta mu '_ {y} & = {frac {sin heta (mu _ {y} mu _ {z} cos varphi + mu _ {x } sin varphi)} {sqrt {1-mu _ {z} ^ {2}}}} + mu _ {y} cos heta mu '_ {z} & = - {sqrt {1-mu _ {z} ^ {2}}} sin heta cos varphi + mu _ {z} cos heta end {zarovnáno}}}$

Pro zvláštní případ

${displaystyle {egin {aligned} mu _ {z} = 1end {aligned}}}$

použití

${displaystyle {egin {aligned} mu '_ {x} & = sin heta cos varphi mu' _ {y} & = sin heta sin varphi mu '_ {z} & = cos heta end {align}}}$

nebo

${displaystyle {egin {aligned} mu _ {z} = - 1end {aligned}}}$

použití

${displaystyle {egin {aligned} mu '_ {x} & = sin heta cos varphi mu' _ {y} & = - sin heta sin varphi mu '_ {z} & = - cos heta end {zarovnáno} }}$

C-kód:

/ ************************ Indicatrix ********************** Nový směr kosiny po rozptylu podle úhlu theta, fi. * mux new = (sin (theta) * (mux * muz * cos (fi) -muy * sin (fi))) / sqrt (1-muz ^ 2) + mux * cos ( theta) * muy new = (sin (theta) * (muy * muz * cos (fi) + mux * sin (fi))) / sqrt (1-muz ^ 2) + muy * cos (theta) * muz new = - sqrt (1-muz ^ 2) * sin (theta) * cos (fi) + muz * cos (theta) * ----------------------- ---------------------------------- * Vstup: * muxy, muys, tlamy - směr kosinu před srážkou * mutheta , fi - kosinus polárního úhlu a azimutálního úhlu * -------------------------------------- ------------------- * Výstup: * muxd, muyd, muzd - směr kosinus po kolizi * ---------------- ----------------------------------------- * / void Indicatrix (double muxs, double muys, double muz, double mutheta, double fi, double * muxd, double * muyd, double * muzd) {double costheta = mutheta; double sintheta = sqrt (1,0-costheta * costheta); // sin (theta) double sinfi = sin (fi); dvojitý cosfi = cos (fi); if (muzs == 1,0) {* muxd = sintheta * cosfi; * muyd = sintheta * sinfi; * muzd = costheta; } elseif (muzs == -1.0) {* muxd = sintheta * cosfi; * muyd = -sintheta * sinfi; * muzd = -costheta; } else {double denom = sqrt (1,0-tlamy * tlamy); dvojitý muzcosfi = tlamy * cosfi; * muxd = sintheta * (muxy * muzcosfi-muys * sinfi) / denom + muxy * costheta; * muyd = sintheta * (muys * muzcosfi + muxs * sinfi) / denom + muys * costheta; * muzd = -denom * sintheta * cosfi + tlamy * costheta; }}

Krok 4: Ukončení fotonu

Pokud fotonový paket zažil mnoho interakcí, pro většinu aplikací má zbývající váha v paketu malý dopad. Ve výsledku je nutné určit prostředky pro zakončení fotonových paketů dostatečně malé hmotnosti. Jednoduchá metoda by použila prahovou hodnotu a pokud je hmotnost fotonového paketu pod prahovou hodnotou, paket je považován za mrtvý. Výše uvedená metoda je omezená, protože nešetří energií. Aby byla celková energie konstantní, a Ruská ruleta technika se často používá pro fotony pod určitou prahovou hodnotou. Tato technika používá ruletovou konstantu m určit, zda foton přežije či nikoli. Fotonový paket má jednu šanci m přežít, v takovém případě bude mít novou váhu mW kde Ž je počáteční hmotnost (tato nová hmotnost v průměru šetří energii). Jindy je hmotnost fotonu nastavena na 0 a foton je ukončen. To je matematicky vyjádřeno níže:

${displaystyle W = {egin {cases} mW & xi leq 1 / m 0 & xi> 1 / mend {cases}}}$

Jednotky grafického zpracování (GPU) a rychlé simulace přenosu fotonů v Monte Carlu

Simulace migrace fotonů v zakaleném médiu v Monte Carlu je vysoce paralelizovatelným problémem, kdy se velké množství fotonů šíří nezávisle, ale podle stejných pravidel a různých sekvencí náhodných čísel. Díky paralelní povaze tohoto speciálního typu simulace Monte Carlo je velmi vhodný pro provádění na grafické procesorové jednotce (GPU). Vydání programovatelných GPU zahájilo takový vývoj a od roku 2008 existuje několik zpráv o použití GPU pro vysokorychlostní simulaci migrace fotonů v Monte Carlu.^[1][2][3][4]

Tento základní přístup lze sám paralelizovat použitím více GPU propojených dohromady. Jedním z příkladů je „GPU Cluster MCML“, který lze stáhnout z webu autorů (Monte Carlo Simulation of Light Transport in Multi-layered Turbid Media Based on GPU Clusters):http://bmp.hust.edu.cn/GPU_Cluster/GPU_Cluster_MCML.HTM

Viz také

• Rovnice radiačního přenosu a teorie difúze pro transport fotonů v biologické tkáni
• Metoda Monte Carlo
• Konvoluce pro optické širokopásmové reakce v rozptylovém médiu
• Metody Monte Carlo pro elektronový transport

Odkazy na další zdroje Monte Carlo

• Optical Imaging Laboratory na Washingtonské univerzitě v St. Louis (MCML)
• Oregonské lékařské laserové centrum
• Fotonová migrace Výzkum v Monte Carlu na Lund University ve Švédsku GPU akcelerace simulací Monte Carlo a škálovatelné Monte Carlo. Otevřete zdrojový kód ke stažení.
• Cloudové Monte Carlo pro lehký transport v zakaleném rozptylovacím médiu Tento nástroj je zdarma k použití ve výzkumu a nekomerčních činnostech.
• Light Transport in Tissue as an Example of Monte Carlo Simulation (with C ++ source code).

Reference

• Wang, L-H; Wu Hsin-I (2007). Biomedicínská optika: principy a zobrazování. Wiley.
• L.-H. Wang; S. L. Jacques; L.-Q. Zheng (1995). "MCML - Monte Carlo modelování přenosu světla ve vícevrstvých tkáních". Počítačové metody a programy v biomedicíně. 47 (2): 131–146. doi:10.1016 / 0169-2607 (95) 01640-F.
• L.-H. Wang; S. L. Jacques; L.-Q. Zheng (1997). „Konvoluce - konvoluce pro reakce na fotonový paprsek konečného průměru dopadající na vícevrstvé tkáně“ (PDF). Počítačové metody a programy v biomedicíně. 54 (3): 141–150. doi:10.1016 / S0169-2607 (97) 00021-7.
• S. L. Jacques; L.-H. Wang (1995). „Monte Carlo modelování přenosu světla v tkáních“ (PDF). V A. J. Welch; M. J. C. van Gemert (eds.). Optická tepelná odezva laserem ozářené tkáně. New York: Plenum Press. str. 73–100.
• L.-H. Wang; S. L. Jacques (1994). „Optimalizované radiální a úhlové polohy v modelování Monte Carlo“ (PDF). Lékařská fyzika. 21 (7): 1081–1083. Bibcode:1994MedPh..21.1081W. doi:10.1118/1.597351. PMID  7968840.

Vložené odkazy