C++ Triangulation::is_infinite方法代码示例

void Foam::DelaunayMeshTools::writeProcessorInterface
(
    const fileName& fName,
    const Triangulation& t,
    const faceList& faces
)
{
    OFstream str(fName);

    pointField points(t.number_of_finite_cells(), point::max);

    for
    (
        typename Triangulation::Finite_cells_iterator cit =
            t.finite_cells_begin();
        cit != t.finite_cells_end();
        ++cit
    )
    {
        if (!cit->hasFarPoint() && !t.is_infinite(cit))
        {
            points[cit->cellIndex()] = cit->dual();
        }
    }

    meshTools::writeOBJ(str, faces, points);
}

void test_case(int n) {

	// cout << "---- " << n << " ----" << endl;

	// Read all infected people
	std::vector<K::Point_2> infected;
	infected.reserve(n);

	for(int i=0; i<n; i++) {
		// cin >> x[i] >> y[i];
		K::Point_2 p;
		cin >> p;
		infected.push_back(p);
		cout << "Read in point " << p << endl;
	}

	// Construct Delauney triangulation
	Triangulation t;
	t.insert(infected.begin(), infected.end());

	// Read all healthy people
	int m;
	cin >> m;

	for(int i=0; i<m; i++) {
		K::Point_2 escaper;
		long d;
		cin >> escaper >> d;

		// --- Find an escape path for this person ---
		// Find out at which face we are
		Face_handle current_face = t.locate(escaper);

		// Check if we are already outside
		if(t.is_infinite(current_face)) {
			cout << "y";
			continue;
		}
		// Check if we are already getting infected
		/*K::Point_2 nearest_infected = t.nearest_vertex(escaper, current_face)->point();
		cout << "Nearest infected person: " << nearest_infected << endl;
		int dx = nearest_infected.x() - escaper.x();
		int dy = nearest_infected.y() - escaper.y();
		long nearest_sqd = dx * dx + dy * dy;
		if(nearest_sqd < d) {
			cout << "n";
			continue;
		}*/

		// Recurse
		vector<Face_handle> visited;
		bool result = recurse(current_face, d, visited, t);

		if(result)
			cout << "POSSIBLE TO ESCAPE" << endl;
		else
			cout << "COULDN'T ESCAPE :(" << endl;
	}

}

// ---------------------------------------------------------
// initialize
// ----------
// Initialize some of the attributes of the triangulation.
// ---------------------------------------------------------
void
initialize(Triangulation &triang)
{
   // set vertex id.
   int id = 0;
   for(FVI vit = triang.finite_vertices_begin();
	vit != triang.finite_vertices_end(); vit ++)
   {
	   vit->id = id++;
	   vit->visited = false;
    	   vit->bad = false;
    	   vit->bad_neighbor = false;
   }

   // set cell id.
   id = 0;
   for(ACI cit = triang.all_cells_begin();
	cit != triang.all_cells_end(); cit ++)
   {
	   cit->id = id++;
    	   cit->visited = false;
    	   cit->outside = false;
    	   cit->transp = false;

	   for(int id = 0 ; id < 4; id++)
	   {
      		cit->set_cocone_flag(id,false);
      		cit->neighbor(id)->set_cocone_flag(cit->neighbor(id)->index(cit),false);
      		cit->bdy[id] = false;
      		cit->opaque[id] = false;
      		for(int k = 0; k < 4; k ++)
	    		cit->umbrella_member[id][k] = -1;
    	   }
    	   
	   // set the convex hull points.
    	   if(! triang.is_infinite(cit)) continue;

    	   for(int i = 0; i < 4; i ++)
    	   {
	    	   if(! triang.is_infinite(cit->vertex(i))) continue;
	    	   cit->vertex((i+1)%4)->set_convex_hull(true);
	    	   cit->vertex((i+2)%4)->set_convex_hull(true);
	    	   cit->vertex((i+3)%4)->set_convex_hull(true);
    	   }
   }
}

bool          
is_inf_VF(const Triangulation& triang,
          const Cell_handle& c, const int uid, const int vid)
{
   Facet_circulator fcirc = triang.incident_facets(Edge(c,uid,vid));
   Facet_circulator begin = fcirc;
   do{
      Cell_handle cur_c = (*fcirc).first;
      if( triang.is_infinite( cur_c ) ) return true;
      fcirc ++;
   } while(fcirc != begin);
   return false;
}

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

int
main(int,char*[])
{
  Triangulation t;

  t.insert(Point(0.1,0,1));
  t.insert(Point(1,0,1));
  t.insert(Point(0.2,0.2, 2));
  t.insert(Point(0,1,2));
  t.insert(Point(0,2,3));

  vertex_iterator vit, ve;
  // Associate indices to the vertices
  int index = 0;
  // boost::tie assigns the first and second element of the std::pair
  // returned by boost::vertices to the variables vit and ve
  for(boost::tie(vit,ve) = vertices(t); vit!=ve; ++vit ){
    vertex_descriptor  vd = *vit;
    if(! t.is_infinite(vd)){
      vertex_id_map[vd]= index++;
    }
  }

  std::cerr << index << " vertices" << std::endl;
  index = 0;
  face_iterator fit,fe;
  for(boost::tie(fit,fe) = faces(t); fit!= fe; ++fit){
    face_descriptor fd = *fit;
    halfedge_descriptor hd = halfedge(fd,t);
    halfedge_descriptor n = next(hd,t);
    
    halfedge_descriptor nn = next(n,t);
    if(next(nn,t) != hd){
      std::cerr << "the face is not a triangle" << std::endl;
    }
    
    ++index;
  }
  
  std::cerr << index << " faces" << std::endl;
  index = 0;

  edge_iterator eit,ee;
  for(boost::tie(eit,ee) = edges(t); eit!= ee; ++eit){
    edge_descriptor ed = *eit;
    vertex_descriptor vd = source(ed,t);
    CGAL_USE(vd);
    ++index;
  }

  std::cerr << index << " edges" << std::endl;
  index = 0;

  halfedge_iterator hit,he;
  for(boost::tie(hit,he) = halfedges(t); hit!= he; ++hit){
    halfedge_descriptor hd = *hit;
    vertex_descriptor vd = source(hd,t);
    CGAL_USE(vd);
    ++index;
  }
  std::cerr << index << " halfedges" << std::endl;

  std::cerr << num_vertices(t) << " " << num_edges(t) << " " << num_halfedges(t) << " " << num_faces(t) << std::endl;

  typedef boost::property_map<Triangulation, boost::vertex_point_t>::type Ppmap;
  Ppmap ppmap = get(boost::vertex_point, t);
 

  for(vertex_descriptor vd : vertices_around_target(*vertices(t).first, t)){
    std::cout <<  ppmap[vd] << std::endl;
  }


  ppmap[*(++vertices(t).first)] = Point(78,1,2);
  std::cout << " changed point of vertex " << ppmap[*(++vertices(t).first)] << std::endl;

  return 0;
}

bool is_degenerate_VF(const Triangulation& triang, const Cell_handle& c, const int& fid, const int& uid, const int& vid, const Point& d, const char* prefix)
{
	char degen_op_filename[200];
	strcat(strcpy(degen_op_filename, prefix), ".degen_VF");
	// an extra check - probably not needed.
	if (triang.is_infinite(c) || triang.is_infinite(c->neighbor(fid)) || triang.is_infinite(c->neighbor(6 - fid - uid - vid)))
	{
		return true;
	}
	vector<Cell_handle> VF;
	Facet_circulator fcirc = triang.incident_facets(Edge(c,uid,vid));
	Facet_circulator begin = fcirc;
	do
	{
		if (triang.is_infinite((*fcirc).first))
		{
			cerr << "< Inf VF >";
			return true; // by check-1 it is degenerate.
		}
		Cell_handle cur_c = (*fcirc).first;
		int cur_fid = (*fcirc).second;
		// check if cur_c and its cur_fid neighbors are cospherical.
		if (is_cospherical_pair(triang, Facet(cur_c,cur_fid)))
		{
			cerr << "< Cosph VF >";
			return true; // by check-2 it is degenerate.
		}
		fcirc ++;
	}
	while (fcirc != begin);
	// check-3
	Point vv[3];
	vv[0] = c->voronoi();
	vv[1] = c->neighbor(fid)->voronoi();
	vv[2] = c->neighbor(6 - fid - uid - vid)->voronoi();
	Vector v[3];
	v[0] = vv[0] - d;
	v[1] = vv[1] - d;
	v[2] = vv[2] - d;
	Vector v1xv0 = CGAL::cross_product(v[1], v[0]);
	Vector v0xv2 = CGAL::cross_product(v[0], v[2]);
	if (CGAL::to_double(v1xv0 * v0xv2) < 0)
	{
		ofstream fout;
		fout.open(degen_op_filename, ofstream::app);
		fout << "# prob : v1xv0 * v0xv2 = " <<
			 CGAL::to_double(v1xv0 * v0xv2) << endl;
		fout << "{LIST " << endl;
		fout << "# VF - color yellow " << endl;
		draw_VF(triang, Edge(c, uid, vid), 1, 1, 0, 1, fout);
		fout << "# v0 : segment(driver, voronoi(c)) - color red " << endl;
		draw_segment(Segment(d, vv[0]), 1, 0, 0, 1, fout);
		fout << "# v1 : segment(driver, voronoi(c->neighbor1)) - color green " << endl;
		draw_segment(Segment(d, vv[1]), 0, 1, 0, 1, fout);
		fout << "# v2 : segment(driver, voronoi(c->neighbor2)) - color blue " << endl;
		draw_segment(Segment(d, vv[2]), 0, 0, 1, 1, fout);
		fout << "}" << endl;
		fout.close();
		cerr << "< - v1xv0 * v0xv2 < 0 - >";
		return true;
	}
	return false;
}

void
grow_maximum(Cell_handle c_max, 
             Triangulation& triang, 
             map<int, cell_cluster> &cluster_set)
{
   // mark it visited.
   c_max->visited = true;
   cluster_set[c_max->id].in_cluster = true;
   cluster_set[c_max->id].outside = c_max->outside;

   // Now grow the maximum through the other tetrahedra.
   vector<Facet> bdy_stack;
   for(int i = 0; i < 4; i ++)
      bdy_stack.push_back(Facet(c_max, i));
   while(! bdy_stack.empty())
   {
      Cell_handle old_c = bdy_stack.back().first;
      int old_id = bdy_stack.back().second;
      bdy_stack.pop_back();
      CGAL_assertion(old_c->visited);
      CGAL_assertion(cluster_set[old_c->id].in_cluster);

      #ifndef __INSIDE__
      CGAL_assertion( old_c->outside );
      #endif

      #ifndef __OUTSIDE__
      CGAL_assertion( ! old_c->outside );
      #endif

      Cell_handle new_c = old_c->neighbor(old_id);
      int new_id = new_c->index(old_c);

      // If the new_c is infinite then no point in checking
      // the flow.
      if(triang.is_infinite(new_c))
         continue;
      // if the flow hits the surface continue.
      // if( old_c->outside != new_c->outside)
         // continue;
      // If new_c is already visited continue.
      if(new_c->visited) 
         continue;

      // collect new_c only if new_c flows into old_c via the
      // facet in between.
      if( is_outflow(Facet(new_c, new_id)) )
      {
         // CGAL_assertion( !is_outflow(Facet(old_c, old_id)));
         if(is_outflow(Facet(old_c, old_id))) continue;
         
         // new_c has to satisfy the following.
         CGAL_assertion( !is_maxima(new_c) && !new_c->visited &&
                         !triang.is_infinite(new_c));
         
         new_c->visited = true;
         cluster_set[new_c->id].in_cluster = true;
         cluster_set[new_c->id].outside = c_max->outside;

         // put the cells accessible via new_c into bdy_stack.
         for(int i = 1; i < 4; i ++)
         {
            if(new_c->neighbor((new_id+i)%4)->visited) 
	       continue;
	    bdy_stack.push_back(Facet(new_c, (new_id+i)%4));
         }
         // put new_c into the current cluster.
         // In other words merge the current cluster and the cluster owned by
         // new_c (as it is initialized).
         add_cell_to_cluster(cluster_set, c_max->id, new_c->id);
      }
   }
}

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

		// add constrainedge to the triangulation
		t.insert_constraint( v1, v2 );

		// add edge to the list of created/existing edges
		edges.push_back( e );
	}

	// add additional and optional interior points
	for( std::vector<std::pair<math::Vec3f, math::Vec2f> >::iterator it = interiorPoints.begin(); it != interiorPoints.end(); ++it )
	{
		// update triangulation
		Vertex_handle v = t.insert( Point( it->second.x, it->second.y ) );

		// insertion may return a vertex which already exists (when the position is the same)
		if( vertexMap.find( v ) == vertexMap.end() )
			// create according MeshEx::Vertex and keep mapping to the CGAL vertices
			vertexMap[v] = createVertex( it->first );
	}


	// extract results and create triangles ----------------------------------------

	// now we have the triangulation of the convex hull of the whole problem, now we
	// have to find the faces which are inside the polygon - we mark each face with a
	// segment(inside or outside) property and by finding a face which is adjacent to
	// a infinite face, we find the segment which is outside


	// we employ some floodfilling scheme
	for( Triangulation::Finite_faces_iterator it = t.finite_faces_begin(); it != t.finite_faces_end(); ++it )
		// reset info to -1
		it->info() = -1;

	int outsideSegment = -1;
	findOutsideSegment( t, t.finite_faces_begin(), 0, -1, outsideSegment );

	if( outsideSegment == -1 )
		printf( "error : outsideSegment not found during triangulation\n" );

	for( Triangulation::Finite_faces_iterator it = t.finite_faces_begin(); it != t.finite_faces_end(); ++it )
		if( (it->info() == -1) && (!t.is_infinite(it)) )
			printf( "triangle not touched!\n" );

	// iterate over all faces of the triangulation and create edges/triangles
	for( Triangulation::Finite_faces_iterator it = t.finite_faces_begin(); it != t.finite_faces_end(); ++it )
	{
		Face_handle fh = it;

		// we are only interested in interior triangles
		if( fh->info() == outsideSegment )
			continue;

		MeshEx::Vertex *v0, *v1, *v2;

		v0 = vertexMap[ fh->vertex(0) ];
		v1 = vertexMap[ fh->vertex(1) ];
		v2 = vertexMap[ fh->vertex(2) ];

		MeshEx::Edge *e0, *e1, *e2;

		e0 = e1 = e2 = 0;

		// look for the edges in the edge vector
		for( std::vector<MeshEx::Edge*>::iterator eit = edges.begin(); eit != edges.end(); ++eit )
		{
			MeshEx::Edge *e = *eit;

			if( e->contains(v0) && e->contains(v1) )
				e0 = e;
			else
			if( e->contains(v1) && e->contains(v2) )
				e1 = e;
			else
			if( e->contains(v2) && e->contains(v0) )
				e2 = e;
		}

		// create the edges which could not be found
		if( !e0 )
		{
			e0 = createEdge( v0, v1 );
			edges.push_back(e0);
		}
		if( !e1 )
		{
			e1 = createEdge( v1, v2 );
			edges.push_back(e1);
		}
		if( !e2 )
		{
			e2 = createEdge( v2, v0 );
			edges.push_back(e2);
		}

		// create triangle
		MeshEx::Triangle *tri = createTriangle( v0, v1, v2, e0, e1, e2 );
	}
}