Harris角点检测学习

时间：2017-05-24 10:06:46 阅读：281 评论：0 收藏：0 [点我收藏+]

标签：img 取值 names core name 跟踪计算 stream 不同

1.角点的定义与性质

角点是一种局部特征，具有旋转不变性和不随光照条件变化而变化的特点，一般将图像中曲率足够高或者曲率变化明显的点作为角点。检测得到的角点特征通常用于图像匹配、目标跟踪、运动估计等方面。

2.Harris角点

1）定性描述

该算法中，将图像分为平坦区域、边缘、角点三部分。平坦区域中像素灰度在各个方向上变换都很小，边缘上的像素灰度在某个方向变化很大，但是在另一些方向变化很小；位于角点上的像素灰度则在各个方向上的变化都比较大。这是通过人眼观察得到的直观感受。

2）量化的数学表达

不同区域像素的变化与相似性是相反的关系：变化越大，相似性越小，反之亦然。所以我们可以使用相似性度量其变化程度。图像I(x,y)在中心(x,y)处区域W(x,y)与相对中心(x,y)移动了△x/△y之后的区域相关性可以通过自相关函数给出

w(u,v)是加权函数，可选取常数或者高斯加权函数（一般选取后者）

看到delta标志就会想到泰勒展开式简化，简化流程如下：

技术分享

所以之前的相似性度量可以近似表达为一个二次项函数

技术分享

其本质是一个椭圆函数，椭圆的主方向就是图像变化率最大的主要方向。而椭圆的大小又是由M(x,y)的特征值决定的。

技术分享

套用这套理论到平坦区域、边缘、角点上，则表示为：

（1）边缘；椭圆狭长，一个特征值大、另一个特征值小

（2）平坦区域；两个特征值都很小

（3）角点；两个特征值都很大

技术分享

所以求出矩阵M(x,y)的两个特征值可以用于检测角点，但是计算量太大。Harris角点检测中避免直接求取特征值，而是利用M(x,y)的性质来计算两个特征值的相对大小

技术分享

只有两个特征值都比较大的情况下R取值才比较大

3.算法流程

根据前面的描述，只要计算出图像中各个像素对应的R，并根据一个设定的阈值二值化处理即可得到角点。为了滤除局部区域的非极大值，可以使用非极大值抑制法做进一步的处理。算法整理流程如下：

技术分享

4.opencv实现

个人使用的开发环境是是opencv_2.4.13+vs2012，现将代码贴出

 1 #include <iostream>
 2 #include <core/core.hpp>
 3 #include <highgui/highgui.hpp>
 4 #include <features2d/features2d.hpp>
 5 #include <imgproc/imgproc.hpp>
 6 
 7 using namespace std;
 8 using namespace cv;
 9 
10 int main(int argc, char* argv[])
11 {
12     /* 1.以bmp格式读取图片，注意该图片此时应该放置于.sln所在目录 */
13     Mat image = imread("../church01.jpg", 0);
14     if(!image.data)
15         return 0;
16 
17     namedWindow("originalImage");
18     imshow("originalImage", image);
19     
20     /* 2.计算图像沿x/y方向一阶导数 */
21     Mat xKernel = (Mat_<double>(1,3) << -1, 0, 1);
22     Mat yKernel = xKernel.t();
23 
24     Mat Ix, Iy;
25     filter2D(image, Ix, CV_64F, xKernel);
26     filter2D(image, Iy, CV_64F, yKernel);
27     
28     /* 3.计算矩阵M中的各项 */
29     Mat Ix2, Iy2, Ixy;
30     Ix2 = Ix.mul(Ix);
31     Iy2 = Iy.mul(Iy);
32     Ixy = Ix.mul(Iy);
33 
34     /* 注意这里使用的是一阶高斯滤波而非二阶，尚未
35      * 理解其原因
36      */
37     Mat gaussKernel = getGaussianKernel(7,1);
38     filter2D(Ix2, Ix2, CV_64F, gaussKernel);
39     filter2D(Iy2, Iy2, CV_64F, gaussKernel);
40     filter2D(Ixy, Ixy, CV_64F, gaussKernel);
41 
42     /* 4.根据公式计算R矩阵 */
43     Mat cornerStrength(image.size(), CV_64F);
44     float alpha = 0.1f;
45     int rows, cols;
46     rows = image.rows;
47     cols = image.cols;
48     for(int i = 0; i < rows; i++)
49     {
50         for(int j = 0; j < cols; j++)
51         {
52             double det_m = Ix2.at<double>(i,j) * Iy2.at<double>(i,j) - Ixy.at<double>(i,j) * Ixy.at<double>(i,j);
53             double trace_m = Ix2.at<double>(i,j) + Iy2.at<double>(i,j);
54             cornerStrength.at<double>(i,j) = det_m - alpha*trace_m*trace_m;
55         }
56     }
57 
58     /* 5.阈值化并进行非极大值抑制 */
59     double maxStrength, minStrength;
60     minMaxLoc(cornerStrength, &minStrength, &maxStrength);
61     Mat dilated;
62     Mat localMax;
63     dilate(cornerStrength, dilated, Mat());
64     compare(cornerStrength, dilated, localMax, CMP_EQ);
65 
66     double threshold = 0.01 * maxStrength;
67     Mat cornerMap;
68     cornerMap = cornerStrength > threshold;
69     bitwise_and(cornerMap, localMax, cornerMap);
70     namedWindow("cornerMap");
71     imshow("cornerMap", cornerMap);
72 
73     /* 6.在原图上绘制角点 */
74     image = imread("../church01.jpg", 0);
75     for(int x = 0;x < cornerMap.cols; x++)
76     {
77         for(int y = 0; y < cornerMap.rows; y++)
78         {
79             if(cornerMap.at<uchar>(y,x))
80             {
81                 circle(image, Point(x,y), 3, Scalar(255), 1);
82 
83             }
84         }
85     }
86 
87     namedWindow("Harris Corner");
88     imshow("Harris Corner", image);
89     waitKey();
90     return 0;
91 }