以下6个函数是opencv有关角点检测的函数 ConerHarris, cornoerMinEigenVal,CornorEigenValsAndVecs, preConerDetect, conerSubPix, goodFeaturesToTracks, 其中， 前三个都调用静态函数cornerEigenValsVecs
1、静态函数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 ) { #ifdef HAVE_TEGRA_OPTIMIZATION if (tegra::useTegra() && tegra::cornerEigenValsVecs(src, eigenv, block_size, aperture_size, op_type, k, borderType)) return; #endif #if CV_TRY_AVX bool haveAvx = CV_CPU_HAS_SUPPORT_AVX; #endif #if CV_SIMD128 bool haveSimd = hasSIMD128(); #endif int depth = src.depth(); double scale = (double)(1 << ((aperture_size > 0 ? aperture_size : 3) - 1)) * block_size; if( aperture_size < 0 ) scale *= 2.0; if( depth == CV_8U ) scale *= 255.0; scale = 1.0/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 = cov.ptr<float>(i); const float* dxdata = Dx.ptr<float>(i); const float* dydata = Dy.ptr<float>(i); #if CV_TRY_AVX if( haveAvx ) j = cornerEigenValsVecsLine_AVX(dxdata, dydata, cov_data, size.width); else #endif // CV_TRY_AVX j = 0; #if CV_SIMD128 if( haveSimd ) { for( ; j <= size.width - v_float32x4::nlanes; j += v_float32x4::nlanes ) { v_float32x4 v_dx = v_load(dxdata + j); v_float32x4 v_dy = v_load(dydata + j); v_float32x4 v_dst0, v_dst1, v_dst2; v_dst0 = v_dx * v_dx; v_dst1 = v_dx * v_dy; v_dst2 = v_dy * v_dy; v_store_interleave(cov_data + j * 3, v_dst0, v_dst1, v_dst2); } } #endif // CV_SIMD128 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 ); }
2、preConerDetect函数分析；
preConerDetect： 角点检测的特征图， 利用二阶导数求解角点，二阶导数为零表示角点；
计算Dx^2Dyy + Dy^2Dxx - 2DxDyDxy
void cv::preCornerDetect( InputArray _src, OutputArray _dst, int ksize, int borderType ) { CV_INSTRUMENT_REGION() int type = _src.type(); CV_Assert( type == CV_8UC1 || type == CV_32FC1 ); CV_OCL_RUN( _src.dims() <= 2 && _dst.isUMat(), ocl_preCornerDetect(_src, _dst, ksize, borderType, CV_MAT_DEPTH(type))) Mat Dx, Dy, D2x, D2y, Dxy, src = _src.getMat(); _dst.create( src.size(), CV_32FC1 ); Mat dst = _dst.getMat(); Sobel( src, Dx, CV_32F, 1, 0, ksize, 1, 0, borderType ); Sobel( src, Dy, CV_32F, 0, 1, ksize, 1, 0, borderType ); Sobel( src, D2x, CV_32F, 2, 0, ksize, 1, 0, borderType ); Sobel( src, D2y, CV_32F, 0, 2, ksize, 1, 0, borderType ); Sobel( src, Dxy, CV_32F, 1, 1, ksize, 1, 0, borderType ); double factor = 1 << (ksize - 1); if( src.depth() == CV_8U ) factor *= 255; factor = 1./(factor * factor * factor); #if CV_SIMD128 float factor_f = (float)factor; bool haveSimd = hasSIMD128(); v_float32x4 v_factor = v_setall_f32(factor_f), v_m2 = v_setall_f32(-2.0f); #endif Size size = src.size(); int i, j; for( i = 0; i < size.height; i++ ) { float* dstdata = dst.ptr<float>(i); const float* dxdata = Dx.ptr<float>(i); const float* dydata = Dy.ptr<float>(i); const float* d2xdata = D2x.ptr<float>(i); const float* d2ydata = D2y.ptr<float>(i); const float* dxydata = Dxy.ptr<float>(i); j = 0; #if CV_SIMD128 if (haveSimd) { for( ; j <= size.width - v_float32x4::nlanes; j += v_float32x4::nlanes ) { v_float32x4 v_dx = v_load(dxdata + j); v_float32x4 v_dy = v_load(dydata + j); v_float32x4 v_s1 = (v_dx * v_dx) * v_load(d2ydata + j); v_float32x4 v_s2 = v_muladd((v_dy * v_dy), v_load(d2xdata + j), v_s1); v_float32x4 v_s3 = v_muladd((v_dy * v_dx) * v_load(dxydata + j), v_m2, v_s2); v_store(dstdata + j, v_s3 * v_factor); } } #endif for( ; j < size.width; j++ ) { float dx = dxdata[j]; float dy = dydata[j]; dstdata[j] = (float)(factor*(dx*dx*d2ydata[j] + dy*dy*d2xdata[j] - 2*dx*dy*dxydata[j])); //计算特征图； } } }
3、cornorSubPix函数分析；
角点的真实位置并不一定在整数像素位置，因而为了获取更为精确的角点位置坐标，需要角点坐标达到亚像素级精度；
原理： https://xueyayang.github.io/pdf_posts/%E4%BA%9A%E5%83%8F%E7%B4%A0%E8%A7%92%E7%82%B9%E7%9A%84%E6%B1%82%E6%B3%95.pdf
//亚像素级角点检测； void cv::cornerSubPix( InputArray _image, InputOutputArray _corners, Size win, Size zeroZone, TermCriteria criteria ) { CV_INSTRUMENT_REGION() const int MAX_ITERS = 100; int win_w = win.width * 2 + 1, win_h = win.height * 2 + 1; int i, j, k; int max_iters = (criteria.type & CV_TERMCRIT_ITER) ? MIN(MAX(criteria.maxCount, 1), MAX_ITERS) : MAX_ITERS; double eps = (criteria.type & CV_TERMCRIT_EPS) ? MAX(criteria.epsilon, 0.) : 0; eps *= eps; // use square of error in comparsion operations cv::Mat src = _image.getMat(), cornersmat = _corners.getMat(); int count = cornersmat.checkVector(2, CV_32F); CV_Assert( count >= 0 ); Point2f* corners = cornersmat.ptr<Point2f>(); if( count == 0 ) return; CV_Assert( win.width > 0 && win.height > 0 ); CV_Assert( src.cols >= win.width*2 + 5 && src.rows >= win.height*2 + 5 ); CV_Assert( src.channels() == 1 ); Mat maskm(win_h, win_w, CV_32F), subpix_buf(win_h+2, win_w+2, CV_32F); float* mask = maskm.ptr<float>(); for( i = 0; i < win_h; i++ ) { float y = (float)(i - win.height)/win.height; float vy = std::exp(-y*y); for( j = 0; j < win_w; j++ ) { float x = (float)(j - win.width)/win.width; mask[i * win_w + j] = (float)(vy*std::exp(-x*x)); } } // make zero_zone if( zeroZone.width >= 0 && zeroZone.height >= 0 && zeroZone.width * 2 + 1 < win_w && zeroZone.height * 2 + 1 < win_h ) { for( i = win.height - zeroZone.height; i <= win.height + zeroZone.height; i++ ) { for( j = win.width - zeroZone.width; j <= win.width + zeroZone.width; j++ ) { mask[i * win_w + j] = 0; } } } // do optimization loop for all the points for( int pt_i = 0; pt_i < count; pt_i++ ) { Point2f cT = corners[pt_i], cI = cT; int iter = 0; double err = 0; do { Point2f cI2; double a = 0, b = 0, c = 0, bb1 = 0, bb2 = 0; getRectSubPix(src, Size(win_w+2, win_h+2), cI, subpix_buf, subpix_buf.type()); const float* subpix = &subpix_buf.at<float>(1,1); // process gradient for( i = 0, k = 0; i < win_h; i++, subpix += win_w + 2 ) { double py = i - win.height; for( j = 0; j < win_w; j++, k++ ) { double m = mask[k]; double tgx = subpix[j+1] - subpix[j-1]; double tgy = subpix[j+win_w+2] - subpix[j-win_w-2]; double gxx = tgx * tgx * m; double gxy = tgx * tgy * m; double gyy = tgy * tgy * m; double px = j - win.width; a += gxx; b += gxy; c += gyy; bb1 += gxx * px + gxy * py; bb2 += gxy * px + gyy * py; } } double det=a*c-b*b; if( fabs( det ) <= DBL_EPSILON*DBL_EPSILON ) break; // 2x2 matrix inversion double scale=1.0/det; cI2.x = (float)(cI.x + c*scale*bb1 - b*scale*bb2); cI2.y = (float)(cI.y - b*scale*bb1 + a*scale*bb2); err = (cI2.x - cI.x) * (cI2.x - cI.x) + (cI2.y - cI.y) * (cI2.y - cI.y); cI = cI2; if( cI.x < 0 || cI.x >= src.cols || cI.y < 0 || cI.y >= src.rows ) break; } while( ++iter < max_iters && err > eps ); // if new point is too far from initial, it means poor convergence. // leave initial point as the result if( fabs( cI.x - cT.x ) > win.width || fabs( cI.y - cT.y ) > win.height ) cI = cT; corners[pt_i] = cI; } }
4、goodFeaturesToTrack函数分析；
goodFeaturesToTrack是对hariis的一种改进算法---Shi-Tomasi角点检测算子。
void cv::goodFeaturesToTrack( InputArray _image, OutputArray _corners, int maxCorners, double qualityLevel, double minDistance, InputArray _mask, int blockSize, bool useHarrisDetector, double harrisK ) { CV_INSTRUMENT_REGION() CV_Assert( qualityLevel > 0 && minDistance >= 0 && maxCorners >= 0 ); CV_Assert( _mask.empty() || (_mask.type() == CV_8UC1 && _mask.sameSize(_image)) ); CV_OCL_RUN(_image.dims() <= 2 && _image.isUMat(), ocl_goodFeaturesToTrack(_image, _corners, maxCorners, qualityLevel, minDistance, _mask, blockSize, useHarrisDetector, harrisK)) Mat image = _image.getMat(), eig, tmp; if (image.empty()) { _corners.release(); return; } // Disabled due to bad accuracy CV_OVX_RUN(false && useHarrisDetector && _mask.empty() && !ovx::skipSmallImages<VX_KERNEL_HARRIS_CORNERS>(image.cols, image.rows), openvx_harris(image, _corners, maxCorners, qualityLevel, minDistance, blockSize, harrisK)) if( useHarrisDetector ) cornerHarris( image, eig, blockSize, 3, harrisK ); else cornerMinEigenVal( image, eig, blockSize, 3 ); double maxVal = 0; minMaxLoc( eig, 0, &maxVal, 0, 0, _mask ); threshold( eig, eig, maxVal*qualityLevel, 0, THRESH_TOZERO ); dilate( eig, tmp, Mat()); Size imgsize = image.size(); std::vector<const float*> tmpCorners; // collect list of pointers to features - put them into temporary image Mat mask = _mask.getMat(); for( int y = 1; y < imgsize.height - 1; y++ ) { const float* eig_data = (const float*)eig.ptr(y); const float* tmp_data = (const float*)tmp.ptr(y); const uchar* mask_data = mask.data ? mask.ptr(y) : 0; for( int x = 1; x < imgsize.width - 1; x++ ) { float val = eig_data[x]; if( val != 0 && val == tmp_data[x] && (!mask_data || mask_data[x]) ) tmpCorners.push_back(eig_data + x); } } std::vector<Point2f> corners; size_t i, j, total = tmpCorners.size(), ncorners = 0; if (total == 0) { _corners.release(); return; } std::sort( tmpCorners.begin(), tmpCorners.end(), greaterThanPtr() ); if (minDistance >= 1) //邻域处理； { // Partition the image into larger grids int w = image.cols; int h = image.rows; const int cell_size = cvRound(minDistance); const int grid_width = (w + cell_size - 1) / cell_size; const int grid_height = (h + cell_size - 1) / cell_size; std::vector<std::vector<Point2f> > grid(grid_width*grid_height); minDistance *= minDistance; for( i = 0; i < total; i++ ) { int ofs = (int)((const uchar*)tmpCorners[i] - eig.ptr()); int y = (int)(ofs / eig.step); int x = (int)((ofs - y*eig.step)/sizeof(float)); bool good = true; int x_cell = x / cell_size; int y_cell = y / cell_size; int x1 = x_cell - 1; int y1 = y_cell - 1; int x2 = x_cell + 1; int y2 = y_cell + 1; // boundary check x1 = std::max(0, x1); y1 = std::max(0, y1); x2 = std::min(grid_width-1, x2); y2 = std::min(grid_height-1, y2); for( int yy = y1; yy <= y2; yy++ ) //移动； { for( int xx = x1; xx <= x2; xx++ ) { std::vector <Point2f> &m = grid[yy*grid_width + xx]; if( m.size() ) { for(j = 0; j < m.size(); j++) { float dx = x - m[j].x; float dy = y - m[j].y; if( dx*dx + dy*dy < minDistance ) //邻域内有角点的话， 就剔除； { good = false; goto break_out; } } } } } break_out: if (good) { grid[y_cell*grid_width + x_cell].push_back(Point2f((float)x, (float)y)); corners.push_back(Point2f((float)x, (float)y)); ++ncorners; if( maxCorners > 0 && (int)ncorners == maxCorners ) break; } } } else { for( i = 0; i < total; i++ ) { int ofs = (int)((const uchar*)tmpCorners[i] - eig.ptr()); int y = (int)(ofs / eig.step); int x = (int)((ofs - y*eig.step)/sizeof(float)); corners.push_back(Point2f((float)x, (float)y)); ++ncorners; if( maxCorners > 0 && (int)ncorners == maxCorners ) break; } } Mat(corners).convertTo(_corners, _corners.fixedType() ? _corners.type() : CV_32F); }
注： 该博文为扩展型；