OpenCV中常用的角点检测为Harris角点和ShiTomasi角点。

以OpenCV源代码文件 .\opencv\sources\samples\cpp\tutorial_code\TrackingMotion\cornerDetector_Demo.cpp为例,主要分析其中的这两种角点检测源代码。角点检测数学原理请参考我之前转载的一篇博客 http://www.cnblogs.com/riddick/p/7645904.html,分析的很详细,不再赘述。本文主要分析其源代码:

1. Harris角点检测  

  根据数学上的推导,可以根据图像中某一像素点邻域内构建的协方差矩阵获取特征值和特征向量,根据特征值建立特征表达式,如下:                

(αβ) - k(α+β)^2

  可以根据上式的值得大小来判断该像素点是平坦区域内点、边界点还是角点。下面说一下怎么在原图像中建立协方差矩阵并求取特征值α和β和特征向量t1, t2。

  该例程代码中调用cornerEigenValsAndVecs()函数计算特征值和特征向量。函数原型如下:

void cv::cornerEigenValsAndVecs( InputArray _src, OutputArray _dst, int blockSize, int ksize, int borderType )

  src为输入灰度图像,dst为输出(6通道 CV_32FC(6),依次保存的是α, t1,  β, t2),blockSize为邻域大小,ksize为sobel求取微分时的窗口大小。  

  该函数内部调用cornerEigenValsVecs()函数,原型如下: 

static void  cornerEigenValsVecs( const Mat& src, Mat& eigenv, int block_size,int aperture_size, int op_type, double k=0.,int borderType=BORDER_DEFAULT )

  主要介绍一下op_type这个参数,该参数是一个枚举值,有三个值可以选择(MINEIGENVAL, HARRIS, EIGENVALSVECS)

  ①MINEIGENVAL用于ShiTomasi角点检测中获取两个特征值中较小的那个值,用以获取强角点,随后介绍;

  ②HARRIS在cornerHarris()函数中用到,用于直接利用协方差矩阵获取特征表达式值的大小,k值在此时会被设置,通常为0.04,其他情况下设置为0;

  ③EIGENVALSVECS就是本例程中设置的,求取两个特征值和特征向量。

  在cornerEigenValsVecs()函数中,先利用sobel算子求水平方向和竖直方向的微分,窗口大小为前述,如下代码:         

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 );} 

  然后初始化协方差矩阵cov(三通道,依次保存dx*dx, dx*dy, dy*dy),如下:

for( ; 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 );

  然后就是利用滤波后的协方差矩阵求取特征值和特征向量了,根据设定不同的op_type调用不同的函数计算,本例程中为调用最后一个calcEigenValsVecs()函数,该函数如下:

static void calcEigenValsVecs( const Mat& _cov, Mat& _dst )
{Size size = _cov.size();if( _cov.isContinuous() && _dst.isContinuous() ){size.width *= size.height;size.height = 1;}for( int i = 0; i < size.height; i++ ){const float* cov = _cov.ptr<float>(i);float* dst = _dst.ptr<float>(i);//调用该函数计算2x2协方差矩阵的特征值和特征向量eigen2x2(cov, dst, size.width);}
}

  该函数中调用eigen2x2()函数计算每个像素点处协方差矩阵的2个特征值和2个特征向量,协方差矩阵为如下形式,数据都保存在cov的三个通道中:

   eigen2x2()函数如下:2x2矩阵特征值和特征向量的计算,有线性代数基础的都学过,就不再赘述

static void eigen2x2( const float* cov, float* dst, int n )
{for( int j = 0; j < n; j++ ){double a = cov[j*3];double b = cov[j*3+1];double c = cov[j*3+2];double u = (a + c)*0.5;double v = std::sqrt((a - c)*(a - c)*0.25 + b*b);          //计算两个特征值l1,l2double l1 = u + v;double l2 = u - v;//计算特征值l1对应的特征向量double x = b;double y = l1 - a;double e = fabs(x);if( e + fabs(y) < 1e-4 ){y = b;x = l1 - c;e = fabs(x);if( e + fabs(y) < 1e-4 ){e = 1./(e + fabs(y) + FLT_EPSILON);x *= e, y *= e;}}double d = 1./std::sqrt(x*x + y*y + DBL_EPSILON);//保存特征值l1及其对应的特征向量dst[6*j] = (float)l1;dst[6*j + 2] = (float)(x*d);dst[6*j + 3] = (float)(y*d); //计算特征值l2对应的特征向量x = b;y = l2 - a;e = fabs(x);if( e + fabs(y) < 1e-4 ){y = b;x = l2 - c;e = fabs(x);if( e + fabs(y) < 1e-4 ){e = 1./(e + fabs(y) + FLT_EPSILON);x *= e, y *= e;}}d = 1./std::sqrt(x*x + y*y + DBL_EPSILON);//保存特征值l2及其对应的特征向量dst[6*j + 1] = (float)l2;dst[6*j + 4] = (float)(x*d);dst[6*j + 5] = (float)(y*d);}
}

  求得2个特征值α、β和2个特征向量之后,就是要利用特征值构建特征表达式,通过表达式的值(  (αβ) - k(α+β)^2   )来区分角点,k的值通常设置为0.04:

/* calculate Mc */for( int j = 0; j < src_gray.rows; j++ ){ for( int i = 0; i < src_gray.cols; i++ ){float lambda_1 = myHarris_dst.at<Vec6f>(j, i)[0];float lambda_2 = myHarris_dst.at<Vec6f>(j, i)[1];Mc.at<float>(j,i) = lambda_1*lambda_2 - 0.04f*pow( ( lambda_1 + lambda_2 ), 2 );}}

  代码中利用 minMaxLoc( Mc, &myHarris_minVal, &myHarris_maxVal, 0, 0, Mat() ); 函数获取特征表达式的最大值min和最小值max,通过选取不同的阈值min<=thresh<=max,来指定大于阈值thresh的表达式值对应的点为检测出的角点。并利用circle()函数显示出来。

circle( myHarris_copy, Point(i,j), 4, Scalar( rng.uniform(0,255), rng.uniform(0,255), rng.uniform(0,255) ), -1, 8, 0 ); 

至此,Harris角点检测完成!

 2. ShiTomasi角点检测

  ShiTomasi角点提取是获取harris角点中的强角点,怎么获取强角点呢,那就是只选取两个特征值中较小的那个特征值构建特征表达式,如果较小的特征值都能够满足设定的阈值条件,那么该角点就视为强角点。

  调用  cornerMinEigenVal( src_gray, myShiTomasi_dst, blockSize, apertureSize, BORDER_DEFAULT ); 函数来获取较小的特征值,其实该函数内部依然调用上面所述的函数 cornerEigenValsVecs( src, dst, blockSize, ksize, MINEIGENVAL, 0, borderType ); ,然后将op_type设置为MINEIGENVAL枚举值,进而调用 static void calcMinEigenVal( const Mat& _cov, Mat& _dst ) 函数计算较小的特征值。该函数代码如下:

static void calcMinEigenVal( const Mat& _cov, Mat& _dst )
{int i, j;Size size = _cov.size();
#if CV_TRY_AVXbool haveAvx = CV_CPU_HAS_SUPPORT_AVX;
#endif
#if CV_SIMD128bool haveSimd = hasSIMD128();
#endifif( _cov.isContinuous() && _dst.isContinuous() ){size.width *= size.height;size.height = 1;}for( i = 0; i < size.height; i++ ){const float* cov = _cov.ptr<float>(i);float* dst = _dst.ptr<float>(i);
#if CV_TRY_AVXif( haveAvx )j = calcMinEigenValLine_AVX(cov, dst, size.width);else
#endif // CV_TRY_AVXj = 0;#if CV_SIMD128if( haveSimd ){v_float32x4 half = v_setall_f32(0.5f);for( ; j <= size.width - v_float32x4::nlanes; j += v_float32x4::nlanes ){v_float32x4 v_a, v_b, v_c, v_t;v_load_deinterleave(cov + j*3, v_a, v_b, v_c);v_a *= half;v_c *= half;v_t = v_a - v_c;v_t = v_muladd(v_b, v_b, (v_t * v_t));v_store(dst + j, (v_a + v_c) - v_sqrt(v_t));}}
#endif // CV_SIMD128for( ; j < size.width; j++ ){float a = cov[j*3]*0.5f;float b = cov[j*3+1];float c = cov[j*3+2]*0.5f;

//求根公式计算较小的根,即为较小的特征值dst[j] = (float)((a + c) - std::sqrt((a - c)*(a - c) + b*b));}}
}

  所有像素点处较小的特征值求出后,直接将该特征值作为特征表达式的值。利用 minMaxLoc( Mc, &myHarris_minVal, &myHarris_maxVal, 0, 0, Mat() ); 函数选取最小的min和最大的max,通过调整阈值thresh来设定大于阈值thresh的为显示出来的强角点。

  至此,ShiTomasi角点检测完成!

  #自己写了一个简单的Harris和ShiTomasi角点检测的代码,如下,仅供参考:

  1 #include <opencv2\opencv.hpp>
  2 #include <iostream>
  3 #include <string>
  4
  5 using namespace std;
  6
  7 #define HARRIS
  8
  9 cv::RNG rng(12345);
 10
 11 void calEigen2x2(cv::Mat cov,cv::Mat &eigenValue)
 12 {
 13     int height = cov.rows;
 14     int width = cov.cols;
 15
 16     float *pCov = (float*)cov.data;
 17     float *pEigenValue = (float*)eigenValue.data;
 18     for (int i = 0; i < height; i++)
 19     {
 20         for (int j = 0; j < width; j++)
 21         {
 22             double a = pCov[(i*width + j) * 3 + 0];
 23             double b = pCov[(i*width + j) * 3 + 1];
 24             double c = pCov[(i*width + j) * 3 + 2];
 25
 26             double tmp1 = (a + c) / 2.;
 27             double tmp2 = sqrtf(b*b + (a - c)*(a - c) / 4.);
 28
 29             double alpha = tmp1 - tmp2;
 30             double beta = tmp1 + tmp2;
 31
 32             pEigenValue[(i*width + j) * 2 + 0] =(float) alpha;
 33             pEigenValue[(i*width + j) * 2 + 1] =(float) beta;
 34         }
 35     }
 36 }
 37
 38 void myCalEigenValues(cv::Mat srcImg, cv::Mat &eigenValue, int covWin, int sobelWin)
 39 {
 40     //求微分
 41     cv::Mat sobelx, sobely;
 42     cv::Sobel(srcImg, sobelx, CV_32FC1, 1, 0, sobelWin, 1. / (255*12), 0, 4);
 43     cv::Sobel(srcImg, sobely, CV_32FC1, 0, 1, sobelWin, 1. / (255*12), 0, 4);
 44
 45     cv::Mat cov = cv::Mat::zeros(srcImg.size(), CV_32FC3);
 46     int height = srcImg.rows;
 47     int width = srcImg.cols;
 48     float *pSobelX = (float*)sobelx.data;
 49     float *pSobelY = (float*)sobely.data;
 50     float *pCov = (float*)cov.data;
 51     for (int i = 0; i < height; i++)
 52     {
 53         for (int j = 0; j < width; j++)
 54         {
 55             float dx = pSobelX[i*width + j];
 56             float dy = pSobelY[i*width + j];
 57
 58             pCov[(i*width + j) * 3 + 0] = dx*dx;
 59             pCov[(i*width + j) * 3 + 1] = dx*dy;
 60             pCov[(i*width + j) * 3 + 2] = dy*dy;
 61         }
 62     }
 63
 64     cv::boxFilter(cov, cov, cov.depth(), cv::Size(covWin, covWin), cv::Point(-1, -1), false, 4);
 65
 66     calEigen2x2(cov, eigenValue);
 67 }
 68
 69 void main()
 70 {
 71     string imgPath = "data/srcImg/0.png";
 72     cv::Mat srcImg = cv::imread(imgPath, 1);
 73     cv::Mat grayImg;
 74     cv::cvtColor(srcImg, grayImg, CV_BGR2GRAY);
 75     int height = srcImg.rows;
 76     int width = srcImg.cols;
 77
 78     cv::Mat eigenValue = cv::Mat::zeros(grayImg.size(), CV_32FC2);
 79     int covWin = 3, sobelWin = 3;
 80     myCalEigenValues(grayImg, eigenValue, covWin, sobelWin);
 81
 82     cv::Mat Mc = cv::Mat::zeros(grayImg.size(), CV_32FC1);
 83 #ifndef HARRIS
 84     //计算特征表达式值
 85     for (int i = 0; i<height; i++)
 86     {
 87         for (int j = 0; j < width; j++)
 88         {
 89             float alpha = eigenValue.at<float>(i, j*2+0);
 90             float beta = eigenValue.at<float>(i, j*2+1);
 91
 92             Mc.at<float>(i, j) = alpha*beta - 0.04f*pow((alpha + beta), 2);
 93         }
 94     }
 95 #else //ShiTomasi
 96     for (int i = 0; i<height; i++)
 97     {
 98         for (int j = 0; j < width; j++)
 99         {
100             float alpha = eigenValue.at<float>(i, j * 2 + 0);
101             float beta = eigenValue.at<float>(i, j * 2 + 1);
102
103             float minEigenValue = (alpha > beta) ? (beta) : (alpha);
104
105             Mc.at<float>(i, j) = minEigenValue;
106         }
107     }
108 #endif
109
110     double minVal, maxVal;
111     cv::minMaxLoc(Mc, &minVal, &maxVal, 0, 0, cv::Mat());
112
113     double thresh = (maxVal + minVal) / 2.;
114     for (int i = 0; i < height; i++)
115     {
116         for (int j = 0; j < width; j++)
117         {
118             double value = (double)Mc.at<float>(i, j);
119             if (value > thresh)
120             {
121                 cv::circle(srcImg, cv::Point(j, i), 4, cv::Scalar(rng.uniform(0, 255), rng.uniform(0, 255), rng.uniform(0, 255)), -1, 8, 0);
122             }
123         }
124     }
125
126     cv::namedWindow("show", 0);
127     cv::imshow("show", srcImg);
128     cv::waitKey(0);
129 }

View Code


3. cornerHarris()函数详解

  前面讲述cornerEigenValsVecs()这个函数是提到op_type这个枚举类型,有三个枚举值可以设置。其中MINEIGENVAL 和  EIGENVALSVECS都在前面介绍过。而 HARRIS则在cornerHarris()函数中使用,这是一个公开的OpenCV API函数,函数原型如下:

void cv::cornerHarris( InputArray _src, OutputArray _dst, int blockSize, int ksize, double k, int borderType )

  k值为上述特征表达式中的常数项。该函数内部其实还是调用cornerEigenValsVecs()函数,只不过调用时将op_type设置为枚举值HARRIS。意思就是提取HARRIS角点,然后调用内部的 static void calcHarris( const Mat& _cov, Mat& _dst, double k ) 函数。只不过还内部函数不再计算特征值和特征向量,而是直接计算特征表达式的值。而特征表达式用下式表示:

  其中矩阵M就是前面说的协方差矩阵,det(M)为M的行列式,Tr(M)为M的迹。在程序中代码如下:

for( ; j < size.width; j++ ){float a = cov[j*3];float b = cov[j*3+1];float c = cov[j*3+2];

        //求特征表达式的值dst[j] = (float)(a*c - b*b - k*(a + c)*(a + c));}

  通常算出特征表达式的值后将其归一化到(0-255),然后可以直接设置阈值Thresh,特征表达式的值>Thresh对应的点视为角点。具体可以参见OpenCV例程源代码:.\opencv\sources\samples\cpp\tutorial_code\TrackingMotion\cornerHarris_Demo.cpp

 


转载于:https://www.cnblogs.com/riddick/p/8463763.html

OpenCV角点检测源代码分析(Harris和ShiTomasi角点)相关推荐

  1. OpenCV——角点检测原理分析(Harris,Shi-Tomasi、亚像素级角点检测)

    一.角点(corner) 角点通常被定义为两条边的交点,或者说,角点的局部邻域应该具有两个不同区域的不同方向的边界.角点检测(Corner Detection)是计算机视觉系统中获取图像特征的一种方法 ...

  2. OpenCV之feature2d 模块. 2D特征框架(1)Harris 角点检测子 Shi-Tomasi角点检测子 定制化创建角点检测子 亚像素级的角点检测 特征点检测

    Harris 角点检测子 目标 本教程中我们将涉及: 有哪些特征?它们有什么用? 使用函数 cornerHarris 通过 Harris-Stephens方法检测角点. 理论 有哪些特征? 在计算机视 ...

  3. 【OpenCV 学习笔记】第二十章: 角点检测之:harris算法以及Shi-Tomasi算法

    第二十章: 角点检测之:harris算法以及Shi-Tomasi算法 一张图像,我们可以用很多方法去处理它,就会得到很多不同的特征.比如基于梯度方法我们就能得到图像的边缘特征:比如基于直方图我们就得到 ...

  4. harris角点检测c语言,Harris角点检测原理及实现

    为便于理解,先简要介绍角点的概念和角点检测背景 一.角点及角点检测背景: 角点概念: 角点,通常可理解为两条边的角点,也可理解为像素值在多个方向有显著变化的点或局部区域内某个属性明显的点.如多个轮廓的 ...

  5. 特征提取与检测(一)---Harris与Shi-Tomasi角点检测原理

    一.Harris角点检测原理 1. 何为角点? 下面有两幅不同视角的图像,通过找出对应的角点进行匹配. 我们可以直观的概括下角点所具有的特征: >轮廓之间的交点: >对于同一场景,即使视角 ...

  6. 角点检测汇总:Harris角点及Shi-Tomasi角点检测

    一.角点定义 有定义角点的几段话: 1.角点检测(Corner Detection)是计算机视觉系统中用来获得图像特征的一种方法,广泛应用于运动检测.图像匹配.视频跟踪.三维建模和目标识别等领域中.也 ...

  7. 计算机视觉(角点检测)- 1 - Harris角点检测

    计算机视觉(角点检测)- 1 - Harris角点检测 学习前言 一.Harris角点检测   1.什么是角点?   2.Harris角点检测的基本原理&基本思想   3.Harris角点检测 ...

  8. 《OpenCV3编程入门》学习笔记10 角点检测(一)Harris角点检测

    第10章 角点检测 10.1 Harris角点检测 10.1.1 角点 1.图像特征类型: (1)边缘 (2)角点(感兴趣点) (3)斑点(感兴趣区域) 2.角点定义: (1) 一阶导数(灰度的梯度) ...

  9. python+OpenCV笔记(二十四):Shi-Tomasi角点检测

    Shi-Tomasi角点检测 原理 python+OpenCV笔记(二十二):角点检测原理(Harris角点检测原理.Shi-Tomasi角点检测原理)https://blog.csdn.net/qq ...

最新文章

  1. rails的一些问题
  2. Linuxtone命令一句话
  3. 百度飞浆行人多目标跟踪笔记
  4. 编程入门python语言是多大孩子学的-什么是少儿Python编程?这一篇就够啦!
  5. 数据结构与算法 / 栈(stack)
  6. 【数据结构与算法】之“接雨水”的算法求解
  7. 用ajax下载字节流形式的excel文件
  8. 用python写一个简单的爬虫_用Python从零开始写一个简单爬虫
  9. 20191115英文每日一句
  10. 经典 55道 MySQL面试题及答案
  11. 深入理解JVM虚拟机
  12. WIN7 Activation
  13. ATX结合Maxim实现多设备并行执行压力测试(AUI自动化测试框架)
  14. BeautyGAN图片的高精度美颜
  15. 论文笔记 | Conducting research in marketing with quasi-experiments
  16. 【Go语言入门指南】零基础入门 go 语言 | Golang 入门指南
  17. odoo本地文档功能开发记录
  18. elasticsearch-head设置登录用户和密码
  19. React(二)react脚手架的搭建
  20. 软考考c语言还是java,计算机程序设计工程师技术水平(java)证书就是计算机技术与软件专业技术资格考试的程序员证书么?...

热门文章

  1. 移动端 H5 概念术语(一)
  2. 51单片机温度传感器DS18B20
  3. iSCSI网络SCSI接口
  4. ESP8266深度睡眠计时器唤醒
  5. DBLP数据集用weka数据挖掘 xml转csv格式文件
  6. Google登录授权详细过程
  7. 虚幻基础之Gameplay游戏框架之Level和World
  8. 华硕fl5600l装固态并重装系统到固态
  9. 时间序列--残差分析
  10. java模仿滴滴的程序_小程序模仿滴滴打车