C++ Point_set类代码示例

  void advancing_front (const Point_set& points, Scene_polyhedron_item* new_item, double size,
                        double radius_ratio_bound = 5., double beta = 0.52)
  {

    // TODO: build DT with indices
    
    Polyhedron& P = * const_cast<Polyhedron*>(new_item->polyhedron());
    Radius filter (size);

    typedef CGAL::Advancing_front_surface_reconstruction_vertex_base_3<Kernel> LVb;
    typedef CGAL::Advancing_front_surface_reconstruction_cell_base_3<Kernel> LCb;

    typedef CGAL::Triangulation_data_structure_3<LVb,LCb> Tds;
    typedef CGAL::Delaunay_triangulation_3<Kernel,Tds> Triangulation_3;

    typedef CGAL::Advancing_front_surface_reconstruction<Triangulation_3, Radius> Reconstruction;

    Triangulation_3 dt( boost::make_transform_iterator(points.begin_or_selection_begin(), Point_set_make_pair_point_index(points)),
                        boost::make_transform_iterator(points.end(), Point_set_make_pair_point_index(points)) );

    Reconstruction R(dt, filter);
    R.run(radius_ratio_bound, beta);
    CGAL::AFSR::construct_polyhedron(P, R);
						  
  }

  void advancing_front (const Point_set& points, Scene_polyhedron_item* new_item, double size)
  {
    Polyhedron& P = * const_cast<Polyhedron*>(new_item->polyhedron());
    Radius filter(10 * size);

    CGAL::advancing_front_surface_reconstruction (points.begin (), points.end (), P, filter);
						  
  }

void print_point_set (const Point_set& point_set)
{
  std::cerr << "Content of point set:" << std::endl;
  for (Point_set::const_iterator it = point_set.begin();
       it != point_set.end(); ++ it)
    std::cerr << "* Point " << *it
              << ": " << point_set.point(*it) // or point_set[it]
              << " with normal " << point_set.normal(*it)
              << std::endl;
}

  void advancing_front (const Point_set& points, Scene_polyhedron_item* new_item, double size,
                        double radius_ratio_bound = 5., double beta = 0.52)
  {
    Polyhedron& P = * const_cast<Polyhedron*>(new_item->polyhedron());
    Radius filter (size);

    CGAL::advancing_front_surface_reconstruction (points.begin (), points.end (), P, filter,
                                                  radius_ratio_bound, beta);
						  
  }

  void compute_normals (Point_set& points, unsigned int neighbors)
  {
    CGAL::jet_estimate_normals<Concurrency_tag>(points.begin_or_selection_begin(), points.end(),
                                                points.point_map(),
                                                points.normal_map(),
                                                2 * neighbors);

    points.set_first_selected (CGAL::mst_orient_normals (points.begin(), points.end(),
                                                         points.point_map(),
                                                         points.normal_map(),
                                                         2 * neighbors));
    points.delete_selection();
  }

Geostat_grid* Pointset_geometry_xml_io::
  read_grid_geometry(QDir dir,const QDomElement& elem, std::string* errors) const {

	std::string grid_name = elem.attribute("name").toStdString();

	QDomElement elemGeom = elem.firstChildElement("Geometry");
	QString grid_size_str = elemGeom.attribute("size");

	bool ok_size;
	int size = grid_size_str.toInt(&ok_size);
	if(!ok_size) {
		errors->append("Failed to get the size of the grid");
		return 0;
	}


  SmartPtr<Named_interface> ni =
    Root::instance()->new_interface( "point_set://" + grid_name+"::"+grid_size_str.toStdString(),
                                     gridModels_manager + "/" + grid_name);
  Point_set* grid = dynamic_cast<Point_set*>( ni.raw_ptr() );
  if(grid==0) return 0;

  QFile file( dir.absolutePath()+"/coordinates.sgems" );
  if( !file.open( QIODevice::ReadOnly ) ) {
  	errors->append("Could not open file coordinates.sgems");
  	return grid;
  }

	QDataStream stream( &file );
#if QT_VERSION >= 0x040600
	stream.setFloatingPointPrecision(QDataStream::SinglePrecision);
#endif
 /*
  // read the coordinates of the points
  std::string coordsfile = dir.absolutePath().toStdString()+"/coordinates.sgems";
  std::fstream stream(coordsfile.c_str(),std::ios::in | std::ios::binary );
  if(stream.bad()) {
  	errors->append("Could not open file "+coordsfile);
  	return grid;
  }
*/
  std::vector<Point_set::location_type > coords;
  for( unsigned int k = 0; k < size; k ++ ) {
  	GsTLCoord x,y,z;
    stream >> x >> y >> z;
    coords.push_back( Point_set::location_type( x,y,z) );
  }
  grid->point_locations( coords );

  return grid;
}

  void compute_normals (Point_set& points, unsigned int neighbors)
  {
    CGAL::jet_estimate_normals<Concurrency_tag>(points.begin(), points.end(),
                                                CGAL::make_normal_of_point_with_normal_pmap(Point_set::value_type()),
                                                2 * neighbors);
    
    points.erase (CGAL::mst_orient_normals (points.begin(), points.end(),
					    CGAL::make_normal_of_point_with_normal_pmap(Point_set::value_type()),
					    2 * neighbors),
		  points.end ());
  }

 void smooth_point_set (Point_set& points, unsigned int scale)
 {
   CGAL::jet_smooth_point_set<Concurrency_tag>(points.begin(), points.end(),
                                               scale);
 }

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

  void simplify_point_set (Point_set& points, double size)
  {
    points.erase (CGAL::grid_simplify_point_set (points.begin (), points.end (), size),
		  points.end ());
  }

  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)
//.........这里部分代码省略.........

Geostat_grid* Csv_poinset_infilter::read( std::ifstream& infile ) {

	QByteArray tmp = dialog_->name().simplified().toLatin1();
	std::string name( tmp.constData() );

  int X_col_id = dialog_->X_column_index();
  int Y_col_id = dialog_->Y_column_index();
  int Z_col_id = dialog_->Z_column_index();

  if( X_col_id == Y_col_id || X_col_id == Z_col_id || Y_col_id == Z_col_id ) {
    GsTLcerr << "The same column was selected for multiple coordinates \n" << gstlIO::end;
    return 0;
  }

  bool use_no_data_value = dialog_->use_no_data_value();
  float no_data_value = GsTLGridProperty::no_data_value;
  if( dialog_->use_no_data_value() ) {
    no_data_value = dialog_->no_data_value();
  }


  std::string str;
  std::getline(infile, str);
  QString qstr(str.c_str());
  QStringList property_names = qstr.split(",");

  bool is_x_provided = dialog_->X_column_name() != "None";
  bool is_y_provided = dialog_->Y_column_name() != "None";
  bool is_z_provided = dialog_->Z_column_name() != "None";
  if(is_x_provided) property_names.removeOne(dialog_->X_column_name());
  if(is_y_provided) property_names.removeOne(dialog_->Y_column_name());
  if(is_z_provided) property_names.removeOne(dialog_->Z_column_name());


  std::vector< std::vector< QString > > property_values( property_names.size() );
  std::vector< Point_set::location_type > point_locations;

//  For a csv file no data value is indicated by an empty field e.g. {34,,5.5}
  while( std::getline(infile, str) ) {
    qstr = str.c_str();
    QStringList fields = qstr.split(",");

    Point_set::location_type loc;
    if(is_x_provided) loc[0] = fields[X_col_id].toDouble();
    if(is_y_provided) loc[1] = fields[Y_col_id].toDouble();
    if(is_z_provided) loc[2] = fields[Z_col_id].toDouble();
    point_locations.push_back(loc);

    unsigned int i=0;
   for(unsigned int j=0;j<fields.size();j++) {
      if(j==0) i=0;
      if(j != X_col_id && j != Y_col_id && j != Z_col_id) {
        property_values[i].push_back(fields[j]);
        i++;
      }
    }
  }
  
//   done reading file
//----------------------------

  int point_set_size = point_locations.size();
  appli_message( "read " << point_set_size << " points" );

  // We now have a vector containing all the locations and another one with
  // all the property values.
  // Create a pointset, initialize it with the data we collected, and we're done
 
  // ask manager to get a new pointset and initialize it
  SmartPtr<Named_interface> ni =
    Root::instance()->interface( gridModels_manager + "/" + name );

  if( ni.raw_ptr() != 0 ) {
    GsTLcerr << "object " << name << " already exists\n" << gstlIO::end;
    return 0;
  }

  std::string size_str = String_Op::to_string( point_set_size );
  ni = Root::instance()->new_interface( "point_set://" + size_str, 
				  gridModels_manager + "/" + name );
  Point_set* pset = dynamic_cast<Point_set*>( ni.raw_ptr() );
  appli_assert( pset != 0 );

  pset->point_locations( point_locations );


  for( unsigned int k= 0; k < property_names.size(); k++ ) {
  // Need to find out if property is categorical
    unsigned int check_size = std::min(30,static_cast<int>(property_values[k].size()));
    bool is_categ = false;
    for(unsigned int i=0; i<check_size ; i++ ) {
      bool ok;
      property_values[k][i].toFloat(&ok);
      if(!ok)  {
        is_categ = true;
        break;
      }
    }

    if(!is_categ) {
//.........这里部分代码省略.........

int main() {

  const int grid_size = 10;
  // Build grid with locally varying mean
  Cartesian_grid* lvm_grid = new Cartesian_grid( grid_size, grid_size, 1 );
  GsTLGridProperty* prop = lvm_grid->add_property( "trend", typeid( float ) );
  
  for( int i=0; i< grid_size*grid_size/2 ; i++ ) {
    prop->set_value( 0.0, i );
  }
  for( int i=grid_size*grid_size/2; i< grid_size*grid_size ; i++ ) {
    prop->set_value( 10.0, i );
  }
  Colocated_neighborhood* coloc_neigh = 
    dynamic_cast<Colocated_neighborhood*>( 
	  lvm_grid->colocated_neighborhood( "trend" )
       );


  // Build kriging grid 
  Cartesian_grid* krig_grid = new Cartesian_grid( grid_size, grid_size, 1 );
  GsTLGridProperty* krig_prop = 
    krig_grid->add_property( string("krig"), typeid( float )  );
  krig_grid->select_property( "krig");

  // Build harddata grid
  const int pointset_size = 4;
  Point_set* harddata = new Point_set( pointset_size );
  std::vector<GsTLPoint> locations;
  locations.push_back( GsTLPoint( 0,0,0 ) );
  locations.push_back( GsTLPoint( 1,5,0 ) );
  locations.push_back( GsTLPoint( 8,8,0 ) );
  locations.push_back( GsTLPoint( 5,2,0 ) );
  harddata->point_locations( locations );
  GsTLGridProperty* hard_prop = harddata->add_property( "poro" );
  
  for( int i=0; i<pointset_size; i++ ) {
    hard_prop->set_value( i, i );
  }

  harddata->select_property( "poro" );


  // Set up covariance
  Covariance<GsTLPoint> cov;
  cov.nugget(0.1);
  cov.add_structure( "Spherical" );
  cov.sill( 0, 0.9 );
  cov.set_geometry( 0, 10,10,10, 0,0,0 );

  Grid_initializer initializer;
  initializer.assign( krig_grid,
		      krig_prop,
		      harddata,
		      "poro" );
	
  for( int i=0; i< krig_prop->size(); i++ ) {
    if( krig_prop->is_harddata( i ) )
      cout << "value at " << i << ": " 
	   << krig_prop->get_value( i ) << endl;
  }

  krig_grid->select_property( "krig");

  Neighborhood* neighborhood = krig_grid->neighborhood( 20, 20, 1, 0,0,0,
							&cov, true );

  typedef GsTLPoint Location;
  typedef std::vector<double>::const_iterator weight_iterator;
  typedef SKConstraints_impl< Neighborhood, Location > SKConstraints;
  typedef SK_local_mean_combiner<weight_iterator, Neighborhood,
                                 Colocated_neighborhood> LVM_combiner;
  typedef Kriging_constraints< Neighborhood, Location > KrigingConstraints;
  typedef Kriging_combiner< weight_iterator, Neighborhood > KrigingCombiner;


  LVM_combiner combiner( *coloc_neigh );
  SKConstraints constraints;


  // initialize the algo
  Kriging algo;
  algo.simul_grid_ = krig_grid;
  algo.property_name_ = "krig";
  algo.harddata_grid_ = 0;
  algo.neighborhood_ = neighborhood;
  algo.covar_ = cov;
  algo.combiner_ = new KrigingCombiner( &combiner );
  algo.Kconstraints_ = new KrigingConstraints( &constraints );


  // Run and output the results
  algo.execute();

  ofstream ofile( "result.out" );
  if( !ofile ) {
    cerr << "can't create file result.out" << endl;
    return 1;
  }

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