第二章 理理論論基礎
2.2 蒙地卡羅羅演算法
蒙地卡羅羅演算法[22]為一種利利用大量量隨機取樣以模擬複雜系統現象或結果的 方法。此方法由 John von Neumann，Stanislaw Ulam 和 Nicholas Metropolis 等人在 1940 年年代核子武器計畫中為研究中子與原子核的交互作用所發明。在預 測 交 互 作 用 複 雜 的 系 統 時 ， 時 常 不不 存 在 確 定 性 演 算 法 (Deterministic algorithms)，透過蒙地卡羅羅演算法，配合系統中所有可能交互作用的適當機率率率 模型，其綜合性的結果將可以被準確預測。在光學的研究中，蒙地卡羅羅法在光子
7
傳遞行行為的模擬上是一個具彈性卻嚴謹的演算法，可同時考慮光子與組織交互作 用中各種可能的變數數，並且可避免傳統輻輻射傳導方程式 (radiative transfer equation)難以獲得解析解而必須以近似形式計算的問題，以獲得更更加精精準的預 測結果。
在 2.1 中描述光子與介質的作用包含吸收與散射等光學事件。當一束光子進 入介質，它會經歷歷一連連串串散射與反射事件，不不斷改變其行行徑方向，最後以被吸收、
離離開介質表面或是消逝在介質深處三種結果收場。吾人所觀察到的漫反射現象，
即由所有離離開介質表面的光子所構成。在蒙地卡羅羅的模擬中，我們把單一光子路路 徑看作一連連串串隨機長度度步及其隨之而來來的事件所組成。
A. 隨機步的長度度
在一特定介質中，光子在兩兩光學事件之間的平均行行徑距離離與其散射係數數(μs) 及吸收係數數(μa)有關。其單一路路徑事件長由機率率率決定，我們稱之為隨機步。若若用 s 表示隨機步長度度，其機率率率密度度p(s)可被描述為：
(2.1)
其中 μt=μs +μa，是為此介質的衰減係數數。
在 軟 體 的 實 作 上 ， 我 們 需 要 將 系 統 提 供 的 均 勻 隨 機 變 數數 (uniform distribution)映射至此機率率率分佈函數數以產生在此機率率率分佈函數數下的隨機亂亂數數。
在此我們使用其機率率率累累積分佈函數數(cumulative distribution function)：
(2.2)
對於任何機率率率分佈函數數，其機率率率累累積分佈函數數值域皆為[0,1]，正好可以映射至 系統中範圍[0,1]的均勻隨機變數數。設在系統中取得的隨機亂亂數數值為 ξ，則可令令：
(2.3)
8
經移項後可得其隨機路路徑長s 為：
(2.4)
在每一段隨機路路徑的終點，光子可能面臨臨吸收或散射兩兩種事件，抑或在行行徑中遭 遇介質交界面而造成折射或反射。這些事件，影響著該光子的存亡及行行進方向。
接下來來將敘述蒙地卡羅羅模擬中對於此些光子狀狀態的描述。
B. 光子傳播方向
光子在介質中的傳遞俱有三個自由度度。在以入射光源為原點的三維空間介質 模擬中，我們定義其方向向量量[μx, μy, μz ]。在光子傳遞的過程中，可能造成 其傳遞方向改變的因素包含散射、反射及折射。
I. 散射事件
當光子的散射事件發生時，其散射方位角及偏折角由機率率率決定。散射方位角 ψ 在各方向機率率率密度度相等，可將其描述為值域在[0 ,2π]的均勻散佈函數數，由系 統均勻隨機亂亂數數 ξ 線性轉換而得：
(2.5) 描 述 散射 偏 折角機 率率率 分佈 的 函數數 稱作 散 射 相位 函 數數 (scattering `phase function)。類類似於光子隨機步長度度的亂亂數數產生方法，在蒙地卡羅羅的實作上使用 相位函數數的累累積機率率率分佈
(2.6)
對應到系統中值域[0,1]的均勻隨機亂亂數數 ξ：
(2.7)
9
則可得均勻隨機變數數 ξ 與偏折角隨機變數數 θ 的轉換函數數F(ξ)：
(2.8)
在本研究中使用 Henyey-Greenstein 相位函數數及由 Kortun 等人利利用時域有限差分 法(finite-difference time-domain)配合細胞模型所推得的組織相位函數數[23]。
當隨機方位角及偏折角產生後，可計算出新的光子傳遞方向向量量[μ’x,μ’y,μ’z]：
(2.9)
其中[μx, μy, μz ]代表原光子傳遞方向向量量。
II. 折射與反射
除跟隨於隨機步後的散射事件，若若一光子在其行行進隨機步過程中遭遇介質交 界面，則會發生折射或反射，造成光子行行進方向的改變。當一光子由折射率率率 n1介 質入射其與折射率率率 n2介質的交界面，若若 n1 > n2，且其入射角大於臨臨界角 φ：
(2.10)
則發生全反射。光子的新方向向量量[μ’x,μ’y,μ’z ] = [μx, μy, -μz ]。
在其他狀狀況下，光子可能發生折射或反射。根據 Fresnel Equations，垂直 (Rs) 與平行行(Rp)入射平面的兩兩種極化光發生反射的光子比例例分別為：
(2.11)
10
其中，θ1為該折射事件的入射角，n1 及 n2 分別為兩兩介質的折射率率率。針對為非極 化光，其反射率率率可視為兩兩者的平均：
(2.12)
在此情況下，一光子的反射與否可由系統中值域為[0,1]的均勻隨機亂亂數數 ξ 決定：
若若 ξ < R，則該光子遭到反射，其反射後之方向向量量同於全反射；若若 ξ > R，
則該光子進入介質 2 且發生折射，其折射後的方向向量量與原方向向量量的關係式為：
(2.13)
C. 光子權重
一顆光子在介質的傳遞過程中，有一定的機率率率因為被吸收而消失。在一個隨 機步結束後，一光子因被吸收而消失的機率率率為 μa/(μs +μa)。在蒙地卡羅羅的模 擬中採用了了替代性的做法來來描述這種現象：光子在隨機步結束之後並不不因單次的 吸收而整顆消失，而是在每一顆光子初始時賦予其一生命權重 w，並在往後每一 個隨機步固定因吸收而損失部分的權重 ΔW：
(2.14)
若若光子在此衰減過程中權重低於某一特定量量值Wth，則僅有 1/m 的機率率率存活。其命 運交由在此處被稱作俄羅羅斯輪輪盤(Russian roulette)的均勻隨機亂亂數數 ξ 決定：
(2.15)
其中為維持系統中光子整體權重不不變，在使用俄羅羅斯輪輪盤後存活的光子，其權重 被重新設為mW。這樣的過程持續至光子離離開該介質而終止。
11