CN110298803B - 一种基于压缩感知的图像条带噪声去除方法 - Google Patents

Abstract

本发明涉及的是一种基于压缩感知的图像条带噪声去除方法，它包括：步骤一、建立包含条带噪声图像的观测矩阵；步骤二、建立包含条带噪声图像Curvelet域的稀疏表示方法；步骤三、建立条带噪声去除正则化问题最优化的求解算法，去除条带噪声，得到去除条带噪声的图像。本发明在条带噪声去除过程，提出的凸投影和滤波操作交替进行的指数阈值收缩迭代算法，在迭代初期快速阈值收缩，加速噪声的去除；在迭代后期减缓阈值收缩，加强图像细节信息保留，能充分利用Curvelet域稀疏表示性能与压缩感知的优势，提高去噪效果。

Description

一种基于压缩感知的图像条带噪声去除方法

技术领域：

本发明涉及的是图像处理中条带噪声去除技术，具体涉及的是一种基于压缩感知的图像条带噪声去除方法。

背景技术：

图像中的条带噪声具有一定的随机性、方向性、突变性，传统的图像条带噪声去除方法(统计学方法、辐射均衡化方法和滤波方法)存在一定的局限性。统计学方法适用于图像灰度分布均匀的情况，反之效果不理想。辐射均衡化方法要求图像范围内的目标具有同一性或近似同一性，当图像中的目标物不具备同一性时，效果会受到较大的影响。滤波法主要有傅立叶变换滤波法和小波变换滤波法等。该类方法的问题在于很难找到一个合适的频率将噪声和信号完全分离，导致在有效信号和噪声之间的频率重叠区域，不能彻底去除条带噪声，甚至会出现局部灰度退化现象。

发明内容：

本发明的目的是提供一种基于压缩感知的图像条带噪声去除方法，这种基于压缩感知的图像条带噪声去除方法用于解决现有图像条带噪声去除方法不能彻底去除条带噪声，甚至会出现局部灰度退化的问题。

本发明解决其技术问题所采用的技术方案是：这种基于压缩感知的图像条带噪声去除方法:

步骤一、建立包含条带噪声图像的观测矩阵；

设不含条带噪声的图像为x，大小为N_r×N_s，包含条带噪声的图像为y(N_r×N_s)，观测矩阵构造方法描述为，

对于水平方向条带噪声，通过y＝R*x得到包含水平方向条带噪声的观测图像；

对于竖直方向的条带噪声，通过y^T＝R*x^T得到包含竖直方向条带噪声的观测图像，T为矩阵转置运算；

(1)对于给定图像的条带噪声所在图像中行序号的集合Ω，

(2)设置观测矩阵R为一个大小为N_r×N_r的单位矩阵；

(3)处理任意一个行序号ω，

循环：

设置R(ω,ω)＝0

(4)得到包含条带噪声图像的观测矩阵R；

步骤二、建立包含条带噪声图像Curvelet域的稀疏表示方法，其包括：

(1)Curvelet域双变量阈值函数：

其中γ_i表示当前系数，γ_i-1表示当前系数的父系数，(g)₊运算表示为

σ_ε为条带噪声信号系数均方差，σ为图像有效信号系数均方差，τ为收敛控制因子；

Curvelet域双变量阈值函数中的参数估计：采用最高尺度子带的中值估计作为σ_ε的估计值，以当前系数γ_i为中心的3*3邻域窗口的均方差作为σ的估计值，通过迭代过程中逐渐收缩控制因子τ，逼近原始图像，提高去噪质量；

(2)Curvelet域相邻尺度系数γ_i和γ_i-1对应关系确定方法：

建立变换系数父子对应关系，设当前系数γ_i所在的子带矩阵

尺寸为m_i行，n_i列，l为方向，i为尺度，k为位置；设父系数γ_i-1所在子带矩阵为/>

的尺寸为m_i-1行，n_i-1列；m_i与m_i-1的关系存在m_i＞m_i-1和m_i≤m_i-1两种情况，n_i和n_i-1的关系也存在n_i＞n_i-1和n_i≤n_i-1两种情况；

当m_i＞m_i-1时需要行扩展

矩阵，使得m_i＝m_i-1的处理方法：

(1)令

(2)如果p＞1

对于任意一个c∈{1,2…p}循环：

插值扩展

(3)令q＝mod(m_i,m_i-1)

(4)如果q＞1

对于任意一个c∈{1,2…q}循环：

对称扩展

(5)得到扩展之后的矩阵

截取1:m_i行/>

得到新的

同样的方式处理当m_i≤m_i-1；当n_i＞n_i-1时；当n_i≤n_i-1时三种情况；

同理，在处理完n_i和n_i-1的情况后，得到的矩阵

和/>

尺寸相同，系数之间存在1:1的对应关系；

步骤三、建立条带噪声去除正则化问题最优化的求解算法；

(1)给定迭代停止参数δ，收敛控制因子τ，最大迭代次数K，包含条带噪声图像观测数据y；

(2)设置迭代计数器k＝0，Curvelet变换算子为W，Curvelet变换级数为z，迭代初值为x⁽⁰⁾＝y，Curvelet变换系数为α,变换尺度为i的Curvelet变换系数为α_i；

(3)条带噪声去除处理过程：

循环：

对x^k进行3*3的维纳窗口滤波

凸投影操作x^(k+1)＝x^k+R^T(y-R^T(x^k))

Curvelet正变换：α^(k+1)＝W(x^(k+1))

对于任意一个变换尺度i∈{1,2…z}循环：

双变量阈值函数滤波：

Curvelet反变换：x^k+1＝W^T(α^(k+1))

凸投影操作x^(k+1)＝x^(k+1)+R^T(y-R^T(x^(k+1)))

指数模型收缩τ

k＝k+1

直到：||x^(k+1)-x^(k)||₂＜＝δ或者k≥K循环结束；

x＝x^(k+1)；

(4)得到去除条带噪声的图像x。

本发明具有以下有益效果：

(1)本发明通过信号处理、稀疏表示及最优化求解等技术，将受条带噪声干扰的原始图像内容恢复出来，以达到条带噪声去除的目的。该方法是基于数据驱动的，具有广泛的适应性。本发明提出的观测矩阵建立方法，生成的观测矩阵与图像数据维数的数量级相同，节省了内存空间；提出的Curvelet域相邻尺度间系数父子关系的对应是通过同一方向，相邻尺度间子带矩阵的插值扩展、对称扩展方法实现，提高了准确性；提出的Curvelet域双变量阈值收缩公式在去噪过程加强了图像细节的保留，提出的系统设计与实现的方法采用面向对象的技术，具有较高的实用性。

(2)本发明在条带噪声去除过程，提出的凸投影和滤波操作交替进行的指数阈值收缩迭代算法，在迭代初期快速阈值收缩，加速噪声的去除；在迭代后期减缓阈值收缩，加强图像细节信息保留，能充分利用Curvelet域稀疏表示性能与压缩感知的优势，提高去噪效果。

(3)本发明通过该相邻尺度系数γ_i和γ_i-1对应关系的确定方法，使得Curvelet域父子系数间相关性更加明显，提高稀疏表示效果。

附图说明

图1是本发明根据包含条带噪声图像(左下)，构造的观测矩阵结构(右下)示意图，该观测矩阵反映出图像中的第1，3……等行包含水平方向条带噪声的情况。

图2是本发明根据同方向、较高尺度(子系数)子带矩阵的尺寸，插值扩展较低尺度(父系数)子带矩阵的示意图，假设

的情况。

图3是包含竖直随机条带噪声的barbara图像；

图4是利用本发明方法迭代100次去噪处理后的图3；

图5是包含竖直随机条带噪声的lenna图像；

图6是利用本发明方法迭代100次去噪处理后的图5；

图7是本发明去除包含竖直随机条带噪声barbara图像和lenna图像两个实验PSNR随迭代次数的变化曲线；

图8是包含水平随机条带噪声的barbara图像；

图9是利用本发明方法迭代100次去噪处理后的图8；

图10是包含水平随机条带噪声的lenna图像；

图11是利用本发明方法迭代100次去噪处理后的图10；

图12是本发明去除包含水平随机条带噪声barbara图像和lenna图像两个实验PSNR随迭代次数的变化曲线。

具体实施方式

下面结合附图对本发明做进一步的说明：

结合图1、图2所示，这种基于压缩感知的图像条带噪声去除方法和系统，将图像中的条带噪声认为是以条带形式、随机缺失的图像信息，这种随机缺失的观测结果等价于通过服从某种随机分布的矩阵跟原始图像做乘积运算得到的观测结果，因此条带噪声的去除转化为在压缩感知的思想下恢复出原始图像的问题。本发明在压缩感知理论框架下，利用Curvelet变换具有良好的图像稀疏表示性能，及L₁范数求解优化方法，首先，建立图像条带噪声观测矩阵；然后，建立Curvelet域条带噪声图像的稀疏表示方式，主要包括针对包含条带噪声图像Curvelet域多尺度间系数联系对应，双变量阈值收缩公式建立两部分内容；然后，求解算法采用凸投影和滤波操作交替进行的指数阈值收缩迭代算法，提高了条带噪声去除效果。

具体如下：

一、建立包含条带噪声图像的观测矩阵

设不含条带噪声的图像为x，(大小为N_r×N_s)，包含条带噪声的图像为y(N_r×N_s)，观测矩阵构造方法描述为(以水平方向条带噪声为例)：

(1)对于给定图像的条带噪声所在图像中行序号的集合Ω，

(2)设置观测矩阵R为一个大小为N_r×N_r的单位矩阵。

(3)处理任意一个行序号ω，

循环：设置R(ω,ω)＝0

(4)得到包含条带噪声图像的观测矩阵R。

本发明提出通过y＝R*x即可得到包含水平方向条带噪声的观测图像，同理，对于竖直方向的条带噪声可以转化为由y^T＝R*x^T得到。

二、建立包含条带噪声图像Curvelet域的稀疏表示方法，主要包括针对含条带噪声图像提出新的Curvelet域双变量阈值收缩公式、新的多尺度间系数对应关系两部分内容。

1.Curvelet域双变量阈值收缩公式描述如下：

上式为Curvelet域基于贝叶斯估计的包含条带噪声图像双变量阈值收缩公式，其中γ_i表示当前系数，γ_i-1表示当前系数的父系数，(g)₊运算表示为

σ_ε为条带噪声信号系数均方差，σ为图像有效信号系数均方差，τ为收敛控制因子。采用最高尺度子带的中值估计作为σ_ε的估计值，以当前系数γ_i为中心的3*3邻域窗口的均方差作为σ的估计值，本发明通过迭代过程中逐渐收缩控制因子τ，逼近原始图像，从而提高去噪质量。

2.Curvelet域相邻尺度系数γ_i和γ_i-1对应关系确定方法：

设当前系数γ_i所在的子带矩阵

尺寸为m_i行，n_i列，l为方向，i为尺度，k为位置；设父系数γ_i-1所在子带矩阵为/>

的尺寸为m_i-1行，n_i-1列。由于Curvelet域

的尺寸与/>

的尺寸没有统一的规律，二者的大小跟方向与尺度关系很大，因此m_i与m_i-1的关系存在m_i＞m_i-1和m_i≤m_i-1两种情况，类似的n_i和n_i-1的关系也存在两种情况。下面给出当m_i＞m_i-1时需要行扩展/>

矩阵，使得m_i＝m_i-1的处理方法描述，其他3种情况处理方法类似：

(1)令

(2)如果p＞1

对于任意一个c∈{1,2…p}循环：

插值扩展

(3)令q＝mod(m_i,m_i-1)

(4)如果q＞1

对于任意一个c∈{1,2…q}循环：

对称扩展

(5)得到扩展之后的矩阵

截取1:m_i行/>

得到新的

同理，在处理完n_i和n_i-1的情况后，得到的矩阵

和/>

尺寸相同，系数之间存在1:1的对应关系，本发明通过该相邻尺度系数γ_i和γ_i-1对应关系的确定方法，使得Curvelet域父子系数间相关性更加明显，提高稀疏表示效果。

三、设计条带噪声去除正则化问题最优化的求解算法描述为：

(1)给定迭代停止参数δ，收敛控制因子τ，最大迭代次数K，包含条带噪声图像观测数据y。

(2)设置迭代计数器k＝0，Curvelet变换算子为W，Curvelet变换级数为z，迭代初值为x⁽⁰⁾＝y。

(3)条带噪声去除处理过程：

循环：

对x^k进行3*3的维纳窗口滤波

凸投影操作x^(k+1)＝x^k+R^T(y-R^T(x^k))

Curvelet正变换：α^(k+1)＝W(x^(k+1))

对于任意一个变换尺度i∈{1,2…z}循环：

双变量阈值函数滤波：

Curvelet反变换：x^k+1＝W^T(α^(k+1))

凸投影操作x^(k+1)＝x^(k+1)+R^T(y-R^T(x^(k+1)))

指数模型收缩τ

k＝k+1

直到：||x^(k+1)-x^(k)||₂＜＝δ或者k≥K循环结束。

x＝x^(k+1)。

(4)得到去除条带噪声的图像x。

本发明在条带噪声去除过程，提出的凸投影和滤波操作交替进行的指数阈值收缩迭代算法，在迭代初期快速阈值收缩，加速噪声的去除；在迭代后期减缓阈值收缩，加强图像细节信息保留，能充分利用Curvelet域稀疏表示性能与压缩感知的优势，提高去噪效果。

实验证明实例:

以标准测试图像中的barbara和lenna图像为代表，由于lenna图像中平滑区域较多，barbara图像纹理区域较多，因此该实验用例具有一定的代表性，分别有40％的列包含竖直随机条带噪声如图3、图5所示，利用本发明方法迭代100次去噪结果如图4、图6所示，PSNR随迭代次数的变化曲线如图7所示，PSNR分别提高17dB和24dB左右。

为检验本发明对于水平条带噪声的去除效果，分别在barbara图像和lenna图像中50％的行加入水平随机条带噪声如图8、图10所示，图9、图11为在本发明方法迭代100次后的效果，PSNR分别提高19dB和24dB左右，PSNR随迭代次数的变化曲线如图12所示。

综上所述，与其他的图像条带噪声去除方法相比，本发明具有较好的适应性与噪声去除效果，在纹理区域和平滑区域的条带噪声去除结果都较为理想，并且本发明方法与系统运行具有良好的收敛性、稳定性。

本发明旨在利用压缩感知理论，提出一套方法组合和相应系统实现方法，提高条带噪声去除的适应性，改善条带噪声去除的效果、从技术上保证条带噪声去除的准确性与科学性。因此，本发明对于改进当前图像条带噪声去除方法，提高图像条带噪声去除效果具有重要意义。

以上所述为本发明的较佳实施例而已，但本发明不应该局限于该实施例和附图所公开的内容。所以凡是不脱离本发明所公开的技术方案下完成的等效或修改，都落入本发明保护的范围。

Claims (1)

1.一种基于压缩感知的图像条带噪声去除方法，其特征在于包括如下步骤：

步骤一、建立包含条带噪声图像的观测矩阵；

设不含条带噪声的图像为x，大小为N_r×N_s，包含条带噪声的图像为y(N_r×N_s)，观测矩阵构造方法描述为，

对于水平方向条带噪声，通过y＝R*x得到包含水平方向条带噪声的观测图像；

对于竖直方向的条带噪声，通过y^T＝R*x^T得到包含竖直方向条带噪声的观测图像，T为矩阵转置运算；

(1)对于给定图像的条带噪声所在图像中行序号的集合Ω，

(2)设置观测矩阵R为一个大小为N_r×N_r的单位矩阵；

(3)处理任意一个行序号ω，

循环：

设置R(ω,ω)＝0

(4)得到包含条带噪声图像的观测矩阵R；

步骤二、建立包含条带噪声图像Curvelet域的稀疏表示方法，其包括：

(1)Curvelet域双变量阈值函数：

其中γ_i表示当前系数，γ_i-1表示当前系数的父系数，(g)₊运算表示为

σ_ε为条带噪声信号系数均方差，σ为图像有效信号系数均方差，τ为收敛控制因子；

Curvelet域双变量阈值函数中的参数估计：采用最高尺度子带的中值估计作为σ_ε的估计值，以当前系数γ_i为中心的3*3邻域窗口的均方差作为σ的估计值，通过迭代过程中逐渐收缩控制因子τ，逼近原始图像，提高去噪质量；

(2)Curvelet域相邻尺度系数γ_i和γ_i-1对应关系确定方法：

建立变换系数父子对应关系，设当前系数γ_i所在的子带矩阵

尺寸为m_i行，n_i列，l为方向，i为尺度，k为位置；设父系数γ_i-1所在子带矩阵为

的尺寸为m_i-1行，n_i-1列；m_i与m_i-1的关系存在m_i＞m_i-1和m_i≤m_i-1两种情况，n_i和n_i-1的关系也存在n_i＞n_i-1和n_i≤n_i-1两种情况；

当m_i＞m_i-1时需要行扩展

矩阵，使得m_i＝m_i-1的处理方法：

(1)令

(2)如果p＞1

对于任意一个c∈{1,2…p}循环：

插值扩展

(3)令q＝mod(m_i,m_i-1)

(4)如果q＞1

对于任意一个c∈{1,2…q}循环：

对称扩展

(5)得到扩展之后的矩阵

截取1:m_i行

得到新的

同样的方式处理当m_i≤m_i-1；当n_i＞n_i-1时；当n_i≤n_i-1时三种情况；

同理，在处理完n_i和n_i-1的情况后，得到的矩阵

和

尺寸相同，系数之间存在1:1的对应关系；

步骤三、建立条带噪声去除正则化问题最优化的求解算法；

(1)给定迭代停止参数δ，收敛控制因子τ，最大迭代次数K，包含条带噪声图像观测数据y；

(2)设置迭代计数器k＝0，Curvelet变换算子为W，Curvelet变换级数为z，迭代初值为x⁽⁰⁾＝y，Curvelet变换系数为α,变换尺度为i的Curvelet变换系数为α_i；

(3)条带噪声去除处理过程：

循环：

对x^k进行3*3的维纳窗口滤波

凸投影操作x^(k+1)＝x^k+R^T(y-R^T(x^k))

Curvelet正变换：α^(k+1)＝W(x^(k+1))

对于任意一个变换尺度i∈{1,2…z}循环：

双变量阈值函数滤波：

Curvelet反变换：x^k+1＝W^T(α^(k+1))

凸投影操作x^(k+1)＝x^(k+1)+R^T(y-R^T(x^(k+1)))

指数模型收缩τ

k＝k+1

直到：||x^(k+1)-x^(k)||₂＜＝δ或者k≥K循环结束；

x＝x^(k+1)；

(4)得到去除条带噪声的图像x。

Images