本发明属于数字图像处理技术领域,涉及一种图像去噪方法,特别涉及一种基于超像素聚类和稀疏表示的图像去噪方法,可应用于图像分类、目标识别、边缘检测等要求对图像进行去噪预处理的场合。
背景技术:
由于受到成像设备和成像环境的限制,数字图像在采集、转换或传输的过程中不可避免地受到噪声的污染。噪声的存在使得图像质量下降,并影响到后续图像处理。为了获得高质量的图像,就必须对图像进行去噪处理。因此,图像去噪在图像处理领域占据着重要的地位。
随着国内外图像去噪技术不断发展,研究人员相继提出了许多图像去噪方法。目前图像去噪方法主要分为三类:空间域去噪方法,频率域去噪方法,稀疏变换域去噪方法。空间域去噪方法主要是利用局部窗口内像素灰度值的连续性来对当前像素点的灰度值进行调整,达到去噪的目的。该类去噪方法主要包括均值滤波,中值滤波、非局部均值滤波(non-localmeans,nlm)等,其中最经典的是nlm算法。nlm算法通过对相似图像块做加权平均来估计参考块的中心点,从而降低噪声,虽然nlm算法相比其它空间域去噪方法,取得了较好的去噪效果,但是峰值信噪比仍然较低,同时去噪后的图像边缘、纹理区域模糊。
频率域去噪方法主要是将图像从空间域变换到频率域,再对频率域系数进行处理,最后将频率域系数反变换到空间域,得到去噪后的图像,该类去噪方法主要包括小波变换去噪方法和多尺度几何分析。小波变换去噪方法缺少方向选择性,不适宜表示图像边缘、轮廓等线性奇异性的结构特征,且过于依赖阈值的选择,导致其去噪效果差。多尺度几何分析缺乏灵活性,对不同的结构特征需要选择不同的变换,而一幅图像含有多种不同的结构。
稀疏变换域去噪方法主要通过对含噪图像进行学习,得到能够反映图像特征的字典,然后利用得到的字典对图像进行重构,从而达到去噪的目的。这类去噪方法中比较经典的方法有k-svd算法。k-svd算法在提取的图像块中随机选取若干图像块作为训练样本,训练得到具有数据自适应性的字典,但由于随机选取若干图像块作为训练样本的操作忽视了图像的结构特征、边缘特征和纹理特征,导致得到的字典不能对图像的这些特征进行很好地描述,并且由于训练得到的字典存在噪声,导致稀疏分解得到的稀疏系数图像信息的描述不准确,最终导致去噪图像峰值信噪比较低,边缘、纹理等细节信息丢失,图像去噪效果差。
技术实现要素:
本发明的目的在于克服上述现有技术存在的缺陷,提出了一种基于超像素聚类和稀疏表示的图像去噪方法,用以解决现有图像去噪方法中存在的去噪图像峰值信噪比低和细节信息丢失的技术问题。
为实现上述目的,本发明采取的技术方案包括如下步骤:
步骤1,输入一幅含有标准方差为δ的高斯白噪声的图像in;
步骤2,首先设定图像in的超像素数目为r,并对图像in进行超像素分割,得到超像素集合{spi|i=1,2,...,r},其次定义一个空的相似矩阵s,计算超像素集合{spi|i=1,2,...,r}中每两个超像素
步骤3,设定类的个数为k,并利用相似矩阵s,对超像素集合{spi|i=1,2,...,r}中的超像素进行聚类,得到相似超像素集合{crk|k=1,2,...,k},其中k是相似超像素集合{crk|k=1,2,...,k}中相似超像素的序号,crk是相似超像素集合{crk|k=1,2,...,k}中第k簇相似超像素;
步骤4,对相似超像素集合{crk|k=1,2,...,k}中每簇相似超像素分别进行重叠取块,得到k个图像块子集合,再以该k个图像块子集合中的每个图像块子集合为元素组成图像块子集集合{{blkt|t=1,2,...,tk}|k=1,2,...,k},并将该k个图像块子集合进行合并,得到图像块集合{blkt|k=1,2,...,k;t=1,2,...,tk},其中,{blkt|t=1,2,...,tk}是图像块子集集合{{blkt|t=1,2,...,tk}|k=1,2,...,k}中第k个图像块子集合,t是从相似超像素集合{crk|k=1,2,...,k}中第k簇相似超像素crk中提取的图像块的序号,blkt是从相似超像素集合{crk|k=1,2,...,k}中第k簇相似超像素crk中提取的第t个图像块,tk是相似超像素集合{crk|k=1,2,...,k}中第k簇相似超像素crk中提取的图像块的数目;
步骤5,对图像块子集集合{{blkt|t=1,2,...,tk}|k=1,2,...,k}中的每个图像块子集合分别进行字典训练,得到字典集合{dk|k=1,2,...,k},其中,dk是字典集合{dk|k=1,2,...,k}中第k个字典;
步骤6,设迭代变量为
步骤7,设定选取相似图像块的数目l,为图像块集合{blkt|k=1,2,...,k;t=1,2,...,tk}中的每个图像块选取l个相似图像块,并计算图像块集合{blkt|k=1,2,...,k;t=1,2,...,tk}中每个图像块的l个相似图像块的稀疏系数加权和,得到加权稀疏系数集合
步骤7a,计算图像块子集合{blkt|k=1,2,...,k;t=1,2,...,tk}中图像块blkt与图像块子集合{blkt|t=1,2,...,tk}中除图像块blkt以外的其它图像块之间的相似度,再对得到的相似度按从大到小的顺序进行排序,从图像块子集合{blkt|t=1,2,...,tk}中选取前l个相似度对应的图像块作为图像块blkt的相似图像块,并对图像块集合{blkt|k=1,2,...,k;t=1,2,...,tk}中除图像块blkt以外的其它图像块进行相同的操作,得到相似度集合
步骤7b,利用相似度集合
步骤8,利用加权稀疏系数集合
其中,ykt表示将图像块blkt的灰度值矩阵进行列化得到的灰度值向量,γ是用以平衡图像块blkt重构误差和稀疏度的归一化参数;
步骤9,设定迭代变量阈值λ,并判断迭代变量
步骤10,利用字典集合{dk|k=1,2...,k}和稀疏系数集合
本发明与现有的技术相比,具有以下优点:
1.本发明由于在获取字典的过程中,通过对超像素进行聚类,并对相似超像素进行字典学习,有效地学习并利用了图像的结构特征、边缘特征、纹理特征以及非局部相似性,能够得到对图像的结构特征、边缘特征和纹理特征描述更加有效的字典,与现有技术相比,有效地提高了去噪图像的峰值信噪比,同时更好地保留了去噪图像的边缘、纹理等细节信息。
2.本发明由于利用相似图像块的加权稀疏系数对图像块稀疏分解过程进行约束,降低了字典中噪声对稀疏系数的影响,能够得到对图像信息描述更加准确的稀疏系数,与现有技术相比,进一步提高去噪图像的峰值信噪比,同时更加完整地保留了去噪图像的边缘、纹理等细节信息。
附图说明
图1是本发明的实现流程图;
图2是本发明仿真实验使用的十幅标准测试图像;
图3是本发明与现有技术对monarch图像的去噪效果对比图;
图4是本发明与现有技术对house图像的去噪效果对比图。
具体实施方式:
以下结合附图和具体实施例,对本发明作进一步详细说明:
参照图1,一种基于超像素聚类和稀疏表示的图像去噪方法,包括如下步骤:
步骤1,输入一幅含有标准方差为δ的高斯白噪声的图像in。
在本实施例中,采用的是分辨率为的512×512灰度图像。
步骤2,首先设定图像in的超像素数目为r,并对图像in进行超像素分割,得到超像素集合{spi|i=1,2,...,r},其次定义一个空的相似矩阵s,计算超像素集合{spi|i=1,2,...,r}中每两个超像素
步骤2a,超像素的个数r的设置不是固定的值,在本实施例中,超像素的个数设置为r=500,且对图像in进行超像素分割,可采用多种算法如简单线性迭代聚类算法(simplelineriteratorclustering,slic)算法,normalizedcut算法,mean-shift算法,quick-shift算法。本实例采用简单线性迭代聚类算法,相比其它超像素分割算法,该算法在运行速度、生成超像素的紧凑度、轮廓保持方面都比较理想,其实现步骤为:
步骤2a1,计算分割完成以后的每一个超像素的像素点数目估计值pn与边长估计值st,其中,
步骤2a2,在图像in平面内,以像素点为基本单位,以step为垂直方向和水平方向的步长,从第rw行像素点开始,均匀地选取r个聚类中心,得到聚类中心集合{cq|q=1,2,...,r},其中,step=st,
步骤2a3,在聚类中心集合{cq|q=1,2,...,r}中的聚类中心cq的ns×ns邻域内,计算每个像素点的梯度值,选取梯度值最小的像素点替换聚类中心集合{cq|q=1,2,...,r}中的聚类中心cq,对聚类中心集合{cq|q=1,2,...,r}中除聚类中心cq以外的其它聚类中心进行相同的操作,得到新的聚类中心集合{cq|q=1,2,...,r};
步骤2a4,设定迭代变量θ,并初始化为0,并在2st×2st的搜索窗口内,将像素点分配给与其距离最小的聚类中心,得到r簇相似像素点,其中计算任意像素点px=[g,x,y]t到任意聚类中心cx=[gc,xc,yc]的距离ds的公式为:
其中,g是像素点px的灰度值,x是像素点px在x轴方向的位置坐标值,y是像素点px在y轴方向的位置坐标值,gc是聚类中心cx的灰度值,xc是聚类中心cx在x轴方向的位置坐标值,yc是聚类中心cx在y轴方向的位置坐标值,
步骤2a5,分别计算每簇相似像素点的均值,作为每簇相似像素点的新的聚类中心,并更新聚类中心集合{cq|q=1,2,...,r};
步骤2a6,设定迭代变量阈值ω,并判断迭代变量
经验数据表明,只需迭代10次即可实现连续两次的聚类中心误差不超过5%,因此,本实例将迭代次数设置为10次。
步骤2b,上述的计算超像素集合{spi|i=1,2,...,r}中每两个超像素
步骤2b1,计算超像素集合{spi|i=1,2,...,r}中每个超像素的特征向量,得到超像素特征向量集合{ui|i=1,2,...,r},计算公式为:
其中,ui是超像素特征向量集合{ui|i=1,2,...,r}中的第i个特征向量,γi是超像素spi中包含的像素点的数目,j表示超像素spi中像素点的序号,且j=1,2,...,γi,fj表示超像素spi中第j个像素点的特征向量,且fj=[g,ix,iy,ixx,iyy,β×x,β×y]t,g表示超像素spi中第j个像素点的灰度值,ix,iy,ixx,iyy分别表示超像素spi中第j个像素点在x轴方向和y轴方向的一阶导数与二阶导数;x和y分别表示超像素spi中第j个像素点在x轴方向的坐标值与y轴方向的坐标值,β是位置特征与其它特征之间的平衡因子,其取值范围是(0,1];
本实例将β设置为0.5,且在求取位置坐标值时,在图像in平面上以中心像素点为原点,水平方向为x轴方向,垂直方向为y轴方向,建立坐标系。
步骤2b2,计算超像素集合{spi|i=1,2,...,r}中每个超像素的协方差矩阵,得到协方差矩阵集合{mi|i=1,2,...,r},计算公式为:
其中,mi是超像素spi的协方差矩阵,a和b分别是协方差矩阵mi中元素的行号和列序号,mi(a,b)是矩阵mi中第a行第b列的元素,且a=1,2,...,7,b=1,2,...,7,a′和b′是超像素spi中第j个像素点的特征向量fj中两个元素的序号,且a′=a,b′=b,fj(a′)是超像素spi中第j个像素点的特征向量fj中序号为a′的元素,fj(b′)是超像素spi中第j个像素点的特征向量fj中序号为b′的元素,a″和b″是超像素特征向量集合{ui|i=1,2,...,r}中的第i个特征向量中两个元素的序号,且a″=a′=a,b″=b′=b,ui(a″)是超像素特征向量集合{ui|i=1,2,...,r}中的第i个特征向量中序号为a″的元素,ui(b″)是超像素特征向量集合{ui|i=1,2,...,r}中的第i个特征向量中序号为b″的元素;
步骤2b3,计算超像素集合{spi|i=1,2,...,r}中任意两个超像素
其中,i1和i2是超像素集合{spi|i=1,2,...,r}中任意两个超像素的序号,且i1=1,2,...,r,i2=1,2,...,r,i1≠i2,
步骤2b4,将超像素相似度集合
其中,r1和r2是相似矩阵s中元素的行序号和列序号,r1=1,2,...,r,r2=1,2,...,r,s(r1,r2)是相似矩阵s中第r1行第r2列元素,
步骤3,设定类的个数为k,并利用相似矩阵s,对超像素集合{spi|i=1,2,...,r}中的超像素进行聚类,得到相似超像素集合{crk|k=1,2,...,k},其中k是相似超像素集合{crk|k=1,2,...,k}中相似超像素的序号,crk是相似超像素集合{crk|k=1,2,...,k}中第k簇相似超像素。
类的个数k的设置不是固定的值,在本实施例中,类的个数为k=40,且上述的对超像素进行聚类,可采用多种算法如近邻传播算法,k-means算法,谱聚类,稀疏子空间聚类算法,拉普拉斯稀疏子空间聚类算法,本实例采用拉普拉斯稀疏子空间聚类算法,该算法具有对噪声具有鲁棒性,对含噪数据聚类效果好的优点,实现步骤为:
步骤3a,利用相似矩阵s,计算对角矩阵e:
其中,r1′是对角矩阵e中对角元素的行序号和列序号,r1′=1,2,...,r,e(r1′,r1′)是对角矩阵e中第r1′行第r1′列的元素,r1和r2是相似矩阵s中元素的行序号和列序号,且r1=r1′,r2=1,2,...,r,s(r1,r2)是相似矩阵s中第r1行第r2列的元素;
步骤3b,利用相似矩阵s和对角矩阵e,计算拉普拉斯矩阵l,计算公式为:
l=e-s;
步骤3c,利用超像素特征向量集合{ui|i=1,2,...,r}和拉普拉斯矩阵l,计算超像素集合{spi|i=1,2,...,r}中每个超像素的稀疏系数,得到稀疏系数矩阵c,计算超像素稀疏系数的公式为:
其中,矩阵u={u1,u2,...,ur},矩阵
上述计算超像素稀疏系数的公式将拉普拉斯矩阵l引入到稀疏子空间聚类公式中,是为了使相似的超像素具有相似的稀疏系数,利用稀疏域的非局部相似性来提高稀疏系数的准确性,以达到更好的超像素聚类效果。
在本实施例中,λ=0.01,η=0.2。
步骤3d,对稀疏系数矩阵c进行更新,得到对称矩阵
其中,k1′和k2′是对称矩阵
步骤3e,建立无向图g,将超像素特征向量集合{ui|i=1,2,...,r}中的每个超像素特征向量作为无向图g的顶点,得到顶点集合{vi|i=1,2,...,r},且将对称矩阵
谱聚类算法对无向图g进行划分时,无向图g的邻接矩阵是对称矩阵
步骤4,对相似超像素集合{crk|k=1,2,...,k}中每簇相似超像素分别进行重叠取块,得到k个图像块子集合,再以该k个图像块子集合中的每个图像块子集合为元素组成图像块子集集合{{blkt|t=1,2,...,tk}|k=1,2,...,k},并将该k个图像块子集合进行合并,得到图像块集合{blkt|k=1,2,...,k;t=1,2,...,tk},其中,{blkt|t=1,2,...,tk}是图像块子集集合{{blkt|t=1,2,...,tk}|k=1,2,...,k}中第k个图像块子集合,t是从相似超像素集合{crk|k=1,2,...,k}中第k簇相似超像素crk中提取的图像块的序号,blkt是从相似超像素集合{crk|k=1,2,...,k}中第k簇相似超像素crk中提取的第t个图像块,tk是相似超像素集合{crk|k=1,2,...,k}中第k簇相似超像素crk中提取的图像块的数目。
图像块的边长p的设置不是固定的值,但是p应该是奇数,在本实施例中,图像块的边长p设置为7,且上述的对相似超像素集合{crk|k=1,2,...,k}中每簇相似超像素分别进行重叠取块,实现步骤为:
步骤4a,设定图像块边长p,在图像in平面内,以图像in的边界像素点为中心,镜像复制p′个像素点,得到图像i'n,其中,
步骤4b,在图像i'n平面内,以相似超像素集合{crk|k=1,2,...,k}中每簇相似超像中的像素点为中心,提取p×p大小的图像块,得到k个图像块子集合。
步骤5,对图像块子集集合{{blkt|t=1,2,...,tk}|k=1,2,...,k}中的每个图像块子集合分别进行字典训练,得到字典集合{dk|k=1,2,...,k},其中,dk是字典集合{dk|k=1,2,...,k}中第k个字典。
上述的对图像块子集集合{{blkt|t=1,2,...,tk}|k=1,2,...,k}中的每个图像块子集合分别进行字典训练,可采用多种字典学习算法,如小波基字典,k-svd算法,主成分分析算法等,本实例采用主成分分析算法,具有计算速度快,得到字典对数据具有自适应性的优点,实现步骤为:
步骤5a,计算图像块子集集合{{blkt|t=1,2,...,tk}|k=1,2,...,k}中的每个图像块子集合{blkt|t=1,2,...,tk}的特征矩阵pk,计算公式为:
其中,bk是相似超像素集合{crk|k=1,2,...,k}中第k簇相似超像素crk对应的图像块子集合{blkt|t=1,2,...,tk}的灰度值矩阵,
步骤5b,计算图像块子集集合{{blkt|t=1,2,...,tk}|k=1,2,...,k}中的每个图像块子集合{blkt|t=1,2,...,tk}对应的字典,得到字典集合{dk|k=1,2,...,k},计算公式为:
其中,
步骤6,设迭代变量为
对图像块集合{blkt|k=1,2,...,k;t=1,2,...,tk}中所有图像块进行稀疏分解,采用广义正交匹配追踪算法。
步骤7,设定选取相似图像块的数目l,为图像块集合{blkt|k=1,2,...,k;t=1,2,...,tk}中的每个图像块选取l个相似图像块,并计算图像块集合{blkt|k=1,2,...,k;t=1,2,...,tk}中每个图像块的l个相似图像块的稀疏系数加权和,得到加权稀疏系数集合
步骤7a,计算图像块子集合{blkt|k=1,2,...,k;t=1,2,...,tk}中图像块blkt与图像块子集合{blkt|t=1,2,...,tk}中除图像块blkt以外的其它图像块之间的相似度,再对得到的相似度按从大到小的顺序进行排序,从图像块子集合{blkt|t=1,2,...,tk}中选取前l个相似度对应的图像块作为图像块blkt的相似图像块,并对图像块集合{blkt|k=1,2,...,k;t=1,2,...,tk}中除图像块blkt以外的其它图像块进行相同的操作,得到相似度集合
相似图像块的数目l的设置不是固定的值,在本实施例中,相似图像块的数目l=10,过多的相似图像块的数目可能导致去噪后的图像边缘模糊,过少的相似图像块的数目导致相似图像块的稀疏系数加权和对图像块的稀疏分解过程影响小,且上述的计算图像块集合{blkt|k=1,2,...,k;t=1,2,...,tk}中图像块blkt与图像块子集合{blkt|t=1,2,...,tk}中除图像块blkt以外的其它图像块之间的相似度,计算公式为:
其中,τ是图像块集合{blkt|t=1,2,...,tk}中除图像块blkt以外的任意图像块的序号,且τ=1,2,...,tk,τ≠t,blkτ是图像块集合{blkt|t=1,2,...,tk}中除图像块blkt以外的任意图像块,ykt与ykτ分别是图像块blkτ与图像块blkτ对应的灰度值列向量,
步骤7b,利用相似度集合
上述的计算图像块集合{blkt|k=1,2,...,k;t=1,2,...,tk}中每个图像块的l个相似图像块的稀疏系数加权和,计算公式为:
其中,
步骤8,利用加权稀疏系数集合
其中,ykt表示将图像块blkt的灰度值矩阵进行列化得到的灰度值向量,γ是用以平衡图像块blkt重构误差和稀疏度的归一化参数。
在本实施例中,
步骤9,设定迭代变量阈值λ,并判断迭代变量
最大迭代次数t的设置不是固定的值,在本实施例中,最大迭代次数λ=10。
步骤10,利用字典集合{dk|k=1,2...,k}和稀疏系数集合
其中,
以下结合仿真实验,对本发明的技术效果作进一步描述:
1.仿真条件和内容:
采用matlabr2010a软件在配置为corei3-21203.30ghz,内存4g,windows764位操作系统的计算机上,使用本发明和现有技术对十幅标准测试图像进行去噪仿真实验,其中,本发明与现有技术对monarch图像和house图像的去噪效果对比结果如图3和图4所示。
2.仿真结果分析:
图2是本发明仿真实验的十幅标准测试图像,从左至右,从上到下,图像的名字依次lena,monarch,house,parrot,barbara,pepper,couple,cameraman,straw,man。本发明仿真实验分别对十幅标准测试图像添加高斯白噪声,得到人工合成的待去噪图像,并且使用峰值信噪比(peaksignaltonoiseratio,psnr)以及图像细节保留程度作为衡量去噪效果的指标。
参照图3,图3(a)是原始monarch图像,图3(b)是含有标准方差为20的高斯白噪声的待去噪monarch图像,图3(c)是nlm方法的去噪效果图,图3(d)是k-svd方法的去噪效果图,图3(e)是bm3d方法的去噪效果图,图3(f)是本发明的去噪效果图,且图3中每个图像的矩形框区域是图像的局部放大图。
由图3可以看出:与其它对比方法相比,本发明对图像的细节信息保留地更加完整,如局部放大图所示,本发明对蝴蝶的触角和蝴蝶翅膀上纹络的边缘保留地更加完整,更加清晰,证明了本发明方法可以实现更好的去噪效果。
参照图4,图4(a)是原始house图像,图4(b)含有标准方差为20的高斯白噪声的待去噪house图像,4(c)是nlm方法的去噪效果图,4(d)是k-svd方法的去噪效果图,图4(e)是bm3d方法的去噪效果图,图4(f)是本发明的去噪效果图,且图4中每个图像的矩形框区域是图像的局部放大图。
由图4可以看出:与其它对比方法相比,本发明对图像的细节信息保留地更加完整,如局部放大图所示,本发明对烟囱的边缘和排气管与房顶交界处细节信息保留地更加完整,更加清晰,证明了本发明方法可以实现更好的去噪效果。
为了进一步分析本发明与其它对比方法的去噪效果,表1给出了本发明方法与其它对比方法的对含有不同标准方差的高斯白噪声的标准测试图像进行去噪的psnr值的对比。表1中,第一行单元格中的数值是仿真实验添加的高斯白噪声的标准方差δ,第一列单元格中的内容是仿真实验的图像名称,且含有四个数值的单元格中,左上角的数值是nlm方法的psnr值,右上角的数值是k-svd方法的psnr值,左下角的数值是bm3d方法的psnr值,右下角的数值是本发明的psnr值。
表1本发明与现有技术对含有不同标准方差的高斯白噪声的标准测试图像进行去噪的psnr值
从表1可以看出对于含有不同标准方差的高斯白噪声的不同的标准测试图像进行去噪时,与其它对比算法相比,本发明的psnr值明显高于nlm算法和k-svd算法,且高于或者接近于bm3d算法。
从图3、图4以及表1可以看出证明了本发明方法可以实现比nlm算法和k-svd算法更好的去噪效果,且比bm3d算法更好或者相近的去噪效果。