C++ Facet_iterator::facet_begin方法代码示例

  void CGALPolyhedronToVTKPolydata_converter(Polyhedron *polyhedron, vtkSmartPointer<vtkPolyData> polydata)
  {
	  //convert from cgal polyhedron to vtk poly data
	  //first the points
	  vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
	  vtkSmartPointer<vtkCellArray> faces = vtkSmartPointer<vtkCellArray>::New();

	  //iterator over all vertices in  polyhedron, add them as points
	  std::map<Vertex_handle, vtkIdType> V;
	  vtkIdType inum = 0;
	  for ( Vertex_iterator v = polyhedron->vertices_begin(); v != polyhedron->vertices_end(); ++v)
	  {
		   points->InsertNextPoint(v->point()[0],v->point()[1],v->point()[2]);
		   V[v] = inum++;
	  }
	  //now iterate over all faces in polyhedron
	  for ( Facet_iterator i = polyhedron->facets_begin(); i != polyhedron->facets_end(); ++i)
	  {
		  Halfedge_around_facet_circulator j = i->facet_begin();
		  faces->InsertNextCell(CGAL::circulator_size(j));
		  do
		  {
			  //get indice of vertex, insert new cell point into faces
			  faces->InsertCellPoint(V[j->vertex()]);

		  } while(++j != i->facet_begin());
	  }

	  // Set the points and faces of the polydata
	  polydata->SetPoints(points);
	  polydata->SetPolys(faces);
  }

int main()
{
    std::vector<Point_3> points;

    points.push_back(Point_3(2.0f, 3.535533905932738f, 3.535533905932737f));
    points.push_back(Point_3(4.0f, 2.0f, 0.0f));
    points.push_back(Point_3(0.0f, 2.0f, 0.0f));
    points.push_back(Point_3(1.0f, 0.0f, 0.0f));
    points.push_back(Point_3(4.0f, 1.414213562373095f, 1.414213562373095f));
    points.push_back(Point_3(0.0f, 1.414213562373095f, 1.414213562373095f));
    points.push_back(Point_3(3.0f, 0.0f, 0.0f));
    points.push_back(Point_3(2.0f, 5.0f, 0.0f));

    Polyhedron P;

    CGAL::convex_hull_3(points.begin(), points.end(), P);

    std::cout << "- Number of vertices  = " << P.size_of_vertices()    << std::endl;
    std::cout << "- Number of edges     = " << P.size_of_halfedges()/2 << std::endl;
    std::cout << "- Number of faces     = " << P.size_of_facets()      << std::endl;

    for ( Facet_iterator i = P.facets_begin(); i != P.facets_end(); ++i)
    {
        Halfedge_facet_circulator j = i->facet_begin();
        CGAL_assertion( CGAL::circulator_size(j) >= 3);
        std::cout << CGAL::circulator_size(j) << ' ';
        do{
            //std::cout << ' ' << std::distance(P.vertices_begin(), j->vertex());
            std::cout << " (" << j->vertex()->point().x() << ' ' << j->vertex()->point().y() << ' ' << j->vertex()->point().z() << ')' << ", ";
        } while ( ++j != i->facet_begin());

        std::cout << std::endl;
    }
    return 0;
}

void renewAttrs( Polyhedron & poly )
{
    // assign IDs to vertices
    for ( Vertex_iterator vi = poly.vertices_begin(); vi != poly.vertices_end(); ++vi ) {
	vi->nCuts()	= 0;
    }

    // assign IDs to faces: the dual vertex and the dual coordinates-the center coordinates
    for ( Facet_iterator fi = poly.facets_begin(); fi != poly.facets_end(); ++fi ) {
	fi->piece()	= NO_INDEX;
	Vector3 sum( 0.0, 0.0, 0.0 );
	Halfedge_facet_circulator hfc = fi->facet_begin();
	do {
	    sum = sum + ( hfc->vertex()->point() - CGAL::ORIGIN );
	} while ( ++hfc != fi->facet_begin() );
	sum = sum / ( double )CGAL::circulator_size( hfc );
	fi->center() = CGAL::ORIGIN + sum;
    }

    // assign the same IDs to the identical/ opposite halfedges
    for ( Halfedge_iterator hi = poly.halfedges_begin(); hi != poly.halfedges_end(); ++hi ) {
	hi->label()		= DEFAULT_LABEL;
	hi->match()		= DEFAULT_LABEL;
	hi->cycle()		= NO_INDEX;
	hi->connect()		= true;
	hi->visit()		= false;
	hi->path()		= NO_INDEX;
	hi->orient()		= true;
	// We individually handle hi->fixed
    }
}

//------------------------------------------------------------------------------
//	Initialize the dual center 
//------------------------------------------------------------------------------
void initAttrs( Polyhedron & poly )
{
    // assign IDs to vertices
    int vID = 0;
    for ( Vertex_iterator vi = poly.vertices_begin(); vi != poly.vertices_end(); ++vi ) {
	vi->id()	= vID++;
	vi->label()	= DEFAULT_LABEL;
	vi->nCuts()	= 0;
    }
    cerr << "Total number of vertices = " << vID << endl;

    // assign IDs to faces: the dual vertex and the dual coordinates-the center coordinates
    int fID = 0;
    for ( Facet_iterator fi = poly.facets_begin(); fi != poly.facets_end(); ++fi ) {
	fi->id()	= fID++;
	// fi->label()	= DEFAULT_LABEL;
	fi->piece()	= NO_INDEX;
	Vector3 sum( 0.0, 0.0, 0.0 );
	Halfedge_facet_circulator hfc = fi->facet_begin();
	do {
	    sum = sum + ( hfc->vertex()->point() - CGAL::ORIGIN );
	} while ( ++hfc != fi->facet_begin() );
	sum = sum / ( double )CGAL::circulator_size( hfc );
	fi->center() = CGAL::ORIGIN + sum;
	// cerr << " Barycenter No. " << fid-1 << " : " << bary[ fid-1 ] << endl;
    }
    cerr << "Total number of facets = " << fID << endl;

    // initialize the halfedge IDs
    for ( Halfedge_iterator hi = poly.halfedges_begin(); hi != poly.halfedges_end(); ++hi ) {
	hi->id() = NO_INDEX;
    }

    // assign the same IDs to the identical/ opposite halfedges
    int eID = 0;
    for ( Halfedge_iterator hi = poly.halfedges_begin(); hi != poly.halfedges_end(); ++hi ) {
	if ( hi->id() == NO_INDEX ) {
	    assert( hi->opposite()->id() == NO_INDEX );
	    hi->id() = eID;
	    hi->opposite()->id() = eID;
	    eID++;
	    hi->label()		= hi->opposite()->label()	= DEFAULT_LABEL;
	    hi->match()		= hi->opposite()->match()	= DEFAULT_LABEL;
	    hi->cycle()		= hi->opposite()->cycle()	= NO_INDEX;
	    hi->weight()	= hi->opposite()->weight()	= 0.0;
	    hi->connect()	= hi->opposite()->connect()	= true;
	    hi->visit()		= hi->opposite()->visit()	= false;
	    // We individually handle hi->fixed
	}
    }
    cerr << "Total number of edges = " << eID << endl;
}

void geometryUtils::subdivide_create_center_vertex(Polyhedron& P, Facet_iterator f) {
    Vector_3 vec(0.0, 0.0, 0.0);
    std::size_t order = 0;
    HF_circulator h = f->facet_begin();
    do {
        vec = vec + (h->vertex()->point() - CGAL::ORIGIN);
        ++order;
    } while (++h != f->facet_begin());
    CGAL_assertion(order >= 3); // guaranteed by definition of polyhedron
    Point_3 center = CGAL::ORIGIN + (vec / static_cast<double> (order));
    Halfedge_handle new_center = P.create_center_vertex(f->halfedge());
    new_center->vertex()->point() = center;
}

F& apply_function_object_polyhedron (Polyhedron &P,
                                     F &f) {
    typedef typename Polyhedron::Halfedge_around_facet_circulator Hafc;
    typedef typename Polyhedron::Facet_iterator Facet_iterator;

    f.reset();

    for (Facet_iterator fit = P.facets_begin();
         fit != P.facets_end();
         fit++) {
        Hafc h0 = fit->facet_begin(), hf = h0--, hs = hf;
        hs ++;

        while (1) {
            // Apply 'f' on each triangle of the polyhedron's facet
            f ( h0->vertex()->point(),
                hf->vertex()->point(),
                hs->vertex()->point() );

            if (hs == h0)
                break;

            hs++; hf++;
        }
    }

    f.end();

    return f;
}

 int main(int argc, const char **argv )
 {
   std::vector<Point> points;
   CGAL::Random_points_on_sphere_3<Point> g;

   size_t N = 0;
   if (argc > 1)
     N = atof(argv[1]);
   N = std::max(size_t(100), N);

   for (size_t i = 0; i < N; ++i)
     points.push_back(rescale(*g++));

   for (size_t n = 0; n < 100; ++n)
     {
       std::cerr << "step " << n << ":\n\t";
       lloyd_step(points);
     }

   Polyhedron P;
   CGAL::convex_hull_3(points.begin(), points.end(), P);

   CGAL::set_ascii_mode( std::cout);
   std::cout << "OFF" << std::endl << P.size_of_vertices() << ' '
	     << P.size_of_facets() << " 0" << std::endl;
   std::copy( P.points_begin(), P.points_end(),
	      std::ostream_iterator<Point>( std::cout, "\n"));
   for (  Facet_iterator i = P.facets_begin(); i != P.facets_end(); ++i) {
     Halfedge_facet_circulator j = i->facet_begin();
     // Facets in polyhedral surfaces are at least triangles.
     CGAL_assertion( CGAL::circulator_size(j) >= 3);
     std::cout << CGAL::circulator_size(j) << ' ';
     do {
       std::cout << ' ' << std::distance(P.vertices_begin(), j->vertex());
     } while ( ++j != i->facet_begin());
     std::cout << std::endl;
   }

   std::ofstream os ("test.cloud");
   std::copy(points.begin(), points.end(),
	     std::ostream_iterator<Point>(os, "\n"));
}

void geometryUtils::computeGaussCurvature(Polyhedron* P) {
    for (Facet_iterator i = P->facets_begin(); i != P->facets_end(); i++) {

        int e = 0;
        Halfedge_around_facet_circulator edge = i->facet_begin();
        do {
            i->kappa[e] = computeLocalGaussCurvature(edge->vertex());
            int H = floor(i->kappa[e] * geometryUtils::kappaMax + 180);

            //            std::cout << H;
            //            std::cout << " ";

            i->color[e] = HSVtoRGB(H, 1.0, 1.0);

            e++;
        } while (++edge != i->facet_begin());
        //        std::cout << edge->kappa;
        //        std::cout << " ";
    }
};

void Scene_polyhedron_item::compute_elements()
{

    v_poly.resize(0);
    v_edge.resize(0);
    normal_flat.resize(0);
    normal_smooth.resize(0);
    Polyhedron& polyhedron =*poly;
    //FACETS

    typedef Polyhedron::Traits	    Kernel;
    typedef Kernel::Point_3	    Point;
    typedef Kernel::Vector_3	    Vector;
    typedef Polyhedron::Facet_iterator Facet_iterator;
    typedef Polyhedron::Halfedge_around_facet_circulator HF_circulator;



    Facet_iterator f;
    for(f = polyhedron.facets_begin();
      f != polyhedron.facets_end();
      f++)
    {

      // If Flat shading: 1 normal per polygon

        Vector n = CGAL::Polygon_mesh_processing::compute_face_normal(f, polyhedron);

        normal_flat.push_back(n.x()); normal_flat.push_back(n.y()); normal_flat.push_back(n.z());
        normal_flat.push_back(n.x()); normal_flat.push_back(n.y()); normal_flat.push_back(n.z());
        normal_flat.push_back(n.x()); normal_flat.push_back(n.y()); normal_flat.push_back(n.z());


      // revolve around current face to get vertices
      HF_circulator he = f->facet_begin();
      HF_circulator end = he;
      CGAL_For_all(he,end)
      {

          Vector n = CGAL::Polygon_mesh_processing::compute_vertex_normal(he->vertex(), polyhedron);
          normal_smooth.push_back(n.x()); normal_smooth.push_back(n.y()); normal_smooth.push_back(n.z());

        const Point& p = he->vertex()->point();
        v_poly.push_back(p.x()); v_poly.push_back(p.y()); v_poly.push_back(p.z());
      }
    }

  double AreaFacetTriangleSeg(Facet_iterator &f)
	{
		Halfedge_around_facet_circulator pHalfedge = f->facet_begin();
		Point3d P = pHalfedge->vertex()->point();
		Point3d Q = pHalfedge->next()->vertex()->point();
		Point3d R = pHalfedge->next()->next()->vertex()->point();

		Vector PQ=Q-P;
                //Vector PR=R-P; // MT
		Vector QR=R-Q;


		Vector normal	=	CGAL::cross_product(PQ,QR);
		double area=0.5*sqrt(normal*normal);

		return area;

	}

void geometryUtils::renderPolyhedron(Polyhedron * pmesh) {

    glBegin(GL_TRIANGLES);
    for (Facet_iterator i = pmesh->facets_begin(); i != pmesh->facets_end(); i++) {


        Halfedge_around_facet_circulator j = i->facet_begin();
        glColor3d(i->color[0].R, i->color[0].G, i->color[0].B);
        glVertex3d(CGAL::to_double(j->vertex()->point().x()), CGAL::to_double(j->vertex()->point().y()), CGAL::to_double(j->vertex()->point().z()));

        glColor3d(i->color[1].R, i->color[1].G, i->color[1].B);
        glVertex3d(CGAL::to_double(j->next()->vertex()->point().x()), CGAL::to_double(j->next()->vertex()->point().y()), CGAL::to_double(j->next()->vertex()->point().z()));

        glColor3d(i->color[2].R, i->color[2].G, i->color[2].B);
        glVertex3d(CGAL::to_double(j->next()->next()->vertex()->point().x()), CGAL::to_double(j->next()->next()->vertex()->point().y()), CGAL::to_double(j->next()->next()->vertex()->point().z()));
    }
    glEnd();
}

//.........这里部分代码省略.........
	E = p.Draw();
	P = n.Draw();
	Morph Mor(E,50);
	Mor.Do(P);*/

	//Noise No(5,0.3,0,s.Center,Z_ax);
	//No.Do(P);

	//Skew Sk(30,s.Center,Z_ax,false,20,-20);
	//Sk.Do(P);

	//Smooth Sm(1);
	//Sm.Do(P);

	//Spherify Sph(50);
	//Sph.Do(P);

	//Squeeze Sq(-30,s.Center,Z_ax,false,10,0);
	//Sq.Do(P);

	//Stretch St(-20,s.Center,Z_ax,true,50,-50);
	//St.Do(P);
	
	//Taper Ta(3,s.Center,X_ax,false,20,-20);
	//Ta.Do(P);

	//Twist Tw(270,s.Center,Z_ax,true,-5,15);
	//Tw.Do(P);

	/*Box_3 B(20, 30, 60, 20, 30, 30); 
	Polyhedron P = B.Draw();
	Twist Tw1(270, B.Center, Z_ax, true, 30, 10);
	Twist Tw2(-270, B.Center, Z_ax, true, -10, -30);
	Stretch St(30, B.Center, Z_ax, true, 20,-20);
	Squeeze Sq(15, B.Center, Z_ax);
	Tw1.Do(P);
	Tw2.Do(P);
	Sq.Do(P);
	St.Do(P);*/

	std::ofstream of("C:\\123.off");

	Box_3 B(20, 30, 60, 20, 30, 30); 
Polyhedron P = B.Draw();
Twist Tw1(270, B.Center, Z_ax, true, 30, 10);
Twist Tw2(-270, B.Center, Z_ax, true, -10, -30);
Stretch St(30, B.Center, Z_ax, true, 20,-20);
Squeeze Sq(15, B.Center, Z_ax);
Tw1.Do(P);
Tw2.Do(P);
Sq.Do(P);
St.Do(P);



	//// Write polyhedron in Object File Format (OFF).
	CGAL::set_ascii_mode( of );
	of << "OFF" << std::endl << P.size_of_vertices() << ' ' << P.size_of_facets() << " 0" << std::endl;
	std::copy( P.points_begin(), P.points_end(), std::ostream_iterator<Point_3>( of, "\n"));

	for (  Facet_iterator i = P.facets_begin(); i != P.facets_end(); ++i) 
	{
		Halfedge_facet_circulator j = i->facet_begin();
		// Facets in polyhedral surfaces are at least triangles.
		CGAL_assertion( CGAL::circulator_size(j) >= 3);
		of << CGAL::circulator_size(j) << ' ';
		
		do 
		{
			of << ' ' << std::distance(P.vertices_begin(), j->vertex());
		} while ( ++j != i->facet_begin());
		
		of << std::endl;
	}

	/* Write polyhedron in (OBJ).
	CGAL::set_ascii_mode( oof );
	oof << "# " << P.size_of_vertices() << ' ' << std::endl <<"# "<< P.size_of_facets() << std::endl;
	oof<<"v ";
	std::copy( P.points_begin(), P.points_end(), std::ostream_iterator<Point_3>( oof, "\nv "));
	oof<<"_ _ _\n";
	for (  Facet_iterator i = P.facets_begin(); i != P.facets_end(); ++i)
	{
		Halfedge_facet_circulator j = i->facet_begin();
		// Facets in polyhedral surfaces are at least triangles.
		//CGAL_assertion( CGAL::circulator_size(j) >= 3);
       
		oof << 'f' << ' ';
       
		do
		{
			oof << ' ' << std::distance(P.vertices_begin(), j->vertex())+1;
		} while ( ++j != i->facet_begin());
       
		oof << std::endl;
	}
	*/

	return 0;
}

void
Scene_textured_polyhedron_item::compute_normals_and_vertices(void)
{
    positions_facets.resize(0);
    positions_lines.resize(0);
    textures_map_facets.resize(0);
    textures_map_lines.resize(0);
    normals.resize(0);

    typedef ::EPIC_kernel Kernel;
    typedef CGAL::Textured_items Items;
    typedef Kernel::Point_3 Point;
    typedef Kernel::Vector_3 Vector;
    typedef CGAL::Polyhedron_3<Kernel,Items> Base;

    typedef Base::Halfedge_around_facet_circulator Halfedge_around_facet_circulator;
    typedef Base::Edge_iterator Edge_iterator;
    typedef Base::Facet Facet;
    typedef Base::Facet_iterator Facet_iterator;

    //Facets
    Facet_iterator f = poly->facets_begin();

    for(f = poly->facets_begin();
        f != poly->facets_end();
        f++)
    {

        Halfedge_around_facet_circulator he = f->facet_begin();
        Halfedge_around_facet_circulator end = he;
        CGAL_For_all(he,end)
        {

            // If Flat shading:1 normal per polygon added once per vertex
            if (cur_shading == Flat || cur_shading == FlatPlusEdges)
            {

                Vector n = CGAL::Polygon_mesh_processing::
                  compute_face_normal(f, static_cast<Base&>(*poly));
                normals.push_back(n[0]);
                normals.push_back(n[1]);
                normals.push_back(n[2]);
            }

            // If Gouraud shading: 1 normal per vertex
            else if(cur_shading == Gouraud)
            {

                const Facet::Normal_3& n = he->vertex()->normal();
                normals.push_back(n[0]);
                normals.push_back(n[1]);
                normals.push_back(n[2]);
            }

            //position
            const Point& p = he->vertex()->point();
            positions_facets.push_back(p.x());
            positions_facets.push_back(p.y());
            positions_facets.push_back(p.z());
            positions_facets.push_back(1.0);

            const double u = he->vertex()->u();
            const double v = he->vertex()->v();
            textures_map_facets.push_back(u);
            textures_map_facets.push_back(v);
        }


    }
    //Lines
    typedef Kernel::Point_3		Point;
    typedef Base::Edge_iterator	Edge_iterator;

    Edge_iterator he;

    for(he = poly->edges_begin();
        he != poly->edges_end();
        he++)
    {

        const Point& a = he->vertex()->point();
        const Point& b = he->opposite()->vertex()->point();
        positions_lines.push_back(a.x());
        positions_lines.push_back(a.y());
        positions_lines.push_back(a.z());
        positions_lines.push_back(1.0);

        const double u = he->vertex()->u();
        const double v = he->vertex()->v();
        textures_map_lines.push_back(u);
        textures_map_lines.push_back(v);

        positions_lines.push_back(b.x());
        positions_lines.push_back(b.y());
        positions_lines.push_back(b.z());
        positions_lines.push_back(1.0);

        const double ou = he->opposite()->vertex()->u();
        const double ov = he->opposite()->vertex()->v();
        textures_map_lines.push_back(ou);
//.........这里部分代码省略.........

void ElPoly::DefineSkin(int NSample){
  std::list<Weighted_point> l;
  FT shrinkfactor = 0.5;
  double *Plot  = new double[pNType()*CUBE(NSample)];
  double *Count = new double[CUBE(NSample)];
  double Thre = 10.;
  double Radius = pEdge(0)/(double)NSample;
  for(int p=0;p<pNPart();p++){
    int t = pType(p);
    int vx = (int)(pPos(p,0)/pEdge(0)*NSample);
    int vy = (int)(pPos(p,1)/pEdge(1)*NSample);
    int vz = (int)(pPos(p,2)/pEdge(2)*NSample);
    int vTot = (vz*NSample+vy)*NSample+vx;
    Plot[vTot*pNType()+t] += 1.;
  }
  double *Norm = (double *)calloc(pNType(),sizeof(double));
  for(int t=0;t<pNType();t++){
    for(int v=0;v<CUBE(NSample);v++){
      if(Norm[t] < Plot[v*pNType()+t])
	Norm[t] = Plot[v*pNType()+t];
    }
    Norm[t] = Norm[t] <= 0. ? 1. : Norm[t];
  }
  for(int vx=0;vx<NSample;vx++){
    double x = vx*pEdge(0)/(double)NSample;
    for(int vy=0;vy<NSample;vy++){
      double y = vy*pEdge(1)/(double)NSample;
      for(int vz=0;vz<NSample;vz++){
	double z = vz*pEdge(2)/(double)NSample;
	int vTot = (vz*NSample+vy)*NSample+vx;
	if(Plot[vTot*pNType()] > Thre){
	  l.push_front(Weighted_point(Bare_point(x,y,z),Radius));
	}
      }
    }
  }

  Polyhedron Polyhe;

  Skin_surface_3 skin_surface(l.begin(), l.end(), shrinkfactor);
  CGAL::mesh_skin_surface_3(skin_surface, Polyhe);

  //  CGAL::subdivide_skin_surface_mesh_3(skin_surface, Polyhe);

  // std::ofstream out("mesh.off");
  // out << Polyhe;

  glDeleteLists(Dr->Particles,1);
  Dr->Particles = glGenLists(1);
  glNewList(Dr->Particles,GL_COMPILE);

  // Polyhedron::Facet_iterator fcUp = Polyhe.facets_begin();
  // for(;fcUp != Polyhe.facets_end(); ++fcUp){
  //   Polyhedron::Supports_facet_halfedge = fcUp.halfedge();
  //   //Halfedge_around_facet_circulator heUp = fcUp.halfedge();
  // }

  // for (Vertex_iterator vit = Polyhe.vertices_begin();vit != Polyhe.vertices_end(); vit++){
  //   // Vector n = policy.normal(vit);
  //   // n = n/sqrt(n*n);
  //   cout << vit->point() << std::endl;
  //   Halfedge_iterator heUp = Polyhe.halfedges_begin();
  //   for(;heUp != Polyhe.halfedges_end(); ++heUp){
  //     //Polyhedron::Halfedge_handle Half = *heUp;
  //     Vertex_handle veUp = heUp->vertex();
  //     K::Point_3 pf1 = vit->point();
  //   }
  // }
  CGAL::Inverse_index<Vertex_handle> index(Polyhe.vertices_begin(),
  					   Polyhe.vertices_end());

  for(Facet_iterator fi = Polyhe.facets_begin();fi != Polyhe.facets_end(); ++fi) {
    HFC hc = fi->facet_begin();
    HFC hc_end = hc;
    Polyhedron::Vertex_handle vf1 = (*hc).vertex();
    hc++;
    Polyhedron::Vertex_handle vf2 = (*hc).vertex();
    hc++;
    Polyhedron::Vertex_handle vf3 = (*hc).vertex();
    hc++;
    K::Point_3 pf1 = vf1->point();
    K::Point_3 pf2 = vf2->point();
    K::Point_3 pf3 = vf3->point();
    Vettore v1(pf1.x(),pf1.y(),pf1.z());
    Vettore v2(pf2.x(),pf2.y(),pf2.z());
    Vettore v3(pf3.x(),pf3.y(),pf3.z());
    Vettore vN(3);
    v1.Mult(InvScaleUn);
    v2.Mult(InvScaleUn);
    v3.Mult(InvScaleUn);
    vN = (v1-v2) ^ (v3-v2);
    //if(vN.Norm() > 2.*AreaMean) continue;
    double Sfumatura = .3*Mat->Casuale();
    glColor4f(0.1,.4+Sfumatura,0.2,1.);
    //glColor4f(HueUp[p].r,HueUp[p].g,HueUp[p].b,HueUp[p].a);
    DrTria(&v1,&v2,&v3,&vN);
    glColor4f(1.,.0,0.,1.);
    DrTriaContour(&v1,&v2,&v3);


//.........这里部分代码省略.........

//.........这里部分代码省略.........
  mwSize nrowsTri = mxGetM(inTRI->pm);
  mwSize nrowsX = mxGetM(inX->pm);

#ifdef DEBUG  
  std::cout << "Number of facets read = " << mesh.size_of_facets() << std::endl;
  std::cout << "Number of vertices read = " << mesh.size_of_vertices() << std::endl;
#endif

  if (nrowsTri != mesh.size_of_facets()) {
    mexErrMsgTxt(("Input " + inTRI->name + ": Number of triangles read into mesh different from triangles provided at the input").c_str());
  }
  if (nrowsX != mesh.size_of_vertices()) {
    mexErrMsgTxt(("Input " + inX->name + ": Number of vertices read into mesh different from vertices provided at the input").c_str());
  }

  // sort halfedges such that the non-border edges precede the
  // border edges. We need to do this before any halfedge iterator
  // operations are valid
  mesh.normalize_border();
  
#ifdef DEBUG
  std::cout << "Number of border halfedges = " << mesh.size_of_border_halfedges() << std::endl;
#endif
  
  // number of holes we have filled
  mwIndex n = 0;

  // a closed mesh with no holes will have no border edges. What we do
  // is grab a border halfedge and close the associated hole. This
  // makes the rest of the iterators invalid, so we have to normalize
  // the mesh again. Then we iterate, looking for a new border
  // halfedge, filling the hole, etc.
  //
  // Note that confusingly, mesh.border_halfedges_begin() gives a
  // pointer to the halfedge that is NOT a border in a border
  // edge. The border halfedge is instead
  // mesh.border_halfedges_begin()->opposite()
  while (!mesh.is_closed()) {
    
    // exit if user pressed Ctrl+C
    ctrlcCheckPoint(__FILE__, __LINE__);

    // get the first hole we can find, and close it
    mesh.fill_hole(mesh.border_halfedges_begin()->opposite());

    // increase the counter of number of holes we have filled
    n++;

    // renormalize mesh so that halfedge iterators are again valid
    mesh.normalize_border();

  }

  // split all facets to triangles
  CGAL::triangulate_polyhedron<Polyhedron>(mesh);

  // copy output with number of holes filled
  std::vector<double> nout(1, n);
  matlabExport->CopyVectorOfScalarsToMatlab<double, std::vector<double> >(outN, nout, 1);

  // allocate memory for Matlab outputs
  double *tri = matlabExport->AllocateMatrixInMatlab<double>(outTRI, mesh.size_of_facets(), 3);

   // extract the triangles of the solution
  // snippet adapted from CgalMeshSegmentation.cpp

  // vertices coordinates. Assign indices to the vertices by defining
  // a map between their handles and the index
  std::map<Vertex_handle, int> V;
  int inum = 0;
  for(Vertex_iterator vit = mesh.vertices_begin(); vit != mesh.vertices_end(); ++vit) {

    // save to internal list of vertices
    V[vit] = inum++;

  }  

  // triangles given as (i,j,k), where each index corresponds to a vertex in x
  mwIndex row = 0;
  for (Facet_iterator fit = mesh.facets_begin(); fit != mesh.facets_end(); ++fit, ++row) {

    if (fit->facet_degree() != 3) {
      std::cerr << "Facet has " << fit->facet_degree() << " edges" << std::endl;
      mexErrMsgTxt("Facet does not have 3 edges");
    }

    // go around the half-edges of the facet, to extract the vertices
    Halfedge_around_facet_circulator heit = fit->facet_begin();
    int idx = 0;
    do {
      
      // extract triangle indices and save to Matlab output
      // note that Matlab indices go like 1, 2, 3..., while C++ indices go like 0, 1, 2...
      tri[row + idx * mesh.size_of_facets()] = 1 + V[heit->vertex()];
      idx++;

    } while (++heit != fit->facet_begin());
    
  }
}