CN103927730B - 基于Primal Sketch修正及矩阵填充的图像降噪方法 - Google Patents

Abstract

本发明涉及基于PrimalSketch修正及矩阵填充的图像降噪方法，至少包括：步骤101：采用PrimalSketch模型获取加噪图像I₁的素描图P；步骤102：提出两阶段素描线修正规则并修正素描图P，对修正后的素描图P₂中每个素描线段上的素描点沿其所在素描线段的方向设计方向窗口，得到区域图。将图像I₁按照区域图映射为结构区域E₁和非结构区域E₂；步骤103：同时利用图像块方差统计的方法将图像I₁划分为光滑区域E₃与非光滑区域E₄；步骤104：将与非结构区域E₂重叠的光滑区域和非光滑区域分别记为光滑区域E₂₃和非光滑区域E₂₄；步骤105：对E₁和E₂₄采用基于矩阵填充的图像降噪方法进行降噪处理，而E₂₃采用非局部均值方法降噪处理；步骤106：合并E₁、E₂₃及E₂₄的降噪结果，得到图像I₁最终的降噪结果图。

Description

基于Primal Sketch修正及矩阵填充的图像降噪方法

技术领域

本发明属于图像处理技术领域，特别是基于Primal Sketch修正及矩阵填充的图像降噪方法，涉及Primal Sketch模型、矩阵填充方法及自然图像的降噪方法，可应用于自然图像的降噪处理。

背景技术

现存的图像降噪算法有许多种，虽然它们可以起到降噪的效果，但几乎每种方法都存在着自身的优势与缺陷，例如非局部均值降噪方法相比传统图像降噪方法具有更好的降噪效果，但其算法复杂度较大，并且许多不相似像素参与加权平均会引入不可忽视的累加误差；而三维块匹配降噪方法(BM3D)利用图像自身相似性能够较好的保持边缘与纹理等重要信息，但却在图像光滑区域的人工痕迹比较明显。通常情况下，一幅图像存在着若干不同的区域，并且这些区域有着不尽相同的特性，因此研究不同区域的特性对于图像处理具有非常重要的意义。大量的实践与理论表明：对于图像的不同区域采用不同的降噪方法进行处理可以达到更好的降噪效果。

目前已有一类基于区域划分的图像降噪方法，其主要思路如下：给定一幅加噪图像，使用Primal Sketch结构提取的方法将其划分为结构区域与非结构区域，同时采用图像块方差统计的方法将其划分为光滑区域与非光滑区域，将与非结构区域重叠的光滑区域和非光滑区域分别作为最终划分出的光滑区域和非光滑区域，并保持结构区域的初始划分结果，则该加噪图像被划分为结构区域、光滑区域及非光滑区域。对于光滑区域采用基于均值块的非局部均值方法，而对于结构区域、非光滑区域均采用三维块匹配降噪方法(BM3D)，最后将三个区域的降噪结果合并起来，得到最终的降噪图像。该方法的提出者在相关论文中已表明其降噪的良好效果，然而，这类方法在现阶段仍存在如下所述的不足之处：

(1)该类方法使用Primal Sketch标记图像的结构区域，而结构区域标记的准确度将会影响光滑区域、非光滑区域的划分结果，进一步将会影响最终合并而成的降噪结果。然而Primal Sketch对于噪声的反应较为敏感，虽然该类方法针对这一问题采用基于小波变换的图像重构进行预降噪处理，但结构区域标记的结果仍不理想，通过与实际情况对比不难发现，某些标记出现不同程度的偏差，准确度有所欠缺。再者，针对噪声图像首先进行预降噪处理虽然在一定程度上抑制噪声对Primal Sketch的影响，但亦会导致结构区域近邻区域的模糊化，使得最终的降噪效果不尽人意。由此可见，利用Primal Sketch进行结构区域标记的准确度有待提高，且预降噪处理带来的模糊化问题有待解决。

(2)在实际应用中，对于某一降噪算法的评价不仅需要考虑其降噪效果，其计算效率与运算耗时也是非常重要的评价指标。该类方法采用三维块匹配降噪方法(BM3D)处理结构区域、非光滑区域，BM3D利用图像自身的相似性，不仅使得降噪结果具有较高的信噪比，而且对纹理和边缘的保持也取得较佳的效果。但与此同时，BM3D存在着步骤繁琐、计算量较大及运算耗时较长等缺点。因此，BM3D在执行效率方面存在的问题将会影响到该类方法的整体执行效率，这也将促使我们寻求某种不仅在降噪效果方面可以与BM3D相媲美而且在执行效率方面优于BM3D的降噪方法。

发明内容

本发明的目的在于克服现有的基于区域划分的图像降噪方法的两处不足，首先针对上述不足之一，提出一种针对噪声图像结构区域标记的两阶段素描线修正规则，并将该修正规则应用于基于区域划分方法中结构区域标记的修正，从而提高结构区域标记的准确度，进而改善最终的降噪结果。再者，针对上述不足之二，提出一种基于矩阵填充的图像降噪方法，采用该方法处理已划分出的结构区域、非光滑区域，仿真实验的结果表明该方法不仅能够同BM3D降噪方法那样较好的保持边缘与纹理等信息，且执行效率较BM3D得到很大的提升。

为实现上述目的，本发明的主要思路是，基于Primal Sketch修正及矩阵填充的图像降噪方法，具体步骤包括如下：

步骤101：输入一幅自然图像I，添加强度为δ的高斯白噪声N，得到加噪图像I₁，当δ不大于30时，对图像I₁直接施以Primal Sketch得到其初始素描图；而当δ的取值大于30时，首先采用基于小波变换的图像重构进行预降噪处理得到重构图像I₂，然后对I₂施以PrimalSketch得到其初始素描图，初始素描图记为P；

步骤102：对初始素描图P施以两阶段素描线修正规则：阶段一、在初始素描图P中找出长度大于阈值Thr1、且包含的线段个数大于2的弯曲素描线{C_i,i＝1,2,...,n}，其中C_i表示某一弯曲素描线，修正这些素描线中若干线段交汇处突变比较尖锐的部分，生成中间结果P₁；阶段二、在中间结果P₁的基础上进行进一步处理，即修正那些长度大于阈值Thr2、整体接近直线的素描线，最终生成的结果即为修正后的素描图P₂；对修正后的素描图P₂中每个素描线段上的素描点沿其所在素描线段的方向设计大小为7*7的方向窗口，得到图像I₁的区域图R。将图像I₁按照区域图R映射为结构区域E₁和非结构区域E₂；

步骤103：将图像I₁分割为若干大小为w*w(w通常取值为11)的正方形窗口{W_h,h＝1,2,...,n}，其中W_h表示某一正方形窗口，使用方差统计的方法依据下述公式判断某一窗口W_h是否为光滑区域：

其中表示窗口W_h矩阵的标准差，参数θ取值范围为5～10。

依据该公式将图像I₁划分为光滑区域E₃与非光滑区域E₄；

步骤104：将光滑区域E₃、非光滑区域E₄与非结构区域E₂的重叠区域分别记为光滑区域E₂₃、非光滑区域E₂₄，并保持结构区域的初始划分结果E₁，则图像I₁最终被划分为结构区域E₁、光滑区域E₂₃及非光滑区域E₂₄；

步骤105：对结构区域E₁及非光滑区域E₂₄采用基于矩阵填充的图像降噪方法进行降噪处理，降噪结果分别记为D₁、D₂；而光滑区域E₂₃采用非局部均值方法降噪处理，降噪结果记为D₃；

步骤106：合并结构区域E₁、非光滑区域E₂₄及光滑区域E₂₃的降噪结果D₁、D₂、D₃，从而得到图像I₁最终的降噪图像D₄。

所述的步骤102，包括如下步骤：

步骤201：素描线修正阶段一的具体过程为：首先求解阈值Thr1，统计出由PrimalSketch生成的素描图中所有素描线的长度，然后从1开始每隔5个像素作为一个分组项，统计出分组集合{L_l,l＝1,2,...,n}，并求出该分组中素描线总个数最多的分组，记该分组的长度区间为[L,L+4]，取该区间范围的平均值作为阈值Thr1，即取值L+2。阈值Thr1的求解式如公式(1)所示：

其中，max{L₁,L₂,...,L_n}表示分组集合{L_l,l＝1,2,...,n}中素描线总个数最多的分组项，表示长度区间的左边界。

本方法将素描图中弯曲的素描线划分为两类，一类呈现为多重波动的曲线，另一类则呈现为或凸或凹的简单弧线；第一类素描线中线段的方向具有多变性，若整体考虑则很难提取出用于修正的潜在规则，将此类素描线划分为若干第二类或凸或凹的简单弧线，然后按照简单弧线的修正规则进行逐一处理。

对于某一素描线，统计其所有素描线段的方向{D_d,d＝1,2,...,n}(方向范围为[0，π])，取其中的最小值、最大值，分别记作DL、DU，如果DL与DU相差小于45度且方向值集中在区间[DL，DU]某一极小的范围时(据统计，该范围取值为[(DL+DU)/2-π/12，(DL+DU)/2+π/12]时出现直线误判现象的概率较小)，则判定该素描线接近于直线，并将这一类素描线记录下来，留至下一步骤中使用。若DL与DU相差不小于45度且方向值大致均匀的分布在区间[DL，DU]，即当该素描线的整体方向呈现出非严格递增的趋势时，则判定其属于第二类曲线。某一素描线是否接近于直线以及是否为第二类曲线的判定条件分别如下列公式(2)、(3)所示：

D_t∈[(DL+DU)/2-π/12，(DL+DU)/2+π/12]，t＝1,2,...,n且DU-DL|＜π/4 (2)

其中DL＝min{D_d,d＝1,2,...,n}，DU＝max{D_d,d＝1,2,...,n}。

对于被判定为弯曲且整体方向呈现出非严格递增的趋势的素描线，首先找出素描线中的线段交汇处突变比较尖锐的部分，求解过程为：计算出某一素描线中所有线段之间的夹角{A_a,a＝1,2,...,n}，另外，弯曲且整体方向呈现出非严格递增的趋势的素描线中所有线段之间的夹角中的多数应该是大于90度的钝角，而夹角较小的线段交汇处的突变比较尖锐，本方法规定夹角小于30度的线段交汇处的突变比较尖锐，并称其为尖锐夹角，统计出所有满足这一条件的尖锐夹角SA，其中SA的求解式如下：

SA＝{A_b||A_b|≤π/6,b＝1,2,...,n} (4)

对于尖锐夹角A_b逐一进行以下操作：获取构成尖锐夹角A_b的两条线段S₁、S₂，然后判断S₁、S₂是否为该素描线的首线段或尾线段：若S₁、S₂均非该素描线的首线段或尾线段，进一步获取分别与S₁、S₂相邻的线段S₃、S₄，线段S₃、S₄的方向角分别记为α、β，然后将线段S₁、S₂合并为线段S′₁，其中以线段S₄的终点作为S′₁的起点ps，以方向角α、β的平均值γ作为S′₁的方向角，通过线段栅格化使其不断逼近线段S₃的起点，最终使得线段S₃的起点作为线段S′₁的终点。若S₁、S₂存在其一是该素描线的首线段或尾线段，则将线段S₁、线段S₂合并为线段S′₁，其中S′₁的起点为线段S₁的起点ps，终点为线段S₂的终点pe，并以ps与pe的连线作为方向角。

至此完成第一阶段修正处理，生成Primal Sketch的中间结果P₁；

步骤202：素描线修正阶段二的具体过程为：首先求解阈值Thr2，其过程与步骤201中求解阈值Thr1的过程相似，在此从略。第一阶段中已记录所有接近直线的素描线，针对这类接近直线的素描线，逐一进行以下修正过程：

统计其所有线段的主导方向值{DD_z,z＝1,2,...,n}，根据主导方向值将接近直线的素描线划分为三类：接近水平的素描线(主导方向值紧邻0度)、接近垂直的素描线(主导方向值紧邻180度)、倾斜但接近直线的素描线(除去主导方向值紧邻0度、180度之外的)。

对接近水平的素描线，找出那些沿着非主导方向(即次要方向)的线段，逐一进行如下操作：从该线段出发，直至遇到第一个沿着主导方向的线段停止，将这些连续的、沿着非主导方向的线段一并处理，并将其称为非主导线段组。若某一素描线包含线段{S_k,k＝1,2,...,n},而正在处理的非主导方向的直线为S_g,遇到的第一个沿着主导方向的线段为S_j,则非主导线段组T表示为：

T＝{S_g,S_g+1,...,S_j-1} (5)

如果T中的第一条线段为S_g为该素描线的第一条线段，则取下一条线段的起始端点的纵坐标作为非主导线段组的修正基准，记为Y；其他情况下取上一线段的终止节点的纵坐标作为非主导线段组的修正基准。

对非主导线段组进行以下修正操作：保持这些非主导线段的横坐标不变，在保证整个素描线的形状接近直线的情况下，根据上述得到的修正基准Y，将非主导线段的纵坐标设置为Y+1或Y-1，这样就完成了单一非主导线段组的修正。

另外两类素描线的修正过程与上述过程类似，在此从略。最终生成的结果即为修正后的素描图；

所述的步骤105，包括如下步骤：

步骤301：对于光滑区域E₂₃采用结合结构信息的非局部均值图像去噪方法(ZL201110091450.2)中提出的非局部均值降噪方法，首先以光滑区域内的每一个像素点p为中心作一个大小为5×5的邻域均值，得到一幅均值图v，然后采用下面的公式进行降噪处理：

其中，p_i+k,j+l为第i+k行第j+l列的待去噪的像素点，k,l＝-f,...,0,...,+f，f为控制窗口大小的参数，为待去噪点p_i+k,j+l的估计值，S_p'为以像素p'为中心的搜索窗，q'为S_p'中的像素点，v(N_p')表示以p'为中心点的矩阵，v(N_q')表示以q'为中心点的矩阵，β为衰减参数，v(q'_i+k,j+l)为q'_i+k,j+l的像素值，Z(p')为归一化系数，由下述公式所表示

光滑区域E₂₃降噪结果记为D₃；

步骤302：非光滑区域E₂₄采用基于矩阵填充的图像降噪方法，该方法在整体上可以分为三个步骤：块匹配、矩阵填充及合成结果图像，具体步骤参见下文；

步骤303：对于结构区域E₁，首先沿着像素所在的边缘方向选取相似块，将这些相似块并为合成矩阵，然后采用基于矩阵填充的图像降噪方法对合成矩阵进行降噪处理，降噪方法的具体步骤参照步骤302，结构区域E₁降噪结果记为D₁。

所述的步骤302，包括如下步骤：

步骤401：首先将非光滑区域E₂₄划分为若干8*8的块矩阵，且块矩阵之间的步长为4，然后逐块使用一种四步骤搜索算法(4SS)对其进行相似块匹配操作，将当前正在处理的块矩阵称作参考块，表示为p_i,j，其中i、j为参考块的左上角坐标，以p_i,j左上角坐标作为中心点在Ns*Ns的局部区域进行相似块匹配操作。在此通过计算块矩阵间的平均绝对误差MAD(Mean Absolute Difference)判断两个块矩阵是否具有相似性，其中

C_xy表示当前正在进行相似块判断的块矩阵，R_xy表示参考块，τ_match表示判断块之间是否相似的最大阈值，本方法中取值为300。若MAD<τ_match，则断定这两个块矩阵相似，否则断定为不相似。假设与参照块p_x,y相似的块矩阵有M个，表示为将每个块逐列排放成一列向量，然后合并这些列向量，表示为一n²×m的合成矩阵P_x,y：

P_x,y＝(p_1,x,y,p_2,x,y,...,p_M,x,y)

其中p_1,x,y表示参考块p_x,y，而参考块p_x,y的相似块匹配结果即为合成矩阵P_i,j；

步骤402：对于合成矩阵P_x,y，仅保留那些非常可信的元素而除掉那些不可信的元素，不可信元素的判断准则为：若某一元素的数值远远偏离其所在行的所有元素的平均值，则该元素不可信；不可信元素源自两方面：其一是那些受到较大幅度的高斯噪声干扰的像素，其二来自于那些误匹配的块矩阵。从P_x,y的不完整版本P_x,y|Ω中恢复出Q_x,y由以下公式的最小化问题进行解决：

其中，n₁×n₂为块矩阵的大小n₁＝n²，n₂＝M，Ω表示所有可信元素的集合，p表示可信集Ω的大小与块矩阵中所有像素个数的比值；

步骤403：合并降噪后的块矩阵作为非光滑区域的最终降噪结果D₂，由于块矩阵间相互重叠，每个像素点可能被若干降噪后的块矩阵包含，取平均值作为其最终降噪后的像素值。

本发明与现有技术相比，具有如下优点：

1、条件性去除降噪预处理，区域划分的准确度有较大提升

本发明将一种针对噪声图像结构区域标记的两阶段素描线修正规则应用于现有方法中结构区域标记的修正，较大的提升了结构区域标记的准确度，进而使得区域划分的结果更为准确，最终的降噪结果也因此得以改善。另外，相对现有方法，当噪声强度不大于30时，本发明无需对噪声图像采用基于小波变换的图像重构进行预降噪处理，而是直接对其施以Primal Sketch，通过素描线修正的方式降低噪声对结构区域标记的影响，同时亦消除预降噪处理带来的结构区域标记近邻区域的模糊化问题。

2、执行效率较为高效，且降噪效果显著

本发明采用基于矩阵填充的图像降噪方法处理非光滑区域、结构区域，该方法通过将相似块合成矩阵的降噪问题转化为低秩矩阵的填充问题以达到降噪的目的，较大程度上改善了现有的基于区域划分的图像降噪方法中采用BM3D所存在的步骤繁琐、计算量大及计算耗时长等不足之处。该降噪方法不仅在降噪效果方面很好的保持边缘与纹理细节，而且在执行效率方面优于BM3D。

下面结合附图对本发明实施例做进一步描述。

附图说明

图1是本发明的流程图；

图2是本发明中两阶段素描线修正的流程图；

图3是本发明中的Lena图像；

图4是本发明中对图3进行Primal Sketch生成的初始素描图；

图5是本发明中对图3添加δ为25的高斯白噪声得到的Lena加噪图像；

图6是本发明中由图5进行Primal Sketch生成的初始素描图；

图7是本发明中由图7进行两阶段素描线修正后的素描图；

图8是本发明中图7结构区域的最终标记示意图；

图9是本发明中由图5划分出的光滑区域、非光滑区域与结构区域；

图10是本发明中图5最终的降噪结果图；

图11是本发明中的人工合成图像一、二及其初始素描图；

图12是本发明中的人工合成图像一、二及其初始素描图修正结果；

图13是本发明中素描线长度统计直方图；

图14是本发明中某一素描线所有线段的方向统计直方图；

图15是本发明中尖锐夹角的示意图；

图16是本发明中对于图15中尖锐夹角的修正结果示意图；

图17是本发明中δ为25的Lena加噪图像的降噪对比图；

图18是本发明中δ为45的Lena加噪图像的降噪对比图；

图19是本发明中图17降噪对比图对应的局部放大图；

图20是本发明中图18降噪对比图对应的局部放大图；

图21是本发明中δ为45的barbara加噪图像的降噪对比图；

图22是本发明中图21降噪对比图对应的局部放大图；

具体实施方式

本发明的主要思路是，对于一幅加噪图像，当其噪声强度δ的取值不大于30时，对该图像直接施以Primal Sketch得到其初始素描图；而当δ的取值大于30时，首先采用基于小波变换的图像重构进行预降噪处理然后再施以Primal Sketch得到其初始素描图。将本发明提出的两阶段素描线修正规则应用于上述初始素描图的修正，对修正后的素描图中每个素描线段上的素描点沿其所在素描线段的方向设计方向窗口，得到该图像的区域图。将该图像按照区域图映射为结构区域和非结构区域，同时采用图像块方差统计的方法将该图像划分为光滑区域与非光滑区域。将与非结构区域重叠的光滑区域和非光滑区域分别作为最终划分出的光滑区域和非光滑区域，并保持结构区域的初始划分结果，则该加噪图像被划分为结构区域、光滑区域及非光滑区域。对于光滑区域采用基于均值块的非局部均值方法，而对于结构区域、非光滑区域均采用本发明提出的基于矩阵填充的降噪方法来取代现有方法所采用的三维块匹配降噪方法(BM3D)，该降噪算法的核心思想在于将。最后将这三个区域的降噪结果合并起来，从而得到最终的降噪图像。

参照图1，本发明的具体实施步骤如下：

步骤101，对于图3所示的Lena图像施以Primal Sketch得到其初始素描图，如图4所示。对于图3所示的Lena图像施以δ为25的高斯白噪声(为简化叙述，在此略去施以δ为45的高斯白噪声时所对应的具体实施步骤，仅在下文中阐明其降噪结果)，结果如图5所示。对图5所示的Lean加噪图像直接施以Primal Sketch得到其初始素描图，如图6所示。

其中提及的Primal Sketch，参见Cheng-en Guo等人于2007年发表在计算机视觉与图形理解杂志上的文章《Primal Sketch：Integrating Texture and Structure》；

步骤102，参照图2，对图6所示的初始素描图施以下述两阶段素描线修正规则，得到修正后的素描图，如图7所示。对图7所示的修正素描图中每个素描线段上的素描点沿其所在素描线段的方向设计大小为7*7的方向窗口，得到与图4相对应的区域图，如图8所示，然后将图5按照图8所示的区域图映射为结构区域和非结构区域；

其中，两阶段素描线修正规则的制定思路及其实施的具体步骤如下所述：人工合成两幅大小为512*512的图像，其一含有宽度分别为1、3、5、7、9、11像素的较短水平线、较短垂直线、较长水平线、较长垂直线、倾角为45度直线以及倾角为135度直线，如图11(a)所示；其二含有宽度分别为1、3、5、7、9、11像素的较短曲线、较长曲线，如图11(b)所示。对于图11(a)、图11(b)，分别施以δ为10、15、25、30、35、40、45的高斯白噪声，然后在这些噪声图像上直接施以Primal Sketch，通过这些仿真实验与结果分析发现，由Primal Sketch生成的初始素描图存在如下问题：噪声越大，素描图的效果就越差，接近直线的素描线的偏离度(统计出同一素描线中所有线段的主导方向，在此把其他次要方向的线段相对主导方向的线段的偏离程度简称为偏离度)越明显；较长且弯曲的素描线中的若干线段交汇处的突变比较尖锐。这样的问题将会影响到Primal Sketch对噪声图像进行结构区域标记的准确度，进一步影响最终的降噪效果。针对上述问题，我们对素描图中受到噪声影响较显著地素描线进行如下两阶段修正：

1a)找出长度(即像素个数)大于阈值Thr1、所包含的线段个数大于2的弯曲素描线，修正这些素描线中若干线段交汇处的突变比较尖锐的部分，生成中间结果。

结合图13进一步阐明这一过程，该图中所示的直方图的横坐标表示素描线的长度，取值范围为1～512(在此假设单一素描线的长度不超过图像2的宽度或高度)，纵坐标表示某一相应素描线长度出现的次数，据统计，长度位于区间[1，5]之间的素描线总个数最多，因此，在此选取平均值3作为阈值Thr1。

结合图14进一步阐明：其中左图表示某一素描图中的其中一条弯曲素描线，右图表示该素描线中所有线段方向的统计直方图，其中DL为45度，DU为135度，方向值均匀的分布在区间[DL，DU]上，因此可以判定该素描线是第二类型的曲线。

结合图15、图16阐明弯曲素描线的修正过程：如图14所示，线段S₁、S₂之间的夹角可以看做尖锐夹角，分别求取相邻线段S₃、S₄的方向角α、β，将线段S₁、S₂合并为S′₁，图16即为图15所示的尖锐夹角修正后的结果示意图。

至此完成第一阶段修正处理，生成Primal Sketch的中间结果。

1b)在第一阶段生成的中间结果的基础上进行后处理，即修正长度大于阈值Thr2、整体接近直线的部分，生成的结果即为最终的素描图。

以图12为例进一步说明上述两阶段修正规则：图12(a)、图12(b)分别为11(a)、图11(b)施以δ为25的高斯白噪声得到的加噪图像，图12(c)、图12(e)分别为图12(a)在未施以及施以修正规则时的结果图，将12(c)、图12(e)对比可知，第二阶段的修正规则取得较好的修正效果；图12(d)、图12(e)分别为图12(b)在未施以及施以修正规则时的结果图，将图12(d)、图12(e)对比可知，第一阶段的修正规则取得较好的修正效果。

步骤103，对图5所示的Lena加噪图像采用图像块方差统计的方法将其划分为光滑区域、非光滑区域。

步骤104，将与非结构区域重叠的光滑区域和非光滑区域分别作为最终划分出的光滑区域和非光滑区域，并保持结构区域的初始划分结果。至此，图5所示的Lena加噪图像已完成光滑区域、结构区域及非光滑区域的划分，如图9所示，其中白色标记的区域为结构区域，黑色标记的部分为光滑区域，剩余的其他区域作为非光滑区域。

步骤105，对于图5所示的Lena加噪图像已划分出的光滑区域、结构区域及非光滑区域分别采用与其相应的降噪方法：

5a)图5的光滑区域采用基于均值块的非局部均值降噪方法；

5b)图5的非光滑区域采用基于矩阵填充的图像降噪方法；

5c)图5的结构区域采用基于矩阵填充的图像降噪方法；

步骤106，合并图5中光滑区域、结构区域及非光滑区域的降噪结果，得到最终的降噪图像，如图10所示。

本发明的优点由以下仿真的条件、内容及结果进一步说明：

1、仿真条件

条件1，用于制定两阶段素描线修正规则的仿真实验中使用两幅大小为512*512的人工合成图像，其一含有宽度分别为1、3、5、7、9、11像素的较短水平线、较短垂直线、较长水平线、较长垂直线、倾角为45度直线以及倾角为135度直线，如图11(a)所示；其二含有宽度分别为1、3、5、7、9、11像素的较短曲线、较长曲线，如图11(b)所示。

条件2，在验证两阶段素描线修正规则及基于矩阵填充的图像降噪方法有效性的仿真实验中使用Lena、barbara等图像。在本仿真实验中研究不同强度的高斯白噪声下的降噪效果，为了叙述的简便，在此仅以图3所示的Lena图像为例，对其添加δ＝25及δ＝45的高斯白噪声，Lena加噪图像如图4所示。

条件3，上述仿真实验均在单一噪声类型(高斯白噪声)情况下进行，且噪声强度δ的取值以30为分界，在不大于30的情况下无需对噪声图像进行降噪预处理。

2、仿真内容与结果

仿真1，通过对图12(a)、图12(b)进行Primal Sketch仿真实验与结果分析，制定出两阶段的素描线修正规则，由图12(f)可看出第一阶段的修正规则对于图11(a)中弯曲素描线的尖锐夹角的修正效果较好，由图12(e)亦可观察到第二阶段的修正规则对于图11(b)中接近直线的素描线中那些偏离其主导方向的线段已进行修正。

仿真2，对如图3所示的Lena图像施以δ为25的高斯白噪声，结果如图5所示。对图5所示的Lean加噪图像直接施以Primal Sketch得到其初始素描图，如图6所示。然后将两阶段的素描线修正规则应用于图6中初始素描图的修正，结果如图7所示。图8是图5所示的Lena加噪图像未使用两阶段的素描线修正规则得到的结构区域标记。然后，对图3所示的Lena图像施以δ为45的高斯白噪声，除首先采用基于小波变换的图像重构进行预降噪处理外，其他流程同上所述，在此从略。

仿真3，采用方差统计的方法将图5所示的Lena加噪图像划分为光滑区域、非光滑区域。

仿真4，将图5中与非结构区域重叠的光滑区域和非光滑区域分别作为最终划分出的光滑区域和非光滑区域，并保持结构区域的初始划分结果，则该加噪图像被划分为结构区域、光滑区域及非光滑区域。图9是图5所示的Lena加噪图像的最终区域划分结果，其中白色标记的区域为结构区域，黑色标记的部分为光滑区域，剩余的其他区域作为非光滑区域。

仿真5，对于图5所示的Lena加噪图像已划分出的光滑区域、结构区域及非光滑区域分别采用与其相应的降噪方法并保存其降噪结果。

仿真6，合并图5所示的Lena加噪图像中不同区域的降噪结果，得到其最终的降噪结果，如图10所示。

仿真7，为了与本发明方法做对比，分别采用三维块匹配(BM3D)、基于区域划分的降噪方法(RS)、非局部均值(NL)、基于矩阵填充的降噪方法(MC)、本发明方法五种降噪方法对Lena、barbara加噪图像进行降噪处理，其降噪结果的PSNR、SSIM如表1所示。结合表1数据及降噪结果附图可知，本发明方法较三维块匹配(BM3D)能够更好的保持边缘与纹理细节，改善非局部均值(NL)降噪结果存在的伪纹理现象，在基于区域划分的降噪方法(RS)基础上施以两阶段的素描线修正规则可起到提高边缘信息的准确性、改善降噪效果的作用；结合表1、表2，对比基于矩阵填充的降噪方法(MC)、三维块匹配(BM3D)的降噪结果及执行时间可知，基于矩阵填充的降噪方法(MC)不仅在降噪效果方面与三维块匹配(BM3D)相接近，且较大幅度的缩短降噪处理的执行时间、提升整体流程的执行效率。另外，本发明中采用基于矩阵填充的降噪方法替代现有方法中采用的三维块匹配(BM3D)，使得整体流程的执行效率得以提升。

表1五种不同方法降噪结果图的PSNR、SSIM比较

表2本发明与BM3D、MC在执行时间方面的比较

图17所示为δ等于25的Lena加噪图像的降噪对比图，其中图17(a)～(e)分别对应δ等于25的Lena噪声图像及本发明方法、BM3D、RS、NL、MC降噪结果。图18所示为δ等于45的Lena加噪图像的降噪对比图，其中图18(a)～(e)分别对应δ等于45的Lena噪声图像及本发明方法、BM3D、RS、NL、MC降噪结果。如图19、图20中各种降噪方法得到的降噪结果的局部放大图所示，通过对比本发明方法、基于区域划分的降噪方法(RS)、三维块匹配(BM3D)降噪方法得出的降噪结果的局部放大图(b)、图(c)及图(d)可知，两阶段的素描线修正规则可起到提高边缘信息的准确性、改善降噪效果的作用；对比局部放大图(b)、图(e)可知本发明可改善非局部均值(NL)降噪结果存在的伪纹理现象；而对比局部放大图(c)、图(f)可知基于矩阵填充降噪方法(MC)在降噪效果方面与三维块匹配(BM3D)相接近。图21为δ＝45的barbara加噪图像的降噪结果对比图，图22为图21中对比图的局部放大图。

另外，仿真实验的结果表明，本发明在同等强度的高斯白噪声下，相比基于区域划分的图像降噪方法，其降噪结果图的PSNR、SSIM通常情况下要大一些，这表明本发明提出的两阶段素描图修正规则能够起到提升区域划分准确度、改善最终降噪效果的作用。

综上所述，本发明是一种基于Primal Sketch修正及矩阵填充的图像降噪方法，与现有方法对比，不仅有效的改善区域划分的准确度，而且较大幅度的提升了执行效率，使降噪效果优于基于区域划分的降噪方法。

Claims (3)

1.基于Primal Sketch修正及矩阵填充的图像降噪方法，至少包括如下步骤：

步骤101：输入一幅自然图像I，添加强度为δ的高斯白噪声N，得到加噪图像I₁，当δ不大于30时，对图像I₁直接施以Primal Sketch得到其初始素描图；而当δ的取值大于30时，首先采用基于小波变换的图像重构进行预降噪处理得到重构图像I₂，然后对I₂施以PrimalSketch得到其初始素描图，初始素描图记为P；

步骤102：对初始素描图P施以两阶段素描线修正规则：阶段一、在初始素描图P中找出长度大于阈值Thr1、且包含的线段个数大于2的弯曲素描线{C_i,i＝1,2,...,n₃}，其中C_i表示某一弯曲素描线，修正这些素描线中若干线段交汇处突变比较尖锐的部分，生成中间结果P₁；阶段二、在中间结果P₁的基础上进行进一步处理，即修正那些长度大于阈值Thr2、整体接近直线的素描线，最终生成的结果即为修正后的素描图P₂；对修正后的素描图P₂中每个素描线段上的素描点沿其所在素描线段的方向设计大小为7*7的方向窗口，得到图像I₁的区域图R，将图像I₁按照区域图R映射为结构区域E₁和非结构区域E₂；

步骤103：将图像I₁分割为若干大小为w*w的正方形窗口{W_h,h＝1,2,...,n₄}，其中W_h表示某一正方形窗口，使用方差统计的方法依据下述公式判断某一窗口W_h是否为光滑区域：

$\sqrt{\mathrm{var}\left(W_h\right)}-q<\mathrm{\&theta;}$

其中表示窗口W_h矩阵的标准差，参数θ取值范围为5～10，依据该公式将图像I₁划分为光滑区域E₃与非光滑区域E₄；

步骤104：将光滑区域E₃、非光滑区域E₄与非结构区域E₂的重叠区域分别记为光滑区域E₂₃、非光滑区域E₂₄，并保持结构区域的初始划分结果E₁，则图像I₁最终被划分为结构区域E₁、光滑区域E₂₃及非光滑区域E₂₄；

步骤105：对结构区域E₁及非光滑区域E₂₄采用基于矩阵填充的图像降噪方法进行降噪处理，降噪结果分别记为D₁、D₂；而光滑区域E₂₃采用非局部均值方法降噪处理，降噪结果记为D₃；

步骤106：合并结构区域E₁、非光滑区域E₂₄及光滑区域E₂₃的降噪结果D₁、D₂、D₃，从而得到图像I₁最终的降噪图像D₄；

所述的步骤102，包括如下步骤：

步骤201：素描线修正阶段一的具体过程为：首先求解阈值Thr1，统计出由PrimalSketch生成的素描图中所有素描线的长度，然后从1开始每隔5个像素作为一个分组项，统计出分组集合{L_l,l＝1,2,...,n₅}，并求出该分组中素描线总个数最多的分组，记该分组的长度区间为[L,L+4]，取该区间范围的平均值作为阈值Thr1，即取值L+2，阈值Thr1的求解式如公式(1)所示：

其中，max{L₁,L₂,...,L_n6}表示分组集合{L_l,l＝1,2,...,n₅}中素描线总个数最多的分组项，表示长度区间的左边界；

将素描图中弯曲的素描线划分为两类，一类呈现为多重波动的曲线，另一类则呈现为或凸或凹的简单弧线；第一类素描线中线段的方向具有多变性，若整体考虑则很难提取出用于修正的潜在规则，将此类素描线划分为若干第二类或凸或凹的简单弧线，然后按照简单弧线的修正规则进行逐一处理；

对于某一素描线，统计其所有素描线段的方向{D_d,d＝1,2,...,n₇}，取其中的最小值、最大值，分别记作DL、DU，如果DL与DU相差小于45度且方向值集中在区间[DL，DU]某一极小的范围时，则判定该素描线接近于直线，并将这一类素描线记录下来，留至下一步骤中使用，若DL与DU相差不小于45度且方向值大致均匀的分布在区间[DL，DU]，即当该素描线的整体方向呈现出非严格递增的趋势时，则判定其属于第二类曲线；某一素描线是否接近于直线以及是否为第二类曲线的判定条件分别如下列公式(2)、(3)所示：

D_t∈[(DL+DU)/2-π/12，(DL+DU)/2+π/12]，t＝1,2,...,n₈且|DU-DL|＜π/4 (2)

且|DU-DL|≥π/4 (3)

其中DL＝min{D_d,d＝1,2,...,n₇}，DU＝max{D_d,d＝1,2,...,n₇}；

对于被判定为弯曲且整体方向呈现出非严格递增的趋势的素描线，首先找出素描线中的线段交汇处突变比较尖锐的部分，求解过程为：计算出某一素描线中所有线段之间的夹角{A_a,a＝1,2,...,n₉}，另外，弯曲且整体方向呈现出非严格递增的趋势的素描线中所有线段之间的夹角中的多数应该是大于90度的钝角，而夹角较小的线段交汇处的突变比较尖锐，夹角小于30度的线段交汇处的突变比较尖锐，并称其为尖锐夹角，统计出所有满足这一条件的尖锐夹角SA，其中SA的求解式如下：

SA＝{A_b||A_b|≤π/6,b＝1,2,...,n₁₀} (4)

对于尖锐夹角A_b逐一进行以下操作：获取构成尖锐夹角A_b的两条线段S₁、S₂，然后判断S₁、S₂是否为该素描线的首线段或尾线段：若S₁、S₂均非该素描线的首线段或尾线段，进一步获取分别与S₁、S₂相邻的线段S₃、S₄，线段S₃、S₄的方向角分别记为α、β，然后将线段S₁、S₂合并为线段S′₁，其中以线段S₄的终点作为S′₁的起点ps，以方向角α、β的平均值γ作为S′₁的方向角，通过线段栅格化使其不断逼近线段S₃的起点，最终使得线段S₃的起点作为线段S′₁的终点；若S₁、S₂存在其一是该素描线的首线段或尾线段，则将线段S₁、线段S₂合并为线段S′₁，其中S′₁的起点为线段S₁的起点ps，终点为线段S₂的终点pe，并以ps与pe的连线作为方向角；

至此完成第一阶段修正处理，生成Primal Sketch的中间结果P₁；

步骤202：素描线修正阶段二的具体过程为：首先用与步骤201求解阈值Thr1过程相似的步骤求解阈值Thr2，用接近直线的素描线，逐一进行以下修正过程：

统计其所有线段的主导方向值{DD_z,z＝1,2,...,n₁₁}，根据主导方向值将接近直线的素描线划分为三类：接近水平的素描线、接近垂直的素描线、倾斜但接近直线的素描线；对接近水平的素描线，找出那些沿着非主导方向的线段，逐一进行如下操作：从该线段出发，直至遇到第一个沿着主导方向的线段停止，将这些连续的、沿着非主导方向的线段一并处理，并将其称为非主导线段组；若某一素描线包含线段{S_k,k＝1,2,...,n₁₂},而正在处理的非主导方向的直线为S_g,遇到的第一个沿着主导方向的线段为S_j,则非主导线段组T表示为：

T＝{S_g,S_g+1,...,S_j-1} (5)

如果T中的第一条线段S_g为该素描线的第一条线段，则取下一条线段的起始端点的纵坐标作为非主导线段组的修正基准，记为Y；其他情况下取上一线段的终止节点的纵坐标作为非主导线段组的修正基准；

对非主导线段组进行以下修正操作：保持这些非主导线段的横坐标不变，在保证整个素描线的形状接近直线的情况下，根据上述得到的修正基准Y，将非主导线段的纵坐标设置为Y+1或Y-1，这样就完成了单一非主导线段组的修正；最终生成的结果即为修正后的素描图。

2.根据权利要求1所述的基于Primal Sketch修正及矩阵填充的图像降噪方法，所述的步骤105，包括如下步骤：

步骤301：对于光滑区域E₂₃采用结合结构信息的非局部均值图像去噪方法中提出的非局部均值降噪方法，首先以光滑区域内的每一个像素点p为中心作一个大小为5×5的邻域均值，得到一幅均值图v，然后采用下面的公式进行降噪处理：

其中，p_i+k,j+l为第i+k行第j+l列的待降噪像素点，k,l＝-f,...,0,...,+f，f为控制窗口大小的参数，为待去噪点p_i+k,j+l的估计值，S_p'为以像素p'为中心的搜索窗，q'为S_p'中的像素点，v(N_p')表示以p'为中心点的矩阵，v(N_q')表示以q'为中心点的矩阵，β为衰减参数，v(q'_i+k,j+l)为q'_i+k,j+l的像素值，Z(p')为归一化系数，由下述公式所表示

$Z\left(p^{\&prime;}\right)=\underset{q^{\&prime;}\&Element;S_{p^{\&prime;}}}{\&Sigma;}{e}^{-\left|\right|v\left(N_{p^{\&prime;}}\right)-v\left(N_{q^{\&prime;}}\right)|{|}^2/\mathrm{\&beta;}}$

光滑区域E₂₃降噪结果记为D₃；

步骤302：非光滑区域E₂₄采用基于矩阵填充的图像降噪方法，该方法在整体上可以分为三个步骤：块匹配、矩阵填充及合成结果图像；

步骤303：对于结构区域E₁，首先沿着像素所在的边缘方向选取相似块，将这些相似块并为合成矩阵，然后采用基于矩阵填充的图像降噪方法对合成矩阵进行降噪处理，结构区域E₁降噪结果记为D₁。

3.根据权利要求2所述的基于Primal Sketch修正及矩阵填充的图像降噪方法，所述的步骤302，包括如下步骤：

步骤401：首先将非光滑区域E₂₄划分为若干8*8的块矩阵，且块矩阵之间的步长为4，然后逐块使用一种四步骤搜索算法对其进行相似块匹配操作，将当前正在处理的块矩阵称作参考块，表示为p_x,y，其中x、y为参考块的左上角坐标，以左上角坐标作为中心点在Ns*Ns的局部区域进行相似块匹配操作；在此通过计算块矩阵间的平均绝对误差MAD来判断两个块矩阵是否具有相似性，其中

$MAD=\frac{1}{N^2}\underset{x=0}{\overset{N-1}{\&Sigma;}}\underset{y=0}{\overset{N-1}{\&Sigma;}}\left|C_{xy}-R_{xy}\right|$

C_xy表示当前正在进行相似块判断的块矩阵，R_xy表示参考块，τ_match表示判断块之间是否相似的最大阈值，本方法中取值为300；若MAD<τ_match，则断定这两个块矩阵相似，否则断定为不相似；假设与参照块p_x,y相似的块矩阵有M个，表示为将每个块逐列排放成一列向量，然后合并这些列向量，表示为一n²×M的合成矩阵P_x,y：

P_x,y＝(p_1,x,y,p_2,x,y,...,p_M,x,y)

其中p_1,x,y表示参考块p_x,y，而参考块p_x,y的相似块匹配结果即为合成矩阵P_x,y；

步骤402：对于合成矩阵P_x,y，仅保留那些非常可信的元素而除掉那些不可信的元素，不可信元素的判断准则为：若某一元素的数值远远偏离其所在行的所有元素的平均值，则该元素不可信；不可信元素源自两方面：其一是那些受到较大幅度的高斯噪声干扰的像素，其二来自于那些误匹配的块矩阵；从P_x,y的不完整版本P_x,y|Ω中恢复出Q_x,y由以下公式的最小化问题进行解决：

$\underset{Q_{x,y}}{min}\frac{1}{2}\left|\right|Q_{x,y}{|}_{\mathrm{\&Omega;}}-P_{x,y}{|}_{\mathrm{\&Omega;}}|{|}_F^2+\mathrm{\&mu;}|\left|Q_{x,y}\right|{|}_{*}$

其中，n₁×n₂为块矩阵的大小n₁＝n²，n₂＝M，Ω表示所有可信元素的集合，p表示可信集Ω的大小与块矩阵中所有像素个数的比值；

步骤403：合并降噪后的块矩阵作为非光滑区域的最终降噪结果D₂，由于块矩阵间相互重叠，每个像素点可能被若干降噪后的块矩阵包含，取平均值作为其最终降噪后的像素值。