prof.cpp

这是个人脸识别程序

					}
				default:
					SysErr("WindowType %x", WindowType);
					break;
				}
	}
return &Window;
}

//-----------------------------------------------------------------------------
// Get the profiles from Grads
// The Grads arrays must be initialized before calling this function

static void Get2dProfFromGradMats (Vec &Prof, 												   // out
				 	unsigned ProfSpec, int iSubProf, const MatVec &Grads, const SHAPE &Shape,  // in
					int iPoint, int ixOffset, int iyOffset, int nProfWidth, bool fDumpRawProf) // in
{
DASSERT(fOdd(nProfWidth));
DASSERT(Grads[iSubProf].nrows() && Grads[iSubProf].ncols());
	// ensure gradient matrices are inited, done previously by InitGradsIfNeeded...

unsigned SubProfSpec = GetSubProfSpec(ProfSpec, iSubProf);

int nProfWidthEachSide = (nProfWidth - 1) / 2;
int nrows = Grads[iSubProf].nrows();
int ncols = Grads[iSubProf].ncols();

// convert shape coords to mat coords, offset by ixOffset and iyOffset and nProfWidthEachSide

int iy = iRound(Shape(iPoint, VY)) - iyOffset;
int ix = iRound(Shape(iPoint, VX)) - ixOffset;

int iyMin = nrows/2 - 1 - iy - nProfWidthEachSide;
int iyMax = nrows/2 - 1 - iy + nProfWidthEachSide;
int ixMin = ncols/2     + ix - nProfWidthEachSide;
int ixMax = ncols/2     + ix + nProfWidthEachSide;

#if CONF_fRawProfs
if (fDumpRawProf)
	{
	if (!pgRawProfFile)
		pgRawProfFile = Fopen(CONF_sRawProfFile, "w");
	Fprintf(pgRawProfFile, "# iPoint %d ixOffset %d iyOffset %d\n", iPoint,ixOffset,iyOffset);
	Fprintf(pgRawProfFile, "2D prof, not saved\n");	// TODO for now
	}
#endif

Mat *pWindow = pGetProfWindow(SubProfSpec, nProfWidth);
double *pData = Prof.m->data;						// for efficiency, acces mat buf directly
													// TODO make this a func in mat.hpp
ASSERT(Prof.nelems() == nProfWidth * nProfWidth);	// make sure we don't overrun Prof
const double *pGrad = Grads[iSubProf].m->data;		// ditto
int Tda = Grads[iSubProf].m->tda;
int iProf = 0;

// For efficiency there are two loops: one to hande the case where the profile is entirely in the
// image boundaries and one for the case where the profile goes off the edge the edge of the image

if (ixMin >= 0 && ixMax < ncols && iyMin >= 0 && iyMax < nrows)	// profile is in image boundary?
	{
	// Profile is in image boundary
	// A further efficiency measure is to once again have two loops,
	// since we can ignore equally weighted windows

	if ((SubProfSpec & PROF_WindowField) == PROF_WindowEquallyWeighted)
		{
		for (int iy = iyMin; iy <= iyMax; iy++)					// no weights
			{
			const double *pGrad1 = pGrad + iy * Tda + ixMin;
			for (int ix = ixMin; ix <= ixMax; ix++)
				*pData++ = *pGrad1++;
			}
		}
	else														// weights
		for (int iy = iyMin; iy <= iyMax; iy++)
			{
			const double *pGrad1 = pGrad + iy * Tda + ixMin;
			for (int ix = ixMin; ix <= ixMax; ix++)
				*pData++ = (*pWindow)(iy-iyMin, ix-ixMin) * *pGrad1++;
			}
	}
else															// profile crosses the image boundary
	{
	for (int iy = iyMin; iy <= iyMax; iy++)
		{
		const double *pGrad1 = pGrad + iy * Tda + ixMin;
		for (int ix = ixMin; ix <= ixMax; ix++)
			{
			if (ix >= 0 && ix < ncols && iy >= 0 && iy < nrows)
				*pData++ = (*pWindow)(iy-iyMin, ix-ixMin) * *pGrad1;
			else
				*pData++ = CONF_UnusedGradVal;
			pGrad1++;
			}
		}
	}
#if 0
if (iPoint == 15)
	Prof.print("Prof ", "%g ");
#endif
}

//-----------------------------------------------------------------------------
// Depending on SubProfSpec, this gets edgeness, cornerness, or harris-stephens value
// Edgeness or cornerness calculated as per ScottCootesTaylor "Improving Appearance Model..."

static double GetEdgeness (int ix, int iy, const Image &Img, const unsigned SubProfSpec)
{
#define CONF_nEWidth 	0
#define CONF_Harris_k	0.2

#define CONF_EGrad		1
#define CONF_EVxlA		0
#define CONF_EVxlB		0

double A = 0, B = 0, C = 0;

for (int j = -CONF_nEWidth; j <= CONF_nEWidth; j++)	// for all pixels in nEWidth window around ix,iy
	for (int i = -CONF_nEWidth; i <= CONF_nEWidth; i++)
		{
		int iy1 = iy+j;
		int ix1 = ix+i;
		if (ix1 > 1 && ix1 < Img.width-2 && iy1 > 1 && iy1 < Img.height-2)
			{
			#if CONF_EGrad
				double dx = Img(ix1, iy1) - Img(ix1+1, iy1);
				double dy = Img(ix1, iy1) - Img(ix1,   iy1+1);
				A += dx * dx;
				B += dy * dy;
				C += dx * dy;
			#elif CONF_EVxlA
				// o1 o2 o3
				// o4    o5
				// o6 o7 o8
				double o1 = Img(ix1-1, iy1-1);
				double o2 = Img(ix1,   iy1-1);
				double o3 = Img(ix1+1, iy1-1);
				double o4 = Img(ix1-1, iy1);

				double o5 = Img(ix1+1, iy1);
				double o6 = Img(ix1-1, iy1+1);
				double o7 = Img(ix1,   iy1+1);
				double o8 = Img(ix1+1, iy1+1);
 				A += (o3+o8 - 4*(o1+o6)) + 2*(o5-o4);
				B += (o1+o3 - 4*(o6+o8)) + 2*(o2-o7);
				C += (o1+o3 - 4*(o6+o8)) + 2*(o2-o7);
			#elif CONF_EVxlB	// based on vil_corners.vxx
				// 1 2 3
				// 4 0 5
				// 6 7 8

				// x gradients
				double xg1 = Img(ix1-1, iy1-1) - Img(ix1,   iy1-1);
				double xg2 = Img(ix1,   iy1-1) - Img(ix1+1, iy1-1);
				double xg3 = Img(ix1+1, iy1-1) - Img(ix1+2, iy1-1);
				double xg4 = Img(ix1-1, iy1)   - Img(ix1,   iy1);
				double xg0 = Img(ix1,   iy1)   - Img(ix1+1, iy1);
				double xg5 = Img(ix1+1, iy1)   - Img(ix1+2, iy1);
				double xg6 = Img(ix1-1, iy1+1) - Img(ix1,   iy1+1);
				double xg7 = Img(ix1,   iy1+1) - Img(ix1+1, iy1+1);
				double xg8 = Img(ix1+1, iy1+1) - Img(ix1+2, iy1+1);

				// y gradients
				double yg1 = Img(ix1-1, iy1-1) - Img(ix1-1, iy1);
				double yg2 = Img(ix1,   iy1-1) - Img(ix1,   iy1);
				double yg3 = Img(ix1+1, iy1-1) - Img(ix1+1, iy1);
				double yg4 = Img(ix1-1, iy1)   - Img(ix1-1, iy1+1);
				double yg0 = Img(ix1,   iy1)   - Img(ix1,   iy1+1);
				double yg5 = Img(ix1+1, iy1)   - Img(ix1+1, iy1+1);
				double yg6 = Img(ix1-1, iy1+1) - Img(ix1-1, iy1+2);
				double yg7 = Img(ix1,   iy1+1) - Img(ix1,   iy1+2);
				double yg8 = Img(ix1+1, iy1+1) - Img(ix1+1, iy1+2);

 				// A += (xg3+xg8 - (xg1+xg6)) + 2*(xg5-xg4);
				// B += (yg1+yg3 - (yg6+yg8)) + 2*(yg2-yg7);
				// C += (xg1+xg3 - (xg6+xg8)) + 2*(xg2-xg7);
 				A += xg5-xg4;
				B += yg2-yg7;
				C += xg1+yg1;
			#else	// sobel
				// o1 o2 o3
				// o4    o5
				// o6 o7 o8
				double o1 = Img(ix1-1, iy1-1);
				double o2 = Img(ix1,   iy1-1);
				double o3 = Img(ix1+1, iy1-1);
				double o4 = Img(ix1-1, iy1);
				double o5 = Img(ix1+1, iy1);
				double o6 = Img(ix1-1, iy1+1);
				double o7 = Img(ix1,   iy1+1);
				double o8 = Img(ix1+1, iy1+1);
				double dx = -o1 - 2*o4 - o6 + o3 + 2*o5 + o8;
				double dy = -o1 - 2*o2 - o3 + o6 + 2*o7 + o8;
				A += dx * dx;
				B += dy * dy;
				C += dx * dy;
			#endif
			}
		}

double Ret = 0;
switch (SubProfSpec & PROF_TBits)
	{
	case PROF_Edgeness:
		Ret = (A + B) * sqrt((A - B) * (A - B) + 4 * C * C);
		break;
	case PROF_Cornerness:
		Ret = 2 * (A * B - C * C);
		break;
	case PROF_HarrisStephens:
		Ret = A * B - C * C - CONF_Harris_k * (A+B) * (A+B);
		break;
	default:
		SysErr("GetEdgeness");
	}
return Ret;
}

//-----------------------------------------------------------------------------
static void Get2dProfFromImage (Vec &Prof, 											// out
				 	unsigned SubProfSpec, const Image &Img, const SHAPE &Shape, 	// in
					int iPoint, int ixOffset, int iyOffset, int nProfWidth)			// in
{
ASSERT(fOdd(nProfWidth));

ASSERT((SubProfSpec & PROF_FBit) == 0);

ASSERT((SubProfSpec & PROF_TBits) == PROF_Edgeness   ||
	     (SubProfSpec & PROF_TBits) == PROF_Cornerness ||
	     (SubProfSpec & PROF_TBits) == PROF_HarrisStephens);

int nProfWidthEachSide = (nProfWidth - 1) / 2;

int iy = iRound(Shape(iPoint, VY)) - iyOffset;
int ix = iRound(Shape(iPoint, VX)) - ixOffset;

int iyMin = iy - nProfWidthEachSide;
int iyMax = iy + nProfWidthEachSide;
int ixMin = ix - nProfWidthEachSide;
int ixMax = ix + nProfWidthEachSide;

#if 0
if (iPoint == 15)
	{
	RgbImage RgbImg(Img);
	DrawShape(RgbImg, Shape, C_RED, 1.0);

	DrawRectInBmp(RgbImg,
		xShapeToImg(ixMin, Img), Img.height/2-iyMax, nProfWidth, nProfWidth, C_YELLOW);

	SHAPE Shape(3,2);
	Shape(0,VX) = ixMin; Shape(0,VY) = iyMin;
	Shape(1,VX) = ixMax; Shape(1,VY) = iyMax;
	Shape(2,VX) = ix;    Shape(2,VY) = iy;
	DrawShape(RgbImg, Shape, C_BLUE, 1.0);

	static int i;
	char sPath[SLEN]; sprintf(sPath, "out/%2.2d", i);
	WriteBmp(RgbImg, sPath, VERBOSE);
	i++;
	}
#endif
double *pData = Prof.m->data;						// for efficiency, acces mat buf directly
													// TODO make this a func in mat.hpp
ASSERT(Prof.nelems() == nProfWidth * nProfWidth);	// make sure we don't overrun Prof

iyMin = Img.height/2  - iyMin - 1;					// shape to image coords
iyMax = Img.height/2  - iyMax - 1;
ixMin = Img.width/2  + ixMin;
ixMax = Img.width/2  + ixMax;

for (iy = iyMax; iy <= iyMin; iy++)
	for (ix = ixMin; ix <= ixMax; ix++)
		*pData++ = GetEdgeness(ix, iy, Img, SubProfSpec);

#if 0
if (iPoint == 15)
	Prof.print("Prof ", "%g ");
#endif
}

//-----------------------------------------------------------------------------
// Return true if got the prof, false if not (because iPoint is off the image)

void Get2dProf (Vec &Prof, 														// out
		 	const unsigned ProfSpec, const int iSubProf,						// in
			const Image &Img, const MatVec &Grads, 								// in
			const SHAPE &Shape, const int iPoint, const int ixOffset, 			// in
			const int iyOffset, const int nProfWidth, const bool fDumpRawProf)	// in
{
unsigned SubProfSpec = GetSubProfSpec(ProfSpec, iSubProf);
bool fFromGradMats = true;

if ((SubProfSpec & PROF_FBit) == 0)
	{
	switch (SubProfSpec & PROF_TBits)
		{
		case PROF_Edgeness:
		case PROF_Cornerness:
		case PROF_HarrisStephens:
			fFromGradMats = false;
		}
	}
if (fFromGradMats)
	Get2dProfFromGradMats(Prof, ProfSpec, iSubProf, Grads, Shape, iPoint,
		ixOffset, iyOffset, nProfWidth, fDumpRawProf);
else
	Get2dProfFromImage(Prof, SubProfSpec, Img, Shape, iPoint, ixOffset, iyOffset, nProfWidth);

Normalize2d(Prof, SubProfSpec);

#if CONF_fDebug2dProfs
if (iPoint == 30 && ixOffset == 0 && iyOffset == 0)
	{
	lprintf("iPoint %d ixOffset %d iyOffset %d\n", iPoint, ixOffset, iyOffset);
	// show the profile matrix
	static int iBmp = 0; char s[SLEN]; sprintf(s, "out/Prof%2.2d.bmp", iBmp++);
	static int iGrad = 0;
	Mat SquareProf(Prof);
	SquareProf.reshapeMe(nProfWidth, nProfWidth, false, 0, 0, 1, Prof.nelems());
	char s1[SLEN];
	sprintf(s1, "\nProf iPoint %d ixOffset %d iyOffset %d\n", iPoint, ixOffset, iyOffset);
	SquareProf.print(s1, "%14.8f", NULL, NULL, 0, true);
	WriteMatAsBmp(SquareProf, s, true, true);
	// Err("early");
	}
#endif
}

//-----------------------------------------------------------------------------
void WriteProfSpecs (int iLev, int nPoints, FILE *pFile)
{
// for ProfSpecs we use the same layout as a (double) Mat (but this is unsigned hex data)

Fprintf(pFile, "{ %d 1\n", nPoints);
for (int iPoint = 0; iPoint < nPoints; iPoint++)
	Fprintf(pFile, "%x\n", GetGenProfSpec(iLev, iPoint));
Fprintf(pFile, "}\n");
}

//-----------------------------------------------------------------------------
// This duplicates the functionality of Mat.sread but for cvec<unsigned>
// One day replace this function and ReadIntVec TODO

void ReadProfSpecs (cvec<unsigned> &ProfSpecs, 								// out
					  int nPoints, const char sAsmFile[], FILE *pAsmFile)	// in
{
char s[SLEN];

// header "{ 68 1"

if (!Fgets(s, SLEN-1, pAsmFile))
	SysErr("Can't read ProfSpecs header from %s", sAsmFile);
char sBrace[SLEN];
int nrows = -1, ncols = -1;
if (sscanf(s, "%s %d %d", sBrace, &nrows, &ncols) != 3 || sBrace[0] != '{' || sBrace[1] != 0)
	SysErr("Bad format ProfSpecs header: %s", s);

if (nrows != nPoints || ncols != 1)
	SysErr("Bad number of rows or cols in ProfSpecs header, got %d %d expected %d 1",
		nrows, ncols, nPoints, s);

// body, lines like "1000000"

ProfSpecs.resize(nPoints);
for (int iPoint = 0; iPoint < nPoints; iPoint++)
	{
	if (!Fgets(s, SLEN-1, pAsmFile))
		SysErr("Can't read ProfSpecs from %s", sAsmFile);
	unsigned ProfSpec = -1;
	if (sscanf(s, "%x", &ProfSpec) != 1)
		SysErr("Can't read ProfSpecs from %s", sAsmFile);
	ProfSpecs[iPoint] = ProfSpec;
	}
// footer "}"

if (!Fgets(s, SLEN-1, pAsmFile))
	SysErr("Can't read ProfSpecs footer %s", sAsmFile);
if (s[0] != '}')
	SysErr("Bad format ProfSpecs footer, expected closing brace, got %s", s);
}

//-----------------------------------------------------------------------------
// After calling this function, m(i,j) corresponds to Img(j,i).  Note the reversal of i,j.
//
// This means m.ncols() == Img.width
// 			  m.nrows() == Img.height
//
// And m.print() prints a matrix that must be rotated anticlockwise
// by 90 degrees to "look like" the corresponding image.

static void ImageToMat (Mat &m, const Image &Img)
{
m.dim(Img.height, Img.width);

const int Tda = m.m->tda;

for (int iy = 0; iy < Img.height; iy++)
	{
	double *pData = m.m->data + iy * Tda;		// for efficiency, acces mat buf directly
	for (int ix = 0; ix < Img.width; ix++)
		*pData++ = Img(ix, iy);
	}
}

//-----------------------------------------------------------------------------