opencv-角点检测之Harris角点检测

Posted 2020-12-31 tang-tangt

 

转自：https://blog.csdn.net/poem_qianmo/article/details/29356187

 

先看看程序运行截图：

 

 

一、引言：关于兴趣点（interest points）

 

在图像处理和与计算机视觉领域，兴趣点（interest points），或称作关键点（keypoints）、特征点（feature points） 被大量用于解决物体识别，图像识别、图像匹配、视觉跟踪、三维重建等一系列的问题。我们不再观察整幅图，而是选择某些特殊的点，然后对他们进行局部有的放矢的分析。如果能检测到足够多的这种点，同时他们的区分度很高，并且可以精确定位稳定的特征，那么这个方法就有使用价值。

 

图像特征类型可以被分为如下三种：

• <1>边缘
• <2>角点 (感兴趣关键点)
• <3>斑点(Blobs)(感兴趣区域)

 

其中，角点是个很特殊的存在。他们在图像中可以轻易地定位，同时，他们在人造物体场景，比如门、窗、桌等出随处可见。因为角点位于两条边缘的交点处，代表了两个边缘变化的方向上的点,，所以他们是可以精确定位的二维特征，甚至可以达到亚像素的精度。且其图像梯度有很高的变化，这种变化是可以用来帮助检测角点的。需要注意的是，角点与位于相同强度区域上的点不同，与物体轮廓上的点也不同，因为轮廓点难以在相同的其他物体上精确定位。

 

二、角点检测算法的分类

在当前的图像处理领域，角点检测算法可归纳为三类：

• <1>基于灰度图像的角点检测
• <2>基于二值图像的角点检测
• <3>基于轮廓曲线的角点检测

而基于灰度图像的角点检测又可分为基于梯度、基于模板和基于模板梯度组合三类方法，其中基于模板的方法主要考虑像素领域点的灰度变化，即图像亮度的变化，将与邻点亮度对比足够大的点定义为角点。常见的基于模板的角点检测算法有Kitchen-Rosenfeld角点检测算法，Harris角点检测算法、KLT角点检测算法及SUSAN角点检测算法。和其他角点检测算法相比，SUSAN角点检测算法具有算法简单、位置准确、抗噪声能力强等特点。

 

三、角点的定义

“如果某一点在任意方向的一个微小变动都会引起灰度很大的变化，那么我们就把它称之为角点”

角点检测(Corner Detection)是计算机视觉系统中用来获得图像特征的一种方法，广泛应用于运动检测、图像匹配、视频跟踪、三维建模和目标识别等领域中。也称为特征点检测。

角点通常被定义为两条边的交点，更严格的说，角点的局部邻域应该具有两个不同区域的不同方向的边界。而实际应用中，大多数所谓的角点检测方法检测的是拥有特定特征的图像点，而不仅仅是“角点”。这些特征点在图像中有具体的坐标，并具有某些数学特征，如局部最大或最小灰度、某些梯度特征等。

现有的角点检测算法并不是都十分的健壮。很多方法都要求有大量的训练集和冗余数据来防止或减少错误特征的出现。另外，角点检测方法的一个很重要的评价标准是其对多幅图像中相同或相似特征的检测能力，并且能够应对光照变化、图像旋转等图像变化。

在我们解决问题时，往往希望找到特征点，“特征”顾名思义，指能描述物体本质的东西，还有一种解释就是这个特征微小的变化都会对物体的某一属性产生重大的影响。而角点就是这样的特征。

观察日常生活中的“角落”就会发现，“角落”可以视为所有平面的交汇处，或者说是所有表面的发起处。假设我们要改变一个墙角的位置，那么由它而出发的平面势必都要有很大的变化。所以，这就引出了图像角点的定义。

 

我们知道，特征检测与匹配是计算机视觉应用中非常重要的一部分，这需要寻找图像之间的特征建立对应关系。图像中的点作为图像的特殊位置，是很常用的一类特征，点的局部特征也可以叫做“关键特征点”(keypoint feature)，或“兴趣点”(interest point)，或“角点”(conrner)。

 

另外，关于角点的具体描述可以有几种：

• 一阶导数(即灰度的梯度)的局部最大所对应的像素点；
• 两条及两条以上边缘的交点；
• 图像中梯度值和梯度方向的变化速率都很高的点；
• 角点处的一阶导数最大，二阶导数为零，指示物体边缘变化不连续的方向。

 

四、cornerHarris函数详解

 

cornerHarris 函数用于在OpenCV中运行Harris角点检测算子处理图像。和cornerMinEigenVal( )以及cornerEigenValsAndVecs( )函数类似，cornerHarris 函数对于每一个像素(x,y)在blockSize × blockSize 邻域内，计算2x2梯度的协方差矩阵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 type)

• 第一个参数，InputArray类型的src，输入数组，填单通道 , 8或32位浮点类型的Mat即可。
• 第二个参数，OutputArray类型的dst，函数调用后的运算结果存在这里，即这个参数用于存放输出结果，且和第一个参数中的Mat变量有一样的尺寸和类型。
• 第三个参数，double类型的thresh，阈值的具体值。
• 第四个参数，double类型的maxval，当第五个参数阈值类型type取 CV_THRESH_BINARY 或CV_THRESH_BINARY_INV 阈值类型时的最大值.
• 第五个参数，int类型的type，阈值类型,。threshold( )函数支持的对图像取阈值的方法由其确定，具体用法如下图：

 

 

 

 

而图形化的阈值描述如下图：

 

  

讲解完这两个函数，让我们看一个调用示例程序：

//-----------------------------------【头文件包含部分】---------------------------------------
 
// 描述：包含程序所依赖的头文件
 
//----------------------------------------------------------------------------------------------
 
#include <opencv2/opencv.hpp>
 
#include <opencv2/imgproc/imgproc.hpp>
 
//-----------------------------------【命名空间声明部分】---------------------------------------
 
// 描述：包含程序所使用的命名空间
 
//-----------------------------------------------------------------------------------------------
 
using namespace cv;
 
 
int main()
 
{
 
//以灰度模式载入图像并显示
 
Mat srcImage = imread("1.jpg", 0);
 
imshow("原始图", srcImage);
 
 
 
//进行Harris角点检测找出角点
 
Mat cornerStrength;
 
cornerHarris(srcImage, cornerStrength, 2, 3, 0.01);
 
 
 
//对灰度图进行阈值操作，得到二值图并显示
 
Mat harrisCorner;
 
threshold(cornerStrength, harrisCorner, 0.00001, 255, THRESH_BINARY);
 
imshow("角点检测后的二值效果图", harrisCorner);
 
 
waitKey(0);
 
return 0;
 
}

 

运行截图：

   

 

 

 

六、本文相关核心函数在OpenCV中的实现源代码

 

这个部分贴出OpenCV中本文相关函数的源码实现细节，来给想了解实现细节的小伙伴们参考。浅墨暂时不在源码的细节上挖深作详细注释。

 

6.1 OpenCV2.X中cornerHarris函数源代码

源码路径： …opencvsourcesmodulesimgprocsrccorner.cpp

 

void cv::cornerHarris( InputArray _src,OutputArray _dst, int blockSize, int ksize, double k, int borderType )
 
{
 
Mat src = _src.getMat();
 
_dst.create( src.size(), CV_32F );
 
Mat dst = _dst.getMat();
 
cornerEigenValsVecs( src, dst, blockSize, ksize, HARRIS, k, borderType);
 
}
 

 

可见cornerHarris内部其实是调用了cornerEigenValsVecs函数，我们看看其实现源码：

 

 
static void
 
cornerEigenValsVecs( const Mat& src,Mat& eigenv, int block_size,
 
int aperture_size, intop_type, double k=0.,
 
intborderType=BORDER_DEFAULT )
 
{
 
#ifdef HAVE_TEGRA_OPTIMIZATION
 
if (tegra::cornerEigenValsVecs(src, eigenv, block_size, aperture_size,op_type, k, borderType))
 
return;
 
#endif
 
  
int depth = src.depth();
 
double scale = (double)(1 << ((aperture_size > 0 ?aperture_size : 3) - 1)) * block_size;
 
if( aperture_size < 0 )
 
scale *= 2.;
 
if( depth == CV_8U )
 
scale *= 255.;
 
scale = 1./scale;
 
  
CV_Assert( src.type() == CV_8UC1 || src.type() == CV_32FC1 );
 
 
Mat Dx, Dy;
 
if( aperture_size > 0 )
 
{
 
Sobel( src, Dx, CV_32F, 1, 0, aperture_size, scale, 0, borderType );
 
Sobel( src, Dy, CV_32F, 0, 1, aperture_size, scale, 0, borderType );
 
}
 
else
 
{
 
Scharr( src, Dx, CV_32F, 1, 0, scale, 0, borderType );
 
Scharr( src, Dy, CV_32F, 0, 1, scale, 0, borderType );
 
}
 
  
Size size = src.size();
 
Mat cov( size, CV_32FC3 );
 
int i, j;
 
 
for( i = 0; i < size.height; i++ )
 
{
 
float* cov_data = (float*)(cov.data + i*cov.step);
 
const float* dxdata = (const float*)(Dx.data + i*Dx.step);
 
const float* dydata = (const float*)(Dy.data + i*Dy.step);
 
 
for( j = 0; j < size.width; j++ )
 
{
 
float dx = dxdata[j];
 
float dy = dydata[j];
 
 
 
cov_data[j*3] = dx*dx;
 
cov_data[j*3+1] = dx*dy;
 
cov_data[j*3+2] = dy*dy;
 
}
 
}
 
 
boxFilter(cov, cov, cov.depth(), Size(block_size, block_size),
 
Point(-1,-1), false, borderType );
 
 
if( op_type == MINEIGENVAL )
 
calcMinEigenVal( cov, eigenv );
 
else if( op_type == HARRIS )
 
calcHarris( cov, eigenv, k );
 
else if( op_type == EIGENVALSVECS )
 
calcEigenValsVecs( cov, eigenv );
 
} 
 
}

 

6.1 OpenCV2.X中Threshold函数源代码

 

路径：…opencvsourcesmodulesimgprocsrc hresh.cpp

 

double cv::threshold( InputArray _src,OutputArray _dst, double thresh, double maxval, int type )
 
{
 
Mat src = _src.getMat();
 
bool use_otsu = (type & THRESH_OTSU) != 0;
 
type &= THRESH_MASK;
 
 
 
if( use_otsu )
 
{
 
CV_Assert( src.type() == CV_8UC1 );
 
thresh = getThreshVal_Otsu_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;
 
}

 

 

另外在贴上与之相关的自适应阈值操作函数的源码adaptiveThreshold：

 

void cv::adaptiveThreshold( InputArray_src, OutputArray _dst, double maxValue,
 
int method, inttype, int blockSize, double delta )
 
{
 
Mat src = _src.getMat();
 
CV_Assert( src.type() == CV_8UC1 );
 
CV_Assert( blockSize % 2 == 1 && blockSize > 1 );
 
Size size = src.size();
 
 
_dst.create( size, src.type() );
 
Mat dst = _dst.getMat();
 
 
 
if( maxValue < 0 )
 
{
 
dst = Scalar(0);
 
return;
 
}
 
 
Mat mean;
 
 
if( src.data != dst.data )
 
mean = dst;
 
 
if( method == ADAPTIVE_THRESH_MEAN_C )
 
boxFilter( src, mean, src.type(), Size(blockSize, blockSize),
 
Point(-1,-1), true,BORDER_REPLICATE );
 
else if( method == ADAPTIVE_THRESH_GAUSSIAN_C )
 
GaussianBlur( src, mean, Size(blockSize, blockSize), 0, 0,BORDER_REPLICATE );
 
else
 
CV_Error( CV_StsBadFlag, "Unknown/unsupported adaptive thresholdmethod" );
 
int i, j;
 
uchar imaxval = saturate_cast<uchar>(maxValue);
 
int idelta = type == THRESH_BINARY ? cvCeil(delta) : cvFloor(delta);
 
uchar tab[768];
 
 
if( type == CV_THRESH_BINARY )
 
for( i = 0; i < 768; i++ )
 
tab[i] = (uchar)(i - 255 > -idelta ? imaxval : 0);
 
else if( type == CV_THRESH_BINARY_INV )
 
for( i = 0; i < 768; i++ )
 
tab[i] = (uchar)(i - 255 <= -idelta ? imaxval : 0);
 
else
 
CV_Error( CV_StsBadFlag, "Unknown/unsupported threshold type");
 
if( src.isContinuous() && mean.isContinuous() &&dst.isContinuous() )
 
{
 
size.width *= size.height;
 
size.height = 1;
 
}
 
for( i = 0; i < size.height; i++ )
 
{
 
const uchar* sdata = src.data + src.step*i;
 
const uchar* mdata = mean.data + mean.step*i;
 
uchar* ddata = dst.data + dst.step*i; 
 
for( j = 0; j < size.width; j++ )
 
ddata[j] = tab[sdata[j] - mdata[j] + 255];
 
}
 
}