C++ Polyhedron::size_of_facets方法代码示例

void two_tetrahedrons()
{
  Polyhedron a;

  make_tetrahedron(a, 
                   Point(1.0, 0.0, 0.0),
                   Point(2.0, 0.0, 0.0),
                   Point(1.5, 1.0, 0.0),
                   Point(1.5, .5, 10.0));

  Polyhedron b;
  make_tetrahedron(b,
                   Point(0.0, 0., .5),
                   Point(0.0, 0.0, 1.5),
                   Point(0.0, 1.0, 1.0),
                   Point(10.0, .5, 1.0));

  if (a.is_pure_triangle())
    std::cout << "a is pure triangle" << std::endl;

  if (b.is_pure_triangle())
    std::cout << "b is pure triangle" << std::endl;

  Polyhedron &biggest = a.size_of_facets() > b.size_of_facets() ? a : b;
  Polyhedron &smallest = a.size_of_facets() > b.size_of_facets() ? b : a;

  std::list<std::list<boost::tuple<Facet_handle, Facet_handle, Segment> > > polylines;
  {
    std::list<boost::tuple<Facet_handle, Facet_handle, Segment> > intersections;
    compute_intersections(biggest, smallest, std::back_inserter(intersections));

    for (std::list<boost::tuple<Facet_handle, Facet_handle, Segment> >::iterator it = intersections.begin();
         it != intersections.end(); it++)
    {
      {
        Halfedge_handle h = it->get<0>()->halfedge();
        Triangle t(h->vertex()->point(), h->next()->vertex()->point(), h->next()->next()->vertex()->point());
        assert(t.has_on(it->get<2>().source()));
        assert(t.has_on(it->get<2>().target()));
      }
      {
        Halfedge_handle h = it->get<1>()->halfedge();
        Triangle t(h->vertex()->point(), h->next()->vertex()->point(), h->next()->next()->vertex()->point());
        assert(t.has_on(it->get<2>().source()));
        assert(t.has_on(it->get<2>().target()));
      }
    }
    sort_polylines<Polyhedron>(biggest, smallest, intersections, polylines);
  }

  std::list<std::vector<typename Polyhedron::Halfedge_handle> > intersection_list;

  split_facets<Polyhedron, 0>(biggest,  polylines, intersection_list);
  //split_facets<Polyhedron, 1>(smallest, polylines);

}

void mergeCoplanar(Polyhedron& p,bool step2) {
	int facetsBefore =  p.size_of_facets();
	
	p.normalize_border();
	if(!p.size_of_border_halfedges()) {
		// Calculate normals only once in advace! so all tris should be coplanar with the original
		std::transform( p.facets_begin(), p.facets_end(),p.planes_begin(),Plane_equation());			// Calculate plane equations (only works on tri<- = bs)=true
		bool coplanarFound = true;
		std::vector<Polyhedron::Halfedge_handle> skipHEs;
		while (coplanarFound) {			
			coplanarFound = false;																		// Set coplanarFound false
			int percCount = 1;
			for (Polyhedron::Halfedge_iterator hit = p.halfedges_begin(); hit != p.halfedges_end(); ++hit,++percCount){ // Loop through all halfedges
				if (is_coplanar(hit,true)){																// If coplanar and equals semantics
					Polyhedron::Halfedge_handle removeMe = hit;
					while (CGAL::circulator_size(removeMe->vertex_begin()) < 3)							// Mover handle to beginning of linestring
						removeMe = removeMe->next();
					
					bool jcnh = false;
					if (!step2) jcnh = joinCreatesNoHole (hit);
					else		jcnh = joinCreatesNoHole2(hit);
					if (jcnh){														// If no holes will be created

						std::cout << "\rFacets   before/after: "<<facetsBefore<<" -> "<< p.size_of_facets()<<". ("<<100*percCount/p.size_of_halfedges()<<"%)";

						while (CGAL::circulator_size(removeMe->opposite()->vertex_begin()) < 3)			// Join vertexes until at the other end of linestring
							if (removeMe->facet_degree()>3 && removeMe->opposite()->facet_degree()>3)
								removeMe = (p.join_vertex(removeMe))->next()->opposite();
							else																		// One of the faces turned into a triangle ->remove center vertex
								break;																	
						if (CGAL::circulator_size(removeMe->opposite()->vertex_begin()) < 3)			// Remove remained of the border
							p.erase_center_vertex(removeMe->opposite());								// if two segments remain remove center point
						else
							p.join_facet(removeMe);														// if one segment remains join facets
						coplanarFound = true;
						break;
					} else { // simplify border, but how to do this safely? not optimal solution implemented. Should do: add inward offseted point of intersection etc.
						if (std::find(skipHEs.begin(), skipHEs.end(),hit)!=skipHEs.end()) {					// Skip if hit in skipList
							while (CGAL::circulator_size(removeMe->opposite()->vertex_begin()) < 3)	{		// Join vertexes until at the other end of linestring
								if (removeMe->facet_degree()>3 && removeMe->opposite()->facet_degree()>3)
									if (triDoesNotIntersectFacet(removeMe))									// if tri reMe,reME->prev does not intersect left or right facet
										removeMe = (p.join_vertex(removeMe))->next()->opposite();			// remove removeME
									else {
										skipHEs.push_back(removeMe);
										skipHEs.push_back(removeMe->opposite());
										removeMe = removeMe->prev();										// move removeME one halfedge back
									}
								else break;																	// stop if only a triangle remains or at other end																
							}	
							skipHEs.push_back(removeMe);
							skipHEs.push_back(removeMe->opposite());
	}	}	}	}	}	}
	//if (!step2) mergeCoplanar(p,true);
	//else 
		std::cout << "\rFacets   before/after: "<<facetsBefore<<" -> "<< p.size_of_facets()<<". (100%)"<<std::endl;
}

int main()
{
  std::stringstream ss;
  ss << "\
OFF\n6 4 0\n\
0 1 0\n\
0 0 0\n\
1 0 0\n\
1 1 0\n\
2 1 0\n\
2 0 0\n\
3 0 1 2\n\
3 0 2 3\n\
3 2 5 3\n\
3 5 4 3\n";

  Polyhedron P;
  ss >> P;

  assert( P.size_of_vertices() == 6);
  assert( P.size_of_facets() == 4);
  assert( P.is_valid() );

  //consider vertex 3 and set its halfedge to be on the border
  Polyhedron::Vertex_iterator vit=P.vertices_begin();
  std::advance(vit, 3);
  assert( vit->point() == K::Point_3(1, 1, 0) );
  Polyhedron::Halfedge_handle h=vit->halfedge();
  while ( !h->is_border() )
    h=h->next()->opposite();
  vit->VBase::set_halfedge(h);
  assert( vit->halfedge()->vertex() == vit );

  //consider vertex 4 and set its halfedge to be on the border
  ++vit;
  assert( vit->point() == K::Point_3(2, 1, 0) );
  h=vit->halfedge();
  while ( !h->is_border() )
    h=h->next()->opposite();
  vit->VBase::set_halfedge(h);
  assert( vit->halfedge()->vertex() == vit );

  //try to append a facet
  Appender modifier;
  P.delegate(modifier);

  assert( P.size_of_vertices() == 7);
  assert( P.size_of_facets() == 5);
  assert( P.is_valid() );
}

void geometryUtils::subdivide(Polyhedron& P) {
    if (P.size_of_facets() == 0)
        return;
    // We use that new vertices/halfedges/facets are appended at the end.
    std::size_t nv = P.size_of_vertices();
    Vertex_iterator last_v = P.vertices_end();
    --last_v; // the last of the old vertices
    Edge_iterator last_e = P.edges_end();
    --last_e; // the last of the old edges
    Facet_iterator last_f = P.facets_end();
    --last_f; // the last of the old facets
    Facet_iterator f = P.facets_begin(); // create new center vertices
    do {
        geometryUtils::subdivide_create_center_vertex(P, f);
    } while (f++ != last_f);
    std::vector<Point_3> pts; // smooth the old vertices
    pts.reserve(nv); // get intermediate space for the new points
    ++last_v; // make it the past-the-end position again
    std::transform(P.vertices_begin(), last_v, std::back_inserter(pts),
            Smooth_old_vertex());
    std::copy(pts.begin(), pts.end(), P.points_begin());
    Edge_iterator e = P.edges_begin(); // flip the old edges
    ++last_e; // make it the past-the-end position again
    while (e != last_e) {
        Halfedge_handle h = e;
        ++e; // careful, incr. before flip since flip destroys current edge
        geometryUtils::subdivide_flip_edge(P, h);
    };
    CGAL_postcondition(P.is_valid());
};

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 two_boxes()
{
  Polyhedron a;
  make_box(0,0,0, 4, 5, 2, a);

  Polyhedron b;
  make_box(1, 1, -1, 2, 2, 1, b);

  if (a.is_pure_triangle())
    std::cout << "a is pure triangle" << std::endl;

  if (b.is_pure_triangle())
    std::cout << "b is pure triangle" << std::endl;

  Polyhedron &biggest = a.size_of_facets() > b.size_of_facets() ? a : b;
  Polyhedron &smallest = a.size_of_facets() > b.size_of_facets() ? b : a;

  std::list<std::list<boost::tuple<Facet_handle, Facet_handle, Segment> > > polylines;
  {
    std::list<boost::tuple<Facet_handle, Facet_handle, Segment> > intersections;
    compute_intersections(biggest, smallest, std::back_inserter(intersections));

    for (std::list<boost::tuple<Facet_handle, Facet_handle, Segment> >::iterator it = intersections.begin();
         it != intersections.end(); it++)
    {
      {
        Halfedge_handle h = it->get<0>()->halfedge();
        Triangle t(h->vertex()->point(), h->next()->vertex()->point(), h->next()->next()->vertex()->point());
        assert(t.has_on(it->get<2>().source()));
        assert(t.has_on(it->get<2>().target()));
      }
      {
        Halfedge_handle h = it->get<1>()->halfedge();
        Triangle t(h->vertex()->point(), h->next()->vertex()->point(), h->next()->next()->vertex()->point());
        assert(t.has_on(it->get<2>().source()));
        assert(t.has_on(it->get<2>().target()));
      }
    }
    sort_polylines<Polyhedron>(biggest, smallest, intersections, polylines);
  }

  std::list<std::vector<Halfedge_handle> > a_edges;
  split_facets<Polyhedron, 0>(biggest, polylines, a_edges);
  check_splitting<Polyhedron, 0>(biggest, polylines, a_edges);
  //split_facets<Polyhedron, 1>(smallest, /* smallest, */ polylines);
}

void intersection( const Polyhedron& P) {
    std::vector<Box> boxes;
    boxes.reserve( P.size_of_facets());
    for ( Facet_const_iterator i = P.facets_begin(); i != P.facets_end(); ++i){
        boxes.push_back(
            Box( i->halfedge()->vertex()->point().bbox()
               + i->halfedge()->next()->vertex()->point().bbox()
               + i->halfedge()->next()->next()->vertex()->point().bbox(),
                 i));
    }
    std::vector<const Box*> box_ptr;
    box_ptr.reserve( P.size_of_facets());
    for ( std::vector<Box>::iterator j = boxes.begin(); j != boxes.end(); ++j){
        box_ptr.push_back( &*j);
    }
    CGAL::box_self_intersection_d( box_ptr.begin(), box_ptr.end(),
                                   Intersect_facets(), std::ptrdiff_t(2000));
}

 Polyhedron generate_polyhedron(
         const MatrixFr& vertices, const MatrixIr& faces) {
     Polyhedron P;
     PolyhedronBuilder<HalfedgeDS> triangle(vertices, faces);
     P.delegate(triangle);
     assert(vertices.rows() == P.size_of_vertices());
     assert(faces.rows() == P.size_of_facets());
     return P;
 }

//------------------------------------------------------------------------------
//	Compute and store the results of dual graphs into appropriate variables
//
//  computeDual( mesh, whole, bary, mid ); 
//------------------------------------------------------------------------------
void computeDual( Polyhedron & poly, Graph & dual,
		  vector< Point3 > & bary, vector< Point3 > & mid )
{
    // initialize the array of face barycenters: refer to the point at which the gravitational forces exerted by 2 objects are equal

    bary.clear();
    bary.resize( poly.size_of_facets() );

    // initialize the dual graph with dual vertices
    dual.clear();
    int fid = 0;
    for ( Facet_iterator fi = poly.facets_begin(); fi != poly.facets_end(); ++fi ) {
	// add a vertex: each face equals to each dual vertex
	add_vertex( dual );
	bary[ fid++ ] = fi->center(); 
	// cerr << " Barycenter No. " << fid-1 << " : " << bary[ fid-1 ] << endl;
    }

    int size_of_edges = poly.size_of_halfedges() / 2; //redundant

    // initialize the array of midpoints
    mid.clear();
    mid.resize( size_of_edges );

    // initialize the array of lengths
    // length.clear();
    // length.resize( size_of_edges );

    // construct the connecitivity of the dual graph
    for ( Halfedge_iterator hi = poly.halfedges_begin(); hi != poly.halfedges_end(); ++hi ) {
	int origID = hi->facet()->id(); // the face
	int destID = hi->opposite()->facet()->id(); // the neighbor face
	Point3 origCoord = hi->vertex()->point();
	Point3 destCoord = hi->opposite()->vertex()->point();
	Vector3 dispVec = origCoord - destCoord;
	Point3 midCoord = destCoord + 0.5 * dispVec;
	// double curLength = sqrt( dispVec.squared_length() ); // the length between two connected vertices
#ifdef DEBUG
	cerr << origID << " -- " << destID << endl;
#endif	// DEBUG	
	
	if ( origID < destID )
		add_edge( origID, destID, hi->id(), dual );  
		
	mid[ hi->id() ] = midCoord;
	hi->mid() = hi->opposite()->mid() = midCoord;
	hi->weight() = hi->opposite()->weight() = 1.0;
    }
}

IGL_INLINE void igl::copyleft::cgal::polyhedron_to_mesh(
  const Polyhedron & poly,
  Eigen::MatrixXd & V,
  Eigen::MatrixXi & F)
{
  using namespace std;
  V.resize(poly.size_of_vertices(),3);
  F.resize(poly.size_of_facets(),3);
  typedef typename Polyhedron::Vertex_const_iterator Vertex_iterator;
  std::map<Vertex_iterator,size_t> vertex_to_index;
  {
    size_t v = 0;
    for(
      typename Polyhedron::Vertex_const_iterator p = poly.vertices_begin();
      p != poly.vertices_end();
      p++)
    {
      V(v,0) = p->point().x();
      V(v,1) = p->point().y();
      V(v,2) = p->point().z();
      vertex_to_index[p] = v;
      v++;
    }
  }
  {
    size_t f = 0;
    for(
      typename Polyhedron::Facet_const_iterator facet = poly.facets_begin();
      facet != poly.facets_end();
      ++facet)
    {
      typename Polyhedron::Halfedge_around_facet_const_circulator he = 
        facet->facet_begin();
      // Facets in polyhedral surfaces are at least triangles.
      assert(CGAL::circulator_size(he) == 3 && "Facets should be triangles");
      size_t c = 0;
      do {
        //// This is stooopidly slow
        // F(f,c) = std::distance(poly.vertices_begin(), he->vertex());
        F(f,c) = vertex_to_index[he->vertex()];
        c++;
      } while ( ++he != facet->facet_begin());
      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 subdiv_border( Polyhedron& P) {
    if ( P.size_of_facets() == 0)
        return;
    // We use that new halfedges are appended at the end.
    Edge_iterator last_e = P.edges_end();
    -- last_e;  // the last of the old edges
    Edge_iterator e = P.edges_begin(); // create trisected border edges
    do {
        if ( e->opposite()->is_border())
            trisect_border_halfedge( P, e->opposite());
        else if ( e->is_border())
            trisect_border_halfedge( P, e);
    } while ( e++ != last_e);
    e = P.edges_begin();               // smooth points on border edges
    std::vector<Point> pts;  // store new smoothed points temporarily
    do {
        if ( e->opposite()->is_border())
            smooth_border_vertices( e->opposite(), std::back_inserter(pts));
        else if ( e->is_border())
            smooth_border_vertices( e, std::back_inserter(pts));
    } while ( e++ != last_e);
    e = P.edges_begin(); // copy smoothed points back
    std::vector<Point>::iterator i = pts.begin();
    do {
        if ( e->opposite()->is_border()) {
            e->vertex()->point() = *i++;
            e->opposite()->vertex()->point() = *i++;
            e->opposite()->next()->vertex()->point() = *i++;
        } else if ( e->is_border()) {
            e->opposite()->vertex()->point() = *i++;
            e->vertex()->point() = *i++;
            e->next()->vertex()->point() = *i++;
        }
    } while ( e++ != last_e);
    CGAL_assertion( i == pts.end());
    CGAL_postcondition( P.is_valid());
}

// a helper method for running different iterators
void running_iterators( Polyhedron& P) {
    if ( P.size_of_facets() == 0)
        return;

    std::size_t nv = P.size_of_vertices();

    std::cout << "The number of vertices in the Polyhedron: " << nv << std::endl;
    std::cout << "The number of facets in the Polyhedron: " << P.size_of_facets() << std::endl;
    std::cout << "The number of half edges in the Polyhedron: " << P.size_of_halfedges() << std::endl;

    std::cout << std:: endl;

    Polyhedron::Vertex_iterator last_v = P.vertices_end();
    -- last_v;  // the last of the old vertices

    Polyhedron::Edge_iterator last_e = P.edges_end();
    -- last_e;  // the last of the old edges

    Polyhedron::Facet_iterator last_f = P.facets_end();
    -- last_f;  // the last of the old facets

    int k = 0;
    Polyhedron::Facet_iterator f = P.facets_begin();

    do {
    	std::cout << "Printing a facet index: " << k++ <<  std::endl;

    	f->halfedge();

    } while ( f++ != last_f);

    std::cout  << std::endl;

    // -------------------------------------------------
    // traverse the vertices
    // -------------------------------------------------

    std::cout << "Printing the vertex indices: " << std::endl;

     int n=0;
     for (Polyhedron::Vertex_iterator vi = P.vertices_begin(); vi != P.vertices_end(); ++vi)
     {
    	 Kernel::Point_3 p;
    	 p = vi->point();
    	 std::cout << "Vertex index: "  << n++ << std::endl;
    	 std::cout << "p.x() = "  << p.x() << std::endl;
    	 std::cout << "p.y() = "  << p.y() << std::endl;
    	 std::cout << "p.z() = "  << p.z() << std::endl;

     }

     std::cout  << std::endl;

     // -------------------------------------------------
     // traverse the edges
     // -------------------------------------------------

     std::cout << "Iterating over the edges.... " << std::endl;

     n=0;
     for (Polyhedron::Edge_iterator ei = P.edges_begin(); ei != P.edges_end(); ++ei)
     {
    	 ei->next();
    	 Kernel::Point_3 p;
    	 p =  ei->vertex()->point();
    	 std::cout << "For edge index: " << n++ << std::endl;
    	 std::cout << "p.x() = "  << p.x() << std::endl;
		 std::cout << "p.y() = "  << p.y() << std::endl;
		 std::cout << "p.z() = "  << p.z() << std::endl;

     }
     std::cout  << std::endl;

	 // -----------------------------------------------
	 // Do something else with the edge iterators
	 // -----------------------------------------------

    Polyhedron::Edge_iterator e = P.edges_begin();
    ++ last_e; // make it the past-the-end position again

    while ( e != last_e) {
    	Polyhedron::Halfedge_handle h = e;
        ++e;
    };

    CGAL_postcondition( P.is_valid());
}

int main()
#endif
{
    CGAL::Timer total_timer;
    total_timer.start();

    std::cerr << "PARAMETERIZATION" << std::endl;

    //***************************************
    // Read options on the command line
    //***************************************

    std::string type;               // default: Floater param
    std::string border;             // default: circular border param.
    std::string solver;             // default: OpenNL solver
    std::string input;              // required
    std::string output;             // default: out.eps
    try
    {
#ifdef CGAL_USE_BOOST_PROGRAM_OPTIONS
        po::options_description desc("Allowed options");
        desc.add_options()
            ("help,h", "prints this help message")
            ("type,t", po::value<std::string>(&type)->default_value("floater"),
            "parameterization method: floater, conformal, barycentric, authalic or lscm")
            ("border,b", po::value<std::string>(&border)->default_value("circle"),
            "border shape: circle, square or 2pts (lscm only)")
            ("solver,s", po::value<std::string>(&solver)->default_value("opennl"),
            "solver: opennl")
            ("input,i", po::value<std::string>(&input)->default_value(""),
            "input mesh (OFF)")
            ("output,o", po::value<std::string>(&output)->default_value("out.eps"),
            "output file (EPS or OBJ)")
            ;

        po::positional_options_description p;
        p.add("input", 1);
        p.add("output", 1);

        po::variables_map vm;
        po::store(po::command_line_parser(argc, argv).options(desc).positional(p).run(), vm);
        po::notify(vm);

        if (vm.count("help")) {
            std::cout << desc << "\n";
            return 1;
        }
#else
        std::cerr << "Command-line options require Boost.ProgramOptions" << std::endl;
        std::cerr << "Use hard-coded options" << std::endl;
        border = "square";
        type = "floater";
        solver = "opennl";
        input = "data/rotor.off";
        output = "rotor_floater_square_opennl_parameterized.obj";
#endif
    }
    catch(std::exception& e) {
      std::cerr << "error: " << e.what() << "\n";
      return 1;
    }
    catch(...) {
      std::cerr << "Exception of unknown type!\n";
      throw;
    }

    //***************************************
    // Read the mesh
    //***************************************

    CGAL::Timer task_timer;
    task_timer.start();

    // Read the mesh
    std::ifstream stream(input.c_str());
    Polyhedron mesh;
    stream >> mesh;
    if(!stream || !mesh.is_valid() || mesh.empty())
    {
        std::cerr << "Error: cannot read OFF file " << input << std::endl;
        return EXIT_FAILURE;
    }

    std::cerr << "Read file " << input << ": "
              << task_timer.time() << " seconds "
              << "(" << mesh.size_of_facets() << " facets, "
              << mesh.size_of_vertices() << " vertices)" << std::endl;
    task_timer.reset();

    //***************************************
    // Create mesh adaptor
    //***************************************

    // The Surface_mesh_parameterization package needs an adaptor to handle Polyhedron_ex meshes
    Parameterization_polyhedron_adaptor mesh_adaptor(mesh);

    // The parameterization methods support only meshes that
    // are topological disks => we need to compute a cutting path
    // that makes the mesh a "virtual" topological disk
    //
//.........这里部分代码省略.........

/*
 * mexFunction(): entry point for the mex function
 */
void mexFunction(int nlhs, mxArray *plhs[], 
		 int nrhs, const mxArray *prhs[]) {

  // interface to deal with input arguments from Matlab
  enum InputIndexType {IN_TRI, IN_X, IN_METHOD, IN_ITER, InputIndexType_MAX};
  MatlabImportFilter::Pointer matlabImport = MatlabImportFilter::New();
  matlabImport->ConnectToMatlabFunctionInput(nrhs, prhs);

  // check that we have all input arguments
  matlabImport->CheckNumberOfArguments(2, InputIndexType_MAX);

  // register the inputs for this function at the import filter
  MatlabInputPointer inTRI =        matlabImport->RegisterInput(IN_TRI, "TRI");
  MatlabInputPointer inX =          matlabImport->RegisterInput(IN_X, "X");
  MatlabInputPointer inMETHOD =     matlabImport->RegisterInput(IN_METHOD, "METHOD");
  MatlabInputPointer inITER =       matlabImport->RegisterInput(IN_ITER, "ITER");

  // interface to deal with outputs to Matlab
  enum OutputIndexType {OUT_TRI, OUT_N, OutputIndexType_MAX};
  MatlabExportFilter::Pointer matlabExport = MatlabExportFilter::New();
  matlabExport->ConnectToMatlabFunctionOutput(nlhs, plhs);

  // check number of outputs the user is asking for
  matlabExport->CheckNumberOfArguments(0, OutputIndexType_MAX);

  // register the outputs for this function at the export filter
  typedef MatlabExportFilter::MatlabOutputPointer MatlabOutputPointer;
  MatlabOutputPointer outTRI = matlabExport->RegisterOutput(OUT_TRI, "TRI");
  MatlabOutputPointer outN = matlabExport->RegisterOutput(OUT_N, "N");

  // if any of the inputs is empty, the output is empty too
  if (mxIsEmpty(prhs[IN_TRI]) || mxIsEmpty(prhs[IN_X])) {
    matlabExport->CopyEmptyArrayToMatlab(outTRI);
    matlabExport->CopyEmptyArrayToMatlab(outN);
    return;
  }

  // polyhedron to contain the input mesh
  Polyhedron mesh;
  PolyhedronBuilder<Polyhedron> builder(matlabImport, inTRI, inX);
  mesh.delegate(builder);

  // get size of input matrix with the points
  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
//.........这里部分代码省略.........