Copying Tom's ideas

£bla

Author: GabijaBern

raycloudtools/rayextract/rayextract.cpp

@@ -5,10 +5,10 @@
 // Author: Thomas Lowe
 #include "raylib/extraction/rayclusters.h"
 #include "raylib/extraction/rayforest.h"
+#include "raylib/extraction/rayleaves.h"
 #include "raylib/extraction/rayterrain.h"
 #include "raylib/extraction/raytrees.h"
 #include "raylib/extraction/raytrunks.h"
-#include "raylib/extraction/rayleaves.h"
 #include "raylib/raycloud.h"
 #include "raylib/rayforestgen.h"
 #include "raylib/rayforeststructure.h"
@@ -25,7 +25,8 @@ static std::string extract_type;
 
 void usage(int exit_code = 1)
 {
-  const bool none = extract_type != "terrain" && extract_type != "trunks" && extract_type != "forest" && extract_type != "trees" && extract_type != "leaves";
+  const bool none = extract_type != "terrain" && extract_type != "trunks" && extract_type != "forest" &&
+                    extract_type != "trees" && extract_type != "leaves";
   // clang-format off
   std::cout << "Extract natural features into a text file or mesh file" << std::endl;
   std::cout << "usage:" << std::endl;
@@ -54,21 +55,22 @@ void usage(int exit_code = 1)
     std::cout << "                            --max_diameter 0.9   - (-m) maximum trunk diameter in segmenting trees" << std::endl;
     std::cout << "                            --crop_length 1.0    - (-p) crops small branches to this distance from end" << std::endl;
     std::cout << "                            --distance_limit 1   - (-d) maximum distance between neighbour points in a tree" << std::endl;
-    std::cout << "                            --height_min 2       - (-h) minimum height counted as a tree" << std::endl;
-    std::cout << "                            --min_radius 0.01    - (-r) minimum tree radius to consider" << std::endl;
+    std::cout << "                            --height_min 2       - (-h) minimum height tree to reconstruct" << std::endl;
+    std::cout << "                            --radius_min 0       - (-r) minimum radius tree to reconstruct" << std::endl;
     std::cout << "                            --girth_height_ratio 0.12 - (-i) the amount up tree's height to estimate trunk girth" << std::endl;
     std::cout << "                            --global_taper 0.024 - (-a) force a taper value (diameter per length) for trees under global_taper_factor of max tree height. Use 0 to estimate global taper from the data" << std::endl;
     std::cout << "                            --global_taper_factor 0.3- (-o) 1 estimates same taper for whole scan, 0 is per-tree tapering. Like a soft cutoff at this amount of max tree height" << std::endl;
     std::cout << "                            --gravity_factor 0.3 - (-f) larger values preference vertical trees" << std::endl;
     std::cout << "                            --branch_segmentation- (-b) _segmented.ply is per branch segment" << std::endl;
     std::cout << "                            --grid_width 10      - (-w) crops results assuming cloud has been gridded with given width" << std::endl;
-    std::cout << "                            --grid_origin 0,0,0 - location of origin within grid cell that overlaps it. Defaults to a cell-centre origin (at grid_width/2 in each axis) matching raysplit grid. 0,0,0 is for a vertex origin." << std::endl;
     std::cout << "                            --use_rays           - (-u) use rays to reduce trunk radius overestimation in noisy cloud data" << std::endl;
-    std::cout << "                            (for internal constants -c -g -s see source file rayextract)" << std::endl;
+    std::cout << "                            (for internal constants -c -g -s -d see source file rayextract)" << std::endl;
   // These are the internal parameters that I don't expose as they are 'advanced' only, you shouldn't need to adjust them
   //  std::cout << "                            --cylinder_length_to_width 4- (-c) how slender the cylinders are" << std::endl;
   //  std::cout << "                            --gap_ratio 0.016    - (-g) will split for lateral gaps at this multiple of branch length" << std::endl;
   //  std::cout << "                            --span_ratio 4.5     - (-s) will split when branch width spans this multiple of radius" << std::endl;
+  //  std::cout << "                            --grid_origin 0,0    - (-d) location of grid corner (any of them) when grid_width used, use 0,0 for grid with vertex at 0,0. 
+  //                                                                   Default is -grid_width/2,-grid_width/2 to match the grid in raysplit grid" << std::endl;
   }
   if (extract_type == "leaves" || none)
   {
@@ -94,33 +96,33 @@ int rayExtract(int argc, char *argv[])
   }
   ray::FileArgument cloud_file, mesh_file, trunks_file, trees_file, leaf_file;
   ray::TextArgument forest("forest"), trees("trees"), trunks("trunks"), terrain("terrain"), leaves("leaves");
-
-  ray::Vector3dArgument grid_origin;
-  ray::OptionalKeyValueArgument grid_origin_option("grid_origin", &grid_origin);
   ray::OptionalKeyValueArgument groundmesh_option("ground", 'g', &mesh_file);
   ray::OptionalKeyValueArgument trunks_option("trunks", 't', &trunks_file);
   ray::DoubleArgument gradient(0.001, 1000.0, 1.0), global_taper(0.0, 1.0), global_taper_factor(0.0, 1.0);
   ray::OptionalKeyValueArgument gradient_option("gradient", 'g', &gradient);
-  ray::OptionalFlagArgument exclude_rays("exclude_rays", 'e'), segment_branches("branch_segmentation", 'b'), stalks("stalks", 's'), use_rays("use_rays", 'u');
+  ray::OptionalFlagArgument exclude_rays("exclude_rays", 'e'), segment_branches("branch_segmentation", 'b'),
+    stalks("stalks", 's'), use_rays("use_rays", 'u');
   ray::DoubleArgument width(0.01, 10.0, 0.25), drop(0.001, 1.0), max_gradient(0.01, 5.0), min_gradient(0.01, 5.0);
 
   ray::DoubleArgument max_diameter(0.01, 100.0), distance_limit(0.01, 10.0), height_min(0.01, 1000.0),
-    min_diameter(0.01, 100.0), leaf_area(0.00001, 1.0, 0.002), leaf_droop(0.0, 10.0, 0.1), crop_length(0.01, 100.0);;
-  ray::DoubleArgument girth_height_ratio(0.001, 0.5), length_to_radius(0.01, 10000.0), cylinder_length_to_width(0.1, 20.0), gap_ratio(0.01, 10.0),
-    span_ratio(0.01, 10.0), min_radius(0.01, 100.0);
-  ray::DoubleArgument gravity_factor(0.0, 100.0), grid_width(1.0, 100000.0),
-    grid_overlap(0.0, 0.9);
+    radius_min(0.0, 1000.0), min_diameter(0.01, 100.0), leaf_area(0.00001, 1.0, 0.002), leaf_droop(0.0, 10.0, 0.1),
+    crop_length(0.01, 100.0);
+  ;
+  ray::DoubleArgument girth_height_ratio(0.001, 0.5), length_to_radius(0.01, 10000.0),
+    cylinder_length_to_width(0.1, 20.0), gap_ratio(0.01, 10.0), span_ratio(0.01, 10.0);
+  ray::DoubleArgument gravity_factor(0.0, 100.0), grid_width(1.0, 100000.0), grid_overlap(0.0, 0.9);
+  ray::Vector2dArgument grid_origin(-1e10, 1e10);
   ray::OptionalKeyValueArgument max_diameter_option("max_diameter", 'm', &max_diameter);
   ray::OptionalKeyValueArgument crop_length_option("crop_length", 'n', &crop_length);
   ray::OptionalKeyValueArgument distance_limit_option("distance_limit", 'd', &distance_limit);
   ray::OptionalKeyValueArgument height_min_option("height_min", 'h', &height_min);
+  ray::OptionalKeyValueArgument radius_min_option("radius_min", 'r', &radius_min);
   ray::OptionalKeyValueArgument girth_height_ratio_option("girth_height_ratio", 'i', &girth_height_ratio);
   ray::OptionalKeyValueArgument cylinder_length_to_width_option("cylinder_length_to_width", 'c',
                                                                 &cylinder_length_to_width);
-
-  ray::OptionalKeyValueArgument min_radius_option("min_radius", 'r', &min_radius);
   ray::OptionalKeyValueArgument gap_ratio_option("gap_ratio", 'g', &gap_ratio);
   ray::OptionalKeyValueArgument span_ratio_option("span_ratio", 's', &span_ratio);
+  ray::OptionalKeyValueArgument grid_origin_option("grid_origin", 'd', &grid_origin);
   ray::OptionalKeyValueArgument gravity_factor_option("gravity_factor", 'f', &gravity_factor);
   ray::OptionalKeyValueArgument grid_width_option("grid_width", 'w', &grid_width);
   ray::OptionalKeyValueArgument global_taper_option("global_taper", 'a', &global_taper);
@@ -142,11 +144,12 @@ int rayExtract(int argc, char *argv[])
     { &groundmesh_option, &trunks_option, &width_option, &smooth_option, &drop_option, &verbose });
   bool extract_trees = ray::parseCommandLine(
     argc, argv, { &trees, &cloud_file, &mesh_file },
-    { &max_diameter_option, &distance_limit_option, &height_min_option, &crop_length_option, &girth_height_ratio_option,
-      &cylinder_length_to_width_option, &gap_ratio_option, &span_ratio_option, &gravity_factor_option,
-      &segment_branches, &grid_width_option, &global_taper_option, &global_taper_factor_option, &use_rays, &verbose,
-      &grid_origin_option, &min_radius_option });
-  bool extract_leaves = ray::parseCommandLine(argc, argv, { &leaves, &cloud_file, &trees_file }, { &leaf_option, &leaf_area_option, &leaf_droop_option, &stalks });
+    { &max_diameter_option, &distance_limit_option, &height_min_option, &radius_min_option, &crop_length_option,
+      &girth_height_ratio_option, &cylinder_length_to_width_option, &gap_ratio_option, &span_ratio_option,
+      &grid_origin_option, &gravity_factor_option, &segment_branches, &grid_width_option, &global_taper_option,
+      &global_taper_factor_option, &use_rays, &verbose });
+  bool extract_leaves = ray::parseCommandLine(argc, argv, { &leaves, &cloud_file, &trees_file },
+                                              { &leaf_option, &leaf_area_option, &leaf_droop_option, &stalks });
 
 
   if (!extract_trunks && !extract_forest && !extract_terrain && !extract_trees && !extract_leaves)
@@ -199,6 +202,10 @@ int rayExtract(int argc, char *argv[])
     {
       params.height_min = height_min.value();
     }
+    if (radius_min_option.isSet())
+    {
+      params.radius_min = radius_min.value();
+    }
     if (crop_length_option.isSet())
     {
       params.crop_length = crop_length.value();
@@ -226,6 +233,12 @@ int rayExtract(int argc, char *argv[])
     if (grid_width_option.isSet())
     {
       params.grid_width = grid_width.value();
+      params.grid_origin =
+        -Eigen::Vector2d(grid_width.value(), grid_width.value()) / 2.0;  // the centred grid used by raysplit grid
+    }
+    if (grid_origin_option.isSet())
+    {
+      params.grid_origin = grid_origin.value();
     }
     if (global_taper_option.isSet())
     {
@@ -235,14 +248,6 @@ int rayExtract(int argc, char *argv[])
     {
       params.global_taper_factor = global_taper_factor.value();
     }
-    if (grid_origin_option.isSet())
-    {
-      params.grid_origin = grid_origin.value();
-    }
-    if (min_radius_option.isSet())
-    {
-      params.min_radius = min_radius.value();
-    }
     params.use_rays = use_rays.isSet();
     params.segment_branches = segment_branches.isSet();
 
@@ -259,11 +264,15 @@ int rayExtract(int argc, char *argv[])
     ray::ForestStructure forest;
     if (!forest.load(cloud_file.nameStub() + "_trees.txt"))
     {
-      usage();
+      std::cerr << "Unable to load "
+                << cloud_file.nameStub() +
+                     "_trees.txt to generate tree mesh file, this could mean that there were no trees output"
+                << std::endl;
+      exit(true);
     }
     ray::Mesh tree_mesh;
     forest.generateSmoothMesh(tree_mesh, -1, 1, 1, 1);
-    ray::writePlyMesh(cloud_file.nameStub() + "_trees_mesh.ply", tree_mesh, true);    
+    ray::writePlyMesh(cloud_file.nameStub() + "_trees_mesh.ply", tree_mesh, true);
   }
   // extract the tree locations from a larger, aerial view of a forest
   else if (extract_forest)
@@ -321,8 +330,8 @@ int rayExtract(int argc, char *argv[])
   }
   else if (extract_leaves)
   {
-    ray::generateLeaves(cloud_file.nameStub(), trees_file.name(), leaf_file.name(), 
-      leaf_area.value(), leaf_droop.value(), stalks.isSet());
+    ray::generateLeaves(cloud_file.nameStub(), trees_file.name(), leaf_file.name(), leaf_area.value(),
+                        leaf_droop.value(), stalks.isSet());
   }
   else
   {

raylib/extraction/raytrees.cpp

@@ -7,13 +7,16 @@
 #include <nabo/nabo.h>
 #include "rayclusters.h"
 
+#include <unordered_set>
+
 namespace ray
 {
 TreesParams::TreesParams()
   : max_diameter(0.9)
   , crop_length(1.0)
   , distance_limit(1.0)
   , height_min(2.0)
+  , radius_min(0.0)
   , girth_height_ratio(0.12)
   , cylinder_length_to_width(4.0)
   , gap_ratio(0.016)
@@ -23,7 +26,7 @@ TreesParams::TreesParams()
   , segment_branches(false)
   , global_taper(0.012)
   , global_taper_factor(0.3)
-  , grid_origin(0, 0, 0)
+  , grid_origin(0, 0)
 {}
 
 /// The main reconstruction algorithm
@@ -148,7 +151,7 @@ Trees::Trees(Cloud &cloud, const Eigen::Vector3d &offset, const Mesh &mesh, cons
       continue;
     }
 
-    if (params_->min_radius && estimated_radius < params_->min_radius)
+    if (params_->radius_min && estimated_radius < params_->radius_min)
     {
       sections_[sec_].tip = base + Eigen::Vector3d(0, 0, 0.01);
       sections_[sec_].total_weight = 1e-10;
@@ -300,9 +303,9 @@ Trees::Trees(Cloud &cloud, const Eigen::Vector3d &offset, const Mesh &mesh, cons
     removeOutOfBoundSections(cloud, min_bound, max_bound, offset);
   }
 
-  if (params_->min_radius)
+  if (params_->radius_min)
   {
-    removeSmallRadiusTrees();
+    removeSmallRadiusTrees(verbose, offset, children);
   }
 
   std::vector<int> root_segs(cloud.ends.size(), -1);
@@ -1082,14 +1085,19 @@ void Trees::removeOutOfBoundSections(const Cloud &cloud, Eigen::Vector3d &min_bo
 {
   const double width = params_->grid_width;
   cloud.calcBounds(&min_bound, &max_bound);
-  const Eigen::Vector3d mid = (min_bound + max_bound) / 2.0 + offset + params_->grid_origin;
-  const Eigen::Vector2d inds(std::round(mid[0] / width), std::round(mid[1] / width));
-  min_bound[0] = width * (inds[0] - 0.5) - offset[0];
-  min_bound[1] = width * (inds[1] - 0.5) - offset[1];
-  max_bound[0] = width * (inds[0] + 0.5) - offset[0];
-  max_bound[1] = width * (inds[1] + 0.5) - offset[1];
-  std::cerr << "min bound: " << (min_bound + offset).transpose() << ", max bound: " << (max_bound + offset).transpose()
-            << " offset " << offset << std::endl;
+  const Eigen::Vector3d mid = (min_bound + max_bound) / 2.0 + offset;
+  Eigen::Vector2d mid_local(mid[0], mid[1]);
+  mid_local -= params_->grid_origin;
+  const Eigen::Vector2d inds(std::floor(mid_local[0] / width), std::floor(mid_local[1] / width));
+  min_bound[0] = width * (inds[0]) + params_->grid_origin[0] - offset[0];
+  min_bound[1] = width * (inds[1]) + params_->grid_origin[1] - offset[1];
+  max_bound[0] = width * (inds[0] + 1.0) + params_->grid_origin[0] - offset[0];
+  max_bound[1] = width * (inds[1] + 1.0) + params_->grid_origin[1] - offset[1];
+
+  std::cout << "min bound: " << (min_bound + offset).transpose() << ", max bound: " << (max_bound + offset).transpose()
+            << std::endl;
+
+  int number_of_trees_outside = 0;
 
   // disable trees out of bounds
   for (auto &section : sections_)
@@ -1102,8 +1110,10 @@ void Trees::removeOutOfBoundSections(const Cloud &cloud, Eigen::Vector3d &min_bo
     if (pos[0] < min_bound[0] || pos[0] > max_bound[0] || pos[1] < min_bound[1] || pos[1] > max_bound[1])
     {
       section.children.clear();  // make it a non-tree
+      number_of_trees_outside++;
     }
   }
+  std::cerr << "Number of trees outside of section: " << number_of_trees_outside << std::endl;
 }
 
 // colour the cloud by tree id, or by branch segment id
@@ -1269,11 +1279,19 @@ bool Trees::save(const std::string &filename, const Eigen::Vector3d &offset, boo
   return true;
 }
 
-void Trees::removeSmallRadiusTrees()
+void Trees::removeSmallRadiusTrees(bool verbose, const Eigen::Vector3d &offset,
+                                   const std::vector<std::vector<int>> &children)
 {
+  std::cerr << "Number of section " << sections_.size() << std::endl;
+  ray::Cloud debug_cloud;
   for (sec_ = 0; sec_ < (int)sections_.size(); sec_++)
   {
-    std::cerr << "------- processing " << sec_ << "th section" << std::endl;
+    if (sections_[sec_].parent >= 0 || sections_[sec_].children.empty())
+    {
+      continue;
+    }
+    if (verbose)
+      std::cerr << "------- processing " << sec_ << "th section" << std::endl;
     double best_accuracy = -1.0;
     std::vector<int> nodes;                    // used
     Eigen::Vector3d base = getRootPosition();  // used
@@ -1301,7 +1319,7 @@ void Trees::removeSmallRadiusTrees()
       if (verbose)
       {
         // choose random color for the section
-
+        ray::RGBA section_color(0, 0, 0, 255);
         for (auto &node : nodes)
         {
           debug_cloud.addRay(Eigen::Vector3d(0, 0, 0), points_[node].pos, 0.0, section_color);
@@ -1324,9 +1342,11 @@ void Trees::removeSmallRadiusTrees()
       double circle_coverage = 0;
       Eigen::Vector3d circle_centre;
       double radius = estimateCylinderRadiusUsingRANSAC(nodes, up, accuracy, circle_coverage, circle_centre);
-
-      std::cerr << sec_ << "Estimated radius  " << radius << " and accuracy " << accuracy << " and coverage "
-                << circle_coverage << " j " << j << std::endl;
+      if (verbose)
+      {
+        std::cerr << sec_ << "Estimated radius  " << radius << " and accuracy " << accuracy << " and coverage "
+                  << circle_coverage << " j " << j << std::endl;
+      }
       if (accuracy > 0.4 && circle_coverage > 0.5)
       {
         radius_estimates.push_back(radius);
@@ -1351,20 +1371,27 @@ void Trees::removeSmallRadiusTrees()
         std::accumulate(radius_estimates.begin(), radius_estimates.end(), 0.0) / radius_estimates.size();
       double accuracy_average =
         std::accumulate(accuracy_estimates.begin(), accuracy_estimates.end(), 0.0) / accuracy_estimates.size();
-
-      std::cerr << sec_ << "radius_average  " << radius_average << std::endl;
+      if (verbose)
+        std::cerr << sec_ << "radius_average  " << radius_average << std::endl;
       if (radius_average < 0.05 && accuracy_average > 0.4)
       {
-        std::cerr << "Removing tree with average radius " << radius_average << " sec " << sec_ << std::endl;
+        if (verbose)
+          std::cerr << "Removing tree with average radius " << radius_average << " sec " << sec_ << std::endl;
         sections_[sec_].children.clear();
         continue;
       }
     }
   }
+  if (verbose)
+  {
+    debug_cloud.translate(offset);
+    debug_cloud.save("../data/debug.ply");
+  }
 }
 
 double Trees::estimateCylinderRadiusUsingRANSAC(const std::vector<int> &nodes, const Eigen::Vector3d &dir,
-                                                double &accuracy, double &circle_coverage)
+                                                double &accuracy, double &circle_coverage,
+                                                Eigen::Vector3d &circle_centre)
 {
   double best_rad = 0.0;
   double best_accuracy = 0.0;
@@ -1456,6 +1483,7 @@ double Trees::estimateCylinderRadiusUsingRANSAC(const std::vector<int> &nodes, c
       best_rad = avg_radius;
       best_accuracy = static_cast<double>(inliers.size()) / nodes.size();
       circle_coverage = coverage;
+      circle_centre = centroid;
     }
   }