Harris角点检测

opencv官网教程
参考１
参考２

理论

在图像中取一个滑动窗口$W_{(xy)}$，不断移动窗口，观察窗口内的像素变化情况：

• 任何一个方向上变化都很大：　角点
• 一个方向变化大，另一个方向没有变化：边缘
• 所有方向都没什么变化：平坦区域

数学表达：

1. 对于灰度图像，滑动窗口，计算窗口的像素变化：
   $$E(u,v)=\sum _{x,y}{w(x,y)[I(x+u,y+v)-I(x,y)]}^2$$
   $w(x,y)$：窗口的在(x,y)的位置
   $I(x,y)$：在(x,y) 的像素值
   
   

2. 对 $$\sum _{x,y}{[I(x+u,y+v)-I(x,y)]}^2$$进行泰勒级数展开：
   $$\approx \sum _{x,y}[I(x,y)+u{I}_x+v{I}_y-I(x,y){]}^2$$
   展开多项式：
   $$\approx \sum _{x,y}{u}^2{I}_x^2+2uv{I}_x{I}_y+v^2{I}_y^2$$
   
   

3. 对$E(u,v)$写成矩阵的形式：$$E(u,v)\approx \left[\begin{array}{c}u,v\end{array}\right]\left(\sum _{x,y}w(x,y)\left[\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right]\right)\left[\begin{array}{c}u\\ v\end{array}\right]$$

4. 我们让中间的部分记作$M$(协方差矩阵)，则有：
   $$E(u,v)=\left[\begin{array}{c}u,v\end{array}\right]M\left[\begin{array}{c}u\\ v\end{array}\right]$$

5. 我们规定$R$代表每个窗口的得分，根据$R$来判断窗口内是否有角点：$$R=det(M)-k(Trace(M){)}^2$$
   $det(M)={\lambda }_1{\lambda }_2$
   $Trace(M)={\lambda }_1+{\lambda }_2$
   
   

6. 结论：
   $R$ > threashold ：这是一个角点
   $R$ < 0 并且　数值较大：这是边缘
   $|R|=$：小值：这是平坦区域

---

代码

#include <iostream>
#include <opencv2/opencv.hpp>


using namespace cv;
using namespace std;


Mat src, src_gray;
int thresh = 200;
int max_thresh = 255;

const char* source_window = "Source Image";
const char* corners_window = "Corners detected";

void cornerHarris_Demo(int, void*);


int main(void)
{

    src = imread("../res/image.png");
    if(src.empty())
    {
        cout << "can't load the image, please confirm the path" << endl;
        return -1;
    }


    cvtColor(src, src_gray, COLOR_BGR2GRAY);
    namedWindow(source_window);
    createTrackbar("Thresh", source_window, &thresh, max_thresh, cornerHarris_Demo);
    imshow(source_window ,src);
    cornerHarris_Demo(0, 0);

    waitKey();
    return 0;
}

void cornerHarris_Demo(int, void*)
{
    int blockSize = 2;
    int apertureSize = 3;
    double k = 0.04;

    Mat dst = Mat::zeros(src.size(), CV_32FC1);
    cornerHarris(src_gray, dst, blockSize, apertureSize, k);

    Mat dst_norm, dst_norm_scaled;
    normalize(dst, dst_norm, 0, 255, NORM_MINMAX, CV_32FC1, Mat()); // normalize dst to dst_norm in the range of 0-255
    convertScaleAbs(dst_norm, dst_norm_scaled);   // dst(x,y) = saturate_cast<uchar>(src(x,y)*alpha+beta),defult:  alpha:1  beta:0

    for(int i=0; i<dst_norm.rows; i++)
    {
        for(int j=0; j<dst_norm.cols; j++)
        {
            if((int)dst_norm.at<float>(i,j) > thresh)
            {
                circle(dst_norm_scaled, Point(i, j), 5, Scalar(0), 2); 
                circle(src, Point(i, j), 5, Scalar(255,0,0), 2);
            }
        }
    }

    namedWindow(corners_window);
    imshow(corners_window, dst_norm_scaled);
    imshow(source_window ,src);
}

结果：

---

OpenCV API

1. 在图像上进行Harris角点检测，对于每个像素$(x,y)$ ,计算2x2的梯度的协方差矩阵$M$ over blockSize x blockSize neighborhood, 然后计算:$$dst(x,y)=det(M)-k(Trace(M){)}^2$$

cornerHarris
(
InputArray src,     // 单通道, 8bits 或者 float 图
OutputArray dst,     // type CV_32FC1 and the same size as src，
int blockSize,     //理论里面的窗口$w$的大小
int ksize,     // 计算梯度用的Sobel核的大小
double k,     // 上面公式的k ，一般取0.04-0.06
int borderType = BORDER_DEFAULT
)