/*------------------------------------------------------------------------------*/
/** 
 *  \file   GW_VoronoiMesh.cpp
 *  \brief  Definition of class \c GW_VoronoiMesh
 *  \author Gabriel Peyré
 *  \date   4-12-2003
 */ 
/*------------------------------------------------------------------------------*/


#ifdef GW_SCCSID
    static const char* sccsid = "@(#) GW_VoronoiMesh.cpp(c) Gabriel Peyré2003";
#endif // GW_SCCSID

#include "stdafx.h"
#include "GW_VoronoiMesh.h"

#ifndef GW_USE_INLINE
    #include "GW_VoronoiMesh.inl"
#endif

using namespace GW;

T_VoronoiVertexList GW_VoronoiMesh::CurrentTargetVertex_;
GW_VoronoiVertex* GW_VoronoiMesh::pCurrentVoronoiVertex_ = NULL;
T_VoronoiVertexMap GW_VoronoiMesh::VoronoiVertexMap_;
T_VertexPathMultiMap GW_VoronoiMesh::BoundaryEdgeMap_;

GW_U32 GW_VoronoiMesh::nNbrBaseVertex_RD_	= 0;
T_FloatMap* GW_VoronoiMesh::pCurWeights_	= NULL;

/*------------------------------------------------------------------------------*/
// Name : GW_VoronoiMesh::PerformFastMarching
/**
 *  \param  Mesh [GW_GeodesicMesh&] The mesh.
 *  \author Gabriel Peyré
 *  \date   4-13-2003
 * 
 *  Just an helper method that launch a fire from each base point.
 */
/*------------------------------------------------------------------------------*/
void GW_VoronoiMesh::PerformFastMarching( GW_GeodesicMesh& Mesh, T_GeodesicVertexList& VertList )
{
	/* add the points as starting points */
	Mesh.ResetGeodesicMesh();
	for( IT_GeodesicVertexList it = VertList.begin(); it!=VertList.end(); ++it )
	{
		GW_GeodesicVertex* pVert = *it;
		GW_ASSERT( pVert!=NULL );
		Mesh.AddStartVertex( *pVert );
	}
	/* perform a firestart */
	Mesh.PerformFastMarching(); 
}

/*------------------------------------------------------------------------------*/
// Name : GW_VoronoiMesh::FastMarchingCallbackFunction_VertexInsersion
/**
 *  \param  CurVert [GW_GeodesicVertex&] The current vertex.
 *  \param  rNewDist [GW_Float] The new distance.
 *  \return [GW_Bool] Insert this vertex to active list ?
 *  \author Gabriel Peyré
 *  \date   5-13-2003
 * 
 *  Add the vertex only if the distance decreases.
 */
/*------------------------------------------------------------------------------*/
GW_Bool GW_VoronoiMesh::FastMarchingCallbackFunction_VertexInsersion( GW_GeodesicVertex& CurVert, GW_Float rNewDist )
{
	return CurVert.GetDistance() > rNewDist;
}

/*------------------------------------------------------------------------------*/
// Name : GW_VoronoiMesh::AddFurthestPoint
/**
*  \param  Mesh [GW_GeodesicMesh&] Fine mesh.
*  \author Gabriel Peyré
*  \date   4-14-2003
* 
*  Add to the list of point the furthest from all nodes.
*/
/*------------------------------------------------------------------------------*/
GW_U32 GW_VoronoiMesh::AddFurthestPoint( T_GeodesicVertexList& VertList,GW_GeodesicMesh& Mesh, GW_Bool bUseRandomStartVertex )
{
	if( VertList.empty() )
	{
		Mesh.ResetGeodesicMesh();
		/* choose 1 point at random */
		GW_GeodesicVertex* pStartVertex = (GW_GeodesicVertex*) Mesh.GetRandomVertex();
		GW_ASSERT( pStartVertex!=NULL );
		if( pStartVertex==NULL )
			return 0;		
		VertList.push_back( pStartVertex );
		if( !bUseRandomStartVertex )
		{
			/* find furthest point from this one */
			GW_VoronoiMesh::PerformFastMarching( Mesh, VertList );
			GW_GeodesicVertex* pSelectedVert = GW_VoronoiMesh::FindMaxVertex( Mesh );
			GW_ASSERT( pSelectedVert!=NULL );
			VertList.push_back( pSelectedVert );
			VertList.pop_front();	// remove random choosen vertex
			Mesh.ResetGeodesicMesh();
		}
	}
	else
	{
		Mesh.RegisterVertexInsersionCallbackFunction( GW_VoronoiMesh::FastMarchingCallbackFunction_VertexInsersion );

		/* perform marching */
		GW_VoronoiMesh::ResetOnlyVertexState( Mesh );
		GW_GeodesicVertex* pStartVert = VertList.back();
		GW_ASSERT( pStartVert!=NULL );
		Mesh.PerformFastMarching( pStartVert  );

		/* find intersection points */
		GW_GeodesicVertex* pSelectedVert = GW_VoronoiMesh::FindMaxVertex( Mesh );
		GW_ASSERT( pSelectedVert!=NULL );
		/* find the triangle where the maximum occurs */
#if 0
		GW_Face* pMaxFace = GW_VoronoiMesh::FindMaxFace( *pSelectedVert );
		GW_ASSERT( pMaxFace!=NULL );
		/* this is a critical point */
		GW_Float x=1.0f/3, y=1.0f/3, z=1.0f/3;
		GW_GeodesicVertex* pNewVert = (GW_GeodesicVertex*) Mesh.InsertVertexInFace( *pMaxFace, x,y,z );
		pNewVert->SetFront( pSelectedVert->GetFront() );
		pNewVert->SetDistance( pSelectedVert->GetDistance() );
#else
		GW_GeodesicVertex* pNewVert = pSelectedVert;
#endif
		/* compute the real point */
		VertList.push_back( pNewVert );

		Mesh.RegisterVertexInsersionCallbackFunction( NULL );
	}

	return 1;
}

/*------------------------------------------------------------------------------*/
// Name : GW_VoronoiMesh::FindMaxVertex
/**
 *  \param  Mesh [GW_GeodesicMesh&] The mesh.
 *  \return [GW_GeodesicVertex*] The max distance vertex.
 *  \author Gabriel Peyré
 *  \date   5-14-2003
 * 
 *  Find the vertex with maximum distance.
 */
/*------------------------------------------------------------------------------*/
GW_GeodesicVertex* GW_VoronoiMesh::FindMaxVertex( GW_GeodesicMesh& Mesh )
{
	GW_GeodesicVertex* pSelectedVert = NULL;
	GW_Float rMaxDist = -GW_INFINITE;
	for( GW_U32 i=0; i<Mesh.GetNbrVertex(); ++i )
	{
		GW_GeodesicVertex* pVert = (GW_GeodesicVertex*) Mesh.GetVertex( i );
		GW_ASSERT( pVert!=NULL );
		if( pVert->GetFace()!=NULL && pVert->GetDistance()>rMaxDist )
		{
			rMaxDist = pVert->GetDistance();
			pSelectedVert = pVert;
		}
	}
	return pSelectedVert;
}


/*------------------------------------------------------------------------------*/
// Name : GW_VoronoiMesh::FindMaxFace
/**
 *  \param  Vert [GW_GeodesicVertex&] The current vertex.
 *  \return [GW_Face*] The face.
 *  \author Gabriel Peyré
 *  \date   4-15-2003
 * 
 *  Find the face of maximum distance.
 */
/*------------------------------------------------------------------------------*/
GW_Face* GW_VoronoiMesh::FindMaxFace( GW_GeodesicVertex& Vert )
{
	GW_Float rBestDistance = -GW_INFINITE;
	GW_Face* pCurFace_ = NULL;
	for( GW_VertexIterator it = Vert.BeginVertexIterator(); it!=Vert.EndVertexIterator(); ++it )
	{
		GW_GeodesicVertex* pVert = (GW_GeodesicVertex*)  *it;
		if( pVert->GetDistance()>rBestDistance )
		{
			rBestDistance = pVert->GetDistance();
			GW_GeodesicVertex* pVert1 = (GW_GeodesicVertex*) it.GetLeftVertex();
			GW_GeodesicVertex* pVert2 = (GW_GeodesicVertex*) it.GetRightVertex();
			if( pVert1!=NULL && pVert2!=NULL )
			{
				if( pVert1->GetDistance()<pVert2->GetDistance() )
					pCurFace_ = it.GetLeftFace();
				else
					pCurFace_ = it.GetRightFace();
			}
			else if( pVert1!=NULL )
			{
				pCurFace_ = it.GetLeftFace();
			}
			else
			{
				GW_ASSERT( pVert2!=NULL );
				pCurFace_ = it.GetRightFace();
			}
		}
	}	
	return pCurFace_;
}


/*------------------------------------------------------------------------------*/
// Name : GW_VoronoiMesh::AddFurthestPointsIterate
/**
 *  \param  Mesh [GW_GeodesicMesh&] The mesh.
 *  \param  nNbrIterations [GW_U32] Number of points.
 *  \author Gabriel Peyré
 *  \date   4-14-2003
 * 
 *  Add a given number of furthest points.
 */
/*------------------------------------------------------------------------------*/
GW_U32 GW_VoronoiMesh::AddFurthestPointsIterate( T_GeodesicVertexList& VertList,GW_GeodesicMesh& Mesh, GW_U32 nNbrIterations, GW_Bool bUseRandomStartVertex, GW_Bool bUseProgressBar )
{
	GW_U32 nNbrPoints = 0;
	GW_ProgressBar pb;
	if( bUseProgressBar )
		pb.Begin();
	for( GW_U32 i=0; i<nNbrIterations; ++i )
	{
		nNbrPoints += GW_VoronoiMesh::AddFurthestPoint( VertList, Mesh, bUseRandomStartVertex );
		if( bUseProgressBar )
			pb.Update(((GW_Float) i+1)/((GW_Float) nNbrIterations));
	}
	if( bUseProgressBar )
	{
		pb.End();
		cout << endl;
	}
	return nNbrPoints;
}




/*------------------------------------------------------------------------------*/
// Name : GW_VoronoiMesh::FastMarchingCallbackFunction_MeshBuilding
/**
*  \param  CurVert [GW_GeodesicVertex&] The current vertex.
*  \return [GW_Bool] The new dead vertex.
*  \author Gabriel Peyré
*  \date   5-13-2003
* 
*  Test if the vertex is a saddle point. 
*/
/*------------------------------------------------------------------------------*/
void GW_VoronoiMesh::FastMarchingCallbackFunction_MeshBuilding( GW_GeodesicVertex& CurVert )
{
	GW_GeodesicVertex* pFront = CurVert.GetFront();
	GW_ASSERT( pFront!=NULL );
	/* retrieve the voronoi vertex corresponding to the front */
	GW_VoronoiVertex* pVoronoiVert0 = GW_VoronoiMesh::GetVoronoiFromGeodesic( *pFront );
	GW_ASSERT( pVoronoiVert0!=NULL );
	/* test if this point is a saddle point */
	for( GW_VertexIterator it=CurVert.BeginVertexIterator(); it!=CurVert.EndVertexIterator(); ++it )
	{
		GW_GeodesicVertex* pNeighborVert = (GW_GeodesicVertex*) *it;
		GW_GeodesicVertex* pNeighborFront = pNeighborVert->GetFront();
		if( pNeighborFront!=NULL && pNeighborFront!=pFront )
		{
			/* that's it ! */
			GW_VoronoiVertex* pVoronoiVert1 = GW_VoronoiMesh::GetVoronoiFromGeodesic( *pNeighborFront );
			GW_ASSERT( pVoronoiVert1!=NULL );
			if( !pVoronoiVert0->IsNeighbor(*pVoronoiVert1) )
			{
				GW_ASSERT( ! pVoronoiVert1->IsNeighbor(*pVoronoiVert0) );
				/* test for manifold structure for the future edge */
				GW_U32 nNumTriangle = 0;
				for( IT_VoronoiVertexList itVoronoi=pVoronoiVert1->BeginNeighborIterator(); itVoronoi!=pVoronoiVert1->EndNeighborIterator(); ++itVoronoi )
				{
					GW_VoronoiVertex* pVoronoiVert2 = *itVoronoi;
					GW_ASSERT( pVoronoiVert2!=NULL );
					if( pVoronoiVert2!=pVoronoiVert0 && pVoronoiVert0->IsNeighbor(*pVoronoiVert2) )
						nNumTriangle++;
				}
				if( nNumTriangle<=2 )
				{
					pVoronoiVert0->AddNeighbor( *pVoronoiVert1 );
					pVoronoiVert1->AddNeighbor( *pVoronoiVert0 );
				}
				/* test for manifold structure on the newly created edges */
				for( IT_VoronoiVertexList itVoronoi=pVoronoiVert1->BeginNeighborIterator(); itVoronoi!=pVoronoiVert1->EndNeighborIterator(); ++itVoronoi )
				{
					GW_VoronoiVertex* pVoronoiVert2 = *itVoronoi;
					GW_ASSERT( pVoronoiVert2!=NULL );
					if( pVoronoiVert2!=pVoronoiVert0 && pVoronoiVert0->IsNeighbor(*pVoronoiVert2) )
					{
						/* we are on a newly created face : test for structure integrity */
						if( !GW_VoronoiMesh::TestManifoldStructure(*pVoronoiVert0, *pVoronoiVert2) )
						{
							pVoronoiVert0->RemoveNeighbor( *pVoronoiVert1 );
							pVoronoiVert1->RemoveNeighbor( *pVoronoiVert0 );
							break;
						}
						if( !GW_VoronoiMesh::TestManifoldStructure(*pVoronoiVert1, *pVoronoiVert2) )
						{
							pVoronoiVert0->RemoveNeighbor( *pVoronoiVert1 );
							pVoronoiVert1->RemoveNeighbor( *pVoronoiVert0 );
							break;
						}
					}
				}
			}
			else
			{
				/* the graph must be undirected */
				GW_ASSERT( pVoronoiVert1->IsNeighbor(*pVoronoiVert0) );
			}
		}
	}
}


/*------------------------------------------------------------------------------*/
// Name : GW_VoronoiMesh::TestManifoldStructure
/**
 *  \param  Vert0 [GW_VoronoiVertex&] 1st vertex.
 *  \param  Vert1 [GW_VoronoiVertex&] 2nd vertex.
 *  \return [GW_Bool] Answer.
 *  \author Gabriel Peyré
 *  \date   5-15-2003
 * 
 *  Test if the edge has only 2 neighbors.
 */
/*------------------------------------------------------------------------------*/
GW_Bool GW_VoronoiMesh::TestManifoldStructure( GW_VoronoiVertex& Vert0, GW_VoronoiVertex& Vert1 )
{
	GW_U32 nNumTriangle = 0;
	for( IT_VoronoiVertexList itVoronoi=Vert1.BeginNeighborIterator(); itVoronoi!=Vert1.EndNeighborIterator(); ++itVoronoi )
	{
		GW_VoronoiVertex* pVoronoiVert2 = *itVoronoi;
		GW_ASSERT( pVoronoiVert2!=NULL );
		if( pVoronoiVert2!=&Vert0 && Vert0.IsNeighbor(*pVoronoiVert2) )
			nNumTriangle++;
	}
	return nNumTriangle<=2;
}

/*------------------------------------------------------------------------------*/
// Name : GW_VoronoiMesh::CreateVoronoiVertex
/**
 *  \author Gabriel Peyré
 *  \date   6-1-2003
 * 
 *  Just create the voronoi vertex.
 */
/*------------------------------------------------------------------------------*/
void GW_VoronoiMesh::CreateVoronoiVertex()
{
	/* Create Vornoi vertex and make the inverse map GeodesicVertex->VoronoiVertex */
	this->SetNbrVertex( (GW_U32) BaseVertexList_.size() );
	GW_U32 nNum = 0;
	for( IT_GeodesicVertexList it = BaseVertexList_.begin(); it!=BaseVertexList_.end(); ++it )
	{
		GW_GeodesicVertex* pVert = *it;
		GW_ASSERT( pVert!=NULL );
		/* create a new voronoi vertex */
		GW_VoronoiVertex* pVornoiVert = new GW_VoronoiVertex;
		pVornoiVert->SetBaseVertex( *pVert );
		this->SetVertex( nNum, pVornoiVert );
		nNum++;
		/* build the inverse map */
		GW_ASSERT( VoronoiVertexMap_.find( pVert->GetID() )==VoronoiVertexMap_.end() );
		VoronoiVertexMap_[ pVert->GetID() ] = pVornoiVert;
	}
}


/*------------------------------------------------------------------------------*/
// Name : GW_VoronoiMesh::BuildMesh
/**
 *  \param  Mesh [GW_Mesh&] The mesh that will be build.
 *  \author Gabriel Peyré
 *  \date   4-13-2003
 * 
 *  Build a mesh from scratch.
 */
/*------------------------------------------------------------------------------*/
void GW_VoronoiMesh::BuildMesh( GW_GeodesicMesh& Mesh, GW_Bool bFixHole )
{
	/* Create Vornoi vertex and make the inverse map GeodesicVertex->VoronoiVertex */
	this->CreateVoronoiVertex();

#if 1	// simple method

	GW_OutputComment("Recomputing the whole Voronoi diagram.");
	GW_VoronoiMesh::PerformFastMarching( Mesh, BaseVertexList_ );

	/* find the faces */
	GW_OutputComment("Building the faces.");
	T_FaceMap FaceMap;	// to store the faces already built.
	for( GW_U32 i=0; i<Mesh.GetNbrFace(); ++i )
	{
		GW_Face* pFace = Mesh.GetFace(i);	GW_ASSERT( pFace!=NULL );
		GW_GeodesicVertex* pGeo[3];
		GW_GeodesicVertex* pFront[3];
		for( GW_U32 j=0; j<3; ++j )
		{
			pGeo[j] = (GW_GeodesicVertex*) pFace->GetVertex(j); GW_ASSERT( pGeo[j]!=NULL );
			pFront[j] = pGeo[j]->GetFront();
		}
		if( pFront[0]!=pFront[1] && pFront[1]!=pFront[2] && pFront[2]!=pFront[0] )
		{
			GW_U32 nID = GW_Vertex::ComputeUniqueId( *pFront[0], *pFront[1], *pFront[2] );
			if( FaceMap.find(nID)==FaceMap.end() )
			{
				/* create the face */
				GW_Face& Face = this->CreateNewFace();
				FaceMap[nID] = &Face;
				for( GW_U32 j=0; j<3; ++j )		// assign the vertices
				{
					GW_U32 nID = pFront[j]->GetID();
					GW_ASSERT( VoronoiVertexMap_.find(nID)!=VoronoiVertexMap_.end() );
					GW_VoronoiVertex* pVorVert = VoronoiVertexMap_[nID];	GW_ASSERT( pVorVert!=NULL );
					Face.SetVertex( *pVorVert, j );
				}
			}
		}
	}

#else

	GW_OutputComment("Computing voronoi diagrams.");
	/* perform once more a firestart to set up connectivity */		
	Mesh.RegisterNewDeadVertexCallbackFunction( GW_VoronoiMesh::FastMarchingCallbackFunction_MeshBuilding );
	Mesh.ResetGeodesicMesh();
	GW_VoronoiMesh::PerformFastMarching( Mesh, BaseVertexList_ );
	Mesh.RegisterNewDeadVertexCallbackFunction( NULL );

	/* build the faces */
	T_FaceMap FaceMap;	// to store the faces already built.

	GW_OutputComment("Building voronoi mesh faces.");
	for( IT_GeodesicVertexList it = BaseVertexList_.begin(); it!=BaseVertexList_.end(); ++it )
	{
		GW_GeodesicVertex* pVert0 = *it;
		GW_ASSERT( pVert0!=NULL );
		/* retrive the corresponding voronoi vertex */
		GW_VoronoiVertex* pVoronoiVert0 = GW_VoronoiMesh::GetVoronoiFromGeodesic( *pVert0 );
		GW_ASSERT( pVoronoiVert0!=NULL );
		for( IT_VoronoiVertexList itVoronoi1=pVoronoiVert0->BeginNeighborIterator(); itVoronoi1!=pVoronoiVert0->EndNeighborIterator(); ++itVoronoi1 )
		{
			GW_VoronoiVertex* pVoronoiVert1 = *itVoronoi1;
			GW_ASSERT( pVoronoiVert1!=NULL );
			GW_U32 nNumTriangle = 0;
			for( IT_VoronoiVertexList itVoronoi2=pVoronoiVert1->BeginNeighborIterator(); itVoronoi2!=pVoronoiVert1->EndNeighborIterator(); ++itVoronoi2 )
			{
				GW_VoronoiVertex* pVoronoiVert2 = *itVoronoi2;
				GW_ASSERT( pVoronoiVert2!=NULL );
				if( pVoronoiVert2!=pVoronoiVert0 && pVoronoiVert0->IsNeighbor(*pVoronoiVert2) )
				{
					/* yes, we find a triangle ! Test if it wasn't already constructed */
					GW_U32 nUniqueId = GW_Vertex::ComputeUniqueId( *pVoronoiVert0, *pVoronoiVert1, *pVoronoiVert2 );
					if( FaceMap.find(nUniqueId)==FaceMap.end() )
					{
						nNumTriangle++;
						GW_ASSERT( nNumTriangle<=2 );	// assert manifold structure
						/* this is the 1st time we encounter this face. */
						GW_Face* pFace = &Mesh.CreateNewFace();
						/* set up the face */
						pFace->SetVertex( *pVoronoiVert0, *pVoronoiVert1, *pVoronoiVert2 );
						FaceMap[nUniqueId] = pFace;
					}
				}
			}
		}
	}

#endif

	
	/* assign the faces */
	this->SetNbrFace( (GW_U32) FaceMap.size() );
	GW_U32 nNum = 0;
	for( IT_FaceMap it = FaceMap.begin(); it!=FaceMap.end(); ++it )
	{
		this->SetFace( nNum, it->second );
		nNum++;
	}
	/* rebuild connectivity */
	GW_OutputComment("Building connectivity.");
	this->BuildConnectivity();
	/* try to fill the holes */
	if( bFixHole )
	{
		GW_OutputComment("Fixing holes.");
		this->FixHole();
		/* re-rebuild connectivity */
		this->BuildConnectivity();
	}

}

/*------------------------------------------------------------------------------*/
// Name : GW_VoronoiMesh::FixHole
/**
 *  \author Gabriel Peyré
 *  \date   5-16-2003
 * 
 *  Try to fix the hole in the boundary.
 */
/*------------------------------------------------------------------------------*/
void GW_VoronoiMesh::FixHole()
{
	typedef std::pair<GW_VoronoiVertex*, GW_VoronoiVertex*> T_VertexPair;
	typedef std::list<T_VertexPair>	T_VertexPairList;
	typedef T_VertexPairList::iterator	IT_VertexPairList;
	T_VertexPairList VertexPairList;
	for( GW_U32 i=0; i<this->GetNbrFace(); ++i )
	{
		GW_Face* pFace = this->GetFace(i);
		GW_ASSERT( pFace!=NULL );
		for( GW_U32 nV = 0; nV<3; ++nV )
		{
			if( pFace->GetFaceNeighbor(nV)==NULL )
			{
				GW_VoronoiVertex* pVert1 = (GW_VoronoiVertex*) pFace->GetVertex( (nV+1)%3 );
				GW_VoronoiVertex* pVert2 = (GW_VoronoiVertex*) pFace->GetVertex( (nV+2)%3 );
				VertexPairList.push_back( T_VertexPair(pVert1,pVert2) );
			}
		}
	}
	char str[50];
	sprintf( str, "%d boundary edges detected.", VertexPairList.size() );
	GW_OutputComment( str );
	while( !VertexPairList.empty() )
	{
		T_VertexPairList HoleBorder;
		T_VertexPair StartEdge = VertexPairList.front();
		T_VertexPair CurEdge = StartEdge;
		VertexPairList.pop_front();
		HoleBorder.push_back( CurEdge );
		/* try to find the hole border */
		GW_Bool bNextEdgeFound = GW_False;
		while( true )
		{
			bNextEdgeFound = GW_False;
			for( IT_VertexPairList it = VertexPairList.begin(); it!=VertexPairList.end(); ++it )
			{
				T_VertexPair NewEdge = *it;
				if( (NewEdge.first==CurEdge.second) && (NewEdge.second!=CurEdge.first) )
				{
					CurEdge = NewEdge;
					HoleBorder.push_back( CurEdge );
					bNextEdgeFound = GW_True;
					VertexPairList.erase( it );
					break;
				}
				if( (NewEdge.second==CurEdge.second)  && (NewEdge.first!=CurEdge.first) )
				{
					CurEdge = T_VertexPair( NewEdge.second, NewEdge.first);
					HoleBorder.push_back( CurEdge );
					bNextEdgeFound = GW_True;
					VertexPairList.erase( it );
					break;
				}
			}
			if( !bNextEdgeFound )
				break;		// the hole cannot be completed
			if( StartEdge.first == CurEdge.second )
				break;		// the hole is completed
		}
		if( bNextEdgeFound && HoleBorder.size()>2 )	// that means we have a full hole
		{
			char str[50];
			sprintf( str, "Filing a hole of %d vertex.", HoleBorder.size() );
			GW_OutputComment( str );
			IT_VertexPairList it = HoleBorder.begin();
			GW_VoronoiVertex* pVert0 = it->first;		GW_ASSERT( pVert0!=NULL );
			it++;
			GW_VoronoiVertex* pVert1 = it->first;		GW_ASSERT( pVert1!=NULL );
			it++;
			for( ; it!=HoleBorder.end(); ++it )
			{
				GW_VoronoiVertex* pVert2 = it->first;	GW_ASSERT( pVert2!=NULL );
				/* test for manifold structure before creating a new edge [v0,v2] */
				GW_Bool bManifold = GW_True;
				GW_Face* pFace1, *pFace2;
				pVert0->GetFaces( *pVert1, pFace1, pFace2 );
				if( pFace1!=NULL && pFace2!=NULL )
					bManifold = GW_False;
				pVert1->GetFaces( *pVert2, pFace1, pFace2 );
				if( pFace1!=NULL && pFace2!=NULL )
					bManifold = GW_False;
				pVert0->GetFaces( *pVert2, pFace1, pFace2 );
				if( pFace1!=NULL && pFace2!=NULL )
					bManifold = GW_False;
				if( bManifold )
				{
					GW_Face& Face = this->CreateNewFace();
					Face.SetVertex( *pVert0, *pVert1, *pVert2 );
					this->AddFace( Face );
				}				
				pVert1 = pVert2;
			}
		}
	}
}


/*------------------------------------------------------------------------------*/
// Name : GW_VoronoiMesh::CreateNewVertex
/**
*  \return [GW_Vertex&] The newly created vertex.
*  \author Gabriel Peyré
*  \date   4-9-2003
* 
*  Allocate memory for a new vertex. 
*/
/*------------------------------------------------------------------------------*/
GW_Vertex& GW_VoronoiMesh::CreateNewVertex()
{
	return *(new GW_VoronoiVertex);
}

/*------------------------------------------------------------------------------*/
// Name : GW_VoronoiMesh::CreateNewFace
/**
*  \return [GW_Face&] The newly created face.
*  \author Gabriel Peyré
*  \date   4-9-2003
* 
*  Allocate memory for a new face.
*/
/*------------------------------------------------------------------------------*/
GW_Face& GW_VoronoiMesh::CreateNewFace()
{
	return *(new GW_Face);
}

/*------------------------------------------------------------------------------*/
// Name : FastMarchingCallbackFunction_Boundaries
/**
*  \author Gabriel Peyré
*  \date   4-12-2003
* 
*  A callback function for geodesic computations.
*/
/*------------------------------------------------------------------------------*/
GW_Bool GW_VoronoiMesh::FastMarchingCallbackFunction_Boundaries( GW_GeodesicVertex& CurVert )
{
	// \todo for the moment propagate on the full neighborhood
//	for( IT_VoronoiVertexList it=CurrentTargetVertex_.begin(); it!=CurrentTargetVertex_.end(); ++it )
//	{
	GW_ASSERT( pCurrentVoronoiVertex_!=NULL );
	for( GW_VertexIterator it=pCurrentVoronoiVertex_->BeginVertexIterator(); it!=pCurrentVoronoiVertex_->EndVertexIterator(); ++it )
	{
		GW_VoronoiVertex* pVornoiVert = (GW_VoronoiVertex*) *it;
		GW_ASSERT( pVornoiVert!=NULL );
		GW_GeodesicVertex* pVert = pVornoiVert->GetBaseVertex();
		GW_ASSERT( pVert!=NULL );
		if( pVert->GetState()!=GW_GeodesicVertex::kDead )
			return GW_False;
	}
	return GW_True;
}

/*------------------------------------------------------------------------------*/
// Name : GW_VoronoiMesh::BuildGeodesicBoundaries
/**
 *  \author Gabriel Peyré
 *  \date   4-21-2003
 * 
 *  Compute the geodesic path corresponding to each base mesh vertex.
 */
/*------------------------------------------------------------------------------*/
void GW_VoronoiMesh::BuildGeodesicBoundaries( GW_GeodesicMesh& Mesh )
{
	Mesh.RegisterForceStopCallbackFunction( FastMarchingCallbackFunction_Boundaries );

	VertexPathMap_.clear();
	GeodesicDistanceMap_.clear();
	BoundaryEdgeMap_.clear();

	for( GW_U32 i=0; i<this->GetNbrVertex(); ++i )
	{
		CurrentTargetVertex_.clear();
		pCurrentVoronoiVertex_ = (GW_VoronoiVertex*) this->GetVertex(i);
		GW_ASSERT( pCurrentVoronoiVertex_!=NULL );

		/* set the list of vertex to which the boundary needs to be computed */
		for( GW_VertexIterator it=pCurrentVoronoiVertex_->BeginVertexIterator(); it!=pCurrentVoronoiVertex_->EndVertexIterator(); ++it )
		{
			GW_VoronoiVertex* pVert = (GW_VoronoiVertex*) *it;
			GW_ASSERT( pVert!=NULL );
			GW_U32 nID = GW_Vertex::ComputeUniqueId( *pCurrentVoronoiVertex_, *pVert );
			if( VertexPathMap_.find( nID )==VertexPathMap_.end() )
				CurrentTargetVertex_.push_back( pVert );
		}

		/* compute fast marching */
		GW_GeodesicVertex* pGeodesicVert = pCurrentVoronoiVertex_->GetBaseVertex();
		GW_ASSERT( pGeodesicVert!=NULL );
		Mesh.ResetGeodesicMesh();
		GW_OutputComment("Performing Fast Marching.");
		Mesh.PerformFastMarching( pGeodesicVert );

		/* track back and compute geodesic path */
		for( IT_VoronoiVertexList it=CurrentTargetVertex_.begin(); it!=CurrentTargetVertex_.end(); ++it )
		{
			GW_VoronoiVertex* pVert = *it;
			GW_ASSERT( pVert!=NULL );
			GW_GeodesicVertex* pGeodesicVert = pVert->GetBaseVertex();
			GW_ASSERT( pGeodesicVert!=NULL );
			GW_ASSERT( pGeodesicVert->GetState()==GW_GeodesicVertex::kDead );
			GW_U32 nID = GW_Vertex::ComputeUniqueId( *pCurrentVoronoiVertex_, *pVert );
			GW_ASSERT( VertexPathMap_.find( nID ) ==VertexPathMap_.end() );
			T_GeodesicVertexList* pVertexGeodesicEdge = new T_GeodesicVertexList;
			VertexPathMap_[nID] = pVertexGeodesicEdge;
			/* track back and compute the path */
			GW_OutputComment("Computing geodesic path.");
			GW_GeodesicPath GeodesicPath;
			GeodesicPath.ComputePath( *pGeodesicVert, 5000 );
			/* convert the path into vertex */
			this->AddPathToMeshVertex( Mesh, GeodesicPath, *pVertexGeodesicEdge );
		}			
	}

	/* \todo This is just a debug test, remove it */
	for( IT_VertexPathMap it = VertexPathMap_.begin(); it!=VertexPathMap_.end(); ++it )
	{
		T_GeodesicVertexList* pPath = it->second;
		GW_GeodesicVertex* pPrevVert = NULL;
		GW_U32 nNum = 0;
		for( IT_GeodesicVertexList itVert = pPath->begin(); itVert!=pPath->end(); ++itVert )
		{
			GW_GeodesicVertex* pVert = *itVert;
			if( pPrevVert==NULL )
				pPrevVert = pVert;
			else
			{
				GW_Face* pFace1, *pFace2;
				pPrevVert->GetFaces( *pVert, pFace1, pFace2 );
				// GW_ASSERT( pFace1!=NULL || pFace2!=NULL )	// \todo FIX THIS BUG
				if( pFace1!=NULL || pFace2!=NULL )
				{
//					cout << "Gap in a geodesic boundary." << endl;
				}					
				pPrevVert = pVert;
			}
			nNum++;
		}
	}

	Mesh.RegisterForceStopCallbackFunction( NULL );
	Mesh.ResetGeodesicMesh();
}


/*------------------------------------------------------------------------------*/
// Name : GW_VoronoiMesh::AddPathToMeshVertex
/**
 *  \param  CurPath [T_GeodesicVertexList&] The path to test.
 *  \author Gabriel Peyré
 *  \date   5-14-2003
 * 
 *  Test for path intersection and fix cracks.
 */
/*------------------------------------------------------------------------------*/
void GW_VoronoiMesh::AddPathToMeshVertex( GW_GeodesicMesh& Mesh, GW_GeodesicPath& GeodesicPath, T_GeodesicVertexList& VertexPath )
{
	GW_GeodesicVertex* pPrevVert = NULL;
	T_GeodesicPointList& PointList = GeodesicPath.GetPointList();
	for( IT_GeodesicPointList it=PointList.begin(); it!=PointList.end(); ++it )
	{
		GW_GeodesicPoint* pPoint = *it;						GW_ASSERT( pPoint!=NULL );
		GW_GeodesicVertex* pVert1 = pPoint->GetVertex1();	GW_ASSERT( pVert1!=NULL );
		GW_GeodesicVertex* pVert2 = pPoint->GetVertex2();	GW_ASSERT( pVert2!=NULL );
		GW_Float rCoord = pPoint->GetCoord();
		GW_Bool bIsNewVertCreated;
		GW_GeodesicVertex* pNewVert = (GW_GeodesicVertex*) Mesh.InsertVertexInEdge( *pPoint->GetVertex1(), *pPoint->GetVertex2(), pPoint->GetCoord(), bIsNewVertCreated );
		if( bIsNewVertCreated )
		{
			pNewVert->SetDistance( pPoint->GetVertex1()->GetDistance()*rCoord + 
								   pPoint->GetVertex2()->GetDistance()*(1-rCoord) );
			pNewVert->SetState( pPoint->GetVertex1()->GetState() );
			pNewVert->SetFront( pPoint->GetVertex1()->GetFront() );
		}
		/* add the vertex to the path */
		VertexPath.push_back(pNewVert);

		/* test for self intersection */
		IT_GeodesicPointList itSelf = it;
		itSelf++;
		if( bIsNewVertCreated )		// fix only if new vertex is added
		while( itSelf!=PointList.end() )
		{
			GW_GeodesicPoint* pPointSelf = *itSelf;	GW_ASSERT( pPointSelf!=NULL );
			GW_GeodesicVertex* pVertSelf1 = pPointSelf->GetVertex1();	GW_ASSERT( pVertSelf1!=NULL );
			GW_GeodesicVertex* pVertSelf2 = pPointSelf->GetVertex2();	GW_ASSERT( pVertSelf1!=NULL );
			GW_Float rCoordSelf = pPointSelf->GetCoord();
			if( pVertSelf2==pVert1 && pVertSelf1==pVert2 )
			{
				/* we should swap the two vertex */
				GW_GeodesicVertex* pTemp = pVertSelf1;
				pVertSelf1 = pVertSelf2;
				pVertSelf2 = pTemp;
				rCoordSelf = 1-rCoordSelf;
			}
			if( pVertSelf2==pVert2 && pVertSelf1==pVert1 )
			{
				/* there is a self intersection */
				if( rCoordSelf>rCoord )
				{
					pPointSelf->SetVertex1( *pVert1 );
					pPointSelf->SetVertex2( *pNewVert );
					pPointSelf->SetCoord( (rCoordSelf-rCoord)/(1-rCoord) );
				}
				else
				{
					pPointSelf->SetVertex1( *pNewVert );
					pPointSelf->SetVertex2( *pVert2 );
					if( rCoord>GW_EPSILON )
						pPointSelf->SetCoord( rCoordSelf/rCoord );
					else
						pPointSelf->SetCoord(0);
				}
			}
			itSelf++;
		}

		/* now test for intersection with previous path */
		GW_U32 nID = GW_Vertex::ComputeUniqueId( *pVert1, *pVert2 );

		/* to avoid to reprocess already processed lists */
		std::list<T_GeodesicVertexList*> PathToProccess;
		IT_VertexPathMultiMultiMap itIntersectedPath = BoundaryEdgeMap_.find(nID);
		if( bIsNewVertCreated )		// fix only if new vertex is added
		while( itIntersectedPath!=BoundaryEdgeMap_.end() &&  itIntersectedPath->first==nID )
		{
			PathToProccess.push_back( itIntersectedPath->second );
			itIntersectedPath++;
		}
		for( std::list<T_GeodesicVertexList*>::iterator itIntersectedPath=PathToProccess.begin(); itIntersectedPath!=PathToProccess.end(); ++itIntersectedPath )
		{
			/* oups, intersection, we should fix it */
			T_GeodesicVertexList* pIntersectedPath = *itIntersectedPath;
			GW_ASSERT( pIntersectedPath!=NULL );
			/* find where the crossing occurs */
			IT_GeodesicVertexList itIntersec = pIntersectedPath->begin();
			GW_GeodesicVertex* pIntersect1 = NULL;
			GW_GeodesicVertex* pIntersect2 = NULL;
			while( !(itIntersec==pIntersectedPath->end()) )
			{
				pIntersect2 = *itIntersec;
				if( (pIntersect1==pVert1 && pIntersect2==pVert2) ||
					(pIntersect1==pVert2 && pIntersect2==pVert1) )	
					break;		// we have found the position in the path
				pIntersect1 = pIntersect2;
				itIntersec++;
			}
			if( pIntersect1!=NULL && pIntersect2!=NULL && itIntersec!=pIntersectedPath->end() )
			{
				/* insert new vertex in path */
				pIntersectedPath->insert( itIntersec, pNewVert );
				/* divide previous edge in two */
				// BoundaryEdgeMap_.erase( nID );	// \todo find the correct iterator and remove it
				nID = GW_Vertex::ComputeUniqueId( *pVert1, *pNewVert );
				BoundaryEdgeMap_.insert( std::pair<GW_U32, T_GeodesicVertexList*>( nID, pIntersectedPath ) );
				nID = GW_Vertex::ComputeUniqueId( *pVert2, *pNewVert );
				BoundaryEdgeMap_.insert( std::pair<GW_U32, T_GeodesicVertexList*>( nID, pIntersectedPath ) );
			}
			else
				GW_ASSERT( GW_False );
		}
		/* add the small newly created edge in the map */
		if( pPrevVert==NULL )
			pPrevVert = pNewVert;
		else
		{	
			nID = GW_Vertex::ComputeUniqueId( *pPrevVert, *pNewVert );
			pPrevVert = pNewVert;
			BoundaryEdgeMap_.insert( std::pair<GW_U32, T_GeodesicVertexList*>( nID, &VertexPath ) );
		}
	}	
}

/*------------------------------------------------------------------------------*/
// Name : FastMarchingCallbackFunction_Parametrization
/**
*  \author Gabriel Peyré
*  \date   4-12-2003
* 
*  A callback function for geodesic computations.
*/
/*------------------------------------------------------------------------------*/
GW_Bool GW_VoronoiMesh::FastMarchingCallbackFunction_Parametrization( GW_GeodesicVertex& CurVert )
{
	GW_ASSERT( pCurrentVoronoiVertex_!=NULL );

	if( !CurVert.GetIsStoppingVertex() || !CurVert.GetBoundaryReached() )
	{
		/* add the parameter value */
		CurVert.AddParameterVertex( *pCurrentVoronoiVertex_, CurVert.GetDistance() );
		if( CurVert.GetIsStoppingVertex() )
			CurVert.SetBoundaryReached( GW_True );
	}
	/* always continue the computations */
	return GW_False;
}


/*------------------------------------------------------------------------------*/
// Name : GW_VoronoiMesh::PerformLocalFastMarching
/**
 *  \param  Mesh [GW_GeodesicMesh&] The vertex around which we perform fast marching.
 *  \param  Vert [GW_VoronoiVertex&] The base mesh.
 *  \author Gabriel Peyré
 *  \date   4-28-2003
 * 
 *  Perform a fast marching on the voronoi faces around a given 
 *  vertex.
 */
/*------------------------------------------------------------------------------*/
void GW_VoronoiMesh::PerformLocalFastMarching( GW_GeodesicMesh& Mesh, GW_VoronoiVertex& Vert )
{
	Mesh.ResetGeodesicMesh();
	/* set up a boundary to stop the front propagation */
	for( GW_FaceIterator it=Vert.BeginFaceIterator(); it!=Vert.EndFaceIterator(); ++it )
	{
		GW_Face* pFace = *it;
		GW_ASSERT( pFace!=NULL );
		GW_VoronoiVertex* pVert1 = (GW_VoronoiVertex*) pFace->GetNextVertex( Vert );
		GW_VoronoiVertex* pVert2 = (GW_VoronoiVertex*) pFace->GetNextVertex( *pVert1 );
		GW_ASSERT( pVert1!=NULL );
		GW_ASSERT( pVert2!=NULL );
		GW_U32 nID = GW_Vertex::ComputeUniqueId( *pVert1, *pVert2 );
		/* get the geodesic edge */
		GW_ASSERT( VertexPathMap_.find(nID)!=VertexPathMap_.end() );
		T_GeodesicVertexList* pVertexGeodesicEdge = VertexPathMap_[nID];
		GW_ASSERT( pVertexGeodesicEdge!=NULL );
		if( pVertexGeodesicEdge!=NULL )
		for( IT_GeodesicVertexList VertIt = pVertexGeodesicEdge->begin(); VertIt!=pVertexGeodesicEdge->end(); ++VertIt )
		{
			GW_GeodesicVertex* pGeodesicVert = *VertIt;
			pGeodesicVert->SetStoppingVertex( GW_True );
		}
	}
	/* now compute the fast marching */
	GW_ASSERT( Vert.GetBaseVertex()!=NULL );
	Mesh.PerformFastMarching( Vert.GetBaseVertex() );
}


/*------------------------------------------------------------------------------*/
// Name : GW_VoronoiMesh::BuildGeodesicParametrization
/**
 *  \author Gabriel Peyré
 *  \date   4-26-2003
 * 
 *  Compute the parameter of each vertex of the base mesh.
 */
/*------------------------------------------------------------------------------*/
void GW_VoronoiMesh::BuildGeodesicParametrization( GW_GeodesicMesh& Mesh )
{
	Mesh.RegisterForceStopCallbackFunction( FastMarchingCallbackFunction_Parametrization );
	Mesh.ResetParametrizationData();
	/* for each base vertex, perform front propagation on the surrounding faces */
	for( GW_U32 i=0; i<this->GetNbrVertex(); ++i )
	{
		CurrentTargetVertex_.clear();
		pCurrentVoronoiVertex_ = (GW_VoronoiVertex*) this->GetVertex(i);	// set up global data for callback function
		GW_ASSERT( pCurrentVoronoiVertex_!=NULL );

		GW_OutputComment("Performing local fast marching.");
		this->PerformLocalFastMarching( Mesh, *pCurrentVoronoiVertex_ );

		/* register geodesic boundaries length */
		for( GW_VertexIterator it=pCurrentVoronoiVertex_->BeginVertexIterator(); it!=pCurrentVoronoiVertex_->EndVertexIterator(); ++it )
		{
			GW_VoronoiVertex* pNeighbor = (GW_VoronoiVertex*) *it;
			GW_ASSERT( pNeighbor!=NULL );
			GW_GeodesicVertex* pGeodesicVert = pNeighbor->GetBaseVertex();
			GW_ASSERT( pGeodesicVert!=NULL );
			GW_U32 nID = GW_Vertex::ComputeUniqueId( *pCurrentVoronoiVertex_, *pNeighbor );
			if( GeodesicDistanceMap_.find(nID)==GeodesicDistanceMap_.end() )
			{
				GW_ASSERT( pGeodesicVert->GetState()==GW_GeodesicVertex::kDead );
				/* this is the first time the distance is computed */
				if( pGeodesicVert->GetState()==GW_GeodesicVertex::kDead )
					GeodesicDistanceMap_[nID]	 = pGeodesicVert->GetDistance();
			}
			else
			{
				GW_ASSERT( pGeodesicVert->GetState()==GW_GeodesicVertex::kDead );
				if( pGeodesicVert->GetState()==GW_GeodesicVertex::kDead )
				{
					/* distance alreay computed in the reverse direction : take the average */
					GW_Float rPrevDist = GeodesicDistanceMap_[nID];
					GW_Float rNewDist = pGeodesicVert->GetDistance();
					GeodesicDistanceMap_[nID] = GW_MIN(rPrevDist, rNewDist); //(rPrevDist+rNewDist)*0.5;
				}
			}
		}

	}

	/* now really compute the 3 parameters for each vertex */
	GW_OutputComment("Computing geodesic parametrisation.");
	this->ComputeVertexParameters( Mesh );

	pCurrentVoronoiVertex_ = NULL;
	Mesh.RegisterForceStopCallbackFunction( NULL );

}

/*------------------------------------------------------------------------------*/
// Name : GW_VoronoiMesh::ComputeVertexParameters
/**
 *  \param  Mesh [GW_GeodesicMesh&] The mesh.
 *  \author Gabriel Peyré
 *  \date   4-27-2003
 * 
 *  Using geodesic distance from 3 base vertex, compute the 
 *  parameter for each vertex.
 */
/*------------------------------------------------------------------------------*/
void GW_VoronoiMesh::ComputeVertexParameters( GW_GeodesicMesh& Mesh )
{
	for( GW_U32 i=0; i<Mesh.GetNbrVertex(); ++i )
	{
		GW_GeodesicVertex* pVert = (GW_GeodesicVertex*) Mesh.GetVertex( i );
		GW_ASSERT( pVert!=NULL );
		GW_Float d0, d1, d2;
		GW_VoronoiVertex* pVornoiVert0 = pVert->GetParameterVertex( 0, d0 );
		GW_VoronoiVertex* pVornoiVert1 = pVert->GetParameterVertex( 1, d1 );
		GW_VoronoiVertex* pVornoiVert2 = pVert->GetParameterVertex( 2, d2 );
		if( pVornoiVert0!=NULL && pVornoiVert1!=NULL && pVornoiVert2!=NULL )
		{
			GW_ASSERT( pVornoiVert0!=NULL );
			GW_ASSERT( pVornoiVert1!=NULL );
			GW_ASSERT( pVornoiVert2!=NULL );
			GW_GeodesicVertex* pVert0 = pVornoiVert0->GetBaseVertex();
			GW_GeodesicVertex* pVert1 = pVornoiVert1->GetBaseVertex();
			GW_GeodesicVertex* pVert2 = pVornoiVert2->GetBaseVertex();
			GW_ASSERT( pVert0!=NULL );
			GW_ASSERT( pVert1!=NULL );
			GW_ASSERT( pVert2!=NULL );
			GW_U32 nID0 = GW_Vertex::ComputeUniqueId( *pVornoiVert1, *pVornoiVert2 );
			GW_U32 nID1 = GW_Vertex::ComputeUniqueId( *pVornoiVert0, *pVornoiVert2 );
			GW_U32 nID2 = GW_Vertex::ComputeUniqueId( *pVornoiVert0, *pVornoiVert1 );
			if( GeodesicDistanceMap_.find(nID0)!=GeodesicDistanceMap_.end() &&
				GeodesicDistanceMap_.find(nID1)!=GeodesicDistanceMap_.end() &&
				GeodesicDistanceMap_.find(nID2)!=GeodesicDistanceMap_.end() )
			{
				/* length of each triangle segment. */
				GW_Float l0 = GeodesicDistanceMap_[nID0];
				GW_Float l1 = GeodesicDistanceMap_[nID1];
				GW_Float l2 = GeodesicDistanceMap_[nID2];

				/* now use Heron rule to compute parameter */
				GW_Float a0 = GW_Maths::ComputeTriangleArea( l0, d1, d2 );
				GW_Float a1 = GW_Maths::ComputeTriangleArea( d0, l1, d2 );
				GW_Float a2 = GW_Maths::ComputeTriangleArea( d0, d1, l2 );
				GW_Float t = GW_Maths::ComputeTriangleArea( l0, l1, l2 );	// area of the big triangle
				GW_Float a = a0 + a1 + a2;
				GW_ASSERT( a>GW_EPSILON );
				/* test if the point is outside */
				if( a1+a2>t )
					pVert->SetParameterVertex( 0, a1/(a1+a2), a2/(a1+a2) );
				else if( a0+a2>t )
					pVert->SetParameterVertex( a0/(a0+a2), 0, a2/(a0+a2) );
				else if( a0+a1>t )
					pVert->SetParameterVertex( a0/(a0+a1), a1/(a0+a1), 0 );
				else
				{
					/* the point is inside the triangle */
					pVert->SetParameterVertex( a0/a, a1/a, a2/a );
				}
				pVert->SetParameterVertex( a0/a, a1/a, a2/a );
			}
			else
			{
				pVert->SetParameterVertex( 0, 0, 0 );
			}
		}
	}

	/* fix parameter for voronoi vertex */
	for( GW_U32 i=0; i<this->GetNbrVertex(); ++i )
	{
		GW_VoronoiVertex* pVoronoiVert0 = (GW_VoronoiVertex*) this->GetVertex( i );
		GW_ASSERT( pVoronoiVert0!=NULL );
		GW_GeodesicVertex* pVert = pVoronoiVert0->GetBaseVertex();
		GW_ASSERT( pVert!=NULL );
		GW_Face* pFace = pVoronoiVert0->GetFace();
		GW_ASSERT( pFace!=NULL );
		GW_VoronoiVertex* pVoronoiVert1 = (GW_VoronoiVertex*) pFace->GetNextVertex( *pVoronoiVert0 );
		GW_ASSERT( pVoronoiVert1!=NULL );
		GW_VoronoiVertex* pVoronoiVert2 = (GW_VoronoiVertex*) pFace->GetNextVertex( *pVoronoiVert1 );
		GW_ASSERT( pVoronoiVert2!=NULL );
		pVert->ResetParametrizationData();
		pVert->AddParameterVertex( *pVoronoiVert0, 1 );
		pVert->AddParameterVertex( *pVoronoiVert1, 0 );
		pVert->AddParameterVertex( *pVoronoiVert2, 0 );
	}

	/* find the center of each voronoi face */
	CentralParameterMap_.clear();
	T_FloatMap BestVertexValue;
	for( GW_U32 i=0; i<Mesh.GetNbrVertex(); ++i )
	{
		GW_GeodesicVertex* pVert = (GW_GeodesicVertex*) Mesh.GetVertex( i );
		GW_ASSERT( pVert!=NULL );
		GW_Float a,b,c;
		GW_VoronoiVertex* pVert0 = pVert->GetParameterVertex( 0, a );
		GW_ASSERT( pVert0!=NULL );
		GW_VoronoiVertex* pVert1 = pVert->GetParameterVertex( 1, b );
		GW_ASSERT( pVert1!=NULL );
		GW_VoronoiVertex* pVert2 = pVert->GetParameterVertex( 2, c );
		if( pVert2!=NULL )
		{
			GW_Float rBestDist = GW_INFINITE;
			GW_U32 nID = GW_Vertex::ComputeUniqueId( *pVert0, *pVert1, *pVert2 );
			if( BestVertexValue.find(nID)!=BestVertexValue.end() )
				rBestDist = BestVertexValue[nID];
			GW_Float rDist = GW_ABS(a - 1.0/3.0) + GW_ABS(b - 1.0/3.0) + GW_ABS(c - 1.0/3.0);
			if( rDist<rBestDist )
			{
				BestVertexValue[nID] = rDist;
				CentralParameterMap_[nID] = pVert;
			}
		}
	}
}

/*------------------------------------------------------------------------------*/
// Name : GW_VoronoiMesh::ResetOnlyVertexState
/**
 *  \param  Mesh [GW_GeodesicMesh&] The mesh.
 *  \author Gabriel Peyré
 *  \date   5-13-2003
 * 
 *  Resst only the state of each vertex of the mesh.
 */
/*------------------------------------------------------------------------------*/
void GW_VoronoiMesh::ResetOnlyVertexState( GW_GeodesicMesh& Mesh )
{
	for( GW_U32 i = 0; i<Mesh.GetNbrVertex(); ++i  )
	{
		GW_GeodesicVertex* pVert = (GW_GeodesicVertex*) Mesh.GetVertex(i);
		GW_ASSERT( pVert!=NULL );
		pVert->SetState( GW_GeodesicVertex::kFar );
	}
}



void GW_VoronoiMesh::PrepareInterpolation( GW_GeodesicMesh& Mesh )
{
	/* duplicate the geodesic data for each of the vertex */
	for( GW_U32 i=0; i<Mesh.GetNbrVertex(); ++i )
	{
		GW_GeodesicVertex* pVert = (GW_GeodesicVertex*) Mesh.GetVertex(i);
		GW_ASSERT( pVert!=NULL );
		GW_GeodesicInformationDuplicata* pDuplicata = new GW_GeodesicInformationDuplicata( *pVert );
	}
}


/*------------------------------------------------------------------------------*/
// Name : GW_VoronoiMesh::FastMarchingCallbackFunction_VertexInsersionNN
/**
*  \param  CurVert [GW_GeodesicVertex&] The current vertex.
*  \param  rNewDist [GW_Float] The new distance.
*  \return [GW_Bool] Insert this vertex to active list ?
*  \author Gabriel Peyré
*  \date   5-13-2003
* 
*  Add the vertex only if the distance decreases.
*/
/*------------------------------------------------------------------------------*/
GW_Bool GW_VoronoiMesh::FastMarchingCallbackFunction_VertexInsersionNN( GW_GeodesicVertex& CurVert, GW_Float rNewDist )
{
	GW_GeodesicInformationDuplicata* pDuplicata = (GW_GeodesicInformationDuplicata*) CurVert.GetUserData();
	GW_ASSERT( pDuplicata!=NULL );
	return pDuplicata->rDistance_+GW_EPSILON >= rNewDist;
}


/*------------------------------------------------------------------------------*/
// Name : GW_VoronoiMesh::FlattenBasePoints
/**
 *  \param  Mesh [GW_GeodesicMesh&] The original mesh.
 *  \param  FlatteningMap [T_Vector2DMap&] A link between the basis geodesic vertex and their 2D positions.
 *  \author Gabriel Peyré
 *  \date   6-1-2003
 * 
 *  Compute a 2D position for each base point, using multidimensional scaling.
 */
/*------------------------------------------------------------------------------*/
void GW_VoronoiMesh::FlattenBasePoints( GW_GeodesicMesh& Mesh, T_Vector2DMap& FlatteningMap )
{
	char str[100];
	/* create a voronoi vertex for each basis vertex */
	this->CreateVoronoiVertex();

	GW_U32 n = this->GetNbrVertex();
	GW_MatrixNxP dist_matrix(n,n,-1);


	/* compute the distance matrix */
	for( GW_U32 i=0; i<n; ++i )
	{
		GW_VoronoiVertex* pVoronoiVerti = (GW_VoronoiVertex*) this->GetVertex(i);					GW_ASSERT( pVoronoiVerti!=NULL );
		GW_GeodesicVertex* pGeodesicVerti = (GW_GeodesicVertex*) pVoronoiVerti->GetBaseVertex();	GW_ASSERT( pGeodesicVerti!=NULL );
		
		sprintf( str, "Computing distance information for vertex %d.", i );
		GW_OutputComment( str );

		/* compute the distance map */
		Mesh.ResetGeodesicMesh();
		Mesh.PerformFastMarching( pGeodesicVerti );

		for( GW_U32 j=0; j<n; ++j )
		{
			GW_VoronoiVertex* pVoronoiVertj = (GW_VoronoiVertex*) this->GetVertex(j);					GW_ASSERT( pVoronoiVerti!=NULL );
			GW_GeodesicVertex* pGeodesicVertj = (GW_GeodesicVertex*) pVoronoiVertj->GetBaseVertex();	GW_ASSERT( pGeodesicVerti!=NULL );
			GW_Float rDist = pGeodesicVertj->GetDistance();
			if( dist_matrix[i][j]<0 )
				dist_matrix[i][j] = rDist*rDist;
			else
				dist_matrix[i][j] = (dist_matrix[i][j]+rDist*rDist)*0.5;

			if( dist_matrix[j][i]<0 )
				dist_matrix[j][i] = rDist*rDist;
			else
				dist_matrix[j][i] = (dist_matrix[j][i]+rDist*rDist)*0.5;
		}
	}

	/* center the matrix */
	sprintf( str, "Performing the multidimensional scaling.", i );
	GW_OutputComment( str );
	GW_MatrixNxP J(n,n,-1.0/n);
	for( GW_U32 i=0; i<n; ++i )
		J[i][i] += 1;

	/* center the matrix */
	dist_matrix = J*dist_matrix*J*(-0.5);

	/* perform eigen-decomposition */
	GW_MatrixNxP v(n,n);
	GW_VectorND RealEig(n);
	dist_matrix.Eigenvalue( v, NULL, &RealEig, NULL );

	GW_Float rScaleX = 0;
	GW_Float rScaleY = 0;
	GW_I32 nEigenvectorX = -1;
	GW_I32 nEigenvectorY = -1;
	/* find the two maximum eigenvalues */
	for( GW_U32 i=0; i<n; ++i )
	{
		if( RealEig[i]>=rScaleX )
		{
			rScaleY = rScaleX;
			nEigenvectorY = nEigenvectorX;
			rScaleX = RealEig[i];
			nEigenvectorX = i;
		}
		else if( RealEig[i]>=rScaleY )
		{
			rScaleY = RealEig[i]; 
			nEigenvectorY = i;
		}
	}
	rScaleX = sqrt(rScaleX);
	rScaleY = sqrt(rScaleY);
	GW_ASSERT( nEigenvectorX>=0 );
	GW_ASSERT( nEigenvectorY>=0 );

	/* compute the 2D position */
	FlatteningMap.clear();
	if( nEigenvectorX>=0 && nEigenvectorY>=0 )
	for( GW_U32 i=0; i<n; ++i )
	{
		GW_VoronoiVertex* pVoronoiVerti = (GW_VoronoiVertex*) this->GetVertex(i);					GW_ASSERT( pVoronoiVerti!=NULL );
		GW_GeodesicVertex* pGeodesicVerti = (GW_GeodesicVertex*) pVoronoiVerti->GetBaseVertex();	GW_ASSERT( pGeodesicVerti!=NULL );

		GW_Vector2D pos;
		pos[0] = rScaleX*v[i][nEigenvectorX];
		pos[1] = rScaleY*v[i][nEigenvectorY];
		FlatteningMap[ pGeodesicVerti->GetID() ] = pos;
	}
}


/*------------------------------------------------------------------------------*/
// Name : GW_VoronoiMesh::FastMarchingCallbackFunction_VertexInsersionRD1
/**
*  \param  CurVert [GW_GeodesicVertex&] The current vertex.
*  \param  rNewDist [GW_Float] The new distance.
*  \return [GW_Bool] Insert this vertex to active list ?
*  \author Gabriel Peyré
*  \date   5-13-2003
* 
*  Here we only record the natural neighbors in the weights map.
*/
/*------------------------------------------------------------------------------*/
GW_Bool GW_VoronoiMesh::FastMarchingCallbackFunction_VertexInsersionRD1( GW_GeodesicVertex& CurVert, GW_Float rNewDist )
{
	GW_ASSERT( pCurWeights_!=NULL );
	GW_GeodesicInformationDuplicata* pDuplicata = (GW_GeodesicInformationDuplicata*) CurVert.GetUserData();
	GW_ASSERT( pDuplicata!=NULL );
	GW_U32 nId = pDuplicata->pFront_->GetID();
	(*pCurWeights_)[ nId ] = -1;	// this is one of our natural neighbor, yeah.
	return pDuplicata->rDistance_+GW_EPSILON >= rNewDist;
}


/*------------------------------------------------------------------------------*/
// Name : GW_VoronoiMesh::FastMarchingCallbackFunction_VertexInsersionRD2
/**
*  \param  CurVert [GW_GeodesicVertex&] The current vertex.
*  \param  rNewDist [GW_Float] The new distance.
*  \return [GW_Bool] Insert this vertex to active list ?
*  \author Gabriel Peyré
*  \date   5-13-2003
* 
*  Here we record the reciprocical of the distance.
*/
/*------------------------------------------------------------------------------*/
GW_Bool GW_VoronoiMesh::FastMarchingCallbackFunction_VertexInsersionRD2( GW_GeodesicVertex& CurVert, GW_Float rNewDist )
{
	GW_ASSERT( pCurWeights_!=NULL );
	GW_U32 nId = CurVert.GetID();
	if( pCurWeights_->find(nId)!=pCurWeights_->end() )	// this is one of our natural neighbor, yeah.
	{
		if( (*pCurWeights_)[nId]<0 )
		{
			/* first time we encounter this vertex, sounds good */
			if( rNewDist>0 )
				(*pCurWeights_)[nId] = 1.0/(rNewDist);
			else
				(*pCurWeights_)[nId] = GW_INFINITE;
			nNbrBaseVertex_RD_++;
		}
		else
			GW_ASSERT( GW_False );
	}
	return GW_True;
}


GW_Bool GW_VoronoiMesh::FastMarchingCallbackFunction_ForceStopRD( GW_GeodesicVertex& Vert )
{
	return nNbrBaseVertex_RD_>=pCurWeights_->size();
}





///////////////////////////////////////////////////////////////////////////////
//  Copyright (c) Gabriel Peyré
///////////////////////////////////////////////////////////////////////////////
//                               END OF FILE                                 //
///////////////////////////////////////////////////////////////////////////////
