CN105205840A - 一种电子断层图像中胶体金的直径估计及自动识别方法 - Google Patents

Abstract

本发明提供一种电子断层图像中胶体金的直径估计方法，包括：1)在原始胶体金直径附近的数值范围内得到一系列估计直径，对于每个估计直径执行步骤2)；2)对于每个估计直径，计算该估计直径下的人工模板与平均模板的互相关系数作为模板匹配互相关系数，其中，所述人工模板是根据估计直径直接生成的胶体金模板，所述平均模板是基于当前估计直径提取胶体金后，对所提取的胶体金区域做图像平均得到的胶体金模板；3)根据各个估计直径所对应的所述模板匹配互相关系数确定胶体金的精确直径。本发明还提供了相应的胶体金的自动识别方法。本发明对胶体金直径和数目等参数的依赖性较低，具有更高的识别率。

Description

一种电子断层图像中胶体金的直径估计及自动识别方法

技术领域

本发明涉及结构生物学电子断层成像技术领域，具体地说，本发明涉及一种电子断层图像中胶体金的直径估计及自动识别方法。

背景技术

电子显微镜(简称电镜)三维重构技术，利用电子显微镜拍摄的生物大分子投影图像序列进行三维重构得到生物大分子的三维密度结构，是结构生物学研究中的一种主流的技术支持手段。

根据投影的不同特点和分子样品的适用范围，电镜三维重构技术可分为：电子晶体学、单颗粒分析和电子断层成像技术。其中电子断层成像技术能够重构出不具有全同性且不需要结晶的细胞及亚细胞超微结构，具有其它技术不可替代的优势。受一系列机器及人工误差的影响，投影图像序列可能出现一系列的偏移旋转，且样品在空间中也会发生单纯的投影图像序列操作无法修复的运动，这极大地影响了重构结果精度，为了得到更高分辨率的三维重构结果，需要在进行三维重构之前对投影图像序列进行对位及空间参数标定。

基于胶体金的对位方法是目前最常用且对位精度最高的方法，该方法通过在样品中植入胶体金作为标记点对每幅样品投影图像进行对位。然而，为了获得高质量的对位结果，基于胶体金的对位一般需要较多的人工干预，这使得对位的效率较低，且由于对人工的依赖大其对位结果的质量并不稳定。

目前有很多关于胶体金自动对位的研究，一般存在两个难题：第一个是胶体金的提取，第二个是胶体金的追踪。本发明主要关注的是胶体金的提取。现有技术中，通常用模板匹配的方法提取胶体金的位置，模板匹配的方法可以快速有效的确定胶体金位置。对于不同的数据集或者同一个数据集的不同图片，需要构建相应的模板。而现有的胶体金自动提取算法(例如：RAPTOR算法，可参考文献RAPTOR：Markovrandomfieldbasedautomaticimagealignmentforelectrontomography；Automaticacquisitionoffiducialmarkersandalignmentofimagesintiltseriesforelectrontomography)十分依赖于人工输入的胶体金直径和数目等参数，如果人为估计的胶体金直径和数目偏差较大，会导致胶体金的提取结果存在较大偏差，甚至发生错误。

发明内容

本发明的任务是克服现有技术的不足，提供一种电子断层图像中胶体金的直径估计和相应的胶体金的自动识别解决方案。

与现有技术相比，本发明具有下列技术效果：

1、本发明对胶体金直径和数目等参数的依赖性较低。

2、本发明具有更高的识别率。

3、本发明能对所识别的胶体金进行一定筛选，排除质量较低的识别结果，从而方便后续的胶体金跟踪及电子断层图像序列的对齐。

附图说明

以下，结合附图来详细说明本发明的实施例，其中：

图1示出了根据本发明一个实施例的电子断层图像中胶体金的自动识别与定位方法的流程图；

图2示出了本发明一个实施例中的根据所给出的胶体金直径进行胶体金提取和定位的流程示意图；

图3示出了测试数据A在0度角下的投影的电镜图像；

图4示出了测试数据B在0度角下的投影的电镜图像；

图5示出了本发明一个实施例的胶体金提取和定位方法对测试数据A在0角度下的投影的电镜图像上所精确定位的胶体金点；

图6示出了本发明一个实施例的胶体金提取和定位方法对测试数据A在高角度下(65度)的投影的电镜图像上所精确定位的胶体金点；

图7示出了本发明一个实施例的胶体金提取和定位方法在测试数据B在0角度下的投影的电镜图像上所精确定位的胶体金点；

图8示出了本发明一个实施例的胶体金提取和定位方法在测试数据B在高角度下(75度)的投影的电镜图像上所精确定位的胶体金点；

图9示出了本发明一个实施例的胶体金提取和定位方法对测试数据A的所有图像做胶体金识别最后所精确定位的点的数量与最后没有定位到的点的数量的比较。

图10示出了本发明一个实施例的胶体金提取和定位方法对测试数据B的所有图像做胶体金识别最后所精确定位的点的数量与最后没有定位到的点的数量的比较。

图11示出了本发明一个实施例的胶体金提取和定位方法在测试数据B在高角度下(-80度)的投影的电镜图像上所精确定位的胶体金点。

具体实施方式

图1示出了根据本发明一个实施例的电子断层图像中胶体金的自动识别与定位方法的流程图，它包括下列步骤：

步骤1：用户输入胶体金直径和电子断层图像。电子断层图像通常由多个不同角度下拍摄的电镜图像组成。用户输入的胶体金直径通常为预估的胶体金直径。本文中也将这个直径称为原始胶体金直径。

步骤2：基于用户输入的胶体金直径，在一定范围内对胶体金的精确直径进行搜索，得到胶体金的精确直径。发明人在研究中发现，人工输入的胶体金直径往往与实际电镜图像中的胶体金直接差别较大，直接基于人工输入的胶体金直径创建胶体金的模板进而进行模板匹配，将会影响胶体金提取结果准确度。发明人针对这一问题进行了深入研究，发现胶体金平均化模板可以取得更好的效果，并进而提出了胶体金直径精确估计方案，从而有效地提高了胶体金提取和定位的准确度。

本步骤中，基于用于输入的原始胶体金直径进行模糊识别，得出胶体金的精确直径的方案包括下列子步骤：

步骤201：在原始胶体金直径附近的数值范围内建立一个搜索空间，以一定步长遍历这个搜索空间，其中遍历到的每个直径称为估计直径。

步骤202：对于每个估计直径，计算该估计直径下的人工模板与平均模板的互相关系数作为模板匹配互相关系数。

其中，人工模板是根据已知胶体金直径所创建的人工模板，平均模板是通过对已经提取的胶体金做图像加和再平均化所得的胶体金平均化模板。下面介绍两种不同模板的创建方法。

人工模板：依据胶体金在投影图像中呈现的固有特性(胶体金具有大小相当，形状类似圆形，与周围图像的对比度高等固有特性)生成胶体金的模板。创建模板的实现比较简单，即根据已知的直径，画一个黑色实心圆，背景为白色，构成一个模板图像。

平均模板：在有一个已经优化的胶体金直径参数所创建的人工模板的前提下提取一张典型的电镜图像中的胶体金，将所提取的胶体金区域的图像叠加求平均得到胶体金平均化模板，这个胶体金平均化模板就能够用于提取其它电镜图像的胶体金。

步骤203：然后根据各个估计直径所对应的所述模板匹配互相关系数，确定胶体金的精确直径。人工模板与平均模板的互相关系数越大，即模板匹配互相关系数越大，表明人工模板与平均模板的匹配程度越高，其对应的估计直径就越接近实际的胶体金直径。因此可以根据各个估计直径所对应的模板匹配互相关系数来寻找胶体金的精确直径。

步骤3：基于精确直径创建胶体金模板，得出胶体金提取和定位结果。

步骤4：输出所有胶体金的位置。

前述实施例的步骤2中，通过预先建立的搜索空间和步长来搜索精确直径，这种做法通常能够获得一个比原始胶体金直径更优的搜索结果，并且搜索方案比较简单。然而，由于搜索方案的局限性，例如只能依次搜索预先建立的搜索空间中的有限个点，因此它所得到的模板匹配互相关系数信息是有限的，在某些情况下，这种简单搜索方案所得的直径可能仍然与实际胶体金直径存在一定偏差。因此，发明人提出了另一个实施例，该实施例基于最大期望(EM)思想对输入的胶体金直径进行优化的策略，来实现胶体金直径的精确估计。

该实施例中，所述步骤201包括：根据用户输入的胶体金直径值初始化直径d。初始化搜索范围ξ和精确范围ξ₀。初始化直径的搜索空间：DIA＝{d_i|d-ξ/2<d_i<d+ξ/2；i＝1,2…m}。遍历搜索空间DIA，对于其中每个直径d_i执行步骤202。

所述步骤202包括下列子步骤：

步骤2021：基于搜索空间的某个直径d_i创建胶体金人工模板。

步骤2022：基于直径d_i，用已知确切直径下的胶体金提取法提取胶体金位置。

步骤2023：将所有提取出的胶体金区域做平均，生成一个平均模板。此处平均是指图像平均：首先提取出胶体金区域，然后将胶体金区域做图像叠加处理，再除以所提取的胶体金区域数目，就得到了所述平均模板，这个平均模板就是胶体金平均化模板，也可称为胶体金平均模板。

步骤2024：计算胶体金人工模板与胶体金平均模板的互相关值k_i。

步骤203包括下列子步骤：

步骤2031：找出互相关值k_i最大的和第二大的直径：d_max,d_sub。

步骤2032：计算：d＝(d_max+d_sub)/2；ξ＝(d_max-d_sub)。更新搜索空间DIA＝{d_i|d-ξ/2<d_i<d+ξ/2；i＝1,2…m}。

步骤2033：遍历更新后的搜索空间，进一步搜索精确直径，直到ξ≤ξ₀。本步骤中，遍历更新后的搜索空间进一步搜索精确直径的方法与第一轮搜索的方法相同，即重复执行前文所述的步骤202以及执行步骤2031和步骤2032即可。

该实施例中，最后得到的直径d即为优化后的胶体金直径。搜索范围ξ和精确范围ξ₀的选取具有一定的经验性。搜索范围ξ代表在用户输入的初始直径d的基础上，±ξ/2的范围内搜索最优半径，所以ξ值越大，搜索范围越广，对用户输入的容错度越高，但是搜索时间也越长。优选地，在一个例子中设置ξ的值为10。精确范围ξ₀代表直径的精确程度，ξ₀值越小，得到的直径越精确，但是也需要更多的计算时间。优选地，在一个例子中ξ₀的取值为0.5。

图2示出了本发明一个实施例中的根据所给出的胶体金直径进行胶体金提取和定位的流程示意图，该流程可以用于上述步骤3中，也可以用于其它任何一个需要在已知胶体金直径前提下提取胶体金的步骤或子步骤中。参考图2，该实施例中，根据所给出的胶体金直径创建胶体金模板后，胶体金提取和定位的流程包括下列步骤：

1.对输入的电镜投影图像(下文中简称电镜图像)做归一化。将胶体金模板与电镜图像做互相关计算并归一化，得到互相关系数图像。其中胶体金模板即可以是人工模板，也可以是平均模板。

2.对互相关系数图像进行网格式的搜索。电镜图像中分布良好的胶体金是不应该重叠的，因此，在已知准确的胶体金直径的前提下，一张图像中的胶体金的容量也是固定的。本实施例所进行的网格式的搜索如下：将互相关系数图像等分成多个正方形区域，每个区域长与宽都为dia(dia即胶体金直径)；在每个小正方形区域内，搜索得到其中的最大的互相关系数值及其所对应的电镜图像上的像素的坐标。对于一张长为w宽为h的电镜图像，理论上一共可以得到的胶体金位置候选点个数为w×h/dia²。对于电镜图像上的每个像素，它都拥有两个标签：该像素在电镜图像上的值(通常为亮度值)，该像素所对应的互相关系数的值。在此，设立两个阈值：thre1及thre2，其中thre1为互相关系数图像的像素平均值，thre2为电镜图像的像素平均值。当一个正方形区域内的一个极值点的坐标所对应的互相关系数的值大于thre1且其所对应的电镜图像上的像素值大于thre2，则认为所对应的电镜图像上的像素是一个胶体金候选点。

3.分析胶体金候选点集中两个标签互相关数值和像素值的分布情况，计算数值集合{所有候选点的互相关数值×像素值}的均值μ和标准差σ。利用候选点的互相关数值与电镜图像像素值的乘积做为一标的数，候选点集进行筛选。在此设立thre3，在一个例子中，通过大量实验确定阈值设为thre3＝μ+2.5σ。对于所对应的互相关数值与像素值的乘积大于thre3的候选点，将其加入集合FID，以方便后续分析。

4.基于归一化后的电镜图像，对FID中的点进行区域生长。区域生长所指的是搜索原胶体金候选点相临的像素点，若其相临的像素点的值大于一定的阈值，则认定该相临的像素点与的原胶体金候选点是相似的。若一个像素点被判断为相似，则进一步搜索该像素点的相临点，若被判断为并不相似，则停止在该点的相临点的搜索，最终这些相似像素将组成为一个区域，称为胶体金区域。胶体金区域的大小与两个因素有关，其一为电镜图像本身的像素值的变化，其二为所设置的阈值thre3。本步骤的目标并不是根据区域生成来进一步精确定位胶体金，而是通过区域生长进一步排除一点分布及背影并不理想的胶体金。因此，所设置的阈值为以侯先点为中心，周围胶体金半径区域内的所有像素的均值。最终，如果区域生长所得的区域的直径大于1.4×dia(胶体金直径)，即认为当前的胶体金候选点对应的区域并非真实的胶体金，将该候选点排除，如果区域生长所得的区域的直径不大于1.4×dia(胶体金直径)，则保留该候选点。

5.用重心法确定FID中的胶体金的精确位置。即对于每个保留的胶体金候选点所对应的胶体金区域，用重心法确定该处胶体金的精确位置。

上述步骤1至5中，互相关计算(NCC)是用来确定投影图像某个区域与胶体金模板在形状上的相似性。对投影图像和互相关矩阵的归一化(归一化并末截断数值，即末改变数值的相对分布)是为了方便进一步的分析。然后图像被划分为w×h/dia²个区域，对于每一个区域，互相关值最大的点被选为胶体金位置的候选点。在实现时，本实施例用了非常宽松的阈值(投影图像与互相关结果图的像素平均值)来确定候选点集，这样在最大程度上避免了阈值选取对最终结果的干扰，可最大限度地选中所有胶体金。候选点集中包含了几乎所有的胶体金位置，但也包含大量的非胶体金点，所以必须对候选点集做进一步的筛选。通过大量实验，发现候选点集中坐标点所对应的互相关值×像素值的直方图分布满足高斯分布。因此把聚类分析运用到这个问题中。在实现时，把坐标点所对应的互相关值与像素值的乘积的分布拟合到一个高斯模型中，计算其均值μ和标准差σ。通过大量实验，发现阈值取μ+2.5σ可以较准确的筛选出胶体金位置。最后用区域生长的方法获取胶体金的区域，然后根据形状再对胶体金做一轮筛选，然后用重心法精确胶体金位置，这样可以实现亚像素精确度。

发明人使用了两组数据进行检测，以验证本发明的实际效果。第一组数据为大小1024×1024的线粒体数据(测试数据A)；第二组数据为大小2048×2048的聚合分子高倍数显微数据(测试数据B)。图3和图4分别示出了测试数据A和B在0度角下的投影的电镜图像。图5和图6分别示出了本发明一个实施例的胶体金提取和定位方法对测试数据A在低角度下(0度)下和高角度下(65度)的投影的电镜图像上所精确定位的胶体金点，其中黑色圆圈中的胶体金为本发明一个实施例的方法定位得到的胶体金，黑色圆圈的中心即为所定位的胶体金的中心。图7和图8分别示出了本发明一个实施例的胶体金提取和定位方法对测试数据B在低角度下(0度)下和高角度下(75度)的投影的电镜图像上所精确定位的胶体金点。通过肉眼观察，可以发现本发明的方法在这几张示例图上的识别率是非常高的，都在所有胶体金的数量的95％以上。图9示出了本发明一个实施例的胶体金提取和定位方法对测试数据A的所有图像做胶体金识别最后所精确定位的点的数量与最后没有定位到的点的数量的比较；图10示出了本发明一个实施例的胶体金提取和定位方法对测试数据B的所有图像做胶体金识别最后所精确定位的点的数量与最后没有定位到的点的数量的比较。可以发现，除了仅有的前后几张的高角度上的图像外，本发明的方法的胶体金识别率都在95％以上，有很多图像上的识别率近乎100％。另外，需要说明的是，根据后继的算法上的分析可以看出，本发明在高角度上并不是无法识别出足够多的胶体金，而是由于形状不好或是背景噪声太强，所识别出的胶体金并不适合用做后继的对位，因此所识别出的胶体金被舍弃了一部分。在此，特别选取测试数据B第1张投影图像作为例子，图11示出了本发明一个实施例的胶体金提取和定位方法在测试数据B在高角度下(-80度)的投影的电镜图像上所精确定位的胶体金点，其中，测试数据B在高角度下(-80度)的投影的电镜图像为测试数据B第1张投影图像。从图11的黑色区域的胶体金可以看出，这些胶体金或是形状不好或是背景噪声太强，是不适合用作电镜图像的胶体金对位的，因此在区域生长阶段抛弃了这些不适合的点，以提高后续对位处理的准确度。

最后应说明的是，以上实施例仅用以描述本发明的技术方案而不是对本技术方法进行限制，本发明在应用上可以延伸为其它的修改、变化、应用和实施例，并且认为所有这样的修改、变化、应用、实施例都在本发明的精神和教导范围内。

Claims (13)

1.一种电子断层图像中胶体金的直径估计方法，包括下列步骤：

1)在原始胶体金直径附近的数值范围内得到一系列估计直径，对于每个估计直径执行步骤2)；

2)对于每个估计直径，计算该估计直径下的人工模板与平均模板的互相关系数作为模板匹配互相关系数，其中，所述人工模板是根据估计直径直接生成的胶体金模板，所述平均模板是基于当前估计直径提取胶体金后，对所提取的胶体金区域做图像平均得到的胶体金模板；

3)根据各个估计直径所对应的所述模板匹配互相关系数确定胶体金的精确直径。

2.根据权利要求1所述的电子断层图像中胶体金的自动识别与定位方法，其特征在于，所述步骤1)包括下列子步骤：

11)根据原始胶体金直径值，初始化直径d；

12)初始化搜索范围ξ；

13)初始化直径的搜索空间：DIA＝{d_i|d-ξ/2＜d_i＜d+ξ/2；i＝1,2…m}。

3.根据权利要求2所述的电子断层图像中胶体金的自动识别与定位方法，其特征在于，所述步骤2)还包括：对于搜索空间中的每个估计直径d_i，基于估计直径d_i直接创建人工模板，基于该人工模板从电镜图像识别出胶体金区域，将所识别出的所有胶体金区域做图像平均，生成该估计直径d_i对应的平均模板，计算每个估计直径d_i对应的人工模板与胶体金平均模板的互相关值k_i。

4.根据权利要求3所述的电子断层图像中胶体金的自动识别与定位方法，所述步骤3)包括下列子步骤：

31)找出使互相关值k_i最大的和第二大的估计直径：d_max,d_sub；

32)计算d＝(d_max+d_sub)/2；ξ＝(d_max-d_sub)；

33)更新直径搜索的搜索空间：DIA＝{d_i|d-ξ/2＜d_i＜d+ξ/2；i＝1,2…m}，遍历更新后的搜索空间重新执行步骤2)，再重新执行步骤31)、32)，直至ξ≤ξ₀，ξ₀为预设的精确范围阈值。

5.根据权利要求4所述的电子断层图像中胶体金的自动识别与定位方法，所述步骤33)还包括：当ξ≤ξ₀时，将此时的d作为精确直径。

6.一种电子断层图像中胶体金的自动识别方法，包括下列步骤：

a)根据权利要求1～5中任意一项所述的电子断层图像中胶体金的直径估计方法得出胶体金的精确直径；

b)基于所述精确直径对应的胶体金模板，得出胶体金识别结果。

7.根据权利要求6所述的电子断层图像中胶体金的自动识别方法，其特征在于，所述步骤b)包括下列子步骤：

b1)对输入的电镜图像做归一化，将胶体金模板与电镜图像做互相关计算并归一化，得到互相关系数图像；

b2)将互相关系数图像等分成多个正方形区域，每个区域长与宽都为当前估计直径dia；在每个正方形区域内，搜索得到其中的最大的互相关系数值及其所对应的电镜图像上的像素的坐标，根据该像素的坐标所对应的电镜图像上的值以及所对应的互相关系数的值的大小，判别该像素是否为胶体金候选点；

b3)根据互相关数值与电镜图像像素值的乘积，进一步筛选胶体金候选点。

8.根据权利要求7所述的电子断层图像中胶体金的自动识别方法，其特征在于，所述步骤b2)中，当一个正方形区域内的一个极值点的坐标所对应的互相关系数的值大于thre1且其所对应的电镜图像上的像素值大于thre2，则认为所对应的电镜图像上的像素是一个胶体金候选点；其中，thre1为互相关系数图像的像素平均值，thre2为电镜图像的像素平均值。

9.根据权利要求8所述的电子断层图像中胶体金的自动识别方法，其特征在于，所述步骤b3)还包括：对于每个胶体金候选点，如果其互相关数值与电镜图像像素值的乘积未达到预设的阈值thre3时，则将该胶体金候选点筛除。

10.根据权利要求9所述的电子断层图像中胶体金的自动识别方法，其特征在于，所述步骤b3)还进一步包括：所述预设的阈值thre3根据步骤b2)所得的所有胶体金候选点的互相关数值与电镜图像像素值的乘积的均值μ和标准差σ确定。

11.根据权利要求6所述的电子断层图像中胶体金的自动识别方法，其特征在于，所述步骤b)包括进一步包括下述子步骤：

b4)基于归一化后的电镜图像，对步骤b3)筛选后保留的每个胶体金候选点进行区域生长得到对应的胶体金区域，排除生长后的区域直径过大的胶体金候选点。

12.根据权利要求11所述的电子断层图像中胶体金的自动识别方法，其特征在于，所述步骤b)包括进一步包括下述子步骤：

b5)基于步骤b4)所得到的各个胶体金候选点所对应的生长后区域，根据重心法确定各个胶体金的精确位置。

13.根据权利要求11所述的电子断层图像中胶体金的自动识别方法，其特征在于，所述步骤b4)中，对每个胶体金候选点进行区域生长的方法为：搜索原胶体金候选点相临的像素点，若其相临的像素点的值大于一定的阈值，则认定该相临的像素点与的原胶体金候选点是相似的，并进一步搜索该像素点的相临点，否则，停止在该点的相临点的搜索，最终这些相似像素将组成为一个与原胶体金候选点对应的生长后的胶体金区域。