add search radius plugin

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)
 

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 <ranges>` 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 <ranges>` 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 <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]
+

environment.yml

@@ -8,6 +8,9 @@ dependencies:
   - black
   - isort
   - shapely
+  - gdal
+  - geopanda
+  - cmake
   - pip:
     - ign-pdal-tools
 

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

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, "" );
 
 }

src/filter_radius_search/CMakeLists.txt

@@ -0,0 +1,17 @@
+
+file( GLOB_RECURSE GD_SRCS 
+	${CMAKE_SOURCE_DIR}/src/filter_radius_search/*.hpp
+	${CMAKE_SOURCE_DIR}/src/filter_radius_search/*.cpp
+	${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
+)
+

src/filter_radius_search/radius_searchFilter.cpp

@@ -0,0 +1,169 @@
+#include "radius_searchFilter.hpp"
+
+#include <pdal/PipelineManager.hpp>
+#include <pdal/StageFactory.hpp>
+#include <pdal/util/ProgramArgs.hpp>
+
+#include <pdal/dimension.hpp>
+
+#include <iostream>
+#include <utility>
+
+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_srcDomainSpec);
+    args.add("reference_domain", "Selects which points will be considered as potential neighbors", m_args->m_referenceDomainSpec);
+    args.add("radius", "Distance of neighbors to consult", m_args->m_radius);
+    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::initializeDomain(StringList domainSpec, std::vector<DimRange> &domain)
+{
+    for (auto const& r : domainSpec)
+    {
+        try
+        {
+            DimRange range;
+            range.parse(r);
+            domain.push_back(range);
+        }
+        catch (const DimRange::error& err)
+        {
+            throwError("Invalid 'domain' option: '" + r + "': " + err.what());
+        }
+    }
+}
+
+void RadiusSearchFilter::addDimensions(PointLayoutPtr layout)
+{
+    m_args->m_dim = layout->registerOrAssignDim(m_args->m_nameAddAttribute, Dimension::Type::Double);
+}
+
+void RadiusSearchFilter::initialize()
+{
+    this->initializeDomain(m_args->m_referenceDomainSpec, m_args->m_referenceDomain);
+    this->initializeDomain(m_args->m_srcDomainSpec, m_args->m_srcDomain);
+
+    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::preparedDomain(std::vector<DimRange> &domain, PointLayoutPtr layout)
+{
+    for (auto& r : domain)
+    {
+        r.m_id = layout->findDim(r.m_name);
+        if (r.m_id == Dimension::Id::Unknown)
+            throwError("Invalid dimension name in 'srcDomain' option: '" + r.m_name + "'.");
+    }
+    std::sort(domain.begin(), domain.end());
+}
+
+void RadiusSearchFilter::prepared(PointTableRef table)
+{
+    PointLayoutPtr layout(table.layout());
+    this->preparedDomain(m_args->m_srcDomain, layout);
+    this->preparedDomain(m_args->m_referenceDomain, 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);
+
+    for (DimRange& r : m_args->m_srcDomain)
+    {   // process only points that satisfy a domain condition
+        if (r.valuePasses(point.getFieldAs<double>(r.m_id)))
+        {
+            doOneNoDomain(point);
+            break;
+        }
+    }
+    return true;
+}
+
+void RadiusSearchFilter::filter(PointView& view)
+{
+    PointRef point_src(view, 0);
+    // Create a kd tree only with the points in the reference domain (to make the search faster)
+    PointRef temp(view, 0);
+
+    refView = view.makeNew();
+    for (PointId id = 0; id < view.size(); ++id)
+    {
+        temp.setPointId(id);
+
+        // initialisation by default
+        temp.setField(m_args->m_dim, int64_t(0));
+
+        for (DimRange& r : m_args->m_referenceDomain)
+        {
+            // process only points that satisfy a domain condition
+            if (r.valuePasses(temp.getFieldAs<double>(r.m_id)))
+            {
+                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
+

src/filter_radius_search/radius_searchFilter.hpp

@@ -0,0 +1,58 @@
+#pragma once
+
+#include <pdal/Filter.hpp>
+#include <pdal/KDIndex.hpp>
+#include <unordered_map>
+
+#include "../utils/DimRange.hpp"
+
+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
+    {
+        StringList m_referenceDomainSpec;
+        std::vector<DimRange> m_referenceDomain;
+        StringList m_srcDomainSpec;
+        std::vector<DimRange> m_srcDomain;
+        double m_radius;
+        PointIdList m_ptsToUpdate;
+        std::string m_nameAddAttribute;
+        Dimension::Id m_dim;
+        bool search3d;
+    };
+    std::unique_ptr<RadiusSearchArgs> m_args;
+    PointViewPtr refView;
+
+    virtual void addArgs(ProgramArgs& args);
+    virtual void preparedDomain(std::vector<DimRange> &domain, PointLayoutPtr layout);
+    virtual void prepared(PointTableRef table);
+    virtual void filter(PointView& view);
+    virtual void initializeDomain(StringList domainSpec, std::vector<DimRange> &domain);
+    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