Fix/no split single tree

raylib/extraction/raysegment.cpp

@@ -166,8 +166,10 @@ std::vector<std::vector<int>> getRootsAndSegment(std::vector<Vertex> &points, co
     Eigen::ArrayXXd::Constant(static_cast<int>(lowfield.rows()), static_cast<int>(lowfield.cols()), std::numeric_limits<double>::lowest());
   for (const auto &point : points)
   {
-    Eigen::Vector3i index = ((point.pos - box_min) / pixel_width).cast<int>();
-    heightfield(index[0], index[1]) = std::max(heightfield(index[0], index[1]), point.pos[2]);
+    const Eigen::Vector3i raw = ((point.pos - box_min) / pixel_width).cast<int>();
+    const int cx = std::max(0, std::min(raw[0], static_cast<int>(heightfield.rows()) - 1));
+    const int cy = std::max(0, std::min(raw[1], static_cast<int>(heightfield.cols()) - 1));
+    heightfield(cx, cy) = std::max(heightfield(cx, cy), point.pos[2]);
   }
   // make heightfield relative to the ground
   for (int i = 0; i < heightfield.rows(); i++)
@@ -185,8 +187,10 @@ std::vector<std::vector<int>> getRootsAndSegment(std::vector<Vertex> &points, co
     points[ind].distance_to_ground = 0.0;
     points[ind].score = 0.0;
     points[ind].root = ind;
-    const Eigen::Vector3i index = ((points[ind].pos - box_min) / pixel_width).cast<int>();
-    closest_node.push(QueueNode(0, 0, heightfield(index[0], index[1]), ind, ind));
+    const Eigen::Vector3i raw = ((points[ind].pos - box_min) / pixel_width).cast<int>();
+    const int cx = std::max(0, std::min(raw[0], static_cast<int>(heightfield.rows()) - 1));
+    const int cy = std::max(0, std::min(raw[1], static_cast<int>(heightfield.cols()) - 1));
+    closest_node.push(QueueNode(0, 0, heightfield(cx, cy), ind, ind));
   }
 
   // perform Djikstra's shortest path to ground algorithm to fill in the parent indices in 'points'
@@ -204,9 +208,11 @@ std::vector<std::vector<int>> getRootsAndSegment(std::vector<Vertex> &points, co
     {
       continue;
     }
-    const Eigen::Vector3i index = ((points[point.root].pos - box_min) / pixel_width).cast<int>();
-    counts(index[0], index[1])++;
-    heights(index[0], index[1]) = std::max(heights(index[0], index[1]), point.pos[2] - lowfield(index[0], index[1]));
+    const Eigen::Vector3i raw = ((points[point.root].pos - box_min) / pixel_width).cast<int>();
+    const int cx = std::max(0, std::min(raw[0], static_cast<int>(counts.rows()) - 1));
+    const int cy = std::max(0, std::min(raw[1], static_cast<int>(counts.cols()) - 1));
+    counts(cx, cy)++;
+    heights(cx, cy) = std::max(heights(cx, cy), point.pos[2] - lowfield(cx, cy));
   }
 
   // in order to avoid boundary artefacts, we create a 2x2 summed array:
@@ -270,8 +276,10 @@ std::vector<std::vector<int>> getRootsAndSegment(std::vector<Vertex> &points, co
   std::vector<std::vector<int>> roots_lists(sums.rows() * sums.cols());
   for (int i = roots_start; i < static_cast<int>(points.size()); i++)
   {
-    const Eigen::Vector3i index = ((points[i].pos - box_min) / pixel_width).cast<int>();
-    const Eigen::Vector2i best_index = bests[index[0] + static_cast<int>(sums.rows()) * index[1]];
+    const Eigen::Vector3i raw = ((points[i].pos - box_min) / pixel_width).cast<int>();
+    const int cx = std::max(0, std::min(raw[0], static_cast<int>(sums.rows()) - 1));
+    const int cy = std::max(0, std::min(raw[1], static_cast<int>(sums.cols()) - 1));
+    const Eigen::Vector2i best_index = bests[cx + static_cast<int>(sums.rows()) * cy];
     const double max_height = max_heights(best_index[0], best_index[1]);
     if (max_height >= height_min)
     {
@@ -295,6 +303,83 @@ std::vector<std::vector<int>> getRootsAndSegment(std::vector<Vertex> &points, co
     }
   }
 
+  // Merge all root clusters into one when the scene looks like a single isolated tree
+  // that the NMS has over-segmented. Two conditions must both hold:
+  //   1. All canopy cells form a single connected region (no inter-tree gaps).
+  //   2. The total canopy area is small relative to the canopy height — consistent
+  //      with a single tree crown rather than a multi-tree stand.
+  if (roots_set.size() > 1)
+  {
+    const int rows = static_cast<int>(heightfield.rows());
+    const int cols = static_cast<int>(heightfield.cols());
+
+    // BFS connected-component labeling on cells with canopy >= height_min
+    std::vector<int> cell_label(rows * cols, -1);
+    int num_labels = 0;
+    for (int x = 0; x < rows; x++)
+    {
+      for (int y = 0; y < cols; y++)
+      {
+        if (cell_label[x + rows * y] >= 0 || heightfield(x, y) < height_min)
+        {
+          continue;
+        }
+        const int label = num_labels++;
+        std::queue<Eigen::Vector2i> q;
+        q.push(Eigen::Vector2i(x, y));
+        while (!q.empty())
+        {
+          Eigen::Vector2i c = q.front();
+          q.pop();
+          const int idx = c[0] + rows * c[1];
+          if (cell_label[idx] >= 0)
+          {
+            continue;
+          }
+          cell_label[idx] = label;
+          for (int dx = -1; dx <= 1; dx++)
+          {
+            for (int dy = -1; dy <= 1; dy++)
+            {
+              if (dx == 0 && dy == 0)
+              {
+                continue;
+              }
+              const int nx = c[0] + dx, ny = c[1] + dy;
+              if (nx >= 0 && nx < rows && ny >= 0 && ny < cols &&
+                  cell_label[nx + rows * ny] < 0 && heightfield(nx, ny) >= height_min)
+              {
+                q.push(Eigen::Vector2i(nx, ny));
+              }
+            }
+          }
+        }
+      }
+    }
+
+    // estimate how many trees the canopy area could plausibly contain:
+    // one tree crown typically occupies ~height^2 of horizontal area
+    int canopy_cells = 0;
+    for (int x = 0; x < rows; x++)
+      for (int y = 0; y < cols; y++)
+        if (heightfield(x, y) >= height_min)
+          canopy_cells++;
+    const double max_h = heightfield.maxCoeff();
+    const double expected_trees = canopy_cells * pixel_width * pixel_width / (max_h * max_h);
+
+    // merge only when there is exactly one canopy region and the total canopy
+    // footprint is consistent with a single tree (expected_trees < 1.5)
+    if (num_labels == 1 && expected_trees < 1.5)
+    {
+      std::vector<int> combined;
+      for (auto &cluster : roots_set)
+      {
+        combined.insert(combined.end(), cluster.begin(), cluster.end());
+      }
+      roots_set = {std::move(combined)};
+    }
+  }
+
   return roots_set;
 }
 

raylib/extraction/raytrees.cpp

@@ -622,6 +622,9 @@ void Trees::bifurcate(const std::vector<std::vector<int>> &clusters, double thic
       new_node.radius_scale *= max_distances[i] / std::sqrt(total_area);
       if (par == -1)
       {
+        // make this a branch of the trunk rather than a new separate root/tree,
+        // so that a merged single-tree root cluster stays as one tree
+        new_node.parent = sec_;
         new_node.radius_scale = 1.0;
       }
       if (new_node.max_distance_to_end > params_->crop_length)  // but only add if they are large enough
@@ -691,9 +694,10 @@ void Trees::bifurcate(const std::vector<std::vector<int>> &clusters, double thic
         {
           sections_[par].children.push_back(static_cast<int>(sections_.size()));
         }
-        if (par == -1)
+        else
         {
-          new_node.root = (int)sections_.size();
+          // par == -1: new_node.parent was set to sec_, so register as a child of this trunk
+          sections_[sec_].children.push_back(static_cast<int>(sections_.size()));
         }
         sections_.push_back(new_node);
       }

raylib/raylaz.cpp

@@ -4,6 +4,9 @@
 //
 // Author: Thomas Lowe
 #include "raylaz.h"
+#include <algorithm>
+#include <fstream>
+#include <limits>
 #include "raylib/rayprogress.h"
 #include "raylib/rayprogressthread.h"
 #include "rayunused.h"
@@ -262,6 +265,7 @@ LasWriter::LasWriter(const std::string &file_name)
   , writer_handle_(nullptr)
   , header_(nullptr)
   , point_(nullptr)
+  , points_written_(0)
 {
   if (laszip_create(&writer_handle_))
   {
@@ -316,6 +320,19 @@ LasWriter::~LasWriter()
     laszip_update_inventory(writer_handle_);
     laszip_close_writer(writer_handle_);
     laszip_destroy(writer_handle_);
+    // laszip_close_writer clobbers number_of_point_records for streaming writes (our count
+    // wasn't known at open_writer). Patch it on disk: LAS 1.2 legacy count at offset 107.
+    if (points_written_ > 0)
+    {
+      std::fstream f(file_name_, std::ios::in | std::ios::out | std::ios::binary);
+      if (f.is_open())
+      {
+        const laszip_U32 count =
+          static_cast<laszip_U32>(std::min<uint64_t>(points_written_, std::numeric_limits<laszip_U32>::max()));
+        f.seekp(107);
+        f.write(reinterpret_cast<const char *>(&count), sizeof(count));
+      }
+    }
   }
 #else
   std::cerr << "writeLas: cannot write file as WITHLAS not enabled. Enable using: cmake .. -DWITH_LAS=true"
@@ -345,6 +362,7 @@ bool LasWriter::writeChunk(const std::vector<Eigen::Vector3d> &points, const std
       point_->gps_time = times[i];
     laszip_write_point(writer_handle_);
   }
+  points_written_ += points.size();
   return true;
 #else   // RAYLIB_WITH_LAS
   RAYLIB_UNUSED(points);

raylib/raylaz.h

@@ -52,6 +52,7 @@ class RAYLIB_EXPORT LasWriter
   laszip_POINTER writer_handle_;
   laszip_header_struct *header_;
   laszip_point_struct *point_;
+  uint64_t points_written_;
 #endif  // RAYLIB_WITH_LAS
 };
 }  // namespace ray