一种神经元形状的数字化重建方法及其系统

1.本发明涉及生物医学图形处理领域，尤其涉及一种神经元形状的数字化重建方法及其系统。


背景技术：

2.大脑是人的身体中最重要的器官，机体各项生命活动的正常运转都离不开大脑的控制和调节；同时，大脑也是人体最复杂的器官，大脑中包含有大量的神经元、神经胶质细胞，这些细胞又与大脑中丰富的血管形成相互连接的错综复杂的血管神经网络。神经元是组成大脑中神经系统结构和产生大脑功能的基本单元，它由胞体和细胞突起组成，其中细胞突起又根据形态和功能的差异被分为轴突(axon)和树突(dendrite)两种类型。大脑的复杂功能依赖于大量的神经元以及神经元之间相互连接所形成的神经环路。神经元突起的形态变化影响着神经元的活动与功能，不同神经元的轴突、树突的形态存在着明显差异，因此，研究神经元形态具有重要的意义。
3.神经元形态数据的获取需要依赖多种技术的结合。通过稀疏标记技术和光学成像技术对脑组织进行成像，可以获取亚微米级别分辨率的神经元图像数据。从这些图像数据中获取矢量化的神经元形态数据，则需要神经元形态重建技术。神经元形态重建主要包含三部分内容：神经元胞体的图像分割、神经元突起的重建、神经元突起内部精细结构的识别(如轴突扣结(bouton)和树突棘(spine)等)。其中，神经元突起状重建是神经元形态重建中基础却困难的内容。
4.目前，神经元突起重建的方法主要是通过分析神经元图像的图像特性，然后通过这些图像的特性设计或者选用相应的数学模型分析图像，进而通过点和线相连的方式提取出神经元突起的拓扑结构，称之为“神经元突起路径”，这些提取出的路径即为神经元突起的重建结果，神经元突起重建结果通常采用多段圆柱状模型相连的具有三维空间坐标的和圆柱半径数据的矢量文件表示。
5.目前这些神经元突起重建方法的关注点主要是神经元突起路径的完整性和路径之间连接的正确性，即这些方法主要是将神经元图像中人可以看到的同一个神经元突起的路径，全部完整的描绘出来，以及确保神经元图像中不同神经元突起不会被识别成同一个神经元突起，分支之间的连接关系和人通过图像看到的连接关系一致，却鲜少有人关注神经元突起的形状重建及重建结果的准确性。
6.神经元的形状重建是指，在重建出神经元突起路径的基础之上结合神经元半径等形态学参数重建出能够描述神经元形状的矢量化文件。而重建结果的准确性包含两个方面，第一个是重建出来的神经元突起路径和图像中神经元突起中心线的贴合程度，第二个是结合神经纤维图像和神经元突起路径估算出的神经元半径的精确度。相较于神经突起路径重建，在确保重建结果的准确性的状态下，神经元的形状重建更能够为后续的神经元形态学分析和神经元模拟提供更多不可或缺的关键数据。
7.针对神经元的形状重建结果的准确性中涉及到的神经元突起路径和图像中神经
元突起中心线的贴合程度的问题，由于图像中神经元突起信号和背景信号复杂多变，因此神经元突起信号的完整性和连接的正确性对于没有人工干预的机器自动识别算法来说是非常困难的。然而对于有人工辅助的神经元突起的重建过程来说，虽然神经元突起的完整性和正确性得到了保证，但是重建出来的神经元突起路径的精确度却有待提高。为了提高重建的神经元突起路径的精确度，现有的方法多采用从图像中神经元突起的灰度特性出发校正神经元突起路径。然而，由于样本标记及成像条件等多种因素的影响，图像中神经元突起信号和背景信号并不是一层不变的。这种复杂的信号，最终会导致神经元突起路径出现锯齿化的现象。而这种锯齿化的神经元突起路径和图像中神经元突起的光滑形态有较大差异。
8.针对神经元的形状重建结果的准确性中涉及到神经元突起半径估计的精确度的问题，现有技术的主要思路是通过提取神经元突起路径节点附近的图像灰度特征，然后对图像进行分割并估计其半径。采用这种思路的方法，提取的神经元突起的半径的精确度，依赖于图像质量的优劣。由于生物标记和成像等多种因素的影响，全脑图像数据中的神经元突起的信号强度和背景信号强度不一定能始终如一的保持一致。变化的信号对于严重依赖于图像特征的算法来说，会导致估算出的神经元突起半径的精确度很不稳定。
9.综上所述,现有的神经元形状重建技术较少,而且目前常用的自动及半自动神经元突起重建算法无法保证重建出来的神经元突起路径的完整性和正确性，目前的神经元突起重建算法对神经元突起路径的半径估算较为粗略，在应用于全脑图像数据时，易受到图像质量的影响，导致计算出来的神经元突起路径节点的半径数据精确度不高。


技术实现要素：

10.因此，为解决上述问题，本发明提供了一种神经元形状的数字化重建方法及其系统。
11.本发明是通过以下技术方案实现的：
12.一种神经元形状的数字化重建方法，包括以下步骤：
13.s1：通过线性插值法在神经元突起路径中插入一系列用于描述神经元突起路径的中心线的等间距的节点p0；
14.s2：依次计算以p0为中心点的局部范围内信号强度最强的像素节点p1，并以p1作为p0的校正点；
15.s3：通过角度约束条件和信号强度约束条件对节点p1进行循环检验和校正，得到精细校正后的节点p2，连接各节点p2，得到精细校正后的神经元突起路径；
16.s4：提取神经元突起路径上各节点的正交平面；
17.s5：分割各正交平面，提取出前景信号；
18.s6：估算各正交平面的半径；
19.s7：循环判断计算结果的准确性后，输出计算结果。
20.优选的，步骤s2中，分别在三轴方向上通过确定以p0为中心点的局部范围内所有像素点和高斯核函数的卷积，从而得出以p0为中心点的局部范围内信号强度最强的像素节点p1的三轴坐标，p1坐标的计算公式为：
21.22.其中，gau表示高斯卷积核，g(p0)为p0在三维坐标轴上的函数。
23.优选的，所述角度约束条件为：神经元突起路径节点p1与其父、子节点连线形成的角度大于等于120
°
；所述信号约束条件为：判断神经元突起路径节点p1的信号强度大于其周边邻域以及其父、子节点的信号强度。
24.优选的，步骤s3中，若节点p1同时符合角度约束条件和信号约束条件，则保留结果；若节点p1不符合其中一个或两个条件，则重复步骤s2、s3。
25.优选的，所述步骤s4包括：
26.s41：根据节点p2与其父节点pf、子节点pc的位置关系，得到圆弧c1；
27.s42：根据圆弧c1、节点p2、pf以及pc的拓扑关系，计算出圆弧c1在节点p2处的切向量v1；
28.s43：通过切向量v1得到以v1为法向量的平面a1，所述a1为神经元突起路径节点p2的正交平面。
29.优选的，所述步骤s5包括：
30.s51：提取出平面a1上的所有像素点，得到神经元突起路径节点p2的正交平面图像；
31.s52：增强所述正交平面图像中神经元突起的前景信号和背景信号之间的差距；
32.s53：将所述正交平面图像上的像素点二值化，提取出前景信号；
33.s54：利用图像形态学的开运算分离出干扰信号；
34.s55：使用连通域分析，标记出分离干扰信号后的正交平面a2。
35.优选的，所述步骤s6包括：
36.s61：计算出正交平面a2中前景信号的像素点个数n2；
37.s62：利用图像分辨率和像素点个数n2计算出正交平面a2中所有像素点所占面积s2；
38.s63：通过圆形面积计算公式得到正交平面a2的半径。
39.优选的，所述步骤s7中，循环判断计算结果的准确性，若计算结果准确，则输出计算结果；若计算结果不准确，则重复步骤s2-步骤s6，直至判断计算结果准确后输出。
40.优选的，所述步骤s7中，基于如下公式判断计算结果的准确性：
[0041][0042][0043]
α＝a*cosγo b(1，4)
[0044]
其中，p为使度量值f取得最大值的点，γ为搜索路径上的输入节点p2与其父子节点所形成的两个向量之间的夹角，gp表示路径上的点p的信号强度，β为代表信号强度的权重因子，α为代表角度大小的权重因子，dis表示搜索点距离校正起始点的距离，a、b、c为设定值，γo表示输入的节点p2与其父、子节点所形成的两个向量之间的夹角。
[0045]
本发明还揭示了一种神经元形状的数字化重建系统，包括：
[0046]
神经元突起路径线性插值单元，用于通过线性插值法在神经元突起路径中插入一系列用于描述神经元突起路径的中心线的等间距的节点p0；
[0047]
神经元突起路径初步校正单元，用于依次计算以p0为中心点的局部范围内信号强
度最强的像素节点p1，并以p1作为p0的校正点；
[0048]
神经元突起路径精细校正单元，用于通过角度约束条件和信号强度约束条件对节点p1进行循环检验和校正，用于得到精细校正后的节点p2，连接各节点p2，得到精细校正后的神经元突起路径；
[0049]
路径节点正交面提取单元，用于提取神经元突起路径上各节点的正交平面；
[0050]
路径节点正交面分割单元，用于分割各正交平面，提取出前景信号；
[0051]
半径估算单元，用于估算各正交平面的半径；
[0052]
循环判断单元，用于循环判断计算结果的准确性后，输出计算结果。
[0053]
本发明技术方案的有益效果主要体现在：
[0054]
1、通过初步校正找到神经元突起路径中信号最强的像素点，再使用角度约束条件和信号约束条件对神经元突起路径中的信号最强的像素点进行循环判断、精确校正，利用神经元突起光滑的特性作为约束条件，从而有效避免神经元突起路径出现锯齿化的现象，同时提高了神经元突起路径的精确度。
[0055]
2、利用代表神经元突起路径的中心线的一系列节点实现神经元突起路径的校正和神经元形状的重建，节点的间距可根据重建精度调整，设置灵活，相对于现有方法中采用多段圆柱状模型相连的重建方式，重建结果更加精确。
[0056]
3、以神经元突起路径节点的拓扑结构为基础进行路径节点的半径估算，对图像质量依赖性较小，估算得到的神经元突起路径节点的半径精确度稳定，能够有效提高半径估算的精确度，在全脑三维神经元突起的图像数据中，避免了由于标记和成像等多种因素的影响而导致三维图像的质量波动的问题，可适用于全脑图像数据，解决了传统方法严重依赖全脑图像质量和图像特征的难题。
[0057]
4、使用神经元突起路径节点局部灰度特征估计半径信息重建神经元突起路径的方法，利用灰度线性变化适当增强图像中神经元突起的信号和背景信号之间的差距，再利用二值化提取出前景信号，然后利用图像形态学的开运算分离出干扰信号，降低了背景信号的干扰。
[0058]
5、在计算结果输出前，循环判断校正后的神经元突起路径以及正交平面半径结果，进一步对计算结果进行检验，提高了神经元形状重建结果的准确性。
[0059]
6、在神经元突起路径重建的基础之上，提高神经元突起路径和图像中神经元突起中心线的贴合程度，并结合神经元路径节点的半径信息，达到神经元形状的精准重建。
[0060]
8、全自动运行，无需人工干预，可自行对全部神经元的形态学文件进行处理，在合适的计算资源和图像数据的支撑下，能够极大的节省人力和时间，提高神经元形态数据的质量，可适用的计算平台类型多。
附图说明
[0061]
图1：是本发明所述一种神经元形状的数字化重建方法的流程图；
[0062]
图2：是本发明所述神经元形状的数字化重建系统的程序流程图；
[0063]
图3：是本发明小鼠试验中神经元突起路径示意图；
[0064]
图4：是图3中a部分的放大图；
[0065]
图5：是对图4进行步骤s1操作后的示意图；
[0066]
图6：是对图5进行步骤s2初步校正后的示意图；
[0067]
图7：是对图6进行步骤s3精细校正后的示意图；
[0068]
图8：是确定图6中各节点p2处切向量v1的示意图；
[0069]
图9：是节点p2的正交平面a1示意图；
[0070]
图10：是正交平面a1灰度线性变化后的示意图；
[0071]
图11：是对图10进行二值化操作后的示意图；
[0072]
图12：是对图11进行图像形态学的开运算操作后的示意图；
[0073]
图13：是对图12中的正交平面的半径估算操作示意图。
具体实施方式
[0074]
为使本发明的目的、优点和特点能够更加清楚、详细地展示，将通过下面优选实施例的非限制性说明进行图示和解释。该实施例仅是应用本发明技术方案的典型范例，凡采取等同替换或者等效变换而形成的技术方案，均落在本发明要求保护的范围之内。
[0075]
同时声明，在方案的描述中，需要说明的是，术语“中心”、“上”、“下”、“左”、“右”、“前”、“后”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。
[0076]
此外，本方案中的术语“第一”、“第二”仅用于描述目的，而不能理解为指示或者暗示对重要性的排序，或者隐含指明所示的技术特征的数量。因此，限定有“第一”、“第二”的特征可以明示或者隐含地包括一个或者更多个该特征。
[0077]
在本发明中，“多个”的含义是两个或者两个以上，除非另有明确具体的限定。
[0078]
本发明揭示了一种神经元形状的数字化重建方法，将该方法应用在fmost成像，gfp荧光标记染色，图像分辨率为0.2
×
0.2
×
1μm3的小鼠全脑三维图像数据集上，结合gtree软件半自动重建单个长程投射神经元神经元突起路径为例得到附图3-图13。
[0079]
如图1所示，本发明的一种神经元形状的数字化重建方法，包括以下步骤：
[0080]
s1：通过线性插值法在神经元突起路径中插入一系列用于描述神经元突起路径的中心线的等间距的节点p0；
[0081]
如图3-图5所示，以图3中的a部分为例，在神经元突起路径中插入一系列用于描述神经元突起路径的中心线的等间距的节点p0，为后续校正提供多个计算样本。
[0082]
具体地，在本实施例中，所述线性插值法是指：先在神经元突起路径中预设一最小距离d，当某一对父子节点之间距离大于d时，在两者之间等间距地插入若干节点p0，且扩充节点p0的类型默认与两骨架点中的父节点一致；为避免相邻的骨架点在原始数据中对应同一体素，因此最小距离d的选取主要与原始数据的体素分辨率有关，例如，若使用原始数据的分辨率为0.2μm
×
0.2μm
×
1μm，基于上述原则，则最小距离d应设为2μm，为轴向分辨率的两倍。
[0083]
s2：依次计算以p0为中心点的局部范围内信号强度最强的像素节点p1，并以p1作为p0的校正点；
[0084]
如图6所示，在一实施例中，通过均值漂移法计算以p0为中心点的局部范围内信号强度最强的像素节点p1，对神经元突出路径进行初步校正，对比图5、图6可以明显看出，代
表神经元突出路径的中心线的变化，具体地，分别在三轴方向上通过确定以p0为中心点的局部范围内所有像素点和高斯核函数的卷积，从而得出以p0为中心点的局部范围内信号强度最强的像素节点p1的三轴坐标，p1坐标的计算公式为：
[0085][0086]
其中，gau表示高斯卷积核，g(p0)为p0在三维坐标轴上的函数。
[0087]
具体地，g(p)表示以点p为中心的局部图像的函数，p1为函数g(p0)中能够得到最大输出值的参数，利用上述公式计算出以p0为中心点的局部范围内，分别在x轴方向、y轴方向和z轴方向上的灰度值最大值(在计算其中一个方向的灰度最大值时，另外两个方向的坐标不变)，在三轴方向上取得灰度值最大值的相应坐标值分别为信号强度最强的像素节点p1的三轴坐标值，在确定节点p1的三轴坐标后，将p1设置为p0的校正点。
[0088]
其中，考虑到神经纤维的直径和轴向分辨率等因素，本发明中轴向坐标移动区间半宽初始值选为2，即假设点p的初始坐标为(x,y,z)，那么在初始条件下轴向的坐标区间为[z-2,z 2]
[0089]
s3：通过角度约束条件和信号强度约束条件对节点p1进行循环检验和校正，得到精细校正后的节点p2，连接各节点p2，得到精细校正后的神经元突起路径；
[0090]
其中，当节点p1与其父、子节点连线形成的角度小于120
°
时，这一骨架点有可能偏离了正确的位置，因此所述角度约束条件为：神经元突起路径节点p1与其父、子节点连线形成的角度大于等于120
°
；所述信号约束条件为：判断神经元突起路径节点p1的信号强度大于其周边邻域以及其父、子节点的信号强度；需注意，起始点和末端点分别缺少父节点和子节点，所以不对起始点和末端点进行校正。
[0091]
如图7所示，在本实施例中，若节点p1同时符合角度约束条件和信号约束条件，则保留结果；若节点p1不符合其中一个或两个条件，则重复步骤s2、s3，直至获得同时符合角度约束条件和信号约束条件的节点p2，并将节点p2设置为p1的校正点，连接各节点p2，得到与图像中神经元突起中心线高度贴合的神经元突起路径，对比图6、图7，代表神经元突出路径的中心线随着各节点p2对节点p1的精细校正，进一步发生变化，循环校正中心线上的各节点p1，使中心线与神经元突起路径图像的贴合程度提高。
[0092]
s4：提取神经元突起路径上各节点的正交平面；
[0093]
在一实施例中，结合精细校正后的神经元突起路径节点p2及其父子节点所在的所在位置的拓扑关系，计算得到在p2点处沿图像中神经元纤维中心线的切向量，以该切向量作为法向量的平面即为所提取的神经元突起路径节点p2的正交平面a1。
[0094]
其中包括以下步骤：
[0095]
s41：根据节点p2与其父节点pf、子节点pc的位置关系，得到圆弧c1；
[0096]
其中，对于校正之后的神经元突起路径节点p2，其父节点pf和子节点pc，可以确定一个半径非常大的内接圆，使得这三点同时存在于该圆弧c1上。
[0097]
s42：根据圆弧c1、节点p2、pf以及pc的拓扑关系，计算出圆弧c1在节点p2处的切向量v1；
[0098]
如图8所示，在确定节点p2的父节点pf和子节点pc以及圆弧c1后，在圆弧c1上节点p2处找到切向量v1。
[0099]
s43：通过切向量v1得到以v1为法向量的平面a1，所述a1为神经元突起路径节点p2的
正交平面。
[0100]
如图9所示，为节点p2处的正交平面a1，此时正交平面a1的图像较为模糊，不能明显分辨前景信号和背景信号。
[0101]
s5：分割各正交平面，提取出前景信号；
[0102]
在一实施例中，首先对正交平面a1进行图像的灰度线性变换处理，对神经元突起路径节点p2的像素信号增强处理，然后通过自适应图像二值化方法，提取出前景信号，接着使用图像形态学的开运算，去除干扰噪声及信号，最后使用连通域分析，分割出神经元突起的圆柱状结构像素点所构成的神经元突起路径节点正交平面a2。
[0103]
其中包括以下步骤：
[0104]
s51：提取出平面a1上的所有像素点，得到神经元突起路径节点p2的正交平面图像；
[0105]
根据计算得到的a1的代数式，基于全脑三维图像数据，可以提取每一个神经元突起路径节点p2的正交平面a1上的所有像素点，得到每一个神经元突起路径节点p2点处的正交平面图像。
[0106]
s52：增强所述正交平面图像中神经元突起的前景信号和背景信号之间的差距；如图10所示，使用图像的灰度线性变换，对神经元突起路径节点p2的像素信号增强处理，对比图9、图10，此时图10中神经元突起的信号和背景信号之间的差距与图9相比明显增强。
[0107]
s53：将所述正交平面图像上的像素点二值化，提取出前景信号；
[0108]
如图11所示，使用自适应图像二值化方法，将图像二值化，整个图像呈现出明显的只有黑和白的视觉效果，进一步突出并提取出前景信号。
[0109]
s54：利用图像形态学的开运算分离出干扰信号；
[0110]
其中，图像形态学中的开运算相当于对图像先进行腐蚀运算再进行膨胀运算，可以消除离散点和"毛刺"，可以将两个物体分开，在本实施例中，如图12所示，利用开运算操作将突起纤维圆柱状结构所代表的前景信号和其外围的干扰信号分离。
[0111]
s55：使用连通域分析，标记出分离干扰信号后的正交平面a2。
[0112]
s6：估算各正交平面的半径；
[0113]
如图13所示，在一实施例中，计算出正交平面a2的前景信号的像素点个数n2，根据图像分辨率的相关参数，计算得到n2个像素点对应的面积s2，使用圆形近似突起纤维的正交平面形状，根据圆形面积计算公式，可以推导出，其半径的估算值为:
[0114]
其中包括以下步骤：
[0115]
s61：计算出正交平面a2中前景信号的像素点个数n2；
[0116]
s62：利用图像分辨率和像素点个数n2计算出正交平面a2中所有像素点所占面积s2；
[0117]
s63：通过圆形面积计算公式得到正交平面a2的半径。
[0118]
s7：循环判断计算结果的准确性后，输出计算结果。
[0119]
所述步骤s7中，循环判断计算结果的准确性，若计算结果准确，则输出计算结果；若计算结果不准确，则重复步骤s2-步骤s6，直至判断计算结果准确后输出。
[0120]
所述步骤s7中，基于如下公式判断计算结果的准确性：
[0121]
[0122][0123]
α＝a*cosγo b(1，4)
[0124]
其中，p为使度量值f取得最大值的点，γ为搜索路径上的输入节点p2与其父子节点所形成的两个向量之间的夹角，gp表示路径上的点p的信号强度，β为代表信号强度的权重因子，α为代表角度大小的权重因子，dis表示搜索点距离校正起始点的距离，a、b、c为设定值，γo表示输入的节点p2与其父、子节点所形成的两个向量之间的夹角。
[0125]
具体地，将精细校正后的神经元突起路径节点p2输入公式(1，2)，若节点p2为搜索路径上使度量值f取得最大值的点，则判断节点p2为最终校正所得的节点，即节点p2的计算结果准确。
[0126]
具体地，γ表示搜索路径上的点p与父子节点形成的两个向量之间的夹角，γ越大则所搜索的路径越平滑。
[0127]
具体地，在搜索过程中，β越大，代表信号强度对度量值f的影响越大，当go小于基准信号强度时，意味着输入的节点p2很可能偏离了前景信号，在校正搜索过程中，需要更多地考虑信号强度变化对f的影响，反之若go大于基准信号强度，说明该骨架点已经位于前景信号内，此时信号强度变化对f的影响小；γo表示输入的节点p2与其父、子节点所形成的两个向量之间的夹角，当γo越小，α越大，此时角度对度量值f的影响越大，说明该骨架点可能在之前的校正过程中受到其他纤维信号的影响，偏离了正确位置，所以在搜索过程中需要更多地考虑角度的影响。
[0128]
具体地，a、b是两个经验因子，a一般取值为0.4，b一般取值为0.7，c是小于1的参数，一般选择为0.98。
[0129]
如图2所示，本发明还揭示了一种神经元形状的数字化重建系统，包括：
[0130]
神经元突起路径线性插值单元，用于通过线性插值法在神经元突起路径中插入一系列用于描述神经元突起路径的中心线的等间距的节点p0；
[0131]
神经元突起路径初步校正单元，用于依次计算以p0为中心点的局部范围内信号强度最强的像素节点p1，并以p1作为p0的校正点；
[0132]
神经元突起路径精细校正单元，用于通过角度约束条件和信号强度约束条件对节点p1进行循环检验和校正，用于得到精细校正后的节点p2，连接各节点p2，得到精细校正后的神经元突起路径；
[0133]
路径节点正交面提取单元，用于提取神经元突起路径上各节点的正交平面；
[0134]
路径节点正交面分割单元，用于分割各正交平面，提取出前景信号；
[0135]
半径估算单元，用于估算各正交平面的半径；
[0136]
循环判断单元，用于循环判断计算结果的准确性后，输出计算结果。
[0137]
本发明尚有多种实施方式，凡采用等同变换或者等效变换而形成的所有技术方案，均落在本发明的保护范围之内。