shapemodel.cpp

这是个人脸识别程序

// $masm\shapeModel.cpp 1.5 milbo$ routines for using the shape model during search
// Warning: this is raw research code -- expect it to be quite messy.
// milbo durban jun 06
//-----------------------------------------------------------------------------
// 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"

//-----------------------------------------------------------------------------
// Utilities for SortVec (which is defined below)

static const Vec *pgVec;

static int __cdecl CompareForSortVec (const void *pArg1, const void *pArg2)
{
int i1 = *(int *)pArg1;
int i2 = *(int *)pArg2;
double  diff = (*pgVec)(i1) - (*pgVec)(i2);
if (diff > 0)
	return -1;
else if (diff < 0)
	return 1;
else
	return 0;
}

//-----------------------------------------------------------------------------
// SortVec(v) returns the INDICES of v, sorted such that v(Sorted[i]) >= v(Sorted[i+1])
//
// v is unchanged
//
// Caller must free the returned array of indices

static void SortVec (int Sorted[],			// out
						const Vec &v)		// in
{
for (int i = 0; i < v.nelems(); i++)
	Sorted[i] = i;

pgVec = &v;

qsort(Sorted, v.nelems(), sizeof(int), CompareForSortVec);
}

//-----------------------------------------------------------------------------
// Set unused points to 0,0
// TODO look for a more efficient way of doing this

void SetUnusedPointsToZero (SHAPE &Shape, const SHAPE &RefShape)
{
for (int iPoint = 0; iPoint < RefShape.nrows(); iPoint++)
	if (!fPointUsed(RefShape, iPoint))
		{
		Shape(iPoint, VX) = 0;
		Shape(iPoint, VY) = 0;
		}
}

//-----------------------------------------------------------------------------
// Limit the values of b to make sure the generated shape is plausible.
// See CootesTaylor section 4.6 www.isbe.man.ac.uk/~bim/Models/app_models.pdf

static void nLimitB (Vec &b, 												// io
			 		const Vec EigVals, int nPrincipalEigs, double BMax)		// in
{
#if CONF_fClipB //---------------------------------------

// Method 1: apply hard limits i.e. clip each b[i] to BMax * sqrt(lambda_i)

for (int iEig = 0; iEig < nPrincipalEigs; iEig++)
	{
	double Limit = BMax * sqrt(EigVals(iEig));
	if (b(iEig) < -Limit)
		b(iEig) = -Limit;
	else if (b(iEig) > Limit)
		b(iEig) = Limit;
	}

#else //------------------------------------------------

// Method 2: constrain b to be in a hyperellipsoid
// i.e. scale each b[i] so sqrt(total) < BMax

double Limit = 0;
for (int iEig = 0; iEig < nPrincipalEigs; iEig++)
	Limit += (b(iEig) * b(iEig)) / EigVals(iEig);

Limit = sqrt(Limit);

if (Limit > BMax)
	for (iEig = 0; iEig < nPrincipalEigs; iEig++)
		b(iEig) *= BMax / Limit;		// pro-rated reduction

#if _DEBUG	// sanity check: recalculate with the new b and check that total < BMax
Limit = 0;
for (int i = 0; i < nPrincipalEigs; i++)
	Limit += (b(i) * b(i)) / EigVals(i);
if (sqrt(Limit) > BMax + 0.001)
	SysErr("LimitB %g %g", sqrt(Limit), BMax);
#endif

#endif //-----------------------------------------------

// discard unused higher components

while (iEig < EigVals.nrows())
	b(iEig++) = 0.0;
}

//-----------------------------------------------------------------------------
static SHAPE ConformShapeToModel_Simple (Vec &b,										// io
								const SHAPE &InShape, const tAsmModel &Model, int iLev,	// in
								const Vec *pWeights)									// in
{
Vec	MeanShape(Model.AsmLevs[iLev].MeanShape);
SetUnusedPointsToZero(MeanShape, InShape);

// For some calculations we need to see shapes (nrows x 2) as vectors (1 x 2*nrows).
// This is an alias to do this for MeanShape.
// Note that this is a "view" so if you change MeanShapeAsVec you are changing MeanShape too, and vice versa.

VecView	MeanShapeAsVec(MeanShape.viewAsCol());

// Find the Model b that best fits InShape to the Model.
// CootesTaylor 2004 Sec 7.2 says just one iteration is best, so there is no loop here.
// In fact, slightly improved results (especially for AR shapes) can be got
// by iterating, see ConformShapeToModel_Fancy

SHAPE x(InShape.nrows(), 2);
x.viewAsCol() = MeanShapeAsVec + Model.EigVecs * b;			// generate a model shape (this actually updates x)
Mat Pose(AlignShape(x, InShape, pWeights));					// find pose to best map the model shape x to InShape
SHAPE y(TransformShape(InShape, Pose.inverse()));			// project InShape into the model space
b = Model.EigInverse * (y.viewAsCol() - MeanShapeAsVec);	// update model params to match y
nLimitB(b, Model.AsmLevs[iLev].EigVals, Model.AsmLevs[iLev].nEigs, Model.AsmLevs[iLev].BMax);						// make sure we stay within model limits

// We now have b
// Generate OutShape from the model using our calculated b, and align OutShape to InShape

SHAPE OutShape(InShape.nrows(), 2);
OutShape.viewAsCol() = Model.EigVecs * b;

OutShape = TransformShape(Model.AsmLevs[iLev].MeanShape + OutShape, Pose);
SetUnusedPointsToZero(OutShape, InShape);
return OutShape;
}

//-----------------------------------------------------------------------------
// Like ConformShapeToModel_Simple but iterates looking for the best fit
// TODO need to add code to deal with unused points?

static SHAPE ConformShapeToModel_Fancy (Vec &b,												// io
								const SHAPE &InShape, const tAsmModel &Model, int iLev,		// in
								const Vec *pWeights, 										// in
								bool fShapeModelFinalIter)									// in
{
Vec	MeanShape(Model.AsmLevs[iLev].MeanShape);
SetUnusedPointsToZero(MeanShape, InShape);

// For some calculations we need to see shapes (nrows x 2) as vectors (1 x 2*nrows).
// This is an alias to do this for MeanShape.
// Note that this is a "view" so if you change MeanShapeAsVec you are changing MeanShape too, and vice versa.

VecView	MeanShapeAsVec(MeanShape.viewAsCol());

// find the Model b that best fits InShape to the Model

SHAPE Shape(InShape);						// working shape

#if CONF_fRetainOrgB
Vec OrgB(b);								// remember b passed in to this function
#endif

int nEigs = Model.AsmLevs[iLev].nEigs;
double BMax = Model.AsmLevs[iLev].BMax;
if (fShapeModelFinalIter)										// if final iter in main ASM search loop then loosen up the model
	{
	nEigs = Model.AsmLevs[iLev].nEigsFinal;
	BMax = Model.AsmLevs[iLev].BMaxFinal;
	}
for (int iter = 0; iter < Model.nShapeModelIters; iter++)
	{
	SHAPE x(InShape.nrows(), 2);
	x.viewAsCol() = MeanShapeAsVec + Model.EigVecs * b;			// generate a model shape (this actually updates x)
#if CONF_fAffineShapeModelAlign
	DASSERT(pWeights == NULL);
	Mat Pose(AlignShape_Affine(x, InShape));					// find pose to best map the model shape x to InShape
#elif CONF_fIterativeShapeModelAlign
	Mat Pose(AlignShape_Iteration(x, InShape, pWeights));
#else
	Mat Pose(AlignShape(x, InShape, pWeights));
#endif
	SHAPE y(TransformShape(InShape, Pose.inverse()));			// project InShape into the model space
	b = Model.EigInverse * (y.viewAsCol() - MeanShapeAsVec);	// update model params to match y
	nLimitB(b, Model.AsmLevs[iLev].EigVals, nEigs, BMax);		// make sure we stay within model limits

	// We now have b.
	// Generate Shape from the model using our calculated b, and align Shape to InShape.

	Shape.viewAsCol() = Model.EigVecs * b;
	Shape = TransformShape(Model.AsmLevs[iLev].MeanShape + Shape, Pose);
#if 0	//TODO put this into report
	#define CONF_fUseInShapeInFinalIter 0
	#define CONF_MinDist 				0		// 0.0 to use only model points, 1.0 to use InShape points whose dist < 1.0

	if (fShapeModelFinalIter && (CONF_MinDist != 0) && (CONF_fUseInShapeInFinalIter || iter < Model.nShapeModelIters-1))
		{
		// Replace points in Shape with the original InShape points if the InShape-model distance is small.
		// The idea is to prevent bad points from pulling the whole model off.

		Vec Dist(ShapeDist(Shape, InShape));	// distances between corresponding points in Shape and InShape
		int *iSortedIndices = (int *)malloc(Dist.nelems() * sizeof(int));
		SortVec(iSortedIndices, Dist);
		for (int i = 0; i < Shape.nrows(); i++)
			{
			if (Dist(i) < CONF_MinDist / pow(2, iLev))
				{
				int iRow = iSortedIndices[i];
				Shape(iRow, VX) = InShape(iRow, VX);
				Shape(iRow, VY) = InShape(iRow, VY);
				}
			}
		free(iSortedIndices);
		}
#endif
#if CONF_fRetainOrgB
	if (iter < Model.nShapeModelIters-1)
		b = OrgB;
#endif
	}
return Shape;
}

//-----------------------------------------------------------------------------
// Copy InShape and return it conformed to the model.
// In other words, generate a model shape that is as close as possible to InShape.
// InShape itself isn't changed.
//
// To match the model to Shape we need to find Pose and b in the following equation:
//
// Shape = Pose * (MeanShape + EigVecs * b)		(where = actually means approx equal)
//
// Note: The returned shape has all the points in the model, even if InShape
// has unused points (i.e. equal to 0,0)

SHAPE ConformShapeToModel (Vec &b,									// io
								const SHAPE &InShape, 				// in
								const tAsmModel &Model, int iLev,	// in
								const Vec *pWeights,				// in
								bool fShapeModelFinalIter)			// in
{
if (Model.fFancyShapeModel)
	return ConformShapeToModel_Fancy(b, InShape, Model, iLev, pWeights, fShapeModelFinalIter);
else
	return ConformShapeToModel_Simple(b, InShape, Model, iLev, pWeights);
}

//-----------------------------------------------------------------------------
void PrintEigInfo (bool fgBrief, const Vec &EigVals, int nPrincipalEigs)
{
double Total = 0;
for (int i = 0; i < nPrincipalEigs; i++)
	Total += EigVals(i);

bprintf(fgBrief, "Fraction of shape variation %.3f Ratio of first eig to last used eig %.3f\n",
	Total / EigVals.sum(), EigVals(0) / EigVals(nPrincipalEigs-1));
}