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

void mergeColinear(Polyhedron& p) {
	int vertBefore =  p.size_of_vertices();
	bool colinearFound = true;
	while (colinearFound) {
		colinearFound = false; // Set temporarily to false, if merge degments then set true

		int percCount = 1;
		for (Polyhedron::Halfedge_iterator hit = p.halfedges_begin(); hit != p.halfedges_end(); ++hit,++percCount){
			if (CGAL::circulator_size(hit->vertex_begin()) == 1) {
				std::cerr << "WARNING: Loose edge, but how??"<<std::endl;
				while (CGAL::circulator_size(hit->vertex_begin()) == 1)
					hit = p.join_vertex(hit->opposite());
				break;
			}

			if ((CGAL::circulator_size(hit->vertex_begin()) == 2) &&			// if only two he connected to vertex
				(hit->facet_degree()>3 && hit->opposite()->facet_degree()>3)) { // if faces are not triangles // prob faster

				Vector_3 cur(hit->prev()->vertex()->point(),hit->vertex()->point());
				Vector_3 nex(hit->vertex()->point(),hit->next()->vertex()->point());
				if ( is_colinear(cur,nex)) {									// check if colinear
					std::cout << "\rVertices before/after: "<<vertBefore<<" -> "<< p.size_of_vertices()<<". ("<<100*percCount/p.size_of_halfedges()<<"%)";
					p.join_vertex(hit->opposite());								// move cur to prev point
					colinearFound = true;
					break;
	}	}	}	}
	std::cout << "\rVertices before/after: "<<vertBefore<<" -> "<< p.size_of_vertices()<<". (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 smoothMesh( Polyhedron & poly, unsigned int nTimes )
{
    int nV = poly.size_of_vertices();
    const double lambda = SMOOTHING_LAMBDA;
    const double mu	= SMOOTHING_MU;

    vector< Point3 > shrink  ( nV );
    vector< Point3 > expand  ( nV );

    for ( unsigned int k = 0; k < nTimes; ++k ) {
	// copy the vertex coordinates
	for ( Vertex_iterator vi = poly.vertices_begin(); vi != poly.vertices_end(); ++vi ) {
	    shrink  [ vi->id() ] = vi->point();
	}
	// shrinking stage
	for ( Vertex_iterator vi = poly.vertices_begin(); vi != poly.vertices_end(); ++vi ) {
	    moveVertex( vi, shrink[ vi->id() ], lambda );
	}
	// copy back the vertex coordinates
	for ( Vertex_iterator vi = poly.vertices_begin(); vi != poly.vertices_end(); ++vi ) {
	    vi->point()		= shrink[ vi->id() ];
	    expand[ vi->id() ]	= shrink[ vi->id() ];
	}
	// expanding stage
	for ( Vertex_iterator vi = poly.vertices_begin(); vi != poly.vertices_end(); ++vi ) {
	    moveVertex( vi, expand[ vi->id() ], mu );
	}
	// copy back the vertex coordinates
	for ( Vertex_iterator vi = poly.vertices_begin(); vi != poly.vertices_end(); ++vi ) {
	    vi->point()		= expand[ vi->id() ];
	}
    }
}

 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;
 }

int main()
{
  typedef CGAL::Polyhedron_3<Kernel> Polyhedron;
  typedef CGAL::Polyhedron_corefinement<Polyhedron> Corefinement;

  std::stringstream fpa(cube_a);
  std::stringstream fpb(cube_b);
  Polyhedron pa;
  Polyhedron pb;

  fpa >> pa;
  fpb >> pb;

  assert( pa.size_of_vertices() == 8);
  assert( pb.size_of_vertices() == 8);
  
  {
  Corefinement coref;
  std::list<std::vector<Kernel::Point_3> > polylines;
  std::vector<std::pair<Polyhedron*, int> > result;
  coref( pa, pb, std::back_inserter(polylines), std::back_inserter(result), Corefinement::Intersection_tag );
  
  assert( polylines.size() == 1 );  
  assert( polylines.begin()->size() == 1 );  
  assert( result.size() == 0 );
  }

  pb.clear();
  std::stringstream fpc(inside_a);
  fpc >> pb;
  assert( pb.size_of_vertices() == 4);

  {
  Corefinement coref;
  std::list<std::vector<Kernel::Point_3> > polylines;
  std::vector<std::pair<Polyhedron*, int> > result;
  coref( pa, pb, std::back_inserter(polylines), std::back_inserter(result), Corefinement::Intersection_tag );
  
  assert( polylines.size() == 4 ); 
  assert( polylines.begin()->size() == 1 );  
  assert( result.size() == 1 );
  assert( result[0].first->size_of_vertices() == 4);
  }
}

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"));
}

//------------------------------------------------------------------------------
//	Normalize the 3D triangulated mesh with the display window
//------------------------------------------------------------------------------
void normalizeMesh( Polyhedron & poly, vector< Segment3 > & bone )
{
    Vector3 sum, ave;

    for ( Vertex_iterator vi = poly.vertices_begin(); vi != poly.vertices_end(); ++vi ) {
	sum = sum + ( vi->point() - CGAL::ORIGIN );
    }
    ave = sum / ( double )poly.size_of_vertices();

    for ( Vertex_iterator vi = poly.vertices_begin(); vi != poly.vertices_end(); ++vi ) {
	vi->point() = vi->point() - ave;
    }
    
    cerr << " ave = " << ave << endl;
    Transformation3 translate( CGAL::TRANSLATION, -ave );

    // fabs:absolute values-no negative values
    double sideMax = 0.0;
    for ( Vertex_iterator vi = poly.vertices_begin(); vi != poly.vertices_end(); ++vi ) {
	if ( fabs( vi->point().x() ) > sideMax ) sideMax = fabs( vi->point().x() );
	if ( fabs( vi->point().y() ) > sideMax ) sideMax = fabs( vi->point().y() );
	if ( fabs( vi->point().z() ) > sideMax ) sideMax = fabs( vi->point().z() );
    }
    
	// sideMax: the largest number
    for ( Vertex_iterator vi = poly.vertices_begin(); vi != poly.vertices_end(); ++vi ) {
	vi->point() = CGAL::ORIGIN + ( vi->point() - CGAL::ORIGIN ) / sideMax;
    }

    Transformation3 scale( CGAL::SCALING, 1.0/sideMax );

    Transformation3 composite = scale * translate;
	    
    for ( unsigned int k = 0; k < bone.size(); ++k ) {
	bone[ k ] = bone[ k ].transform( composite );
    }
}

void alignMesh( Polyhedron & poly )
{
    int num = poly.size_of_vertices();

    // initialization here is very important!!
    Point3 ave( 0.0, 0.0, 0.0 );
    for ( Vertex_iterator vi = poly.vertices_begin(); vi != poly.vertices_end(); ++vi ) {
	ave = ave + ( vi->point() - CGAL::ORIGIN );
    }
    ave = CGAL::ORIGIN + ( ave - CGAL::ORIGIN )/( double )num;

    unsigned int dim	= 3;
    double * data	= new double [ num*dim ];

    int nPoints = 0;
    for ( Vertex_iterator vi = poly.vertices_begin(); vi != poly.vertices_end(); ++vi ) {
	data[ nPoints * dim + 0 ] = vi->point().x() - ave.x();
	data[ nPoints * dim + 1 ] = vi->point().y() - ave.y();
	data[ nPoints * dim + 2 ] = vi->point().z() - ave.z();
	nPoints++;
    }

    assert( nPoints == ( int )num );

    /*****************************************
      analyze the engenstructure of X^T X
     *****************************************/
    /* define the matrix X */
    gsl_matrix_view X = gsl_matrix_view_array( data, num, dim );

    /* memory reallocation */
    // proj = ( double * )realloc( proj, sizeof( double ) * PDIM * num );

    /* calculate the covariance matrix B */
    gsl_matrix * B              = gsl_matrix_alloc( dim, dim );
    gsl_blas_dgemm( CblasTrans, CblasNoTrans, 1.0,
                    &(X.matrix),
                    &(X.matrix), 0.0,
                    B );
       
    /* divided by the number of samples */
    gsl_matrix_scale( B, 1.0/(double)num );

    gsl_vector * eVal       = gsl_vector_alloc( dim );
    gsl_matrix * eVec       = gsl_matrix_alloc( dim, dim );
    gsl_eigen_symmv_workspace * w
                                = gsl_eigen_symmv_alloc( dim );

    // eigenanalysis of the matrix B
    gsl_eigen_symmv( B, eVal, eVec, w );
    // release the memory of w
    gsl_eigen_symmv_free( w );
    // sort eigenvalues in a descending order
    gsl_eigen_symmv_sort( eVal, eVec, GSL_EIGEN_SORT_VAL_DESC );

    // #ifdef MYDEBUG
    for ( unsigned int i = 0; i < dim; ++i ) {
	cerr << "Eigenvalue No. " << i << " = " << gsl_vector_get( eVal, i ) << endl;
	cerr << "Eigenvector No. " << i << endl;
	double length = 0.0;
	for ( unsigned int j = 0; j < dim; ++j ) {
	    cerr << gsl_matrix_get( eVec, i, j ) << " ";
	    length += gsl_matrix_get( eVec, i, j )*gsl_matrix_get( eVec, i, j );
	}
	cerr << " length = " << length << endl;
    }
    // #endif	// MYDEBUG

    // 0 1 2, 0 2 1, 
    

    Transformation3 map;
    map = 
	Transformation3( gsl_matrix_get(eVec,1,0), gsl_matrix_get(eVec,0,0), gsl_matrix_get(eVec,2,0),
			 gsl_matrix_get(eVec,1,1), gsl_matrix_get(eVec,0,1), gsl_matrix_get(eVec,2,1),
			 gsl_matrix_get(eVec,1,2), gsl_matrix_get(eVec,0,2), gsl_matrix_get(eVec,2,2) );
    if ( map.is_odd() ) {
	cerr << " Transformation matrix reflected" << endl;
	map = 
	    Transformation3( gsl_matrix_get(eVec,1,0), gsl_matrix_get(eVec,0,0), -gsl_matrix_get(eVec,2,0),
			     gsl_matrix_get(eVec,1,1), gsl_matrix_get(eVec,0,1), -gsl_matrix_get(eVec,2,1),
			     gsl_matrix_get(eVec,1,2), gsl_matrix_get(eVec,0,2), -gsl_matrix_get(eVec,2,2) );
    }

    for ( unsigned int i = 0; i < dim; ++i ) {
	cerr << "| ";
	for ( unsigned int j = 0; j < dim; ++j ) {
	    cerr << map.cartesian( i, j ) << " ";
	}
	cerr << "|" << endl;
    }

    transformMesh( poly, map );

    return;
}

//------------------------------------------------------------------------------
//	Align the mesh
//------------------------------------------------------------------------------
Vector3 principalAxis( Polyhedron & poly )
{
    int num = poly.size_of_vertices();

    // initialization here is very important!!
    Point3 ave( 0.0, 0.0, 0.0 );
    for ( Vertex_iterator vi = poly.vertices_begin(); vi != poly.vertices_end(); ++vi ) {
	ave = ave + ( vi->point() - CGAL::ORIGIN );
    }
    ave = CGAL::ORIGIN + ( ave - CGAL::ORIGIN )/( double )num;

    unsigned int dim	= 3;
    double * data	= new double [ num*dim ];

    int nPoints = 0;
    for ( Vertex_iterator vi = poly.vertices_begin(); vi != poly.vertices_end(); ++vi ) {
	data[ nPoints * dim + 0 ] = vi->point().x() - ave.x();
	data[ nPoints * dim + 1 ] = vi->point().y() - ave.y();
	data[ nPoints * dim + 2 ] = vi->point().z() - ave.z();
	nPoints++;
    }

    assert( nPoints == ( int )num );

    /*****************************************
      analyze the engenstructure of X^T X
     *****************************************/
    /* define the matrix X */
    gsl_matrix_view X = gsl_matrix_view_array( data, num, dim );

    /* memory reallocation */
    // proj = ( double * )realloc( proj, sizeof( double ) * PDIM * num );

    /* calculate the covariance matrix B */
    gsl_matrix * B              = gsl_matrix_alloc( dim, dim );
    gsl_blas_dgemm( CblasTrans, CblasNoTrans, 1.0,
                    &(X.matrix),
                    &(X.matrix), 0.0,
                    B );
       
    /* divided by the number of samples */
    gsl_matrix_scale( B, 1.0/(double)num );

    gsl_vector * eVal       = gsl_vector_alloc( dim );
    gsl_matrix * eVec       = gsl_matrix_alloc( dim, dim );
    gsl_eigen_symmv_workspace * w
                                = gsl_eigen_symmv_alloc( dim );

    // eigenanalysis of the matrix B
    gsl_eigen_symmv( B, eVal, eVec, w );
    // release the memory of w
    gsl_eigen_symmv_free( w );
    // sort eigenvalues in a descending order
    gsl_eigen_symmv_sort( eVal, eVec, GSL_EIGEN_SORT_VAL_DESC );

#ifdef MYDEBUG
    for ( unsigned int i = 0; i < dim; ++i ) {
	cerr << "Eigenvalue No. " << i << " = " << gsl_vector_get( eVal, i ) << endl;
	cerr << "Eigenvector No. " << i << endl;
	for ( unsigned int j = 0; j < dim; ++j ) {
	    length += gsl_matrix_get( eVec, i, j )*gsl_matrix_get( eVec, i, j );
	}
	cerr << " length = " << length << endl;
    }
#endif	// MYDEBUG

    Vector3 ref( gsl_matrix_get( eVec, 0, 0 ),
		 gsl_matrix_get( eVec, 0, 1 ),
		 gsl_matrix_get( eVec, 0, 2 ) );
    return ref;

#ifdef DEBUG
    gsl_vector_view eachVec = gsl_matrix_column( eigenVec, 0 );
    double cosRot = gsl_matrix_get( eigenVec, 0, 1 );
    double sinRot = gsl_matrix_get( eigenVec, 1, 1 );
#ifdef DEBUG
    cerr << " 2nd axis : " << cosRot << " , " << sinRot << endl;
#endif	// DEBUG

    Transformation2 rotate( CGAL::ROTATION, -sinRot, cosRot );
    for ( unsigned int i = 0; i < subpatch.size(); ++i ) {
	subpatch[ i ]->triangle() = subpatch[ i ]->triangle().transform( rotate );
    }
#endif	// DEBUG
}

int main(int argc, char* argv[])
{
  std::cerr.precision(17);

  std::ifstream points_file(argv[2]);

  std::vector<Point> points;
  std::copy(std::istream_iterator<Point>(points_file),
            std::istream_iterator<Point>(),
            std::back_inserter(points));


  int gridsize = 10;
  if(argc>3){
    gridsize = boost::lexical_cast<int>(argv[3]);
  }
  std::cerr << "gridsize = " << gridsize << std::endl;
  
  int nb_points=points.size();
  
  std::vector<bool> ray_res(nb_points);
  std::vector<bool> grid_res(nb_points);
  
  //using ray
  {
    Polyhedron polyhedron;
    std::ifstream polyhedron_file(argv[1]);
    polyhedron_file >> polyhedron;
    std::cerr << "|V| = " << polyhedron.size_of_vertices() << std::endl;
    
    CGAL::Timer timer;
    timer.start();
    CGAL::Point_inside_polyhedron_3<Polyhedron,K> inside_with_ray(polyhedron,0);
    timer.stop();
    std::cerr <<"Using ray"<< std::endl;
    std::cerr << "  Preprocessing took " << timer.time() << " sec." << std::endl;
    timer.reset();  
    int n_inside = 0;
    
    timer.start();
    for(int k=0;k<nb_points;++k){
      ray_res[k]=inside_with_ray(points[k]);
      if(ray_res[k]){
        ++n_inside;
      }
    }
    timer.stop();
    std::cerr << "  " << n_inside << " points inside " << std::endl;
    std::cerr << "  " << points.size() - n_inside << " points outside "  << std::endl;
    std::cerr << " Queries took " << timer.time() << " sec." << std::endl; 

  }
  
  //using grid
  {
    Polyhedron polyhedron;
    std::ifstream polyhedron_file(argv[1]);
    polyhedron_file >> polyhedron;
    std::cerr << "|V| = " << polyhedron.size_of_vertices() << std::endl;
    
    CGAL::Timer timer;
    timer.start();
    CGAL::Point_inside_polyhedron_3<Polyhedron,K> inside_with_grid(polyhedron, gridsize);
    timer.stop();
    std::cerr <<"Using grid"<< std::endl;
    std::cerr << "  Preprocessing took " << timer.time() << " sec." << std::endl;
    timer.reset();

    if(argc>5){
      random_points(argv[4],inside_with_grid.bbox() ,boost::lexical_cast<int>(argv[5]) );
    }
    
    int n_inside = 0;
    timer.start();
    for(int k=0;k<nb_points;++k){
      grid_res[k]=inside_with_grid(points[k]);
      if(grid_res[k]){
        ++n_inside;
      }
    }
    timer.stop();
    std::cerr << "  " << n_inside << " points inside " << std::endl;
    std::cerr << "  " << points.size() - n_inside << " points outside "  << std::endl;
    std::cerr << "  Queries took " << timer.time() << " sec." << std::endl; 
  }

  for(int k=0;k<nb_points;++k){
    if(ray_res[k]!=grid_res[k]){
      std::cerr << "WARNING: Result is different for point " << k << std::endl;
    }
  }
  
  //using original code
  {
    Polyhedron polyhedron;
    std::ifstream polyhedron_file(argv[1]);
    polyhedron_file >> polyhedron;
    std::cerr << "|V| = " << polyhedron.size_of_vertices() << std::endl;

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

// 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());
}

/*
 * 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
//.........这里部分代码省略.........