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

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

// Compute all the finite voronoi vertices and the circumradius of all the finite cells.
void compute_voronoi_vertex_and_cell_radius(Triangulation& triang)
{
	bool is_there_any_problem_in_VV_computation = false;
	for (FCI cit = triang.finite_cells_begin();
			cit != triang.finite_cells_end(); cit ++)
	{
		//we tell CGAL to call our function if there is a problem
		//we also tell it not to die if things go haywire
		//CGAL::Failure_function old_ff = CGAL::set_error_handler(failure_func);
		//CGAL::Failure_behaviour old_fb = CGAL::set_error_behaviour(CGAL::CONTINUE);
		// be optimistic :-)
		//this is a global
		cgal_failed = false;
		cit->set_voronoi(triang.dual(cit));
		bool is_correct_computation = !cgal_failed;
		is_there_any_problem_in_VV_computation |= !is_correct_computation;
		if (cgal_failed)
		{
			// set cc the centroid of the cell.
			Vector cc = CGAL::NULL_VECTOR;
			for (int i = 0; i < 4; i ++)
			{
				cc = cc + (cit->vertex(i)->point() - CGAL::ORIGIN);
			}
			cc = (1./4.)*cc;
			cit->set_voronoi(CGAL::ORIGIN + cc);
		}
		//put everything back the way we found it,
		//CGAL::set_error_handler(old_ff);
		//CGAL::set_error_behaviour(old_fb);
		// set the cell radius.
		cit->set_cell_radius(CGAL::to_double((cit->vertex(0)->point()-cit->voronoi()) *(cit->vertex(0)->point()-cit->voronoi())));
	}
	return;
}

vector<int>
compute_smax(Triangulation& triang, 
             map<int, cell_cluster> &cluster_set,
             const double& mr)
{
   for(ACI cit = triang.all_cells_begin();
      cit != triang.all_cells_end(); cit ++)
   {
      cluster_set[cit->id] = cell_cluster(cit->id);
      cit->visited = false;
   }

   int max_cnt = 0;
   for(FCI cit = triang.finite_cells_begin();
      cit != triang.finite_cells_end(); cit ++)
   {
      if( ! is_maxima(cit) ) continue;

      #ifndef __OUTSIDE__
      if( cit->outside ) continue;
      #endif

      #ifndef __INSIDE__
      if( ! cit->outside ) continue;
      #endif
      
      if(max_cnt++%1000 == 0) cerr << "+";
      grow_maximum(cit, triang, cluster_set);
   }
   cerr << ".";

   // club_segment(triang, cluster_set, mr );
   // cerr << ".";
   club_contiguous_segment(triang, cluster_set );
   cerr << ".";

   // Compute the volume of each cluster. Remember after merging the 
   // 'rep' field is more useful than cluster_id.
   vector<int> cluster_ids;
   vector<double> cluster_vol;
   cluster_ids.clear();
   cluster_vol.clear();
   calc_cluster_volume_and_store_with_cluster_rep(triang, 
	  					  cluster_set,
			                          cluster_vol, 
						  cluster_ids);
   cerr << ".";

   // Sort the clusters with respect to the volumes.
   vector<int> sorted_indices;
   sorted_indices.clear();
   sort_cluster_wrt_volume(cluster_vol, cluster_ids, sorted_indices);
   cerr << ".";
   return sorted_indices;
}

int main()
{
  Triangulation tr;
  int a, b, d;
  for (a=0;a!=4;a++)
    for (b=0;b!=4;b++)
      for (d=0;d!=4;d++)
	tr.insert(Point((a*b-d*a)*10 +a ,(a-b+d +5*b)*100,
			    a*a-d*d-b));
  Triangulation::Finite_cells_iterator cit=tr.finite_cells_begin();
  for(; cit != tr.finite_cells_end(); ++cit) {
    Point circum = tr.dual(cit);
    CGAL_USE(circum);
  }
  return 0;
}

int main()
{
	Triangulation triangulation;

	boost::filesystem::path input_pathname = "C:/Carleton/CGAL-4.4/demo/Polyhedron/data/elephant.off";

	// create_cubes(triangulation, 2, 2, 1 );
	read_off( triangulation, input_pathname.string() );

	// read_off(triangulation, "C:/Carleton/Meshes/holmes_off/geometry/octahedron.off"); 
	// read_off(triangulation, "C:/Carleton/CGAL-4.4/demo/Polyhedron/data/cube.off");
	// read_off(triangulation, "C:/Carleton/CGAL-4.4/demo/Polyhedron/data/ellipsoid.off");

#if 0
	for (auto cell = triangulation.finite_cells_begin(); cell != triangulation.finite_cells_end(); ++cell)
	{
		for (int i = 0; i < 4; ++i)
		{
			Point p = cell->vertex(i)->point();
			assert(-10.0 < p.x() && p.x() < +10.0);
			assert(-10.0 < p.y() && p.y() < +10.0);
			assert(-10.0 < p.z() && p.z() < +10.0);
		}
	}
#endif

	set_cell_and_vertex_ids(triangulation);
	set_random_weights(triangulation);
	propagate_weights(triangulation);

	std::cout << "Number of finite vertices : " << triangulation.number_of_vertices() << std::endl;
	std::cout << "Number of finite edges    : " << triangulation.number_of_finite_edges() << std::endl;
	std::cout << "Number of finite facets   : " << triangulation.number_of_finite_facets() << std::endl;
	std::cout << "Number of finite cells    : " << triangulation.number_of_finite_cells() << std::endl;

	std::string filename = input_pathname.filename().stem().string() + "_tet.vtk";
	write_vtk( triangulation, filename );

	if (triangulation.number_of_finite_cells() < 100)
	{
		dump_triangulation(triangulation);
	}

	Graph graph;

	create_steiner_points(graph,triangulation);

	// the distances are temporary, so we choose an external property for that
	std::vector<double> distances(num_vertices(graph));
	std::vector<GraphNode_descriptor> predecessors(num_vertices(graph));

	boost::dijkstra_shortest_paths(
		graph, 
		*vertices(graph).first, 
		boost::weight_map(get(&GraphEdge::weight, graph)).
		distance_map(boost::make_iterator_property_map(distances.begin(), get(boost::vertex_index, graph))).
		predecessor_map(boost::make_iterator_property_map(predecessors.begin(), get(boost::vertex_index, graph)))
	);
	
	filename = input_pathname.filename().stem().string() + "_wsp.vtk";
	write_shortest_path_vtk( graph, predecessors, distances, filename );

	// write_graph_dot("graph.dot", graph);

	std::cout << "This is the end..." << std::endl;

	return EXIT_SUCCESS;
}

void
curate_tr(Triangulation& triang,
          map<int, cell_cluster> &cluster_set,
          const vector<int> &sorted_cluster_index_vector,
          const int output_pocket_count,
          const int output_tunnel_count)
{
   int tun_void_cnt = 0;
   for(int i = 0; i < (int)sorted_cluster_index_vector.size(); i++)
   {
      if(i >= (int)sorted_cluster_index_vector.size())
         break;
      int cl_id = sorted_cluster_index_vector[i];
      if( cluster_set[cl_id].mouth_cnt < 1 )
         tun_void_cnt++;
      else if( cluster_set[cl_id].mouth_cnt > 1 )
         tun_void_cnt++;
      else
         continue;
      if( tun_void_cnt <= output_tunnel_count )
      {
         cerr << "Tunnel/Void number " << tun_void_cnt << " is not to be curated." << endl;
         continue;
      }
      // currently we do not curate more than 100 tunnel/void.
      if( tun_void_cnt > 100 ) break;
      // curate this tunnel/void.
      for(FCI cit = triang.finite_cells_begin();
         cit != triang.finite_cells_end(); cit ++)
      {
	 if(cluster_set[cit->id].find() != cl_id ) continue;
         cit->outside = false;
         // tag this cell.
         cit->c_tag = true;
      }
   }

   int pocket_cnt = 0;
   for(int i = 0; i < (int)sorted_cluster_index_vector.size(); i++)
   {
      if(i >= (int)sorted_cluster_index_vector.size())
         break;
      int cl_id = sorted_cluster_index_vector[i];
      if( cluster_set[cl_id].mouth_cnt < 1 )
         continue;
      else if( cluster_set[cl_id].mouth_cnt > 1 )
         continue;
      else
         pocket_cnt++;

      if( pocket_cnt <= output_pocket_count )
      {
         cerr << "Pocket number " << pocket_cnt << " is not to be curated." << endl;
         continue;
      }
      // currently we do not curate more than 100 tunnel/void.
      if( pocket_cnt > 100 ) break;

      // curate this pocket.
      for(FCI cit = triang.finite_cells_begin();
         cit != triang.finite_cells_end(); cit ++)
      {
	 if(cluster_set[cit->id].find() != cl_id ) continue;
         cit->outside = false;
         // tag this cell.
         cit->c_tag = true;
      }
   }
   return;
}

// ---------------------------------------------------------
// compute_voronoi_vertex_and_cell_radius
// --------------------------------------
// Compute all the finite voronoi vertices and the circumradius
// of all the finite cells.
// ---------------------------------------------------------
void
compute_voronoi_vertex_and_cell_radius(Triangulation &triang)
{
  
   bool is_there_any_problem_in_VV_computation = false;
       
   for(FCI cit = triang.finite_cells_begin();
        cit != triang.finite_cells_end(); cit ++)
   {
           // For each cell call the circumcenter computation.
           // The function will set a boolean variable passed
           // as a parameter showing if the computation is correct.
           bool is_correct_computation = true;
           cit->set_voronoi(nondg_voronoi_point( cit->vertex(0)->point(),
                                                 cit->vertex(1)->point(),
                                                 cit->vertex(2)->point(),
                                                 cit->vertex(3)->point(),
                                                 is_correct_computation
                                               )
                           );
           is_there_any_problem_in_VV_computation |= !is_correct_computation;

           // The problem arises when the computation is not correct
           // that is, when is_correct_computation = false.

           // If is_correct_computation = true then we need to check
           // the radius. Sometimes, even if the circumcenter computation
           // is correct, radius of the cell is too big to fit in double
           // and that creates a problem. Sometimes, (possibly due to 
           // overflow) it is <= 0. In the next loop we will check for
           // the validity of radius and if that also works well we will
           // return. Otherwise, we will go the next part where the 
           // degeneracies are taken care of.
           
           if( is_correct_computation )
           {
                bool is_correct_radius = true;
                
                // check if the radius fits good in double.
          
                double r = CGAL::to_double((cit->voronoi() - cit->vertex(0)->point()) *
                                           (cit->voronoi() - cit->vertex(0)->point())
                                          );

                if( isnan(r) || isinf(r) )
                      is_correct_radius = false;
                if( r <= 0 )
                      is_correct_radius = false;

                // if it does, go back and collect the next cell
                // and do the same computation.
                if( is_correct_radius )
                      continue;
                is_there_any_problem_in_VV_computation |= !is_correct_radius;
           }

           cerr << " < bad > ";
           cit->set_dirty(true);

           // The flow comes here means either the voronoi computation 
           // is incorrect or the cell radius is junk.
           // Either case, we need to take measures.
           // Our assumption is 
           // This happens when the four points of the cell are coplanar
           // atleast that is reflected by the choice of arithmatic.
           // Our approach is to approximate the circumcenter of the cell
           // by the circumcenter of one of the triangular facets.

           cit->set_voronoi(dg_voronoi_point( cit->vertex(0)->point(),
                                              cit->vertex(1)->point(),
                                              cit->vertex(2)->point(),
                                              cit->vertex(3)->point(),
                                              is_correct_computation
                                            )
                           );

           // debug
           
           // Certain checks to make sure we are not setting anything bad
           // about the voronoi information in any cell.
           
           double cx = CGAL::to_double(cit->voronoi().x());
           double cy = CGAL::to_double(cit->voronoi().y());
           double cz = CGAL::to_double(cit->voronoi().z());

           // check 1 : the coordinates of the voronoi point fits well in double.
           CGAL_assertion( ( ( ! isnan(cx)) && ( ! isinf(cx)) ) &&
                           ( ( ! isnan(cy)) && ( ! isinf(cy)) ) &&
                           ( ( ! isnan(cz)) && ( ! isinf(cz)) ) 
                         );
           // check 2 : the radius of the cell fits well in double and is non-negative.

           double r = CGAL::to_double((cit->voronoi() - cit->vertex(0)->point()) *
                                      (cit->voronoi() - cit->vertex(0)->point())
//.........这里部分代码省略.........