Harris角点检测

            Harris角点检测的原理，在网上有很多了，我根据它的原理实现了相应的检测方法。

            

/***********************************************
***File Name:main.cpp
***@Function:
***@Author: zhangjie
***@Mail:zhjkobe2013@live.com
***@Create time:2015-04-16
***@Modify time:2015-
***********************************************/

#include <cstdlib>
#include <iostream>
#include <cmath>
#include "FTImage.h"
const int Max_Corner=800;
using namespace std;


struct Point3
{
	int x;
	int y;
	double dval;
	Point3()
	{
		x=0;
		y=0;
		dval=0.0f;
	}
	void setval(int mx,int my,double r)
	{
         x=mx;
         y=my;
         dval=r;
    }
    
};
struct Feature
{
    double val1;
    double val2;
    double val3;
};
const int offset[16][2]=
{
	{0,-3},{1,-3},{2,-2},{3,-1},
	{3,0},{3,1},{2,2},{1,3},
	{0,3},{-1,3},{-2,2},{-3,1},
    {-3,0},{-3,-1},{-2,-2},{-1,-3}
};
Feature Compute_Feature(float Ixx,float Ixy,float Iyy)
{ 
	Feature res;
	float det=(Ixx-Iyy)*(Ixx-Iyy)+4*(Ixy*Ixy);
	det=sqrtf(det);
	//if(Ixx+Iyy==0) printf("error\n");
	res.val1=(Ixx+Iyy-det)/2;
	//if(res.val1<=0) printf("error\n");
	res.val2=(Ixx+Iyy+det)/2;
	res.val3=res.val1+res.val2;
	return res;
}
static Point3 Feature_Points[Max_Corner];
float GaussWeight(int m,int n)
{
   int ty=abs(m);
   int tx=abs(n);

   float sigma=1.0f;
   float res=exp(-(tx*tx+ty*ty)/(2*sigma*sigma));
   return res;
}
void DrawCircle(uchar* ptr,int x,int y,int val,int stride)
{
	for(int i=0;i<16;i++)
	{
		int ox=x+offset[i][0];
		int oy=y+offset[i][1];
	    ptr[oy*stride+ox*3]=255;
	    ptr[oy*stride+ox*3+1]=0;
	    ptr[oy*stride+ox*3+2]=0;
	}
}

void sobel_x(uchar* ptr,int row,int col,double* ptrx)
{
  const int sobelx[3][3]={{-1,0,1},{-2,0,2},{-1,0,1}};
  double sum;
  
  
  for(int i=0;i<row;i++)
   for(int j=0;j<col;j++)
   {  
     if(i<1||j<1||i>row-2||j>col-2)
           ptrx[i*col+j]=ptr[i*col+j];
     else  
     {
        sum=0.0f;
        for(int m=i-1;m<=i+1;m++)
         for(int n=j-1;n<=j+1;n++)
         {
           sum=sum+((int)ptr[m*col+n])*sobelx[m-i+1][n-j+1]*1.0;
         }   
         ptrx[i*col+j]=abs(sum)/255;
     }
     
      
   }
    
    
  /*const int sobelx[3]={-1,0,1};
  for(int i=0;i<row;i++)
   for(int j=0;j<col;j++)
   {  
     if(j<1||j>col-2)
           ptrx[i*col+j]=0.0f;
     else  
     {
           ptrx[i*col+j]=(double)(ptr[i*col+j+1]-ptr[i*col+j-1]);
     }
     
   }*/
}
void sobel_y(uchar* ptr,int row,int col,double* ptry)
{
   const int sobely[3][3]={{-1,-2,-1},{0,0,0},{1,2,1}};
   double sum;
   
   for(int i=0;i<row;i++)
    for(int j=0;j<col;j++)
    {  
      if(i<1||j<1||i>row-2||j>col-2)
           ptry[i*col+j]=ptr[i*col+j];
      else  
      {
        sum=0.0f;
        for(int m=i-1;m<=i+1;m++)
         for(int n=j-1;n<=j+1;n++)
         {
           sum=sum+((int)ptr[m*col+n])*sobely[m-i+1][n-j+1]*1.0;
         }  
         ptry[i*col+j]=abs(sum)/255; 
      }
     
    }
    
  /*const int sobely[3]={-1,0,1};
  for(int i=0;i<row;i++)
   for(int j=0;j<col;j++)
   {  
     if(i<1||i>row-2)
           ptry[i*col+j]=0.0f;
     else  
     {
           ptry[i*col+j]=(double)(ptr[(i+1)*col+j]-ptr[(i-1)*col+j]);
     }
     
   }*/
   
}
class Harris
{
public:
    Harris();
    //Harris(int nCorner,int thresh)
    ~Harris();
    bool Harris_Corner(uchar* ptr,int row,int col,uchar* rgbptr,int stride);
    //void Draw_Corner(uchar* ptr,int stride);
    //void Print_Corner();
private:
	double rate;
	int skipe_Pixel;
	bool Use_GaussWeight;
	double thresh;
	int radius;
};

Harris::Harris()
{
	rate=0.04;
	skipe_Pixel=5;
	Use_GaussWeight=true;
	thresh=50;
	radius=5;
	//Feature_Points=(Point3*)malloc(Max_Corner*sizeof(Point3));
}

Harris::~Harris()
{
	//if(Feature_Points) free(Feature_Points);
}

bool Harris::Harris_Corner(uchar* ptr,int row,int col,uchar* rgbptr,int stride)
{
    double* ptrx=(double*)malloc(row*col*sizeof(double));
    double* ptry=(double*)malloc(row*col*sizeof(double));
    int* Feature_Map=(int*)malloc(row*col*sizeof(int));
    memset((void*)Feature_Map,0,row*col*sizeof(int));
    
	sobel_x(ptr,row,col,ptrx);
	sobel_y(ptr,row,col,ptry);
    
    int count=0;
	for(int i=radius;i<row-radius;i+=skipe_Pixel)
	{
		for(int j=radius;j<col-radius;j+=skipe_Pixel)
		{
			double Ixx=0;
			double Ixy=0;
			double Iyy=0;

			for(int m=-radius;m<=radius;m++)
			   for(int n=-radius;n<=radius;n++)
			   {
			   	 Ixx+=ptrx[(i+m)*col+(j+n)]*ptrx[(i+m)*col+(j+n)];//*GaussWeight(m,n);
			   	 Ixy+=ptrx[(i+m)*col+(j+n)]*ptry[(i+m)*col+(j+n)];//*GaussWeight(m,n);
			   	 Iyy+=ptry[(i+m)*col+(j+n)]*ptry[(i+m)*col+(j+n)];//*GaussWeight(m,n);
			   }
            Feature Fval=Compute_Feature(Ixx,Ixy,Iyy);
            //double R=Fval.val1*Fval.val2-rate*(Fval.val3*Fval.val3);
            double r1=Ixx*Iyy-Ixy*Ixy;
            double r2=Ixx+Iyy;
            double R=r1-rate*r2*r2;
            if(R>thresh&&count<Max_Corner&&Feature_Map[i*col+j]==0) 
            { 
                   
            		count=count+1; 
            		Feature_Points[count].setval(j,i,R);
            		int x=Feature_Points[count].x;
                    int y=Feature_Points[count].y;
                    //double val=Feature_Points[count].dval;
           		    printf("count=%d,x=%d,y=%d,val=%f\n",count,x,y,R);
           		    DrawCircle(rgbptr,x,y,R,stride);
           		    
           		    for(int p=-radius;p<=radius;p++)
			          for(int q=-radius;q<=radius;q++)
			          {
                         Feature_Map[(i+p)*col+j+q]=1;
                      }
			   
            }

		}
	}
	
	/*for(int i=0;i<Max_Corner;i++)
	{
	  if(Feature_Points[i].dval>0.0) DrawCircle(rgbptr,Feature_Points[i],stride);
	  else
	    {printf("corner num is %d\n",i); break;}	
	}
	
	for(int i=0;i<Max_Corner;i++)
	{
      printf("i=%d,x=%d,y=%d,R=%f\n",i,Feature_Points[i].x,Feature_Points[i].y,Feature_Points[i].dval);
	  if(Feature_Points[i].dval==0.0f) break;
	  else
		cout<<i<<":"<<"["<<Feature_Points[i].x<<","<<Feature_Points[i].y<<"]"<<endl;
	}*/

   free(ptrx);
   free(ptry);
   free(Feature_Map);
   return true;
} 


/*void Harris::Draw_Corner(uchar* ptr,int stride)
{
	

}


void Harris::Print_Corner()
{
	
}*/




int main(int argc, char *argv[])
{
    printf("start\n");
    FTImage ImgGray,ImgRGB;
    ImgGray.FTLoadImageGray("church01.bmp",480,740);
    ImgRGB.FTLoadImageRGB("church01.bmp",480,740);
    //Harris* m_Harris=new Harris();
     Harris m_Harris;
    m_Harris.Harris_Corner(ImgGray.pixels,740,480,ImgRGB.pixels,ImgRGB.GetWidthStep());
    //m_Harris.Draw_Corner();
    //m_Harris.Print_Corner();
    
    printf("end\n");
    ImgRGB.FTSaveImage("res.bmp");
    //if(m_Harris) delete m_Harris;
    getchar();
    system("PAUSE");
    return EXIT_SUCCESS;
    
  
}