add search radius plugin

.gitignore

@@ -1,4 +1,5 @@
 xcode
 install
 __pycache__
-test/__pycache__
+test/__pycache_
+test/.idea

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)
 

Dockerfile

@@ -9,6 +9,7 @@ RUN apt-get update && apt-get install --no-install-recommends -y cmake make buil
 
 COPY src src
 COPY CMakeLists.txt CMakeLists.txt
+COPY macro macro
 
 RUN cmake -G"Unix Makefiles" -DCONDA_PREFIX=$CONDA_PREFIX -DCMAKE_BUILD_TYPE=Release  
 RUN make -j4 install
@@ -17,7 +18,9 @@ FROM debian:bullseye-slim
 
 COPY --from=build /opt/conda/envs/pdal_ign_plugin /opt/conda/envs/pdal_ign_plugin
 COPY --from=build /tmp/install/lib /tmp/install/lib
-
+COPY --from=build /tmp/macro /macro 
+
 ENV PATH=$PATH:/opt/conda/envs/pdal_ign_plugin/bin/
 ENV PROJ_LIB=/opt/conda/envs/pdal_ign_plugin/share/proj/
 ENV PDAL_DRIVER_PATH=/tmp/install/lib
+

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

doc/grid_radius_search.md

@@ -0,0 +1,51 @@
+# 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]
+
+**3d_search**: Search in 3d (as a ball). [Default: false]
+
+**2d_search_above**: If 3d_search is false, take points only if exists points between the referent point and a distance above. [Default: 0.]
+
+**2d_search_bellow**: If 3d_search is false, take points only if exists points between the referent point and a distance bellow. [Default: 0.]
+

environment.yml

@@ -1,13 +1,16 @@
 name: pdal_ign_plugin
 channels:
   - conda-forge
+  - anaconda
 dependencies:
   - pdal
   - python-pdal
   - pytest
   - black
   - isort
   - shapely
+  - gdal
+  - cmake
   - pip:
     - ign-pdal-tools
 

macro/ex_filtering_points.py

@@ -0,0 +1,55 @@
+import argparse
+import pdal
+from macro import *
+
+"""
+This tool shows how to use functions of macro in a pdal pipeline
+"""
+
+def parse_args():
+    parser = argparse.ArgumentParser("Tools for apply pdal pipelines")
+    parser.add_argument("--input", "-i", type=str, required=True, help="Input las file")
+    parser.add_argument("--output", "-o", type=str, required=True, help="Output las file")
+    return parser.parse_args()
+
+
+if __name__ == "__main__":
+    args = parse_args()
+
+    pipeline = pdal.Reader.las(args.input)
+
+    # step 1 a 9
+    pipeline = add_radius_search(pipeline, 1, False, "Classification==2", build_condition("Classification", [4,5]), "Classification=102")
+    pipeline = add_radius_search(pipeline, 1, False, "Classification==102", "Classification==2", "Classification=2")
+    pipeline = add_radius_search(pipeline, 1, False, "Classification==3", "Classification==5", "Classification=103")
+    pipeline = add_grid_decimation(pipeline, 0.75, "max", build_condition("Classification", [4,5,102,103]), "Classification=100")
+    pipeline |= pdal.Filter.assign(value="Classification=2", where="Classification==102")
+    pipeline |= pdal.Filter.assign(value="Classification=3", where="Classification==103")
+    pipeline = add_grid_decimation(pipeline, 0.5, "max", "Classification==2", "Classification=102")
+    pipeline = add_grid_decimation(pipeline, 0.5, "max", build_condition("Classification", [2,3,4,5,6,9,17,64,100]), "Classification=200")
+    pipeline = add_radius_search(pipeline, 1.5, False, "Classification==102", build_condition("Classification", [4,5,6,9,17,64,100]), "Classification=100")
+    pipeline |= pdal.Filter.assign(value="Classification=2", where="Classification==102")
+
+    # step 10
+    pipeline = add_radius_search(pipeline, 1.5, False, "Classification==2", "Classification==17", "Classification=102")
+    pipeline = add_radius_search(pipeline, 1.5, False, "Classification==102", build_condition("Classification", [2,3,4,5]), "Classification=2")
+
+    # step 11
+    pipeline = add_radius_search(pipeline, 1.5, False, "Classification==3", "Classification==17", "Classification=103")
+    pipeline = add_radius_search(pipeline, 1.5, False, "Classification==103", build_condition("Classification", [2,3,4,5]), "Classification=3")
+
+    # step 12
+    pipeline = add_radius_search(pipeline, 1.5, False, "Classification==4", "Classification==17", "Classification=104")
+    pipeline = add_radius_search(pipeline, 1.5, False, "Classification==104", build_condition("Classification", [2,3,4,5]), "Classification=4")
+
+    # step 13
+    pipeline = add_radius_search(pipeline, 1.5, False, "Classification==5", "Classification==17", "Classification=105")
+    pipeline = add_radius_search(pipeline, 1.5, False, "Classification==105", build_condition("Classification", [2,3,4,5]), "Classification=5")
+
+    # step 14
+    pipeline = add_radius_search(pipeline, 1.5, False, "Classification==9", "Classification==17", "Classification=109")
+    pipeline = add_radius_search(pipeline, 1.5, False, "Classification==109", "Classification==9", "Classification=9")
+
+    pipeline |= pdal.Writer.las(extra_dims="all",minor_version=4,dataformat_id=6,filename=args.output)
+    pipeline.execute()
+

macro/macro.py

@@ -0,0 +1,81 @@
+import argparse
+import pdal
+
+"""
+Some useful filters combinations for complete pdal pipeline
+"""
+
+
+def add_radius_search(pipeline, radius, search_3d, condition_src, condition_ref, condition_out ):
+    """
+    search points from "condition_src" closed from "condition_ref", and reassign them to "condition_out"
+    This combination is equivalent to the CloseBy macro of TerraScan
+    radius : the search distance
+    search_3d : the distance reseach is in 3d if True
+    condition_src, condition_ref, condition_out : a pdal condition as "Classification==2"
+    """
+    pipeline |= pdal.Filter.ferry(dimensions=f"=>REF_DOMAIN, =>SRS_DOMAIN, =>radius_search")
+    pipeline |= pdal.Filter.assign(value=["SRS_DOMAIN = 0", f"SRS_DOMAIN = 1 WHERE {condition_src}",
+                                          "REF_DOMAIN = 0", f"REF_DOMAIN = 1 WHERE {condition_ref}",
+                                          "radius_search = 0"])
+    pipeline |= pdal.Filter.radius_search(radius=radius, src_domain="SRS_DOMAIN",reference_domain="REF_DOMAIN",
+                                          output_name_attribute="radius_search", search_3d=search_3d)
+    pipeline |= pdal.Filter.assign(value=condition_out,where="radius_search==1")
+    return pipeline
+
+
+
+def add_grid_decimation(pipeline, grid_resolution, output_type, condition, condition_out):
+    """
+    Select a points in a grid  from "condition"; points not selected are reassign to "condition_out"
+    This combination is equivalent to the Thin Points macro of TerraScan
+    grid_resolution : resolution of the grid
+    output_type : "max" or "min" (the highest or lower points of the grid)
+    condition, condition_out : a pdal condition as "Classification==2"
+    """
+    pipeline |= pdal.Filter.ferry(dimensions=f"=>grid,")
+    pipeline |= pdal.Filter.assign(value="grid = 0")
+    pipeline |= pdal.Filter.grid_decimation(resolution=grid_resolution, output_name_attribut="grid",
+                                            output_type=output_type, where=condition)
+    pipeline |= pdal.Filter.assign(value=condition_out,where=f"grid==0 && ({condition})")
+    return pipeline
+
+
+
+def classify_hgt_ground(pipeline, hmin, hmax, condition, condition_out):
+    """
+       reassign points from "condition" between "hmin" and "hmax" of the ground to "condition_out"
+       This combination is equivalent to the ClassifyHgtGrd macro of TerraScan
+       condition, condition_out : a pdal condition as "Classification==2"
+    """
+    pipeline |= pdal.Filter.hag_delaunay(allow_extrapolation=True)
+    condition_h = f"HeightAboveGround>{hmin} && HeightAboveGround<={hmax}"
+    condition_h += " && " + condition
+    pipeline |= pdal.Filter.assign(value=condition_out, where=condition_h)
+    return pipeline
+
+
+
+def keep_non_planar_pts(pipeline, condition, condition_out):
+    """
+       reassign points from "condition" who are planar to "condition_out"
+       This combination is equivalent to the ClassifyModelKey macro of TerraScan
+       condition, condition_out : a pdal condition as "Classification==2"
+    """
+    pipeline |= pdal.Filter.approximatecoplanar(knn=8,thresh1=25,thresh2=6,where=condition)
+    pipeline |= pdal.Filter.assign(value=condition_out,where=f"Coplanar==0 && ({condition})")
+    return pipeline
+
+
+
+
+def build_condition(key, values):
+    """
+       build 'key==values[0] || key==values[1] ...'
+    """
+    condition = ""
+    for v in values:
+        condition += key+"=="+str(v)
+        if v!=values[-1]:condition += " || "
+    return condition
+

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, "" );
 
 }
@@ -81,8 +81,8 @@ void GridDecimationFilter::processOne(BOX2D bounds, PointRef& point, PointViewPt
     double y = point.getFieldAs<double>(Dimension::Id::Y);
     int id = point.getFieldAs<double>(Dimension::Id::PointId);
 
-    double d_width_pt = std::floor((x - bounds.minx) / m_args->m_edgeLength) + 1;
-    double d_height_pt = std::floor((y - bounds.miny) / m_args->m_edgeLength) + 1;
+    double d_width_pt = std::floor((x - bounds.minx) / m_args->m_edgeLength);
+    double d_height_pt = std::floor((y - bounds.miny) / m_args->m_edgeLength);
 
     int width = static_cast<int>(d_width_pt);
     int height = static_cast<int>(d_height_pt);

src/filter_radius_search/CMakeLists.txt

@@ -0,0 +1,16 @@
+
+file( GLOB_RECURSE GD_SRCS 
+	${CMAKE_SOURCE_DIR}/src/filter_radius_search/*.hpp
+	${CMAKE_SOURCE_DIR}/src/filter_radius_search/*.cpp)
+
+PDAL_CREATE_PLUGIN(
+    TYPE filter
+    NAME radius_search
+    VERSION 1.0
+    SOURCES ${GD_SRCS}
+)
+
+install(TARGETS
+	pdal_plugin_filter_radius_search
+)
+