Skip to content

Commit

Permalink
Enable radisAssignFilter with max2d_above or max2d_below = 0
Browse files Browse the repository at this point in the history
  • Loading branch information
@leavauchier
leavauchier committed Jun 6, 2024
1 parent 5eed8c2 commit 892aeca
Show file tree
Hide file tree
Showing 4 changed files with 138 additions and 114 deletions.
4 changes: 2 additions & 2 deletions doc/radius_assign.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.]

214 changes: 111 additions & 103 deletions src/filter_radius_assign/RadiusAssignFilter.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<double>(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<double>(Dimension::Id::Z);
if (m_args->m_max2d_below >= 0 || m_args->m_max2d_above >= 0)
{
double Zpt = refView->point(ptId).getFieldAs<double>(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<double>(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<int8_t>(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<int8_t>(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<int8_t>(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<int8_t>(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

30 changes: 23 additions & 7 deletions test/test_radius_assign.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@

import numpy as np
import pdal
import pytest

pt_x = 1639825.1
pt_y = 1454924.6
Expand All @@ -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)
Expand Down Expand Up @@ -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,
},
]

Expand Down Expand Up @@ -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
4 changes: 2 additions & 2 deletions test/utils.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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`.")

0 comments on commit 892aeca

Please sign in to comment.