基于ICA和NLTV的图像去噪方法与流程

本发明属于图像处理技术领域，涉及了去除加性噪声的图像去噪领域。

背景技术：

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

独立分量分析(Independent Component Analysis，ICA)是20世纪90年代发展起来的一种新的信号处理技术。基本的ICA是指从多个源信号的线性混合信号中分离出各源信号的技术。除了需要已知源信号是统计独立外，无其他先验知识，ICA是伴随着盲信源问题而发展起来的，故又称盲分离。在时间信号处理方面由于ICA可以在n个原信号和线性混合系数都未知的情况下可以把这n个原信号从他们线性混合成的n个混合信号中分离出来，并且它的假设条件仅是这n个原信号之间是相互独立的，所以ICA在时间信号分离方面具有良好的应用优势。在图像信号处理方面，加性噪声图像可以看成是噪声与纯净图像的线性相加，而实际中噪声和原始纯净图像一般是相互独立的，符合ICA应用的假设条件，所以ICA可以应用在含有加性噪声情况下的图像去噪处理。

但ICA方法在进行图像去噪时，需要至少2幅带噪声的图像或者一幅纯净无噪图像和一幅噪声图像，并且图像之间相互独立。在只能提供一幅含噪图像的情况，ICA方法不具备图像去噪应用条件。本发明将ICA和NLTV相结合进行图像去噪。在对含加性噪声的输入含噪图像进行NLTV初步去噪后，得到初步去噪图像，这个初步去噪图像在一定程度上可以认为是另一张带噪声图像。再利用这二幅图像，应用ICA方法进行图像去噪。

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模型具有在去除图像噪声的同时能较好地保留图像的纹理细节的特点。本文中所提到的NLTV模型是指GUY GILBOA和STANLEY OSHER提出的非局部TV模型。

技术实现要素：

本发明要解决ICA技术只能应用于只有单一输入含噪图像情况下的图像去噪问题的缺点，提出一种基于ICA和NLTV的图像去噪方法，提高了ICA在图像去噪领域的应用价值。

本发明方法利用NLTV方法对输入的含噪图像u₀进行初步去噪，得到初步去噪图像u₁，在获得的u₁、u₀基础上应用ICA方法对输入图像u₀去噪，得到ICA方法去噪后的图像u₂。为了获得更好的去噪效果，对u₂再一次进行NLTV去噪，输出的去噪结果图像作为最后去噪结果图像。本发明方法和NLTV相比具有更好的图像去噪效果。去噪后的图像有更高的峰值信噪比。

本发明所述的基于ICA和NLTV的图像去噪方法，包括如下步骤：

(1)首先输入N×N大小的含噪图像u₀。

(2)设置本发明去噪方法中的相关参数，包括ICA迭代收敛阈值e，第一轮NLTV的搜索窗口大小N₁×N₁、邻域窗口大小N₂×N₂、保真参数λ₁、像素相似度权重函数ω₁的参数h₁、像素相似度权重函数ω₁中高斯核的标准差σ₁，分裂的Bregman迭代辅助变量b₁^k的初始值b₁⁰、平滑参数θ₁。第二轮NLTV的搜索窗口大小N₃×N₃、邻域窗口大小N₄×N₄、保真参数λ₂、像素相似度权重函数ω₂的权重函数参数h₂、像素相似度权重函数ω₂中高斯核的标准差σ₂，分裂的Bregman迭代辅助变量b₂^k的初始值b₂⁰、平滑参数θ₂。

(3)对噪声图像u₀用NLTV方法进行第一轮去噪。在本步骤中，对NLTV模型采用分裂的Bregman算法进行数值迭代求解，迭代一次后获得初步去噪图像u₁。用于基于ICA方法的图像去噪[步骤(4)-(10)]。

(4)将u₀与u₁看成是纯噪图像与无噪图像线性混合成的两个信号，基于这两个图像信号构建混合矩阵S，对S进行中心化后得到矩阵S_e，对S_e进行白化处理，得到处理结果Z。

(5)建立初始解混矩阵B。

(6)建立初始随机列向量L，L各向量元素取值范围是[0,1]。

(7)对L进行迭代。目的是通过L的迭代运算结果来逼近解混矩阵B的一个列向量的真值。

(8)若满足L迭代的停止条件，则转向步骤(9)，否则返回步骤(7)。

(9)用迭代结果L替换B的一个列向量，如果B的列向量全被替换过一次，则转向步骤(10)，否则返回步骤(6)。

(10)进行W＝B^TS运算。W的第一行行向量经转换即得到去噪后的图像u₂，转换规则是：W的第一个行向量的第(α-1)×N+β列元素转换为u₂的第α列第β行元素，其中α和β取值均为1,2,3,...,N。

(11)为了对u₂进一步去噪，建立第二轮NLTV模型。

(12)对步骤(11)建立的NLTV模型，用分裂的Bregman算法求解，得到三步式数值迭代计算公式。并设迭代计数变量初值k＝0。

(13)用步骤(12)中得到的分裂的Bregman算法三步式数值迭代计算公式进行顺序迭代运算，得到本次迭代的输出图像u₃^k+1。

(14)计算步骤(13)中的迭代输出图像u₃^k+1的峰值信噪比PSNR，如果本次迭代后输出图像u₃^k+1的峰值信噪比PSNR小于等于上一次迭代输出图像u₃^k的峰值信噪比PSNR，则迭代结束，将上一次迭代输出图像u₃^k作为最优值u_final输出，即令u_final＝u₃^k，并转到步骤(15)；如果本次迭代后输出图像u₃^k+1的峰值信噪比PSNR大于上一次迭代输出图像u₃^k的峰值信噪比，则更新k值：k＝k+1,并返回到步骤(13)，继续迭代运算。

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

本发明的优点是：将具有只要满足n个原信号之间相互独立，不需要其他的先验条件之优势的ICA应用于单幅图像去噪；通过NLTV去噪方法对输入的含噪声图像u₀进行初步去噪，利用初步去噪后的图像u₁与原噪声图像u₀进行ICA迭代，去除图像u₀的噪声，得到去噪后图像u₂，然后再对u₂再一次进行NLTV去噪，输出的去噪结果图像作为最后去噪结果图像，以获得更好的去噪效果；本发明采用了分裂Bregman迭代进行NLTV模型的数值实现，它将迭代过程分成三步，降低了迭代的复杂性，加快去噪运算的速度。

附图说明

附图1是本发明方法的流程图。

附图2是本发明仿真实验所用的原始图像，图2a是Lena图像，图2b是Peppers图像。

附图3是本发明对原始Peppers图像加均值为零标准差为15的高斯噪声后的图像进行去噪的仿真结果，其中，图3a是待去噪的含噪图像(所含噪声的标准差为15)，图3b是NLTV去噪结果，图3c是本发明方法的去噪结果。

具体实施方式

下面结合附图进一步说明本发明。

本发明所述的基于ICA和NLTV的图像去噪方法，步骤如下：

(1)首先输入N×N大小的含噪图像u₀。

(2)设置本发明方法中的相关参数，包括ICA迭代收敛阈值e，第一轮NLTV的搜索窗口大小N₁×N₁、邻域窗口大小N₂×N₂、保真参数λ₁、像素相似度权重函数ω₁的参数h₁、像素相似度权重函数ω₁中高斯核的标准差σ₁，分裂的Bregman迭代辅助变量b₁^k的初始值b₁⁰、平滑参数θ₁。第二轮NLTV的搜索窗口大小N₃×N₃、邻域窗口大小N₄×N₄、保真参数λ₂、像素相似度权重函数ω₂的参数h₂、像素相似度权重函数ω₂中高斯核的标准差σ₂，分裂的Bregman迭代辅助变量b₂^k的初始值b₂⁰、平滑参数θ₂。

(3)本发明中为了应用ICA方法进行图像去噪，除了输入的含噪图像u₀，还需另一幅含噪图像u₁，本发明通过用NLTV方法对含噪图像u₀进行初步去噪来获得u₁(称为第一轮NLTV去噪)。首先建立NLTV模型(称为第一轮NLTV模型)：其中J(u₁)是第一轮NLTV模型的目标函数，Ω是u₀、u₁的图像空间，λ₁是保真参数，u₀是输入的含噪声的待去噪图像，u₁是去噪后的图像；是GUY GILBOA和STANLEY OSHER提出的非局部梯度算子，其中变量p₁表示当前像素点，变量q₁表示以p₁为中心的搜索窗口内的一点，u₁(p₁)、u₁(q₁)分别是图像u₁上的点p₁和q₁的像素灰度值；是图像u₀中的两个像素点p₁和q₁的相似度权重函数，表示图像u₀中以p₁为中心的大小为N₂×N₂的图像块和以q₁为中心的大小为N₂×N₂的图像块之间的高斯加权距离，表示求和范围是以p₁或q₁为中心的N₂×N₂邻域内的每一像素点(不包括p₁或q₁自身)，共N₂×N₂-1项,是标准差为σ₁的高斯核，h₁是权重函数ω₁(p₁,q₁)的常数参数，在步骤(2)中对h₁、σ₁进行赋值。对该NLTV模型，采用分裂的Bregman算法进行迭代逼近求解。引入辅助函数w₁^k和辅助变量b₁^k，构造如下三步迭代格式：

其中，k的取值是0，1，2，…,等非负整数，迭代初始值u₁⁰＝u₀，b^k、w^k分别表示分裂的Bregman迭代的辅助变量和辅助函数，θ₁是控制迭代结果的平滑参数，b₁^k的初始值b₁⁰、以及λ₁和θ₁的赋值在步骤(2)中进行预设。

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

顺序地应用公式(3‐4)、(3‐5)、(3‐6)进行迭代，

其中表示求和范围是以p₁为中心的搜索窗口N₁×N₁内除p₁以外的每一像素点，求和项数共N₁×N₁-1项，ω₁(p₁,q₁)是u₀的两个像素点p₁和q₁的相似度权重函数。

在本步骤中，设初值k＝0，顺序地按式(3-4)、(3-5)、(3-6)进行迭代运算一次，得到初步去噪后的图像u₁＝u₁¹；

(4)使用u₀与u₁进行中心化和白化处理。具体方法是：首先将大小N×N的图像u₀转换成1×N²的行向量X，转换规则是u₀的第α(α＝1,2,3,...,N)列第β(β＝1,2,3,...,N)行元素转为X的第(α-1)×N+β列元素；用同样方法将大小N×N的图像u₁转换成1×N²的行向量Y。用X和Y构建混合矩阵并进行中心化和白化。中心化是将原始数据减去平均数。白化也称为球化，它的本质是去相关。如果均值为零的随机向量O＝[o₁,o₂,...,o_n]^T满足E{OO^T}＝I，其中I是单位矩阵，那么随机向量O＝[o₁,o₂,...,o_n]^T是白化向量。

对S进行中心化后得到矩阵S_e：

其中x_i和y_i分别是X和Y的第i个元素，i＝1,2,…,N×N，是X的所有元素的均值，是Y的所有元素的均值，

对S_e进行白化处理，即白化矩阵W₀与S_e相乘，得到Z：

Z＝W₀S_e

其中白化矩阵W₀＝Λ^-1/2U^T，Λ是S_e^T的协方差矩阵的特征值矩阵，U是S_e^T的协方差矩阵的特征向量矩阵。Z将参与后续的步骤。这一处理可以降低后续步骤的计算复杂度。

(5)设置解混矩阵B的初始值为2×2的零矩阵。

(6)建立一个大小为2×1，的随机列向量L，L各元素取值范围是[0～1]。

(7)对L进行迭代。迭代公式如下：

L＝E{Zg(L^TZ)}-E{g₁(L^TZ)}L，

L＝L-BB^TL,

L＝L/||L||,

其中E{·}是均值运算，g(·)为任意二次函数,在本发明中令g(·)为g(t)＝t²。g₁(·)是g(·)的一阶导数。

(8)如果L满足||L^TL|-1|<e，则转到步骤(9)；如果L不满足||L^TL|-1|<e，则返回到步骤(7)。其中e是ICA方法中的收敛阈值，为常数，其值在步骤(2)预先设置。

(9)用L来替换更新B中的一个列向量。如果B的所有列向量都被替换更新一次，则转到步骤(10)。如果B的所有列向量没有被替换更新完，则返回到步骤(6)。

(10)利用得到的解混矩阵B分离出噪声，获得去噪后的图像。步骤如下：先计算解混结果W＝B^TS，再将W中的两个1×N²大小行向量都转换成N×N大小的矩阵，则可以得到分离后的去噪图像u₂和噪声图像。转换规则是：W的第一个行向量的第(α-1)×N+β列元素转为u₂的第α列第β行元素，其中α和β取值均为1,2,3,...,N。

(11)为了再次提升去噪效果，将u₂进行第二轮NLTV去噪。建立NLTV模型其中J(u₃)是第二轮NLTV模型的目标函数，Ω表示u₂、u₃的图像空间,λ₂是保真参数，u₃是去噪后的图像；是GUY GILBOA和STANLEY OSHER提出的非局部梯度算子，其中p₂表法当前像素点，q₂给示以p₂为中心的搜索窗口内的一点，u₃(p₂)、u₃(q₂)分别是图像u₃上的点p₂和q₂的像素灰度值；

是图像u₂的两个像素点p₂和q₂的相似度权重函数，表示图像u₂中以p₂为中心的大小为N₄×N₄的图像块和以q₂为中心的大小为N₄×N₄的图像块之间的高斯加权距离，表示求和范围是以p₂或q₂为中心的N₄×N₄邻域内的每一像素点(不包括p₂或q₂自身)，共N₄×N₄-1项，是标准差为σ₂的高斯核，h₂是权重函数ω₂(p₂,q₂)的常数参数，在步骤(2)中对h₂、σ₂进行赋值。

(12)对步骤(11)建立的NLTV模型，采用分裂的Bregman迭代实现，迭代过程分成三步，如式(12-1)、(12-2)和式(12-3)所示，降低了迭代的复杂性，加快了迭代速度。

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

求解式(12-1)和式(12-2)，并数值化，式(12-3)也数值化，由此得到数值化后的迭代格式，如式(13-1)、(13-2)和式(13-3)所示。第一次迭代前，设k＝0。

(13)顺序地按式(13-1)、(13-2)和式(13-3)进行迭代计算。

其中ω₂(p₂,q₂)是u₂的像素点p₂和q₂间的相似度权重函数，表示求和范围是以p₂为中心的搜索窗口N₃×N₃内除p₂以外的每一像素点，求和项数共N₃×N₃-1项；

(14)计算步骤(13)的输出图像u₃^k+1的峰值信噪比PSNR，如果本次迭代后输出图像u₃^k+1的峰值信噪比PSNR小于等于上一次迭代输出图像u₃^k的峰值信噪比PSNR，则满足迭代停止条件，将上一次迭代输出图像u₃^k作为最优值u_final输出，即令u_final＝u₃^k，并转到步骤(15)；如果本次迭代后输出图像u₃^k+1的峰值信噪比PSNR大于上一次迭代输出图像u₃^k的峰值信噪比PSNR，则不满足迭代停止条件，更新k值：k＝k+1，并返回到步骤(13)，继续迭代运算。

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

本发明的效果可以通过以下实验进一步证实。

1、实验条件：

使用Matlab软件对256×256尺寸的Lena、Peppers和House图像进行仿真实验。在图2所示的原始干净的图像上分别叠加标准差为15和20的高斯噪声(均值为零)后作为待去噪的输入图像。本实验的参数为：

对于标准差为15的Lena噪声图：e＝0.0001、N₁×N₁＝5×5、N₂×N₂＝5×5、h₁＝9、σ₁＝6、λ₁＝1、θ₁＝11.5、b₁＝0、w₁＝0、N₃×N₃＝5×5、N₄×N₄＝5×5、h₂＝13、σ₂＝6、λ₂＝1、θ₂＝11.5、b₂＝0、w₂＝0。

对于标准差为20的Lena噪声图：e＝0.0001、N₁×N₁＝3×3、N₂×N₂＝3×3、h₁＝9、σ₁＝6、λ₁＝1、θ₁＝11.5、b₁＝0、w₁＝0、N₃×N₃＝5×5、N₄×N₄＝5×5、h₂＝13、σ₂＝6、λ₂＝1、θ₂＝11.5、b₂＝0、w₂＝0。

对于标准差为15的Peppers噪声图：e＝0.0001、N₁×N₁＝5×5、N₂×N₂＝5×5、h₁＝9、σ₁＝6、λ₁＝1、θ₁＝11.5、b₁＝0、w₁＝0、N₃×N₃＝5×5、N₄×N₄＝5×5、h₂＝13、σ₂＝8、λ₂＝1、θ₂＝11.5、b₂＝0、w₂＝0。

对于标准差为20的Peppers噪声图：e＝0.0001、N₁×N₁＝3×3、N₂×N₂＝3×3、h₁＝14、σ₁＝6、λ₁＝1、θ₁＝11.5、b₁＝0、w₁＝0、N₃×N₃＝5×5、N₄×N₄＝5×5、h₂＝14、σ₂＝6、λ₂＝1、θ₂＝11.5、b₂＝0、w₂＝0。

2、实验内容

按照上面所述的实验步骤进行Matlab实验仿真,并将本发明方法与NLTV方法(采用分裂的Bregman算法迭代逼近来求解)进行比较。

3、实验结果

实验结果见表一以及附图3。实验结果表明本发明方法比NLTV方法有更高的峰值信噪比(PSNR)，去噪效果更好。

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

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