openCV2馬拉松第19圈——Harris角點檢測（自己實現）

電腦視覺討論群162501053轉載請註明：http://blog.csdn.net/abcd1992719g/article/details/26824529


收入囊中

• 使用OpenCV的connerHarris實現角點檢測
• 自己實現Harris演算法

下面是自己實現的一個
因為閥值設定比較高，所以房屋周圍沒有找出來


葵花寶典在此之前，我們講過邊緣的檢測，邊緣檢測的基本原理就是x方向或者y方向梯度變化很大，角點，顧名思義，就是兩個方向的梯度變化都很大。
左1，平滑地區，沒有邊緣和角點，視窗在任何方向移動都無變化左2，邊緣地區，在邊緣方向移動沒有變化左3，角點地區，在任何方向移動都有顯著的變化
下面我們定義

E(u,v)=∑(x,y)∈W [I(x+u,y+v)−I(x,y)]2（w是我們的視窗,[u,v]是我們的shift，也就是移動）

我們想瞭解微小移動對E到底有何影響

於是我們對I進行一階泰勒展開

I(x+u,y+v)=I(x,y)+Ixu+Iyv+higherorder terms

≈I(x,y)+Ixu+Iyv 

=I(x,y)+[Ix Iy ][u v]T

再代入下面的公式

E(u,v)=∑[I(x+u,y+v)−I(x,y)]2(x,y)∈W 

我們得到

因此，我們可以寫成

其中，M是一個二階矩矩陣，可以由我們原始圖片的差分得到

我們來想像一下，假使梯度是沿著x或者y方向的，也就是說Ix = 0或者是Iy = 0,那麼我們的M就是
如果a或者b很接近0，也就是說他們很小，那麼這就不是個角點，角點的地方a,b肯定都很大
我們發現，E是一個關於U,V的二次曲線，令E(u,v) = CONST,那麼E就是一個橢圓而M,由矩陣對角化可以寫成左邊的格式，我們可以得到兩個特徵值
長軸短軸由特徵值決定，橢圓方向由矩陣R決定，下面沒有用到R,所以你可以忽略它


在這裡，我們的特徵值就派上用場了，當兩個特徵值都很大很相近，說明我的橢圓很小，整個E變化很劇烈，那麼我就找到了我的角點。下面是 Harris給出的公式

• det(M) = 
• trace(M) = 

那麼很多人就在糾結要怎麼計算了。

det = Ix2*Iy2 - Ixy^2trace = Ix2 + Iy2

［還有另一種經驗公式，，我自己的實現就是採用了這一種］
說了這麼多理論，實現起來其實沒有那麼困難，下面我們就來看看演算法的步驟吧。

1. 用sobel運算元分別計算出水平梯度和垂直梯度

-1 0 1

-2 0 2

-1 0 1

-1 -2 -1 

0  0  0

1  2  1

2. 計算高斯二階矩陣，就是一個視窗裡所有差分和

3.計算響應函數R

4.設定閥值

5.非極大值抑制，判斷這個點的響應函數R是不是周圍最大的

強調一下，harris對仿射變換隻有部分不變性質。平移和旋轉具有covariant,但是scaling不具有，如

初識API

C++: void cornerHarris(InputArray  src, OutputArray  dst, int  blockSize, int  ksize, double  k, int  borderType=BORDER_DEFAULT )

	src – 8位元或者是浮點數矩陣. dst – 存放harris響應函數，類型是CV_32FC1因為計算出的值比較大，size和src一樣 . blockSize – 視窗大小. ksize – Sobel運算元大小. k – 一般是0.04-0.06. borderType  .

對每個像素 計算  共變數矩陣  over a  neighborhood. 用的是harris響應函數

荷槍實彈下面是官方sample,我就不多解釋了，純粹是API的調用，要注意的就是得到響應函數後進行了歸一化到［0，255］

#include "opencv2/highgui/highgui.hpp"#include "opencv2/imgproc/imgproc.hpp"#include <iostream>#include <stdio.h>#include <stdlib.h>using namespace cv;using namespace std;/// Global variablesMat src, src_gray;int thresh = 200;int max_thresh = 255;char* source_window = "Source image";char* corners_window = "Corners detected";/// Function headervoid cornerHarris_demo( int, void* );/** @function main */int main( int argc, char** argv ){  /// Load source image and convert it to gray  src = imread( argv[1], 1 );  cvtColor( src, src_gray, CV_BGR2GRAY );  /// Create a window and a trackbar  namedWindow( source_window, CV_WINDOW_AUTOSIZE );  createTrackbar( "Threshold: ", source_window, &thresh, max_thresh, cornerHarris_demo );  imshow( source_window, src );  cornerHarris_demo( 0, 0 );  waitKey(0);  return(0);}/** @function cornerHarris_demo */void cornerHarris_demo( int, void* ){  Mat dst, dst_norm, dst_norm_scaled;  dst = Mat::zeros( src.size(), CV_32FC1 );  /// Detector parameters  int blockSize = 2;  int apertureSize = 3;  double k = 0.04;  /// Detecting corners  cornerHarris( src_gray, dst, blockSize, apertureSize, k, BORDER_DEFAULT );  /// Normalizing  normalize( dst, dst_norm, 0, 255, NORM_MINMAX, CV_32FC1, Mat() );  convertScaleAbs( dst_norm, dst_norm_scaled );  /// Drawing a circle around corners  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 );              }          }     }  /// Showing the result  namedWindow( corners_window, CV_WINDOW_AUTOSIZE );  imshow( corners_window, dst_norm_scaled );}

舉一反三下面是我自己的實現，我找了好多地方都沒找到別人的原始碼...

#include "opencv2/highgui/highgui.hpp"  #include "opencv2/imgproc/imgproc.hpp"  #include <cmath>#include <iostream>using namespace cv;   Mat harris(Mat &im, double sigma, int thresh, int radius){    Mat dx,dy,Ix,Iy,Ix2,Iy2,Ixy,cim;    Sobel( im, Ix, CV_64F, 1, 0, 3);   //演算法第一步，計算水平垂直差分    Sobel( im, Iy, CV_64F, 0, 1, 3);    int ksize = max(1, (int)(6*sigma));    if(ksize % 2 == 0)    ksize++;GaussianBlur(Ix.mul(Ix), Ix2, Size(ksize, ksize), sigma);   //演算法第二步，計算二階高斯差分矩陣   GaussianBlur(Iy.mul(Iy), Iy2, Size(ksize, ksize), sigma);GaussianBlur(Ix.mul(Iy), Ixy, Size(ksize, ksize), sigma);//Harris corner measure//cim = (Ix2.*Iy2 - Ixy.^2)./(Ix2 + Iy2);cim = (Ix2.mul(Iy2) - Ixy.mul(Ixy)) / (Ix2+Iy2);            //演算法第三步，計算響應函數，我使用了另外一種    Mat structedElement(radius, radius, CV_8U, Scalar(1));    Mat mx,norm_cim;    normalize( cim, norm_cim, 0, 255, NORM_MINMAX, CV_8U, Mat() );dilate(norm_cim, mx, structedElement);norm_cim = ( norm_cim == mx) & (norm_cim>thresh);           //演算法第4第5步融合，非極大值抑制和閥值檢測return norm_cim;}int main( int, char** argv )  {  Mat src,gray;    src = imread( argv[1] );    cvtColor( src, gray, CV_RGB2GRAY );    Mat corners = harris(gray, 1.5, 30, 2);for( int j = 0; j < corners.rows ; j++ ) { for( int i = 0; i < corners.cols; i++ ) {            if( corners.at<unsigned char>(j,i) > 0){circle( gray, Point( i, j ), 3,  Scalar(0), 2, 8, 0 );}}}        namedWindow("result", 1);     imshow("result", gray);    waitKey();      return 0;  }

另外附加福利harris matlab版本 http://www.cs.illinois.edu/~slazebni/spring13/harris.mharris python版本 http://www.janeriksolem.net/2009/01/harris-corner-detector-in-python.html