CN107798663B - 一种基于偏微分方程和bm3d的无参数化图像恢复方法 - Google Patents

Abstract

本发明提出了一种基于偏微分方程和BM3D的无参数化图像恢复方法，属于图像处理技术领域。传统的BM3D改进方法是直接利用现有的噪声参数估计方法得到一个参数估计，再直接带入原始BM3D方法中，但结果不稳定，本发明所述图像恢复方法利用偏微分方程对传统BM3D方法的去噪过程和噪声方差估计过程进行了改进，特点是跟偏微分方程预处理紧密融合，能使本发明所述的改进BM3D方法具有高精度和高稳定性。本发明所述一种基于偏微分方程和BM3D的无参数化图像恢复方法能够应用各种图像处理领域。

Description

一种基于偏微分方程和BM3D的无参数化图像恢复方法

技术领域

本发明涉及一种基于偏微分方程和BM3D的无参数化图像恢复方法，属于图像处理技术领域。

背景技术

BM3D(Block Matching and 3D Filtering，图像块匹配算法)方法是图像去噪领域中的一个优秀算法，其结果具有PSNR值高，视觉效果好，图像细节信息保持完整等优点。但其也有明显缺点，主要体现在BM3D方法的处理效果对于噪声图像中的噪声水平参数敏感度极高。具体体现在以下几个方面：

1.当输入的噪声标准差在40附近时，BM3D处理结果的PSNR(峰值信噪比)值会产生突变；

2.当输入的噪声标准差非整数时，BM3D处理结果的PSNR值会产生剧烈震荡；

3.准确的噪声标准差并非BM3D方法的最佳输入参数。并且由于噪声的随机性，相同噪声水平下的两幅不同噪声图像经过BM3D方法处理后，最佳的输入参数值也可能有所不同。

BM3D方法需要同时输入噪声图像和噪声水平，才能有效进行去噪。对于实际应用问题，图像的真实噪声水平是未知的，首先需要对噪声水平进行估计，再将估计值输入到BM3D算法中，完成后续处理。但是由于上述1缺陷，若真实噪声标准差在40附近，那么输入估计的噪声参数可能会严重降低恢复质量。由于上述2缺陷，噪声参数估计的结果普遍为非整数，直接带入BM3D算法也会对图像恢复结果产生影响。由于上述3缺陷，即使利用一些高精度的噪声水平估计方法得到了精确的噪声水平值，直接带入BM3D方法中也不一定得到最优的处理结果。

发明内容

本发明为了解决现有BM3D方法存在的上述技术问题，提出了一种基于偏微分方程和BM3D的无参数化图像恢复方法，所采取的技术方案如下：

一种基于偏微分方程和BM3D的无参数化图像恢复方法，所述图像恢复方法：首先，利用偏微分方程对噪声图像进行去噪预处理获得基本估计图，然后将所述基本估计图进行分块而获得基本估计图块，并利用距离匹配方法将所述基本估计图块匹配为基本估计图组；

其次，将所述基本估计图进行图像边界探测与去除处理后获得平滑无边界图像，通过灰度分割将所述平滑无边界图像分割出单灰度平滑区域，依次计算所述单灰度平滑区域中每个区域的区域样本方差，对获得的所述区域样本方差进行加权平均处理后获得噪声方差估计值；

在此，将噪声图像进行分块并获得噪声图块，将噪声图块按基本估计图块的距离匹配结果进行分组进而获得噪声图组；

最后通过协同维纳滤波方法，利用所述基本估计图组和所述噪声方差估计对所述噪声图组进行处理并获得图组估计，将所述图组估计以所述噪声方差估计为权重参数进行加权合成处理后获得最终估计图像；所述最终估计图像即为恢复后图像。

进一步地，所述图像修复的具体过程为：

步骤一：输入待处理的噪声图像f；

步骤二：利用基于各项异性扩散方程的图像去噪方法对噪声图像f做预处理，并得到得到基本估计图像u_PDE；所述基于各项异性扩散方程的图像去噪方法使用的偏微分方程具有以下形式：

u(x，0)＝f(x)，x∈Ω

步骤三：利用步骤二所述偏微分方程中的探测基本估计图像u_PDE的边界区域，并用阈值法去除边界，得到平滑无边界图像u_flat：

其中，α表示所述阈值法中的阈值；

步骤四：对步骤三所述平滑无边界图像u_flat按灰度值0,1,…,255进行分割，得到一组单灰度平滑区域

步骤五：分别计算步骤四所述单灰度平滑区域在噪声图像f中对应的区域的灰度值样本方差

其中，n_i表示中像素点个数，为样本均值，i＝0，1，2，...，255；

步骤六：对步骤五所述像素点个数n_i大于图像像素点总数1％的区域的灰度值样本方差取平均值，得到最终噪声方差估计其中P,Q分别表示图像的长宽像素点总数；

步骤七：按分块规则分别对步骤二所述噪声图像f和基本估计图像u_PDE进行分块处理，分别获得噪声图块{F_x}_x∈Ω和基本估计图块

所述分块规则为：以像素点x为左上角顶点，取8×8像素大小的正方形图像块；

步骤八：将步骤七所述基本估计图块进行分组，并获得基本估计图组

其中，x_R表示参考图块的左上角顶点坐标，||·||₂表示l²距离，再按照这种分组规则，将噪声图像和基本估计图像中的图像块分别堆叠成一系列三维图组：

噪声图块三维图组：

基本估计图块三维图组：

步骤九：对步骤八所述噪声图块三维图组做协同维纳滤波，得到图组估计：

其中，包含的图像块为T_3D是三维的正交变换，由二维离散余弦变换(2D-DCT)和一维Haar变换构成，为其逆变换；“＊”表示点对点的乘法运算；是维纳收缩系数，具体定义为

步骤十：将步骤九所述图组估计中的图像块还原到其原始位置，对重叠部分做加权平均处理，得到最终估计图：

其中，权值示性函数

进一步地，步骤二所述基本估计图像u_PDE的具体获得过程为：

第一步：对图像函数u的空间变量进行离散化处理并获得像素点x；离散后的像素点x的行、列坐标分别由p，q表示；图像的尺寸为P×Q像素，则p，q的取值范围分别为{0，1，2，...，P}和{0，1，2，...，Q}；

第二步：对图像函数u的时间变量进行离散化处理，取固定的时间步长τ，以迭代计算的方式实现方程在时间上的演化，方程中时间变量定义域的上界为T＝τN；N为总计迭代次数，并且N值取范围为100——200次；

第三步：定义方向导数的离散规则为

其中N，S，W，E为四个方向的简称，x可以替换为任何所需变量。

根据初值条件u(x，0)＝f(x)，x∈Ω，可获得图像函数u的迭代初始化条件为u⁰＝f；

根据边值条件可获得图像函数u的边界离散规则

u_p，0＝u_p，1，

u_p，Q＝u_p，Q-1，

u_0，q＝u_1，q，

u_P，q＝u_P-1，q.

结合边界处处理规则对n＝0，1，2，...，N-1依次按方程

的数值格式进行迭代计算：

其中，n表示当前的迭代次数，n的取值范围为{0，1，2，...，N-1}，经过N次迭代后，最终获得的u^N即为基本估计图像u_PDE。

进一步地，步骤二中所述偏微分方程中的函数结构采用其中，NLPM为该方程的名称标记，表示“非局部PM模型”，卷积尺度参数σ取值范围为0.5—1，参数K的取值范围为1—3。

进一步地，步骤三所述阈值α的取值范围为0.1—0.3。

本发明有益效果：

本发明所述一种基于偏微分方程和BM3D的无参数化图像恢复方法通过利用偏微分方程的处理以及省略BM3D参数处理的设计，使图像恢复方法在输入噪声标准差在40附近的情况下，以及当输入的噪声标准差非整数时，使PSNR完全不会产生震荡和突变，并且有效提高PSNR值和运行效率，极大程度上提高了图像恢复的成像质量。同时，大大提高了图像处理方法对图像恢复的处理的稳定性。同时，由于本发明所述图像恢复方法采用偏微分方程进行预处理，能够提高本发明图像恢复方法与噪声水平估计方法应用的稳定性，能够保证相同噪声水平下的不同噪声图像输入到本发明所述方法中进行处理，其图像恢复质量和处理结果皆为相同的最优结果。

附图说明

图1为本发明所述图像恢复方法的整体流程框图。

图2为传统BM3D方法的结构框图。

图3为本发明所述图像处理方法改进后的BM3D结构框图。

图4为本发明所述图像处理方法的噪声参数估计方法框架图。

图5为不同扩散方程图像去噪方法实验结果图比较图。

图6为扩散系数函数进行边缘检测的有效性效果图。

图7为本发明所述图像处理方法实验效果图一。

图8为本发明所述图像处理方法实验效果图二。

图9为本发明所述图像处理方法实验效果图三。

图10为本发明所述图像处理方法实验效果图四。

图11为本发明所述图像处理方法实验效果图五。

图12为本发明所述图像处理方法实验效果图六。

图13为本发明所述图像处理方法实验效果图七。

图14为本发明所述图像处理方法实验效果图八。

图15为本发明所述图像处理方法实验效果图九。

图16为本发明所述图像处理方法实验效果图十。

具体实施方式

下面结合具体实施例对本发明做进一步说明，但本发明不受实施例的限制。

实施例1：

一种基于偏微分方程和BM3D的无参数化图像恢复方法，所述图像恢复方法：首先，利用偏微分方程对噪声图像进行去噪预处理获得基本估计图，然后将所述基本估计图进行分块而获得基本估计图块，并利用距离匹配方法将所述基本估计图块匹配为基本估计图组；

其次，将所述基本估计图进行图像边界探测与去除处理后获得平滑无边界图像，通过灰度分割将所述平滑无边界图像分割出单灰度平滑区域，依次计算所述单灰度平滑区域中每个区域的区域样本方差，对获得的所述区域样本方差进行加权平均处理后获得噪声方差估计值；

在此，将噪声图像进行分块并获得噪声图块，将噪声图块按基本估计图块的距离匹配结果进行分组进而获得噪声图组；

最后通过协同维纳滤波方法，利用所述基本估计图组和所述噪声方差估计对所述噪声图组进行处理并获得图组估计，将所述图组估计以所述噪声方差估计为权重参数进行加权合成处理后获得最终估计图像；所述最终估计图像即为恢复后图像。

进一步地，所述图像修复的具体过程为：

步骤一：输入待处理的噪声图像f；

步骤二：利用基于各项异性扩散方程的图像去噪方法对噪声图像f做预处理，并得到得到基本估计图像u_PDE；所述基于各项异性扩散方程的图像去噪方法使用的偏微分方程具有以下形式：

u(x，0)＝f(x)，x∈Ω (2)

由于扩散方程图像去噪方法的核心思想就是将图像中的噪声视为“杂质”，从而将图像去噪过程视为“杂质扩散”的过程，因此以下将“去噪”和“扩散”视为同一过程。将图像视为“在整个去噪过程中，随时间变化的连续函数”，初始时刻为噪声图像，而结束时刻为去噪后的图像。

上述三个函数方程中，u＝u(x，t)即为连续变化的图像，其自变量的定义域为(x，t)∈Ω×(0，T)，Ω为图像所占的空间范围，T为整个图像去噪过程所用的时间。(1)式为用作图像去噪的扩散方程，其中div是散度符号，u_σ表示高斯卷积后的图像函数，具体为卷积尺度参数σ＞0，取值范围为0.5—1。(2)式为方程的初值条件，即初始时刻的图像函数u(x，t)等于噪声图像f，相当于为扩散方程(1)指明了处理对象；(3)式为方程的边值条件，它对扩散过程进行了限制，使得最终的扩散结果具有较好性质。

是一个关于和的函数，表示扩散方程的扩散系数，写成这种形式基本可以包括大多数的扩散方程图像去噪方法，不同方法的差异也基本是在这里体现。例如：

其中，K为扩散过程的相关参数，取值范围为1—3，(4)式为经典PM模型(“原始PM模型”)，(5)式为“正则化PM模型”，(6)式为本发明所述图像恢复方法中的创新模型，称为“非局部PM模型”。

如图5所示，(a)为噪声图像，原始PM方法处理图(b)的灰度平坦区域存在较多的黑白斑点，正则化PM方法处理图(c)在边界处发生模糊现象，而本发明提出的非局部PM方法可以有效避免以上两种缺陷。同时，在表1所示的处理结果PSNR指标上，采用新型偏微分方程也有效提高了PSNR，即比原始PM方法最大提高1.59dB，比正则化PM方法最大提高1.54dB。

表1不同扩散方程图像去噪方法处理结果PSNR数据比较

步骤三：利用步骤二所述偏微分方程中的探测基本估计图像u_PDE的边界区域，并用阈值法去除边界，得到平滑无边界图像u_flat：利用偏微分方程进行边界区域检测的有效性如图6所示，

其中，α表示所述阈值法中的阈值；所述阈值α的取值范围为0.1—0.3。

步骤四：对步骤三所述平滑无边界图像u_flat按灰度值0,1,…,255进行分割，得到一组单灰度平滑区域

步骤五：分别计算步骤四所述单灰度平滑区域在噪声图像f中对应的区域的灰度值样本方差

其中，n_i表示中像素点个数，为样本均值，i＝0，1，2，...，255；

步骤六：对步骤五所述像素点个数n_i大于图像像素点总数1％的区域的灰度值样本方差取平均值，得到最终噪声方差估计其中P,Q分别表示图像的长宽像素点总数；

步骤七：按分块规则分别对步骤二所述噪声图像f和基本估计图像u_PDE进行分块处理，分别获得噪声图块{F_x}_x∈Ω和基本估计图块

所述分块规则为：以像素点x为左上角顶点，取8×8像素大小的正方形图像块；

步骤八：将步骤七所述基本估计图块进行分组，并获得基本估计图组

其中，x_R表示参考图块的左上角顶点坐标，||·||₂表示l²距离，再按照这种分组规则，将噪声图像和基本估计图像中的图像块分别堆叠成一系列三维图组：

噪声图块三维图组：

基本估计图块三维图组：

步骤九：对步骤八所述噪声图块三维图组做协同维纳滤波，得到图组估计：

其中，包含的图像块为T_3D是三维的正交变换，由二维离散余弦变换(2D-DCT)和一维Haar变换构成，为其逆变换；“★”表示点对点的乘法运算；是维纳收缩系数，具体定义为

步骤十：将步骤九所述图组估计中的图像块还原到其原始位置，对重叠部分做加权平均处理，得到最终估计图：

其中，权值示性函数

其中，步骤二所述基本估计图像u_PDE的具体获得过程为：

第一步：对图像函数u的空间变量进行离散化处理并获得像素点x；离散后的像素点x的行、列坐标分别由p,q表示；图像的尺寸为P×Q像素，则p,q的取值范围分别为{0，1，2，...，P}和{0，1，2，...，Q}；

第二步：对图像函数u的时间变量进行离散化处理，取固定的时间步长τ，以迭代计算的方式实现方程在时间上的演化，方程中时间变量定义域的上界为T＝τN；N为总计迭代次数N取值因具体问题而异，在本实施例中，偏微分方程仅作为预处理操作，因此，N值取范围为100—200次。

第三步：定义方向导数的离散规则为

其中N,S,W,E为四个方向的简称，x可以替换为任何所需变量。

根据初值条件u(x，0)＝f(x)，x∈Ω，可获得图像函数u的迭代初始化条件为u⁰＝f；

根据边值条件可获得图像函数u的边界离散规则

u_p，0＝u_p，1，

u_p，Q＝u_p，Q-1，

u_0，q＝u_1，q，

u_P，q＝u_P-1,q.

结合边界处处理规则对n＝0，1，2，...，N-1依次按方程

的数值格式进行迭代计算：

其中，n表示当前的迭代次数，n的取值范围为{0，1，2，...，N-1}，经过N次迭代后，最终获得的u^N即为基本估计图像u_PDE。

传统的BM3D算法流程如图2所示：第一步的目的是为第二步提供预处理图像，便于第二步中图像块匹配得到更精确的结果，也便于维纳滤波过程中的参数计算。但图像中的噪声会影响第一步的处理结果，进而会对后续操作产生不利影响，导致技术问题中列举的多种问题。

本实施例所述一种基于偏微分方程和BM3D的无参数化图像恢复方法对传统BM3D算法框架进行了改进，利用偏微分方程预处理对噪声图像进行图像去噪的方式替代传统BM3D算法第一步硬阈值处理，具体改进后的BM3D算法框架如图2所示。将传统的BM3D算法第一步用偏微分方程预处理进行替代能够完全消除原始BM3D方法对噪声参数敏感的问题。

在对传统BM3D算法框架中第一步进行改进后，本实施例结合第一步的改进方法提出原创性的噪声水平估计方法，即基于偏微分方程的噪声参数估计方法(而非传统BM3D改进算法中直接利用某种现有的噪声水平估计方法从噪声图像中获得噪声参数估计，再输入原始BM3D方法的操作)，这样，能够实现噪声参数估计与图像去噪处理的完全融合，进一步提高图像恢复质量。改进后的BM3D方法框架如图3所示。

实验结果对比：

传统BM3D方法的图像处理实验结果如图7-9所示：

图7表示传统BM3D方法输入参数为40附近时，处理结果的PSNR值产生突变。图8表示传统BM3D方法处理结果PSNR在噪声水平输入值为非整数时产生震荡现象。图9表示当准确的噪声标准差并非BM3D方法的最佳输入参数时，真实噪声标准差均为20的两幅不同噪声图像经过传统BM3D方法处理后输出结果PSNR与噪声参数输入值的关系，可以发现，两组结果的最佳输入参数均与真实噪声水平不同，且彼此不同。说明传统方法的稳定性很差。

本发明所述图像恢复方法的实验结果如图10-16所示：

图10表示本发明所述图像恢复方法处理结果与噪声参数输入值的关系，

根据图10所示的本发明所述图像恢复方法处理结果的PSNR值随输入噪声参数值的变化情况，可以发现PSNR变化曲线具有更高的连续性，原始方法中出现的突变和震荡现象均消失。说明本发明所述图像恢复方法具有更高的稳定性。

图11展示了本发明所述无参数化图像恢复方法与估计参数的传统BM3D方法的处理结果，可以发现本发明所述无参数化图像恢复方法无论是PSNR值上，还是在运行效率上，新方法都具有显著优势，即有效的提高了是PSNR值和运行效率。

图12至图16为本发明所述方法实际应用的效果对比图，几组图像中，左侧图像均为单反相机拍摄得到的自然图像(带噪声)，右侧图像为无参数化图像恢复方法的恢复结果(去噪图像)。由于单反相机拍摄的自然图像中噪声信息未知，几组实验均无任何噪声参数的输入，完全由本发明所述方法进行自动处理，真正实现了无参数化图像恢复。实验结果均显示了本发明所述方法可以有效去除图像中的未知水平噪声；图像边界、纹理信息得到完好的保护与恢复；图像的色彩过渡区域仍保持较好的连续性，并未出现色块分离、阶梯效应等现象；同时恢复图像整体具有较好的视觉效果。

图15和图16组实验表明新方法同样可以处理彩色图像。实验结果表明，方法在处理彩色图像时，除具有与处理灰度图像相同的优点外(如噪声去除彻底、边界纹理保护完整、过渡区域无阶梯效应)，还具有色彩保持(即恢复图像与原始图像无色差)、无人为痕迹、视觉效果极佳等优点。

虽然本发明已以较佳的实施例公开如上，但其并非用以限定本发明，任何熟悉此技术的人，在不脱离本发明的精神和范围内，都可以做各种改动和修饰，因此本发明的保护范围应该以权利要求书所界定的为准。

Claims (3)

1.一种基于偏微分方程和BM3D的无参数化图像恢复方法，其特征在于，所述图像恢复方法：

首先，利用偏微分方程对噪声图像进行去噪预处理获得基本估计图，然后将所述基本估计图进行分块而获得基本估计图块，并利用距离匹配方法将所述基本估计图块匹配为基本估计图组；

其次，将所述基本估计图进行图像边界探测与去除处理后获得平滑无边界图像，通过灰度分割将所述平滑无边界图像分割出单灰度平滑区域，依次计算所述单灰度平滑区域中每个区域的区域样本方差，对获得的所述区域样本方差进行加权平均处理后获得噪声方差估计值；

在此，将噪声图像进行分块并获得噪声图块，将噪声图块按基本估计图块的距离匹配结果进行分组进而获得噪声图组；

最后通过协同维纳滤波方法，利用所述基本估计图组和所述噪声方差估计对所述噪声图组进行处理并获得图组估计，将所述图组估计以所述噪声方差估计为权重参数进行加权合成处理后获得最终估计图像；所述最终估计图像即为恢复后图像；

所述图像恢复的具体过程为：

步骤一：输入待处理的噪声图像f；

步骤二：利用基于各项异性扩散方程的图像去噪方法对噪声图像f做预处理，并得到得到基本估计图像u_PDE；所述基于各项异性扩散方程的图像去噪方法使用的偏微分方程具有以下形式：

u(x，0)＝f(x)，x∈Ω

其中，T为整个图像去噪过程所用的时间；

步骤三：利用步骤二所述偏微分方程中的探测基本估计图像u_PDE的边界区域，并用阈值法去除边界，得到平滑无边界图像u_flat：

其中，α表示所述阈值法中的阈值；

步骤四：对步骤三所述平滑无边界图像u_flat按灰度值0,1,…,255进行分割，得到一组单灰度平滑区域

步骤五：分别计算步骤四所述单灰度平滑区域在噪声图像f中对应的区域的灰度值样本方差

其中，n_i表示中像素点个数，为样本均值，i＝0，1，2，...，255；

步骤六：对步骤五所述像素点个数n_i大于图像像素点总数1％的区域的灰度值样本方差取平均值，得到最终噪声方差估计其中P,Q分别表示图像的长宽像素点总数；

步骤七：按分块规则分别对步骤二所述噪声图像f和基本估计图像u_PDE进行分块处理，分别获得噪声图块{F_x}_x∈Ω和基本估计图块

所述分块规则为：以像素点x为左上角顶点，取8×8像素大小的正方形图像块；

步骤八：将步骤七所述基本估计图块进行分组，并获得基本估计图组

其中，x_R表示参考图块的左上角顶点坐标，||·||₂表示l²距离，再按照这种分组规则，将噪声图像和基本估计图像中的图像块分别堆叠成一系列三维图组：

噪声图块三维图组：

基本估计图块三维图组：

步骤九：对步骤八所述噪声图块三维图组做协同维纳滤波，得到图组估计：

其中，包含的图像块为T_3D是三维的正交变换，由二维离散余弦变换(2D-DCT)和一维Haar变换构成，为其逆变换；“＊”表示点对点的乘法运算；是维纳收缩系数，具体定义为

步骤十：将步骤九所述图组估计中的图像块还原到其原始位置，对重叠部分做加权平均处理，得到最终估计图

其中，权值示性函数

步骤二所述基本估计图像u_PDE的具体获得过程为：

第一步：对图像函数u的空间变量进行离散化处理并获得像素点x；离散后的像素点x的行、列坐标分别由p,q表示；图像的尺寸为P×Q像素，则p,q的取值范围分别为{0，1，2，...，P}和{0，1，2，...，Q}；

第二步：对图像函数u的时间变量进行离散化处理，取固定的时间步长τ，以迭代计算的方式实现偏微分方程在时间上的演化，偏微分方程中时间变量定义域的上界为T＝τN；N为总计迭代次数，并且N值取范围为100—200次；

第三步：定义方向导数的离散规则为:

其中N,S,W,E为四个方向的简称，x可以替换为任何所需变量；

根据初值条件u(x，0)＝f(x)，x∈Ω，可获得图像函数u的迭代初始化条件为u⁰＝f；

根据边值条件可获得图像函数u的边界离散规则

u_p，0＝u_p，1，

u_p，Q＝u_p，Q-1，

u_0，q＝u_1，q，

u_P，q＝u_P-1，q。

结合边界处处理规则对n＝0，1，2，...，N-1依次按偏微分方程

的数值格式进行迭代计算：

其中，n表示当前的迭代次数，n的取值范围为{0，1，2，...，N-1}，经过N次迭代后，最终获得的u^N即为基本估计图像u_PDE。

2.根据权利要求1所述图像恢复方法，其特征在于，步骤二中所述偏微分方程中的函数结构采用卷积尺度参数σ取值范围为0.5-1，参数K的取值范围为1-3。

3.根据权利要求2所述图像恢复方法，其特征在于，步骤三所述阈值α的取值范围为0.1-0.3。