CN113658317B - 电子显微镜连拍图像处理方法和装置 - Google Patents

Abstract

本公开提供一种电子显微镜连拍图像处理方法和装置，该方法包括以下步骤：获取电子显微镜连拍系列图像；由连拍系列图像中的各整张图像按照时序叠加形成三维矩阵；将三维矩阵中的各个图像时域‑频域变换来变换到频域，得到三维频谱；将三维频谱经三维层叠滤波器进行滤波处理，得到滤波后的三维频谱；将滤波后的三维频谱进行所述时域‑频域变换的反变换，得到经滤波后的输出图像。本发明实施例的方法能够快速、有效地去除图像中的噪声。

Description

电子显微镜连拍图像处理方法和装置

技术领域

本发明图像处理技术，尤其涉及一种电子显微镜连拍图像处理方法和装置，特别是涉及一种透射电子显微镜连拍图像处理方法和装置。

背景技术

透射式电子显微镜(简称透射电镜)是在原子尺度研究材料结构的重要观测工具。透射电镜的高分辨成像模式下，以近似于平行的电子束照射到样品上，样品后表面的出射波经由电磁透镜的偏转后，可以获得放大的像，由位于像面的CCD、CMOS等照相机采集图像。近几年图像采集系统迅速发展，目前最新的Gatan K2、K3等系列高速相机能够以每秒百张的速度记录图像，能够快速记录下样品在电子辐照下的动态变化，其优异的可探测量子效率(DQE，detective quantum efficiency)和灵敏度，甚至让易受辐照的MOF样品(有机骨架样品)原子结构得以呈现。

连拍技术在拍摄易受辐射样品时尤为常见。例如石墨烯、单层氮化硼等单层样品，可以使用传统CCD相机拍摄多张图像后，叠加以提高信噪比。电镜硬件条件更好、成像稳定时，还可以得到原子结构随时间变化的图像。另一方面，某些样品不能接受大剂量的电子辐照，所以常以较少的电子剂量多次对样品拍照。极低剂量的拍照，需要使用Gatan K2、K3系列相机，以5％甚至更低的曝光率对样品进行连拍，可以对MOF样品、生物蛋白质等极易受电子辐照损失的样品进行结构研究。

连拍技术广泛应用于高分辨电镜成像使用中，由于研究的样品多为不耐电子辐射的样品，所以连拍后得到的系列图像通常需要校正漂移，因此，现有技术中对连拍图像常用的一种图像处理技术是对连拍得到的系列图像校正漂移之后直接叠加图像来提高图像的信噪比。但现有技术中，图像直接叠加后得到的信噪比的提升是比较有限的。例如，如果图像的信噪比为snr，则N张图像叠加后的图像信噪比为

并且，叠加后的图像无法反映各图像之间的变化。图像叠加后，图像与图像之间的差异被完全抹掉。

在不考虑调制传递函数(MTF，Modulation Transfer Function)的作用时，图像噪声污染和图像复原的一般过程还可以用如图1所示的过程来描述：原始图像f(x)被加性噪声n(x)污染后得到g(x)，g(x)通常就是记录的实验图像；将复原算子R作用到g(x)进行用于通过滤波去除噪声污染而得到复原图像的复原过程，即R(g)，得到复原图像f′(x)。由于复原处理是在频域进行的，因此复原算子包括对采集图像的时域-频域变换处理和滤波处理。其中，滤波处理可利用维纳滤波法对图像进行处理，维纳滤波是图像处理领域应用较广泛的方法，其对图像的处理在频域进行，维纳滤波器W(k)的一般性公式写作：

f′(x)＝F^-1{W(k)G(k)} (1)

其中，f′(x)表示滤波后得到的图像；G(k)为实验采集的图像g(x)的傅里叶变换形式；F^-1{}表示反傅里叶变换；|F(k)|²为期望图像f(x)的功率谱，|N(k)|²为噪声n(x)的功率谱。

尽管维纳滤波法会大大提高图像信噪比，但在图像处理过程中，会出现伪影，所谓伪影，即滤波像上出现一些不存在的假象。通过滤波方法提高图像信噪比，容易出现假象。当透射电镜研究的是具有周期性结构的晶体时，滤波产生的伪影可能会让晶格特征延伸到真空。而对每张图单独处理产生的伪影，将十分复杂。

现有技术中还有一种对常规图像进行滤波的方法，即3D块匹配算法(BM3D算法)，BM3D算法是2007年提出的滤波方法。该方法使用了块匹配的方方式，其是将图像分割成若干小的像素块，在当前帧和前后若干帧以参考块为中心的相似区域内挑选出与参考块最匹配的多个相近像素块，将相近像素块叠加三维像素数组，然后再对这个三维像素数组进行三维离散余弦变换、硬阈值滤波和三维离散余弦反变换，将反变换结果进行重构处理，即将滤波后的图像塞回原图位置，从而完成整张图像的滤波。图2所示为现有BM3D滤波过程的示意图，如图2所示，该BM3D算法的实现流程有两大步骤，第一步为基础估计，第二步为最终估计，而这两大步之中，又分别包括三个部分：相似块分组(如图2中的通过块匹配分组)，阈值协同滤波和聚合，其中，第一步中，相似块分组包括通过块匹配进行分组，并将匹配的块组合在一起构成三维矩阵；阈值协同滤波包括对三维矩阵进行从时域到频域的三维变换(如三维离散余弦变换)，变换后通过对三维矩阵设置阈值的方式将小于阈值的系数置0，以减少噪声，随后进行频域到时域的反变换(如三维离散余弦反变换)，并逐块估计，得到对图像块去噪之后的预计结果；聚合是将估计值返回初始位置，并通过对目标块和对应位置的块采用加权平均的方法进行处理，得到第一步基础估计的结果。第二步中的相似块分析是在结果图像和原始图像上分别进行逐块估计，得到两个三维矩阵；第二步中的协同滤波步骤中对两个三维矩阵进行三维变换，然后借助第一步中图像三维变换后的能量谱对原始图像组合得到的三维数组进行维纳滤波，滤波完之后进行数组反变换，最终得到对分解后图像的估计结果。第二步中的聚合过程与第一步中基本相同，不再详述。BM3D算法能够有效地减少图像上的噪声，提高图像的质量，因此BM3D算法被认为是目前除了深度学习之外滤波效果最好的滤波方法。但是，BM3D针对的都是常规图片的处理，且BM3D算法的块匹配计算量相当大，由于每次待处理的区域很小，一般是很小的二幂次平方(如8*8或16*16像素)，因此针对整个图像，BM3D的计算速度非常慢，计算过程复杂。BM3D算法也有针对视频的处理的，在视频上可以选择不同帧的相似区域进行块匹配，之后再进行三维层叠的滤波。因为视频是动态变化的，所以也是需要分割并选择相似区域，再匹配相似的其他帧区域进行块匹配。由于BM3D每次进行匹配的区域非常小，因此计算过程复杂，计算非常耗时。

若将BM3D方法应用于对电镜图像的去噪处理，则不仅会存在由于块匹配的区域很小而导致计算非常耗时，而且，由于电镜的晶格像具有周期重复性，还存在小区域块匹配时，所选择的不同帧的相似区域实际上未必对应着图像的正确区域的问题，即很容易出现块匹配错误。

如何针对电镜连拍得到的系列图像进行快速有效的去噪处理，是一个有待解决的问题。

发明内容

鉴于此，本发明实施例提供了一种电子显微镜连拍图像处理方法和装置，以消除现有技术中存在的一个或更多个缺陷。

本发明的技术方案如下：

一种电子显微镜连拍图像处理方法，该方法包括以下步骤：

获取电子显微镜连拍系列完整图像；

对按照拍摄时序排列的连拍系列完整图像进行从时域到频域的三维变换，得到三维频谱图像；

将所述三维频谱图像经三维层叠滤波器进行滤波处理，得到滤波后的三维频谱图像；

将滤波后的三维频谱图像进行从频域到时域的三维反变换，得到经滤波后的输出图像。

可选地，所述方法还包括：在形成三维矩阵之前对所述连拍系列图像进行图像漂移去除处理；或者在形成三维矩阵之前对基于各连拍系列图像的偏移量对各连拍系列图像进行平移处理。

可选地，所述方法还包括：在进行从时域到频域的三维变换之前对所述连拍系列图像进行图像漂移去除处理。

可选地，在确定不需要对按照拍摄时序排列的连拍系列完整图像进行平移处理的情况下，所述对按照拍摄时序排列的连拍系列完整图像进行从时域到频域的三维变换，得到三维频谱图像包括：将按照拍摄时序排列的连拍系列完整图像叠加成三维矩阵，将叠加后的三维矩阵进行从时域到频域的三维变换，得到三维频谱图像；

在确定要对按照拍摄时序排列的连拍系列完整图像进行平移处理的情况下，所述对按照拍摄时序排列的连拍系列完整图像进行从时域到频域的三维变换，得到三维频谱图像包括：对各完整图像的行和列进行从时域到频域的二维变换，得到二维频谱图像；将得到的二维频谱图像进行校正平移，将校正平移后的二维频谱图像按照拍摄时序进行叠加，形成三维矩阵；以及将三维矩阵进行第三维从时域到频域的变换，得到三维频谱图像。

可选地，所述时域-频域变换为三维傅里叶变换，所述频域-时域变换为三维反傅里叶变换；所述三维傅里叶变换包括：对系列图像的所有行和列做二维傅里叶变换，该二维傅里叶变换包括分别对系列图像的所有行和所有列做一维傅里叶变换，基于一维傅里叶变换结果得到二维傅里叶变换后的图像；将二维傅里叶变换后的图像沿着与行和列正交的第三维方向进行一维傅里叶变换，得到三维傅里叶变换结果；以及

所述三维反傅里叶变换包括：对经滤波后频谱图像的所有行和列做二维反傅里叶变换，该二维反傅里叶变换包括分别对经滤波后频谱图像的所有行和所有列做一维反傅里叶变换，基于一维反傅里叶变换结果得到二维反傅里叶变换后的图像；将二维反傅里叶变换后的图像沿着与行和列正交的第三维方向进行一维反傅里叶变换，得到三维反傅里叶变换结果。

可选地，所述三维傅里叶变换的函数满足以下公式：

其中，I(k_x，k_y，k_z)为图像的频域表示，I(x，y，z)为图像的时域表示；x，y，z表示时域的三个方向，k_x，k_y，k_z表示连拍图像在频域的三个方向；

所述三维反傅里叶变换的函数满足以下公式：

可选地，所述方法还包括：将三维滤波处理后，基于频谱图像在第三维方向的频率不为0的分量获得图像间差异。

可选地，所述滤波处理采用的滤波函数满足如下公式：

f_3D′(x)＝F^-1{W_3D(k)G_3D(k)}；

其中，f_3D′(x)表示滤波后得到的图像三维矩阵；W_3D(k)表示滤波器滤波函数；G_3D(k)为连拍系列图像的三维矩阵的傅里叶变换形式；F^-1{}表示反傅里叶变换；|G_3D(k)|²为期望图像的功率谱，σ²为噪声的功率谱，σ表示噪声功率；或者

所述滤波处理采用的滤波函数满足如下公式：

f_3D′(x)＝F^-1{W_3D(k)G_3D(k)}；

其中，f_3D′(x)表示滤波后得到的图像三维矩阵；W_3D(k)表示滤波器滤波函数；G_3D(k)为连拍系列图像的三维矩阵的傅里叶变换形式；F^-1{}表示反傅里叶变换；|G_3D(k)|²-σ²表示期望图像的功率谱，σ表示噪声功率。

可选地，所述方法还包括：将滤波后得到的图像数据集训练建立的深度学习模型，以利用训练的深度学习模型进行去噪处理；

所述深度学习模型为DnCNN网络模型；或者，

所述深度学习模型为DnCNN网络模型和GAN去噪网络模型。

可选地，所述电子显微镜为透射电子显微镜；所述电子显微镜连拍图像为透射电子显微镜连拍晶格图像。

本发明的另一方面，还提供一种电子显微镜连拍图像处理装置，该装置包括处理器和存储器，所述存储器中存储有计算机指令，所述处理器用于执行所述存储器中存储的计算机指令，当所述计算机指令被处理器执行时该装置实现如前所述方法的步骤。

本发明的另一方面，还提供一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时实现如前所述方法的步骤。

本发明的电子显微镜连拍图像处理方法和装置能够快速、有效去除噪声，且滤波像中伪影较少。并且，本发明能够以较高的信噪比复原连拍图像中每张图像的细节。本发明实施例还适用于低剂量甚至极低剂量、高噪声的连拍图像的快速处理。

本发明的附加优点、目的，以及特征将在下面的描述中将部分地加以阐述，且将对于本领域普通技术人员在研究下文后部分地变得明显，或者可以根据本发明的实践而获知。本发明的目的和其它优点可以通过在书面说明及其权利要求书以及附图中具体指出的结构实现到并获得。

本领域技术人员将会理解的是，能够用本发明实现的目的和优点不限于以上具体所述，并且根据以下详细说明将更清楚地理解本发明能够实现的上述和其他目的。

附图说明

此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，并不构成对本发明的限定。附图中的部件不是成比例绘制的，而只是为了示出本发明的原理。为了便于示出和描述本发明的一些部分，附图中对应部分可能被放大，即，相对于依据本发明实际制造的示例性装置中的其它部件可能变得更大。在附图中：

图1示出了现有技术中噪声污染及图像复原的过程示意图。

图2所示为现有BM3D滤波过程的示意图。

图3所示为本发明一实施例中透射电镜连拍图像的去噪方法流程示意图。

图4所示为本发明一实施例中三维图像层叠的示意图。

图5所示为本发明一实施例中去噪处理算法的示意图。

图6所示为本发明一实施例中的算法流程示意。

图7A至图7C分别为未经处理的电镜图像、采用本发明的3D层叠滤波方法和经传统维纳滤波方法处理后的电镜图像。

图8所示为本发明一实施例中沿着图7B和图7C中的剖线的像素强度示意图。

图9A至图9C分别为未经处理的电镜图像、采用本发明的3D层叠滤波方法和经传统维纳滤波方法处理的电镜图像处理后的电镜图像。

图10中的(a)、(b)、(c)分别示出了未经处理的电镜图像、采用本发明的3D层叠滤波方法处理后的电镜图像和经现有的BM3D方法处理后的电镜图像。

图11为本发明实施例中利用三维层叠滤波法生成的一些训练集的示意图。

图12为本发明实施例中将三维层叠滤波方法应用于深度学习网络模型的示意图。

图13为使用DnCNN模型来进行电镜图像去噪的示意图。

图14为DnCNN网络模型的去噪效果与3D层叠去噪效果的比较图。

图15为将GAN去噪模型中鉴别器的处理流程示意图。

图16为GAN去噪模型的去噪效果与3D层叠去噪效果的比较图。

图17示出了本发明的各种去噪方法与现有方法的效果比较示意图，(a)、(b)、(c)、(d)、(e)分别对应原始图像、传统维纳滤波去噪图像、三维层叠维纳滤波图像、DnCNN模型去噪图像、和GAN模型去噪图像。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，下面结合实施方式和附图，对本发明做进一步详细说明。在此，本发明的示意性实施方式及其说明用于解释本发明，但并不作为对本发明的限定。

在此，还需要说明的是，为了避免因不必要的细节而模糊了本发明，在附图中仅仅示出了与根据本发明的方案密切相关的结构和/或处理步骤，而省略了与本发明关系不大的其他细节。

应该强调，术语“包括/包含”在本文使用时指特征、要素、步骤或组件的存在，但并不排除一个或更多个其它特征、要素、步骤或组件的存在或附加。

对于现有的图像去噪处理要么去噪效果不够好，要么去噪过于耗时的问题，本公开针对透射电镜连拍图像，提出了一种电镜连拍图像的去噪方法，该方法在BM3D的基础上进行了改进，使得改进后的方法消除了现有技术中的一个或更多个缺陷，使得既提高了滤波速度，又提高了滤波效果。

图3所示为本发明一实施例中透射电镜连拍图像的去噪方法流程示意图。如图3所示，本方法包括如下步骤：

步骤S110，获取透射电镜连拍系列图像。

本步骤中，采用Gatan K2、K3等系列高速相机对样品进行连拍可获得透射电镜连拍系列图像，以快速记录下样品在电子辐照下的动态变化。例如可从连拍得到的实验像中获取N张实验像，g₁(x,y),g₂(x,y),…,g_N(x,y)，该N张实验像可以是从一次连拍中得到的所有系列图像中抽取的一部分，也可以是一次连拍得到的全部连拍图像。

由于通常存在样品漂移，因此连拍系列图像(实验像)之间往往存在图像漂移，在存在图像漂移的情况下，在获得连拍系列图像之后，本发明的去噪方法优选地还包括对获得的连拍系列图像进行消除漂移的预处理。互相关运算常用来计算两图之间的漂移，在具体实现时，可将两图的频谱相乘之后再做反傅里叶变换，就可以快速地得到互相关图。在频谱相乘时，可只取两图各自的相位部分，这样得到的互相关图在求解两图平移更为精确。此外，还可参照申请号为201710576072.4的中国专利申请，通过在各帧图像中确定特征点，并基于特征点来进行对齐方法进行图像漂移消除处理。这些方法或其他的现有消除漂移的方法都可用来尽量地消除系列图像的漂移。

由于消除图像漂移的这些方法均为现有技术，因此，在此不进行赘述。

步骤S120，对按照拍摄时序排列的连拍系列完整图像进行从时域到频域的三维变换，得到三维频谱图像。

由于是高速连拍，连拍系列电镜图像本身具有极大相似性，因此，在本发明实施例中，不需要BM3D这么麻烦的处理过程，即无需对图像分割成小块，可直接将连拍获得的系列图像中各个整张的图像按照拍摄时序进行层叠。

在本发明一个实施例中，在不存在图像平移或无需考虑图像平移的情况下，可以将得到的N张实验像，各实验像g_n(x，y)(n＝1，2，...，N)直接叠加成三维矩阵，其中，该N张实验像可以是经漂移消除处理之后的实验像，也可以是无需图像漂移消除处理而未经漂移消除处理的实验像。图4所示为多张连拍整图的三维层叠，X、Y方向分别为图像像素的行方向和列方向，Z方向为层叠的第三维，它可以认为是时间序列的排列。

由于图像的复原(去噪处理)是在频率域处理的，由此，本发明实施例中，在将N张连拍系列图像(实验像)按照拍摄时序叠加而得到三维矩阵之后，进一步对得到的三维矩阵做从时域到频域的3D变换，得到三维频谱图像G′(k_x，k_y，k_z)。通过3D变换将系列图像变换到频率域之后，便可以采滤波方法进行去噪处理。

在本发明一实施例中，3D变换采用的是三维傅里叶变换。三维傅里叶变换的具体可表示为：

其中，I(k_x，k_y，k_z)为三维矩阵的频域表示，I(x，y，z)为三维矩阵的时域表示；x，y，z表示连拍图像在时域的三个方向，k_x，k_y，k_z表示连拍图像在频域的三个方向。该三维傅里叶变化实际上是由两个一维傅里叶变换组成的二维傅里叶变换扩展而来，该二维傅里叶变换的计算过程是：对连拍图像的所有行做一维傅里叶变换，之后对连拍图像的所有列做一维傅里叶变换。然后基于二维傅里叶变换进一步扩展地理解到三维傅里叶变换：在得到所有各自图像的二维傅里叶变换后，沿第三个方向做一维傅里叶变换，如果图像的尺寸是U行*V行，那么第三个方向包含了U*V个一维傅里叶变换。

在本发明另一实施例中，在存在图像平移且需考虑图像平移的情况下，对于N张实验像中的每张实验像g_n(x，y)，如果图像存在平移，先将每张图做时域到频域的二维变换得到二维频谱图像G_n(k_x，k_y)，并乘以平移矩阵后得到平移后的二维频谱图像G′_n(k_x，k_y)，然后将平移后的二维频谱图像G′_n(k_x，k_y)按照拍摄时序在第三维方向上叠加成三维矩阵G′(k_x，k_y，z)后，对三维矩阵的第三维做从时域到频域的一维变换得到三维频谱图像G′(k_x，k_y，k_z)。通过时域到频域的变换将系列图像变换到频率域之后，便可以采滤波方法进行去噪处理。本实施例中，也是按照拍摄时序排列的连拍系列完整图像进行了从时域到频域的三维变换，也即，先进行了时域到频域的二维变换，平移后在第三维方向上进行了时域到频域的变换。作为示例，时域到频域的三维变换可以是傅里叶变换，即首先对实验像进行二维傅里叶变换(包括对行和列进行的两个一维傅里叶变换)，对变换后的二维频谱图像平移并叠加得到三维矩阵后，对第三维进行一维傅里叶变换，从而得到三维频谱图像。

傅里叶变换的本质可以理解为将信号分解为不同频率大小的三角函数exp(-ikx)，其中k为频率。当k＝0时，得到的信号恰为该维方向的所有数值的平均值；随着k的增加，三角函数的频率增大，所包含的图像细节也更多。

由此，运用三维傅里叶变换对图像进行变换后，得到频谱在第三维方向上，可以认为是所有图像在某个频谱(k_x，k_y)的傅里叶展开。由于z方向包括的图像总数并不多，以N作为图像总数，那么，本发明实施例把三维傅里叶变换后的频谱写作I(k_x，k_y，k_z＝-(N-1)/2N)...、I(k_x，k_y，k_z＝-2/N)、I(k_x，k_y，k_z＝-1/N)、I(k_x，k_y，k_z＝0)、I(k_x，k_y，k_z＝1/N)，I(k_x，k_y，k_z＝2/N)、...I(k_x，k_y，k_z＝(N-1)/2N)。频谱的中心分量I(k_x，k_y，k_z＝0)为所有图像的频谱的平均。k_z增大，恰对应着图像间频谱差异的高频分量。

基于上述三维傅里叶变换的分析，可以了解到，第三维的傅里叶变换，兼顾了系列图像的“同”与“异”。“同”是指I(k_x，k_y，k_z＝0)分量，是所有图像的频谱平均；“异”指的是I(k_x，k_y，k_z≠0)，可以分解系列图像的频谱差异，并随着k_z增大，图像间频谱差异的高频分量越大。

本发明实施例中，时域到频域的傅里叶变换仅为举例，但本发明并不限于此，还可以采用其他变换方法。

步骤S130，将三维频谱图像经三维层叠滤波器进行滤波处理，得到滤波后的三维频谱图像。

大部分的二维滤波方法都可以基于二维到三维的扩展方式推广到三维频谱空间进行滤波，且图像间的差异要保留多少可以通过控制I(k_x，k_y，k_z≠0)项的数目2(N-1)予以调整，也可以说通过控制z方向的图像总数来予以调整。

在本发明实施例中，可以采用多种滤波算法来实现对三维频谱图像的滤波处理。为了便于描述，本发明实施例中，将对连拍系列图像进行步骤S120的变换处理后进行的滤波算法称为三维层叠滤波算法，相应的滤波器可称为三维层叠滤波器。整个的去噪方法流程可称为三维层叠滤波去噪方法。

如图5所示，三维频谱图像G′(k_x，k_y，k_z)经滤波器滤波后得到G″(k_x，k_y，k_z)。

步骤S140，将三维层叠滤波器处理后的各个图像进行从频域到时域的三维反变换，得到估计的输出图像，即去噪后的复原图像。

作为示例，从频域到时域的三维反变换可以是三维反傅里叶变换。三维反傅里叶变换的具体可表示为：

二维反傅里叶变换的计算过程是，对三维频谱图像中各图像的所有行做一维反傅里叶变换，之后对各图像的所有列做一维反傅里叶变换，即对图像的行和列进行二维反傅里叶变换。在得到所有各自图像的二维反傅里叶变换后，在沿第三个方向(k_z方向)做一维反傅里叶变换，如果图像的尺寸是U行*V行，那么第三个方向包含了U*V个一维反傅里叶变换。

如图5所示，G″(k_x，k_y，k_z)做三维反傅里叶变换后，得到的即为滤波后的复原图像f′_n(x，y)。

将三维滤波处理后，基于频谱图像在第三维方向的频率不为0的分量获得的复原图像可作为体现出图像间的差异性的图像。

针对透射电镜图像，噪声的功率谱|N(k)|²可以用真空区域的图像进行估计，期望图像的功率谱|F(k)|²可以用实验像的功率谱|G(k)|²和噪声功率谱估计得到，即|F(k)|²＝|G(k)|²-|N(k)|²。进一步地，如果将|N(k)|²认为是均匀噪声，可以简化为σ²；当然，也可以对每张图像的噪声功率谱进行估计，而不做均匀噪声的近似。另一种更为简单的处理方法就是将|F(k)|²直接用|G(k)|²来估计。

因此，在本发明实施例中，三维层叠滤波器可采用多种实现方法，下面列举两种。

(1)使用第一种简单的3D层叠滤波方法来比较层叠后的滤波效果与传统方法的差别。

在该滤波方法中，假设噪声功率为σ，三维傅里叶变换后的频谱图像为G_3D(k)。期望像的噪声功率谱以G_3D(k)近似，滤波器可写作：

f_3D′(x)＝F^-1{W_3D(k)G_3D(k)} (5)

图7A是连拍图像中的一帧，图像大小为256*256。原图的信噪比为7.4。本示例中是提取一次连拍图像中的9张连拍图像对一帧图像进行分析，图7B给出本发明实施例采用的3D层叠滤波方法对该帧图像的滤波结果，图7B中左图和右图都是同一张，而右图给出了剖线用于进一步分析了真空区域是否有伪影出现，采用本发明的滤波器滤波后图像信噪比估计为75.25。图7C是采用传统维纳滤波方法得到的复原图像，信噪比为60.53。值得注意的是，滤波后传统维纳滤波方法的信噪比低于本发明中三维层叠滤波法的信噪比，但在真空区域，传统维纳滤波方法得到的复原图却出现了更严重的伪影，见图7C中箭头所指部分。

图8示出了沿着图7B和图7C中的剖线的像素强度示意图。如图8所示，箭头示意地方为真空，但维纳滤波法得到的图像出现了更明显的假象。

(2)使用第二种三维层叠滤波算法

假设噪声功率为σ，三维傅里叶变换的图像为G_3D(k)。期望像的噪声功率谱以|G_3D(k)|²-σ²近似，层叠滤波器可写作：

f_3D′(x)＝F^-1{W_3D(k)G_3D(k)} (7)

图9A为原始实验像，图9B为采用本发明的三维层叠滤波器得到的结果，图9B为普通维纳滤波方法得到的去噪图像的结果。采用本发明的三维层叠滤波的信噪比为96.65，采用传统维纳滤波方法的信噪比为82.08。滤波后特征应该尽可能的反映原图特征，虽然传统方法处理结果的原子特征更明显，但传统维纳滤波方法的伪影更明显，例如箭头所指的那排原子。

如上所述，在本发明实施例的滤波方法，通过输入若干张连拍图像(优选为经漂移处理的图像)，将图像层叠为三维矩阵g(x,y,z)后，做三维傅里叶变换，再将得到的三维频谱乘以特定的3D层叠滤波器后，做反傅里叶三维变换，即得到估计的输出图像f’(x,y,z)，如图6的简化示图所示意。在连拍图像需要平移处理的情况下，可先将各连拍图像进行二维傅里叶变换，变换后进行图像平移，并在平移后进行第三维傅里叶变换得到的三维频谱乘，然后将得到的三维频谱乘以特定的3D层叠滤波器后，做反傅里叶三维变换，得到估计的输出图像f’(x,y,z)。

本发明还给出了采用现有的BM3D方法与本发明的3D层叠滤波方法进行去噪处理后得到的图像的对比图，如图10所示。图10中的(a)行、(b)行、(c)行分别表示未经处理的3幅原始电镜图像、采用本发明的3D层叠滤波方法得到的去噪后的3幅电镜图像和经现有的BM3D方法处理后得到的3幅电镜图像。图10中可以清晰地看出，利用本发明的3D层叠滤波方法的去噪效果远远大于传统的BM3D去噪方法。

此外，若采用现有的BM3D的方法来对连拍图像进行去噪处理，对于9张图像的去噪可能需27秒左右的时间才能获得去噪图像，而采用本发明的方法，1.5秒左右的时间就能够获得同样9张图像的去噪效果很好的去噪图像。

本发明的如上电镜连拍图像处理方法(去噪方法)能够快速、有效去除噪声，实现透射电镜图像的快速处理，且滤波像中伪影较少。在相近的信噪比结果的情况下，本发明的去噪方法得到的伪影比传统方法小许多，结果更可信。并且，通过选取频谱图像在第三维方向的频率不为0的分量在滤波后进行反傅里叶变换，可以获得图像间差异，即能够以较高的信噪比复原连拍图像中每张图像的细节。本发明的方法还尤其适用于低剂量图像甚至极低剂量(如约5％左右及更低的曝光率)、高噪声的连拍图像的快速处理。经过测试，传统方法需要叠加数十张图像才能看清的细节，现在只需要几张就能恢复图像。

更进一步地，在本发实施例中，还将上述连拍系列图像的去噪结果用作深度学习模型的训练数据，以进一步利用深度学习模型进行连拍图像的去噪处理。

训练深度学习模型之前需要准备大量的训练数据，用以提高网络模型的泛化能力。对于本发明实施例的三维层叠滤波方法，可以通过增大参数σ来获得更高的图像质量，但是随着参数的增大，图像边缘的伪影会越来越明显，因此在本发明一些实施例中，可首先将原始图像进行三维层叠滤波，先去除一些噪声，然后进行深度学习模型训练，这样可以大大提高网络的学习效率，减轻网络的训练难度，然后将预处理图像与原始图像成对生成数据集。

在一示例中，选取了414张实验图像，即取46次9张系列连拍式电子显微像进行三维层叠滤波处理后生成的预处理图像作为深度学习网络的无噪输入，原始图像作为有噪输入，将二者整合为数据集，然后从中选取了400组图像作为训练集，5组图像作为验证集，9组图像作为测试集来训练深度学习网络模型，图11所示为生成的一些训练集，图11中，a行为原始电子显微像，b行为三维层叠滤波后的图像。

本发明一实施例中，深度学习网络模型可采用DnCNN模型。图12所示为基于本发明的三维层叠滤波方法对原始图像进行预处理，先去除一部分噪声，然后生成有噪和无噪的成对图像作为数据集，并用该数据集训练深度学习去噪模型DnCNN以进一步去除图像中的噪声的流程图。

DnCNN模型是于2017年提出的一种深度学习去噪算法，主要是基于残差学习(Residual Learning)结合批量归一化(BN，Batch Normalization)技术。由于深度学习网络的输入具有很大的数据量，在使用批量梯度下降法训练网络的情况下需要花费大量的时间来计算每个输入的梯度值，导致训练速度十分缓慢，因此本发明实施例中采用小批量梯度下降法来训练网络、学习图像中的特征。

其网络结构如图13所示，首先输入一张待处理图像，通过神经网络学习该图像与理想输出图像之间的残差图像即噪声图像，然后再通过将噪声图像从输入图像中去除来达到去噪的效果。

图13中，第1层采用卷积块Conv+ReLU，即先对输入图像以零填充(no padding)方式与多个(如64个)3x3卷积核执行卷积运算，然后使用线性整流单元(ReLU，RectifiedLinear Unit)函数进行激活运算得到下一层的输入图像。第2-16层采用卷积块Conv+BN+ReLU，即先将上一层的输出图像以零填充(no padding)方式与多个(如64个)3x3卷积核执行卷积运算后进行批量归一化处理，然后使用ReLU函数进行激活运算得到下一层的输入图像。第17层即最后一层只采用了卷积运算(Conv)，即将第16层的输出图像以零填充(nopadding)方式与多个(如64个)3x3卷积核执行卷积运算后得到通过残差学习而获得的噪声图像，然后再用第一层的输入图像减去噪声图像得到输出图像即去噪后的图像。

图14为DnCNN的去噪效果与3D层叠去噪效果的比较图，其中(a)行为原始实验图像，(b)行为采用本发明实施例的3D层叠去噪方法得到的去噪图像，(c)为采用DnCNN得到的去噪图像。用肉眼可以清晰的看出，DnCNN模型在测试集上的表现，在有原子的范围内，训练后的DnCNN模型的去噪效果仍然比较好，但伪影有增多。

为了进一步降低DnCNN模型中出现的伪影，本发明一实施例中，还提出了GAN去噪(GANDenoising)网络模型，通过将DnCNN模型用于生成器来生成一个无噪声的图像，同时与一个基于理想图像参考的鉴别器进行竞争，使生成的图像进一步得到改善，从而使生成的图像达到更好的去噪目的，且达到更逼真的效果。本发明实施例中，GAN去噪模型是基于生成对抗网络(GAN)的图像去噪GAN(ImageDenoisingGAN)模型改进而来的。

在生成器方面，本发明实施例使用的是上述的DnCNN网络结构，具体内容不再赘述。在鉴别器方面，GANDenoising模型使用了六个具有批量归一化块和泄露线性整流单元(LReLU，Leaky Rectified Linear Unit)激活函数的卷积网络，其网络结构如图15所示，第1-6层，采用卷积块Conv+BN+LReLU，即将上一层的输出图像与3x3卷积核执行卷积运算后进行批量归一化处理，然后使用LReLU函数进行激活运算得到下一层的输入图像。第7层，最后采用Conv+BN和一个Sigmoid函数将输出映射到规范化为[0,1]的概率分数，表示的是图像是“真实的”的概率。

图16为采用GAN去噪模型的去噪效果与3D层叠去噪效果的比较图，其中(a)行为原始实验图像，(b)行为采用本发明实施例的3D层叠去噪方法得到的去噪图像，(c)为采用GAN去噪模型得到的去噪图像。用肉眼可以清晰的看出，GAN去噪模型在测试集上的表现，在有原子的范围内，训练后的GAN去噪模型的去噪效果仍然比较好，且伪影减少。

图17示出了本发明的各种去噪方法与现有方法的效果比较示意图，(a)、(b)、(c)、(d)、(e)分别对应原始图像、传统维纳滤波去噪图像、三维层叠维纳滤波图像、DnCNN模型去噪图像、和GAN模型去噪图像。图17可以看出，三维层叠维纳滤波去噪图像具有最优的去噪效果，并且数据处理速度与传统BM3D相比快了数十倍。DnCNN模型去噪图像、和GAN模型去噪图像在保持了好的去噪效果的同时，可以实现更加快速地滤波，只不过防止伪影的效果略差。

如上描述的电子显微镜连拍图像去噪算法，不仅适用于透射电子显微镜连拍图像的去噪处理，还适用于扫描电子显微镜等图像的去噪处理。

与本发明的方法相应地，本发明实施例还提供了一种能够执行电子显微镜连拍图像处理的计算机设备，该计算机设备为通用数据处理装置，其包括通用的计算机硬件结构，其至少包括处理器和存储器。存储器适于存储处理器可执行的一条或多条指令或程序。该一条或多条指令或程序被处理器执行以实现上述电子显微镜连拍图像处理方法中的步骤。

本公开还涉及存储介质，其上可以存储有计算机程序代码，当程序代码被执行时可以实现本发明的方法的各种实施例，该存储介质可以是有形存储介质，诸如光盘、U盘、软盘、硬盘等。

需要明确的是，本发明并不局限于上文所描述并在图中示出的特定配置和处理。为了简明起见，这里省略了对已知方法的详细描述。在上述实施例中，描述和示出了若干具体的步骤作为示例。但是，本发明的方法过程并不限于所描述和示出的具体步骤，本领域的技术人员可以在领会本发明的精神后，作出各种改变、修改和添加，或者改变步骤之间的顺序。

本领域普通技术人员应该可以明白，结合本文中所公开的实施方式描述的各示例性的组成部分、系统和方法，能够以硬件、软件或者二者的结合来实现。具体究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本发明的范围。当以硬件方式实现时，其可以例如是电子电路、专用集成电路(ASIC)、适当的固件、插件、功能卡等等。当以软件方式实现时，本发明的元素是被用于执行所需任务的程序或者代码段。程序或者代码段可以存储在机器可读介质中，或者通过载波中携带的数据信号在传输介质或者通信链路上传送。“机器可读介质”可以包括能够存储或传输信息的任何介质。机器可读介质的例子包括电子电路、半导体存储器设备、ROM、闪存、可擦除ROM(EROM)、软盘、CD-ROM、光盘、硬盘、光纤介质、射频(RF)链路，等等。代码段可以经由诸如因特网、内联网等的计算机网络被下载。

还需要说明的是，本发明中提及的示例性实施例，基于一系列的步骤或者装置描述一些方法或系统。但是，本发明不局限于上述步骤的顺序，也就是说，可以按照实施例中提及的顺序执行步骤，也可以不同于实施例中的顺序，或者若干步骤同时执行。

本发明中，针对一个实施方式描述和/或例示的特征，可以在一个或更多个其它实施方式中以相同方式或以类似方式使用，和/或与其他实施方式的特征相结合或代替其他实施方式的特征。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明实施例可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (9)

1.一种电子显微镜连拍图像处理方法，其特征在于，该方法包括以下步骤：

获取电子显微镜连拍系列完整图像；

对按照拍摄时序排列的连拍系列完整图像进行从时域到频域的三维变换，得到三维频谱图像；

将所述三维频谱图像经三维层叠滤波器进行滤波处理，得到滤波后的三维频谱图像；

将滤波后的三维频谱图像进行从频域到时域的三维反变换，得到经滤波后输出的复原图像；

所述三维变换和二维变换分别为三维傅里叶变换和二维傅里叶变换，三维反变换为三维反傅里叶变换；

所述三维傅里叶变换包括：

对每张图像的所有行和列做二维傅里叶变换，该二维傅里叶变换包括分别对系列图像的所有行和所有列做一维傅里叶变换，基于一维傅里叶变换结果得到二维傅里叶变换后的图像；

将二维傅里叶变换后的图像沿着与行和列正交的第三维方向进行一维傅里叶变换，得到三维傅里叶变换结果；以及

所述三维反傅里叶变换包括：

对经滤波后频谱图像的所有行和列做二维反傅里叶变换，该二维反傅里叶变换包括分别对经滤波后频谱图像的所有行和所有列做一维反傅里叶变换，基于一维反傅里叶变换结果得到二维反傅里叶变换后的图像；

将二维反傅里叶变换后的图像沿着与行和列正交的第三维方向进行一维反傅里叶变换，得到三维反傅里叶变换结果。

2.根据权利要求1所述的方法，其特征在于，所述方法还包括：在进行从时域到频域的三维变换之前对所述连拍系列完整图像进行图像漂移去除处理。

3.根据权利要求1或2所述的方法，其特征在于，

在确定不需要对按照拍摄时序排列的连拍系列完整图像进行平移处理的情况下，所述对按照拍摄时序排列的连拍系列完整图像进行从时域到频域的三维变换，得到三维频谱图像包括：将按照拍摄时序排列的连拍系列完整图像叠加成三维矩阵，将叠加后的三维矩阵进行从时域到频域的三维变换，得到三维频谱图像；

在确定要对按照拍摄时序排列的连拍系列完整图像进行平移处理的情况下，所述对按照拍摄时序排列的连拍系列完整图像进行从时域到频域的三维变换，得到三维频谱图像包括：对各完整图像的行和列进行从时域到频域的二维变换，得到二维频谱图像；将得到的二维频谱图像进行校正平移，将校正平移后的二维频谱图像按照拍摄时序进行叠加，形成三维矩阵；以及将三维矩阵进行第三维从时域到频域的变换，得到三维频谱图像。

4.根据权利要求1所述的方法，其特征在于，

所述三维傅里叶变换的函数满足以下公式：

其中，I(k_x，k_y，k_z)为图像的频域表示，I(x，y，z)为图像的时域表示；x，y，z表示时域的三个方向，k_x，k_y，k_z表示连拍图像在频域的三个方向；

所述三维反傅里叶变换的函数满足以下公式：

5.根据权利要求1所述的方法，其特征在于，所述方法还包括：

将三维滤波处理后，基于频谱图像在第三维方向的频率不为0的分量获得的复原图像作为体现差异性的图像。

6.根据权利要求1所述的方法，其特征在于，

所述滤波处理采用的滤波函数满足如下公式：

f_3D′(x)＝F^-1{W_3D(k)G_3D(k)}；

其中，f_3D′(x)表示滤波后得到的图像三维矩阵；W_3D(k)表示滤波器滤波函数；G_3D(k)为连拍系列图像的三维矩阵的傅里叶变换形式；F^-1{}表示反傅里叶变换；|G_3D(k)|²为期望图像的功率谱，σ²为噪声的功率谱，σ表示噪声功率；或者

所述滤波处理采用的滤波函数满足如下公式：

f_3D′(x)＝F^-1{W_3D(k)G_3D(k)}；

其中，f_3D′(x)表示滤波后得到的图像三维矩阵；W_3D(k)表示滤波器滤波函数；G_3D(k)为连拍系列图像的三维矩阵的傅里叶变换形式；F^-1{}表示反傅里叶变换；|G_3D(k)|²-σ²表示期望图像的功率谱，σ表示噪声功率。

7.根据权利要求1所述的方法，其特征在于，所述电子显微镜为透射电子显微镜；所述电子显微镜连拍图像为透射电子显微镜连拍晶格图像；

所述方法还包括：将滤波后得到的图像数据集训练建立的深度学习模型，以利用训练的深度学习模型进行去噪处理；

所述深度学习模型为DnCNN网络模型；或者，

所述深度学习模型为DnCNN网络模型和GAN去噪网络模型。

8.一种电子显微镜连拍图像处理装置，包括处理器和存储器，其特征在于，所述存储器中存储有计算机指令，所述处理器用于执行所述存储器中存储的计算机指令，当所述计算机指令被处理器执行时该装置实现如权利要求1至7中任意一项所述方法的步骤。

9.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，该程序被处理器执行时实现如权利要求1至7中任意一项所述方法的步骤。

Images