基于边缘梯度的模板匹配算法

#include "GeoMatch.h" using namespace cv; GeoMatch::GeoMatch(void) { noOfCordinates = 0; // 初始化 modelDefined = false; } int GeoMatch::CreateGeoMatchModel( Mat templateArr, double maxContrast, double minContrast) { Mat srcstub, src =templateArr; Size size; size.width = src.cols; size.height = src.rows; modelHeight = src.rows; //保存模板图像的尺寸 modelWidth = src.cols; noOfCordinates = 0; //初始化 cordinates = new CvPoint[modelWidth *modelHeight]; //用于模板图像中选定点的坐标//Allocate memory for coorinates of selected points in template image edgeMagnitude = new double[modelWidth *modelHeight]; //用于保存选定点边缘量值 //Allocate memory for edge magnitude for selected points edgeDerivativeX = new double[modelWidth *modelHeight]; //用于保存选定点边缘x方向量值 //Allocate memory for edge X derivative for selected points edgeDerivativeY = new double[modelWidth *modelHeight]; //用于保存选定点边缘y方向量值 ////Allocate memory for edge Y derivative for selected points Mat gx, gy; Sobel(templateArr, gx, CV_16S, 1, 0, 3, 1, 0, BORDER_DEFAULT); Sobel(templateArr, gy, CV_16S, 0, 1, 3, 1, 0, BORDER_DEFAULT); Mat nmsEdges = Mat::zeros(templateArr.rows, templateArr.cols, CV_8UC1); const short* _sdx; const short* _sdy; double fdx, fdy; double MagG, DirG; double MaxGradient = -99999.99; double direction; int *orients = new int[size.height *size.width]; int count = 0, i, j; // count variable; double **magMat; CreateDoubleMatrix(magMat, size); for (i = 1; i < size.height - 1; i++) { for (j = 1; j < size.width - 1; j++) { _sdx = (short*)(gx.ptr<uchar>(i)); _sdy = (short*)(gy.ptr<uchar>(i)); fdx = _sdx[j]; fdy = _sdy[j]; //读入sobel的x，y梯度 // read x, y derivatives MagG = sqrt((float)(fdx*fdx) + (float)(fdy*fdy)); //边界强度 direction = cvFastArctan((float)fdy, (float)fdx); //边界方向 magMat[i][j] = MagG; if (MagG>MaxGradient) MaxGradient = MagG; //正则化 get maximum gradient value for normalizing. // get closest angle from 0, 45, 90, 135 set if ((direction>0 && direction < 22.5) || (direction >157.5 && direction < 202.5) || (direction>337.5 && direction<360)) direction = 0; else if ((direction>22.5 && direction < 67.5) || (direction >202.5 && direction <247.5)) direction = 45; else if ((direction >67.5 && direction < 112.5) || (direction>247.5 && direction<292.5)) direction = 90; else if ((direction >112.5 && direction < 157.5) || (direction>292.5 && direction<337.5)) direction = 135; else direction = 0; orients[count] = (int)direction; count++; } } count = 0; // init count // 非最大抑制算法 non maximum suppression double leftPixel, rightPixel; for (i = 1; i < size.height - 1; i++) { for (j = 1; j < size.width - 1; j++) { switch (orients[count]) { case 0: leftPixel = magMat[i][j - 1]; rightPixel = magMat[i][j + 1]; break; case 45: leftPixel = magMat[i - 1][j + 1]; rightPixel = magMat[i + 1][j - 1]; break; case 90: leftPixel = magMat[i - 1][j]; rightPixel = magMat[i + 1][j]; break; case 135: leftPixel = magMat[i - 1][j - 1]; rightPixel = magMat[i + 1][j + 1]; break; } // 将当前像素值与相邻像素进行比较compare current pixels value with adjacent pixels if ((magMat[i][j] < leftPixel) || (magMat[i][j] < rightPixel)) nmsEdges.at<uchar>(i, j) = 0; else nmsEdges.at<uchar>(i, j) = (uchar)(magMat[i][j] / MaxGradient * 255); count++; } } int RSum = 0, CSum = 0; int curX, curY; int flag = 1; //滞后阈值Hysterisis threshold for (i = 1; i < size.height - 1; i++) { for (j = 1; j < size.width; j++) { _sdx = (short*)(gx.ptr<uchar>(i)); _sdy = (short*)(gy.ptr<uchar>(i)); fdx = _sdx[j]; fdy = _sdy[j]; MagG = sqrt(fdx*fdx + fdy*fdy); //梯度强度Magnitude = Sqrt(gx^2 +gy^2) DirG = cvFastArctan((float)fdy, (float)fdx); //梯度方向Direction = tan(y/x) flag = 1; if (((double)(nmsEdges.at<uchar>(i, j))) < maxContrast) { if ((double)(nmsEdges.at<uchar>(i, j))< minContrast) { (nmsEdges.at<uchar>(i, j)) = 0; flag = 0; // 不是边缘，remove from edge } else { // 如果8个相邻像素中的任何一个不大于最大值，则从边缘移除 //if any of 8 neighboring pixel is not greater than max contraxt remove from edge if ((((double)(nmsEdges.at<uchar>(i-1, j-1))) < maxContrast) && (((double)(nmsEdges.at<uchar>(i - 1, j))) < maxContrast) && (((double)((nmsEdges.at<uchar>(i - 1, j + 1)))) < maxContrast) && (((double)((nmsEdges.at<uchar>(i , j - 1)))) < maxContrast) && (((double)((nmsEdges.at<uchar>(i , j + 1)))) < maxContrast) && (((double)((nmsEdges.at<uchar>(i + 1, j - 1)))) < maxContrast) && (((double)((nmsEdges.at<uchar>(i + 1, j )))) < maxContrast) && (((double)((nmsEdges.at<uchar>(i + 1, j + 1)))) < maxContrast)) { (nmsEdges.at<uchar>(i , j)) = 0; flag = 0; } } } // 保存提取到的边缘信息save selected edge information curX = i; curY = j; if (flag != 0) { if (fdx != 0 || fdy != 0) { RSum = RSum + curX; CSum = CSum + curY; // Row sum and column sum for center of gravity cordinates[noOfCordinates].x = curX; cordinates[noOfCordinates].y = curY; edgeDerivativeX[noOfCordinates] = fdx; edgeDerivativeY[noOfCordinates] = fdy; //handle divide by zero if (MagG != 0) edgeMagnitude[noOfCordinates] = 1 / MagG; // gradient magnitude else edgeMagnitude[noOfCordinates] = 0; noOfCordinates++; } } } } centerOfGravity.x = RSum / noOfCordinates; // center of gravity centerOfGravity.y = CSum / noOfCordinates; // center of gravity // change coordinates to reflect center of gravity for (int m = 0; m<noOfCordinates; m++) { int temp; temp = cordinates[m].x; cordinates[m].x = temp - centerOfGravity.x; temp = cordinates[m].y; cordinates[m].y = temp - centerOfGravity.y; } delete[] orients; modelDefined = true; return 1; } double GeoMatch::FindGeoMatchModel(Mat srcarr, double minScore, double greediness, Point *resultPoint) { double resultScore = 0; double partialSum = 0; double sumOfCoords = 0; double partialScore; const short* _Sdx; const short* _Sdy; int i, j, m; // count variables double iTx, iTy, iSx, iSy; double gradMag; int curX, curY; //pyrDown(srcarr, srcarr, Size(srcarr.cols * 0.5, srcarr.rows * 0.5)); double **matGradMag; //梯度强度矩阵 Mat src = srcarr; Size size; size.width = src.cols; size.height = src.rows; CreateDoubleMatrix(matGradMag, size); // 保存梯度强度值 Mat gx, gy; Sobel(src, gx, CV_16S, 1, 0, 3, 1, 0, BORDER_DEFAULT); Sobel(src, gy, CV_16S, 0, 1, 3, 1, 0, BORDER_DEFAULT); // stoping criterias to search for model double normMinScore = minScore / noOfCordinates; // precompute minumum score double normGreediness = ((1 - greediness * minScore) / (1 - greediness)) / noOfCordinates; // precompute greedniness for (i = 0; i < size.height; i++) { _Sdx = (short*)(gx.ptr<uchar>(i)); _Sdy = (short*)(gy.ptr<uchar>(i)); for (j = 0; j < size.width; j++) { iSx = _Sdx[j]; // 源图像的X梯度X derivative of Source image iSy = _Sdy[j]; // 源图像的Y梯度 Y derivative of Source image gradMag = sqrt((iSx*iSx) + (iSy*iSy)); //梯度强度Magnitude = Sqrt(dx^2 +dy^2) if (gradMag != 0) // 避免0做分母 hande divide by zero matGradMag[i][j] = 1 / gradMag; // 1/Sqrt(dx^2 +dy^2) else matGradMag[i][j] = 0; } } for (i = 0; i < size.height; i++) { for (j = 0; j < size.width; j++) { partialSum = 0; // initilize partialSum measure for (m = 0; m<noOfCordinates; m++) { curX = i + cordinates[m].x