算法基础

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

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

算法原理

通过这个图可以看到,大概就是在某一个像素点周围有一个滑动窗口,当这个滑动窗口在各个方向进行小范围移动的时候,我们看一下像素的平均灰度值的变化情况。公式可以表示为:
E ( u , v ) = x , y w ( x , y ) [ I ( x + u , y + v ) I ( x , y ) ] 2 E(u, v)=\sum_{x,y}w(x,y)[I(x+u, y+v)-I(x,y)]^2 E(u,v)=x,yw(x,y)[I(x+u,y+v)I(x,y)]2,其中 u , v u,v u,v是窗口在水平,竖直方向的偏移,然后,这里要对在(x+u,y+v)在(x,y)点进行二阶泰勒级数展开,然后取一阶近似,这个公式的展开,可以百度一下二阶泰勒级数的近视,得到 E ( u , v ) E(u,v) E(u,v)约等于 x , y w ( x , y ) [ I ( x + y ) + u I x + v I y I ( x , y ) ] 2 = x , y [ u I x + v I y ] 2 = ( u , v ) x , y w ( x , y ) [ <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> ] [ <mstyle displaystyle="false" scriptlevel="0"> u </mstyle> <mstyle displaystyle="false" scriptlevel="0"> v </mstyle> ] \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} x,yw(x,y)[I(x+y)+uIx+vIyI(x,y)]2=x,y[uIx+vIy]2=(u,v)x,yw(x,y)[IxIy,IxIyIxIy,IyIy][uv],其中 x , y w ( x , y ) [ <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> ] \sum_{x,y}w(x,y)\begin{bmatrix}I_xI_y , I_xI_y\\ I_xI_y, I_yI_y\end{bmatrix} x,yw(x,y)[IxIy,IxIyIxIy,IyIy]被我们叫做结构张量。我们需要寻找一个这样的点,无论我们如何取(u,v),E(u,v)都是变化较大的,这个像素点就是我们需要找的角点。接下来的推导,大家可以去看:https://blog.csdn.net/linqianbi/article/details/78930239 这位博主的,我的理解就是结构张量的 I x I_x Ix I y I_y Iy的变化情况可以反应像素是边(edge)还是角点(corner),还是面(flat)。然后通过两个变量来判断角点不诗很方便,所以论文提出了一个角点响应函数R(corner response function), R = d e t M k ( t r a c e M ) 2 R=detM-k(trace M)^2 R=detMk(traceM)2,其中 d e t M = λ 1 λ 2 detM = \lambda_1\lambda_2 detM=λ1λ2 t r a c e M = λ 1 + λ 2 trace M=\lambda1+\lambda_2 traceM=λ1+λ2,k是一个 [ 0.04 , 0.06 ] [0.04,0.06] [0.04,0.06]的常数,现在的判断逻辑为:

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

算法步骤

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

效果

原图
效果图