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

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

Geostat_grid*
Simulacre_input_filter::read_pointset( QDataStream& stream, 
                                      std::string* errors ) {
  
  char* name;
  stream >> name;
  std::string grid_name( name );

  Q_INT32 version;
  stream >> version;
  if( version < 100 ) {
    errors->append( "file too old" );
    return 0;
  }

  Q_UINT32 size;
  stream >> size;
  
  // create a new point set of the correct size
  std::ostringstream ostr;
  ostr << size;
  std::string final_grid_name;
  SmartPtr<Named_interface> ni = 
    Root::instance()->new_interface( "point_set://" + ostr.str(),
                                     gridModels_manager + "/" + grid_name,
                                     &final_grid_name );
  Point_set* grid = dynamic_cast<Point_set*>( ni.raw_ptr() );

  // get the property names
  Q_UINT32 properties_count;
  stream >> properties_count;

  std::vector< char* > prop_names( properties_count );
  for( unsigned int i = 0; i < properties_count; i++ ) 
    stream >> prop_names[i];
  
  // read the coordinates of the points
  std::vector<Point_set::location_type > coords;
  for( unsigned int k = 0; k < size; k ++ ) {
    float x,y,z;
    stream >> x >> y >> z;
    coords.push_back( Point_set::location_type( x,y,z) );
  }
  grid->point_locations( coords );


  // read all the properties
  for( unsigned int j = 0; j < properties_count; j++ ) {
    std::string prop_name( prop_names[j] );
    GsTLGridProperty* prop = grid->add_property( prop_name );
    for( GsTLInt k = 0; k < size ; k++ ) {
      float val;
      stream >> val;
      prop->set_value( val, k );
    }
  }

  for( unsigned int l = 0; l < properties_count; l++ ) {
    delete [] prop_names[l];
  }
  delete [] name;

  return grid;
}

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

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

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