一种基于图像金字塔与NLM的冷冻电镜图像增强方法

一种基于图像金字塔与nlm的冷冻电镜图像增强方法
技术领域
1.本发明属于结构生物学分析技术领域，具体涉及一种基于图像金字塔与 nlm的冷冻电镜图像增强方法。


背景技术：

2.冷冻电子显微镜(cryogenic electron microscopy，cryo-em)技术是目前测定生物大分子结构的强有力手段及主流方法之一。使用冷冻电镜技术重构生物分子三维结构的总体步骤大致可分为4步：样本纯化与制备、数据采集与分析、二维颗粒图像挑选与分类、三维结构重构。由于生物样本的制备技术和电子显微镜硬件系统的限制，使得冷冻电镜采集到的图像信噪比极低，且存在部分非样本颗粒干扰，导致自动挑选冷冻电镜图像中的样本颗粒较为困难。
3.目前常见的颗粒挑选方法包括模板匹配法、图像分割法、深度学习的方法等。模板匹配法通过计算模板与图像待检测窗口之间的相关分数，来判断区域是否为颗粒或背景噪声，其中模板可以是人工挑选或者模拟的颗粒和噪声。图像分割法一般使用图像分割技术从经过图像增强处理的冷冻电镜图像中分割背景与颗粒。由于近年来深度学习的卓越表现，大量深度学习方法被用于颗粒挑选，深度学习方法通过训练深度学习网络或模型，使其具有区分背景与颗粒的能力，实现颗粒挑选。
4.由于冷冻电镜图像信噪比极低的问题，导致图像特征难以提取，除了图像分割法对图像有预处理外，其他方法一般没有或者只有简单的预处理，图像中颗粒区域特征没有得到足够的增强，使得主流的颗粒挑选方法流程非常复杂，适应场景有限，颗粒挑选的整体精确度不够高。


技术实现要素：

5.本发明针对现有技术的不足，提供一种基于图像金字塔与nlm的冷冻电镜图像增强方法，用于增强电镜图像的质量，降低噪声对颗粒区域的影响，使得边缘更加明晰，从而使后续颗粒挑选流程可以简化，并且可以提高挑选的整体效率。
6.为了达到上述目的，本发明提供的技术方案是一种基于图像金字塔与nlm 的冷冻电镜图像增强方法，包括以下步骤：
7.步骤1，输入冷冻电镜图像数据；
8.步骤2，对冷冻电镜图像使用非局部均值去噪算法去除图像噪声；
9.步骤3，对去除噪声的图像使用基于离散余弦变换的背景矫正算法，消减图像亮度不均现象；
10.步骤4，设定金字塔图像最小尺寸minsize，复制背景矫正后的图像构建图像金字塔；
11.步骤5，判断背景矫正后图像尺寸是否小于设定值minsize，若小于，则图像金字塔构建完成，并执行步骤6，否则对背景矫正后图像使用下采样技术，重新执行步骤2-步骤4；
12.步骤6，对图像金字塔每一层使用图像直方图变换和上采样技术，使金字塔各级图像尺寸相同；
13.步骤7，将图像金字塔各级图像形成的相同尺寸图像累加取均值，获取融合图像。
14.而且，所述步骤2中非局部均值去噪算法计算公式如下：
[0015][0016]
式中，v表示含噪声的图像，x为噪声图像v中的像素点，ω
x
为像素点x的搜索窗口，w(x,y)为邻域内像素点x与像素点y之间的相似度权重，为去噪后的图像中像素点x的像素值。
[0017]
w(x,y)的计算方式如下：
[0018][0019]
式中，d2为两个相似度计算块的欧氏距离，σ为噪声的标准差，h为一个基于σ的滤波平滑参数。
[0020]
d2的计算方式如下：
[0021]
d2＝||v(x)-v(y)||2ꢀꢀ
(3)
[0022]
式中，v(x)，v(y)分别为像素点x与像素点y的相似度计算块。
[0023]
而且，所述步骤3中对于尺寸为m
×
n图像f(x,y)，二维离散余弦变换 dct的过程如下：
[0024][0025][0026][0027]
式中，(x,y)、(u,v)均为像素坐标，c(u,v)为离散余弦变换后的dct系数矩阵。
[0028]
对系数矩阵c(u,v)使用一个低通滤波器，再经过离散余弦变换的逆变换过程得到图像的近似背景。
[0029]
低通滤波器公式如下：
[0030][0031]
式中，u
th
,v
th
分别为u,v方向的阈值，c
lp
(u,v)为低通滤波后的dct系数矩阵。
[0032]
使用二维离散余弦变换的逆变换idct获得图像背景的过程如下：
[0033]
[0034][0035][0036]
式中，(x,y)、(u,v)均为像素坐标，fb(x,y)表示近似的背景图像， c
lp
(u,v)为低通滤波后的dct系数矩阵。
[0037]
根据原始图像f(x,y)与背景图像fb(x,y)进行图像矫正的公式如下：
[0038]fsub
(x,y)＝f(x,y)-fb(x,y)
ꢀꢀ
(11)
[0039][0040]
式中，f
sub
(x,y)为原始图像减去背景后获得的矫正图像，f
norm
(x,y)是对矫正图像进行归一化处理后的图像。
[0041]
而且，所述步骤6中直方图变换公式如下：
[0042][0043]
式中，a为直方图变换后亮度的最高值，b为变换前背景像素分布的均值， k为控制变换函数图像的坡度的系数，f(x,y)表示直方图变换前的金字塔图像， f
ht
(x,y)表示直方图变换后的金字塔图像。
[0044]
而且，所述步骤步骤7中图像最终融合方式如下：
[0045][0046]
式中，upsampling(ik)表示对图像ik的上采样，k为金字塔图像的层数。
[0047]
与现有技术相比，本发明具有如下优点：
[0048]
1)使用了非局部均值去噪的方法，能够更有效地抑制图像中噪声，减小噪声对后续图像增强步骤和颗粒挑选步骤的负面影响；
[0049]
2)使用了图像金字塔、图像背景矫正与图像直方图变换的方法，使得增强图像中颗粒的特征更为明显，颗粒与背景之间的区分度更高，更有利于对颗粒图像的提取与分割。
附图说明
[0050]
图1为本发明实施例的技术流程图。
[0051]
图2为本发明与其它两种图像增强方法的结果对比图，其中图2(a)为模拟电镜图像与三种图像增强方法得到的增强图，图2(b)为图2(a)的二值图，图2(c)为真实电镜图像与三种图像增强方法得到的增强图，图2(d)为图2 (c)的二值图。
[0052]
图3为带噪图像(实线)与增强图像(虚线)的二值图像psnr与ssim 结果对比图，其中图3(a)为二值图像psnr的对比图，图3(b)为二值图像 ssim的对比图。
具体实施方式
[0053]
本发明提供一种基于图像金字塔与nlm的冷冻电镜图像增强方法，下面结合附图和实施例对本发明的技术方案作进一步说明。
[0054]
如图1所示，本发明实施例的流程包括以下步骤：
[0055]
步骤1，输入冷冻电镜图像数据。
[0056]
步骤2，对冷冻电镜图像使用非局部均值(nlm)去噪算法去除图像噪声。
[0057]
非局部均值去噪算法计算公式如下：
[0058][0059]
式中，v表示含噪声的图像，x为噪声图像v中的像素点，ω
x
为像素点x的搜索窗口，w(x,y)为邻域内像素点x与像素点y之间的相似度权重，为去噪后的图像中像素点x的像素值。
[0060]
w(x,y)的计算方式如下：
[0061][0062]
式中，d2为两个相似度计算块的欧氏距离，σ为噪声的标准差，h为一个基于σ的滤波平滑参数。
[0063]
d2的计算方式如下：
[0064]
d2＝||v(x)-v(y)||2ꢀꢀ
(3)
[0065]
式中，v(x)，v(y)分别为像素点x与像素点y的相似度计算块(相似块)。
[0066]
本实施例中，除了使用nlm去噪算法，还可以使用高斯滤波、中值滤波、维纳滤波、导向滤波、双边滤波、k-svd、bm3d、knnm等常规去噪方法，也可以使用dncnn、noise2noise等基于深度学习的去噪方法。考虑算法模型的复杂度、算法的计算时间与去噪效果等，使用nlm去噪算法能够以较高的效率达到较好的去噪效果，在nlm算法中选用较小的邻域可以减少大量计算、提高图像处理速度，本实施例中每个像素点的搜索窗口尺寸为21
×
21，相似块图像尺寸选定为7
×
7像素。
[0067]
步骤3，对去除噪声的图像使用基于离散余弦变换的背景矫正算法，消减图像亮度不均现象。
[0068]
对于尺寸为m
×
n图像f(x,y)，二维离散余弦变换(dct)的过程如下：
[0069][0070][0071][0072]
式中，(x,y)、(u,v)均为像素坐标，c(u,v)为离散余弦变换后的dct系数矩阵。
[0073]
在dct的系数矩阵c(u,v)中，低频部分为图像能量集中的部分，图像比较平坦，一般为背景；高频部分为图像变化比较剧烈的区域，一般为边界。因此，对系数矩阵c(u,v)使用一个低通滤波器，再经过离散余弦变换的逆变换过程即可得到图像的近似背景。
[0074]
低通滤波器公式如下：
[0075][0076]
式中，u
th
,v
th
分别为u,v方向的阈值，c
lp
(u,v)为低通滤波后的dct系数矩阵。
[0077]
通过对去噪后冷冻电镜图像亮度不均情况的分析，发现整体背景亮度不均区域分布在对角线上的两个角落区域，或者图像中间区域，尺寸大概占据图像的四分之一到二分之一的区域。在使用dct做背景矫正的实验中，也发现两个阈值u
th
,v
th
均取3时，可以兼顾多种类型的亮度不均图像，且可以取得满足后续处理的效果。
[0078]
使用二维离散余弦变换的逆变换(inversediscretecosinetransform,idct) 获得图像背景的过程如下：
[0079][0080][0081][0082]
式中，(x,y)、(u,v)均为像素坐标，fb(x,y)表示近似的背景图像， c
lp
(u,v)为低通滤波后的dct系数矩阵。
[0083]
根据原始图像f(x,y)与背景图像fb(x,y)进行图像矫正的公式如下：
[0084]fsub
(x,y)＝f(x,y)-fb(x,y)
ꢀꢀ
(11)
[0085][0086]
式中，f
sub
(x,y)为原始图像减去背景后获得的矫正图像，f
norm
(x,y)是对矫正图像进行归一化处理后的图像。
[0087]
本实施例中，对图像背景的矫正方法，除了使用基于dct变换获取图像背景外，还可以使用图像高斯模糊、图像顶帽底帽变换、傅里叶变换等方式来获得图像背景，但存在矫正效果较差或者计算时间较长等问题，使用dct变换可以达到实时处理的速度，也可以取得较好的矫正效果。
[0088]
步骤4，设定金字塔图像最小尺寸minsize，复制背景矫正后的图像构建图像金字塔。
[0089]
通过对原尺寸为4096
×
4096像素的冷冻电镜图像进行分析发现：对小尺寸颗粒(80
×
80像素及以下)，当图像尺寸缩小至200
×
200像素时，大部分颗粒基本不可与背景区分，为了兼顾对不同尺寸颗粒的处理，本实施例中金字塔图像最小尺寸设定为200
×
200像素。
[0090]
步骤5，判断背景矫正后图像尺寸是否小于设定值minsize，若小于，则图像金字塔构建完成，并执行步骤6；否则对背景矫正后图像使用下采样技术，重新执行步骤2-步骤4。
[0091]
步骤6，对图像金字塔每一层使用图像直方图变换和上采样技术，使金字塔各级图像尺寸相同。
[0092]
直方图变换公式如下：
[0093][0094]
式中，a为直方图变换后亮度的最高值，b为变换前背景像素分布的均值，k为控制变换函数图像的坡度的系数，f(x,y)表示直方图变换前的金字塔图像， f
ht
(x,y)表示直方图变换后的金字塔图像。
[0095]
对冷冻电镜图像经过图像非局部均值去噪和背景矫正后的图像进行分析，发现背景中仍然存在部分区域亮度过高，或者部分背景区域与颗粒区域区分度较差的情况。本实施例经过多次实验尝试，发现将a设置为230，b设置为110， k设置为20时，可以使亮度过高区域最高亮度降至230，避免后续处理可能出现的过曝问题，整体背景亮度得到不同程度的提高，使颗粒区域的亮度降低。对比直方图变换前后的图像直方图可发现，处理后颗粒区域像素集中在50以下，背景像素集中在200以上，增加了颗粒与背景之间的区分度。
[0096]
步骤7，将图像金字塔各级图像形成的相同尺寸图像累加取均值，获取融合图像。
[0097]
对金字塔图像的每一层进行上采样，上采样结果图像与原始图像具有相同的尺寸，图像最终融合方式如下：
[0098][0099]
式中，ii表示金字塔图像的第i层，upsampling(ii)表示对金字塔图像的第i 层进行上采样，k为金字塔图像的总层数。
[0100]
下面通过对比实验来对本发明的效果做进一步的阐述。
[0101]
一、仿真条件：
[0102]
实验图像数据有两类，一种是利用生物大分子三维结构投影得到的模拟电镜图像，一种是冷冻电镜拍摄的真实电镜图像。本文使用的三维结构文件为生物大分子的断层扫描图像，所有的数据文件都来自于存放冷冻电镜图像与断层扫描图像的公开存储库emdb，使用的四个代表性断层扫描图像文件编号为emd-2660、emd-2824、emd-3137、emd-3645。使用的真实冷冻电镜图像数据集编号为empiar-10028、empiar-10017、empiar-10033、empiar
‑ꢀ
10089，断层扫描图像文件是由上述真实电镜像数据集经过颗粒挑选和三维重构形成的，真实电镜编号顺序与上述断层扫描图像文件编号顺序对应。实验中，所有方法使用的语言都由python实现，两个对比方法(autocryopicker, supercryoempicker中的图像增强方法)的运行的软件为matlab。
[0103]
二、实验内容和结果
[0104]
分别将电镜图像数据集原始电镜图像作为本发明方法与两种对比方法的输入数据。在本发明方法处理流程中，对电镜图像使用非局部均值去噪的方法去除图像噪声，接着使用离散余弦变换矫正图像背景，然后对图像下采样获得下级图像，对下级图像再次使用
上述的去噪、背景矫正和下采样处理获得下级图像，循环此过程直至图像尺寸达到设定值。上述循环过程中得到的尺度逐次减小的背景矫正后图像，形成图像金字塔中分辨率逐步降低的每一级。对金字塔每一级图像进行直方图变换和上采样后获得相同尺寸图像，最后将上采样结果图像融合，即获得最终的增强图像。对比方法处理流程中，原始电镜图像作为输入数据，经过对比方法中的一系列图像处理步骤，得到最终增强图像。
[0105]
对模拟电镜图像和真实电镜图像经过三种图像增强方法处理，从主观和客观两方面进行评价，主观上对比三种方法的得到的增强图像和二值图像，客观上使用峰值信噪比psnr(peaksignal-to-noiseratio)和结构相似ssim (structuralsimilarity)来衡量模拟电镜图像经过图像增强得到的图像二值图与干净无噪声图像二值图之间的相似度。
[0106]
本发明与对比方法对图像增强结果如图2所示，对模拟图像增强结果二值图计算psnr和ssim的实验数据结果如表1所示：
[0107]
表1三种图像增强方法得到的psnr和ssim
[0108][0109]
由表1可知，本发明方法得到的图像平均峰值信噪比psnr为13.1687，比autocryopicker中的增强方法8.8393提高约4.3db，比supercryoempicker中的增强方法11.4388提高约1.7db。从平均结构相似ssim上来看，本方法的 ssim为0.8576，比autocryopicker中的增强方法0.6729提高约0.1847，比 supercryoempicker中的增强方法0.8428提高约0.0148。由结果对比可知，本发明使用的基于图像金字塔与nlm的冷冻电镜图像增强方法，能在冷冻电镜图像增强阶段中获得更优的结果。由图2可知，从主观上对比增强结果图像，本发明的方法能够得到更好的图像增强效果。
[0110]
对二值图像计算的psnr与ssim结果对比如图3所示，实线为处理前图像的结果，虚线为处理后图像的结果。由图3可知，在噪声比例较大(0.5～0.9) 时，处理后的图像psnr和ssim均有明显地提升。由于使用了图像金字塔，在多次下采样后图像丢失部分细节结构，在噪声比例较小(0.1～0.3)时，处理后图像的psnr和ssim结果比处理前结果更低或相近。但真实冷冻电镜图像的特点就是信噪比极低，低噪声比例的情况并不影响本发明方法对真实冷冻电镜图像的处理。
[0111]
具体实施时，以上流程可采用计算机软件技术实现自动运行流程。
[0112]
本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替
代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。