Harris角点检测

转载于 2018-05-01 10:14:33 发布 · 344 阅读

文章标签：

#Harris

计算机视觉专栏收录该内容

7 篇文章

订阅专栏

1、主要思想：
如果像素周围显示存在多于一个方向的边，我们认为该点为兴趣点。该点就称为角点。

我们把图像域中点 $x$ 上的对称半正定矩阵 $M_I=M_I(x)$ 定义为：
$M_I=\nabla I \nabla I^T = \left[ \begin{matrix} I_x \\ I_y \end{matrix} \right] \left[ \begin{matrix} I_x & I_y \end{matrix} \right]=\left[ \begin{matrix} I_x^2 & I_xI_y \\ I_xI_y & I_y^2 \end{matrix} \right]$

其中 $\nabla I$ 为包含导数 $I_x$ 和 $I_y$ 的图像梯度， $M_I$ 的秩为1，特征值为 $\lambda_1=|\nabla I|^2$ 和 $\lambda_2=0$
现在对于图像的每一个像素，我们可以计算出该矩阵
选择权重矩阵 $W$ （通常为高斯滤波器 $G_{\sigma}$ ），我们可以得到卷积
$\overline{M_I}=W*M_I$
该卷积的目的是得到 $M_I$ 在周围像素上的局部平均。计算出的矩阵 $\overline{M_I}$ 称为Harris矩阵。 $W$ 的宽度决定了在像素 $x$ 周围的感兴趣区域。
（1） $\lambda_1$ 和 $\lambda_2$ 都是很大的整数，则该 $x$ 点为角点
（2） $\lambda_1$ 很大， $\lambda_2 \approx 0$ ，则该区域存在一个边
（3）如果 $\lambda_1 \approx \lambda_2 \approx 0$ ，该区域内为空。

指示函数： $det(\overline{M_I})-k\;trace(\overline{M_I})^2=\lambda_1\lambda_2-k(\lambda_1+\lambda_2)^2$
为了去除加权常数 $k$ ，我们通常使用商数：
$\dfrac{det(\overline{M_I})}{trace(\overline{M_I})^2}$

通常，两个（相同大小）像素块 $I_1(x)$ 和 $I_2(x)$ 的相关矩阵定义为：
$c(I_1,I_2)=\sum\limits_xf(I_1(x),I_2(x))$
对于互相关矩阵，函数 $f(I_1,I_2)=I_1I_2$ ， $c(I_1,I_2)$ 的值越大，像素块 $I_1$ 和 $I_2$ 的相似度越高。

归一化的互相关矩阵是互相关矩阵的一种变形，可以定义为：
$ncc(I_1,I_2)=\dfrac{1}{n-1}\sum\limits_x\dfrac{(I_1(x)-\mu_1)}{\sigma_1} \cdot \dfrac{(I_2(x)-\mu_2)}{\sigma_2}$
该方法对图像亮度变化具有稳健性

2、实现步骤
（1）利用水平与竖直差分算子对图像进行卷积操作，计算得到相应的 $f_x$ 、 $f_y$ ，根据实对称矩阵 $M$ 的组成，计算对应元素的值。
（2）利用高斯函数对矩阵 $M$ 进行平滑操作，得到新的 $M$ 矩阵，步骤1和2可以改变顺序，也可以先对图像进行高斯滤波，再求相应方向上梯度大小。
（3）对每一像素和给定的邻域窗口，计算局部特征结果矩阵 $M$ 的特征值和响应函数 $C(i,j)=det(M)-k(trace(M))^2$ 。
（4）选取响应函数的阈值，根据非极大值抑制原理，同时满足阈值及某邻域内的局部极大值为候选角点。

#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
using namespace cv;
using namespace std;

void CornerHarris(const Mat& srcImage, Mat& result,
                  int blockSize, int kSize, double k)
{
    Mat src;
    srcImage.copyTo(src);
    result.create(src.size(), CV_32F);
    int depth = src.depth();
    // 检测掩膜尺寸
    double scale = (double)(1 << ((kSize > 0 ?
kSize : 3) - 1)) * blockSize;
    if (depth == CV_8U)
        scale *= 255.;
    scale = 1. / scale;
    // sobel滤波
    Mat dx, dy;
    Sobel(src, dx, CV_32F, 1, 0, kSize, scale, 0);
    Sobel(src, dy, CV_32F, 0, 1, kSize, scale, 0);
    Size size = src.size();
    cv::Mat cov(size, CV_32FC3);
    int i, j;
    // 求解水平与竖直梯度
    for (i = 0; i < size.height; i++){
        float *covData = (float*)(cov.data + i*cov.step);
        const float *dxData = (const float*)(dx.data + i*dx.step);
        const float *dyData = (const float*)(dy.data + i*dy.step);
        for (j = 0; j < size.width; j++)
        {
            float dx_ = dxData[j];
            float dy_ = dyData[j];
            covData[3 * j] = dx_*dx_;
            covData[3 * j + 1] = dx_*dy_;
            covData[3 * j + 2] = dy_*dy_;
        }
    }
    // 计算窗口内求和
    boxFilter(cov, cov, cov.depth(),
        Size(blockSize, blockSize), Point(-1, -1), false);
    // 判断图像连续性
    if (cov.isContinuous() && result.isContinuous())
    {
        size.width *= size.height;
        size.height = 1;
    }
    else
        size = result.size();
    // 计算响应函数
    for (i = 0; i < size.height; i++)
    {
        // 获取图像矩阵指针
        float *resultData = (float*)(result.data + i*result.step);
        const float *covData = (const float*)(cov.data + i*cov.step);
        for (j = 0; j < size.width; j++)
        {
            // 焦点响应生成
            float a = covData[3 * j];
            float b = covData[3 * j + 1];
            float c = covData[3 * j + 2];
            resultData[j] = a*c - b*b - k*(a + c)*(a + c);
        }
    }
}

int main()
{
    cv::Mat srcImage = cv::imread("sea.jpg");
    if (!srcImage.data)
        return -1;
    cv::imshow("srcImage", srcImage);
    cv::Mat srcGray, result;
    cvtColor(srcImage, srcGray, COLOR_BGR2GRAY);
    result = Mat::zeros(srcImage.size(), CV_32FC1);
    // 角点检测参数
    int blockSize = 2;
    int apertureSize = 3;
    double k = 0.04;
    // 角点检测
    // cornerHarris( srcGray, result, blockSize, apertureSize, k, BORDER_DEFAULT );
    CornerHarris(srcGray, result, blockSize, apertureSize, k);
    // 矩阵归一化
    normalize(result, result, 0, 255, NORM_MINMAX, CV_32FC1, Mat());
    convertScaleAbs(result, result);
    // 绘图角点检测结果
    for (int j = 0; j < result.rows; j++)
    {
        for (int i = 0; i < result.cols; i++)
        {
            if ((int)(result.at<uchar>(j, i)) > 150)
            {
                circle(srcImage, Point(i, j), 5, Scalar(0), 2, 8, 0);
            }
        }
    }
    cv::imshow("result", srcImage);
    cv::waitKey(0);
    return 0;
}

#include "opencv2/highgui.hpp"
#include "opencv2/imgproc.hpp"
#include <iostream>
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( int argc, char** argv )
{
  src = cv::imread("images/object.jpg");
  if ( src.empty() )
  {
    cout << "Could not open or find the image!\n" << endl;
    cout << "Usage: " << argv[0] << " <Input image>" << endl;
    return -1;
  }
  cvtColor( src, src_gray, COLOR_BGR2GRAY );
  namedWindow( source_window, WINDOW_AUTOSIZE );
  createTrackbar( "Threshold: ", source_window, &thresh, max_thresh, cornerHarris_demo );
  imshow( source_window, src );
  cornerHarris_demo( 0, 0 );
  waitKey(0);
  return(0);
}

void cornerHarris_demo( int, void* )
{
  Mat dst, dst_norm, dst_norm_scaled;
  dst = Mat::zeros( src.size(), CV_32FC1 );
  int blockSize = 2;
  int apertureSize = 3;
  double k = 0.04;
  cornerHarris( src_gray, dst, blockSize, apertureSize, k, BORDER_DEFAULT );
  normalize( dst, dst_norm, 0, 255, NORM_MINMAX, CV_32FC1, Mat() );
  convertScaleAbs( dst_norm, dst_norm_scaled );
  for( int j = 0; j < dst_norm.rows ; j++ )
     { for( int i = 0; i < dst_norm.cols; i++ )
          {
            if( (int) dst_norm.at<float>(j,i) > thresh )
              {
               circle( dst_norm_scaled, Point( i, j ), 5,  Scalar(0), 2, 8, 0 );
              }
          }
     }
  namedWindow( corners_window, WINDOW_AUTOSIZE );
  imshow( corners_window, dst_norm_scaled );
}

#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
using namespace cv;
using namespace std;
Mat heatmap;

enum MyShape{MyCIRCLE=0,MyRECTANGLE,MyELLIPSE};
struct ParamColorMap {
    int iColormap;
    Mat img;
};
String winName="False color";
static const String ColorMaps[] = { "Autumn", "Bone", "Jet", "Winter", "Rainbow", "Ocean", "Summer",
                                    "Spring", "Cool", "HSV", "Pink", "Hot", "Parula", "User defined (random)"};
static void TrackColorMap(int x, void *r)
{
    ParamColorMap *p = (ParamColorMap*)r;
    Mat dst;
    p->iColormap= x;
    if (x == COLORMAP_PARULA + 1)
    {
        Mat lutRND(256, 1, CV_8UC3);
        randu(lutRND, Scalar(0, 0, 0), Scalar(255, 255, 255));
        applyColorMap(p->img, dst, lutRND);
    }
    else
        applyColorMap(p->img,dst,p->iColormap);
    putText(dst, "Colormap : "+ColorMaps[p->iColormap], Point(10, 20), FONT_HERSHEY_SIMPLEX, 0.8, Scalar(255, 255, 255),2);
    imshow(winName, dst);
}

void CornerHarris(const Mat& srcImage, Mat& result,
                  int blockSize, int kSize, double k)
{
    Mat src;
    srcImage.copyTo(src);
    result.create(src.size(), CV_32F);
    int depth = src.depth();
    // 检测掩膜尺寸
    double scale = (double)(1 << ((kSize > 0 ?
kSize : 3) - 1)) * blockSize;
    if (depth == CV_8U)
        scale *= 255.;
    scale = 1. / scale;
    GaussianBlur(src, src, Size(blockSize, blockSize), 0, 0, BORDER_DEFAULT);
    // sobel滤波
    Mat dx, dy;
    Sobel(src, dx, CV_32F, 1, 0, kSize, scale, 0);
    Sobel(src, dy, CV_32F, 0, 1, kSize, scale, 0);
    Size size = src.size();
    cv::Mat cov(size, CV_32FC3);
    int i, j;
    // 求解水平与竖直梯度
    for (i = 0; i < size.height; i++){
        float *covData = (float*)(cov.data + i*cov.step);
        const float *dxData = (const float*)(dx.data + i*dx.step);
        const float *dyData = (const float*)(dy.data + i*dy.step);
        for (j = 0; j < size.width; j++)
        {
            float dx_ = dxData[j];
            float dy_ = dyData[j];
            covData[3 * j] = dx_*dx_;
            covData[3 * j + 1] = dx_*dy_;
            covData[3 * j + 2] = dy_*dy_;
        }
    }
    // 计算窗口内求和
    //boxFilter(cov, cov, cov.depth(), Size(blockSize, blockSize), Point(-1, -1), false);

    GaussianBlur(cov, cov, Size(blockSize, blockSize), 0, 0, BORDER_DEFAULT);
    // 判断图像连续性
    if (cov.isContinuous() && result.isContinuous())
    {
        size.width *= size.height;
        size.height = 1;
    }
    else
        size = result.size();

    for (i = 0; i < result.rows; i++)
    {
        // 获取图像矩阵指针
        float *resultData = (float*)(result.data + i*result.step);
        const float *covData = (const float*)(cov.data + i*cov.step);
        for (j = 0; j < result.cols; j++)
        {
            // 焦点响应生成
            float a = covData[3 * j];
            float b = covData[3 * j + 1];
            float c = covData[3 * j + 2];
            resultData[j] = a*c - b*b - k*(a + c)*(a + c);
        }
    }
    /*
    // 计算响应函数
    for (i = 0; i < size.height; i++)
    {
        // 获取图像矩阵指针
        float *resultData = (float*)(result.data + i*result.step);
        const float *covData = (const float*)(cov.data + i*cov.step);
        for (j = 0; j < size.width; j++)
        {
            // 焦点响应生成
            float a = covData[3 * j];
            float b = covData[3 * j + 1];
            float c = covData[3 * j + 2];
            resultData[j] = a*c - b*b - k*(a + c)*(a + c);
        }
    }*/
}

cv::Mat getHeatMap(cv::Mat input) // input is of type CV_8UC1, return is of type CV_8UC3
{
    cv::Mat result(input.rows, input.cols, CV_8UC3);
    int row = input.rows;
    int col = input.cols;
    for(int i=0;i<row;i++){
        for(int j=0;j<col;j++){
            int val = input.at<uchar>(i,j);
            result.at<cv::Vec3b>(i,j)[0]=val;
            result.at<cv::Vec3b>(i,j)[1]=val;
            result.at<cv::Vec3b>(i,j)[1]=255-val;
        }
    }
    return result;
}

int main()
{
    cv::Mat srcImage = cv::imread("images/hourse.PNG");
    if (!srcImage.data)
        return -1;
    cv::imshow("srcImage", srcImage);
    cv::Mat srcGray, result;
    cvtColor(srcImage, srcGray, COLOR_BGR2GRAY);
    result = Mat::zeros(srcImage.size(), CV_32FC1);
    // 角点检测参数
    int blockSize = 3;//2 3 4 5 6 7
    int apertureSize = 3;//3 5 7 9
    double k = 0.04;
    // 角点检测
    // cornerHarris( srcGray, result, blockSize, apertureSize, k, BORDER_DEFAULT );
    CornerHarris(srcGray, result, blockSize, apertureSize, k);
    //cornerHarris(srcGray, result, blockSize, apertureSize, k);
    // 矩阵归一化
    normalize(result, result, 0, 255, NORM_MINMAX, CV_32FC1, Mat());
    convertScaleAbs(result, result);

    Mat tmp = result.clone();
    // 绘图角点检测结果
    for (int j = 0; j < result.rows; j++)
    {
        for (int i = 0; i < result.cols; i++)
        {
            cout<< (int)(result.at<uchar>(j, i)) <<" ";
            if ((int)(result.at<uchar>(j, i)) > 100)
            {
                circle(srcImage, Point(i, j), 5, Scalar(0), 2, 8, 0);
            }
        }
        cout<<endl;

    }
    cv::imshow("result", srcImage);
    tmp.convertTo(tmp,CV_8UC1);
    heatmap = getHeatMap(result);
    imshow("heatmap",heatmap);
    ParamColorMap p;
    p.img = tmp;
    p.iColormap = 0;
    imshow("Gray image",tmp);
    namedWindow(winName);
    createTrackbar("colormap", winName,&p.iColormap,1,TrackColorMap,(void*)&p);
    setTrackbarMin("colormap", winName, COLORMAP_AUTUMN);
    setTrackbarMax("colormap", winName, COLORMAP_PARULA+1);
    setTrackbarPos("colormap", winName, -1);
    TrackColorMap(0, (void*)&p);
    cout << "Press a key to exit" << endl;

    cv::waitKey(0);
    return 0;
}
/*
cv::Mat getHeatMap(cv::Mat input) // input is of type CV_8UC1, return is of type CV_8UC3
{
    cv::Mat result(input.rows, input.cols, CV_8UC3);
    for (int yy = 0; yy < input.rows; ++yy)
    {
        for (int xx = 0; xx < input.cols; ++xx)
        {
            int pixelValue = input.at<uchar>(yy, xx);
            if (pixelValue < 128) {
                result.at<cv::Vec3b>(yy, xx) = cv::Vec3b(0, 0 + 2*pixelValue, 255 - 2 * pixelValue);
            } else {
                result.at<cv::Vec3b>(yy, xx) = cv::Vec3b(0 + 2*pixelValue, 255 - 2 * pixelValue, 0);
            }
        }
    }
    return result;
}
*/

#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
using namespace cv;
using namespace std;
Mat heatmap;

enum MyShape{MyCIRCLE=0,MyRECTANGLE,MyELLIPSE};
struct ParamColorMap {
    int iColormap;
    Mat img;
};

String winName="False color";
static const String ColorMaps[] = { "Autumn", "Bone", "Jet", "Winter", "Rainbow", "Ocean", "Summer",
                                    "Spring", "Cool", "HSV", "Pink", "Hot", "Parula", "User defined (random)"};

static void TrackColorMap(int x, void *r)
{
    ParamColorMap *p = (ParamColorMap*)r;
    Mat dst;
    p->iColormap= x;
    if (x == COLORMAP_PARULA + 1)
    {
        Mat lutRND(256, 1, CV_8UC3);
        randu(lutRND, Scalar(0, 0, 0), Scalar(255, 255, 255));
        applyColorMap(p->img, dst, lutRND);
    }
    else
        applyColorMap(p->img,dst,p->iColormap);
    putText(dst, "Colormap : "+ColorMaps[p->iColormap], Point(10, 20), FONT_HERSHEY_SIMPLEX, 0.8, Scalar(255, 255, 255),2);
    imshow(winName, dst);
}

void CornerHarris(const Mat& srcImage, Mat& result,
                  int blockSize, int kSize, double k)
{
    Mat src;
    srcImage.copyTo(src);
    result.create(src.size(), CV_32F);
    int depth = src.depth();
    // 检测掩膜尺寸
    double scale = (double)(1 << ((kSize > 0 ?
kSize : 3) - 1)) * blockSize;
    if (depth == CV_8U)
        scale *= 255.;
    scale = 1. / scale;
    GaussianBlur(src, src, Size(blockSize, blockSize), 0, 0, BORDER_DEFAULT);
    // sobel滤波
    Mat dx, dy;
    Sobel(src, dx, CV_32F, 1, 0, kSize, scale, 0);
    Sobel(src, dy, CV_32F, 0, 1, kSize, scale, 0);
    Size size = src.size();
    cv::Mat cov(size, CV_32FC3);
    int i, j;
    // 求解水平与竖直梯度
    for (i = 0; i < size.height; i++){
        float *covData = (float*)(cov.data + i*cov.step);
        const float *dxData = (const float*)(dx.data + i*dx.step);
        const float *dyData = (const float*)(dy.data + i*dy.step);
        for (j = 0; j < size.width; j++)
        {
            float dx_ = dxData[j];
            float dy_ = dyData[j];
            covData[3 * j] = dx_*dx_;
            covData[3 * j + 1] = dx_*dy_;
            covData[3 * j + 2] = dy_*dy_;
        }
    }
    // 计算窗口内求和
    //boxFilter(cov, cov, cov.depth(), Size(blockSize, blockSize), Point(-1, -1), false);

    GaussianBlur(cov, cov, Size(blockSize, blockSize), 0, 0, BORDER_DEFAULT);
    // 判断图像连续性
    if (cov.isContinuous() && result.isContinuous())
    {
        size.width *= size.height;
        size.height = 1;
    }
    else
        size = result.size();

    for (i = 0; i < result.rows; i++)
    {
        // 获取图像矩阵指针
        float *resultData = (float*)(result.data + i*result.step);
        const float *covData = (const float*)(cov.data + i*cov.step);
        for (j = 0; j < result.cols; j++)
        {
            // 焦点响应生成
            float a = covData[3 * j];
            float b = covData[3 * j + 1];
            float c = covData[3 * j + 2];
            resultData[j] = a*c - b*b - k*(a + c)*(a + c);
        }
    }
}

cv::Mat getHeatMap(cv::Mat input) // input is of type CV_8UC1, return is of type CV_8UC3
{
    cv::Mat result(input.rows, input.cols, CV_8UC3);
    int row = input.rows;
    int col = input.cols;
    for(int i=0;i<row;i++){
        for(int j=0;j<col;j++){
            int val = input.at<uchar>(i,j);
            result.at<cv::Vec3b>(i,j)[0]=val;
            result.at<cv::Vec3b>(i,j)[1]=val;
            result.at<cv::Vec3b>(i,j)[1]=255-val;
        }
    }
    return result;
}

int main()
{
    cv::Mat srcImage = cv::imread("images/object.jpg");
    if (!srcImage.data)
        return -1;
    cv::imshow("srcImage", srcImage);
    cv::Mat srcGray, result;
    cvtColor(srcImage, srcGray, COLOR_BGR2GRAY);
    result = Mat::zeros(srcImage.size(), CV_32FC1);
    // 角点检测参数
    int blockSize = 3;//2 3 4 5 6 7
    int apertureSize = 3;//3 5 7 9
    double k = 0.04;
    // 角点检测
    // cornerHarris( srcGray, result, blockSize, apertureSize, k, BORDER_DEFAULT );
    CornerHarris(srcGray, result, blockSize, apertureSize, k);
    //cornerHarris(srcGray, result, blockSize, apertureSize, k);
    // 矩阵归一化
    normalize(result, result, 0, 255, NORM_MINMAX, CV_32FC1, Mat());
    convertScaleAbs(result, result);

    Mat tmp = result.clone();
    // 绘图角点检测结果
    for (int j = 0; j < result.rows; j++)
    {
        for (int i = 0; i < result.cols; i++)
        {
            cout<< (int)(result.at<uchar>(j, i)) <<" ";
            if ((int)(result.at<uchar>(j, i)) > 120)
            {
                circle(srcImage, Point(i, j), 5, Scalar(0), 2, 8, 0);
            }
        }
        cout<<endl;

    }
    cv::imshow("result", srcImage);
    tmp.convertTo(tmp,CV_8UC1);
    ParamColorMap p;
    p.img = tmp;
    p.iColormap = 0;
    imshow("Gray image",tmp);
    namedWindow(winName);
    createTrackbar("colormap", winName,&p.iColormap,1,TrackColorMap,(void*)&p);
    setTrackbarMin("colormap", winName, COLORMAP_AUTUMN);
    setTrackbarMax("colormap", winName, COLORMAP_PARULA+1);
    setTrackbarPos("colormap", winName, -1);
    TrackColorMap(0, (void*)&p);
    cout << "Press a key to exit" << endl;

    cv::waitKey(0);
    return 0;
}