CN103824294A - 一种电子断层图像序列对位方法 - Google Patents

Abstract

本发明提供一种电子断层图像序列对位方法，包括下列步骤：1）采用能够保留图像仿射变换信息的特征提取算法，从每幅电子断层投影图像中提取特征点；2）对不同电子断层投影图像的特征点进行匹配，获得多幅电子断层投影图像的特征点之间的对应关系；3）基于仿射模型，根据步骤2）所获得的多幅电子断层投影图像的特征点之间的对应关系，估计各电子断层投影图像所对应的对位参数。本发明不依赖胶体金标记，且能够提高对位准确度。

Description

一种电子断层图像序列对位方法

技术领域

本发明涉及电子断层成像技术领域，具体地说，本发明涉及一种电子断层图像序列对位方法。

背景技术

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

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

传统电子断层图像对位技术中，通过在样品中植入胶体金作为标记点对每幅样品投影图像进行对位。然而，胶体金的植入增加了噪声，对最终重构结果的分辨率有不利影响。另外，某些样品无法植入胶体金，因此传统电子断层图像对位技术的适用范围受限。

近年来，人们提出了一种无胶体金标记对位技术，目前主要可分为二大类。第一类是通过互相关方法求出图像之间偏移旋转参数后进行对位;第二类是通过虚拟特征点方法求出图像空间参数后对位。第一类方法大部分只能在二维图像平面内纠正图像平移和旋转，而基于样品三维重构-再投影的互相关改进方法则能求出少量空间参数，但非常耗时。第二类方法计算量适中，且通过引入计算机视觉的方法，理论上可估计样品的三维空间参数，恢复更多的对位参数，但实践上该方法估计的对位参数的准确度往往难以达到应用要求。

发明内容

因此，本发明的任务是克服现有技术的不足，提供一种准确度高的不依赖胶体金标记的电子断层图像对位解决方案。

为实现上述发明目的，本发明提供了一种电子断层图像序列对位方法，包括下列步骤：

1）采用能够保留图像仿射变换信息的特征提取算法，从每幅电子断层投影图像中提取特征点；

2）对不同电子断层投影图像的特征点进行匹配，获得多幅电子断层投影图像的特征点之间的对应关系；

3）基于仿射模型，根据步骤2）所获得的多幅电子断层投影图像的特征点之间的对应关系，估计各电子断层投影图像所对应的对位参数。

其中，所述步骤1）中，采用SIFT算子或者SURF算法从每幅电子断层投影图像中提取特征点。

其中，所述步骤2）包括下列子步骤：

21）采用局部搜索法，对每幅电子断层投影图像，将其与其它电子断层投影图像进行两幅电子断层投影图像之间的特征点匹配；

22）将步骤21）中所匹配的特征点匹配对串联，得到对应于被投影样品相同特征区块的各个不同角度的电子断层投影图像的特征点组成的特征链。

其中，所述步骤21）中，采用随机抽样一致原则，基于二极几何去除错误匹配对。

其中，所述步骤21）中，每幅所述电子断层投影图像仅与邻近的电子断层投影图像进行特征点匹配。

其中，所述步骤3）中，各电子断层投影图像所对应的对位参数包括：每幅电子断层投影图像所对应的相机空间旋转参数、每幅电子断层投影图像所对应的样品在空间中的偏移参数以及各特征点所对应的样品原始特征点空间坐标。

其中，所述步骤3）包括下列子步骤：

31）根据穿过图像的特征链数目选取两张图像做为初始待估计序列；

32）设置初始待估计序列的对位参数的初值，所述对位参数包括图像所对应的相机空间旋转矩阵R、样品在空间中的偏移t、以及对应于图像中匹配特征点的样品空间点(X,Y,Z,1)^T，其中，R的初值通过初始旋转角度确定，t的初值通过投影图像对应的相机旋转角度和相机深度共同确定，根据待估计序列中各特征点的值，基于R、t的初值进行反投影，获得所对应的每一个空间点(X,Y,Z,1)^T的初值；

对当前待估计序列中的每幅图像，寻找使得效用函数E最小的参数R、t，

$E=\underset{i,j}{\mathrm{\Σ}}{({\mathrm{Proj}}_i\left({\tilde{X}}_j\right)-x_{i,j})}^2{\mathrm{\δ}}_{i,j}$

是第j个空间点的估计值，

是对第j个空间点的估计值进行第i像面的重投影的输出值，x_i,j为第j个空间点在第i像面的投影的测量值，δ_i,j表明点j在投影i中是否可见，如可见则δ_i,j取1，如不可见则δ_i,j取0；

33）将一幅新电子断层图像加入，得到新的待估计序列；

34）以当前已知或已估计出的电子断层投影图像所对应的对位参数作为初值，通过位于低角度的投影点的反投影获得种子序列中新出现的空间点的初值；

35）将当前种子序列中的每幅图像，寻找使得所述效用函数E最小的参数R、t，并据此更新当前待估计序列中每幅图像所对应的参数R、t；

36）重复执行步骤33）至35）直至所有的图像均已加入待估计序列。

其中，所述步骤32）、35）中，采用稀疏捆绑调整技术来求解使得所述效用函数E最小的参数R、t。

其中，所述步骤35）包括下列子步骤：

351）输入当前待估计序列中所有测量得到的特征点x_i,j及初始化的对位参数变量；

352）求解使得所述效用函数E最小的每幅图像的对位参数并用所求得的对位参数替换对位参数变量的初值；

353）根据所得的对位参数进行重投影；

354）计算重投影的估计值与实际测量值的残差其中表示第j个空间点在第i像面的投影估计值，x_i,j表示所测量的投影点的值；

355）判断当前是否有估计点的残差大于预设阀值σ，如果是，则执行步骤356），否则执行步骤357）；

356）将特征链上的残差大于预设阀值σ的估计点移除，回到步骤352）；

357）输出使得效用函数E最小的对位参数。

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

1、本发明不依赖胶体金标记，且能够提高对位准确度。

2、本发明能够在计算量相对较小的前提下，实现不依赖胶体金标记的高准确度对位。

附图说明

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

图1示出了获取样品的电子断层投影图像序列的示意图；

图2示出了本发明一个实施例的电子断层图像对位方法流程图；

图3示出了本发明一个实施例中所使用的增量式捆绑调整算法流程图；

图4示出了线粒体投影序列的互相关方法对齐后的三维重构结果中三维样品中心部分的x-y截面；

图5示出了线粒体投影序列的互相关方法对齐后的三维重构结果中三维样品中心部分的y-z截面；

图6示出了线粒体投影序列的本发明实施例方法对齐后的三维重构结果中三维样品中心部分的x-y截面；

图7示出了线粒体投影序列的本发明实施例方法对齐后的三维重构结果中三维样品中心部分的y-z截面；

图8示出了囊泡组织投影序列的互相关方法对齐后的三维重构结果中三维样品中心部分的x-y截面；

图9示出了囊泡组织投影序列的互相关方法对齐后的三维重构结果中三维样品中心部分的y-z截面；

图10示出了囊泡组织投影序列的本发明实施例的方法对齐后的三维重构结果中三维样品中心部分的x-y截面；

图11示出了囊泡组织投影序列的本发明实施例的方法对齐后的三维重构结果中三维样品中心部分的y-z截面；

图12示出了中心粒投影序列的互相关方法对齐后的三维重构结果中三维样品中心部分的x-y截面；

图13示出了中心粒投影序列的互相关方法对齐后的三维重构结果中三维样品偏上部分的x-y截面；

图14示出了中心粒投影序列的互相关方法对齐后的三维重构结果中三维样品中心部分的y-z截面；

图15示出了中心粒投影序列的互相关方法对齐后的三维重构结果中三维样品中心部分的x-y截面；

图16示出了中心粒投影序列的互相关方法对齐后的三维重构结果中三维样品偏上部分的x-y截面；

图17示出了中心粒投影序列的互相关方法对齐后的三维重构结果中三维样品中心部分的y-z截面；

图18示出了中心粒投影序列的有胶体金对位方法对齐后的三维重构结果中三维样品中心部分的x-y截面；

图19示出了中心粒投影序列的有胶体金对位方法对齐后的三维重构结果中三维样品偏上部分的x-y截面；

图20示出了中心粒投影序列的有胶体金对位方法对齐后的三维重构结果中三维样品中心部分的y-z截面；

图21示出了基于小鼠线粒体投影序列的两种三维重构的再投影图像与原始序列的相关系数曲线的对比；

图22示出了示出了基于猪动脉内皮细胞囊泡组织投影序列的两种三维重构的再投影图像与原始序列的相关系数曲线的对比；

图23示出了基于中心粒投影序列的多种三维重构的再投影图像与原始序列的相关系数曲线的对比。

具体实施方式

发明人经过深入研究，发现传统的基于特征的图像对位方案具有以下问题：1.传统的基于特征的图像对位方案未能充分地利用图像信息。它们往往只是简单的将图像进行小区域的划分，或是通过计算梯度极值寻找特征区域，这类利用简单的局部自相关方法抽取图像信息的做法会造成信息量损失，导致抽取的有效特征点偏少,无法有效恢复仿射旋转等信息,而发明人研究发现，仿射旋转等信息的损失，是造成传统的基于特征的图像对位方案准确度不佳的重要原因之一。2.传统的基于特征的图像对位方案中，错误特征点无法有效去除。其原因一是在抽取图像信息时的信息损失，另一重要原因是错误点去除方法不够鲁棒。3.传统的基于特征的图像对位方案中，往往采用传统的回归方法或是牛顿法来进行参数估计，而三维参数估计及投影数据的拟合是一个非线性的过程，特征点偏少可导致三维参数估计置信度降低，而错误的特征点更是可能将最终估计值引入一个错误的局部值，因此传统的回归方法或是牛顿法很难有效估计样品的三维参数。

基于以上分析，本发明提出了一种改进的无胶体金标记的电子断层图像对位方法。为了使本领域的技术人员更好地理解本发明的技术方案，下面结合附图和具体实施例对本发明作进一步的详细说明。

为便于理解，首先简要介绍图像序列的拍摄过程。图1示出了获取样品的电子断层投影图像序列的示意图，

如图1所示，将样品1（例如生物组织切片）置于Transmission electronmicroscopy(TEM)成像系统中，相机（图1中未示出）围绕样品旋转，从多个不同角度对该样品进行成像，从而获得多个不同角度的电子断层投影图像（即图1（a）部分下方所示出的一系列图像），这些不同角度的电子断层投影图像就组成了本发明所要处理的图像序列。图1（b）部分则示出了同一样品空间点与各投影图像中的点的对应关系。对图像序列进行分析，能够重构原样品的三维结构。而电子断层图像对位就是重构原样品的三维结构的关键技术之一。

图2示出了本发明的一个实施例的无胶体金标记的电子断层图像对位方法的流程图。如图2所示，该无胶体金标记的电子断层图像对位方法包括以下四个步骤：

步骤一：投影图像序列特征点提取。即接收输入的投影图像序列，从每张投影图像序列中抽取出具有一系列稳定特征的近似于点的微小区域，这些区域的特征在图像序列中稳定存在，这些被抽取出的近似于点的微小区域就是特征点。特征点提取就是抽取图像中的特征点并记录特征点坐标和与其相对应的特征描述向量的集合。

本实施例中，用SIFT算子对每幅图像提取特征点，每个特征点生成一个128维的向量。SIFT算子是被美国专利文献US6,711,293所公开的一种提取图像特征的算子。在过去传统的基于特征的对位方法中，使用过很多其他算子。但是，它们有一个共同的特点，就是对特征点的特征描述仍是基于互相关。而互相关对于仿射变换、图像旋转等鲁棒性不高。发明人研究发现，SIFT在特征向量的生成过程中采用了一种循环编码策略，很好的克服了对位图像的仿射变换带来的误差及电镜图像的低信噪比、带来了更好的旋转识别、并且有效地利用差分金字塔的多尺度信息，非常适合用于无胶体金标记的电子断层图像对位。当然，在其它实施例中，除了SIFT算子以外，还可以采用其它的能够保留仿射变换信息的图像特征提取算法来提取投影图像序列的特征点，例如SURF(Speeded Up Robust Features)法。

步骤二：基于所记录的各特征点的特征描述向量，对投影图像序列的特征点进行匹配。即对从步骤一中得到的各个图像中的特征点进行关系映射，以确定来自序列中不同图像的特征点是否代表同一个区域，从而得到不同图像的特征点之间的对应关系。

根据本发明的一个实施例，该步骤包括：两幅图像间的匹配和图像序列的匹配两个过程，下面分别介绍。

（1）两幅图像间的匹配算法

特征点匹配是生成特征链的基础。根据本发明的一个实施例，对于SIFT的向量特征，使用欧式距离来判定两个特征点是否相似。对于电子断层图像序列，前后两张图像中相对应的点的偏移区域具有一定的先验知识，并且对于一些具有对称结构或自相似的生物样品图像，进行全局搜索本身会造成一定程度的错误匹配。因此，根据本发明的一个实施例，提出了一个基于二维数据结构的局部搜索方法完成两幅图像间的匹配，包括下列步骤：

1、有输入A图与B图的特征点集，对于A图中的每个向量，在B图的某对应位置取一窗口函数，输出窗口内的点集。

2、计算A中向量与窗口向量的欧式距离，并计录是否为当前最小或次小距离。

3、若最小距离小于次小距离乘以阈值，则匹配成功，计录当前匹配。

4、重复1～3过程，直到所有特征点被遍历。

5、对于所匹配的特征点对，使用随机抽样一致（RANSAC）原则基于二极几何去除错误匹配对。

上述两幅图像间的匹配中，窗口函数利用图像序列的先验信息做为映射依据并由一二维数据结构支撑。阈值的选取具有一定的经验性，通过大量的实验，发现阈值的优选区间为0.4～0.6。事实上，即便阈值处于优选区间，错误匹配及一配多这样的情况也是难免的，因此，根据本发明的一个实施例，进一步地使用随机抽样一致方法来去除所得的匹配集合中错误的匹配对。

来自不同图像的特征点集的两两匹配(即每两副图像都完成一次匹配)是一个比较耗时的工作。特征点对应于实际样品的某一特征区块（特征区块是近似于空间点的微小区块，下文不再赘述），在多个不同角度对其进行投影，所得的投影结果（即图像）中均可能含有该特征区块的一系列特征，因此特征点通常不会孤立存在，而是在多个角度所对应的投影图像中存在相互关联的对应于同一实际样品特征区块的投影特征点（即文中所称特征点）。所以特征点是稳定可靠的并且在投影图像之间具有传递性；并且，图像序列是渐变的，一张投影图像中出现的一个被匹配的特征点，其相对应的特征点最有可能在其附近的邻居图像中出现，即，对于一个在nth图像（表示第n幅图像）的特征点，其对等点将很有可能出现在n-1th、n+1th、n+2th等图像（n-1th、n+1th、n+2th分别表示第n-1、n+1、n+2幅图像）上。基于上述分析，若要寻找nth与n+2th的匹配对，则匹配nth与n+1th、n+1th与n+2th即可。

因噪声或形变，匹配对的传递方式不会像上述分析那么理想，所以在一个实施例中，提出了以下的匹配策略：

1、初始化step＝1，当step≤MAX_STEP时；

2、对于序列中的每幅图像，用前文中描述的两幅图像间的匹配方法来匹配nth与(n+step)th的特征集；

3、step值加1；

4、重复1～3步骤，直到不再满足step≤MAX_STEP的条件。

理论上，MAX_STEP应该越大越好。但是，考虑到计算量和计算速度的要求，MAX_STEP应该适中。在一个优选实施例中，取MAX_STEP为3，如前文所述，图像序列是渐变的，一张投影图像中出现的一个被匹配的特征点，其相对应的特征点最有可能在其附近的邻居图像中出现，因此，在对于一幅图像，在临近的3幅邻居图像中寻找可匹配到大部分的特征点。并且由于特征点具有传递性，在对下一幅图像进行匹配时，对应于实际样品同一特征区块的特征点的匹配关系够自然地传递，因此能够获得足够的信息以进行后续的计算。

（2）投影图像序列特征链生成

投影图像序列特征链生成就是对所得到的特征点对应关系进行跟踪，以最终形成一个链状的传递关系，为便于描述，本文中将这种链状的传递关系称为特征链。

在求得两两图像之间匹配的特征点，也即匹配对。本步骤中将这些匹配对中拥有相同特征点的匹配对串联起来，形成特征链。例如有(1,2)-(2,3)和(1,2)-(3,4)两个匹配对，则可串联成(1,2)-(2,3)-(3,4)，意为此特征链包括第1幅图像的第1个特征点、第2幅图像的第3个特征点和第3幅图像的第4个特征点。

本实施例在具体实现上将两两匹配的特征点通过填入一个三维结构生成可快速检索的特征链。具体地说，在存储特征点时，在一个实施例中，使用存储序列化的坐标信息的平衡二叉树作为同一投影图像的特征点的存储集合，其中平衡二叉树中的每个结点代表该投影图像的一个特征点并存储相应的数据。另一方面，基于特征点的匹配关系，每个平衡二叉树的每个结点还存储与其它平衡二叉树的相应结点的链接关系，这样就形成了三维结构。本实施例的三维结构能够在匹配对跨度有限（不在间隔过大的投影图像之间进行特征点匹配）的情况下，利用特征点的传递性保证特征链长度，因此，能够在保证对位质量的前提下较好地控制计算量，提高图像对位的运算速度。

步骤三：基于步骤二获得的特征链，估计对位参数及其它空间参数，估计对位参数及其它空间参数就是通过已经得到的图像特征链还原电镜的成像过程。在获得图像序列中各图像所对应的对位参数后，即可得到对位后的图像序列。

下面结合实施例对本步骤做进一步地描述。

（1）投影过程及最优化目标建模

对于TEM成像，因其场景的深度大小变化远小于场景距照像机的距离，所以选用仿射模型是极为合适的。在此引入如下的射影模型：

$\left(\begin{array}{c}u\\ v\\ z_c\end{array}\right)=K(I,0)\left(\begin{array}{cc}R& t\\ 0& 1\end{array}\right)\left(\begin{array}{c}X\\ Y\\ Z\\ 1\end{array}\right)---\left(1\right)$

其中K为相机的内参矩阵，(I,0)为投影变换矩阵，R是3×3正交矩阵，表示相机在空间的旋转，可以写为欧拉角形式R_αR_βR_γ其中：

$R_{\mathrm{\α}}=\left(\begin{array}{ccc}1& 0& 0\\ 0& \cos \mathrm{\α}& \sin \mathrm{\α}\\ 0& -\sin \mathrm{\α}& \cos \mathrm{\α}\end{array}\right)R_{\mathrm{\β}}=\left(\begin{array}{ccc}\cos \mathrm{\β}& 0& -\sin \mathrm{\β}\\ 0& 1& 0\\ \sin \mathrm{\β}& 0& \cos \mathrm{\β}\end{array}\right)R_{\mathrm{\γ}}=\left(\begin{array}{ccc}\cos \mathrm{\γ}& \sin \mathrm{\γ}& 0\\ -\sin \mathrm{\γ}& \cos \mathrm{\γ}& 0\\ 0& 0& 1\end{array}\right)$

t为3×1向量，近似表示样品在空间中的偏移（或

为相机的空间方位）。(u,v,z_c)^T对于空间点的齐次形式的投影坐标，即真实图像上点坐标可由(u/z_c,v/z_c)表示。(X,Y,Z,1)^T为空间点的齐次坐标。通常，默认K在整个成像过程中保持不变，在整个参数优化过程中，需要估计的对位参数和其它空间参数包括：旋转矩阵R，空间偏移t，样品空间点(X,Y,Z,1)^T。

在参数估计的过程中，初始值按下述方法设置：将α设为读角器读出的角度值，β和λ分别设定为零，即默认样品在操作过程中除了样品平面在旋转轴的旋转外是无其它旋转的。同理，对于t向量，选取0度倾斜角相机为坐标原点，其它相机依次跟据旋转角度及相机深度确定t的初值。而K取固定值。对于每一个空间点(X,Y,Z,1)^T，其初值由其对应的位于低角度的投影点的反投影获得。

为估计R、t、(u,v,z_c)^T、(X,Y,Z,1)^T四组参数，采用如下形式的效用函数：

$E=\underset{i,j}{\mathrm{\Σ}}{({\mathrm{Proj}}_i\left({\tilde{X}}_j\right)-x_{i,j})}^2{\mathrm{\δ}}_{i,j}---\left(2\right)$

是空间特征点的坐标，

是上述全投影过程的输出值。x_i,j为测量的特征点空间坐标。为方便计算，上述公式中，

及x_i,j的值都以欧式坐标表示。δ_i,j表明点j在投影i中是否可见，如可见则δ_i,j取1，如不可见则δ_i,j取0。对位参数估计的过程就是寻找使得E最小的参数R、t。

上面的投影变换中，投影参数（即对位参数）及特征点的位置都是未知的，并且它们相互影响。因此，若同时估计上述未知量，参数估计问题就成为了非线性最优化问题。本实施例中有大量的投影点做为参数输入，而相对应的，投影参数较少，在一个实施例中，采用增量式的捆绑调整的方法进行参数估计，从而最大程度地利用已知的投影点，在保证对位精度的同时大幅提升计算速度。下面结合实施例对增量式的捆绑调整进行详细说明。

（2）增量式的捆绑调整

增量式的捆绑调整解决了三个重要的问题：1、输入变量中存在大量特征点及空间点，如何进行求解；2、每次参数估计过程的鲁棒性保证；3、特征链过短情况下的捆绑调整。

对于第1个问题，采用稀疏捆绑调整（sparse bundle adjustment）技术来求解参数，采用稀疏捆绑调整是一种求解非线性模型的算法,基于Levenberg-Marquardt算法（简称为LM算法）。使得f为一个映射向量p∈R^m至估计向量的函数，即

LM算法最基本的思想是在f的邻近区域p内，对于一个小的距离||δ_p||(||·||表示L2正则)，进行仿射近似。f可近似地估计为f(p+δ_p)≈f(p)+Jδ_p,其中J是f的Jacobian矩阵。当表示映射的函数f的Jacobian矩阵中包含大量零元素时，稀疏捆绑调整通过稀疏化版本的LM算法存储及计算解决待估计值过多的问题。通过将效用函数表示为

其中ε为残差向量,通过交换等式，将所有待评估参数表示为向量P∈R^M以适应LM算法形式。

关于第2个问题，即使在特征点匹配阶段去除了大多数的两两不一致，但仍会有一部分异常点没有被排除，此时，为了保证方法的鲁棒性，可在每次捆绑调整后，对所得空间参数进行反投影，计算其残差

其中

表示第j个空间点在第i像面的投影估计值，x_i,j表示所测量的投影点。如果一条链上的某估计点的残差大于预设阀值σ，在下次计算中链上该点将被移除：σ_thre＝min{σ_min,σ_med,σ_max},σ_max及σ_min为预设的上下限,σ_med是计算所得残差的中值。捆绑调整中，每次循环执行一次参数估计及异常值检测，直到再无异常值被检测出来。

对于第3个问题，因特征提取所得到的链长比较短，无法一次对所有的投影进行参数标定，此时，可以用一种增量式的方法从低角度开始进行计算以解决链过于短的问题。具体如下：首先，选取有足够Baseline(基线)的，且通过图像的链数足够多的两张图像做为种子序列并对其进行参数标定；之后，再次地，在原先估计中加入一些新的图像（通常为已估计的图像的邻居，它们可以保证有足够的已估点的对等点），先每次加入一个相机（即一个投影角度），保证原先参数固定（即已计算所得的空间点及相机参数固定）并估计新加入图像对应的相机的参数（因捆绑调整算法只保证达到局部最优，因此先对单相机进行参数调整可优化初值），再对所有图像进行捆绑调整并更新所有参数；重复上述过程直到再无可估计的图像。这样，就克服了在电镜数据中特征链长较短的问题。

下面根据本发明的一个实施例，给出一个具体的增量式的捆绑调整估计方法，该估计方法基于已知的图像序列及其特征链的前提下，估计每幅图像所对应的相机空间旋转参数，及其图像序列中各特征点所对应的样品的各个空间点的位置坐标，具体流程如图3所示。这个流程主要包含两个部分，一个主流程和一个附流程。图3中标记a)b)c)主要对应了上文中需要解决的三个问题的解决方法在流程中的位置。图3中主流程即增量式的捆绑调整估计方法的主要流程，而子流程则用于解释主流程中的捆绑调整的具体过程。

增量式的捆绑调整估计方法的主流程包括下列步骤：

1）选择初始图像对：选取通过图像的特征链数目足够多的两张图像做为种子序列，该种子序列为初始的待估计序列；

2）对图像进行捆绑调整：设置种子序列中低角度的图像所对应的相机空间旋转矩阵R、样品在空间中的偏移t、以及每一个空间点(X,Y,Z,1)^T的初值。其中，R的初值通过初始旋转角度确定，t的初值通过投影图像对应的相机旋转角度和相机深度共同确定，然后根据图像序列中各特征点的值，基于R、t的初值进行反投影，获得所对应的每一个空间点(X,Y,Z,1)^T的初值。将图像匹配的特征点代入前文的公式(2)，对序列中的每幅图像，寻找使得效用函数E最小的参数R、t。

3）选择合适图像加入调整流程：选择与已估计序列共有较多特征链的图像加入当前的待估计序列（即将新的图像加入调整流程），得到新的待估计序列。

4）初始化捆绑调整中各参数：将待估计序列中其它图像已知或已估计出的相机空间旋转矩阵R、样品在空间中的偏移t作为本次估计的初值。新出现的每一个空间点(X,Y,Z,1)^T初值由其对应的位于低角度的投影点的反投影获得，反投影依据本步骤已知或已估计出的R,t的值进行。

5）捆绑调整:将待估计序列中的各图像的参数及空间点值的初始估计代入前文的公式(2)，对序列中的每幅图像，寻找使得效用函数E最小的参数R,t。在一个优选实施例中，本步骤的捆绑调整中引入了抑制异常点的机制，将在下文的子流程中做进一步地说明。

6）重复执行步骤3）至5）直至所有的图像均已加入序列。这样就估计出了所有图像所对应的相机空间旋转矩阵R和样品在空间中的偏移t，以及图像序列中所有特征点的位置(X,Y,Z,1)^T。

根据一个优选实施例，所述步骤5）中，抑制异常点的捆绑调整的子流程包括下列步骤：

51）输入所有的测量得到的特征点x_i,j及初始化的参数变量。

52）求解更新每幅图像的平移变量t空间旋转变量R及每一个空间点(X,Y,Z,1)^T。其中，采用前文中的稀疏捆绑调整技术来降低计算量。

53）对所得空间参数进行重投影。

54）计算重投影的估计值与实际测量值的残差

其中表示第j个空间点在第i像面的投影估计值，x_i,j表示所测量的投影点。

55）判断当前是否有估计点的残差大于预设阀值σ，如果是，则执行步骤56），否则执行步骤57）。

56）将特征链上的残差大于预设阀值σ的估计点移除，回到步骤52）。

57）输出使得效用函数E最小的参数R,t。

采用三组数据分别基于两种对位方法进行三维重构测试。其中，第一组数据是小鼠线粒体投影序列；第二组数据是猪动脉内皮细胞囊泡组织投影序列；第三组数据是中心粒投影序列。两种对位方法分别是本发明一个实施例的对位方法和现有技术中典型的互相关对位方法。在对位完成后，基于两种方法的对位后图像序列均采用加权背投影（Weighted Back-Projection）进行三维重构，这样三维重构的精度就反映了对位的精度。

上述三组数据的三维重构结果显示，与传统的互相关对齐方法相比，基于本发明的一个实施例的对位方案的噪声及人工尾迹显著降低。具体地，图4示出了线粒体投影序列的互相关方法对齐后的三维重构结果中三维样品中心部分的x-y截面；图5示出了线粒体投影序列的互相关方法对齐后的三维重构结果中三维样品中心部分的y-z截面；图6示出了线粒体投影序列的本发明实施例方法对齐后的三维重构结果中三维样品中心部分的x-y截面；图7示出了线粒体投影序列的本发明实施例方法对齐后的三维重构结果中三维样品中心部分的y-z截面；观察上述图中方框部分及其他背景部分，可以明显地发现本发明实施例的方法所得到的线粒体投影序列的重构结果细节更清晰，且人工尾迹更少。

图8示出了囊泡组织投影序列的互相关方法对齐后的三维重构结果中三维样品中心部分的x-y截面；图9示出了囊泡组织投影序列的互相关方法对齐后的三维重构结果中三维样品中心部分的y-z截面；图10示出了囊泡组织投影序列的本发明实施例的方法对齐后的三维重构结果中三维样品中心部分的x-y截面；图11示出了囊泡组织投影序列的本发明实施例的方法对齐后的三维重构结果中三维样品中心部分的y-z截面；观察上述图中方框部分及其他背景部分，可以明显地发现本发明实施例的方法所得到囊泡组织投影序列的重构结果细节更清晰，且人工尾迹更少。

图12示出了中心粒投影序列的互相关方法对齐后的三维重构结果中三维样品中心部分的x-y截面；图13示出了中心粒投影序列的互相关方法对齐后的三维重构结果中三维样品偏上部分的x-y截面；图14示出了中心粒投影序列的互相关方法对齐后的三维重构结果中三维样品中心部分的y-z截面；图15示出了中心粒投影序列的互相关方法对齐后的三维重构结果中三维样品中心部分的x-y截面；图16示出了中心粒投影序列的互相关方法对齐后的三维重构结果中三维样品偏上部分的x-y截面；图17示出了中心粒投影序列的互相关方法对齐后的三维重构结果中三维样品中心部分的y-z截面；另外，对于中心粒投影序列，增加了一组基于有胶体金对位的重构数据。18示出了中心粒投影序列的有胶体金对位方法对齐后的三维重构结果中三维样品中心部分的x-y截面；图19示出了中心粒投影序列的有胶体金对位方法对齐后的三维重构结果中三维样品偏上部分的x-y截面；图20示出了中心粒投影序列的有胶体金对位方法对齐后的三维重构结果中三维样品中心部分的y-z截面；中心粒投影序列这组样品数据比较小且采集的比较良好，但还是可以明显地发现本发明实施例的方法以及有胶体金方法所得到的结果细节更清晰，且人工尾迹更少（胶体金越扭曲,说明效果越差）。有胶体金的对位方法是现在ET图像对位中最为良好的方法，但适用面较窄。例如上面的线粒体和囊泡组织的数据，就无法做有胶体金方法的对位。

进一步地，对三维重构结果进行重投影获得虚拟图像序列，计算该虚拟图像序列与原始图像序列的相关系数曲线，相关度越高，说明三维重构结果越准确，即对原始图像的对位精度越高。图4示出了基于小鼠线粒体投影序列的两种三维重构的再投影图像与原始序列的相关系数曲线的对比；图5示出了示出了基于猪动脉内皮细胞囊泡组织投影序列的两种三维重构的再投影图像与原始序列的相关系数曲线的对比；图6示出了基于中心粒投影序列的多种三维重构的再投影图像与原始序列的相关系数曲线的对比。可以看出，参与测试的本发明实施例的方案的相关系数显著高于互相关对齐方法。在基于中心粒投影序列的三维重构中，参与测试的本发明实施例的方案的相关系数与有胶体金标记的对位方法基本相当。

对位时间计算速度视拍摄所得的数据大小而定。现今ET图像大小一般为1024×1024和2048×2048两种。在测试中，基于当前主流的PC平台（CPU为英特尔酷睿i5，内存4G），本发明实施例的方法在1024×1024的数据集上的运行时间为10分钟左右；在2048×2048的数据集上的运行时间约为30分钟左右。可以看出，本发明能够在实现不依赖胶体金标记的高准确度对位的同时，将计算量控制在较小范围，从而达到较快的计算速度。

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

Claims (11)

1.一种电子断层图像序列对位方法，包括下列步骤：

1）从每幅电子断层投影图像中提取特征点，所述特征点能够保留图像仿射变换信息；

2）对不同电子断层投影图像的特征点进行匹配，获得多幅电子断层投影图像的特征点之间的对应关系；所述对应关系是对应于同一实际样品特征区块的不同角度电子断层投影图像中的特征点之间的关联关系；

3）基于仿射模型，根据步骤2）所获得的多幅电子断层投影图像的特征点之间的对应关系，估计各电子断层投影图像所对应的对位参数。

2.根据权利要求1所述的电子断层图像序列对位方法，其特征在于，所述步骤1）中，采用SIFT算子或者SURF算法从每幅电子断层投影图像中提取特征点。

3.根据权利要求1所述的电子断层图像序列对位方法，其特征在于，所述步骤2）包括下列子步骤：

21）采用局部搜索法，对每幅电子断层投影图像，将其与其它电子断层投影图像进行两幅电子断层投影图像之间的特征点匹配；

22）将步骤21）中所匹配的特征点匹配对串联，得到对应于被投影样品相同特征区块的各个不同角度的电子断层投影图像的特征点组成的特征链。

4.根据权利要求3所述的电子断层图像序列对位方法，其特征在于，所述步骤21）中，采用随机抽样一致原则，基于二极几何去除错误匹配对。

5.根据权利要求3所述的电子断层图像序列对位方法，其特征在于，所述步骤21）中，每幅所述电子断层投影图像仅与邻近的电子断层投影图像进行特征点匹配。

6.根据权利要求4所述的电子断层图像序列对位方法，其特征在于，所述步骤3）中，各电子断层投影图像所对应的对位参数包括：每幅电子断层投影图像所对应的相机空间旋转参数、每幅电子断层投影图像所对应的样品在空间中的偏移参数以及各特征点所对应的样品原始特征点空间坐标。

7.根据权利要求6所述的电子断层图像序列对位方法，其特征在于，所述步骤3）中，首先选取电子断层投影图像序列中的部分图像，根据已知的电子断层投影图像的特征点之间的对应关系，基于误差最小原则对当前所选取的图像的所述对位参数进行估计；然后再加入新的图像，根据前一步所估计的所述对位参数以及已知的电子断层投影图像的特征点之间的对应关系，基于误差最小原则对加入新的图像后的当前所选取的图像的所述对位参数进行估计和更新；不断重复直至所有图像均被加入且获得所有图像的所述对位参数。

8.根据权利要求7所述的电子断层图像序列对位方法，其特征在于，所述步骤3）中采用稀疏捆绑调整技术估计所述对位参数。

9.根据权利要求6所述的电子断层图像序列对位方法，其特征在于，所述步骤3）包括下列子步骤：

31）根据穿过图像的特征链数目选取两张图像做为初始待估计序列；

32）设置初始待估计序列的对位参数的初值，所述对位参数包括图像所对应的相机空间旋转矩阵R、样品在空间中的偏移t、以及对应于图像中匹配特征点的样品空间点(X,Y,Z,1)^T，其中，R的初值通过初始旋转角度确定，t的初值通过投影图像对应的相机旋转角度和相机深度共同确定，根据待估计序列中各特征点的值，基于R、t的初值进行反投影，获得所对应的每一个空间点(X,Y,Z,1)^T的初值；

对当前待估计序列中的每幅图像，寻找使得效用函数E最小的参数R、t，

$E=\underset{i,j}{\mathrm{\Σ}}{({\mathrm{Proj}}_i\left({\tilde{X}}_j\right)-x_{i,j})}^2{\mathrm{\δ}}_{i,j}$

是第j个空间点的估计值，

是对第j个空间点的估计值进行第i像面的重投影的输出值，x_i,j为第j个空间点在第i像面的投影的测量值，δ_i,j表明点j在投影i中是否可见，如可见则δ_i,j取1，如不可见则δ_i,j取0；

33）将一幅新电子断层图像加入，得到新的待估计序列；

34）以当前已知或已估计出的电子断层投影图像所对应的对位参数作为初值，通过位于低角度的投影点的反投影获得种子序列中新出现的空间点的初值；

35）将当前种子序列中的每幅图像，寻找使得所述效用函数E最小的参数R、t，并据此更新当前待估计序列中每幅图像所对应的参数R、t；

36）重复执行步骤33）至35）直至所有的图像均已加入待估计序列。

10.根据权利要求9所述的电子断层图像序列对位方法，其特征在于，所述步骤32）、35）中，采用稀疏捆绑调整技术来求解使得所述效用函数E最小的参数R、t。

11.根据权利要求10所述的电子断层图像序列对位方法，其特征在于，所述步骤35）包括下列子步骤：

351）输入当前待估计序列中所有测量得到的特征点x_i,j及初始化的对位参数变量；

352）求解使得所述效用函数E最小的每幅图像的对位参数并用所求得的对位参数替换对位参数变量的初值；

353）根据所得的对位参数进行重投影；

354）计算重投影的估计值与实际测量值的残差

其中

表示第j个空间点在第i像面的投影估计值，x_i,j表示所测量的投影点的值；

355）判断当前是否有估计点的残差大于预设阀值σ，如果是，则执行步骤356），否则执行步骤357）；

356）将特征链上的残差大于预设阀值σ的估计点移除，回到步骤352）；

357）输出使得效用函数E最小的对位参数。

Images