OGS
Functions
Raster2PointCloud.cpp File Reference

Detailed Description

Copyright
Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org) Distributed under a Modified BSD License. See accompanying file LICENSE.txt or http://www.opengeosys.org/project/license

Definition in file Raster2PointCloud.cpp.

#include <algorithm>
#include <memory>
#include <string>
#include <tclap/CmdLine.h>
#include <vtkPolyDataAlgorithm.h>
#include <vtkSmartPointer.h>
#include <vtkXMLPolyDataWriter.h>
#include "Applications/DataExplorer/VtkVis/VtkGeoImageSource.h"
#include "Applications/DataExplorer/VtkVis/VtkImageDataToPointCloudFilter.h"
#include "Applications/FileIO/AsciiRasterInterface.h"
#include "BaseLib/FileTools.h"
#include "InfoLib/GitInfo.h"
Include dependency graph for Raster2PointCloud.cpp:

Go to the source code of this file.

Functions

int main (int argc, char *argv[])
 

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)

Definition at line 26 of file Raster2PointCloud.cpp.

27 {
28  TCLAP::CmdLine cmd(
29  "Batch conversion of raster files into VTK point cloud data.\n"
30  "A series of raster files is converted into point clouds by creating a "
31  "number of random points based on the intensity value in each pixel "
32  "(the higher the value, the more points. The [x,y]-extend of these "
33  "random points is limited by the extent of the pixel they represent, "
34  "the elevation is limited by the max-elevation parameter. Likewise, "
35  "the "
36  "maximum number of points per pixel is limited by the max-points "
37  "parameter. The increase of the number of points in relation to the "
38  "pixel value is indicated by gamma, where 1 is a linear increase and "
39  "values smaller or larger than 1 indicate a logarithmic or exponential "
40  "increase, respectively.\n\n"
41  "OpenGeoSys-6 software, version " +
42  GitInfoLib::GitInfo::ogs_version +
43  ".\n"
44  "Copyright (c) 2012-2021, OpenGeoSys Community "
45  "(http://www.opengeosys.org)",
46  ' ', GitInfoLib::GitInfo::ogs_version);
47  TCLAP::ValueArg<double> gamma_arg(
48  "g", "gamma",
49  "Exponental increase coefficient. A value of '1' indicates a linear "
50  "increase, smaller or larger values indicate a logarithmic or "
51  "exponential increase, respectively (default: 1).",
52  false, 0.0, "exponential increase");
53  cmd.add(gamma_arg);
54  TCLAP::ValueArg<std::size_t> max_points_arg(
55  "p", "max-points",
56  "Maximum number of points per pixel (default: 1000).", false, 0,
57  "maximum points per pixel");
58  cmd.add(max_points_arg);
59  TCLAP::ValueArg<double> max_height_arg(
60  "e", "max-elevation",
61  "Maximum elevation of randomly created points (default: 5000).", false,
62  0.0, "maximum point elevation");
63  cmd.add(max_height_arg);
64  TCLAP::ValueArg<std::string> output_arg(
65  "o", "output",
66  "Path to output polydata files (*.vtp). Specifying 'file.vtp' will "
67  "make the programme write a series of files called 'file0.vtp', "
68  "'file1.vtp', etc.",
69  true, "", "output file name");
70  cmd.add(output_arg);
71  TCLAP::ValueArg<std::string> input_arg(
72  "i", "input",
73  "Path to input raster files (*.asc). Specifying 'file.asc' will make "
74  "the programme look for a series of files called 'file0.asc', "
75  "'file1.asc', etc.",
76  true, "", "input file name");
77  cmd.add(input_arg);
78  cmd.parse(argc, argv);
79 
80  std::string const input_name = input_arg.getValue().c_str();
81  std::string const output_name = output_arg.getValue().c_str();
82 
83  double const max_height =
84  (max_height_arg.isSet()) ? max_height_arg.getValue() : 5000;
85  std::size_t const max_points =
86  (max_points_arg.isSet()) ? max_points_arg.getValue() : 1000;
87  double const gamma = (gamma_arg.isSet()) ? gamma_arg.getValue() : 1;
88 
89  std::string const ibase_name = BaseLib::dropFileExtension(input_name);
90  std::string const ifile_ext = BaseLib::getFileExtension(input_name);
91 
92  std::string const obase_name = BaseLib::dropFileExtension(output_name);
93  std::string const ofile_ext = BaseLib::getFileExtension(output_name);
94 
95  // Rainevents using point cloud filter
96  std::cout << "Starting export... ";
97  double global_min = std::numeric_limits<double>::max();
98  double global_max = 0;
99 
100  std::size_t n_rasters = 0;
101  for (std::size_t i = 0;; ++i)
102  {
103  std::string const file_name =
104  ibase_name + std::to_string(i) + ifile_ext;
105  std::unique_ptr<GeoLib::Raster> r(
106  FileIO::AsciiRasterInterface::getRasterFromASCFile(file_name));
107  if (r == nullptr)
108  {
109  n_rasters = i;
110  break;
111  }
112 
113  for (auto it = r->begin(); it != r->end(); ++it)
114  {
115  if (*it != r->getHeader().no_data)
116  {
117  global_min = std::min(*it, global_min);
118  global_max = std::max(*it, global_max);
119  }
120  }
121  }
122 
123  for (std::size_t i = 0; i < n_rasters; ++i)
124  {
125  std::string const file_name =
126  ibase_name + std::to_string(i) + ifile_ext;
127  vtkSmartPointer<VtkGeoImageSource> geo_image =
128  vtkSmartPointer<VtkGeoImageSource>::New();
129  geo_image->readImage(QString::fromStdString(file_name));
130  VtkImageDataToPointCloudFilter* const point_cloud_filter =
131  VtkImageDataToPointCloudFilter::New();
132  point_cloud_filter->SetInputConnection(geo_image->GetOutputPort());
133  geo_image->Update();
134  point_cloud_filter->SetMinValueRange(global_min);
135  point_cloud_filter->SetMaxValueRange(global_max);
136  point_cloud_filter->SetMinHeight(0);
137  point_cloud_filter->SetMaxHeight(max_height);
138  point_cloud_filter->SetIsLinear(false);
139  point_cloud_filter->SetMinNumberOfPointsPerCell(0);
140  point_cloud_filter->SetMaxNumberOfPointsPerCell(max_points);
141  point_cloud_filter->SetGamma(gamma);
142  point_cloud_filter->Update();
143  vtkSmartPointer<vtkPolyDataAlgorithm> const pd_alg =
144  dynamic_cast<vtkPolyDataAlgorithm*>(point_cloud_filter);
145  if (pd_alg)
146  {
147  std::string const output =
148  obase_name + std::to_string(i) + ofile_ext;
149  vtkSmartPointer<vtkXMLPolyDataWriter> pd_writer =
150  vtkSmartPointer<vtkXMLPolyDataWriter>::New();
151  pd_writer->SetInputData(pd_alg->GetOutputDataObject(0));
152  pd_writer->SetFileName(output.c_str());
153  pd_writer->Write();
154  }
155  }
156  std::cout << "done." << std::endl;
157  return EXIT_SUCCESS;
158 }
FileIO::AsciiRasterInterface::getRasterFromASCFile
static GeoLib::Raster * getRasterFromASCFile(std::string const &fname)
Reads an ArcGis ASC raster file.
Definition: AsciiRasterInterface.cpp:125
VtkImageDataToPointCloudFilter::New
static VtkImageDataToPointCloudFilter * New()
Create a new objects (required because of VTKs reference counting)
BaseLib::getFileExtension
std::string getFileExtension(const std::string &path)
Definition: FileTools.cpp:186
BaseLib::dropFileExtension
std::string dropFileExtension(std::string const &filename)
Definition: FileTools.cpp:169
GitInfoLib::GitInfo::ogs_version
GITINFOLIB_EXPORT const std::string ogs_version
MathLib::r
static const double r
Definition: GaussLegendreTet.cpp:64

References BaseLib::dropFileExtension(), BaseLib::getFileExtension(), FileIO::AsciiRasterInterface::getRasterFromASCFile(), VtkImageDataToPointCloudFilter::New(), GitInfoLib::GitInfo::ogs_version, and MathLib::r.