一种基于独立分量分析和奇异值分解的非局部TV图像去噪方法，为了减少噪声对像素相似度权重函数的干扰，提高像素相似度权重赋值的准确性，构建了一种基于待去噪图像中各领域图像块最大奇异值的新型像素相似度权重函数，用于NLTV模型，从而得到新的NLTV模型。通过新的NLTV方法对输入的唯一含噪图像u0进行初步去噪，得到另一含噪图像u1，把初步去噪图像u1认为是另一幅含噪输入图像，在获得的u1、u0基础上应用ICA方法对输入图像u0去噪，得到ICA方法去噪后的图像u2，提高了ICA在图像去噪领域的应用价值。为了得到更好的去噪效果，对图像u2再一次应用新的NLTV方法去噪，得到最后去噪图像ufinal。
1.基于独立分量分析和奇异值分解的非局部TV图像去噪方法，步骤如下：
(2)设置方法中的相关参数，包括ICA迭代收敛阈值e，第一轮NLTV的搜索窗口大小N1×N1、邻域窗口大小N2×N2、保真参数λ1、像素相似度权重函数ω1的参数h1和j1、像素相似度权重函数ω1中高斯核的标准差σ1，分裂的Bregman迭代辅助变量b1k的初始值、平滑参数θ1；第二轮NLTV的搜索窗口大小N3×N3、邻域窗口大小N4×N4、保真参数λ2、像素相似度权重函数ω2的参数h2和j2、像素相似度权重函数ω2中高斯核的标准差σ2，分裂的Bregman迭代辅助变k量b2的初始值、平滑参数θ2；
(3)设Mx是以步骤(1)输入的含噪图像u0中像素点x∈Ω为中心的大小为N2×N2的邻域图像块的像素灰度值矩阵，Ω为u0的图像空间；对图像块Mx进行奇异值分解：Mx＝UxΛxVxT；式中Ux、Vx分别是Mx的左奇异矩阵和右奇异矩阵，大小都是N2×N2；Λx是Mx的奇异值矩阵，大小也是N2×N2，Λx的对角元素不为零，其他元素全为零，Λx的对角元素一共有N2个：按从大到小排列为： Λx的对角元素就是Mx的奇
异值，包含了Mx的全部特征；奇异值的定义：对于m×n阶矩阵C，CTC的n个特征值的非负平方根叫作C的奇异值；
(4)利用步骤(3)得到的图像块的奇异值构建新的像素相似度权重函数；由于奇异值包含了图像块的主要特征，所以相似的图像块之间的奇异值是相近的；在图像块Mx的奇异值中，大的奇异值包含了图像块的主要特征，小的奇异值包含了图像块的次要特征；并且在含有噪声的图像中，噪声不是图像的主要特征，所以在构建新型相似度权重函数时，只用图像块的最大奇异值，通过图像块的主要特征来判断图像块之间的相似性，这样就排除了噪声的干扰；构建图像u0中两个像素点p1和q1的相似度权重函数：其中p1为当前像素点，q1是以p1为中心的搜索窗口内的一点，ap1是以p1为中心的、大小为N2×N2的图像块的像素灰度值矩阵的最大奇异值，aq1是以q1为中心的、大小为N2×N2的图像块的像素灰度值矩阵的最大奇异值； 表示以p1为中心的大小为N2×N2的图像块与以q1为中心的大小为N2×N2的图像块之间的高斯加权距离， 表示求和范围是以p1或q1为中心的的大小为N2×N2的邻域内的每一像素点，不包括p1或q1自身，共N2×N2-1项求和项， 是标准差为σ1的高斯核，h1和j1是常数，h1、j1通过干预指数函数的衰减速度来控制权重函数的大小,h1和j1的取值越大，权重函数的值越接近1，算法收敛速度快，但是难以达到最优值，h1和j1的取值越小，权重函数的值越接近0，算法经过多次迭代可以收敛到最优值，但是耗费时间多，h1和j1的取值原则上要综合以上两点因素，取大小合适的值；
(5)为了应用ICA方法进行图像去噪，除了输入的含噪图像u0，还需另一幅含噪图像u1，通过用NLTV方法对含噪图像u0进行初步去噪来获得u1，称为第一轮NLTV去噪；具体见步骤(5)至(6)；首先建立基于步骤(4)得到的相似度权重函数的新的NLTV模型，称为第一轮新NLTV模型： 其中， 是第一轮NLTV模型的目标函数，λ1是保真参数，Ω是u0和u1的图像空间，u0是输入的含噪图像，u1是去噪后的得到的图像， 是GUY GILBOA和STANLEY OSHER提出的非局部梯度算子，其中变量
p1表示当前像素点，变量q1表示以p1为中心的大小为N1×N1的搜索窗口内的一点，u1(p1),u1(q1)分别是像素点p1、q1的灰度值；本步骤中使用的像素相似度权重函数ω1(p1,q1)在步骤(4)中建立；
(6)对于步骤(5)建立的模型，采用分裂的Bregman迭代算法进行数值实现；引入辅助函数w1k和辅助变量b1k，构造如下迭代格式：其中，k的取值是0，1，2，…,等非负整数，迭代初始值u10＝u0，b1k、w1k分别表示分裂的Bregman迭代的辅助变量和辅助函数，λ1就是步骤(5)建立的NLTV模型中的保真参数，θ1是控制迭代结果的平滑参数，b1k的初始值b10、以及λ1和θ1的赋值在步骤(2)中进行预设；
求解式(1)和式(2)，并数值化，式(3)也数值化，由此得到数值化后的迭代格式如式(4)、式(5)和式(6)所示；
其中， 表示求和范围是以p1为中心的搜索窗口N1×N1内除p1以外的每一像素点，求和项数共N1×N1-1项，ω1(p1,q1)是u0的相似度权重函数；
在本步骤中，设初值k＝0，顺序地按式(4)、式(5)和式(6)进行迭代运算一次，得到初步去噪后的图像u1＝u11；
(7)使用u0与u1构建混合矩阵S，对S进行中心化后得到矩阵Se，对Se进行白化处理，得到处理结果Z；具体方法是：首先将大小N×N的图像u0转换成1×N2的行向量X，转换规则是u0的第α(α＝1,2,3,...,N)列第β(β＝1,2,3,...,N)行元素转为X的第(α-1)×N+β列元素；用同样方法将大小N×N的图像u1转换成1×N2的行向量Y；用X和Y构建混合矩阵 并进行中心化和白化；中心化是将将原始数据减去平均数；白化也称为球化，它的本质是去相关；
如果均值为零的随机向量O＝[o1,o2,...,on]T满足E(OOT)＝I，其中I是单位矩阵，那么随机向量O＝[o1,o2,...,on]T是白化向量，公式如下：对S进行中心化后得到矩阵Se：
其中xi和yi分别是X和Y的第i个元素，i＝1,2,…,N×N， 是X的所有元素的均值， 是Y的所有元素的均值，对Se进行白化处理，即白化矩阵W0与Se相乘，得到Z：
其中白化矩阵W0＝Λ-1/2UT，Λ是SeT的协方差矩阵的特征值矩阵，U是SeT的协方差矩阵的特征向量矩阵；Z将参与后续的步骤；这一处理可以降低后续步骤的计算复杂度；
(9)建立一个2×1的随机列向量L，L各元素取值范围是[0～1]；
(10)对L进行迭代，目的是通过L的迭代算法的运算结果来逼近解混矩阵B的一个列向量的真值；迭代公式如下：L＝E{Zg(LTZ)}-E{g1(LTZ)}L
其中E{·}是均值运算；g(·)为任意二次函数；令g(·)为g(x)＝x2；g1(·)是g(·)的一阶导数；
(11)如果L满足||LTL|-1|＜e，则转到步骤(12)；如果L不满足||LTL|-1|＜e，则返回到步骤(10)；其中e是ICA方法中的收敛阈值，为常数，其值在步骤(2)预先设置；
(12)用L来替换更新B中的一个列向量；如果B的所有列向量都被替换更新一次，则转到步骤(13)；如果B的所有列向量没有被替换更新完，则返回到步骤(9)；
(13)利用得到的解混矩阵B可以分离出噪声，获得去噪后的图像；步骤如下：先计算解混结果W＝BTS；W是解混后得到的结果，再将W中的两个1×N2大小行向量都转换成N×N大小的矩阵，则可以得到分离后的去噪图像u2和噪声图像；转换规则是：W的第一个行向量的第(α-1)×N+β列元素转换为u2的第α列第β行元素，其中α和β取值均为1,2,3,...,N；
(14)按照步骤(3)的方法获得图像u2中以图像各点为中心大小为N4×N4的各图像块的奇异值和最大，和步骤(3)不同的是本步骤中每个图像块的奇异值共有N4个；
(15)基于步骤(14)得到的图像块的奇异值，构建u2的两个像素点p2和q2的相似度权重函数：其中p2为当前像素点，q2是以p2为中心的搜索窗口内的一点，ap2是以p2为中心的、大小为N4×N4的图像块的像素灰度值矩阵的最大奇异值，aq2是以q2为中心的、大小为N2×N2的图像块的像素灰度值矩阵的最大奇异值；
表示以p2为中心的大小为N4×N4的图像块与以q2为中
心的大小为N4×N4的图像块之间的高斯加权距离， 表示求和范围是以p2或q2为中心的大小为N4×N4的邻域内的每一像素点，不包括p2或q2自身，共N4×N4-1项求和项， 是标准差为σ2的高斯核；h2和j2是常数，h2、j2通过干预指数函数的衰减速度来控制权重函数的大小,h2和j2的取值越大，权重函数的值越接近1，算法收敛速度快，但是难以达到最优值，h2和j2的取值越小，权重函数的值越接近0，算法经过多次迭代可以收敛到最优值，但是耗费时间多，h2和j2的取值原则上要综合以上两点因素，取大小合适的值；
(16)为了再次提升去噪效果，对图像u2进行第二轮NLTV去噪；首先建立NLTV模型：其中， 是第二轮NLTV模型的目标函数，λ2是
保真参数，Ω是u2和u3的图像空间，u3是去噪后的图像， 是GUY GILBOA和STANLEY OSHER提出的非局部梯度算子，其中p2为第二轮NLTV去噪的当前像素点，q2是以p2为中心的搜索窗口内的一点，u3(p2)、u3(q2)分别是图像上的点p2和q2的像素灰度值；本步骤中使用的权重函数ω2(p2,q2)已在步骤(15)中建立；
(17)对步骤(16)建立的模型，采用分裂的Bregman迭代方法逼近求解，迭代过程分成如下三步：其中，k的取值是0，1，2，…,等非负整数，迭代初始值u30＝u2，b2k、w2k分别表示分裂的Bregman迭代的辅助变量和辅助函数，λ2就是步骤(16)建立的非局部TV模型中的保真参数，θ2是控制迭代结果的平滑参数，b2k的初始值b20、以及λ2和θ2的赋值在步骤(2)中进行预设；
求解式(7)和式(8)，并数值化，式(9)也数值化，由此得到数值化后的迭代格式，如式(10)、式(11)和式(12)所示；第一次迭代运算前，令k＝0；
(18)顺序地应用公式(10)、(11)、(12)进行迭代运算，
其中ω2(p2,q2)是图像u2中像素p2和q2间的相似度权重函数， 表示求和范围是以p2为中心的搜索窗口N3×N3内除p2以外的每一像素点，求和项数共N3×N3-1项；
(19)计算步骤(18)中的本次迭代输出图像u3k+1的峰值信噪比PSNR，如果本次迭代输出图像u3k+1的峰值信噪比PSNR小于等于上一次迭代输出图像u3k的峰值信噪比PSNR，则将上一k k次迭代输出图像u3作为最优值ufinal输出，即令ufinal＝u3 ，并转到步骤(20)；如果本次迭代后输出图像的峰值信噪比PSNR大于上一次迭代输出图像的峰值信噪比PSNR，则令k＝k+1，并返回到步骤(18)，继续迭代运算；
基于独立分量分析和奇异值分解的非局部TV图像去噪方法
技术领域
[0001] 本发明属于图像处理技术领域，涉及了去除加性噪声的图像去噪领域。
背景技术
[0002] 图像去噪旨在通过对被噪声污染的图像进行某种处理，以降低噪声对原始有用信息的影响，尽可能的还原出被噪声污染前的原始图像。
[0003] GUY GILBOA和STANLEY OSHER提出的非局部TV模型(Nonlocal Total Variation,NLTV)是将GUY GILBOA和STANLEY OSHER提出的非局部算子(见文献：NONLOCAL OPERATORS WITH APPLICATIONS TO IMAGE PROCESSING.SIAM Multiscale Modeling and Simulation.Vol.7,No.3,pp.1005–1028)引入到Rudin-Osher-Fatemi提出的总变分(Total Variation，TV)模型中，该模型包含保真项和正则项。保真项表示观测图像(即待去噪图像)与去噪后图像的接近程度，正则项表达去噪后图像光滑程度的先验知识。NLTV模型具有在去除图像噪声的同时能较好地保留图像的纹理细节的特点。
[0004] 在正则项中，相似度权重函数具有很重要的作用。其表达式与含噪图像的图像块中各像素灰度值直接相关，。但是由于噪声的存在，噪声会对像素值产生干扰，进而影响像素相似度权重。奇异值分解是一种基于特征向量的矩阵变换方法，在信号处理、模式识别、数字水印技术等方面都得到了应用。由于图像具有矩阵结构，因此将其运用在图像处理中，用来提取图像的主要特征，本发明用奇异值分解方法提取待去噪图像中各领域图像块的奇异值，基于提取的图像块奇异值，构建新的像素相似度权重函数，以降低噪声对像素相似度权重的影响，再将构建的新的像素相似度权重函数应用到NLTV模型，从而得到新的NLTV模型。
[0005] 独立分量分析(independent component analysis，ICA)是20世纪90年代发展起来的一种新的信号处理技术。基本的ICA是指从多个源信号的线性混合信号中分离出源信号的技术。除了需要已知源信号是统计独立外，无其他先验知识，ICA是伴随着盲信源问题而发展起来的，故又称盲分离。在时间信号处理方面由于ICA可以在n个原信号和线性混合系数都未知的情况下可以把这n个原信号从他们线性混合成的n个混合信号中分离出来，并且它的假设条件仅是这n个原信号之间是相互独立的，所以ICA在时间信号分离方面具有良好的应用及优势。
[0006] ICA作为盲信号分离的一种有效工具，在假设各原信号相互独立的条件下可以将各原信号从它们的混合信号中分离出来。图像可以看成是一种信号，含噪声的图像是噪声与无噪声图像这两种独立信号的混合信号，这符合ICA的假设条件，所以ICA可以将噪声图像分离成噪声与无噪声图像这两种原信号，达到图像去噪的目的。但是ICA应用于图像去噪需要至少两幅含噪输入图像，在只有一幅含噪输入图像时，ICA无法应用于图像去噪，本发明方法解决ICA在图像去噪领域应用的这一不足。
发明内容
[0007] 本发明要克服现有技术的噪声对像素相似度权重函数的干扰大、像素相似度权重赋值的准确性低的缺点，提出一种基于独立分量分析和奇异值分解的非局部TV图像去噪方法。
[0008] 本发明构建了一种基于待去噪图像中各领域图像块最大奇异值的新型像素相似度权重函数，用于NLTV模型，从而得到新的NLTV模型。通过新的NLTV方法对输入的唯一含噪图像u0进行初步去噪，得到另一含噪图像u1，把初步去噪图像u1认为是另一幅含噪输入图像，在获得的u1、u0基础上应用ICA方法对输入图像u0去噪，得到ICA方法去噪后的图像u2，提高了ICA在图像去噪领域的应用价值。为了得到更好的去噪效果，对图像u2再一次应用新的NLTV方法去噪，得到本发明方法的最后去噪图像ufinal。
[0009] 本文中所提到的NLTV模型是指GUY GILBOA和STANLEY OSHER提出的非局部TV模型。本文中所提到的新的NLTV模型是指NLTV模型中的像素相似度权重函数用本发明构建的权重函数代替后的改进模型。
[0010] 本发明所述的基于独立分量分析和奇异值分解的非局部TV图像去噪方法图，步骤如下：
[0011] (1)首先输入大小为N×N的含噪图像u0。
[0012] (2)设置算法中的相关参数，包括ICA方法中迭代收敛阈值e，第一轮NLTV的搜索窗口大小N1×N1、邻域窗口大小N2×N2、保真参数λ1、像素相似度权重函数ω1的参数h1和j1、像素相似度权重函数ω1中高斯核的标准差σ1，分裂的Bregman迭代辅助变量b1k的初始值b10、平滑参数θ1。第二轮NLTV的搜索窗口大小N3×N3、邻域窗口大小N4×N4、保真参数λ2、像素相似度权重函数ω2的参数h2和j2、像素相似度权重函数ω2中高斯核的标准差σ2，分裂的Bregman迭代辅助变量b2k的初始值b20、平滑参数θ2。
[0013] (3)通过奇异值分解方法获得图像u0中各邻域图像块(图像块大小为N2×N2)的奇异值矩阵和最大奇异值。
[0014] (4)构建基于步骤(3)得到的最大奇异值的图像u0的像素间相似度权重函数。
[0015] (5)应用步骤(4)构建的权重函数，建立用于对u0进行第一轮NLTV方法初步去噪的新的非局部TV模型。
[0016] (6)对步骤(5)建立的新的非局部TV模型对噪声图像u0进行NLTV方法第一轮去噪，得到去噪图像u1，用于以下的ICA方法去噪[步骤(7)-(13)]。
[0017] (7)将u0与u1看成是纯噪图像与无噪图像线性混合成的两个信号，基于这两个图像信号构建混合矩阵S，对S进行中心化后得到矩阵Se，对Se进行白化处理，得到处理结果Z。
[0018] (8)建立初始解混矩阵B。
[0019] (9)建立随机列向量L，L各向量元素取值范围是[0,1]。
[0020] (10)对L进行迭代。目的是通过L的迭代运算结果来逼近解混矩阵B的一个列向量的真值。
[0021] (11)若满足L迭代的停止条件，则转向步骤(12)，否则返回步骤(10)。
[0022] (12)用迭代结果L来替换更新B中的一个列向量。如果B的所有列向量都被替换更新一次，则转到步骤(13)。如果B的所有列向量没有被替换更新完，则返回到步骤(9)；
[0023] (13)通过W＝BTS运算，分离出噪声，从而得到去噪后的图像u2。W的第一行行向量经变换即可得到去噪后的图像u2，转换规则是：W的第一个行向量的第(α-1)×N+β列元素转换为u2的第α列第β行元素，其中α和β取值均为1,2,3,...,N。
[0024] (14)为了对u2进一步去噪，进行第二轮NLTV去噪，首先在本步骤，通过奇异值分解方法获得图像u2中各领域图像块(图像块大小为N4×N4)的奇异值矩阵和最大奇异值。
[0025] (15)构建基于步骤(14)得到的最大奇异值的图像u2中二个像素间相似度权重函数。
[0026] (16)应用步骤(15)构建的权重函数，建立用于对u2进行第二轮NLTV方法进一步去噪的新的非局部TV模型。
[0027] (17)对步骤(16)建立的新的非局部TV模型,用分裂的Bregman算法求解，得到三步式数值迭代计算公式。并设迭代计数变量初值k＝0。
[0028] (18)用步骤(17)中得到的分裂的Bregman算法三步式数值迭代计算公式进行顺序k+1迭代运算，得到本次迭代的输出图像u3 。
[0029] (19)计算步骤(18)中的本次迭代输出图像u3k+1的峰值信噪比，如果本次迭代输出图像u3k+1的峰值信噪比小于等于上一次迭代输出图像u3k的峰值信噪比，则将上一次迭代输出图像u3k作为最优值ufinal输出，即令ufinal＝u3k，并转到步骤(20)；如果本次迭代后输出图像的峰值信噪比大于上一次迭代输出图像的峰值信噪比PSNR，则令k＝k+1，并返回到步骤(18)，继续迭代运算。
[0030] (20)将ufinal作为最后去噪结果图像输出。
[0031] 本发明与现有的技术相比具有以下优点：
[0032] (1)建立新型的权重函数。通过奇异值分解得到每个领域图像块的奇异值矩阵和最大奇异值，选取包含图像块主要信息的最大奇异值来构建新的像素相似度权重函数，并将该权重函数用于NLTV模型进行图像去噪。减少了噪声对于权重分配干扰，提高了权重分配的准确性。
[0033] (2)通过对输入的含噪图像u0初步去噪，得到另一幅含噪图像u1，再应用独立分量分析方法进行图像去噪，克服了独立分量分析方法在只能提供一幅含噪图像情况下不能用于图像去噪的局限性。
附图说明
[0034] 图1是本发明方法的流程图。
[0035] 图2是本发明仿真实验所用的原始图像，图2a是Lena图像，图2b是Peppers图像。
[0036] 图3是本发明对原始Peppers图像加均值为零标准差为15的高斯噪声后的图像进行去噪的仿真结果，其中，图3a是待去噪的含噪图像(所含噪声的标准差为15)，图3b是NLTV去噪结果，图3c是本发明方法的去噪结果。
具体实施方式
[0037] 下面对本发明做进一步说明。
[0038] 本发明方法所述的基于独立分量分析和奇异值分解的非局部TV图像去噪方法图，步骤如下：
[0039] (1)首先输入N×N大小的含噪图像u0；
[0040] (2)设置算法中的相关参数，包括ICA迭代收敛阈值e，第一轮NLTV的搜索窗口大小N1×N1、邻域窗口大小N2×N2、保真参数λ1、像素相似度权重函数ω1的参数h1和j1、像素相似度权重函数ω1中高斯核的标准差σ1，分裂的Bregman迭代辅助变量b1k的初始值、平滑参数θ1。第二轮NLTV的搜索窗口大小N3×N3、邻域窗口大小N4×N4、保真参数λ2、像素相似度权重函数ω2的参数h2和j2、像素相似度权重函数ω2中高斯核的标准差σ2，分裂的Bregman迭代辅助变量b2k的初始值、平滑参数θ2；
[0041] (3)设Mx是以步骤(1)输入的含噪图像u0中像素点x∈Ω为中心的大小为N2×N2的T图像块的像素灰度值矩阵，Ω为u0的图像空间。对图像块Mx进行奇异值分解：Mx＝UxΛxVx 。
式中Ux、Vx分别是Mx的左奇异矩阵和右奇异矩阵，大小都是N2×N2。Λx是Mx的奇异值矩阵，大小也是N2×N2，Λx的对角元素不为零，其他元素全为零，Λx的对角元素一共有N2个：
按照从大到小的顺序排列： Λx的对角元素就是
Mx的奇异值，包含了Mx的全部特征。奇异值的定义：对于m×n阶矩阵C，CTC的n个特征值的非负平方根叫作C的奇异值；
[0042] (4)利用步骤(3)得到的图像块的奇异值构建新的像素相似度权重函数；由于奇异值包含了图像块的主要特征，所以相似的图像块之间的奇异值是相近的；在图像块Mx的奇异值ax1、ax2、…、 中，大的奇异值包含了图像块的主要特征，小的奇异值包含了图像块的次要特征；并且在含有噪声的图像中，噪声不是图像的主要特征，所以本发明在构建新型相似度权重函数时，只用图像块的最大奇异值，通过图像块的主要特征来判断图像块之间的相似性，这样就排除了噪声的干扰；构建图像u0中两个像素点p1和q1的相似度权重函数：
[0043]
[0044] 其中p1为当前像素点，q1是以p1为中心的搜索窗口内的一点，ap1是以p1为中心的、大小为N2×N2的图像块的像素灰度值矩阵的最大奇异值，aq1是以q1为中心的、大小为N2×N2的图像块的像素灰度值矩阵的最大奇异值； 表示以p1为中心的大小为N2×N2的图像块与以q1为中心的大小为N2×N2的图像块之间的高斯加权距离，表示求和范围是以p1或q1为中心的的大小为N2×N2的邻域内的每一像素点(不包括p1或q1自身)，共N2×N2-1项求和项， 是标准差为σ1的高斯核，h1和j1是常数，h1、j1通过干预指数函数的衰减速度来控制权重函数的大小,h1和j1的取值越大，权重函数的值越接近1，算法收敛速度快，但是难以达到最优值，h1和j1的取值越小，权重函数的值越接近0，算法经过多次迭代可以收敛到最优值，但是耗费时间多，h1和j1的取值原则上要综合以上两点因素，取大小合适的值；
[0045] (5)本发明中为了应用ICA方法进行图像去噪，除了输入的含噪图像u0，还需另一幅含噪图像u1，本发明通过用NLTV方法对含噪图像u0进行初步去噪来获得u1(称为第一轮NLTV去噪)。具体见步骤(5)至(6)。首先建立基于步骤(4)得到的相似度权重函数的新的NLTV模型(称为第一轮NLTV模型)： 其中，是第一轮NLTV模型的目标函数，λ1是保真参数，Ω是u0和u1的图像空间，u0是输入的含噪图像，u1是去噪后的得到的图像， 是GUY GILBOA和STANLEY OSHER提出的非局部梯度算子， 其
中变量p1表示当前像素点，变量q1表示以p1为中心的大小为N1×N1的搜索窗口内的一点，u1(p1),u1(q1)分别是像素点p1、q1的灰度值。本步骤中使用的像素相似度权重函数ω1(p1,q1)在步骤(4)中建立；
[0046] (6)对于步骤(5)建立的模型，采用分裂的Bregman迭代算法进行数值实现；引入辅助函数w1k和辅助变量b1k，构造如下迭代格式：
[0047]
[0048]
[0049]
[0050] 其中，k的取值是0，1，2，…,等非负整数，迭代初始值u10＝u0，b1k、w1k分别表示分裂的Bregman迭代的辅助变量和辅助函数，λ1就是步骤(5)建立的NLTV模型中的保真参数，θ1是控制迭代结果的平滑参数，b1k的初始值b10、以及λ1和θ1的赋值在步骤(2)中进行预设。
[0051] 求解式(1)和式(2)，并数值化，式(3)也数值化，由此得到数值化后的迭代格式如式(4)、式(5)和式(6)所示。
[0052]
[0053]
[0054]
[0055] 其中， 表示求和范围是以p1为中心的搜索窗口N1×N1内除p1以外的每一像素点，求和项数共N1×N1-1项，ω1(p1,q1)是u0的相似度权重函数。
[0056] 在本步骤中，设初值k＝0，顺序地按式(4)、式(5)和式(6)进行迭代运算一次，得到1
初步去噪后的图像u1＝u1。
[0057] (7)使用u0与u1进行中心化和白化处理。具体方法是：首先将大小N×N的图像u0转换成1×N2的行向量X，转换规则是u0的第α(α＝1,2,3,...,N)列第β(β＝1,2,3,...,N)行元素转为X的第(α-1)×N+β列元素；用同样方法将大小N×N的图像u1转换成1×N2的行向量Y。用X和Y构建混合矩阵 并进行中心化和白化。中心化是将将原始数据减去平均数。
白化也称为球化，它的本质是去相关。如果均值为零的随机向量O＝[o1,o2,...,on]T满足E(OOT)＝I，其中I是单位矩阵，那么随机向量O＝[o1,o2,...,on]T是白化向量。公式如下：
[0058] 对S进行中心化后得到矩阵Se：
[0059]
[0060] 其中xi和yi分别是X和Y的第i个元素，i＝1,2,…,N×N， 是X的所有元素的均值，是Y的所有元素的均值，
[0061] 对Se进行白化处理，即白化矩阵W0与Se相乘，得到Z：
[0062] Z＝W0Se
[0063] 其中白化矩阵W0＝Λ-1/2UT，Λ是SeT的协方差矩阵的特征值矩阵，U是SeT的协方差矩阵的特征向量矩阵。Z将参与后续的步骤。这一处理可以降低后续步骤的计算复杂度。
[0064] (8)设置解混矩阵B的初始值为2×2的零矩阵；
[0065] (9)建立一个大小为2×1，的随机列向量L，L各元素取值范围是[0～1]；
[0066] (10)对L进行迭代，目的是通过L的迭代算法的运算结果来逼近解混矩阵B的一个列向量的真值。迭代公式如下：
[0067] L＝E{Zg(LTZ)}-E{g1(LTZ)}L
[0068] L＝L-BBTL；
[0069] L＝L/||L||；
[0070] 其中E{·}是均值运算，g(·)为任意二次函数。在本发明中令g(·)为g(x)＝x2。g1(·)是g(·)的一阶导数；
[0071] (11)如果L满足||LTL|-1|＜e，则转到步骤(12)；如果L不满足||LTL|-1|＜e，则返回到步骤(10)。其中e是ICA方法中的收敛阈值，为常数，其值在步骤(2)预先设置；
[0072] (12)用L来替换更新B中的一个列向量。如果B的所有列向量都被替换更新一次，则转到步骤(13)。如果B的所有列向量没有被替换更新完，则返回到步骤(9)；
[0073] (13)利用得到的解混矩阵B可以分离出噪声，获得去噪后的图像。步骤如下：先计T 2算解混结果W＝BS。W是解混后得到的结果，再将W中的两个1×N大小行向量都转换成N×N大小的矩阵，则可以得到分离后的去噪图像u2和噪声图像。转换规则是：W的第一个行向量的第(α-1)×N+β列元素转换为u2的第α列第β行元素，其中α和β取值均为1,2,3,...,N；
[0074] (14)按照步骤(3)的方法得到图像u2中以图像各点为中心大小为N4×N4的各图像块的奇异值，和步骤(3)不同的是本步骤中每个图像块的奇异值共有N4个；
[0075] (15)基于步骤(14)得到的图像块的奇异值，构建u2的两个像素点p2和q2的相似度权重函数：
[0076]
[0077] 其中p2为当前像素点，q2是以p2为中心的搜索窗口内的一点，ap2是以p2为中心的、大小为N4×N4的图像块的像素灰度值矩阵的最大奇异值，aq2是以q2为中心的、大小为N2×N2的图像块的像素灰度值矩阵的最大奇异值； 表示以p2为中心的大小为N4×N4的图像块与以q2为中心的大小为N4×N4的图像块之间的高斯加权距离，表示求和范围是以p2或q2为中心的大小为N4×N4的邻域内的每一像素点(不包括p2或q2自身)，共N4×N4-1项求和项， 是标准差为σ2的高斯核；h2和j2是常数，h2、j2通过干预指数函数的衰减速度来控制权重函数的大小,h2和j2的取值越大，权重函数的值越接近1，算法收敛速度快，但是难以达到最优值，h2和j2的取值越小，权重函数的值越接近0，算法经过多次迭代可以收敛到最优值，但是耗费时间多，h2和j2的取值原则上要综合以上两点因素，取大小合适的值；
[0078] (16)为了再次提升去噪效果，对图像u2进行第二轮NLTV去噪。首先建立NLTV模型：其中， 是第二轮NLTV模型
的目标函数，λ2是保真参数，Ω是u2和u3的图像空间，u3是去噪后的图像， 是GUY GILBOA和STANLEY OSHER提出的非局部梯度算子，
其中p2为第二轮NLTV去噪的当前像素点，q2是
以p2为中心的搜索窗口内的一点，u3(p2)、u3(q2)分别是图像上的点p2和q2的像素灰度值。本步骤中使用的权重函数ω2(p2,q2)已在步骤(15)中建立；
[0079] (17)对步骤(16)建立的模型，采用分裂的Bregman迭代实现，迭代过程分成如下三步：
[0080]
[0081]
[0082]
[0083] 其中，k的取值是0，1，2，…,等非负整数，迭代初始值u30＝u2，b2k、w2k分别表示分裂的Bregman迭代的辅助变量和辅助函数，λ2就是步骤(16)建立的非局部TV模型中的保真参k 0数，θ2是控制迭代结果的平滑参数，b2的初始值b2、以及λ2和θ2的赋值在步骤(2)中进行预设。
[0084] 求解式(7)和式(8)，并数值化，式(9)也数值化，由此得到数值化后的迭代格式，如式(10)、式(11)和式(12)所示；第一次迭代运算前，令k＝0；
[0085] (18)顺序地应用公式(10)、(11)、(12)进行迭代运算，
[0086]
[0087]
[0088]
[0089] 其中ω2(p2,q2)是图像u2中像素p2和q2间的相似度权重函数， 表示求和范围是以p2为中心的搜索窗口N3×N3内除p2以外的每一像素点，求和项数共N3×N3-1项；
[0090] (19)计算步骤(18)中的本次迭代输出图像u3k+1的峰值信噪比PSNR，如果本次迭代输出图像u3k+1的峰值信噪比PSNR小于等于上一次迭代输出图像u3k的峰值信噪比PSNR，则将上一次迭代输出图像u3k作为最优值ufinal输出，即令ufinal＝u3k，并转到步骤(20)；如果本次迭代后输出图像的峰值信噪比PSNR大于上一次迭代输出图像的峰值信噪比PSNR，则令k＝k+1，并返回到步骤(18)，继续迭代运算。
[0091] (20)将ufinal作为最后去噪结果图像输出。
[0092] 本发明效果可以通过以下实验进一步证实：
[0093] (一)实验条件
[0094] 使用Matlab软件对Lena图和Peppers图进行测试,输入图像所含高斯白噪声的标准差为15和20。本发明的参数为：
[0095] 标准差为15的Lena噪声图：e＝0.0001、N1×N1＝5×5、N2×N2＝5×5、h1＝9、j1＝31.62、σ1＝5、λ1＝1、θ1＝11.5、b1＝0、w1＝0、N3×N3＝5×5、N4×N4＝5×5、h2＝15、j2＝
31.62、σ2＝5、λ2＝1、θ2＝11.5、b2＝0、w2＝0。
[0096] 标准差为20的Lena噪声图：e＝0.0001、N1×N1＝3×3、N2×N2＝3×3、h1＝14、j1＝31.62、σ1＝5、λ1＝1、θ1＝11.5、b1＝0、w1＝0、N3×N3＝7×7、N4×N4＝5×5、h2＝17.6、j2＝
31.62、σ2＝5、λ2＝1、θ2＝11.5、b2＝0、w2＝0。标准差为15的Peppers噪声图：e＝0.0001、N1×N1＝5×5、N2×N2＝5×5、h1＝9、j1＝31.62、σ1＝5、λ1＝1、θ1＝11.5、b1＝0、w1＝0、N3×N3＝5×5、N4×N4＝5×5、h2＝15、j2＝31.62、σ2＝5、λ2＝1、θ2＝11.5、b2＝0、w2＝0。
[0097] 标准差为20的Peppers噪声图：e＝0.0001、N1×N1＝3×3、N2×N2＝3×3、h1＝14、j1＝31.62、σ1＝5、λ1＝1、θ1＝11.5、b1＝0、w1＝0、N3×N3＝7×7、N4×N4＝5×5、h2＝17.6、j2＝31.62、σ2＝5、λ2＝1、θ2＝11.5、b2＝0、w2＝0。
[0098] (二)实验内容
[0099] 按照上面所述的实验步骤进行Matlab实验仿真,并将本发明算法与NLTV算法进行比较。本发明去噪算法迭代结束条件是一旦本次迭代后输出图像的峰值信噪比PSNR小于上一次迭代的峰值信噪比PSNR，则迭代停止，并将上一次迭代的结果作为最优值输出；如果未满足停止条件，则继续迭代。这样就得到了去噪效果最好的图像。用来对比的NLTV去噪算法迭代结束条件是一旦本次迭代后输出图像的峰值信噪比PSNR小于上一次迭代的峰值信噪比PSNR，则迭代停止，并将上一次迭代的结果作为最优值输出；如果未满足停止条件，则继续迭代。这样就得到了去噪效果最好的图像。
[0100] (三)实验结果
[0101] 实验结果见表一和附图3。实验结果表明本发明算法去噪的峰值信噪比(PSNR)比NLTV模型更高，去噪效果更好。
[0102] 以上所述,仅是本发明的较佳实施例而已,并不对本发明做形式上的限制,凡是依据本发明对以上实例所做的简单修改，等同变化与修饰,均仍属本发明技术方案的范围内。
[0103] 表一 用本发明方法和NLTV方法去噪前后图像的峰值信噪比
[0104]
