角点检测（2）harris算子

#include "stdafx.h"
#include <iostream>
#include <stdlib.h>
#include <cv.h>
#include <cxcore.h>
#include <highgui.h>
#include <math.h>
using namespace cv;  
using namespace std;
int main(int argc, char** argv)
{
cout << "Corner Detection OpenCV!"<<endl;
char* filename="2.jpg";
IplImage* imgRGB = cvLoadImage("2.jpg");
IplImage* imgRGB2 = cvLoadImage("1.jpg");
IplImage* imgGrey = cvLoadImage("1.jpg",CV_LOAD_IMAGE_GRAYSCALE);

if (imgGrey==NULL){//image validation
cout << "No valid image input."<<endl;
char c=getchar();
return 1;
}
int w=imgGrey->width;
int h=imgGrey->height;

IplImage* eig_image = cvCreateImage(cvSize(w, h),IPL_DEPTH_32F, 1);
IplImage* temp_image = cvCreateImage(cvSize(w, h),IPL_DEPTH_32F, 1);

const int MAX_CORNERS = 140;//estimate a corner number
CvPoint2D32f corners[MAX_CORNERS] = {0};// coordinates of corners
//CvPoint2D32f* corners = new CvPoint2D32f[ MAX_CORNERS ]; //another method of declaring an array
int corner_count = MAX_CORNERS;
double quality_level = 0.1;//threshold for the eigenvalues
double min_distance = 5;//minimum distance between two corners
int eig_block_size = 3;//window size
int use_harris = false;//use 'harris method' or not

//----------initial guess by cvGoodFeaturesToTrack---------------
cvGoodFeaturesToTrack(imgGrey,
eig_image, // output
temp_image,
corners,
&corner_count,
quality_level,
min_distance,
NULL,
eig_block_size,
use_harris);


int r=2; //rectangle size
int lineWidth=1; // rectangle line width
//-----draw good feature corners on the original RGB image---------
for (int i=0;i<corner_count;i++){
cvRectangle(imgRGB2, cvPoint(corners[i].x-r,corners[i].y-r),
cvPoint(corners[i].x+r,corners[i].y+r), cvScalar(255,0,0),lineWidth);
}

int half_win_size=3;//the window size will be 3+1+3=7
int iteration=20;
double epislon=0.1;
cvFindCornerSubPix(
imgGrey,
corners,
corner_count,
cvSize(half_win_size,half_win_size),
cvSize(-1,-1),//no ignoring the neighbours of the center corner
cvTermCriteria(CV_TERMCRIT_ITER|CV_TERMCRIT_EPS,iteration,epislon)
);

//------draw subpix corners on another original RGB image------------
for (int i=0;i<corner_count;i++){
cvRectangle(imgRGB, cvPoint(corners[i].x-r,corners[i].y-r),
cvPoint(corners[i].x+r,corners[i].y+r), cvScalar(0,0,255),lineWidth);
}

//to display a coordinate of the third corner
cout<<"x="<<corners[2].x;
cout<<",y="<<corners[2].y<<endl;

cvNamedWindow("cvFindCornerSubPix", CV_WINDOW_AUTOSIZE );
cvShowImage( "cvFindCornerSubPix", imgRGB );
cvNamedWindow("cvGoodFeaturesToTrack", CV_WINDOW_AUTOSIZE );
cvShowImage( "cvGoodFeaturesToTrack", imgRGB2 );


cvWaitKey(0);
cvReleaseImage(&imgGrey);
cvReleaseImage(&imgRGB);
cvReleaseImage(&imgRGB2);
cvDestroyWindow("cvGoodFeaturesToTrack");
cvDestroyWindow("cvFindCornerSubPix");

//char c=getchar();
return 0;
}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%   Harris角点提取算法                                          % 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
clear; 
filename='1.jpg'; 
X= imread(filename);     % 读取图像 
Info=imfinfo(filename); 
if Info.BitDepth>8 
    f=rgb2gray(X); 
end 
% 
% fx = [5 0 -5;8 0 -8;5 0 -5];          % 高斯函数一阶微分，x方向(用于改进的Harris角点提取算法) 
ori_im=double(f)/255;                   %unit8转化为64为双精度double64 
fx = [-2 -1 0 1 2];                     % x方向梯度算子(用于Harris角点提取算法) 
Ix = filter2(fx,ori_im);                % x方向滤波 
% fy = [5 8 5;0 0 0;-5 -8 -5];          % 高斯函数一阶微分，y方向(用于改进的Harris角点提取算法) 
fy = [-2;-1;0;1;2];                     % y方向梯度算子(用于Harris角点提取算法) 
Iy = filter2(fy,ori_im);                % y方向滤波 
Ix2 = Ix.^2; 
Iy2 = Iy.^2; 
Ixy = Ix.*Iy; 
clear Ix; 
clear Iy; 
 
h= fspecial('gaussian',[7 7],2);        % 产生7*7的高斯窗函数，sigma=2 
 
Ix2 = filter2(h,Ix2); 
Iy2 = filter2(h,Iy2); 
Ixy = filter2(h,Ixy); 
 
height = size(ori_im,1); 
width = size(ori_im,2); 
result = zeros(height,width);           % 纪录角点位置，角点处值为1 
 
R = zeros(height,width); 
 
Rmax = 0;                              % 图像中最大的R值 
for i = 1:height 
    for j = 1:width 
        M = [Ix2(i,j) Ixy(i,j);Ixy(i,j) Iy2(i,j)];             % auto correlation matrix 
        R(i,j) = det(M)-0.06*(trace(M))^2;                     % 计算R 
        if R(i,j) > Rmax 
            Rmax = R(i,j); 
        end; 
    end; 
end; 
 
cnt = 0; 
for i = 2:height-1 
    for j = 2:width-1 
        % 进行非极大抑制，窗口大小3*3 
        if R(i,j) > 0.01*Rmax && R(i,j) > R(i-1,j-1) && R(i,j) > R(i-1,j) && R(i,j) > R(i-1,j+1) && R(i,j) > R(i,j-1) && R(i,j) > R(i,j+1) && R(i,j) > R(i+1,j-1) && R(i,j) > R(i+1,j) && R(i,j) > R(i+1,j+1) 
            result(i,j) = 1; 
            cnt = cnt+1; 
        end; 
    end; 
end; 
[posc, posr] = find(result == 1); 
cnt                                      % 角点个数 
imshow(ori_im) 
hold on; 
plot(posr,posc,'r+');