diff --git a/.gitignore b/.gitignore index 0518d62..d8c4899 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,5 @@ xcode install __pycache__ -test/__pycache__ +test/__pycache_ +test/.idea diff --git a/CMakeLists.txt b/CMakeLists.txt index 47fe925..bd23b83 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -3,6 +3,7 @@ cmake_minimum_required( VERSION 3.5 ) project(MY_READER LANGUAGES CXX) set(CMAKE_PREFIX_PATH ${CONDA_PREFIX}) +set(CMAKE_XCODE_ATTRIBUTE_OTHER_CODE_SIGN_FLAGS "-o linker-signed") find_package(PDAL REQUIRED) @@ -11,4 +12,5 @@ set(CMAKE_DEBUG_POSTFIX d) ## add plugin add_subdirectory(src/filter_grid_decimation) +add_subdirectory(src/filter_radius_search) diff --git a/ci/build.sh b/ci/build.sh index d745403..b226faf 100755 --- a/ci/build.sh +++ b/ci/build.sh @@ -24,6 +24,7 @@ fi conda activate pdal_ign_plugin export CONDA_PREFIX=$CONDA_PREFIX +echo conda is $CONDA_PREFIX mkdir build cd build diff --git a/doc/grid_radius_search.md b/doc/grid_radius_search.md new file mode 100755 index 0000000..09b385b --- /dev/null +++ b/doc/grid_radius_search.md @@ -0,0 +1,45 @@ +# filter grid decimation + +Purpose +--------------------------------------------------------------------------------------------------------- + +The **radius search filter** add a new attribut where the value depends on their neighbors in a given radius: For each point in the domain src_domain_, if it has any neighbor with a distance lower than radius_ that belongs to the domain reference_domain_, it is updated. + + +Example +--------------------------------------------------------------------------------------------------------- + +This pipeline updates the Keypoint dimension of all points with classification 1 to 2 (unclassified and ground) that are closer than 1 meter from a point with classification 6 (building) + + +``` + [ + "file-input.las", + { + "type" : "filters.radiussearch", + "src_domain" : "Classification[1:2]", + "reference_domain" : "Classification[6:6]", + "radius" : 1, + "output_name_attribut": "radius" + }, + "output.las" + ] +``` + +Options +--------------------------------------------------------------------------------------------------------------------------------------------------------------------- + +**src_domain** : + A :ref:`range ` which selects points to be processed by the filter. Can be specified multiple times. Points satisfying any range will be processed + +**reference_domain** : + A :ref:`range ` which selects points that can are considered as potential neighbors. Can be specified multiple times. + +**radius** : + An positive float which specifies the radius for the neighbors search. + +**update_expression** : + A list of :ref:`assignment expressions ` to be applied to the points that satisfy the radius search. The list of values is evaluated in order. + +**output_name_attribut**: The name of the new attribut. [Default: radius] + diff --git a/environment.yml b/environment.yml index fdc1eaf..781039d 100755 --- a/environment.yml +++ b/environment.yml @@ -1,6 +1,7 @@ name: pdal_ign_plugin channels: - conda-forge + - anaconda dependencies: - pdal - python-pdal @@ -8,6 +9,8 @@ dependencies: - black - isort - shapely + - gdal + - cmake - pip: - ign-pdal-tools diff --git a/src/filter_grid_decimation/CMakeLists.txt b/src/filter_grid_decimation/CMakeLists.txt index 2e358ca..40ffda7 100755 --- a/src/filter_grid_decimation/CMakeLists.txt +++ b/src/filter_grid_decimation/CMakeLists.txt @@ -1,7 +1,7 @@ -file( GLOB_RECURSE GD_SRCS ${CMAKE_SOURCE_DIR} *) - -message("GD_SRCS ${GD_SRCS}") +file( GLOB_RECURSE GD_SRCS + ${CMAKE_SOURCE_DIR}/src/filter_grid_decimation/*.hpp + ${CMAKE_SOURCE_DIR}/src/filter_grid_decimation/*.cpp) PDAL_CREATE_PLUGIN( TYPE filter diff --git a/src/filter_grid_decimation/grid_decimationFilter.cpp b/src/filter_grid_decimation/grid_decimationFilter.cpp index 60f03a9..69a9e28 100755 --- a/src/filter_grid_decimation/grid_decimationFilter.cpp +++ b/src/filter_grid_decimation/grid_decimationFilter.cpp @@ -39,7 +39,7 @@ void GridDecimationFilter::addArgs(ProgramArgs& args) { args.add("resolution", "Cell edge size, in units of X/Y",m_args->m_edgeLength, 1.); args.add("output_type", "Point keept into the cells ('min', 'max')", m_args->m_methodKeep, "max" ); - args.add("output_name_attribut", "Name of the add attribut", m_args->m_nameAddAttribut, "grid" ); + args.add("output_name_attribut", "Name of the added attribut", m_args->m_nameAddAttribut, "grid" ); args.add("output_wkt", "Export the grid as wkt", m_args->m_nameWktgrid, "" ); } diff --git a/src/filter_radius_search/CMakeLists.txt b/src/filter_radius_search/CMakeLists.txt new file mode 100755 index 0000000..7f63b22 --- /dev/null +++ b/src/filter_radius_search/CMakeLists.txt @@ -0,0 +1,15 @@ + +file( GLOB_RECURSE GD_SRCS + ${CMAKE_SOURCE_DIR}/src/utils/DimRange.*) + +PDAL_CREATE_PLUGIN( + TYPE filter + NAME radius_search + VERSION 1.0 + SOURCES ${GD_SRCS} +) + +install(TARGETS + pdal_plugin_filter_radius_search +) + diff --git a/src/filter_radius_search/radius_searchFilter.cpp b/src/filter_radius_search/radius_searchFilter.cpp new file mode 100644 index 0000000..82fb455 --- /dev/null +++ b/src/filter_radius_search/radius_searchFilter.cpp @@ -0,0 +1,126 @@ +#include "radius_searchFilter.hpp" + +#include +#include +#include + +#include + +#include +#include + +namespace pdal +{ + +static PluginInfo const s_info = PluginInfo( + "filters.radius_search", + "Re-assign some point attributes based KNN voting", + "" ); + +CREATE_SHARED_STAGE(RadiusSearchFilter, s_info) + +std::string RadiusSearchFilter::getName() const { return s_info.name; } + +RadiusSearchFilter::RadiusSearchFilter() : +m_args(new RadiusSearchFilter::RadiusSearchArgs) +{} + + +RadiusSearchFilter::~RadiusSearchFilter() +{} + + +void RadiusSearchFilter::addArgs(ProgramArgs& args) +{ + args.add("src_domain", "Selects which points will be subject to radius-based neighbors search", m_args->m_srcDomain, "SRS_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_name_attribute", "Name of the added attribut", m_args->m_nameAddAttribute, "radius" ); + args.add("3d_search", "Search in 3d", m_args->search3d, false ); +} + +void RadiusSearchFilter::addDimensions(PointLayoutPtr layout) +{ + m_args->m_dim = layout->registerOrAssignDim(m_args->m_nameAddAttribute, Dimension::Type::Double); + 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 RadiusSearchFilter::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_nameAddAttribute.empty()) + throwError("The output_name_attribut must be given."); +} + +void RadiusSearchFilter::prepared(PointTableRef table) +{ + PointLayoutPtr layout(table.layout()); +} + +void RadiusSearchFilter::ready(PointTableRef) +{ + m_args->m_ptsToUpdate.clear(); +} + +void RadiusSearchFilter::doOneNoDomain(PointRef &point) +{ + // build3dIndex and build2dIndex are buuild 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; + + m_args->m_ptsToUpdate.push_back(point.pointId()); +} + + +// update point. kdi and temp both reference the NN point cloud +bool RadiusSearchFilter::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 RadiusSearchFilter::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)); + + // process only points that satisfy a domain condition + //if (r.valuePasses(temp.getFieldAs(r.m_id))) + if (temp.getFieldAs(m_args->m_dim_ref)>0) + { + refView->appendPoint(view, id); + break; + } + } + + 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/src/filter_radius_search/radius_searchFilter.hpp b/src/filter_radius_search/radius_searchFilter.hpp new file mode 100644 index 0000000..515c155 --- /dev/null +++ b/src/filter_radius_search/radius_searchFilter.hpp @@ -0,0 +1,53 @@ +#pragma once + +#include +#include +#include + +extern "C" int32_t RadiusSearchFilter_ExitFunc(); +extern "C" PF_ExitFunc RadiusSearchFilter_InitPlugin(); + +namespace pdal +{ + +class PDAL_DLL RadiusSearchFilter : public Filter +{ +public: + RadiusSearchFilter(); + ~RadiusSearchFilter(); + + static void * create(); + static int32_t destroy(void *); + std::string getName() const; + +private: + + struct RadiusSearchArgs + { + std::string m_referenceDomain; + std::string m_srcDomain; + double m_radius; + PointIdList m_ptsToUpdate; + std::string m_nameAddAttribute; + Dimension::Id m_dim; + bool search3d; + Dimension::Id m_dim_ref, m_dim_src; + }; + std::unique_ptr m_args; + PointViewPtr refView; + + virtual void addArgs(ProgramArgs& args); + virtual void prepared(PointTableRef table); + virtual void filter(PointView& view); + virtual void initialize(); + virtual void addDimensions(PointLayoutPtr layout); + virtual void ready(PointTableRef); + + bool doOne(PointRef& point); + void doOneNoDomain(PointRef &point); + + RadiusSearchFilter& operator=(const RadiusSearchFilter&) = delete; + RadiusSearchFilter(const RadiusSearchFilter&) = delete; +}; + +} // namespace pdal diff --git a/test/test_grid_decimation.py b/test/test_grid_decimation.py index 9491d5a..ab7f8d9 100755 --- a/test/test_grid_decimation.py +++ b/test/test_grid_decimation.py @@ -1,17 +1,15 @@ import csv import json import math -import os import tempfile from test import utils import pdal import pdaltools.las_info as li import pytest -import shapely - def test_grid_decimation(): + ini_las = "test/data/4_6.las" resolution = 10 @@ -26,7 +24,6 @@ def test_grid_decimation(): d_width = math.floor((bounds[0][1] - bounds[0][0]) / resolution) + 1 d_height = math.floor((bounds[1][1] - bounds[1][0]) / resolution) + 1 nb_dalle = d_width * d_height - print("size of the grid", nb_dalle) PIPELINE = [ {"type": "readers.las", "filename": ini_las}, diff --git a/test/test_radius_search.py b/test/test_radius_search.py new file mode 100755 index 0000000..2e077e2 --- /dev/null +++ b/test/test_radius_search.py @@ -0,0 +1,117 @@ +import json +import tempfile +from test import utils + +from shapely.geometry import Point + +from random import * +import numpy as np +import pdal +import random as rand +from math import * + +import pytest + +def test_radius_search(): + + distance_radius = 1 + + #shapely is 2d library + def distance2d(pt1, pt2): + return pt1.distance(pt2) + def distance3d(pt1, pt2): + return sqrt( (pt1.x-pt2.x)**2 + (pt1.y-pt2.y)**2 + (pt1.z-pt2.z)**2 ) + + pt_x = 1639825.15 + pt_y = 1454924.63 + pt_z = 7072.17 + pt_ini = (pt_x, pt_y, pt_z, 1) + pt1 = Point(pt_x, pt_y, pt_z) + + dtype = [('X', 'SRS_DOMAIN" + }, + { + "type": "filters.ferry", + "dimensions": "=>REF_DOMAIN" + }, + { + "type": "filters.assign", + "value": [ + "SRS_DOMAIN = 1 WHERE Classification==2", + "SRS_DOMAIN = 0 WHERE Classification!=2", + "REF_DOMAIN = 1 WHERE Classification==1", + "REF_DOMAIN = 0 WHERE Classification!=1", + ], + }, + { + "type": filter, + "radius": "1.", + "src_domain": "SRS_DOMAIN", + "reference_domain": "REF_DOMAIN", + "output_name_attribute": "radius_2D", + "3d_search":False + }, + { + "type": filter, + "radius": "1.", + "src_domain": "SRS_DOMAIN", + "reference_domain": "REF_DOMAIN", + "output_name_attribute": "radius_3D", + "3d_search": True + } + ] + + pipeline = pdal.Pipeline(json.dumps(PIPELINE)) + + # execute the pipeline + pipeline.execute() + arrays = pipeline.arrays + array = arrays[0] + + nb_pts_radius_2d = 0 + nb_pts_radius_3d = 0 + for pt in array: + if pt["radius_2D"] > 0: + nb_pts_radius_2d += 1 + if pt["radius_3D"] > 0: + nb_pts_radius_3d += 1 + + assert nb_pts_radius_2d == nb_points_take_2d + assert nb_pts_radius_3d == nb_points_take_3d \ No newline at end of file