Spatial Indexing

apps/CMakeLists.txt

@@ -19,8 +19,23 @@ add_executable(${LAS2LAS++_EXE_NAME} las2las++.cpp)
 
 target_link_libraries(${LAS2LAS++_EXE_NAME} ${LIBRARY_NAME} OpenMP::OpenMP_CXX)
 
+set(INSPECT_VLRS_EXE_NAME las++-inspect-vlrs)
+
+add_executable(${INSPECT_VLRS_EXE_NAME} inspect_vlrs.cpp)
+
+target_link_libraries(${INSPECT_VLRS_EXE_NAME} ${LIBRARY_NAME}
+                      OpenMP::OpenMP_CXX)
+
+set(VALIDATE_SPATIAL_INDEX_EXE_NAME las++-validate-spatial-index)
+
+add_executable(${VALIDATE_SPATIAL_INDEX_EXE_NAME} validate_spatial_index.cpp)
+
+target_link_libraries(${VALIDATE_SPATIAL_INDEX_EXE_NAME} ${LIBRARY_NAME}
+                      OpenMP::OpenMP_CXX)
+
 install(
-  TARGETS ${LAS2LAS++_EXE_NAME}
+  TARGETS ${LAS2LAS++_EXE_NAME} ${INSPECT_VLRS_EXE_NAME}
+          ${VALIDATE_SPATIAL_INDEX_EXE_NAME}
   DESTINATION bin
   COMPONENT applications)
 

apps/inspect_vlrs.cpp

@@ -0,0 +1,103 @@
+/*
+ * SPDX-FileCopyrightText: (c) 2025 Trailblaze Software, all rights reserved
+ * SPDX-License-Identifier: LGPL-2.1-or-later
+ *
+ * This library is free software; you can redistribute it and/or modify it under
+ * the terms of the GNU Lesser General Public License as published by the Free
+ * Software Foundation; version 2.1.
+ *
+ * This library is distributed in the hope that it will be useful, but WITHOUT
+ * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ * FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+ * details.
+ *
+ * For LGPL2 incompatible licensing or development requests, please contact
+ * trailblaze.software@gmail.com
+ */
+
+#include <cmath>
+#include <fstream>
+#include <iomanip>
+#include <iostream>
+#include <string>
+
+#include "las_reader.hpp"
+#include "spatial_index.hpp"
+#include "utilities/assert.hpp"
+#include "utilities/printing.hpp"
+#include "vlr.hpp"
+
+using namespace laspp;
+
+int main(int argc, char* argv[]) {
+  if (argc != 2) {
+    std::cerr << "Usage: " << argv[0] << " <laz_file>" << std::endl;
+    return 1;
+  }
+
+  std::string filename = argv[1];
+  std::filesystem::path file_path(filename);
+  LASReader reader(file_path);
+
+  std::cout << "=== VLRs ===" << std::endl;
+  const auto& vlrs = reader.vlr_headers();
+  for (size_t i = 0; i < vlrs.size(); ++i) {
+    const auto& vlr = vlrs[i];
+    std::cout << "VLR " << i << ":" << std::endl;
+    std::cout << indented(vlr, "  ");
+    std::cout << "  Is LAZ VLR: " << (vlr.is_laz_vlr() ? "yes" : "no") << std::endl;
+    std::cout << "  Is LAStools spatial index VLR: "
+              << (vlr.is_lastools_spatial_index_vlr() ? "yes" : "no") << std::endl;
+    std::cout << std::endl;
+  }
+
+  std::cout << "=== EVLRs ===" << std::endl;
+  const auto& evlrs = reader.evlr_headers();
+  for (size_t i = 0; i < evlrs.size(); ++i) {
+    const auto& evlr = evlrs[i];
+    std::cout << "EVLR " << i << ":" << std::endl;
+    std::cout << indented(evlr, "  ");
+    std::cout << "  Is LAZ VLR: " << (evlr.is_laz_vlr() ? "yes" : "no") << std::endl;
+    std::cout << "  Is LAStools spatial index EVLR: "
+              << (evlr.is_lastools_spatial_index_evlr() ? "yes" : "no") << std::endl;
+
+    // If it's a LAStools spatial index, show more details
+    if (evlr.is_lastools_spatial_index_evlr()) {
+      LASPP_ASSERT(reader.has_lastools_spatial_index(),
+                   "Failed to load LAStools spatial index from EVLR.");
+      std::cout << "Spatial index: " << std::endl;
+      const auto& index = reader.lastools_spatial_index();
+      const auto& quadtree = index.quadtree_header();
+
+      std::cout << "  Spatial index header:" << std::endl;
+      std::cout << indented(quadtree, "    ");
+      std::cout << "  Number of cells: " << index.num_cells() << std::endl;
+
+      // Show cell details
+      size_t total_intervals = 0;
+      for (const auto& [cell_index, cell] : index.cells()) {
+        total_intervals += cell.intervals.size();
+      }
+      std::cout << "  Total intervals: " << total_intervals << std::endl;
+      std::cout << "  Average intervals per cell: " << std::fixed << std::setprecision(1)
+                << (index.num_cells() > 0 ? static_cast<double>(total_intervals) /
+                                                static_cast<double>(index.num_cells())
+                                          : 0.0)
+                << std::endl;
+
+      // Show first few cells as examples
+      std::cout << indented(limited_map(index.cells(), 3), "    ") << std::endl;
+    }
+    std::cout << std::endl;
+  }
+
+  if (reader.has_lastools_spatial_index()) {
+    const auto& spatial_index = reader.lastools_spatial_index();
+    std::cout << "Loaded LAStools spatial index with " << spatial_index.num_cells() << " cells."
+              << std::endl;
+    std::cout << "Spatial index details:" << std::endl;
+    std::cout << indented(spatial_index, "  ") << std::endl;
+  }
+
+  return 0;
+}

apps/las2las++.cpp

@@ -15,89 +15,43 @@
  * trailblaze.software@gmail.com
  */
 
-#include <cstring>
 #include <filesystem>
 #include <fstream>
 #include <iostream>
+#include <string>
 
-#include "las_point.hpp"
 #include "las_reader.hpp"
 #include "las_writer.hpp"
+#include "utilities/assert.hpp"
 
-class LASPoint {
-  std::array<int32_t, 3> position;
-  double gps_time;
-
- public:
-  void operator=(const laspp::LASPointFormat0& point) {
-    position[0] = point.x;
-    position[1] = point.y;
-    position[2] = point.z;
-  }
-
-  void operator=(const laspp::GPSTime& point) { gps_time = point; }
-
-  operator laspp::LASPointFormat0() const {
-    laspp::LASPointFormat0 point;
-    point.x = position[0];
-    point.y = position[1];
-    point.z = position[2];
-    point.intensity = 0;
-    point.bit_byte = 0;
-    point.classification_byte = 0;
-    point.scan_angle_rank = 0;
-    point.user_data = 0;
-    point.point_source_id = 0;
-    return point;
-  }
-
-  operator laspp::GPSTime() const { return laspp::GPSTime(gps_time); }
+int main(int argc, char* argv[]) {
+  bool add_spatial_index_flag = false;
+  int file_arg_start = 1;
 
-  friend std::ostream& operator<<(std::ostream& os, const LASPoint& point) {
-    os << "Position: (" << point.position[0] << ", " << point.position[1] << ", "
-       << point.position[2] << ")";
-    os << " GPS Time: " << point.gps_time;
-    return os;
+  for (int i = 1; i < argc; ++i) {
+    if (std::string(argv[i]) == "--add-spatial-index" || std::string(argv[i]) == "-s") {
+      add_spatial_index_flag = true;
+      file_arg_start++;
+    }
   }
-};
-
-template <typename PointType>
-void read_and_write_points(laspp::LASReader& reader, laspp::LASWriter& writer) {
-  std::vector<PointType> points(reader.num_points());
-  reader.read_chunks<PointType>(points, {0, reader.num_chunks()});
-  std::cout << points[0] << std::endl;
-  std::cout << points[points.size() - 1] << std::endl;
-  writer.write_points<PointType>(points);
-}
 
-int main(int argc, char* argv[]) {
-  if (argc != 3) {
-    std::cerr << "Usage: " << argv[0] << " <in_file> <out_file>" << std::endl;
+  if (argc - file_arg_start != 2) {
+    std::cerr << "Usage: " << argv[0] << " [--add-spatial-index|-s] <in_file> <out_file>"
+              << std::endl;
+    std::cerr
+        << "  --add-spatial-index, -s: Add spatial index to output file (points will be reordered)"
+        << std::endl;
     return 1;
   }
 
-  std::filesystem::path in_file(argv[1]);
-  std::ifstream ifs(in_file, std::ios::binary);
-  LASPP_ASSERT(ifs.is_open(), "Failed to open ", in_file);
-  laspp::LASReader reader(ifs);
-  std::cout << reader.header() << std::endl;
-
-  if (reader.geo_keys()) {
-    std::cout << "GeoTIFF Projection Info:" << std::endl;
-    std::cout << *reader.geo_keys() << std::endl;
-  }
-
-  auto vlrs = reader.vlr_headers();
-  for (const auto& vlr : vlrs) {
-    std::cout << vlr << std::endl;
-  }
+  std::string in_file_str = argv[file_arg_start];
+  std::string out_file_str = argv[file_arg_start + 1];
 
-  auto evlrs = reader.evlr_headers();
-  for (const auto& evlr : evlrs) {
-    std::cout << evlr << std::endl;
-  }
+  std::filesystem::path in_file(in_file_str);
+  laspp::LASReader reader(in_file);
+  std::cout << reader.header() << std::endl;
 
-  std::filesystem::path out_file(argv[2]);
+  std::filesystem::path out_file(out_file_str);
   std::fstream ofs;
   ofs.open(out_file, std::ios::binary | std::ios::in | std::ios::out | std::ios::trunc);
   LASPP_ASSERT(ofs.is_open(), "Failed to open ", out_file);
@@ -121,28 +75,32 @@ int main(int argc, char* argv[]) {
     }
     laspp::LASWriter writer(ofs, point_format);
 
-    for (const auto& vlr : reader.vlr_headers()) {
-      if (vlr.is_laz_vlr()) {
-        continue;
-      }
-      writer.write_vlr(vlr, reader.read_vlr_data(vlr));
-    }
-
-    writer.header().transform() = reader.header().transform();
-
-    {
-      std::vector<LASPoint> points(reader.num_points());
-
-      // reader.read_chunks<LASPoint>(points, {0, reader.num_chunks()});
-      // std::cout << points[0] << std::endl;
-      // std::cout << points[1] << std::endl;
-      // std::cout << points[1000] << std::endl;
-      // std::cout << points[points.size() - 1] << std::endl;
-
-      // writer.write_points<LASPoint>(points);
-      LASPP_SWITCH_OVER_POINT_TYPE(reader.header().point_format(), read_and_write_points, reader,
-                                   writer);
-    }
+    // Copy everything from reader to writer
+    writer.copy_from_reader(reader, add_spatial_index_flag);
+
+    LASPP_ASSERT_EQ(reader.num_points(), writer.header().num_points(),
+                    "Number of points in output file does not match input file");
+    LASPP_ASSERT_EQ(reader.header().num_points_by_return(), writer.header().num_points_by_return(),
+                    "Number of points by return in output file does not match input file");
+    LASPP_ASSERT_LE(reader.header().bounds().min_x(), writer.header().bounds().min_x(),
+                    "Output file min x is less than input file min x");
+    LASPP_ASSERT_GE(reader.header().bounds().max_x(), writer.header().bounds().max_x(),
+                    "Output file max x is greater than input file max x");
+    LASPP_ASSERT_LE(reader.header().bounds().min_y(), writer.header().bounds().min_y(),
+                    "Output file min y is less than input file min y");
+    LASPP_ASSERT_GE(reader.header().bounds().max_y(), writer.header().bounds().max_y(),
+                    "Output file max y is greater than input file max y");
+    LASPP_ASSERT_LE(reader.header().bounds().min_z(), writer.header().bounds().min_z(),
+                    "Output file min z is less than input file min z");
+    LASPP_ASSERT_GE(reader.header().bounds().max_z(), writer.header().bounds().max_z(),
+                    "Output file max z is greater than input file max z");
+    LASPP_ASSERT_EQ(reader.header().transform().scale_factors(),
+                    writer.header().transform().scale_factors(),
+                    "Scale factors in output file do not match input file");
+    LASPP_ASSERT_EQ(reader.header().transform().offsets(), writer.header().transform().offsets(),
+                    "Offsets in output file do not match input file");
+
+    std::cout << writer.header() << std::endl;
   }
 
   return 0;