prof.cpp

这是个人脸识别程序

//		i.  call PrepareProf1D once with nProfWidth set to the length of the entire profile
//		ii. call Get1dProf multiple times, incrementing iOffset each time
//
// The two-step approach is for speed since we only have to prepare the profile
// once for multiple iOffset's.
// (We don't need the two-set approach for 2D profiles because we don't use whiskers
// for 2D profs and thus don't have to sample along a whisker).
//
// Shape is used for figuring out the direction of the whisker
// The whisker center point is x and y which may not be on the shape if there is
// orthogonal offset
// iOrthOffset is offset from the point orthogonally to the whisker

void PrepareProf1D (const Image &Img, const Mat &Shape, unsigned SubProfSpec, 	// in
					const tLand Lands[], 										// in
					const SHAPE &AlignedAvShape,								// in
					int iPoint, int nProfWidth, double ProfAngle,				// in
					int iOrthOffset)											// in
{
DASSERT((SubProfSpec & PROF_2d) == 0);
DASSERT(nProfWidth < CONF_nMaxProfWidth1D);

ngProfWidth = nProfWidth;
igProfPoint = iPoint;		// for sanity checking in Get1dProf later
pgProfImage = Img.buf;		// ditto
gSubProfSpec = SubProfSpec;	// ditto

int width = Img.width, height = Img.height;
int	nSamplePoints = (nProfWidth - 1) / 2; // number of profile sample points
										  // (nProfWidth is +-nSamplePoints and middle point)

double DeltaX;		// steps we move away along whisker from iPoint when sampling image
double DeltaY;		// DeltaX is along whisker, DeltaY is orthogonal to whisker
GetProfStepSize(&DeltaX, &DeltaY, Shape, iPoint, ProfAngle, Lands, AlignedAvShape);

double x = GetX(Shape(iPoint, VX), 0, iOrthOffset, DeltaX, DeltaY);
double y = GetY(Shape(iPoint, VY), 0, iOrthOffset, DeltaX, DeltaY);

// if whisker is vertical, PrevPix is below this pixel

double PrevPix = iGetPixel(Img, iRound(GetX(x, -nSamplePoints-1, 0, DeltaX, DeltaY)),
								iRound(GetY(y, -nSamplePoints-1, 0, DeltaX, DeltaY)));

if (iPoint == CONF_iProfTestPoint)
	lprintf("\nPrepareProf1D%d SubProfSpec %x width %d height %d ",
		CONF_iProfTestPoint, SubProfSpec, width, height);

unsigned TypeField = (SubProfSpec & PROF_TBits);

#if CONF_fDebug1dProfs
if (iPoint == 69 || iPoint == 77)
	lprintf("PrepareProf1D rawProf iPoint %2d:                ", iPoint);
#endif
for (int iSample = -nSamplePoints; iSample <= nSamplePoints; iSample++)
	{
	int ix = iRound(GetX(x, iSample, 0, DeltaX, DeltaY));
	int iy = iRound(GetY(y, iSample, 0, DeltaX, DeltaY));
	double Pix = iGetPixel(Img, ix, iy);
#if CONF_fDebug1dProfs
	if (iPoint == 69 || iPoint == 77)
		lprintf("%7.0f ", Pix);
#endif
	if (iPoint == CONF_iProfTestPoint && iSample > -3 && iSample < 3)
		lprintf("%s%d%+d=%.0f%s ", (iSample==0? "[":""), ix, iy, Pix, (iSample==0? "]":""));
#if CONF_fRawProfs
	gRawProf(0, iSample + nSamplePoints) = Pix;
#endif
	switch (TypeField)
		{
		case PROF_Grad:
			gProf(0, iSample + nSamplePoints) = Pix - PrevPix;
			break;
		case PROF_GradMagBelowRight:	// actually for 1d profs this is just the pixel on
										// the right, not pixels below-right
			gProf(0, iSample + nSamplePoints) = ABS(Pix - PrevPix);
			break;
		default:
			SysErr("PrepareProf1D TypeField %x", TypeField);
			break;
		}
	PrevPix = Pix;
	}
#if CONF_fDebug1dProfs
if (iPoint == 69 || iPoint == 77)
	lprintf("\n");
#endif
if (iPoint == CONF_iProfTestPoint)
	{
	GetProfStepSize(&DeltaX, &DeltaY, Shape, iPoint, ProfAngle, Lands, AlignedAvShape);
	lprintf("\nPrepareProf1D%d width %d XY %g %g DeltaX %.4f DeltaY %.4f\n",
		CONF_iProfTestPoint, width, Shape(iPoint, VX), Shape(iPoint, VY), DeltaX, DeltaY);
	}
}

//-----------------------------------------------------------------------------
// Sum==0 happens on the shirt of BioId B1274.jpg for example -- the shirt
// is uniform so all deltas are 0

static inline bool fDoubleIsZero (double Double, Vec &Prof)
{
if (fEqual(Double, 0))
	{
#if CONF_fPrintProfileWarn
	Warn("Double==0 %s", pgCurImageName);
	Prof.print("Prof ");
#endif
	return true;
	}
return false;
}

//-----------------------------------------------------------------------------
// Return true if profile normalization causes linear dependence between profiles.
//
// This is the case if we normalize all profiles to the same length i.e. divide by sum(Prof)
// This is NOT the case (nearly always) if we normalize by dividing by absSum(Prof).
//
// See Tatsuoka "Multivariate Analysis" p148
// Using linearly dependent profiles causes a singular covar matrix, which is a problem.
// We solve this problem using CONF_fDropEigToPreventSingCovar or
// CONF_fDropProfToStopColinearity

bool fProfileNormalizationCausesLinearDependence (unsigned SubProfSpec)
{
unsigned NormalizationField = (SubProfSpec & PROF_NormalizationField);

return NormalizationField == PROF_Flat;
	//TODO this falsely returns true for certain types of 1D profs,
	//   but it doesn't seem to matter
}

//-----------------------------------------------------------------------------
// This supports the following normalization method:
//		PROF_Flat with TBits
//			PROF_Grad 				(normalize using abs values)
//			PROF_GradMagBelowRight	(normalize using actual values - save time but using abs vals)
//
// It once supported PROF_SigmAbsSum but this gave worse results so I removed it.

static void Normalize1d (Vec &Prof, 				// io
						unsigned SubProfSpec)		// in
{
DASSERT((SubProfSpec & PROF_2d) == 0);

switch (SubProfSpec & PROF_NormalizationField)		// normalization field
	{
	case PROF_Flat:
		{
		double Sum;
		switch (SubProfSpec & PROF_TBits)	// type field
			{
			case PROF_Grad:					// profile values can be negative or positive
											// so must use absSum()
				Sum = Prof.absSum();
				break;
			case PROF_GradMagBelowRight: 	// profile values are all positive
										 	//  so faster to use sum()
				Sum = Prof.sum();
				break;
			default:
				SysErr("Get1dProf %x", SubProfSpec);
			}
		Sum /= Prof.ncols() * gNormalizedProfLen;
		if (!fDoubleIsZero(Sum, Prof))
			Prof /= Sum;
		break;
		}
	default:
		SysErr("Get1dProf %x", SubProfSpec);
	}
}

//-----------------------------------------------------------------------------
// This supports the following normalization methods:
//
//		PROF_Flat
//		PROF_SigmAbsSum
//		PROF_SigmLen
//
// This routine needs to be fast.

static void Normalize2d (Vec &Prof, 						// io
						 const unsigned SubProfSpec)		// in
{
extern double SigmoidScaleGlobal;
DASSERT(SubProfSpec & PROF_2d);
const int nelems = Prof.nelems();
static bool fPrinted = false;
int iCol = nelems;
double * const pData = Prof.m->data;			// for efficiency, access mat buf directly

#if 0	// for debugging, to look at prof before and after
if (!fPrinted)
	{ Prof.print("\nx <- c(", "%g,"); lprintf(")\n"); }
#endif

switch (SubProfSpec & PROF_NormalizationField)	// normalization field
	{
	case PROF_Flat:
		{
		const double Len = Prof.len();
		if (Len <= 0)							// prevent div by zero
			{
			if (!fPrinted)
				{
				Warn("Normalize2d: Len <= 0");
				fPrinted = true;
				}
			}
		else
			{
			const double Scale = nelems * gNormalizedProfLen / Len;
			while (--iCol >= 0)
				pData[iCol] *= Scale;
			}
		break;
		}
	case PROF_SigmAbsSum:
		{
		// Note that we scale but don't center here because
		// we don't want to center an all positive profile

		const double AbsSum = Prof.absSum();
		if (AbsSum <= 0)
			{
			if (!fPrinted)
				{
				Warn("Normalize2d: AbsSum <= 0");
				fPrinted = true;
				}
			}
		else
			{
			const double Scale = nelems * gNormalizedProfLen / AbsSum;
			while (--iCol >= 0)
				{
				pData[iCol] *= Scale;
				pData[iCol] /= (fabs(pData[iCol]) + SigmoidScaleGlobal / 10);
				}
			}
 		break;
		}
	case PROF_SigmLen:
		{
		// Once again, scale but not center

		const double Len = Prof.len();
		if (Len <= 0)
			Warn("Normalize2d: Len <= 0");
		else
			{
			const double Scale = nelems * gNormalizedProfLen / Len;
			while (--iCol >= 0)
				{
				pData[iCol] *= Scale;
				pData[iCol] /= (fabs(pData[iCol]) + SigmoidScaleGlobal);
				}
			}
 		break;
		}
	default:
		SysErr("Normalize2d %x", SubProfSpec);
		break;
	}
#if 0
if (!fPrinted)
	{ Prof.print("y <- c(", "%g,"); lprintf(")\n"); fPrinted = true; }
#endif
}

//-----------------------------------------------------------------------------
// Use this after preparing global vars by calling PrepareProf1D
// iOffset is offset along whisker from point, can be +ve or -ve

void Get1dProf (Vec &Prof,												// out
				 unsigned SubProfSpec, const Image &Img, int iPoint, 	// in
				 int iOffset, int iRowInGlobal, bool fDumpRawProf)		// in
{
DASSERT((SubProfSpec & PROF_2d) == 0);
DASSERT((SubProfSpec & PROF_FBit) == 0);
DASSERT(iPoint == igProfPoint);
DASSERT(Img.buf == pgProfImage);
DASSERT(SubProfSpec == gSubProfSpec);

int nProfWidth = Prof.ncols();		// +-nSamplePoints and middle point
Prof = gProf.view(iRowInGlobal, ngProfWidth/2 + iOffset - nProfWidth/2, 1, nProfWidth);

#if CONF_fDebug1dProfs && 0
if (iPoint == 69 || iPoint == 77)
	{
	char s[SLEN]; sprintf(s, "Get1dProf iPoint %2d:                              ", iPoint);
	Prof.print(s);
	}
#endif
#if CONF_fRawProfs
if (fDumpRawProf)
	{
	if (!pgRawProfFile)
		pgRawProfFile = Fopen(CONF_sRawProfFile, "w");
	Fprintf(pgRawProfFile,
		"# iPoint %d iOffset %d iRowInGlobal %d\n", iPoint, iOffset, iRowInGlobal);
	gRawProf.view(iRowInGlobal, ngProfWidth/2 + iOffset - nProfWidth/2,
		1, nProfWidth).write(CONF_sRawProfFile, pgRawProfFile);
	}
#endif

Normalize1d(Prof, SubProfSpec);
}

//-----------------------------------------------------------------------------
// Do sigmoid normalization on x and y gradients
// as per ScottCootesTaylor "Improving Appearance Model..."

#if 0 //TODO not yet used because we don't have a mechanism for
	  //   passing more than one grad to the normalization routines
static void NormalizeXYGrads (VecView &xGrad, VecView &yGrad)
{
DASSERT(xGrad.ncols() == yGrad.ncols());

double xLen = xGrad.sum();
double yLen = yGrad.sum();
const double TotalLen = sqrt(xLen * xLen + yLen * yLen);
const double NormFactor = xGrad.nelems() * gNormalizedProfLen;

for (int iCol = 0; iCol < xGrad.ncols(); iCol++)
	{
	double Len = sqrt(xGrad(0, iCol) * xGrad(0, iCol) + yGrad(0, iCol) * yGrad(0, iCol));
	xGrad(0, iCol) *= NormFactor / (Len + TotalLen);
	yGrad(0, iCol) *= NormFactor / (Len + TotalLen);
	}
}
#endif

//-----------------------------------------------------------------------------
static Mat *pGetProfWindow (unsigned SubProfSpec, int nProfWidth)
{
static Mat Window;	// cache for speed: use last window if possible, since
					// window type seldom changes
static int nLastProfWidth = 0;
static unsigned LastWindowType = unsigned(-1);

unsigned WindowType = (SubProfSpec & PROF_WindowField);

if (nProfWidth != nLastProfWidth || WindowType != LastWindowType)
	{
	// first time for this nProfWidth and WindowType, so initialize Window

	nLastProfWidth = nProfWidth;
	LastWindowType = WindowType;
	int nProfWidthEachSide = (nProfWidth - 1) / 2;
	Window.dim(nProfWidth, nProfWidth);
	for (int iy = -nProfWidthEachSide; iy <= nProfWidthEachSide; iy++)
		for (int ix = -nProfWidthEachSide; ix <= nProfWidthEachSide; ix++)
			switch (WindowType)
				{
				case PROF_WindowEquallyWeighted:
					Window(iy+nProfWidthEachSide, ix+nProfWidthEachSide) = 1;
					break;
				case PROF_WindowGaussian:
					{
					// Window will look like this (but divide the numbers shown here by 100):
					//
					//    24  30  35  41  45  48  49  48  45  41  35  30  24
					//    30  37  44  51  56  59  61  59  56  51  44  37  30
					//    35  44  53  61  67  71  73  71  67  61  53  44  35
					//    41  51  61  70  77  82  84  82  77  70  61  51  41
					//    45  56  67  77  85  90  92  90  85  77  67  56  45
					//    48  59  71  82  90  96  98  96  90  82  71  59  48
					//    49  61  73  84  92  98 100  98  92  84  73  61  49
					//    48  59  71  82  90  96  98  96  90  82  71  59  48
					//    45  56  67  77  85  90  92  90  85  77  67  56  45
					//    41  51  61  70  77  82  84  82  77  70  61  51  41
					//    35  44  53  61  67  71  73  71  67  61  53  44  35
					//    30  37  44  51  56  59  61  59  56  51  44  37  30
					//    24  30  35  41  45  48  49  48  45  41  35  30  24

					const double Width2 = (nProfWidthEachSide-1) * (nProfWidthEachSide-1);
					Window(iy+nProfWidthEachSide, ix+nProfWidthEachSide) =
						exp(-((ix * ix) + (iy * iy))/(2 * Width2));
					break;
					}
				case PROF_WindowCircle:
					{
					// Window will look like this:
					//
					// 0 0 0 0 0 1 1 1 0 0 0 0 0
					// 0 0 0 1 1 1 1 1 1 1 0 0 0
					// 0 0 1 1 1 1 1 1 1 1 1 0 0
					// 0 1 1 1 1 1 1 1 1 1 1 1 0
					// 0 1 1 1 1 1 1 1 1 1 1 1 0
					// 1 1 1 1 1 1 1 1 1 1 1 1 1
					// 1 1 1 1 1 1 1 1 1 1 1 1 1
					// 1 1 1 1 1 1 1 1 1 1 1 1 1
					// 0 1 1 1 1 1 1 1 1 1 1 1 0
					// 0 1 1 1 1 1 1 1 1 1 1 1 0
					// 0 0 1 1 1 1 1 1 1 1 1 0 0
					// 0 0 0 1 1 1 1 1 1 1 0 0 0
					// 0 0 0 0 0 1 1 1 0 0 0 0 0

					const double Width2 = nProfWidthEachSide * nProfWidthEachSide;
					if (ix * ix + iy * iy <= Width2 + 1)
						Window(iy+nProfWidthEachSide, ix+nProfWidthEachSide) = 1;	// TODO
					else
						Window(iy+nProfWidthEachSide, ix+nProfWidthEachSide) = 0.001;
									//TODO can't actually use 0 else covar matrix will be singular
					break;