CN104605823A - 一种基于pfMC模型的荧光扩散光学断层图像重建方法 - Google Patents

Abstract

本发明公开了一种基于微扰荧光蒙特卡罗(pfMC)模型的荧光扩散光学断层图像重建方法，属于生物医学工程技术领域，包括：确定探测区域，在探测区域选择若干扫描点，获取探测器上相对应的荧光光强分布；建立描述组织光学参数空间结构的三维数字模型，根据光源扫描的位置和方向进行激发光子的前向白蒙特卡罗(white Monte Carlo)模拟，跟踪激发光子并记录到达探测器的光子路径；依据pfMC方法计算荧光光子权重，并计算出荧光雅可比矩阵；最后通过GPU集群迭代重建计算出组织内的荧光团位置及其荧光系数；本发明的优点是，基于体素化的微扰MC荧光模型，基于GPU集群加速的重建迭代过程，为三维荧光断层成像系统提供了一种快速的、准确的重建方法。

Description

一种基于pfMC模型的荧光扩散光学断层图像重建方法

技术领域

本发明属于生物医学工程技术领域，特别是涉及一种基于微扰荧光蒙特卡罗(pfMC)模型的荧光扩散光学断层图像重建方法。

背景技术

荧光扩散光学断层成像利用激光光束激发生物体内荧光，采用高灵敏度探测器在生物表面探测激发出的荧光信号。利用探测器对小动物表面进行多角度拍摄，根据探测器上拍摄的信号结合数学物理模型和优化算法逆向推导生物体内荧光参数的分布，即对荧光标记物进行三维定位和荧光标记物浓度测定。荧光扩散光学断层成像是一种使用荧光探针的生物医学光学功能成像技术，通过收集从组织表面透射出来的荧光分布信息，重建组织内部的荧光系数分布。这种技术具有无创,对机体无害等优点，因而被广泛应用于癌症诊断，药物研发以及基因表达可视化的研究。

为了重建荧光系数分布，即求解荧光扩散光学断层成像的逆问题，首先需要建立一种能够准确模拟组织中激发光和荧光光子的产生与传播的正向模型。最早建立的模型之一是Grunbaum^[1]于1992提出的用一步或两步马尔可夫随机过程来描述光子在组织中的传播。在此统计模型的基础上后来又发展出了随机行走模型和蒙特卡罗模型。另一类重要的模型是麦克斯韦辐射传输模型及其近似模型，如扩散模型。麦克斯韦辐射传输模型能够准确的描述光子在各种组织内传播，但是其在复杂边界的组织中的求解是异常复杂的，因此各种近似模型被发展出来替代辐射传输模型。扩散模型是一种易于求解的近似模型(如用有限元求解)，并且在一定组织条件下能够比较准确的模拟光子在组织中的传播，但是在低散射系数或高吸收系数的组织中，却存在很大的误差。蒙特卡洛方法能够准确描述光子在组织内的输运行为，适用于各种复杂边界及各种复杂异构情形。然而这种统计模拟十分耗费计算机资源，一般需要耗费大量时间。

其次是在正向模型的基础上建立雅可比矩阵，然后再迭代重建荧光系数分布。然而各种正向模型中，适应于各种复杂条件三维组织模型，精确度高且具有实际可执行性的只有蒙特卡罗模型。但是直接的荧光蒙特卡罗模型是不适应于建立荧光雅可比矩阵的。为了快速建立雅可比矩阵，微扰法被应用到荧光蒙特卡罗(pMC)方法中。在微扰荧光蒙特卡罗方法中，由white MC获取光子在体素内的路径信息，再由pMC方法对光子路径进行处理，建立雅可比矩阵。这种微扰法它不受成像系统的限制，基于体素的路径信息存贮能够大为减少数据量，具有速度和精度的优势。

参考文献

[1]F.A.Grunbaum,“Diffuse tomography:Computational aspects of the isotropic case,”Inverse Problems 8,421-433(1992).

发明内容

为了克服上述现有技术的不足，本发明提供了一种基于微扰荧光蒙特卡罗(pfMC)模型的荧光扩散光学断层图像重建方法，本发明的目的是在保持MC模拟“金标准”的精度，同时改善MC模拟的耗时特性，实现快速荧光断层图像重建。

本发明所采用的技术方案是：

一种基于pfMC模型的荧光扩散光学断层图像重建方法，其特征在于包括下列步骤：

(1)确定探测区域，选择光源扫描位置，并获取相对应的探测器上的荧光强度分布；

(2)数字化待探测对象，建立描述组织光学参数三维分布的体素化模型，在不同光源位置与方向进行白蒙特卡罗模拟，记录在探测器收集到的激发光子的权重及在每个体素内的路径信息和相应的物理量；

(3)依据pfMC模型，计算荧光光子权重；

(4)由激发光子的路径信息构建荧光雅可比矩阵；

(5)根据计算得到荧光光子权重和荧光雅可比矩阵，迭代重建计算出组织内的荧光团位置及荧光系数。

步骤(2)中，采用GPU集群加速整个白蒙特卡罗过程，记录被探测器采集到的激发光子的相关物理量，被记录在体素v内的相关物理量L_v为：

$L_v=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\tilde{l}}_i=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{l}_i\·e^{{\mathrm{\μ}}_{\mathrm{av}}{l}_i}$

式中：n为体素v内的散射次数，μ_av为体素v内的吸收系数，l_i为体素v内的第i次散射路径长度。

GPU集群加速整个白蒙特卡罗过程具体如下：

2.1在集群上，依据光源数目及集群节点数目分配各个节点上进行MC模拟的次数，并在集群启动后把相应的光源位置和方向信息输送到各个节点上；

2.2在各个节点上，根据要模拟的光子数目分配内存，将组织模型信息分配到GPU内的纹理寄存器，光学参数信息分配到GPU内的常数寄存器来提高读取速度，并根据GPU的技术信息确定最优的CUDA运行中并行块和线程的数目；

2.3在GPU内的各个线程上，每个线程上模拟一个光子包，并跟踪光子发射后在组织内散射及吸收的过程，记录光子在每个体素内的路径相关物理量，直到光子溢出组织后终止光子模拟。

步骤(3)，计算荧光权重w^em(p)计算公式如下：

$w^{\mathrm{em}}\left(p\right)=w^{\mathrm{ex}}\left(p\right)\mathrm{\η}\underset{v=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\μ}}_{\mathrm{afv}}{L}_v$

式中：p为位置及单位方向的六维矢量，k为光子在组织内的路径所经过的体素数目，μ_afv是体素v内荧光系数，η是量子产率，w^ex为激发光子权重。

步骤(4)中，荧光雅可比矩阵的计算公式如下：

$\begin{array}{c}{J}_{\mathrm{jd}}^v=\frac{\∂w^{\mathrm{em}}\left(p\right){|}_{\det \mathrm{ectord}}}{\∂{\mathrm{\μ}}_{\mathrm{afv}}}\\ =w^{\mathrm{ex}}\left(p\right){|}_{\det \mathrm{ectord}}\mathrm{\η}{L}_v\end{array}$

式中：j为节点i上的光源编号，d为探测器编号，v为需要重建的模型体素编号，J为荧光雅可比矩阵，为对应体素v、光源j、探测器d的荧光雅可比矩阵元。

步骤(5)整个荧光重建过程均通过GPU集群加速，对各个节点上计算得到的荧光雅可比矩阵矩阵块进行处理，节点i上矩阵块计算公式如下：

D_i＝[D₁₁…D_1t；…；D_s1…D_st]^T

式中：s为节点i上的光源数目，t为探测器总数，n为需要重建的模型体素数目，D_i为节点i上所有光源下的探测器探测到的荧光强度分布与pfMC模拟得到的荧光分布的差值，D_st为节点i上第s光源下的探测器中第t个探测单元探测到的荧光强度分布与相对应的dfMC模拟得到的荧光分布的差值。

步骤(5)，由矩阵块计算的中间量J_iJ_i ^TX与J_iD_i被传输到主节上再进行共轭梯度迭代求解荧光分布，这些中间量比荧光雅可比矩阵块占用的存贮空间少,迭代计算公式如下：

$\underset{X}{\min }f\left(X\right)=\left|\right|J^T\mathrm{JX}-J^{T}D\left|\right|=\left|\right|\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(J_i{J_i}^{T}X-J_i{D}_i)\left|\right|$

式中：i为节点编号，N为总节点数目，X为荧光团吸收系数分布矩阵，J_i ^T为矩阵J_i的转置矩阵，f(X)为目标函数。

与现有技术相比，本发明的有益效果为：能够实现对生物组织内荧光团的三维分布与荧光系数的重建；具有重建精度高，速度快的优势。

附图说明

图1为为本发明GPU集群下并行加速下白蒙特卡罗模拟流程图；

图2为本发明GPU集群并行加速下pfMC三维荧光重建流程图；

图3a、图3b分别为本发明重建体模实验结果图。

具体实施方式

下面结合附图对本发明进一步说明。

本发明的实施步骤如下：

1、确定探测区域，选择光源扫描位置，并在不同扫描位置下得到探测器上的荧光强度分布D，以供荧光反演求解；

2、数字化探测组织，建立体素模型，在不同光源位置与方向下进行GPU集群加速下白蒙特卡罗模拟，获取被探测器收集到的激发光子在每个体素内的路径相关物理量Lv，GPU集群加速下的白蒙特卡罗过程如图1所示，步骤如下：

2.1在集群上，依据光源数目及集群节点数目分配各个节点上进行MC模拟的次数，并在集群启动后把相应的光源位置和方向信息输送到各个节点上；

2.2在各个节点上，根据要模拟的光子数目分配内存，可以将组织模型信息分配到GPU内的纹理寄存器，光学参数信息分配到GPU内的常数寄存器来提高读取速度，并根据GPU的技术信息确定最优的CUDA运行中并行块和线程的数目；

2.3在GPU内的各个线程上，每个线程上模拟一个光子包，并跟踪光子发射后在组织内散射及吸收的过程，记录光子在每个体素内的路径相关物理量，直到光子溢出组织后终止光子模拟。

3、GPU集群下，在每个光源下探测到的荧光强度分布D及MC模拟得到的大量沿光子路径记录的物理信息进行荧光分布反演求解，具体步骤如图2所示：

3.1在集群上，根据白蒙特卡罗模拟时的集群分配信息进行重建过程中的集群分配，使得每个节点上对应有一部分光源的激发光子路径信息，如果总共有n个光源，那么每个节点上有1/n个光源的MC模拟路径信息；

3.2在每个节点上，载入相对应的激发光子路径信息，并依据pfMC方法计算相对应节点i的荧光荧光雅可比矩阵块J_i，计算公式如下：

$\begin{array}{c}{J}_{\mathrm{jd}}^v=\frac{\∂w^{\mathrm{em}}\left(p\right){|}_{\det \mathrm{ectord}}}{\∂{\mathrm{\μ}}_{\mathrm{afv}}}\\ =w^{\mathrm{ex}}\left(p\right){|}_{\det \mathrm{ectord}}\mathrm{\η}{L}_v\end{array}$

式中：j为节点i上的光源编号，d为探测器编号，v为需要重建的模型体素编号，J为荧光雅可比矩阵，为对应体素v，光源j，探测器d的荧光雅可比矩阵元。

在每个节点上，依据计算得到的荧光雅可比矩阵块，计算共轭梯度法所需要的中间量，并将这些中间量传输到主节点上，这种中间量要比荧光雅可比矩阵块占的内存小，因此传输更快，J_i和D_i计算公式如下：

D_i＝[D₁₁…D_1t；…；D_s1…D_st]^T

式中：i为节点编号，t为探测器总数，n为需要重建的模型体素数目，D_i为节点i上所有光源下的探测器探测到的荧光强度分布与pfMC模拟得到的荧光分布的差值，D_st为节点i上第s光源下的探测器中第t个探测单元探测到的荧光强度分布与相对应的dfMC模拟得到的荧光分布的差值。

3.3在主节点上，由上述的中间量J_iJ_i ^TX与J_iD_i进行共轭梯度法迭代，求解荧光团位置及其吸收系数,具体计算公式如下：

$\underset{X}{\min }f\left(X\right)=\left|\right|J^T\mathrm{JX}-J^{T}D\left|\right|=\left|\right|\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(J_i{J_i}^{T}X-J_i{D}_i)\left|\right|$

式中：i为节点编号，N为总节点数目，X为荧光团吸收系数分布矩阵，J_i ^T为矩阵J_i的转置矩阵，f(X)为目标函数。

下面以一个简单的仿真体模实验来验证本发明实施例提供的基于微扰荧光蒙特卡罗模型的荧光扩散光学断层成像方法的可行性，详见下文描述：

如图3a所示，设计了一个圆柱形仿真模型，在模型内假定填充有脂肪水溶液，并垂直插入1根细荧光棒以及荧光球，参数如表1所示。在MC模拟中，圆柱形仿真模型被剖分成161X150X161个体素。在不同光源位置分别进行MC模拟，并记录下所有到达探测器的光子的路径信息及出射组织时的状态信息(权重、方向、位置)及体素内的路径相关物理量。这里我们模拟的到达探测器的光子数为300万。同时我们使用荧光MC模拟直接得到探测器上所需的实验数据，并且加上5％的高斯随机噪声。根据模拟光子的路径信息及权重信息可以计算出雅克比矩阵，仿真模型的信息如下表所示：

表1：组织光学参数

这里我们对模型进行360o全空间扫描了，共计480个光源位置，如图3a、图3b所示。同时总共采集了70000多个有效探测器数据。在CPU/GPU集群框架下，进行白蒙特卡罗模拟，得到了大量光子路径相关物理量，并通过pMC方法计算得到荧光雅可比矩阵，最后采用共轭梯度法迭代重建荧光系数分布。如图3a、图3b所示是迭代20次的重建结果。

Claims (7)

1.一种基于pfMC模型的荧光扩散光学断层图像重建方法，其特征在于包括下列步骤：

(1)确定探测区域，选择光源扫描位置，并获取相对应的探测器上的荧光强度分布；

(2)数字化待探测对象，建立描述组织光学参数三维分布的体素化模型，在不同光源位置与方向进行白蒙特卡罗模拟，记录在探测器收集到的激发光子的权重及在每个体素内的路径信息和相应的物理量；

(3)依据pfMC模型，计算荧光光子权重；

(4)由激发光子的路径信息构建荧光雅可比矩阵；

(5)根据计算得到荧光光子权重和荧光雅可比矩阵，迭代重建计算出组织内的荧光团位置及荧光系数。

2.根据权利要求1所述的基于pfMC模型的荧光扩散光学断层图像重建方法，其特征在于步骤(2)中，采用GPU集群加速整个白蒙特卡罗过程，记录被探测器采集到的激发光子的相关物理量，被记录在体素v内的路径相关物理量L_v：为：

$L_v=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\tilde{l}}_i=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{l}_i\·e^{{\mathrm{\μ}}_{\mathrm{av}}{l}_i}$

式中：n为体素v内的散射次数，μ_av为体素v内的吸收系数，l_i为体素v内的第i次散射路径长度。

3.根据权利要求2所述的基于pfMC模型的荧光扩散光学断层图像重建方法，其特征在于步骤(2)中GPU集群加速整个白蒙特卡罗过程具体如下：

2.1在集群上，依据光源数目及集群节点数目分配各个节点上进行MC模拟的次数，并在集群启动后把相应的光源位置和方向信息输送到各个节点上；

2.2在各个节点上，根据要模拟的光子数目分配内存，将组织模型信息分配到GPU内的纹理寄存器，光学参数信息分配到GPU内的常数寄存器来提高读取速度，并根据GPU的技术信息确定最优的CUDA运行中并行块和线程的数目；

2.3在GPU内的各个线程上，每个线程上模拟一个光子包，并跟踪光子发射后在组织内散射及吸收的过程，记录光子在每个体素内的路径相关物理量，直到光子溢出组织后终止光子模拟。

4.根据权利要求1所述的基于pfMC模型的荧光扩散光学断层图像重建方法，其特征在于步骤(3)，计算荧光权重w^em(p)计算公式如下：

$w^{\mathrm{em}}\left(p\right)=w^{\mathrm{ex}}\left(p\right)\mathrm{\η}\underset{v=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\μ}}_{\mathrm{afv}}{L}_v$

式中：p为位置及单位方向的六维矢量，k为光子在组织内的路径所经过的体素数目，μ_afv是体素v内荧光系数，η是量子产率，w^ex为激发光子权重。

5.根据权利要求1所述的基于pfMC模型的荧光扩散光学断层图像重建方法，其特征在于步骤(4)中，荧光雅可比矩阵的计算公式如下：

$\begin{array}{c}{J}_{\mathrm{jd}}^v=\frac{{{{\∂}_w}^{\mathrm{em}}\left(p\right)|}_{\det \mathrm{e\; ctor\; d}}}{\∂{\mathrm{\μ}}_{\mathrm{afv}}}\\ =w^{\mathrm{ex}}{\left(p\right)|}_{\det \mathrm{ector\; d}}\mathrm{\η}{L}_v\end{array}$

式中：j为节点i上的光源编号，d为探测器编号，v为需要重建的模型体素编号，J为荧光雅可比矩阵，为对应体素v、光源j、探测器d的荧光雅可比矩阵元。

6.根据权利要求1所述的基于pfMC模型的快速荧光扩散光学断层图像重建方法，其特征在于步骤(5)，整个荧光重建过程均通过GPU集群加速，对各个节点上计算得到的荧光雅可比矩阵矩阵块进行处理，节点i上矩阵块计算公式如下：

D_i＝[D₁₁…D_1t；…；D_s1…D_st]^T

式中：s为节点i上的光源数目，t为探测器总数，n为需要重建的模型体素数目，D_i为节点i上所有光源下的探测器探测到的荧光强度分布与pfMC模拟得到的荧光分布的差值，D_st为节点i上第s光源下的探测器中第t个探测单元探测到的荧光强度分布与相对应的dfMC模拟得到的荧光分布的差值。

7.根据权利要求6所述的基于pfMC模型的快速荧光扩散光学断层图像重建方法，其特征在于步骤(5)，矩阵块计算的中间量J_iJ_i ^TX与J_iD_i被传输到主节上再进行共轭梯度迭代求解荧光分布，这些中间量比荧光雅可比矩阵块占用的存贮空间少，迭代计算公式如下：

$\underset{X}{\min }f\left(X\right)=\left|\right|J^T\mathrm{JX}-J^{T}D\left|\right|=\left|\right|\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(J_i{J_i}^{T}X-J_i{D}_i)\left|\right|$

式中：i为节点编号，N为总节点数目，X为荧光团吸收系数分布矩阵，J_i ^T为矩阵J_i的转置矩阵，f(X)为目标函数。