6.1 Harris 角点检测
6.1.1角点检测概述及原理
关于兴趣点（interest points）
在图像处理和与计算机视觉领域，兴趣点（interest points），或称作关键点（keypoints）、特征点（feature points） 被大量用于解决物体识别，图像识别、图像匹配、视觉跟踪、三维重建等一系列的问题。我们不再观察整幅图，而是选择某些特殊的点，然后对他们进行局部有的放矢的分析。如果能检测到足够多的这种点，同时他们的区分度很高，并且可以精确定位稳定的特征，那么这个方法就有使用价值。
图像特征类型可以被分为如下三种：
- 边缘；
- 角点 (感兴趣关键点)；
- 斑点(Blobs)(感兴趣区域)；
其中，角点是个很特殊的存在。他们在图像中可以轻易地定位，同时，他们在人造物体场景，比如门、窗、桌等出随处可见。因为角点位于两条边缘的交点处，代表了两个边缘变化的方向上的点,，所以他们是可以精确定位的二维特征，甚至可以达到亚像素的精度。且其图像梯度有很高的变化，这种变化是可以用来帮助检测角点的。需要注意的是，角点与位于相同强度区域上的点不同，与物体轮廓上的点也不同，因为轮廓点难以在相同的其他物体上精确定位。
角点检测算法的分类
在当前的图像处理领域，角点检测算法可归纳为三类：
- 基于灰度图像的角点检测；
- 基于二值图像的角点检测；
- 基于轮廓曲线的角点检测。
而基于灰度图像的角点检测又可分为基于梯度、基于模板和基于模板梯度组合三类方法，其中基于模板的方法主要考虑像素领域点的灰度变化，即图像亮度的变化，将与邻点亮度对比足够大的点定义为角点。常见的基于模板的角点检测算法有Kitchen-Rosenfeld角点检测算法，Harris角点检测算法、KLT角点检测算法及SUSAN角点检测算法。和其他角点检测算法相比，SUSAN角点检测算法具有算法简单、位置准确、抗噪声能力强等特点。
角点的定义
“如果某一点在任意方向的一个微小变动都会引起灰度很大的变化，那么我们就把它称之为角点”。
角点检测(Corner Detection)是计算机视觉系统中用来获得图像特征的一种方法，广泛应用于运动检测、图像匹配、视频跟踪、三维建模和目标识别等领域中。也称为特征点检测。
角点通常被定义为两条边的交点，更严格的说，角点的局部邻域应该具有两个不同区域的不同方向的边界。而实际应用中，大多数所谓的角点检测方法检测的是拥有特定特征的图像点，而不仅仅是“角点”。这些特征点在图像中有具体的坐标，并具有某些数学特征，如局部最大或最小灰度、某些梯度特征等。
现有的角点检测算法并不是都十分的健壮。很多方法都要求有大量的训练集和冗余数据来防止或减少错误特征的出现。另外，角点检测方法的一个很重要的评价标准是其对多幅图像中相同或相似特征的检测能力，并且能够应对光照变化、图像旋转等图像变化。
在我们解决问题时，往往希望找到特征点，“特征”顾名思义，指能描述物体本质的东西，还有一种解释就是这个特征微小的变化都会对物体的某一属性产生重大的影响。而角点就是这样的特征。
观察日常生活中的“角落”就会发现，“角落”可以视为所有平面的交汇处，或者说是所有表面的发起处。假设我们要改变一个墙角的位置，那么由它而出发的平面势必都要有很大的变化。所以，这就引出了图像角点的定义。
我们知道，特征检测与匹配是计算机视觉应用中非常重要的一部分，这需要寻找图像之间的特征建立对应关系。图像中的点作为图像的特殊位置，是很常用的一类特征，点的局部特征也可以叫做“关键特征点”(keypoint feature)，或“兴趣点”(interest point)，或“角点”(conrner)。
另外，关于角点的具体描述可以有几种：
- 一阶导数(即灰度的梯度)的局部最大所对应的像素点；
- 两条及两条以上边缘的交点；
- 图像中梯度值和梯度方向的变化速率都很高的点；
- 角点处的一阶导数最大，二阶导数为零，指示物体边缘变化不连续的方向。
角点检测的工作原理
由于角点代表了图像像素梯度变化，我们将寻找这个”变化”。考虑到一个灰度图像
I
I
。划动窗口 w(x,y)w(x,y) (with displacements
u
u
在 xx方向和
v
v
方向) 计算像素灰度变化。
其中:
w(x,y)w(x,y) is the window at position
(x,y)
(
x
,
y
)
。
I(x,y)
I
(
x
,
y
)
is the intensity at
(x,y)
(
x
,
y
)
。
i(x+u,y+v)
i
(
x
+
u
,
y
+
v
)
is the intensity at the moved window
(x+u,y+v)
(
x
+
u
,
y
+
v
)
。
为了寻找带角点的窗口，我们搜索像素灰度变化较大的窗口。于是, 我们期望最大化以下式子：
使用泰勒(Taylor)展开式：
式子可以展开为：
一个矩阵表达式可以写为:
表示为:
因此我们有等式:
每个窗口中计算得到一个值。这个值决定了这个窗口中是否包含了角点：
其中：
一个窗口，它的分数 大于一个特定值，这个窗口就可以被认为是”角点”
6.1.2 Harris 角点检测相关API
6.1.2.1 Harris角点检测
当一个窗口在图像上移动，在平滑区域如图1-(a)，窗口在各个方向上没有变化。在边缘上如图1-(b)，窗口在边缘的方向上没有变化。在角点处如图1-(c)，窗口在各个方向上具有变化。Harris角点检测正是利用了这个直观的物理现象，通过窗口在各个方向上的变化程度，决定是否为角点。
图1
将图像窗口平移
[u,v]
[
u
,
v
]
产生灰度变化
E(u,v)
E
(
u
,
v
)
由：
， 得到：
对于局部微小的移动量 [u,v]，近似表达为：
其中M是 2*2 矩阵，可由图像的导数求得：
E(u,v)的椭圆形式如下图：
图2
定义角点响应函数
R
R
为：
Harris角点检测算法就是对角点响应函数R进行阈值处理：R > threshold，即提取R的局部极大值。
参考论文：A Combined Corner and Edge Detector-1988. （见附件）
6.1.2.2 Harris 角点检测相关API函数讲解
角点检测：cornerHarris函数详解
cornerHarris 函数用于在OpenCV中运行Harris角点检测算子处理图像。和cornerMinEigenVal( )以及cornerEigenValsAndVecs( )函数类似，cornerHarris 函数对于每一个像素(x,y)(x,y)在邻域blockSize×blockSize内，计算2x2梯度的协方差矩阵
M((x,y))
M
(
(
x
,
y
)
)
，接着它计算如下式子：
即可以找出输出图中的局部最大值，即找出了角点。
C++: void cornerHarris( InputArray src,
OutputArray dst,
int blockSize,
int ksize,
double k,
intborderType=BORDER_DEFAULT ) 【参数】
第一个参数，InputArray类型的src，输入图像，即源图像，填Mat类的对象即可，且需为单通道8位或者浮点型图像。
第二个参数，OutputArray类型的dst，函数调用后的运算结果存在这里，即这个参数用于存放Harris角点检测的输出结果，和源图片有一样的尺寸和类型。
第三个参数，int类型的blockSize，表示邻域的大小，更多的详细信息在cornerEigenValsAndVecs（）中有讲到。
第四个参数，int类型的ksize，表示Sobel()算子的孔径大小。
第五个参数，double类型的k，Harris参数。
第六个参数，int类型的borderType，图像像素的边界模式，注意它有默认值BORDER_DEFAULT。更详细的解释，参考borderInterpolate( )函数。
接着我们一起过一遍稍后需要用到的Threshold函数的解析，然后看一个以cornerHarris为核心的示例程序。
阈值操作：Threshold()函数详解
函数Threshold( ) 对单通道数组应用固定阈值操作。该函数的典型应用是对灰度图像进行阈值操作得到二值图像。(另外，compare( )函数也可以达到此目的) 或者是去掉噪声，例如过滤很小或很大象素值的图像点。
C++: double threshold( InputArray src,
OutputArray dst,
double thresh,
double maxVal,
int thresholdType) 【参数】
第一个参数，InputArray类型的src，输入数组，填单通道 , 8或32位浮点类型的Mat即可。
第二个参数，OutputArray类型的dst，函数调用后的运算结果存在这里，即这个参数用于存放输出结果，且和第一个参数中的Mat变量有一样的尺寸和类型。
第三个参数，double类型的thresh，阈值的具体值。
第四个参数，double类型的maxval，当第五个参数阈值类型type取 CV_THRESH_BINARY 或CV_THRESH_BINARY_INV 阈值类型时的最大值.
第五个参数，int类型的type，阈值类型,。threshold( )函数支持的对图像取阈值的方法由其确定，具体用法如下。
而图形化的阈值描述如下图：
图3
6.1.2.3 Harris 角点检测相关API函数源代码
角点检测：cornerHarris函数源代码
/*【cornerHarris ( )源代码】***********************************************************
* @Version：OpenCV 3.0.0（Opnencv2和Opnencv3差别不大，Linux和PC的对应版本源码完全一样，均在对应的安装目录下）
* @源码路径：…\opencv\sources\modules\imgproc\src\ corner.cpp
* @起始行数：597行
********************************************************************************/
void cv::cornerHarris( InputArray _src, OutputArray _dst, int blockSize, int ksize, double k, int borderType )
{
CV_OCL_RUN(_src.dims() <= 2 && _dst.isUMat(),
ocl_cornerMinEigenValVecs(_src, _dst, blockSize, ksize, k, borderType, HARRIS))
Mat src = _src.getMat();
_dst.create( src.size(), CV_32FC1 );
Mat dst = _dst.getMat();
#if IPP_VERSION_X100 >= 801 && 0
CV_IPP_CHECK()
{
int type = src.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);
int borderTypeNI = borderType & ~BORDER_ISOLATED;
bool isolated = (borderType & BORDER_ISOLATED) != 0;
if ( (ksize == 3 || ksize == 5) && (type == CV_8UC1 || type == CV_32FC1) &&
(borderTypeNI == BORDER_CONSTANT || borderTypeNI == BORDER_REPLICATE) && cn == 1 && (!src.isSubmatrix() || isolated) )
{
IppiSize roisize = { src.cols, src.rows };
IppiMaskSize masksize = ksize == 5 ? ippMskSize5x5 : ippMskSize3x3;
IppDataType datatype = type == CV_8UC1 ? ipp8u : ipp32f;
Ipp32s bufsize = 0;
double scale = (double)(1 << ((ksize > 0 ? ksize : 3) - 1)) * blockSize;
if (ksize < 0)
scale *= 2.0;
if (depth == CV_8U)
scale *= 255.0;
scale = std::pow(scale, -4.0);
if (ippiHarrisCornerGetBufferSize(roisize, masksize, blockSize, datatype, cn, &bufsize) >= 0)
{
Ipp8u * buffer = ippsMalloc_8u(bufsize);
IppiDifferentialKernel filterType = ksize > 0 ? ippFilterSobel : ippFilterScharr;
IppiBorderType borderTypeIpp = borderTypeNI == BORDER_CONSTANT ? ippBorderConst : ippBorderRepl;
IppStatus status = (IppStatus)-1;
if (depth == CV_8U)
status = ippiHarrisCorner_8u32f_C1R((const Ipp8u *)src.data, (int)src.step, (Ipp32f *)dst.data, (int)dst.step, roisize,
filterType, masksize, blockSize, (Ipp32f)k, (Ipp32f)scale, borderTypeIpp, 0, buffer);
else if (depth == CV_32F)
status = ippiHarrisCorner_32f_C1R((const Ipp32f *)src.data, (int)src.step, (Ipp32f *)dst.data, (int)dst.step, roisize,
filterType, masksize, blockSize, (Ipp32f)k, (Ipp32f)scale, borderTypeIpp, 0, buffer);
ippsFree(buffer);
if (status >= 0)
{
CV_IMPL_ADD(CV_IMPL_IPP);
return;
}
}
setIppErrorStatus();
}
}
#endif
cornerEigenValsVecs( src, dst, blockSize, ksize, HARRIS, k, borderType );
} 阈值操作：Threshold()函数详解
/*【threshold ( )函数源代码】*********************************************************
* @Version：OpenCV 3.0.0（Opnencv2和Opnencv3差别不大，Linux和PC的对应版本源码完全一样，均在对应的安装目录下）
* @源码路径：…\opencv\sources\modules\imgproc\src\ thresh.cpp
* @起始行数：1186行
********************************************************************************/
double cv::threshold( InputArray _src, OutputArray _dst, double thresh, double maxval, int type )
{
CV_OCL_RUN_(_src.dims() <= 2 && _dst.isUMat(),
ocl_threshold(_src, _dst, thresh, maxval, type), thresh)
Mat src = _src.getMat();
int automatic_thresh = (type & ~CV_THRESH_MASK);
type &= THRESH_MASK;
CV_Assert( automatic_thresh != (CV_THRESH_OTSU | CV_THRESH_TRIANGLE) );
if( automatic_thresh == CV_THRESH_OTSU )
{
CV_Assert( src.type() == CV_8UC1 );
thresh = getThreshVal_Otsu_8u( src );
}
else if( automatic_thresh == CV_THRESH_TRIANGLE )
{
CV_Assert( src.type() == CV_8UC1 );
thresh = getThreshVal_Triangle_8u( src );
}
_dst.create( src.size(), src.type() );
Mat dst = _dst.getMat();
if( src.depth() == CV_8U )
{
int ithresh = cvFloor(thresh);
thresh = ithresh;
int imaxval = cvRound(maxval);
if( type == THRESH_TRUNC )
imaxval = ithresh;
imaxval = saturate_cast<uchar>(imaxval);
if( ithresh < 0 || ithresh >= 255 )
{
if( type == THRESH_BINARY || type == THRESH_BINARY_INV ||
((type == THRESH_TRUNC || type == THRESH_TOZERO_INV) && ithresh < 0) ||
(type == THRESH_TOZERO && ithresh >= 255) )
{
int v = type == THRESH_BINARY ? (ithresh >= 255 ? 0 : imaxval) :
type == THRESH_BINARY_INV ? (ithresh >= 255 ? imaxval : 0) :
/*type == THRESH_TRUNC ? imaxval :*/ 0;
dst.setTo(v);
}
else
src.copyTo(dst);
return thresh;
}
thresh = ithresh;
maxval = imaxval;
}
else if( src.depth() == CV_16S )
{
int ithresh = cvFloor(thresh);
thresh = ithresh;
int imaxval = cvRound(maxval);
if( type == THRESH_TRUNC )
imaxval = ithresh;
imaxval = saturate_cast<short>(imaxval);
if( ithresh < SHRT_MIN || ithresh >= SHRT_MAX )
{
if( type == THRESH_BINARY || type == THRESH_BINARY_INV ||
((type == THRESH_TRUNC || type == THRESH_TOZERO_INV) && ithresh < SHRT_MIN) ||
(type == THRESH_TOZERO && ithresh >= SHRT_MAX) )
{
int v = type == THRESH_BINARY ? (ithresh >= SHRT_MAX ? 0 : imaxval) :
type == THRESH_BINARY_INV ? (ithresh >= SHRT_MAX ? imaxval : 0) :
/*type == THRESH_TRUNC ? imaxval :*/ 0;
dst.setTo(v);
}
else
src.copyTo(dst);
return thresh;
}
thresh = ithresh;
maxval = imaxval;
}
else if( src.depth() == CV_32F )
;
else
CV_Error( CV_StsUnsupportedFormat, "" );
parallel_for_(Range(0, dst.rows),
ThresholdRunner(src, dst, thresh, maxval, type),
dst.total()/(double)(1<<16));
return thresh;
} 好了，让我们看一个调用示例程序，代码参看附件【demo1】。
图4
6.1.3 Harris 角点检测综合实例
本次综合示例为调节滚动条来控制阈值，以控制的harris检测角点的数量。一共有三个图片窗口，分别为显示原始图的窗口，包含滚动条的彩色效果图窗口，以及灰度图效果图窗口。代码参看附件【demo2】。
图5
图6