代码改变世界

Harris角点

2014-10-09 17:51 ☆Ronny丶 阅读(...) 评论(...) 编辑 收藏

1. 不同类型的角点

在现实世界中,角点对应于物体的拐角,道路的十字路口、丁字路口等。从图像分析的角度来定义角点可以有以下两种定义:

  1. 角点可以是两个边缘的角点;
  2. 角点是邻域内具有两个主方向的特征点;

前者往往需要对图像边缘进行编码,这在很大程度上依赖于图像的分割与边缘提取,具有相当大的难度和计算量,且一旦待检测目标局部发生变化,很可能导致操作的失败。早期主要有Rosenfeld和Freeman等人的方法,后期有CSS等方法。

基于图像灰度的方法通过计算点的曲率及梯度来检测角点,避免了第一类方法存在的缺陷,此类方法主要有Moravec算子、Forstner算子、Harris算子、SUSAN算子等。

image

这篇文章主要介绍的Harris角点检测的算法原理,比较著名的角点检测方法还有jianbo Shi和Carlo Tomasi提出的Shi-Tomasi算法,这个算法开始主要是为了解决跟踪问题,用来衡量两幅图像的相似度,我们也可以把它看为Harris算法的改进。OpenCV中已经对它进行了实现,接口函数名为GoodFeaturesToTrack()。另外还有一个著名的角点检测算子即SUSAN算子,SUSAN是Smallest Univalue Segment Assimilating Nucleus(最小核值相似区)的缩写。SUSAN使用一个圆形模板和一个圆的中心点,通过圆中心点像素与模板圆内其他像素值的比较,统计出与圆中心像素近似的像元数量,当这样的像元数量小于某一个阈值时,就被认为是要检测的角点。我觉得可以把SUSAN算子看为Harris算法的一个简化。这个算法原理非常简单,算法效率也高,所以在OpenCV中,它的接口函数名称为:FAST()

2. Harris角点

2.1 基本原理

人眼对角点的识别通常是在一个局部的小区域或小窗口完成的。如果在各个方向上移动这个特征的小窗口,窗口内区域的灰度发生了较大的变化,那么就认为在窗口内遇到了角点。如果这个特定的窗口在图像各个方向上移动时,窗口内图像的灰度没有发生变化,那么窗口内就不存在角点;如果窗口在某一个方向移动时,窗口内图像的灰度发生了较大的变化,而在另一些方向上没有发生变化,那么,窗口内的图像可能就是一条直线的线段。

image

对于图像$I(x,y)$,当在点$(x,y)$处平移$(\Delta x,\Delta y)$后的自相似性,可以通过自相关函数给出:

$$c(x,y;\Delta x,\Delta y) = \sum_{(u,v)\in W(x,y)}w(u,v)(I(u,v) – I(u+\Delta x,v+\Delta y))^2$$

其中,$W(x,y)$是以点$(x,y)$为中心的窗口,$w(u,v)$为加权函数,它既可是常数,也可以是高斯加权函数。

image

根据泰勒展开,对图像$I(x,y)$在平移$(\Delta x,\Delta y)$后进行一阶近似:

$$I(u+\Delta x,v+\Delta y) = I(u,v)+I_x(u,v)\Delta x+I_y(u,v)\Delta y+O(\Delta x^2,\Delta y^2)\approx I(u,v)+I_x(u,v)\Delta x+I_y(u,v)\Delta y$$

其中,$I_x,I_y$是图像$I(x,y)$的偏导数,这样的话,自相关函数则可以简化为:

$$c(x,y;\Delta x,\Delta y)\approx \sum_w (I_x(u,v)\Delta x+I_y(u,v)\Delta y)^2=[\Delta x,\Delta y]M(x,y)\begin{bmatrix}\Delta x \\ \Delta y\end{bmatrix}$$

其中

$$M(x,y)=\sum_w \begin{bmatrix}I_x(x,y)^2&I_x(x,y)I_y(x,y) \\ I_x(x,y)I_y(x,y)&I_y(x,y)^2\end{bmatrix} = \begin{bmatrix}\sum_w I_x(x,y)^2&\sum_w I_x(x,y)I_y(x,y) \\\sum_w I_x(x,y)I_y(x,y)&\sum_w I_y(x,y)^2\end{bmatrix}=\begin{bmatrix}A&C\\C&B\end{bmatrix} $$

也就是说图像$I(x,y)$在点$(x,y)$处平移$(\Delta x,\Delta y)$后的自相关函数可以近似为二项函数:

$$c(x,y;\Delta x,\Delta y)\approx A\Delta x^2+2C\Delta x\Delta y+B\Delta y^2$$

其中

$$A=\sum_w I_x^2, B=\sum_w I_y^2,C=\sum_w I_x I_y$$

二次项函数本质上就是一个椭圆函数。椭圆的扁率和尺寸是由$M(x,y)$的特征值$\lambda_1、\lambda_2$决定的,椭贺的方向是由$M(x,y)$的特征矢量决定的,如下图所示,椭圆方程为:

$$[\Delta x,\Delta y]M(x,y)\begin{bmatrix}\Delta x \\ \Delta y\end{bmatrix} = 1$$

image

椭圆函数特征值与图像中的角点、直线(边缘)和平面之间的关系如下图所示。共可分为三种情况:

  • 图像中的直线。一个特征值大,另一个特征值小,$\lambda_1\gg \lambda_2$或$\lambda_2\gg \lambda_1$。自相关函数值在某一方向上大,在其他方向上小。
  • 图像中的平面。两个特征值都小,且近似相等;自相关函数数值在各个方向上都小。
  • 图像中的角点。两个特征值都大,且近似相等,自相关函数在所有方向都增大。

    image

根据二次项函数特征值的计算公式,我们可以求$M(x,y)$矩阵的特征值。但是Harris给出的角点差别方法并不需要计算具体的特征值,而是计算一个角点响应值$R$来判断角点。$R$的计算公式为:

$$R=det \boldsymbol{M} - \alpha(trace\boldsymbol{M})^2$$

式中,$det\boldsymbol{M}$为矩阵$\boldsymbol{M}=\begin{bmatrix}A&B\\B&C\end{bmatrix}$的行列式;$trace\boldsymbol{M}$为矩阵$\boldsymbol{M}$的直迹;$\alpha$为经常常数,取值范围为0.04~0.06。事实上,特征是隐含在$det\boldsymbol{M}$和$trace\boldsymbol{M}$中,因为,

$$det\boldsymbol{M} = \lambda_1\lambda_2=AC-B^2$$

$$trace\boldsymbol{M}=\lambda_2+\lambda_2=A+C$$

2.2 Harris角点算法实现

根据上述讨论,可以将Harris图像角点检测算法归纳如下,共分以下五步:

1. 计算图像$I(x,y)$在$X$和$Y$两个方向的梯度$I_x、I_y$。

$$I_x=\frac{\partial I}{\partial x}=I\otimes(-1\ 0\ 1),I_y =\frac{\partial I}{\partial x}=I\otimes(-1\ 0\ 1)^T $$

2. 计算图像两个方向梯度的乘积。

$$I_x^2=I_x\cdot I_y,I_y^2=I_y\cdot I_y,I_{xy}=I_x\cdot I_y$$

3. 使用高斯函数对$I_x^2、I_y^2和I_{xy}$进行高斯加权(取$\sigma=1$),生成矩阵$M$的元素$A、B$和$C$。

$$A=g(I_x^2)=I_x^2\otimes w,C=g(I_y^2)=I_y^2\otimes w,B=g(I_{x,y})=I_{xy}\otimes w$$

4. 计算每个像素的Harris响应值$R$,并对小于某一阈值$t$的$R$置为零。

$$R=\{R:det \boldsymbol{M} - \alpha(trace\boldsymbol{M})^2<t\}$$

5. 在$3\times3$或$5\times5$的邻域内进行非最大值抑制,局部最大值点即为图像中的角点。

Harris角点检测的C++实现代码:https://github.com/RonnyYoung/ImageFeatures/blob/master/source/harris.cpp

2.3 Harris角点的性质

1. 参数$\alpha$对角点检测的影响

假设已经得到了矩阵$\boldsymbol{M}$的特征值$\lambda_1\ge\lambda_2\ge0$,令$\lambda_2=k\lambda_1,0\le k\le 1$。由特征值与矩阵$\boldsymbol{M}$的直迹和行列式的关系可得:

$$det\boldsymbol{M}=\prod_i\lambda_i \ \ \ \ \ \  trace\boldsymbol{M}=\sum_i\lambda_i$$

从而可以得到角点的响应

$$R=\lambda_2\lambda_2=\alpha(\lambda_2+\lambda_2)^2=\lambda^2(k-\alpha(1+k)^2)$$

假设$R\ge0$,则有:

$$0\le \alpha \le\frac{k}{(1+k)^2}\le0.25$$

对于较小的$k$值,$R\approx\lambda^2(k-\alpha),\alpha<k$。

由此,可以得出这样的结论:增大$\alpha$的值,将减小角点响应值$R$,降低角点检测的灵性,减少被检测角点的数量;减小$\alpha$值,将增大角点响应值$R$,增加角点检测的灵敏性,增加被检测角点的数量。

2. Harris角点检测算子对亮度和对比度的变化不敏感

这是因为在进行Harris角点检测时,使用了微分算子对图像进行微分运算,而微分运算对图像密度的拉升或收缩和对亮度的抬高或下降不敏感。换言之,对亮度和对比度的仿射变换并不改变Harris响应的极值点出现的位置,但是,由于阈值的选择,可能会影响角点检测的数量。

image image

3. Harris角点检测算子具有旋转不变性

Harris角点检测算子使用的是角点附近的区域灰度二阶矩矩阵。而二阶矩矩阵可以表示成一个椭圆,椭圆的长短轴正是二阶矩矩阵特征值平方根的倒数。当特征椭圆转动时,特征值并不发生变化,所以判断角点响应值$R$也不发生变化,由此说明Harris角点检测算子具有旋转不变性。

4. Harris角点检测算子不具有尺度不变性

如下图所示,当右图被缩小时,在检测窗口尺寸不变的前提下,在窗口内所包含图像的内容是完全不同的。左侧的图像可能被检测为边缘或曲线,而右侧的图像则可能被检测为一个角点。

image

2.4 Harris的OpenCV接口

OpenCV的Hairrs角点检测的函数为cornerHairrs(),但是它的输出是一幅浮点值图像,浮点值越高,表明越可能是特征角点,我们需要对图像进行阈值化。

C++: void cornerHarris(InputArray src, OutputArray dst, int blockSize, int apertureSize, double k, int borderType = BORDER_DEFAULT);
  • src – 输入的单通道8-bit或浮点图像。
  • dst – 存储着Harris角点响应的图像矩阵,大小与输入图像大小相同,是一个浮点型矩阵。
  • blockSize – 邻域大小。
  • apertureSize – 扩展的微分算子大。
  • k – 响应公式中的,参数$\alpha$。
  • boderType – 边界处理的类型。
int main()
{
    Mat image = imread("../buliding.png");
    Mat gray;
    cvtColor(image, gray, CV_BGR2GRAY);

    Mat cornerStrength;
    cornerHarris(gray, cornerStrength, 3, 3, 0.01);
    threshold(cornerStrength, cornerStrength, 0.001, 255, THRESH_BINARY);
    return 0;
}

 

  image     image   image

从上面上间一幅图像我们可以看到,有很多角点都是粘连在一起的,我们下面通过加入非极大值抑制来进一步去除一些粘在一起的角点。

非极大值抑制原理是,在一个窗口内,如果有多个角点则用值最大的那个角点,其他的角点都删除,窗口大小这里我们用3*3,程序中通过图像的膨胀运算来达到检测极大值的目的,因为默认参数的膨胀运算就是用窗口内的最大值替代当前的灰度值。

int main()
{
    Mat image = imread("buliding.png");
    Mat gray;
    cvtColor(image, gray, CV_BGR2GRAY);

    Mat cornerStrength;
    cornerHarris(gray, cornerStrength, 3, 3, 0.01);

    double maxStrength;
    double minStrength;
    // 找到图像中的最大、最小值
    minMaxLoc(cornerStrength, &minStrength, &maxStrength);

    Mat dilated;
    Mat locaMax;
    // 膨胀图像,最找出图像中全部的局部最大值点
    dilate(cornerStrength, dilated, Mat());
    // compare是一个逻辑比较函数,返回两幅图像中对应点相同的二值图像
    compare(cornerStrength, dilated, locaMax, CMP_EQ);

    Mat cornerMap;
    double qualityLevel = 0.01;
    double th = qualityLevel*maxStrength; // 阈值计算
    threshold(cornerStrength, cornerMap, th, 255, THRESH_BINARY);
    cornerMap.convertTo(cornerMap, CV_8U);
    // 逐点的位运算
    bitwise_and(cornerMap, locaMax, cornerMap);

    drawCornerOnImage(image, cornerMap);
    namedWindow("result");
    imshow("result", image);
    waitKey();

    return 0;
}
void drawCornerOnImage(Mat& image, const Mat&binary)
{
    Mat_<uchar>::const_iterator it = binary.begin<uchar>();
    Mat_<uchar>::const_iterator itd = binary.end<uchar>();
    for (int i = 0; it != itd; it++, i++)
    {
        if (*it)
            circle(image, Point(i%image.cols, i / image.cols), 3, Scalar(0, 255, 0), 1);
    }
}

现在我们得到的效果就比默认的函数得到的结果有相当的改善,如上面最右边效果图。

3. 多尺度Harris角点

3.1 多尺度Harris角点的原理

虽然Harris角点检测算子具有部分图像灰度变化的不变性和旋转不变性,但它不具有尺度不变性。但是尺度不变性对图像特征来说至关重要。人们在使用肉眼识别物体时,不管物体远近,尺寸的变化都能认识物体,这是因为人的眼睛在辨识物体时具有较强的尺度不变性。在图像特征提取:尺度空间理论这篇文章里就已经讲到了高斯尺度空间的概念。下面将Harris角点检测算子与高斯尺度空间表示相结合,使用Harris角点检测算子具有尺度的不变性。

仿照Harris角点检测中二阶矩的表示方法,使用$M=\mu(x,\sigma_I,\sigma_D)$为尺度自适应的二阶矩:

$$\boldsymbol{M}=\mu(x,\sigma_I,\sigma_D)=\sigma_D^2g(\sigma_I)\otimes\begin{bmatrix}L_x^2(x,\sigma_D)&L_xL_y(x,\sigma_D)\\L_xL_y(x,\sigma_D)&L_y^2(x,\sigma_D)\end{bmatrix} $$

其中,$g(\sigma_I)$表示尺度为$sigma_I$的高斯卷积核,$x$表示图像的位置。与高斯测度空间类似,使用$L(x)$表示经过高斯平滑后的图像,符号$\otimes$表示卷积,$L_x(x,\sigma_D)$和$L_y(x,\sigma_D)$表示对图像使用高斯$g(\sigma_D)$函数进行平滑后,在$x$或$y$方向取其微分的结果,即$L_x=\partial_xL$和$L_y=\partial_yL$。通常将$\sigma_I$称为积分尺度,它是决定Harris角点当前尺度的变量,$\sigma_D$为微分尺度或局部尺度,它是决定角点附近微分值变化的变量。显然,积分尺度$\sigma_I$应该大于微分尺度$\sigma_D$。

3.2 多尺度Harris角点实现

首先,检测算法从预先定义的一组尺度中进行积分尺度搜索,这一组尺度定义为

$$\sigma_1\dots\sigma_n = \sigma_0\dots k^n\sigma_0$$

一般情况下使用k=1.4。为了减少搜索的复杂性,对于微分尺度$\sigma_D$的选择,我们采用在积分尺度的基础上,乘以一个比例常数,即$\sigma_D=s\sigma_I$,一般取s=0.7。这样,通常使用积分和微分的尺度,便可以生成$\mu(x,\sigma_I,\sigma_D)$,再利用Harris角点判断准则,对角点进行搜索,具体可以分两步进行。

1. 与Harris角点搜索类似,对于给定的尺度空间值$\sigma_D$,进行如下角点响应值计算和判断:

$$cornerness = det(\mu(x,\sigma_n)-\alpha trace^2(\mu(x,\sigma_n)))>threshold_H$$

2. 对于满足1中条件的点,在点的8邻域内进行角点响应最大值搜索(即非最大值抑制)出在8邻域内角点响应最大值的点。对于每个尺度$\sigma_n(1,2,\dots,n)$都进行如上搜索。

由于位置空间的候选点并不一定在尺度空间上也能成为候选点,所以,我们还要在尺度空间上进行搜索,找到该点的所谓特征尺度值。搜索特征尺度值也分两步。

1. 对于位置空间搜索到的每个候选点,进行拉普拉斯响应计算,并满足其绝对值大于给定的阈值条件:

$$F(x,\sigma_n) = \sigma_n^2|L_{xx}(x,\sigma_n)+L_{yy}(x,\sigma_n)| \ge threshold_L$$

2. 与邻近的两个尺度空间的拉普拉斯响应值进行比较,使其满足:

$$F(x,\sigma_n)>F(x,\sigma_l), \ \ \ l\in \{n-1.n+1\}$$

满足上述条件的尺度值就是该点的特征尺度值。这样,我们就找到了在位置空间和尺度空间都满足条件的Harris角点。

多尺度Harris角点检测C++实现:https://github.com/RonnyYoung/ImageFeatures/blob/master/source/harrisLaplace.cpp

4. 更多的讨论

在上面描述的Harris角点具有光照不变性、旋转不变性、尺度不变性,但是严格意义上来说并不具备仿射不变性。Harris-Affine是一种新颖的检测仿射不变特征点的方法,可以处理明显的仿射变换,包括大尺度变化和明显的视角变化。Harris-Affine主要是依据了以下三个思路:

  1. 特征点周围的二阶矩的计算对区域进行的归一化,具有仿射不变性;
  2. 通过在尺度空间上归一化微分的局部极大值求解来精化对应尺度;
  3. 自适应仿射Harris检测器能够精确定位牲点;

这篇文章不对Harris-Affine作进一步的描述,有时间会对这一算法做专门的分析,有兴趣的可以参考原论文:Scale & Affine Invariant Interest Point Detectors.

5. 参考资料

[1] 《图像局部不变特征与描述》王永明,王贵锦。

[2] Harris角点及Shi-Tomasi角点检测

[3] 图像特征提取PPT

[4] Harris角点检测算法 1

[5] OpenCV Harris角点检测

[6] Opencv学习笔记(五)Harris角点检测