CN106875345A - 基于奇异值权重函数的非局部tv模型图像去噪方法 - Google Patents

Abstract

一种基于奇异值权重函数的非局部TV模型图像去噪方法。步骤如下：(1)首先输入噪声图像。(2)设置算法的相关参数，包括非局部搜索窗口大小N₁×N₁、邻域窗口大小N₂×N₂、像素相似度权重函数的参数h、j、高斯核的标准差σ，分裂的Bregman迭代辅助变量的初始值b⁰、保真参数λ，以及平滑参数θ。(3)通过奇异值分解方法获得图像块的最大奇异值。(4)构建基于该最大奇异值的新的像素相似度权重函数。(5)应用步骤(4)构建的权重函数，建立非局部TV模型。(6)对步骤(5)建立的非局部TV模型采用分裂的Bregman算法进行求解。(7)通过分裂的Bregman算法数值迭代运算获得去噪图像。(8)如果迭代满足停止条件，输出迭代最优结果图像并转向步骤(9)，如果不满足停止条件，则返回步骤(7)继续迭代。(9)将步骤(8)的迭代最优结果图像作为最后去噪结果图像。

Description

基于奇异值权重函数的非局部TV模型图像去噪方法

一、技术领域

本发明属于图像处理技术领域，具体涉及去除加性噪声的图像去噪领域和改进的非局部TV模型的图像去噪方法。

二、背景技术

图像去噪旨在通过对被噪声污染的图像进行某种处理，以降低噪声对原始有用信息的影响，尽可能地还原出被噪声污染前的图像。

GUY GILBOA和STANLEY OSHER提出的非局部TV模型(Nonlocal Total Variation,NLTV)是将GUY GILBOA和STANLEY OSHER提出的非局部算子(见文献：NONLOCAL OPERATORSWITH APPLICATIONS TO IMAGE PROCESSING.SIAM Multiscale Modeling andSimulation.Vol.7,No.3,pp.1005–1028)引入到Rudin-Osher-Fatemi提出的总变分(TotalVariation，TV)模型中而得到的，NLTV模型具有在去除图像噪声的同时能较好地保留图像的纹理细节的特点，该模型包含保真项和正则项。保真项表示观测图像(即待去噪的图像)与去噪后图像的接近程度；正则项表达去噪后图像光滑程度的先验知识。在正则项中，像素相似度权重函数具有很重要的作用，其表达式与含噪图像的图像块中各像素灰度值直接相关，但是由于噪声的存在，噪声会对像素灰度值产生干扰，进而影响像素相似度权重。

奇异值分解是一种基于特征向量的矩阵变换方法，在信号处理、模式识别、数字水印技术等方面都得到了应用。由于图像具有矩阵结构，因此本发明将其用于提取图像块的主要特征，将提取的图像主要特征用于构建新的像素相似度权重函数，以降低噪声对像素相似度权重的影响，再将构建的新的像素相似度权重函数应用到NLTV模型，得到新的NLTV模型，用于图像去噪。

本文中所提到的NLTV模型是指GUY GILBOA和STANLEY OSHER提出的非局部TV模型。

三、发明内容

本发明的目的是克服目前现有基于非局部理论的图像去噪技术中，在图像受到噪声干扰后，图像的像素值被噪声污染而改变，像素相似度权重赋值会受到噪声干扰的缺陷。本发明通过构建一种新型的像素相似度权重函数来减少噪声的干扰，提高像素相似度权重赋值的准确性。本发明所述的基于奇异值权重函数的非局部TV模型图像去噪方法包括如下步骤：

(1)首先输入N₀×N₀大小的待去噪的图像f；

(2)设置本发明方法的相关参数，包括非局部搜索窗口大小N₁×N₁、邻域窗口大小N₂×N₂、像素相似度权重函数的参数h和j，高斯核的标准差σ，分裂的Bregman迭代辅助变量b^k的初始值b⁰、平滑参数θ，以及保真参数λ；

(3)通过奇异值分解方法获得图像f中各领域图像块(图像块大小为N₂×N₂)的奇异值矩阵和最大奇异值；

(4)构建基于步骤(3)得到的最大奇异值的图像f的像素相似度权重函数；

(5)应用步骤(4)构建的权重函数，建立基于该权重函数的用于对图像f进行去噪的新的非局部TV模型；

(6)对步骤(5)建立的新的非局部TV模型，采用分裂的Bregman算法，通过逐次迭代逼近来求解；并设迭代计数变量初值k＝0；

(7)进行分裂的Bregman算法数值迭代运算，获得本次迭代后的输出图像u^k+1；

(8)每次完成步骤(7)的迭代后，计算出迭代结果u^k+1的评价值，即u^k+1的峰值信噪比(PSNR)；如果本次迭代后的输出图像u^k+1的峰值信噪比大于上一次迭代输出图像u^k的峰值信噪比，则不满足迭代停止条件，令k＝k+1，并返回步骤(7)，继续迭代运算；如果本次迭代结果u^k+1的峰值信噪比小于等于上一次迭代输出图像u^k的峰值信噪比，则满足迭代停止条件，将上一次迭代输出图像u^k作为最优值u_final输出，即令u_final＝u^k，并转到步骤(9)；

(9)将结果u_final作为最后去噪结果图像输出。

本发明的优点是建立新的像素相似度权重函数。将图像块进行奇异值分解并选取最大奇异值，将选取的最大奇异值用于像素相似度权重函数的构建，减少了噪声对权重函数的干扰，能更准确的对像素相似度权重进行赋值。将由此构建的基于图像块的像素相似度权重函数应用到NLTV模型，用于图像去噪，提高图像去噪的效果。本发明方法适合用于所含噪声模型为加性噪声的含噪图像的去噪。

四、附图说明

图1为本发明方法的流程图。

图2是原始图像未加噪的图像，其中，图2(11)是Peppers图，图2(12)是Cameraman图，图2(13)是Lena图。

图3为本发明对Lena图去噪仿真结果，其中，图3(11)是待去噪的图像(是在原始图像上叠加了均值为零标准差为10的高斯噪声后的结果)，图3(12)非局部TV方法去噪，图3(13)本发明方法去噪。

五、具体实施方式

本发明所述的基于奇异值权重函数的非局部TV模型图像去噪方法，具体步骤如下：

(1)首先输入N₀×N₀大小的待去噪图像f；

(2)设置本发明方法的相关参数，包括非局部搜索窗口大小N₁×N₁、邻域窗口大小N₂×N₂、像素相似度权重函数的参数h和j，高斯核的标准差σ，分裂的Bregman迭代辅助变量b^k的初始值b⁰、平滑参数θ和保真参数λ；

(3)设M_x是以步骤(1)输入的含噪图像f中像素点x∈Ω为中心的大小为N₂×N₂的图像块像素灰度值矩阵,Ω为f的图像空间，不同的像素点对应不同的图像块。将各图像块M_x进行奇异值分解：M_x＝U_xΛ_xV_x ^T。式中U_x、V_x分别是M_x的左奇异矩阵和右奇异矩阵，大小都是N₂×N₂。Λ_x是M_x的奇异值矩阵，它的大小也是N₂×N₂，它的对角元素不为零，其他元素全为零，Λ_x对角元素一共有N₂个：按从大到小排列为：Λ_x的对角元素就是M_x的奇异值，它包含了M_x的全部特征。奇异值的定义：对于m×n阶矩阵A，A^TA的n个特征值的非负平方根叫作A的奇异值；

(4)利用步骤(3)得到的图像块的奇异值构建新的像素相似度权重函数。由于奇异值包含了图像块的主要特征，所以相似的图像块之间的奇异值是相近的；在图像块的奇异值 中，大的奇异值包含了图像块的主要特征，小的奇异值包含了图像块的次要特征；并且在含有噪声的图像中，噪声不是图像的主要特征，所以在构建新的像素相似度权重函数时，只选用图像块的最大奇异值，通过图像块的主要特征来判断图像块之间的相似性，这样就减少了噪声的干扰；构建图像f中两个像素点x和y的相似度权重函数：

其中x为当前像素点，y是以x为中心的搜索窗口内的一点，a_x是以x为中心、大小为N₂×N₂的邻域内的像素灰度值矩阵M_x的最大奇异值，a_y是以y为中心、大小为N₂×N₂的邻域的像素灰度值矩阵M_y的最大奇异值；表示在图像f中以x为中心的大小为N₂×N₂的图像块与以y为中心的大小为N₂×N₂的图像块之间的高斯加权距离，表示求和范围是以x或y为中心的大小为N₂×N₂的邻域内的每一像素点(不包括x或y自身)，共N₂×N₂-1项，G_σ(c)是标准差为σ的高斯核函数；h和j是常数，通过干预指数函数的衰减速度来控制权重函数ω的大小,h和j的取值越大，权重函数ω的值越接近1，算法收敛速度快，但是难以达到最优值，h和j的取值越小，权重函数ω的值越接近0，算法经过多次迭代可以收敛到最优值，但是耗费时间多，h和j的取值原则上要综合以上两点因素，取大小合适的值；

(5)建立非局部TV模型：其中J(u)是模型的目标函数，λ是保真参数，f是输入的含噪图像，u是去噪后得到的图像，Ω是f、u的图像空间，是GUY GILBOA和STANLEY OSHER提出的非局部梯度算子，u(y)和u(x)分别是图像u上的点x和y的像素灰度值，

(6)对步骤(5)建立的NLTV模型，采用分裂的Bregman算法迭代逼近来求解，迭代过程分成三步，降低了迭代的复杂性，加快了迭代速度。迭代格式如下：

其中，k的取值是0，1，2，…,等非负整数，迭代初始值u⁰＝f，b^k、w^k分别表示分裂的Bregman迭代的辅助变量和辅助函数，λ就是步骤(5)建立的非局部TV模型中的保真参数，θ是控制迭代结果的平滑参数；b^k的初始值b⁰、以及λ和θ的赋值已在步骤(2)中进行预设。

求解式(6-1)和式(6-2)，并数值化，式(6-3)也数值化，由此得到数值化后的三步迭代格式，如步骤(7)中的式(7-1)、(7-2)和式(7-3)所示。第一次迭代运算前，令k＝0；

(7)顺序地应用公式(7-1)、(7-2)、(7-3)进行迭代，

其中表示求和范围是以x为中心的搜索窗口N₁×N₁内的每一像素点y(不包括x自身)，共N₁×N₁-1项，ω(x,y)是f中像素x和y间的相似度权重函数，由步骤(4)计算得到；

(8)每次完成步骤(7)的迭代后，计算出迭代结果u^k+1的评价值，即u^k+1的峰值信噪比(PSNR)；如果本次迭代后的输出图像u^k+1的峰值信噪比大于上一次迭代输出图像u^k的峰值信噪比，则不满足迭代停止条件，令k＝k+1，并返回步骤(7)，继续迭代运算；如果本次迭代结果u^k+1的峰值信噪比小于等于上一次迭代输出图像u^k的峰值信噪比，则满足迭代停止条件，将上一次迭代输出图像u^k作为最优值u_final输出，即令u_final＝u^k，并转到步骤(9)；

(9)将结果u_final作为最后去噪结果图像输出。

本发明去噪效果可以通过以下实验进一步证实：

(一)实验条件

使用Matlab软件对如附2所示的256×256尺寸的Peppers、Cameraman、Lena灰度图像进行测试。在原图像中中分别添加平均值为零标准差为10,15,20的高斯噪声后作为待去噪的图像输入。本发明的参数为：对于标准差为10的含噪图像，去噪参数为搜索窗口5×5、邻域窗口5×5、h＝9.4、j＝31.62、σ＝5、b⁰＝0、w⁰＝0、λ＝1、θ＝12。对于标准差为15的含噪声图像，去噪参数为搜索窗口5×5、邻域窗口5×5、h＝14、j＝31.62、σ＝6、b⁰＝0、w⁰＝0、λ＝1、θ＝11.5。对于标准差为20的含噪图像，去噪参数为搜索窗口5×5、邻域窗口5×5、h＝17、j＝31.62、σ＝5、b⁰＝0、w⁰＝0、λ＝1、θ＝9。迭代停止的条件是每迭代一次，都计算出评价值，即迭代结果u^k+1的峰值信噪比PSNR，如果本次迭代的评价值大于上一次的评价值，则迭代继续，如果本次迭代的评价值小于等于上一次的评价值，则迭代停止，并将上一次迭代的结果u^k作为最优值输出。这样就得到了去噪效果最好的图像。

用来对比的NLTV去噪方法参数：对于标准差为10的含噪图像，去噪参数为搜索窗口5×5、邻域窗口5×5、h＝9、σ＝5、b⁰＝0、w⁰＝0、λ＝1、θ＝12。对于标准差为15的噪声图像去噪参数为搜索窗口5×5、邻域窗口5×5、h＝13、σ＝5、b⁰＝0、w⁰＝0、λ＝1、θ＝11.5。对于标准差为20的含噪图像，去噪参数为搜索窗口5×5、邻域窗口5×5、h＝17、σ＝5、b⁰＝0、w⁰＝0、λ＝1、θ＝9。迭代停止的条件是每迭代一次，都计算出评价值，即峰值信噪比PSNR，如果本次迭代的评价值大于上一次的评价值，则返回到步骤(6)继续迭代，如果本次迭代的评价值小于等于上一次的评价值，则迭代停止，并将上一次迭代的结果作为最优值输出。这样就得到了去噪效果最好的图像。

(二)实验内容

按照上面所述的实验步骤进行实验仿真,并将本发明方法与非局部TV方法进行比较。

(三)实验结果

实验结果见表一和附图3。实验结果表明本发明去噪方法与NLTV模型去噪方法(也采用分裂的Bregman算法迭代逼近来求解)相比，本发明方法去噪后图像的峰值信噪比更高，去噪效果更好。

以上所述,仅是本发明的较佳实施例而已,并不对本发明做形式上的限制,凡是依据本发明对以上实例所做的简单修改，等同变化与修饰,均仍属本发明技术方案的范围内。

表一用本发明方法和NLTV方法去噪前后图像的峰值信噪比。

Claims (1)

1.基于奇异值权重函数的非局部TV模型图像去噪方法，包括如下步骤：

(1)首先输入N₀×N₀大小的待去噪图像f；

(2)设置本发明方法的相关参数，包括非局部搜索窗口大小N₁×N₁、邻域窗口大小N₂×N₂、像素相似度权重函数的参数h和j，高斯核的标准差σ，分裂的Bregman迭代辅助变量b^k的初始值b⁰、平滑参数θ和保真参数λ；

(3)设M_x是以步骤(1)输入的含噪图像f中像素点x∈Ω为中心的大小为N₂×N₂的图像块像素灰度值矩阵,Ω为f的图像空间，不同的像素点对应不同的图像块。将各图像块M_x进行奇异值分解：M_x＝U_xΛ_xV_x ^T。式中U_x、V_x分别是M_x的左奇异矩阵和右奇异矩阵，大小都是N₂×N₂。Λ_x是M_x的奇异值矩阵，它的大小也是N₂×N₂，它的对角元素不为零，其他元素全为零，Λ_x对角元素一共有N₂个：按从大到小排列为：Λ_x的对角元素就是M_x的奇异值，它包含了M_x的全部特征。奇异值的定义：对于m×n阶矩阵A，A^TA的n个特征值的非负平方根叫作A的奇异值；

(4)利用步骤(3)得到的图像块的奇异值构建新的像素相似度权重函数。由于奇异值包含了图像块的主要特征，所以相似的图像块之间的奇异值是相近的；在图像块的奇异值 中，大的奇异值包含了图像块的主要特征，小的奇异值包含了图像块的次要特征；并且在含有噪声的图像中，噪声不是图像的主要特征，所以在构建新的像素相似度权重函数时，只选用图像块的最大奇异值，通过图像块的主要特征来判断图像块之间的相似性，这样就减少了噪声的干扰；构建图像f中两个像素点x和y的相似度权重函数：

$\mathrm{\ω}(x,y)=e^{-\{\frac{\underset{c}{\Σ}\[G_{\mathrm{\σ}}\left(c\right)|f(x+c)-f(y+c){|}^2\]}{h^2}+\frac{{(a_x-a_y)}^2}{j^2}\}},$

其中x为当前像素点，y是以x为中心的搜索窗口内的一点，a_x是以x为中心、大小为N₂×N₂的邻域内的像素灰度值矩阵M_x的最大奇异值，a_y是以y为中心、大小为N₂×N₂的邻域的像素灰度值矩阵M_y的最大奇异值；表示在图像f中以x为中心的大小为N₂×N₂的图像块与以y为中心的大小为N₂×N₂的图像块之间的高斯加权距离，表示求和范围是以x或y为中心的大小为N₂×N₂的邻域内的每一像素点(不包括x或y自身)，共N₂×N₂-1项，G_σ(c)是标准差为σ的高斯核函数；h和j是常数，通过干预指数函数的衰减速度来控制权重函数ω的大小,h和j的取值越大，权重函数ω的值越接近1，算法收敛速度快，但是难以达到最优值，h和j的取值越小，权重函数ω的值越接近0，算法经过多次迭代可以收敛到最优值，但是耗费时间多，h和j的取值原则上要综合以上两点因素，取大小合适的值；

(5)建立非局部TV模型：其中J(u)是模型的目标函数，λ是保真参数，f是输入的含噪图像，u是去噪后得到的图像，Ω是f、u的图像空间，是GUY GILBOA和STANLEY OSHER提出的非局部梯度算子，u(y)和u(x)分别是图像u上的点x和y的像素灰度值，

(6)对步骤(5)建立的NLTV模型，采用分裂的Bregman算法迭代逼近来求解，迭代过程分成三步，降低了迭代的复杂性，加快了迭代速度。迭代格式如下：

$w^{k+1}=\arg \underset{w}{\min }\{\∫|w|dx+\frac{\mathrm{\θ}}{2}\∫{(w-{\▿}_{NL}{u}^k-b^k)}^{2}dx\}---(6-1)$

$u^{k+1}=\arg \underset{u}{\min }\{\mathrm{\λ}\∫{(u-f)}^{2}dx+\frac{\mathrm{\θ}}{2}\∫{(w^{k+1}-{\▿}_{NL}u-b^k)}^{2}dx\}---(6-2)$

$b^{k+1}=b^k+{\▿}_{NL}{u}^{k+1}-w^{k+1}---(6-3)$

其中，k的取值是0，1，2，…,等非负整数，迭代初始值u⁰＝f，b^k、w^k分别表示分裂的Bregman迭代的辅助变量和辅助函数，λ就是步骤(5)建立的非局部TV模型中的保真参数，θ是控制迭代结果的平滑参数；b^k的初始值b⁰、以及λ和θ的赋值已在步骤(2)中进行预设。

求解式(6-1)和式(6-2)，并数值化，式(6-3)也数值化，由此得到数值化后的三步迭代格式，如步骤(7)中的式(7-1)、(7-2)和式(7-3)所示。第一次迭代运算前，令k＝0；

(7)顺序地应用公式(7-1)、(7-2)、(7-3)进行迭代，

$w^{k+1}(x,y)=\max (\sqrt{\underset{y}{\Σ}{\{\sqrt{\mathrm{\ω}(x,y)}\[u^k\left(y\right)-u^k\left(x\right)\]+b^k(x,y)\}}^2}-\frac{1}{\mathrm{\θ}},0)\frac{\sqrt{\mathrm{\ω}(x,y)}\[u^k\left(y\right)-u^k\left(x\right)\]+b^k(x,y)}{\sqrt{\underset{y}{\Σ}{\{\sqrt{\mathrm{\ω}(x,y)}\[u^k\left(y\right)-u^k(x,y)\]+b^k(x,y)\}}^2}}---(7-1)$

$u^{k+1}\left(x\right)=\frac{{\mathrm{\λu}}^k\left(x\right)+\frac{\mathrm{\θ}}{2}\underset{y}{\Σ}\{2u^k\left(y\right)\mathrm{\ω}(x,y)-\sqrt{\mathrm{\ω}(x,y)}\[w^{k+1}(x,y)-w^{k+1}(y,x)+b^k(y,x)-b^k(x,y)\]\}}{\mathrm{\λ}+\mathrm{\θ}\underset{y}{\Σ}\mathrm{\ω}(x,y)}---(7-2)$

$b^{k+1}(x,y)=b^k(x,y)+\sqrt{\mathrm{\ω}(x,y)}\[u^{k+1}\left(y\right)-u^{k+1}\left(x\right)\]-w^{k+1}(x,y)---(7-3)$

其中表示求和范围是以x为中心的搜索窗口N₁×N₁内的每一像素点y(不包括x自身)，共N₁×N₁-1项，ω(x,y)是f中像素x和y间的相似度权重函数，由步骤(4)计算得到；

(8)每次完成步骤(7)的迭代后，计算出迭代结果u^k+1的评价值，即u^k+1的峰值信噪比(PSNR)；如果本次迭代后的输出图像u^k+1的峰值信噪比大于上一次迭代输出图像u^k的峰值信噪比，则不满足迭代停止条件，令k＝k+1，并返回步骤(7)，继续迭代运算；如果本次迭代结果u^k+1的峰值信噪比小于等于上一次迭代输出图像u^k的峰值信噪比，则满足迭代停止条件，将上一次迭代输出图像u^k作为最优值u_final输出，即令u_final＝u^k，并转到步骤(9)；

(9)将结果u_final作为最后去噪结果图像输出。