CN102324043B - 基于dct的特征描述算子及优化空间量化的图像匹配方法 - Google Patents

Abstract

本发明涉及一种基于DCT的特征描述算子及优化空间量化的图像匹配方法，包括以下步骤：(1)对给定的N×N大小的图像块进行DCT空间的线性分解，得到基图像块的线性组合；(2)对图像块进行仿射扭曲并对姿态空间进行优化，得到量化姿态数N；(3)计算模板图像N个量化姿态DCT基图像块在不同投影系数下的组合，得到均值图像块

(4)将模板图像和检测到的特征点周围的局部图像块进行图像匹配。本发明设计合理，减少了离线进行计算的采样次数，节约了特征描绘算子的建立时间，减小了图像匹配中的计算量和数据存储量，降低了量化姿态级数，进而了提高图像匹配的速度和精度，可在OWD算法的基础上大幅提高目标匹配的实时性能。

Description

基于DCT的特征描述算子及优化空间量化的图像匹配方法

技术领域

本发明属于图像处理领域，尤其是一种基于DCT的特征描述算子及优化空间量化的图像匹配方法。

背景技术

图像特征匹配技术是计算机视觉、模式识别领域的基本问题。图像间共同部分匹配的基本目标就是找到在不同视角姿态下，不同投影参数图像中共同拥有的局部区域，或者三维物体上相似点之间的对应关系。根据采集图像的摄像机间的视点位置关系，在特征匹配过程中采集图像的设备之间的视点变化幅度较大，因此，图像场景内容、尺度、亮度等均发生了明显的变化，此时只能利用针对这些变化具有不变性的局部特征来完成图像的匹配任务。

图像特征匹配过程中，最关键的环节是提取图像块中具有不变性特征点描绘算子。国外一些专家学者提出了一些描述算子的建立方法。在David Lowe提出的尺度不变特征变换(Scale Invariant Feature Transform，SIFT)描绘算子中，通过在关键点为中心的邻域窗口内采样，并用直方图统计邻域像素的梯度方向，直方图的峰值代表该关键点处邻域梯度的主方向，即作为该关键点的方向，然后根据主方向进行调整来消除仿射变换中的旋转造成的差异，然而对于摄像机视角变化较大时，SFIT描绘算子的匹配性能下降很大，特征点的建立也需要较大的计算量。Matas通过提取图像中具有不变特性的最稳定极值区域(Maximally Stable Extremal Regions，MSER)来进行图像匹配，其原理是通过选取适当的阈值得到连通分量，并对这些连通分量的平稳性进行检测获得最终的平稳区域，MESR描述算子对大尺度变换的适应性不强，由于图像像素值的邻接性不再保持，得到的区域集合也和模板有所不同，其区分度不强，且计算复杂度较高。Stefan提出了一种基于统计学习建立特征描绘算子的单向描述符(One Way Descriptor，OWD)，此方法对特征点周围图像进行所有姿态下的扭曲，得到模板图像块的均值图像，该方法能够实时估计当前图像和模板图像之间的投影变换姿态参数，使得目标匹配可以在线实时进行图像块的添加，且较其他特征描绘算子更加快速和稳健。

在OWD特征描绘算子建立过程中，模板图像块的均值图像计算中需要利用多种采样姿态对图像块进行投影变换，由于这一操作的计算量巨大，很难满足实时匹配的需求，因此需要通过少量离线计算来提高在线的计算效率。OWD方法中采用了基于主成分分析(Principal Component Analysis，PCA)的方法，得到一组标准正交基来表示图像中的特征点，这种方法可以较好的将模板图像进行线性分解，但是需要的运算量较大，且由于PCA没有固定的变换矩阵，为了得到图像的PCA空间投影需存储大量的PCA基。

发明内容

本发明的目的在于克服现有技术的不足，提供一种基于DCT的特征描述算子及优化空间量化的图像匹配方法，该方法使用基于离散余弦变换(DCT)图像线性分解方法和姿态空间量化方法，减少了离线进行计算的采样次数，降低了图像匹配中的计算量和数据存储量，进而了提高图像匹配的速度。

本发明解决其技术问题是采取以下技术方案实现的：

一种基于DCT的特征描述算子及优化空间量化的图像匹配方法，包括以下步骤：

(1)对给定的N×N大小的图像块进行DCT空间的线性分解，得到基图像块的线性组合；

(2)对图像块进行仿射扭曲并对姿态空间进行优化，得到量化姿态数N；

(3)计算模板图像N个量化姿态DCT基图像块在不同投影系数下的组合，得到均值图像块

(4)将模板图像和检测到的特征点周围的局部图像块进行图像匹配。

而且，所述的步骤(1)线性分解后的图像块形式为：

$I=\underset{n=0}{\overset{N\×N-1}{\mathrm{\Σ}}}{C}_n{f}_n$

式中，C_n为DCT变换系数，f_n为DCT变换后的频域函数。

而且，所述步骤(2)对姿态空间进行优化的方法包括以下步骤：

(1)在仿射摄像机模型下得到镜头拉升参数λ、摄像机的轴向旋转角ψ、摄像机距离目标物体等距离球面上的纬度角θ、经度角Φ；

(2)针对不同的纬度θ采用不同的量化精度，在低纬度采用小的经度步长，高纬度采用较大的经度步长；

(3)仿射量化点采用图像平面正方向的四分之一球面分布覆盖姿态变换，总量化级数为80，得到量化姿态数N。

而且，所述均值图像块

的计算方法为：

$\stackrel{\‾}{p_{i,h}}\∝\frac{1}{M}\underset{j=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=0}{\overset{N\×N-1}{\mathrm{\Σ}}}w(C_n{f}_n,H_{h,j})=\underset{n=0}{\overset{N\×N-1}{\mathrm{\Σ}}}{C}_n\stackrel{\‾}{f_n}$

式中：C_n为DCT变换系数，

为基图像块，N为量化姿态数。

而且，所述的基图像块

采用如下公式进行计算：

$\stackrel{\‾}{f_n}=\frac{1}{M}\underset{j=1}{\overset{M}{\mathrm{\Σ}}}w(f_n,H_{h,j})$

该基图像块

采用线下预先计算并存储。

而且，所述步骤(4)的图像匹配包括以下步骤：

(1)将模板图像和检测到的特征点周围的局部图像块在初始姿态h下进行配准，使新图像块和进行姿态参数h扭曲后的图像之间差值的平方和最小，然后计算图像块和模板各个姿态的均值图像块的最小极值，并得到最接近的姿态h；

(2)添加增量姿态参数Δh，将姿态参数h进行更新，判断其二阶范数是否低于一个设定的阈值ε，如果低于ε，则认为h已经收敛，即为预期的姿态。

本发明的优点和积极效果是：

本发明首先对模板图像块进行基于DCT的线性分解，然后计算图像块在仿射摄像机模型优化姿态空间中的均值图像块，将新加入的图像和均值图像块进行匹配得到粗略姿态，最后进行迭代得到目标图像的精确姿态。该方法减少了离线进行计算的姿态空间采样次数，节约了特征描绘算子的建立时间，减小了图像匹配中的计算量和数据存储量，降低了量化姿态级数，进而了提高图像匹配的速度和精度，可在OWD算法的基础上大幅提高目标匹配的实时性能。

附图说明

图1是图像进行DCT线性分解的基图像块示意图；

图2是局部图像经过扭曲函数仿射变换后得到的不同姿态的图像示意图；

图3是仿射摄像机参数化模型示意图；

图4是图像仿射姿态参数量化示意图。

图5是图像匹配速度比较示意图；

图6是图像匹配精度比较示意图；

图7是图像匹配过程存储空间比较示意图。

具体实施方式

以下结合附图对本发明实施例做进一步详述：

一种基于DCT的特征描述算子及优化空间量化的图像匹配方法，从两个方面提高提高图像处理中的目标匹配的速度。一方面，在目标图像的均值图像块计算中，提出了基于DCT(Discrete Cosine Transform，离散余弦变换)空间的线性分解方法，用来提高特征描述算子的建立速度，降低图像仿射扭曲计算均值图像块过程中的计算量；另一方面在在仿射姿态空间量化中，引入了仿射摄像机模型对描述算子的姿态空间进行最佳的参数化和量化，以进一步降低量化姿态级数。下面对本发明进行详细说明，本图像匹配方法包括以下步骤：

步骤1：对给定的N×N大小的图像块进行DCT空间的线性分解，得到基图像块的线性组合；

在本步骤中，通过DCT变换对图像进行分解，将图像表示为基图像的线性组合。对于一个给定的N×N大小的图像块I(x，y)，其二维离散余弦变换及逆变换可以由如下定义：

$C(u,v)=\underset{x=0}{\overset{N-1}{\mathrm{\Σ}}}\underset{y=0}{\overset{N-1}{\mathrm{\Σ}}}I(x,y)\mathrm{\α}\left(u\right)\mathrm{\α}\left(v\right)\cos \left[\frac{\mathrm{\π}(2x+1)u}{2N}\right]\cos \left[\frac{\mathrm{\π}(2y+1)v}{2N}\right]$

$=\underset{x=0}{\overset{N-1}{\mathrm{\Σ}}}\underset{y=0}{\overset{N-1}{\mathrm{\Σ}}}I(x,y)f(x,y,u,v)---\left(1\right)$

$I(u,v)=\underset{u=0}{\overset{N-1}{\mathrm{\Σ}}}\underset{v=0}{\overset{N-1}{\mathrm{\Σ}}}C(u,v)\mathrm{\α}\left(u\right)\mathrm{\α}\left(v\right)\cos \left[\frac{\mathrm{\π}(2x+1)u}{2N}\right]\cos \left[\frac{\mathrm{\π}(2y+1)v}{2N}\right]$

$=\underset{u=0}{\overset{N-1}{\mathrm{\Σ}}}\underset{v=0}{\overset{N-1}{\mathrm{\Σ}}}C(u,v)f(x,y,u,v)---\left(2\right)$

其中f(x，y，u，v)是和原始图像等大小的DCT基图像，在固定图像中坐标点(x，y)时，其对应基图像块坐标由(u，v)来确定。其中u，v＝0，1，2，…，N-1，当u＝0时 $\mathrm{\α}(*)=\sqrt{1/N},$ 否则 $\mathrm{\α}(*)=\sqrt{2/N}:$

$f(x,y,u,v)=\mathrm{\α}\left(u\right)\mathrm{\α}\left(v\right)\cos \left[\frac{\mathrm{\π}(2x+1)u}{2N}\right]\cos \left[\frac{\mathrm{\π}(2y+1)v}{2N}\right]---\left(3\right)$

分解后的DCT基图像块如图1所示。

在对图像进行分解，需要计算一维DCT系数。具体方法为：

对DCT系数采用自左上角到右下角的之字形扫描，将其排列为一个一维序列，

因此可以将图像分解为如下形式：

$I=\underset{n=0}{\overset{N\×N-1}{\mathrm{\Σ}}}{C}_n{f}_n---\left(4\right)$

其中，C_n为DCT变换系数，f_n为DCT变换后的频域函数。由上式可见，任一个图像块均可以在DCT空间进行投影分解后表示成基图像的线性组合。

步骤2：对图像块进行仿射扭曲并对姿态空间进行优化，得到量化姿态数N；

在本步骤中，假设扭曲变换可以用原始图像和扭曲后图像的简单线性置换来表示，对扭曲函数w(x；H)进行线性近似，扭曲后的图像如图2所示。模板图像由一组均值图像块

构成，其中N为量化姿态数：

$\stackrel{\‾}{p_{i,h}}=\frac{1}{M}\underset{j=1}{\overset{M}{\mathrm{\Σ}}}w(p_i,H_{h,j})---\left(5\right)$

将图像块变换为基图像的线性组合后，则可以得到：

$\stackrel{\‾}{p_{i,h}}=\frac{1}{M}\underset{j=1}{\overset{M}{\mathrm{\Σ}}}w(\underset{n=0}{\overset{N\×N-1}{\mathrm{\Σ}}}{C}_n{f}_n,H_{h,j})---\left(6\right)$

在扭曲函数线性的假设下可以将上式写为：

$\stackrel{\‾}{p_{i,h}}\∝\frac{1}{M}\underset{j=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=0}{\overset{N\×N-1}{\mathrm{\Σ}}}w(C_n{f}_n,H_{h,j})$

$=\frac{1}{M}\underset{n=0}{\overset{N\×N-1}{\mathrm{\Σ}}}\underset{j=1}{\overset{M}{\mathrm{\Σ}}}w(C_n{f}_n,H_{h,j})---\left(7\right)$

$=\underset{n=0}{\overset{N\×N-1}{\mathrm{\Σ}}}{C}_n\frac{1}{M}\underset{j=1}{\overset{M}{\mathrm{\Σ}}}w(f_n,H_{h,j})$

$=\underset{n=0}{\overset{N\×N-1}{\mathrm{\Σ}}}{C}_n\stackrel{\‾}{f_n}$

上式中，

M是对量化姿态的随机采样个数，可以在线下预先对

进行计算并存储，其均值图像块

为

的线性组合，与量化姿态的采样个数M无关。

本步骤采用仿射摄像机模型对图像块进行量化姿态，如图3所示，由于图片中的特征点多来自相对平坦的目标物体，因此近似仿射摄像机获取的图像和真实的透视图像间的差异就会足够小。为减小姿态空间的量化级数，本方法提出了仿射摄像机参数化模型，来覆盖所有的姿态，具体包括如下步骤：

1、根据仿射变换的分解形式，对图像块A可以分解为：

其中：λ表示摄像机的镜头拉升参数，ψ表示摄像机的在原位以摄像机光轴为轴向的旋转角，θ为摄像机在距离目标物体等距离球面上的纬度角，Φ为摄像机在距离目标物体等距离球面上的经度角。其中λ＞0，令t＝1/cosθ则λt是矩阵A的行列式，R_i是旋转矩阵Φ∈[0，π)，T_t是对角矩阵。

2、针对不同的纬度θ对经度采用不同的量化经度，在低纬度采用小的经度步长，高纬度采用较大的经度步长。量化的纬度角范围θ∈[0，π/2)，经度角范围为Φ∈[0，π)，所有点均分布在图像像平面正前方的四分之一球面。

3、由图4量化示意图所示，做进一步优化，实验表明，当量化级数为80时，即可覆盖所有的姿态，从而得到量化姿态数N。

步骤3：计算模板图像N个量化姿态DCT基图像块在不同投影系数下的组合，得到均值图像块

所述均值图像块

的计算方法为：

$\stackrel{\‾}{p_{i,h}}\∝\frac{1}{M}\underset{j=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=0}{\overset{N\×N-1}{\mathrm{\Σ}}}w(C_n{f}_n,H_{h,j})=\underset{n=0}{\overset{N\×N-1}{\mathrm{\Σ}}}{C}_n\stackrel{\‾}{f_n}$

式中：C_n为DCT变换系数，

为基图像块，N为量化姿态数。

通过上述三个步骤即可完成单向描述符的特征描绘算子的建立过程。

步骤4：对图像进行匹配。

在本步骤中包括图像的粗略匹配和精确匹配两个步骤：

1、图像的粗略匹配步骤为：将模板图像和检测到的特征点周围的局部图像块在初始姿态h下进行配准，使新图像块和进行姿态参数h扭曲后的图像之间差值的平方和最小，然后计算图像块和模板各个姿态的均值图像块的最小极值，并得到最接近的姿态h；

我们可以非常简单的方法来实现特征的的匹配：

$h=\underset{i=i\hat{d},h}{\arg \min }{\left|\right|n\left(p\right)-n\left(\stackrel{\‾}{p_{i,h}}\right)\left|\right|}^2---\left(9\right)$

其中n(·)是归一化函数，将图像块p减去均值并除以图像块的标准差，进行归一化能够使描绘算子能够对光照变化具有较好的鲁棒性。在匹配时，我们并不对整个图像的所有像素值进行式(9)的计算，而是在图像中划分网格，并取每个网格内图像像素灰度值的局部极大值进行计算，从而大大的减少了计算量。

2、图像的精确匹配步骤为：添加增量姿态参数Δh，将姿态参数h进行更新，判断其二阶范数是否低于一个设定的阈值ε，如果低于ε，则认为h已经收敛，即为预期的姿态。

本步骤的具体过程为：在得到了粗略的姿态估计之后，我们通过一个仅需少量运算的最优化过程得到一个最佳的精确姿态参数。将我们的问题归纳到图像配准的范围内，我们的目标是将一个模板图像T(x)配准对齐到一幅输入图像I(x)，其中x＝(x，y)^T是包含图像像素坐标的一维列向量。也即对模板图像和新检测到的一个特征点周围的局部图像块在初始姿态参数为h下进行配准，以使得模板图像和将新图像采用姿态参数h扭曲到标准模板后的图像之间差值的平方和最小。

$\underset{x}{\mathrm{\Σ}}{[I\left(w(x;h)\right)-T\left(x\right)]}^2---\left(10\right)$

为了对上式进行最优化，我们此时已经获得了初步的粗略姿态h，在此基础上通过对其增加一个增量Δh来进行迭代优化。下面的表达式便是迭代最小化的一步：

$\underset{x}{\mathrm{\Σ}}{[I\left(w(x;h+\mathrm{\Δh})\right)-T\left(x\right)]}^2---\left(11\right)$

根据参数增量Δh，我们可以将参数h进行更新：

h←h+Δh            (12)

这两步构成了迭代优化过程，通常我们通过判断参数更新向量Δh的二阶范数是否已经低于一个设定的阈值ε，即||Δh||≤ε时便认为该迭代过程已经收敛，并且参数h即我们预期的精确姿态。

至此，完成了本发明的图像匹配方法。

实例验证中采用的仿真环境为Visual C++2005，仿真计算机仿的硬件配置为

Processor Core^TM 2Duo 2.53GHz处理器，4GB内存。结果分析的对照组为1)采用本发明方法的目标匹配；2)采用PCA的单向描述符目标匹配方法。采用的Graffiti标准测试图在各个仿射参数添加随机的扰动之后进行仿射变换合成了一系列4000幅图像进行仿射变换不变性的对比实验，我们还采用目前公认比较出色的MSER进行特征检测，并建立SIFT描绘算子进行了对比。同时，也对算法所需的存储空间做了对比。

图5说明了本发明提出的算法和基于PCA的单向描述符算法的匹配速度对比。结果表明采用本算法进行匹配时，匹配帧速率为后者的两倍。

图6说明了不同算法在视角变化的情况下的匹配精度。在采用80个仿射姿态的情况下，采用DCT和仿射摄像机的描绘算子在视角变化较大的情况下匹配准确率要高于基于PCA的单向描述符。其中的Gepard即为采用PCA的单向描述符算法。

图7说明了不同算法所需的线下计算存储空间。实验数据可得，PCA描绘算子和DCT和仿射摄像机模型描绘算子所需存储空间和保留基个数直接的关系成线性，但是PCA描绘算子对应的斜率要明显大于DCT描绘算子。在实验中我们保留210个基即可获得非常鲁棒的匹配精度，在此情况下，本发明比PCA描绘算子可以节省22.4MB内存空间。

需要强调的是，本发明所述的实施例是说明性的，而不是限定性的，因此本发明并不限于具体实施方式中所述的实施例，凡是由本领域技术人员根据本发明的技术方案得出的其他实施方式，同样属于本发明保护的范围。

Claims (2)

1.一种基于DCT的特征描述算子及优化空间量化的图像匹配方法，其特征在于：包括以下步骤：

⑴对给定的N×N大小的图像块进行DCT空间的线性分解，得到基图像块的线性组合，该基图像块形式为：

$I=\underset{n=0}{\overset{N\×N-1}{\mathrm{\Σ}}}{C}_n{f}_n$

式中，C_n为DCT变换系数，f_n为DCT变换后的频域函数；

⑵对图像块进行仿射扭曲并对姿态空间进行优化，得到量化姿态数M；该优化方法包括以下步骤：

①在仿射摄像机模型下得到镜头拉升参数λ、摄像机的轴向旋转角ψ、摄像机距离目标物体等距离球面上的纬度角θ、经度角Φ；

②针对不同的纬度θ采用不同的量化精度，在低纬度采用小的经度步长，高纬度采用较大的经度步长；

③仿射量化点采用图像平面正方向的四分之一球面分布覆盖姿态变换，得到量化姿态数M，我们取量化级数为80；

⑶计算模板图像M个量化姿态DCT基图像块在不同投影系数下的组合，得到均值图像块

该均值图像块

的计算方法为：

$\stackrel{\‾}{p_{i,h}}\∝\frac{1}{M}\underset{j=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=0}{\overset{N\×N-1}{\mathrm{\Σ}}}w(C_n,f_n,H_{h,j})=\underset{n=0}{\overset{N\×N-1}{\mathrm{\Σ}}}{C}_n\stackrel{\‾}{f_n}$

式中：C_n为DCT变换系数，

为基图像块，M为量化姿态数；

⑷将模板图像和检测到的特征点周围的局部图像块进行图像匹配，包括以下步骤：

①将模板图像和检测到的特征点周围的局部图像块在初始姿态h下进行配准，使新图像块和进行姿态参数h扭曲后的图像之间差值的平方和最小，然后计算图像块和模板各个姿态的均值图像块的最小极值，并得到最接近的姿态h；

②添加增量姿态参数△h，将姿态参数h进行更新，判断其二阶范数是否低于一个设定的阈值ε，如果低于ε，则认为h已经收敛，即为预期的姿态。

2.根据权利要求1所述的基于DCT的特征描述算子及优化空间量化的图像匹配方法，其特征在于：所述的基图像块

采用如下公式进行计算：

$\stackrel{\‾}{f_n}=\frac{1}{M}\underset{j=1}{\overset{M}{\mathrm{\Σ}}}w(f_n,H_{h,j})$

该基图像块

采用线下预先计算并存储。

Images