基于一维快速哈特莱变换和匹配追踪的图像稀疏分解快速算法

摘要

本发明公布了一种基于一维快速哈特莱变换和匹配追踪的图像稀疏分解快速算法，它首先构造了核心原子库，并将原始二维图像转换成一维实信号，稀疏分解原子库中的原子转换成一维原子，然后利用一维快速哈特莱变换实现图像或图像的残差与原子的互相关运算，寻找最佳原子，最终实现图像的分解。其优点是对图像的稀疏分解速度快，且重构图像视觉效果好。

说明书

技术领域

本发明涉及一种用于静止图像去噪、图像识别和图像压缩等方面的匹配追踪图像稀疏分解快速算法。

背景技术

在数字图像处理的实际工程应用中，数字图像的表达或分解是一个关键环节。1994年，Mallat等人提出了图像稀疏分解的匹配追踪(Matching Pursuit，匹配追踪)方法(BERGEAU F，MALLAT S.Matching pursuit of images.Proceedings ofIEEE-SP.Piladel-phia，PA，USA，1994.330-333)。此后，稀疏分解在图像处理领域的应用得到不断发展。在图像稀疏分解中，选择何种类型的原子构造过完备原子库是实现时频搜索的关键。针对过完备原子库，为了更好地表示图像内容，研究者提出了非对称原子库、二维Gabor原子库等，同时新发展的Ridgelet、Curvelet、Bandelet和Contourlet等也可作为原子模型以形成过完备原子库。

目前研究者已提出了多种图像稀疏分解算法，如基于匹配追踪的图像稀疏分解、BP算法、MOF算法和BOB算法等，其中最常用的方法是基于匹配追踪的图像稀疏分解。和其他图像稀疏分解算法一样，基于匹配追踪的图像稀疏分解也存在计算量巨大、运算速度慢的问题。研究发现基于匹配追踪的图像稀疏分解的主要运算量是分解中寻找最佳原子的匹配运算，即图像或图像残差与原子的内积运算。由于这种内积运算是在高维空间的内积运算(空间维数和图像大小相同)，而且要进行很多次(每一步的内积运算次数与原子库中原子的个数相同)。因此，寻求图像稀疏分解中内积运算的快速算法是提高基于匹配追踪的图像稀疏分解速度的关键。为此，研究者提出采用遗传算法(李恒建，尹忠科，王建英.基于量子遗传算法的图像稀疏分解[J].西南交通大学学报.2007，42(1)：19-23)来提高其计算速度，然而由于遗传算法是一局部搜索的优化算法，很难在全部原子中寻找最优时频原子，从而导致重构出的图像效果差；其他类似基于智能计算的方法都存在类似的问题(李恒建，尹忠科，张家树，王建英.基于混沌变异粒子群优化算法的图像稀疏分解[J].西南交通大学学报.2008，43(4)509-513)；此外，有学者提出在过完备原子库上按原子结构的分类特性生成原子的快速算法(华泽玺，尹忠科，黄雄华.信号在过完备库上分解中原子形成的快速算法[J].西南交通大学学报，2005，40(3)：402-405.)，该方法一定程度上提高了运算速度，但是仍然需要大量的高维内积运算.；为实现在全局中找最佳原子，有学者提出利用快速傅里叶变换(Fast Fourier Transform，FFT)来实现互相关运算来寻找最佳原子，但是快速傅里叶变换涉及到复数运算，需将实数信号人为添加虚数变量，增加了计算的复杂度，计算速度提高有限；在快速傅里叶变换算法基础上，针对一维实信号，刘浩等提出利用快速哈特莱变换(FastHartley Transform，FHT)实现匹配追踪的快速算法(刘浩，潘伟.基于FHT的实信号稀疏分解快速算法[J].西南交通大学学报.2009，44(1))，这种算法在计算速度方面优于基于快速傅里叶变换的匹配追踪快速算法，又是一种全局搜索算法，避免了遗传算法只能做到局部最优搜索的局限性，然而它只能实现一维实信号分解。

总之，采用匹配追踪方法进行图像的稀疏分解在静止图像的去噪、图像识别和图像压缩多个方面取得了成功，但目前还很难被推广而产业化，其主要原因是图像稀疏分解的匹配追踪计算量巨大，计算时间在现有条件下令人无法忍受。因此，图像稀疏分解快算法的研究已成为国内外近年来的前沿性研究课题，不仅具有重要的学术价值，也具有广泛的潜在应用前景。

发明内容

本发明的目的是提供一种基于一维快速哈特莱变换和匹配追踪的图像稀疏分解快速算法，该算法对图像的稀疏分解速度快，且重构图像视觉效果好。

本发明解决其技术问题，所采用的技术方案为：基于一维快速哈特莱变换和匹配追踪的图像稀疏分解快速算法，包括如下步骤：

(1)核心原子库的形成：

对大小为M×N的原始图像f，构造稀疏分解用的原子库，原子g_γ由参数组γ＝(θ，u，v，s_x，s_y)定义，其中θ代表原子的旋转角度，u、v分别为原子在x、y方向上的平移，s_x、s_y分别为原子在x、y方向上的伸缩；令u＝M/2和v＝N/2形成二维的核心原子库{g_γ′}，其中参数组γ′按以下方法离散化：

${\gamma }^{\prime }=(\theta ,\frac{M}{2},\frac{N}{2},s_x,s_y)=(\frac{\pi {\theta }^{\prime }}{\max (M,N)},\frac{M}{2},\frac{N}{2},2^{\frac{\mathrm{sty}}{5}},2^{\frac{\mathrm{sty}}{5}})$

其中1≤θ′≤max(M，N)，0≤stx≤5log^M-5，$0\le \mathrm{sty}\le 5{\log }_2^N-5,$θ′、stx、sty均在相应区间上取整数；

(2)原子与图像转化：

将(1)步形成的核心原子库中的所有原子g_γ′转换成一维原子g′_γ′；将大小为M×N的原始图像f转换成一维信号f′；

(3)图像的初始分解：

对(2)步得到的图像的一维信号f′与所有一维原子g′_γ′采用一维快速哈特莱变换算法做互相关运算，得到一维信号f′与过完备原子库{g_γ}中所有一维原子g′_γ的匹配系数p_γ，并记录匹配系数p_γ绝对值最大的一维原子对应的参数组${\gamma }_0=({\theta }_0,u_0,v_0,s_{x_0},s_{y_0})$和对应的匹配系数将图像的一维信号f′减去最匹配原子与匹配系数的乘积，得到图像残差R_f′；

(4)图像分解：

对当前得到的图像残差R_f′与所有一维原子g′_γ′采用一维快速哈特莱变换算法做互相关运算，得到图像残差R_f′与过完备原子库{g_γ}中每个原子的匹配系数p_γ，并记录匹配系数p_γ绝对值最大的原子对应的参数组${\gamma }_i=({\theta }_i,u_i,v_i,s_{x_i},s_{y_i})$和对应的匹配系数其中i为图像分解的当前次数；将图像残差R_f′减去最匹配原子与匹配系数的乘积，得到更新后的图像残差R_f′；

重复以上的图像分解操作，当图像分解结果满足后续处理要求时，结束操作，得到原始图像f匹配追踪稀疏分解的结果$\{(P_{{\gamma }_i},{\gamma }_i)=(P_{{\gamma }_i},{\theta }_i,u_i,v_i,s_{x_i},s_{y_i})|i=\mathrm{0,1,2},...,n\},$其中n为图像分解操作的总次数。

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

一、本发明在每步为图像或图像残差寻找最佳原子的匹配运算中，将其与一组原子的内积运算转化为与一个原子的一次互相关运算，在降低运算次数的同时，还实现了在整个过完备原子库中的全部原子进行逐一比较，然后选取出真正的最佳原子，是一种全局最优算法。因此，本发明的重构图像视觉效果好，而且有效克服了基于遗传算法等智能计算方法的图像稀疏分解由于局部寻优，即在随机选取的一部分原子中寻找最佳原子(近似最优原子)而带来的结果随机性问题。同时，本发明采用一维快速哈特莱变换算法实现图像或图像残差与原子的互相关运算，不仅提高了图像稀疏分解的速度，而且降低了对数据存储容量的要求。离散哈特莱变换(Discrete Hartley Transform，DHT)是一种类似于离散傅里叶变换(Discrete Fourier Transform，DFT)的积分变换，具有DFT变换的大部分特性，但其正反变换的积分核相同。实函数的DHT变换仍是实函数，从而避免了变换过程中的复数运算，因此数据存储容量仅为FFT的一半，而且快速哈特莱变换可以采用FFT的结构形式，能够保障很高的运算速度，其运算量大约是FFT的一半。

二、本发明固定原子位置参u＝M/2和v＝N/2形成核心原子库，核心原子库中的原子个数远小于过完备原子库，降低了稀疏分解对存储量的要求，达到了计算量和存储量之间一个很好的平衡。同时，还能保持稀疏分解效果和时间不受影响。这是因为核心原子库包含了过完备原子库中所有不同形状的原子，取出核心原子库中的一个原子，在进行互相关运算时，通过平移(即参数u和v的连续变化取值)即可得到过完备原子库中具有相同形状，但其中心位置不同的一组原子，在每一组原子中寻找出匹配系数绝对值最大的原子，然后比较各组的匹配系数绝对值最大原子的匹配系数的绝对值，再从中找出匹配系数绝对值最大的原子，即在过完备原子库中的所有原子中寻找出最佳匹配原子，因此不会降低图像分解的质量；过完备原子库中的非核心原子库的原子仅是在互相关运算中临时通过平移产生，比较后的非最佳匹配原子不用储存与记录，其存储要求的容量小；而且，由于原子的平移几乎不需要计算量，所以整个稀疏分解过程的速度几乎不会变化。

三、本发明将原始图像或图像残差按一行一行地抽取首尾相接成转换成一维信号，同时核心原子库中的原子也按一行一行地抽取首尾相接转换成一维原子，既为后面用一维快速哈特莱变换进行分解和互相关运算提供了可能，又节省数据存储空间，同时还不会影响相邻像素点的相关性，图像匹配运算中效果也不会受到影响。

总之，本发明方法作为一种全局搜索最优算法，结合过完备库的结构特性、匹配追踪技术和一维快速哈特莱变换技术，把二维图像转换成一维实信号，在提高图像稀疏分解质量、重构图像视觉效果好的情况下，还提高了稀疏分解的速度。即，实现了在解的质量方面和计算的复杂性方面取得一个较好的平衡。

下面结合附图和具体实施方式对本发明作进一步详细说明。

附图说明

图1为本发明方法对一实际图像进行分解时使用的原始图像及最终分解后的原子图、图像残差和重构图像。其中a分图为从标准Lena图像(128×128)中抽取的一小块大小为32×32的待处理原始图像；b分图为第30次图像分解时得到的最佳匹配原子的示意图；c分图是第30次图像分解后的残差图像；d分图是第30次图像分解后重建的结果，即重构的图像。

图2为基于遗传算法的匹配追踪图像稀疏分解算法对图1的同一实际图像进行分解时使用的原始图像及最终分解后的原子图、图像残差和重构图像。其中a分图与图1的a分图为同一图像；b分图为第30次图像分解时得到的最佳匹配原子的示意图；c分图是第30次图像分解后的残差图像；d分图是第30次图像分解后重建的结果，即重构的图像。

具体实施方式

实施例

基于一维快速哈特莱变换和匹配追踪的图像稀疏分解快速算法，包括如下步骤：

(1)核心原子库的形成：

对大小为M×N的原始图像f，构造稀疏分解用的原子库，非对称原子的基本形式表示如下：

$g(x,y)=(4x^2-2)e^{-(x^2+y^2)}$

通过对基本非对称原子进行旋转、平移和伸缩变换，可以得到一系列原子g_γ，从而形成原子库D＝{g_γ}_γ∈Γ。

$g_{\gamma }=g_{\theta }(\frac{x-u}{s_x},\frac{y-v}{s_y})$

从而，原子g_γ可以由参数组γ＝(θ，u，v，s_x，s_y)定义，其中θ代表原子的旋转角度，u、v分别为原子在x、y方向上的平移，s_x、s_y分别为原子在x、y方向上的伸缩；令u＝M/2和v＝N/2形成核心原子库{g_γ′}，其中参数组γ′按以下方法离散化：

${\gamma }^{\prime }=(\theta ,\frac{M}{2},\frac{N}{2},s_x,s_y)=(\frac{\pi {\theta }^{\prime }}{\max (M,N)},\frac{M}{2},\frac{N}{2},2^{\frac{\mathrm{sty}}{5}},2^{\frac{\mathrm{sty}}{5}})$

其中1≤θ′≤max(M，N)，0≤stx≤5log^M-5，$0\le \mathrm{sty}\le 5{\log }_2^N-5,$θ′、stx、sty均在相应区间上取整数。

在过完备原子库中，参数θ、s_x和s_y定义了一个原子的形状，而参数u和v定义了一个原子的中心位置。如果原子库中的不同原子具有相同的参数θ、s_x、s_y和不同的参数u、v，这样的原子各方面的特征都相同，只有中心点位置不同而已。从理论上讲，这种不同的原子包含在原子库中是完全必要的，但从实际计算的角度讲，原子库中包含众多这样的原子是根本没有必要的。它违背了原子库良好结构定义中的含义(一是原子库中应包含尽可能多的原子个数和种类，以达到稀疏分解的目的，获得稀疏分解的良好效果；二是原子库中应尽可能不包含相近似的原子，以满足存储量和计算量方面的要求)。因此，一个具有良好结构的原子库中，参数u和v应是不变的，且u＝M/2，v＝N/2。这样原子库的大小将大大减小，适合一般计算机内存的要求。在稀疏分解过程中，达到了计算量和存储量之间一个很好的平衡。

(2)原子与图像转化：

将(1)步形成的核心原子库中的所有原子gγ′转换成一维原子g′_γ′；将大小为M×N的原始图像f转换成一维信号f′；

将二维图像数组按一行一行地抽取首尾相接成一维的数组，这样既方便后面用一维快速哈特莱变换进行分解和互相关运算，又节省数据存储空间。同时还不会影响相邻像素点的相关性，在做图像匹配运算中效果会好一些。

(3)图像的初始分解：

对(2)步得到的图像的一维信号f′与所有一维原子g′_γ′采用一维快速哈特莱变换算法做互相关运算，得到一维信号f′与过完备原子库{g_γ}中所有一维原子g′_γ的匹配系数p_γ，并记录匹配系数p_γ绝对值最大的一维原子对应的参数组${\gamma }_0=({\theta }_0,u_0,v_0,s_{x_0},s_{y_0})$和对应的匹配系数将图像的一维信号f′减去最匹配原子与匹配系数的乘积，得到图像残差R_f′；

(4)图像分解：

对当前步得到的图像残差R_f′与所有一维原子g′_γ′采用一维快速哈特莱变换算法做互相关运算，得到图像残差R_f′与过完备原子库{g_γ}中每个原子的匹配系数p_γ，并记录匹配系数p_γ绝对值最大的原子对应的参数组${\gamma }_i=({\theta }_i,u_i,v_i,s_{x_i},s_{y_i})$和对应的匹配系数其中i为图像分解的当前次数；将图像残差R_f′减去最匹配原子与匹配系数的乘积，得到更新后的图像残差R_f′；

重复以上的图像分解操作，当图像分解结果满足后续处理要求时，结束操作，得到原始图像f匹配追踪稀疏分解的结果$\{(P_{{\gamma }_i},{\gamma }_i)=(P_{{\gamma }_i},{\theta }_i,u_i,v_i,s_{x_i},s_{y_i})|i=\mathrm{0,1,2},...,n\},$其中n为图像分解操作的总次数。

在图像稀疏分解中，寻找最佳原子的过程，就是计算图像f或图像的残余R_f′在对应原子g_γ′上的分量，即计算内积<R_f′，g_γ′>，对具有参数θ、s_x、s_y的一个原子g_γ′，u取值[0，M-1]，v取值[0，N-1]，则该原子要和图像或图像残差值作M×N次内积<R_f′，g_γ′>计算，所以M×N次内积<R_f′，g_γ′>计算可以转换成一次R_f′和g_γ′的互相关运算。虽然从理论上讲，这样的转换，计算量几乎没有变化，但是计算效率却大大提高。由于利用快速哈特莱变换可以快速实现互相关运算，所以我们利用快速哈特莱变换对以上算法作进一步的改进，这样的改进，丝毫不会影响图像稀疏分解的效果，但却可以大大提高稀疏分解的速度。因为一般而言，用快速哈特莱变换算法实现互相关运算，可以使计算速度提高2～3个数量级。

本发明的一维快速哈特莱变换实现互相关运算为现有算法，其原理如下：长度为L的有限长实序列x(n)的离散Hartley(DHT)变换定义为：

$X^h\left(k\right)=\mathrm{DHT}\left[x\left(n\right)\right]=\underset{n=0}{\overset{L-1}{\Sigma }}x\left(n\right)\mathrm{cas}\frac{2\mathrm{\pi nk}}{N},$其中k＝0，1，...，L-1

对信号奇偶抽取，并考虑周期性得到DHT的快速运算方法，快速哈特莱变换：

$X^h\left(k\right)=X_0^h\left(k\right)+\cos \left(\frac{2\pi }{L}k\right)X_1^h\left(k\right)+\sin \left(\frac{2\pi }{L}k\right)X_1^h(\frac{L}{2}-k)$

可见，L＝2^K(K＝1，2，…)长度的序列的DHT可以转换为2个L/2点的奇偶采样序列的DHT，而每个L/2长度序列的DHT又能转换为2个L/4长度序列的DHT，逐次分解下去，最终可分解到长度为2的序列的DHT。就像一维快速傅里叶变换能大大提高离散傅里叶变换速度一样，一维快速哈特莱变换大大地提高了离散哈特莱变换的速度。在快速哈特莱变换的每个蝶形单元，需要3次实数加法运算和2次实数乘法运算，而一维快速傅里叶的每个蝶形单元，由于有复数运算，共需6次实数加法运算和4次实数乘法运算，可见一维快速哈特莱变换比快速傅里叶变换运算速度更快。

设r(n)为x(n)和y(n)的循环互相关，令R^h(k)为r(n)的DHT，则有下述DHT循环互相关定理：

$R^h\left(k\right)=X^h\left(k\right)Y_e^h\left(k\right)-Y_o^h\left(k\right)X^h(N-k)$

$r\left(n\right)=\mathrm{IDHT}\left[R^h\left(k\right)\right]=\frac{1}{N}\mathrm{DHT}\left[R^h\left(k\right)\right],$其中n＝0，1，...，L-1

通过以上DHT循环互相关定理，可以快速实现互相关运算，从而把图像或图像的残差与一类原子的内积运算<R_f′，g_γ′>转化为图像或图像残差与一个原子的一次互相关运算(即一次R_f′和g_γ′的互相关运算)，用一维快速哈特莱变换技术来实现这个互相关运算。

仿真实验

本发明的效果可以通过以下仿真实验结果进行验证：

在相同的环境(计算机软，硬件配置)下进行本发明方法与遗传算法的仿真实验。鉴于遗传算法具有一定的随机性，为了实验数据具有参考价值，以下所有的遗传算法的数据，均是20次试验的平均值。

图1为本发明方法对一实际图像进行仿真实验时使用的原始图像及最终分解后的原子图、图像残差和重构图像。其中a分图为从标准Lena图像(128×128)中抽取的一小块大小为32×32的待处理原始图像；b分图为第30次图像分解时得到的最佳匹配原子的示意图；c分图是第30次图像分解后的残差图像；d分图是第30次图像分解后重建的结果，即重构的图像。

图2为基于遗传算法的匹配追踪图像稀疏分解算法对图1的同一实际图像进行仿真实验时使用的原始图像及最终分解后的原子图、图像残差和重构图像。其中a分图与图1的a分图为同一图像；b分图为第30次图像分解时得到的最佳匹配原子的示意图；c分图是第30次图像分解后的残差图像；d分图是第30次图像分解后重建的结果，即重构的图像。

对比图1的c分图和图2的c分图可以看出，本发明算法最后得到的图像残差视觉效果不好，相比于基于遗传算法的匹配追踪图像稀疏分解算法要模糊许多，即图像残差的能量做到了最小化。说明本发明算法在相同的图像分解次数下，更好的逼近或近似的表示了原始图像。

对比图1的d分图和图2的d分图可以看出，在相同的图像分解次数下，基于遗传算法的匹配追踪图像稀疏分解算法的重构的图像视觉效果不是太好，有一些模糊不清。本发明方法的重构图像的视觉效果优于基于遗传算法的匹配追踪图像稀疏分解算法的重构图像。

以下的表1为不同图像分解次数下本发明算法与遗传算法的匹配追踪算法重构图像的PSNR(遗传算法试验的代数和个体数分别选为40、21)

表1

从表1可以看出：在同样的图像分解次数情况下，本发明方法重构图像的PSNR均比遗传算法的高。在图像分解次数为10次时，本发明方法重构图像的PSNR值要比遗传算法的高2.2161dB，当图像分解次数为90次时，本发明方法重构图像的PSNR值要比遗传算法的高5.7381dB。并且随着图像分解次数的增加，本发明方法较遗传算法的PSNR提高越明显。

以下的表2为不同PSNR情况下图像分解的完成时间。

表2

本发明方法随着图像分解次数增加PSNR值也增加，表2中两种算法的各5个PSNR值依次对应的图像分解次数为10，30，50，70，90；而遗传算法在与本发明方法的图像分解次数相同的情况下，通过增加代数和个体数实现上表相同的PSNR值。

从表2可看出，本发明方法较遗传算法的速度更快，当PSNR值在20.0314时，本发明方法比遗传算法在图像分解速度方面提高15％左右，当PSNR值在32.0138时，本发明方法比遗传算法在图像分解速度方面提高1倍以上。并且随着PSNR值的提高本发明方法的速度提高越明显，

以下的表3为不同分解速度下重构图像的PSNR。

表3

本发明方法随着图像分解次数的增加分解完成时间也相应的增加，表2中两种算法的各5个分解完成时间值依次对应的图像分解次数为10，30，50，70，90；而遗传算法在做到图像分解次数与本发明方法相同的情况下，试验中通过调整代数和个体数也实现上表相同的分解完成时间值。

从表3看出，在图像分解完成时间相同情况下，本发明方法重构图像的PSNR值较遗传算法更高，当分解时间为54.25325秒时，本发明方法较遗传算法在重构图像的PSNR方面高0.3019dB，当分解时间为561.0184秒时，本发明方法较遗传算法在重构图像的PSNR方面高1.0991dB。并且随着图像分解时间的增加，本发明方法的重构图像的PSNR值提高更明显。

本发明方法其图像分解的重复操作次数n，取决于后续处理对其重构图像的峰值信噪比(PSNR)的要求，当PSNR要求高时，分解的重复操作次数n多，反之则少。对于M×N的原始图像，若后续处理为图像压缩，其重复操作次数n一般小于2M或2N；若后续处理为图像去噪，其重复操作次数n一般大于M或N，但一般小于2M或2N。