本文的内容基本来源于以下博客,主要是为了方便以后再次学习。
http://blog.csdn.net/newthinker_wei/article/details/45603583
http://www.360doc.com/content/15/1212/23/20007814_519967668.shtml
https://github.com/RonnyYoung/ImageFeatures/blob/master/source/harris.cpp
Harris原理分析:
人眼对角点的识别通常是在一个局部的小区域或小窗口完成的。如果这个特定的小窗口在各个方向上移动时,窗口内图像的灰度发生了较大的变化,那么就认为在窗口内遇到了角点。如果这个特定的窗口在图像各个方向上移动时,窗口内图像的灰度没有发生变化,那么窗口内就不存在角点;如果窗口在某一个方向移动时,窗口内图像的灰度发生了较大的变化,而在另一些方向上没有发生变化,那么,窗口内的图像可能就是一条线段。如下图:
首先,对于图像,在点(x, y)处平移(u, v)后产生的灰度变化的自相关函数如下:
其中,窗口函数可以是平坦的,也可以是高斯的,如下图:
将 I(x+u, y+v) 进行泰勒展开如下:
则有
由于是局部微小的移动量(u, v),所以可以近似得到下面忽略余项之后的表达式,是一个二项式函数:
所以,
其中,M的表达式如下,可由图像的导数求得:
上面说到,忽略余项之后的表达式为一个二项式函数。本质上,二项式函数就是一个椭圆函数,椭圆的扁率和尺寸是由M(x,y)的特征值λ1、λ2决定的,椭圆的方向是由M(x,y)的特征矢量决定的,如下图所示,椭圆方程为:
椭圆函数的特征值与图像中的角点、直线(边缘)和平面之间的关系如下图所示, 可分为三种情况:
a. 图像中的直线:一个特征值大,另一个特征值小,λ1>λ2或λ2>λ1;自相关函数值在某一方向上大,在其他方向上小。
b. 图像中的平面:两个特征值都小,且近似相等;自相关函数数值在各个方向上都小。
c. 图像中的角点:两个特征值都大,且近似相等;自相关函数在所有方向都增大。
我们是通过M的两个特征值的大小对图像进行分类,所以,定义角点相应函数R:
其中k为经验常数,一般是一个远小于1的数, 在opencv中的默认值是0.04。
对于边缘部分,λ1>>λ2且 λ2< k*λ1,角点相应函数R: λ1λ2 - k(λ1+λ2)^2 ≈ λ1λ2 - k(λ1)^2 < λ1 k λ1 - k(λ1)^2 = 0,因此,在边缘处的角点响应值小于0。
对于非边缘部分,λ1、λ2相差不大,可以认为(λ1+λ2)^2 ≈ 4λ1λ2,角点相应函数R: λ1λ2 - k(λ1+λ2)^2 ≈ λ1λ2- 4kλ1λ2 = (1- 4k)λ1λ2,由于k远小于1,因此, (1- 4k) ≈ 1,所以λ1λ2 - k(λ1+λ2)^2 ≈ λ1λ2
在角点处,λ1、λ2都比较大,故角点响应值也比较大;在平坦区域λ1、λ2都很小,故角点响应值也很小。
所以,上图可以转化为:
其中:
• R 只与M的特征值有关
• 角点:R 为大数值正数
• 边缘:R 为大数值负数
• 平坦区:R 为小数值
在判断角点的时候,对角点响应函数R进行阈值处理:R > threshold,并且提取R的局部极大值。
Harris角点的缺点:
由Harris角点的定义可知,该角点具有旋转不变性。
Harris角点不具有尺度不变性。
Harris代码实现:
void detectHarrisCorners(const Mat& imgSrc, Mat& imgDst, double alpha)
{
Mat gray;
if (imgSrc.channels() == 3)
{
cvtColor(imgSrc, gray, CV_BGR2GRAY);
}
else
{
gray = imgSrc.clone();
}
gray.convertTo(gray, CV_64F);
Mat xKernel = (Mat_<double>(1, 3) << -1, 0, 1);
Mat yKernel = xKernel.t();
Mat Ix, Iy;
filter2D(gray, Ix, CV_64F, xKernel);
filter2D(gray, Iy, CV_64F, yKernel);
Mat Ix2, Iy2, Ixy;
Ix2 = Ix.mul(Ix);
Iy2 = Iy.mul(Iy);
Ixy = Ix.mul(Iy);
Mat gaussKernel = getGaussianKernel(7, 1);
filter2D(Ix2, Ix2, CV_64F, gaussKernel);
filter2D(Iy2, Iy2, CV_64F, gaussKernel);
filter2D(Ixy, Ixy, CV_64F, gaussKernel);
// 计算每个像素点的角点响应函数
Mat cornerStrength(gray.size(), gray.type());
for (int i = 0; i < gray.rows; i++)
{
for (int j = 0; j < gray.cols; j++)
{
double det_m = Ix2.at<double>(i, j) * Iy2.at<double>(i, j) - Ixy.at<double>(i, j) * Ixy.at<double>(i, j);
double trace_m = Ix2.at<double>(i, j) + Iy2.at<double>(i, j);
cornerStrength.at<double>(i, j) = det_m - alpha * trace_m *trace_m;
}
}
double maxStrength;
minMaxLoc(cornerStrength, NULL, &maxStrength, NULL, NULL);
Mat dilated;
Mat localMax;
// 默认采用3*3核进行膨胀,膨胀之后,除了局部最大值点和原来相同,其它非局部最大值点被3*3邻域内的最大值取代
dilate(cornerStrength, dilated, Mat());
// 与原图相比,只剩下和原图值相同的点,这些点都是局部最大值点,保存到localMax
compare(cornerStrength, dilated, localMax, CMP_EQ);
Mat cornerMap;
double qualityLevel = 0.01;
// 根据角点响应最大值计算阈值
double thresh = qualityLevel * maxStrength;
// 根据阈值得到角点图
cornerMap = cornerStrength > thresh;
// 和局部最大值图进行与操作,完成非最大值抑制
bitwise_and(cornerMap, localMax, cornerMap);
imgDst = cornerMap.clone();
}
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);
}
}