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

摘要

本发明涉及一种基于集群式GPU加速的荧光蒙特卡罗方法，能同时模拟多个源，大大节省了模拟光子在生物组织中传输的时间，考虑了组织中荧光团对激发光的吸收作用，分别追踪了生物组织内激发光和荧光的光子传输。该方法精度高，能够获取真实生物组织中的光传输信息，这些丰富的信息为光学成像系统的优化提供了依据，为光诊断和光治疗提供了精确的指导信息。

主权项

一种基于集群式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)计算光子移动一步的步长ss＝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)，其表达式为：然后，根据该矢量表达式计算光子与界面发生作用的入射角θ和折射角θt,其计算式为：cosθ＝a·ux+b·uy+c·uzcosθt＝(1‑(1‑cos2θ)·n12/n22)1/2这里，kx,ky,kz表示光子遇到边界之前的方向余弦值，n1和n2表示遇到界面前后所在组织的折射系数；最后分四种情况计算光子遇到界面之后的方向：如果cosθ＝0,设置光子遇到边界之后方向不变；如果cosθ＝1，当ε＞(n2‑n1)2/(n2+n1)2,设置光子遇到边界之后方向不变，当ε≤(n2‑n1)2/(n2+n1)2，设置光子反射，其方向变为：kxi=kx-2cosθ*akyi=ky-2cosθ*bkzi=kz-2cosθ*c]]>如果θ＞sin‑1(n2/n1)，设置光子发生全反射，其方向变为：kxi=kx-2cosθ*akyi=ky-2cosθ*bkzi=kz-2cosθ*c]]>在其他情况下，设置光子与界面发生作用后的方向为：其中，r的表达式为：(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。