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 个特征向量, 该函数如下: 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,l2
- double 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_AVX
- bool haveAvx = CV_CPU_HAS_SUPPORT_AVX;
- #endif
- #if CV_SIMD128
- bool haveSimd = hasSIMD128();
- #endif
- if( _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_AVX
- if( haveAvx )
- j = calcMinEigenValLine_AVX(cov, dst, size.width);
- else
- #endif // CV_TRY_AVX
- j = 0;
- #if CV_SIMD128
- if( 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_SIMD128
- for( ; 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 角点检测完成!
来源: https://www.cnblogs.com/riddick/p/8463763.html