masm.cpp

这是个人脸识别程序

	SubShape(4, VX) = x - Half; SubShape(4, VY) = y - Half;
	RgbImage RgbImg(ScaledImg);
	DrawShape(RgbImg, SubShape, C_YELLOW, 1.0, IM_CONNECT_DOTS);
	WriteBmp(RgbImg, s, VERBOSE);
#endif
	int nSubProfs = nSubProfsForProfSpec(ProfSpec);
	for (int iSub = 0; iSub < nSubProfs; iSub++)
		{
		Mat *pProfs = &Stats.Profs[iLev][iPoint][iSub];		// profiles of current point
		Get2dProf(pProfs->row(iProf), ProfSpec, iSub, ScaledImg, Grads,
			ScaledShape, iPoint, iOffset, iOrthOffset, CONF_nProfWidth2d, true);
		}
	}
else										// one dimensional prof
	{
	DASSERT((ProfSpec & 0x00ffff00) == 0); 					// this code supports only 1 one-dimensional profile per point
	Mat *pProfs = &Stats.Profs[iLev][iPoint][0];			// profiles of current point
	Get1dProf(pProfs->row(iProf), ProfSpec, ScaledImg, iPoint, iOffset, 0, true);
	}
Stats.nProfs[iLev][iPoint] = ++iProf;
}

//-----------------------------------------------------------------------------
static void SynthesizeProf (Mat &NewProf, 		// out
							const Mat &OldProf,	// in
							unsigned ProfSpec)	// in
{
CheckSameDim(NewProf, OldProf, "SynthesizeProf");

if (ProfSpec & PROF_2d)						// two dimensional prof?
	{
	// copy OldProf into NewProf but reverse the order of each row

	MatView OldSquare(OldProf.viewAsSquare());
	MatView NewSquare(NewProf.viewAsSquare());
	int ncols = OldSquare.ncols();
	for (int iRow = 0; iRow < OldSquare.nrows(); iRow++)
		for (int iCol = 0; iCol < ncols; iCol++)
			NewSquare(iRow, iCol) = OldSquare(iRow, ncols - iCol - 1);	// reverse each row
	}
else
	{
	// Copy OldProf into NewProf (no reversal is needed because
	// whisker directions chosen so that we don't need to mirror 1D profiles).
	// Whisker directions are determined by iPrev and iNext in atland.hpp:gLandTabAll
	// TODO For a picture see BothEyesWithNewLandmarksAndWhiskers.jpg

	NewProf = OldProf;
	}
}

//-----------------------------------------------------------------------------
// This is the brother function to CollectProfForOneLandmark but instead of collecting
// the profile for iPoint from the image, we get the profile from iPartner's profile.

static void SynthesizeProfFromPartner (tStats &Stats,									// io, Stats for iPoint are updated
					int iShape, int iLev, int iPoint, int iPartner, const tShapes &Shapes)	// in
{
ASSERT(CONF_nSynthesizedEyePoints == 8);				// for now we can only synthesize exactly 8 new eye points 76_REye0...83_REye7

unsigned Bits = Shapes.pLand[iPoint].Bits;
ASSERT(Bits & FA_Synthesize);
unsigned PartnerBits = Shapes.pLand[iPartner].Bits;
ASSERT(Bits == PartnerBits);

// Using symmetry to generate landmarks only works if whiskers are symmetrical, so fgExplicitPrevNext must be true
// This only applies to 1D profiles because 2D profiles don't use whiskers.
// So check here that this condition if true

if ((fgExplicitPrevNext == false) && !(GetGenProfSpec(iLev, iPoint) & PROF_2d))
	SysErr("SynthesizeProfFromPartner: 1D synthesized profile yet fgExplicitPrevNext==false, iLev %d iPoint %d", iLev, iPoint);

int iProf = Stats.nProfs[iLev][iPoint];							// start where we left off last time

unsigned ProfSpec = GetGenProfSpec(iLev, iPoint);
ASSERT(ProfSpec == GetGenProfSpec(iLev, iPartner));
int nSubProfs = nSubProfsForProfSpec(ProfSpec);

for (int iSub = 0; iSub < nSubProfs; iSub++)
	{
	Mat *pProfs = &Stats.Profs[iLev][iPoint][iSub];				// profiles of current point
	Mat *pPartnerProfs = &Stats.Profs[iLev][iPartner][iSub];	// profiles of corresponding partner
	SynthesizeProf(pProfs->row(iProf), pPartnerProfs->row(iProf), ProfSpec);
	}
Stats.nProfs[iLev][iPoint] = ++iProf;

if (Stats.nProfs[iLev][iPoint] != Stats.nProfs[iLev][iPartner])
	SysErr("A Stats.nProfs[%d][iPoint %d] %d != Stats.nProfs[%d][iPartner %d] %d\n",
		iLev, iPoint, Stats.nProfs[iLev][iPoint], iLev, iPartner, Stats.nProfs[iLev][iPartner]);
}

//-----------------------------------------------------------------------------
static void CollectProfs (tStats &Stats, 		// io
							tShapes &Shapes)	// in: but will be modified if CONF_nGenPrescaleFaceWidth != 0
{
if (CONF_nGenPrescaleFaceWidth != 0)
	{
	ASSERT(!CONF_fAlignImageToAlignedShape);
	int nrows = Shapes.AvShape.nrows();
	if (MRJaw2 >= nrows || MLJaw2 >= nrows)
		SysErr("CONF_nGenPrescaleFaceWidth %d != 0 but can't scale shapes "
			"because Shapes.AvShape.nrows=%d so can't get jaw width",
			CONF_nGenPrescaleFaceWidth, nrows);
	else
		{
		double Scale = (double)CONF_nGenPrescaleFaceWidth / (Shapes.AvShape(MRJaw2, VX) - Shapes.AvShape(MLJaw2, VY));
		Shapes.AvShape *= Scale;
		}
	}
InitPacifyUser(Shapes.nShapes);

for (int iShape = 0; iShape < Shapes.nShapes; iShape++)	// outer loop is image-by-image so we only have to read each image once
	{
#if CONF_fAlignImageToAlignedShape
	SHAPE *pShape = &Shapes.Shapes[iShape];
#else
	SHAPE *pShape = &Shapes.OrgShapes[iShape];
#endif
	PacifyUser(iShape);
	Image Img;
	sLoadImageGivenDirs(Img, Shapes.sImageDirs, &Shapes.TagStrings[iShape][FNAME_OFFSET]);

	if (CONF_fAlignImageToAlignedShape)
		{
		ASSERT(!CONF_nGenPrescaleFaceWidth);
		AlignImageToAlignedShape(Img, Shapes.Shapes[iShape], Shapes.OrgShapes[iShape]);
		}
	else if (CONF_nGenPrescaleFaceWidth != 0)
		{
		// convert all face widths to CONF_nGenPrescaleFaceWidth

		double Scale = (double)CONF_nGenPrescaleFaceWidth / ((*pShape)(MRJaw2, VX) - (*pShape)(MLJaw2, VX));
		ScaleImage(Img, Scale * Img.width, Scale * Img.height, QUIET, CONF_fGenPrescaleBilinear);
		*pShape *= Scale;
		}
	if (CONF_fEqualizeImages)
		{
		int nxMin, nxMax, nyMin, nyMax;
		GetEqualizationMask(nxMin, nxMax, nyMin, nyMax, *pShape, Img, CONF_EqualizationMaskMargin);
		EqualizeAndCorrectImage(Img, nxMin, nxMax, nyMin, nyMax, QUIET);
		}
	// if fgExplicitPrevNext is true:
	//   AlignedAvShape is needed (in PrepareProf1D) only for shapes with missing landmarks.
	//   This is because AlignedAvShape is used only for calculating prev and next landmarks
	//   when a landmark is missing. So we optimize here by only calculating AlignedAvShape when
	//   it is needed i.e. if landmarks are missing for this shape, or if we want
	//   the aligned shape for a debugging image.
	// if fgExplicitPrevNext is false: AlignedAvShape isn't used

	SHAPE AlignedAvShape(Shapes.AvShape);
	if (fgExplicitPrevNext && (Shapes.nUsedLandmarks[iShape] != Shapes.nPoints || CONF_fShowSampledImageWithShapes))
		AlignShape(AlignedAvShape, *pShape);

	pgCurImageName = &Shapes.TagStrings[iShape][FNAME_OFFSET];	// for debugging
	for (int iLev = 0; iLev < CONF_nLevs; iLev++)				// loop level-by-level so we only have to reduce each image once
		{
		double Scale = GetPyrScale(iLev, CONF_PyrRatio);
		SHAPE ScaledShape(*pShape / Scale);
		SHAPE ScaledAlignedAvShape(AlignedAvShape / Scale);
		Image ScaledImg;
		ReduceImageAssign(ScaledImg, Img, Scale, CONF_PyrReduceMethod);
#if CONF_fShowSampledImageWithShapes
		ShowSampledImageWithShapes(iLev, iShape, ScaledImg, ScaledShape, ScaledAlignedAvShape, Shapes);
#endif CONF_fShowSampledImageWithShapes
		MatVec Grads;
		InitGradsIfNeededForModelGeneration(Grads, ScaledImg, iLev, Shapes.nPoints);
		for (int iPoint = 0; iPoint <  Shapes.nPoints; iPoint++)
			{
			bool fCollect = fPointUsed(ScaledShape, iPoint);			// skip landmarks that are unused i.e. 0,0

			#if CONF_fUseOnlyUnobscuredFeats							// skip landmarks that are obscured by glasses etc.
			if (Shapes.TagInts[iShape] & Shapes.pLand[iPoint].Bits & (FA_Glasses|FA_Beard|FA_Mustache))
				fCollect = false;
			#endif

			#if CONF_fUseOnlyUnobscuredFeats && CONF_nSynthesizedEyePoints
			else
			#endif

			#if CONF_nSynthesizedEyePoints
			if (iPoint >= MLEye0 && iPoint <= MLEye7)		// can we collect left eye profile?
				{
				if (!Shapes.fExtraLEyeLandsUsed[iShape])		// not if left eye point is unused
					fCollect = false;
				}
			else if (iPoint >= MREye0 && iPoint <= MREye7)		// can we collect right eye profile?
				{
				if (!Shapes.fExtraREyeLandsUsed[iShape])		// not if right eye point is unused
					fCollect = false;
				}
			#endif

			if (fCollect)
				CollectProfForOneLandmark(Stats, iShape, iLev, iPoint, Shapes.pLand, ScaledImg, ScaledShape, ScaledAlignedAvShape, Grads);
			}
#if CONF_nSynthesizedEyePoints
		for (iPoint = 0; iPoint <  Shapes.nPoints; iPoint++)
			{
			bool fSynthesize = false;		// becomes true if we want to synthesize iPoint from iPartner

			if (iPoint >= MLEye0 && iPoint <= MLEye7 &&										// can we synthesize left eye points from right eye?
				!Shapes.fExtraLEyeLandsUsed[iShape] && Shapes.fExtraREyeLandsUsed[iShape])	// only if left eye point is unused and right eye point is avaliable
				{
				fSynthesize = true;
				}
			if (iPoint >= MREye0 && iPoint <= MREye7 &&										// can we synthesize right eye points from left eye?
				!Shapes.fExtraREyeLandsUsed[iShape] && Shapes.fExtraLEyeLandsUsed[iShape])	// only if right eye point is unused and left eye point is avaliable
				{
				fSynthesize = true;
				}
			#if CONF_fUseOnlyUnobscuredFeats														// skip landmarks that are obscured by glasses etc.
			if (Shapes.TagInts[iShape] & Shapes.pLand[iPoint].Bits & (FA_Glasses|FA_Beard|FA_Mustache))
				fSynthesize = false;
			#endif

			if (fSynthesize)
				{
				int iPartner = Shapes.pLand[iPoint].iPartner;
				ASSERT(iPartner > 0);
				SynthesizeProfFromPartner(Stats, iShape, iLev, iPoint, iPartner, Shapes);			// synthesize iPoint from iPartner
				}
			}
#endif
		}
	}
lprintf("0 %s\n", sSecsElapsed());
}

//-----------------------------------------------------------------------------
// This does the following:
//   Do spectral decomposition: A = Q D Q'  where D is a diagonal matrix of (sorted) eig values
//   Generate Dinv by taking the reciprocal of elements of D that are non-zero, leave zero elements at zero
//   Generate inverted matrix = Q Dinv Q
//
// If A is not singular this returns the same as Mat::inverse() except for small numerical errors
// arising during eigenvector generation.  This routine is also much slower than Mat::inverse().

static void InvertMatRemovingSmallEigs (Mat &A, 			// io: returns inverted A
								 		int &nColinear)		// out: number of colinear colums
{
Vec EigVals(A.nrows());
Vec EigVecs = GetEigsForSymMat(A, EigVals);			// returns sorted eigs
Mat Diag(A.nrows(), A.nrows());

nColinear = 0;
double Min = EigVals(0) / 1000;

for (int i = 0; i < EigVals.nelems(); i++)
	{
	if (ABS(EigVals(i)) < Min)
		{
		Diag(i, i) = 0;
		nColinear++;
		}
	else
		Diag(i, i) = 1 / EigVals(i);
	}
A = EigVecs * Diag * EigVecs.t();
}

//-----------------------------------------------------------------------------
static void InvertMatWithRidgeRegresssion (Mat &A) 	// io: returns inverted A
{
A = CONF_RidgeLambdaForSingCovars * IdentMat(A.nrows()) + (1 - CONF_RidgeLambdaForSingCovars) * A;

A.invertMe();
}

//-----------------------------------------------------------------------------
// This does the following:
//   Do spectral decomposition: A = Q D Q'  where D is a diagonal matrix of (sorted) eig values
//   Set small or negative elements of D to zero
//   Reconstruct A using the new D in A = Q D Q'

static void RemoveSmallEigs (Mat &A, 			// io:
							 int &nDiscarded,	// out:
							 double MaxRatio)	// in:
{
Vec EigVals(A.nrows());
Vec EigVecs = GetEigsForSymMat(A, EigVals);			// returns sorted eigs
Mat Diag(A.nrows(), A.nrows());

double Min = A.trace() / MaxRatio;

nDiscarded = 0;
for (int i = 0; i < EigVals.nelems(); i++)
	if (EigVals(i) < Min)
		nDiscarded++;
	else
		Diag(i, i) = EigVals(i);

A = EigVecs * Diag * EigVecs.t();
}

//-----------------------------------------------------------------------------
// Return true if everything ok with the matrix Covar, with description of Covar in s
// Return false with reason in s if something not good about Covar, writes reason to log file

static bool fCheckCovarMat (const Mat &Covar, int nAllowedZeroEigs)
{
bool f = false;

if (fEqual(Covar(0, 0), 0))						// quick check for first variance
	lprintf("\nCovar first elem is 0 ");
else if (!Covar.fIeeeNormal())
	lprintf("\nCovar not IEEE normal ");
else if (Covar.maxAbsElem() > FLT_MAX)			// use FLT_MAX as a convenient big number (for checking valid inversion of matrix)
	lprintf("\nCovar elem too big ");
else if (!fPosDef1(Covar, nAllowedZeroEigs))	// will test for symmetry too
	lprintf("\nCovar not pos def ");
else
	f = true;

return f;
}

//-----------------------------------------------------------------------------
static void GenAvProf (Mat &AvProf, 					// out
					   const Mat &Profs, int nProfs) 	// in
{
for (int iProf = 0; iProf < nProfs; iProf++)
	AvProf += Profs.row(iProf);				// add current row of Profs

AvProf /= nProfs;
}

//-----------------------------------------------------------------------------
// This generates the inverted covariance matrix Covar from Profs and AvProfs.
//
// This routine does various processing (depending on #defines) to ensure
// that the covariance matrix is actually invertible and positive definite.
// This is needed, for example, if we previously normalized all profiles to the same length
// -- this makes the profiles colinear resulting in a singular covariance matrix.
//
// If the inverted matrix is still not ok after the above processing, we convert
// it to an  identity matrix (so all covariances are 0 and all variances are 1)