算法基础

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

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

算法原理

通过这个图可以看到，大概就是在某一个像素点周围有一个滑动窗口，当这个滑动窗口在各个方向进行小范围移动的时候，我们看一下像素的平均灰度值的变化情况。公式可以表示为：
$E (u, 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) + u I_{x} + v I_{y} - I (x, y)]^{2} = \sum_{x, y} [u I_{x} + v I_{y}]^{2} = (u, v) \sum_{x, y} w (x, y) [\begin{matrix} <mstyle displaystyle="false" scriptlevel="0"> I_{x} I_{y}, I_{x} I_{y} </mstyle> \\ <mstyle displaystyle="false" scriptlevel="0"> I_{x} I_{y}, I_{y} I_{y} </mstyle> \end{matrix}] [\begin{matrix} <mstyle displaystyle="false" scriptlevel="0"> u </mstyle> \\ <mstyle displaystyle="false" scriptlevel="0"> v </mstyle> \end{matrix}]$ ，其中 $\sum_{x, y} w (x, y) [\begin{matrix} <mstyle displaystyle="false" scriptlevel="0"> I_{x} I_{y}, I_{x} I_{y} </mstyle> \\ <mstyle displaystyle="false" scriptlevel="0"> I_{x} I_{y}, I_{y} I_{y} </mstyle> \end{matrix}]$ 被我们叫做结构张量。我们需要寻找一个这样的点，无论我们如何取（u,v），E(u,v)都是变化较大的，这个像素点就是我们需要找的角点。接下来的推导，大家可以去看：https://blog.csdn.net/linqianbi/article/details/78930239 这位博主的，我的理解就是结构张量的 $I_{x}$ 和 $I_{y}$ 的变化情况可以反应像素是边（edge）还是角点(corner)，还是面(flat)。然后通过两个变量来判断角点不诗很方便，所以论文提出了一个角点响应函数R(corner response function), $R = d e t M - k (t r a c e M)^{2}$ ,其中 $d e t M = λ_{1} λ_{2}$ ， $t r a c e M = λ 1 + λ_{2}$ ，k是一个 $[0.04, 0.06]$ 的常数，现在的判断逻辑为:

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

算法步骤

1利用Sobel算子计算除XY方向的梯度值
2计算出 $l x^{2}$ ， $l y^{2}$ , $l x * l y$
3利用高斯函数对 $l x^{2}, l y^{2}, l x * l y$ 进行滤波
4计算局部特征结果矩阵M的特征值和响应函数 $C (i, j) = D e t (M) - k (t r a c e (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);
            }
        }
    }
}

效果

原图
效果图

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

算法基础

算法原理

算法步骤

代码实现

效果

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