diff --git a/raycloudtools/rayextract/rayextract.cpp b/raycloudtools/rayextract/rayextract.cpp index 5f0d25c..010c2cd 100644 --- a/raycloudtools/rayextract/rayextract.cpp +++ b/raycloudtools/rayextract/rayextract.cpp @@ -54,7 +54,8 @@ 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 << " --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; @@ -62,11 +63,13 @@ void usage(int exit_code = 1) 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 << " --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 -j 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 - (-j) 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) { @@ -99,21 +102,24 @@ int rayExtract(int argc, char *argv[]) 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), + ray::DoubleArgument max_diameter(0.01, 100.0), distance_limit(0.01, 10.0), height_min(0.01, 1000.0), 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 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", 'j', &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); @@ -135,8 +141,8 @@ 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, + { &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 }); @@ -191,6 +197,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(); @@ -218,7 +228,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()) { params.global_taper = 0.5 * global_taper.value(); @@ -243,7 +258,8 @@ 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); diff --git a/raylib/extraction/raytrees.cpp b/raylib/extraction/raytrees.cpp index 99233ce..84dc4a9 100644 --- a/raylib/extraction/raytrees.cpp +++ b/raylib/extraction/raytrees.cpp @@ -14,6 +14,7 @@ TreesParams::TreesParams() , 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,6 +24,7 @@ TreesParams::TreesParams() , segment_branches(false) , global_taper(0.012) , global_taper_factor(0.3) + , grid_origin(0,0) {} /// The main reconstruction algorithm @@ -283,6 +285,10 @@ Trees::Trees(Cloud &cloud, const Eigen::Vector3d &offset, const Mesh &mesh, cons removeOutOfBoundSections(cloud, min_bound, max_bound, offset); } + if (params_->radius_min > 0.0) + { + removeSmallRadiusTrees(); + } std::vector root_segs(cloud.ends.size(), -1); // now colour the ray cloud based on the segmentation segmentCloud(cloud, root_segs, section_ids); @@ -1051,11 +1057,13 @@ 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; - 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]; + 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; // disable trees out of bounds @@ -1073,6 +1081,22 @@ void Trees::removeOutOfBoundSections(const Cloud &cloud, Eigen::Vector3d &min_bo } } +// filtering out these smaller trees must be done at the end, as that is when we have a final estimation of global taper (and therefore radius) +void Trees::removeSmallRadiusTrees() +{ + for (auto §ion : sections_) + { + if (section.parent >= 0 || section.children.empty()) + { + continue; + } + if (radius(section) < params_->radius_min) + { + section.children.clear(); // make it a non-tree + } + } +} + // colour the cloud by tree id, or by branch segment id void Trees::segmentCloud(Cloud &cloud, std::vector &root_segs, const std::vector §ion_ids) { diff --git a/raylib/extraction/raytrees.h b/raylib/extraction/raytrees.h index 523101e..621a5be 100644 --- a/raylib/extraction/raytrees.h +++ b/raylib/extraction/raytrees.h @@ -22,16 +22,18 @@ struct RAYLIB_EXPORT TreesParams double crop_length; // distance to end of branch where it crops the branch, and doesn't generate further geometry double distance_limit; // maximum distance between points that can count as connected double height_min; // minimum height for a tree. Lower values are considered undergrowth and excluded + double radius_min; // minimum radius for a tree. Note that all trees above height_min are processed before being filtered through this double girth_height_ratio; // how far up tree to measure girth double cylinder_length_to_width; // the slenderness of the branch segment cylinders double gap_ratio; // max gap per branch length double span_ratio; // points that span a larger width determine that a branch has become two - double gravity_factor; // preferences branches that are less lateral, so penalises implausable horizontal branches - double grid_width; // used on a grid cell with overlap, to remove trees with a base in the overlap zone - bool segment_branches; // flag to output the ray cloud coloured by branch segment index rather than by tree index - double global_taper; // forced global taper, uses global_taper_factor to define how much it is applied - double global_taper_factor; // 0 estimates per-tree tapering, 1 uses per-scan tapering, 0.5 is mid-way on mid-weight trees - bool use_rays; // use the full rays in order to estimate a smaller radius when points are not all on the real branch + double gravity_factor; // preferences branches that are less lateral, so penalises implausable horizontal branches + double grid_width; // used on a grid cell with overlap, to remove trees with a base in the overlap zone + bool segment_branches; // flag to output the ray cloud coloured by branch segment index rather than by tree index + double global_taper; // forced global taper, uses global_taper_factor to define how much it is applied + double global_taper_factor; // 0 estimates per-tree tapering, 1 uses per-scan tapering, 0.5 is mid-way on mid-weight trees + bool use_rays; // use the full rays in order to estimate a smaller radius when points are not all on the real branch + Eigen::Vector2d grid_origin; // used when grid_width is set, for cropping purposes }; struct BranchSection; // forwards declaration @@ -99,6 +101,8 @@ class RAYLIB_EXPORT Trees double radius(const BranchSection §ion) const; /// remove elements of nodes that are too distant to the set of end points bool removeDistantPoints(std::vector &nodes); + /// remove trees with radius under the specified minimum + void removeSmallRadiusTrees(); // cached data that is used throughout the processing method int sec_; diff --git a/raylib/raylaz.cpp b/raylib/raylaz.cpp index 46fa9e6..130f61a 100644 --- a/raylib/raylaz.cpp +++ b/raylib/raylaz.cpp @@ -29,7 +29,7 @@ bool readLas(const std::string &file_name, if (ifs.fail()) { - std::cerr << "readLas: failed to open stream" << std::endl; + std::cerr << "readLas: failed to open file " << file_name << ", error: " << std::strerror(errno) << std::endl; return false; }