C++ Point_set::size方法代码示例

  void scale_space (const Point_set& points, ItemsInserter items,
		    unsigned int scale, bool generate_smooth = false,
                    bool separate_shells = false, bool force_manifold = true,
                    unsigned int samples = 300, unsigned int iterations = 4)
  {
    ScaleSpace reconstruct (scale, samples);
    reconstruct.reconstruct_surface(points.begin (), points.end (), iterations,
                                    separate_shells, force_manifold);

    for( unsigned int sh = 0; sh < reconstruct.number_of_shells(); ++sh )
      {
        Scene_polygon_soup_item* new_item
          = new Scene_polygon_soup_item ();
        new_item->setColor(Qt::magenta);
        new_item->setRenderingMode(FlatPlusEdges);
        new_item->init_polygon_soup(points.size(), reconstruct.number_of_triangles ());

        Scene_polygon_soup_item* smooth_item = NULL;
        if (generate_smooth)
          {
            smooth_item = new Scene_polygon_soup_item ();
            smooth_item->setColor(Qt::magenta);
            smooth_item->setRenderingMode(FlatPlusEdges);
            smooth_item->init_polygon_soup(points.size(), reconstruct.number_of_triangles ());
          }

        std::map<unsigned int, unsigned int> map_i2i;
        unsigned int current_index = 0;
    
        for (ScaleSpace::Triple_iterator it = reconstruct.shell_begin (sh);
             it != reconstruct.shell_end (sh); ++ it)
          {
            for (unsigned int ind = 0; ind < 3; ++ ind)
              {
                if (map_i2i.find ((*it)[ind]) == map_i2i.end ())
                  {
                    map_i2i.insert (std::make_pair ((*it)[ind], current_index ++));
                    Point p = points[(*it)[ind]].position();
                    new_item->new_vertex (p.x (), p.y (), p.z ());
                    
                    if (generate_smooth)
                      {
                        p = *(reconstruct.points_begin() + (*it)[ind]);
                        smooth_item->new_vertex (p.x (), p.y (), p.z ());
                      }
                  }
              }
            new_item->new_triangle( map_i2i[(*it)[0]],
                                    map_i2i[(*it)[1]],
                                    map_i2i[(*it)[2]] );
            if (generate_smooth)
              smooth_item->new_triangle( map_i2i[(*it)[0]],
                                         map_i2i[(*it)[1]],
                                         map_i2i[(*it)[2]] );
              
          }

        *(items ++) = new_item;
        if (generate_smooth)
          *(items ++) = smooth_item;
      }

    if (force_manifold)
      {
        std::ptrdiff_t num = std::distance( reconstruct.garbage_begin(  ),
                                            reconstruct.garbage_end(  ) );

        Scene_polygon_soup_item* new_item
          = new Scene_polygon_soup_item ();
        new_item->setColor(Qt::blue);
        new_item->setRenderingMode(FlatPlusEdges);
        new_item->init_polygon_soup(points.size(), num);

        Scene_polygon_soup_item* smooth_item = NULL;
        if (generate_smooth)
          {
            smooth_item = new Scene_polygon_soup_item ();
            smooth_item->setColor(Qt::blue);
            smooth_item->setRenderingMode(FlatPlusEdges);
            smooth_item->init_polygon_soup(points.size(), num);
          }

        std::map<unsigned int, unsigned int> map_i2i;

        unsigned int current_index = 0;
        for (ScaleSpace::Triple_iterator it=reconstruct.garbage_begin(),
               end=reconstruct.garbage_end();it!=end;++it)
          {
            for (unsigned int ind = 0; ind < 3; ++ ind)
              {
                if (map_i2i.find ((*it)[ind]) == map_i2i.end ())
                  {
                    map_i2i.insert (std::make_pair ((*it)[ind], current_index ++));
                    Point p = points[(*it)[ind]].position();
                    new_item->new_vertex (p.x (), p.y (), p.z ());
                    
                    if (generate_smooth)
                      {
                        p = *(reconstruct.points_begin() + (*it)[ind]);
                        smooth_item->new_vertex (p.x (), p.y (), p.z ());
//.........这里部分代码省略.........

  unsigned int scale_of_noise (const Point_set& points, double& size)
  {
    Tree tree(points.begin(), points.end());
    
    double ratio_kept = (points.size() < 1000)
      ? 1. : 1000. / (points.size());
    
    std::vector<Point> subset;
    for (std::size_t i = 0; i < points.size (); ++ i)
      if (rand() / (double)RAND_MAX < ratio_kept)
    	subset.push_back (points[i]);
    
    std::vector<unsigned int> scales;
    generate_scales (std::back_inserter (scales));

    std::vector<unsigned int> chosen;
    
    for (std::size_t i = 0; i < subset.size (); ++ i)
      {
    	Neighbor_search search(tree, subset[i],scales.back());
	double current = 0.;
    	unsigned int nb = 0;
    	std::size_t index = 0;
	double minimum = (std::numeric_limits<double>::max)();
	unsigned int c = 0;
	
    	for (Search_iterator search_iterator = search.begin();
    	     search_iterator != search.end (); ++ search_iterator, ++ nb)
    	  {
	    current += search_iterator->second;

    	    if (nb + 1 == scales[index])
    	      {
		double score = std::sqrt (current / scales[index])
		  / std::pow (scales[index], 0.375); // NB ^ (5/12)

		if (score < minimum)
		  {
		    minimum = score;
		    c = scales[index];
		  }

    		++ index;
    		if (index == scales.size ())
    		  break;
    	      }
    	  }
	chosen.push_back (c);
      }

    std::sort (chosen.begin (), chosen.end());

    unsigned int noise_scale = chosen[chosen.size() / 2];
    size = 0.;
    for (std::size_t i = 0; i < subset.size (); ++ i)
      {
    	Neighbor_search search(tree, subset[i], noise_scale);
	size += std::sqrt ((-- search.end())->second);
      }
    size /= subset.size();

    
    return noise_scale;
  }

int main (int, char**)
{
  Point_set pts;

  pts.add_normal_map();
  
  bool map_added = false;
  Size_t_map echo_map;
  Color_map color_map;

  boost::tie (echo_map, map_added) = pts.add_property_map<std::size_t> ("echo");
  CGAL_assertion (map_added);
  boost::tie (color_map, map_added) = pts.add_property_map<Classification::RGB_Color> ("color");
  CGAL_assertion (map_added);

  for (std::size_t i = 0; i < 1000; ++ i)
  {
    Point_set::iterator it
      = pts.insert (Point (CGAL::get_default_random().get_double(),
                           CGAL::get_default_random().get_double(),
                           CGAL::get_default_random().get_double()),
                    Vector (CGAL::get_default_random().get_double(),
                            CGAL::get_default_random().get_double(),
                            CGAL::get_default_random().get_double()));
    echo_map[*it] = std::size_t(CGAL::get_default_random().get_int(0, 4));
    color_map[*it] = CGAL::make_array ((unsigned char)(CGAL::get_default_random().get_int(0, 255)),
                                       (unsigned char)(CGAL::get_default_random().get_int(0, 255)),
                                       (unsigned char)(CGAL::get_default_random().get_int(0, 255)));
  }

  Feature_set features;
  Feature_generator generator (features, pts, pts.point_map(),
                               5,  // using 5 scales
                               pts.normal_map(),
                               color_map, echo_map);

  CGAL_assertion (generator.number_of_scales() == 5);
  CGAL_assertion (features.size() == 80);

  Label_set labels;

  std::vector<int> training_set (pts.size(), -1);
  for (std::size_t i = 0; i < 20; ++ i)
  {
    std::ostringstream oss;
    oss << "label_" << i;
    Label_handle lh = labels.add(oss.str().c_str());

    for (std::size_t j = 0; j < 10; ++ j)
      training_set[std::size_t(CGAL::get_default_random().get_int(0, int(training_set.size())))] = int(i);
  }
  CGAL_assertion (labels.size() == 20);
  
  Classifier classifier (labels, features);
  
  classifier.train<CGAL::Sequential_tag> (training_set, 800);
#ifdef CGAL_LINKED_WITH_TBB
  classifier.train<CGAL::Parallel_tag> (training_set, 800);
#endif

  std::vector<int> label_indices(pts.size(), -1);
  
  Classification::classify<CGAL::Sequential_tag>
    (pts, labels, classifier, label_indices);

  Classification::classify_with_local_smoothing<CGAL::Sequential_tag>
    (pts, pts.point_map(), labels, classifier,
     generator.neighborhood().sphere_neighbor_query(0.01f),
     label_indices);

  Classification::classify_with_graphcut<CGAL::Sequential_tag>
    (pts, pts.point_map(), labels, classifier,
     generator.neighborhood().k_neighbor_query(12),
     0.2f, 10, label_indices);

#ifdef CGAL_LINKED_WITH_TBB
  Classification::classify<CGAL::Sequential_tag>
    (pts, labels, classifier, label_indices);

  Classification::classify_with_local_smoothing<CGAL::Sequential_tag>
    (pts, pts.point_map(), labels, classifier,
     generator.neighborhood().sphere_neighbor_query(0.01f),
     label_indices);

  Classification::classify_with_graphcut<CGAL::Sequential_tag>
    (pts, pts.point_map(), labels, classifier,
     generator.neighborhood().k_neighbor_query(12),
     0.2f, 10, label_indices);
#endif

  Classification::Evaluation evaluation (labels, training_set, label_indices);
  
  return EXIT_SUCCESS;
}

  unsigned int scale_of_anisotropy (const Point_set& points, double& size)
  {
    Tree tree(points.begin(), points.end());
    
    double ratio_kept = (points.size() < 1000)
      ? 1. : 1000. / (points.size());
    
    std::vector<Point> subset;
    for (std::size_t i = 0; i < points.size (); ++ i)
      if (rand() / (double)RAND_MAX < ratio_kept)
    	subset.push_back (points[i]);
    
    std::vector<unsigned int> scales;
    generate_scales (std::back_inserter (scales));

    std::vector<unsigned int> chosen;

    for (std::size_t i = 0; i < subset.size (); ++ i)
      {
    	Neighbor_search search(tree, subset[i],scales.back());
	double current = 0.;
    	unsigned int nb = 0;
    	std::size_t index = 0;
	double maximum = 0.;
	unsigned int c = 0;
	
    	for (Search_iterator search_iterator = search.begin();
    	     search_iterator != search.end (); ++ search_iterator, ++ nb)
    	  {
	    current += search_iterator->second;

    	    if (nb + 1 == scales[index])
    	      {
		double score = std::sqrt (current / scales[index])
		  / std::pow (scales[index], 0.75); // NB ^ (3/4)

		if (score > maximum)
		  {
		    maximum = score;
		    c = scales[index];
		  }

    		++ index;
    		if (index == scales.size ())
    		  break;
    	      }
    	  }
	chosen.push_back (c);
      }

    double mean = 0.;
    for (std::size_t i = 0; i < chosen.size(); ++ i)
      mean += chosen[i];
    mean /= chosen.size();
    unsigned int aniso_scale = static_cast<unsigned int>(mean);

    size = 0.;
    for (std::size_t i = 0; i < subset.size (); ++ i)
      {
    	Neighbor_search search(tree, subset[i], aniso_scale);
	size += std::sqrt ((-- search.end())->second);
      }
    size /= subset.size();
    
    return aniso_scale;
  }

bool poisson_reconstruct(FaceGraph* graph,
                         Point_set& points,
                         typename Traits::FT sm_angle, // Min triangle angle (degrees).
                         typename Traits::FT sm_radius, // Max triangle size w.r.t. point set average spacing.
                         typename Traits::FT sm_distance, // Approximation error w.r.t. point set average spacing.
                         const QString& solver_name, // solver name
                         bool use_two_passes,
                         bool do_not_fill_holes)
{
  // Poisson implicit function
  typedef CGAL::Poisson_reconstruction_function<Traits> Poisson_reconstruction_function;

  // Surface mesher
  typedef CGAL::Surface_mesh_default_triangulation_3 STr;
  typedef CGAL::Surface_mesh_complex_2_in_triangulation_3<STr> C2t3;
  typedef CGAL::Implicit_surface_3<Traits, Poisson_reconstruction_function> Surface_3;

  // AABB tree
  typedef CGAL::AABB_face_graph_triangle_primitive<FaceGraph> Primitive;
  typedef CGAL::AABB_traits<Traits, Primitive> AABB_traits;
  typedef CGAL::AABB_tree<AABB_traits> AABB_tree;
    CGAL::Timer task_timer; task_timer.start();

    //***************************************
    // Checks requirements
    //***************************************

    if (points.size() == 0)
    {
      std::cerr << "Error: empty point set" << std::endl;
     return false;
    }

    bool points_have_normals = points.has_normal_map();
    if ( ! points_have_normals )
    {
      std::cerr << "Input point set not supported: this reconstruction method requires oriented normals" << std::endl;
     return false;
    }

    CGAL::Timer reconstruction_timer; reconstruction_timer.start();

    //***************************************
    // Computes implicit function
    //***************************************

 
    std::cerr << "Computes Poisson implicit function "
              << "using " << solver_name.toLatin1().data() << " solver...\n";
              
    
    // Creates implicit function from the point set.
    // Note: this method requires an iterator over points
    // + property maps to access each point's position and normal.
    Poisson_reconstruction_function function(points.begin_or_selection_begin(), points.end(),
                                             points.point_map(), points.normal_map());

    bool ok = false;    
    #ifdef CGAL_EIGEN3_ENABLED
    if(solver_name=="Eigen - built-in simplicial LDLt")
    {
      CGAL::Eigen_solver_traits<Eigen::SimplicialCholesky<CGAL::Eigen_sparse_matrix<double>::EigenType> > solver;
      ok = function.compute_implicit_function(solver, use_two_passes);
    }
    if(solver_name=="Eigen - built-in CG")
    {
      CGAL::Eigen_solver_traits<Eigen::ConjugateGradient<CGAL::Eigen_sparse_matrix<double>::EigenType> > solver;
      solver.solver().setTolerance(1e-6);
      solver.solver().setMaxIterations(1000);
      ok = function.compute_implicit_function(solver, use_two_passes);
    }
    #endif

    // Computes the Poisson indicator function f()
    // at each vertex of the triangulation.
    if ( ! ok )
    {
      std::cerr << "Error: cannot compute implicit function" << std::endl;
     return false;
    }

    // Prints status
    std::cerr << "Total implicit function (triangulation+refinement+solver): " << task_timer.time() << " seconds\n";
    task_timer.reset();

    //***************************************
    // Surface mesh generation
    //***************************************

    std::cerr << "Surface meshing...\n";

    // Computes average spacing
    Kernel::FT average_spacing = CGAL::compute_average_spacing<Concurrency_tag>(points.all_or_selection_if_not_empty(),
                                                                                6 /* knn = 1 ring */,
                                                                                points.parameters());

    // Gets one point inside the implicit surface
    Kernel::Point_3 inner_point = function.get_inner_point();
    Kernel::FT inner_point_value = function(inner_point);
    if(inner_point_value >= 0.0)
//.........这里部分代码省略.........