CN105581779B - 一种直接融合结构成像的生物发光断层成像重建的方法 - Google Patents

Abstract

一种直接融合结构成像的生物发光断层成像重建的方法，属于医学图像处理领域，本方法针对现有的生物发光断层成像重建算法中存在的问题，提出了一种融合其他成像模态进行BLT重建的方法，它利用结构成像方式获得的重建图像作为先验信息计算得到自适应的正则化矩阵，进而实现BLT重建。首先基于扩散方程模拟光在生物组织中的传输规律，利用有限元方法进行数值求解；为实现融合结构成像进行BLT重建，并考虑到荧光光源的稀疏特性，提出了一种不需对结构图像组织分割便可融合的方法。为实现对上述方法的求解，提出了基于Split‑Bregman的迭代方法进行求解。

Description

一种直接融合结构成像的生物发光断层成像重建的方法

技术领域

本发明属于医学图像处理领域，涉及一种直接融合结构成像的生物发光断层成像重建算法。

背景技术

随着人类对生命奥秘的不断探索以及现代计算机技术的飞速发展，各种医学影像技术和医学成像设备也进入了快速发展期。从传统的计算机断层成像(ComputedTomography，CT)、超生成像(Ultrasound Imaging，USI)、核磁共振成像(MagneticResonance Imaging，MRI)、正电子成像(Positron Emission Tomography，PET)、单光子发射断层成像(Single-Photon Emission Computed Tomography，SPECT)发展到激发荧光断层成像(Fluorescent Molecular Tomography，FMT)和生物发光断层成像(Bioluminescence Tomography，BLT)等多种成像模态，这些成像技术极大推动了生命科学的进步。

生物发光断层成像(以下简称BLT)作为一种重要的成像模态，利用其高灵敏度、背景噪声低等特性，可以在生物体表获得荧光图像重建生物体内荧光源的分布以及细胞分子水平变化。BLT由于能够反映信号深度的信息，已经成为光学分子影像的一个重要分支。生物发光断层成像技术不需要外在光源的激发，而是通过一种生物化学发光反应在体内发光。体内产生的荧光在生物组织内部以某种规律传播并不断地与生物组织发生相互作用，并到达体表。最后，在生物组织体表利用高灵敏度的探测器获得的荧光图像就可以重建出荧光光源在小动物体内的分布情况，从而在本质上揭示在体分子的活动规律。

然而生物体组织具有强散射、弱吸收的特点，荧光光子无法在生物组织中沿直线传输，而是经历了大量的散射过程，导致BLT逆问题在数学上是一个高度病态的问题，外界微小的测量扰动都会给重建结果带来很大的变化。因此，降低BLT的病态性，如何唯一、准确地重建荧光光源是BLT研究的重点及热点。

近年来，国内外研究学者为了解决BLT成像的准确性和可靠性等难题做了各种尝试。国内外研究人员为降低其病态性做了很多工作，一般都是从不同角度为该问题提供各种不同的先验知识。其中，融合结构成像的方法能够起到很好的BLT重建效果。但在融合结构成像时往往需要人为对结构成像图片进行基于不同生物组织的分割，这样往往会消耗过多的时间，同时人为带有主观性的分割方式也会为BLT重建带来不可预知的误差。

在BLT重建时大都是基于正则化方法将光源重建问题转变成一个非线性的最优化问题。正则化方法最早在上世纪60年代由苏联科学院院士Tikhonov提出，用来求解形如Fx＝y方程的不适定问题。正则化的核心思想是在求解问题过程中引入先验信息，从而得到对原问题有意义的解。对于基于正则化的重建方法，正则化参数对结果有较大影响。参数的选取是否适当，不仅影响着BLT重建算法的收敛性以及收敛速度，而且影响着所求得的解是否收敛于问题的真实解。同时，正则化项中的正则化矩阵的选择也会影响到解。在传统的重建模型中，正则化矩阵一般被设定为单位矩阵，其能够为成像区域内的所有节点提供相同的权重。如果能够基于其他成像模态所得到的先验信息求得自适应的正则化矩阵，就可以为成像区域内的点赋予不同权重，从而提升重建效果。现有的一些基于正则化的BLT重建方法也会融合结构成像作为先验信息，而这些方法又分为硬先验与软先验。硬先验方法是将同一解剖结构的光学特性参数看作是一个重建参数以大量减少重建参数的数量，其结果高度地依赖于先验信息即结构成像的准确性。软先验方法则是通过在目标函数的惩罚项中引入正则化矩阵来实现的，而正则矩阵是根据先验信息构造的，由于正则化系数可调整，这使得软先验具有更强的适应性。

本发明提出一种直接融合结构成像进行BLT重建的方法。首先基于扩散方程模拟光在生物组织中的传输规律，利用有限元方法进行数值求解；为实现融合结构成像进行BLT重建，并考虑到荧光光源的稀疏特性，提出了一种不需对结构图像组织分割便可融合的方法。为实现对上述方法的求解，提出了基于Split-Bregman的迭代方法进行求解。

发明内容

针对现有的生物发光断层成像重建算法中存在的问题，本发明提出了一种融合其他成像模态进行BLT重建的方法，它利用结构成像方式获得的重建图像作为先验信息计算得到自适应的正则化矩阵，进而实现BLT重建。

本发明采用的技术方案为：基于扩散方程并利用有限元方法及正则化思想将BLT重建问题转为最小二乘问题得到表达式，融合结构成像计算得到正则化矩阵，通过Split-Bregman的迭代方法求解方程，得到BLT重建结果。

在BLT实验中，认为生物发光的光源是稳定的，即工作在连续波的状态。由于BLT的实验是在全黑密闭的暗室下进行的，因此用稳态的扩散方程和Robin边界条件来描述光在生物组织中的传播过程，然而在光源位置较深或光源强度较弱的情况下，使用单光谱测量值很难反演出光源的正确分布，因此，引入多光谱进行光源重建，将扩散近似方程和Robin边界条件进行改变，增加λ参数来描述谱段。即有如下方程：

其中，r表示坐标向量，Ω表示生物组织区域，代表组织边界，D(r,λ)和μ_a(r,λ)分别表示与波长λ有关的扩散系数和吸收系数，Φ(r,λ)和S(r,λ)分别表示波长为λ时的光子密度和光源密度，v(r)表示边界上的单位外法线，A(r；n,n′)表示为

其中n,n′表示区域内、外的折射率，R(r)由n计算得出，R(r)≈-1.4399n^-2+0.7099n^-1+0.6881+0.0636n，定义荧光图像是在一组互不连接的表面片段Υ_i上进行采集，片段Υ_i光滑且连续，并满足那么得边界测量的光流密度Q(r,λ)表示为：

在实验中，使用滤波片将到达生物体体表的光分成M个波段τ₁,…,τ_M，τ_i＝[λ_i-1,λ_i)。基于有限元理论，得到波段τ_i上的Φ(r,τ_i)∈H¹(Ω)的弱解形式：

Ψ(r,τ_i)为从一个试探空间中选取的试探函数，选取要求为在Ω上，该函数本身和其所有的一阶偏导数可积，并且它的边界条件与Φ(r,τ_i)的边界条件相同。

下面对生物组织求解区域进行离散化，根据有限元方法，将Ω离散为含N_T个单元和N_P个节点。基于有限元方法，将公式(5)在单个谱段τ_l整理为如下方程：

(K(τ_l)+C(τ_l)+B(τ_l))Φ(τ_l)＝M(τ_l)Φ(τ_l)＝F(τ_l)S(τ_l) (6)

矩阵K、C、B和F的元素分别为：

公式(6)和(7)建立了未知光源密度和边界测量值之间的关系。由公式(6)可知，M(τ_l)＝K(τ_l)+C(τ_l)+B(τ_l)，由于M(τ_l)是稀疏的正定矩阵，因此有：

Φ(τ_l)＝M^-1(τ_l)F(τ_l)S(τ_l) (8)

将上式进行简化整理，边界测量值和未知光源密度之间的关系即表示为：

Φ(τ_l)＝A(τ_l)S(τ_l) (9)

在多光谱情况下，每个波段τ_l的能量比ω(τ_l)通过预先的谱分析测得，即有S(τ_l)＝ω(τ_l)S，S代表总光子密度。将所有谱段的光源及边界测量值进行整合，可得：

Φ＝AS (10)

其中，A和Φ可表示为：

基于正则化思想，通常将该问题转为最小二乘问题，即表现为如下形式：

其中，L表示正则化矩阵。在传统的模型中，正则化矩阵一般设定为单位矩阵，在本方法中L通过结构成像方法获得，具体表达式如下：

其中γ_i、γ_j为其他图像的像素归一化灰度值；r_i、r_j表示有限元节点的坐标；M_i满足(n表示有限元节点数)；σ_g和σ_d为尺度因子，其中σ_g表示灰度特征差异，σ_d用于衡量节点的距离差异；θ(x)表示阶跃函数：

本方法采用Split-Bregman迭代方法对公式(13)进行求解，首先‖Lx‖₁可以被转换为基于有限元的矩阵形式|Lx|，实用u来表示Lx，将问题转换为如下形式：

通过添加惩罚项的方式将公式(16)变为一个无约束问题，如下式：

其中μ是惩罚参数。最终，通过Split-Bregman迭代转换为两个子问题：

b_k+1＝b_k+u_k+1-Lx_k+1 (18)

为了使问题的求解更为简化，将第一个子问题拆分得到两个子问题，这样就得到了如下三个步骤：

步骤3：b_k+1＝b_k+u_k+1-Lx_k+1 (19)

这样，通过Split-Bregman迭代方法即可解得x。

本方法针对现有的生物发光断层成像重建算法中存在的问题，提出一种直接融合其他结构成像模态进行BLT重建的方法。首先基于扩散方程模拟光在生物组织中的传输规律，并利用有限元方法进行数值求解；考虑到荧光光源的稀疏特性，提出了一种不需对结构图像组织分割便可融合的方法以实现融合结构成像进行BLT重建。为实现对上述方法的求解，提出了基于Split-Bregman的迭代方法进行求解，最终得到重建结果。

附图说明

图1为实验用到的两区域仿体示意图，黑色区域中面积较小的白色圆形区域(异常区域)与黑色区域(正常区域)的光吸收与散射系数不同；

图2为一个两区域的仿体结构成像示意图，图中区域1的灰度值被设为80，区域2的灰度值为50；

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

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

具体实施方式

下面根据实施示例与附图对本发明进行说明。

首先，通过nirfast建立仿体。本实验用到的是一个两区域仿体，仿体结构如图1所示。仿体的前向有限元节点数设置为3508个。由于BLT逆问题求解中基于单光谱的重建结果可能存在不唯一性，因此实验中采用了两个波段。不同光波段在仿体中拥有不同的光吸收系数与光散射系数，实验所用到的波段在仿体中的光学特性参数如表1所示。

表1不同波段的仿体光学特性参数

通过添加光源得到的仿体的初始光源分布如图3所示。

一般情况下，为了保证重建时计算量不会过大，重建模型的有限元节点数往往小于前向的有限元节点数。在本实验中，重建的有限元节点数设为1309个。计算得到重建的有限元模型边界结点并与前向的节点进行对应，根据仿体的初始光源分布通过计算得到边界光子数密度测量值Φ。

根据如图2所示两区域的仿体结构成像示意图计算出正则化矩阵L，计算公式如公式(14)所示。在本实例中尺度因子σ_g和σ_d的取值分别设定为0.001与0.1。

根据已经得到的Φ，L与通过计算得到的A，求解公式(13)即可得到BLT重建结果。在求解公式(13)时，本发明使用公式(19)所示的Split-Bregman迭代方法，计算出如图4所示的光源重建结果。

Claims (1)

1.一种直接融合结构成像的生物发光断层成像重建的方法，其特征在于：

该方法利用结构成像方式获得的重建图像作为先验信息计算得到自适应的正则化矩阵，进而实现BLT重建；

基于扩散方程并利用有限元方法及正则化思想将BLT重建问题转为最小二乘问题得到表达式，融合结构成像计算得到正则化矩阵，通过Split-Bregman的迭代方法求解方程，得到BLT重建结果；

在BLT实验中，认为生物发光的光源是稳定的，即工作在连续波的状态；由于BLT的实验是在全黑密闭的暗室下进行的，因此用稳态的扩散方程和Robin边界条件来描述光在生物组织中的传播过程，然而在光源位置较深或光源强度较弱的情况下，使用单光谱测量值很难反演出光源的正确分布，因此，引入多光谱进行光源重建，将扩散近似方程和Robin边界条件进行改变，增加λ参数来描述谱段；即有如下方程：

其中，r表示坐标向量，Ω表示生物组织区域，代表组织边界，D(r,λ)和μ_a(r,λ)分别表示与波长λ有关的扩散系数和吸收系数，Φ(r,λ)和S(r,λ)分别表示波长为λ时的光子密度和光源密度，v(r)表示边界上的单位外法线，A(r；n,n′)表示为

其中n,n′表示区域内、外的折射率，R(r)由n计算得出，R(r)≈-1.4399n^-2+0.7099n^-1+0.6881+0.0636n，定义荧光图像是在一组互不连接的表面片段Υ_i上进行采集，片段Υ_i光滑且连续，并满足那么得边界测量的光流密度Q(r,λ)表示为：

在实验中，使用滤波片将到达生物体体表的光分成M个波段τ₁,…,τ_M，τ_i＝[λ_i-1,λ_i)；基于有限元理论，得到波段τ_i上的Φ(r,τ_i)∈H¹(Ω)的弱解形式：

Ψ(r,τ_i)为从一个试探空间中选取的试探函数，选取要求为在Ω上，该函数本身和其所有的一阶偏导数可积，并且它的边界条件与Φ(r,τ_i)的边界条件相同；

下面对生物组织求解区域进行离散化，根据有限元方法，将Ω离散为含N_T个单元和N_P个节点；基于有限元方法，将公式(5)在单个谱段τ_l整理为如下方程：

(K(τ_l)+C(τ_l)+B(τ_l))Φ(τ_l)＝M(τ_l)Φ(τ_l)＝F(τ_l)S(τ_l) (6)

矩阵K、C、B和F的元素分别为：

公式(6)和(7)建立了未知光源密度和边界测量值之间的关系；由公式(6)可知，M(τ_l)＝K(τ_l)+C(τ_l)+B(τ_l)，由于M(τ_l)是稀疏的正定矩阵，因此有：

Φ(τ_l)＝M^-1(τ_l)F(τ_l)S(T_l) (8)

将上式进行简化整理，边界测量值和未知光源密度之间的关系即表示为：

Φ(τ_l)＝A(τ_l)S(τ_l) (9)

在多光谱情况下，每个波段τ_l的能量比ω(τ_l)通过预先的谱分析测得，即有S(τ_l)＝ω(τ_l)S，S代表总光子密度；将所有谱段的光源及边界测量值进行整合，可得：

Φ＝AS (10)

其中，A和Φ可表示为：

基于正则化思想，通常将该问题转为最小二乘问题，即表现为如下形式：

其中，φ(S,Φ)代表数据拟合项，λ为正则化参数，代表正则化项；在本题中， 将S由变量x表示，得到：

其中，L表示正则化矩阵；在传统的模型中，正则化矩阵一般设定为单位矩阵，在本方法中L通过结构成像方法获得，具体表达式如下：

其中γ_i、γ_j为其他图像的像素归一化灰度值；r_i、r_j表示有限元节点的坐标；M_i满足(n表示有限元节点数)；σ_g和σ_d为尺度因子，其中σ_g表示灰度特征差异，σ_d用于衡量节点的距离差异；θ(x)表示阶跃函数：

本方法采用Split-Bregman迭代方法对公式(13)进行求解，首先‖Lx‖₁可以被转换为基于有限元的矩阵形式|Lx|，使 用u来表示Lx，将问题转换为如下形式：

通过添加惩罚项的方式将公式(16)变为一个无约束问题，如下式：

其中μ是惩罚参数；最终，通过Split-Bregman迭代转换为两个子问题：

为了使问题的求解更为简化，将第一个子问题拆分得到两个子问题，这样就得到了如下三个步骤：

步骤1：

步骤2：

步骤3：b_k+1＝b_k+u_k+1-Lx_k+1 (19)

这样，通过Split-Bregman迭代方法即可解得x。