一种基于多光源分辨的生物发光断层成像重建方法

摘要

本发明公开了一种基于多光源分辨的生物发光断层成像重建方法，利用光学特性参数与解剖结构信息作为先验知识，结合多光谱荧光数据，提高了BLT光源重建的准确性以及重建图像的分辨率；采用自适应可行区域迭代收缩策略，减少了BLT重建问题的病态性，降低了重建结果对于重建算法的依赖性；提出了一种结合AP聚类和K-means聚类的混合聚类方法，避免了单独采用各自聚类算法的局限性与不足。本方法不仅直接使用了网格节点的能量分布信息，而且充分考虑了重建节点所属四面体能量、中心坐标、体积等信息，在最大程度上提升了聚类算法的稳定性，减少了算法计算时间复杂度。

说明书

技术领域

本发明属于分子影像领域，涉及一种基于多光源分辨的生物发光断层成像重建方法。

背景技术

生物发光断层成像(以下简称BLT)是光学分子影像技术中的一种重要成像模态，它利用荧光素酶与荧光素底物在活性分子或细胞内的生化反应进行成像；通过发光光源的精确定位，定性、定量地反映生物体分子细胞水平的变化。

BLT重建是一个不适定问题，其解不唯一、并且易受测量噪声的影响。在目前研究中主要通过设置光源可行区域，自适应有限元、自适应网格、有效的正则化方法等策略来提高重建结果的准确性。冯金超在BLT研究中提出了一种依据重建光源能量密度的可行区域自动确定方法，避免了人工标定可行区域的不稳定性；自动可行区域的区域选择很大程度上依赖于重建的光源，由于一般不能保证每次重建光源的准确性，因此这种自动可行区域重建方法并不普适；吕玉杰在BLT研究中依据一种参数估计方法对有限元网格进行了自适应细分，不仅提高了定量、定位精度，而且使重建算法的鲁棒性得到了提高；自适应有限元和自适应网格虽然有很高的稳定性，然而网格细分必须重新计算系统矩阵，数学计算复杂，计算代价较大。在BLT重建算法研究方面常用的算法包括：基于L₂范数的重建算法(Tikhonov正则化)；全变差重建算法；基于L₁范数的重建算法(IVTCG重建算法，Homtopy算法等)，基于贪婪思想的重建算法；基于L_P(0<P<1)范数的重建算法。然而，这些数学优化算法都面临参数的确定，不适当的参数将严重影响算法的重建结果。除了单光源重建之外，在BLT预临床研究中多个病变区域的标定、中心定位、大小分析对于指导临床治疗有着更加重要的意义。然而，多光源重建相比单光源重建其重建难度更大，目前对于多目标的分辨，在BLT研究中并没有一个鲁棒、高效的重建方法。因此，多光源重建一直作为BLT研究的一个主要方面。

参考文献：

【1】VirostkoJ,PowersAC,JansenED.Validationofluminescentsourcereconstructionusingsingle-viewspectrallyresolvedbioluminescenceimages[J].Appliedoptics,2007,46(13):2540-2547.

【2】FreyBJ,DueckD.Clusteringbypassingmessagesbetweendatapoints[J].science,2007,315(5814):972-976.

发明内容

针对上述现有技术中存在的问题和缺陷，本发明的目的在于，提供一种基于多光源分辨的生物发光断层成像重建方法。

为了实现上述目的，本发明采用如下技术方案：

一种基于多光源分辨的生物发光断层成像重建方法，具体包括以下步骤：

步骤1：获得成像目标的解剖结构信息、光学特异性参数以及表面多光谱荧光数据；

步骤2：对步骤1得到的表面多光谱荧光数据进行归一化处理，得到归一化后的表面荧光数据；

步骤3：基于光传输模型和有限元方法构建表面荧光数据与成像目标内部光源分布的线性关系；

步骤4：将步骤3得到的线性关系转化为最优化问题；

步骤5：对步骤4中的最优化问题，采用自适应可行区域迭代收缩策略进行代数迭代求解，得到最优解X^*，完成重建。

进一步地，所述的基于多光源分辨的生物发光断层成像重建方法，还包括以下步骤：

步骤6：对步骤5得到的最优解X^*进行AP中心聚类，计算得到每一类的聚类中心m；

步骤7：将步骤6中得到的聚类中心，作为k-means聚类算法的初始聚类中心，进行k-means聚类，计算得到每一类的聚类中心C；

步骤8：将步骤5得到的最优解X^*进行几何3D展示，实现成像目标内部光源的重建结果的显示；对步骤7得到的聚类中心C进行几何3D展示，实现成像目标内部光源中心位置的确定。

具体地，所述步骤3的具体过程包括：

步骤3.1：根据步骤1得到的解剖结构信息和光学特异性参数得到光传输模型，所述的光传输模型为扩散近似方程，利用扩散近似方程来描述光在成像目标内的传输过程；

步骤3.2：利用nirfast工具对成像目标进行网格离散，获得成像目标的离散网格数据，包括网格节点和四面体，得到网格节点和四面体之间的对应关系；

步骤3.3：根据离散网格数据和有限元理论，对步骤3.1中的扩散近似方程进行离散，进而构建表面荧光数据与成像目标内部光源分布的线性关系：

Φ(λ)＝M(λ)X

其中，Φ(λ)为波长为λ的表面荧光数据，M(λ)表示波长λ下的荧光数据的系统矩阵，X是要求解的成像目标内部光源能量分布，是非负的；

步骤3.4：对步骤3.3中的系统矩阵M(λ)进行以下处理：

则步骤3.3中的线性关系变为其中为处理后的系统矩阵。

具体地，所述步骤4中将得到的线性关系转化为最优问题的方法为：

将步骤3.4中的进行以下处理：

具体地，所述步骤5的具体过程包括：

步骤5.1：将初始光源可行区域R设置为整个成像目标，并且设定区域收缩因子其中N_I为整个成像目标内包含的节点，N_F为迭代终止时非0解的数目，L_max为区域迭代收缩的次数；

步骤5.2：结合可行区域R，步骤4的最优化问题变为其中，为可行区域R下的系统矩阵，采用代数迭代算法ART，得到可行区域R下的最优解X(R)；

步骤5.3：根据步骤5.2求得的可行区域R下的最优解X(R)，以及此可行区域下的系统矩阵计算成像目标表面的光通量Φ^c(λ)。计算穿透成像目标的多光谱荧光数据与光通量Φ^c(λ)之差，即其中，i表示迭代次数(0＜i≤L_max)；

步骤5.4：将步骤5.3得到的f(i)与min(f(1),…,f(i-1))进行比较，如果f(i)小于min(f(1),…,f(i-1))，则将全局最优解X^*(R)用X(R)进行替换，否则，不作替换；

步骤5.5：将步骤5.2求得的最优解X(R)从大到小进行排序，选取前length(R)/β个值对应的节点来构建新的可行区域R；

步骤5.6：判断迭代次数i是否达到最大值L_max，如果达到，终止区域迭代收缩，输出全局最优解X^*(R)；如果没有达到，转至步骤5.2；

步骤5.7：由全局最优解X^*(R)得到整个成像目标内的最优解X^*，若节点位于可行区域R内，X^*＝X^*(R)；若节点位于可行区域R外部，X^*＝0。

具体地，所述步骤6的过程包括：

步骤6.1：选取X^*中的非0元素对应的节点构建相似度矩阵D，其中，D为对称矩阵D(i,k)＝D(k,i)，即

其中其中，第i个节点的坐标为(x_i,y_i,z_i)，第k个节点的坐标为(x_k,y_k,z_k)，X^*_i为第i个节点的能量，X^*_k为第k个节点的能量；

步骤6.2：记初始参考度矩阵P为0，初始适选矩阵A为0，初始代表矩阵R为0，设定聚类迭代次数；

步骤6.3：计算代表矩阵R^ok1和适选矩阵A^ok1，计算公式如下：

其中，R^ok1(i,k)表示点k适合作为点i的类代表点的合适程度，A^ok1(i,k)表示点i适合作为点k的类代表点的合适程度，k'满足k'≠k，i'满足i'≠k；

步骤6.4：计算更新后的代表矩阵R和适选矩阵A，交替更新过程如下：

其中，R_n+1表示第n+1次迭代得到的代表矩阵R_n+1^ok1更新后的代表矩阵；R_n表示第n次迭代得到的代表矩阵R_n^ok1更新后的代表矩阵；A_n+1表示第n+1次迭代得到的代表矩阵A_n+1^ok1更新后的代表矩阵；A_n表示第n次迭代得到的代表矩阵A_n^ok1更新后的代表矩阵；λ表示阻尼系数；

步骤6.5：判断聚类迭代次数是否达到设定的聚类迭代次数，如果达到，终止迭代，输出代表矩阵R和适选矩阵A，否则，转至步骤6.3；

步骤6.6：根据代表矩阵R和适选矩阵A对各个节点进行分类，并确定每一类的聚类中心，具体方法如下：

对于任意节点i，如果满足arg_kmax{A(i,k)+R(i,k)}，则节点i和节点k属于同一类；对于每一类内的任意节点m，计算A(m,m)+R(m,m)，选取A(m,m)+R(m,m)中的最大值对应的节点m，作为此类的聚类中心。

具体地，所述步骤7的过程包括：

步骤7.1：根据步骤3得到的网格节点和四面体之间的对应关系，计算所有四面体的能量值，选取能量值非0的四面体作为聚类对象，计算此聚类对象的几何中心；

步骤7.2：依次计算每个能量值非0的四面体的几何中心与各初始聚类中心之间的距离{l₁,…l_c}，其中，c为初始聚类中心的数目，根据最小距离l_i将该四面体划分为第i类；

步骤7.3：重新计算每一类的聚类中心C：

其中p(i)、a(i)、c(i)为每一类第i个非0四面体的能量、体积以及几何中心。

与现有技术相比，本发明具有以下技术效果：

1、本发明利用光学特性参数与解剖结构信息作为先验知识，结合多光谱荧光数据，提高了BLT光源重建的准确性以及重建图像的分辨率。

2、本发明采用自适应可行区域迭代收缩策略，减少了BLT重建问题的病态性，降低了重建结果对于重建算法的依赖性。

3、本发明将重建光源中心的确定问题转换为聚类分析问题，采用AP聚类和k-means聚类算法，有效地利用了AP聚类无需事先指定聚类数目的优点，以及k-means聚类算法的相对可伸缩性，同时又避免了AP聚类的不稳定性以及k-means聚类必须已知簇均值的局限性；在混合聚类过程不仅直接使用了网格节点的能量分布信息，而且充分考虑了节点所属四面体的重建能量、中心坐标、体积体等信息，在最大程度上提升了聚类算法的稳定性，减少了计算时间复杂度。

附图说明

图1为本发明的方法流程图；

图2(a)为仿真实验模型及中心相距10mm的光源分布图；

图2(b)为仿真实验模型及中心相距5mm的光源分布图；

图3(a)为中心相距10mm的探测器测量生物表面光强分布图；

图3(b)为中心相距5mm的探测器测量生物表面光强分布图；

图4(a)为本发明在光源中心相距为10mm情况下的重建光源分布图；

图4(b)为本发明在光源中心相距为10mm情况下的光源聚类中心分布图；

图5(a)为本发明在光源中心相距为5mm情况下的重建光源分布图；

图5(b)为本发明在光源中心相距为5mm情况下的光源聚类中心分布图；

下面结合附图和具体实施方式对本发明的方案做进一步详细地解释和说明。

具体实施方式

实施例

遵从上述技术方案，参见图1，本实施例的基于多光源分辨的生物发光断层成像重建方法，具体包括以下步骤：

步骤1：获得成像目标的解剖结构信息、光学特异性参数以及表面多光谱荧光数据

针对皮下病变部位进行研究，成像目标定义为不考虑边界，且具有一定厚度的皮肤组织，参照文献【1】，生成半径为42mm，高度为24mm的圆柱体进行模拟。

对成像目标进行CT扫描，得到其解剖结构信息，如形状、尺寸、位置等信息。

将上述成像目标固定在电控旋转台上，将光学检测仪器放置在成像目标的上方，成像目标内部生成的光源发出荧光，其波长由生物探针材料本身决定，其波长范围为590～650nm，上述荧光在皮肤中的传输的光学特异性参数，参照文献【1】予以设置。

荧光穿透成像目标，光学检测仪器检测穿透成像目标的多光谱荧光数据，上述荧光数据并未采集成像目标侧面及底面的荧光数据，因此称为单角度数据采集。

所述的光学检测仪器采用高灵敏度的CCD相机。

步骤2：对步骤1得到的表面多光谱荧光数据进行归一化处理

由于不同光谱的荧光吸收和散射性质不同，因此透射到成像目标表面的荧光光强有很大的差异，为了消除这种差异对重建结果的影响，对每一波段的表面多光谱荧光数据进行归一化处理：

其中，Φ(λ)为波长为λ的表面荧光数据，为归一化后的波长为λ的表面荧光数据。

步骤3：基于光传输模型和有限元方法构建表面荧光数据与成像目标内部光源分布的线性关系。具体实现方式如下：

步骤3.1：根据步骤1得到的解剖结构信息和光学特异性参数得到光传输模型，所述的光传输模型为扩散近似方程，利用扩散近似方程来描述光在成像目标内的传输过程。

步骤3.2：利用nirfast工具对成像目标进行网格离散，获得成像目标的离散网格数据，包括网格节点和四面体，得到网格节点和四面体之间的对应关系。

步骤3.3：根据离散网格数据和有限元理论，对步骤3.1中的扩散近似方程进行离散，进而构建表面荧光数据与成像目标内部光源分布的线性关系：

Φ(λ)＝M(λ)X

其中，M(λ)是系统矩阵，表示成像目标内所有单元对光的吸收和散射作用，X是要求解的成像目标内部光源能量分布，是非负的。

步骤3.4：由于步骤2对不同光谱下的表面荧光数据做了归一化处理，因此对步骤3.3中的系统矩阵M做如下处理:

则步骤3.3中的线性关系也变为其中为处理后的系统矩阵。

步骤4：将步骤3.4得到的线性关系转化为最优化问题：

步骤5：对步骤4中最优化问题，采用自适应可行区域迭代收缩策略进行代数迭代求解，得到最优解X^*。具体实现方法如下：

步骤5.1：初始光源可行区域R设置为整个成像目标，并且设定区域收缩因子其中N_I为整个成像目标内包含的节点，N_F为迭代终止时非0解的数目，设为1，L_max为区域迭代收缩的次数，设为50。

步骤5.2：结合可行区域R，步骤4的最优化问题变为其中，为可行区域R下的系统矩阵，具体求解采用代数迭代算法(ART)，得到可行区域R下的最优解X(R)，其反应可行区域R下的每个节点的能量大小。

步骤5.3：根据步骤5.2求得的可行区域R下的最优解X(R)，以及此可行区域下的系统矩阵计算成像目标表面的光通量Φ^c(λ)。计算光学检测仪器检测到的穿透成像目标的多光谱荧光数据与光通量之差，即其中，i表示迭代次数(0＜i≤L_max)。

步骤5.4：将步骤5.3得到的f(i)与min(f(1),…,f(i-1))进行比较，如果f(i)小于min(f(1),…,f(i-1))，则将全局最优解X^*(R)用X(R)进行替换，否则，不作替换。。

步骤5.5：将步骤5.2求得的最优解X(R)从大到小进行排序，选取前length(R)/β个值对应的节点来构建新的可行区域R。

步骤5.6：判断迭代次数i是否达到最大值L_max，如果达到，终止区域迭代收缩，输出全局最优解X^*(R)；如果没有达到，转至步骤5.2。

步骤5.7：由全局最优解X^*(R)得到整个成像目标内的最优解X^*，若节点位于可行区域R内，X^*＝X^*(R)；若节点位于可行区域R外部，X^*＝0。

最优解X^*表示整个成像目标内各个节点的能量分布，其展示在空间分布图上，即实现成像目标内部生成的光源的重建。

为了准确获取成像目标内部生成的光源的中心位置，本发明的方法还包括以下步骤：

步骤6：将步骤5得到的最优解X^*进行AP中心聚类，得到每一类的聚类中心。由于X^*每个元素对应的是网格节点的能量，在中心聚类时，选取X^*中的非0元素对应的节点，按照其空间分布进行聚类，由于AP聚类可以自动确定聚类数目和聚类中心，因此采用AP聚类算法在确定光源数目的基础上，初步对各个光源中心定位。

具体实现方法如下：

步骤6.1：选取X^*中的非0元素对应的节点构建相似度矩阵D，其中，D为对称矩阵D(i,k)＝D(k,i)，即

其中其中，第i个节点的坐标为(x_i,y_i,z_i)，第k个节点的坐标为(x_k,y_k,z_k)，X^*_i为第i个节点的能量，X^*_k为第k个节点的能量。

步骤6.2：记初始参考度矩阵P为0，初始适选矩阵A为0，初始代表矩阵R为0，聚类迭代次数为1000。

步骤6.3：AP聚类算法传递两种类型的消息矩阵：代表矩阵R^ok1和适选矩阵A^ok1，参考文献【2】，计算公式如下：

其中，R^ok1(i,k)表示点k适合作为点i的类代表点的合适程度。A^ok1(i,k)表示点i适合作为点k的类代表点的合适程度。k'满足k'≠k，i'满足i'≠k。

步骤6.4：AP聚类算法的迭代过程就是这两个消息矩阵R和A交替更新的过程，设阻尼系数λ为0.8，计算更新后的代表矩阵R和适选矩阵A，交替更新过程如下：

其中，R_n+1表示第n+1次迭代得到的代表矩阵R_n+1^ok1更新后的代表矩阵；R_n表示第n次迭代得到的代表矩阵R_n^ok1更新后的代表矩阵；A_n+1表示第n+1次迭代得到的代表矩阵A_n+1^ok1更新后的代表矩阵；A_n表示第n次迭代得到的代表矩阵A_n^ok1更新后的代表矩阵；

步骤6.5：判断聚类迭代次数是否达到设定的聚类迭代次数，如果达到，终止迭代，输出代表矩阵R和适选矩阵A，否则，转至步骤6.3。

步骤6.6：根据代表矩阵R和适选矩阵A确定聚类中心m

对于任意节点i，如果满足arg_kmax{A(i,k)+R(i,k)}，则节点i和节点k属于同一类；对于每一类内的任意节点m，计算A(m,m)+R(m,m)，选取A(m,m)+R(m,m)中的最大值对应的节点m，作为此类的聚类中心。

步骤7：对步骤6中AP聚类得到的聚类中心，作为k-means聚类算法的初始聚类中心，进行k-means聚类，得到每一类的聚类中心C。由于最优解的稀疏性导致AP聚类结果的稳定性较差，因此k-means聚类依据能量非0四面体的几何中心进行聚类。如具体实现方法如下：

步骤7.1：根据步骤3得到的网格节点和四面体之间的对应关系，计算所有四面体的能量值，选取能量值非0的四面体作为聚类对象，计算此聚类对象的几何中心；上述设计可充分利用网格节点和四面体之间的相关性，获取更多的非零数据进行聚类，最大程度上降低了由于最优解X^*的稀疏性导致的偶然误差，采用四面体空间分布进行聚类，增加了聚类数据源的可靠性与稳定性，从而提高了聚类分析的准确度进而降低了试验误差。

步骤7.2：依次计算每个能量值非0的四面体的几何中心与各初始聚类中心之间的距离{l₁,…l_c}，其中，c为初始聚类中心的数目，根据最小距离l_i将该四面体划分为第i类。

步骤7.3：重新计算每一类的聚类中心C：

其中p(i)、a(i)、c(i)为每一类第i个非0四面体的能量、体积以及几何中心。

步骤8：利用matlab相关工具包，以及nirfast、iso2mesh相关函数将步骤5得到的最优解X^*进行几何3D展示，实现成像目标内部光源的重建结果的显示；对步骤7得到的聚类中心C进行几何3D展示，实现成像目标内部光源中心位置的确定。

实验结果分析

结合附图2、附图3、附图4对本发明的重建结果和聚类中心做进一步的描述。

用于实验的仿体如附图2所示。其为半径为42mm，高为24mm的圆柱，主要模拟深层皮下组织的病变。在试验中设置了中心相距10mm和5mm的两组仿体实验，其中图2(a)模拟5个中心相距为10mm的光源，图2(b)模拟5个中心相距为5mm的光源。

重建过程所需的表面光强分布图如附图3所示。由于模拟皮肤组织，所以图3的表面光强分布均在仿体上表面采集。图3(a)是5个中心相距为10mm的光源表面测量数据。图3(b)是5个中心相距为5mm的光源表面测量数据。对比可知中心相距距离的减少对于重建难度的影响。

最终的多光源重建结果如附图4和附图5所示。其中，附图4为光源中心相距10mm情况下，基于本发明的实验重建结果和聚类结果。附图5为光源中心相距5mm情况下，基于本发明的实验重建结果和聚类结果。其中图4(a)和图5(a)为重建结果展示图，红色区域为重建的光源分布，蓝色的圆圈代表真实的光源中心，蓝色的十字代表聚类结果中心。图4(b)和图5(b)展示了最终的聚类结果图，中心点为最终聚类中心。在临床应用中，由图4(a)、图5(a)和图4(b)、图5(b)可以清楚的得到病变区域，并准确确定区域中心。基于本发明的重建方法，重建光源位置误差小，算法稳定，多目标中心定位准确，是一种有效的多光源生物发光断层成像重建方法。