Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Adjustments to https://github.com/GabijaBern/raycloudtools/tree/new_command_line_params #24

Open
wants to merge 3 commits into
base: main
Choose a base branch
from
Open

Adjustments to https://github.com/GabijaBern/raycloudtools/tree/new_command_line_params #24

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
ff1f85f
first commit of changes, now testing
thomas-lowe Nov 15, 2024
33da1c0
better error reporting on no trees output
thomas-lowe Nov 15, 2024
ac44c1e
changed single letter for grid origin param
thomas-lowe Nov 22, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
28 changes: 22 additions & 6 deletions raycloudtools/rayextract/rayextract.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -54,19 +54,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 << " --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 << " --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)
{
Expand Down Expand Up @@ -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", '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);
Expand All @@ -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 });

Expand Down Expand Up @@ -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();
Expand Down Expand Up @@ -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();
Expand All @@ -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);
Expand Down
35 changes: 30 additions & 5 deletions raylib/extraction/raytrees.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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
Expand Down Expand Up @@ -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)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am curious to understand the code a little bit better. If I run this approach, will it break the algorithm? What does bifurcate() method do?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

removeSmallRadiusTrees() is called towards the end when the radius() function will be accurate. We can't filter based on estimated_radius because this is not the final estimation, just an initial estimate.

bifurcate() is a function for splitting a segment at a tree branching point, so it adds to segments (cylinders) to the list, the key line being: sections_.push_back(new_node); Usually it is just a split in two, but it may also split into multiple branches.

{
removeSmallRadiusTrees();
}
std::vector<int> root_segs(cloud.ends.size(), -1);
// now colour the ray cloud based on the segmentation
segmentCloud(cloud, root_segs, section_ids);
Expand Down Expand Up @@ -1051,11 +1057,14 @@ 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;
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

And shouldn't this be mid_local +=?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Basically it is: grid_origin + snap_to_grid(mid - grid_origin)
which means snapping to an offset grid.
I've tested it on simple cases and seems to work.

The '-d' treads on the distance_limit param, so I had to change it to -j in my PR. One reason why I like to minimise command line parameters.

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];
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If I am understanding your implementation correctly, you don't need to add the params_->grid_origin[0], for subsequent calls as well, otherwise, you remove the grid origin shift?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Basically it is: grid_origin + snap_to_grid(mid - grid_origin)
which means snapping to an offset grid.
I've tested it on simple cases and seems to work.

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
Expand All @@ -1073,6 +1082,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 &section : 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<int> &root_segs, const std::vector<int> &section_ids)
{
Expand Down
16 changes: 10 additions & 6 deletions raylib/extraction/raytrees.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -99,6 +101,8 @@ class RAYLIB_EXPORT Trees
double radius(const BranchSection &section) const;
/// remove elements of nodes that are too distant to the set of end points
bool removeDistantPoints(std::vector<int> &nodes);
/// remove trees with radius under the specified minimum
void removeSmallRadiusTrees();

// cached data that is used throughout the processing method
int sec_;
Expand Down
2 changes: 1 addition & 1 deletion raylib/raylaz.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am not certain this error will inform if file does not exist.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It does. Just tested it.
"readLas: failed to open file doesnottexist.blh, error: No such file or directory"

return false;
}

Expand Down