CN104331892A - 一种基于形态学的神经元识别与分析方法 - Google Patents

Abstract

本发明提出一种基于形态学的神经元识别与分析方法。采集神经细胞的显微图像数据；对显微图像数据进行预处理得到去噪声后的图像；对去除噪声的图像进行阈值分割得到去除背景后的神经元图像；从神经元图像中提取单个神经元；提取单个神经元的基础参数，基础参数包括骨架、细胞体、细胞体面积轴突或树突的数量、轴突或树突的长度。本发明可以自动识别和分析显微图像中的神经元，且识别速度快、效果理想，不仅能直观看出神经元形状，还能输出有效的神经元基础参数，大大减轻了研究者阅片的劳动强度，同时利于后续数据分析。

Description

一种基于形态学的神经元识别与分析方法

技术领域

本发明属于医学图像处理技术领域，具体涉及一种基于形态学的神经元识别与分析方法。

背景技术

神经元是组成神经系统的基本结构和功能单位，也称神经细胞，包括神经细胞体和神经纤维。神经纤维是由神经元的轴突或树突、髓鞘和神经膜组成。髓鞘是由髓磷脂和蛋白质组成，包在轴突或树突的外面，有绝缘作用，神经膜是一种神经胶质细胞，呈薄膜状，包在神经纤维外面，具有保护和再生的作用。轴突或树突是神经元的信号输出通道，轴突通常较树突细，是从细胞体发出的一根较长的呈圆柱形的细长突起，每个神经元只有一个轴突，一般自细胞体发出，与一个或多个目标神经元发生连接。

在当前有关人和动物的神经系统研究中，神经显微形态分析被广泛应用于正常的或病理状态下的神经科学研究或诊断中。神经元及轴突的形态分析，对于对正确理解神经发育、损伤程度及规律，掌握神经形态特征规律及分布具有重要意义。在神经元形态研究中，单项形态参数并不具有实质性的研究参考价值，需进行多项基础参数及其互相关性的研究。这些基础参数包括:纤维或轴突数量、周长和神经细胞体面积等。在以上获得的基础参数的基础上，进而可以计算出其它有较高研究价值的参数:髓鞘厚度与轴突直径之比；轴突直径与纤维直径之比(g-ratio)。

复旦大学博士学位论文《视神经纤维自动识别与分析关键技术研究》提出了一种视神经纤维的识别与分析方法，包括图像预处理、纤维轮廓识别以及形态数据分析这三个步骤。该方法中采用区域生长法粗分割图像，采用snake精确分割神经纤维，虽能自动识别神经纤维，但是识别速度较慢，并且输出的形态数据局限于纤维直径。

发明内容

本发明提出一种基于形态学的神经元识别与分析方法，可以自动识别和分析显微图像中的神经元，且识别速度快、效果理想，不仅能直观看出神经元形状，还能输出有效的神经元基础参数，大大减轻了研究者阅片的劳动强度，同时利于后续数据分析。

为了解决上述技术问题，本发明提供一种基于形态学的神经元的识别与分析方法，步骤如下：

第一步，采集神经细胞的显微图像数据；

第二步，对第一步得到的显微图像数据进行预处理，得到去噪声后的图像；

第三步，对第二步去除噪声的图像进行阈值分割，得到去除背景后的神经元图像；

第四步，从神经元图像中提取单个神经元；

第五步，提取单个神经元的基础参数，基础参数包括骨架、细胞体、细胞体面积轴突或树突的数量、轴突或树突的长度。

较佳地，第二步中采用中值滤波方法对显微图像进行去噪处理。

较佳地，第三步中，采用自适应阈值函数T(x,y)进行阈值分割，方法如公式(1)所示：

$g(x,y)=\left\{\begin{array}{c}g(x,y),g(x,y)>T(x,y)\\ 0,g(x,y)\≤T(x,y)\end{array}\right.---\left(1\right)$

公式(1)中，自适应阈值函数T(x,y)＝g₀(x,y)+T₀；图像g₀(x,y)是对灰度值g(x,y)的局部补偿，常数T₀是图像g₀(x,y)的全局阈值，采用Otsu方法确定最佳全局阈值T₀。

较佳地，第四步中，采用八邻域算法识别并标记各个连通域，所识别的连通域即为神经元。

进一步，第五步中，所述骨架提取方法为：设单个神经元连通域为集合A，则集合A的骨架S(A)用腐蚀和开操作表达，表达式如公式(2)所示，

$S\left(A\right)=\underset{k=0}{\overset{K}{\∪}}{S}_k\left(A\right)---\left(2\right)$

式(2)中，腐蚀和开操作k次后的骨架S_k(A)＝(AΘkB)-(AΘkB)οB，B是一个结构元，(AΘkB)表示对集合A的连续k次腐蚀，(AΘkB)如公式(3)所示，

(AΘkB)＝((…(AΘB)ΘB)Θ…)ΘB)          (3)

式(2)中，K是集合A被腐蚀成空集前的最后一次迭代步骤，且

进一步，第五步中，所述细胞体提取方法为：

首先，提取连通域的质心，设单个神经元图像中像素点(x,y)的灰度值为g(x,y)，则连通域质心如公式(4)所示，

$\hat{x}=\frac{\mathrm{\Σ\Σxg}(x,y)}{\mathrm{\Σ\Σg}(x,y)},\hat{y}=\frac{\mathrm{\Σ\Σyg}(x,y)}{\mathrm{\Σ\Σg}(x,y)}---\left(4\right)$

其次，对连通域质心进行单位膨胀；

假设C为神经元图像中的集合，D为神经元图像中一个圆形结构元，D对C进行膨胀的方式如公式(5)所示，

$C\⊕D=\left\{z\right|[{\left(\hat{D}\right)}_z\∩C]\⊆C\}.---\left(5\right)$

公式(5)中，为膨胀后的图像，D的圆心为质心半径为一个像素，初次膨胀时集合C为质心

再次，质心被结构元单位膨胀后，比较单位膨胀操作获得的集合与神经元连通域为集合A的重叠性，若两者之间没有重叠点，则用集合替代单位膨胀操作中的集合C并重复单位膨胀，直至集合与集合A恰好有一个重叠点；

最后，以质心作为圆心，以质心到所述重叠点的距离作为半径，构建圆形结构元E，用结构元E对单个神经元连通域集合A进行开操作，如公式(6)所示，开操作结果AοE即为神经元的细胞体图像，

进一步，第五步中，所述细胞体面积提取方法为：将所述细胞体图像中的像素数量乘以显微镜的放大系数获得细胞体的面积。

进一步，第五步中，所述轴突或树突数量的提取方法为：将所述获得的骨架图像与细胞体图像重叠，骨架和细胞体的交点记为轴突或树突的起点，并将交点的数量作为轴突或树突的数量参数输出。

进一步，所述轴突或树突长度的提取方法为：

首先，以一个轴突或树突的起点为起点，沿骨架轨迹逐邻域判断像素数量标定属于该轴突或树突的支点，直至骨架的终点；当八邻域内出现三个及三个以上像素不为0的像素点时，则邻域中心为该轴突或树突的支点；当八邻域内只有一个像素不为0的像素点时，则邻域中心为该轴突或树突的终点；

然后，以一个轴突或树突的起点为起点，沿骨架轨迹，遍历属于该轴突或树突的支点，直至相应的骨架终点，使用递推函数计算起点到支点、支点到支点、支点到终点的像素总数，再乘以显微镜的放大系数作为该轴突或树突的长度参数输出。

本发明与现有技术相比，其显著优点在于，(1)本发明使用结合了Otsu确定最佳阈值的自适应阈值分割方法，能够快速有效地识别神经元；(2)本发明结合形态学算法直接获得细胞体图像，便于直直观观察，同时还可以获得包含胞体面积、轴突或树突数量、轴突或树突长度的基础参数，有利于后续参数分析。

附图说明

图1是本发明基于形态学的神经元识别与分析方法的流程示意图。

图2是本发明中单个神经元基础参数的提取流程图。

图3是本发明中阈值分割后多个神经元的图像。

图4是本发明中单个神经元的骨架图像。

图5是本发明中单个神经元的细胞体图像。

图6是本发明中包含骨架的单个神经元图像，其中s1、s2、s3、s4是轴突或树突的起点，e1、e2、e3、e4是对应轴突或树突的终点，由于该图像中轴突或树突无支点，则s1-e1为一个完整的轴突或树突。

具体实施方式

如图1所示，本发明基于形态学的神经元的识别与分析方法，步骤如下：

第一步，采集神经细胞的显微图像数据；

第二步，对第一步得到的显微图像数据进行预处理，得到去噪后的图像。

电噪声、外界的轻微震动等干扰会影响显微图像质量，因此图像的滤波去噪是必要的，常用的空间滤波方法有均值滤波、中值滤波和自适应滤波等。根据髓鞘和轴突的灰度特性，本发明采用中值滤波就能达到比较理想的效果。中值滤波使用一个像素领域中的灰度级中值来替代该像素值，常用的二维的领域窗口如：方形、圆形、十字形、菱形等，通常领域内像素个数为奇数。

第三步，对第二步去除噪声的图像进行阈值分割，得到去除背景后的神经元图像。

阈值处理通过设定不同的特征阈值，把图像像素分为若干类，去除背景突显目标，常用的特征包括：直接来自原始图像的灰度或彩色特征、由原始灰度或彩色值变换得到的特征。在许多情况下，采用的特征随图像位置的变化而变化，因此，本发明考虑到实际处理图像，采用自适应阈值函数T(x,y)对像素点(x，y)灰度值为g(x,y)的图像进行阈值处理，方法如公式(1)所示：

$g(x,y)=\left\{\begin{array}{c}g(x,y),g(x,y)>T(x,y)\\ 0,g(x,y)\≤T(x,y)\end{array}\right.---\left(1\right)$

公式(1)中，T(x,y)＝g₀(x,y)+T₀；图像g₀(x,y)是对g(x,y)的局部补偿，常数T₀是图像g₀(x,y)的全局阈值。本发明采用Otsu方法确定最佳全局阈值T₀，原理详见电子工业出版社出版的《数字图像处理(第三版)》。阈值分割后的图像如图3所示。

第四步，提取单个神经元。

第三步阈值分割后的图像中出现众多包含神经纤维和神经细胞体的连通域，根据灰度特性，采用八邻域算法识别并标记各个连通域，所识别的连通域即为本发明所需的神经元。

第五步，提取单个神经元的基础参数。下面结合图2详细说明针对单个神经元可提取的参数以及方法：

(1)提取骨架。

设单个神经元连通域为集合A，则集合A的骨架S(A)用腐蚀和开操作表达，表达式如公式(2)所示，

$S\left(A\right)=\underset{k=0}{\overset{K}{\∪}}{S}_k\left(A\right)---\left(2\right)$

公式(2)中，腐蚀和开操作k次后的骨架S_k(A)＝(AΘkB)-(AΘkB)οB，B是一个结构元，(AΘkB)表示对集合A的连续k次腐蚀，(AΘkB)如公式(3)所示，

(AΘkB)＝((…(AΘB)ΘB)Θ…)ΘB)           (3)

K是集合A被腐蚀成空集前的最后一次迭代步骤，即：提取骨架后的图像如图4所示。

(2)提取细胞体。

首先，提取连通域的质心。

设单个神经元图像中像素点(x,y)的灰度值为g(x,y)，则连通域质心如公式(4)所示，

$\hat{x}=\frac{\mathrm{\Σ\Σxg}(x,y)}{\mathrm{\Σ\Σg}(x,y)},\hat{y}=\frac{\mathrm{\Σ\Σyg}(x,y)}{\mathrm{\Σ\Σg}(x,y)}---\left(4\right)$

其次，对连通域质心进行单位膨胀。

假设C为神经元图像中的集合，D为神经元图像中一个圆形结构元，则D对C进行膨胀的方式如公式(5)所示，

$C\⊕D=\left\{z\right|[{\left(\hat{D}\right)}_z\∩C]\⊆C\}.---\left(5\right)$

公式(5)中，为膨胀后的图像，D的圆心为质心半径为一个像素，初次膨胀时集合C为质心

再次，确定重叠点。

质心被结构元单位膨胀后，比较单位膨胀操作获得的集合与神经元连通域为集合A的重叠性，若两者之间没有重叠点，则用集合替代单位膨胀操作中的集合C并重复单位膨胀，直至集合与集合A恰好有一个重叠点。

最后，对连通域进行开操作。

以质心作为圆心，以质心到所述重叠点的距离作为半径，构建圆形结构元E，则结构元E对单个神经元连通域集合A的开操作结果AοE如公式(6)所示，开操作结果AοE即为神经元的细胞体图像，如图5所示。

(3)提取细胞体的面积参数。

将所述细胞体图像中的像素数量乘以显微镜的放大系数，即可获得细胞体的面积参数。

(4)提取轴突或树突的数量参数。

将前述获得的骨架图像与细胞体图像重叠，骨架和细胞体的交点记为轴突或树突的起点(如图6中的s1、s2、s3、s4)，并将交点的数量作为轴突或树突的数量参数输出。

(5)提取轴突或树突的长度参数。

首先，确定轴突或树突的支点和终点。

以第五步的(4)中的一个轴突或树突的起点为起点，沿第五步的(1)图像中的骨架轨迹，逐邻域判断像素数量标定属于该轴突或树突的支点，直至骨架的终点(如图6中的e1、e2、e3、e4)。当八邻域内出现三个及三个以上像素不为0的像素点时，则邻域中心为该轴突或树突的支点；当八邻域内只有一个像素不为0的像素点时，则邻域中心为该轴突或树突的终点，即骨架的终点。

其次，计算轴突或树突长度。

以第五步的(4)中的一个轴突或树突的起点为起点，沿骨架轨迹，遍历属于该轴突或树突的支点，直至相应的骨架终点，使用递推函数计算起点到支点、支点到支点、支点到终点的像素总数，再乘以显微镜的放大系数作为该轴突或树突的长度参数输出。特例，若一个轴突或树突无支点，如图6所示，则计算起点s1到终点e1的像素数量，再乘以显微镜的放大系数作为该轴突或树突s1-e1的长度参数输出。

第六步，待第四步提取的全部神经元均按照第五步提取完基础参数后，输出神经元的基础参数。

Claims (9)

1.一种基于形态学的神经元的识别与分析方法，其特征在于，步骤如下：

第一步，采集神经细胞的显微图像数据；

第二步，对第一步得到的显微图像数据进行预处理，得到去噪声后的图像；

第三步，对第二步去除噪声的图像进行阈值分割，得到去除背景后的神经元图像；

第四步，从神经元图像中提取单个神经元；

第五步，提取单个神经元的基础参数，基础参数包括骨架、细胞体、细胞体面积轴突或树突的数量、轴突或树突的长度。

2.如权利要求1所述基于形态学的神经元的识别与分析方法，其特征在于，第二步中采用中值滤波方法对显微图像进行去噪处理。

3.如权利要求1所述基于形态学的神经元的识别与分析方法，其特征在于，第三步中，采用自适应阈值函数T(x,y)进行阈值分割，方法如公式(1)所示：

$g(x,y)=\left\{\begin{array}{c}g(x,y),g(x,y)>T(x,y)\\ 0,g(x,y)\≤T(x,y)\end{array}\right.---\left(1\right)$

公式(1)中，自适应阈值函数T(x,y)＝g₀(x,y)+T₀；图像g₀(x,y)是对灰度值g(x,y)的局部补偿，常数T₀是图像g₀(x,y)的全局阈值，采用Otsu方法确定最佳全局阈值T₀。

4.如权利要求1所述基于形态学的神经元的识别与分析方法，其特征在于，第四步中，采用八邻域算法识别并标记各个连通域，所识别的连通域即为神经元。

5.如权利要求1所述基于形态学的神经元的识别与分析方法，其特征在于，第五步中，所述骨架提取方法为：设单个神经元连通域为集合A，则集合A的骨架S(A)用腐蚀和开操作表达，表达式如公式(2)所示，

$S\left(A\right)=\underset{k=0}{\overset{K}{\∪}}{S}_k\left(A\right)---\left(2\right)$

式(2)中，腐蚀和开操作k次后的骨架S_k(A)＝(AΘkB)-(AΘkB)оB，B是一个结构元，(AΘkB)表示对集合A的连续k次腐蚀，(AΘkB)如公式(3)所示，

(AΘkB)＝((…(AΘB)ΘB)Θ…)ΘB)     (3)

式(2)中，K是集合A被腐蚀成空集前的最后一次迭代步骤，且

6.如权利要求5所述基于形态学的神经元的识别与分析方法，其特征在于，第五步中，所述细胞体提取方法为：

6.1提取连通域的质心，设单个神经元图像中像素点(x,y)的灰度值为g(x,y)，则连通域质心如公式(4)所示，

$\hat{x}=\frac{\mathrm{\Σ\Σxg}(x,y)}{\mathrm{\Σ\Σg}(x,y)},\hat{y}=\frac{\mathrm{\Σ\Σyg}(x,y)}{\mathrm{\Σ\Σg}(x,y)}---\left(4\right)$

6.2对连通域质心进行单位膨胀；

假设C为神经元图像中的集合，D为神经元图像中一个圆形结构元，D对C进行膨胀的方式如公式(5)所示，

$C\⊕D=\left\{z\right|[{\left(\hat{D}\right)}_z\∩C]\⊇C\}.---\left(5\right)$

公式(5)中，为膨胀后的图像，D的圆心为质心半径为一个像素，初次膨胀时集合C为质心

6.3质心被结构元单位膨胀后，比较单位膨胀操作获得的集合与神经元连通域为集合A的重叠性，若两者之间没有重叠点，则用集合替代单位膨胀操作中的集合C并重复单位膨胀，直至集合与集合A恰好有一个重叠点；

6.4以质心作为圆心，以质心到所述重叠点的距离作为半径，构建圆形结构元E，用结构元E对单个神经元连通域集合A进行开操作，如公式(6)所示，开操作结果AоE即为神经元的细胞体图像，

7.如权利要求5所述基于形态学的神经元的识别与分析方法，其特征在于，第五步中，

所述细胞体面积提取方法为：将所述细胞体图像中的像素数量乘以显微镜的放大系数获得细胞体的面积。

8.如权利要求5所述基于形态学的神经元的识别与分析方法，其特征在于，第五步中，所述轴突或树突数量的提取方法为：将所述获得的骨架图像与细胞体图像重叠，骨架和细胞体的交点记为轴突或树突的起点，并将交点的数量作为轴突或树突的数量参数输出。

9.如权利要求5所述基于形态学的神经元的识别与分析方法，其特征在于，第五步中，所述轴突或树突长度的提取方法为：

9.1以一个轴突或树突的起点为起点，沿骨架轨迹逐邻域判断像素数量标定属于该轴突或树突的支点，直至骨架的终点；当八邻域内出现三个及三个以上像素不为0的像素点时，则邻域中心为该轴突或树突的支点；当八邻域内只有一个像素不为0的像素点时，则邻域中心为该轴突或树突的终点；

9.2以一个轴突或树突的起点为起点，沿骨架轨迹，遍历属于该轴突或树突的支点，直至相应的骨架终点，使用递推函数计算起点到支点、支点到支点、支点到终点的像素总数，再乘以显微镜的放大系数作为该轴突或树突的长度参数输出。