基于局部多原子匹配追踪的图像稀疏分解快速方法

摘要

本发明公开一种基于局部多原子匹配追踪的图像稀疏分解快速方法，主要解决现有图像稀疏分解方法全局搜索和/或单原子选择导致的计算复杂度高、运算速度慢的问题。其实现步骤为：(1)生成并存储核原子；(2)稀疏分解初始化；(3)多原子匹配搜索；(4)图像残差更新；(5)分解结束判断；(6)局部搜索核原子集更新。其优点是将局部搜索和多原子选择两种匹配追踪方式相结合，最大限度地降低了运算复杂度；图像残差更新采用逐原子依次更新的方式，保证了稀疏分解的逼近质量。与传统的匹配追踪算法相比，本发明方法运算复杂度低5个数量级，与当今其他先进方法相比，本发明方法在运算速度和逼近性能上都有较明显的优势。

说明书

技术领域

本发明属于模式识别与图像处理领域，涉及一种基于局部多原子匹配追踪的图像稀疏分解快速方法，可应用于图像处理和计算机视觉中。

背景技术

以小波分析为基础，Mallat和Zhang提出了信号在超完备冗余字典中的稀疏表示和稀疏分解思想，此后该思想被推广到二维图像，并成功应用于图像处理的许多方面，如图像去噪、人脸识别、超分辨率重建、图像压缩编码等。图像稀疏分解问题是一个NP-hard组合搜索问题，直接求解十分困难，现有的方法可分为三类：(1)贪婪算法，其特点是每次迭代选取一个最能匹配残差图像的原子，如匹配追踪(Matching Pursuit, MP)算法、正交匹配追踪(Orthogonal Matching Pursuit, OMP)算法等；(2) 松弛算法，其特点是将稀疏分解对应的l₀优化问题松弛为易求解的l_p(p≥1)优化问题，如基追踪(Basis Pursuit, BP)算法；(3) 函数逼近法，其特点是利用构造的特殊函数序列来逼近l₀范数，如光滑L0范数(Smoothed L0 Norm, SL0)系列算法。上述方法的一个共同特点是计算量十分巨大，这也是图像稀疏分解实际应用发展缓慢的主要原因之一。

国外学者Gribonval和Vandergheynst对MP算法在理论上做出的收敛证明，使MP算法有了充分的理论依据，MP算法相对较低的计算复杂度也使得它成为是稀疏分解方法中最常用的一种。为降低图像稀疏分解方法的运算复杂度，研究人员对MP算法提出了众多改进方法。如全局匹配追踪(Full Search Matching Pursuit, FSMP)算法(Figueras i Ventura R M, Vandergheynst P, Frossard P. Low-rate and flexible image coding with redundant representations[J]. IEEE Transactions on Image Processing, 2006, 15(3):726-739)将快速傅里叶变换(Fast Fourier Transform，FFT)引入到MP算法中用于内积的批量计算，大为提高了稀疏分解速度。由于FFT实际上是复数运算，对于虚部为0的实信号来说，增加了不必要的计算开销，为此，基于快速哈特莱变换(Fast Hartley Transform，FHT)的实信号稀疏分解算法被提出(刘浩, 潘伟. 基于FHT的实信号稀疏分解快速算法. 西南交通大学学报, 2009, 44(1):45-48)，并被应用于二维图像(尹忠科, 王在磊, 和红杰, 王建英. 基于一维快速哈特莱变换和匹配追踪的图像稀疏分解快速方法:中国, 200910216797.8[P]. 2011-12-21)，其运算速度比MP算法提高1至2个数量级。以上算法由于每次迭代只搜索一个最佳原子，速度提升有限，且将图像转为一维信号后再作互相关运算存在较大误差。M项追踪(M-Term Pursuit，MTP)算法(Rahmoune A, Vandergheynst P, Frossard P. Sparse approximation using M-term pursuit and application in image and video coding[J]. IEEE Transactions on Image Processing, 2012, 21(4):1950-1962)将冗余字典划分成若干非相干子字典，然后从子字典中并行一次选择多个匹配原子，极大地提高了稀疏分解速度。但子字典的划分对稀疏逼近质量有较大影响。AMMP(Approximate M-fold Matching Pursuit)算法(Gan Tao, He Yan-Min, Zhu Wei-Le. Fast M-fold matching pursuit algorithm for image approximation[J]. Journal of Systems Engineering and Electronics, 2009, 20(4): 883-888)基于对原子间互相关信息的估计，并行一次选择多个原子，在微小精度损失的情况下大幅降低了稀疏分解运算复杂度。然而，每次迭代都对所有原子不加鉴别地全部作互相关比较，仍然浪费了很多运算时间。

总之，目前已有的基于匹配追踪的图像稀疏分解方法，包括改进的方法，由于计算量大、运算速度慢的原因还很难被推广产业化。因此，进一步研究图像稀疏分解快速方法不仅具有重要的学术价值，也具有广泛的潜在应用前景。

发明内容

本发明的目的是为了解决全局搜索和/或单原子选择导致的图像稀疏分解方法运算速度慢的问题，而提出一种基于局部多原子匹配追踪的图像稀疏分解快速方法。

本发明的目的是通过下述技术方案实现的。

基于局部多原子匹配追踪的图像稀疏分解快速方法，其特征在于步骤如下：

(1) 生成并存储核原子：对大小为A×B的原始图像F，采用合适的原子生成函数，令原子参数组中的中心位置参数为(0,0)、其余参数用适当方法离散化后生成K个大小为(2A－1)×(2B－1)的核原子，对所有核原子截取其能量大于其总能量99%的最小区域存储于计算机内存；

(2) 稀疏分解初始化：设定匹配原子总数T、最优性因子α、相干性阈值δ和局部搜索比例ρ；令图像残差R=F，用二维快速哈特莱变换实现R与步骤(1)存储的K个核原子的离散互相关运算，得到K个投影系数矩阵，然后求出每个投影系数矩阵中元素绝对值的最大值，得到K个核原子的匹配能力值，将匹配能力值从大到小排序，取前L个能力值对应的核原子构成下一轮迭代的局部搜索核原子集Θ，同时取前L个能力值对应的投影系数矩阵构成本轮迭代的匹配搜索系数矩阵集Φ，其中L为不大于ρK的最大整数；置多原子匹配集G_Λ为空集；

(3) 多原子匹配搜索：求出匹配搜索系数矩阵集Φ中全部元素绝对值的最大值，得到最大匹配能力值Y及其对应的普通原子g_γ1，将g_γ1作为第一个原子放入到多原子匹配集G_Λ中；从Φ中依次取出绝对值不小于αY的元素对应的普通原子，判断该普通原子与G_Λ中所有原子的相干系数是否都小于δ，若是，则将该普通原子加入到G_Λ中，若不是，则G_Λ保持不变；反复筛选直至Φ中绝对值不小于αY的全部元素对应的普通原子判断完毕，得到本轮迭代的多原子匹配集G_Λ；记录G_Λ中每个原子在Φ中对应的元素值c_i及其对应的参数组γ_i，其中i表示g_γi为G_Λ中第i个原子；

(4) 图像残差更新：将图像残差R减去多原子匹配集G_Λ中第一个原子g_γ1与c₁的乘积，得到本轮迭代第一次更新的图像残差R’，令R=R’；再将R减去G_Λ中第二个原子g_γ2与c₂的乘积，得到本轮迭代第二次更新的图像残差R’，再令R=R’；反复更新，直到G_Λ中所有原子都完成残差更新操作；

(5) 分解结束判断：如果搜索到的匹配原子个数超过T时，则结束操作，得到原始图像F的稀疏分解结果{(c_j,γ_j) | j=1,2,…,T}；否则转步骤(6)进行新一轮迭代；

(6) 局部搜索核原子集更新：用二维快速哈特莱变换计算当前图像残差R与局部搜索核原子集Θ中所有核原子的离散互相关，得到L个新的投影系数矩阵，将这些新矩阵构成新的匹配搜索系数矩阵集Φ，求出Φ中每个系数矩阵元素绝对值的最大值，得到L个新的匹配能力值，将新得到的匹配能力值与其他未更新的K－L个匹配能力值合在一起并从大到小排序，取前L个能力值对应的核原子构成新的局部搜索核原子集Θ；将多原子匹配集G_Λ清空，转步骤(3)。

上述步骤(2)中的最优性因子α的取值区间为[0.7,1]，相干性阈值δ的取值区间为[0,0.1]，局部搜索比例ρ的取值区间为[0.1,0.3]。

与现有技术相比，本发明的有益效果在于：

1．将局部搜索和多原子选择两种匹配追踪方式相结合，最大限度地降低了运算复杂度，从而极大地提高了稀疏分解速度。

现有的基于匹配追踪(MP)算法的图像稀疏分解方法都属于全局搜索方法，即使是引入快速傅里叶变换(FFT)或快速哈特莱变换(FHT)用于内积的批量计算，由于每次迭代只搜索一个最佳匹配原子，其运算速度的提升仍比较有限。已有的多原子匹配追踪算法采用一次迭代搜索多个匹配原子的逼近方式，为图像稀疏分解方法的研究提供了新的方向，但这类方法依然属于全局搜索方法，因为每次迭代都对字典中的所有原子不加鉴别地全部作互相关比较，浪费了很多不必要的运算时间。为克服以上技术的不足，本发明提供一种局部搜索方式，局部搜索是指每轮迭代只将残差图像与小部分而不是全部核原子作互相关运算，然后在互相关结果中搜索匹配原子，因为残差图像在原子空间上的投影值在相邻代变化较小，故在前一轮迭代对核原子的匹配能力作降序排序后取前一小部分核原子，可以保证这小部分核原子对应的普通原子是最好的一批；局部搜索方式为图像稀疏分解方法的研究提供了另一种新方向。为融合现有的先进技术，本发明在分析MP算法的特性和挖掘相邻代的核原子排序规律的基础上，引入二维FHT用于批量计算内积，再将局部搜索和多原子选择两种匹配追踪方式有机结合并应用于图像稀疏分解，最大限度地提高了运算速度。

2．图像残差更新采用逐原子依次更新的方式，保证了稀疏分解的逼近质量。

本发明迭代一轮选择多个原子，其图像残差更新方式不是将这多个原子与对应投影值做乘积再累加后一次更新，也不是将残差在这多个原子所张成的子空间上作正交投影后再更新，而是将这些原子逐个依次地进行图像残差更新，其过程与传统MP算法的残差更新方法一致，因此本发明方法所得稀疏分解结果图的逼近质量更高，与MP算法的结果很接近，相比于正交投影方法，本发明的更新方法节省了正交投影计算时间。

3．核原子大小取为(2A－1)×(2B－1)，能精确对应各平移位置的普通原子。

已有的稀疏分解技术对核原子大小未做说明，如果核原子大小与原始图像大小相同，会造成中心位置平移后部分区域缺失，导致对应的普通原子表达误差较大。对于大小为A×B的原始图像，本发明将核原子大小取为(2A－1)×(2B－1)，那么不论核原子中心位置平移到何处，截取相应区域后都能精确得到对应的普通原子。

附图说明

图1是用来对128×128图像进行稀疏分解的核原子与普通原子之间的转化示意图，其中(a)为255×255核原子，虚线框内区域表示同参数的128×128核原子；(b)表示255×255核原子平移后得到的中心位置为(80,120)的128×128普通原子；(c)表示128×128核原子平移后得到的中心位置为(80,120)的128×128普通原子。

图2是用本发明方法对Barbara(256×256)图像进行稀疏分解得到的由不同个数匹配原子稀疏表示的结果图，其中(a)是原始图像，(b)是由200个匹配原子稀疏表示的结果图，(c)是由800个匹配原子稀疏表示的结果图，(d)是由3000个匹配原子稀疏表示的结果图。

图3是给定峰值信噪比(PSNR)的条件下，本发明方法与其他三种先进方法：FSMP算法、MTP算法、AMMP算法的稀疏分解耗时比较。

具体实施方式

以下结合附图和实施例对本发明中的技术方案作详细说明。

基于局部多原子匹配追踪的图像稀疏分解快速方法，其实施步骤如下：

(1) 生成并存储核原子：对大小为A×B的原始图像F，原子生成函数采用如下的二维Gabor函数：

式中，γ={a,b,σ,η,θ,ω,φ}为原子参数组，(a,b)表示原子的中心位置，σ、η、θ、λ和φ分别表示原子的尺度、纵横比、旋转角度、波长和相位。用频率带宽ω替换波长λ以便原子在视觉上的理解更直观，转换关系式为：

令h={σ,η,θ,ω,φ}，对参数组h中的参数离散化：σ={2ⁿ | n取[0,log₂N－1]区间的整数, N=min(A,B)}、η={0.3,0.5,0.7,1}、θ={nπ/12 | n=0,…,11}、ω={0.5,0.7,1,2}、φ={0,π/2}；设x取[－A＋1,A－1]上的整数值、y取[－B＋1,B－1]上的整数值，并令(a,b)=(0,0)，则由二维Gabor函数可生成K个大小为(2A－1)×(2B－1)的核原子g_h；这里的核原子个数K由σ、η、θ、ω和φ的离散值组合个数决定，核原子大小(2A－1)×(2B－1)由x、y的取值范围决定。如果参数组γ中a取[0,A－1]区间的整数、b取[0,B－1]区间的整数、其他参数取值同h，设x取[0,A－1]上的整数值、y取[0,B－1]上的整数值，则可得到KAB个大小为A×B的普通原子g_γ。所有普通原子构成冗余字典D={g_γ}，相应地，所有核原子构成核字典H={g_h}。

由上可知，每个核原子g_h都对应着AB个普通原子g_γ，只要知道中心位置(a,b)，截取核原子的相应区域就能得到对应的普通原子。如图1所示，是用来对128×128图像进行稀疏分解的核原子与普通原子之间的转化示意图，其中(a)为255×255核原子，虚线框内区域表示同参数的128×128核原子；(b)表示255×255核原子平移后得到的中心位置为(80,120)的128×128普通原子；(c)表示128×128核原子平移后得到的中心位置为(80,120)的128×128普通原子。可以看出，(b)精确表达了相应的普通原子，而(c)表达的普通原子误差较大。图1解释了将核原子大小取为(2A－1)×(2B－1)的必要性。

将核原子存于计算机内存直接读取，可减省稀疏分解过程中生成原子所花的时间。但核原子太大且个数太多，在普通计算机上存储并非易事，为了即节省存储量又不影响稀疏分解的逼近质量，对所有核原子截取其能量大于其总能量99%的最小区域存储于计算机内存，故存储于内存的核原子其大小低于(2A－1)×(2B－1)，尺度σ越小，存于内存的核原子越小，从而所占内存越少。

(2) 稀疏分解初始化：设定匹配原子总数T=3000、最优性因子α=0.8、相干性阈值δ=0.01和局部搜索比例ρ=0.2；令图像残差R=F，用表示R与核原子g_h的离散互相关运算结果，包含着g_h对应的AB个普通原子g_γ与R的内积值(即投影系数)，对应关系式如下：

和都称为投影系数矩阵。用二维快速哈特莱变换实现离散互相关运算可提高计算速度。

投影系数矩阵中元素绝对值的最大值或代表核原子g_h的匹配能力，其中符号表示求矩阵的无穷大范数得到最大的元素绝对值。对核字典H中所有K个核原子的匹配能力值从大到小排序，取前L个能力值对应的核原子构成下一轮迭代的局部搜索核原子集Θ，即

其中，L为不大于ρK的最大整数，是g_h的匹配能力值在所有核原子能力值中的排位；同时取前L个能力值对应的投影系数矩阵构成本轮迭代的匹配搜索系数矩阵集Φ，即

置多原子匹配集G_Λ为空集。

(3) 多原子匹配搜索：求出匹配搜索系数矩阵集Φ中全部元素绝对值的最大值，得到最大匹配能力值Y，根据Y在矩阵中的行列位置用步骤(2)中的对应关系式得到对应的普通原子g_γ1，将g_γ1作为第一个原子放入到多原子匹配集G_Λ中。从Φ中依次取出绝对值不小于αY的元素对应的普通原子g_γ，用公式

判断g_γ与G_Λ中所有原子的相干系数是否都小于δ，若是，则将g_γ加入到G_Λ中；若不是，则G_Λ保持不变；反复筛选直至Φ中绝对值不小于αY的全部元素对应的普通原子判断完毕，得到本轮迭代的多原子匹配集G_Λ。记录G_Λ中每个原子在Φ中对应的元素值c_i及其对应的参数组γ_i，其中i表示g_γi为G_Λ中第i个原子。

(4) 图像残差更新：将图像残差R减去多原子匹配集G_Λ中第一个原子g_γ1与c₁的乘积，得到本轮迭代第一次更新的图像残差R’，令R=R’；再将R减去G_Λ中第二个原子g_γ2与c₂的乘积，得到本轮迭代第二次更新的图像残差R’，再令R=R’；反复更新，直到G_Λ中所有原子都完成残差更新操作。

(5) 分解结束判断：如果搜索到的匹配原子个数超过T时，则结束操作，得到原始图像F的稀疏分解结果{(c_j,γ_j) | j=1,2,…,T}；否则转步骤(6)进行新一轮迭代。

(6) 局部搜索核原子集更新：用二维快速哈特莱变换计算当前图像残差R与局部搜索核原子集Θ中所有核原子的离散互相关，得到L个新的投影系数矩阵，将这L个新矩阵构成新的匹配搜索系数矩阵集Φ，即

求出Φ中每个的元素绝对值的最大值，得到L个新的匹配能力值，将新得到的匹配能力值与其他未更新的K－L个匹配能力值合在一起并从大到小排序，取前L个能力值对应的核原子构成新的局部搜索核原子集Θ；将多原子匹配集G_Λ清空，转步骤(3)。

本发明中的二维快速哈特莱变换(FHT)实现互相关运算为现有技术，其基本原理如下：

若设ψ(k,l)为f(k,l)与g(k,l)的二维离散互相关函数，则由离散互相关定理可得到ψ(k,l)的二维离散哈特来变换(Discrete Hartley Transform, DHT)的计算公式为：

这里，Ψ(u,v)为ψ(k,l)的二维DHT，F_e(u,v)、F_o(u,v)分别为f(k,l)的二维DHT的偶分量和奇分量，G_e(u,v)、G_o(u,v)分别为g(k,l)的二维DHT的偶分量和奇分量。对Ψ(u,v)作二维DHT逆变换就能得到互相关结果ψ(k,l)。

用二维DHT实现互相关运算的优势在于：DHT是实数变换，运算效率高；DHT的正、逆变换具有相同的形式，软件实现程序相同；DHT的快速算法FHT能大大提高运算速度。FHT的思想是：将长度为2^M序列的DHT转换为2个2^M-1点的奇偶采样序列的DHT，而每个长度为2^M-1序列的DHT又能转换为2个2^M-2点采样序列的DHT，逐次分解下去，最终分解成长度为2序列的DHT，用M级“哈特莱蝶形”计算DHT，从而减少实数加法和乘法的计算次数，实现快速计算的目的。

本发明的优点由以下的仿真数据和图像进一步说明。

图2是用本发明方法对Barbara(256×256)图像进行稀疏分解得到的由不同个数匹配原子稀疏表示的结果图，其中(a)是原始图像，(b)是由200个匹配原子稀疏表示的结果图，(c)是由800个匹配原子稀疏表示的结果图，(d)是由3000个匹配原子稀疏表示的结果图。由图2看出，用200个原子已能稀疏表示出原始图像的大致轮廓，此时的峰值信噪比(PSNR)为24.68dB；当匹配原子数为800个时，稀疏表示的结果图已经比较清晰，PSNR为29.47dB；匹配原子数达到3000个时，稀疏表示的结果图接近原始图像，此时PSNR为35.92dB。图2说明了本发明方法作图像稀疏分解的过程，验证了本发明方法的有效性，同时也表明图像稀疏分解结果的逼近质量并不是随匹配原子个数增大而均匀增加的。

图3是给定峰值信噪比(PSNR)的条件下，本发明方法与其他三种先进方法：FSMP算法、MTP算法、AMMP算法的稀疏分解耗时比较。由图3看出，随PSNR的增加，FSMP算法的运算时间成幂级数增长，MTP算法和AMMP算法的运算耗时增长比FSMP算法慢，但AMMP的运算速度比MTP快。相比于其他算法，本发明方法的耗时增长很缓慢，随PSNR的增大其运算速度的优势越来越明显，这是由于随着迭代次数的增加，匹配搜索系数矩阵集Φ中绝对值大于αY的元素越来越多，每一轮迭代选取的非相干匹配原子数也更多，从而加快了稀疏分解的速度。

表1列出了本发明方法与FSMP算法、MTP算法、AMMP算法对四幅512×512标准测试图像：Lena、Barbara、Goldhill和Peppers在运行时间和逼近质量上的实验结果，其中的时间单元是在相同软硬件条件下，以FSMP的运行时间为基准折算所得。从表1看出，与FSMP相比，本发明方法在微小精度损失的情况下，具有显著的速度优势，且随着逼近原子数的增多该优势迅速增强；与多原子搜索算法MTP和AMMP相比，本发明方法的运行时间更少、逼近质量更高，这得益于FHT的应用和残差更新方式的不同。当搜索800个匹配原子时，本发明方法的速度提升分别约为FSMP算法、MTP算法、AMMP算法的150倍、27倍、3倍左右，而逼近质量PSNR仅比FSMP算法至多低0.14dB，比MTP算法和AMMP算法的逼近质量略高。

表1  本发明方法与FSMP、MTP、AMMP算法的运行时间和逼近质量的比较