CN101286187A - 光在生物组织中传输特性的定量蒙特卡罗模拟方法 - Google Patents

Abstract

光在生物组织中传输特性的定量蒙特卡罗模拟方法属于光谱技术应用，该方法先将目标生物组织的空间结构描述为一个三维数字矩阵，将入射光源表征为设定数目的光子的集合，将入射光源位置和入射光方向赋给每个光子的初始位置和方向；追踪每个光子的传输过程，然后输出光吸收矩阵和所有逸出光子的信息。本发明通过建立光在组织中传输的蒙特卡罗模型，定量地模拟出任意三维结构的生物组织中的光传输特性。该方法精度高，稳定性好，能够获得真实生物组织中光传输的丰富特性。应用该方法，可以优化生物组织光学检测/成像系统的设计，为光治疗剂量定量化提供精确信息，为光谱技术在生物组织中的检测区域定位及信号定量分析提供准确可靠的工具。

Description

光在生物组织中传输特性的定量蒙特卡罗模拟方法

技术领域

本发明属于光谱技术应用，数学仿真和生物医学工程领域，涉及一种光在组织中传输特性的定量蒙特卡罗(Monte Carlo)模拟方法。

背景技术

将光谱技术应用于生物组织的检测和成像，相对于现有的成熟技术，如PET，fMRI，肌电/心电/脑电等，具有非侵入测量，时间分辨率高，成本低廉等优点，对疾病诊断与监测，光学临床治疗、脑功能机制以及生理病理研究具有重要意义。近年来，光谱技术在生物组织检测成像的研究和临床应用中受到越来越广泛的好评，但是其准确性和可靠性还没有被广泛接受，其原因是光，尤其是最常用的近红外光，在生物组织中的复杂传输特性还不能定量获取，限制了光谱技术的应用发展。定量获取任意波长光在生物组织中的传输特性，将优化生物组织光学检测仪器的设计，促进光谱技术在临床诊断治疗领域的应用，推动医学成像，激光外科手术，激光动力学疗法和辐射剂量定量化等临床应用方面的发展，同时对生物组织的新陈代谢，生理结构状态和功能等科学问题的探索也具有重大理论和现实意义。

目前，获取光在生物组织中传输特性的方法主要分为两类。一类是解析扩散方程方法，该方法需要描述组织结构的边界条件，对于普遍具有三维复杂结构的真实生物组织，边界条件太复杂，难以获得扩散方程的解析解。另一类是蒙特卡罗模拟方法，该方法不需要边界条件，灵活且精度高，可用于获取光在任意生物组织中的传输特性。近年来，越来越多的研究者致力于建立获取光在组织中传输特性的蒙特卡罗方法。其中，‘MCML’方法(S.L.Jacques and L.Wang.″Monte Carlo modeling of light transport intissues，″in 1995 Proc OPTICAL-THERMAL RESPONSE OFLASER-IRRADIATED TISSUE Conf.，PP.73-99.和L.Wang，S.L.Jacques，and L.Zheng，″MCML--Monte Carlo modeling of light transport inmulti-layered tissues，″Comput.Methods.Programs.Biomed.，vol.47，pp.131-146，1995.)是国际上应用最为广泛，精度最高的方法，但是该方法只适用于柱状多层的简化生物组织，不能用于具有复杂三维结构的生物组织。此外，近年来国外研发出的可用于三维结构生物组织的新方法——‘tMCimg’(D.A.Boas，J.P.Culver，J.J.Stott，and A.K.Dunn，″Threedimensional Monte Carlo code for photon migration through complexheterogeneous media including the adult human head，″Opt.Express.，vol.10，no.3，pp.159-170，2002.)，具有组织模型限制和精度不高的缺点。

发明内容

本发明的目的在于提供一种光在生物组织中传输特性的定量蒙特卡罗模拟方法，该方法可应用于任意三维结构的生物组织，模拟精度高。

本发明提供的光在生物组织中传输特性的定量蒙特卡罗模拟方法，其步骤包括：

(1)将目标生物组织的空间结构描述为一个三维数字矩阵，即组织模型，矩阵中的元素对应目标生物组织的体素，每个元素的数值为标识组织类型的数字；设定目标生物组织的各类组织光学特性参数，包括吸收系数，散射系数，折射系数和各向异性因子；设定一个与组织模型大小相同的三维空矩阵，其中每个元素用于记录其对应生物组织中一个体素对光的吸收量；

(2)将入射光源表征为设定数目的光子的集合，将入射光源位置和入射光方向赋给每个光子的初始位置和方向；

(3)追踪每个光子的传输过程，首先投放光子，然后以由随机数确定的步长移动，在遇到边界的情况下发生反射或折射，一旦光子折射到生物组织外就停止追踪，否则在光子移动一个步长后发生光吸收而衰减能量，发生散射而改变方向，继续下一步移动的追踪直到光子能量衰减到足够小，应用俄罗斯轮回盘算法，如果光子在轮回盘中复活，继续追踪光子下一步的移动，如果在轮回盘中死亡，停止对该光子的追踪；

(4)追踪完所有光子后输出光吸收矩阵和所有逸出光子的信息。

本发明建立了一个新的光子在生物组织中传输特性的蒙特卡罗方法，可以高精度地定量模拟计算出任意波长的光在任意生物组织中的传输特性，包括多种参数值，例如吸收分布三维图、光通量分布三维图、漫射光能量分布三维图、透射光能量分布三维图、漫射光的方向分辨特性曲线，透射光的方向分辨特性曲线、总漫射量、总吸收量、总透射量或漫射/透射光在各层组织中路径长度等等。这些丰富的信息为生物组织光学检测/成像系统设计的优化提供了方便，为光治疗方法的剂量定量化提供精确的指导信息，为光谱技术在生物组织中检测区域定位以及信号定量分析提供了准确可靠的工具，为促进光谱技术的应用发展提供平台。

附图说明

图1为本发明提出的光在生物组织中传输特性的定量蒙特卡罗方法基本流程图。

图2为实施例中示意用的人脑组织结构剖面图。

图3为本发明的模拟方法得到的光在人前额脑中传输的光通量分布，相邻两子图间隔为4mm，其中子图g与图2位置对应。

图4为现有方法‘MCML’得到的光在人前额脑中传输的光通量分布。

图5为本发明的模拟方法与现有精度最高的方法‘MCML’，现有可用于三维结构组织的新方法‘tMCimg’得到的光在层状简化生物组织中的光通量分布和吸收特性曲线的比较结果。

具体实施方式

本发明提供的光在生物组织内传输特性的定量蒙特卡罗模拟方法，建立了光子在多体素生物组织中传输的蒙特卡罗模型，包括随机变量取样，组织结构的三维空间系统，光子表征，光子入射，光子步长，光子吸收，光子散射，不同组织之间的界面对光子迁移的作用，光子迁移终止，光子传输特性记录以及整个蒙特卡罗模型的流程。下面结合附图和实例对本发明作进一步详细的说明。

如图1所示，本发明的实施步骤如下：

(1)将目标生物组织的空间结构描述为一个三维数字矩阵，即组织模型。矩阵中的元素对应目标生物组织的体素，每个元素的数值为标识组织类型的数字。体素越小，描述组织结构的数学组织模型就越逼近于真实生物组织结构，模拟得到的光子传输特性精度就越高。实际模拟过程中，可以选取实际生物组织中最薄或截面积最小的那层组织的1/10厚度或1/10截面最小边长作为体素边长。对于脑活动检测或肌肉中有氧代谢活动检测中的模拟，可选择体素边长为0.5mm。

(2)设定目标生物组织的各类组织光学特性参数，包括吸收系数，散射系数，折射系数和各向异性因子。这些参数值的精度也会影响模拟精度，查阅已经发表的文献来获取这些参数是一较佳的实施例。值得注意的是，这些参数的选定要求与目标光源的波长一致。

(3)设定一个与组织模型大小相同的三维空矩阵，其中每个元素用来记录其对应生物组织中的一个体素对光的吸收量。

步骤(1)、(2)和(3)之间没有严格的时序，可以并列或次序互换。

(4)将入射光源表征为设定数目的光子的集合，将入射光源位置和入射光方向赋给每个光子的初始位置和方向。模拟使用的光子数越高，模拟计算的精度越高。其数值可设定为光源能量与探测光能量比值的10倍，一般取值10⁵以上。

(5)开始追踪每个光子的传输过程。对于每个光子，具体追踪过程如下：

(5.1)投放光子。如果光子在生物组织外，使用迭代算法将光子按照初始方向移动到生物组织表面，然后用菲涅尔定律计算出光子在组织表面发生折射后的光子能量和方向。如果光子在生物组织中，不改变光子位置，方向和能量。

(5.2)计算光子的一个备用参量-剩余步长Sleft：Sleft＝-ln(ξ)。

(5.3)根据光子当前位置在组织模型中提取对应元素的数值，然后根据该数值提取其表示的一种组织的光学参数，包括各项异性因子g，吸收系数μ_a，散射系数μ_s和折射系数n。

(5.4)计算光子移动一步的步长s

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

这里，ξ是由计算机伪随机数产生器产生的均匀分布在(0，1)的伪随机数，dx，dy，dz依次为组织模型中每个体素的长，宽，高。

(5.5)按照光子当前方向和(5.4)中计算出的步长伪移动光子一步，判断该过程中光子是否穿过不同组织之间的界面。如果是，进入步骤(5.6)；如果否，将光子移动一个步长，更新光子的当前位置，然后进入步骤(5.9)。

(5.6)确定光子在第一个界面上的撞击点。方法为：首先找出如下表达式中的最小项；

(max([x/dx]，[x₀/dx])-x/dx)/ux

(max([y/dy]，[y₀/dy])-y/dy)/uy

(max([z/dz]，[z₀/dz])-z/dy)/uz

这里，x，y，z表示光子当前位置，x₀，y₀，z₀是伪移动一个步长后的位置。如果第一项最小，那么光子在界面上的作用点表达式为：

x0＝max([x]，[x₀])

y0＝(x0-x)/ux·uy+y

z0＝(x0-x)/ux·uz+z

其他情况类推。

(5.7)计算光子在(5.6)中确定的撞击点与界面发生相互作用后的方向。具体方法为：首先计算界面法向矢量(a，b，c)，其表达式为：

$(a,b,c)=\frac{\left(\right[x]-[x_0],[y]-[y_0],[z]-[z_0\left]\right)}{{({\left(\right[x]-[x_0\left]\right)}^2+{\left(\right[y]-[y_0\left]\right)}^2+{\left(\right[z]-[z_0\left]\right)}^2)}^{1/2}}$

然后，根据该矢量表达式计算光子与界面发生作用的入射角θ和折射角θ_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₂表示遇到界面前后所在组织的折射系数。

最后分四种情况计算光子遇到界面之后的方向：

如果cosθ＝0，设置光子遇到界面之后方向不变；

如果cosθ＝1，当ε＞(n₂-n₁)²/(n₂+n₁)²，设置光子遇到界面之后方向不变；当ε≤(n₂-n₁)²/(n₂+n₁)²，设置光子发生反射，其方向变为：

$\left\{\begin{array}{c}{k}_{\mathrm{xi}}=k_x-2\cos \mathrm{\θ}\·a\\ k_{\mathrm{yi}}=k_y-2\cos \mathrm{\θ}\·b\\ k_{\mathrm{zi}}=k_z-2\cos \mathrm{\θ}\·c\end{array}\right.$

如果θ＞sin^-1(n₂/n₁)，设置光子发生全反射，方向变为：

$\left\{\begin{array}{c}{k}_{\mathrm{xi}}=k_x-2\cos \mathrm{\θ}\·a\\ k_{\mathrm{yi}}=k_y-2\cos \mathrm{\θ}\·b\\ k_{\mathrm{zi}}=k_z-2\cos \mathrm{\θ}\·c\end{array}\right.$

在其他情况下，设置光子与界面发生作用后的方向为：

其中，r的表达式为：

$r=\frac{1}{2}[{\left(\frac{n_1\·\cos {\mathrm{\θ}}_t-n_2\·\cos \mathrm{\θ}}{n_1\·\cos {\mathrm{\θ}}_t+n_2\·\cos \mathrm{\θ}}\right)}^2+{\left(\frac{n_1\·\cos \mathrm{\θ}-n_2\·\cos {\mathrm{\θ}}_t}{n_1\·\cos \mathrm{\θ}+n_2\·\cos {\mathrm{\θ}}_t}\right)}^2]$

(5.8)更新已移动的步长s为光子当前位置至步骤(5.7)中与界面发生作用的撞击点的距离，然后将光子当前位置更新到上述撞击点。

(5.9)将光子在当前已移动的步长中的吸收累加到步骤(3)中三维空矩阵的相应元素中。该元素与该步长所穿过的最后一个体素对应。光吸收量A的计算式为：A＝w(1-exp(-μ_a·s))。

如果光子在(5.7)中折射到空气中，即为光子逸出组织，设定光子死亡且更新Sleft＝0，并记录该光子的当前属性信息，包括位置、方向、能量、经过的路径长等；其它情况下更新Sleft：Sleft＝Sleft-s·μ_s。

如果Sleft＞0，进入步骤(5.3)；否则进入步骤(5.10)。

(5.10)如果光子已经死亡，停止对该光子的追踪；否则进入步骤(5.11)。

(5.11)光子发生散射。散射后光子方向(μ′_x，μ′_y，μ′_z)更新为：

${{\mathrm{\μ}}^{\′}}_x=\sin \mathrm{\α}({\mathrm{\μ}}_x{\mathrm{\μ}}_z\cos \mathrm{\ψ}-{\mathrm{\μ}}_y\sin \mathrm{\ψ})/\sqrt{1-{{\mathrm{\μ}}_z}^2}+{\mathrm{\μ}}_x\cos \mathrm{\α}$

${{\mathrm{\μ}}^{\′}}_y=\sin \mathrm{\α}({\mathrm{\μ}}_y{\mathrm{\μ}}_z\cos \mathrm{\ψ}-{\mathrm{\μ}}_x\sin \mathrm{\ψ})/\sqrt{1-{{\mathrm{\μ}}_z}^2}+{\mathrm{\μ}}_y\cos \mathrm{\α}$

${{\mathrm{\μ}}^{\′}}_z=-\sin \mathrm{\α}\cos \mathrm{\ψ}\sqrt{1-{{\mathrm{\μ}}_z}^2}+{\mathrm{\μ}}_z\cos \mathrm{\α}$

这里，(μ_x，μ_y，μ_z)表示光子散射前的方向，α和ψ分别表示光子散射前的偏转角和方位角，其表达式为：

$\cos \mathrm{\α}=\left\{\begin{array}{cc}[1+g^2-{(1-g^2/1+g+2\mathrm{g\ξ})}^2]/2g& \mathrm{ifg}\≠0\\ 2\mathrm{\ξ}-1& \mathrm{ifg}=0\end{array}\right.$

ψ＝2πξ

(5.12)比较光子当前能量与设定阈值的大小关系。该阈值表示可继续追踪光子传输过程的最小能量值，一般情况下为10^-4。如果比设定阈值大，进入步骤(5.2)，否则进入步骤(5.13)。

(5.13)应用俄罗斯轮回盘算法，如果光子在轮回盘中死亡，停止对该光子的追踪，否则进入步骤(5.2)。

(6)重复如上过程追踪完所有光子。然后，输出光吸收矩阵和所有逸出光子的信息。

(7)基于步骤(6)中的输出，进行统计计算，转换为其它所需获取的传输特性。例如：将光吸收量矩阵中各元素除以其对应体素的吸收系数，得到通量分布；将光吸收量矩阵所有元素的值相加，然后除以光子数目得到总吸收量；将所有在光源同侧逸出的光子(即漫射光子)当前能量相加，然后除以光子数目得到总漫射量；将所有在光源对侧逸出的光子(即透射光子)当前能量相加，然后除以光子数目得到总透射量；将光吸收矩阵中对应同一深度体素的所有元素相加，即得到光吸收随深度变化的曲线；将光通量矩阵中对应同一深度体素的所有元素相加，即得到光通量随深度变化的曲线；将漫射光子中与光子入射方向同一角度范围内的光子能量相加，即得到漫射光能量与漫射角度之间的关系曲线；将透射光子中与光子入射方向同一角度范围内的光子能量相加，即得到透射光能量与漫射角度之间的关系曲线；将漫射光子中与光子入射位置距离在同一范围内的光子能量相加，即得到漫射光能量与距离之间的关系曲线；将透射光子中与光子入射位置距离在同一范围内的光子能量相加，即得到透射光能量与距离之间的关系曲线；将所有漫射光子在一类组织中经过的路径加到一起，即得到漫射光在不同组织中的路径长度；将所有透射光子在一类组织中经过的路径加到一起，即得到透射光在不同组织中的路径长度。

对于使用本发明的方法对光在目标生物组织中传输特性进行模拟得到的输出文件，采用matlab进行读取、绘图、及参数转换是一较佳实施例。

下面通过实例进一步阐述本发明。

实施例1：

以体素边长为0.4mm，光子数为10⁷和下表中的组织光学参数使用本发明方法模拟光在如图2剖面图所示的数字人前额脑中的传输，所得的光通量分布如图3所示。图3中的子图g对应图2中的白色线框内区域。使用现有的新方法‘tMCimg’无法对该组织模型进行模拟，原因是该方法不能应用于表面弯曲的生物组织。以相同参数使用现有的受到广泛使用的‘MCML’方法所得的任意纵截面的光通量分布如图4所示。图3中不同横截面的分布图不同，与脑组织的三维组织结构吻合。图4中各纵截面中的光通量分布相同，只能反映光通量随深度的变化关系，不能反映脑组织结构对光通量分布的影响。比较图3和图4，表明本发明方法较之现有方法能更真实更稳定地计算光在生物组织中的传输特性。

实施例2

以柱状多层简化组织模型检验本发明方法的精度，并与现有精度最高的‘MCML’方法得到结果进行比较。体素边长为0.5mm，光子数为10⁵。光学参数如下表所示。

下表为本发明方法与现有精度最高的方法‘MCML’，以及可用于三维结构组织的新方法‘tMCimg’得到的光在层状简化生物组织中传输的漫射系数的比较结果。图5的上子图为光通量分布图，下子图为光吸收随生物组织深度变化的曲线。下表和图5均表明本发明在精度上不低于现有被认为精度最高的方法‘MCML’的水平。

表中数值为均值±方差

Claims (4)

1、 一种光在生物组织中传输特性的定量蒙特卡罗模拟方法，其步骤包括：

(1)将目标生物组织的空间结构描述为一个三维数字矩阵，即组织模型，矩阵中的元素对应目标生物组织的体素，每个元素的数值为标识组织类型的数字；设定目标生物组织的各类组织光学特性参数，包括吸收系数，散射系数，折射系数和各向异性因子；设定一个与组织模型大小相同的三维空矩阵，其中每个元素用于记录其对应生物组织中一个体素对光的吸收量；

(2)将入射光源表征为设定数目的光子的集合，将入射光源位置和入射光方向赋给每个光子的初始位置和方向；

(3)追踪每个光子的传输过程，首先投放光子，然后以由随机数确定的步长移动，在遇到边界的情况下发生反射或折射，一旦光子折射到生物组织外就停止追踪，否则在光子移动一个步长后发生光吸收而衰减能量，发生散射而改变方向，继续下一步移动的追踪直到光子能量衰减到足够小，应用俄罗斯轮回盘算法，如果光子在轮回盘中复活，继续追踪光子下一步的移动，如果在轮回盘中死亡，停止对该光子的追踪；

(4)追踪完所有光子后输出光吸收矩阵和所有逸出光子的信息。

2、 根据权利要求1所述的定量蒙特卡罗模拟方法，其特征在于：步骤(3)包括下述过程

(3.1)按照下述方式投放光子：如果光子在生物组织外，使用迭代算法将光子按照初始方向移动到生物组织表面，然后用菲涅尔定律计算出光子在组织表面发生折射后的光子能量和方向；如果光子在生物组织中，不改变光子位置，方向和能量；

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

(3.3)根据光子当前位置在组织模型中提取对应元素的数值，然后根据该数值提取其表示的一种组织的光学参数，包括各项异性因子g，吸收系数μ_a，散射系数μ_s和折射系数n；

(3.4)计算光子移动一步的步长s

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

其中，ξ是由计算机伪随机数产生器产生的均匀分布在(0，1)的伪随机数，dx，dy，dz依次为组织模型中每个体素的长，宽，高；

(3.5)按照光子当前方向和步骤(3.4)中计算出的步长伪移动光子一步，判断该过程中光子是否穿过不同组织之间的界面；如果是，进入步骤(3.6)；如果否，将光子移动一个步长，更新光子的当前位置，然后进入步骤(3.9)；

(3.6)确定光子在第一个界面上的撞击点；

(3.7)计算光子在步骤(3.6)中确定的撞击点与界面发生相互作用后的方向；具体方法为：首先计算界面法向矢量(a，b，c)，其表达式为：

其中，x，y，z表示光子当前位置，x₀，y₀，z₀是伪移动一个步长后的位置；然后，根据该矢量表达式计算光子与界面发生作用的入射角θ和折射角θ_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₂分别表示遇到界面前、后所在组织的折射系数；

最后分四种情况计算光子遇到界面之后的方向：

如果cosθ＝0，设置光子遇到界面之后方向不变；

如果cosθ＝1，当ε＞(n₂-n₁)²/(n₂+n₁)²，设置光子遇到界面之后方向不变；当ε≤(n₂-n₁)²/(n₂+n₁)²，设置光子发生反射，其方向变为：

$\left\{\begin{array}{c}{k}_{\mathrm{xi}}=k_x-2\cos \mathrm{\θ}\·a\\ k_{\mathrm{yi}}=k_y-2\cos \mathrm{\θ}\·b\\ k_{\mathrm{zi}}=k_z-2\cos \mathrm{\θ}\·c\end{array}\right.$

如果θ＞sin^-1(n₂/n₁)，设置光子发生全反射，方向为：

$\left\{\begin{array}{c}{k}_{\mathrm{xi}}=k_x-2\cos \mathrm{\θ}\·a\\ k_{\mathrm{yi}}=k_y-2\cos \mathrm{\θ}\·b\\ k_{\mathrm{zi}}=k_z-2\cos \mathrm{\θ}\·c\end{array}\right.$

在其他情况下，设置光子与界面发生作用后的方向分别为：

其中，r的表达式为：

$r=\frac{1}{2}[{\left(\frac{n_1\·\cos {\mathrm{\θ}}_t-n_2\·\cos \mathrm{\θ}}{n_1\·\cos {\mathrm{\θ}}_t+n_2\·\cos \mathrm{\θ}}\right)}^2+{\left(\frac{n_1\·\cos \mathrm{\θ}-n_2\·\cos {\mathrm{\θ}}_t}{n_1\·\cos \mathrm{\θ}+n_2\·\cos {\mathrm{\θ}}_t}\right)}^2]$

(3.8)更新已移动的步长s为光子当前位置至步骤(3.7)中与界面发生作用的撞击点的距离，然后将光子当前位置更新到上述撞击点；

(3.9)将光子在当前已移动的步长中的吸收累加到步骤(1)中三维空矩阵的相应元素中；该元素与该步长所穿过的最后一个体素对应；光吸收量A的计算式为：A＝w(1-exp(-μ_a·s))；

如果光子在步骤(3.7)中折射到空气中，即为光子逸出组织，设定光子死亡且更新Sleft＝0，并记录该光子的当前属性信息，包括位置、方向、能量、经过的路径长等；其它情况下更新Sleft：Sleft＝Sleft-s·μ_s；

如果Sleft＞0，进入步骤(3.3)；否则进入步骤(3.10)；

(3.10)如果光子已经死亡，停止对该光子的追踪，转入步骤(4)；否则进入步骤(3.11)；

(3.11)光子发生散射，散射后光子方向(μ′_x，μ′_y，μ′_z)更新为：

${{\mathrm{\μ}}^{\′}}_x=\sin \mathrm{\α}({\mathrm{\μ}}_x{\mathrm{\μ}}_z\cos \mathrm{\ψ}-{\mathrm{\μ}}_y\sin \mathrm{\ψ})/\sqrt{1-{{\mathrm{\μ}}_z}^2}+{\mathrm{\μ}}_x\cos \mathrm{\α}$

${{\mathrm{\μ}}^{\′}}_y=\sin \mathrm{\α}({\mathrm{\μ}}_y{\mathrm{\μ}}_z\cos \mathrm{\ψ}-{\mathrm{\μ}}_x\sin \mathrm{\ψ})/\sqrt{1-{{\mathrm{\μ}}_z}^2}+{\mathrm{\μ}}_y\cos \mathrm{\α}$

${{\mathrm{\μ}}^{\′}}_z=-\sin \mathrm{\α}\cos \mathrm{\ψ}\sqrt{1-{{\mathrm{\μ}}_z}^2}+{\mathrm{\μ}}_z\cos \mathrm{\α}$

这里，(μ_x，μ_y，μ_z)表示光子散射前的方向，α和ψ分别表示光子散射前的偏转角和方位角，其表达式为：

$\cos \mathrm{\α}=\left\{\begin{array}{ccc}[1+g^2-{(1-g^2/1+g+2\mathrm{g\ξ})}^2]/2g& \mathrm{if}& g\≠0\\ 2\mathrm{\ξ}-1& \mathrm{if}& g=0\end{array}\right.$

ψ＝2πξ

(3.12)比较光子当前能量与设定阈值的大小关系；如果比设定阈值大，进入步骤(3.2)，否则进入步骤(3.13)；

(3.13)应用俄罗斯轮回盘算法，如果光子在轮回盘中死亡，停止对该光子的追踪，转入步骤(4)，否则进入步骤(3.2)。

3、 根据权利要求2所述的定量蒙特卡罗模拟方法，其特征在于：步骤(3.6)按照下述方法确定光子在第一个界面上的撞击点：

首先找出如下表达式中的最小项；

(max([x/dx]，[x₀/dx])-x/dx)/ux

(max([y/dy]，[y₀/dy])-y/dy)/uy

(max([z/dz]，[z₀/dz])-z/dy)/uz

如果第一项最小，那么光子在界面上的撞击点表达式为：

x0＝max([x]，[x₀])

y0＝(x0-x)/ux·uy+y

z0＝(x0-x)/ux·uz+z

其他情况类推。

4、 根据权利要求1、2或3所述的定量蒙特卡罗模拟方法，其特征在于：该方法还包括步骤(5)，其过程为：基于步骤(4)中的输出，进行统计计算并转换为下述所需获取的传输特性中的一个或任意几个：

将光吸收量矩阵中各元素除以其对应体素的吸收系数，得到通量分布；

将光吸收量矩阵所有元素的值相加，然后除以光子数目得到总吸收量；

将所有在光源同侧逸出的光子当前能量相加，然后除以光子数目得到总漫射量；

将所有在光源对侧逸出的光子当前能量相加，然后除以光子数目得到总透射量；

将光吸收矩阵中对应同一深度体素的所有元素相加，即得到光吸收随深度变化的曲线；

将光通量矩阵中对应同一深度体素的所有元素相加，即得到光通量随深度变化的曲线；

将漫射光子中与光子入射方向同一角度范围内的光子能量相加，即得到漫射光能量与漫射角度之间的关系曲线；

将透射光子中与光子入射方向同一角度范围内的光子能量相加，即得到透射光能量与漫射角度之间的关系曲线；

将漫射光子中与光子入射位置距离在同一范围内的光子能量相加，即得到漫射光能量与距离之间的关系曲线；

将透射光子中与光子入射位置距离在同一范围内的光子能量相加，即得到透射光能量与距离之间的关系曲线；

将所有漫射光子在一类组织中经过的路径加到一起，即得到漫射光在不同组织中的路径长度；

将所有透射光子在一类组织中经过的路径加到一起，即得到透射光在不同组织中的路径长度。

Images