一种基于多光源分辨的生物发光断层成像重建方法
【技术领域】
[0001] 本发明属于分子影像领域，涉及一种基于多光源分辨的生物发光断层成像重建方 法。
【背景技术】
[0002] 生物发光断层成像（以下简称BLT)是光学分子影像技术中的一种重要成像模态， 它利用荧光素酶与荧光素底物在活性分子或细胞内的生化反应进行成像；通过发光光源的 精确定位，定性、定量地反映生物体分子细胞水平的变化。
[0003] BLT重建是一个不适定问题，其解不唯一、并且易受测量噪声的影响。在目前研 究中主要通过设置光源可行区域，自适应有限元、自适应网格、有效的正则化方法等策略来 提高重建结果的准确性。冯金超在BLT研究中提出了一种依据重建光源能量密度的可行 区域自动确定方法，避免了人工标定可行区域的不稳定性；自动可行区域的区域选择很大 程度上依赖于重建的光源，由于一般不能保证每次重建光源的准确性，因此这种自动可行 区域重建方法并不普适；吕玉杰在BLT研究中依据一种参数估计方法对有限元网格进行了 自适应细分，不仅提高了定量、定位精度，而且使重建算法的鲁棒性得到了提高；自适应有 限元和自适应网格虽然有很高的稳定性，然而网格细分必须重新计算系统矩阵，数学计算 复杂，计算代价较大。在BLT重建算法研究方面常用的算法包括：基于1^范数的重建算法 (Tikhonov正则化）；全变差重建算法；基于1^范数的重建算法（IVTCG重建算法，Homtopy 算法等），基于贪婪思想的重建算法；基于LP(0〈P〈1)范数的重建算法。然而，这些数学优化 算法都面临参数的确定，不适当的参数将严重影响算法的重建结果。除了单光源重建之外， 在BLT预临床研究中多个病变区域的标定、中心定位、大小分析对于指导临床治疗有着更 加重要的意义。然而，多光源重建相比单光源重建其重建难度更大，目前对于多目标的分 辨，在BLT研究中并没有一个鲁棒、高效的重建方法。因此，多光源重建一直作为BLT研究 的一个主要方面。
[0004] 参考文献：
[0005] 【I 】Virostko J，Powers A C，Jansen E D. Val idation of luminescent source reconstruction using single-view spectrally resolved bioluminescence images[J] · Applied optics，2007, 46 (13):2540-2547.
[0006] 【2】Frey B J，Dueck D.Clustering by passing messages between data points [J]· science，2007, 315(5814) :972-976.
【发明内容】
[0007] 针对上述现有技术中存在的问题和缺陷，本发明的目的在于，提供一种基于多光 源分辨的生物发光断层成像重建方法。
[0008] 为了实现上述目的，本发明采用如下技术方案：
[0009] -种基于多光源分辨的生物发光断层成像重建方法，具体包括以下步骤：
[0010] 步骤I :获得成像目标的解剖结构信息、光学特异性参数以及表面多光谱荧光数 据；
[0011] 步骤2 :对步骤1得到的表面多光谱荧光数据进行归一化处理，得到归一化后的表 面荧光数据；
[0012] 步骤3 :基于光传输模型和有限元方法构建表面荧光数据与成像目标内部光源分 布的线性关系；
[0013] 步骤4 :将步骤3得到的线性关系转化为最优化问题；
[0014] 步骤5 :对步骤4中的最优化问题，采用自适应可行区域迭代收缩策略进行代数迭 代求解，得到最优解X%完成重建。
[0015] 进一步地，所述的基于多光源分辨的生物发光断层成像重建方法，还包括以下步 骤：
[0016] 步骤6 :对步骤5得到的最优解矿进行AP中心聚类，计算得到每一类的聚类中心 m ；
[0017] 步骤7 :将步骤6中得到的聚类中心，作为k-means聚类算法的初始聚类中心，进 行k-means聚类，计算得到每一类的聚类中心C ;
[0018] 步骤8 :将步骤5得到的最优解Xl行几何3D展示，实现成像目标内部光源的重 建结果的显示；对步骤7得到的聚类中心C进行几何3D展示，实现成像目标内部光源中心 位置的确定。
[0019] 具体地，所述步骤3的具体过程包括：
[0020] 步骤3. 1 :根据步骤1得到的解剖结构信息和光学特异性参数得到光传输模型，所 述的光传输模型为扩散近似方程，利用扩散近似方程来描述光在成像目标内的传输过程；
[0021] 步骤3. 2 :利用nirfast工具对成像目标进行网格离散，获得成像目标的离散网格 数据，包括网格节点和四面体，得到网格节点和四面体之间的对应关系；
[0022] 步骤3. 3 :根据离散网格数据和有限元理论，对步骤3. 1中的扩散近似方程进行离 散，进而构建表面焚光数据与成像目标内部光源分布的线性关系：
[0024] 其中，Φ (λ)为波长为λ的表面荧光数据，Μ( λ)表示波长λ下的荧光数据的系 统矩阵，X是要求解的成像目标内部光源能量分布，是非负的；
[0025] 步骤3. 4 :对步骤3. 3中的系统矩阵Μ( λ )进行以下处理：
[0027] 则步骤3. 3中的线性关系变为
为处理后的系统矩阵。
[0028] 具体地，所述步骤4中将得到的线性关系转化为最优问题的方法为：
[0029] 将步骤3. 4中的Q p )进行以下处理：
[0031 ] 具体地，所述步骤5的具体过程包括：
[0032] 步骤5. I :将初始光源可行区域R设置为整个成像目标，并且设定区域收缩因子
其中N1为整个成像目标内包含的节点，N $迭代终止时非0解的数目， Lniax为区域迭代收缩的次数；
[0033] 步骤5. 2 :结合可行区域R，步骤4的最优化问题变为
其 中，
为可行区域R下的系统矩阵，采用代数迭代算法ART，得到可行区域R下的最优 解 X(R);
[0034] 步骤5. 3 :根据步骤5. 2求得的可行区域R下的最优解X (R)，以及此可行区域下的 系统矩阵AU,计算成像目标表面的光通量Φ?λ)。计算穿透成像目标的多光谱荧光 数据Φ?Α)与光通量Φε(λ)之差，即
其 中，i表示迭代次数（0 < i彡Lmax);
[0035] 步骤5. 4 :将步骤5. 3得到的f⑴与min(f (1)，…，f (i-Ι))进行比较，如果f(i) 小于min (f (1)，…，f (i-1))，则将全局最优解矿(R)用X (R)进行替换，否则，不作替换；
[0036] 步骤5. 5 :将步骤5. 2求得的最优解X (R)从大到小进行排序，选取前length (R) / β个值对应的节点来构建新的可行区域R ;
[0037] 步骤5. 6 :判断迭代次数i是否达到最大值Lniax，如果达到，终止区域迭代收缩，输 出全局最优解Xi(R);如果没有达到，转至步骤5. 2 ;
[0038] 步骤5.7 :由全局最优解Xi(R)得到整个成像目标内的最优解X%若节点位于可行 区域R内，Xi= Xi(R);若节点位于可行区域R外部，Xi= 0。
[0039] 具体地，所述步骤6的过程包括：
[0040] 步骤6. 1 :选取Xi中的非0元素对应的节点构建相似度矩阵D，其中，D为对称矩 阵 D(i,