CN117557623B - 一种冷冻电镜图像序列的精准快速对齐方法 - Google Patents

Abstract

本发明公开了一种冷冻电镜图像序列的精准快速对齐方法，属于图像处理技术领域。本发明提出了一种在非线性最小二乘中加入迭代式的权重估计算法，该算法基于L₂范数进行优化，通过对每次优化后标记点残差的计算赋予每个标记点不同的权值，为残差较小的点赋予较低的权重，为残差较高的点赋予较低的权重或不赋予权重，从而减小离群点的影响，提高了算法的鲁棒性。本发明避免了离群点对投影参数优化造成的影响，保证最终对齐的准确性，为对齐之后的重建提供保障。本发明在冷冻电镜图像序列对齐的捆绑调节过程中，在存在离群点的情况下仍可实现快速精准的投影参数优化，实现图像序列对齐。

Description

一种冷冻电镜图像序列的精准快速对齐方法

技术领域

本发明涉及图像处理技术领域，特别涉及一种冷冻电镜图像序列的精准快速对齐方法。

背景技术

冷冻电子显微镜技术在当前结构生物学领域中发挥了重要的作用，可以通过与三维重建技术结合，还原生物大分子的三维结构。在结构生物学的三维结构重建领域中，通常需要获得原物体在一定旋转范围内的离散角度的投影图像序列(即倾斜序列)，并通过该序列还原出物体原本的三维结构。在获得投影图像序列到反求解获得三维结构的过程中，一个核心的步骤是对投影图像序列进行对齐，其中基于胶体金标记的图像序列对齐是最精准的，而胶体金的图像对齐需要构造完整的点链并进行投影参数标定。点链的完整度以及图像序列的投影参数标定对样品三维重构结果精度的好坏有着直接影响，因此是至关重要的环节。

针对上述问题，图像序列的投影参数标定最常用的方法是捆绑调整，捆绑调整实际上是一个大型的稀疏几何参数估计问题，其中参数是捆绑在一起的三维特征坐标和相机位姿。捆绑调整解决的问题可以归约成非线性最小二乘问题，目标是最小化重投影误差的L₂范数，在L₂范数中参数是相机参数以及三维点坐标。尽管捆绑调整在近二十年内迅速发展并逐步优化，但该算法的鲁棒性依旧是待优化的部分，即使数据集中仅有一小部分离群点，也会对最终的结果产生较大的影响。传统的数据处理流程在高角度的投影图像的胶体金跟踪中会出现很大一部分因跟踪误差过大而被错误抛弃的胶体金，导致胶体金点链缺失的情况，造成参数估计中的偏差。将捆绑调节算法中最小化重投影误差的L₂范数换成L₁范数可以有效的减小离群点对投影参数标定的影响，但快速且完备地生成的胶体金点链仍是当前领域的一大难点。

虽然L₁范数可以减少离群点所占权重值，但是求解最小化L₁范数的重投影误差问题最困难的就是其并不连续可导，因此无法用常用的最小二乘法直接计算。为了最小化L₁范数，有人提出采用对数障碍内点法，首先将L₁范数问题转化为近似等价的线性规划的形式，然后将该线性规划的限制条件放置到对数障碍函数中，在线性规划的目标函数前添加近似精度的参数t，组合这两部分得到一个新的连续可导、且与原函数近似等价的函数，最后通过更新策略，在不断获得牛顿步长的循环中，选择可以使得原目标函数最小的投影参数，获得在误差适合的范围内最精确的参数解。但此方法速度过慢，即便是面对几百个基准点的投影参数优化都需要十几分钟，常见的冷冻电镜图像序列投影参数优化都是成千上万个标记点，因此基于L₁范数的无法满足时间上的要求。

发明内容

本发明的目的在于提供一种冷冻电镜图像序列的精准快速对齐方法，以弥补现有技术的不足。

在基于胶体金的冷冻电镜图像对齐时，可能会因为基准点标记算法的误差造成偏移较大的基准点存在甚至发生丢失现象造成点链不完整，特别是高角度的图像，基准点的坐标误差更大，并且在优化过程中如果投影参数的初始化过差可能无法进行有效的优化求解，此时会对最后的投影参数标定造成较大的影响。

为达到上述目的，本发明是通过以下技术方案实现的：

一种冷冻电镜图像序列的精准快速对齐方法，该方法包括如下步骤：

S1：获取冷冻电镜图像数据，对投影图像序列上的胶体金进行识别；

S2：对低角度(-50°～+50°)图像上的胶体金进行匹配跟踪；

S3：梳理投影图像序列中的胶体金匹配对，生成贯穿整个投影图像序列的胶体金链，即二维投影坐标对应关系；

S4：基于正交投影模型进行投影参数的标定，首先进行初始化投影参数，之后通过捆绑调整算法求解得到初步的低角度下的胶体金的三维坐标预估计，通过三维坐标的预估计来进一步估算高角度投影图像上的胶体金虚拟坐标分布；

S5：通过与未加入的高角度投影图像相匹配，从低角度投影向高角度投影进行迭代式的序列跟踪，利用赋权的迭代式捆绑调节算法实现投影参数再标定，最终实现对整个序列的可用胶体金跟踪；

S6：投影参数标定完成后，即得到优化后的投影参数后，将投影参数与未对齐的图像序列进行计算，得到对齐后的冷冻电镜图像序列，完成对齐。

进一步的，所述S4中，使用初始的正交投影模型，其数学建模如下：

s表示投影图像的尺度的变化，γ投影图像的片内旋转，α表示投影图像在成像过程中的绕转轴的倾转角度，β表示投影图像在成像过程中转轴的倾转角度，(u,v)表示测量所得的胶体金坐标；

进一步的，所述S4中提到的投影参数初始化，其计算流程如下：

对s，α，β，γ以及平移变量t进行初始化；因为图像为刚性变化，因此s为1，α是投影图像在成像过程中的绕转轴的倾转角度，接近为0，因此初始化为0；主要对平移变量tx，ty以及γ进行计算初始化计算；其计算步骤如下：

(1)把图像的中心确定，初始化度数最低图像的平移变量tx,ty为0；利用如下的仿射变换公式计算每两张图片之间的仿射变换关系：

之后从中心开始向左向右，分别初始化，从而计算每张图片的tx，ty；

(2)在得到初步的平移变量之后，对γ进行初始化，从-90°到90°寻找最佳的γ，以5°为一个间隔进行寻找；

以当前的角度进行平移变量的计算：

之后对根据不同角度图像上的标记点利用γ的数值进行片内旋转计算，并与之前计算的tx，ty进行偏移更新计算；将相同角度图像上的标记点利用最小二乘进行优化；之后计算每个点的优化残差。寻找残差最小的角度作为γ；

(3)利用之前得到的γ进行更新平移变量tx，ty；

(4)将γ的搜索空间缩小为-5°到5°，以1°为以间隔再次利用(2)步骤进行寻找最优的γ值，并更新tx，ty的值，之后在此取值下，对ty再次利用最小二乘进行更新，并记录当前的残差值；

(5)迭代进行步骤(4)直到残差不在减小停止，停止后利用当前的tx，ty以及γ值对整个序列进行最小二乘优化，得到最终解。

进一步的，所述赋权的迭代式捆绑调节算法，包括以下步骤：

(1)将相机参数与基准点的二维坐标作为输入；

(2)将标记点和相机参数进行分组，分组后将L₂范数作为目标函数对输入值进行优化；

(3)计算优化后的值的残差，利用残差计算每个基准点的归一化后的中位数的绝对偏差(MADN)，通过MADN的值为标记点进行分配权重；

(4)将带权的基准点再次利用L₂范数的目标函数进行迭代式的优化；

(5)当每个基准点的权重基本不再发生改变的时候优化停止，得到最终的相机参数与基准点的二维坐标。

更进一步的，所述捆绑调节目标是实现胶体金的估计值与测量值最小化：

其中表示估计所得到的空间三维胶体金坐标，Proj(·)为投影过程，x_i,j是测量得到的胶体金坐标δ_i,j说明第j个胶体金是否在第i张投影中可见。

更进一步的，所述步骤(3)中权重的计算具体如下：

d_i,j＝e_i,j/(k*MADN)

MADN＝median(|e_i,j-median(e)|)/0.6475

其中，e_i,j表示每个胶体金的残差，k是标记权重为零的阈值，w_i,j是每个胶体金的权重。

与现有技术相比，本发明的优点和有益效果是：

本发明提出了全新的图像对齐流程，实现了对高角度图像上标记点的补充保证了点链的完整性减少投影参数的标定偏差，并在投影参数标定时利用迭代求解的方式提供了可靠的初始化投影参数，之后利用一种在非线性最小二乘中加入迭代式的权重估计算法，该算法基于L₂范数进行优化，通过对每次优化后标记点残差的计算赋予每个标记点不同的权值，为残差较小的点赋予较低的权重，为残差较高的点赋予较低的权重或不赋予权重，从而减小离群点的影响，提高了算法的鲁棒性。相较于基于L₁范数为目标参数的优化，速度可以大大提高，保障了在具备准确性的同时达到L₂范数优化的速度。

本发明可以产生尽可能完整和胶体金坐标链，并避免了离群点对投影参数优化造成的影响，保证最终对齐的准确性，为对齐之后的重建提供保障。本发明在冷冻电镜图像序列对齐过程中实现了高角度图像序列的标记点补充，并在捆绑调节过程中，提供了准确的初始化估计并在存在离群点的情况下仍可实现快速精准的投影参数优化，实现图像序列对齐。

附图说明

图1是本发明迭代式基于胶体金标记的冷冻电镜图像序列对齐流程。

图2是基于L₂范数的迭代式的权重估计算法流程。

图3是Nitrosop3和Vibrio数据集对齐前后基准点轨迹的叠加在图像空间中的分布图；其中:(a)Nitrosop3原始基准点轨迹的叠加(像素数中的x-y坐标)；(b)Nitrosop3对齐基准点后的轨迹叠加；(c)Vibrio原始基准点轨迹的叠加(像素数中的x-y坐标)；(d)Vibrio对齐基准点后的轨迹叠加。

具体实施方式

下面结合实施例对本发明所述的技术方案作进一步地描述说明。

实施例1

一种冷冻电镜图像序列的精准快速对齐方法，如图1所示，该方法包括如下步骤：

S1：获取冷冻电镜图像数据，对投影图像序列上的胶体金进行识别；

S2：对低角度(-50°～+50°)图像上的胶体金进行匹配跟踪；

S3：梳理投影图像序列中的胶体金匹配对，生成贯穿整个投影图像序列的胶体金链，即二维投影坐标对应关系；

S4：基于正交投影模型进行投影参数的标定，首先进行初始化投影参数，之后通过捆绑调整算法求解得到初步的低角度下的胶体金的三维坐标预估计，通过三维坐标的预估计来进一步估算高角度投影图像上的胶体金虚拟坐标分布；

S5：通过与未加入的高角度投影图像相匹配，从低角度投影向高角度投影进行迭代式的序列跟踪，利用赋权的迭代式捆绑调节算法实现投影参数再标定，最终实现对整个序列的可用胶体金跟踪；

S6：投影参数标定完成后，即得到优化后的投影参数后，将投影参数与未对齐的图像序列进行计算，得到对齐后的冷冻电镜图像序列，完成对齐。

在一实施例，所述S4中，使用初始的正交投影模型，其数学建模如下：

s表示投影图像的尺度的变化，γ投影图像的片内旋转，α表示投影图像在成像过程中的绕转轴的倾转角度，β表示投影图像在成像过程中转轴的倾转角度，(u,v)表示测量所得的胶体金坐标。

在一个实施例中的投影参数初始化，其计算流程如下：

对s，α，β，γ以及平移变量t进行初始化。因为图像为刚性变化，因此s为1。因此α是投影图像在成像过程中的绕转轴的倾转角度，接近为0，因此初始化为0。主要对平移变量tx，ty以及γ进行计算初始化计算。其计算步骤如下：

(1)把图像的中心确定，初始化度数最低图像的平移变量tx,ty为0。利用如下的仿射变

换公式计算每两张图片之间的仿射变换关系：

之后从中心开始向左向右，分别初始化，从而计算每张图片的tx，ty。

(2)在得到初步的平移变量之后，对γ进行初始化。从-90°到90°寻找最佳的γ，以5°

为一个间隔进行寻找。因为

所以以当前的角度进行平移变量的计算：

之后对根据不同角度图像上的标记点利用γ的数值进行片内旋转计算，并与之前计算的tx，ty进行偏移更新计算。将相同角度图像上的标记点利用最小二乘进行优化。之后计算每个点的优化残差。寻找残差最小的角度作为γ。

(3)利用之前得到的γ进行更新平移变量tx，ty。

(4)将γ的搜索空间缩小为-5°到5°，以1°为以间隔再次利用(2)步骤进行寻找最优的γ值，并更新tx，ty的值，之后在此取值下，对ty再次利用最小二乘进行更新，并记录当前的残差值。

(5)迭代进行步骤(4)直到残差不在减小停止。停止后利用当前的tx，ty以及γ值对整个序列进行最小二乘优化，得到最终解。

在一个实施例中，所述赋权的迭代式捆绑调节算法，包括以下步骤：

(1)将相机参数与基准点的二维坐标作为输入；

(2)将标记点和相机参数进行分组，分组后将L₂范数作为目标函数对输入值进行优化；

(3)计算优化后的值的残差，利用残差计算每个基准点的归一化后的中位数的绝对偏差(MADN)，通过MADN的值为标记点进行分配权重；

(4)将带权的基准点再次利用L₂范数的目标函数进行迭代式的优化；

(5)当每个基准点的权重基本不再发生改变的时候优化停止，得到最终的相机参数与基准点的二维坐标。图2展示了权重迭代算法的思路，通过每次的优化使离群点的权重逐渐减小，从而避免离群点的影响。

在一个实施例中，所述捆绑调节目标是实现胶体金的估计值与测量值最小化：

其中表示估计所得到的空间三维胶体金坐标，Proj(·)为投影过程，x_i,j是测量得到的胶体金坐标δ_i,j说明第j个胶体金是否在第i张投影中可见。

在一个实施例中，所述步骤(3)中权重的计算具体如下：

d_i,j＝e_i,j/(k*MADN)

MADN＝median(|e_i,j-median(e)|)/0.6475

其中，e_i,j表示每个基准点的残差，k是标记权重为零的阈值，默认值为4.685。这个值通常用于稳健拟合，因为在误差呈正态分布的情况下，它的渐近效率可达95％。

实施例2

以实施例1为基础进行具体实验验证，验证数据集选取：Nitrosop3和Vibrio数据集。

在此展示了这两套数据集分别利用L₂和本发明赋权迭代算法优化后的残差，从表1可以看到本发明赋权迭代算法运行得到的残差更小，并且残差的方差更小，说明得到的投影参数更加稳定并且更精确。

表1

为了更加直观的展示对齐效果，图3展示了这两套数据集对齐前后基准点轨迹在空间中的分布状态。从图3可以看出，两个数据集经过本发明处理后，由离散程度高、杂乱的数据形态，最终得到了清晰排序的数据集，即最终完成了图像序列的对齐。

在上述实施例的基础上，本发明继续对其中涉及到的技术特征及该技术特征在本发明中所起到的功能、作用进行详细的描述，以帮助本领域的技术人员充分理解本发明的技术方案并且予以重现。

最后，虽然本说明书按照实施方式加以描述，但并非每个实施方式仅包含一个独立的技术方案，说明书的这种叙述方式仅仅是为清楚起见，本领域技术人员应当将说明书作为一个整体，各实施例中的技术方案也可以经适当组合，形成本领域技术人员可以理解的其他实施方式。

Claims (5)

1.一种冷冻电镜图像序列的精准快速对齐方法，其特征在于，该方法包括如下步骤：

S1：获取冷冻电镜图像数据，对投影图像序列上的胶体金进行识别；

S2：对低角度-50°～+50°图像上的胶体金进行匹配跟踪；

S3：梳理投影图像序列中的胶体金匹配对，生成贯穿整个投影图像序列的胶体金链，即二维投影坐标对应关系；

S4：基于正交投影模型进行投影参数的标定，首先进行初始化投影参数，之后通过捆绑调整算法求解得到初步的低角度下的胶体金的三维坐标预估计，通过三维坐标的预估计来进一步估算高角度投影图像上的胶体金虚拟坐标分布；

S5：通过与未加入的高角度投影图像相匹配，从低角度投影向高角度投影进行迭代式的序列跟踪，利用赋权的迭代式捆绑调节算法实现投影参数再标定，最终实现对整个序列的可用胶体金跟踪；所述赋权的迭代式捆绑调节算法，包括以下步骤：

(1)将相机参数与基准点的二维坐标作为输入；

(2)将标记点和相机参数进行分组，分组后将L₂范数作为目标函数对输入值进行优化；

(3)计算优化后的值的残差，利用残差计算每个基准点的归一化后的中位数的绝对偏差MADN，通过MADN的值为标记点进行分配权重；

(4)将带权的基准点再次利用L₂范数的目标函数进行迭代式的优化；

(5)当每个基准点的权重基本不再发生改变的时候优化停止，得到最终的相机参数与基准点的二维坐标；

S6：投影参数标定完成后，即得到优化后的投影参数后，将投影参数与未对齐的图像序列进行计算，得到对齐后的冷冻电镜图像序列，完成对齐。

2.如权利要求1所述的精准快速对齐方法，其特征在于，所述S4中，使用初始的正交投影模型，其数学建模如下：

s表示投影图像的尺度的变化，γ投影图像的片内旋转，α表示投影图像在成像过程中的绕转轴的倾转角度，β表示投影图像在成像过程中转轴的倾转角度，(u,v)表示测量所得的胶体金坐标，t表示平移变量；

3.如权利要求2所述的精准快速对齐方法，其特征在于，所述S4中，所述S4中提到的投影参数初始化，其计算流程如下：

对s，α，β，γ以及平移变量t进行初始化；因为图像为刚性变化，因此s为1，α是投影图像在成像过程中的绕转轴的倾转角度，接近为0，因此初始化为0；主要对平移变量tx，ty以及γ进行计算初始化计算；其计算步骤如下：

S4-1：把图像的中心确定，初始化度数最低图像的平移变量tx,ty为0；利用如下的仿射变换公式计算每两张图片之间的仿射变换关系：

之后从中心开始向左向右，分别初始化，从而计算每张图片的tx，ty；

S4-2：在得到初步的平移变量之后，对γ进行初始化，从-90°到90°寻找最佳的γ，以5°为一个间隔进行寻找；

以当前的角度进行平移变量的计算：

之后对根据不同角度图像上的标记点利用γ的数值进行片内旋转计算，并与之前计算的tx，ty进行偏移更新计算；将相同角度图像上的标记点利用最小二乘进行优化；之后计算每个点的优化残差；寻找残差最小的角度作为γ；

S4-3：利用之前得到的γ进行更新平移变量tx，ty；

S4-4：将γ的搜索空间缩小为-5°到5°，以1°为以间隔再次利用S4-2步骤进行寻找最优的γ值，并更新tx，ty的值，之后在此取值下，对ty再次利用最小二乘进行更新，并记录当前的残差值；

S4-5：迭代进行步骤S4-4直到残差不再减小停止，停止后利用当前的tx，ty以及γ值对整个序列进行最小二乘优化，得到最终解。

4.如权利要求1所述的精准快速对齐方法，其特征在于，所述S5中，所述目标函数是实现胶体金的估计值与测量值最小化：

其中表示估计所得到的空间三维胶体金坐标，Proj(·)为投影过程，x_i,j是测量得到的胶体金坐标，δ_i,j说明第j个胶体金是否在第i张投影中可见。

5.如权利要求4所述的精准快速对齐方法，其特征在于，所述S5中，所述步骤(3)中权重的计算具体如下：

d_i,j＝e_i,j/(k*MADN)

MADN＝median(|e_i,j-median(e_i,j)|)/0.6475

其中，e_i,j表示每个胶体金的残差，k是标记权重为零的阈值，w_i,j是每个胶体金的权重。