[發(fā)明專利]光在生物組織中傳輸特性的定量蒙特卡羅模擬方法

有效

• 申請?zhí)枺?/td>: 200810047962.7
• 申請日: 2008-06-10
• 公開（公告）號: CN101286187A
• 公開（公告）日: 2008-10-15
• 發(fā)明（設(shè)計）人: 駱清銘;龔輝;李婷
• 申請（專利權(quán)）人: 華中科技大學(xué)
• 主分類號: G06F19/00
• 分類號: G06F19/00;A61B5/00;A61N5/06
• 代理公司: 華中科技大學(xué)專利中心
• 代理人: 曹葆青

1、?一種光在生物組織中傳輸特性的定量蒙特卡羅模擬方法，其步驟包括：

(1)將目標(biāo)生物組織的空間結(jié)構(gòu)描述為一個三維數(shù)字矩陣，即組織模型，矩陣中的元素對應(yīng)目標(biāo)生物組織的體素，每個元素的數(shù)值為標(biāo)識組織類型的數(shù)字；設(shè)定目標(biāo)生物組織的各類組織光學(xué)特性參數(shù)，包括吸收系數(shù)，散射系數(shù)，折射系數(shù)和各向異性因子；設(shè)定一個與組織模型大小相同的三維空矩陣，其中每個元素用于記錄其對應(yīng)生物組織中一個體素對光的吸收量；

(2)將入射光源表征為設(shè)定數(shù)目的光子的集合，將入射光源位置和入射光方向賦給每個光子的初始位置和方向；

(3)追蹤每個光子的傳輸過程，首先投放光子，然后以由隨機(jī)數(shù)確定的步長移動，在遇到邊界的情況下發(fā)生反射或折射，一旦光子折射到生物組織外就停止追蹤，否則在光子移動一個步長后發(fā)生光吸收而衰減能量，發(fā)生散射而改變方向，繼續(xù)下一步移動的追蹤直到光子能量衰減到足夠小，應(yīng)用俄羅斯輪回盤算法，如果光子在輪回盤中復(fù)活，繼續(xù)追蹤光子下一步的移動，如果在輪回盤中死亡，停止對該光子的追蹤；

(4)追蹤完所有光子后輸出光吸收矩陣和所有逸出光子的信息。

2、?根據(jù)權(quán)利要求1所述的定量蒙特卡羅模擬方法，其特征在于：步驟(3)包括下述過程

(3.1)按照下述方式投放光子：如果光子在生物組織外，使用迭代算法將光子按照初始方向移動到生物組織表面，然后用菲涅爾定律計算出光子在組織表面發(fā)生折射后的光子能量和方向；如果光子在生物組織中，不改變光子位置，方向和能量；

(3.2)計算光子的剩余步長Sleft：Sleft＝-ln(ξ)；

(3.3)根據(jù)光子當(dāng)前位置在組織模型中提取對應(yīng)元素的數(shù)值，然后根據(jù)該數(shù)值提取其表示的一種組織的光學(xué)參數(shù)，包括各項異性因子g，吸收系數(shù)μ_a，散射系數(shù)μ_s和折射系數(shù)n；

(3.4)計算光子移動一步的步長s

s＝min(-ln(ξ)/μ_s，min(dx，dy，dz))

其中，ξ是由計算機(jī)偽隨機(jī)數(shù)產(chǎn)生器產(chǎn)生的均勻分布在(0，1)的偽隨機(jī)數(shù)，dx，dy，dz依次為組織模型中每個體素的長，寬，高；

(3.5)按照光子當(dāng)前方向和步驟(3.4)中計算出的步長偽移動光子一步，判斷該過程中光子是否穿過不同組織之間的界面；如果是，進(jìn)入步驟(3.6)；如果否，將光子移動一個步長，更新光子的當(dāng)前位置，然后進(jìn)入步驟(3.9)；

(3.6)確定光子在第一個界面上的撞擊點；

(3.7)計算光子在步驟(3.6)中確定的撞擊點與界面發(fā)生相互作用后的方向；具體方法為：首先計算界面法向矢量(a，b，c)，其表達(dá)式為：

其中，x，y，z表示光子當(dāng)前位置，x₀，y₀，z₀是偽移動一個步長后的位置；然后，根據(jù)該矢量表達(dá)式計算光子與界面發(fā)生作用的入射角θ和折射角θ_t，其計算式為：

cosθ＝a·k_x+b·k_y+c·k_z

cosθ_t＝(1-(1-cos²θ)·n₁²/n₂²)^1/2

其中，k_x，k_y，k_z表示光子遇到界面之前的方向余弦值，n₁和n₂分別表示遇到界面前、后所在組織的折射系數(shù)；

最后分四種情況計算光子遇到界面之后的方向：

如果cosθ＝0，設(shè)置光子遇到界面之后方向不變；

如果cosθ＝1，當(dāng)ε＞(n₂-n₁)²/(n₂+n₁)²，設(shè)置光子遇到界面之后方向不變；當(dāng)ε≤(n₂-n₁)²/(n₂+n₁)²，設(shè)置光子發(fā)生反射，其方向變?yōu)椋?/p>

kxi=kx-2cosθ·akyi=ky-2cosθ·bkzi=kz-2cosθ·c]]>

如果θ＞sin^-1(n₂/n₁)，設(shè)置光子發(fā)生全反射，方向為：

kxi=kx-2cosθ·akyi=ky-2cosθ·bkzi=kz-2cosθ·c]]>

在其他情況下，設(shè)置光子與界面發(fā)生作用后的方向分別為：

其中，r的表達(dá)式為：

r=12[(n1·cosθt-n2·cosθn1·cosθt+n2·cosθ)2+(n1·cosθ-n2·cosθtn1·cosθ+n2·cosθt)2]]]>

(3.8)更新已移動的步長s為光子當(dāng)前位置至步驟(3.7)中與界面發(fā)生作用的撞擊點的距離，然后將光子當(dāng)前位置更新到上述撞擊點；

(3.9)將光子在當(dāng)前已移動的步長中的吸收累加到步驟(1)中三維空矩陣的相應(yīng)元素中；該元素與該步長所穿過的最后一個體素對應(yīng)；光吸收量A的計算式為：A＝w(1-exp(-μ_a·s))；

如果光子在步驟(3.7)中折射到空氣中，即為光子逸出組織，設(shè)定光子死亡且更新Sleft＝0，并記錄該光子的當(dāng)前屬性信息，包括位置、方向、能量、經(jīng)過的路徑長等；其它情況下更新Sleft：Sleft＝Sleft-s·μ_s；

如果Sleft＞0，進(jìn)入步驟(3.3)；否則進(jìn)入步驟(3.10)；

(3.10)如果光子已經(jīng)死亡，停止對該光子的追蹤，轉(zhuǎn)入步驟(4)；否則進(jìn)入步驟(3.11)；

(3.11)光子發(fā)生散射，散射后光子方向(μ′_x，μ′_y，μ′_z)更新為：

μ′x=sinα(μxμzcosψ-μysinψ)/1-μz2+μxcosα]]>

μ′y=sinα(μyμzcosψ-μxsinψ)/1-μz2+μycosα]]>

μ′z=-sinαcosψ1-μz2+μzcosα]]>

這里，(μ_x，μ_y，μ_z)表示光子散射前的方向，α和ψ分別表示光子散射前的偏轉(zhuǎn)角和方位角，其表達(dá)式為：

cosα=[1+g2-(1-g2/1+g+2gξ)2]/2gifg≠02ξ-1ifg=0]]>

ψ＝2πξ

(3.12)比較光子當(dāng)前能量與設(shè)定閾值的大小關(guān)系；如果比設(shè)定閾值大，進(jìn)入步驟(3.2)，否則進(jìn)入步驟(3.13)；

(3.13)應(yīng)用俄羅斯輪回盤算法，如果光子在輪回盤中死亡，停止對該光子的追蹤，轉(zhuǎn)入步驟(4)，否則進(jìn)入步驟(3.2)。