CN103077506A - 结合局部和非局部的自适应图像去噪方法 - Google Patents

Abstract

本发明公开了一种结合局部和非局部的自适应图像去噪方法，主要解决现有去噪方法去噪效果不佳问题。其实现步骤：①输入含噪图像；②估计含噪图像的噪声标准差；③以含噪图像的任意像素点为中心提取像素矢量，并计算像素矢量非局部均值；④对含噪图像的所有像素点进行步骤3操作；⑤对所有的像素矢量进行聚类，对每类像素矢量子集训练字典；⑥对所有的像素矢量分别进行自适应去噪；⑦将所有的去噪像素矢量拉成图像块，并聚集得到去噪图像；⑧判断迭代是否完成，若完成则输出去噪图像，否则将去噪图像作为含噪图像转到步骤2进入下一步迭代。本发明能有效去除含有高斯白噪声的自然图像中的噪声，可用于医学影像、视频多媒体等领域的数字图像预处理。

Description

结合局部和非局部的自适应图像去噪方法

技术领域

本发明属于图像处理技术领域，具体的说是一种结合局部和非局部的自适应图像去噪方法，可用于医学影像、天文学影像、视频多媒体等领域的数字图像预处理。

背景技术

图像去噪技术解决了图像在获取、编码及传输等过程中受到的各种噪声干扰所导致的图像质量下降问题，提高了图像质量，是图像处理中的重要环节和研究内容。

图像去噪技术大致是从空域和变换域两方面来进行的，目前效果较好的空域滤波方法有非局部均值滤波去噪方法NLM、稀疏表示下的图像去噪方法KSVD等，较好的变换域的滤波方法有三维块匹配去噪方法BM3D等。

非局部均值滤波去噪方法NLM，是Buades等人在文献“A non local algorithm forimage denoising.IEEE Conference on Computer Vision and Pattern Recognition，2005,vol.2,pp:60-65.”中提出来的。该方法通过计算图像中两个像素点邻域的相似度来确定该点对所要求的点的信息补偿程度，中心点的灰度值为邻域内像素点灰度值的加权平均。该方法存在的不足是，过于强调邻域内像素点的非局部信息而忽略了中心像素点自身的局部信息，导致图像边缘细节信息丢失严重。

稀疏表示下的图像去噪方法KSVD，是Elad等人在文献“Image denoising via spareand redundant representations over learned dictionaries.IEEE Transactions on ImageProcessing，2006,vol.15no.12,pp.3736-3745.”中提出来的。该方法是用图像块在冗余字典上的稀疏近似来实现图像去噪，图像块在冗余字典上越稀疏，则逼近原图像块的信息就越准确。该方法存由于只是对图像块自身进行稀疏逼近，利用的只是图像块自身的局部信息，而忽视了邻域图像块的非局部信息，因而对图像的平滑区域去噪效果不理想。

三维块匹配去噪方法BM3D，是Dabov等人在文献“Image denoising by sparse3-Dtransform-domain collaborative filtering.IEEE Transactions on Image Processing，2007,vol.16no.8,pp.2080-2095.”中提出来的。该方法不仅利用了图像的结构信息，而且结合了变换域的阈值方法，是一种结合局部、非局部、多尺度稀疏的自适应非参数化的滤波技术。但该方法仍然存在的不足是，对图像块稀疏逼近采用的是离散余弦变换DCT字典，由于该字典是一个固定字典，不能有效地逼近原图像的边缘和细节信息，易造成原图中部分边缘和纹理细节丢失。

发明内容

本发明的目的在于针对上述去噪方法的不足，提出一种结合局部和非局部的自适应图像去噪方法，以保持图像的边缘和纹理细节，提高图像的去噪效果。

实现本发明的技术思路是：利用基于非局部均值NLM的自适应去噪算法对图像块进行去噪，并通过BM3D的聚集方法对去噪后的图像块进行重构获得去噪图像，其具体实现步骤包括如下：

（1）输入一幅N行M列的含噪图像Y，设置最大迭次数γ和停止参数δ；

（2）采用下式估计含噪图像Y的噪声标准差σ_n：

${\mathrm{\σ}}_n=\frac{\mathrm{median}\left(\mathrm{abs}\right|W\left|\right)}{0.6745},$

其中，W为含噪图像Y进行小波分解得到的第一层高频系数，abs|·|是取绝对值操作，median(·)是取中值操作；

（3）以含噪图像Y中的任意像素点为中心，分别确定一个l×l大小的位置区域Γ_i和一个s×s大小的搜索窗Ω_i，提取含噪图像Y的位置区域Γ_i的像素得到图像块y_i，其中i∈{1，...,N×M}，l、s的取值范围分别为5~11和21~41奇数个像素；

（4）将l×l大小的图像块y_i按列拉成l²×1大小的像素矢量v_i，按照下式分别计算像素矢量v_i和搜索窗Ω_i内所有像素矢量集合{v_ip}之间的相似度，得到相似度集合{w_ip}：

$w_{\mathrm{ip}}=\exp \left(\frac{-{\left|\right|v_i-v_{\mathrm{ip}}\left|\right|}_2^2}{h}\right),$

其中，p=1，...,s²，

表示二范数，h为衰减因子，h的取值范围10σ_n~15σ_n；

（5）对相似度集合{w_ip}中的s²个相似度按照从大到小排序，选取前k个相似度作为像素矢量v_i的相似度集合{w_ij}，将该相似度集合{w_ij}对应的像素矢量作为像素矢量v_i的相似像素矢量集合{v_ij}，其中j=1，...,k，k的取值范围为10~30；

（6）计算相似像素矢量集合{v_ij}的像素矢量非局部均值

${\stackrel{\_}{v}}_i=\frac{\underset{j=1}{\overset{k}{\mathrm{\Σ}}}{w}_{\mathrm{ij}}{v}_{\mathrm{ij}}}{\underset{j=1}{\overset{k}{\mathrm{\Σ}}}{w}_{\mathrm{ij}}};$

（7）在含噪图像Y中，对所有的像素点分别进行步骤（3）至步骤（6）操作，得到N×M个像素矢量v_i，N×M个相似像素矢量集合{v_ij}，N×M个像素矢量非局部均值

（8）采用K均值k-means聚类方法将N×M个像素矢量v_i分成f类，得到f个像素矢量子集U_t和f个聚类中心c_t，其中t∈{1，...,f}，f的取值范围为15~60；

（9）采用主成分分析PCA方法分别对f个像素矢量子集U_t训练一个l²×l²大小的字典D_t；

（10）利用f个字典D_t对N×M个像素矢量v_i分别进行自适应去噪，得到N×M个去噪像素矢量

（11）将N×M个l²×1大小的去噪像素矢量

分别按列拉成N×M个l×l大小的去噪图像块

（12）采用基于三维块匹配去噪BM3D的图像块聚集方法对N×M个去噪图像块

进行聚集，得到去噪图像X；

（13）计算含噪图像Y和去噪图像X之间的均方误差e(X,Y)；

（14）判断均方误差e(X,Y)是否小于停止参数δ或者本次迭代次数是否达到最大迭代次数γ，若满足两个条件中的一个，则终止迭代，输出去噪图像X，否则将去噪图像X作为含噪图像Y，重复步骤（2）到步骤（13）。

本发明与现有的技术相比具有以下优点：

a）本发明采用了基于聚类的字典学习的稀疏表示方式，与现有技术采用的固定字典稀疏表示方式相比，具有更强的稀疏能力。

b）本发明是基于非局部均值NLM的自适应图像去噪，结合了图像的局部结构特性和非局部相似性特点，与现有技术相比，能更有效地保持图像的边缘和纹理细节。

c）本发明采用了一种基于噪声重新估计的迭代去噪方式，有效地降低了前一步去噪结果对后一步去噪效果的影响。

附图说明

图1是本发明的实现流程图；

图2是本发明仿真使用的未加噪图像；

图3是本发明仿真使用的加噪图像；

图4是用现有方法和本发明方法对图3的去噪效果比较图。

具体实施方式

以下参照附图对本发明的具体实现及效果做进一步的详细说明。

参照图1，本发明的实现步骤如下：

步骤1，输入一幅N行M列的含噪图像Y，设置最大迭次数γ和停止参数δ，其中γ、δ的取值范围分别为9~15和0.01~0.03，本实例中γ、δ的取值分别为12和0.02。

步骤2，采用下式估计含噪图像Y的噪声标准差σ_n：

${\mathrm{\σ}}_n=\frac{\mathrm{median}\left(\mathrm{abs}\right|W\left|\right)}{0.6745},$

其中，W为含噪图像Y进行小波分解得到的第一层高频系数，abs|·|是取绝对值操作，median(·)是取中值操作。

步骤3，以含噪图像Y中的任意像素点为中心,分别确定一个l×l大小的位置区域Γ_i和一个s×s大小的搜索窗Ω_i，提取含噪图像Y的位置区域Γ_i的像素得到图像块y_i，其中i∈{1,...,N×M}，l、s的取值范围分别为5~11和21~41奇数个像素，本实例中l、s的取值分别为12、31。

步骤4，将l×l大小的图像块y_i按列拉成l²×1大小的像素矢量v_i，按照下式分别计算像素矢量v_i和搜索窗Ω_i内所有像素矢量集合{v_ip}之间的相似度，得到相似度集合{w_ip}：

$w_{\mathrm{ip}}=\exp \left(\frac{-{\left|\right|v_i-v_{\mathrm{ip}}\left|\right|}_2^2}{h}\right),$

其中，p=1，...,s²，h为衰减因子，

表示二范数，h的取值范围10σ_n~15σ_n，本实例中h的取值为12σ_n。

步骤5，对相似度集合{w_ip}中的s²个相似度按照从大到小排序，选取前k个相似度作为像素矢量v_i的相似度集合{w_ij}，将该相似度集合{w_ij}对应的像素矢量作为像素矢量v_i的相似像素矢量集合{v_ij}，其中j=1，...,k，k的取值范围为10~30，本实例中k的取值为16。

步骤6，计算相似像素矢量集合{v_ij}的像素矢量非局部均值

${\stackrel{\_}{v}}_i=\frac{\underset{j=1}{\overset{k}{\mathrm{\Σ}}}{w}_{\mathrm{ij}}{v}_{\mathrm{ij}}}{\underset{j=1}{\overset{k}{\mathrm{\Σ}}}{w}_{\mathrm{ij}}}.$

步骤7，在含噪图像Y中，对所有的像素点分别进行步骤3至步骤6操作，得到N×M个像素矢量v_i，N×M个相似像素矢量集合{v_ij}，N×M个像素非局部矢量均值

步骤8，采用K均值k-means聚类方法将N×M个像素矢量v_i分成f类，得到f个像素矢量子集U_t和f个聚类中心c_t，其中t∈{1,...,f}，f的取值范围为15~60本实例中f的取值为40。

步骤9，采用主成分分析PCA方法分别对f个像素矢量子集U_t训练一个l²×l²大小的字典D_t。

步骤10，利用步骤9中得到的f个字典D_t，对N×M个像素矢量v_i分别进行自适应去噪，得到N×M个去噪像素矢量

10a）分别计算像素矢量v_i与f个聚类中心c_t之间的欧式距离d(v_i,c_t)：

$d(v_i,c_t)={\left|\right|v_i-c_t\left|\right|}_2^2;$

10b）在f个欧式距离d(v_i,c_t)中选取欧式距离d(v_i,c_t)最小值时c_t对应的下角标g，其中g∈{1,...,f}；

10c）从步骤9中得到的f个PCA字典D_t中选取下角标为g的字典，得到稀疏字典D_g；

10d）利用上述稀疏字典D_g以及步骤7中得到的像素矢量非局部均值

和相似像素矢量集合{v_ij}，计算自适应参数矢量λ：

$\mathrm{\λ}=\frac{{\left[\frac{1}{k}\underset{j=1}{\overset{k}{\mathrm{\Σ}}}\right[\mathrm{diag}\left(D_g^T({\stackrel{\‾}{v}}_i-v_{\mathrm{ij}})\right){]}^2-{\mathrm{\σ}}_n^{2}I]}_{+}}{{\mathrm{\σ}}_n^{2}I},$

其中，I为单位矩阵，σ_n为步骤2中计算得到的图像的噪声标准差，diag(·)为将矢量转成对角化矩阵的操作，

表示稀疏字典D_g的转置，[·]₊为矩阵取正运算；

10e）利用稀疏字典D_g和像素矢量均值

计算像素矢量v_i的去噪像素矢量

${\hat{v}}_i=D_g\frac{\left(D_g^T{\stackrel{\‾}{v}}_i\right)+\mathrm{\λ}*\left(D_g^T{v}_i\right)}{I+\mathrm{\λ}};$

10f）对N×M个像素矢量分别重复步骤10a)至步骤10e)，得到N×M个去噪像素矢量。

步骤11，将N×M个l²×1大小的去噪像素矢量分别按列拉成N×M个l×l大小的去噪图像块

步骤12，采用基于三维块匹配去噪BM3D的图像块聚集方法对N×M个去噪图像块

进行聚集，得到去噪图像X，其中，采用的图像块聚集方法来源于文献“Imagedenoising by sparse3-D transform-domain collaborative filtering.IEEE Transactions onImage Processing，2007,vol.16no.8,pp.2080-2095.”。

12a）初始化两个大小为N×M的零矩阵，得到标记矩阵R和重构矩阵Q；

12b）采用形状参数为2的凯瑟Kaiser窗，初始化一个大小为l×l的聚集矩阵Θ；

12c）采用下式分别对N×M个去噪图像块

进行聚集加权，得到N×M个加权去噪图像块z_i：

$z_i=<\mathrm{\&Theta;},{\hat{y}}_i>,$

其中，i∈{1，...,N×M}，<·,·>表示矩阵的内积运算；

12d）根据步骤3中得到的图像块y_i所在含噪图像Y中的位置区域Γ_i，分别在标记矩阵R和重构矩阵Q的相同位置确定一个相同大小的区域，得到标记区域Φ_i和重构区域ψ_i；

12e）对标记矩阵R中区域Φ_i的元素值加上聚集矩阵Θ的元素值，得到标记更新矩阵

12f）对重构矩阵Q中区域ψ_i的元素值加上加权去噪图像块z_i的元素值，得到重构更新矩阵

12g）对标记矩阵R的所有位置进行步骤12d）至步骤12f）操作，得到最终标记更新矩阵

和最终重构更新矩阵

12h）利用上述的最终标记更新矩阵

和最终重构更新矩阵

计算去噪图像X：

$X=\mathrm{div}(\hat{Q},\hat{R}),$

其中，div(·,·)表示矩阵的点除运算，即

表示矩阵

中的每个元素分别除以矩阵

相应位置上的元素。

步骤13，计算含噪图像Y和去噪图像X之间的均方误差e(X,Y)：

$e(X,Y)=\frac{1}{256*\mathrm{NM}}{\left|\right|X-Y\left|\right|}_2^2.$

步骤14，判断均方误差e(X,Y)是否小于停止参数δ或者本次迭代次数是否达到最大迭代次数γ，若满足两个条件中的一个，则终止迭代，输出去噪图像X；否则将去噪图像X作为含噪图像Y，重复步骤2到步骤13。

本发明的效果可以进一步通过以下内容进行说明：

1．仿真条件：

在CPU为pentium4处理器，主频为1.86GHZ，内存2G，操作系统为WINDOWSXP SP3，仿真平台为Matlab7.10的计算机上进行。

2．评价指标：

图像去噪效果的评价分为主观和客观两个方面。在主观上评价一幅图像去噪效果的优劣主要是通过人眼的视觉特性来衡量，图像质量好，感觉清晰则去噪效果好，反之效果则差；在客观上评价一幅图像的去噪效果，本发明采用峰值信噪比PSNR来衡量。令真实图像为X，去噪图像为Y，X(t)为真实图像X中第t个像素的灰度值，X′(t)为去噪图像X′中第t个像素的幅值，Ω为图像中的像素集合，令X_max=max{X(t),t∈Ω}，则峰值信噪比按如下公式计算：

$\mathrm{RSNR}=10\log \frac{X_{\max }^2}{\frac{1}{\left|\mathrm{\&Omega;}\right|}\underset{t\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}[X\&prime;\left(t\right)-X\left(t\right){]}^2},$

其中，|Ω|表示图像像素集合Ω中的像素个数。

3．仿真图像：

仿真图像选择图2所示的四幅原始测试图像，其中，图2（a）为Lena图像，图2（b）为House图像，图2（c）为Barbara图像，图2（d）为Man图像，四幅图像均为灰度图像，灰度级为256。图3为图2（c）Barbara图像加入噪声方差为35的高斯随机白噪声后的图像。

4.仿真内容：

对比方法1：Buades等在文献“A non local algorithm for image denoising.IEEEConference on Computer Vision and Pattern  Recognition，2005,vol.2,60-65.”中提出的图像非局部均值滤波去噪方法，简称NLM；

对比方法2：Elad等在文献“Image denoising via spare and redundant representationsover learned dictionaries.IEEE Transactions on Image Processing，2006,vol.15no.12,pp.3736-3745.”中提出的基于稀疏表示下的图像去噪方法，简称KSVD；

对比方法3：Dabov等在文献“Image denoising by sparse3-D transform-domaincollaborative filtering.IEEE Transactions on Image Processing，2007,vol.16no.8,pp.2080-2095.”中提出的基于三维块匹配去噪方法，简称BM3D。

仿真1，利用所述的三种现有方法和本发明方法对图2中的四幅图分别在噪声方差为15、20、25和35四种情况下进行了去噪仿真。四幅图像在四种情况下的峰值信噪比结果列于下表1中，其中sigma为噪声方差，PSNR的数据最大值用粗体突出显示。

表1不同的去噪方法对不同图像在不同噪声水平Sigma下的峰值信噪比PSNR

从表1的四幅图像在四种噪声水平下的去噪后的PSNR可以看出，本发明方法在四种噪声水平下的PSNR均优于其他三种方法。从单幅图像的相同噪声水平下，本发明方法的PSNR都比其它方法高，尤其是Barbara在噪声水平为35的情况下，本发明方法的PSNR较NLM高出2.27db，较KSVD高出1.5db，较BM3D高出0.22db。从平均值上来看，本发明方法的PSNR较NLM高出1.76db，较KSVD高出1.15db，较BM3D高出0.1db，所以，可以充分地验证本发明方法较其它三种方法对自像去噪的更具有有效性。

仿真2，利用所述的三种现有方法和本发明方法对图3进行去噪仿真，结果如图4所示，其中，图4(a)是利用NLM方法对图3进行去噪的结果图，图4(b)是利用KSVD方法对图3进行去噪的结果图，图4(c)是利用BM3D方法对图3进行去噪的结果图，图4(d)是本发明方法对图3进行去噪的结果图。

从图4可以看出，本发明的图像边缘和细节保持能力优于NLM方法，图像轮廓纹理保持效果优于KSVD方法和BM3D方法。

综上所述，本发明的去噪效果优于现有的三种去噪方法。

Claims (4)

1.一种结合局部和非局部的自适应图像去噪方法，包括如下步骤：

（1）输入一幅N行M列的含噪图像Y，设置最大迭次数γ和停止参数δ；

（2）采用下式估计含噪图像Y的噪声标准差σ_n：

${\mathrm{\&sigma;}}_n=\frac{\mathrm{median}\left(\mathrm{abs}\right|W\left|\right)}{0.6745},$

其中，W为含噪图像Y进行小波分解得到的第一层高频系数，abs|·|是取绝对值操作，median(·)是取中值操作；

（3）以含噪图像Y中的任意像素点为中心，分别确定一个l×l大小的位置区域Γ_i和一个s×s大小的搜索窗Ω_i，提取含噪图像Y的位置区域Γ_i的像素得到图像块y_i，其中i∈{1,...,N×M}，l、s的取值范围分别为5~11和21~41奇数个像素；

（4）将l×l大小的图像块y_i按列拉成l²×1大小的像素矢量v_i，按照下式分别计算像素矢量v_i和搜索窗Ω_i内所有像素矢量集合{v_ip}之间的相似度，得到相似度集合{w_ip}：

$w_{\mathrm{ip}}=\exp \left(\frac{-{\left|\right|v_i-v_{\mathrm{ip}}\left|\right|}_2^2}{h}\right),$

其中，p=1,...,s2，

表示二范数，h为衰减因子，h的取值范围10σ_n~15σ_n；

（5）对相似度集合{w_ip}中的s²个相似度按照从大到小排序，选取前k个相似度作为像素矢量v_i的相似度集合{w_ij}，将该相似度集合{w_ij}对应的像素矢量作为像素矢量v_i的相似像素矢量集合{v_ij}，其中j=1，...,k，k的取值范围为10~30；

（6）计算相似像素矢量集合{v_ij}的像素矢量非局部均值

${\stackrel{\_}{v}}_i=\frac{\underset{j=1}{\overset{k}{\mathrm{\&Sigma;}}}{w}_{\mathrm{ij}}{v}_{\mathrm{ij}}}{\underset{j=1}{\overset{k}{\mathrm{\&Sigma;}}}{w}_{\mathrm{ij}}};$

（7）在含噪图像Y中，对所有的像素点分别进行步骤（3）至步骤（6）操作，得到N×M个像素矢量v_i，N×M个相似像素矢量集合{v_ij}，N×M个像素矢量非局部均值

（8）采用K均值k-means聚类方法将N×M个像素矢量v_i分成f类，得到f个像素矢量子集U_t和f个聚类中心c_t，其中t∈{1,...,f}，f的取值范围为15~60；

（9）采用主成分分析PCA方法分别对f个像素矢量子集U_t训练一个l²×l²大小的字典D_t；

（10）利用f个字典D_t对N×M个像素矢量v_i分别进行自适应去噪，得到N×M个去噪像素矢量

（11）将N×M个l²×1大小的去噪像素矢量分别按列拉成N×M个l×l大小的去噪图像块

（12）采用基于三维块匹配去噪BM3D的图像块聚集方法对N×M个去噪图像块

进行聚集，得到去噪图像X；

（13）计算含噪图像Y和去噪图像X之间的均方误差e(X,Y)；

（14）判断均方误差e(X,Y)是否小于停止参数δ或者本次迭代次数是否达到最大迭代次数γ，若满足两个条件中的一个，则终止迭代，输出去噪图像X，否则将去噪图像X作为含噪图像Y，重复步骤（2）到步骤（13）。

2.根据权利要求1所述的图像去噪方法，其中步骤（10）所述的利用f个字典D_t对N×M个像素矢量v_i分别进行自适应去噪，按如下步骤进行：

10a）分别计算像素矢量v_i与f个聚类中心c_t之间的欧式距离d(v_i,c_t)：

$d(v_i,c_t)={\left|\right|v_i-c_t\left|\right|}_2^2,$

其中，i∈{1，...,N×M}，t=1，...,f，f为聚类的类别数，

表示二范数；

10b）在f个欧式距离d(v_i,c_t)中选取欧式距离d(v_i,c_t)最小值时c_t对应的下角标g，其中g∈{1，...,f}；

10c）从步骤（9）中得到的f个PCA字典D_t中选取下角标为g的字典，得到稀疏字典D_g；

10d）利用上述稀疏字典D_g以及步骤（7）中得到的像素矢量非局部均值

和相似像素矢量集合{v_ij}，计算自适应参数矢量λ：

$\mathrm{\&lambda;}=\frac{{\left[\frac{1}{k}\underset{j=1}{\overset{k}{\mathrm{\&Sigma;}}}\right[\mathrm{diag}\left(D_g^T({\stackrel{\_}{v}}_i-v_{\mathrm{ij}})\right){]}^2-{\mathrm{\&sigma;}}_n^{2}I]}_{+}}{{\mathrm{\&sigma;}}_n^{2}I},$

其中，j=1，...,k，k为相似像素矢量的个数，I为单位矩阵，σ_n为步骤（2）中计算得到的图像的噪声标准差，diag(·)为将矢量转成对角化矩阵的操作，

表示稀疏字典D_g的转置，[·]₊为矩阵取正运算；

10e）利用稀疏字典D_g和像素矢量均值计算像素矢量v_i的去噪像素矢量

${\hat{v}}_i=D_g\frac{\left(D_g^T{\stackrel{\_}{v}}_i\right)+\mathrm{\&lambda;}*\left(D_g^T{v}_i\right)}{I+\mathrm{\&lambda;}};$

10f）对N×M个像素矢量分别重复步骤10a)至步骤10e)，得到N×M个去噪像素矢量。

3.根据权利要求1所述的图像去噪方法，其中步骤（12）所述的采用基于三维块匹配去噪BM3D的图像块聚集方法对N×M个去噪图像块进行聚集，按如下步骤进行：

12a）初始化两个大小为N×M的零矩阵，得到标记矩阵R和重构矩阵Q；

12b）采用形状参数为2的凯瑟Kaiser窗，初始化一个大小为l×l的聚集矩阵Θ，其中l为图像块的边长；

12c）采用下式分别对N×M个去噪图像块

进行聚集加权，得到N×M个加权去噪图像块z_i：

$z_i=<\mathrm{\&Theta;},{\hat{y}}_i>,$

其中，i∈{1，...,N×M}，<·,·>表示矩阵的内积运算；

12d）根据步骤（3）中得到的图像块y_i所在含噪图像Y中的位置区域Γ_i，分别在标记矩阵R和重构矩阵Q的相同位置确定一个相同大小的区域，得到标记区域Φ_i和重构区域ψ_i；

12e）对标记矩阵R中区域Φ_i的元素值加上聚集矩阵Θ的元素值，得到标记更新矩阵

12f）对重构矩阵Q中区域ψ_i的元素值加上加权去噪图像块z_i的元素值，得到重构更新矩阵

12g）对标记矩阵R的所有位置进行步骤12d）至步骤12f）操作，得到最终标记更新矩阵

和最终重构更新矩阵

12h）利用上述的最终标记更新矩阵

和最终重构更新矩阵

计算去噪图像X：

$X=\mathrm{div}(\hat{Q},\hat{R}),$

其中，div(·,·)表示矩阵的点除运算，即

表示矩阵

中的每个元素分别除以矩阵相应位置上的元素。

4.根据权利要求1所述的图像去噪方法，其中步骤（13）所述的计算含噪图像Y和去噪图像X之间的均方误差e(X,Y)，按如下公式进行：

$e(X,Y)=\frac{1}{256*\mathrm{NM}}{\left|\right|X-Y\left|\right|}_2^2,$

其中，N、M表示含噪图像Y的行数、列数，

表示二范数。

Images