CN109886888A - 一种基于l1范数最小化模型的神经纤维骨架校正方法 - Google Patents
一种基于l1范数最小化模型的神经纤维骨架校正方法 Download PDFInfo
- Publication number
- CN109886888A CN109886888A CN201910081471.2A CN201910081471A CN109886888A CN 109886888 A CN109886888 A CN 109886888A CN 201910081471 A CN201910081471 A CN 201910081471A CN 109886888 A CN109886888 A CN 109886888A
- Authority
- CN
- China
- Prior art keywords
- skeleton
- nerve fibre
- point
- fiber
- endpoints
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Landscapes
- Magnetic Resonance Imaging Apparatus (AREA)
Abstract
本发明公开了一种基于L1范数最小化模型的神经纤维骨架校正方法,包括:获取存在扭曲神经纤维结构的原始图像;对原始图像进行重建,获得待校正的神经纤维初始骨架;固定初始骨架上的两个端点,基于L1范数最小化模型,校正两个端点之间所有节点位置。本发明通过L1范数最小化模型校正神经纤维骨架,构建的L1范数最小化模型包括两部分,一部分用于测算节点的信号值,保证纤维骨架点尽可能向其在真实图像中局部信号最强处聚集,另一部分运用骨架点之间的二阶差分反映纤维骨架的光滑性,在尽可能维持纤维平滑的特征中,保留真实纤维的扭曲结构稀疏分布的特点,相比现有以L2范数模型为基础的校正方法,更适用于真实图像中存在扭曲结构的纤维骨架校正。
Description
技术领域
本发明属于图像处理领域,更具体地,涉及一种基于L1范数最小化模型的神经纤维骨架校正方法。
背景技术
神经元是由一系列管状结构的纤维,通过一定次序连接而成的树形结构,其形态重建是为了将这类树形结构精确量化,从而为神经科学研究提供定量精准的研究材料。现有的先进成像和标记方法已能提供近乎完整的精细神经纤维图像。
然而,神经图像中复杂的纤维结构阻碍着神经元纤维骨架的准确获取。具体地,神经纤维中存在许多纤维方向突变的扭曲结构,现有自动重建方法(例如基于L2范数模型的骨架校正方法)在刻画这类结构上尚存在困难,这将直接导致后续依赖于重建结果的形态统计信息发生错误,影响相关的神经科学研究。
发明内容
针对现有技术的缺陷,本发明的目的在于解决现有技术重建结果偏离真实神经图像的技术问题。
为实现上述目的,第一方面,本发明实施例提供了一种基于L1范数最小化模型的神经纤维骨架校正方法,所述方法包括以下步骤:
步骤S1.获取存在扭曲神经纤维结构的原始图像;
步骤S2.对原始图像进行重建,获得待校正的神经纤维初始骨架;
步骤S3.固定初始骨架上的两个端点,基于L1范数最小化模型,校正两个端点之间所有节点位置。
具体地,固定初始骨架上的两个端点p1、pn,校正两个端点之间所有节点p2、…、pn-1位置,所述基于L1范数最小化模型如下:
其中,pi为骨架上第i个节点位置,为校正后第i个节点位置,i=2,…,n-1,n为骨架上节点总数,为L1范数,λ为常数,取值为0.5;Λi是三维图像内点pi的一个体素范围的邻域,点p表示Λi中的一个体素,s(p)是点p处的图像信号强度,|| ||2表示2范数;σ为常数,取值范围为
具体地,步骤S3包括以下步骤:
S301.将所述基于L1范数最小化模型转化为增广拉格朗日方程的形式,具体如下:
其中,di=2pi-pi-1-pi+1;
其中,ri为依据线性约束di=2pi-pi-1-pi+1引入的对偶变量r中的第i个分量,μ为常数,取值0.5,<ri,(2pi-pi-1-pi+1)-di>表示ri与(2pi-pi-1-pi+1)-di的内积;
S302.固定dk和rk的值,基于增广拉格朗日方程,更新纤维骨架上除去两个端点外所有骨架点的位置,位置集合记为pk+1,k为迭代次数;
S303.固定pk+1和rk的值,使用软阈值法求解dk+1;
S304.依据已更新的pk+1和dk+1,求解rk+1;
S305.对步骤S301~S304获得的纤维校正结果,进行N次迭代,得到校正后神经纤维最终骨架,N≥10,1≤k≤N。
具体地,步骤S302中,第(k+1)次迭代出的骨架上除去两个端点外所有骨架点pk+1的计算公式如下:
其中,为第k次迭代时计算得到的第i个骨架点的位置,ri为矩阵r中的第i个列向量,di为矩阵d中的第i个列向量,n为待校正纤维骨架点,rk与dk相等,矩阵d中的第i个元素k为迭代次数。
具体地,步骤S302具体为:利用牛顿法求解第一项用梯度法分段求解后两项,第一项和后两项之和满足i=2,…,n-1,S为后两项之和。
具体地,所述使用软阈值法求解dk+1,具体如下:
其中,dk+1为的集合,为第(k+1)次迭代时计算得到的第i个骨架点的位置,T为软阈值操作算子,并满足Tλ(ω)=[tλ(ω1),tλ(ω2),...]T,tλ(ωi)=sgn(ωi)max{0,|ωi-λ|}。
具体地,rk+1的计算公式如下:
其中,rk+1为的集合,δ为步长。
具体地,δ取值0.5。
具体地,进行前10次迭代计算时,计算pk+1和dk+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.对原始图像进行重建,获得待校正的神经纤维初始骨架。
使用自动重建算法对原始图像进行重建,获得待校正的神经纤维重建骨架。骨架由一系列节点表示p1,p2,…,pn-1,pn,本发明实施例n=80。如图3所示,圆框内初始重建骨架明显偏离真实图像。
步骤S3.固定初始骨架上的两个端点,基于L1范数最小化模型,校正两个端点之间所有节点位置。
固定初始骨架上的两个端点(p1、pn),基于L1范数最小化模型的纤维骨架校正算法,校正两个端点之间所有节点(p2、…、pn-1)位置。所述校正模型如下:
其中,pi为骨架上第i个节点位置,为校正后第i个节点位置,n为骨架上节点总数,为L1范数,λ为常数,取值为0.5;Λi是三维图像内点pi的一个体素范围的邻域,点p表示Λi中的一个体素,s(p)是点p处的图像信号强度,|| ||2表示2范数;σ为常数,取值范围为
参数λ用于控制纤维的光滑性,第一项函数g(pi)用于测算节点的信号值,保证纤维骨架点尽可能向其在真实图像中局部信号最强处聚集,本发明实施例中采用均值漂移(mean shift)算法实现。通过最小化g(pi)可令当前骨架点pi移动至图像局部信号峰值处。参数σ具体取值是根据真实纤维的半径决定。
第二项函数运用骨架点之间的二阶差分反映纤维骨架的光滑性。其作用是,在尽可能维持纤维平滑的特征中,保留真实纤维的扭曲结构稀疏分布的特点。本发明实施例中采用L1范数最小化模型实现。
步骤S3具体包括以下步骤:
S301.将校正模型的优化问题转化为增广拉格朗日方程的形式。
其中,di=2pi-pi-1-pi+1。
其中,ri为依据线性约束di=2pi-pi-1-pi+1引入的对偶变量r中的第i个分量,μ为常数,取值0.5,<ri,(2pi-pi-1-pi+1)-di>表示ri与(2pi-pi-1-pi+1)-di的内积。
S302.固定dk和rk的值,基于增广拉格朗日方程,更新纤维骨架上除去两个端点外所有骨架点的位置,位置集合记为pk+1,k为迭代次数。
第(k+1)次迭代出的骨架上除去两个端点外所有骨架点pk+1的计算公式如下:
其中,为第k次迭代时计算得到的第i个骨架点的位置,ri为矩阵r中的第i个列向量,di为矩阵d中的第i个列向量,n为待校正纤维骨架点,rk与dk相等,矩阵d中的第i个元素k为迭代次数。
具体地,首先利用牛顿法求解第一项具体如下:
其中,(p-pi)2可表示为(p-pi)2=(p-pi)TSigmaM(p-pi),SigmaM为引入的方向矩阵。
假如,当前骨架点pi相关的和两方向近乎一致,这表明当前方向Dircx对骨架点pi的校正影响力很小,因此方向矩阵对应修正为SigmaM=inv(0.1*DircxDircxT+DircyDircyT+DirczDirczT),其中,Dircx为和构成的方向,Dircy和Dircz为与Dircx垂直的两个方向。反之,SigmaM为单位矩阵。取当前点pi的邻域范围为9×9×7体素,上述邻域选取范围在保证考虑邻域信号对当前信号影响力的同时,又不会过度增大计算量。
再用梯度法分段求解后两项,最后第一项和后两项之和满足式子可得到与骨架点位置相关的pk+1值的解:
S303.固定pk+1和rk的值,使用软阈值法求解dk+1。
其中,dk+1为的集合,为第(k+1)次迭代时计算得到的第i个骨架点的位置,T为软阈值操作算子,并满足Tλ(ω)=[tλ(ω1),tλ(ω2),...]T,tλ(ωi)=sgn(ωi)max{0,|ωi-λ|}。
S304.依据已更新的pk+1和dk+1,求解rk+1。
其中,rk+1为的集合,δ为步长。
由于拉格朗日增广方程是关于r的线性函数,因此对偶问题可由梯度下降法求解,每次下降的步长δ为预先设定的值,本发明实施例中选取0.5。
S305.对于步骤S301~S304获得的纤维校正结果,进行N次迭代,可得到较为贴合真实图像的纤维重建骨架,N≥10。
对神经纤维重建骨架进行多次迭代校正(迭代次数为10~20次),直至校正结果能贴合真实图像中的纤维。进行前10次迭代计算时,S302和S303中计算pk+1和dk+1所涉及到的常数项和常数矩阵μ-1λ(初值为单位矩阵,维度与dk一致)在每次迭代后均按照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所述的神经纤维骨架校正方法,其特征在于,固定初始骨架上的两个端点p1、pn,校正两个端点之间所有节点p2、…、pn-1位置,所述基于L1范数最小化模型如下:
其中,pi为骨架上第i个节点位置,为校正后第i个节点位置,i=2,…,n-1,n为骨架上节点总数,为L1范数,λ为常数,取值为0.5;Λi是三维图像内点pi的一个体素范围的邻域,点p表示Λi中的一个体素,s(p)是点p处的图像信号强度,||||2表示2范数;σ为常数,取值范围为
3.如权利要求1所述的神经纤维骨架校正方法,其特征在于,步骤S3包括以下步骤:
S301.将所述基于L1范数最小化模型转化为增广拉格朗日方程的形式,具体如下:
其中,di=2pi-pi-1-pi+1;
其中,ri为依据线性约束di=2pi-pi-1-pi+1引入的对偶变量r中的第i个分量,μ为常数,取值0.5,<ri,(2pi-pi-1-pi+1)-di>表示ri与(2pi-pi-1-pi+1)-di的内积;
S302.固定dk和rk的值,基于增广拉格朗日方程,更新纤维骨架上除去两个端点外所有骨架点的位置,位置集合记为pk+1,k为迭代次数;
S303.固定pk+1和rk的值,使用软阈值法求解dk+1;
S304.依据已更新的pk+1和dk+1,求解rk+1;
S305.对步骤S301~S304获得的纤维校正结果,进行N次迭代,得到校正后神经纤维最终骨架,N≥10,1≤k≤N。
4.如权利要求3所述的神经纤维骨架校正方法,其特征在于,步骤S302中,第(k+1)次迭代出的骨架上除去两个端点外所有骨架点pk+1的计算公式如下:
其中,为第k次迭代时计算得到的第i个骨架点的位置,ri为矩阵r中的第i个列向量,di为矩阵d中的第i个列向量,n为待校正纤维骨架点,rk与dk相等,矩阵d中的第i个元素k为迭代次数。
5.如权利要求4所述的神经纤维骨架校正方法,其特征在于,步骤S302具体为:利用牛顿法求解第一项用梯度法分段求解后两项,第一项和后两项之和满足i=2,…,n-1,S为后两项之和。
6.如权利要求3所述的神经纤维骨架校正方法,其特征在于,所述使用软阈值法求解dk +1,具体如下:
其中,dk+1为的集合,为第(k+1)次迭代时计算得到的第i个骨架点的位置,T为软阈值操作算子,并满足Tλ(ω)=[tλ(ω1),tλ(ω2),...]T,tλ(ωi)=sgn(ωi)max{0,|ωi-λ|}。
7.如权利要求3所述的神经纤维骨架校正方法,其特征在于,rk+1的计算公式如下:
其中,rk+1为的集合,δ为步长。
8.如权利要求7所述的神经纤维骨架校正方法,其特征在于,δ取值0.5。
9.如权利要求3所述的神经纤维骨架校正方法,其特征在于,进行前10次迭代计算时,计算pk+1和dk+1所涉及到的常数项和常数矩阵μ-1λ在每次迭代后,均按照1.5倍线性增大;10次迭代后常数项和常数矩阵维持不变。
10.一种计算机可读存储介质,其特征在于,所述计算机可读存储介质上存储有计算机程序,所述计算机程序被处理器执行时实现如权利要求1至9任一项所述的神经纤维骨架校正方法。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910081471.2A CN109886888B (zh) | 2019-01-28 | 2019-01-28 | 一种基于l1范数最小化模型的神经纤维骨架校正方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910081471.2A CN109886888B (zh) | 2019-01-28 | 2019-01-28 | 一种基于l1范数最小化模型的神经纤维骨架校正方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN109886888A true CN109886888A (zh) | 2019-06-14 |
CN109886888B CN109886888B (zh) | 2021-05-18 |
Family
ID=66927122
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201910081471.2A Active CN109886888B (zh) | 2019-01-28 | 2019-01-28 | 一种基于l1范数最小化模型的神经纤维骨架校正方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN109886888B (zh) |
Cited By (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113591616A (zh) * | 2021-07-14 | 2021-11-02 | 华中科技大学 | 一种基于前景点聚类的神经纤维重建方法和系统 |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
EP2766851A4 (en) * | 2011-10-12 | 2016-05-04 | Seno Medical Instr Inc | SYSTEM AND METHOD FOR COLLECTING OPTOA TACTICAL DATA AND FOR PRODUCING PARAMETER CARDS THEREFOR |
CN106296610A (zh) * | 2016-08-05 | 2017-01-04 | 天津大学 | 基于低秩矩阵分析的三维骨架修复方法 |
CN106683178A (zh) * | 2016-12-30 | 2017-05-17 | 天津大学 | 基于图论的低秩矩阵恢复三维骨架方法 |
-
2019
- 2019-01-28 CN CN201910081471.2A patent/CN109886888B/zh active Active
Patent Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
EP2766851A4 (en) * | 2011-10-12 | 2016-05-04 | Seno Medical Instr Inc | SYSTEM AND METHOD FOR COLLECTING OPTOA TACTICAL DATA AND FOR PRODUCING PARAMETER CARDS THEREFOR |
CN106296610A (zh) * | 2016-08-05 | 2017-01-04 | 天津大学 | 基于低秩矩阵分析的三维骨架修复方法 |
CN106683178A (zh) * | 2016-12-30 | 2017-05-17 | 天津大学 | 基于图论的低秩矩阵恢复三维骨架方法 |
Non-Patent Citations (2)
Title |
---|
HANG ZHOU 等: "Fast Quantifying Discrepancies Between Brain-wide Neuron", 《14TH INTERNATIONAL CONFERENCE ON PHOTONICS AND IMAGING IN》 * |
YOUSEF AL-KOFAHI 等: "Improved Automatic Detection and Segmentation", 《IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING》 * |
Cited By (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113591616A (zh) * | 2021-07-14 | 2021-11-02 | 华中科技大学 | 一种基于前景点聚类的神经纤维重建方法和系统 |
CN113591616B (zh) * | 2021-07-14 | 2024-02-13 | 华中科技大学 | 一种基于前景点聚类的神经纤维重建方法和系统 |
Also Published As
Publication number | Publication date |
---|---|
CN109886888B (zh) | 2021-05-18 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN109727195B (zh) | 一种图像超分辨率重构方法 | |
Escobar et al. | Simultaneous untangling and smoothing of tetrahedral meshes | |
WO2021022929A1 (zh) | 一种单帧图像超分辨率重建方法 | |
CN106204447A (zh) | 基于总变差分和卷积神经网络的超分辨率重建方法 | |
CN106600538A (zh) | 一种基于区域深度卷积神经网络的人脸超分辨率算法 | |
CN110852393A (zh) | 一种遥感图像的分割方法及系统 | |
CN108550108A (zh) | 一种基于相位迭代最小化的傅里叶叠层成像图像重建方法 | |
CN106169095B (zh) | 主动学习大数据标注方法和系统 | |
CN106127689A (zh) | 图像视频超分辨率方法和装置 | |
CN106526839B (zh) | 一种基于模式的同步无波前自适应光学系统 | |
CN109785279B (zh) | 一种基于深度学习的图像融合重建方法 | |
CN113096020B (zh) | 基于平均模式生成对抗网络的书法字体创作方法 | |
CN109598676A (zh) | 一种基于哈达玛变换的单幅图像超分辨率方法 | |
CN109242771A (zh) | 一种超分辨率图像重建方法及装置、计算机可读存储介质和计算机设备 | |
CN112184547B (zh) | 红外图像的超分辨率方法及计算机可读存储介质 | |
CN109447897B (zh) | 一种真实场景图像合成方法及系统 | |
DE112016007098T5 (de) | Indexierung von voxeln für das 3d-drucken | |
KR102543690B1 (ko) | 사전 정보 학습 기반 영상 업스케일링 장치 및 방법 | |
CN113744136A (zh) | 基于通道约束多特征融合的图像超分辨率重建方法和系统 | |
CN113592711A (zh) | 点云数据不均匀的三维重建方法、系统、设备及存储介质 | |
CN109886888A (zh) | 一种基于l1范数最小化模型的神经纤维骨架校正方法 | |
CN107424119A (zh) | 一种单图像的超分辨率方法 | |
CN111105354A (zh) | 基于多源深度残差网络的深度图像超分辨率方法及装置 | |
US20230260083A1 (en) | Computer-implemented method, computer program product and system for processing images | |
CN113674156A (zh) | 一种图像超分辨率重建的方法及系统 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |