OmniSciDB  b28c0d5765
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PdalLoader.cpp
Go to the documentation of this file.
1 /*
2  * Copyright 2022 HEAVY.AI, Inc.
3  *
4  * Licensed under the Apache License, Version 2.0 (the "License");
5  * you may not use this file except in compliance with the License.
6  * You may obtain a copy of the License at
7  *
8  * http://www.apache.org/licenses/LICENSE-2.0
9  *
10  * Unless required by applicable law or agreed to in writing, software
11  * distributed under the License is distributed on an "AS IS" BASIS,
12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13  * See the License for the specific language governing permissions and
14  * limitations under the License.
15  */
16 
17 #ifdef HAVE_PDAL
18 #ifndef __CUDACC__
19 
20 #include <algorithm>
21 #include <cstring>
22 #include <filesystem>
23 #include <memory>
24 #include <string>
25 #include <unordered_map>
26 #include <vector>
27 
28 #include <tbb/parallel_for.h>
29 
30 #include <pdal/Options.hpp>
31 #include <pdal/PointTable.hpp>
32 #include <pdal/PointView.hpp>
33 #include <pdal/SpatialReference.hpp>
34 #include <pdal/filters/ReprojectionFilter.hpp>
35 #include <pdal/io/LasHeader.hpp>
36 #include <pdal/io/LasReader.hpp>
37 #include <pdal/private/SrsTransform.hpp>
38 
39 #include "PdalLoader.h"
40 
41 namespace PdalLoader {
42 
43 #ifdef _WIN32
44 #pragma comment(linker "/INCLUDE:get_metadata_for_file")
45 #else
46 __attribute__((__used__))
47 #endif
48 std::shared_ptr<LidarMetadata> get_metadata_for_file(const std::string& filename,
49  const std::string& out_srs) {
50  std::shared_ptr<LidarMetadata> lidar_metadata = std::make_shared<LidarMetadata>();
51 
52  try {
53  pdal::Options las_opts;
54  pdal::Option las_opt("filename", filename);
55  las_opts.add(las_opt);
56  pdal::PointTable table;
57  pdal::LasReader las_reader;
58  las_reader.setOptions(las_opts);
59  las_reader.prepare(table);
60  pdal::LasHeader las_header = las_reader.header();
61  lidar_metadata->filename = filename;
62  lidar_metadata->file_source_id = las_header.fileSourceId();
63  lidar_metadata->version_major = las_header.versionMajor();
64  lidar_metadata->version_minor = las_header.versionMinor();
65  lidar_metadata->creation_year = las_header.creationYear();
66  lidar_metadata->is_compressed = las_header.compressed();
67  lidar_metadata->num_points = las_header.pointCount();
68  const auto& dims = las_header.usedDims();
69  lidar_metadata->num_dims = dims.size();
70  lidar_metadata->has_time = las_header.hasTime();
71  lidar_metadata->has_color = las_header.hasColor();
72  lidar_metadata->has_wave = las_header.hasWave();
73  lidar_metadata->has_infrared = las_header.hasInfrared();
74  lidar_metadata->has_14_point_format = las_header.has14PointFormat();
75 
76  lidar_metadata->source_spatial_reference = las_header.srs();
77  if (out_srs.empty()) {
78  lidar_metadata->is_transformed = false;
79  lidar_metadata->transformed_spatial_reference = las_header.srs();
80  } else {
81  lidar_metadata->is_transformed = true;
82  lidar_metadata->transformed_spatial_reference = pdal::SpatialReference(out_srs);
83  }
84 
85  lidar_metadata->specified_utm_zone =
86  lidar_metadata->source_spatial_reference.getUTMZone();
87 
88  // todo: transform these to out_srs
89  lidar_metadata->source_x_bounds =
90  std::make_pair(las_header.minX(), las_header.maxX());
91  lidar_metadata->source_y_bounds =
92  std::make_pair(las_header.minY(), las_header.maxY());
93  lidar_metadata->source_z_bounds =
94  std::make_pair(las_header.minZ(), las_header.maxZ());
95 
96  if (lidar_metadata->is_transformed) {
97  pdal::SrsTransform transform =
98  pdal::SrsTransform(lidar_metadata->source_spatial_reference,
99  lidar_metadata->transformed_spatial_reference);
100  // Need to copy bounds as transform changes them in place
101  std::pair<double, double> x_bounds = lidar_metadata->source_x_bounds;
102  std::pair<double, double> y_bounds = lidar_metadata->source_y_bounds;
103  std::pair<double, double> z_bounds = lidar_metadata->source_z_bounds;
104  bool transform_ok =
105  transform.transform(x_bounds.first, y_bounds.first, z_bounds.first);
106  if (!transform_ok) {
107  throw std::runtime_error("Could not transform bounding box coordinate");
108  }
109  transform_ok =
110  transform.transform(x_bounds.second, y_bounds.second, z_bounds.second);
111  if (!transform_ok) {
112  throw std::runtime_error("Could not transform bounding box coordinate");
113  }
114  lidar_metadata->transformed_x_bounds = x_bounds;
115  lidar_metadata->transformed_y_bounds = y_bounds;
116  lidar_metadata->transformed_z_bounds = z_bounds;
117  } else {
118  lidar_metadata->transformed_x_bounds = lidar_metadata->source_x_bounds;
119  lidar_metadata->transformed_y_bounds = lidar_metadata->source_y_bounds;
120  lidar_metadata->transformed_z_bounds = lidar_metadata->source_z_bounds;
121  }
122  } catch (std::exception& e) {
123  // std::cout << "Could not read from file: " << filename << " with error: " <<
124  // e.what()
125  // << std::endl;
126  }
127  return lidar_metadata;
128 }
129 
130 #ifdef _WIN32
131 #pragma comment(linker "/INCLUDE:get_metadata_for_files")
132 #else
133 __attribute__((__used__))
134 #endif
135 std::pair<int64_t, std::vector<std::shared_ptr<LidarMetadata>>> get_metadata_for_files(
136  const std::vector<std::filesystem::path>& file_paths,
137  const std::string& out_srs) {
138  auto timer = DEBUG_TIMER(__func__);
139  const size_t num_files(file_paths.size());
140  std::vector<std::shared_ptr<LidarMetadata>> lidar_metadata_vec(num_files);
141  std::atomic<uint32_t> file_counter = 0;
142  std::atomic<int64_t> num_points = 0;
144  tbb::blocked_range<size_t>(0, num_files), [&](const tbb::blocked_range<size_t>& r) {
145  const size_t start_file_idx = r.begin();
146  const size_t end_file_idx = r.end();
147  for (size_t file_idx = start_file_idx; file_idx != end_file_idx; ++file_idx) {
148  const auto& file_path_string = file_paths[file_idx].string();
149  std::string cache_key = make_cache_key(file_path_string, out_srs, "metadata");
150  std::shared_ptr<LidarMetadata> lidar_metadata;
151  if (metadata_cache.isKeyCached(cache_key)) {
152  try {
153  lidar_metadata = metadata_cache.getDataForKey(cache_key);
154  } catch (std::exception& e) {
155  lidar_metadata = get_metadata_for_file(file_path_string, out_srs);
156  }
157  } else {
158  lidar_metadata = get_metadata_for_file(file_path_string, out_srs);
159  metadata_cache.putDataForKey(cache_key, lidar_metadata);
160  }
161  if (lidar_metadata->num_points > 0) {
162  const auto out_idx = file_counter.fetch_add(1, std::memory_order_relaxed);
163  num_points.fetch_add(lidar_metadata->num_points, std::memory_order_relaxed);
164  lidar_metadata_vec[out_idx] = lidar_metadata;
165  }
166  }
167  });
168  lidar_metadata_vec.resize(file_counter);
169  return std::make_pair(num_points.load(), lidar_metadata_vec);
170 }
171 
172 #ifdef _WIN32
173 #pragma comment(linker "/INCLUDE:filter_files_by_bounds")
174 #else
175 __attribute__((__used__))
176 #endif
177 std::pair<int64_t, std::vector<std::shared_ptr<LidarMetadata>>> filter_files_by_bounds(
178  const std::vector<std::shared_ptr<LidarMetadata>>& lidar_metadata_vec,
179  const double x_min,
180  const double x_max,
181  const double y_min,
182  const double y_max) {
183  std::vector<std::shared_ptr<LidarMetadata>> filtered_lidar_metadata_vec;
184  int64_t unfiltered_num_points = 0;
185  int64_t filtered_num_points = 0;
186 
187  for (const auto& lidar_metadata : lidar_metadata_vec) {
188  unfiltered_num_points += lidar_metadata->num_points;
189 
190  if (x_max < lidar_metadata->transformed_x_bounds.first ||
191  x_min > lidar_metadata->transformed_x_bounds.second ||
192  y_max < lidar_metadata->transformed_y_bounds.first ||
193  y_min > lidar_metadata->transformed_y_bounds.second) {
194  // doesn't pass metadata filter - discard
195  continue;
196  }
197  filtered_num_points += lidar_metadata->num_points;
198  filtered_lidar_metadata_vec.emplace_back(lidar_metadata);
199  }
200  return std::make_pair(filtered_num_points, filtered_lidar_metadata_vec);
201 }
202 
203 #ifdef _WIN32
204 #pragma comment(linker "/INCLUDE:keys_are_cached")
205 #else
206 __attribute__((__used__))
207 #endif
208 bool keys_are_cached(
209  const std::string& filename,
210  const std::string& out_srs,
211  const size_t num_records,
212  const std::vector<std::pair<std::string, size_t>>& sub_keys_and_byte_sizes) {
213  for (const auto& sub_key_and_byte_size : sub_keys_and_byte_sizes) {
214  const std::string cache_key(
215  make_cache_key(filename, out_srs, sub_key_and_byte_size.first));
216  const size_t num_bytes = num_records * sub_key_and_byte_size.second;
217  if (!data_cache.isKeyCachedAndSameLength(cache_key, num_bytes)) {
218  return false;
219  }
220  }
221  return true;
222 }
223 
224 #ifdef _WIN32
225 #pragma comment(linker "/INCLUDE:openAndPrepareFile")
226 #else
227 __attribute__((__used__))
228 #endif
229 void PdalFileMgr::openAndPrepareFile() {
230  auto timer = DEBUG_TIMER(__func__);
231  if (is_ready_) {
232  return;
233  }
234  pdal::Option las_reader_filename_opt("filename", lidar_metadata_->filename);
235  pdal::Options las_reader_opts;
236  las_reader_opts.add(las_reader_filename_opt);
237  las_reader_.setOptions(las_reader_opts);
238 
239  if (!out_srs_.empty() && out_srs_ != "none") {
240  pdal::Options reprojection_opts;
241  try {
242  reprojection_opts.add("out_srs",
243  lidar_metadata_->transformed_spatial_reference.getWKT());
244  pdal::ReprojectionFilter reprojection_filter;
245  reprojection_filter.setOptions(reprojection_opts);
246  reprojection_filter.setInput(las_reader_);
247  reprojection_filter.prepare(point_table_);
248  point_view_set_ = reprojection_filter.execute(point_table_);
249  } catch (pdal::pdal_error& e) {
250  // std::cout << "Encountered error in reprojection: " << e.what() << std::endl;
251  }
252  } else {
253  las_reader_.prepare(point_table_);
254  point_view_set_ = las_reader_.execute(point_table_);
255  }
256 
257  point_view_ = *point_view_set_.begin();
258  if (static_cast<int64_t>(las_reader_.header().pointCount()) !=
259  lidar_metadata_->num_points) {
260  // std::cout << "Encountered mismatch between file header point count and "
261  // "metadata."
262  // << std::endl;
263  return;
264  }
265  is_ready_ = true;
266 }
267 
268 #ifdef _WIN32
269 #pragma comment(linker "/INCLUDE:readData")
270 #else
271 __attribute__((__used__))
272 #endif
273 void PdalFileMgr::readData(double* x,
274  double* y,
275  double* z,
276  int32_t* intensity,
277  int8_t* return_num,
278  int8_t* num_returns,
279  int8_t* scan_direction_flag,
280  int8_t* edge_of_flight_line_flag,
281  int16_t* classification,
282  int8_t* scan_angle_rank,
283  const int64_t num_points) {
284  auto timer = DEBUG_TIMER(__func__);
285  for (int64_t p = 0; p < num_points; ++p) {
286  x[p] = point_view_->getFieldAs<double>(pdal::Dimension::Id::X, p);
287  y[p] = point_view_->getFieldAs<double>(pdal::Dimension::Id::Y, p);
288  z[p] = point_view_->getFieldAs<double>(pdal::Dimension::Id::Z, p);
289  intensity[p] = point_view_->getFieldAs<int32_t>(pdal::Dimension::Id::Intensity, p);
290  return_num[p] = point_view_->getFieldAs<int8_t>(pdal::Dimension::Id::ReturnNumber, p);
291  num_returns[p] =
292  point_view_->getFieldAs<int8_t>(pdal::Dimension::Id::NumberOfReturns, p);
293  scan_direction_flag[p] =
294  point_view_->getFieldAs<int8_t>(pdal::Dimension::Id::ScanDirectionFlag, p);
295  edge_of_flight_line_flag[p] =
296  point_view_->getFieldAs<int8_t>(pdal::Dimension::Id::EdgeOfFlightLine, p);
297  classification[p] =
298  point_view_->getFieldAs<int16_t>(pdal::Dimension::Id::Classification, p);
299  scan_angle_rank[p] =
300  point_view_->getFieldAs<int8_t>(pdal::Dimension::Id::ScanAngleRank, p);
301  }
302 }
303 
304 } // namespace PdalLoader
305 
306 #endif // HAVE_PDAL
307 #endif // __CUDACC__
FORCE_INLINE T __attribute__((__may_alias__))*may_alias_ptr(T *ptr)
Definition: TypePunning.h:32
OUTPUT transform(INPUT const &input, FUNC const &func)
Definition: misc.h:320
void parallel_for(const blocked_range< Int > &range, const Body &body, const Partitioner &p=Partitioner())
#define DEBUG_TIMER(name)
Definition: Logger.h:374