text: #include "stdafx.h"
  #include "GW_Mesh.h"
  #include "GW_VertexIterator.h"
  
  #ifndef GW_USE_INLINE
      #include "GW_Mesh.inl"
  #endif
  
  using namespace GW;
  
  GW_Mesh* GW_Mesh::pStaticThis_ = NULL;
  
  
  /*------------------------------------------------------------------------------*/
  // Name : GW_Mesh::operator=
  /**
   *  \param  v [GW_Mesh&] Mesh.
   *  \return [GW_Mesh&] *this
   *
   *  Copy operator.
   */
  /*------------------------------------------------------------------------------*/
  GW_Mesh& GW_Mesh::operator=(const GW_Mesh& Mesh)
  {
      this->SetNbrVertex( Mesh.GetNbrVertex() );
      this->SetNbrFace( Mesh.GetNbrFace() );
  
      for( GW_U32 i=0; i<this->GetNbrVertex(); ++i )
      {
          if( this->GetVertex(i)==NULL )
          {
              GW_Vertex& NewVert = this->CreateNewVertex();
              this->SetVertex(i, &NewVert);
          }
          GW_Vertex& NewVert = *this->GetVertex(i);
          const GW_Vertex& OriginalVert = *Mesh.GetVertex(i);
          NewVert = OriginalVert;
          /* resolve face attachement */
          const GW_Face* pFace = OriginalVert.GetFace();
          if( pFace!=NULL )
              NewVert.SetFace( *this->GetFace( pFace->GetID() ) );
      }
      for( GW_U32 i=0; i<this->GetNbrFace(); ++i )
      {
          if( this->GetFace(i)==NULL )
          {
              GW_Face& NewFace = this->CreateNewFace();
              this->SetFace(i, &NewFace);
          }
          GW_Face& NewFace = *this->GetFace(i);
          const GW_Face& OriginalFace = *Mesh.GetFace(i);
          NewFace = OriginalFace;
          /* resolve Vertex and neighbor */
          for( GW_U32 j=0; j<3; ++j )
          {
              GW_U32 VertID = OriginalFace.GetVertex(j)->GetID();
              NewFace.SetVertex( *this->GetVertex(VertID), j );
              const GW_Face* pNeigh = OriginalFace.GetFaceNeighbor(j);
              if( pNeigh==NULL )
                  NewFace.SetFaceNeighbor(NULL, j);
              else
                  NewFace.SetFaceNeighbor( this->GetFace( pNeigh->GetID() ), j);
          }
      }
  
      return *this;
  }
  
  
  /*------------------------------------------------------------------------------*/
  // Name : GW_Mesh::SetNbrFace
  /**
   *  \param  nNum [GW_U32] New number of faces.
   *  Resize the mesh.
   */
  /*------------------------------------------------------------------------------*/
  void GW_Mesh::SetNbrFace( GW_U32 nNum )
  {
      GW_U32 nOldSize = this->GetNbrFace();
      if( nNum<nOldSize )
      {
          /* check if the vertex at the end should be deleted */
          for( GW_U32 i=nNum; i<nOldSize; ++i )
              GW_SmartCounter::CheckAndDelete( this->GetFace( i ) );
          FaceVector_.resize( nNum );
      }
      if( nNum>nOldSize )
      {
          FaceVector_.resize( nNum );
          /* set to NULL newly appended pointers */
          for( GW_U32 i=nOldSize; i<nNum; ++i )
              this->SetFace( i, NULL );
      }
  }
  
  /*------------------------------------------------------------------------------*/
  // Name : GW_Mesh::SetNbrVertex
  /**
   *  \param  nNum [GW_U32] New number of vertex.
   *
   *  Resize the mesh.
   */
  /*------------------------------------------------------------------------------*/
  void GW_Mesh::SetNbrVertex( GW_U32 nNum )
  {
      GW_U32 nOldSize = this->GetNbrVertex();
      if( nNum<nOldSize )
      {
          /* check if the vertex at the end should be deleted */
          for( GW_U32 i=nNum; i<nOldSize; ++i )
              GW_SmartCounter::CheckAndDelete( this->GetVertex( i ) );
          //!!VertexVector_.resize( nNum );
          VertexVector_size=nNum;
      }
      if( nNum>nOldSize )
      {
          //!!VertexVector_.resize( nNum );
          GW_Vertex** nv = new GW_Vertex*[VertexVector_size=nNum];            //make a larger vector for vertices
          for(GW_U32 i=0;i<nOldSize;++i) nv[i] = VertexVector_[i];                //copy the old vertices (with whatever use counter they have)
          //if(nOldSize>0)
          delete[] VertexVector_;
          VertexVector_ = nv;
  
          //for( GW_U32 i=0; i<nNum; ++i )
          //    VertexVector_[i] = new GW_Vertex(); //LUK:
  
          /* set to NULL newly appended pointers */
          for( GW_U32 i=nOldSize; i<nNum; ++i )
              VertexVector_[i] = NULL;
              //this->SetVertex( i, NULL );
      }
  }
  
  
  
  
  /*------------------------------------------------------------------------------*/
  // Name : GW_Mesh::GetBoundingRadius
  /**
   *  \return [GW_Float] Radius.
   *
   *  Get the maximum distance from (0,0,0)
   */
  /*------------------------------------------------------------------------------*/
  GW_Float GW_Mesh::GetBoundingRadius()
  {
      GW_Float rRadius = 0;
      for( GW_U32 i=0; i<this->GetNbrVertex(); ++i )
      {
          if( this->GetVertex(i)!=NULL )
              rRadius = GW_MAX(this->GetVertex(i)->GetPosition().SquareNorm(), rRadius);
      }
  
      return (GW_Float) sqrt(rRadius);
  }
  
  
  void GW_Mesh::GetBoundingBox( GW_Vector3D& min, GW_Vector3D& max )
  {
      min.SetValue( GW_INFINITE );
      max.SetValue( -GW_INFINITE );
      for( GW_U32 i=0; i<this->GetNbrVertex(); ++i )
      {
          if( this->GetVertex(i)!=NULL )
          {
              GW_Vector3D& pos = this->GetVertex(i)->GetPosition();
              min[0] = GW_MIN( min[0], pos[0] );
              min[1] = GW_MIN( min[1], pos[1] );
              min[2] = GW_MIN( min[2], pos[2] );
              max[0] = GW_MAX( max[0], pos[0] );
              max[1] = GW_MAX( max[1], pos[1] );
              max[2] = GW_MAX( max[2], pos[2] );
          }
      }
  }
  
  /*------------------------------------------------------------------------------*/
  // Name : GW_Mesh::GetBarycenter
  /**
  *  \return [GW_Vector3D] The iso-barycenter of the mesh.
  *
  *  Get the maximum distance from (0,0,0).
  */
  /*------------------------------------------------------------------------------*/
  GW_Vector3D GW_Mesh::GetBarycenter()
  {
      GW_Vector3D Barycenter;
      for( GW_U32 i=0; i<this->GetNbrVertex(); ++i )
      {
          if( this->GetVertex(i)!=NULL )
              Barycenter += this->GetVertex(i)->GetPosition();
      }
  
      if( this->GetNbrVertex()>0 )
          Barycenter /= (GW_Float) this->GetNbrVertex();
      return Barycenter;
  }
  
  /*------------------------------------------------------------------------------*/
  // Name : GW_Mesh::ScaleVertex
  /**
   *  \param  rScale [GW_Float] Scaling factor
   *
   *  Scale each vertex.
   */
  /*------------------------------------------------------------------------------*/
  void GW_Mesh::ScaleVertex( GW_Float rScale )
  {
      for( GW_U32 i=0; i<this->GetNbrVertex(); ++i )
      {
          if( this->GetVertex(i)!=NULL )
              this->GetVertex(i)->GetPosition() *= rScale;
      }
  }
  
  /*------------------------------------------------------------------------------*/
  // Name : GW_Mesh::TranslateVertex
  /**
  *  \param  Vect [GW_Float] Translation vector.
  *
  *  Scale each vertex.
  */
  /*------------------------------------------------------------------------------*/
  void GW_Mesh::TranslateVertex( const GW_Vector3D& Vect )
  {
      for( GW_U32 i=0; i<this->GetNbrVertex(); ++i )
      {
          if( this->GetVertex(i)!=NULL )
              this->GetVertex(i)->GetPosition() += Vect;
      }
  }
  
  /*------------------------------------------------------------------------------*/
  // Name : GW_Mesh::BuildNormal
  /**
   *  Compute vertex normals form faces.
   */
  /*------------------------------------------------------------------------------*/
  void GW_Mesh::BuildRawNormal()
  {
      //!!for( IT_VertexVector it=VertexVector_.begin(); it!=VertexVector_.end(); ++it )
      for(GW_Vertex** it = VertexVector_; it-VertexVector_<VertexVector_size; ++it)
      {
          GW_Vertex* pVert = *it;
          GW_ASSERT( pVert!=NULL );
          pVert->BuildRawNormal();
      }
  }
  
  /*------------------------------------------------------------------------------*/
  // Name : GW_Mesh::BuildCurvatureData
  /**
   *
   *  Compute all curvature data using comlex schemes. This includes
   *  normal, curvatures and curvature directions.
   */
  /*------------------------------------------------------------------------------*/
  void GW_Mesh::BuildCurvatureData()
  {
      //!!for( IT_VertexVector it=VertexVector_.begin(); it!=VertexVector_.end(); ++it )
      for(GW_Vertex** it = VertexVector_; it-VertexVector_<VertexVector_size; ++it)
      {
          GW_Vertex* pVert = *it;
          GW_ASSERT( pVert!=NULL );
          pVert->BuildCurvatureData();
      }
  }
  
  
  /*------------------------------------------------------------------------------*/
  // Name : GW_Mesh::BuildConnectivity
  /**
  *  Call this method when you have set the vertex and the face.
  *    This will set up the neighboorhood for each face.
  */
  /*------------------------------------------------------------------------------*/
  void GW_Mesh::BuildConnectivity()
  {
      T_FaceList* VertexToFaceMap = new T_FaceList[this->GetNbrVertex()];
  
      // build the inverse map vertex->face
      for(IT_FaceVector it = FaceVector_.begin(); it!=FaceVector_.end(); ++it )
      {
          GW_Face* pFace = *it;
          GW_ASSERT( pFace!=NULL );
          for( GW_U32 i=0; i<3; ++i )
          {
              GW_Vertex* pVert = pFace->GetVertex(i);
              GW_ASSERT(pVert!=NULL);
              GW_ASSERT( pVert->GetID() <= this->GetNbrVertex() );
              VertexToFaceMap[pVert->GetID()].push_back( pFace );
          }
      }
  
      // now we can set up connectivity
      for( IT_FaceVector it=FaceVector_.begin(); it!=FaceVector_.end(); ++it )
      {
          GW_Face* pFace = *it;
          GW_ASSERT( pFace!=NULL );
  
          /* set up the neigbooring faces of the 3 vertices */
          T_FaceList* pFaceLists[3];
          for( GW_U32 i=0; i<3; ++i )
          {
              GW_Vertex* pVert = pFace->GetVertex(i);
              pFaceLists[i] = &VertexToFaceMap[pVert->GetID()];
          }
  
          /* compute neighbor in the 3 directions */
          for( GW_U32 i=0; i<3; ++i )
          {
              GW_Face* pNeighbor = NULL;
              GW_U32 i1 = (i+1)%3;
              GW_U32 i2 = (i+2)%3;
              /* we must find the intersection of the surrounding faces of these 2 vertex */
              GW_Bool bFind = GW_False;
              for( IT_FaceList it1 = pFaceLists[i1]->begin(); it1!=pFaceLists[i1]->end() && bFind!=GW_True; ++it1 )
              {
                  GW_Face* pFace1 = *it1;
                  for( IT_FaceList it2 = pFaceLists[i2]->begin(); it2!=pFaceLists[i2]->end() && bFind!=GW_True; ++it2 )
                  {
                      GW_Face* pFace2 = *it2;
                      if( pFace1==pFace2 && pFace1!=pFace )
                      {
                          pNeighbor = pFace1;
                          bFind=GW_True;
                      }
                  }
              }
              pFace->SetFaceNeighbor( pNeighbor, i );
              /* make some test on the neighbor to assure symetry in the connectivity relationship */
              if( pNeighbor!=NULL )
              {
                  GW_I32 nEdgeNumber = pNeighbor->GetEdgeNumber( *pFace->GetVertex(i1),*pFace->GetVertex(i2) );
                  GW_ASSERT( nEdgeNumber>=0 );
                  pNeighbor->SetFaceNeighbor( pFace, nEdgeNumber );
              }
          }
      }
  
      GW_DELETEARRAY( VertexToFaceMap );
  }
  
  
  /*------------------------------------------------------------------------------*/
  // Name : GW_Mesh::GetArea
  /**
   *  \return [GW_Float] Total area.
   *
   *  Compute the total area of the triangulated mesh.
   */
  /*------------------------------------------------------------------------------*/
  GW_Float GW_Mesh::GetArea()
  {
      GW_Float rArea = 0;
  
      for( IT_FaceVector it=FaceVector_.begin(); it!=FaceVector_.end(); ++it)
      {
          GW_Face* pFace = *it;
          GW_ASSERT( pFace!=NULL );
          GW_Vertex* v0 = pFace->GetVertex(0);
          GW_Vertex* v1 = pFace->GetVertex(1);
          GW_Vertex* v2 = pFace->GetVertex(2);
          if( v0!=NULL && v1!=NULL && v2!=NULL )
          {
              GW_Vector3D e1 = v1->GetPosition() - v0->GetPosition();
              GW_Vector3D e2 = v2->GetPosition() - v0->GetPosition();
              rArea += ~(e1 ^ e2);
          }
      }
  
      return (GW_Float) 0.5*rArea;
  }
  
  /*------------------------------------------------------------------------------*/
  // Name : GW_Mesh::GetRandomVertex
  /**
   *  \return [GW_Vertex*] The vertex. NULL if it was impossible.
   *
   *  Get Return a vertex at random.
   */
  /*------------------------------------------------------------------------------*/
  GW_Vertex* GW_Mesh::GetRandomVertex()
  {
      GW_U32 nNumber = 0;
      GW_Vertex* pStartVertex = NULL;
      while( pStartVertex==NULL )
      {
          if( nNumber>=this->GetNbrVertex()/10 )
              return NULL;
          GW_U32 nNumVert = (GW_U32) floor(GW_RAND*this->GetNbrVertex());
          pStartVertex = this->GetVertex( nNumVert );
          if( pStartVertex->GetFace()==NULL )
              pStartVertex = NULL;
          nNumber++;
      }
      return pStartVertex;
  }
  
  /*------------------------------------------------------------------------------*/
  // Name : GW_Mesh::InsertVertexInFace
  /**
   *  \param  Face [GW_Face&] The face.
   *  \param  x [GW_Float] x barycentric coordinate.
   *  \param  y [GW_Float] y barycentric coordinate.
   *  \param  z [GW_Float] z barycentric coordinate.
   *    \return The newly created point.
   *
   *  Insert a vertex in a face.
   */
  /*------------------------------------------------------------------------------*/
  GW_Vertex* GW_Mesh::InsertVertexInFace( GW_Face& Face, GW_Float x, GW_Float y, GW_Float z )
  {
      GW_Vertex* pVert0 = Face.GetVertex(0);
      GW_Vertex* pVert1 = Face.GetVertex(1);
      GW_Vertex* pVert2 = Face.GetVertex(2);
      GW_ASSERT( pVert0!=NULL );
      GW_ASSERT( pVert1!=NULL );
      GW_ASSERT( pVert2!=NULL );
      /* create two new-faces */
      GW_Face* pFace1 = &this->CreateNewFace();
      GW_Face* pFace2 = &this->CreateNewFace();
      this->SetNbrFace( this->GetNbrFace()+2 );
      this->SetFace( this->GetNbrFace()-2, pFace1 );
      this->SetFace( this->GetNbrFace()-1, pFace2 );
      /* create one new vertex */
      GW_Vertex* pNewVert = &this->CreateNewVertex();
      GW_Vector3D pos = pVert0->GetPosition()*x + pVert1->GetPosition()*y + pVert2->GetPosition()*z;
      pNewVert->SetPosition( pos );
      pNewVert->BuildRawNormal();
      this->SetNbrVertex( this->GetNbrVertex()+1 );
      this->SetVertex( this->GetNbrVertex()-1, pNewVert );
      /* assign vertex to faces */
      pFace1->SetVertex( *pVert0, *pVert1, *pNewVert );
      pFace2->SetVertex( *pNewVert, *pVert1, *pVert2 );
      Face.SetVertex( *pVert0, *pNewVert, *pVert2 );
      /* assign dependance vertex->mother face */
      pNewVert->SetFace( Face );
      pVert0->SetFace( Face );
      pVert1->SetFace( *pFace1 );
      pVert2->SetFace( Face );
      /* outer faces */
      if( Face.GetFaceNeighbor(2)!=NULL )
      {
          GW_I32 nEdgeNumber = Face.GetFaceNeighbor(2)->GetEdgeNumber( Face );
          GW_ASSERT( nEdgeNumber>=0 );
          Face.GetFaceNeighbor(2)->SetFaceNeighbor( pFace1, nEdgeNumber );
      }
      if( Face.GetFaceNeighbor(0)!=NULL )
      {
          GW_I32 nEdgeNumber = Face.GetFaceNeighbor(0)->GetEdgeNumber( Face );
          GW_ASSERT( nEdgeNumber>=0 );
          Face.GetFaceNeighbor(0)->SetFaceNeighbor( pFace2, nEdgeNumber );
      }
      /* build connectivity of inner faces */
      pFace1->SetFaceNeighbor( pFace2, &Face, Face.GetFaceNeighbor(2) );
      pFace2->SetFaceNeighbor( Face.GetFaceNeighbor(0), &Face,  pFace1 );
      GW_Face* pTempFace = Face.GetFaceNeighbor(1);
      Face.SetFaceNeighbor(  pFace2, pTempFace, pFace1 );
  
      return pNewVert;
  }
  
  
  /*------------------------------------------------------------------------------*/
  // Name : GW_Mesh::InsertVertexInEdge
  /**
   *  \param  Vert1 [GW_Vertex&] 1st vertex.
   *  \param  Vert2 [GW_Vertex&] 2nd vertex.
   *  \param  x [GW_Float] Barycentric coords in [0,1].
   *  \return [GW_Vertex&] The vertex. This is a newly created one only if 0<<x<<1.
   *
   *  Insert a vertex in the middle of an edge.
   */
  /*------------------------------------------------------------------------------*/
  GW_Vertex* GW_Mesh::InsertVertexInEdge( GW_Vertex& Vert1, GW_Vertex& Vert2, GW_Float x, GW_Bool& bIsNewVertCreated )
  {
      if( x<GW_EPSILON )
      {
          bIsNewVertCreated = GW_False;
          return &Vert2;
      }
      if( x>1-GW_EPSILON )
      {
          bIsNewVertCreated = GW_False;
          return &Vert1;
      }
      bIsNewVertCreated = GW_True;
      /* create the new vertex */
      GW_Vertex* pNewVert = &this->CreateNewVertex();
      this->SetNbrVertex( this->GetNbrVertex()+1 );
      this->SetVertex( this->GetNbrVertex()-1, pNewVert );
      /* set position */
      pNewVert->SetPosition( Vert1.GetPosition()*x + Vert2.GetPosition()*(1-x) );
      /* retrieve the neigbor faces face */
      GW_Face* pFace1 = NULL;
      GW_Face* pFace2 = NULL;
      Vert1.GetFaces( Vert2, pFace1, pFace2 );
      GW_ASSERT( pFace1!=NULL || pFace2!=NULL );
      /* assign the face of new vertex */
      if( pFace1!=NULL )
          pNewVert->SetFace( *pFace1 );
      else if( pFace2!=NULL )
          pNewVert->SetFace( *pFace2 );
      GW_I32 nVert1Num1, nVert1Num2, nVert1Num3, nVert2Num1, nVert2Num2, nVert2Num3;
      GW_Face* pNewFace1 = NULL;
      GW_Face* pNewFace2 = NULL;
      if( pFace1!=NULL )
      {
          nVert1Num1 = pFace1->GetEdgeNumber( Vert1 );
          GW_ASSERT( nVert1Num1>=0 );
          nVert1Num2 = pFace1->GetEdgeNumber( Vert2 );
          GW_ASSERT( nVert1Num2>=0 );
          nVert1Num3 = 3-nVert1Num2-nVert1Num1;
          pNewFace1 = &this->CreateNewFace();
          this->SetNbrFace( this->GetNbrFace()+1 );
          this->SetFace( this->GetNbrFace()-1, pNewFace1 );
          pNewFace1->SetVertex( *pFace1->GetVertex(nVert1Num3), nVert1Num3 );
          pNewFace1->SetVertex( Vert2, nVert1Num2 );
          pNewFace1->SetVertex( *pNewVert, nVert1Num1 );
          /* connectivity between Face1 and new face */
          GW_Face* pFaceNeighbor = pFace1->GetFaceNeighbor(nVert1Num1);
          pNewFace1->SetFaceNeighbor( pFaceNeighbor, nVert1Num1 );
          pNewFace1->SetFaceNeighbor( pFace1, nVert1Num2 );
          if( pFaceNeighbor!=NULL )
          {
              GW_I32 nNum = pFaceNeighbor->GetEdgeNumber( Vert2, *pFace1->GetVertex(nVert1Num3) );
              GW_ASSERT( nNum>=0 );
              pFaceNeighbor->SetFaceNeighbor( pNewFace1, nNum );
          }
          /* connectivity for face 1 */
          pFace1->SetFaceNeighbor( pNewFace1, nVert1Num1 );
          pFace1->SetVertex( *pNewVert, nVert1Num2 );
          /* reassign vertex 2 */
          Vert2.SetFace( *pNewFace1 );
      }
      if( pFace2!=NULL )
      {
          nVert2Num1 = pFace2->GetEdgeNumber( Vert1 );
          GW_ASSERT( nVert2Num1>=0 );
          nVert2Num2 = pFace2->GetEdgeNumber( Vert2 );
          GW_ASSERT( nVert2Num2>=0 );
          nVert2Num3 = 3-nVert2Num2-nVert2Num1;
          pNewFace2 = &this->CreateNewFace();
          this->SetNbrFace( this->GetNbrFace()+1 );
          this->SetFace( this->GetNbrFace()-1, pNewFace2 );
          pNewFace2->SetVertex( *pFace2->GetVertex(nVert2Num3), nVert2Num3 );
          pNewFace2->SetVertex( Vert2, nVert2Num2 );
          pNewFace2->SetVertex( *pNewVert, nVert2Num1 );
          /* connectivity between Face2 and new face */
          GW_Face* pFaceNeighbor = pFace2->GetFaceNeighbor(nVert2Num1);
          pNewFace2->SetFaceNeighbor( pFaceNeighbor, nVert2Num1 );
          pNewFace2->SetFaceNeighbor( pFace2, nVert2Num2 );
          if( pFaceNeighbor!=NULL )
          {
              GW_I32 nNum = pFaceNeighbor->GetEdgeNumber( Vert2, *pFace2->GetVertex(nVert2Num3) );
              GW_ASSERT( nNum>=0 );
              pFaceNeighbor->SetFaceNeighbor( pNewFace2, nNum );
          }
          /* connectivity for face 1 */
          pFace2->SetFaceNeighbor( pNewFace2, nVert2Num1 );
          pFace2->SetVertex( *pNewVert, nVert2Num2 );
          /* reassign vertex 2 */
          Vert2.SetFace( *pNewFace2 );
      }
      /* set inter connectivity */
      if( pNewFace1!=NULL )
          pNewFace1->SetFaceNeighbor( pNewFace2, nVert1Num3 );
      if( pNewFace2!=NULL )
          pNewFace2->SetFaceNeighbor( pNewFace1, nVert2Num3 );
  
      /* last thing : compute normal */
      pNewVert->BuildRawNormal();
  
      return pNewVert;
  }
  
  /*------------------------------------------------------------------------------*/
  // Name : GW_Mesh::IterateConnectedComponent_Vertex
  /**
  *  \param  pCallback [VertexIterate_Callback] The function.
  *
  *  Iterate a callback on each vertex of a connected component.
  */
  /*------------------------------------------------------------------------------*/
  void GW_Mesh::IterateConnectedComponent_Vertex( GW_Vertex& start_vert, VertexIterate_Callback pCallback )
  {
      /* march on the voronoi diagram */
      T_VertexList VertexToProceed;
      VertexToProceed.push_back( &start_vert );
      T_VertexMap VertexDone;
      VertexDone[ start_vert.GetID() ] = &start_vert;
  
  
      while( !VertexToProceed.empty() )
      {
          GW_Vertex* pVert = VertexToProceed.front();
          GW_ASSERT( pVert!=NULL );
          VertexToProceed.pop_front();
  
          /* cut the face */
          pCallback( *pVert );
  
          /* add neighbors */
          for( GW_VertexIterator it = pVert->BeginVertexIterator(); it!=pVert->EndVertexIterator(); ++it )
          {
              GW_Vertex* pNewVert = (GW_Vertex*) *it;
              if( pNewVert==NULL )
                  break;
              GW_ASSERT( pNewVert!=NULL );
              if( VertexDone.find(pNewVert->GetID())==VertexDone.end() )
              {
                  VertexToProceed.push_back( pNewVert );
                  VertexDone[ pNewVert->GetID() ] = pNewVert;    // so that it won't be added anymore
              }
          }
      }
  }
  
  /*------------------------------------------------------------------------------*/
  // Name : GW_Mesh::IterateConnectedComponent_Face
  /**
  *  \param  pCallback [VertexIterate_Callback] The function.
  *
  *  Iterate a callback on each vertex of a connected component.
  */
  /*------------------------------------------------------------------------------*/
  void GW_Mesh::IterateConnectedComponent_Face( GW_Face& start_face, FaceIterate_Callback pCallback )
  {
      /* march on the voronoi diagram */
      T_FaceList FaceToProceed;
      FaceToProceed.push_back( &start_face );
      T_FaceMap FaceDone;
      FaceDone[ start_face.GetID() ] = &start_face;
  
  
      while( !FaceToProceed.empty() )
      {
          GW_Face* pFace = FaceToProceed.front();
          GW_ASSERT( pFace!=NULL );
          FaceToProceed.pop_front();
  
          /* cut the face */
          pCallback( *pFace );
  
          /* add neighbors */
          for( GW_U32 i=0; i<3; ++i )
          {
              GW_Face* pNewFace = pFace->GetFaceNeighbor(i);
              if( pNewFace!=NULL && FaceDone.find(pNewFace->GetID())==FaceDone.end() )
              {
                  FaceToProceed.push_back( pNewFace );
                  FaceDone[ pNewFace->GetID() ] = pNewFace;    // so that it won't be added anymore
              }
          }
      }
  }
  
  /*------------------------------------------------------------------------------*/
  // Name : GW_Mesh::ReOrientMesh
  /**
   *  \param  start_face [GW_Face&] The seed face.
   *
   *  Reorient the mesh using this start vertex.
   */
  /*------------------------------------------------------------------------------*/
  void GW_Mesh::ReOrientMesh( GW_Face& start_face )
  {
      /* march on the voronoi diagram */
      T_FaceList FaceToProceed;
      FaceToProceed.push_back( &start_face );
      T_FaceMap FaceDone;
      FaceDone[ start_face.GetID() ] = &start_face;
  
  
      while( !FaceToProceed.empty() )
      {
          GW_Face* pFace = FaceToProceed.front();
          GW_ASSERT( pFace!=NULL );
          FaceToProceed.pop_front();
  
          /* add neighbors */
          for( GW_U32 i=0; i<3; ++i )
          {
              GW_Vertex* pVertDir = pFace->GetVertex(i);    GW_ASSERT( pVertDir!=NULL );
              GW_Face* pNewFace = pFace->GetFaceNeighbor(*pVertDir);
              if( pNewFace!=NULL && FaceDone.find(pNewFace->GetID())==FaceDone.end() )
              {
                  /* find the two other vertices */
                  GW_U32 i1 = (i+1)%3;
                  GW_U32 i2 = (i+2)%3;
                  GW_Vertex* pNewVert[3];
                  pNewVert[0] = pFace->GetVertex(i2);                        GW_ASSERT( pNewVert[0]!=NULL );
                  pNewVert[1] = pFace->GetVertex(i1);                        GW_ASSERT( pNewVert[1]!=NULL );
                  pNewVert[2] = pNewFace->GetVertex(*pNewVert[0], *pNewVert[1]);    GW_ASSERT( pNewVert[2]!=NULL );
                  GW_Face* pNeigh[3];
                  pNeigh[0] = pNewFace->GetFaceNeighbor( *pNewVert[0] );
                  pNeigh[1] = pNewFace->GetFaceNeighbor( *pNewVert[1] );
                  pNeigh[2] = pNewFace->GetFaceNeighbor( *pNewVert[2] );
                  /* reorient the face */
                  pNewFace->SetVertex( *pNewVert[0], *pNewVert[1], *pNewVert[2] );
                  pNewFace->SetFaceNeighbor( pNeigh[0], pNeigh[1], pNeigh[2] );
                  FaceToProceed.push_back( pNewFace );
                  FaceDone[ pNewFace->GetID() ] = pNewFace;    // so that it won't be added anymore
              }
          }
      }
  
      /* check for global orientation (just an heuristic) */
      GW_Face* pFace = this->GetFace(0);    GW_ASSERT( pFace!=NULL );
      GW_Vector3D v = pFace->GetVertex(0)->GetPosition() +
                      pFace->GetVertex(1)->GetPosition() +
                      pFace->GetVertex(2)->GetPosition();
      GW_Vector3D n = pFace->ComputeNormal();
      if( n*v<0 )
          this->FlipOrientation();
  }
  
  
  /*------------------------------------------------------------------------------*/
  // Name : GW_Mesh::FlipOrientation
  /**
   *  Flip each face.
   */
  /*------------------------------------------------------------------------------*/
  void GW_Mesh::FlipOrientation()
  {
      for( GW_U32 i=0; i<this->GetNbrFace(); ++i )
      {
          GW_Face* pFace = this->GetFace(i);    GW_ASSERT( pFace!=NULL );
          pFace->SetVertex( *pFace->GetVertex(1), *pFace->GetVertex(0), *pFace->GetVertex(2) );
          pFace->SetFaceNeighbor( pFace->GetFaceNeighbor(1), pFace->GetFaceNeighbor(0), pFace->GetFaceNeighbor(2) );
      }
  }
  
  
  /*------------------------------------------------------------------------------*/
  // Name : GW_Mesh::ReOrientNormals
  /**
   *  Align normals so that they are in the same direction with face
   *  normals.
   */
  /*------------------------------------------------------------------------------*/
  void GW_Mesh::ReOrientNormals()
  {
      for( GW_U32 i=0; i<this->GetNbrFace(); ++i )
      {
          GW_Face* pFace = this->GetFace(i);    GW_ASSERT( pFace!=NULL );
          GW_Vector3D n = pFace->ComputeNormal();
          for( GW_U32 k=0; k<3; ++k )
          {
              GW_Vector3D& nv = pFace->GetVertex(k)->GetNormal();
              if( nv*n < 0 )
                  nv = -nv;
          }
      }
  }
  
  
  /*------------------------------------------------------------------------------*/
  // Name : GW_Mesh::FlipNormals
  /**
   *  Flip each normal of the mesh.
   */
  /*------------------------------------------------------------------------------*/
  void GW_Mesh::FlipNormals()
  {
      for( GW_U32 i=0; i<this->GetNbrVertex(); ++i )
      {
          GW_Vertex* pVert = this->GetVertex(i);    GW_ASSERT( pVert!=NULL );
          GW_Vector3D& n = pVert->GetNormal();
          n = -n;
      }
  }
  
  
  void GW_Mesh::CheckIntegrity()
  {
      for( GW_U32 i=0; i<this->GetNbrVertex(); ++i )
      {
          GW_Vertex* pVert = this->GetVertex(i); GW_ASSERT( pVert!=NULL );
          GW_Face* pFace = pVert->GetFace();    GW_ASSERT( pFace!=NULL );
          if( pFace!=NULL && pFace->GetVertex(0)!=pVert &&
              pFace->GetVertex(1)!=pVert &&
              pFace->GetVertex(2)!=pVert )
              GW_ASSERT( GW_False );
      }
      for( GW_U32 i=0; i<this->GetNbrFace(); ++i )
      {
          GW_Face* pFace = this->GetFace(i);    GW_ASSERT( pFace!=NULL );
          for( GW_U32 k=0; k<3; ++k )
          {
              GW_U32 k1 = (k+1)%3;
              GW_U32 k2 = (k+2)%3;
              GW_Face* pNeighFace = pFace->GetFaceNeighbor(k);
              GW_Vertex* pV1 = pFace->GetVertex(k1);    GW_ASSERT( pV1!=NULL );
              GW_Vertex* pV2 = pFace->GetVertex(k2);    GW_ASSERT( pV2!=NULL );
              if( pNeighFace!=NULL )
              {
                  GW_ASSERT( pNeighFace->GetFaceNeighbor(*pV1, *pV2)==pFace );
                  GW_ASSERT( pFace->GetFaceNeighbor(*pV1, *pV2)==pNeighFace);
              }
          }
      }
  }
  
  void GW_Mesh::ExtractBoundary( GW_Vertex& seed, T_VertexList& boundary, T_VertexMap* pExtracted )
  {
      GW_ASSERT( seed.IsBoundaryVertex() );
      GW_Vertex* pPrev = NULL;
      GW_Vertex* pCur = &seed;
      GW_Vertex* pNext = NULL;
      GW_U32 num = 0;
      do
      {
          num++;
          boundary.push_back(pCur);
          if( pExtracted!=NULL )
              (*pExtracted)[ pCur->GetID() ] = pCur;
          pNext = NULL;
          for( GW_VertexIterator it = pCur->BeginVertexIterator();
                  (it!=pCur->EndVertexIterator()) && (pNext==NULL); ++it )
          {
              GW_Vertex* pVert = *it;
              if( pVert->IsBoundaryVertex() && pVert!=pPrev )
                  pNext = pVert;
          }
          GW_ASSERT( pNext!=NULL );
          pPrev = pCur;
          pCur = pNext;
      }
      while( pCur!=&seed && pNext!=NULL && num<this->GetNbrVertex() );
  }
  
  
  void GW_Mesh::ExtractAllBoundaries( std::list<T_VertexList>& boundary_list )
  {
      T_VertexMap AlreadyExtracted;
      for( GW_U32 i=0; i<this->GetNbrVertex(); ++i )
      {
          GW_Vertex* pVert = this->GetVertex(i);    GW_ASSERT( pVert!=NULL );
          if( pVert->IsBoundaryVertex() && AlreadyExtracted.find(i)==AlreadyExtracted.end() )
          {
              T_VertexList boundary;
              this->ExtractBoundary( *pVert, boundary, &AlreadyExtracted );
              boundary_list.push_back( boundary );
          }
      }
  }
  
  GW_Float GW_Mesh::GetPerimeter( GW_U32* nNbrBoundaries )
  {
      GW_Float rPerimeter = 0;
      std::list<T_VertexList> boundary_list;
      this->ExtractAllBoundaries( boundary_list );
      if( nNbrBoundaries!=NULL )
          *nNbrBoundaries = (GW_U32) boundary_list.size();
      for( std::list<T_VertexList>::iterator it = boundary_list.begin(); it!=boundary_list.end(); ++it )
      {
          T_VertexList& boundary = *it;
          rPerimeter += GW_Mesh::GetPerimeter( boundary );
      }
  
      return rPerimeter;
  }
  
  GW_Float GW_Mesh::GetPerimeter( T_VertexList& boundary, GW_Bool bCyclic )
  {
      GW_Float rPerimeter = 0;
      GW_Vertex* pPrev = NULL;
      GW_Vertex* pVert = NULL;
      for( IT_VertexList it = boundary.begin(); it!=boundary.end(); ++it )
      {
          pVert = *it;
          if( pPrev!=NULL )
          {
              rPerimeter += (pPrev->GetPosition()-pVert->GetPosition()).Norm();
          }
          pPrev = pVert;
      }
      if( boundary.size()>1 && bCyclic )
      {
          rPerimeter += (boundary.front()->GetPosition()-pVert->GetPosition()).Norm();
      }
      return rPerimeter;
  }
  