Harris角点检测算法原理及C++实现

关注

原创最新推荐文章于 2026-05-31 00:04:03 发布 4.8k 阅读 AI 写同款· GEO 优化 ›

【图像预处理】 Harris角点检测器原理及C++实现在图像中，点的位置很多时候是很有用的。我们经常会在图像处理的过程中提取我们的兴趣点，而角点是很好的兴趣点。在解决许多实际问题时，使用角点特征的效果优于线（或者边缘）。这是因为观察一条移动直线时，只有其垂直于线段的运动分量能被观测到，而语线段共线的分量则不可见。角点则不同，他们提供唯一匹配。 Harris角点检测器是一种比较常用的检测器，它对二维平移和旋转、少量光照变化和视角变化具有鲁棒性，并且计算量阅读详情

算法基础

首先要知道什么是角点？角点指的是两条边的交点，在图像处理中是指那些可以表示像素变化，例如局部最大最小的灰度的点。图像的特征类型一般为：

边缘
角点（感兴趣的关键点）
斑点（感兴趣区域）
角点在保留图像的重要特征的同时，可以有效的减少信息的计算量，使其信息含量很高，可以有效的提高计算速度，和有利于图像的匹配，使得实时处理成为可能。

算法原理

在这里插入图片描述通过这个图可以看到，大概就是在某一个像素点周围有一个滑动窗口，当这个滑动窗口在各个方向进行小范围移动的时候，我们看一下像素的平均灰度值的变化情况。公式可以表示为：
$v)=\sum_{x,y}w(x,y)[I(x+u, y+v)-I(x,y)]^2$ ，其中 $u, v$ 是窗口在水平，竖直方向的偏移，然后，这里要对在(x+u,y+v)在(x,y)点进行二阶泰勒级数展开，然后取一阶近似，这个公式的展开，可以百度一下二阶泰勒级数的近视，得到 $E (u, v)$ 约等于 $\sum_{x,y}w(x,y)[I(x+y)+uI_x+vI_y-I(x,y)]^2=\sum_{x,y}[uI_x+vI_y]^2=(u ,v)\sum_{x,y}w(x,y)\begin{bmatrix}I_xI_y , I_xI_y\\ I_xI_y, I_yI_y\end{bmatrix}\begin{bmatrix}u\\v \end{bmatrix}$ ，其中 $\sum_{x,y}w(x,y)\begin{bmatrix}I_xI_y , I_xI_y\\ I_xI_y, I_yI_y\end{bmatrix}$ 被我们叫做结构张量。我们需要寻找一个这样的点，无论我们如何取（u,v），E(u,v)都是变化较大的，这个像素点就是我们需要找的角点。接下来的推导，大家可以去看：https://blog.youkuaiyun.com/linqianbi/article/details/78930239 这位博主的，我的理解就是结构张量的 $I_x$ 和 $I_y$ 的变化情况可以反应像素是边（edge）还是角点(corner)，还是面(flat)。然后通过两个变量来判断角点不诗很方便，所以论文提出了一个角点响应函数R(corner response function), $R=detM-k(trace M)^2$ ,其中 $\lambda_1\lambda_2$ ， $M=\lambda1+\lambda_2$ ，k是一个 $[0.04, 0.06]$ 的常数，现在的判断逻辑为:

角点：R为大数值整数
边缘：R为大数值负数
平坦区：绝对值R是小数值

算法步骤

1利用Sobel算子计算除XY方向的梯度值
2计算出 $lx^2$ ， $ly^2$ , $l x * l y$
3利用高斯函数对 $lx^2,ly^2,lx*ly$ 进行滤波
4计算局部特征结果矩阵M的特征值和响应函数 $C(i,j)=Det(M)-k(trace(M))^2(0.04<=k<=0.06)$
5将计算出响应函数的值C进行非极大值抑制，滤除一些不是角点的点，同时满足大于设定的阈值

代码实现

Mat RGB2GRAY(const Mat &src){
    if(!src.data || src.channels()!=3){
        exit(1);
    }
    int row = src.rows;
    int col = src.cols;
    Mat gray = Mat(row, col, CV_8UC1);
    for(int i = 0; i < row; i++){
        for(int j = 0; j < col; j++){
            gray.at<uchar>(i, j) = (uchar)(0.114 * src.at<Vec3b>(i, j)[0] + 0.587 * src.at<Vec3b>(i, j)[1] + 0.299 * src.at<Vec3b>(i, j)[2]);
        }
    }
    return gray;
}

void SobelGradDirction(Mat &src, Mat &sobelX, Mat &sobelY){
    int row = src.rows;
    int col = src.cols;
    sobelX = Mat::zeros(src.size(), CV_32SC1);
    sobelY = Mat::zeros(src.size(), CV_32SC1);

    for(int i = 0; i < row-1; i++){
        for(int j = 0; j < col-1; j++){
            double gradY = src.at<uchar>(i+1, j-1) + src.at<uchar>(i+1, j) * 2 + src.at<uchar>(i+1, j+1) - src.at<uchar>(i-1, j-1) - src.at<uchar>(i-1, j) * 2 - src.at<uchar>(i-1, j+1);
            sobelY.at<float>(i, j) = abs(gradY);
            double gradX = src.at<uchar>(i-1, j+1) + src.at<uchar>(i, j+1) * 2 + src.at<uchar>(i+1, j+1) - src.at<uchar>(i-1, j-1) - src.at<uchar>(i, j-1) * 2 - src.at<uchar>(i+1, j-1);
            sobelX.at<float>(i, j) = abs(gradX);
        }
    }
    //将梯度数组转换为8位无符号整形
    convertScaleAbs(sobelX, sobelX);
    convertScaleAbs(sobelY, sobelY);
}

Mat SobelXX(const Mat src){
    int row = src.rows;
    int col = src.cols;
    Mat_<float> dst(row, col, CV_32FC1);
    for(int i = 0; i < row; i++){
        for(int j = 0; j < col; j++){
            dst.at<float>(i, j) = src.at<uchar>(i, j) * src.at<uchar>(i, j);
        }
    }
    return dst;
}

Mat SobelYY(const Mat src){
    int row = src.rows;
    int col = src.cols;
    Mat_<float> dst(row, col, CV_32FC1);
    for(int i = 0; i < row; i++){
        for(int j = 0; j < col; j++){
            dst.at<float>(i, j) = src.at<uchar>(i, j) * src.at<uchar>(i, j);
        }
    }
    return dst;
}

Mat SobelXY(const Mat src1, const Mat &src2){
    int row = src1.rows;
    int col = src1.cols;
    Mat_<float> dst(row, col, CV_32FC1);
    for(int i = 0; i < row; i++){
        for(int j = 0; j < col; j++){
            dst.at<float>(i, j) = src1.at<uchar>(i, j) * src2.at<uchar>(i, j);
        }
    }
    return dst;
}

void separateGaussianFilter(Mat_<float> &src, Mat_<float> &dst, int ksize, double sigma){
    CV_Assert(src.channels()==1 || src.channels() == 3); //只处理单通道或者三通道图像
    //生成一维的
    double *matrix = new double[ksize];
    double sum = 0;
    int origin = ksize / 2;
    for(int i = 0; i < ksize; i++){
        double g = exp(-(i-origin) * (i-origin) / (2 * sigma * sigma));
        sum += g;
        matrix[i] = g;
    }
    for(int i = 0; i < ksize; i++) matrix[i] /= sum;
    int border = ksize / 2;
    copyMakeBorder(src, dst, border, border, border, border, BORDER_CONSTANT);
    int channels = dst.channels();
    int rows = dst.rows - border;
    int cols = dst.cols - border;
    //水平方向
    for(int i = border; i < rows - border; i++){
        for(int j = border; j < cols - border; j++){
            float sum[3] = {0};
            for(int k = -border; k<=border; k++){
                if(channels == 1){
                    sum[0] += matrix[border + k] * dst.at<float>(i, j+k);
                }else if(channels == 3){
                    Vec3f rgb = dst.at<Vec3f>(i, j+k);
                    sum[0] += matrix[border+k] * rgb[0];
                    sum[1] += matrix[border+k] * rgb[1];
                    sum[2] += matrix[border+k] * rgb[2];
                }
            }
            for(int k = 0; k < channels; k++){
                if(sum[k] < 0) sum[k] = 0;
                else if(sum[k] > 255) sum[k] = 255;
            }
            if(channels == 1)
                dst.at<float>(i, j) = static_cast<float>(sum[0]);
            else if(channels == 3){
                Vec3f rgb = {static_cast<float>(sum[0]), static_cast<float>(sum[1]), static_cast<float>(sum[2])};
                dst.at<Vec3f>(i, j) = rgb;
            }
        }
    }
    //竖直方向
    for(int i = border; i < rows - border; i++){
        for(int j = border; j < cols - border; j++){
            float sum[3] = {0};
            for(int k = -border; k<=border; k++){
                if(channels == 1){
                    sum[0] += matrix[border + k] * dst.at<float>(i+k, j);
                }else if(channels == 3){
                    Vec3f rgb = dst.at<Vec3f>(i+k, j);
                    sum[0] += matrix[border+k] * rgb[0];
                    sum[1] += matrix[border+k] * rgb[1];
                    sum[2] += matrix[border+k] * rgb[2];
                }
            }
            for(int k = 0; k < channels; k++){
                if(sum[k] < 0) sum[k] = 0;
                else if(sum[k] > 255) sum[k] = 255;
            }
            if(channels == 1)
                dst.at<float>(i, j) = static_cast<float>(sum[0]);
            else if(channels == 3){
                Vec3f rgb = {static_cast<float>(sum[0]), static_cast<float>(sum[1]), static_cast<float>(sum[2])};
                dst.at<Vec3f>(i, j) = rgb;
            }
        }
    }
    delete [] matrix;
}

Mat harrisResponse(Mat_<float> GaussXX, Mat_<float> GaussYY, Mat_<float> GaussXY, float k){
    int row = GaussXX.rows;
    int col = GaussXX.cols;
    Mat_<float> dst(row, col, CV_32FC1);
    for(int i = 0; i < row; i++){
        for(int j = 0; j < col; j++){
            float a = GaussXX.at<float>(i, j);
            float b = GaussYY.at<float>(i, j);
            float c = GaussXY.at<float>(i, j);
            dst.at<float>(i, j) = a * b - c * c - k * (a + b) * (a + b);
        }
    }
    return dst;
}

int dir[8][2] = {0, 1, 0, -1, 1, 0, -1, 0, 1, 1, 1, -1, -1, 1, -1, -1};

void MaxLocalValue(Mat_<float>&resultData, Mat srcGray, Mat &resultImage, int kSize){
    int r = kSize / 2;
    resultImage = srcGray.clone();
    int row = resultImage.rows;
    int col = resultImage.cols;
    for(int i = r; i < row - r; i++){
        for(int j = r; j < col - r; j++){
            bool flag = true;
            for(int k = 0; k < 8; k++){
                int tx = i + dir[k][0];
                int ty = j + dir[k][1];
                if(resultData.at<float>(i, j) < resultData.at<float>(tx, ty)){
                    flag = false;
                }
            }
            if(flag && (int)resultData.at<float>(i, j) > 18000){
                circle(resultImage, Point(i, j), 5, Scalar(0, 0, 255), 2, 8, 0);
            }
        }
    }
}