www.pudn.com > Face3DModel.zip > moraveccornerdetector.cpp
// MoravecCornerDetector.cpp
//
//////////////////////////////////////////////////////////////////////
#include "stdafx.h"
#include "global.h"
#include "memory.h"
#include "math.h"
#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif
#define GAUSS_FUNC_VAL 100
unsigned char GetGray(unsigned char x1,unsigned char x2,unsigned char x3)
{
return (float)x1*0.299+x2*0.587+x3*0.114;
}
//
void CornerDetectOnGrayImage(int nWidth, int nHeight, BYTE *pbGray, BYTE *pbEdge)
{
int row, col;
int i;
double Gauss[GAUSS_FUNC_VAL];
//use notation in paper 'A Combined Corner and Edge Detector',
//Harrris,C. and Stephens,M.(1988).
int *piPartialDerivativeX = new int[nWidth * nHeight];
int *piPartialDerivativeY = new int[nWidth * nHeight];
float *piXX = new float[nWidth * nHeight];
float *piYY = new float[nWidth * nHeight];
float *piXY = new float[nWidth * nHeight];
double *pdR = new double[nWidth * nHeight];
double A,B,C;
double RMax;
PartialDerivativeX(nWidth, nHeight, pbGray, piPartialDerivativeX);
PartialDerivativeY(nWidth, nHeight, pbGray, piPartialDerivativeY);
ScalarMultiple(nWidth, nHeight, piPartialDerivativeX,
piPartialDerivativeX, piXX);
ScalarMultiple(nWidth, nHeight, piPartialDerivativeY,
piPartialDerivativeY, piYY);
ScalarMultiple(nWidth, nHeight, piPartialDerivativeX,
piPartialDerivativeY, piXY);
// Parameters can be adjusted
int nStandardDev = 2;
double k = 0.06; //0.06
double theta = 0.00001; //0.000005
for (i = 0; i < GAUSS_FUNC_VAL; i++)
{
Gauss[i] = gaussian(i, nStandardDev);
}
//start from the top row
for (row = nHeight - 1 - 1; row > 0; row--)
{
//start from the left column
for (col = 1; col < nWidth - 1; col++)
{
A=ConvolutionProductWithW(nWidth, nHeight, row, col, piXX, Gauss);
B=ConvolutionProductWithW(nWidth, nHeight, row, col, piYY, Gauss);
C=ConvolutionProductWithW(nWidth, nHeight, row, col, piXY, Gauss);
pdR[row*nWidth+col] = A*B - C*C - k*(A+B)*(A+B);
}
}
RMax = FindMax(nWidth * nHeight, pdR);
Threshhold(nWidth, nHeight, pdR, pbEdge, RMax, theta);
delete [] piPartialDerivativeX;
delete [] piPartialDerivativeY;
delete [] piXX;
delete [] piYY;
delete [] piXY;
delete [] pdR;
}
void PartialDerivativeX(int nWidth, int nHeight, BYTE *pbGray, int *piPartialDerivativeX)
{
int row, col;
//start from the top row
for (row = nHeight - 1 - 1; row > 0; row--)
{
//start from the left column
for (col = 1; col < nWidth - 1; col++)
{
piPartialDerivativeX[row*nWidth+col] = pbGray[row*nWidth+col+1] - pbGray[row*nWidth+col-1];
}
}
//margin derivative is set to 0
for (row = nHeight - 1; row >= 0; row--)
{
piPartialDerivativeX[row*nWidth] = 0;
piPartialDerivativeX[row*nWidth+nWidth-1] = 0;
}
for (col = 0; col < nWidth; col++)
{
piPartialDerivativeX[col] = 0;
piPartialDerivativeX[nHeight*nWidth-1-col] = 0;
}
}
void PartialDerivativeY(int nWidth, int nHeight, BYTE *pbGray, int *piPartialDerivativeY)
{
int row, col;
//start from the top row
for (row = nHeight - 1 - 1; row > 0; row--)
{
//start from the left column
for (col = 1; col < nWidth - 1; col++)
{
piPartialDerivativeY[row*nWidth+col] = pbGray[(row+1)*nWidth+col] - pbGray[(row-1)*nWidth+col];
}
}
//margin derivative is set to 0
for (row = nHeight - 1; row >= 0; row--)
{
piPartialDerivativeY[row*nWidth] = 0;
piPartialDerivativeY[row*nWidth+nWidth-1] = 0;
}
for (col = 0; col < nWidth; col++)
{
piPartialDerivativeY[col] = 0;
piPartialDerivativeY[nHeight*nWidth-1-col] = 0;
}
}
///////////////////////////////////////////////////////////////////////
// gaussian function (centred at the origin and ignoring the factor
// of 1/(s*sqrt(2*PI)) )
///////////////////////////////////////////////////////////////////////
double gaussian(double x, double s)
{
return(exp(-x/(2*s*s)));
}
void ScalarMultiple(int nWidth, int nHeight, int *piSrc1, int *piSrc2, float *piDst)
{
int i;
int length = nWidth*nHeight;
for (i = 0; i < length; i++)
{
piDst[i] = piSrc1[i] * piSrc2[i];
}
}
//convolution product of piX and w(Gaussian Function)
double ConvolutionProductWithW(int nWidth, int nHeight, int row, int col, float *piX, double *Gauss)
{
int i, j;
double sum = 0;
//the parameter which can be adjusted
//when the value is higher, the result is better,
//but the running speed is much slower
int nWindowScale = 5;
//start from the top row
for (i = min(nHeight - 1 - 1, row + nWindowScale);
i >= max(1, row - nWindowScale); i--)
{
//start from the left column
for (j = max(1, col - nWindowScale);
j <= min(nWidth - 1 - 1, col + nWindowScale); j++)
{
sum += piX[i*nWidth+j] * Gauss[(i-row)*(i-row)+(j-col)*(j-col)];
}
}
return sum;
}
double FindMax(int nLength, double *pdR)
{
int i;
double max = 0;
for (i = 0; i < nLength; i++)
if (pdR[i] > max) max = pdR[i];
return max;
}
void Threshhold(int nWidth, int nHeight, double *pdR, BYTE *pbEdge, double RMax, double theta)
{
int i,j;
int row, col;
int maxIndex;
double maxValue;
double LowerLimit = RMax * theta;
int neighbor = 5; //this parameter can be adjusted
// clear data which is less than the LowerLimit
for (i = 0; i < nWidth * nHeight; i++)
{
if (pdR[i] < LowerLimit) pdR[i] = 0;
}
// make the interest points seperate
//start from the top row
for (row = nHeight - 1 - 1; row > 0; row--)
{
//start from the left column
for (col = 1; col < nWidth - 1; col++)
{
if (pdR[row*nWidth + col] > 0)
{
// check the adjacent elements
maxIndex = min(nHeight - 1 - 1, row + neighbor)*nWidth
+ max(1, col - neighbor);
maxValue = pdR[maxIndex];
for (i = min(nHeight - 1 - 1, row + neighbor);
i >= max(1, row - neighbor); i--)
{
for(j = max(1, col - neighbor);
j <= min(nWidth - 1 - 1, col + neighbor); j++)
{
if (pdR[i*nWidth + j] > maxValue)
{
pdR[maxIndex] = 0;
maxIndex = i*nWidth + j;
maxValue = pdR[maxIndex];
}
else if(i*nWidth + j != maxIndex)
pdR[i*nWidth + j] = 0;
}
}
}
}
}
for (i = 0; i < nWidth * nHeight; i++)
{
if (pdR[i] > 0) pbEdge[i] = 255;
}
}