一种基于集群式gpu加速的荧光蒙特卡罗模拟方法
【专利摘要】本发明涉及一种基于集群式GPU加速的荧光蒙特卡罗方法，能同时模拟多个源，大大节省了模拟光子在生物组织中传输的时间，考虑了组织中荧光团对激发光的吸收作用，分别追踪了生物组织内激发光和荧光的光子传输。该方法精度高，能够获取真实生物组织中的光传输信息，这些丰富的信息为光学成像系统的优化提供了依据，为光诊断和光治疗提供了精确的指导信息。
【专利说明】一种基于集群式GPU加速的焚光蒙特卡罗模拟方法
【技术领域】
[0001] 本发明涉及一种基于集群式GPU加速的荧光蒙特卡罗模拟方法，属于计算机科 学、数学仿真和生物医学工程领域。
【背景技术】
[0002] 在组织光学领域，光传输过程建模的目标是发展任意组织内光子分布特性的准确 快速的数学定量方法，它包括两个层次的含义：其一是建立模型作为评价其它特定应用模 型的"金标准";其二是针对组织光学具体应用中的特定条件发展快速有效地实用模型。蒙 特卡罗方法的基本特征是对随机性问题进行仿真，能有效解决随机性问题，甚至对许多确 定性方法所难以解决的随机性问题都能较方便地解决。蒙特卡罗法被该领域研究者奉为金 标准。
[0003] 现有的蒙特卡罗代码只考虑了激发光在生物组织中的传播，而忽视了生物组织中 荧光团的存在。激发光被荧光团吸收后可产生荧光，可在生物组织中继续传播。由于生物组 织的不同组成成份会发出不同特征的荧光，通过追踪这些荧光并获取其信息，有利于生物 医学的诊断和治疗。并且，随着模拟的光子数目增加，模拟源的数目不断增加，组织结构越 来越复杂，计算机模拟的时间也成倍在增加，现有的代码已经不能满足对模拟耗时的需求。 发明一种优化程序结构，采用集群和GPU加速实现并行高速运算的荧光蒙特卡罗模拟的代 码可以满足此需求。
【发明内容】
[0004] 本发明目的在于提供了一种基于集群式GPU加速的荧光蒙特卡罗模拟方法，该技 术可以同时模拟多个源，大大节省了模拟光子在生物组织中传输的时间，并分别追踪了生 物组织内激发光和荧光的光子传输。
[0005] -种基于集群式GPU加速的荧光蒙特卡罗模拟方法，其特征在于包括以下步骤：
[0006] (1)将目标生物组织的空间结构分割成一个三维体素模型，设定一个与三维体素 模型大小相同的三维数字矩阵，矩阵中每个元素的数值为标识的组织类型；设定目标生物 组织的各类组织光学特性参数，组织光学特性参数为激发光的吸收系数、激发光的散射系 数、荧光的吸收系数、荧光的散射系数、荧光团的吸收系数、荧光团的折射系数和荧光团的 各向异性因子；
[0007] (2)通过MPI这种消息传递通信协议，主节点获取各子节点的GPU设备的数量，将 需要计算的源的数量平均分配给每个子节点的GPU设备，各子节点确定GPU设备上网格和 块的维度和尺寸，并根据该节点GPU设备的数目开辟相应数目的子进程，将计算任务分配 给各个子进程，启动并行编程与计算平台CUDA，各子节点CPU分配内存空间、显存空间，将 要计算的数据从内存复制到显存上；
[0008] (3)将入射光源表征为设定数目的光子的集合，将入射光源位置和入射光方向赋 给每个光子的初始位置和方向；
[0009] (4)在目标生物组织内追踪每个激发光光子的传输过程，将其中耗费时间长的大 规模数据并行、高计算密度的步骤安排在GPU上并行执行，其余步骤安排在CPU上串行执 行；CPU分配内存，用于存放GPU输出数据，将计算后显存上的数据复制到内存上，各子节点 返回运行信息给主节点；
[0010] (5)如激发光光子被目标生物组织的荧光团吸收并产生荧光，追踪每个荧光光子 的传输过程；将其中耗费时间长的大规模数据并行、高计算密度的步骤安排在GPU上并行 执行，其余步骤安排在CPU上串行执行；CPU分配内存，用于存放GPU输出数据，将计算后显 存上的数据复制到内存上，各子节点返回运行信息给主节点；
[0011] (6)追踪所有光子后输出所有逸出的激发光光子信息和荧光光子信息，释放内存 和显存空间并退出CUDA。
[0012] 步骤（4)具体包括以下步骤进行：
[0013] (4. 1)投放激发光光子；
[0014] (4. 2)激发光光子移动一个由随机数确定的步长，判定激发光光子在此过程中是 否遇到外界边界，如果遇到边界，则转到步骤（4. 3)，否则转到步骤（4.4);
[0015] (4. 3)如果激发光光子发生反射，则更新激发光光子的步长和方向余弦信息，如果 激发光光子发生折射到组织外，则设定光子死亡，停止对光子追踪；
[0016] (4. 4)通过设定随机数判定激发光光子是被吸收或散射，如果激发光光子被吸收， 则转到步骤（4. 5)，如果激发光光子被散射，则更新光子方向余弦信息，转到步骤（4.2);
[0017] (4. 5)如果激发光光子被吸收不产生荧光，则设定光子死亡，停止对光子追踪；
[0018] 步骤（5)具体包括以下步骤进行：
[0019] (5. 1)如果激发光光子被荧光团吸收并产生荧光，荧光的初始方向由各向同性散 射的偏转角和方位角决定，继续追踪荧光光子下一步的移动；
[0020] (5. 2)如果荧光光子散射，则更新光子方向余弦信息，转到步骤（5. 3);如果荧光 光子被吸收，则设定光子死亡，停止对光子追踪；
[0021] (5. 3)荧光光子移动一个由随机数确定的步长，判定荧光光子在此过程中是否遇 到外界边界，如果遇到边界，则转到步骤（5. 4)，否则转到步骤（5.2);
[0022] (5. 4)如果荧光光子发生反射，则更新荧光光子的步长和方向余弦信息，转到步骤 (5. 2);如果荧光光子发生折射到组织外，则设定光子死亡，停止对光子追踪。
[0023] 基于步骤（6)中的输出，进行统计计算并转化为下述所需获取的传输特性中的一 个或任意几个：
[0024] 将同一入射方向内，所有逸出的激发光光子能量相加，即可得激发光逸出能量与 漫射角度之间的关系曲线；
[0025] 将同一入射方向内，所有逸出的荧光光光子能量相加，即可得荧光逸出能量与漫 射角度之间的关系曲线。
[0026] 本发明建立了一种超高速的荧光蒙特卡罗模拟方法，可以同时获取逸出的激发光 光子信息和荧光光子信息。这些丰富的信息为光学成像系统的优化提供了依据，为光诊断 和光治疗提供了精确的指导彳目息。
【专利附图】
【附图说明】
[0027] 图1为本发明的基本流程图。
【具体实施方式】
[0028] 结合附图对本发明作进一步的描述。
[0029] 如图1所示，本发明的实施步骤如下：
[0030] (1)将目标生物组织的空间结构描述为一个三维数字矩阵，即组织模型。矩阵中的 元素对应目标生物组织的体素，每个体素的数值为标识组织类型的数字；体素越小，描述组 织结构的组织模型越逼近于真实生物组织结构，模拟得到的光子传输特性精度就越高；
[0031] (2)设定目标生物组织的各类组织特性参数，组织光学特性参数为激发光的吸收 系数、激发光的散射系数、荧光的吸收系数、荧光的散射系数、荧光团的吸收系数、荧光团的 折射系数和荧光团的各向异性因子；这些参数值的精度也会影响模拟精度，可查阅相应文 献来获取这些参数。值得注意是，这些参数的选定要求与目标光源波长一致；
[0032] (3)通过MPI这种消息传递通信协议，主节点获取各子节点的GPU设备的数量，读 入数据，将需要计算的源的数量平均分配给每个子节点的GPU设备；
[0033] (4)各子节点计算该节点各GPU设备的可用资源，确定网格和块的维度和尺寸；
[0034] (5)各子节点CPU的主进程根据该节点GPU设备的数目开辟相应数目的子进程，将 计算任务分配给各个子进程；
[0035] (6)启动CUDA，各子节点CPU分配内存空间、显存空间，将要计算的数据从内存复 制到显存上；
[0036] (7)将入射光源表征为设定数目光子的集合，将入射光源位置和入射光方向赋给 每个光子的初始位置和方向；
[0037] (8)追踪每个激发光光子的传输过程；将其中耗费时间长的大规模数据并行、高 计算密度的步骤安排在GPU上并行执行，其余步骤安排在CPU上串行执行，CPU分配内存， 用于存放GPU输出数据，将计算后显存上的数据复制到内存上，各子节点返回运行信息给 主节点；
[0038] 追踪每个激发光光子的传输过程具体步骤为：
[0039] (8. 1)投放激发光光子，如果光子处在组织表面之上，将该位置设定为光子发射位 置；如果光子处在组织表面之外，将光子自动移向组织表面，其具体操作使用迭代方法；如 果光子处在组织内部，则将依照此位置直接发射光子，即光子不发生与表面的相互作用，而 直接开始光子在组织内部的传输；
[0040] (8. 2)计算光子的剩余步长sleft:Sleft = In ;
[0041] (8. 3)计算光子移动一步的步长s
[0042] s = min (In () /t, min (dx, dy, dz))
[0043] 这里，是由计算机伪随机数产生器产生的均匀分布在（0, I)的伪随机数， dx，dy，dz依次为组织模型中每个体素的长、宽、高；当光子为激发光光子时：t = a s+afp ;当 光子为荧光光子时：t = af sf。a和s分别是激发光的吸收系数和散射系数；af和sf分别是荧 光的吸收系数和散射系数；afp是荧光团的吸收系数；
[0044] (8. 4)按照光子当前方向和（8. 3)中计算出的步长移动光子一步，判断该过程中 光子是否穿过不同组织之间的界面；如果是，进入步骤（8.5);如果否，将光子移动一个步 长，更新光子当前位置，然后进入步骤（8.8);
[0045] (8. 5)确定光子在界面上的撞击点，方法为：首先找出如下表达式中的最小项：
[0046] x_t = (I [x/dx]-x/dx)/ux
[0047] y_t = (l[y/dy]-y/dy)/uy
[0048] z_t = (I [z/dz]-z/dz)/uz
[0049] 这里，x，y，z表示光子当前位置；如果第一项最小，那么光子在界面上的作用点表 达式为：
[0050] x0 = X ux+x_t
[0051] y0 = y uy+y_t
[0052] z0 = z uz+z_t
[0053] 这里，x0, y0, zO是光子在界面的作用点，其他情况类推；
[0054] (8. 6)计算光子在（8. 5)中确定的撞击点与界面发生相互作用后的方向，具体方 法为：首先计算界面法向矢量（a，b，c)，其表达式为：
【权利要求】
1. 一种基于集群式GPU加速的荧光蒙特卡罗模拟方法，其特征在于包括以下步骤： (1) 将目标生物组织的空间结构分割成一个三维体素模型，设定一个与三维体素模型 大小相同的三维数字矩阵，矩阵中每个元素的数值为标识的组织类型；设定目标生物组织 的各类组织光学特性参数，组织光学特性参数为激发光的吸收系数、激发光的散射系数、荧 光的吸收系数、荧光的散射系数、荧光团的吸收系数、荧光团的折射系数和荧光团的各向异 性因子； (2) 通过MPI这种消息传递通信协议，主节点获取各子节点的GPU设备的数量，将需要 计算的源的数量平均分配给每个子节点的GPU设备，各子节点确定GPU设备上网格和块的 维度和尺寸，并根据该节点GPU设备的数目开辟相应数目的子进程，将计算任务分配给各 个子进程，启动并行编程与计算平台CUDA，各子节点CPU分配内存空间、显存空间，将要计 算的数据从内存复制到显存上； (3) 将入射光源表征为设定数目的光子的集合，将入射光源位置和入射光方向赋给每 个光子的初始位置和方向； (4) 在目标生物组织内追踪每个激发光光子的传输过程，将其中耗费时间长的大规模 数据并行、高计算密度的步骤安排在上并行执行，其余步骤安排在CPU上串行执行；CPU 分配内存，用于存放GPU输出数据，将计算后显存上的数据复制到内存上，各子节点返回运 行信息给主节点； (5) 如激发光光子被目标生物组织的荧光团吸收并产生荧光，追踪每个荧光光子的传 输过程；将其中耗费时间长的大规模数据并行、高计算密度的步骤安排在GPU上并行执行， 其余步骤安排在CPU上串行执行；CPU分配内存，用于存放GPU输出数据，将计算后显存上 的数据复制到内存上，各子节点返回运行信息给主节点； (6) 追踪所有光子后输出所有逸出的激发光光子信息和荧光光子信息，释放内存和显 存空间并退出CUDA。
2. 根据权利要求1所述的基于集群式GPU加速的荧光蒙特卡罗模拟方法，其特征在于 步骤（4)具体包括以下步骤进行： (2. 1)投放激发光光子； (2.2) 激发光光子移动一个由随机数确定的步长，判定激发光光子此过程中是否遇到 外界边界，如果遇到边界，则转到步骤（2. 3)，否则转到步骤（2.4); (2. 3)如果激发光光子发生反射，则更新激发光光子的步长和方向余弦信息，如果激发 光光子发生折射到组织外，则设定光子死亡，停止对光子追踪； (2. 4)通过设定随机数判定激发光光子是被吸收或散射，如果激发光光子被吸收，则转 到步骤（2. 5)，如果激发光光子被散射，则更新光子方向余弦信息，转到步骤（2.2); (2. 5)如果激发光光子被吸收不产生荧光，则设定光子死亡，停止对光子追踪。
3. 根据权利要求1所述的基于集群式GPU加速的荧光蒙特卡罗模拟方法，其特征在于 步骤（5)具体包括以下步骤进行： (3. 1)如果激发光光子被荧光团吸收并产生荧光，荧光的初始方向由各向同性散射的 偏转角和方位角决定，继续追踪荧光光子下一步的移动； (3.2) 如果荧光光子散射，则更新光子方向余弦信息，转到步骤（3.3);如果荧光光子 被吸收，则设定光子死亡，停止对光子追踪； (3. 3)荧光光子移动一个由随机数确定的步长，判定荧光光子在此过程中是否遇到外 界边界，如果遇到边界，则转到步骤（3. 4)，否则转到步骤（3.2); (3.4)如果荧光光子发生反射，则更新荧光光子的步长和方向余弦信息，转到步骤 (3. 2);如果荧光光子发生折射到组织外，则设定光子死亡，停止对光子追踪。
4.根据权利1所述的基于集群式GPU加速的荧光蒙特卡罗模拟方法，其特征在于：基 于步骤￠)中的输出，进行统计计算并转化为下述所需获取的传输特性中的一个或任意几 个： 将同一入射方向内，所有逸出的激发光光子能量相加，即可得激发光逸出能量与漫射 角度之间的关系曲线； 将同一入射方向内，所有逸出的荧光光光子能量相加，即可得荧光逸出能量与漫射角 度之间的关系曲线。
【文档编号】G06F19/12GK104331641SQ201410534214
【公开日】2015年2月4日 申请日期:2014年10月11日 优先权日:2014年10月11日
【发明者】骆清铭, 邓勇, 罗召洋, 江旭 申请人:华中科技大学