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

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

int main(int argc, char* argv[]) {
    CGAL::Timer user_time;
    cerr << "Loading OFF file ... " << endl;
    user_time.start();
    Polyhedron P;
    std::ifstream in1((argc>1)?argv[1]:"data/tetra_intersected_by_triangle.off");
    in1 >> P;
    cerr << "Loading OFF file   : " << user_time.time() << " seconds." << endl;
    if ( ! P.is_pure_triangle()) {
        cerr << "The input object is not triangulated. Cannot intersect."
                  << endl;
        exit(1);
    }
    user_time.reset();
    cerr << "Intersection ... " << endl;
    intersection( P);
    cerr << "Intersection       : " << user_time.time() << " seconds." << endl;
    write_off();

    return 0;
}

void set_semantic_AABB_C2V(Polyhedron& exteriorPolyhe,PolVector& polyVec) {

	if (exteriorPolyhe.is_pure_triangle()) {
		std::transform( exteriorPolyhe.facets_begin(), exteriorPolyhe.facets_end(),exteriorPolyhe.planes_begin(),Plane_equation());
		std::vector<std::string>	semList;
		std::vector<std::shared_ptr<AAbbTree>>		treeList;

		// Build Trees. One for each semantic
		for(PolVector::iterator pvIt = polyVec.begin();pvIt!=polyVec.end();++pvIt) {// Get AABB trees of all semantics
			if (pvIt->is_pure_triangle()) {
				std::string semP = pvIt->facets_begin()->semanticBLA;
				std::vector<std::string>::iterator  strIt = std::find(semList.begin(), semList.end(),semP);
				if (strIt==semList.end()) {											// If new sematic
					semList.push_back(semP);										// Add sem
					std::shared_ptr<AAbbTree> tree = std::make_shared<AAbbTree>(pvIt->facets_begin(),pvIt->facets_end());		// Create tree
					tree->accelerate_distance_queries();							// accelerate
					treeList.push_back(tree);										// Add tree
				} else																					// If not new
					treeList[strIt-semList.begin()]->insert(pvIt->facets_begin(),pvIt->facets_end());	// Append to tree
			} else std::cerr << "ERROR: Not pure triangle (set_semantic_AABB2C2V)" << std::endl;
		}
		

		// For each facet calculate the least distance to each tree
		std::string semListStr = boost::algorithm::join((semList), " ");
		int percCount = 1;
		Polyhedron::Facet_iterator exfIt;						// Iterate over exterior faces
		for (exfIt = exteriorPolyhe.facets_begin(); exfIt != exteriorPolyhe.facets_end(); ++exfIt,++percCount) {	
			
			std::cout << "\r"<<semListStr<<". ("<<100*percCount/exteriorPolyhe.size_of_facets()<<"%)";

			Vector_3 orthVec = exfIt->plane().orthogonal_vector();
			normalizeVector(orthVec);	//if (!normalizeVector(ortVec)) continue;

			std::vector<distSemFace> dsfList(semList.size());										
			Point_3 centerPoint = comp_facetCentroid(exfIt);						// Compute centroid
			std::vector<Kernel::FT> leastSemDistances;
			for (int intIt=0;intIt<(int)treeList.size();++intIt) {					// Loop over all trees
				AAbbTree::Point_and_primitive_id pp = treeList[intIt]->closest_point_and_primitive(centerPoint);
				dsfList[intIt].dist	= CGAL::squared_distance(centerPoint,pp.first);	// Store distance semantic and facet for each tree
				dsfList[intIt].sem	= semList[intIt];
				dsfList[intIt].fh	= pp.second;
			}

			std::sort(dsfList.begin(),dsfList.end(),by_dist());

			exfIt->leastSqDistance = dsfList[0].dist;					// least sqrt distance
			if (exfIt->isMinkFacet = dsfList[0].dist > SEMANTIC_DISTANCE_THRESHOLD) {
				exfIt->semanticBLA = TO_DIST_SEMANTIC;						// Default semantic if too distant
				continue;
			} else
				exfIt->semanticBLA = dsfList[0].sem;					// Semantics of closest

			Vector_3 faceNormal;
			Kernel::FT faceSqArea;
			double minAngle = 10;
			Kernel::FT maxArea= 0;

			for (std::vector<distSemFace>::iterator slIt = dsfList.begin();slIt != dsfList.end();++slIt)// HANDLE ANYTHING AS LESS IMPORTANT
				if (slIt->dist < dsfList[0].dist+OVERLAP_DIST_THRESHOLD) {	// Check if Equidistant
					pointVector facetPoints = comp_facetPoints(exfIt);
					CGAL::normal_vector_newell_3(facetPoints.begin(),facetPoints.end(),faceNormal); // Calculate normal vector, ortVec set to zero in newell
					double angle = comp_angle(orthVec,faceNormal);
					if (angle!=-1 && angle < minAngle+OVERLAP_ANGLE_THRESHOLD) {
						if (minAngle >= angle+OVERLAP_ANGLE_THRESHOLD)		exfIt->equidistSems.clear();
						if (angle < minAngle)								minAngle = angle;
						faceSqArea = comp_facetSquaredArea(facetPoints);
						if (faceSqArea>maxArea-OVERLAP_AREA_THRESHOLD) {
							if (maxArea<=faceSqArea-OVERLAP_AREA_THRESHOLD)	exfIt->equidistSems.clear();
							if (faceSqArea>maxArea)							maxArea = faceSqArea;
							exfIt->equidistSems.push_back(slIt->sem);				// Add equidist semantics
						}
					}
				}
		}
		std::cout << "\r"<<semListStr<<". (100%)" << std::endl;
	}else std::cerr << "ERROR: Not pure triangle (set_semantic_AABB2C2V)" << std::endl;
}

TrianglesList meshSimplification(TrianglesList &triangles, int stopPredicate) {

	#ifdef MESHSIMPLIFICATION_LOG
	CGAL::Timer timer;
	timer.start();
	#endif

	TrianglesList result;

	try
	{
		Polyhedron P;

		#ifdef MESHSIMPLIFICATION_LOG
		std::cout << "Start Building Polyhedron surface... " << std::endl;
		#endif

		Build_triangle_mesh_coherent_surface<HalfedgeDS> triangle(triangles);
		P.delegate(triangle);
		P.normalize_border();

		#ifdef MESHSIMPLIFICATION_LOG
		std::cout << "Completed Building Polyhedron surface:" << std::endl;
		std::cout << "Polyhedron is_pure_triangle: " << P.is_pure_triangle() << std::endl;
		std::cout << "Polyhedron is_closed: " << P.is_closed() << std::endl;
		std::cout << "Polyhedron is_pure_bivalent : " << P.is_pure_bivalent () << std::endl;
		std::cout << "Polyhedron is_pure_trivalent: " << P.is_pure_trivalent() << std::endl;
		std::cout << "Polyhedron is_valid 0: " << P.is_valid(false, 0) << std::endl;
		std::cout << "Polyhedron is_valid 1: " << P.is_valid(false, 1) << std::endl;
		std::cout << "Polyhedron is_valid 2: " << P.is_valid(false, 2) << std::endl;
		std::cout << "Polyhedron is_valid 3: " << P.is_valid(false, 3) << std::endl;
		std::cout << "Polyhedron is_valid 4: " << P.is_valid(false, 4) << std::endl;
		std::cout << "Polyhedron normalized_border_is_valid : " << P.normalized_border_is_valid(false) << std::endl;
		#endif

		#ifdef MESHSIMPLIFICATION_LOG
		std::cout << "Start edge_collapse... " << std::endl;
		#endif

		SMS::Count_stop_predicate<Polyhedron> stop(stopPredicate);

		int removedEdges = SMS::edge_collapse(P, stop,
				CGAL::vertex_index_map(boost::get(CGAL::vertex_external_index, P)).edge_index_map(boost::get(CGAL::edge_external_index ,P))
		);

		#ifdef MESHSIMPLIFICATION_LOG
		std::cout << "Completed edge_collapse:" << std::endl;
		std::cout << "Finished with: " << removedEdges << " edges removed and "  << (P.size_of_halfedges()/2) << " final edges." << std::endl;
		#endif

		//Build output result
		for ( Polyhedron::Facet_iterator fit( P.facets_begin() ), fend( P.facets_end() ); fit != fend; ++fit )
		{
			if ( fit->is_triangle() )
			{
				PointCGAL verts[3];
				int tick = 0;

				Polyhedron::Halfedge_around_facet_circulator hit( fit->facet_begin() ), hend( hit );
				do
				{
					if ( tick < 3 )
					{
						verts[tick++] = PointCGAL( hit->vertex()->point().x(), hit->vertex()->point().y(), hit->vertex()->point().z() );
					}
					else
					{
						std::cout << "meshSimplification: We've got facets with more than 3 vertices even though the facet reported to be triangular..." << std::endl;
					}

				} while( ++hit != hend );

				result.push_back( Triangle(verts[0], verts[1], verts[2]) );
			}
			else
			{
				std::cout << "meshSimplification: Skipping non-triangular facet" << std::endl;
			}

		}

	}
	catch (CGAL::Assertion_exception e)
	{
		std::cout << "ERROR: meshSimplification CGAL::Assertion_exception" << e.message() << std::endl;
	}

	#ifdef MESHSIMPLIFICATION_LOG
	timer.stop();
	std::cout << "meshSimplification result with: " << result.size() << " triangles." << std::endl;
	std::cout << "Total meshSimplification time: " << timer.time() << std::endl;
	#endif

	return result;
}