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

polyhedron cgal_to_polyhedron( const Nef_polyhedron &NP ){
    Polyhedron P;
    polyhedron ret;
    
    if( NP.is_simple() ){
        NP.convert_to_polyhedron(P);
        std::vector<double> coords;
        std::vector<int> tris;
        int next_id = 0;
        std::map< Polyhedron::Vertex*, int > vid;
        for( Polyhedron::Vertex_iterator iter=P.vertices_begin(); iter!=P.vertices_end(); iter++ ){
            coords.push_back( CGAL::to_double( (*iter).point().x() ) );
            coords.push_back( CGAL::to_double( (*iter).point().y() ) );
            coords.push_back( CGAL::to_double( (*iter).point().z() ) );
            vid[ &(*iter) ] = next_id++;
        }
        
        for( Polyhedron::Facet_iterator iter=P.facets_begin(); iter!=P.facets_end(); iter++ ){
            Polyhedron::Halfedge_around_facet_circulator j = iter->facet_begin();
            tris.push_back( CGAL::circulator_size(j) );
            do {
                tris.push_back( std::distance(P.vertices_begin(), j->vertex()) );
            } while ( ++j != iter->facet_begin());
        }
        
        ret.initialize_load_from_mesh( coords, tris );
    } else {
        std::cout << "resulting polyhedron is not simple!" << std::endl;
    }
    return ret;
}

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

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

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

void transformMesh( Polyhedron & poly, Transformation3 & map )
{
    for ( Vertex_iterator vi = poly.vertices_begin(); vi != poly.vertices_end(); ++vi ) {
	Point3 store = vi->point().transform( map );
	vi->point() = store;
    }
}

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

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

//------------------------------------------------------------------------------
//	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 modify_vertex_position()
{
	Polyhedron P;
    Halfedge_handle h = P.make_tetrahedron();
    if ( P.is_tetrahedron(h))
	{
		int i(0);
		for(Vertex_iterator vi = P.vertices_begin(); vi != P.vertices_end(); ++vi,++i)
		{
		  std::cout << "before changing vertex " << i << ": " << vi->point().x() << vi->point().y() << vi->point().z() << endl;
		  Point_3 pt(1, 0, 0);
		  vi->point() = pt;
		  std::cout << "after changing vertex " << i << ": " << vi->point().x() << vi->point().y() << vi->point().z() << endl;
		}
	}
}

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

scalarField calcCurvature(const triSurface& surf)
{
    scalarField k(surf.points().size(), 0);

    Polyhedron P;

    buildCGALPolyhedron convert(surf);
    P.delegate(convert);

    // Info<< "Created CGAL Polyhedron with " << label(P.size_of_vertices())
    //     << " vertices and " << label(P.size_of_facets())
    //     << " facets. " << endl;

    // The rest of this function adapted from
    //     CGAL-3.7/examples/Jet_fitting_3/Mesh_estimation.cpp

     //Vertex property map, with std::map
    typedef std::map<Vertex*, int> Vertex2int_map_type;
    typedef boost::associative_property_map< Vertex2int_map_type >
        Vertex_PM_type;
    typedef T_PolyhedralSurf_rings<Polyhedron, Vertex_PM_type > Poly_rings;

    typedef CGAL::Monge_via_jet_fitting<Kernel>         Monge_via_jet_fitting;
    typedef Monge_via_jet_fitting::Monge_form           Monge_form;

    std::vector<Point_3> in_points;  //container for data points

    // default parameter values and global variables
    unsigned int d_fitting = 2;
    unsigned int d_monge = 2;
    unsigned int min_nb_points = (d_fitting + 1)*(d_fitting + 2)/2;

    //initialize the tag of all vertices to -1
    Vertex_iterator vitb = P.vertices_begin();
    Vertex_iterator vite = P.vertices_end();

    Vertex2int_map_type vertex2props;
    Vertex_PM_type vpm(vertex2props);

    CGAL_For_all(vitb, vite)
    {
        put(vpm, &(*vitb), -1);
    }

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

//**********************************************************************************
//returns min coordinates
Vector3r MinCoord(const shared_ptr<Shape>& cm1,const State& state1){
	const Se3r& se3=state1.se3; 
	Polyhedra* A = static_cast<Polyhedra*>(cm1.get());

	//move and rotate CGAL structure Polyhedron
	Matrix3r rot_mat = (se3.orientation).toRotationMatrix();
	Vector3r trans_vec = se3.position;
	Transformation t_rot_trans(rot_mat(0,0),rot_mat(0,1),rot_mat(0,2), trans_vec[0],rot_mat(1,0),rot_mat(1,1),rot_mat(1,2),trans_vec[1],rot_mat(2,0),rot_mat(2,1),rot_mat(2,2),trans_vec[2],1.);
	Polyhedron PA = A->GetPolyhedron();
	std::transform( PA.points_begin(), PA.points_end(), PA.points_begin(), t_rot_trans);
	
	Vector3r minccord = trans_vec;
	for(Polyhedron::Vertex_iterator vi = PA.vertices_begin(); vi != PA.vertices_end(); ++vi){	
		if (vi->point()[0]<minccord[0]) minccord[0]=vi->point()[0];
		if (vi->point()[1]<minccord[1]) minccord[1]=vi->point()[1];
		if (vi->point()[2]<minccord[2]) minccord[2]=vi->point()[2];
	}
	
	return minccord;
}

ofMesh & ofxCGALSkinSurface::makeSkinSurfaceMesh() {
    std::list<Weighted_point> l;
    FT shrinkfactor = shrinkFactor;
    
    for(int i=0; i<points.size(); i++) {
        float px = points[i].point.x;
        float py = points[i].point.y;
        float pz = points[i].point.z;
        float radius = points[i].radius;
        
        l.push_front(Weighted_point(Bare_point(px, py, pz), radius * radius));
    }
    
    Skin_surface_3 skin_surface(l.begin(), l.end(), shrinkfactor);
    
    Polyhedron p;
    CGAL::mesh_skin_surface_3(skin_surface, p);
    
    if(bSubdiv == true) {
        CGAL::subdivide_skin_surface_mesh_3(skin_surface, p);
    }
    
    mesh.clear();
    
    map<Point_3, int> point_indices;
    int count = 0;
    
    for (auto it=p.vertices_begin(); it!=p.vertices_end(); ++it) {
        auto& p = it->point();
        mesh.addVertex(ofVec3f(p.x(), p.y(), p.z()));
        point_indices[p] = count++;
    }
    for (auto it=p.facets_begin(); it!=p.facets_end(); ++it) {
        mesh.addIndex(point_indices[it->halfedge()->vertex()->point()]);
        mesh.addIndex(point_indices[it->halfedge()->next()->vertex()->point()]);
        mesh.addIndex(point_indices[it->halfedge()->prev()->vertex()->point()]);
    }
}

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