asmsearch.cpp

这是个人脸识别程序

// $masm\asmsearch.cpp 1.5 milbo$ active shape routines for searching
// Warning: this is raw research code -- expect it to be quite messy.
// milbo dec05 durban
//-----------------------------------------------------------------------------
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// A copy of the GNU General Public License is available at
// http://www.r-project.org/Licenses/
//-----------------------------------------------------------------------------

#include "all.hpp"

//-----------------------------------------------------------------------------
// This returns the distance between ProfDiff = Prof - AvProf,
// using the inverted covar mat Covar.
// The efficiency of this routine is crucial to ASM search speed.
//
// ProfDiff is a row vector.  In math textbooks, vectors are column vectors
// So we take tranposes differently from the formulae you see in textbooks.
// Also, Covar is actually the inverted covariance matrix.

static double GetProfDist (const Vec &ProfDiff, 		// in: Prof-AvProf
						   const Mat &Covar, const SparseMat *pSparseCovar, const unsigned SubProfSpec, const int nTrimCovar)
{
DASSERT(CONF_fStandardProfDist + CONF_fAbsProfDiff + CONF_fExperimentalProfDist == 1);	// only one of these should be set

double Dist;

#if CONF_fStandardProfDist //------------------------------------------------------

// Following code implements Dist = (Prof - AvProf) * Covar * (Prof - AvProf).t()
// but is faster than just using the standard matrix routines.

if (pSparseCovar != NULL)
	Dist = Sparse_xAx(ProfDiff, *pSparseCovar);
else if (nTrimCovar)													// only needed for old format .asm files
	Dist = Trimmed_xAx(ProfDiff, Covar, SubProfSpec, nTrimCovar);
else
	Dist = xAx(ProfDiff, Covar);

static bool fWarned = false;
if (!fWarned && Dist < 0)
	{
	Warn("NegDist %g %x", Dist, SubProfSpec);	// this happens if the covar mat is not pos def
	fWarned = true;
	}

#elif CONF_fAbsProfDiff  //----------------------------------------------------

// weighted sum of differences (see my Feb06 notebook p37) (TODO doesn't give better results)
// TODO also tried various sigmoid scalings of this result and also sigmoids on QuadProfDist result without improvement

Dist = 0;
for (int i = 0; i < Prof.nelems(); i++)
	for (int j = 0; j < Prof.nelems(); j++)
		{
		double x = Prof(i) - AvProf(i);
		double y = Prof(j) - AvProf(j);
		double w = Covar(i, j); 							// weight
#if 1 // methodA (see notebook)
		double Result = fabs(fabs(w) * x + w * y);
#elif 0 // methodB
		double Result;
		bool fAdd = (fPos(x) == fPos(y));					// add if both elements have same sign
		if (fAdd)
			Result = w * (x + y);
		else
			Result = w * (x - y);
		bool fPosResult = (fAdd? fPos(w): !fPos(w));		// true if sign of result must be positive
		if ((fPosResult && Result < 0) || (!fPosResult && Result > 0))
			Result = -Result;
		DASSERT(fPos(Result) == fPos(w * x * y));			// check that sign of result is same as using xAx
#elif 0 // same as CONF_fStandardProfDist xAx
		double Result = w * x * y;
#endif
		// Result /= (1 + fabs(Result));					// sigmoid scaling
		Dist += Result;
		}

#elif CONF_fExperimentalProfDist  //-------------------------------------------

DASSERT(Covar.nelems());

//TODO this will be faster if we precalculate various elements of this
// (Covar * AvProf.t()).print("Covar * AvProf.t()\n", "%g ");
// (.5 * AvProf  * Covar * AvProf.t()).print("(.5 * AvProf  * Covar * AvProf.t())\n", "%g ");
//Dist = OneElemToDouble(.5 * AvProf  * Covar * AvProf.t() - Prof * Covar * AvProf.t());

// Results for following piece of code m3.asm CONF_nProfWidth1d (5 * 2 + 1)=11
//
// DENOM = 100
// ALL   N=50   A         SDev      B         SDev      T       TProb  ChangeBtoA  BetterA WorseB   Base  ChangeBaseToA
// Fit         0.053     0.019     0.054     0.022     -0.97    0.248      3.4%      44%    44%     0.094   77.8%
// LEyeFit     0.043     0.026     0.046     0.028     -1.33    0.162      5.4%      42%    54%     0.058   35.4%
// REyeFit     0.050     0.034     0.054     0.037     -1.67    0.099      8.5%      42%    58%     0.072   45.0%
// MouthFit    0.045     0.025     0.046     0.025     -0.11    0.394      0.6%      38%    50%     0.105  132.0%
// PointFit35  0.441     0.518     0.441     0.516      0.00    0.397     -0.0%      52%    48%     0.464    5.3%
// Max         0.136               0.146                                   7.8%
//
// DENOM = 50
// ALL   N=50   A         SDev      B         SDev      T       TProb  ChangeBtoA  BetterA WorseB   Base  ChangeBaseToA
// Fit         0.056     0.023     0.054     0.022      1.03    0.232     -2.8%      56%    36%     0.094   67.2%
// LEyeFit     0.045     0.023     0.046     0.028     -0.22    0.387      1.3%      48%    50%     0.058   30.1%
// REyeFit     0.056     0.042     0.054     0.037      0.45    0.358     -2.7%      54%    38%     0.072   30.0%
// MouthFit    0.045     0.024     0.046     0.025     -0.19    0.390      0.7%      52%    42%     0.105  132.2%
// PointFit35  0.435     0.513     0.441     0.516     -0.10    0.395      1.3%      46%    54%     0.464    6.6%
// Max         0.149               0.146                                  -1.7%
//
// DENOM = 1000
// ALL   N=50   A         SDev      B         SDev      T       TProb  ChangeBtoA  BetterA WorseB   Base  ChangeBaseToA
// Fit         0.054     0.020     0.054     0.022     -0.54    0.342      1.6%      50%    44%     0.094   74.7%
// LEyeFit     0.046     0.026     0.046     0.028      0.06    0.396     -0.3%      44%    52%     0.058   28.1%
// REyeFit     0.053     0.037     0.054     0.037     -0.53    0.344      2.1%      50%    44%     0.072   36.4%
// MouthFit    0.042     0.022     0.046     0.025     -1.72    0.092      8.3%      36%    58%     0.105  149.7%
// PointFit35  0.440     0.517     0.441     0.516     -0.01    0.397      0.1%      50%    50%     0.464    5.4%
// Max         0.151               0.146                                  -2.9%
//
// DENOM 500
// ALL   N=50   A         SDev      B         SDev      T       TProb  ChangeBtoA  BetterA WorseB   Base  ChangeBaseToA
// Fit         0.053     0.018     0.054     0.022     -0.89    0.265      2.7%      44%    50%     0.094   76.6%
// LEyeFit     0.044     0.024     0.046     0.028     -0.55    0.340      2.8%      46%    48%     0.058   32.1%
// REyeFit     0.050     0.035     0.054     0.037     -1.60    0.111      7.3%      38%    60%     0.072   43.4%
// MouthFit    0.043     0.023     0.046     0.025     -1.15    0.204      5.5%      36%    58%     0.105  143.3%
// PointFit35  0.436     0.516     0.441     0.516     -0.09    0.395      1.1%      34%    66%     0.464    6.4%
// Max         0.146               0.146                                   0.1%
//
// DENOM 200
// ALL   N=50   A         SDev      B         SDev      T       TProb  ChangeBtoA  BetterA WorseB   Base  ChangeBaseToA
// Fit         0.053     0.019     0.054     0.022     -0.73    0.303      2.3%      50%    44%     0.094   75.9%
// LEyeFit     0.045     0.028     0.046     0.028     -0.47    0.356      1.8%      46%    52%     0.058   30.8%
// REyeFit     0.050     0.035     0.054     0.037     -1.95    0.061      8.5%      40%    56%     0.072   44.9%
// MouthFit    0.044     0.022     0.046     0.025     -0.91    0.261      3.9%      54%    42%     0.105  139.6%
// PointFit35  0.436     0.517     0.441     0.516     -0.09    0.395      1.1%      44%    56%     0.464    6.4%
// Max         0.148               0.146                                  -1.1%

#define DENOM 1000
double Trace = Covar.trace();
Mat Covar1(Covar);
for (int i = 0; i < Covar.nrows(); i++)
	for (int j = 0; j < Covar.ncols(); j++)
		if (i != j && ABS(Covar(i, j)) < Trace / DENOM)
			Covar1(i,j) = 0;

Vec Diff(Prof - AvProf);
Dist = OneElemToDouble(Diff * Covar1 * Diff.t());

static bool fWarned = false;
if (!fWarned  && Dist < 0)
	{
	Warn("NegDist %g", Dist);	// this happens if the covar mat is not pos def
	fWarned = true;
	}

#else   //---------------------------------------------------------------------
SysErr("GetProfDist")			// no distance method specified
#endif   //--------------------------------------------------------------------

return Dist;
}

//-----------------------------------------------------------------------------
// ix and iy are the offset and orthogonal offset wrt iPoint of the desired profile
// Returns the fit
//
// This routine need to be very efficient because it is called a LOT during search.

double GetProfileFit (tSearchImages &SearchImgs, 										// in except SearchImgs.ProgressImg is modified
					  const int iPoint, const int ix, const int iy,
					  const tAsmLev &AsmLev, const SHAPE &Shape, const tLand Lands[],	// in
					  const unsigned ProfSpec, const int nTrimCovar)
{
static Vec Prof;	// it is slightly quicker to allocate this statically and call dimClear() each time
					// because most of the time the dimension is the same as the previous profile
double Fit = 0;

int nSubProfs = nSubProfsForProfSpec(ProfSpec);
for (int iSub = 0; iSub < nSubProfs; iSub++)
	{
	Prof.dimClear(1, AsmLev.Profs[iPoint][iSub].ncols());
	int nProfileWidth = nGetProfWidthFromModel(iPoint, iSub, AsmLev);
	bool fSuccess = true;
	if (AsmLev.ProfSpecs[iPoint] & PROF_2d)		// two dimensional profile?
		Get2dProf(Prof, AsmLev.ProfSpecs[iPoint], iSub, SearchImgs.Img, SearchImgs.Grads, Shape, iPoint, ix, iy, nProfileWidth, false);
	else
		Get1dProf(Prof, AsmLev.ProfSpecs[iPoint], SearchImgs.Img, iPoint, ix, 0, false);

	double Weight = (iSub > 0? CONF_CombineProfWeight: 1.0);
	unsigned SubProfSpec = GetSubProfSpec(ProfSpec, iSub);
	int nTrimCovar1 = ((SubProfSpec & PROF_2d)? nTrimCovar: 0);
	const SparseMat *pSparseCovar = //TODO tidy this up, see ReadAsmLev()
			((AsmLev.Covars[iPoint][iSub].nrows() > 1)? NULL: &AsmLev.SparseCovars[iPoint][iSub]);
	Prof -= AsmLev.Profs[iPoint][iSub];	// for efficiency, use operator-= rather than operator-
	Fit += Weight * GetProfDist(Prof, AsmLev.Covars[iPoint][iSub], pSparseCovar, SubProfSpec, nTrimCovar1);
	AccumulateFeats(ix, iy, iSub, Fit, Prof, ProfSpec);		// only does anything if -G flag was specified
	}
return Fit;
}

//-----------------------------------------------------------------------------
// Read header of .asm file

static void ReadAsmHeader (int &nPoints, int &nLevs,							// out
							double &PyrRatio,									// out
							int &PyrReduceMethod,								// out
							bool &fExplicitPrevNext, double &NormalizedProfLen, // out
							int	 &nScaledFaceWidth, bool &fBilinearRescale,		// out
							int	 &nTrimCovar, double &SigmoidScale,				// out
							const char sAsmFile[], FILE *pAsmFile)				// in
{
// expect something like this: ASM0 [../asm/ms.asm]

char s[SLEN], sFile[SLEN];
int nVersion = -1;
Fgets(s, SLEN-1, pAsmFile);

nPoints = -1;

if (strncmp(s, "ASM", 3) != 0)			// check ASM file magic number
	Err("%s is not an ASM file (first three chars are not \"ASM\")", sAsmFile);

if (sscanf(s+3, "%d %s", &nVersion, sFile) != 2)
	SysErr("Bad format ASM file header: %s", s);

if (nVersion != CONF_nAsmFileVersion)
	SysErr("Bad version number %d in ASM file header", nVersion);

//TODO should create functions for the repeated code below

char sToken[SLEN];
Fgets(s, SLEN-1, pAsmFile);
int Temp = -1;
if (2 != sscanf(s, "%s %d", sToken, &Temp) || 0 != strcmp(sToken, "Points"))
	SysErr("Can't read Points from ASM file");
if (Temp < 1 || Temp > 1000)
	SysErr("Bad number of Points %d in ASM file header", Temp);
nPoints = Temp;

Fgets(s, SLEN-1, pAsmFile);
Temp = -1;
if (2 != sscanf(s, "%s %d", sToken, &Temp) || 0 != strcmp(sToken, "Levs") || Temp < 0)
	SysErr("Can't read Levs from ASM file");
if (Temp < 1 || Temp > CONF_nMaxLevs)
	SysErr("Bad number of levels %d in ASM file header", Temp);
nLevs = Temp;

Fgets(s, SLEN-1, pAsmFile);
float DTemp = -1.0;
if (2 != sscanf(s, "%s %g", sToken, &DTemp) || 0 != strcmp(sToken, "PyramidRatio") || DTemp <= 0 || DTemp > 1000)	// 1000 is arbitrary
	SysErr("Can't read PyramidRatio from ASM file");
PyrRatio = DTemp;

Fgets(s, SLEN-1, pAsmFile);
Temp = -1;
if (2 != sscanf(s, "%s %d", sToken, &Temp) || 0 != strcmp(sToken, "PyramidReduceMethod"))
	SysErr("Can't read PyramidReduceMethod from ASM file");
if (Temp < IM_NEAREST_PIXEL || Temp > IM_AVERAGE_ALL)
	SysErr("Bad PyramidReduceMethod %d in ASM file header", Temp);
PyrReduceMethod = Temp;

Fgets(s, SLEN-1, pAsmFile);
Temp = -1;
if (2 != sscanf(s, "%s %d", sToken, &Temp) || 0 != strcmp(sToken, "ExplicitPrevNext") 
		|| (Temp != 0 && Temp != 1))
	SysErr("Can't read ExplicitPrevNext from ASM file");
fExplicitPrevNext = (Temp == 1);

Fgets(s, SLEN-1, pAsmFile);
DTemp = -1.0;
if (2 != sscanf(s, "%s %g", sToken, &DTemp) || 0 != strcmp(sToken, "NormalizedProfLen") 
		|| DTemp <= 0 || DTemp > 1000)	// 1000 is arbitrary
	SysErr("Can't read NormalizedProfLen from ASM file");
NormalizedProfLen = DTemp;

Fgets(s, SLEN-1, pAsmFile);
Temp = -1;
if (2 != sscanf(s, "%s %d", sToken, &Temp) || 0 != strcmp(sToken, "ScaledFaceWidth") 
		|| Temp < 0 || Temp > 1000)	// 1000 is arbitrary
	SysErr("Can't read ScaledFaceWidth from ASM file");
nScaledFaceWidth = Temp;

Fgets(s, SLEN-1, pAsmFile);
Temp = -1;
if (2 != sscanf(s, "%s %d", sToken, &Temp) || 0 != strcmp(sToken, "BilinearRescale") 
		|| (Temp != 0 && Temp != 1))
	SysErr("Can't read BilinearRescale from ASM file");
fBilinearRescale = (Temp == 1);

Fgets(s, SLEN-1, pAsmFile);
Temp = -1;
if (2 != sscanf(s, "%s %d", sToken, &Temp) || 0 != strcmp(sToken, "TrimCovar")
		|| Temp < 0 || Temp > 1000)	// 1000 is arbitrary
	SysErr("Can't read TrimCovar from ASM file");
nTrimCovar = Temp;

Fgets(s, SLEN-1, pAsmFile);
DTemp = -1.0;
if (2 != sscanf(s, "%s %g", sToken, &DTemp) || 0 != strcmp(sToken, "SigmoidScale")
		|| DTemp <= 0 || DTemp > 1e10)	// 1e10 is arbitrary
	SysErr("Can't read SigmoidScale from ASM file");
SigmoidScale = DTemp;

lprintf("\nnPoints %d nLevs %d PyrRatio %g PyrReduce %d fExplicitPrevNext %d NormalizedProfLen %g\n"
		"nScaledFaceWidth %d fBilinearRescale %d nTrimCovar %d SigmoidScale %g",
	nPoints, nLevs,	PyrRatio, PyrReduceMethod, fExplicitPrevNext, NormalizedProfLen,
	nScaledFaceWidth, fBilinearRescale, nTrimCovar, SigmoidScale);
}

//-----------------------------------------------------------------------------
// Read and check header of a level within an .asm file
// No results are kept, this just checks that the header is good

static void CheckAsmLevHeader (const char sAsmFile[], FILE *pAsmFile, int iLev)
{
// expect something like this: Lev 0 FirstShape xExtent 210 yExtent 243

char s[SLEN], sLev[SLEN], sFirstShape[SLEN];
int nLev;

if (!Fgets(s, SLEN-1, pAsmFile))
	SysErr("Can't read level %d header from %s", iLev, sAsmFile);