CN102222228A - 影像特征点提取方法 - Google Patents

Abstract

本发明公开了一种影像特征点提取方法，该方法首先将数字图像G(m，n)进行多尺度多方向小波变换，分别得到各尺度下各方向影像；然后，分别提取各尺度下影像的特征点；最后选择各尺度下共同的特征点作为数字图像G(m，n)的特征点。该方法的特点是可抑制变换后高频影像中辐射强度变化，很大程度上消除影像不同方向高频信息变化所带来的影响，有利于提取由于影像尺度变化的局部特征。

Description

影像特征点提取方法

技术领域

本发明涉及数字摄影测量和计算机视觉等领域，特别是涉及一种对影像特征点进行提取的方法。

背景技术

影像特征点是由于目标物理与几何特征不同使影像中局部区域的灰度发生明显变化而形成的。特征点是影像灰度曲面不连续，影像灰度信息在局部区域发生明显变化的点集，是影像最基本的特征之一。数字影像中由于点扩展特性，特征表现为一个微小邻域中灰度的急剧变化或分布的均匀性，特征的存在使局部区域中有较大的信息量。影像点特征提取，是实现影像定标、立体匹配、目标描述与识别、运动估计、目标跟踪等重要基础。在数字摄影测量和计算机视觉领域，根据不同应用目标选择有效的点特征提取算子是非常重要的。

影像几何特征分为三类，即点特征、线特征和面特征。在数字影像中，点特征具有旋转不变和不随光照条件变化而变化的特点。在数字摄影测量和计算机视觉等领域应用中，利用点特征作为基本处理基元可以减少计算的数据量，同时又不损失影像的重要灰度信息。在立体影像匹配中，基于点特征的匹配不仅速度高而且匹配可靠性强。在影像目标描述与识别处理中，大量使用点特征作为影像基本处理基元可以简化处理算法提高处理效率。

特征点提取的算法或算子也称为兴趣算子或定位算子，即运用某种算法从影像中提取出感兴趣的点，并确定点的位置。目前常用的主要有

算子、Harris算子、SUSAN算子。

算子、是W.Forstner在1982提出的。该算法在数字摄影测量和机器视觉的立体匹配与目标描述与识别等处理中广泛应用，具有精度高、速度快的优点。

算子是基于误差椭圆理论构建的，首先计算每个像素的Robert梯度，以及以像素点为中心的窗口灰度协方差矩阵，基于误差椭圆的圆度为判定依据。但它的缺点是要确定阈值，因此受图像灰度、对比度变化的影响。Harris算子是C.Harris和M.J.Stephens在1988年提出的一种基于信号的点特征提取算子[36]。算法的思想是给出与自相关函数有关的矩阵，该矩阵的特征值是自相关函数的曲率，如果某点的x方向和y方向的曲率都很高，那么该点就取作特征点。Harris算子只用到一阶差分，计算简单，可以最大限度地提出局部特征点，惟一不足的是精度不高。SUSAN(Smallest Univalve Segment Assimilating Nucleus)算法，最早是由牛津大学的S.M.Smith和J.M.Brady提出的，用于检测图像中的角点。SUSAN算子的基本原理是：与每一图像点对应的局部区域都应具有相同的灰度值。如果某一窗口区域内的每一像素的灰度值与该窗口中心像素的灰度值相同或相近，则这一窗口区域将被称之为USAN(Univalve Segment Assimilating Nucleus)区域。计算图像中每一个像素的USAN，位于边缘上的像素的USAN较小，位于角点上的像素的USAN更小。Susan算法的特点：对角点的检测比对边缘检测的效果好，无需梯度运算，适用于基于角点的影像配准。在抗噪和计算速度方面有较大的改进。存在的缺点是：相似比较函数计算复杂；图像中不同区域处目标与背景的对比程度不一样，取固定阈值t不符合实际情况。

上述的特征点提取算法，在很大程度上受到影像的辐射强度变化、影像尺度变化以及噪声等因素的影响，不能有效、稳定地提取影像特征点。

发明内容

本发明旨在提供一种在辐射强度与尺度变化时，能够提取出稳定的、大量的特征点的影像特征点提取方法。

本发明的影像特征点提取算法的技术方案如下：

本发明的影像特征点提取方法，包括以下步骤：

S01：将数字图像G(m，n)进行多尺度多方向小波变换，分别得到各尺度下各方向影像；

S02：分别提取各尺度下影像的特征点；

S03：选择各尺度下共同的特征点作为数字图像G(m，n)的特征点，并保存。

优选的，还包括步骤S04：将特征点以区别于数字图像G(m，n)的颜色显示在数字图像上。

优选的，所述的区别于数字图像G(m，n)的颜色为红色。

优选的，所述将数字图像G(m，n)进行多尺度多方向小波变换中，多尺度多方向具体为三尺度四方向，所述的四方向为0°45°90°135°，分别得到三尺度四方向的影像。

优选的，进行三尺度四方向的小波变换，通过如下公式实现：

$\left\{\begin{array}{c}{A}_{2^{j+1}}G(l,k)=A_{2^j}G(m,n)*h(m,n)\\ D_{2^{j+1}}^{\mathrm{\α}}G(l,k)=A_{2^j}G(m,n)*g^{\mathrm{\α}}(m,n)\end{array}\right.---\left(1\right)$

式中，

表示尺度j+1的低频影像，

表示尺度j+1，方向为α的高频影像，h、g^α为空间滤波器器响应系数。

优选的，所述S02分别提取各尺度下影像的特征点，通过如下方法实现：

S11：依据尺度2高频影像及下式，计算该尺度影像梯度极大模对应像点位置，作为该尺度待选特征点。

$\sqrt{{\left(D_2^0\right)}^2+{\left(D_2^{90}\right)}^2}=\max ,$ $\sqrt{{\left(D_2^{45}\right)}^2+{\left(D_2^{135}\right)}^2}=\max ---\left(2\right)$

S12：依据尺度1高频影像及下式，计算该尺度影像梯度极大模对应像点位置，作为该尺度待选特征点。

$\sqrt{{\left(D_1^0\right)}^2+{\left(D_1^{90}\right)}^2}=\max ,$ $\sqrt{{\left(D_1^{45}\right)}^2+{\left(D_1^{135}\right)}^2}=\max ---\left(3\right)$

S13：根据以下四个公式，在尺度0低频影像上计算影像兴趣值点

dg₁＝∑|G_o(i，j)-G_o(i，j+1)|  (4)

dg₂＝∑|G_o(i，j)-G_o(i+1，j+1)|(5)

dg₃＝∑|G_o(i，j)-G_o(i+1，j+1)|(6)

dg₄＝∑|G_o(i，j)-G_o(i+1，j-1)|(7)

根据下式对得到的兴趣点进行非极大抑制，对应像点位置即为特征点位置。

v＝min(dg₁，dg₂，dg₃，dg₄)(8)

本发明的有益效果：

本发明的影像特征点提取方法可抑制变换后高频影像中辐射强度的变化，很大程度上消除影像不同方向高频信息变化所带来的影响，有利于提取由于影像尺度变化的局部特征，是一种是影像辐射强度与尺度不变性特征提取算法。解决了传统的点特征提取算法，在很大程度上受到影像的辐射强度变化、影像尺度变化以及噪声等因素的影响，不能有效、稳定地提取影像特征点的问题。

附图说明

图1为本发明的影像特征点提取方法一个实施例的流程图；

图2为3尺度四方向提取特征点时离散方向的小波变换图；

图3为一个具体实施例中尺度2影像特征点的提取结果；

图4为一个具体实施例中尺度1影像特征点的提取结果；

图5为一个具体实施例中尺度0影像特征点的提取结果；

图6为一个具体实施例中最后的提取结果；

图7为采用本发明的特征点提取方法对尺寸缩小30％所得图像的特征点的提取结果与原图提取结果的对比图；

图8为采用本发明的特征点提取方法对亮度下降后图像特征点的提取结果与原图提取结果的对比图；

图9为对比试验中采用本发明的方法的提取结果；

图10为对比试验中采用Harris算法的提取结果；

图11为对比试验中采用SUSAN算法的提取结果；

图12为对比试验中采用

法子的提取结果。

具体实施方式

下面结合附图和具体实施例对发明的技术方案和有益效果进一步进行说明。

本发明的方法可以提取多尺度多方向上的特征点，但是考虑计算效率，3尺度效果为好。

本发明的小波变换通过公式(1)实现：

$\left\{\begin{array}{c}{A}_{2^{j+1}}G(l,k)=A_{2^j}G(m,n)*h(m,n)\\ D_{2^{j+1}}^{\mathrm{\α}}G(l,k)=A_{2^j}G(m,n)*g^{\mathrm{\α}}(m,n)\end{array}\right.---\left(1\right)$

式中，

表示尺度j+1的低频影像，

表示尺度j+1，方向为α的高频影像，h、g^α为空间滤波器器响应系数。

公式(1)的推导过程如下：

首先进行二维方向小波构造

二维光滑函数φ(x，y)沿方向

的一阶偏导数为

${\mathrm{\ψ}}^{\mathrm{\α}}(x,y)=\frac{\∂\mathrm{\φ}(x,y)}{\∂\stackrel{\→}{n}}=\cos \mathrm{\α}\frac{\∂\mathrm{\φ}(x,y)}{\∂x}+\sin \mathrm{\α}\frac{\∂\mathrm{\φ}(x,y)}{\∂y}---\left(9\right)$

函数f(x，y)的方向小波变换为

Wf^α(x，y)＝ψ^α(x，y)*f(x，y)(10)

式中：*表示ψ^α(x，y)与f(x，y)进行卷积运算，α表示小波变换后高频信息方向。

选取高斯核为光滑函数，则(9)式可写为：

ψ^α(x，y)＝cosα·ψ(x)·φ(y)+sinα·φ(x)·ψ(y)

φ(x，y)、ψ^α(x，y)的付里叶变换为：

$\widehat{\mathrm{\φ}}({\mathrm{\ω}}_x,{\mathrm{\ω}}_y,\mathrm{\σ})=\widehat{\mathrm{\φ}}({\mathrm{\ω}}_x,\mathrm{\σ})\·\widehat{\mathrm{\φ}}({\mathrm{\ω}}_y,\mathrm{\σ})$

${\widehat{\mathrm{\ψ}}}^{\mathrm{\α}}({\mathrm{\ω}}_x,{\mathrm{\ω}}_y,\mathrm{\σ})=\cos \mathrm{\α}\·\widehat{\mathrm{\ψ}}({\mathrm{\ω}}_x,\mathrm{\σ})\·\widehat{\mathrm{\φ}}({\mathrm{\ω}}_y,\mathrm{\σ})+\sin \mathrm{\α}\·\widehat{\mathrm{\φ}}({\mathrm{\ω}}_x,\mathrm{\σ})\·\widehat{\mathrm{\ψ}}({\mathrm{\ω}}_y,\mathrm{\σ})$

式中：ω_x，ω_y为频率变量，σ为高斯函数的参数。

存在两尺度方程[-π，π]

$\left\{\begin{array}{c}\widehat{\mathrm{\φ}}(2{\mathrm{\ω}}_x,2{\mathrm{\ω}}_y,\mathrm{\σ})\stackrel{\·}{=}H({\mathrm{\ω}}_x,{\mathrm{\ω}}_y,\mathrm{\σ})\·\widehat{\mathrm{\φ}}({\mathrm{\ω}}_x,{\mathrm{\ω}}_y,\mathrm{\σ})\\ {\widehat{\mathrm{\ψ}}}^{\mathrm{\α}}(2{\mathrm{\ω}}_x,2{\mathrm{\ω}}_y,\mathrm{\σ})\stackrel{\·}{=}G({\mathrm{\ω}}_x,{\mathrm{\ω}}_y,\mathrm{\σ})\·\widehat{\mathrm{\φ}}({\mathrm{\ω}}_x,{\mathrm{\ω}}_y,\mathrm{\σ})\end{array}\right.---\left(11\right)$

由于：

$\widehat{\mathrm{\φ}}(\mathrm{\ω},\mathrm{\σ})=e^{-\frac{{\mathrm{\σ}}^2{\mathrm{\ω}}^2}{2}},$ $\widehat{\mathrm{\ψ}}(\mathrm{\ω},\mathrm{\σ})=\mathrm{j\ω}{e}^{-\frac{{\mathrm{\σ}}^2{\mathrm{\ω}}^2}{2}}$

$H(\mathrm{\ω},\mathrm{\σ})=e^{-\frac{3{\mathrm{\σ}}^2{\mathrm{\ω}}^2}{2}},$ $G(\mathrm{\ω},\mathrm{\σ})=-j2\mathrm{\ω}{e}^{-\frac{3{\mathrm{\σ}}^2{\mathrm{\ω}}^2}{2}}$

可得：

$\left\{\begin{array}{c}H({\mathrm{\ω}}_x,{\mathrm{\ω}}_y,\mathrm{\σ})=H({\mathrm{\ω}}_x,\mathrm{\σ})\·H({\mathrm{\ω}}_y,\mathrm{\σ})\\ G({\mathrm{\ω}}_x,{\mathrm{\ω}}_y,\mathrm{\σ})=\cos \mathrm{\α}\·G({\mathrm{\ω}}_x,\mathrm{\σ})H({\mathrm{\ω}}_y,\mathrm{\σ})+\sin \mathrm{\α}\·H({\mathrm{\ω}}_x,\mathrm{\σ})G({\mathrm{\ω}}_y,\mathrm{\σ})\end{array}\right.---\left(12\right)$

对公式(12)作付里叶逆变换，得到空间滤波器响应系数为

$h(m,n,\mathrm{\σ})=\frac{1}{4{\mathrm{\π}}^2}{\∫}_{-\mathrm{\π}}^{+\mathrm{\π}}({\∫}_{-\mathrm{\π}}^{+\mathrm{\π}}H({\mathrm{\ω}}_x,{\mathrm{\ω}}_y,\mathrm{\σ})\·e^{\mathrm{jm}{\mathrm{\ω}}_x}d{\mathrm{\ω}}_x)\·e^{\mathrm{jn}{\mathrm{\ω}}_y}d{\mathrm{\ω}}_y$

$g^{\mathrm{\α}}(m,n,\mathrm{\σ})=\frac{1}{4{\mathrm{\π}}^2}{\∫}_{-\mathrm{\π}}^{+\mathrm{\π}}({\∫}_{-\mathrm{\π}}^{+\mathrm{\π}}G({\mathrm{\ω}}_x,{\mathrm{\ω}}_y,\mathrm{\σ})\·e^{\mathrm{jm}{\mathrm{\ω}}_x}d{\mathrm{\ω}}_x)\·e^{\mathrm{jn}{\mathrm{\ω}}_y}d{\mathrm{\ω}}_y)$

即：

h(m，n，σ)＝h(m，σ)h(n，σ)

g^α(m，n，σ)＝cosα·g(m，σ)h(n，σ)+sinα·g(n，σ)h(m，σ)

则得影像函数G(m，n)的离散方向小波变换为

$\left\{\begin{array}{c}{A}_{2^{j+1}}G(l,k)=A_{2^j}G(m,n)*h(m,n)\\ D_{2^{j+1}}^{\mathrm{\α}}G(l,k)=A_{2^j}G(m,n)*g^{\mathrm{\α}}(m,n)\end{array}\right.---\left(1\right)$

式中，

表示尺度j+1的低频影像，

表示尺度j+1，方向为α的高频影像。参见附图1，本发明的小波变换提取特征点的方法，首先需要按照公式(1)进行小波变换，变换过程如图2所示，首先从从0尺度低频影像G0分解出1尺度低频影像G1，四个方向高频影像

然后继续从1尺度低频影像G1分解出2尺度低频影像G₂和2尺度4个方向的高频影像

然后分别提取个尺度图像的特征点，具体步骤如下：

S11：依据尺度2高频影像及下式，计算该尺度影像梯度极大模对应像点位置，作为该尺度待选特征点。

$\sqrt{{\left(D_2^0\right)}^2+{\left(D_2^{90}\right)}^2}=\max ,$ $\sqrt{{\left(D_2^{45}\right)}^2+{\left(D_2^{135}\right)}^2}=\max ---\left(2\right)$

S12：依据尺度1高频影像及下式，计算该尺度影像梯度极大模对应像点位置，作为该尺度待选特征点。

$\sqrt{{\left(D_1^0\right)}^2+{\left(D_1^{90}\right)}^2}=\max ,$ $\sqrt{{\left(D_1^{45}\right)}^2+{\left(D_1^{135}\right)}^2}=\max ---\left(3\right)$

S13：根据以下四个公式，在尺度0低频影像上计算影像兴趣值点

dg₁＝∑|G_o(i，j)-G_o(i，j+1)|  (4)

dg₂＝∑|G_o(i，j)-G_o(i+1，j)|  (5)

dg₃＝∑|G_o(i，j)-G_o(i+1，j+1)|(6)

dg₄＝∑|G_o(i，j)-G_o(i+1，j-1)|(7)

根据下式对得到的兴趣点进行非极大抑制，对应像点位置即为特征点位置。

v＝min(dg₁，dg₂，dg₃，dg₄)(8)

附图3、4、5分别为尺度2、尺度1、尺度0的影像特征点的提取结果，图6为最后的提取结果。

图7为采用本发明的特征点提取方法对尺寸缩小30％所得图像的特征点的提取结果与原图提取结果的对比图，其中左图为原图，提取特征点7448个，右图为左图缩小至0.6倍时的图像，提取特征点7473个。

图8为采用本发明的特征点提取方法对亮度下降后图像特征点的提取结果与原图提取结果的对比，其中左图为原图，提取特征点7448个，右图为左图亮度下降后的图像，提取特征点4967个。

图7、8说明了本发明的提取方法具有较强的抗辐射变化和抗尺度变化的特点，在这两种情况下依然能够稳定的提取相当数量的特征点。

综上，本发明的特征点提取方法，首先利用高斯函数作为光滑函数，对高斯光滑函数沿方向矢量求一阶导数，得到方向离散小波。基于离散方向小波对影像进行多尺度变换，构造多尺度影像特征算法。其特点是可抑制变换后高频影像中辐射强度变化。与传统小波变换相比较，本算法另一特点是，除了提取影像水平和竖直方向高频信息，还可提取其它任意方向的高频信息(本算法提取0°45°90°135°四方向高频信息)，这样可在很大程度上消除影像不同方向高频信息变化所带来的影响；第三，本算法引入了多尺度的提取思想，利用小波变换在三尺度下提取影像特征点，有利于提取由于影像尺度变化的局部特征。

为了更好的说明本发明的特征点提取方法的效果，与其它数字摄影测量提取方法的提取结果进行了比较，具体如下：

实验影像为某地无人机平台数字影像，佳能EOS400D像机，影像大小3888像素×2592像素，选取局部影像进行实验(800像素×800像素)。

与传统的Harris算法、法子和SUSAN算法做了对比，图9至12分别显示了采用本发明的方法、Harris算法、

法子和SUSAN算法的提取结果，提取效果和对比结果如下：

表1实验比较

实验采用了小飞机拍摄的某地的影像，该影像有着较为明显的特征变化。从提取的效果来看，Forstner算法较为精准的提取了影像灰度变化的特征点，但数量过少，忽略了较多的可用特征；Harris算法提取的特征点分布较为均匀，也能较好的反应影像中的特征，比较而言其提取的特征点还是过少；SUSAN算法在较短的时间里提取了大量的特征点，但其中的冗余信息过多，其特征集合不适合作为匹配的初始集合，会造成匹配效率的下降；而本文的本方法在5.1秒的时间内提取了近7000个特征点，并且分布非常均匀，同时很好的体现了局部影像特征，在经过影像匹配之后仍然能够保证有足够的匹配点对剩余，相比前面三种方法具有更高的特征点提取效率，且具有尺度与辐射强度不变性。

Claims (6)

1.影像特征点提取方法，其特征在于包括以下步骤：

S01：将数字图像G(m，n)进行多尺度多方向小波变换，分别得到各尺度下各方向影像；

S02：分别提取各尺度下影像的特征点；

S03：选择各尺度下共同的特征点作为数字图像G(m，n)的特征点，并保存。

2.根据权利要求1所述的影像特征点提取方法，其特征在于：还包括步骤S04：将特征点以区别于数字图像G(m，n)的颜色显示在数字图像上。

3.根据权利要求3所述的影像特征点提取方法，其特征在于：所述的区别于数字图像G(m，n)的颜色为红色。

4.根据权利要求1所述的影像特征点提取方法，其特征在于：所述将数字图像G(m，n)进行多尺度多方向小波变换中，多尺度多方向具体为三尺度四方向，所述的四方向为0°45°90°135°，分别得到三尺度四方向的影像。

5.根据权利要求4所述的影像特征点的提取方法，其特征在于：进行三尺度四方向的小波变换，通过如下公式实现：

$\left\{\begin{array}{c}{A}_{2^{j+1}}G(l,k)=A_{2^j}G(m,n)*h(m,n)\\ D_{2^{j+1}}^{\mathrm{\α}}G(l,k)=A_{2^j}G(m,n)*g^{\mathrm{\α}}(m,n)\end{array}\right.$

式中，

表示尺度j+1的低频影像，表示尺度j+1，方向为α的高频影像，h、g^a为空间滤波器响应系数，*代表进行卷积运算。

6.根据权利要求5所述的影像特征点提取方法，其特征在于：所述S02分别提取各尺度下影像的特征点，通过如下方法实现：

S11：依据尺度2高频影像及下式，计算该尺度影像梯度极大模对应像点位置，作为该尺度待选特征点；

$\sqrt{{\left(D_2^0\right)}^2+{\left(D_2^{90}\right)}^2}=\max ,$ $\sqrt{{\left(D_2^{45}\right)}^2+{\left(D_2^{135}\right)}^2}=\max$

S12：依据尺度1高频影像及下式，计算该尺度影像梯度极大模对应像点位置，作为该尺度待选特征点；

$\sqrt{{\left(D_1^0\right)}^2+{\left(D_1^{90}\right)}^2}=\max ,$ $\sqrt{{\left(D_1^{45}\right)}^2+{\left(D_1^{135}\right)}^2}=\max$

S13：根据以下四个公式，在尺度0低频影像上计算影像兴趣值点

dg₁＝∑|G_o(i，j)-G_o(i，j+1)|

dg₂＝∑|G_o(i，j)-G_o(i+1，j)|

dg₃＝∑|G_o(i，j)-G_o(i+1，j+1)|

dg₄＝∑|+G_o(i，j)-G_o(i+1，j-1)|

根据下式对得到的兴趣点进行非极大抑制，对应像点位置即为特征点位置。

v＝min(dg₁，dg₂，dg₃，dg₄)

Images