diff --git a/doc/radius_assign.md b/doc/radius_assign.md index e215895..b35524d 100755 --- a/doc/radius_assign.md +++ b/doc/radius_assign.md @@ -45,7 +45,7 @@ Options **is3d**: Search in 3d (as a ball). [Default: false] -**max2d_above**: If search in 2d : upward maximum distance in Z for potential neighbors (corresponds to a search in a cylinder with a height = max2d_above above the source point). Default (0) = infinite height [Default: 0.] +**max2d_above**: If search in 2d : upward maximum distance in Z for potential neighbors (corresponds to a search in a cylinder with a height = max2d_above above the source point). Values < 0 mean infinite height [Default: -1.] -**max2d_below**: If search in 2d : upward maximum distance in Z for potential neighbors (corresponds to a search in a cylinder with a height = max2d_below below the source point). Default (0) = infinite height [Default: 0.] +**max2d_below**: If search in 2d : upward maximum distance in Z for potential neighbors (corresponds to a search in a cylinder with a height = max2d_below below the source point). Values < 0 mean infinite height [Default: -1.] diff --git a/src/filter_radius_assign/RadiusAssignFilter.cpp b/src/filter_radius_assign/RadiusAssignFilter.cpp index ca54ee0..257051c 100644 --- a/src/filter_radius_assign/RadiusAssignFilter.cpp +++ b/src/filter_radius_assign/RadiusAssignFilter.cpp @@ -12,129 +12,137 @@ namespace pdal { -static PluginInfo const s_info = PluginInfo( - "filters.radius_assign", - "Re-assign some point dimension based on KNN voting", - "" ); + static PluginInfo const s_info = PluginInfo( + "filters.radius_assign", + "Re-assign some point dimension based on KNN voting", + ""); -CREATE_SHARED_STAGE(RadiusAssignFilter, s_info) + CREATE_SHARED_STAGE(RadiusAssignFilter, s_info) -std::string RadiusAssignFilter::getName() const { return s_info.name; } + std::string RadiusAssignFilter::getName() const { return s_info.name; } -RadiusAssignFilter::RadiusAssignFilter() : -m_args(new RadiusAssignFilter::RadiusAssignArgs) -{} + RadiusAssignFilter::RadiusAssignFilter() : m_args(new RadiusAssignFilter::RadiusAssignArgs) + { + } + RadiusAssignFilter::~RadiusAssignFilter() + { + } -RadiusAssignFilter::~RadiusAssignFilter() -{} + void RadiusAssignFilter::addArgs(ProgramArgs &args) + { + args.add("src_domain", "Selects which points will be subject to radius-based neighbors search", m_args->m_srcDomain, "SRC_DOMAIN"); + args.add("reference_domain", "Selects which points will be considered as potential neighbors", m_args->m_referenceDomain, "REF_DOMAIN"); + args.add("radius", "Distance of neighbors to consult", m_args->m_radius, 1.); + args.add("output_dimension", "Name of the added dimension", m_args->m_outputDimension, "radius"); + args.add("is3d", "Search in 3d", m_args->search3d, false); + args.add("max2d_above", "if search in 2d : upward maximum distance in Z for potential neighbors (corresponds to a search in a cylinder with a height = max2d_above above the source point). Values < 0 mean infinite height", m_args->m_max2d_above, -1.); + args.add("max2d_below", "if search in 2d : upward maximum distance in Z for potential neighbors (corresponds to a search in a cylinder with a height = max2d_below below the source point). Values < 0 mean infinite height", m_args->m_max2d_below, -1.); + } + void RadiusAssignFilter::addDimensions(PointLayoutPtr layout) + { + m_args->m_dim = layout->registerOrAssignDim(m_args->m_outputDimension, Dimension::Type::Unsigned8); + m_args->m_dim_ref = layout->registerOrAssignDim(m_args->m_referenceDomain, Dimension::Type::Unsigned8); + m_args->m_dim_src = layout->registerOrAssignDim(m_args->m_srcDomain, Dimension::Type::Unsigned8); + } -void RadiusAssignFilter::addArgs(ProgramArgs& args) -{ - args.add("src_domain", "Selects which points will be subject to radius-based neighbors search", m_args->m_srcDomain, "SRC_DOMAIN"); - args.add("reference_domain", "Selects which points will be considered as potential neighbors", m_args->m_referenceDomain, "REF_DOMAIN"); - args.add("radius", "Distance of neighbors to consult", m_args->m_radius, 1.); - args.add("output_dimension", "Name of the added dimension", m_args->m_outputDimension, "radius"); - args.add("is3d", "Search in 3d", m_args->search3d, false ); - args.add("max2d_above", "if search in 2d : upward maximum distance in Z for potential neighbors (corresponds to a search in a cylinder with a height = max2d_above above the source point). Default (0) = infinite height", m_args->m_max2d_above, 0.); - args.add("max2d_below", "if search in 2d : upward maximum distance in Z for potential neighbors (corresponds to a search in a cylinder with a height = max2d_below below the source point). Default (0) = infinite height", m_args->m_max2d_below, 0.); -} - -void RadiusAssignFilter::addDimensions(PointLayoutPtr layout) -{ - m_args->m_dim = layout->registerOrAssignDim(m_args->m_outputDimension, Dimension::Type::Unsigned8); - m_args->m_dim_ref = layout->registerOrAssignDim(m_args->m_referenceDomain,Dimension::Type::Unsigned8); - m_args->m_dim_src = layout->registerOrAssignDim(m_args->m_srcDomain,Dimension::Type::Unsigned8); -} + void RadiusAssignFilter::initialize() + { + if (m_args->m_referenceDomain.empty()) + throwError("The reference_domain must be given."); + if (m_args->m_radius <= 0) + throwError("Invalid 'radius' option: " + std::to_string(m_args->m_radius) + ", must be > 0"); + if (m_args->m_outputDimension.empty()) + throwError("The output_dimension must be given."); + } -void RadiusAssignFilter::initialize() -{ - if (m_args->m_referenceDomain.empty()) - throwError("The reference_domain must be given."); - if (m_args->m_radius <= 0) - throwError("Invalid 'radius' option: " + std::to_string(m_args->m_radius) + ", must be > 0"); - if (m_args->m_outputDimension.empty()) - throwError("The output_dimension must be given."); -} - -void RadiusAssignFilter::prepared(PointTableRef table) -{ - PointLayoutPtr layout(table.layout()); -} + void RadiusAssignFilter::prepared(PointTableRef table) + { + PointLayoutPtr layout(table.layout()); + } -void RadiusAssignFilter::ready(PointTableRef) -{ - m_args->m_ptsToUpdate.clear(); -} + void RadiusAssignFilter::ready(PointTableRef) + { + m_args->m_ptsToUpdate.clear(); + } -void RadiusAssignFilter::doOneNoDomain(PointRef &point) -{ - // build3dIndex and build2dIndex are build once - PointIdList iNeighbors; - if (m_args->search3d) iNeighbors = refView->build3dIndex().radius(point, m_args->m_radius); - else iNeighbors = refView->build2dIndex().radius(point, m_args->m_radius); - - if (iNeighbors.size() == 0) - return; - - if (!m_args->search3d) + void RadiusAssignFilter::doOneNoDomain(PointRef &point) { - double Zref = point.getFieldAs(Dimension::Id::Z); - if (m_args->m_max2d_below>0 || m_args->m_max2d_above>0) + // build3dIndex and build2dIndex are build once + PointIdList iNeighbors; + if (m_args->search3d) + iNeighbors = refView->build3dIndex().radius(point, m_args->m_radius); + else + iNeighbors = refView->build2dIndex().radius(point, m_args->m_radius); + + if (iNeighbors.size() == 0) + return; + + if (!m_args->search3d) { - bool take (false); - for (PointId ptId : iNeighbors) + double Zref = point.getFieldAs(Dimension::Id::Z); + if (m_args->m_max2d_below >= 0 || m_args->m_max2d_above >= 0) { - double Zpt = refView->point(ptId).getFieldAs(Dimension::Id::Z); - if (m_args->m_max2d_below>0 && Zpt>=Zref && (Zpt-Zref)<=m_args->m_max2d_below) {take=true; break;} - if (m_args->m_max2d_above>0 && Zpt<=Zref && (Zref-Zpt)<=m_args->m_max2d_above) {take=true; break;} + bool take(false); + for (PointId ptId : iNeighbors) + { + double Zpt = refView->point(ptId).getFieldAs(Dimension::Id::Z); + if (m_args->m_max2d_below >= 0 && Zpt >= Zref && (Zpt - Zref) <= m_args->m_max2d_below) + { + take = true; + break; + } + if (m_args->m_max2d_above >= 0 && Zpt <= Zref && (Zref - Zpt) <= m_args->m_max2d_above) + { + take = true; + break; + } + } + if (!take) + return; } - if (!take) return; } - } - - m_args->m_ptsToUpdate.push_back(point.pointId()); -} - -// update point. kdi and temp both reference the NN point cloud -bool RadiusAssignFilter::doOne(PointRef& point) -{ - if (m_args->m_srcDomain.empty()) // No domain, process all points - doOneNoDomain(point); - else if (point.getFieldAs(m_args->m_dim_src)>0) - doOneNoDomain(point); - return true; -} - -void RadiusAssignFilter::filter(PointView& view) -{ - PointRef point_src(view, 0); - PointRef temp(view, 0); - - refView = view.makeNew(); - for (PointId id = 0; id < view.size(); ++id) - { - temp.setPointId(id); - temp.setField(m_args->m_dim, int64_t(0)); // initialisation - - // process only points that satisfy a domain condition - if (temp.getFieldAs(m_args->m_dim_ref)>0) - refView->appendPoint(view, id); + m_args->m_ptsToUpdate.push_back(point.pointId()); } - - for (PointId id = 0; id < view.size(); ++id) + + // update point. kdi and temp both reference the NN point cloud + bool RadiusAssignFilter::doOne(PointRef &point) { - point_src.setPointId(id); - doOne(point_src); + if (m_args->m_srcDomain.empty()) // No domain, process all points + doOneNoDomain(point); + else if (point.getFieldAs(m_args->m_dim_src) > 0) + doOneNoDomain(point); + return true; } - for (auto id: m_args->m_ptsToUpdate) + + void RadiusAssignFilter::filter(PointView &view) { - temp.setPointId(id); - temp.setField(m_args->m_dim, int64_t(1)); + PointRef point_src(view, 0); + PointRef temp(view, 0); + + refView = view.makeNew(); + for (PointId id = 0; id < view.size(); ++id) + { + temp.setPointId(id); + temp.setField(m_args->m_dim, int64_t(0)); // initialisation + + // process only points that satisfy a domain condition + if (temp.getFieldAs(m_args->m_dim_ref) > 0) + refView->appendPoint(view, id); + } + + for (PointId id = 0; id < view.size(); ++id) + { + point_src.setPointId(id); + doOne(point_src); + } + for (auto id : m_args->m_ptsToUpdate) + { + temp.setPointId(id); + temp.setField(m_args->m_dim, int64_t(1)); + } } -} } // namespace pdal - diff --git a/test/test_radius_assign.py b/test/test_radius_assign.py index 3c38244..7f35205 100755 --- a/test/test_radius_assign.py +++ b/test/test_radius_assign.py @@ -6,6 +6,7 @@ import numpy as np import pdal +import pytest pt_x = 1639825.1 pt_y = 1454924.6 @@ -26,7 +27,7 @@ def distance3d(pt1, pt2): ) -def run_filter(arrays_las, distance_radius, search_3d, distance_cylinder=0.0): +def run_filter(arrays_las, distance_radius, search_3d, limit_z_above=-1, limit_w_below=-1): filter = "filters.radius_assign" utils.pdal_has_plugin(filter) @@ -55,8 +56,8 @@ def run_filter(arrays_las, distance_radius, search_3d, distance_cylinder=0.0): "reference_domain": "REF_DOMAIN", "output_dimension": "radius_search", "is3d": search_3d, - "max2d_above": distance_cylinder, - "max2d_below": distance_cylinder, + "max2d_above": limit_z_above, + "max2d_below": limit_w_below, }, ] @@ -144,20 +145,35 @@ def func_test(pt): assert nb_pts_radius_2d == nb_points_take_2d -def test_radius_assign_2d_cylinder(): +@pytest.mark.parametrize( + "limit_z_above, limit_z_below", + [ + (-1, -1), # no limit + (-1, 2), # limit below only + (2, -1), # limit above only + (0, -1), # take all points below only + (-1, 0), # take all points above only + ], +) +def test_radius_assign_2d_cylinder(limit_z_above, limit_z_below): distance_radius = 1 - distance_cylinder = 0.25 + limit_z_above = 0.25 + limit_z_below = 0.25 def func_test(pt): distance_i = distance2d(pt_ini, pt) if distance_i < distance_radius: - if abs(pt_ini[2] - pt[2]) <= distance_cylinder: + if (limit_z_above >= 0) and ((pt[2] - pt_ini[2]) <= limit_z_above): + return 1 + if (limit_z_below >= 0) and ((pt_ini[2] - pt[2]) <= limit_z_below): return 1 return 0 arrays_las, nb_points_take_2d = build_random_points_around_one_point( func_test, distance_radius ) - nb_pts_radius_2d_cylinder = run_filter(arrays_las, distance_radius, False, distance_cylinder) + nb_pts_radius_2d_cylinder = run_filter( + arrays_las, distance_radius, False, limit_z_above, limit_z_below + ) assert nb_pts_radius_2d_cylinder == nb_points_take_2d diff --git a/test/utils.py b/test/utils.py index c84fbe2..7b33c1e 100755 --- a/test/utils.py +++ b/test/utils.py @@ -3,7 +3,7 @@ def pdal_has_plugin(name_filter): - print("init pdal driver : ", os.environ["PDAL_DRIVER_PATH"]) + print("Initialize pdal driver : ", os.environ["PDAL_DRIVER_PATH"]) result = subprocess.run(["pdal", "--drivers"], stdout=subprocess.PIPE) if name_filter not in result.stdout.decode("utf-8"): - raise ValueError("le script " + name_filter + " n'est pas visible") + raise ValueError("Filter " + name_filter + " not found by `pdal --drivers`.")