CN104586366B - 一种激发荧光断层成像重建方法 - Google Patents

Abstract

本发明公开了一种激发荧光断层成像重建方法，包括以下步骤：S1：基于光传输模型和有限元理论，将成像对象的光学特性参数和解剖结构信息作为先验信息，建立表面的测量数据与成像对象内部荧光目标分布的系统方程；S2:对S1建立的系统方程进行预处理，使得预处理后的系统方程中的稀疏矩阵各列相互正交；S3：融合光源空间分布的稀疏特性，从S1中系统方程中选取若干列，构建新的系数矩阵和新的系统方程；S4:将S3中新的系统方程转化为有约束条件的极小化问题，采用联合代数迭代方法对新的系统方程进行求解，从而获得成像对象内部的荧光目标的三维分布与浓度。本发明的有益效果：融合了荧光目标空间分布的稀疏特性，增加了先验信息，可提高重建精度；对系统方程进行预处理，系统方程的系数矩阵各列相互正交，可加快算法收敛，提高重建速度。

Description

一种激发荧光断层成像重建方法

技术领域

本发明涉及激发荧光成像领域，具体涉及一种激发荧光断层成像重建方法。

背景技术

激发荧光断层成像是近年来一种新出现的光学分子影像技术，它利用荧光探针标记包括成像对象(一般为小鼠、裸鼠等生物活体)，在外部激发光源的照射下，荧光探针吸收外部激发光，产生发射光。在成像对象体外，利用高灵敏度的光学检测仪器，如CCD相机，探测透射到成像对象体外的光子。结合光子传输模型，利用数学方法，激发荧光断层成像对成像对象体内的荧光目标进行重建，从而获得荧光探针的位置和浓度信息，在基因病理学、肿瘤诊断学、细胞和分子生物学等学科领域具有广泛的应用前景。

激发荧光断层成像重建属于逆向问题，具有严重的病态性，其本质原因在于生物组织的高散射、地吸收特性，光子在成像对象内部不沿直线传输，而是经过大量的无规则可循的散射过程。同时，由于光学检测仪器探测到的是传输到成像对象边界点的荧光信号，数量有限，而重建区域内部点的数量巨大。由少量测量数据来求解大量未知数，这是一个不适定问题，具有严重的病态性，其解不唯一并且易受测量误差和噪声的影响。如何精确重建成像对象体内的荧光目标，是荧光分子断层成像研究的核心问题。联合代数迭代技术由于运算简单，在断层成像重建中广泛应用。现有基于联合代数迭代技术的激发荧光断层成像重建没有考虑荧光目标空间分布的稀疏特性，同时没有对系统方程进行预处理，重建精度和重建速度有待提高。

发明内容

本发明要解决的技术问题是对基于联合代数迭代技术的激发荧光断层成像重建方法进行改进，提供一种精确度高、运算速度快的激发荧光断层成像重建方法。

本发明的技术解决方案是：一种激发荧光断层成像重建方法，其中：包括以下步骤：

S1：基于光传输模型和有限元理论，将成像对象的光学特性参数和解剖结构信息作为先验信息，建立成像对象的表面测量数据与内部荧光目标分布的系统方程；

S2：对S1建立的系统方程进行预处理，使得预处理后的系统方程的系数矩阵各列相互正交；

S3：融合荧光目标空间分布的稀疏特性，从S2中的系统方程中选取若干列，构建新的系数矩阵和新的系统方程；

S4：加入荧光目标空间分布非负的约束条件，采用联合代数迭代技术对S3中新的系统方程进行求解，从而获得成像对象内部的荧光目标的三维分布。

更为具体的包括以下步骤：

包括以下步骤：

(1)获取测量数据；

(1a)激发光源对固定在电控旋转台上的成像对象进行多角度的透射式断层扫描；

(1b)使用光学检测仪器获得测量数据；

(2)获得成像对象的解剖结构信息以及光学特性参数；

(3)基于光传输模型和有限元理论，将成像对象的解剖结构信息和光学特性参数作为先验信息，建立表面的测量数据与成像对象内部荧光目标分布的系统方程y＝Ax；

(4)对建立的y＝Ax进行如下预处理：

(4a)对A^TA进行特征值分解，得到QΛQ^T＝A^TA，其中T表示转置运算；

(4b)计算权值矩阵W＝Λ^-1Q^TA^T,将权值矩阵分别乘以系统方程的左侧和右侧，得到z＝Bx，式中，z＝Wy＝Λ^-1Q^TA^Ty，B＝Q^T；

(5)从上述矩阵B中挑选k列元素，记为系数矩阵C，C由以下步骤得到：

(5a)初始化，余量r₀＝z，系数矩阵C₀为空矩阵，迭代次数t＝1；

(5b)计算向量α＝|B^Tr_t-1|，取α的最大元素对应的下标，记为m_t，取矩阵B的第m_t列，记为b_mt，b_mt与C_t-1合并，构造系数矩阵C_t，具体为

(5c)计算并将s_t中的小于0的元素设为0；

(5d)计算新的余量r_t＝z-C_ts_t，迭代次数加1，得到t＝t+1；

(5e)重覆步骤5b-5d直到迭代次数满足t＝k；

(6)以步骤(5)得到的s_k作为荧光目标分布初始值，即x₀＝s_k，利用联合代数迭代方法求解荧光目标分布，迭代更新公式为：

x_n＝x_n-1+UC^TV(z-Cx_n-1)

$x_n=\frac{1}{2}(x_n+\mathrm{abs}\left(x_n\right))$

(7)步骤(6)的迭代停止准则为|x_n-x_n-1|/x_n≤10^-6。

本发明的有益效果：

融合了荧光目标空间分布的稀疏特性，增加了先验信息，可提高重建精度；对系统方程进行预处理，系统方程各列相互正交，可加快算法收敛，提高重建速度。

附图说明

图1为本发明方法流程图。

图2为实验所用仿体模型。

图3为联合代数迭代方法得到的重建结果。

图4为本发明方法得到的重建结果。

具体实施方式

实施例：

参阅图1至图4；

步骤1，建立系统方程

根据有限元理论，建立系统方程，该系统方程为激发荧光断层成像对象体表的荧光测量数据与成像对象内部荧光目标分布的线性关系：

y＝Ax (1)

式中，y是已知的荧光测量数据向量，A是已知的系统矩阵，x是要求解的荧光目标分布向量，且x和y的元素均非负；

步骤2，求解系统方程

(1)对建立的系统方程y＝Ax进行如下预处理：

(1a)对A^TA进行特征值分解，得到QΛQ^T＝A^TA，其中T表示转置运算，Q为特征向量构成的矩阵，Λ为特征值构成的对角矩阵；

(1b)计算权值矩阵W＝Λ^-1Q^TA^T，将权值矩阵分别乘以系统方程的左侧和右侧，得到z＝Bx，式中，z＝Wy＝Λ^-1Q^TA^Ty，B＝Q^T；

(2)从上述矩阵B中挑选k列元素，记为系数矩阵C，k取值由荧光目标空间分布的稀疏度决定，C根据以下步骤得到：

(2a)初始化，余量r₀＝z，系数矩阵C₀为空矩阵，迭代次数t＝1；

(2b)计算向量α＝|B^Tr_t-1|，取α的最大元素对应的下标，记为m_t，取矩阵B的第m_t列，记为与C_t-1合并，构造系数矩阵C_t，具体为

(2c)计算并将s_t中的小于0的元素设为0；

(2d)计算新的余量r_t＝z-C_ts_t，迭代次数加1，得到t＝t+1；

(2e)重复步骤2b-2d直到迭代次数满足t＝k；

(3)利用步骤(1)和步骤(2)分布得到的z和C，构建新的系统方程并进行求解，得到成像对象体内的荧光目标分布，新的系统方程为z＝Cx，荧光目标重建可表示为如下最小化问题：

min||z-Cx|| s.t.x≥0 (2)

以步骤(2)得到的s_k作为荧光目标分布初始值x₀＝s_k，采用联合代数迭代方法求解上述最小化问题，迭代更新公式为：

$x_n=x_{n-1}+{\mathrm{UC}}_k^{T}V(z-{\mathrm{Cx}}_{n-1})---\left(3\right)$

$x_n=\frac{1}{2}(x_n+\mathrm{abs}\left(x_n\right))---\left(4\right)$

式(3)中，U为对角矩阵，其第i行、第i列元素为系数矩阵C的第i列元素之和的倒数，即V为对角矩阵，其第i行、第i列元素为系数矩阵C的第i行元素之和的倒数，即式(4)中，abs(x_n)为对x_n的各个元素取绝对值；当|x_n-x_n-1|/x_n≤10^-6时，上述迭代停止。

下面给出一种具体的实施例，结合附图2、附图3和附图4对本发明的重建结果做进一步的描述。

附图2为一边长为20mm的立方形仿体，该仿体内有一个半径为1mm，高度为2mm的圆柱形荧光目标。附图3是利用本发明方法得到的重建结果，显示的是z＝9.5mm截面的荧光目标分布。附图4是利用现有技术中的联合代数迭代方法得到的重建结果，显示的是z＝9.5mm截面的荧光目标分布。该实施例中，荧光目标的真实中心位置为(15,7,9.5)mm，采用本发明方法得到的荧光目标中心位置为(14.8,6.3,9.7)mm，重建误差为0.7mm，迭代500次耗时1.4s。采用已有技术得到的荧光目标中心位置为(13.3,6.0,9.9)mm，重建误差为1.9mm，迭代500次耗时6.3s。基于本发明的方法与已有技术相比，重建精度提高，重建速度加快，是一种有效的激发荧光断层成像重建方法。

上列详细说明是针对本发明可行实施例的具体说明，该实施例并非用以限制本发明的专利范围，凡未脱离本发明所为的等效实施或变更，均应包含于本案的专利范围中。

Claims (2)

1.一种激发荧光断层成像重建方法，其特征在于：包括以下步骤：

S1：基于光传输模型和有限元理论，将成像对象的光学特性参数和解剖结构信息作为先验信息，建立成像对象的表面测量数据与内部荧光目标分布的系统方程；

S2：对S1建立的系统方程进行预处理，使得预处理后的系统方程的系统方程各列相互正交；

S3：融合荧光目标空间分布的稀疏特性，从S2中的系统方程中选取若干列，构建新的系数矩阵和新的系统方程；

S4：加入荧光目标空间分布非负的约束条件，采用联合代数迭代技术对S3中新的系统方程进行求解，从而获得成像对象内部的荧光目标的三维分布。

2.根据权利要求1所述的激发荧光断层成像重建方法，其特征在于：

包括以下步骤：

(1)获取测量数据；

(1a)激发光源对固定在电控旋转台上的成像对象进行多角度的透射式断层扫描；

(1b)使用光学检测仪器获得测量数据；

(2)获得成像对象的解剖结构信息以及光学特性参数；

(3)基于光传输模型和有限元理论，将成像对象的解剖结构信息和光学特性参数作为先验信息，建立表面的测量数据与成像对象内部荧光目标分布的系统方程y＝Ax；式中，y是已知的荧光测量数据向量，A是已知的系统矩阵，x是要求解的荧光目标分布向量，且x和y的元素均非负；

(4)对建立的y＝Ax进行如下预处理：

(4a)对A^TA进行特征值分解，得到QΛQ^T＝A^TA，其中T表示转置运算；

(4b)计算权值矩阵W＝Λ^-1Q^TA^T，将权值矩阵分别乘以系统方程的左侧和右侧，得到z＝Bx，式中，z＝Wy＝Λ^-1Q^TA^Ty，B＝Q^T；

(5)从上述矩阵B中挑选k列元素，记为系数矩阵C，C由以下步骤得到：

(5a)初始化，余量r₀＝z，系数矩阵C₀为空矩阵，迭代次数t＝1；

(5b)计算向量α＝|B^Tr_t-1|，取α的最大元素对应的下标，记为m_t，取矩阵B的第m_t列，记为与C_t-1合并，构造系数矩阵C_t，具体为

(5c)计算并将s_t中的小于0的元素设为0；

(5d)计算新的余量r_t＝z-C_ts_t，迭代次数加1，得到t＝t+1；

(5e)重覆步骤5b-5d直到迭代次数满足t＝k；

(6)以步骤(5)得到的s_k作为荧光目标分布初始值，即x₀＝s_k，利用联合代数迭代方法求解光源分布，迭代更新公式为：

x_n＝x_n-1+UC^TV(z-Cx_n-1)U为对角矩阵，V为对角矩阵，

$x_n=\frac{1}{2}(x_n+abs(x_n\left)\right)$

(7)步骤(6)的迭代停止准则为|x_n-x_n-1|/x_n≤10^-6。