基于双光子图像的突触检测方法与流程

本发明涉及模式识别和神经学技术领域，尤其是涉及一种基于双光子图像的突触检测方法。

背景技术：

大量研究表明，神经元的突触可塑性包括功能可塑性和结构可塑性，与学习和记忆密切相关。突触可塑性调节异常常常被认为是认知、学习和记忆功能障碍的主要原因。许多精神和认知障碍发病机制与突触可塑性调节的异常密切相关,对突触可塑性调节的研究有利于揭示智力发育迟滞、认知功能障碍这类神经精神疾病的病理机制。

近年来，利用正置双光子显微镜系统和荧光探针标记技术方法一直都是研究的热点。现有检测突触方法，就是神经学家通过肉眼一个一个的寻找突触。该方法耗费大量人工时间，而且精度不高，易出现遗漏的情况。

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

以下为与本发明相关的现有技术：

[1]N.Dey,L.Blanc-Feraud,C.Zimmer,“Richardson–Lucy algorithm with total variation regularization for 3D confocal microscope deconvolution,”Microscopy research and technique,vol.69,no.4,pp.260-266,2006.

[2]J.E.Ledoux,“Emotion circuits in the brain,”Annual Review of Neuroscience,vol.23,no.2,pp.155-184,2000.

[3]V.D.Paola,A.Holtmaat,G.Knott,et al.,“Cell Type-Specific Structural Plasticity of Axonal Branches and Boutons in the Adult Neocortex,”Neuron,vol.49,no.6,pp.861-875,2006.

[4]N.Becker,C.J.Wierenga,R.Fonseca,et al.,“LTD induction causes morphological changes of presynaptic boutons and reduces their contacts with spines,”Neuron,vol.60,no.4,pp.590-597,2008.

[5]F.Jing,X.Zhou,J.G.Dy,et al.,“An Automated Pipeline for Dendrite Spine Detection and Tracking of 3D Optical Microscopy Neuron Images of In Vivo Mouse Models,”Neuroinformatics,vol.7,no.2,pp.113-130,2009.

[6]H.Park and M.M.Poo,“Neurotrophin regulation of neural circuit development and function,”Nature Reviews Neuroscience,vol.14,no.1,pp.7-23,2013.

[7]F.W.Grillo,S.Song,L.M.Teles-Grilo,et al.,“Increased axonal bouton dynamics in the aging mouse cortex,”Proceedings of the National Academy of Sciences of the United States of America,vol.110，no.16,pp.1514-1523,2013.

[8]Y.Yang,D.Q.Liu,W.Huang,et al.,“Selective synaptic remodeling of amygdalocortical connections associated with fear memory,”Neuroscience,to be published.

技术实现要素：

本发明实施例提供一种基于双光子图像的突触检测方法，以解决如何提高检测突触的效率和精度的问题。

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

一种基于双光子图像的突触检测方法，所述方法包括：

获取所述双光子图像；

将所述双光子图像分为轴突图像和树突图像；

在所述轴突图像中进行轴突小结的2D检测；

在所述树突图像中进行树突棘的2D检测；

用检测到的轴突小结和树突棘进行2D突触检测；

基于2D突触检测结果，进行3D突触检测。

进一步地，所述在所述轴突图像中进行轴突小结的2D检测，具体包括：

对所述轴突图像进行迭代的3D反卷积处理，得到反卷积图像；

确定所述反卷积图像中的局部极大值点；

基于所述轴突图像，确定增强图像；

对所述增强图像进行二值分割处理；

基于所述局部极大值点以及二值分割处理结果，去除所述二值分割处理结果中不包括所述局部极大值点的区域。

进一步地，所述确定所述反卷积图像中的局部极大值点具体包括：

设定第一阈值；

找出所述反卷积图像中像素值大于第一阈值的疑似局部极大值的点；

判断所述各疑似局部极大值的点是否满足在一预定邻域内的所述像素值是最大的，如果是，则确定所述疑似局部极大值的点为所述局部极大值点。

进一步地，所述基于所述轴突图像，确定增强图像具体包括：

以二维xy坐标系为基准，基于所述轴突图像，根据以下公式确定所述增强图像：

其中，所述g表示所述增强图像；所述fx表示所述轴突图像在x方向上的偏导；所述fy表示所述轴突图像在y方向上的偏导；所述fxx表示所述fx在所述x方向上的偏导；所述fxy表示所述fx在所述y方向上的偏导；所述fyy表示所述fy在所述y方向上的偏导。

进一步地，所述方法包括：所述在所述树突图像中进行树突棘的2D检测具体包括：

对所述树突图像进行归一化和结构增强处理；

骨架化结构增强处理结果中的树突，得到树突骨架二值图；

基于所述树突骨架二值图，确定树突骨架上的分叉点；

基于所述树突骨架二值图，并根据所述分叉点，定位所述树突棘。

进一步地，所述对所述树突图像进行归一化和结构增强处理具体包括：

将所述树突图像归一化到0和1之间；

根据以下公式，设定树突棘界面模型：

其中，所述θ表示所述树突棘的角度；所述σ表示所述树突棘的宽度；所述F(x,y)表示所述树突棘界面模型；所述x、所述y分别表示归一化后的图像中像素点的横、纵坐标；

按照以下模型增强归一化后图像的线结构：

进一步地，所述骨架化结构增强处理结果中的树突，得到树突骨架二值图具体包括：

对结构增强处理结果进行阈值分割，得到树突二值图像；

遍历所述树突二值图像，根据以下条件删除像素点，直至所述树突二值图像不再变化为止，得到第一二值图像：

XH(p)＝1；2≤min{n1(p),n2(p)}≤3；

其中，所述p表示像素点；所述x1,x2,...,x8为从右侧邻点开始按逆时针排序的所述p的8个相邻点像素点；

遍历所述第一二值图像，根据以下条件删除像素点，直至所述树突二值图像不再变化为止，得到所述树突骨架二值图：

XH(q)＝1；2≤min{n1(q),n2(q)}≤3；

其中，所述q表示像素点；所述x1,x2,...,x8为从右侧邻点开始按逆时针排序的所述q的8个相邻点像素点。

进一步地，所述基于所述树突骨架二值图，确定树突骨架上的分叉点具体包括：

构造滤波器；

将所述滤波器与所述树突骨架二值图进行卷积处理，并将卷积图中大于第二阈值的点确定为分叉点。

进一步地，所述基于所述树突骨架二值图，并根据所述分叉点，定位所述树突棘，具体包括：

构造预定半径的结构元素，并对所述树突骨架二值图进行形态学腐蚀；

用所述树突骨架二值图减去腐蚀后的骨架图，得到树突棘备选位置；

基于所述树突骨架二值图，以所述分叉点为中心，提取预定边长的区域；

选取所述分叉点所在区域的骨架部分，并判断是否存在所述树突棘备选位置与该骨架部分相连；如果有，则确定所述树突棘备选位置为树突棘。

进一步地，所述用检测到的轴突小结和树突棘进行2D突触检测具体包括：

用各所述轴突小结的中心点来代表所述轴突小结；

用所述树突棘备选位置的中心点来代表所述树突棘；

计算并记录所述各轴突小结中心点和所有所述树突棘中心点之间距离的最小值；

若所述距离的最小值小于第三阈值，则将所述轴突小结和所述树突棘确定为2D突触。

进一步地，所述基于2D突触检测结果，进行3D突触检测具体包括：

用所述2D突触中轴突小结中心点和树突棘备选位置中心点的中点来代表所述突触的位置；

用所述2D突触中所述轴突小结中心点和所述树突棘备选位置中心点的距离来代表所述突触的真实性；

在所述双光子图像序列切片相邻的前后层中，在所述2D突触位置的邻域内的坐标内查找所述突触；

判断各所述突触的真实性在所述双光子图像序列切片相邻的前后层中是否是最高的，若是，则确定所述突触为3D突触。

本发明实施例提供一种基于双光子图像的突触检测方法。通过将双光子图像分为轴突图像和树突图像，分别检测轴突小结和树突棘来找到疑似突触以供神经学家筛选，用检测到的轴突小结和树突棘进行2D突触检测；然后基于2D突触检测结果，进行3D突触检测，这样大大地提高了检测突触的效率和精度，还具有普适性，且对图像类型不敏感。

附图说明

图1为根据本发明实施例的基于双光子图像的突触检测方法流程示意图；

图2为根据本发明实施例的基于树突二值图像，删除像素点的示意图；

图3为根据本发明实施例的滤波器示意图；

图4为根据本发明实施例的双光子图像序列切片层间关系示意图。

具体实施方式

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

本发明实施例提供一种基于双光子图像的突触检测方法。如图1所示，该方法可以包括：

S100：获取双光子图像。

例如：该图像可以是对鼠脑杏仁核远距离投射过来轴突和听觉皮层的树突进行荧光标记的图像。

S110：将双光子图像分为轴突图像和树突图像。

其中，在双光子图像中轴突和树突是用不同的颜色显示的。通过抽取对应颜色的部分即可从中提取轴突或者树突图像。

S120：在轴突图像中进行轴突小结的2D检测。

本步骤进一步可以包括：

S121：对轴突图像进行迭代的3D反卷积处理，得到反卷积图像。

在实际应用中，本步骤可以通过ImageJ工具中的插件Diffraction PSF 3D和插件Iterative Deconvolution 3D.class来实现。其中，运行插件Diffraction PSF 3D，输入相关参数，生成一个卷积核。优选地，可以按照如下配置来输入参数和选取范围：介质折射率：常用范围(1.3-1.5)；物镜孔径：常用范围(1.0-1.35)；波长：常用范围(460-500nm),最佳地，为488；平面图像分辨率：约200nm；z轴分辨率：约700nm；图像的宽、高、层数依次为512，512，140；纵向球面像差最大孔径：常用范围(120-870)，最佳地，为500-700nm。同时选择原图和生成的卷积核，并输入反卷积次数(优选地为10次)。然后运行插件Iterative Deconvolution 3D.class计算得出反卷积的图像。

S122：确定反卷积图像中的局部极大值点。

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

S1221：预先设定一个第一阈值(例如：1000)。

S1222：找出反卷积图像中像素值大于第一阈值的疑似局部极大值的点。

S1223：判断各疑似局部极大值的点是否满足在一预定邻域内的像素值是最大的，如果是，则执行步骤S1224；否则执行步骤S1225。

S1224：确定该疑似局部极大值的点为局部极大值点。

S1225：确定该疑似局部极大值的点不是局部极大值点。

S123：基于轴突图像，确定增强图像。

其中，本步骤可以包括：

以二维xy坐标系为基准，基于轴突图像，根据以下公式确定增强图像：

其中，g表示增强图像；fx表示轴突图像在x方向上的偏导；fy表示轴突图像在y方向上的偏导；fxx表示fx在x方向上的偏导；fxy表示fx在y方向上的偏导；fyy表示fy在y方向上的偏导。

S124：对增强图像进行二值分割处理。

S125：基于局部极大值点以及二值分割处理结果，去除二值分割处理结果中不包括局部极大值点的区域。

S130：在树突图像中进行树突棘的2D检测。

本步骤可以包括：

S131：对树突图像进行归一化和结构增强处理。

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

S1311：将树突图像归一化到[0,1]。

S1312：根据以下公式，设定树突棘界面模型：

其中，θ表示树突棘的角度；σ表示树突棘的宽度；F(x,y)表示树突棘界面模型；x、y分别表示归一化后的图像中像素点的横、纵坐标。

S1313：按照以下模型增强归一化后图像的线结构：

其中，本步骤中涉及到的参数的含义参见步骤S1312的有关说明，在此不再赘述。

S132：骨架化结构增强处理结果中的树突，得到树突骨架二值图。

本步骤可以包括：

S1321：对结构增强处理结果进行阈值分割，得到树突二值图像。

S1322：遍历树突二值图像，根据以下条件删除像素点，直至树突二值图像不再变化为止，得到第一二值图像：

XH(p)＝1；2≤min{n1(p),n2(p)}≤3；

其中，

其中，p表示像素点；x1,x2,...,x8为从右侧邻点开始按逆时针排序的p的8个相邻点像素点。

本步骤中，XH(p)是一个判别是否删除像素点p的条件。bi决定XH(p)的取值。例如：如果满足x1＝0且(x2＝1或x3＝1)，则b1＝1。进而得出XH(p)的值。

x1,x2,...,x8为从右侧邻点开始按逆时针排序的p的8个相邻点像素点，如图2所示。

S1323：遍历第一二值图像，根据以下条件删除像素点，直至树突二值图像不再变化为止，得到树突骨架二值图：

XH(q)＝1；2≤min{n1(q),n2(q)}≤3；

其中，

其中，q表示像素点；x1,x2,...,x8为从右侧邻点开始按逆时针排序的q的8个相邻点像素点。

有关本步骤中涉及到的参数的说明可以参见步骤S1322的说明，在此不再赘述。

S133：基于树突骨架二值图，确定树突骨架上的分叉点。

本步骤具体可以包括：

S1331：构造滤波器。

例如：构造尺寸为3×3的滤波器，如图3所示。其中，4个顶点取值为0，其余点取值为1。

S1332：将滤波器与树突骨架二值图进行卷积处理，并将卷积图中大于第二阈值的点确定为分叉点。

例如，可以将卷积图中值大于等于4的点确定为分叉点。

S134：基于树突骨架二值图，并根据分叉点，定位树突棘。

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

S1341：构造预定半径的结构元素，并对树突骨架二值图进行形态学腐蚀。

本步骤中，预定半径可以根据树突棘的长度(其占几个像素，例如3)来确定。

S1342：用树突骨架二值图减去腐蚀后的骨架图，得到树突棘备选位置。

S1343：基于树突骨架二值图，以分叉点为中心，提取预定边长的区域。

其中，预定边长可以根据树突棘的长度(其占几个像素，例,31)来确定，优选地，其可以大于树突棘长度的2倍。

S1344：选取分叉点所在区域的骨架部分，并判断是否存在树突棘备选位置与该骨架部分相连；如果有，则执行步骤S1345；否则，执行步骤S1346。

S1345：确定树突棘备选位置为树突棘。

S1346：确定树突棘备选位置不是树突棘。

S140：用检测到的轴突小结和树突棘进行2D突触检测。

本步骤进一步可以包括：

S141：用每一个轴突小结的中心点来代表该轴突小结。

S142：用树突棘备选位置的中心点来代表树突棘。

S143：计算并记录每一个轴突小结中心点和所有树突棘中心点之间距离的最小值d。

S144：若该距离最小值d小于第三阈值，则将轴突小结和树突棘确定为2D突触。

若该距离最小值d大于等于第三阈值，则不将轴突小结和树突棘确定为2D突触。

S150：基于2D突触检测结果，进行3D突触检测。

本步骤可以包括：

S151：用2D突触中轴突小结中心点和树突棘备选位置中心点的中点来代表突触的位置。

S152：用上述2D突触中轴突小结中心点和树突棘备选位置中心点的距离来代表突触的真实性。

轴突小结中心点和所有树突棘中心点之间距离的最小值越小，真实性越高。

S153：在双光子图像序列切片相邻的前后层中，在2D突触位置的邻域内的坐标内查找突触。

S154：判断每一个突触的真实性在双光子图像序列切片相邻的前后层中是否是最高的，若是，则执行步骤S155；否则，执行步骤S156。

S155：确定该突触为3D突触。

S156：确定该突触不是3D突触。

图4示例性地示出了双光子图像序列切片层间关系。其中，r表示双光子图像序列切片；r-1、r+1分别表示双光子图像序列切片相邻的前后层。

在双光子图像序列切片上进行2D分割，然后将2D分割结果在3D上串起来，找到3D突触。

上述实施例中的实施方案可以进一步组合或者替换，且实施例仅仅是对本发明的优选实施例进行描述，并非对本发明的构思和范围进行限定，在不脱离本发明设计思想的前提下，本领域中专业技术人员对本发明的技术方案作出的各种变化和改进，均属于本发明的保护范围。