CN103188989B - 一种动态荧光分子断层图像的重建方法 - Google Patents

Abstract

本发明涉及一种动态荧光分子断层图像的重建方法，包括以下步骤：1）设置一小动物诱发荧光分子成像***；2）采用所述***，在不同时间点全角度、等间隔采集反映荧光探针在成像对象体内代谢分布的荧光投影图像，并将采集图像存储为输入矩阵u；3）根据u生成KL变换矩阵A；4）求解输入矩阵u经KL变换后的数据5）获取经KL变换后需要进行3D重建的L′个KL分量；6）根据光场分布函数Φ_x和格林函数G生成权重矩阵W₁；7）根据公式采用3D FMT重建方法对每一个KL分量分别进行3D重建；8）对KL分量的重建结果进行逆KL变换，得到4D荧光断层序列。本发明可以广泛应用于动态荧光分子图像的重建过程中。

Description

一种动态荧光分子断层图像的重建方法

技术领域

本发明涉及一种荧光分子断层图像的重建方法，特别是关于一种动态荧光分子断层图像的重建方法。

背景技术

目前动态（4D）荧光分子断层图像的重建大都是利用逐帧（frame-by-frame）重建方法完成，即对采集到的每一帧（圈）的荧光投影图像进行逐个重建。此方法的优点是可以直接利用已有的各种诱发荧光分子断层成像（Fluorescence MolecularTomography,FMT)重建方法完成图像重建。然而，在动态FMT成像研究中，例如在药物代谢的研究中，药物（荧光探针）的浓度是随时间变化的，因此所采集到的每一帧测量数据（荧光投影图像）在时间上是高度相关的，如果采用上述逐帧重建方法对荧光投影图像进行4D FMT重建，则忽略了重建过程中荧光探针浓度沿时间轴的相关性，但是如果直接将采集的所有荧光投影图像序列作为一个整体直接进行4D重建，其计算量通常已超出普通计算机的计算能力，例如对一个35帧的荧光投影图像序列（假设在每帧FMT成像过程中，等间隔采集24张不同角度的荧光投影图像）直接进行4D FMT重建，其用于重建的空间-时间权重矩阵将占用大约3234GB的内存，这大大超出了普通计算机的计算能力。

发明内容

针对上述问题，本发明的目的是提供一种充分考虑荧光投影图像序列的时间相关性，且能够有效提高动态荧光分子断层图像重建速度的动态荧光分子断层图像的重建方法。

为实现上述目的，本发明采取以下技术方案：一种动态荧光分子断层图像的重建方法，包括以下步骤：1）设置一包括有计算机、小动物旋转平台装置和荧光成像激发与检测装置的小动物诱发荧光分子成像***；2）采用所述小动物诱发荧光分子成像***，在不同时间点，全角度、等间隔采集反映荧光探针在成像对象体内的代谢分布的荧光投影图像，并将每一时间点所采集的不同角度的荧光投影图像存储为一矩阵序列，将所述矩阵序列标记为输入矩阵u；3）根据输入矩阵u生成KL变换矩阵A；4）基于求解输入矩阵u经KL变换后的数据式中，I_M为M×M的单位矩阵，为克罗内克积；5）获取经KL变换后需要进行3D重建的L′个KL分量；6）根据光场分布函数Φ_x和格林函数G生成权重矩阵W₁，其中，Φ_x描述了对应激发光波长的光子密度，G描述了在荧光谱段，光子在成像对象体内的传播；7）根据公式采用3D FMT重建方法在KL域中对步骤5）中获取的每一个KL分量分别进行3D重建，其中，k＝1,2,...,L′，W^-1 ₁为权重矩阵W₁的逆矩阵；8）对上述步骤7）中获取的KL分量的重建结果进行逆KL变换，得到完整的4D荧光断层图像序列。

所述步骤3）根据输入矩阵u生成KL变换矩阵A的方法为：①构建输入矩阵u的时间协方差矩阵P^time；②通过求解P^timeA^T=A^TD，生成KL变换矩阵A，式中d_l是协方差矩阵P^time的第l个特征值，L是动态成像中扫描的总帧数。

所述步骤3）根据输入矩阵u生成KL变换矩阵A的方法为：通过对输入矩阵u进行奇异值分解求得。

所述步骤5）中获取经KL变换后需要进行3D重建的L′个KL分量采用scree方法、Kaiser方法和基于特征值的累计方差方法中的一种。

所述步骤6）中在检测点r处获取的荧光投影图像采用下式表示：

Φ_m(r)＝∫_r′∈VG(r,r′)n(r′)Φ_x(r′)dr′

其中，W₁(r)＝G(r,r′)Φ_x(r′)ΔV(r)，ΔV是重建区域V被离散成的小网格；V是重建区域；Φ_m是发射荧光波长的光子密度；n(r′)为重建区域V中点r′处的待重建的荧光产额。

所述步骤7）采用3D FMT重建方法在KL域中对步骤5）中获取的每一个KL分量分别进行3D重建中，3D FMT重建方法采用ART方法、LSQR方法和Bayesian方法中的一种。

本发明由于采取以上技术方案，其具有以下优点：1、本发明在不同时间点，全角度、连续采集荧光探针在成像对象体内的代谢分布，且将采集的所有荧光投影图像作为一个整体进行重建，因此在重建过程中充分考虑到了荧光投影图像在时间上的相关性，弥补了现有重建方法对时间相关性考虑的不足。2、本发明通过对测量得到的荧光投影数据进行KL变换，去掉了荧光投影数据间的时间相关性，将复杂的4D FMT重建问题简化为一系列KL分量在KL域中的3D重建问题，因此与直接对测量数据进行4D重建的方法相比降低了4D FMT重建的计算规模。3、本发明采用的KL变换具有很高的数据压缩性，将荧光投影数据进行KL变换后，有用的信号被保存在前几个KL分量中，采用现有的scree方法或其它类似方法提取需要进行FMT重建的KL分量的个数，并对提取的每一KL分量分别采用现有的3D FMT重建方法进行重建，然后将重建后的数据进行逆KL变换则可以恢复出完整的4D断层图像序列，因此有效提高4D FMT重建方法的速度，经本发明的仿真试验表明，采用本发明的重建方法，4D FMT的重建时间大约为传统重建方法的1.5％。本发明可以广泛应用于动态荧光分子图像的重建过程中。

附图说明

图1是本发明基于KL变换的4D FMT重建方法结构框图；

图2是本发明的仿真实验结构示意图，其中，图2（a）是主视示意图，图2（b）是图2（a）的俯视示意图；

图3是吲哚菁绿ICG（Indocyanine Green，ICG）在小鼠肝部的代谢曲线示意图，也是本发明仿真实验中所使用的ICG浓度时间曲线示意图，横坐标为时间，纵坐标为ICG在小鼠内的浓度值；

图4（a）～（d）是采用本发明方法计算所得仿真实验的FMT重建结果的效果示意图，图4（a）是第2帧（2min）荧光投影图像的重建结果效果示意图，图4（b）是第20帧（20min）荧光投影图像重建结果效果示意图；图4（c）是第50帧（50min）荧光投影图像重建结果效果示意图；图4（d）是第120帧（120min）荧光投影图像重建结果效果示意图；

图5是采用本发明的重建方法所获得的ICG在仿真实验中的浓度时间变化曲线示意图，其中，圆形为图3中的实际的ICG浓度，圆点为采用本发明方法计算所得到的ICG浓度，横坐标为采集时间点。

具体实施方式

下面结合附图和实施例对本发明进行详细的描述。

对于FMT成像，基于一阶波恩近似，在检测点r处获取的荧光投影图像可以采用下式表示：

Φ_m(r)＝∫_r′∈VG(r,r′)n(r′)Φ_x(r′)dr′    （1）

式中，V是重建区域；Φ_x为光场分布函数，描述了对应激发光波长的光子密度，Φ_m为发射（荧光）波长的光子密度；n(r′)为重建区域V中点r′处的待重建的荧光产额，通常n(r′)正比于荧光探针的浓度；格林函数G(r,r′)描述在荧光波段，光子在成像对象体内的传播。

对于动态的FMT，重建区域V被离散成小的网格ΔV后，公式（1）可被离散为下述线性方程：

u_k=W₁n_k    （2）

式中，u_k是一个M×1的向量，用以描述第k帧的荧光投影图像；n_k是一个N×1的向量，用以描述待重建的荧光断层图像；W₁(r)＝G(r,r′)Φ_x(r′)ΔV(r)是基于漫射方程生成的M×N权重矩阵，用以描述光子在成像对象体内的传播。

对于4D荧光断层图像重建，现有的方法是采用公式（2）逐帧进行重建，得到每帧的断层图像n_k，然后将所有的断层图像进行组合生成动态的FMT图像序列，上述重建方法忽略了荧光投影图像u_k在时间轴上的相关性。为了解决上述问题，可以将所有测量得到的荧光投影图像序列作为整体进行4D FMT重建，重建公式为：

u＝Wn    （3）

式中，是所采集到的完整荧光投影图像序列，其中描述了第i帧的荧光投影图像；L是动态成像中扫描的总帧数（总次数）；一系列待重建的荧光断层图像，其中描述了第i时刻荧光探针在成像体内的浓度分布；W＝diag[W₁,W₁,…,W₁]是动态FMT成像的空间-时间***矩阵，其大小为ML×NL；但是此种直接对荧光投影图像序列进行4D重建的方法，W的内存占用及随后的计算量都超出了普通计算机能力。

为了解决上述问题，本发明提出了基于KL变换的4D FMT重建方法，即通过对公式（3）进行KL（Karhunen–Loève）变换。在动态FMT的研究中，本发明令λ=[λ¹,λ²,...,λ^L]^T表示像素的时间强度曲线（Time Activity Curve,TAC），KL变换和KL逆变换的公式分别为：

KL变换： $\widetilde{\mathrm{\λ}}=\mathrm{A\λ}---\left(4\right)$

KL逆变换： $\mathrm{\λ}=A^T\widetilde{\mathrm{\λ}}---\left(5\right)$

式中，是经KL变换后的数据；A是KL变换矩阵，A^T是KL变换的逆变换矩阵（在KL变换中KL的逆变换矩阵近似等于KL变换矩阵的转置矩阵）；KL变换矩阵A可以通过以下公式进行求解：

P^timeA^T=A^TD（6）

式中，P^time为所有TAC的协方差矩阵，矩阵大小为L×L：

$\left[P^{\mathrm{time}}\right]=\frac{1}{N-1}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}({\mathrm{\λ}}_i^l-{\stackrel{\‾}{\mathrm{\λ}}}^l)({\mathrm{\λ}}_i^n-{\stackrel{\‾}{\mathrm{\λ}}}^n)---\left(7\right)$

式中， ${\stackrel{\‾}{\mathrm{\λ}}}^l=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\λ}}_i^l,$ ${\stackrel{\‾}{\mathrm{\λ}}}^n=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\λ}}_i^n,$ $D=\mathrm{diag}{\left\{d_1\right\}}_{l=1}^L,$ d_l是协方差矩阵P^time的第l个特征值。在动态FMT重建中，P^time的维数通常较小（L＜200），因此在公式（6）中，特征矢量能够被快速计算。

为了描述荧光投影图像u在时间方向上的KL变换，本发明定义如下KL变换矩阵A_M：

$A_M=A\⊗I_M---\left(8\right)$

式中，I_M是一个M×M的单位矩阵；是克罗内克积（kronecker product）。

将公式（2）的两侧分别乘以A_M，基于Wernick（韦尼克）推导，得到以下公式：

$A_{M}u=A_M\mathrm{Wn}$

$=(A\⊗I_M)(I_L\⊗W_1)n$

$=\left(A{I}_L\right)\⊗\left(I_M{W}_1\right)n---\left(9\right)$

$\left(I_{L}A\right)\⊗\left(W_1{I}_N\right)n$

$=(I_K\⊗W_1)(A\⊗I_N)n$

$={\mathrm{WA}}_{N}n$

令：

$\tilde{u}=A_{M}u---\left(10\right)$

$\tilde{n}=A_{N}n---\left(11\right)$

将公式（10）和公式（11）代入公式（9），公式（9）可以写为：

$\tilde{u}=W\tilde{n}---\left(12\right)$

由于KL变换能够去除荧光投影图像u的时间相关性。因此可以对KL变换后的每个KL分量基于公式（12）在KL域中单独进行重建，即利用KL变换的去相关性，将复杂的4D荧光投影数据的重建转化为一系列KL分量在KL域的3D重建问题，降低了4D FMT重建的计算规模。

此外，由于公式（12）与公式（3）具有相同的成像模型，因此可以采用现有各种的FMT重建方法完成各个KL分量的重建，而不需对其进行修改，对重建后的KL分量再进行逆KL变换，就能够获取完整的4D断层图像序列n：

$n=A_N^T\tilde{n}---\left(13\right)$

式中，I_N是一个N×N的单位矩阵，由于KL变化具有压缩属性，经过KL变换后，信号被保存在前几个KL分量，更高阶的KL分量几乎没有包含有用的信号，通过抛弃高阶KL分量，即仅通过重建前几个KL分量，然后对重建后的KL分量进行逆KL变就可以恢复出完整的4D断层序列，在选择需保留的KL分量个数时可以采用现有的“scree”方法、“Kaiser”方法或基于特征值的累计方差方法。

在实验过程中获得荧光投影图像可以采用现有的非接触、全角度小动物诱发荧光分子成像***，小动物诱发荧光分子成像***包括有计算机、小动物旋转平台装置和荧光成像激发与检测装置。

如图1所示，下面通过具体实施例进一步说明对动态荧光投影图像序列进行4D荧光分子断层图像重建的方法流程，具体包括以下步骤：

1）在不同时间点t₁、t₂…t_L′…t_L，使用全角度、非接触FMT成像***，等时间间隔采集反映荧光探针在成像对象体内的代谢分布的荧光投影图像（Proj.1、Proj.2…Proj.K）。对于每个成像时刻（帧）所采集的不同角度的荧光投影图像，将其转化为一个列向量，通过变换所有帧，得到KL变换的输入矩阵u（图1中纵向虚线框内的方框分别表示各采集时刻所采集到不同角度的荧光投影图像序列1、2…L′…L）。

如图2（a）～（b）所示，本发明以模拟ICG在小鼠肝中的代谢情况（如图3所示）的仿真实验为实施例进一步说明动态荧光投影图像序列的采集和重建过程。

本发明参照小鼠形体大小和组织的平均光学参数，仿真实验设置一直径为3.0cm、高为4.5cm的圆柱作为仿体1，仿体1的光学约化散射系数为10cm^-1，光学吸收系数为0.3cm^-1，将一个直径为0.43cm、高为0.5cm的小圆柱作为荧光棒2放置在仿体1内后，将仿体1放置在小动物旋转平台装置上，开启激光光源照射荧光棒2，利用荧光成像激发与检测装置对成像对象进行全角度的荧光投影图像采集，成像过程中，为了模拟ICG在小鼠肝部的时间代谢曲线，在成像的不同时刻，根据小鼠肝部代谢曲线将不同时刻的ICG浓度赋值给荧光棒（如图3所示），对于每一时间点（帧）成像，等间隔采集24个角度的荧光投影图像，即对每一时间点共生成24个不同角度的荧光投影图像。考虑到由于在在体实验中，一帧的成像时间近似为1min，因此在本发明的仿真实验中，设置每帧图像的采集间隔为1min。具体的动态仿真过程为，首先，根据给定的ICG的浓度时间曲线（如图3所示），在不同的时间点（从1min到120min）为荧光棒2设置不同的ICG浓度。然后，根据公式（2）生成120帧（每一min为一帧）的荧光投影图像序列，用以描述模拟ICG在小鼠肝内不同时刻的浓度分布信息。最后，将所有采集到的荧光投影图像序列（总计120帧，每帧24个不同角度的荧光投影数据）组合成一个整体，形成动态的荧光投影图像序列，作为KL变换的输入矩阵u。

2）根据公式（7），构建输入矩阵u的时间协方差矩阵P^time。

3）通过求解P^timeA^T=A^TD，生成L×L的KL变换矩阵A，在动态FMT成像中，由于总的扫描次数L通常较小，因此变换矩阵A的计算是快速的。

4）基于公式（10），获取经KL变换后的数据（图1中细线框内的 表示经KL变换后得到的KL分量）。

5）采用“scree”方法，保留KL分量中的前L′个KL分量（L′＝2表示取前两个KL分量）。

6）根据光场分布函数Φ_x和格林函数G生成权重矩阵W₁。

7）基于公式k＝1,2,...,L′，使用3D FMT重建方法，对步骤5）中获取的每一个KL分量分别进行3D重建得到重建结果（如图1中粗实线框内得到的L′个KL分量的重建结果）。其中，3D FMT重建方法为现有的方法，可以采用代数重建方法（Algebraic Reconstruction Technique,ART）方法、最小二乘QR分解迭代法（LeastSquares QR，LSQR）方法或Bayesian方法。

8）如图4（a）～（d）所示，基于公式（13），将步骤7）的重建结果进行逆KL变换，得到完整的4D荧光断层图像序列n₁、n₂…n_L（如图1所示）。

上述实施例中，步骤2）和3）根据输入矩阵u生成KL变换矩阵A的方法为还可以采用通过对输入矩阵u进行奇异值分解求得。

上述实施例中，本发明的实施例采用仿真实验，在仿真实验中，改变荧光探针的浓度只需改变浓度设定值即可；在仿体实验中，改变荧光探针的浓度需要将放置荧光探针的玻璃棒取出，清洗，注入新的不同浓度的荧光探针，然后在放回仿体内；在活体实验中，在尾静脉注射荧光探针后，荧光探针浓度随时间自行改变。

如图5所示，将各个时刻重建的ICG的浓度取平均值，生成采用本方法计算得到的不同时刻所对应的ICG浓度的分布图，图中圆点表示120个采集时间点所对应的使用本方法计算得到ICG的浓度，图中的圆形为图3的给定的实际ICG的浓度时间曲线，可以看出两者是基本接近的，充分说明了本发明方法的有效性。

上述各实施例仅用于说明本发明，其中本发明的方法的实施步骤等都是可以有所变化的，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。

Claims (9)

1.一种动态荧光分子断层图像的重建方法，包括以下步骤：

1)设置一包括有计算机、小动物旋转平台装置和荧光成像激发与检测装置的小动物诱发荧光分子成像***；

2)采用所述小动物诱发荧光分子成像***，在不同时间点，全角度、等时间间隔采集反映荧光探针在成像对象体内的代谢分布的荧光投影图像，并将每一时间点所采集的不同角度的荧光投影图像存储为一矩阵序列，将所述矩阵序列标记为输入矩阵u；

3)根据输入矩阵u生成KL变换矩阵A；

4)基于求解输入矩阵u经KL变换后的数据式中，I_M为M×M的单位矩阵，为克罗内克积；

5)获取经KL变换后需要进行3D重建的L'个KL分量；

6)根据光场分布函数Φ_x和格林函数G生成权重矩阵W₁，其中，Φ_x描述了对应激发光波长的光子密度，G描述了在荧光谱段，光子在成像对象体内的传播；

7)根据公式采用3D FMT重建方法在KL域中对步骤5)中获取的每一个KL分量分别进行3D重建，其中，k＝1,2,...,L'，W^-1 ₁为权重矩阵W₁的逆矩阵；

8)对上述步骤7)中获取的KL分量的重建结果进行逆KL变换，得到完整的4D荧光断层图像序列。

2.如权利要求1所述的一种动态荧光分子断层图像的重建方法，其特征在于：所述步骤3)根据输入矩阵u生成KL变换矩阵A的方法为：

①构建输入矩阵u的时间协方差矩阵P^time；

②通过求解P^timeA^T＝A^TD，生成KL变换矩阵A，式中d_l是协方差矩阵P^time的第l个特征值，L是动态成像中扫描的总帧数。

3.如权利要求1所述的一种动态荧光分子断层图像的重建方法，其特征在于：所述步骤3)根据输入矩阵u生成KL变换矩阵A的方法为：通过对输入矩阵u进行奇异值分解求得。

4.如权利要求1或2或3所述的一种动态荧光分子断层图像的重建方法，其特征在于：所述步骤5)中获取经KL变换后需要进行3D重建的L'个KL分量采用基于特征值的累计方差方法。

5.如权利要求1或2或3所述的一种动态荧光分子断层图像的重建方法，其特征在于：所述步骤6)中在检测点r处获取的荧光投影图像采用下式表示：

${\mathrm{\Φ}}_m\left(r\right)={\∫}_{r^{\′}\∈V}^{\mathrm{}}G(r,r^{\′})n\left(r^{\′}\right){\mathrm{\Φ}}_x\left(r^{\′}\right){\mathrm{dr}}^{\′}$

其中，W₁(r)＝G(r,r')Φ_x(r')ΔV(r)，ΔV是重建区域V被离散成的小网格；V是重建区域；Φ_m是发射荧光波长的光子密度；n(r')为重建区域V中点r′处的待重建的荧光产额。

6.如权利要求4所述的一种动态荧光分子断层图像的重建方法，其特征在于：所述步骤6)中在检测点r处获取的荧光投影图像采用下式表示：

${\mathrm{\Φ}}_m\left(r\right)={\∫}_{r^{\′}\∈V}^{\mathrm{}}G(r,r^{\′})n\left(r^{\′}\right){\mathrm{\Φ}}_x\left(r^{\′}\right){\mathrm{dr}}^{\′}$

其中，W₁(r)＝G(r,r')Φ_x(r')ΔV(r)，ΔV是重建区域V被离散成的小网格；V是重建区域；Φ_m是发射荧光波长的光子密度；n(r')为重建区域V中点r′处的待重建的荧光产额。

7.如权利要求1或2或3或6所述的一种动态荧光分子断层图像的重建方法，其特征在于：所述步骤7)采用3D FMT重建方法在KL域中对步骤5)中获取的每一个KL分量分别进行3D重建中，3D FMT重建方法采用ART方法、LSQR方法和Bayesian方法中的一种。

8.如权利要求4所述的一种动态荧光分子断层图像的重建方法，其特征在于：所述步骤7)采用3D FMT重建方法在KL域中对步骤5)中获取的每一个KL分量分别进行3D重建中，3D FMT重建方法采用ART方法、LSQR方法和Bayesian方法中的一种。

9.如权利要求5所述的一种动态荧光分子断层图像的重建方法，其特征在于：所述步骤7)采用3D FMT重建方法在KL域中对步骤5)中获取的每一个KL分量分别进行3D重建中，3D FMT重建方法采用ART方法、LSQR方法和Bayesian方法中的一种。