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

//**********************************************************************************
//test of polyhedron intersection callable from python shell
bool do_Polyhedras_Intersect(const shared_ptr<Shape>& cm1,const shared_ptr<Shape>& cm2,const State& state1,const State& state2){

	const Se3r& se31=state1.se3; 
	const Se3r& se32=state2.se3;
	Polyhedra* A = static_cast<Polyhedra*>(cm1.get());		
	Polyhedra* B = static_cast<Polyhedra*>(cm2.get());

	//move and rotate 1st the CGAL structure Polyhedron
	Matrix3r rot_mat = (se31.orientation).toRotationMatrix();
	Vector3r trans_vec = se31.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);

	//move and rotate 2st the CGAL structure Polyhedron
	rot_mat = (se32.orientation).toRotationMatrix();
	trans_vec = se32.position;
	t_rot_trans = Transformation(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 PB = B->GetPolyhedron();
	std::transform( PB.points_begin(), PB.points_end(), PB.points_begin(), t_rot_trans);

	//calculate plane equations
	std::transform( PA.facets_begin(), PA.facets_end(), PA.planes_begin(),Plane_equation());
	std::transform( PB.facets_begin(), PB.facets_end(), PB.planes_begin(),Plane_equation());


	//call test
	return do_intersect(PA,PB);
}

// Add semantics to Interior room polyhedra
void set_semantic_InteriorLoD4(Polyhedron& polyhe) {
	std::transform( polyhe.facets_begin(), polyhe.facets_end(),polyhe.planes_begin(),Plane_Newel_equation());		
	for (Polyhedron::Facet_iterator fIt = polyhe.facets_begin(); fIt != polyhe.facets_end(); ++fIt) {	// Iterate over faces
		Kernel::FT z = fIt->plane().orthogonal_vector().z();
		if		(z <= -HORIZONTAL_ANGLE_RANGE)	fIt->semanticBLA = "FloorSurface";
		else if (z >=  HORIZONTAL_ANGLE_RANGE)	fIt->semanticBLA = "CeilingSurface";
		else									fIt->semanticBLA = "InteriorWallSurface";
}	}	

bool is_weakly_convex(Polyhedron const& p) {
  for (typename Polyhedron::Edge_const_iterator i = p.edges_begin(); i != p.edges_end(); ++i) {
    typename Polyhedron::Plane_3 p(i->opposite()->vertex()->point(), i->vertex()->point(), i->next()->vertex()->point());
    if (p.has_on_positive_side(i->opposite()->next()->vertex()->point()) &&
        CGAL::squared_distance(p, i->opposite()->next()->vertex()->point()) > 1e-8) {
      return false;
    }
  }
  // Also make sure that there is only one shell:
  boost::unordered_set<typename Polyhedron::Facet_const_handle, typename CGAL::Handle_hash_function> visited;
  // c++11
  // visited.reserve(p.size_of_facets());
  
  std::queue<typename Polyhedron::Facet_const_handle> to_explore;
  to_explore.push(p.facets_begin()); // One arbitrary facet
  visited.insert(to_explore.front());
  
  while (!to_explore.empty()) {
    typename Polyhedron::Facet_const_handle f = to_explore.front();
    to_explore.pop();
    typename Polyhedron::Facet::Halfedge_around_facet_const_circulator he, end;
    end = he = f->facet_begin();
    CGAL_For_all(he,end) {
      typename Polyhedron::Facet_const_handle o = he->opposite()->facet();
      
      if (!visited.count(o)) {
        visited.insert(o);
        to_explore.push(o);
      }
    }
  }

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

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

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

static int number_of_degenerate_facets(Polyhedron& p, const double threshold)
{
  int count = 0;
  for (typename Polyhedron::Facet_iterator facet = p.facets_begin();
       facet != p.facets_end(); facet++)
  {
    dolfin_assert(facet->is_triangle());
    if ( facet_is_degenerate<Polyhedron>(facet, threshold) )
      count++;
  }
  return count;
}

int main()
{
  std::list<Weighted_point> l;
  FT shrinkfactor = 0.5;

  l.push_front(Weighted_point(Bare_point(0, 0, 0), 1));
  l.push_front(Weighted_point(Bare_point(0, 1, 0), 2));
  l.push_front(Weighted_point(Bare_point(0, 0, 2), 1));

  Skin_surface_3 skin_surface(l.begin(), l.end(), shrinkfactor);

  Polyhedron p;
  CGAL::mesh_skin_surface_3(skin_surface, p);

  for (Polyhedron::Facet_iterator fit = p.facets_begin(); fit != p.facets_end(); ++fit) {

    // Retrieve the generating simplex from the regular triangulation
    const Skin_surface_3::Simplex &s = fit->containing_simplex();

    // Get the weighted points from the simplex
    Skin_surface_3::Weighted_point wp[4];
    switch (s.dimension()) {
    case 0: {
      Skin_surface_3::Vertex_handle vh = s;
      wp[0] = vh->point();
      break;
    }
    case 1: {
      Skin_surface_3::Edge e = s;
      wp[0] = e.first->vertex(e.second)->point();
      wp[1] = e.first->vertex(e.third)->point();
      break;
    }
    case 2: {
      Skin_surface_3::Facet f = s;
      wp[0] = f.first->vertex((f.second + 1) & 3)->point();
      wp[1] = f.first->vertex((f.second + 2) & 3)->point();
      wp[2] = f.first->vertex((f.second + 3) & 3)->point();
      break;
    }
    case 3: {
      Skin_surface_3::Cell_handle ch = s;
      wp[0] = ch->vertex(0)->point();
      wp[1] = ch->vertex(1)->point();
      wp[2] = ch->vertex(2)->point();
      wp[3] = ch->vertex(3)->point();
      break;
    }
    }
  }

  return 0;
}

int main(int argc, char** argv)
{
  typedef boost::property_map< 
    Polyhedron,
    CGAL::face_index_t 
    >::const_type Face_index_map;

  std::ifstream in((argc>1)?argv[1]:"cube.off");
  Polyhedron P;
  in >> P ;
  
  // initialize facet indices
  std::size_t i = 0;
  for(Polyhedron::Facet_iterator it = P.facets_begin(); it != P.facets_end(); ++it, ++i) 
  {
    it->id() = i;
  }

  // Ad hoc property_map to store normals. Face_index_map is used to
  // map face_descriptors to a contiguous range of indices. See
  // http://www.boost.org/libs/property_map/doc/vector_property_map.html
  // for details.
  boost::vector_property_map<Vector, Face_index_map> 
    normals(get(CGAL::face_index, P));

  calculate_face_normals(
    P // Graph
    , get(CGAL::vertex_point, P) // map from vertex_descriptor to point
    , normals // map from face_descriptor to Vector_3
    );

  std::cout << "Normals" << std::endl;
  for(Polyhedron::Facet_iterator it = P.facets_begin(); it != P.facets_end(); ++it) { 
    // Facet_iterator is a face_descriptor, so we can use it as the
    // key here
    std::cout << normals[it] << std::endl;
  }

  return 0;
}

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

// Remove small triangles with different semantics
void remove_singularSemantics(Polyhedron& exteriorPolyhe) {
	double SMALL_TRIANGLE_LIMIT = 0.1;
	Polyhedron::Facet_iterator exfIt;						// Iterate over exterior faces

	while (true) {
		bool newSemFound = false;
		for (exfIt = exteriorPolyhe.facets_begin(); exfIt != exteriorPolyhe.facets_end(); ++exfIt) {
			if (exfIt->area < SMALL_TRIANGLE_LIMIT) {
				std::map<std::string,double> semCountList;
				std::map<std::string,double>::iterator semCount;
				Polyhedron::Facet::Halfedge_around_facet_circulator itHDS = exfIt->facet_begin();

				double maxArea = 0;
				std::string newSem;
						
					do {
						if (is_coplanar(itHDS,false)) {
							std::string otherSem = itHDS->opposite()->facet()->semanticBLA;

							// Compute sum of area

							if (is_coplanar(itHDS->opposite()->next(),true) ||
								is_coplanar(itHDS->opposite()->prev(),true)) {
								semCount = semCountList.find(otherSem);
								if (semCount == semCountList.end())
									semCountList[otherSem] = 1; // Might not be needed
								else
									++(semCountList[otherSem]);
							}
						}
					} while (++itHDS != exfIt->facet_begin());
				
					for (semCount = semCountList.begin();semCount!=semCountList.end();++semCount) {
						if (semCount->second > maxArea) {
							maxArea = semCount->second;
							newSem = semCount->first;
						} else if (semCount->second == maxArea &&
							((exfIt->semanticBLA=="Anything"&&semCount->first!="Anything")||semCount->first==exfIt->semanticBLA))
							newSem = semCount->first;
					}
					
				if (maxArea > 0 && newSem!=exfIt->semanticBLA) {
					exfIt->semanticBLA = newSem;
					newSemFound = true;
					break;										// Dafuq? is this correct?
				}
			}
		}
		if (!newSemFound) break;								// Really?
	}
}

 double lloyd_step (std::vector<Point> &points)
 {
   std::vector<Plane> planes;
   for (size_t i = 0; i < points.size(); ++i)
     planes.push_back(tangent_plane(points[i]));

   Polyhedron P;

   CGAL::internal::halfspaces_intersection(planes.begin(),
					   planes.end(), P, K());

   size_t N = points.size();
   points.clear();
   double totarea = 0.0;
   double vararea = 0.0;

   for (Polyhedron::Facet_iterator it = P.facets_begin();
	it != P.facets_end(); ++it)
     {
       Polyhedron::Halfedge_around_facet_circulator
	 h0 = it->facet_begin(), hf = h0--, hs = hf;
       hs ++;

       double farea = 0.0;
       Vector fcentroid (CGAL::NULL_VECTOR);
       while(1)
	 {
	   const Point &a  = h0->vertex()->point(),
	     &b = hf->vertex()->point(),
	     &c = hs->vertex()->point();
	   Vector v = 
	     CGAL::cross_product(b -a, c - a);
	   double tarea = .5 * sqrt(v.squared_length());
	   farea += tarea;
	   fcentroid = fcentroid +  (tarea/3.0) * ((a - CGAL::ORIGIN) + 
						   (b - CGAL::ORIGIN) + 
						   (c - CGAL::ORIGIN));
	   if (hs == h0)
	     break;
	   ++hs; ++hf;
	 }
       totarea += farea;
       vararea += pow(4.0 * M_PI / N - farea, 2.0);
       fcentroid = fcentroid/farea;
       points.push_back(project_back(CGAL::ORIGIN + fcentroid));
     }
   std::cerr << "area = " << totarea << ", "
	     << "vararea = " << vararea << "\n";
   return vararea/totarea;
 }

int main() {
    Point_3 p( 1, 0, 0);
    Point_3 q( 0, 1, 0);
    Point_3 r( 0, 0, 1);
    Point_3 s( 0, 0, 0);
    Polyhedron P;
    P.make_tetrahedron( p, q, r, s);
    std::transform( P.facets_begin(), P.facets_end(), P.planes_begin(),
                    Plane_equation());
    CGAL::set_pretty_mode( std::cout);
    std::copy( P.planes_begin(), P.planes_end(),
               std::ostream_iterator<Plane_3>( std::cout, "\n"));
    return 0;
}