说明：本文为《Monte Carlo modeling of light transport in multi-layered tissues in standard C》翻译，主要介绍Lihong Wang等人于1995编写的MCML程序，研究组织光学中的蒙特卡洛仿真。

原文网站：https://omlc.org/index.html

word文档及相关资料：https://download.csdn.net/download/Turn_on/13527370

0.介绍

蒙特卡洛模拟已用于解决各种物理问题。但是，没有简洁的公认定义。我们希望采用Lux等人的定义（1991）。在蒙特卡洛方法的所有应用中，都建立了一个随机模型，其中某个随机变量（或几个变量的组合）的期望值等于要确定的物理量的值。然后，通过代表上述引入的随机变量的多个独立样本的平均值来估计该期望值。为了构造一系列独立样本，使用跟随待估计变量分布的随机数。

蒙特卡洛对光子传播的模拟提供了一种灵活而严谨的方法来处理混浊组织中的光子传输。该方法描述了光子传播的局部规则，该规则在最简单的情况下表示为概率分布，该概率分布描述了光子与组织相互作用的位点之间的光子运动的步长以及发生散射事件时光子轨迹中的偏转角。该模拟可以同时对多个物理量进行评分。然而，该方法本质上是统计的，并且依赖于计算机计算大量光子的传播。因此，该方法需要大量的计算时间。

所需的光子数量在很大程度上取决于所提出的问题，所需的精度和所需的空间分辨率。例如，为了简单地从具有特定光学特性的组织中了解总漫反射率，通常约3,000个光子可以产生有用的结果。为了绘制圆柱对称问题中的光子φ（r，z）的空间分布，通常至少需要10,000个光子才能产生可??接受的答案。为了在更复杂的三维问题（例如，用直径埋没的血管照射组织的有限直径光束）中绘制空间分布图，所需的光子可能会超过100,000。在这些介绍性说明中要记住的一点是，蒙特卡洛模拟是严格的，但必须是统计的，因此需要大量的计算时间才能达到精度和分辨率。但是，该方法的灵活性使Monte Carlo建模成为功能强大的工具。

本文提出的蒙特卡洛模拟的另一方面值得强调。此处描述的模拟未将光子视为波动现象，并且忽略了诸如相位和极化之类的特征。这些模拟的动机是为了预测浑浊组织中的辐射能传输。光子被大多数组织多重散射，因此相位和极化很快被随机化，并且在能量传输中几乎没有作用。尽管蒙特卡洛模拟可能能够记录相位和极化，并能统计地处理波动现象，但本手册不会考虑这些问题。

蒙特卡洛模拟基于宏观光学特性，假定这些特性在小单位组织体积上均匀扩展。光子-组织相互作用部位之间的平均自由程通常在10-1000 μm范围内，而100 μm是可见光谱中非常典型的值（Cheong等，1990）。举个例子，模拟不处理单元内辐射能量分布的细节。

作为蒙特卡洛模拟的简单示例。我们想在图0.1中给出单个光子包的典型轨迹。光子位置（点）之间的每一步都是可变的，并且等于–ln（ξ）/（μa +μs），其中ξ是随机数，μa和μs分别是吸收系数和散射系数（在此示例中，μa＝0.5cm-1，μs＝ 15 cm-1，g ＝ 0.90。值g是散射的各向异性。随着光子在组织中的移动，其重量从初始值1减小，并在n步后等于a，其中a是反照率（a=μs/（μa+μs））。当光子撞击表面时，一部分光子重量会随反射率逃逸，剩余的重量会内部反射并继续传播。最终，光子重量下降到阈值水平以下，并且该光子的模拟终止。在此示例中，当剩余光子重量的最后一个显着部分在星号（*）指示的位置处的表面逸出时发生终止。通常计算许多光子轨迹（104到106）以产生对介质中光子分布的统计描述。

本手册大致分为两个主要部分。第一部分描述了组织中光子传输的蒙特卡洛模拟的原理，如何在ANSI标准C中实现模拟以及一些计算结果和验证。第二部分为用户提供了使用mcml和修改mcml以适合特殊需要的详细说明，其中mcml代表多层组织的Monte Carlo模拟（Monte Carlo simulations for multi-layered tissues）。附录提供了流程图和mcml的整个源代码，以及一些其他有用的信息。

图0.1. 通过蒙特卡洛模拟计算，一个光子通过均质介质的运动。

1.蒙特卡洛模拟描述

本文中描述的蒙特卡洛模拟处理垂直入射在多层组织上的无限窄光子束的传输。每层都是无限宽的，并由以下参数描述：厚度，折射率，吸收系数μa，散射系数μs和各向异性系数g。还必须给出顶部环境介质（例如空气）和底部环境介质（如果存在）的折射率。尽管真实组织永远不可能无限宽，但是如果它比光子分布的空间范围大得多，则可以这样对其进行处理。吸收系数μa被定义为每单位无限小光程的光子吸收概率，而散射系数μs被定义为每单位无限小光程的光子散射概率。为了简化表示，有时使用总相互作用系数μt，它是吸收系数μa和散射系数μs的总和。相应地，相互作用系数是指单位无限小光程中光子相互作用的概率。各向异性g是偏转角余弦值的平均值（请参见第3.5节）。

光子吸收，能量密度，反射率和透射率是要模拟的物理量。该模拟在三个维度上传播光子，记录由于空间阵列中每个栅格元素的吸收而引起的光子沉积，A（x，y，z）（每J传递能量的J/cm 3或cm-3），最后通过将沉积量除以局部吸收系数μa（以cm -1表示），计算注量φ（x，y，z）（每J传递能量的J/cm2或cm-2）：φ（x，y，z）= A（x，y，z）/μa。由于光子吸收和光子通量可以通过组织的局部吸收系数来回转换，因此我们仅以mcml为单位报告光子吸收（请参见第4.2节和第9.3节）。可以通过在另一个程序转换中转换光子吸收来获得光子通量。模拟还记录了顶部（和底部）表面的光子逸出，以局部反射率（和透射率）（cm–2 sr–1）表示（请参阅第4.1节）。

在mcml的第一个版本中，我们考虑圆柱对称组织模型。因此，我们选择在二维阵列A（r，z）中记录光子沉积，尽管此模拟的光子传播是在三维中进行的。

蒙特卡洛模拟中同时使用了三个坐标系。笛卡尔坐标系用于跟踪光子包。坐标系的原点是组织表面上的光子入射点，z轴始终是指向组织内部的表面的法线，因此xy平面在组织表面上（图1.1）。圆柱坐标系用于对内部光子吸收A（r，z）进行评分，其中r和z分别是圆柱坐标系的径向坐标和z轴坐标。直角坐标系和圆柱坐标系共享原点和z轴。圆柱坐标系的r坐标也用于漫反射率和总透射率。它们分别以Rd（r，α）和Tt（r，α）记录在组织表面上，其中α是光子出射方向与组织表面的法线（反射的–z轴和透射的z轴）之间的夹角。z轴与光子传播方向动态对齐的运动球坐标系用于采样光子包的传播方向变化。在该球坐标系中，首先采样由散射引起的偏转角θ和方位角ψ。然后，根据直角坐标系中的方向余弦更新光子方向（请参见第3.5节）。

为了吸收光子，在z和r方向上建立了二维均匀网格系统。网格线间距分别在z和r方向上为Δz和Δr。z方向和r方向的网格元素总数分别为Nz和Nr。对于漫反射率和透射率，在r和α方向上设置了二维均匀栅格系统。该栅格系统可以与栅格系统共享r方向，以吸收光子。因此，我们只需要为α方向的漫反射率和透射率建立一个额外的一维网格系统。在我们的仿真中，我们总是选择α的范围为[0，π/ 2]，即0≤α≤π/ 2。网格元素的总数为Nα。因此，网格线间隔为Δα=π/（2 Nα）。

这是一个合适的时间，我们在整个模拟过程中始终使用cm作为长度的基本单位，以保持一致性。例如，每个层的厚度以及在r和z方向上的网格线间距以cm为单位。吸收系数和散射系数以cm-1为单位。

图1.1. 在多层组织上建立的笛卡尔坐标系的示意图，y轴指向外部。

在某些讨论中，尽管在程序中使用索引来引用数组元素，数组仍将仅通过网格元素的位置（例如（r，z）或（r，α）），而不是通过网格元素的索引来引用。

蒙特卡洛方法，顾名思义（“掷骰子”），依赖于来自定义明确的概率分布的变量的随机采样。几本书（Cashew等，1959； Lux等，1991；和Kalos等，1986）为蒙特卡洛建模的原理提供了很好的参考。让我们简要回顾一下在蒙特卡洛模拟中采样随机变量的方法。

考虑随机变量χ，这是组织中光子传播的蒙特卡洛模拟所需要的。该变量可以是光子将在光子与组织相互作用的位置之间采取的可变步长，也可以是由于散射事件而使散射的光子可能经历的偏转角。有一个概率密度函数定义了间隔（a，b）上χ的分布。对概率密度函数进行归一化：

(2.1)

为了模拟传播，我们希望能够基于伪随机数生成器反复随机选择χ的值。计算机提供一个随机变量ξ，该变量均匀分布在时间间隔（0，1）上。该均匀分布的随机变量的累积分布函数为：

(2.2)

为了对通常不均匀分布的函数p（χ）进行采样，我们假设存在一个不递减的函数χ= f（ξ）（Kalos等，1986），它将ξ∈（0，1）映射到χ∈（a，b）（图2.1）。变量χ和变量ξ具有一对一的映射关系。随后让以下几率相等：

(2.3a)

或者

(2.3a)

根据累积分布函数的定义，可以将2.3b更改为累积分布函数方程：

(2.4)

根据等式(2.4)左侧的相应概率密度函数扩展累积分布函数 并应用公式 (2.2)的右侧，我们将2.4式转换为：

(2.5)

然后等式(2.5)被用来求解 以获得函数f(ξ1) 。如果假设函数χ=f(ξ) 不增加，则类似的推导将得出等式（2.5）为：

(2.6)

但是，由于（1-ξ1 ）和ξ1 具有相同的分布，因此它们可以互换。因此，（2.5）和等式（2.6）是等效的。在下一章中，（2.5）将被重复调用以采样传播变量。

从图2.1可以了解整个采样过程。使用ξ进行蒙特卡洛选择χ的关键是使ξ在区间[0，ξ1]中的概率等于χ在区间[a，χ1]中的概率。在图2.1中，我们将描绘[0，χ1]上的p（χ）积分的阴影区域等同于描绘[0，ξ1]上的p（ξ）积分的阴影区域。请记住，曲线p（χ）和p（ξ）下的总面积均等于1，这适用于概率密度函数。结果是基于图2.1中阴影区域的相等性，在上边界ξ1和χ1之间进行了一对一映射。换句话说，我们将Fχ（χ1）等于Fξ（ξ1）。（方程式2.4）等同于（方程式2.5）。箭头表示变换过程χ1 = f（ξ1）。对于每个ξ1，选择χ1使得ξ1和χ1的累积分布函数具有相同的值。相应地，阴影线区域相等。也可以在图2.1中看到单调函数f（ξ）始终存在，因为ξ和χ的累积分布函数都是单调的。

图2.1. 基于均匀分布的随机变量ξ对随机变量χ进行采样。

举个例子，考虑光子运动步长s的采样，这将在第3.2节中全面讨论。概率密度函数为：

(2.7)

相互作用系数μt等于μa+μs。在等式2.5中使用此功能。根据随机数ξ得出采样值s1的表达式：

(2.8)

求解值s1：

(2.9a)

如方程式2.6中所述，以上表达式等价于：

(2.9b)

一种简单的减少方差的技术，即隐式光子捕获，用于提高蒙特卡洛模拟的效率。这种技术允许一个人同时沿着一条特定路径等效地传播许多光子作为一个包。每个光子数据包最初分配的权重W等于1。光子从原点正交注入组织，该光子对应于准直的任意狭窄的光子束。

光子的当前位置由笛卡尔坐标（x，y，z）指定。当前的光子方向由单位矢量r指定，单位矢量r可以等效地由方向余弦（μx，μy，μz）描述：

(3.1)

其中x ，y 和z 是沿每个轴的单位向量。光子位置初始化为（0，0，0），方向余弦设置为（0，0，1）。事实证明，这种对笛卡尔坐标系中光子位置和方向的描述（Witt，1977）比圆柱坐标系中的光子位置和方向更简单（Keijzer等人，1989）。

当发射光子时，如果组织表面的边界不匹配，则将发生一些镜面反射。如果外部介质和组织的折射率分别为n1和n2，则指定镜面反射率Rsp（Born等，1986； Hecht，1987）：

(3.2a)

如果第一层是玻璃，位于折射率为n3的一层介质的顶部，则考虑在玻璃层的两个边界上多次反射和透射。然后通过以下方式计算镜面反射率：

(3.2b)

其中r1和r2是玻璃层两个边界上的菲涅耳反射率：

(3.3)

(3.4)

注意，如果镜面反射率定义为光子被反射而不与组织相互作用的概率，则等式3.2a和3.2b并不是严格正确的，尽管它们可能是厚组织的真实镜面反射率的很好的估计。如果要严格区分镜面反射率和漫反射率，则可以跟踪光子数据包经历的交互次数。当我们对反射率进行评分时，如果交互次数不为零，则反射率就是漫反射率。否则，它是镜面反射率。透射率可以类似地区分。

光子重量减少Rsp，镜面反射率Rsp将报告到输出数据文件中。

(3.5)

光子包的步长是根据光子自由路径s∈[0，∞）的概率分布的采样计算得出的，这意味着0≤s

(3.6a)

或者：

(3.6b)

上面的方程 3.6b可以在（0，s1）范围内的s'上积分，并得到指数分布，其中使用P{s≥0}=1：

(3.7)

等式3.7可以重新排列以产生自由路径s的累积分布函数：

(3.8)

可以将累积分布函数分配给第二章中讨论的均匀分布的随机数ξ。可以重新排列方程式，以提供一种选择步长的方法：

(3.9a)

或用ξ代替（1-ξ）：

(3.9b)

等式 3.9b给出了光子-组织相互作用位点之间的平均自由程为1 /μt，因为–ln（ξ）的统计平均值为1，即 =1。还有另一种方法可以得到等式3.9b。用公式 3.8，自由路径s的概率密度函数为：

(3.10)

p（s1）可以代入等式2.5产生式 3.9b，在等式2.5中的积分将恢复到等式3.8中。

在多层混浊介质中，光子包可能会在发生交互作用之前经历多层介质的自由飞行。在这种情况下，等式3.7的对应项变为：

(3.11)

其中，i是层的索引，符号μti是第i层的交互系数，而si是第i层的步长。 总步长s和为：

(3.12)

求和在光子包已行进的所有层上进行。等式 3.11没有考虑光子在边界的反射和透射，因为它们是分开处理的。采样方程式通过等式3.11获得至ξ：

(3.13)

如您所见，等式 3.9b只是3.13的特例。可以将采样解释为总的无量纲步长为–ln（ξ），其中无量纲步长定义为量纲步长si和相互作用系数μti的乘积。由于玻璃层的吸收系数和散射系数为零，因此对方程式3.13的左侧无贡献。等式3.13的详细过程将在3.6和3.7节中讨论。虽然等式3.13看起来很复杂，但其会使3.6和3.7节中的过程的理论基础看起来更简单。

从现在开始，为简单起见，我们将使用步长s而不是s1或ssum。注意，这种采样方法涉及对数函数的计算，这很费时。这在第5.7节中有所反映。可以使用更快的方法来避免对数计算（Ahrens等，1972； Marsaglia，1961； MacLaren等，1964）。

一旦指定了步长s，光子就可以在组织中移动了。光子包的位置通过以下方式更新：

(3.14)

箭头指示数量替换。左侧的变量具有新值，而右侧的变量具有旧值。在实际的c程序中，等号用于此目的。方程式3.14的简单性是使用笛卡尔坐标的主要原因。

一旦光子迈出了一步，就必须计算由于相互作用位点的吸收而导致的光子重量的某种衰减。光子当前重量W的一小部分将沉积在局部网格元素中。计算出的沉积光子重量的量ΔW：

(3.15)

通过添加ΔW来更新沉积在该局部网格元素中的总累积光子重量A（r，z）：

(3.16)

光子重量也必须通过以下方式更新：

(3.17)

具有新权重W的光子数据包将在交互位置处发生散射（稍后讨论）。请注意，整个光子包在步骤结束时会发生相互作用，即吸收或散射。

一旦光子包被移动并且其重量减小，光子包就准备好被散布了。将有一个偏转角θ∈[0，π）和一个方位角ψ∈[0，2π）进行统计采样。 Henyey和Greenstein（1941）最初为银河散射提出的散射函数（请注意，我们在此处定义的散射函数是cosθ的概率密度函数。它与van de Hulst（1980）定义的相位函数的差为常数1/2。）描述了偏转角cosθ的余弦的概率分布：

(3.18)

其中，各向异性g等于，其值在–1到1之间。值0表示各向同性散射，值接近1表示正向散射。雅克等（1987年）实验确定Henyey-Greenstein函数很好地描述了组织中的单个散射。对于组织，g的值范围在0.3到0.98之间，但在可见光谱中，g经常为?0.9。应用公式2.5，cosθ的选择可以表示为随机数ξ的函数：

(3.19)

接下来，采样在0到2π间隔内均匀分布的方位角ψ：

(3.20)

一旦选择了偏转角和方位角，就可以计算出光子包的新方向：

(3.21)

如果光子包的角度与组织表面的法线太近（例如| μz |> 0.99999），则应使用以下公式：

(3.22)

其中当μz为正时，SIGN（μz）返回1；当μz为负时，SIGN（μz）返回–1。最后，更新当前的光子方向：μx = μ'x，μy = μ'y，μz = μ'z。

在两个角度θ和ψ的采样以及方向余弦的更新中，涉及三角函数。由于三角运算需要大量计算，因此我们应尽可能避免使用它们。采样的详细过程可以在文件“ mcmlgo.c”中编写的Spin（）函数中找到（请参阅附录A和附录B.4）。

在一个步骤中，光子包可能撞击组织的边界，该边界位于组织和周围介质之间，其中步长s由等式3.9b计算。例如，光子包可能试图在空气/组织界面处逸出组织。如果是这种情况，则光子包可能以观察到的反射率逃逸（如果还包括后边界，则为透射率）或被边界内部反射。当步长足够大以达到边界时，可以使用不同的方法来解决此问题。让我们首先介绍程序mcml中使用的两种方法之一。

首先，计算缩短的步长s1：

(3.23)

其中z0和z1是当前层的上下边界的z坐标（笛卡尔坐标系见图1.1）。缩短的步长s1是当前光子位置与光子传播方向上的边界之间的距离。当μz为零时，由于光子方向与边界平行，因此光子不会撞击边界。因此，等式3.23不包括μz为零的情况。我们将光子包s1移到边界，并且没有与组织的相互作用（关于移动光子包，请参见第3.3节）。下一步要采取的剩余步长将更新为s←s–s1。如果被内部反射，则光子包将行进剩余的步长。

其次，我们计算光子包被内部反射的概率，这取决于入射角αi到边界上，其中αi = 0表示正交入射。 αi的值计算如下：

(3.24)

斯涅尔定律表明入射角αi，透射角αt和光子从ni入射并传输到nt的介质的折射率之间的关系：

(3.25)

内部反射率R（αi）是根据菲涅耳公式计算的（Born等，1986； Hecht，1987）：

(3.26)

这是两个正交偏振方向的平均反射率。

第三，我们通过产生一个随机数ξ并将该随机数与内部反射率进行比较来确定光子是否被内部反射：

如果ξ≤R（αi），则光子被内部反射；

如果ξ> R（αi），则光子逸出组织 （3.27）

如果光子被内部反射，则光子包会停留在表面上，必须通过反转z分量来更新其方向余弦（μx，μy，μz）：

(3.28)

此时，必须再次检查剩余步长。如果足够大，可以到达另一个边界，则应重复上述过程。如果碰到了组织/组织界面，我们将必须根据以下部分进行处理。否则，如果步长足够小以适合该组织层，则光子包将以较小的步长移动。 在此小步骤结束时，将相应地处理吸收和散射。

另一方面，如果光子包逸出组织，则必须增加特定网格元素（r，αt）的反射率或透射率。反射率Rd（r，αt）或透射率Tt（r，αt）由逸出的光子重量W的大小更新：

(3.29)

由于光子已完全逸出，因此此光子包的跟踪在此结束。可以将新的光子发射到组织中，然后追踪。注意，在我们的仿真中，未散射的透射率（如果有）和漫透射率都被记入Tt（r，αt），没有区别，尽管它们可以像我们在3.1节中讨论的那样加以区别。

值得一提的是对内部反射率建模的另一种方法。并非使光子包的内部反射成为全有或全无的事件，而是每次光子包到达表面边界时都可以使用部分反射方法。 当前光子重量的1–R（αi）分数成功地逃脱了组织，并增加了局部反射率或透射率阵列，例如Rd（r，αt）←R d（r，αt）+ W（ 1–R（αi））。 所有其余的光子权重将被反映，并且光子权重将更新为W←WR（αi）。

两种方法都可以在mcml程序中使用，并且用户可以选择使用哪种方法，具体取决于他们想要得分的物理量。可以更改程序中的标志以在这两种方法之间切换（请参见5.2节）。All-or none方法更快，但部分反射方法应能够减少反射率或透射率的变化。还没有研究通过使用部分反射方法可以减少多少差异。

与第3.5节类似，以等式表示的三角运算次数。为了计算速度，应最小化3.24-3.26。这些计算的详细过程可以在写在文件“ mcmlgo.c”中的RFresnel（）函数中找到（请参阅附录A和附录B.4）。

如果光子步长足够大，足以撞击组织/组织界面，则此步可能会穿过几层组织。考虑一个光子包，它试图在组织1中以μa1，μ s1，n1形成一个大小为s的台阶，但在缩短的步骤s1之后以μa2，μs2，n2到达与组织2的一个交界面。类似于上一节中的讨论，光子包首先被移动到没有交互作用的接口，并且下一步要采取的剩余光子步长被更新为s←s–s1。然后，我们必须通过菲涅尔公式决定统计地反射还是传输光子包。如果反射了光子包，则与上一节中的处理方式相同。但是，如果将光子包传输到组织的下一层，则它必须继续传播而不是被终止。基于等式3.13，必须根据新组织的光学特性转换剩余步长：

(3.30)

其中，μt1和μt2分别是组织1和组织2的相互作用系数。再次检查当前步长s是否有其他边界或界面交叉。重复上述过程，直到步长小到足以容纳一层组织为止。在此小步骤结束时，将相应地处理吸收和散射。

如果光子包在玻璃层中，则光子将移动到玻璃层的边界，而不会更新剩余的光子步长，因为玻璃层中的路径长度不会影响等式3.13的左侧。重要的是要理解，如果光子包穿过多个组织层，则使用等式3.9b的步长和等式3.30的重复使用都基于等式 3.13。

在发射光子包之后，它可以通过反射或透射出组织而自然终止。对于仍在组织内部传播的光子包，如果在许多相互作用步骤后光子重量W已充分减小，使其降至阈值以下（例如Wth = 0.0001），则光子会进一步传播，除非您对光子传播的最后阶段感兴趣，否则产生的信息很少。但是，必须执行适当的端接以确保能量（或光子数量）的守恒，而不会歪曲光子沉积的分布。当W≤Wth时，使用一种称为轮盘赌的技术来终止光子包。轮盘赌技术以mW的重量给光子包以m（例如m = 10）生存的机会。如果光子数据包在轮盘赌中无法幸存，则光子权重将减小为零，并终止光子。

(3.31)

其中ξ是均匀分布的伪随机数（请参见第2章）。该方法节省能量，但以无偏方式终止光子。与轮盘相反，光子轮盘和分裂的组合可以适当地用于减少方差（Hendricks等，1985）。

如前所述，我们在蒙特卡洛模拟过程中记录了光子的反射率，透射率和吸收率。在本章中，我们将详细讨论这些物理量的过程。某些公式的大小显示在公式末尾的方括号中。

z和r方向上的最后一个单元格需要特别注意。因为光子可以传播到网格系统之外，所以当将光子权重记录到漫反射率或透射率阵列或吸收阵列中时，物理位置可能不适合网格系统。在这种情况下，溢出方向上的最后一个单元用于收集光子重量。因此，z和r方向上的最后一个像元在相应位置没有给出实际值。但是，角度α总是在我们为其选择的范围内，即0≤α≤π/ 2，因此在漫反射率和透射率的评分上不会引起问题。

启动光子包后，立即计算镜面反射率。镜面反射后的光子重量被传输到组织中。在模拟过程中，一些光子包可能会离开介质，并且根据光子包的出处将其权重计入漫反射阵列或透射率阵列中。跟踪多个光子包（N）之后，我们分别得到两个计分阵列Rd（r，α）和Tt（r，α）的漫反射率和透射率。它们在程序中分别由Rd-rα[ir，iα]和Tt-rα[ir，iα]内部表示。网格元素中心的坐标通过以下方式计算：

(4.1)

(4.2)

其中ir和iα是r和α的指数。原始数据给出了二维网格系统中每个网格元素中的总光子权重。为了获得二维网格系统每个方向的网格元素中的总光子权重，我们将另一维的二维数组求和：

(4.3)

(4.4)

(4.5)

(4.6)

为了获得总漫反射率和透射率，我们再次对一维阵列求和：

(4.7)

(4.8)

所有这些阵列基于权重为N的N个初始光子数据包，给出每个栅格元素的总光子权重。为了将Rd-rα[ir，iα]和Tt-rα[ir，iα]转换为与每个立体角垂直于光子方向的单位面积的光子概率，将它们除以环形区域在a上的投影。垂直于光子出射方向的平面（Δacosα），围绕环形环的沿α方向的网格线间隔所跨越的立体角（ΔΩ）以及光子包总数（N）：

(4.9)

(4.10)

其中：

(4.11)

(4.12)

其中r和α是根据等式计算的。 4.1和等式 4.2径向分辨的漫反射率Rd-r [ir]和总透射率Tt-r [ir]除以圆环的面积（Δa）和光子包的总数（N），将它们转换为每单位面积的光子概率：

(4.13)

(4.14)

将角度分解的漫反射率Rd-α[iα]和总透射率Tt-α[iα]除以立体角（ΔΩ）和光子包总数（N），将它们转换为每光子概率单位立体角：

(4.15)

(4.16)

将总漫反射率和透射率除以光子包的总数（N），得出概率：

(4.17)

(4.18)

其中[–]表示无量纲单位。

在模拟过程中，将吸收的光子重量计入吸收阵列A（r，z）。 A（r，z）在内部由2D数组Arz [ir，iz]表示，其中ir和iz是r和z方向上网格元素的索引。网格元素中心的坐标可以通过等式4.1及下式计算：

(4.19)

原始数据Arz [ir，iz]给出二维网格系统中每个网格元素中的总光子权重。为了获得每个网格元素在z方向上的总光子权重，我们对r方向上的2D数组求和：

(4.20)

可以从Az [iz]计算每个层A l [layer]中吸收的总光子重量和组织A中吸收的总光子重量：

(4.21)

(4.22)

其中求和范围“层中的iz”包括导致该层中z坐标的所有iz。然后，将这些量适当缩放以得到密度：

(4.23)

(4.24)

(4.25)

(4.26)

数量A给出光子被组织吸收的概率。一维阵列A l [层]给出每一层中光子的吸收概率。此时，Arz [ir，iz]给出吸收的光子概率密度（cm-3），并且可以通过将其除上当前位置所在层的吸收系数μa（cm-1）转换为光子注量φrz（cm-2）：

(4.27)

1D阵列Az [iz]给出z方向上每单位长度的光子概率（cm-1）。也可以将其除以吸收系数μa（cm-1）以得出无量纲的量φz [iz]：

(4.28)

乍一看，这个数量似乎很难理解或多余。但是，等式4.20中原始数据的总和等效于等式7.15中无限宽平面光束的卷积，这将在第7章中讨论。等效性如下所示。 根据等式4.20、4.23和4.24中，最终转换后Az [iz]和Arz [ir，iz]具有以下关系：

(4.29)

其中Δa（ir）在等式4.11中计算，但我们强调它是方程4.29中ir的函数。用等式4.27和4.28，4.29可以转换为：

(4.30)

这是以下积分的数值解：

(4.31)

等式 4.31本质上是等式7.15对于具有常数S差的无限宽扁光束，其中S是无限宽扁光束的功率密度。等式4.31和等式7.15之间的等价关系将在φz（z）替换为F（r，z），将φrz（r，z）替换为G（r''，z），并用r替换r''之后可以看到。因此，φz [iz]给出了一个无限宽的平坦光束的通量，其恒定差就是功率密度S。

程序mcml将仅报告Arz [ir，iz]和Az [iz]而不是φrz [ir，iz]和φz [iz]。将设置程序conv以将Arz [i r，i z]和Az [iz]转换为φrz [ir，iz]和φz [iz]。

在我们的蒙特卡洛模拟中，我们总是建立一个网格系统。计算结果将受到有限网格大小的限制。本节将讨论什么网格大小是最好的。

每个网格元素的平均值位置

该模拟提供每个网格元素中计分的物理量的平均值。现在，问题是应该在哪个位置分配平均值？可以说没有最佳点，因为对物理量的确切答案是未知的。但是，在每个网格元素中物理量的线性近似下，我们可以找到每个网格点的最佳点。在大多数情况下，对于较小的网格尺寸，可以证明线性近似是合理的，因为高阶项远小于线性项。随后将讨论一些特殊场合。

让我们首先讨论r方向上的网格系统，因为r是一个变量，在该变量上将实现有限大小的光子束的卷积（请参见第7章）。r方向上的栅格系统使用Δr作为栅格间隔，共有N个栅格元素。每个网格元素的索引由n表示。每个栅格元素的中心由rn表示：

(4.32)

正如我们已经提到的，蒙特卡洛模拟实际上近似于每个网格元素中物理量Y（r）的平均值，其中Y（r）可以是漫反射率，漫透射率和内部能量密度Arz（r，z）以获得z值。

(4.33)

其中2πrnΔr是当n = 0时环或圆的面积。

如果每个网格元素中的Y（r）线性近似，我们可以证明存在满足的最佳点rb：

(4.34)

其中

(4.35)

证明：Y（r）由关于rb的泰勒级数展开为近似值：

(4.36)

将式4.36代入式4.33得：

即：

(4.37)

如果我们令等式4.37中的方括号为零，求出rb，得到：

(4.38a)

所以等式4.37就变成了：

(4.39)

等式4.38a便可被重新定义为：

(4.38b)

可以将式4.32带入式4.38b的第二项：

(4.38c)

当n=0时，

当n=1时，

当n=2时，

当n=3时，

当n=4时，

…

可以看出，最佳点偏离每个网格元素的中心，并且指向网格框的索引越小，偏差越大。随着索引n变大，最佳点接近网格元素的中心。此现象归因于等式4.33中的2πr因子。对于z方向不存在类似的因子，并且z方向的最佳点应该是每个网格元素的中心。

使用选定网格系统的蒙特卡洛模拟的计算结果始终具有有限的精度，而该精度从根本上受到有限的网格大小的限制。上面的定理仅给出了用模拟结果最能代表函数值的点。

第一次光子与组织相互作用的影响

上面的定理假设Y（r）的可微性，其中Y（r）可以表示特定位置的漫反射率Rd（r），漫透射率Td（r）和在特定z值下内部能量通量φrz（r，z）。该假设应适用于漫反射率和漫透射率。但是，对于内部能量密度，z轴能量密度是脉冲响应（响应无限窄的光子束）的增量函数。因此，您不能假定函数φrz（r，z）在r等于零时，即φrz（r = 0，z）的可微性。Gardner（1992b）等人提供了解决此问题的最佳方法。它们跟踪与介质的第一个光子相互作用，该相互作用始终在z轴上，与其余的相互作用分开。因此，函数φrz（r，z）将不包括产生δ函数的第一光子-组织相互作用。该方法除了具有更好的精度外，还确保了在r等于零时的可微φrz（r，z）。

在mcml的当前版本（1.0版）中，尽管我们打算在以后的版本中添加此方法，但Gardner等人的方法尚未实现。但是，mcml的结果在r方向的网格大小的空间分辨率内仍然是正确的。尽管没有人会使用蒙特卡洛模拟仅吸收半无限介质中的响应，但让我们使用这个简单的示例来说明我们的意思。在这种情况下，Gardner等人的方法将在z轴上产生指数衰减的响应，这是r的增量函数。mcml的1.0版将在第一个网格元素中产生相同的指数衰减响应，这不是增量函数，因为网格元素的体积有限。但是，由于使栅格间距Δr足够小，结果接近Gardner等人的方法。

到目前为止，讨论与有限大小的光子束的卷积无关，这将在第7章中讨论。据加德纳等（1992b），对于半径至少比网格间距Δr大三倍的高斯光束，第一相互作用引起的卷积误差将很小。正如我们将在7.4节中讨论的那样，在卷积程序conv中，扩展的梯形积分而不是在原始网格点上积分，应该进一步减少误差。当然，当分别考虑z轴上的相互作用时，我们的卷积程序conv不受高斯光束半径的限制。相反，如果通过在原始网格点上计算被积分数来近似积分，则光子束的半径仍应足够大以产生合理的积分精度。

由于该定理较晚才被发展起来，因此尚未在程序（mcml和conv）中实现。 我们计划在下一个版本中实现它（请参阅附录H）。

模拟是用ANSI标准C编写的（Plauger等，1989； Kelley等，1990），这使得可以在支持ANSI标准C的任何计算机上执行该程序。到目前为止，该程序已经成功进行了测试在Macintosh，IBM PC兼容，Sun SPARCstation 2，IBM RISC / 6000 POWERstation 320和Silicon Graphics IRIS工作站上运行。本章主要介绍几个规则，程序中使用的重要常量和数据结构，跟踪光子数据包的算法，程序流以及程序的时序。假定具有C的先验知识即可完全理解本章。附录A中详细列出了该程序的流程。附录B逐个文件列出了完整的源代码。在附录C中，我们提供了由UNIX计算机上的make实用程序使用的make文件。附录D和E分别是输入数据文件的模板和样本输出数据文件。附录F为UNIX用户提供了几种C Shell脚本。获取程序所需的信息在附录G中有详细说明。

在编写此程序时，我们遵循以下规则和约定：

符合ANSI标准C，以便可以在各种计算机上执行该程序。

尽可能避免使用全局变量。到目前为止，该程序尚未定义任何全局变量。

避免对程序进行严格限制。例如，我们在运行时通过动态分配消除了对数组元素数量的限制。这意味着只要内存允许，程序就可以接受任意数量的图层或网格线。

预处理程序名称均为大写字母，例如：

#define PREPROCESSORS

对于全局变量，函数名称或数据类型，每个单词的第一个字母为大写，并且单词之间的连接不带下划线，例如：

short GlobalVar;

对于伪变量，每个单词的第一个字母为大写字母，单词之间通过下划线相连，例如：

void NiceFunction(char Dummy_Var);

局部变量都是小写字母，单词之间用下划线连接，例如：

short local_var;

在头文件mcml.h或源文件mcmlmain.c和mcmlgo.c中定义了几个重要的常量。它们在表5.1中列出，并且可以根据您的特殊需要进行更改。程序中使用常量STRLEN定义字符串长度。WEIGHT是阈值光子重量，低于该阈值光子数据包将通过轮盘赌。重量为W的光子小包有机会以新的重量W / CHANCE生存。 COSZERO和COS90D用于使用菲涅耳公式计算反射，并用于计算新的光子方向余弦。|cosα| > COSZERO，α被认为非常接近0或180o。|cosα|

当THINKCPROFILER为1时，可以使用Macintosh上的THINK C编译器的分析器（Symantec，1991年）来监视有关每个功能的程序时序配置文件（请参见5.7节）。当GNUCC为1时，可以使用自由软件基金会（Free Software Foundation）的GNU C编译器（gcc）编译代码，该代码虽然声称符合ANSI C标准，但不支持difftime（）之类的多个功能。因此，当GNUCC为1时，程序中的计时模块将无法工作，尽管否则程序将正常运行。当STANDARDTEST为1时，随机数生成器在喂入固定种子后将生成固定的随机数序列。此功能用于调试程序。通常应将STANDARDTEST设置为0，这会使随机数生成器使用当前时间作为种子。当PARTIALREFLECTION设置为0时，将使用第3.6节中描述的边界（例如，空气/组织边界）处的光子内反射的全有或无模拟机制，其中光子包是全反射或全透射的，取决于菲涅耳反射率与随机数的比较。否则，当PARTIALREFLECTION设置为1时，光子包将被部分反射和透射。通常，我们将PARTIALREFLECTION设置为0，因为全或无的模拟速度更快。

表5.1 程序mcml中的重要常量

常数

文件

数值

含义

WEGHT

mcml.h

1*10-4

阈值weight

CHANCE

mcml.h

0.1

轮盘幸存的机会

STRLEN

mcml.h

256

字符串长度

COSZERO

mcmlgo.c

1-1*10-12

Cosine of ~ 0

COS90D

mcmlgo.c

1*10-6

Cosine of ~ 90

THINKCPROFILER

mcmlmain.c

1/0

在Macintosh上切换THINK C探查器

GNUCC

mcmlmain.c

1/0

GNU C编译器的开关

STANDARDTEST

mcmlgo.c

1/0

切换固定的随机数序列

PARTIAL-

REFLECTION

mcmlgo.c

1/0

边界处部分内部反射的开关

数据结构是程序的重要组成部分。相关参数由结构按逻辑组织，因此程序更易于编写，读取，维护和修改。光子包的参数被分组为单个结构，该结构由以下各项定义：

光子包的位置和行进方向分别由笛卡尔坐标x，y，z和方向余弦（ux，uy，uz）描述。光子包的当前重量由结构构件w表示。成员dead在光子数据包启动时初始化为0，表示光子数据包的状态。如果光子包的重量低于阈值重量，则光子包已退出组织或无法幸免于轮盘赌，则将死成员设置为1。这是用来指示该程序停止跟踪当前的光子数据包。请注意，布尔类型不是ANSI标准C中的内部数据类型。在我们的头文件mcml.h中，它定义为char类型。成员层是光子包所在层的索引。成员类型使用short类型，其值的范围在–32768（2 15）和+32768之间（Plauger等，1989）。尽管可以始终根据光子数据包的笛卡尔坐标和介质的几何结构来识别该层，但是为了提高计算效率而定义了该层。仅当光子包穿过组织/组织界面时才更新成员层。成员s是当前步长的步长，以厘米为单位。当步长足够大以触及边界或界面时，成员sleft用于以无量纲单位存储未完成的步长。例如，如果选择的步长s足够长（如第3.6节所述）以交互系数μt到达当前层的边界，则选择缩短的步长s1作为当前步长，而未完成的步长必须存储。在程序中，实现了以下分配：成员s = s1和成员sleft =（s–s 1）μt。请注意，我们以无量纲单位存储未完成的步长，因此，我们仅需知道当前层的相互作用系数，即可在光子包穿过层时将步长转换回以cm为单位。

描述组织层所需的参数分为一个结构：

顶部和底部边界的笛卡尔坐标分别由z0和z1表示。介质层的折射率，吸收系数，散射系数和各向异性因子分别由成员n，mua，mus和g表示。临界角的余弦分别由成员cos_crit0和cos_crit1表示。用该层相对于两个相邻层的相对折射率来计算它们。

所有输入参数均在以下结构中定义。 该结构的大多数成员由用户在仿真之前提供。

数据输出的文件名是out_fname，其格式（out_fformat）对于ASCII可以是A，对于二进制可以是B。当前，仅支持ASCII格式。成员num_photons表示要模拟的光子包的数量。由于可以模拟更多数量的光子数据包，因此将num_photons成员使用long int类型，其值范围在–2,147,483,648（2 31）和+2,147,483,648（Plauger等，1989）之间。阈值权重由成员Wth表示。重量小于Wth的光子包将经历轮盘赌。网格线间距Δz，Δr和Δα分别由成员dz，dr和da表示。网格元素的数目Nz，Nr和Nα分别由nz，nr和na表示。总层数由成员num_layers表示。成员layerspecs是指向LayerStruct结构的指针。该指针可以动态地分配有结构的阵列，其中每个结构代表一个层，或顶部环境介质或底部环境介质（例如空气）。因此，阵列中有（num_photons + 2）个元素，其中元素0和元素（num_photons + 1）分别存储顶部环境介质和底部环境介质的折射率。动态分配将在后面讨论。

所有的输出数据也被组织成一个结构：

成员Rsp是镜面反射率。指针Rd_ra将动态分配，并且等效地使用，就好像它是r和α上的2D数组一样。Rd_ra是第4.1节中讨论的漫反射率Rd（r，α）的内部表示。成员Rd_r和Rd_a分别是r和α上的一维漫反射率分布。成员Rd是总漫反射率Rd。成员A_rz是2D内部光子分布A（r，z）的表示（请参见第4.2节）。成员A_z和A_1分别是相对于z和层的对应的一维内部光子分布。成员A是整个组织吸收光子的概率。透射率的成员类似于反射率的成员，除了未散射透射率和漫透射率之间没有区别。将所有传输的光子重量计入阵列。

在上面定义的结构中，指针用于表示2D或1D数组。这些指针在运行时根据用户的规范动态分配。因此，用户可以使用不同数量的层数或网格元素数，而无需更改程序的源代码。这提供了程序的灵活性和内存利用率。动态分配程序由Press（1988）修改。这里仅介绍一维阵列分配，矩阵分配可在附录B.5-“ mcmlnr.c”中找到。

图6.4. 比较两种边界条件不匹配的半无限介质内部能量密度与深度z的关系。结果来自蒙特卡洛模拟，每个样本使用mcml封装了100万个光子包。z方向上的网格线间距为0.005厘米，网格元素的数量为200。

我们用指数函数拟合两条曲线的平行部分。曲线的阻尼常数对于匹配边界分别为约1.73cm-1和对于未匹配边界为1.74cm-1。匹配边界的阻尼常数的倒数分别为0.578 cm和不匹配边界的阻尼常数为0.575 cm。它们非常接近根据扩散理论预测的穿透深度（0.57 cm，公式6.2）。

我们使用mcml对具有各种光学特性的半无限介质完成了多个Monte Carlo模拟，并将用户时间与散射系数与吸收系数μs/μa和各向异性g的比值进行拟合。边界匹配（相对折射率为1）和边界不匹配的介质的情况（相对折射率不为1）将单独显示。请注意，用户时间是系统分配给程序运行的时间，而不是实时时间（挂钟时间）。在分时系统中，它们不必相同，并且根据系统的状态，可能不会再现同一运行的实时时间。在mcml中，用户时间被报告到输出数据文件，并且实时时间用于预测仿真过程中仿真何时完成。

在开始针对各种光学特性的多个mcml仿真之前，我们知道完成跟踪光子包所需的时间与光子包终止之前所经历的步骤数成正比。根据第3章中描述的光子传播规则，此步数不应取决于散射系数μs和吸收系数μa的绝对值，而应取决于它们的比率。因此，我们将两个参数之一保持恒定（例如，μs = 100 cm-1），而改变另一个参数（例如，μa）。

匹配的边界

对于边界匹配的介质（相对折射率nrel为1），我们完成了10,000个光子数据包的多个mcml模拟，每个模拟了各种吸收系数μa和各向异性系数，同时将散射系数μs固定为100 cm-1。 结果列于表6.3，其中列出了散射系数和吸收系数之间的比率，而不是分别将两个系数分开列出，用户时间转换为每1000个光子包的秒数，而不是10,000个。随后将讨论“预测的用户时间”和“错误”列。

表6.3. 具有匹配边界的各种媒体的计算时间。列“预测的用户时间(Predicted User Time)”是使用等式的计算值。稍后介绍6.1-6.3。“错误(Error)”列等于（预测的用户时间–用户时间）/（用户时间）* 100。

μ s /μ a

g

User Time

(sec./1000 photons)

Predicted User Time

(sec./1000 photons)

Error

(%)

0.2

0

0.43

0.36

–16.15

1

0

0.78

0.81

3.58

2

0

1.10

1.15

4.42

10

0

2.59

2.56

–0.99

20

0

3.66

3.62

–0.97

100

0

8.46

8.10

–4.23

200

0

12.27

11.46

–6.64

1000

0

29.61

25.61

–13.49

0.2

0.1

0.40

0.38

–5.04

1

0.1

0.80

0.86

7.19

2

0.1

1.10

1.22

10.70

10

0.1

2.80

2.75

–1.82

20

0.1

3.90

3.90

0.10

100

0.1

9.80

8.81

–10.07

200

0.1

13.20

12.52

–5.19

1000

0.1

28.60

28.25

–1.21

0.2

0.5

0.50

0.45

–9.54

1

0.5

1.00

1.08

7.71

2

0.5

1.40

1.57

11.80

10

0.5

3.50

3.73

6.49

20

0.5

5.30

5.42

2.19

100

0.5

12.10

12.90

6.60

200

0.5

17.90

18.74

4.71

1000

0.5

40.60

44.63

9.93

0.2

0.9

0.50

0.49

–1.83

1

0.9

1.20

1.35

12.74

2

0.9

1.90

2.09

10.20

10

0.9

6.30

5.77

–8.39

20

0.9

9.80

8.93

–8.86

100

0.9

25.80

24.62

–4.58

200

0.9

38.10

38.10

0.00

1000

0.9

95.30

105.02

10.20

0.2

0.99

0.50

0.42

–16.18

1

0.99

1.20

1.42

18.68

2

0.99

2.10

2.41

14.85

10

0.99

8.80

8.20

–6.87

20

0.99

16.50

13.88

–15.89

100

0.99

60.60

47.16

–22.18

200

0.99

97.20

79.86

–17.84

1000

0.99

248.80

271.37

9.07

我们以对数-对数标度绘制了每个各向异性系数g的散射系数与吸收系数之比的用户时间（图6.5）。对于每个各向异性因子g，散射系数和吸收系数之间的比率越高，使用时间越长。这是因为与比率较低的介质相比，比率较高的μs/μa介质中的光子包在达到阈值权重之前可以跳更多的步，才有机会被终止。因此，低吸收介质的蒙特卡罗模拟非常慢。如果在低吸收率半无限浑浊介质中作为r的函数的漫反射率是唯一要计算的物理量，则纯Monte Carlo模拟和漫射理论的混合模型（Wang等，1992）要快得多。模型比纯Monte Carlo模拟要好。混合模型的速度对散射系数和吸收系数之间的比率不太敏感。

对于散射系数和吸收系数之间相同的比率，各向异性系数g越大，使用时间越长。这是因为各向异性较大的介质g中的光子包有较少的机会从介质中反射出来，因为散射是更向前的，因此被终止了。

对于每个各向异性因子g，数据点在对数-对数图中均对齐良好。这意味着我们可以使用幂函数拟合每个各向异性因子g的数据点。拟合线如图6.5所示。

两个拟合系数C1（g）和C2（g）取决于各向异性因子g。拟合系数C1（g）相对于各向异性因子g以对数线性标度绘制（图6.6a），并且可以拟合指数函数。同样，拟合系数C2（g）相对于（1-g）以线性对数标度绘制（图6.6b），并且可以拟合对数函数。这三个拟合总结为以下经验公式：

(6.3)

(6.4)

(6.5)

其中t是具有匹配边界的半无限媒体的用户时间，以每1000个光子包的秒数为单位。

图6.5. 用户时间与具有匹配边界的介质的不同各向异性因子g的散射系数μs和吸收系数μa之比。

图6.6. 拟合系数（a）C1（g）和（b）C2（g）与边界匹配的介质的各向异性因子g的关系。

要测试均衡器有多好。在图6.3-6.5中，我们使用它们来计算表6.3中给出的光学性能的用户时间，并在表6.3中的“预测用户时间”列中显示了计算出的用户时间。对于表中所有行（其中之一除外）的相对错误，如表6.3的“错误”列所示，在20％以内。

等式6.3-6.5通常可以用来预测边界匹配的介质的用户时间。但是，必须提及一些限制和注意事项。首先，各向异性因子g限制为小于0.99，而g> 0.99的精度是未知的，因为当g接近1时，拟合系数C2（g）接近无穷大。其次，这些方程式基于Sun SPARCstation 2上的mcml仿真。因此，我们期望方程6.5的比例因子或表6.3中的值在其他计算机系统上。可以通过使用计算机系统上运行的mcml模拟表6.3中具有光学特性的一种或多种介质，并采用计算机上的用户时间与表6.3中的用户时间之比来确定此比例因子。第三，mcml的速度与程序mcml中的阈值权重（表5.1中的WEIGHT）相关，通常为1×10 –4（见表5.1），以及生存轮盘的机会（表5.1中的机会），通常为0.1（请参见表5.1）。如果更改了这两个参数，则等式6.3-6.5不再有效。尚未探索这两个参数将如何影响用户时间。

不匹配的边界

我们对边界不匹配的介质（相对折射率nrel=1）重复上述过程。我们完成了1000个光子包（而不是10,000个匹配边界）的多个mcml模拟，每个光子包具有各种吸收系数和各向异性因子，同时保持散射系数固定为100 cm-1和相对折射率为1.37，这在可见或红外波长对于人体组织是典型的参数。结果列于表6.4。

我们以对数-对数刻度绘制了用户时间与每个各向异性因子g的散射系数与吸收系数之比的函数关系（图6.7）。然后，分别针对g和（1-g）绘制拟合系数C1（g）和C2（g）（图6.8a和b），并分别拟合指数函数和对数函数。拟合给出以下经验公式：

(6.6)

(6.7)

(6.8)

其中t是边界不匹配（nrel = 1.37）的半无限媒体的用户时间，以每1000个光子包的秒数表示。

表6.4. 边界不匹配的各种媒体的计算时间。列“预测的用户时间”是使用等式的计算值。稍后介绍6.4-6.6。 “错误”列等于（预测的用户时间–用户时间）/（用户时间）* 100。

μ s /μ a

g

User Time

(sec./1000 photons)

Predicted User Time

(sec./1000 photons)

Error

(%)

0.2

0

0.43

0.43

–9.74

1

0

0.98

1.05

7.14

2

0

1.42

1.54

8.26

10

0

3.67

3.73

1.51

20

0

5.47

5.45

–0.28

100

0

13.43

13.22

–1.57

200

0

18.67

19.35

3.66

1000

0

50.50

46.90

–7.13

0.2

0.1

0.50

0.44

–11.03

1

0.1

1.00

1.09

8.85

2

0.1

1.50

1.60

6.68

10

0.1

3.83

3.92

2.23

20

0.1

5.90

5.76

–2.44

100

0.1

15.95

14.08

–11.70

200

0.1

22.55

20.71

–8.18

1000

0.1

42.38

50.66

19.54

0.2

0.5

0.53

0.49

–8.10

1

0.5

1.13

1.26

11.25

2

0.5

1.73

1.89

9.31

10

0.5

4.87

4.88

0.22

20

0.5

7.63

7.34

–3.77

100

0.5

19.53

18.95

–2.97

200

0.5

27.08

28.51

5.28

1000

0.5

77.96

73.58

–5.62

0.2

0.9

0.55

0.49

–11.64

1

0.9

1.27

1.45

14.31

2

0.9

2.03

2.33

14.58

10

0.9

7.28

6.95

–4.55

20

0.9

12.10

11.13

–7.99

100

0.9

35.98

33.26

–7.56

200

0.9

55.05

53.28

–3.21

1000

0.9

133.06

159.18

19.63

0.2

0.99

0.55

0.41

–25.96

1

0.99

1.28

1.50

17.16

2

0.99

2.10

2.63

25.19

10

0.99

8.82

9.68

9.77

20

0.99

16.82

16.97

0.92

100

0.99

69.55

62.51

–10.12

200

0.99

120.56

109.60

–9.09

1000

0.99

366.45

403.62

10.14

图6.7. 在相对折射率为1.37的介质中，用户时间与散射系数μs与吸收系数μa对不同各向异性因子g的比之比。

图6.8. 相对系数1.37的介质的拟合系数（a）C 1和（b）C 2与各向异性系数g的关系。

要测试等式.6-6.8是，我们使用它们来计算表6.4中光学特性的用户时间，并在表6.4中的“预测用户时间（predicted user time）”列中显示了计算出的用户时间。对于表中除其中两行（黑体字）以外的所有行，相对误差在20％以内。对于边界匹配的媒体的警告也适用于此。

为了总结匹配边界和不匹配边界的经验公式，我们在表6.5中列出了这些公式。

表6.5. 边界匹配和不匹配的用户时间的经验公式。

Items

Matched(nrel=1)

Mismatched(nrel=1.37)

C1(g)

0.81 exp(0.57 g)

1.05 exp(0.36 g)

C2(g)

0.50-0.13log(1-g)

0.55-0.13log(1-g)

t[sec./1000 photons]

C1(g)(us/ua)C2(g)

C1(g)(us/ua)C2(g)

对于多层组织，我们尚未找到基于其他理论的计算结果与我们的计算进行比较。但是，由于Gardner的合作（Gardner等人，1992），我们将多层组织的Monte Carlo模拟结果与独立编写的Monte Carlo模拟结果进行了比较。比较包括漫反射率与半径Rd（r），其中半径r是光子入射点与观察点之间的距离，透射率与半径Tt（r）以及内部能量密度与z的关系， φz（z）以及内部能量密度与r和z的关系，φrz（r，z）。这种比较至少可以大大减少编程错误的机会。

我们选择了三层组织进行仿真。各层的光学性质示于表6.6。

表6.6. 三层组织的光学性质。

Layer

Refractive

Index n

Absorption

Coeff. (cm–1 )

Scattering

Coeff. (cm–1 )

Anisotropy

Factor g

Thickness

(cm)

1

1.37

1

100

0.9

0.1

2

1.37

1

10

0

0.1

3

1.37

2

10

0.7

0.2

顶部和底部环境介质的折射率都设置为1.0。z和r方向的网格间距均为0.01 cm。在z和r方向上的网格元素数分别为40和50。我们不想解析反射或透射光子的出射角，因此，我们将角度α方向上的网格元素数设置为1。跟踪的光子包数是1000000。程序mcml的实际输入文件如下所示。

克雷格·加德纳（Craig Gardner）为他的模拟代码设置了相同的运行，只是他只使用了100,000个光子（Gardner等，1992）。表6.7比较了这两个模拟的总漫反射率和透射率。

表6.7. 总漫反射率和透射率的比较。

Source

Diffuse Reflectance Rd

Transmittance Tt

Gardner et al., 1992

0.2381

0.0974

Mcml

0.2375

0.0965

漫反射率和透射率之间的比较示于图2和3中。分别为6.9和6.10。对于一些给定的z坐标，作为脉冲响应的注量φrz（r，z）之间的比较是半径r的函数，如图6.11所示。如第4.2节所述，注量φz（z）作为z的函数之间的比较随后在图6.12中给出。所有比较都表明两个独立模拟之间的一致性。

图6.9. 根据Gardner的计算（Gardner等，1992）和mcml模拟，比较漫反射率与半径的函数关系。

图6.10. 根据Gardner的计算（Gardner等人，1992）和mcml模拟，比较透射率与半径的关系。

图6.11. 基于Gardner的计算（Gardner等人，1992年）和mcml仿真，对多个z坐标的注量与半径r的函数之间的比较。

图6.12. 基于Gardner的计算（Gardner等人，1992）和mcml仿真，作为z函数的注量之间的比较。

对于2D数组（例如，作为r和z的能量密度的函数），conv能够以轮廓格式输出数据（请参见10.9和10.10节）。与脉冲响应的r和z有关的能量密度在图6.13中的轮廓线中显示。

图6.13. 基于mcml仿真的注量与脉冲响应的r和z的关系等高线图。

本章将讨论无限窄光子束的卷积蒙特卡洛模拟结果的原理和实现，以产生对有限尺寸的光子束的响应。高斯光束和圆形平面光束被认为是特例。

到目前为止，我们仅处理了对通常入射在多层组织表面上的无限窄光子束的响应。该响应也称为脉冲响??应。但是，实际上所有光子束的大小都是有限的。从理论上讲，我们可以通过分布发射光子数据包的初始位置，直接使用蒙特卡洛模拟来计算对有限尺寸光子束的响应。唯一的问题是，与模拟无限窄的光子束的响应相比，它需要跟踪更多数量的光子包以获得可接受的方差。因此，尽管有时对于某些类型的无法卷积的组织结构（例如带有不规则掩埋物体的组织），它可能是唯一的方法，但该方法并不高效。

幸运的是，我们要处理的系统是线性且不变的。线性意味着如果无限窄的光子束的输入强度乘以一个因子，则响应将乘以相同的因子。这也意味着对两个光子束的响应是对每个光子束的响应之和。不变性是指，当无限窄的光子束沿某个方向水平移动一定距离时，响应也将沿相同方向水平移动相同距离。因此，如果我们假设有限尺寸的光子束是准直的，那么无限窄的光子束的响应将是组织系统的格林函数，并且可以通过格林的卷积来计算有限尺寸的光子束的响应。根据有限尺寸的光子束的轮廓起作用。

注意，上述响应可以是内部吸收分布，也可以是反射率或透射率分布。尽管它可能不是所有三个坐标的三变量函数，但我们通常将响应表示为C（x，y，z）。随后将讨论由于对称性引起的这种简并性。我们将与所考虑的响应类型相对应的格林函数表示为G(x, y, z)。由于光子束通常入射在组织表面上，因此函数G（x，y，z）具有圆柱对称性。如果准直的光子束作为光源具有强度分布S（x，y），则可以通过卷积获得响应（Prahl，1988； Prahl等，1989）：

(7.1a)

或通过x''= x–x'和y''= y–y'进行变量转换：

(7.1b)

在等式中7.1a格林函数是源点（x'，y'）与观察点（x，y）之间的距离的函数，其中距离为：

(7.2)

如果光源的强度轮廓S（x'，y'）也具有圆柱对称性，则（x'，y'）只是相对于坐标系原点的源点（x'，y'）半径的函数，半径为：

(7.3)

因此，等式7.1a可以考虑这些对称性重新制定：

(7.4a)

类似的，等式7.1b也可以用以下对称性来重新表示：

(7.4b)

由于响应C（x，y，z）将具有相同的圆柱对称性，因此可以在圆柱坐标系中更容易地处理该问题，其中7.4a和7.4b可以写成：

(7.5a)

(7.5a)

等式7.5b比7.5a在计算中更具优势，因为θ''上的积分与z无关。因此，对于所有深度z只需计算一次该积分。在随后出现的某些情况下，可以解析地表示θ''上的积分，因此将二维积分转换为一维积分。

变量的转换如图7.1所示。第一坐标系（r'，θ'，z）的原点位于源的中心，如图7.1A和等式7.5a所示。观察点在（r，0）。源的增量区域位于（r'，θ'），其值为S（r'）。源点与观察点之间的距离为d'等于 。

第二个坐标系（r’’，θ’’，z）的原点位于观察点的中心，如图7.1B和等式7.5b所示。源位于（r，0）的中心。源的增量区域位于（r’’，θ’’）处，并且有值S（d’’），其中d’’等于 。源点与观测点之间的距离是r’’。

在以下各节中，我们将高斯光束和圆形平面光束作为圆柱对称的示例，以进一步简化公式7.5b。

图7.1。 两个坐标系之间的变换的图示。点划线的圆是激光源，P是观察点。圆线示意性地表示式7.5a（A）和7.5（B）积分。

在高斯光束的情况下，如果忽略了发散，则可以应用上述卷积。如果高斯光束的1/e2半径用R表示，则光束强度分布为：

(7.6)

其中中心强度（r = 0）S0与总功率P的关系为：

(7.7)

将式7.6带入式7.5b，卷积变成了：

(7.8)

方括号中的积分类似于修改后的Bessel函数的积分表示（Spiegel，1968）：

(7.9a)

可以重新格式化为：

(7.9b)

可以将7.9b带入7.8将其重写：

' (7.10)

其中I0是零阶修改的Bessel函数。

如果光子束在半径R内是均匀的并且已准直，则源函数变为：

(7.11)

其中P是光束的总功率。将式7.11代入式7.5b，卷积变为：

(7.12)

函数Iθ（r，r''）为：

(7.13)

根据等式 7.13和图7.1B，式7.12中整合的极限可以被改为一个有限区域：

(7.14a)

其中

(7.14b)

函数Max接受两个参数中较大的一个。

图7.1B中的圆线说明了等式7.13中的第二种情况，观察点P在源外。当点P在源外时，式7.13的第一种情况永远不会满足。对于点P在源内的情况，读者可以类似地绘制图片。请注意，积分的大部分区域都位于源的外部，因此被积分物中S的值为零。

作为圆形平面光束的一种特殊情况，我们让半径R接近无穷大，这表示无限宽的平面光束。在这种情况下，总功率P也接近无穷大，但是我们可以使用功率密度来描述光束的强度。这种情况下的卷积可以通过简单地让R→∞和P /（πR2）→S来完成，其中S是等式7.14中的功率密度或辐照度（W / cm2）。Iθ（r，r''）恒定为1。等式 7.14变为：

(7.15)

在这两种特殊的光子束情况下，二维积分被转换为一维积分，这在数值计算中明显更快。由于蒙特卡洛模拟将物理量评分到离散的网格点，因此积分算法的最佳选择是扩展梯形法则，该法则由Press等人（1988）用C编写，称为qtrap（）。 “复杂性的提高通常会转化为一种更高阶的方法，其效率仅对于足够光滑的被积物才有效。qtrap是一种选择方法，例如，对于被积物来说是函数的函数，该变量在测量数据点之间线性插值。”

积分的最简单选择是将原始网格点处的被积分值相加，再乘以网格间距。但是，这种方法对集成精度没有任何控制。对于给定的精度，有时此方法提供的精度超出了要求，这浪费了计算时间，而有时却给出的精度更低，这不符合预期。例如，r方向上的原始网格元素的数量为50，我们希望将响应卷积在半径为R的圆形平面光束上，该半径约为5Δr，其中Δr是r方向上的网格间距。要在等式中计算C（0，z）。在7.14a中，从0到R的积分范围仅覆盖5Δr。这意味着将仅完成5个功能评估，这可能会产生无法接受的答案。相反，扩展梯形规则会进行正确的计算，直到达到用户指定的精度为止。

我们对原始函数qtrap（）进行了一些修改，以便将所需的精度作为参数。因此，程序conv的用户可以在运行时更改允许的错误（该程序conv将在下一节中讨论）。

扩展梯形积分中使用的积分值评估顺序如图7.2所示。（Press等，1988）如果我们将f（x）积分到[a，b]上，则在第一步中评估f（a）和f（b），如图7.2中的1和2所示。除非函数为线性，否则此步骤将不会提供足够的精度。为了优化网格，我们在第二步中评估f（（a + b）/ 2），如3所示。我们继续执行此过程，直到积分评估达到指定的精度为止。

请注意，图7.2中第三次评估之后的被积评估顺序类似于一个完美平衡的二叉树。如果我们要存储评估后的函数值，自然将它们存储在二叉树中以进行速度检索（将在后面进行讨论）。

图7.2. 由qtrap（）调用的trapzd（）中的被整数求值序列的图示。 “对例程trapzd（）的顺序调用合并了先前调用的信息，并仅在精炼网格所需的那些新点上评估被积数。最下面的一行显示了第四次调用后的函数求值的总数。”（Press et al，1988）第三次评估后，被积物评估的顺序类似于一个完美平衡的二叉树。

物理量的内插和外推

如图7.2所示，函数qtrap（）将需要计算被积分数，因此，物理量可能不在原始网格点上。线性插值用于介于两个原始网格点之间的那些点。线性外推用于超出原始网格系统的那些点。内插和外推如图7.3所示。实心圆表示网格点处的原始得分值。 实线和虚线分别表示内插和外推。对于r方向上给定数量的网格元素（例如，此图片中的Nr = 8），外推最多可计算到（Nr-0.5），因为线性外推对于超出（Nr- 0.5）。因此，物理量被设置为零，超过（Nr-0.5）。有时，数据可能非常嘈杂，以至于最后一点的功能甚至比第二点到最后一点的功能还要高。在这种情况下，不使用外推法。相反，我们只需将函数值设置为零即可。请注意，正如我们在第4章开头提到的那样，r方向上的最后一个单元用于收集不适合网格系统的光子。因此，最后一个像元中的值通常比第二个到最后一个像元中的值高得多，因此在卷积过程中不使用。一言以蔽之，物理量在间隔[0，rmax]中不为零，其中rmax为：

(7.16)

其中Δr是r方向上的网格间距。

图7.3. 物理量的内插和外推的图示。例如，将r方向Nr上的网格元素的数量设置为8。Δr是r方向上的网格间隔。积分极限为a和b（请参见方程7.21和7.22）。箭头指向要评估被积物的位置。

高斯光束的积分评估

如方程式7.10所示，物理量的评估只是在高斯光束上进行卷积的被整数评估的一部分。尽管由于物理原因积分必须收敛，但方程式7.10的形式由于修改后的Bessel函数会随着参数的增加而迅速增加，因此可能无法直接进行数值计算，因此它可能会超过计算机可以容纳的限制（例如，对于某些计算机为10 +38）。因此，需要适当的公式来计算公式7.10。我们注意到，在参数较大的区域中，修改后的Bessel函数具有以下近似值：

(7.17)

因此，如果我们从I0（）中提取指数项，则可以确保修改的Bessel函数随着参数的增加而减少。我们基于I0（）定义以下新函数：

(7.18a)

或者

(7.18a)

I0e（）应该始终是有界的。请注意， 给出7.17只是为了显示函数I0（）的渐近行为。等式 7.18a绝不带有任何近似值。将式7.6、7.7和7.18b 带入7.10，其变成：

(7.19)

由于指数项和I0e（）项均减少，因此可以计算被积数而不会超出范围。

到目前为止，我们已经解决了如何计算卷积而没有溢出的问题。但是，计算速度是另一个问题。我们发现等式7.19中对exp（）I0e（）的求值是每个积分计算的主要部分，根据要解决的特定问题，该比例最多可以达到90％。对于多变量物理量（例如A（r，z）），在针对不同的z坐标计算积分的同一时间点时，卷积可以重复计算式7.19中的exp（）I0e（）。

因此，如果我们可以保存函数评估，即等式7.19中exp（）I 0e（）的计算，对于一个z坐标，则可以节省很多计算时间。但是，积分将迭代执行，直到达到给定的精度为止。因此，功能评价的数量事先未知。我们只能使用动态数据分配来保存功能评估。此外，由于梯形积分qtrap（）的求值序列类似于图7.2中的二叉树，因此可以使用平衡良好的二叉树来存储用于搜索速度的函数求值。

高斯光束的积分极限

由于积分极限在等式7.14a中对于圆形平面光束是有限的，可以使用函数qtrap（）直接计算积分。相反，等式7.19中的高斯光束积分上限为无穷大。该问题可以通过变量转换和积分来解决，其可以通过常规的midexp（）来计算（Press等，1988）。但是，我们发现这种方法的计算效率不高。通过适当地将等式7.19中的指数项截断，可以将上限减小为有限值，当下式满足时：

(7.20a)

或者

(7.20b)

其中K是可以在卷积程序conv中设置的常数，我们将计算被积，否则我们认为被积可忽略不计。例如，如果我们选择K等于4（在程序中实际使用），则式7.19中的指数项约为1x10e–14，其数量级远大于计分物理量的数量级动态范围。

正如我们在本节开头所讨论的，我们仅计算区间[0，rmax]中的物理量，其中rmax由等式7.16给出。结合这个极限和等式7.20b，式7.19变为：

(7.21a)

(7.21b)

(7.21c)

其中函数Max（）和Min（）分别采用两个参数中的较大值和较小值。

圆扁梁的综合评估

圆扁光束的积分评估比高斯光束的积分评估简单得多。但是，等式7.14a中的Iθ（）的评估非常耗时。类似于对高斯光束的被积评估，使用二叉树存储评估的Iθ（），以加快积分速度（请参阅对高斯光束的讨论）。

圆形扁梁的积分极限

由于积分极限在等式7.14a中对于圆形平面光束是有限的，可以使用函数qtrap（）直接计算积分。正如我们在本节开头所讨论的，我们仅计算区间[0，rmax]中的物理量，其中rmax由等式7.16给出。考虑到这个限制，7.14a变为：

(7.22a)

(7.22b)

(7.22c)

其中函数Max（）和Min（）分别采用两个参数中的较大值和较小值。

卷积错误的来源

在等式7.21c和7.22c中，积分的上限可能受到rmax的限制，rmax是蒙特卡罗模拟过程中r方向的网格极限。在r方向上超出原始网格限制的物理量不会导致卷积，从而导致错误。

让我们首先讨论圆形平面梁的情况，因为它更容易。从等式7.22c，我们知道，当满足：

(7.23a)

或者

(7.23b)

r方向上的有限网格不会影响卷积。否则，卷积将被有限网格沿r方向截断。下一节将介绍这种效果。因此，我们不应该相信r≥rmax–R的卷积数据。换句话说，如果您想观察r在循环中半径为R的平面光束的物理量，在r方向上的网格限制应足够大，以使等式7.23a用mcml执行蒙特卡洛模拟时，将被保留。

对于高斯光束，没有像公式7.23这样的干净地描述了有效范围，因为高斯光束在理论上沿r方向延伸至无穷远。但是，当r >> R时，半径为1 / e2的高斯光束的卷积结果与半径为R的圆形平面光束的卷积结果非常接近（在下一节中显示）。因此，对于高斯光束，我们可以对圆形扁光束使用相同的标准（公式7.23），以达到一定的精度。

另一个错误源是由mcml 1.0进行的蒙特卡洛模拟。该版本的mcml不会像Gardner等人（1992b）那样单独评分第一个交互（请参见第4.3节）。 如果高斯光束的半径小于栅格间距Δr的三倍，则这可能会产生很大的误差。换句话说，应满足以下方程式以获得可靠的卷积：

(7.24)

总而言之，当等式7.23和7.24成立，卷积应该是可靠的。

mcml的蒙特卡罗模拟结果的卷积过程在另一个名为“ conv”的程序中实现。像程序mcml一样，它是用ANSI标准C编写的，因此可以在各种计算机上执行。

高斯光束的卷积结果

我们没有任何标准数据来验证卷积程序。但是，克雷格·加德纳（Craig Gardner）使用他的卷积程序（Gardner等，1992）提供了一些卷积结果。脉冲响应基于第6.6节中讨论的仿真，其中我们使用了相同的混浊介质和网格系统。他根据他的蒙特卡洛模拟结果计算了卷积，然后在6.6节中比较了蒙特卡洛模拟结果之后，我们对其进行了卷积。

入射光子束是高斯束，其总能量为1 J，半径为0.1 cm。卷积的漫反射率和透射率在图7.4和7.5中进行了比较。卷积的注量在图7.6中进行了比较。需要注意的是在图7.4-6中的曲线在等于0.5厘米的r处弯曲得更快。这可以通过等式7.21中的积分来解释。由于网格系统的空间限制（50个径向网格，间距为0.01 cm，或rmax = 0.5cm），随着观察点r接近rmax，积分上限被rmax越来越多地削减。因此，集成会低估真实价值。

图7.4. 根据conv和Gardner的计算比较作为r的函数的漫反射率（Gardner等，1992）。

图7.5. 根据conv和Gardner的计算，将透射率作为r的函数进行比较（Gardner等，1992）。

图7.6. 基于conv和Gardner的计算，对于给定的z坐标，将注量作为r的函数进行比较（Gardner等，1992）。

对于2D数组（例如，作为r和z的能量密度的函数），conv能够以轮廓格式输出数据（请参见10.9和10.10节）。与高斯光束的r和z有关的能量密度如图7.7轮廓线所示。

图7.7. 基于高斯光束的conv的注量[J cm-2]随r和z变化的轮廓图。高斯光束的总能量为1J，1/e2半径为0.1 cm。蒙特卡洛模拟适用于表6.6的三层组织。吸收系数μa和散射系数μs以cm–1为单位。

圆扁光束的卷积结果

对于圆形平面光子束，我们没有其他结果可比较。但是，我们想将圆形平面光子束的结果与高斯束的结果进行比较。高斯光束的结果来自上述计算。圆形平面光束的总能量为1J，半径为0.1cm。在图7.8、7.9和7.10中分别比较了漫反射率，透射率和注量。

图7.8. 使用conv在高斯光束和平面光束上卷积的r的漫反射之间的比较。两个光束的总能量为1J，半径为0.1cm。

图7.9. 使用conv在高斯光束和平面光束上进行卷积的透射率随r的比较。两个光束的总能量为1J，半径为0.1cm。

图7.10. 使用conv比较在高斯光束和平面光束上z等于0.005cm处的注量与r的关系。两个光束的总能量为1J，半径为0.1cm。

可以看出，当r大于约2R时，高斯光束和平面光束的结果几乎相同，其中R是光束的半径。此外，如上一节中所述，当r接近r方向的网格极限0.5 cm时，两种响应都会向下弯曲。

卷积错误

卷积积分是迭代计算的。当积分的新估计值与旧估计值之间的差异仅占新估计值的一小部分时，迭代将停止。用户可以使用命令“e”来控制该小的比例。取值范围是0到1。较小的值可以提供更好的精度，但需要更长的计算时间，反之亦然。通常，建议使用0.001至0.1。有时，高允许误差会导致卷积结果有些不连续。如果发生这种情况，请选择一个较低的允许卷积误差，然后重做卷积。例如，在图7.10中的高斯光束上进行卷积时，允许的卷积误差为0.001。如果选择允许的卷积误差为0.01，则可以看到注量分布的不连续性（图7.11）。

图7.11. 使用z等于0.005cm时在r处的注量与r的函数比较，使用conv在具有不同允许卷积误差的高斯光束上卷积。允许误差为0.01的结果在r=0.15cm附近具有不连续性。使用0.001的允许误差可消除不连续性。高斯光束的总能量为1J，1/e2半径为0.1 cm。

2. 用户使用手册

略