harris1.cpp

VC写的harris角点检测程序。又名plessey

// Harris1.cpp: implementation of the CHarris class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "ImageCodes.h"
#include "Harris1.h"
#include "Afxtempl.h"
#include <math.h>
#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif
extern UINT r0,c0;
extern unsigned int **GrayValue;//存放灰度值，8位 GrayValue[r][c]表示倒数r行c列的灰度值
CArray<CPoint,CPoint&> MaxpointArray;

//CArray<CPoint,CPoint&>MaxpointArray;
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CHarris::CHarris()
{

}

CHarris::~CHarris()
{

}

void CHarris::Harris(double sigma,int thresh,int radius)//从view中调用 CHarris m_harris; m_harris.Harris(0.7,200,3)
{
	// 指向x方向导数的指针
	double * pnGradX ; 
	// 指向y方向导数的指针
	double * pnGradY ;
	double * pnGradXY ;
	double * pMaxValue;
	int y ,i,j;
	int x ; int m_maxvalue;
	CPoint maxpoint;//角点

	MaxpointArray.RemoveAll();
	pnGradX      = new double [r0*c0]         ;
	pnGradY      = new double [r0*c0]         ;
	pnGradXY     = new double [r0*c0]         ;
	pMaxValue    = new double [r0*c0]         ;
	// 计算方向导数
	DirGrad2(pnGradX, pnGradY,pnGradXY) ;
	
	// 对gx gy gxgy 进行高斯滤波
	Gaussian(sigma, pnGradX,pnGradY,pnGradXY) ;
	//求R=(AB-C^2)-k(A+B)^2  cim = (Ix2.*Iy2 - Ixy.^2)./(Ix2 + Iy2 + eps); 
	for(y=0; y<(int)r0; y++)
		for(x=0; x<(int)c0; x++)
		{
		//	pMaxValue[y*c0+x]=((pnGradX[y*c0+x]*pnGradY[y*c0+x]-pnGradXY[y*c0+x]*pnGradXY[y*c0+x]))-0.04*(pnGradX[y*c0+x]+pnGradY[y*c0+x])*(pnGradX[y*c0+x]+pnGradY[y*c0+x]);
			pMaxValue[y*c0+x]=(pnGradX[y*c0+x]*pnGradY[y*c0+x]-pnGradXY[y*c0+x]*pnGradXY[y*c0+x])/(pnGradX[y*c0+x]+pnGradY[y*c0+x]+0.00000001);
		}
	BOOL bfindPoint;//求局部最大值 在小窗口内扫描如果 大于设定阈值 则将该点坐标存入MaxpointArray
	for(y=radius;y<(int)r0-2*radius;y++)
	{
		for(x=radius;x<(int)c0-2*radius;x++)
		{
			bfindPoint=FALSE;
			m_maxvalue=thresh;
			for(i=-radius;i<=radius;i++)
				for(j=-radius;j<=radius;j++)
				{
					if(pMaxValue[(y+i)*c0+x+j]>m_maxvalue)
					{						
						m_maxvalue=(int)pMaxValue[(y+i)*c0+x+j];
						maxpoint.x=x;
						maxpoint.y=r0-y;
						bfindPoint=TRUE;
					}
				}
			if(bfindPoint) MaxpointArray.Add(maxpoint);//存入角点坐标 在view中显示
			x+=2*radius;
		}
		y+=2*radius;

	}
	
	// 释放内存
	delete []pnGradX     ;
	delete []pnGradY     ;
	delete []pnGradXY    ;
	
}

void CHarris::DirGrad2(double *pnGradX , double *pnGradY,double *pnGradXY)
{
	int y ;
	int x ;
	// 计算x方向的方向导数，在边界出进行了处理，防止要访问的象素出界
	for(y=0; y<(int)r0; y++)
	{
		for(x=0; x<(int)c0; x++)
		{//求梯度
			pnGradX[y*c0+x]=(int)(GrayValue[y][min((int)(c0-1),x+1)]-GrayValue[y][max(0,x-1)])
				+(int)(GrayValue[max(0,y-1)][min((int)(c0-1),x+1)]-GrayValue[max(0,y-1)][max(0,x-1)])
				+(int)(GrayValue[min((int)(r0-1),y+1)][min((int)(c0-1),x+1)]-GrayValue[min((int)(r0-1),y+1)][max(0,x-1)]);
		}
	}

	// 计算y方向的方向导数，在边界出进行了处理，防止要访问的象素出界
	for(y=0; y<(int)r0; y++)
	{
		for(x=0; x<(int)c0; x++)
		{
			pnGradY[y*c0+x]=(int)(GrayValue[min((int)(r0-1),y+1)][x]-GrayValue[max(0,y-1)][x])
				+(int)(GrayValue[min((int)(r0-1),y+1)][max(0,x-1)]-GrayValue[max(0,y-1)][max(0,x-1)])
				+(int)(GrayValue[min((int)(r0-1),y+1)][min((int)c0-1,x+1)]-GrayValue[max(0,y-1)][min((int)c0-1,x+1)]);
				//(pUnOrgImg[min((int)(r0-1),y+1)*c0 + x]-pUnOrgImg[max(0,y-1)*c0+x]);
		}
	}
	for(y=0; y<(int)r0; y++)
		for(x=0; x<(int)c0; x++)
			pnGradXY[y*c0+x]=pnGradX[y*c0+x]*pnGradY[y*c0+x];
	for(y=0; y<(int)r0; y++)
		for(x=0; x<(int)c0; x++)
		{
			pnGradX[y*c0+x]=pnGradX[y*c0+x]*pnGradX[y*c0+x];
			pnGradY[y*c0+x]=pnGradY[y*c0+x]*pnGradY[y*c0+x];	
		}
}

void CHarris::Gaussian(double sigma, double *pUnchSmooth1,double *pUnchSmooth2,double *pUnchSmooth3)
{
	UINT y;
	UINT x;
	
	int i,j;
	
	// 高斯滤波器的数组长度
	
	int nWindowSize;
		
	// 数组长度，根据概率论的知识，选取[-3*sigma, 3*sigma]以内的数据。
	// 这些数据会覆盖绝大部分的滤波系数
	nWindowSize = 1 + (int)(2 * ceil(3 * sigma));
	
	//  窗口长度的1/2
	int	nHalfLen;   
	
	// 二维高斯数据滤波器
	double **pdKernel;

	
	// 高斯系数与图象数据的点乘
	double  dDotMul1, dDotMul2,dDotMul3;
	
	// 高斯滤波系数的总和
	double  dWeightSum  ;          
  
	// 中间变量
	double * pdTmp1,* pdTmp2,* pdTmp3 ;
	
	// 分配内存
	pdTmp1 = new double[r0*c0];
	pdTmp2 = new double[r0*c0];
	pdTmp3 = new double[r0*c0];
	pdKernel=new double *[nWindowSize];//动态给数组申请空间
	for(i=0;i<nWindowSize;i++)
		pdKernel[i]=new double[nWindowSize];

	// 产生二维高斯数据滤波器
	MakeGauss2(sigma, pdKernel, nWindowSize) ;
	// MakeGauss返回窗口的长度，利用此变量计算窗口的半长
	nHalfLen = nWindowSize / 2;
	
  // 进行高斯滤波
	for(y=0; y<r0; y++)
	{
		for(x=0; x<c0; x++)
		{
			dDotMul1=dDotMul2=dDotMul3= 0;
			dWeightSum = 0;
			for(j=(-nHalfLen); j<=nHalfLen; j++)
			{
				for(i=(-nHalfLen); i<=nHalfLen; i++)
				{
				// 判断是否在图象内部
					if( (i+x) >= 0  && (i+x) < c0 && (j+y) >= 0  && (j+y) < r0 )
					{
						
						dDotMul1 += pUnchSmooth1[(y+j)*c0+x+i]*pdKernel[nHalfLen+i][nHalfLen+j];
						dDotMul2 += pUnchSmooth2[(y+j)*c0+x+i]*pdKernel[nHalfLen+i][nHalfLen+j];
						dDotMul3 += pUnchSmooth3[(y+j)*c0+x+i]*pdKernel[nHalfLen+i][nHalfLen+j];
						dWeightSum += pdKernel[nHalfLen+i][nHalfLen+j];
					}
				}
			}
			pUnchSmooth1[y*c0 + x] = dDotMul1/dWeightSum;
			pUnchSmooth2[y*c0 + x] = dDotMul2/dWeightSum;
			pUnchSmooth3[y*c0 + x] = dDotMul3/dWeightSum;
		}
	}

	for(i=0;i<nWindowSize;i++)//释放空间
	delete []pdKernel[i];
	delete []pdKernel;
	delete []pdTmp1;
	delete []pdTmp2;
	delete []pdTmp3;
}

void CHarris::MakeGauss2(double sigma, double **pdKernel, int pnWindowSize)
{//二维高斯滤波器
	// 循环控制变量
	int i,j ;
	
	// 数组的中心点
	int nCenter;

	// 数组的某一点到中心点的距离
	double  dDis,dDiy  ; 

	double PI = 3.14159;
	// 中间变量
	double  dValue; 
	double  dSum  ;
	dSum = 0 ; 
	// 中心
	nCenter = (pnWindowSize) / 2;
	for(i=0; i<(pnWindowSize); i++)
	{
		for(j=0; j<(pnWindowSize); j++)
		{
			dDis = (double)(i - nCenter);
			dDiy=  (double)(j - nCenter);
			dValue=exp(-(1/2)*(dDis*dDis+dDiy*dDiy)/(sigma*sigma))/(sqrt(2 * PI) * sigma );
			pdKernel[i][j]=dValue;
			dSum += dValue;
		}
	}
	for(i=0; i<(pnWindowSize); i++)//归一化
	{
		for(j=0; j<(pnWindowSize); j++)
		{
			pdKernel[i][j]/= dSum;
		}
	}
}