prof.cpp

这是个人脸识别程序

// Apply the given Mask to the image
// This assumes that ix and iy are in range i.e ix  min value is 1 and
// max ix is Img.width-2 (not Img.width-1)
// Ditto for iy

static inline double FilterMask33 (const Mat &m,
									const int ix, const int iy, const tImMask33 &Mask)
{
const double *const pData = m.m->data;		// for efficiency, acces mat buf directly
const int Tda = m.m->tda;
double Sum;

const double *pData1 = pData + (iy-1) * Tda;
Sum  = Mask[0][0] * pData1[ix-1];
Sum += Mask[0][1] * pData1[ix];
Sum += Mask[0][2] * pData1[ix+1];

pData1 += Tda;
Sum += Mask[1][0] * pData1[ix-1];
Sum += Mask[1][1] * pData1[ix];
Sum += Mask[1][2] * pData1[ix+1];

pData1 += Tda;
Sum += Mask[2][0] * pData1[ix-1];
Sum += Mask[2][1] * pData1[ix];
Sum += Mask[2][2] * pData1[ix+1];

return Sum;
}

//-----------------------------------------------------------------------------
static void InitGrads_Mask33 (Mat &Grads, 									// out
								const Mat &ImgMat, const tImMask33 &Mask)	// in
{
int ncols = ImgMat.ncols();
int nrows = ImgMat.nrows();

for (int iy = 1; iy < nrows-1; iy++)
	for (int ix = 1; ix < ncols-1; ix++)
		Grads(iy, ix) = FilterMask33(ImgMat, ix, iy, Mask);

// fill edges not covered in above for loops

Grads.fill(CONF_UnusedGradVal, 0,       0,       nrows, 1);
Grads.fill(CONF_UnusedGradVal, 0,       0,       1,     ncols);
Grads.fill(CONF_UnusedGradVal, nrows-1, 0,       nrows, ncols);
Grads.fill(CONF_UnusedGradVal, 0,       ncols-1, nrows, ncols);
}

//-----------------------------------------------------------------------------
// Compare this pixel to the pixel on its right and pixel below it.
// Return sum of squares i.e. always non-negative
// This assumes that ix and iy are in range i.e ix  min value is 0 and
// max ix is Img.width-2 (not Img.width-1)
// Ditto for iy

static void InitGrads_GradMagBelowRight (Mat &Grads, 								// out
										 const Mat &ImgMat, unsigned SubProfSpec)	// in
{
int ncols = ImgMat.ncols();
int nrows = ImgMat.nrows();

for (int iy = 0; iy < nrows-1; iy++)
	for (int ix = 0; ix < ncols-1; ix++)
		{
		const double xDelta = ImgMat(iy,  ix+1) - ImgMat(iy,ix);
		const double yDelta = ImgMat(iy+1,ix)   - ImgMat(iy,ix);
		Grads(iy, ix) = sqrt(xDelta*xDelta + yDelta*yDelta);
		}

// fill edges not covered in above for loops

Grads.fill(CONF_UnusedGradVal, nrows-1, 0,       nrows, ncols);
Grads.fill(CONF_UnusedGradVal, 0,       ncols-1, nrows, ncols);
}

//-----------------------------------------------------------------------------
// Compare this pixel to the pixel on its right and pixel below it.
// Return sum, i.e. can be negative
// This assumes that ix and iy are in range i.e ix  min value is 0 and
// max ix is Img.width-2 (not Img.width-1)
// Ditto for iy

static void InitGrads_GradBelowRight (Mat &Grads, 								// out
									  const Mat &ImgMat, unsigned SubProfSpec)	// in
{
int ncols = ImgMat.ncols();
int nrows = ImgMat.nrows();

for (int iy = 0; iy < nrows-1; iy++)
	for (int ix = 0; ix < ncols-1; ix++)
		{
		const double xDelta = ImgMat(iy,  ix+1) - ImgMat(iy,ix);
		const double yDelta = ImgMat(iy+1,ix)   - ImgMat(iy,ix);
		Grads(iy, ix) = xDelta + yDelta;
		}

// fill edges not covered in above for loops

Grads.fill(CONF_UnusedGradVal, nrows-1, 0,       nrows, ncols);
Grads.fill(CONF_UnusedGradVal, 0,       ncols-1, nrows, ncols);
}

//-----------------------------------------------------------------------------
// Calculate edgeness and cornerness as per ScottCootesTaylor "Improving Appearance Model..."
//
// If iEdgeness >= 0   then init Grads[iEdgeness]   to edgeness
// If iCornerness >= 0 then init Grads[iCornerness] to cornerness
#if 0	// obsolete
static void InitGrads_EdgenessCornerness (MatVec &Grads, 								// out
								const Mat &ImgMat, int iEdgeness, int iCornerness,		// in
								bool fHarrisStephens)
{
#define CONF_EWindowWidth 		5			// a positive integer, use 1 for just current pix
#define CONF_EWindowSigma 		CONF_EWindowWidth
#define CONF_EGaussianWindow 	0
#define CONF_EScale 	  		1e18
#define CONF_HarrisStephens_k	0.2
#define CONF_fEdgenessSobel     1
#define CONF_fNew			    0
#define CONF_fESquare			0

#if CONF_fDebug2dProfs
lprintf("\nInitGrads_EdgenessCornerness iEdgeness=%d iCornerness=%d\nWeights ",
	iEdgeness, iCornerness);
#endif

int ncols = ImgMat.ncols();
int nrows = ImgMat.nrows();

// for every pixel in image except those at the edges
for (int iy = CONF_EWindowWidth; iy < nrows-CONF_EWindowWidth; iy++)
	for (int ix = CONF_EWindowWidth; ix < ncols-CONF_EWindowWidth; ix++)
		{
		double dxdx = 0, dydy =0, dxdy = 0;
		// for all pixels in window around ix,iy
		for (int j = -CONF_EWindowWidth+1; j < CONF_EWindowWidth; j++)
			for (int i = -CONF_EWindowWidth+1; i < CONF_EWindowWidth; i++)
				{
				#if CONF_EGaussianWindow
					double Weight = exp(-double((i * i) + (j * j)) /
							(2 * CONF_EWindowSigma * CONF_EWindowSigma));
				#else // equally weighted
					double Weight = 1;
				#endif
				#if CONF_fNew
					//TODO possibly square these
					double dx  = Weight * (ImgMat(ix, iy) - ImgMat(ix+1, iy));
					double dy  = Weight * (ImgMat(ix, iy) - ImgMat(ix,   iy+1));
					double dxy = Weight * (ImgMat(ix, iy) - ImgMat(ix+1, iy+1));
				#elif CONF_fEdgenessSobel
					// vertical Sobel
					double dx = Weight * FilterMask33(ImgMat, ix, iy, gImMasks[PROF_Mask000]);
					// horizontal Sobel
					double dy = Weight * FilterMask33(ImgMat, ix, iy, gImMasks[PROF_Mask202]);
					double dxy = 0;
				#else // same as PROF_GradMagBelowRight
					double dx = ImgMat(iy+j,  ix+i+1) - ImgMat(iy+j,ix+i);
					double dy = ImgMat(iy+j+1,ix+i)   - ImgMat(iy+j,ix+i);
					double dxy = 0;
				#endif
				#if CONF_fESquare
					dxdx += dx * dx;
					dydy += dy * dy;
					dxdy += dx * dy;
				#else
					dxdx += dx;
					dydy += dy;
					dxdy += dxy;
				#endif
				}
		const double A = dxdx;
		const double B = dydy;
		const double C = dxdy;
		const double C2 = C * C;

		if (fHarrisStephens)
			{
			ASSERT(iEdgeness >= 0);
			ASSERT(iCornerness < 0);
			Grads[iEdgeness](iy, ix) = A * B - C2 - CONF_HarrisStephens_k * (A+B) * (A+B);
			}
		else
			{
			if (iEdgeness >= 0)
				Grads[iEdgeness](iy, ix) =
					(A + B) * sqrt((A - B) * (A - B) + 4 * C2) / CONF_EScale;

			if (iCornerness >= 0)
				Grads[iCornerness](iy, ix) = 2 * (A * B - C2) /  CONF_EScale;
			}
		}

// fill edges not covered in above for loops

if (iEdgeness >= 0)
	{
	Grads[iEdgeness].fill(CONF_UnusedGradVal, 0,  0, nrows, CONF_EWindowWidth);
	Grads[iEdgeness].fill(CONF_UnusedGradVal, 0,  0, CONF_EWindowWidth, ncols);
	Grads[iEdgeness].fill(CONF_UnusedGradVal, nrows-CONF_EWindowWidth, 0, nrows, ncols);
	Grads[iEdgeness].fill(CONF_UnusedGradVal, 0, ncols-CONF_EWindowWidth, nrows, ncols);
	}
if (iCornerness >= 0)
	{
	Grads[iCornerness].fill(CONF_UnusedGradVal, 0, 0, nrows, CONF_EWindowWidth);
	Grads[iCornerness].fill(CONF_UnusedGradVal, 0, 0, CONF_EWindowWidth, ncols);
	Grads[iCornerness].fill(CONF_UnusedGradVal, nrows-CONF_EWindowWidth, 0, nrows, ncols);
	Grads[iCornerness].fill(CONF_UnusedGradVal, 0, ncols-CONF_EWindowWidth, nrows, ncols);
	}
#if CONF_fDebug2dProfs //--------------------------------------------
#define DWIDTH 	5
#define xPIX 	341
#define yPIX 	422
if (iEdgeness >= 0)
	{
	lprintf("\nInitGrads_EdgenessCornerness: Edgeness\n");
	for (iy = yPIX-DWIDTH; iy <= yPIX+DWIDTH; iy++)
		{
		for (int ix = xPIX-DWIDTH; ix <= xPIX+DWIDTH; ix++)
			lprintf("%4.0f ", Grads[iEdgeness](iy,ix));
		lprintf("\n");
		}
	}
if (iCornerness >= 0)
	{
	lprintf("InitGrads_EdgenessCornerness: Cornerness\n");
	for (iy = yPIX-DWIDTH; iy <= yPIX+DWIDTH; iy++)
		{
		for (int ix = xPIX-DWIDTH; ix <= xPIX+DWIDTH; ix++)
			lprintf("%4.0f ", Grads[iCornerness](iy,ix));
		lprintf("\n");
		}
	}
#endif  //-----------------------------------------------------------
}
#endif

//-----------------------------------------------------------------------------
// Given an image Img, generate the gradient arrays Grads for the image.
// ProfSpec tells us what gradients to generate.  We generate grads only for 2D profiles.
// We use the term "gradient" loosely -- the gradient could actually be various
// things that aren't quite the gradient.

void InitGrads (MatVec &Grads, const Image &Img, unsigned ProfSpec)
{
int ncols = Img.width;
int nrows = Img.height;
int nSubProfs = nSubProfsForProfSpec(ProfSpec);
Grads.resize(nSubProfs);

Mat ImgMat;					// convert the image to matrix form so no type
							// conversions needed below (for speed)

ImageToMat(ImgMat, Img);	// note: remember that Img(i,j) == ImgMat(j,i)   (swap indices)

for (int iSub = 0; iSub < nSubProfs; iSub++)		// for each subprofile
	{
	unsigned SubProfSpec = GetSubProfSpec(ProfSpec, iSub);
	if (SubProfSpec & PROF_2d)
		{
		Mat *pMat = &Grads[iSub];
		Grads[iSub].dim(nrows, ncols);
		if (SubProfSpec & PROF_FBit)
			{
			unsigned iFilter = (SubProfSpec & PROF_TBits);
			DASSERT(iFilter < sizeof(gImMasks) / sizeof(gImMasks[0]));
			InitGrads_Mask33(Grads[iSub], ImgMat, gImMasks[iFilter]);
			}
		else switch (SubProfSpec & PROF_TBits)		// commonly used profiles treated
			{										// separately for efficiency
			case PROF_GradMagBelowRight:			// gradient magnitude
				InitGrads_GradMagBelowRight(Grads[iSub], ImgMat, SubProfSpec);
				break;
			case PROF_GradBelowRight:				// gradient
				InitGrads_GradBelowRight(Grads[iSub], ImgMat, SubProfSpec);
				break;
			case PROF_Edgeness:
			case PROF_Cornerness:
			case PROF_HarrisStephens:
				break;
			default:
				SysErr("InitGrads %x", SubProfSpec);	// unknown PROF_TBits
				break;
			}
		}
	}
}

//-----------------------------------------------------------------------------
// Return the 2D profile type for this level -- for generating the model (not using the model)
// Returns 0 if there are no 2D profiles at this level.
//
// This calls GetGenProfSpec() to see what type of gradient to generate
// GetGenProfSpec is part of the model configuration (defined in masconf.cpp)
//
// For now, all the 2D profile types must be the same because we store
// Grads on an image basis, not on a landmark basis (but we do create distinct
// Grads for each level in the image)
// So here we check that they do indeed all match and return their common value.

static unsigned GetGenerated2dProfileTypeForThisLev (int iLev, int nPoints)
{
unsigned ProfSpec = 0;
for (int iPoint = 0; iPoint < nPoints; iPoint++)
	{
	unsigned ProfSpec1 = GetGenProfSpec(iLev, iPoint);
	if (ProfSpec1 & PROF_2d)
		{
		if (ProfSpec == 0)
			ProfSpec = ProfSpec1;
		else if (ProfSpec1 != ProfSpec)
			SysErr("ProfSpec not constant over all 2D profiles");
		}
	}
return ProfSpec;
}

//-----------------------------------------------------------------------------
void InitGradsIfNeededForModelGeneration (MatVec &Grads, 	// out:
			const Image &Img, 	// in: Img already scaled to this pyr lev
			int iLev, int nPoints)
{
DASSERT(CONF_nMaxLevs <= 16);	// following init relies on CONF_nMaxLevs less than 16
DASSERT(iLev < CONF_nMaxLevs);

// Cache the 2d prof specs so we don't have to recalculate every time.
// This code relies on the fact that 2d profiles for all landmarks in
// a given level are the same.
// See also CheckThatAll2dProfsMatch().

static unsigned ProfSpecs2d[16] = {-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1};

if (ProfSpecs2d[iLev] == -1)	// first time so need to initialize?
	ProfSpecs2d[iLev] = GetGenerated2dProfileTypeForThisLev(iLev, nPoints);

if (ProfSpecs2d[iLev] & PROF_2d)
	InitGrads(Grads, Img, ProfSpecs2d[iLev]);
}

//-----------------------------------------------------------------------------
// For now, all the 2D profile types must be the same because we store
// Grads on an image basis, not on a landmark basis.
// So here we check that they do indeed all match.
//
// Returns the first 2d ProfSpec in ProfSpecs, if any

static unsigned CheckThatAll2dProfsMatch (const cvec<unsigned> &ProfSpecs, int nPoints)
{
unsigned ProfSpec = 0;
for (int iPoint = 0; iPoint < nPoints; iPoint++)
	{
	unsigned ProfSpec1 = ProfSpecs[iPoint];
	if (ProfSpec1 & PROF_2d)
		{
		if (ProfSpec == 0)
			ProfSpec = ProfSpec1;
		else if (ProfSpec1 != ProfSpec)
			SysErr("ProfSpec not constant over all 2D profiles");
		}
	}
return ProfSpec;
}

//-----------------------------------------------------------------------------
void InitGradsIfNeededForModelSearch (MatVec &Grads,		// out
				const cvec<unsigned> &ProfSpecs,	// in: prof specs read from ASM file
				const Image &Img, 			// in: Img already scaled to this pyramid level
				int nPoints,
				const char sPath[])	// in: if not "" then write gradient image for debug
{
unsigned ProfSpec = CheckThatAll2dProfsMatch(ProfSpecs, nPoints);

if (ProfSpec & PROF_2d)
	{
	InitGrads(Grads, Img, ProfSpec);
	if (sPath[0])	// need to write gradient images for debugging?
		{
		int nSubProfs = nSubProfsForProfSpec(ProfSpec);
		for (int iSub = 0; iSub < nSubProfs; iSub++)
			{
			char s[SLEN]; sprintf(s, "%s-%d.bmp", sPath, iSub);
			WriteMatAsBmp(Grads[iSub], sPath, PROF_NORMALIZE_BMP, QUIET);
			}
		}
	}
}

//-----------------------------------------------------------------------------
int nGetProfWidth (int ncols, unsigned ProfSpec)
{
if ((ProfSpec & PROF_2d) == 0)	// one dimensional profile?
	return ncols;

static int nLastCols;			// a little bit of caching to avoid calling sqrt
static int nLastSqrt;

if (ncols == nLastCols)
	return nLastSqrt;

nLastCols = ncols;
nLastSqrt = sqrt(ncols);

return nLastSqrt;
}