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

技术特征：

1.基于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，构造如下三步迭代格式：

${w_1}^{k+1}=\arg \underset{w_1}{min}\left\{\∫\left|w_1\right|{\mathrm{dp}}_1+\frac{{\mathrm{\θ}}_1}{2}\∫{(w_1-{\▿}_{NL}{u_1}^k-{b_1}^k)}^2{\mathrm{dp}}_1\right\}---\left(3-1\right)$

${u_1}^{k+1}=\arg \underset{u_1}{min}\{{\mathrm{\λ}}_1\∫{(u_1-u_0)}^2{\mathrm{dp}}_1+\frac{{\mathrm{\θ}}_1}{2}\∫{({w_1}^{k+1}-{\▿}_{NL}{u}_1-{b_1}^k)}^2{\mathrm{dp}}_1\}---\left(3-2\right)$

${b_1}^{k+1}={b_1}^k+{\▿}_{NL}{u_1}^{k+1}-{w_1}^{k+1}---\left(3-3\right)$

其中，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)进行迭代，

${w_1}^{k+1}\left(p_1,q_1\right)=\max \left(\sqrt{\underset{q_1}{\Σ}{\left\{\sqrt{{\mathrm{\ω}}_1\left(p_1,q_1\right)}\[{u_1}^k\left(q_1\right)-{u_1}^k\left(p_1\right)\]+{b_1}^k\left(p_1,q_1\right)\right\}}^2}-\frac{1}{{\mathrm{\θ}}_1},0\right)\frac{\sqrt{{\mathrm{\ω}}_1\left(p_1,q_1\right)}\[{u_1}^k\left(q_1\right)-{u_1}^k\left(p_1\right)\]+{b_1}^k\left(p_1,q_1\right)}{\sqrt{\underset{q_1}{\Σ}{\left\{\sqrt{{\mathrm{\ω}}_1\left(p_1,q_1\right)}\[{u_1}^k\left(q_1\right)-{u_1}^k\left(p_1\right)\]+{b_1}^k\left(p_1,q_1\right)\right\}}^2}}---\left(3-4\right)$

${u_1}^{k+1}\left(p_1\right)=\frac{{\mathrm{\λ}}_1{u_1}^k\left(p_1\right)+\frac{{\mathrm{\θ}}_1}{2}\underset{q_1}{\Σ}\left\{2{u_1}^k\left(q_1\right){\mathrm{\ω}}_1\left(p_1,q_1\right)-\sqrt{{\mathrm{\ω}}_1\left(p_1,q_1\right)}\[{w_1}^{k+1}\left(p_1,q_1\right)-{w_1}^{k+1}\left(q_1,p_1\right)+{b_1}^k\left(q_1,p_1\right)-{b_1}^k\left(p_1,q_1\right)\]\right\}}{{\mathrm{\λ}}_1+{\mathrm{\θ}}_1\underset{q_1}{\Σ}{\mathrm{\ω}}_1\left(p_1,q_1\right)}---\left(3-5\right)$

${b_1}^{k+1}\left(p_1,q_1\right)={b_1}^k\left(p_1,q_1\right)+\sqrt{{\mathrm{\ω}}_1\left(p_1,q_1\right)}\[{u_1}^{k+1}\left(q_1\right)-{u_1}^{k+1}\left(p_1\right)\]-{w_1}^{k+1}\left(p_1,q_1\right)---\left(3-6\right)$

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

在本步骤中，设初值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：

$S_e=\left[\begin{array}{c}{x}_1-\stackrel{\‾}{X},x_2-\stackrel{\‾}{X},\mathrm{...},x_{N\×N}-\stackrel{\‾}{X}\\ y_1-\stackrel{\‾}{Y},y_2-\stackrel{\‾}{Y},\mathrm{...},y_{N\×N}-\stackrel{\‾}{Y}\end{array}\right],$

其中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)所示，降低了迭代的复杂性，加快了迭代速度；

${w_2}^{k+1}=\arg \underset{w_2}{\min }\left\{\&Integral;\left|w_2\right|{\mathrm{dp}}_2+\frac{{\mathrm{\&theta;}}_2}{2}\&Integral;{\left(w_2-{\&dtri;}_{NL}{u_3}^k-{b_2}^k\right)}^2{\mathrm{dp}}_2\right\}---\left(12-1\right)$

${u_3}^{k+1}=\arg \underset{u_3}{\min }\left\{{\mathrm{\&lambda;}}_2\&Integral;{\left(u_3-u_2\right)}^2{\mathrm{dp}}_2+\frac{{\mathrm{\&theta;}}_2}{2}\mathrm{\&lambda;}{\left({w_2}^{k+1}-{\&dtri;}_{NL}{u}_3-{b_2}^k\right)}^2{\mathrm{dp}}_2\right\}---\left(12-2\right)$

${b_2}^{k+1}={b_2}^k+{\&dtri;}_{NL}{u_3}^{k+1}-{w_2}^{k+1}---\left(12-3\right)$

其中，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)进行迭代计算；

${w_2}^{k+1}\left(p_2,q_2\right)=\max \left(\sqrt{\underset{q_2}{\&Sigma;}{\left\{\sqrt{{\mathrm{\&omega;}}_2\left(p_2,q_2\right)}\&lsqb;{u_3}^k\left(q_2\right)-{u_3}^k\left(p_2\right)\&rsqb;+{b_2}^k\left(p_2,q_2\right)\right\}}^2}-\frac{1}{{\mathrm{\&theta;}}_2},0\right)\frac{\sqrt{{\mathrm{\&omega;}}_2\left(p_2,q_2\right)}\&lsqb;{u_3}^k\left(q_2\right)-{u_3}^k\left(p_2\right)\&rsqb;+{b_2}^k\left(p_2,q_2\right)}{\sqrt{\underset{q_2}{\&Sigma;}{\left\{\sqrt{{\mathrm{\&omega;}}_2\left(p_2,q_2\right)}\&lsqb;{u_3}^k\left(q_2\right)-{u_3}^k\left(p_2\right)\&rsqb;+{b_2}^k\left(p_2,q_2\right)\right\}}^2}}---\left(13-1\right)$

${u_3}^{k+1}\left(p_2\right)=\frac{{\mathrm{\&lambda;}}_2{u_3}^k\left(p_2\right)+\frac{{\mathrm{\&theta;}}_2}{2}\underset{q_2}{\&Sigma;}\left\{2{u_3}^k\left(q_2\right){\mathrm{\&omega;}}_2\left(p_2,q_2\right)-\sqrt{{\mathrm{\&omega;}}_2\left(p_2,q_2\right)}{w_2}^{k+1}\left(p_2,q_2\right)-{w_2}^{k+1}\left(q_2,p_2\right)+{b_2}^k\left(q_2,p_2\right)-{b_2}^k\left(p_2,q_2\right)\right\}}{{\mathrm{\&lambda;}}_2+{\mathrm{\&theta;}}_2\underset{q_2}{\&Sigma;}{\mathrm{\&omega;}}_2\left(p_2,q_2\right)}---\left(13-2\right)$

${b_2}^{k+1}\left(p_2,q_2\right)={b_2}^k\left(p_2,q_2\right)+\sqrt{{\mathrm{\&omega;}}_2\left(p_2,q_2\right)}\&lsqb;{u_3}^{k+1}\left(q_2\right)-{u_3}^{k+1}\left(p_2\right)\&rsqb;-{w_2}^{k+1}\left(p_2,q_2\right)---\left(13-3\right)$

其中ω₂(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作为最终去噪结果图像输出。