特征点检测之Harris法

在计算机视觉中，特征点的概念被大量用于解决物体识别、图像匹配、视觉跟踪、三维重建等问题，比如图像中物体的角点，它们是在图像中可被轻易而精确地定位的二维特征。顾名思义，特征点检测的思想是无需观察整幅图像，而是通过选择某些特殊点，然后对它们执行局部分析。如果能检测到足够多的这种点，同时它们的区分度很高，并且可以精确定位稳定的特征，那么这个方法就很有效。这里主要使用Harris特征检测器检测图像角点。开发平台为Qt5.3.2+OpenCV2.4.9。

在此之前，先给出OpenCV中cv::cornerHarris函数的调用方式：

    cv::cornerHarris(image, // 输入图像
                     cornerStrength, // 输出为表示角点强度的32位浮点图像
                     3, // 导数平滑的相邻像素的尺寸
                     3, // 梯度计算的滤波器孔径大小
                     0.01); // Harris参数
 
 
1
2
3
4
5

描述Harris算子的经典论文可参考：

The article by C.Harris, M.J. Stephens, A Combined Corner and Edge Detector, Alvey Vision Conference, pp.147-152, 1988

The article by J. Shi and C. Tomasi, Good features to track, Int. Conference on Computer Vision and Pattern Recognition, pp. 593-600, 1994

The article by K. Mikolajczyk and C. Schmid, Scale and Affine invariant interest point detectors, International Journal of Computer Vision, vol 60, no 1, pp. 63-86, 2004

在第一篇论文中提到，Harris角点检测最直观的解释是在任意两个相互垂直的方向上，都有较大变化的点。为了定义一幅图像中的角点，Harris观察一个假定的特征点周围小窗口内的方向性强度平均变化。

如图所示，假设一个小窗口在图像上移动，在平滑区域如左图所示，窗口在各个方向上均没有变化。对于中间图，小窗口在边缘的方向上移动时也没有变化。而知道小窗口移动到右图的角点处，窗口在各个方向上均有明显的变化。Harris角点检测正是利用了这个直观的物理现象，通过窗口在各个方向上的变化程度，决定是否存在着角点。这里我们考虑偏移量(u,v)，则将图像窗口平移(u,v)产生的E(u,v)可表示为：

由以下公式可得到E(u,v)：

对于局部微小的移动，E(u,v)可近似表达为：

其中M的详细表达式为：

E(u,v)的椭圆形式如下：

E(u,v)是一个协方差矩阵，表现的是所有方向上强度的变化率。该定义涉及图像的一阶导数，这通常是Sobel算子的计算结果。而在OpenCV中cv::cornerHarris函数的第四个参数对应的正是用于计算Sobel滤波器的孔径(aperture)。协方差的两个特征值给出了最大平均强度变化以及垂直方向上的平均强度变化，如果这两个特征值均较低，就认为当前是同质区域；如果其中一个特征值较高，另外一个较低，则认为当前位于边缘上；最后，若两个特征值均较高，则判定当前位于角点处。

因此，定义角点响应函数R，其中k为函数cv::cornerHarris中的最后一个参数；之后，对R进行阈值处理，设定若R大于阈值threshold，则提取出局部极大值：

Harris角点的更多的理论部分可见： 
http://blog.csdn.net/lu597203933/article/details/15088485 
http://blog.csdn.net/xiaowei_cqu/article/details/7805206。

下面记录一下harris角点检测的几种方案。

一、基本的Harris角点检测实现

#include <QCoreApplication>
#include <opencv2/core/core.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/highgui/highgui.hpp>

int main(int argc, char *argv[])
{
    QCoreApplication a(argc, argv);
    // Harris法的简单实现
    cv::Mat image = cv::imread("c:/031.jpg", 0);
    cv::Mat cornerStrength;
    cv::cornerHarris(image,
                     cornerStrength, // 输出为表示角点强度的32位浮点图像
                     3, // 导数平滑的相邻像素的尺寸
                     3, // 梯度计算的滤波器孔径大小
                     0.01); // Harris的相关参数

    // 要在窗口输出，需要转化为CV_8U格式，阈值化即可
    // 角点强度的阈值
    cv::Mat harrisCorner;
    double threshold = 0.0001;
    cv::threshold(cornerStrength,
                  harrisCorner,
                  threshold,
                  255,
                  cv::THRESH_BINARY_INV); // 输出为翻转的二值图像

    cv::namedWindow("Original Image");
    cv::imshow("Original Image", image);
    cv::namedWindow("Harris Corner");
    cv::imshow("Harris Corner", harrisCorner);

    return a.exec();
}
 
 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

得到的结果为二值图像，可以看到图像中角点的位置包含许多圆圈，这与精确定位特征点的目标相悖：

二、改进的Harris角点检测实现(Shi-Tomasi算法)

这里通过封装自定义类来改进角点检测的效果。定义一个类HarrisDetector（其中已封装了Harris参数和相关函数）：

#ifndef HARRISDETECTOR_H
#define HARRISDETECTOR_H
#include <opencv2/core/core.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/highgui/highgui.hpp>

class HarrisDetector
{
private:
    // 表示图像角点强度的32位浮点图像
    cv::Mat cornerStrength;
    // 将输出图像阈值化后的32位浮点图像
    cv::Mat cornerThreshold;
    // 局部极大值图像
    cv::Mat localMax;

    // 以下三个为cornerHarris函数的必要参数
    // 导数平滑的相邻像素的尺寸
    int neighbourhoodPixelSize;
    // 滤波器的孔径大小
    int aperture;
    // Harris参数
    double k;

    // 阈值计算的最大强度
    double maxStrength;
    // 计算得到的阈值
    double threshold;
    // 非极大值抑制的相邻像素的尺寸
    int noneighbourhoodPixelSize;
    // 非极大值抑制的核
    cv::Mat kernel;

public:
    // 初始化参数
    HarrisDetector():neighbourhoodPixelSize(3),
                     aperture(3),
                     k(0.01),
                     maxStrength(0.0),
                     threshold(0.01),
                     noneighbourhoodPixelSize(3)
    {
        setLocalMaxWindowSize(noneighbourhoodPixelSize);
    }

    // 创建非极大值抑制的核
    void setLocalMaxWindowSize(int size);

    // 计算Harris角点
    void detect(const cv::Mat &image);

    // 由Harris的值获取角点图
    cv::Mat getCornerMap(double qualityLevel);

    // 由Harris的值获取特征点
    void getCorners(std::vector<cv::Point> &points, double qualityLevel);

    // 由角点图获取特征点
    void getCorners(std::vector<cv::Point> &points, const cv::Mat& cornerMap);

    // 在特征点的位置绘制图
    void drawOnImage(cv::Mat &image,
                     const std::vector<cv::Point> &points,
                     cv::Scalar color = cv::Scalar(255, 255, 255),
                     int radius = 4, int thickness = 2);
};


#endif // HARRISDETECTOR_H
 
 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69

接着，在harrisdetector.cpp中定义各个函数和初始化的参数：

#include "harrisdetector.h"

// 创建非极大值抑制的核
void HarrisDetector::setLocalMaxWindowSize(int size)
{
    noneighbourhoodPixelSize = size;
    kernel.create(noneighbourhoodPixelSize, noneighbourhoodPixelSize, CV_8U);
}

// 计算Harris角点
void HarrisDetector::detect(const cv::Mat &image)
{
    // Harris计算
    cv::cornerHarris(image,cornerStrength,
                     neighbourhoodPixelSize,
                     aperture,
                     k);
    // 内部阈值计算
    double minStrength; // 未使用
    cv::minMaxLoc(cornerStrength,
                  &minStrength,
                  &maxStrength);
    // 局部极大值检测
    cv::Mat dilate; // 临时图像
    cv::dilate(cornerStrength, dilate, cv::Mat());
    cv::compare(cornerStrength, dilate, localMax, cv::CMP_EQ);
}

// 由Harris的值获取角点图
cv::Mat HarrisDetector::getCornerMap(double qualityLevel)
{
    cv::Mat cornerMap;
    // 对角点图像进行阈值化
    threshold = qualityLevel * maxStrength;
    cv::threshold(cornerStrength, cornerThreshold,
                  threshold,255,cv::THRESH_BINARY);
    // 转换为8位图像
    cornerThreshold.convertTo(cornerMap, CV_8U);
    // 非极大值抑制
    cv::bitwise_and(cornerMap, localMax, cornerMap);
    return cornerMap;
}

// 由Harris的值获取特征点
void HarrisDetector::getCorners(std::vector<cv::Point> &points, double qualityLevel)
{
    // 得到角点图
    cv::Mat cornerMap = getCornerMap(qualityLevel);
    getCorners(points, cornerMap);
}

// 由角点图获取特征点
void HarrisDetector::getCorners(std::vector<cv::Point> &points, const cv::Mat& cornerMap)
{
    // 遍历像素得到所有特征
    for( int y = 0; y < cornerMap.rows; y++ )
    {
        const uchar* rowPtr = cornerMap.ptr<uchar>(y);
        for( int x = 0; x < cornerMap.cols; x++ )
        {
            // 如果是特征点
            if (rowPtr[x])
            {
                points.push_back(cv::Point(x,y));
            }
        }
    }
}

// 在特征点的位置绘制图
void HarrisDetector::drawOnImage(cv::Mat &image,
                 const std::vector<cv::Point> &points,
                 cv::Scalar color,
                 int radius, int thickness)
{
    std::vector<cv::Point>::const_iterator it = points.begin();
    // 对于所有角点，绘制白色圆圈
    while(it != points.end())
    {
        cv::circle(image, *it, radius, color, thickness);
        ++ it;
    }
}
 
 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83

最后，使用该类的步骤如下，直接修改main函数：

#include <QCoreApplication>
#include <opencv2/core/core.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <QDebug>
#include "harrisdetector.h"

int main(int argc, char *argv[])
{
    QCoreApplication a(argc, argv);

    cv::Mat image = cv::imread("c:/031.jpg", 0);
    // 创建Harris对象
    HarrisDetector harris;
    // 计算Harris的值
    harris.detect(image);
    // 检测Harris的角点
    std::vector<cv::Point>pts;
    harris.getCorners(pts, 0.1);
    // 绘制角点图
    harris.drawOnImage(image, pts);
    cv::namedWindow("Harris Corners Final");

    return a.exec();
}
 
 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

生成的图像：

这里为了改进特征点检测结果，添加了额外的非极大值抑制步骤，目的是移除彼此相邻的Harris角点。这就要求Harris角点不只需要得分高于给定阈值，它还必须是局部极大值。在检测中使用了一个技巧，即将Harris得分的图像进行膨胀：

cv::dilate(cornerStrength, dilate, cv::Mat());
 
 
1

这是由于膨胀运算替换每个像素值为相邻范围内的最大值，因此只有局部极大值的点才会保留原样，并通过以下函数进行测试：

cv::compare(cornerStrength, dilate, localMax, cv::CMP_EQ);
 
 
1

其中，localMax矩阵仅在局部极大值的位置为真，因此又可以在getCornerMap函数中用它来抑制所有非极大值的特征（基于cv::bitwise_and函数）。

三、引入适合跟踪的优质特征的Harris检测实现

在浮点处理器的帮助下，为了避免特征值分解而引入的数学上的简化变得微不足道，因此Harris检测可以基于计算而得的特征值。原则上这个修改不会显著影响检测的结果，但是能够避免使用任意的k参数。

以上第二中方法引入了局部极大值的条件，改善了部分效果。然而，特征点倾向于图像中不均匀分布、普遍集中在纹理丰富的部分。这里记录一种新的解决方案：

该方案利用两个特征点之间的最小距离，从Harris得分最高的点开始，仅接受离开有效特征点距离大于特定值的那些点。在OpenCV中提供cv::goodFeaturesToTrack实现这一方法，它检测到的特征能用于视觉跟踪应用中的优质特征集合。其调用方式如下：

    // 计算适合跟踪的优质特征
    std::vector<cv::Point2f> corners;
    cv::goodFeaturesToTrack(image,corners,
        250,    // 返回的最大特征点的数目
        0.01,   // 质量等级
        8); // 两点之间的最小允许距离
 
 
1
2
3
4
5
6

除了质量等级阈值、特征点之间的最小允许距离，该函数还需要指定返回的最大特征点数目，这是因为特征点是按照强度进行排序的。以下给出该方法的实现代码，直接在main函数中添加：

    // 输入图像
    cv::Mat image= cv::imread("c:/031.jpg",0);

    // 计算适合跟踪的优质特征
    std::vector<cv::Point2f> corners;
    cv::goodFeaturesToTrack(image,corners,
        250,    // 返回的最大特征点的数目
        0.01,   // 质量等级
        8); // 两点之间的最小允许距离

    // 遍历所有特征点并画圆圈
    std::vector<cv::Point2f>::const_iterator it= corners.begin();
    while (it!=corners.end())
    {
        cv::circle(image, *it, 3, cv::Scalar(255,255,255), 2);
        ++it;
    }

    // 显示输出结果
    cv::namedWindow("Good Features to Track");
    cv::imshow("Good Features to Track",image);

 
 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

返回生成的结果：

可以看到，该方法显著改进了特征点的分布情况，但是这样也增加了检测的复杂度，因为要求特征点要安装Harris的得分进行排序。该函数也可以指定一个可选的参数，使得按照经典的焦点分数定义进行计算。

其中，cv::goodFeaturesToTrack函数拥有一个封装类cv::GoodFeatureToTrackDetector，它继承自cv::FeaturesDetector类。其用法与以上的Harris类相类似：

    // 特征点向量
    std::vector<cv::KeyPoint> keypoints;
    cv::goodFeaturesToTrackDetector gftt(
        250,    // 返回的最大特征点的数目
        0.01,   // 质量等级
        8); // 两点之间的最小允许距离
    // 使用FeatureDetector的函数进行检测
    gftt.detect(image, keypoints);
 
 
1
2
3
4
5
6
7
8

结果与先前得到的结果是一样的，因为它们调用的是同一个函数。

本节的代码下载地址：http://download.csdn.net/detail/liyuefeilong/8483013