【特徵匹配】Harris及Shi-Tomasi原理及原始碼解析
阿新 • • 發佈:2019-01-01
演算法原理:呼叫cornerMinEigenVal()函式求出每個畫素點自適應矩陣M的較小特徵值,儲存在矩陣eig中,然後找到矩陣eig中最大的畫素值記為maxVal,然後閾值處理,小於qualityLevel*maxVal的特徵值排除掉,最後函式確保所有發現的角點之間具有足夠的距離。
void cv::goodFeaturesToTrack( InputArray _image, OutputArray _corners, int maxCorners, double qualityLevel, double minDistance, InputArray _mask, int blockSize, bool useHarrisDetector, double harrisK ) { Mat image = _image.getMat(), mask = _mask.getMat(); CV_Assert( qualityLevel > 0 && minDistance >= 0 && maxCorners >= 0 ); CV_Assert( mask.empty() || (mask.type() == CV_8UC1 && mask.size() == image.size()) ); Mat eig, tmp; if( useHarrisDetector ) cornerHarris( image, eig, blockSize, 3, harrisK ); //採用Harris角點檢測 else cornerMinEigenVal( image, eig, blockSize, 3 ); //採用Harris改進演算法,eig儲存矩陣M較小的特徵值。見下面演算法實現 double maxVal = 0; minMaxLoc( eig, 0, &maxVal, 0, 0, mask );//儲存eig中最大的值maxVal threshold( eig, eig, maxVal*qualityLevel, 0, THRESH_TOZERO );//閾值處理,小於maxVal*qualityLevel的畫素值歸為0。 dilate( eig, tmp, Mat());//膨脹,3×3的核,為了取區域性極大值 Size imgsize = image.size(); vector<const float*> tmpCorners; // collect list of pointers to features - put them into temporary image 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); } } sort( tmpCorners, greaterThanPtr<float>() ); //按值從大到小排序 vector<Point2f> corners; size_t i, j, total = tmpCorners.size(), ncorners = 0; /* 網格處理,即把影象劃分成正方形網格,每個網格邊長為容忍距離minDistance 以一個角點位置為中心,minDistance為半徑的區域內部不允許出現第二個角點 */ 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.data); 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++ )//檢測角點,minDistance半徑鄰域內,有沒有其他角點出現 { for( int xx = x1; xx <= x2; xx++ ) { 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) { // printf("%d: %d %d -> %d %d, %d, %d -- %d %d %d %d, %d %d, c=%d\n", // i,x, y, x_cell, y_cell, (int)minDistance, cell_size,x1,y1,x2,y2, grid_width,grid_height,c); grid[y_cell*grid_width + x_cell].push_back(Point2f((float)x, (float)y)); corners.push_back(Point2f((float)x, (float)y));//滿足條件的存入corners ++ncorners; if( maxCorners > 0 && (int)ncorners == maxCorners ) break; } } } else //不設定容忍距離 { for( i = 0; i < total; i++ ) { int ofs = (int)((const uchar*)tmpCorners[i] - eig.data); 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); }
求矩陣M最小的特徵值
static void calcMinEigenVal( const Mat& _cov, Mat& _dst ) { int i, j; Size size = _cov.size(); if( _cov.isContinuous() && _dst.isContinuous() ) { size.width *= size.height; size.height = 1; } for( i = 0; i < size.height; i++ )//遍歷所有畫素點 { const float* cov = (const float*)(_cov.data + _cov.step*i); float* dst = (float*)(_dst.data + _dst.step*i); j = 0; for( ; j < size.width; j++ ) { float a = cov[j*3]*0.5f;//cov[j*3]儲存矩陣M左上角元素 float b = cov[j*3+1]; //cov[j*3+1]儲存左下角和右上角元素 float c = cov[j*3+2]*0.5f;//cov[j*3+2]右下角元素 dst[j] = (float)((a + c) - std::sqrt((a - c)*(a - c) + b*b));//求最小特徵值,一元二次方程求根公式 } } }