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