CN105559750A - 组织结构引导的复合正则化生物发光断层成像重建方法 - Google Patents

Abstract

组织结构引导的复合正则化生物发光断层成像重建方法，本方法属于医学图像处理领域；针对现有生物发光断层成像重建算法中存在的上述问题，本方法可以在已知肿瘤所在的生物组织而不知具体位置时很好的实现光源的重建；另外在不知肿瘤所在生物组织时，通过均分权重系数来实现重建。首先采用有限元方法仿真光在生物组织中的传输，并基于生物发光荧光断层成像中光源的稀疏分布特性，采用结构成像引导的复合正则化方法融入更多先验信息以降低BLT重建的病态性，最后在此基础上实现荧光光源的准确重建。在两个波段设置了双光源进行仿真实验，以验证算法的有效性和准确性。本方法不但可以对荧光光源进行准确重建定位，而且可以极大地提高计算效率。

Description

组织结构引导的复合正则化生物发光断层成像重建方法

技术领域

本发明属于医学图像处理领域，涉及一种组织结构引导的复合正则化生物发光断层成像重建算法。

背景技术

光学分子影像是一种快速发展的分子影像技术，它通过将光学过程与一定的分子性质相结合，对目标体内的生物荧光或激发荧光进行分析和处理，并进行定性和定量研究。

光学分子影像技术中最具代表性的成像方式是激发荧光成像技术(FluorescenceImaging)和生物发光成像技术(BioluminescenceImaging,BLI)。它们均是二维的生物发光成像技术，虽然这种成像技术应用起来方便、简单，但这种二维成像方式在应用过程中存在局限性，尤其是对于成像深度的限制，二维的荧光图像不能反映光源深度信息,并且难以定量化。它们仅可以反映生物体内的荧光探针在某一角度的投影信息，而这种投影信息是多个深度的信号叠加而成的。因此，二维成像方法具有很低的分辨率。

在二维生物发光成像技术基础上发展起来的生物发光断层成像技术(BioluminescenceTomography,BLT)由于能够反映信号深度的信息，已经成为光学分子影像的一个重要分支。生物发光断层成像技术不需要外在光源的激发，而是通过一种生物化学发光反应在体内发光。体内产生的荧光在生物组织内部以某种规律传播并不断地与生物组织发生相互作用，并到达体表。最后，在生物组织体表利用高灵敏度的探测器获得的荧光图像就可以重建出荧光光源在小动物体内的分布情况，从而在本质上揭示在体分子的活动规律。

但光子在生物组织中不沿直线传输，而是经历了大量的散射过程，导致BLT逆问题在数学上是一个高度病态的问题，外界微小的测量扰动都会给重建结果带来很大的变化。国内外研究人员为降低其病态性做了很多工作，一般都是从不同角度为该问题提供各种不同的先验知识。

为了克服重建的病态性使重建变得更加稳定，传统的求解方法一般是建立在基于L2范数正则化项的基础上的。但是，使用基于L2范数正则化项会不可避免地使重建结果过于平滑，容易陷入局部最小值，从而降低成像精度。为了较好地解决基于L2范数的正则化方法得到的光源重建结果容易平滑的问题，需要在重建中加入更多的合理的先验信息。作为一种典型的特异性成像，生物发光断层成像可以结合一些结构成像方式，如micro-CT、micro-MRI等进行重建。一般的重建方法中，只是对不同生物组织赋以不同的光学特性参数，结构成像提供的先验信息没有得到充分利用。本发明提供的重建方法可以充分利用结构成像提供的先验信息，当已知肿瘤所在组织而不确定其具体位置时，通过对组织权重赋值的方法来确定其具体位置，使得重建结果更加准确。另外，当不能确定肿瘤所处组织时，组织权重可能需要均分。本发明非常适合用来进行肿瘤的早期检测，并且需要成像的肿瘤分布是稀疏的，因此可以将光源重建问题描述为稀疏约束问题，采用稀疏正则化方法重建。

综合上述问题，本发明提出一种结构成像引导的复合正则化方法来准确重建荧光光源。

发明内容

针对现有生物发光断层成像重建算法中存在的上述问题，本发明提出了一种结构成像引导的复合正则化重建方法，该方法可以在已知肿瘤所在的生物组织而不知具体位置时很好的实现光源的重建；另外在不知肿瘤所在生物组织时，通过均分权重系数来实现重建。

本发明采用的技术方案为：首先采用有限元方法仿真光在生物组织中的传输，并基于生物发光荧光断层成像中光源的稀疏分布特性，采用结构成像引导的复合正则化方法融入更多先验信息以降低BLT重建的病态性，最后在此基础上实现荧光光源的准确重建。

具体包括下述步骤：

步骤一，光在生物组织中的传输规律用扩散方程加以表示。根据变分原理，结合Robin边界条件，建立与扩散近似方程和边值问题等价的弱形式：

${\∫}_{\mathrm{\Ω}}\left(\right(D\left(r\right)\left(\▿\mathrm{\Φ}\left(r\right)\right)\·\▿\mathrm{\Ψ}\left(r\right))+{\mathrm{\μ}}_a\mathrm{\Φ}(r\left)\mathrm{\Ψ}\right(r\left)\right)dr+{\∫}_{\∂\mathrm{\Ω}}\mathrm{\Φ}\left(r\right)\mathrm{\Ψ}\left(r\right)/\left(2I_n\right(r\left)\right)dr{\∫}_{\mathrm{\Ω}}S\left(r\right)\mathrm{\Ψ}\left(r\right)dr---\left(1\right)$ 其中D(r)为扩散系数分布，Φ(r)为光子密度分布，μ_a为光吸收系数，Ψ(r)为任意测试函数，S(r)为光源分布，I_n(r)为使用Robin边界条件时引入的一个量，R_n≈-1.4399n^-2+0.7099n^-1+0.6681+0.636n表示扩算传输内反射系数，n是一个常数，与边界内外的光学折射系数偏差有关。

步骤二，有限元离散。对成像体Ω进行网格剖分，实现离散化。本方法采用二维仿体进行仿真，由于三角网格单元能够实现使用较少数量来近似描述不规则区域，所以，本方法使用三角网格剖分成像物体。根据对近似解精度的要求、单元的几何形状、节点个数和节点的自由度信息来选择单元基函数。由于函数均表示为基函数的线性组合，假设成像物体所划分的网格包含Nn个节点和Nt个三角单元，光子密度分布Φ(r)表示为：

其中，φ_k(r)表示节点i上的函数值，表示节点i处的基函数。同样，光源函数S(r)表示为：

其中，s_k(r)表示节点i上的光源值。

步骤三，合成总体刚度矩阵。将(2)、(3)式代入(1)式中，得到：

该矩阵方程简单表示为：

([K]+[C]+[B])＝MΦ＝FS(5)

其中，M为刚度矩阵，其他各矩阵元素通过下面公式计算得到：

整理(6)式，得到：

M^-1FS＝Φ(7)

进一步地，上式写成如下形式：

AS＝Φ(8)

其中，A＝M^-1F。

对于生物发光断层成像，探测器所能够采集到的测量数据仅仅是成像物体表面的荧光信息，因此在(7)式中，需要将光子密度向量Φ中的非表面元素去掉，并将系统矩阵A中对应的行去掉，从而得到不含表面荧光信息的光子密度分布Φ^m和系统矩阵A，此时的光源重建问题表述为：

AS＝Φ^m(9)

步骤四，利用结构成像得到组织结构。利用micro-CT或micro-MRI等结构成像提供的组织结构信息，通过对生物组织进行分割，得到不同的生物组织，将此作为先验信息融合到后续的重建中。

步骤五，采用复合正则化方法进行重建。目标函数设置如下(10)式所示：

$min\left\{f\left(s\right)\right\}=\left|\right|As-{\mathrm{\Φ}}^m|{|}_2^2+\mathrm{\λ}|\left|s\right|{|}_1+\mathrm{\β}\underset{i=1}{\overset{q}{\Σ}}{\mathrm{\ω}}_i\left|\right|s_G|{|}_2---\left(10\right)$

上式参数的意义如下：

A：利用有限元求解扩散方程得到的系统矩阵；

s：是未知量，待求解；

Φ^m：测量值；

λ：正则化参数1；

β：正则化参数2；

ω_i：不同组织的权重；

s_G:表示其余未知量s分解为的q个不相重合的部分，q可表示将结构图像分割后分成的q个不同组织，ω_i表示对不同组织的权重。

通过求上述最优化问题，可重建出未知量s的分布情况。

附图说明

图1为波长为600nm时的一个两区域仿体结构光学参数分布示意图；

图2为波长为630nm时的一个两区域仿体结构光学参数分布示意图；

图3为仿体初始的光源分布，图中面积较大的黑色圆形区域为仿体区域，黑色区域中面积较小的圆形区域为光源分布区域；

图4为通过本方法重建的光源分布结果。

具体实施方式

下面结合附图和实施例对本发明做更详细的说明。

为了检验本发明提出方法的有效性，通过仿体实验进行验证。在实验中，考虑到光谱分布，选择的两个谱段分别为600nm和630nm，并认为生物组织的光学特性参数已知，详细参数如表1所示,具体结构图如图1和图2所示。

表1生物组织光学特性参数

本发明采用有限元方法产生前向边界测量数据。为有效地避免“逆行为”，在产生前向数据时，前向数据的网格含有3508个节点，重建时所需要的网格含有1309个节点。

本发明中，在两个波段设置了双光源进行仿真实验，以验证算法的有效性和准确性。实验结果表明，本算法不但可以对荧光光源进行准确重建定位，而且可以极大地提高计算效率。

Claims (1)

1.组织结构引导的复合正则化生物发光断层成像重建方法，其特征在于：

首先采用有限元方法仿真光在生物组织中的传输，并基于生物发光荧光断层成像中光源的稀疏分布特性，采用结构成像引导的复合正则化方法融入更多先验信息以降低BLT重建的病态性，最后在此基础上实现荧光光源的准确重建；

具体包括下述步骤：

步骤一，光在生物组织中的传输规律用扩散方程加以表示；根据变分原理，结合Robin边界条件，建立与扩散近似方程和边值问题等价的弱形式：

${\∫}_{\mathrm{\Ω}}((D\left(r\right)(\▿\mathrm{\Φ}\left(r\right))\·\▿\mathrm{\Ψ}\left(r\right))+{\mathrm{\μ}}_a\mathrm{\Φ}\left(r\right)\mathrm{\Ψ}\left(r\right))\mathrm{dr}+{\∫}_{\∂\mathrm{\Ω}}\mathrm{\Φ}\left(r\right)\mathrm{\Ψ}\left(r\right)/\left(2I_n\left(r\right)\right)\mathrm{dr}={\∫}_{\mathrm{\Ω}}S\left(r\right)\mathrm{\Ψ}\left(r\right)\text{dr---}\left(1\right)$

其中D(r)为扩散系数分布，Φ(r)为光子密度分布，μ_a为光吸收系数，Ψ(r)为任意测试函数，S(r)为光源分布，I_n(r)为使用Robin边界条件时引入的一个量，R_n≈-1.4399n^-2+0.7099n^-1+0.6681+0.636n表示扩算传输内反射系数，n是一个常数，与边界内外的光学折射系数偏差有关；

步骤二，有限元离散；对成像体Ω进行网格剖分，实现离散化；本方法采用二维仿体进行仿真，由于三角网格单元能够实现使用较少数量来近似描述不规则区域，所以，本方法使用三角网格剖分成像物体；根据对近似解精度的要求、单元的几何形状、节点个数和节点的自由度信息来选择单元基函数；由于函数均表示为基函数的线性组合，假设成像物体所划分的网格包含N_n个节点和N_t个三角单元，光子密度分布Φ(r)表示为：

其中，φ_k(r)表示节点i上的函数值，表示节点i处的基函数；同样，光源函数S(r)表示为：

其中，s_k(r)表示节点i上的光源值；

步骤三，合成总体刚度矩阵；将(2)、(3)式代入(1)式中，得到：

该矩阵方程简单表示为：

([K]+[C]+[B])＝MΦ＝FS(5)

其中，M为刚度矩阵，其他各矩阵元素通过下面公式计算得到：

整理(6)式，得到：

M^-1FS＝Φ(7)

进一步地，上式写成如下形式：

AS＝Φ(8)

其中，A＝M^-1F；

对于生物发光断层成像，探测器所能够采集到的测量数据仅仅是成像物体表面的荧光信息，因此在(7)式中，需要将光子密度向量Φ中的非表面元素去掉，并将系统矩阵A中对应的行去掉，从而得到不含表面荧光信息的光子密度分布Φ^m和系统矩阵A，此时的光源重建问题表述为：

AS＝Φ^m(9)

步骤四，利用结构成像得到组织结构；利用micro-CT或micro-MRI等结构成像提供的组织结构信息，通过对生物组织进行分割，得到不同的生物组织，将此作为先验信息融合到后续的重建中；

步骤五，采用复合正则化方法进行重建；目标函数设置如下(10)式所示：

$min\left\{f\left(s\right)\right\}=\left|\right|As-{\mathrm{\Φ}}^m|{|}_2^2+\mathrm{\λ}|\left|s\right|{|}_1+\mathrm{\β}\underset{i=1}{\overset{q}{\Σ}}{\mathrm{\ω}}_i\left|\right|s_G|{|}_2---\left(10\right)$

上式参数的意义如下：

A：利用有限元求解扩散方程得到的系统矩阵；

s：是未知量，待求解；

Φ^m：测量值；

λ：正则化参数1；

β：正则化参数2；

ω_i：不同组织的权重；

s_G:表示其余未知量s分解为的q个不相重合的部分，q可表示将结构图像分割后分成的q个不同组织，ω_i表示对不同组织的权重；

通过求上述最优化问题，可重建出未知量s的分布情况。