C++ Face_handle类代码示例

void R_s_k_2::draw_edge(const Edge& edge)
{
  int i = edge.second;
  Face_handle face = edge.first;
  Point a = face->vertex((i+1)%3)->point();
  Point b = face->vertex((i+2)%3)->point();
  draw_segment(a, b);
}

void Mesh::Vertex::computeNormal()
{
	for(int i = 0; i < numFaces(); ++i)
	{
		Face_handle v = face(i);
		m_normal += v->normal();
	}
	m_normal.unitize();
}

bool
emptyLatticeTriangle( const Delaunay & t, const Face_handle & f )
{
  if ( t.is_infinite( f ) ) return false;
  Z2i::Point a( toDGtal(f->vertex(0)->point())),
    b(toDGtal(f->vertex(1)->point())),
    c(toDGtal(f->vertex(2)->point()));
  
  Z2i::Vector ab( b - a ), ac( c - a );
  int d = ab[ 0 ] * ac[ 1 ] - ab[ 1 ] * ac[ 0 ];
  return ( d == 1 ) || (d == -1 );
}

void Foam::CV2D::external_flip(Face_handle& f, int i)
{
    Face_handle n = f->neighbor(i);

    if
    (
        CGAL::ON_POSITIVE_SIDE
     != side_of_oriented_circle(n, f->vertex(i)->point())
    ) return;

    flip(f, i);
    i = n->index(f->vertex(i));
    external_flip(n, i);
}

void Mesh::commitRemove()
{
	// first go over all the points and tell them where they go
	int count = 0;
	for(int i = 0; i < numVtx(); ++i)
	{
		Vertex_handle v = find_vertex(i);
		if (v->m_tran == -1)
			v->m_index = count++;
		else
			v->m_index = -1; // to be removed
		v->m_tran = -1;
	}
	// no go over the faces and transfer their point references
	count = 0;
	for(int i = 0; i < numFaces(); ++i)
	{
		Face_handle f = find_facet(i);
		bool remove = false;
		for(int ii = 0; ii < f->size(); ++ii)
		{
			int newvi = find_vertex(f->m_pi[ii])->m_index;
			f->m_pi[ii] = newvi;
			remove |= (newvi == -1);
		}
		if (!remove)
			f->m_index = count++;
		else 
			f->m_index = -1;
	}
	// actually remove the vertices
	Vertex_iterator vit = vertices_begin();
	while(vit != vertices_end())
	{
		if (vit->m_index == -1)
			vit = m_vtx.erase(vit);
		else
			++vit;
	}

	// remove faces
	Face_iterator fit = faces_begin();
	while(fit != faces_end())
	{
		if (fit->m_index == -1)
			fit = m_face.erase(fit);
		else
			++fit;
	}
}

void Foam::CV2D::fast_restore_Delaunay(Vertex_handle vh)
{
    int i;
    Face_handle f = vh->face(), next, start(f);

    do
    {
        i=f->index(vh);
        if (!is_infinite(f))
        {
            if (!internal_flip(f, cw(i))) external_flip(f, i);
            if (f->neighbor(i) == start) start = f;
        }
        f = f->neighbor(cw(i));
    } while (f != start);
}

 virtual void after_split_face (Face_handle /* old_face */,
                                Face_handle new_face, bool )
 {
   // Assign index to the new face.
   new_face->set_data (n_faces);
   n_faces++;
 }

bool recurse(Face_handle current_face, long d, vector<Face_handle>& visited, Triangulation t) {

	// Remember we visited this node
	visited.push_back(current_face);

	cout << "Recursion at: " << t.dual(current_face) << endl;

	// Base of recursion: check if we are free
	if(t.is_infinite(current_face)) {
		cout << "Infinite face!" << endl;
		return true;
	}

	for(int neighbor_num=0; neighbor_num<3; neighbor_num++) {
		Face_handle neighbor_face = current_face->neighbor(neighbor_num);
		//cout << (current_face == neighbor_face) << endl;
		cout << "\tChecking neighbor of vertex " << current_face->vertex(neighbor_num)->point() << endl;

		cout << "\tPoints: ";
		for(int j=0; j<3; j++) {
			cout << neighbor_face->vertex(j)->point() << ", ";
		}
		cout << endl;

		// If we already visited
		if(find(visited.begin(), visited.end(), current_face) != visited.end()) {
			continue;
		}

		Vertex_handle border_endpoint1 = current_face->vertex((neighbor_num + 1) % 3);
		Vertex_handle border_endpoint2 = current_face->vertex((neighbor_num + 2) % 3);
		K::FT border_length_sq = CGAL::squared_distance(border_endpoint1->point(), border_endpoint2->point());

		if(CGAL::to_double(border_length_sq) >= 4 * d) {	// If we can fit through that edge 
			// cout << "can fit through edge " << neighbor_num << endl;

			// New search starting from neighbor
			if(recurse(neighbor_face, d, visited, t)) {
				return true;
			}
		} else {
			cout << "cannot fit through edge :S" << endl;
		}
	}

	return false;
}

void ApplyObjects(Pmwx& ioMap)
{
	if (gFAAObs.empty()) return;

	Point_2	sw, ne;
	CalcBoundingBox(ioMap, sw, ne);
// 	ioMap.Index();

	int	placed = 0;

//	CGAL::Arr_landmarks_point_location<Arrangement_2>	locator(gMap);
	CGAL::Arr_walk_along_line_point_location<Arrangement_2>	locator(gMap);

	for (FAAObsTable::iterator i = gFAAObs.begin(); i != gFAAObs.end(); ++i)
	{
		if (i->second.kind != NO_VALUE)
		{
			Point_2 loc = Point_2(i->second.lon, i->second.lat);

			DebugAssert(CGAL::is_valid(gMap));
			CGAL::Object obj = locator.locate(loc);
			Face_const_handle ff;
			if(CGAL::assign(ff,obj))
			{
				Face_handle f = ioMap.non_const_handle(ff);
				GISPointFeature_t	feat;
					feat.mFeatType = i->second.kind;
					feat.mLocation = loc;
					if (i->second.agl != DEM_NO_DATA)
						feat.mParams[pf_Height] = i->second.agl;
					feat.mInstantiated = false;
					f->data().mPointFeatures.push_back(feat);
					++placed;
	#if 0
					printf("Placed %s at %lf, %lf\n",
						FetchTokenString(i->second.kind), i->second.lon, i->second.lat);
	#endif
//				if (v.size() > 1)
//					fprintf(stderr,"WARNING (%d,%d): Point feature %lf, %lf matches multiple areas.\n",gMapWest, gMapSouth, CGAL::to_double(loc.x()), CGAL::to_double(loc.y()));

			}
		}
	}
	printf("Placed %d objects.\n", placed);
}

void BuildCGALPolygon<HDS>::operator()( HDS& hds){
	CGAL::Polyhedron_incremental_builder_3<HDS> B( hds, true);
	B.begin_surface( this->pts->GetNumberOfPoints(), this->polys->GetNumberOfCells(), 0);
	typedef typename HDS::Vertex   Vertex;
	typedef typename Vertex::Point Point;
	vtkIdType npts = 0;
	vtkIdType *indx;
	double vertex[3];
	for (int i = 0; i != this->pts->GetNumberOfPoints(); i++ ){
		this->pts->GetPoint(i, vertex);
		B.add_vertex( Point( vertex[0], vertex[1], vertex[2]));

	}
	int j = 0;
	Face_handle face; 
	for (this->polys->InitTraversal(); this->polys->GetNextCell(npts,indx); ){
		if(indx[0] != indx[1] & indx[0] != indx[2] & indx[1] != indx[2]){ 
			//VTK polygons can contain lines where two vertexes are identical. Forget these
			if (B.test_facet(indx, indx+3)){ 
				face = B.begin_facet();
				B.add_vertex_to_facet( indx[0]);
				B.add_vertex_to_facet( indx[1]);
				B.add_vertex_to_facet( indx[2]);
				B.end_facet();
				//cout << this->IoletIdArray->GetValue(j) << endl;
				face->id() = IoletIdArray->GetValue(j) + 2;
				//the face id is size_t i.e. unsigned so we shift this to positive. 1 is wall. 2,3 ... are 
				//the inlets and outlets.  
			}
			else
				cout << "Ignoring Non manifold facet between: " << indx[0] << " " <<  indx[1] << " " << indx[2] << endl;
		}
		else{
			cout << "Eleminated degenerate vertex: " << indx[0] << " " <<  indx[1] << " " << indx[2] << endl;
		}
		++j;
	}
	
	B.end_surface();
	//cout << B.check_unconnected_vertices () << endl;
	B.remove_unconnected_vertices();
	//cout << B.check_unconnected_vertices () << endl;
}

void findOutsideSegment( Triangulation &t, Face_handle fh, int currentSegment, int commingFromIndex, int &outsideSegment )
{

	// if the face is an infinite face
	// then we know the outside segment which is the
	// current segment
	if( t.is_infinite(fh) )
	{
		if( (outsideSegment != -1)&&(outsideSegment != currentSegment) )
			printf( "error : different outsideSegments detected during triangulation (edgeloop not closed?)\n" );
		outsideSegment = currentSegment;
		return;
	}

	// if there is already a segment identifier, then we know that
	// the face has already been visited
	if( fh->info() != -1 )
		return;

	fh->info() = currentSegment;

	for( int i=0; i<3; ++i )
	{
		// get edge associated with the index i
		std::pair<Face_handle, int> edge = std::make_pair( fh, i );

		if( i == commingFromIndex )
			continue;
		

		// if the edge is a constrained edge, then we know this is the border
		if(t.is_constrained(edge))
		{
			//...and we have to pass a madified currentSegment value
			findOutsideSegment( t, fh->neighbor(i), (currentSegment+1)%2, t.mirror_index( fh, i), outsideSegment );
		}else
			//...else we recurse and leave the currentSegment value untouched
			findOutsideSegment( t, fh->neighbor(i), currentSegment, t.mirror_index( fh, i), outsideSegment );
	}
}

void
discoverInfiniteComponent(const CDT & ct)
{
  //when this function is called, all faces are set "in_domain"
  Face_handle start = ct.infinite_face();
  std::list<Face_handle> queue;
  queue.push_back(start);

  while(! queue.empty())
  {
    Face_handle fh = queue.front();
    queue.pop_front();
    fh->set_in_domain(false);
	
    for(int i = 0; i < 3; i++)
    {
      Face_handle fi = fh->neighbor(i);
      if(fi->is_in_domain()
        && !ct.is_constrained(CDT::Edge(fh,i)))
        queue.push_back(fi);
    }
  }
}

void CriticalCurves::setParameters(double radius_1, double radius_2, Arrangements_2 insets_1, Arrangements_2 insets_2)
{
    Arrangement_2_iterator inset_1 = insets_1.begin();
    Arrangement_2_iterator inset_2 = insets_2.begin();

    while (inset_1 != insets_1.end() && inset_2 != insets_2.end())
    {
        Arrangement_2 arrangement;

        // Add the curves of the inset.
        for (Edge_iterator edge = inset_1->edges_begin(); edge != inset_1->edges_end(); ++edge)
        {
            insert(arrangement, edge->curve());
        }

        // Add the critical curves of type I.
        for (Edge_iterator edge = inset_2->edges_begin(); edge != inset_2->edges_end(); ++edge)
        {
            if (CGAL::COLLINEAR == edge->curve().orientation())
            {
                // Displaced a segment.
                Nt_traits nt_traits;
                Algebraic_ft factor = nt_traits.convert(Rational(radius_1) + Rational(radius_2));
                Conic_point_2 source = edge->curve().source();
                Conic_point_2 target = edge->curve().target();
                Algebraic_ft delta_x = target.x() - source.x();
                Algebraic_ft delta_y = target.y() - source.y();
                Algebraic_ft length = nt_traits.sqrt(delta_x * delta_x + delta_y * delta_y);
                Algebraic_ft translation_x = factor * delta_y / length;
                Algebraic_ft translation_y = - factor * delta_x / length;
                Conic_point_2 point_1(source.x() + translation_x, source.y() + translation_y);
                Conic_point_2 point_2(target.x() + translation_x, target.y() + translation_y);
                Algebraic_ft a = - delta_y;
                Algebraic_ft b = delta_x;
                Algebraic_ft c = factor * length - (source.y() * target.x() - source.x() * target.y());
                X_monotone_curve_2 x_monotone_curve(a, b, c, point_1, point_2);
                insert(arrangement, x_monotone_curve);
            }
            else
            {
                // Displaces an arc.
                Rational two(2);
                Rational four(4);

                Rational r = edge->curve().r();
                Rational s = edge->curve().s();
                Rational t = edge->curve().t();
                Rational u = edge->curve().u();
                Rational v = edge->curve().v();
                Rational w = edge->curve().w();

                Nt_traits nt_traits;
                Rational x_center = - u / (two * r);
                Rational y_center = - v / (two * r);
                Rat_point_2 rat_center(x_center, y_center);
                Conic_point_2 center(nt_traits.convert(x_center), nt_traits.convert(y_center));

                Rational radius = Rational(radius_1) + two * Rational(radius_2);

                Algebraic_ft coefficient = nt_traits.convert(radius / Rational(radius_2));

                Conic_point_2 source_1 = edge->curve().source();
                Algebraic_ft x_source_2 = center.x() + coefficient * (source_1.x() - center.x());
                Algebraic_ft y_source_2 = center.y() + coefficient * (source_1.y() - center.y());
                Conic_point_2 source_2(x_source_2, y_source_2);

                Conic_point_2 target_1 = edge->curve().target();
                Algebraic_ft x_target_2 = center.x() + coefficient * (target_1.x() - center.x());
                Algebraic_ft y_target_2 = center.y() + coefficient * (target_1.y() - center.y());
                Conic_point_2 target_2(x_target_2, y_target_2);

                Rat_circle_2 circle(rat_center, radius * radius);

                Conic_arc_2 conic_arc(circle, CGAL::COUNTERCLOCKWISE, source_2, target_2);

                insert(arrangement, conic_arc);
            }
        }

        // Add the critical curves of type II.
        for (Edge_iterator edge = inset_2->edges_begin(); edge != inset_2->edges_end(); ++edge)
        {
            double x = CGAL::to_double(edge->curve().source().x());
            double y = CGAL::to_double(edge->curve().source().y());
            double radius = radius_1 + radius_2;
            Rat_point_2 center(x, y);
            Rat_circle_2 circle(center, radius * radius);
            Conic_arc_2 conic_arc(circle);
            insert(arrangement, conic_arc);
        }

        // Remove the curves which are not include in the inset.
        Objects objects;
        Face_handle face;
        for (Edge_iterator edge = arrangement.edges_begin(); edge != arrangement.edges_end(); ++edge)
        {
            CGAL::zone(*inset_1, edge->curve(), std::back_inserter(objects));
            for (Object_iterator object = objects.begin(); object != objects.end(); ++object)
            {
                if (assign(face, *object))
//.........这里部分代码省略.........

/**
 * Returns an alpha shape but as lines instead of only points.
 */
std::list<std::pair<Shared_Point,Shared_Point> > Hull::alphaHull2D(const std::list<Point> points,double alpha){
	//iHt dt;
	std::list<std::pair<Shared_Point,Shared_Point> > resultingEdges;
	if(points.size() < 3) return resultingEdges;

	#ifdef NO_HULL_CALCULATION
	if(points.size() == 4){
		std::list<Point>::const_iterator it = points.begin();
		Shared_Point p0(new Point(it->p[0], 0, it->p[2]));
		it++;
		Shared_Point p1(new Point(it->p[0], 0, it->p[2]));
		it++;
		Shared_Point p2(new Point(it->p[0], 0, it->p[2]));
		it++;
		Shared_Point p3(new Point(it->p[0], 0, it->p[2]));

		resultingEdges.push_back(std::pair<Shared_Point,Shared_Point>(p0,p1));
		resultingEdges.push_back(std::pair<Shared_Point,Shared_Point>(p1,p2));
		resultingEdges.push_back(std::pair<Shared_Point,Shared_Point>(p2,p3));
		resultingEdges.push_back(std::pair<Shared_Point,Shared_Point>(p3,p0));
		
		return resultingEdges;
	}else{
			fprintf(stderr,"Cannot Compute no hull\n");
	}
	#endif

	std::vector<K::Point_2> dt;
	dt.reserve(points.size());

	for(std::list<Point>::const_iterator it=points.begin(); it!=points.end(); it++){
		K::Point_2 p((*it).p[0],(*it).p[2]);
		dt.push_back(p);
	}

	try{
        Alpha_shape_2 as(dt.begin(),dt.end(),alpha);//Alpha_shape_2::REGULARIZED);
		//Alpha_shape_2 as(dt.begin(),dt.end(),Alpha_shape_2::REGULARIZED);
		//std::cout << "Alpha shape computed in REGULARIZED mode by defaut."<< std::endl;

		//Alpha_iterator opt = as.find_optimal_alpha(1);
		//as.set_alpha(*opt);

		for(Alpha_shape_2::Alpha_shape_edges_iterator it = as.alpha_shape_edges_begin(); it != as.alpha_shape_edges_end(); it++){
			Edge e = (*it);
			Face_handle f = e.first;
			int i = e.second;

			Alpha_shape_2::Vertex_handle vh1 = f->vertex(f->cw(i));
			Alpha_shape_2::Vertex_handle vh2 = f->vertex(f->ccw(i));
			K::Point_2& p1 = vh1->point();
			K::Point_2& p2 = vh2->point();
			
			std::pair<Shared_Point,Shared_Point> pair(Shared_Point(new Point(p1.x(),0,p1.y())),Shared_Point(new Point(p2.x(),0,p2.y())));
			resultingEdges.push_back(pair);
		}

	}catch( ...){
		printf("Catched Unknown Cgal Exception adding plane without generatig shape\n");
	}
	return resultingEdges;

}

int main()
{
  Triangulation t;
  Face_handle fh;

  // Check the empty triangulation
  fh = test_point_location(t, Point(0.5, 0.5), Triangulation::EMPTY);
  CGAL_assertion(fh == Face_handle());

  // Insert the first point
  Point p0(0.5, 0.5);
  Vertex_handle vh0 = t.insert(p0);
  CGAL_assertion(t.is_valid(true));
  CGAL_USE(vh0);

  fh = test_point_location(t, p0, Triangulation::VERTEX);
  CGAL_assertion(fh->has_vertex(vh0));

  fh = test_point_location(t, p0 + Vector(0.1, 0.1), Triangulation::EDGE);
  CGAL_assertion(fh->has_vertex(vh0));

  fh = test_point_location(t, p0 + Vector(-0.1, -0.1), Triangulation::EDGE);
  CGAL_assertion(fh->has_vertex(vh0));

  fh = test_point_location(t, p0 + Vector(-0.2, -0.3), Triangulation::FACE);
  CGAL_assertion(fh->has_vertex(vh0));

  CGAL_assertion(t.is_valid(true));

  // Insert the second point on an edge
  Point p1(0.7, 0.7);
  Vertex_handle vh1 = t.insert(p1);
  CGAL_USE(vh1);
  CGAL_assertion(t.is_valid(true));

  fh = test_point_location(t, p0, Triangulation::VERTEX);
  CGAL_assertion(fh->has_vertex(vh0));

  fh = test_point_location(t, p1, Triangulation::VERTEX);
  CGAL_assertion(fh->has_vertex(vh1));

  fh = test_point_location(t, p0 + Vector(0.1, 0.1), Triangulation::EDGE);
  CGAL_assertion(fh->has_vertex(vh0));
  CGAL_assertion(fh->has_vertex(vh1));

  fh = test_point_location(t, p0 + Vector(-0.1, -0.1), Triangulation::EDGE);
  CGAL_assertion(fh->has_vertex(vh0));
  CGAL_assertion(!fh->has_vertex(vh1));

  fh = test_point_location(t, p1 + Vector(0.1, 0.1), Triangulation::EDGE);
  CGAL_assertion(!fh->has_vertex(vh0));
  CGAL_assertion(fh->has_vertex(vh1));

  fh = test_point_location(t, p0 + Vector(-0.02, -0.03), Triangulation::FACE);
  CGAL_assertion(fh->has_vertex(vh0));

  fh = test_point_location(t, p1 + Vector(-0.02, -0.03), Triangulation::FACE);
  CGAL_assertion(fh->has_vertex(vh1));

  CGAL_assertion(t.is_valid(true));

  // Insert the third point in a face
  Point p2(0.8, 0.6);
  Vertex_handle vh2 = t.insert(p2);
  CGAL_USE(vh2);
  CGAL_assertion(t.is_valid(true));

  fh = test_point_location(t, p0, Triangulation::VERTEX);
  CGAL_assertion(fh->has_vertex(vh0));
  fh = test_point_location(t, p1, Triangulation::VERTEX);
  CGAL_assertion(fh->has_vertex(vh1));
  fh = test_point_location(t, p2, Triangulation::VERTEX);
  CGAL_assertion(fh->has_vertex(vh2));

  fh = test_point_location(t, Point(0.6, 0.6), Triangulation::EDGE);
  CGAL_assertion(fh->has_vertex(vh0));
  CGAL_assertion(fh->has_vertex(vh1));

  test_point_location(t, Point(0.7, 0.6), Triangulation::FACE);
  test_point_location(t, p0 + Vector(-0.02, -0.03), Triangulation::FACE);
  test_point_location(t, p0 + Vector(0.02, -0.03), Triangulation::FACE);
  test_point_location(t, p0 + Vector(-0.02, 0.03), Triangulation::FACE);
  test_point_location(t, p0 + Vector(0.02, 0.03), Triangulation::FACE);

  return 0;
}