基于电镜序列切片图像的突触检测方法与流程

本发明涉及模式识别和神经学领域，具体是涉及一种基于电镜序列切片图像的突触检测方法。

背景技术：

近年来，神经科学是研究的热点，其中突触是重要的研究对象。突触是神经元之间进行信息传递和进行信息整合的关键部位，是神经元之间相互联系的基本单元，由突触前膜、突触间隙和突触后膜组成。突触不仅是研究神经系统正常电生理和功能的基础,同时也是研究神经系统疾病机制的突破口，因此对突触进行研究具有十分重要的意义。

目前基于电镜序列切片大规模获取神经元图像已经不再是难题，而从大体量数据中检测突触很大程度上制约了对突触进行研究。

有鉴于此，特提出本发明。

技术实现要素：

为了解决现有技术中存在的技术问题，本发明实施例提供一种基于电镜序列切片图像的突触检测方法，以实现高效、精准的突触检测。

为了实现上述目的，一方面，提供一下技术方案：

一种基于电镜序列切片图像的突触检测方法，所述方法包括：

获取所述电镜序列切片图像；

对所述电镜序列切片图像进行拼接对准；

针对拼接对准后的单层电镜序列切片图像，利用Adaboost算法进行突触检测；

基于突触检测结果，进行基于突触特征的突触识别；

基于突触识别的结果，进行基于上下文线索的突触过滤；

基于突触过滤的结果，进行突触的三维重构。

进一步地，所述针对拼接对准后的单层电镜序列切片图像，利用Adaboost算法进行突触检测，具体包括：

创建训练样本；

利用所述训练样本的Haar特征来训练分类器，得到一个级联的boosted分类器；

利用所述级联的boosted分类器对所述突触进行识别检测，同时获取所述突触的位置信息。

进一步地，所述基于突触检测结果，进行基于突触特征的突触识别具体包括：

对所述突触的局部图像进行基于形态学的图像分割，得到突触致密带的备选位置；

基于所述突触致密带的备选位置，应用动态规划方法，将所述突触的局部图像分割为突触前、突触致密带、突触后；

根据所述突触前、所述突触后和所述突触致密带分别构造对应的统计特征；

对所述突触前中的囊泡进行基于Adaboost算法的检测，并根据检测的囊泡构造对应的统计信息；

针对所述突触中的所述囊泡建立正负样本库；

使用随机森林方法进行所述突触识别。

进一步地，所述对所述突触的局部图像进行基于形态学的图像分割，得到突触致密带的备选位置，具体包括：

对所述突触的局部图像进行形态学操作，得到所述突触的粗略结构；

利用二次多项式拟合方法将所述突触的粗略结构拟合成第一曲线；

将拟合出的第一曲线的端点确定为起止点，并使用Dijkstra算法，在所述突触的局部图像上确定第一最短路径；

以所述第一最短路径作为种子点，使用基于稀疏标记的区域增长方法获得所述突触致密带的备选位置。

进一步地，所述基于所述突触致密带的备选位置，应用动态规划方法，对突触前和突触后进行分割，具体包括：

利用二次多项式拟合方法将所述突触的粗略结构拟合成贯穿所述突触的局部图像的第二曲线；

将所述第二曲线的端点确定为起止点，并利用Dijkstra算法在所述突触的局部图像上确定第二最短路径；

对所述第二最短路径进行膨胀处理，则将所述突触的局部图像划分为二区域，并计算所述二区域的像素均值，从而将所述突触的局部图像分割为突触前、突触致密带、突触后。

进一步地，所述根据所述突触前、所述突触后和所述突触致密带分别构造对应的统计特征，具体包括：

分别计算所述突触前、所述突触后和所述突触致密区的均值和方差以及K阶矩，从而得到所述突触前、所述突触后和所述突触致密带的特征向量

进一步地，所述对所述突触前中的囊泡进行基于Adaboost算法的检测，并根据检测的囊泡构造对应的统计信息，具体包括：

通过Adaboost算法检测出囊泡信息，并将所述囊泡的数量、均值、方差以及K阶矩添加至所述特征向量，从而构造所述囊泡的统计信息。

进一步地，所述使用随机森林方法进行所述突触识别具体包括：

步骤A：进行自举重采样；

步骤B：随机选择预定数量的所述特征；

步骤C：遍历所述各特征的每一种分割方式，确定最好的分割点，并分割成两个子节点；

步骤D：分别对每一所述子节点，执行步骤C，直至每一节点达到预定纯度为止。

进一步地，所述两个子节点包括左子节点和右子节点；则所述确定最好的分割点包括：

根据以下公式确定信息量：

其中，所述L(R)表示对当前样本进行观察所得到的信息熵；所述p_i表示分到第i类的样本数量/样本总数量；i＝1，……K；所述K表示样本被分的类数；

根据以下公式确定所述特征的信息增益：

其中，所述s表示所述特征；所述|R|表示当前节点处的样本数目；所述|R_ls|表示根据所述特征s分到所述左子节点的样本个数；所述|R_rs|表示根据所述特征s分到所述右子节点的样本个数；所述L(R_ls)表示根据所述特征s分到所述左子节点的样本的信息熵；所述L(R_rs)表示根据所述特征s分到所述右子节点的样本的信息熵；

将使所述信息增益最大的特征确定为所述最好的分割点。

进一步地，所述基于突触识别的结果，进行基于上下文线索的突触过滤，具体包括：

步骤E：确定利用随机森林方法在所述电镜序列切片图像中第i层图像中检测到的第j个突触R_i,j；其中，所述i表示所述电镜序列切片图像中图像的层数；所述j表示所述突触的序号；

步骤F：确定所述突触R_i,j的中心位置和突触R_i+1,k的的中心位置之间的相对距离；其中，所述k取正整数；

步骤G：计算所述突触R_i,j与第i+1层图像上所有突触之间距离的最小值；

步骤H：如果所述最小值小于所述电镜序列切片图像中上下两层图像之间的最大偏移量，则将所述R_i,j和所述R_i+1,k确定为分布在所述电镜序列切片图像中的同一突触；

步骤I：针对所述电镜序列切片图像中所有层图像中检测到的所有突触，执行步骤F至步骤H，直至步骤G得到的距离最小值大于等于所述电镜序列切片图像中上下两层图像之间的最大偏移量；

判断所述突触R_i,j在所述电镜序列切片图像中出现的次数是否大于预定次数，若是，则确定所述突触R_i,j为真正的突触，否则，确定所述突触R_i,j为伪突触。

本发明提供一种基于电镜序列切片图像的突触检测方法，其对电镜序列切片图像进行拼接对准；接着针对拼接对准后的单层电镜序列切片图像，利用Adaboost算法进行突触检测；然后基于突触特征进行突触识别，最后利用上下文线索进行突触过滤，并基于突触过滤的结果，进行突触的三维重构。解决了大规模的突触检测以及统计的问题。实现了高效、精准的突触检测。本发明在各种图像上做了验证，可以在突触检测方面取得了很好的效果。

附图说明

图1为根据本发明实施例的基于电镜序列切片图像的突触检测方法的流程示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合附图和具体实施例，对本发明进一步详细说明。

本发明实施例提供一种基于电镜序列切片图像的突触检测方法。该方法可以包括：

S100：获取电镜序列切片图像。

其中，电镜序列切片图像可以通过AUTM-SEM得到。

S110：对电镜序列切片图像进行拼接对准。

S120：针对拼接对准后的单层电镜序列切片图像，利用Adaboost算法进行突触检测。

本步骤可以进一步包括：

S121：创建训练样本。

本步骤可以根据自建的突触数据库来创建训练集。例如，可以创建具有755个正样本和2755个负样本的训练集。较佳地，可以将训练样本大小一律归一化为24×24像素。

S122：利用训练样本的Haar特征来训练分类器，得到一个级联的boosted分类器。

本步骤还可以进一步包括：

S1221：构造正负样本。

作为示例，可以用(x₁,y₁),…,(x_n,y_n)来表示样本图像。其中，y_i表示正负样本(i＝1,......n)，y_i＝0表示样本图像为负样本，y_i＝1时表示样本图像为正样本，n表示正负样本的总数。

S1222：初始化每个样本对应的权重。

开始时所有正样本或负样本的权重是相同的，当y_i＝0,1时，样本的初始权重可表示为m、l分别表示负样本和正样本的数量；i＝1......n。

S1223：调整每个样本对应的权重，获得不同权重对应的样本集和弱分类器。

本步骤可以进一步包括：

S12231：对权重进行归一化。

例如，通过以下公式对权重进行归一化：其中，n表示正负样本的总数；i＝1......n；ω_t,i和ω_t,j分别表示第t次循环中样本i和样本j的权重。

S12232：根据以下公式计算权重误差：

ε_t＝min_f,p,θ∑_iω_t,i|h(x_i,f,p,θ)-y_i|

其中，ε_t表示权重误差；h(x_i,f,p,θ)表示弱分类器；f表示一种特征，θ表示阈值，p表示调整不等式方向的参数；x表示检测子窗口；y_i表示正负样本，y_i＝0表示样本图像为负样本，y_i＝1时表示样本图像为正样本。

S12233：根据最小权重误差，得到最优弱分类器。

例如：h_t(x)＝h(x,f_t,p_t,θ_t)，其中，f_t、p_t、θ_t分别表示使权重误差达到最小的特征、调整不等式方向的参数、阈值。

S12234：根据以下公式更新各样本的权重：

如果x_i分类正确，则e_i＝0；否则，e_i＝1。

根据正、负样本的标记(正样本标记为1，负样本标记为0)，对于h_t分类错误的样本，加大其对应的权重；而对于分类正确的样本，降低其权重，这样分错的样本就被突出出来，从而得到一个新的样本分布，在新的样本分布下，再次对弱分类器进行训练，得到弱分类器h_t+1。依次类推，经过T次循环，得到T个弱分类器。

S12235：将各弱分类器按一定的权重叠加起来，构成一个强分类器，并将各强分类器组合成级联的boosted分类器。

S123：利用级联的boosted分类器对突触进行识别检测，同时获取突触的位置信息。

本步骤通过移动搜索窗口来检测突触。搜索窗口依次通过每一级分类器，这样就排除了前面几层的检测中大部分的候选区域。全部通过每一级分类器检测的区域即为目标区域。

为了检测整幅单层电镜图像，在该图像中移动搜索窗口，检测目标区域中每一个位置，以此来确定可能的目标突触。

为了搜索不同大小的突触，可以将分类器设计为尺寸改变型分类器，这样比改变待检图像的尺寸大小更为有效。因而，为了在单层电镜图像中检测未知大小的目标物体，通常需要用不同比例大小的搜索窗口对图片进行几次扫描。优选地，本发明实施例根据对突触实际大小进行估计，设定搜索窗口的最小和最大尺寸分别为(150×150)和(400×400)(图像分辨率为2nm/像素)。

S130：基于突触检测结果，进行基于突触特征的突触识别。

本步骤进一步可以包括：

S131：对突触的局部图像进行基于形态学的图像分割，得到突触致密带的备选位置。

突触的局部图像就是检测到突触的搜索窗口内的图像。

本步骤还可以进一步包括：

S1311：对突触的局部图像进行形态学操作，得到突触的粗略结构。

其中，形态学操作可以包括但不限于二值化、腐蚀、开运算、填充。

S1312：利用二次多项式拟合方法将突触的粗略结构拟合成第一曲线。

例如：可以从突触的粗略结构中选取n个点，然后，利用该n个点拟合出如下曲线：

f(x)＝c₀+c₁x+c₂x²

其中，c₀、c₁、c₂分别表示系数；x表示自变量，其可能是横坐标也可能是纵坐标，根据具体情况来判定。

S1313：将拟合出的第一曲线的端点确定为起止点，并使用Dijkstra算法，在突触的局部图像上确定第一最短路径。

具体地，本步骤可以进一步包括：

S13131：确定起点至终点的距离、结点间的距离及结点集合；

作为示例，起点至终点的距离D(z_t)为从起点z₁沿某一路径到终点z_n的所有结点像素值之和其中，I(z_t)表示从起点Z₁到Z_n经过的结点z_t的像素值，z_t表示从起点Z₁到Z_n经过的结点，t＝1，2，…，n。上述描述可以通过下式表示：

对于结点间的距离。例如：结点间的距离d(z_i,z_j)为结点z_i到结点z_j之间的“距离”即结点z_j的像素值I(z_j)。

对于结点集合。例如：P表示结点集合(也即经过的结点)。对P进行初始化，得到P＝{z₁}(表示其中只包括起点)。

S13132：查找结点集合P以外的点对所有不在P中的点有：

S13133：将距离z_t最近的结点加入P中，将D(z_t+1)更新为min[D(z_t+1),D(z_t)+d(z_t+z_t+1)]。

S13134：重复执行步骤S13132和步骤S13133，直至终点z_n包含在P中，从而在突触的局部图像上确定第一最短路径。

S1314：以突触的局部图像上第一最短路径作为种子点，使用基于稀疏标记的区域增长方法获得突触致密带的备选位置。

本步骤可以使用如下区域增长模型来实现：

其中，fu表示人工给定标记的值，例如：在目标处赋值fu＝-∞，背景处赋值fu＝∞，在其他处赋值fu＝0；u表示分割结果，u:Ω→[0,1]，例如：u≥0.5表示目标，u<0.5表示背景；Ω表示整幅单层电镜图像；g:Ω→[0,1]将Ω中较强边界赋予较小的值。

S132：基于突触致密带的备选位置，应用动态规划方法，将突触的局部图像分割为突触前、突触致密带、突触后。

其中，突触前是指轴突(包含囊泡)，突触后是指树突。

本步骤可以进一步包括：

S1321：利用二次多项式拟合方法将突触的粗略结构拟合成贯穿突触的局部图像的第二曲线。

本步骤可以参考拟合第一曲线的方法，在此不再赘述。

S1322：将第二曲线的端点确定为起止点，并利用Dijkstra算法在突触的局部图像上确定第二最短路径。

本步骤在具体实施时可以参考确定第一最短路径的方法，在此不再赘述。

S1323：对第二最短路径进行膨胀处理，则将突触的局部图像划分为二区域，并计算该二区域的像素均值，通过该二区域的像素均值将突触的局部图像分割为突触前、突触致密带、突触后。

因为突触前分布着稠密的囊泡，因此突触前的像素值均值较小，而突触后像素值均值较大；所以，通过计算突触的局部图像分开的两个区域的像素均值，就达到了对突触前和突触后进行分割的目的。

S133：根据突触前、突触后和突触致密带分别构造对应的统计特征。

具体地，本步骤可以通过分别计算突触前、突触后和突触致密区的均值和方差以及K阶矩(例如：三、四、五、六阶矩)，从而得到突触前、突触后和突触致密带的特征矩阵。

下面以突触前为例来对本步骤进行详细说明：

根据以下公式计算突触前的均值：

其中，i表示突触前的像素点序号；A_i表示突触前像素点i的像素值；N表示突触前像素点的数目；μ表示突触前的均值。

根据以下公式计算突触前的方差：

其中，V表示突触前的方差。

根据以下公式计算突触前的K阶矩：

m_k＝E(A-μ)^k

其中，E(x)表示x的期望值。

S134：对突触前中的囊泡进行基于Adaboost算法的检测，并根据检测的囊泡构造对应的统计信息。

本步骤通过Adaboost算法检测出囊泡信息，并将囊泡的数量、均值、方差以及K阶矩(例如：三、四、五、六阶矩)添加至步骤S133得到的特征向量中。

S135：针对突触中的囊泡建立正负样本库。

S136：使用随机森林方法进行突触识别。

其中，建立决策树的步骤可以包括：

S1361：进行自举重采样。

本步骤在N个训练样本中采取有放回的方式抽取N个样本。将重采样的样本作为新的样本看作一个节点。

S1362：随机选择预定数量的特征。

其中，该特征例如可以为均值和方差以及K阶矩。

S1363：遍历每个特征的每一种分割方式，找到最好的分割点，并分割成两个子节点。

其中，分割方式指的是给定一个阈值，若样本对应特征的取值小于该阈值，则将样本分到左子节点N1；否则，将样本分到右子节点N2。

其中，最好的分割点是根据“纯度”来构建的。这里，计算纯度的方法可以为熵。例如：如果样本被分为K类，且每一类的比例P(j)＝分到第j类的样本数量/样本总数量(j＝1，……K)；则，对该事件的结果进行观察所得到的信息量可以用如下度量来定义：

其中，L(R)表示对当前样本进行观察所得到的信息熵；p_i表示分到第i类的样本数量/样本总数量(i＝1，……K)；K表示样本被分的类数。

对于当前节点上的样本，将上述度量称为熵不纯度，它反映了该节点上的特征对样本分类的不纯度。其中每个特征的信息增益定义如下：

其中，s表示特征；|R|表示当前节点处的样本数目；|R_ls|表示根据特征s分到左子节点N1的样本个数；|R_rs|表示根据特征s分到右子节点N2的样本个数；L(R_ls)表示根据特征s分到左子节点N1的样本的信息熵；L(R_rs)表示根据特征s分到右子节点N2的样本的信息熵。

每个特征的信息增益也就是当前节点的不纯度减去子节点不纯度的加权平均数，权重由子节点记录数与当前节点记录数的比例决定。如果特征s的信息增益最大，这表明特征s能够为这个分类系统带来更多的信息，因此选取s作为当前节点，也即作为最好的分割点。

S1364：分别对每个子节点N1和N2执行步骤S1363，直到每个节点达到预定纯度为止。

S140：基于突触识别的结果，进行基于上下文线索的突触过滤。

本步骤中，上下文线索是指电镜序列切片图像中相邻图像之间的上下层关系。

本步骤具体可以包括：

S141：确定利用随机森林方法在电镜序列切片图像中第i层图像中检测到的第j个突触R_i,j。其中，i表示电镜序列切片图像中图像的层数；j表示突触的序号；i和j分别取正整数。

S142：确定突触R_i,j的中心位置和突触R_i+1,k的的中心位置之间的相对距离。其中，k取正整数。

S143：计算突触R_i,j与第i+1层图像上所有突触之间距离的最小值。

S144：如果该最小值小于电镜序列切片图像中上下两层图像之间的最大偏移量，则将R_i,j和R_i+1,k确定分布在电镜序列切片图像中的同一突触。

S145：针对电镜序列切片图像中所有层图像中检测到的所有突触，执行步骤S142至步骤S144，直至步骤S143计算得到的距离最小值大于等于电镜序列切片图像中上下两层图像之间的最大偏移量。

S146：判断突触R_i,j在电镜序列切片图像中出现的次数是否大于预定次数，若是，则确定突触R_i,j为真正的突触；否则，确定突触R_i,j为伪突触。

下面以一具体实施例来对基于突触识别的结果进行基于上下文线索的突触过滤的过程进行详细说明。

步骤A：将利用随机森林方法在电镜序列切片图像中第i层图像中检测到的第j个突触记为R_i,j。将R_i,j所在检测窗口的左上角和右下角坐标分别记为和

步骤B：根据检测到的突触与检测窗口一一对应的关系，将R_i,j公式化为：

步骤C：根据以下公式确定R_i,j的中心位置_Ci,j和R_i+1,k的中心位置：

步骤D：根据以下公式确定R_i,j的中心位置_Ci,j和R_i+1,k的中心位置之间的距离

步骤E：根据以下公式计算突触R_i,j与第i+1层图像上所有突触之间距离的最小值▽X_i,j：

其中，n_i表示所有突触的数量。

步骤F：若▽X_i,j<θ，则R_i,j和R_i+1,k表示分布在相邻序列图像中的同一突触。其中，θ表示电镜序列切片图像中上下两层图像之间的最大偏移量。

步骤G：令i＝i+1，j＝k，并重复步骤C至步骤F，直至▽X_i,j≥θ。

步骤H：若R_i,j在电镜序列切片图像中出现的次数β≥3，则R_i,j为真正的突触；否则，R_i,j为伪突触。

S150：基于突触过滤的结果，进行突触的三维重构。

本步骤进一步可以包括：

S151：提取过滤得到的突触的光滑致密带结构。

本步骤可以进一步包括：

S1511：利用形态学运算，提取过滤得到的突触的致密带。

S1512：使用拟合的方法，对突触致密带的形态进行优化，得到光滑的突触致密带结构。

S152：对致密带进行重构，得到突触的三维结构。

具体地，本步骤可以包括：

S1521：计算第i层第j个突触致密带到第i+1层的第k个突触致密带的相对距离.

S1522：判断上述相对距离是否小于电镜序列切片图像中上下两层图像之间的最大偏移量；若是，则利用FiJi工具将第j个突触致密带和第k个突触致密带进行重构，得到突触致密带三维结构，也即突触的三维结构。

本发明实施例虽然以上述步骤顺序进行的描述，但是本领域技术人员应清楚这些步骤可以并行执行或者以颠倒的顺序执行，在此，本发明实施例对此不作限定，这些以并行或颠倒的顺序执行的技术方案也应包含在本发明的保护范围之内，在此不再赘述。

术语“包括”或者任何其它类似用语旨在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备/装置不仅包括那些要素，而且还包括没有明确列出的其它要素，或者还包括这些过程、方法、物品或者设备/装置所固有的要素。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。