[发明专利]一种基于集群式GPU加速的荧光蒙特卡罗模拟方法

权利要求书

1.一种基于集群式GPU加速的荧光蒙特卡罗模拟方法，其特征在于包括以下步骤：

(1)将目标生物组织的空间结构分割成一个三维体素模型，设定一个与三维体素模型大小相同的三维数字矩阵，矩阵中每个元素的数值为标识的组织类型；设定目标生物组织的各类组织光学特性参数，组织光学特性参数为激发光的吸收系数、激发光的散射系数、荧光的吸收系数、荧光的散射系数、荧光团的吸收系数、荧光团的折射系数和荧光团的各向异性因子；

(2)通过MPI这种消息传递通信协议，主节点获取各子节点的GPU设备的数量，将需要计算的源的数量平均分配给每个子节点的GPU设备，各子节点确定GPU设备上网格和块的维度和尺寸，并根据该节点GPU设备的数目开辟相应数目的子进程，将计算任务分配给各个子进程，启动并行编程与计算平台CUDA，各子节点CPU分配内存空间、显存空间，将要计算的数据从内存复制到显存上；

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

(4)在目标生物组织内追踪每个激发光光子的传输过程，将其中耗费时间长的大规模数据并行、高计算密度的步骤安排在GPU上并行执行，其余步骤安排在CPU上串行执行；CPU分配内存，用于存放GPU输出数据，将计算后显存上的数据复制到内存上，各子节点返回运行信息给主节点；步骤(4)具体包括以下步骤：

(8.1)投放激发光光子，如果光子处在组织表面之上，将该位置设定为光子发射位置；如果光子处在组织表面之外，将光子自动移向组织表面，其具体操作使用迭代方法；如果光子处在组织内部，则将依照此位置直接发射光子，即光子不发生与表面的相互作用，而直接开始光子在组织内部的传输；

(8.2)计算光子的剩余步长sleft:Sleft＝-Inξ；

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

s＝min(-In(ξ)/μ_t,min(dx,dy,dz))

这里ξ是由计算机伪随机数产生器产生的均匀分布在(0,1)的伪随机数，dx,dy,dz依次为组织模型中每个体素的长、宽、高；当光子为激发光光子时：μ_t＝μ_a+μ_s+μ_afp；当光子为荧光光子时：μ_t＝μ_af+μ_sf；μ_a和μ_s分别是激发光的吸收系数和散射系数；μ_af和μ_sf分别是荧光的吸收系数和散射系数；μ_afp是荧光团的吸收系数；

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

(8.5)确定光子在界面上的撞击点，方法为：首先找出如下表达式中的最小项：

x_t=(1+[x/dx]-x/dx)/uxy_t=(1+[y/dy]-y/dy)/uyz_t=(1+[z/dz]-z/dz)/uz]]>

这里，x,y,z表示光子当前位置；如果第一项最小，那么光子在界面上的作用点表达式为：

x0=x+ux*x_ty0=y+uy*y_tz0=z+uz*z_t]]>

这里，x0,y0,z0是光子在界面的作用点；

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

(a,b,c)=([x]-[x0],[y]-[y0],[z]-[z0])(([x]-[x0])2+([y]-[y0])2+([z]-[z0])2)1/2]]>

然后，根据该矢量表达式计算光子与界面发生作用的入射角θ和折射角θ_t,其计算式为：

cosθ＝a·ux+b·uy+c·uz

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₁)²，设置光子反射，其方向变为：

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

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

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

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

其中，r的表达式为：

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

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

(8.8)如果光子在(8.6)中折射到空气中，即为光子逸出组织，设定光子死亡，停止对光子追踪且更新sleft＝0，并记录该光子当前属性信息，包括位置、方向，转到步骤(8.14)；其他情况下Sleft＝Sleft-s*μ_t；如果该光子为激发光光子则μ_t＝μ_a+μ_s+μ_afp；当光子为荧光光子则μ_t＝μ_af+μ_sf；如果Sleft＞0，则进入步骤(8.3)；否则进入步骤(8.9)；

(8.9)通过设定随机数判定光子是被吸收或散射，如果该光子为荧光光子则转到步骤(8.10)，否则则转到步骤(8.11)；

(8.10)当ξ≤μ_af/(μ_af+μ_afp)时，则荧光光子被吸收，设定光子死亡，停止对光子追踪，转到步骤(8.14)；如果ξ＞μ_af/(μ_af+μ_afp)时，则荧光光子散射，转到步骤(8.13)；

(8.11)当ξ≤(μ_a+μ_afp)/(μ_a+μ_afp+μ_s)时，则激发光光子被吸收，进入步骤(8.12)；当ξ＞(μ_a+μ_afp)/(μ_a+μ_afp+μ_s)时，则激发光光子散射，转到步骤(8.13)；

(8.12)通过设定随机数判定光子是否被生物组织中荧光团吸收而产生荧光，荧光的初始方向由各向同性散射的偏转角和方位角决定；当ξ≤μ_afp/(μ_a+μ_afp)时，则产生荧光，更新sleft＝0，转到步骤(8.2)；当ξ＞μ_afp/(μ_a+μ_afp)时，则设定光子死亡，停止对光子追踪，转到步骤(8.14)；

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

μx′=sinθ′1-μz2(μxμzcosψ-μysinψ)+μxcosθ′μy′=sinθ′1-μz2(μyμzcosψ-μxsinψ)+μycosθ′μz′=-sinθ′cosψ1-μz2+μzcosθ′]]>

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

散射方向更新后，转到步骤(8.2)；

(8.14)判定光子是否为最后一个光子；如果是，则继续下一步；如果否，则重复上述过程；

(5)如激发光光子被目标生物组织的荧光团吸收并产生荧光，追踪每个荧光光子的传输过程；将其中耗费时间长的大规模数据并行、高计算密度的步骤安排在GPU上并行执行，其余步骤安排在CPU上串行执行；CPU分配内存，用于存放GPU输出数据，将计算后显存上的数据复制到内存上，各子节点返回运行信息给主节点；

(6)追踪所有光子后输出所有逸出的激发光光子信息和荧光光子信息，释放内存和显存空间并退出CUDA。