一种基于集群式GPU加速的荧光蒙特卡罗模拟方法技术

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

【技术实现步骤摘要】
一种基于集群式GPU加速的荧光蒙特卡罗模拟方法
本专利技术涉及一种基于集群式GPU加速的荧光蒙特卡罗模拟方法，属于计算机科学、数学仿真和生物医学工程领域。
技术介绍
在组织光学领域，光传输过程建模的目标是发展任意组织内光子分布特性的准确快速的数学定量方法，它包括两个层次的含义：其一是建立模型作为评价其它特定应用模型的“金标准”；其二是针对组织光学具体应用中的特定条件发展快速有效地实用模型。蒙特卡罗方法的基本特征是对随机性问题进行仿真，能有效解决随机性问题，甚至对许多确定性方法所难以解决的随机性问题都能较方便地解决。蒙特卡罗法被该领域研究者奉为金标准。现有的蒙特卡罗代码只考虑了激发光在生物组织中的传播，而忽视了生物组织中荧光团的存在。激发光被荧光团吸收后可产生荧光，可在生物组织中继续传播。由于生物组织的不同组成成份会发出不同特征的荧光，通过追踪这些荧光并获取其信息，有利于生物医学的诊断和治疗。并且，随着模拟的光子数目增加，模拟源的数目不断增加，组织结构越来越复杂，计算机模拟的时间也成倍在增加，现有的代码已经不能满足对模拟耗时的需求。专利技术一种优化程序结构，采用集群和GPU加速实现并行高速运算的荧光蒙特卡罗模拟的代码可以满足此需求。
技术实现思路
本专利技术目的在于提供了一种基于集群式GPU加速的荧光蒙特卡罗模拟方法，该技术可以同时模拟多个源，大大节省了模拟光子在生物组织中传输的时间，并分别追踪了生物组织内激发光和荧光的光子传输。一种基于集群式GPU加速的荧光蒙特卡罗模拟方法，其特征在于包括以下步骤：(1)将目标生物组织的空间结构分割成一个三维体素模型，设定一个与三维体素模型大小相同的三维数字矩阵，矩阵中每个元素的数值为标识的组织类型；设定目标生物组织的各类组织光学特性参数，组织光学特性参数为激发光的吸收系数、激发光的散射系数、荧光的吸收系数、荧光的散射系数、荧光团的吸收系数、荧光团的折射系数和荧光团的各向异性因子；(2)通过MPI这种消息传递通信协议，主节点获取各子节点的GPU设备的数量，将需要计算的源的数量平均分配给每个子节点的GPU设备，各子节点确定GPU设备上网格和块的维度和尺寸，并根据该节点GPU设备的数目开辟相应数目的子进程，将计算任务分配给各个子进程，启动并行编程与计算平台CUDA，各子节点CPU分配内存空间、显存空间，将要计算的数据从内存复制到显存上；(3)将入射光源表征为设定数目的光子的集合，将入射光源位置和入射光方向赋给每个光子的初始位置和方向；(4)在目标生物组织内追踪每个激发光光子的传输过程，将其中耗费时间长的大规模数据并行、高计算密度的步骤安排在GPU上并行执行，其余步骤安排在CPU上串行执行；CPU分配内存，用于存放GPU输出数据，将计算后显存上的数据复制到内存上，各子节点返回运行信息给主节点；(5)如激发光光子被目标生物组织的荧光团吸收并产生荧光，追踪每个荧光光子的传输过程；将其中耗费时间长的大规模数据并行、高计算密度的步骤安排在GPU上并行执行，其余步骤安排在CPU上串行执行；CPU分配内存，用于存放GPU输出数据，将计算后显存上的数据复制到内存上，各子节点返回运行信息给主节点；(6)追踪所有光子后输出所有逸出的激发光光子信息和荧光光子信息，释放内存和显存空间并退出CUDA。步骤(4)具体包括以下步骤进行：(4.1)投放激发光光子；(4.2)激发光光子移动一个由随机数确定的步长，判定激发光光子在此过程中是否遇到外界边界，如果遇到边界，则转到步骤(4.3)，否则转到步骤(4.4)；(4.3)如果激发光光子发生反射，则更新激发光光子的步长和方向余弦信息，如果激发光光子发生折射到组织外，则设定光子死亡，停止对光子追踪；(4.4)通过设定随机数判定激发光光子是被吸收或散射，如果激发光光子被吸收，则转到步骤(4.5)，如果激发光光子被散射，则更新光子方向余弦信息，转到步骤(4.2)；(4.5)如果激发光光子被吸收不产生荧光，则设定光子死亡，停止对光子追踪；步骤(5)具体包括以下步骤进行：(5.1)如果激发光光子被荧光团吸收并产生荧光，荧光的初始方向由各向同性散射的偏转角和方位角决定，继续追踪荧光光子下一步的移动；(5.2)如果荧光光子散射，则更新光子方向余弦信息，转到步骤(5.3)；如果荧光光子被吸收，则设定光子死亡，停止对光子追踪；(5.3)荧光光子移动一个由随机数确定的步长，判定荧光光子在此过程中是否遇到外界边界，如果遇到边界，则转到步骤(5.4)，否则转到步骤(5.2)；(5.4)如果荧光光子发生反射，则更新荧光光子的步长和方向余弦信息，转到步骤(5.2)；如果荧光光子发生折射到组织外，则设定光子死亡，停止对光子追踪。基于步骤(6)中的输出，进行统计计算并转化为下述所需获取的传输特性中的一个或任意几个：将同一入射方向内，所有逸出的激发光光子能量相加，即可得激发光逸出能量与漫射角度之间的关系曲线；将同一入射方向内，所有逸出的荧光光光子能量相加，即可得荧光逸出能量与漫射角度之间的关系曲线。本专利技术建立了一种超高速的荧光蒙特卡罗模拟方法，可以同时获取逸出的激发光光子信息和荧光光子信息。这些丰富的信息为光学成像系统的优化提供了依据，为光诊断和光治疗提供了精确的指导信息。附图说明图1为本专利技术的基本流程图。具体实施方式结合附图对本专利技术作进一步的描述。如图1所示，本专利技术的实施步骤如下：(1)将目标生物组织的空间结构描述为一个三维数字矩阵，即组织模型。矩阵中的元素对应目标生物组织的体素，每个体素的数值为标识组织类型的数字；体素越小，描述组织结构的组织模型越逼近于真实生物组织结构，模拟得到的光子传输特性精度就越高；(2)设定目标生物组织的各类组织特性参数，组织光学特性参数为激发光的吸收系数、激发光的散射系数、荧光的吸收系数、荧光的散射系数、荧光团的吸收系数、荧光团的折射系数和荧光团的各向异性因子；这些参数值的精度也会影响模拟精度，可查阅相应文献来获取这些参数。值得注意是，这些参数的选定要求与目标光源波长一致；(3)通过MPI这种消息传递通信协议，主节点获取各子节点的GPU设备的数量，读入数据，将需要计算的源的数量平均分配给每个子节点的GPU设备；(4)各子节点计算该节点各GPU设备的可用资源，确定网格和块的维度和尺寸；(5)各子节点CPU的主进程根据该节点GPU设备的数目开辟相应数目的子进程，将计算任务分配给各个子进程；(6)启动CUDA，各子节点CPU分配内存空间、显存空间，将要计算的数据从内存复制到显存上；(7)将入射光源表征为设定数目光子的集合，将入射光源位置和入射光方向赋给每个光子的初始位置和方向；(8)追踪每个激发光光子的传输过程；将其中耗费时间长的大规模数据并行、高计算密度的步骤安排在GPU上并行执行，其余步骤安排在CPU上串行执行，CPU分配内存，用于存放GPU输出数据，将计算后显存上的数据复制到内存上，各子节点返回运行信息给主节点；追踪每个激发光光子的传输过程具体步骤为：(8.1)投放激发光光子，如果光子处在组织表面之上，将该位置设定为光子发射位置；如果光子处在组织表面之外，将光子自动移向组织表面，其具体操作使用迭代方法；如果光子处在组织内部，则将依照此位置直接发射光子，即光子不发本文档来自技高网...

【技术保护点】
一种基于集群式GPU加速的荧光蒙特卡罗模拟方法，其特征在于包括以下步骤： (1)将目标生物组织的空间结构分割成一个三维体素模型，设定一个与三维体素模型大小相同的三维数字矩阵，矩阵中每个元素的数值为标识的组织类型；设定目标生物组织的各类组织光学特性参数，组织光学特性参数为激发光的吸收系数、激发光的散射系数、荧光的吸收系数、荧光的散射系数、荧光团的吸收系数、荧光团的折射系数和荧光团的各向异性因子； (2)通过MPI这种消息传递通信协议，主节点获取各子节点的GPU设备的数量，将需要计算的源的数量平均分配给每个子节点的GPU设备，各子节点确定GPU设备上网格和块的维度和尺寸，并根据该节点GPU设备的数目开辟相应数目的子进程，将计算任务分配给各个子进程，启动并行编程与计算平台CUDA，各子节点CPU分配内存空间、显存空间，将要计算的数据从内存复制到显存上； (3)将入射光源表征为设定数目的光子的集合，将入射光源位置和入射光方向赋给每个光子的初始位置和方向； (4)在目标生物组织内追踪每个激发光光子的传输过程，将其中耗费时间长的大规模数据并行、高计算密度的步骤安排在GPU上并行执行，其余步骤安排在CPU上串行执行；CPU分配内存，用于存放GPU输出数据，将计算后显存上的数据复制到内存上，各子节点返回运行信息给主节点； (5)如激发光光子被目标生物组织的荧光团吸收并产生荧光，追踪每个荧光光子的传输过程；将其中耗费时间长的大规模数据并行、高计算密度的步骤安排在GPU上并行执行，其余步骤安排在CPU上串行执行；CPU分配内存，用于存放GPU输出数据，将计算后显存上的数据复制到内存上，各子节点返回运行信息给主节点； (6)追踪所有光子后输出所有逸出的激发光光子信息和荧光光子信息，释放内存和显存空间并退出CUDA。...

【技术特征摘要】
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)计算光子移动一步的步长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,y,z表示光子当前位置；如果第一项最小，那么光子在界面上的作用点表达式为：这里，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，设置光子反射，其方向变为：如果θ＞sin-1(n2/n1)，设置光子发生全反射，其方向变为：在其他情况下，设置光子...

【专利技术属性】
技术研发人员：骆清铭，邓勇，罗召洋，江旭，
申请(专利权)人：华中科技大学，
类型：发明
国别省市：湖北;42