基于自适应牛顿硬阈值追踪法生物发光断层成像重建方法

1.本发明属于分子影像技术领域，涉及基于自适应牛顿硬阈值追踪算法的生物发光断层成像的重建方法。


背景技术：

2.分子影像学技术相比于传统的医学成像技术，能够利用特异性的分子探针在分子水平上实时监测肿瘤的发生和发展，为疾病的早发现早诊断早治疗奠定基础。分子影像的主要模态有：正电子发射断层成像(pet)、核磁共振成像(mri)、计算机断层成像(ct)、超声成像(usi)、光学成像(oi)等。光学成像利用特异性分子探针与肿瘤细胞靶向结合，采用ccd相机采集生物体表面的光学信息，分辨率和灵敏度都比较高，价格低廉，并且无电离辐射。生物发光成像(bli)在临床前小动物研究中发挥着越来越重要的作用，作为一种非侵入性和非电离成像方法，bli通常用于肿瘤观察和体内治疗评价。但这种平面光学分子成像技术所得到的二维信号是体内光源在三维空间中穿透不同组织后的叠加信号，受生物组织对光子的吸收和散射及目标距离体表深度等因素影响，它不能准确反映目标在体内的空间位置信息，更不能反映其大小和形状信息，作为bli的一个延伸，生物发光断层成像(blt)是一种更有前景的成像技术，它以无创和动态的方式提供了细胞和分子水平上相互作用过程的三维可视化。
3.不幸的是，blt重建仍然面临一些挑战，包括高不适定性和高计算复杂度。在这些挑战中，最关键的问题是不适定性，这是由生物组织的强散射特性和有限的噪声边界测量引起的，导致了blt重建的不准确性。因此，在blt研究中，如何快速、准确地求解blt逆问题是非常必要的。
4.在过去的几十年里，人们设计了各种方法来获得更精确的解决方案。blt重建通过对解引入不同的约束来缓解不适定问题。对于解的分布约束，常用的约束方法有l1、l2和l
p
(0《p《1)正则化。blt重建的典型方法是基于tikhonov正则化的方法。基于tikhonov正则化的方法在原反问题中加入了l2范数约束。基于tikhonov正则化方法的主要优点是求解l2范数约束问题的比较简单，许多标准的优化方法，如牛顿法和共轭梯度法，都能有效地求解基于tikhonov的方法。然而，l2范数的重建结果过于平滑。另一种常用的方法是稀疏重建方法。blt通常用于早期肿瘤成像，在这一阶段，与受试者全身相比，肿瘤很小并且稀疏。因此，肿瘤可以看做是稀疏源。lasso通常用来恢复稀疏信号，套索在原反问题中加入了l1范数约束。为了有效地求解lasso，人们提出了各种方法。稀疏重建方法可以有效地从有限的测量数据中恢复出信号。因此，稀疏重建方法是blt中最常用的方法。近年来，非凸正则化方法在blt重建中引起了人们的关注。非凸正则化是l
p
(0《p《1)范数正则化。与l
p
(1≤p≤2)范数相比，当p
→
0时，l
p
(0《p《1)范数正则化将更好地逼近l0范数。重建结果表明，非凸正则方法比l1范数方法具有更好的重建效果。
5.值得一提的是，这些方法都是基于凸优化或近似凸优化的理论，一般都面临时间效率低的问题。由于收敛速度慢，blt的规模可达数万维，极大地限制了实际应用，并且正则
化参数选择较为困难。为了优化重建效率，最近的许多工作都采用贪婪算法来代替迭代的近似凸优化方法。它是一种针对欠定线性系统的信号恢复技术，利用解的稀疏性约束，可以从比传统shannon-nyquist采样定理所要求的更少的样本中恢复稀疏信号。正则化正交匹配追踪(romp)通过在每次迭代中选择一个以上的合格原子来减少时间开销。cosamp还旨在通过比romp更好的原子选择方法来优化重建效率。romp在每次迭代中选择的原子被保留，而cosamp在每次迭代中选择的原子可以在下一次迭代中被丢弃。然而，这些算法需要事先知道稀疏性，在实际的重建过程中，稀疏性参数是未知的。稀疏自适应匹配追踪算法(samp)，无需稀疏参数即可完成重建。该算法首先设置一个较小的稀疏度，并在每次迭代中增加固定步长，使估计的稀疏度接近真实稀疏度。这种方法容易导致稀疏值估计不准确，从而影响重建的准确性。
6.综上所述，设计一种重建精度高，速度快，能够自适应计算参数的模型尤为重要。


技术实现要素：

7.针对现有技术中的缺陷和不足，本发明提供了基于自适应牛顿硬阈值追踪法生物发光断层成像重建方法，大大减少了迭代次数，提升了重建速度。
8.本发明是这样实现的：
9.一种基于自适应牛顿硬阈值追踪法生物发光断层成像重建方法，包括以下步骤：
10.步骤一、获取生物体表面光通量、重建目标的解剖结构信息和重建目标的光学特异性参数，利用amira软件对重建目标进行离散以获得网格结构信息；
11.步骤二、基于扩散近似模型和robin边界条件，结合自由空间光传输和有限元理论，把重建目标的解剖结构信息和光学特异性参数作为先验信息，建立生物体表面光通量和生物体内发光目标网格结构信息的线性关系；
12.步骤三、使用自适应牛顿硬阈值追踪算法求解步骤二的线性关系，获得重建结果。
13.可选的，所述的步骤三中使用自适应牛顿硬阈值追踪算法求解步骤二的线性关系包括：
14.通过迭代硬阈值算法选择初始集，定义一个非线性函数自适应调整稀疏度的步长，使估计的稀疏度能够快速接近真实的稀疏度，然后利用最小二乘法对初始集进行修剪得到支撑集，在支撑集上，基于强wolfe线搜索准则，采用修正阻尼牛顿法来获得最优重建结果。
15.可选的，所述的步骤三中使用自适应牛顿硬阈值追踪算法求解步骤二的线性关具体包括：
16.将线性关系转化为基于l0范数的p-稀疏最优化问题：
17.使得||x||0≤p；
18.其中，x表示需要重建的生物体内发光目标网格结构信息，b表示生物体表面光通量，a是一个系统矩阵，f(x)为目标函数，||x||0表示l0范数，即x向量中所有非零元素的个数，p为稀疏度。
19.可选的，所述的步骤三中使用自适应牛顿硬阈值追踪算法求解步骤二的线性关具体为：
20.3.1、采用硬阈值算法进行迭代计算，迭代初始值为零向量：
21.x
k 1
＝hq(xk a
t
(b-axk))；
22.其中a
t
为a的转置，k为迭代次数，hq是一个非线性算子，作用是每次迭代保留q个最大值；
23.3.2、目标函数f(x)在x处的梯度为：
[0024][0025]
代入3.1的计算式得：
[0026][0027]
其中，η是不动点，
“‑”
表示梯度下降方向，将每次迭代的结果取绝对值，再选取最大的q个值作为每次迭代的初始集，记为sk：
[0028][0029]
3.3、令初始支撑集u
k-1
为空集，和初始集合并作为备选集ck：
[0030]ck
＝u
k-1
∪sk；
[0031]
3.4、通过最小二乘法计算当前估计的光源分布
[0032][0033]
3.5、从当前估计的光源分布中选取pk个绝对值最大的作为支撑集：
[0034][0035]
3.6、采用线搜索方法估计最终光源分布并计算最终残差rk：
[0036][0037][0038]
其中是在支撑集上的搜索方向，αk是搜索步长；
[0039]
3.7、若最近一次计算的残差小于或等于阈值，则停止迭代；若最近一次计算的残差大于前一次计算的残差，则自适应调整稀疏度的步长；否则，更新支撑集和残差：
[0040][0041]
p
k 1
＝pk s
k 1
；
[0042]
其中s是步长，v是常数，p是稀疏度。
[0043]
可选的，所述3.6具体包括：
[0044]
1)使用修正阻尼牛顿法计算搜索方向：
[0045][0046]
其中为支撑集的补集，为牛顿方向，定义收敛条件如下：
[0047][0048]
其中γ是常数，若在迭代过程中以上的收敛条件不满足：则采用梯度下降方向为搜索方向：
[0049][0050]
其中为梯度下降的方向；
[0051]
2)使用强wolfe准则计算线搜索步长：
[0052][0053]
其中σ是0到1之间的常数。
[0054]
可选的，所述步骤一具体包括：利用ccd图像传感器探测被测生物体表面的光学信息。
[0055]
可选的，所述步骤二具体包括：
[0056]
a)获取重建目标的解剖结构信息：对计算机断层成像投影数据进行三维重建，采用3dmed软件预处理获得重建目标的三维体数据；采用3dmed软件中的人机交互式半自动化分割方法对体数据进行组织分割，获得重建目标的解剖结构；
[0057]
b)获取光学特异性参数：利用解剖结构信息和应用基于生物组织特异性的光学三维重建方法的基于区域的扩散光学层析成像算法获得重建目标内各个组织的光学特异性参数。
[0058]
可选的，还包括：
[0059]
步骤四、将重建结果和重建目标的解剖结构信息进行图像融合，并进行展示。
[0060]
可选的，所述的展示具体包括：使用tecplot软件进行三维立体和二维截面展示。
[0061]
采用上述方法，本发明的优点如下所示：
[0062]
第一，传统的基于l0范数的贪婪算法需要自行估计稀疏度参数。本发明可以通过迭代自适应调整稀疏度参数：定义一个非线性函数自适应调整稀疏度的步长，使估计的稀疏度快速接近真实的稀疏度，重建位置更加精确，避免了估计带来的误差；
[0063]
第二，使用最小二乘法对初始集进行剪枝得到支撑集，每次迭代将支撑集与初始集合并，增加了选到正确原子的概率；
[0064]
第二，引入强wolfe准则，建立全局收敛，解决了传统牛顿法对初始值的限制，保证了算法的稳定性和鲁棒性；
[0065]
第四，引入修正阻尼牛顿法计算搜索方向，建立二次收敛，大大减少了迭代次数，提升了重建速度。
附图说明
[0066]
附图是用来提供对本公开的进一步理解，并且构成说明书的一部分，与下面的具
体实施方式一起用于解释本公开，但并不构成对本公开的限制。在附图中：
[0067]
图1为数字鼠仿真模型；
[0068]
图2为生物组织表面光分布信息示例，单光源实例；
[0069]
图3为生物组织表面光分布信息示例，双光源示例；
[0070]
图4为本发明实施例提供的单光源数值仿真实验的实验结果，压缩采样匹配追踪(cosamp)算法重建结果；
[0071]
图5为本发明实施例提供的单光源数值仿真实验的实验结果，迭代硬阈值(iht)算法重建结果；
[0072]
图6为本发明实施例提供的单光源数值仿真实验的实验结果，弹性网(elastic-net)算法重建结果；
[0073]
图7为本发明实施例提供的单光源数值仿真实验的实验结果，本发明重建结果；
[0074]
图8为本发明实施例提供的双光源数值仿真实验的实验结果，cosamp算法重建结果；
[0075]
图9为本发明实施例提供的双光源数值仿真实验的实验结果，iht算法重建结果；
[0076]
图10为本发明实施例提供的双光源数值仿真实验的实验结果，elastic-net算法重建结果；
[0077]
图11为本发明实施例提供的双光源数值仿真实验的实验结果，本发明重建结果。
[0078]
具体实施实例
[0079]
为了使本发明更加清楚明白，以下结合说明书附图和具体实施例对本发明的内容做具体说明。应当理解这里所描述的具体实施例仅用于解释本发明不用于限定本发明。
[0080]
本发明中提到的重建目标为生物体，具体还可以指生物体的某组织部分或某器官，比如在实施例中提到的躯干部位，主要脏器肝肾脾胃等。
[0081]
本发明的基于自适应牛顿硬阈值追踪法生物发光断层成像重建方法，包括以下步骤：(需要先在目标生物体内放置发光源)
[0082]
步骤一、获取生物体表面光通量、重建目标的解剖结构信息和光学特异性参数，利用amira软件对重建目标进行离散以获得网格结构信息，即获得需要重建的发光目标的坐标信息；
[0083]
步骤二、基于扩散近似模型和robin边界条件，结合自由空间光传输和有限元理论，把重建目标的解剖结构信息和光学特异性参数作为先验信息，建立生物体表面光通量和生物体内发光目标网格结构信息的线性关系；
[0084]
步骤三、使用自适应牛顿硬阈值追踪算法求解步骤二的线性关系，获得重建结果；
[0085]
步骤四、将重建结果和重建目标的解剖结构信息进行图像融合，并进行展示。
[0086]
进一步，步骤一具体包括：利用ccd(charge coupled device,电荷耦合元件)图像传感器探测被测生物体表面的光学信息。
[0087]
进一步，步骤二具体包括：
[0088]
a)获取重建目标的解剖结构信息
[0089]
对计算机断层成像投影数据进行三维重建，采用3dmed软件预处理获得重建目标的三维体数据；采用3dmed软件中的人机交互式半自动化分割方法对体数据进行组织分割，获得重建目标的解剖结构；
[0090]
b)获取光学特异性参数
[0091]
利用解剖结构信息和应用基于生物组织特异性的光学三维重建方法的基于区域的扩散光学层析成像算法获得重建目标内各个组织的光学特异性参数。
[0092]
进一步，步骤三中使用自适应牛顿硬阈值追踪算法求解步骤二的线性关系包括：
[0093]
通过迭代硬阈值算法选择初始集，定义一个非线性函数自适应调整稀疏度的步长，使估计的稀疏度能够快速接近真实的稀疏度，然后利用最小二乘法对初始集进行修剪得到支撑集，在支撑集上，基于强wolfe线搜索准则，采用修正阻尼牛顿法来获得最优重建结果。
[0094]
进一步，步骤三中使用自适应牛顿硬阈值追踪算法求解步骤二的线性关系具体包括：
[0095]
3.1、将线性关系转化为基于l0范数的p-稀疏最优化问题：
[0096]
使得||x||0≤p；
[0097]
其中，x表示需要重建的生物体内发光目标网格结构信息，b表示生物体表面光通量，a是一个系统矩阵，f(x)为目标函数，||x||0表示l0范数，即x向量中所有非零元素的个数，p为稀疏度；
[0098]
3.2、为了求解上一步骤的最优化问题，首先采用硬阈值算法进行迭代计算，迭代初始值为零向量：
[0099]
x
k 1
＝hq(xk a
t
(b-axk))；
[0100]
其中a
t
为a的转置，k为迭代次数，hq是一个非线性算子，作用是每次迭代保留q个最大值；
[0101]
3.3、目标函数f(x)在x处的梯度为：
[0102][0103]
代入2的计算式得：
[0104][0105]
其中，η是不动点，
“‑”
表示梯度下降方向，将每次迭代的结果取绝对值，再选取最大的q个值作为每次迭代的初始集，记为sk：
[0106][0107]
3.4、令初始支撑集u
k-1
为空集，和初始集合并作为备选集ck：
[0108]ck
＝u
k-1
∪sk；
[0109]
3.5、通过最小二乘法计算当前估计的光源分布
[0110][0111]
3.6、从当前估计的光源分布中选取pk个绝对值最大的作为支撑集：
[0112][0113]
3.7、采用线搜索方法估计最终光源分布并计算最终残差rk：
[0114][0115][0116]
其中是在支撑集上的搜索方向，αk是搜索步长；
[0117]
3.8、若最近一次计算的残差小于或等于阈值(自定义值，如1e-06)，则停止迭代；若最近一次计算的残差大于前一次计算的残差，则自适应调整稀疏度的步长；否则，更新支撑集和残差：
[0118][0119]
p
k 1
＝pk s
k 1
；
[0120]
其中s是步长，v是常数，p是稀疏度；
[0121]
进一步，3.7具体包括：
[0122]
1)使用修正阻尼牛顿法计算搜索方向：
[0123][0124]
其中为支撑集的补集，为牛顿方向。定义收敛条件如下：
[0125][0126]
其中γ是常数。若在迭代过程中以上的收敛条件不满足：则采用梯度下降方向为搜索方向：
[0127][0128]
其中为梯度下降的方向。
[0129]
2)使用强wolfe准则计算线搜索步长：
[0130][0131]
其中σ是0到1之间的常数。
[0132]
进一步，步骤四具体包括：使用tecplot软件进行三维立体和二维截面展示。
[0133]
实施例一：
[0134]
如图1所示，为南加州大学提供的通过pet数据、ct数据和冷冻切片构造的三维数字鼠模型，选取其中高32mm的腹部躯干作为重建目标。包含六个主要器官：肝脏、肺、胃、肾脏、心脏、肌肉。
[0135]
如图2所示，为单光源实例，具体实施步骤如下：
[0136]
1)将一半径为1mm的球形光源植入数字鼠的肝部，中心位置为(18，6.5，14.1mm)；
[0137]
2)采用扩散近似方程描述光在重建目标内的传输过程；
[0138]
3)根据有限元理论，并融合重建目标的解剖信息和光学特异性参数，构建表面测量数据与重建目标内部发光目标分布的线性关系。
[0139]
如图3所示，为双光源实例，具体实施步骤如下：
[0140]
1)将两个半径为1mm的球形光源植入数字鼠的肝部，边与边的间距为6mm，中心位置分别为(13.5，5.7，17.5mm)、(21.5，5.7，17.5mm)；
[0141]
2)采用扩散近似方程描述光在重建目标内的传输过程；
[0142]
3)根据有限元理论，并融合重建目标的解剖信息和光学特异性参数，构建表面测量数据与重建目标内部发光目标分布的线性关系。
[0143]
图4-7是基于本发明的单光源重建结果以及和cosamp、iht、elastic-net三种算法的对比，左边是三维视图，球体为真实光源，球体外的异形体为重建光源；右边是二维视图，切片位置为z＝14.1mm，黑色圆圈为真实光源位置。本发明计算获得的目标中心位置为(18.037，6.325，14.603mm),位置误差le＝||l
r-l0||2＝0.1826mm(lr为重建区域的中心，l0是真实区域的中心)，dice系数(sr是重建区域，s0是真实区域)为0.83，归一化均方误差(x
rec
为重建荧光产额值，x
org
为真实荧光产额值)为0.0120，耗时0.19秒。具体数据见表1。
[0144]
表1单光源对比算法的数值结果
[0145][0146][0147]
图8-11是基于本发明的双光源重建结果以及和cosamp、iht、elastic-net三种算法的对比，左边是三维视图，球体为真实光源，球体外的异形体为重建光源；右边是二维视图，切片位置为z＝17.5mm，黑色圆圈为真实光源位置。本发明计算获得的目标中心位置为(13.266，5.401，17.463mm)和(21.149，5.701，17.261mm),位置误差le分别为0.3818mm和0.4247mm，dice系数分别为为0.83，nrmse为0.0178，耗时0.19秒。双光源的的数值结果如表2所示；
[0148]
表2双光源对比算法的数值结果
[0149][0150]
以上两个实例表明，本发明降低了blt重建的病态性，重建位置误差小，重建速度得到了提高。
[0151]
以上所述为本发明的两个具有代表性的实施例，不用于限制本发明，凡在本发明的精神和原则之内所做的任何修改、等同替换和改进等，均应该包含在本发明的保护范围之内。另外，本发明的具体各个技术特征在一定条件下可以通过任何方式进行组合，为了避免重复，本发明不再另行说明，但仍在本发明的保护范围之内。