25 #include <unordered_map>
28 #include <tbb/parallel_for.h>
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>
41 namespace PdalLoader {
44 #pragma comment(linker "/INCLUDE:get_metadata_for_file")
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>();
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();
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();
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();
81 lidar_metadata->is_transformed =
true;
82 lidar_metadata->transformed_spatial_reference = pdal::SpatialReference(out_srs);
85 lidar_metadata->specified_utm_zone =
86 lidar_metadata->source_spatial_reference.getUTMZone();
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());
96 if (lidar_metadata->is_transformed) {
98 pdal::SrsTransform(lidar_metadata->source_spatial_reference,
99 lidar_metadata->transformed_spatial_reference);
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;
105 transform.transform(x_bounds.first, y_bounds.first, z_bounds.first);
107 throw std::runtime_error(
"Could not transform bounding box coordinate");
110 transform.transform(x_bounds.second, y_bounds.second, z_bounds.second);
112 throw std::runtime_error(
"Could not transform bounding box coordinate");
114 lidar_metadata->transformed_x_bounds = x_bounds;
115 lidar_metadata->transformed_y_bounds = y_bounds;
116 lidar_metadata->transformed_z_bounds = z_bounds;
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;
122 }
catch (std::exception& e) {
127 return lidar_metadata;
131 #pragma comment(linker "/INCLUDE:get_metadata_for_files")
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) {
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)) {
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);
158 lidar_metadata = get_metadata_for_file(file_path_string, out_srs);
159 metadata_cache.putDataForKey(cache_key, lidar_metadata);
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;
168 lidar_metadata_vec.resize(file_counter);
169 return std::make_pair(num_points.load(), lidar_metadata_vec);
173 #pragma comment(linker "/INCLUDE:filter_files_by_bounds")
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,
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;
187 for (
const auto& lidar_metadata : lidar_metadata_vec) {
188 unfiltered_num_points += lidar_metadata->num_points;
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) {
197 filtered_num_points += lidar_metadata->num_points;
198 filtered_lidar_metadata_vec.emplace_back(lidar_metadata);
200 return std::make_pair(filtered_num_points, filtered_lidar_metadata_vec);
204 #pragma comment(linker "/INCLUDE:keys_are_cached")
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)) {
225 #pragma comment(linker "/INCLUDE:openAndPrepareFile")
229 void PdalFileMgr::openAndPrepareFile() {
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);
239 if (!out_srs_.empty() && out_srs_ !=
"none") {
240 pdal::Options reprojection_opts;
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) {
253 las_reader_.prepare(point_table_);
254 point_view_set_ = las_reader_.execute(point_table_);
257 point_view_ = *point_view_set_.begin();
258 if (static_cast<int64_t>(las_reader_.header().pointCount()) !=
259 lidar_metadata_->num_points) {
269 #pragma comment(linker "/INCLUDE:readData")
273 void PdalFileMgr::readData(
double* x,
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) {
285 for (int64_t p = 0; p < num_points; ++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);
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);
298 point_view_->getFieldAs<int16_t>(pdal::Dimension::Id::Classification, p);
300 point_view_->getFieldAs<int8_t>(pdal::Dimension::Id::ScanAngleRank, p);
FORCE_INLINE T __attribute__((__may_alias__))*may_alias_ptr(T *ptr)
OUTPUT transform(INPUT const &input, FUNC const &func)
void parallel_for(const blocked_range< Int > &range, const Body &body, const Partitioner &p=Partitioner())
std::string filename(char const *path)
#define DEBUG_TIMER(name)