CN109886888B - 一种基于l1范数最小化模型的神经纤维骨架校正方法 - Google Patents

Abstract

本发明公开了一种基于L1范数最小化模型的神经纤维骨架校正方法，包括：获取存在扭曲神经纤维结构的原始图像；对原始图像进行重建，获得待校正的神经纤维初始骨架；固定初始骨架上的两个端点，基于L1范数最小化模型，校正两个端点之间所有节点位置。本发明通过L1范数最小化模型校正神经纤维骨架，构建的L1范数最小化模型包括两部分，一部分用于测算节点的信号值，保证纤维骨架点尽可能向其在真实图像中局部信号最强处聚集，另一部分运用骨架点之间的二阶差分反映纤维骨架的光滑性，在尽可能维持纤维平滑的特征中，保留真实纤维的扭曲结构稀疏分布的特点，相比现有以L2范数模型为基础的校正方法，更适用于真实图像中存在扭曲结构的纤维骨架校正。

Description

一种基于L1范数最小化模型的神经纤维骨架校正方法

技术领域

本发明属于图像处理领域，更具体地，涉及一种基于L1范数最小化模型的神经纤维骨架校正方法。

背景技术

神经元是由一系列管状结构的纤维，通过一定次序连接而成的树形结构，其形态重建是为了将这类树形结构精确量化，从而为神经科学研究提供定量精准的研究材料。现有的先进成像和标记方法已能提供近乎完整的精细神经纤维图像。

然而，神经图像中复杂的纤维结构阻碍着神经元纤维骨架的准确获取。具体地，神经纤维中存在许多纤维方向突变的扭曲结构，现有自动重建方法(例如基于L2范数模型的骨架校正方法)在刻画这类结构上尚存在困难，这将直接导致后续依赖于重建结果的形态统计信息发生错误，影响相关的神经科学研究。

发明内容

针对现有技术的缺陷，本发明的目的在于解决现有技术重建结果偏离真实神经图像的技术问题。

为实现上述目的，第一方面，本发明实施例提供了一种基于L1范数最小化模型的神经纤维骨架校正方法，所述方法包括以下步骤：

步骤S1.获取存在扭曲神经纤维结构的原始图像；

步骤S2.对原始图像进行重建，获得待校正的神经纤维初始骨架；

步骤S3.固定初始骨架上的两个端点，基于L1范数最小化模型，校正两个端点之间所有节点位置。

具体地，固定初始骨架上的两个端点p₁、p_n，校正两个端点之间所有节点p₂、…、p_n-1位置，所述基于L1范数最小化模型如下：

其中，pi为骨架上第i个节点位置，

为校正后第i个节点位置，i＝2,…，n-1，n为骨架上节点总数，

为L₁范数，λ为常数，取值为0.5；Λ_i是三维图像内点p_i的一个体素范围的邻域，点p表示Λ_i中的一个体素，s(p)是点p处的图像信号强度，|| ||₂表示2范数；σ为常数，取值范围为

具体地，步骤S3包括以下步骤：

S301.将所述基于L1范数最小化模型转化为增广拉格朗日方程的形式，具体如下：

其中，d_i＝2p_i-p_i-1-p_i+1；

其中，r_i为依据线性约束d_i＝2p_i-p_i-1-p_i+1引入的对偶变量r中的第i个分量，μ为常数，取值0.5，<r_i，(2p_i-p_i-1-p_i+1)-d_i>表示r_i与(2p_i-p_i-1-p_i+1)-d_i的内积；

S302.固定d^k和r^k的值，基于增广拉格朗日方程，更新纤维骨架上除去两个端点外所有骨架点的位置，位置集合记为p^k+1，k为迭代次数；

S303.固定p^k+1和r^k的值，使用软阈值法求解d^k+1；

S304.依据已更新的p^k+1和d^k+1，求解r^k+1；

S305.对步骤S301～S304获得的纤维校正结果，进行N次迭代，得到校正后神经纤维最终骨架，N≥10，1≤k≤N。

具体地，步骤S302中，第(k+1)次迭代出的骨架上除去两个端点外所有骨架点p^k+1的计算公式如下：

其中，

为第k次迭代时计算得到的第i个骨架点的位置，r_i为矩阵r中的第i个列向量，d_i为矩阵d中的第i个列向量，n为待校正纤维骨架点，r^k与d^k相等，矩阵d中的第i个元素

k为迭代次数。

具体地，步骤S302具体为：利用牛顿法求解第一项

用梯度法分段求解后两项，第一项和后两项之和满足

i＝2,…,n-1，S为后两项之和。

具体地，所述使用软阈值法求解d^k+1，具体如下：

其中，d^k+1为

的集合，

为第(k+1)次迭代时计算得到的第i个骨架点的位置，T为软阈值操作算子，并满足T_λ(ω)＝[t_λ(ω₁)，t_λ(ω₂),...]^T，t_λ(ω_i)＝sgn(ω_i)max{0,|ω_i-λ|}。

具体地，r^k+1的计算公式如下：

其中，r^k+1为

的集合，δ为步长。

具体地，δ取值0.5。

具体地，进行前10次迭代计算时，计算p^k+1和d^k+1所涉及到的常数项

和常数矩阵μ^-1λ在每次迭代后，均按照1.5倍线性增大；10次迭代后常数项和常数矩阵维持不变。

第二方面，本发明实施例提供了一种计算机可读存储介质，该计算机可读存储介质上存储有计算机程序，该计算机程序被处理器执行时实现上述第一方面所述的神经纤维骨架校正方法。

总体而言，通过本发明所构思的以上技术方案与现有技术相比，具有以下有益效果：

本发明通过L1范数最小化模型校正神经纤维骨架，构建的L1范数最小化模型包括两部分，一部分用于测算节点的信号值，保证纤维骨架点尽可能向其在真实图像中局部信号最强处聚集，另一部分运用骨架点之间的二阶差分反映纤维骨架的光滑性，在尽可能维持纤维平滑的特征中，保留真实纤维的扭曲结构稀疏分布的特点，相比现有以L2范数模型为基础的校正方法，更适用于真实图像中存在扭曲结构的纤维骨架校正。

附图说明

图1为本发明实施例提供的一种基于L1范数最小化模型的神经纤维骨架校正方法流程图；

图2为本发明实施例提供的神经纤维原始图像；

图3为本发明实施例提供的神经纤维重建骨架的初始结果示意图；

图4为本发明实施例提供的神经纤维重建骨架的校正结果示意图；

图5(a)为本发明实施例提供的存在扭曲结构的纤维及其对应的初始骨架结构示意图；

图5(b)为本发明实施例提供的本发明方法校正结果示意图；

图5(c)为本发明实施例提供的基于L2范数模型的骨架校正方法校正结果示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

如图1所示，一种基于L1范数最小化模型的神经纤维骨架校正方法，所述方法包括以下步骤：

步骤S1.获取存在扭曲神经纤维结构的原始图像；

步骤S2.对原始图像进行重建，获得待校正的神经纤维初始骨架；

步骤S3.固定初始骨架上的两个端点，基于L1范数最小化模型，校正两个端点之间所有节点位置。

步骤S1.获取存在扭曲神经纤维结构的原始图像。

原始图像为管状或柱状图像。一般而言，大部分神经纤维的方向是缓慢变化的，具有平滑性，小部分结构会存在纤维扭曲，方向突变的情况。上述纤维结构特征可抽象为一组信号序列，大部分平缓变化的纤维方向可视为值趋近于零，少部分突变的方向可视为非零值。这类系数在零附近集中的现象即为纤维结构特征的“稀疏性”。如图2所示，原始图像选取荧光显微切片成像系统或功能性双光子共聚焦成像显微镜获取的小鼠大脑切片图像。

步骤S2.对原始图像进行重建，获得待校正的神经纤维初始骨架。

使用自动重建算法对原始图像进行重建，获得待校正的神经纤维重建骨架。骨架由一系列节点表示p₁,p₂,…,p_n-1,p_n，本发明实施例n＝80。如图3所示，圆框内初始重建骨架明显偏离真实图像。

步骤S3.固定初始骨架上的两个端点，基于L1范数最小化模型，校正两个端点之间所有节点位置。

固定初始骨架上的两个端点(p₁、p_n)，基于L1范数最小化模型的纤维骨架校正算法，校正两个端点之间所有节点(p₂、…、p_n-1)位置。所述校正模型如下：

其中，pi为骨架上第i个节点位置，

为校正后第i个节点位置，n为骨架上节点总数，

为L₁范数，λ为常数，取值为0.5；Λ_i是三维图像内点p_i的一个体素范围的邻域，点p表示Λ_i中的一个体素，s(p)是点p处的图像信号强度，|| ||₂表示2范数；σ为常数，取值范围为

参数λ用于控制纤维的光滑性，第一项函数g(p_i)用于测算节点的信号值，保证纤维骨架点尽可能向其在真实图像中局部信号最强处聚集，本发明实施例中采用均值漂移(mean shift)算法实现。通过最小化g(p_i)可令当前骨架点p_i移动至图像局部信号峰值处。参数σ具体取值是根据真实纤维的半径决定。

第二项函数

运用骨架点之间的二阶差分反映纤维骨架的光滑性。其作用是，在尽可能维持纤维平滑的特征中，保留真实纤维的扭曲结构稀疏分布的特点。本发明实施例中采用L1范数最小化模型实现。

步骤S3具体包括以下步骤：

S301.将校正模型的优化问题转化为增广拉格朗日方程的形式。

其中，d_i＝2p_i-p_i-1-p_i+1。

其中，r_i为依据线性约束d_i＝2p_i-p_i-1-p_i+1引入的对偶变量r中的第i个分量，μ为常数，取值0.5，<r_i，(2p_i-p_i-1-p_i+1)-d_i>表示r_i与(2p_i-p_i-1-p_i+1)-d_i的内积。

S302.固定d^k和r^k的值，基于增广拉格朗日方程，更新纤维骨架上除去两个端点外所有骨架点的位置，位置集合记为p^k+1，k为迭代次数。

第(k+1)次迭代出的骨架上除去两个端点外所有骨架点p^k+1的计算公式如下：

其中，

为第k次迭代时计算得到的第i个骨架点的位置，r_i为矩阵r中的第i个列向量，d_i为矩阵d中的第i个列向量，n为待校正纤维骨架点，r^k与d^k相等，矩阵d中的第i个元素

k为迭代次数。

具体地，首先利用牛顿法求解第一项

具体如下：

其中，(p-p_i)²可表示为(p-p_i)²＝(p-p_i)^TSigmaM(p-p_i)，SigmaM为引入的方向矩阵。

假如，当前骨架点p_i相关的

和

两方向近乎一致，这表明当前方向Dircx对骨架点p_i的校正影响力很小，因此方向矩阵对应修正为SigmaM＝inv(0.1*DircxDircx^T+DircyDircy^T+DirczDircz^T)，其中，Dircx为

和

构成的方向，Dircy和Dircz为与Dircx垂直的两个方向。反之，SigmaM为单位矩阵。取当前点p_i的邻域范围为9×9×7体素，上述邻域选取范围在保证考虑邻域信号对当前信号影响力的同时，又不会过度增大计算量。

再用梯度法分段求解后两项，最后第一项和后两项之和满足式子

可得到与骨架点位置相关的p^k+1值的解：

S303.固定p^k+1和r^k的值，使用软阈值法求解d^k+1。

其中，d^k+1为

的集合，

为第(k+1)次迭代时计算得到的第i个骨架点的位置，T为软阈值操作算子，并满足T_λ(ω)＝[t_λ(ω₁),t_λ(ω₂)，...]^T，t_λ(ω_i)＝sgn(ω_i)max{0，|ω_i-λ|}。

S304.依据已更新的p^k+1和d^k+1，求解r^k+1。

其中，r^k+1为

的集合，δ为步长。

由于拉格朗日增广方程是关于r的线性函数，因此对偶问题可由梯度下降法求解，每次下降的步长δ为预先设定的值，本发明实施例中选取0.5。

S305.对于步骤S301～S304获得的纤维校正结果，进行N次迭代，可得到较为贴合真实图像的纤维重建骨架，N≥10。

对神经纤维重建骨架进行多次迭代校正(迭代次数为10～20次)，直至校正结果能贴合真实图像中的纤维。进行前10次迭代计算时，S302和S303中计算p^k+1和d^k+1所涉及到的常数项

和常数矩阵μ^-1λ(初值为单位矩阵，维度与d^k一致)在每次迭代后均按照1.5倍线性增大。10次迭代后常数项和常数矩阵维持不变。经过上述校正后，校正骨架结果见图4，图4中的圆框位置与图3一致，说明本发明能有效校正重建骨架，使之贴合真实纤维内存在的扭曲结构。

如图5(a)所示，存在扭曲结构的纤维及其对应的初始骨架结构，扭曲部分及部分初始重建骨架被放大至右上角处显示。图5(b)和图5(c)分别为利用本发明方法和基于L2范数模型的骨架校正方法校正初始重建骨架后的结果。比较图5(b)和图5(c)中的放大图(图中右上角)可知，本发明方法和基于L2范数模型的骨架校正方法在校正存在扭曲结构的纤维时，存在差异：本发明方法能较好的校正纤维扭曲结构处的真实骨架，而基于L2范数模型的骨架校正方法则会失效。

本发明能有效校正存在扭曲结构的纤维骨架，使之贴合真实图像。本方法可以应用于多类管状纤维结构的校正，并特别适用于具有方向突变的神经轴突纤维的校正。

以上，仅为本申请较佳的具体实施方式，但本申请的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本申请揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本申请的保护范围之内。因此，本申请的保护范围应该以权利要求的保护范围为准。

Claims (10)

1.一种基于L1范数最小化模型的神经纤维骨架校正方法，其特征在于，所述方法包括以下步骤：

步骤S1.获取存在扭曲神经纤维结构的原始图像；

步骤S2.对原始图像进行重建，获得待校正的神经纤维初始骨架；

步骤S3.固定初始骨架上的两个端点，基于L1范数最小化模型，校正两个端点之间所有节点位置。

2.如权利要求1所述的神经纤维骨架校正方法，其特征在于，固定初始骨架上的两个端点p₁、p_n，校正两个端点之间所有节点p₂、…、p_n-1位置，所述基于L1范数最小化模型如下：

其中，p_i为骨架上第i个节点位置，

为校正后第i个节点位置，i＝2,…,n-1，n为骨架上节点总数，g(p_i)用于测算第i个节点的信号值，

为L₁范数，λ为常数，取值为0.5；Λ_i是三维图像内点p_i的一个体素范围的邻域，点p表示Λ_i中的一个体素，s(p)是点p处的图像信号强度，|| ||₂表示2范数；σ为常数，取值范围为

3.如权利要求2所述的神经纤维骨架校正方法，其特征在于，步骤S3包括以下步骤：

S301.将所述基于L1范数最小化模型转化为增广拉格朗日方程的形式，具体如下：

其中，d_i＝2p_i-p_i-1-p_i+1；

其中，r_i为依据线性约束d_i＝2p_i-p_i-1-p_i+1引入的对偶变量r中的第i个分量，μ为常数，取值0.5，<r_i,(2p_i-p_i-1-p_i+1)-d_i>表示r_i与(2p_i-p_i-1-p_i+1)-d_i的内积；

S302.固定d^k和r^k的值，基于增广拉格朗日方程，更新纤维骨架上除去两个端点外所有骨架点的位置，位置集合记为p^k+1，k为迭代次数；

S303.固定p^k+1和r^k的值，使用软阈值法求解d^k+1；

S304.依据已更新的p^k+1和d^k+1，求解r^k+1；

S305.对步骤S301～S304获得的纤维校正结果，进行N次迭代，得到校正后神经纤维最终骨架，N≥10，1≤k≤N。

4.如权利要求3所述的神经纤维骨架校正方法，其特征在于，步骤S302中，第(k+1)次迭代出的骨架上除去两个端点外所有骨架点p^k+1的计算公式如下：

其中，

为第k次迭代时计算得到的第i个骨架点的位置，r_i为矩阵r中的第i个列向量，d_i为矩阵d中的第i个列向量，n为待校正纤维骨架点，r^k与d^k相等，矩阵d中的第i个元素

k为迭代次数。

5.如权利要求4所述的神经纤维骨架校正方法，其特征在于，步骤S302具体为：利用牛顿法求解第一项

用梯度法分段求解后两项，第一项和后两项之和满足

i＝2,…,n-1，S为后两项之和。

6.如权利要求3所述的神经纤维骨架校正方法，其特征在于，所述使用软阈值法求解d^{k +1}，具体如下：

其中，d^k+1为

的集合，

为第(k+1)次迭代时计算得到的第i个骨架点的位置，T为软阈值操作算子，并满足T_λ(ω)＝[t_λ(ω₁),t_λ(ω₂),...]^T，t_λ(ω_i)＝sgn(ω_i)max{0,|ω_i-λ|}。

7.如权利要求3所述的神经纤维骨架校正方法，其特征在于，r^k+1的计算公式如下：

其中，r^k+1为r_i ^k+1的集合，δ为步长。

8.如权利要求7所述的神经纤维骨架校正方法，其特征在于，δ取值0.5。

9.如权利要求3所述的神经纤维骨架校正方法，其特征在于，进行前10次迭代计算时，计算p^k+1和d^k+1所涉及到的常数项

和常数矩阵μ^-1λ在每次迭代后，均按照1.5倍线性增大；10次迭代后常数项和常数矩阵维持不变。

10.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质上存储有计算机程序，所述计算机程序被处理器执行时实现如权利要求1至9任一项所述的神经纤维骨架校正方法。

Images