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/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 0000000..2fe0fe5 --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,16 @@ +repos: + - repo: https://github.com/ambv/black + rev: 23.12.0 + hooks: + - id: black + language_version: python3.11 + - repo: https://github.com/pycqa/flake8 + rev: 6.1.0 + hooks: + - id: flake8 + - repo: https://github.com/pycqa/isort + rev: 5.13.2 + hooks: + - id: isort + name: isort (python) + args: ["--profile", "black"] diff --git a/CMakeLists.txt b/CMakeLists.txt index 47fe925..e5785cb 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_assign) diff --git a/Dockerfile b/Dockerfile index ce43677..3a1cbfd 100755 --- a/Dockerfile +++ b/Dockerfile @@ -1,6 +1,6 @@ FROM mambaorg/micromamba:bullseye-slim as build -COPY environment.yml /environment_docker.yml +COPY environment_docker.yml /environment_docker.yml USER root RUN micromamba env create -f /environment_docker.yml @@ -9,15 +9,19 @@ 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 cmake -G"Unix Makefiles" -DCONDA_PREFIX=$CONDA_PREFIX -DCMAKE_BUILD_TYPE=Release RUN make -j4 install 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 /opt/conda/envs/pdal_ign_plugin /opt/conda/envs/pdal_ign_plugin +RUN mkdir -p /pdal_ign_plugin +COPY --from=build /tmp/install/lib /pdal_ign_plugin/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 +ENV PDAL_DRIVER_PATH=/pdal_ign_plugin/install/lib + diff --git a/README.md b/README.md index 89fda41..c98a2f5 100755 --- a/README.md +++ b/README.md @@ -45,6 +45,8 @@ python -m pytest -s [grid decimation](./doc/grid_decimation.md) +[radius assign](./doc/radius_assign.md) + ## Adding a filter In order to add a filter, you have to add a new folder in the src directory : 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_decimation.md b/doc/grid_decimation.md index 573a46b..6638898 100755 --- a/doc/grid_decimation.md +++ b/doc/grid_decimation.md @@ -1,9 +1,11 @@ # filter grid decimation +**Deprecated** : *better use the gridDecimation filter of pdal > 2.7* + Purpose --------------------------------------------------------------------------------------------------------- -The **grid decimation filter** transform only one point in each cells of a grid calculated from the points cloud and a resolution therm. The transformation is done by the value information. The selected point could be the highest or the lowest point on the cell. It can be used, for exemple, to quickly filter vegetation points in order to keep only the canopy points. A new attribut is created with the value '1' for the grid, and '0' for the other points. +The **grid decimation filter** transform only one point in each cells of a grid calculated from the points cloud and a resolution therm. The transformation is done by the value information. The selected point could be the highest or the lowest point on the cell. It can be used, for exemple, to quickly filter vegetation points in order to keep only the canopy points. A new dimension is created with the value '1' for the grid, and '0' for the other points. Example @@ -18,7 +20,7 @@ This example transform highest points of classification 5 in classification 9, o { "type": "filters.gridDecimation", "output_type":"max", - "output_name_attribut": "grid", + "output_dimension": "grid", "output_wkt":"file-output.wkt" }, { @@ -37,6 +39,6 @@ Options **resolution** : The resolution of the cells in meter. [Default: 1.] -**output_name_attribut**: The name of the new attribut. [Default: grid] +**output_dimension**: The name of the new dimension. [Default: grid] **output_wkt**: the name of the export grid file as wkt polygon. If none, no export [Default:""] diff --git a/doc/radius_assign.md b/doc/radius_assign.md new file mode 100755 index 0000000..e215895 --- /dev/null +++ b/doc/radius_assign.md @@ -0,0 +1,51 @@ +# filter radius assign + +Purpose +--------------------------------------------------------------------------------------------------------- + +The **radius assign filter** overwrites the output_dimension_ dimension with boolean values: +* 1 if the point has any neighbor with a distance lower than radius_ that belongs to the domain reference_domain_ +* 0 otherwise. + + +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.radius_assign", + "src_domain" : "Classification[1:2]", + "reference_domain" : "Classification[6:6]", + "radius" : 1, + "output_dimension": "radius", + "is3d": True + }, + "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. + +**output_dimension**: The name of the new dimension'. [Default: radius] + +**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_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.] + diff --git a/environment.yml b/environment.yml index fdc1eaf..695737c 100755 --- a/environment.yml +++ b/environment.yml @@ -1,14 +1,19 @@ name: pdal_ign_plugin channels: - conda-forge + - anaconda dependencies: - pdal - python-pdal + - gdal +# --------- dev dep --------- # + - cmake + - pre-commit # hooks for applying linters on commit + - black # code formatting + - isort # import sorting + - flake8 # code analysis - pytest - - black - - isort - - shapely +# --------- pip & pip librairies --------- # + - pip - pip: - ign-pdal-tools - - diff --git a/environment_docker.yml b/environment_docker.yml index b49f938..7e23ef6 100755 --- a/environment_docker.yml +++ b/environment_docker.yml @@ -1,6 +1,13 @@ name: pdal_ign_plugin channels: - conda-forge + - anaconda dependencies: - pdal + - python-pdal + - gdal +# --------- pip & pip librairies --------- # + - pip + - pip: + - ign-pdal-tools diff --git a/macro/__init__.py b/macro/__init__.py new file mode 100755 index 0000000..e69de29 diff --git a/macro/ex_filtering_points.py b/macro/ex_filtering_points.py new file mode 100755 index 0000000..f79af6d --- /dev/null +++ b/macro/ex_filtering_points.py @@ -0,0 +1,94 @@ +import argparse +import pdal +import macro + +""" +This tool shows how to use functions of macro in a pdal pipeline +""" + +def parse_args(): + parser = argparse.ArgumentParser("Tool to apply pdal pipelines for DSM and DTM calculation") + parser.add_argument("--input", "-i", type=str, required=True, help="Input las file") + parser.add_argument("--output_las", "-o", type=str, required=True, help="Output cloud las file") + parser.add_argument("--output_dsm", "-s", type=str, required=True, help="Output dsm tiff file") + parser.add_argument("--output_dtm", "-t", type=str, required=True, help="Output dtm tiff file") + return parser.parse_args() + + +if __name__ == "__main__": + args = parse_args() + + pipeline = pdal.Reader.las(args.input) + + ## 1 - recherche des points max de végétation (4,5) sur une grille régulière, avec prise en compte des points sol (2) et basse + ## vegetation (3) proche de la végétation : on les affecte en 100 + + # bouche trou : assigne les points sol en 102 à l'intérieur de la veget (4,5) + pipeline = macro.add_radius_assign(pipeline, 1, False, condition_src="Classification==2", condition_ref=macro.build_condition("Classification", [4,5]), condition_out="Classification=102") + pipeline = macro.add_radius_assign(pipeline, 1, False, condition_src="Classification==102", condition_ref="Classification==2", condition_out="Classification=2") + + # selection des points de veget basse proche de la veget haute : assigne 103 + pipeline = macro.add_radius_assign(pipeline, 1, False, condition_src="Classification==3", condition_ref="Classification==5", condition_out="Classification=103") + + # max des points de veget (et surement veget - 102,103) sur une grille régulière : assigne 100 + pipeline |= pdal.Filter.gridDecimation(resolution=0.75, value="Classification=100", output_type="max", where=macro.build_condition("Classification", [4,5,102,103])) + + # remise à zero des codes 102 et 103 + pipeline |= pdal.Filter.assign(value="Classification=2", where="Classification==102") + pipeline |= pdal.Filter.assign(value="Classification=3", where="Classification==103") + + ## 2 - sélection des points pour DTM et DSM + + # selection de points sol (max) sur une grille régulière + pipeline |= pdal.Filter.gridDecimation(resolution=0.5, value="Classification=102", output_type="max", where="Classification==2") + + # selection de points DSM (max) sur une grille régulière + pipeline |= pdal.Filter.gridDecimation(resolution=0.5, value="Classification=200", output_type="max", where=macro.build_condition("Classification", [2,3,4,5,6,9,17,64,100])) + + # assigne des points sol sélectionnés (102) en 100 : les points proches de la végaétation, des ponts, de l'eau et 64 + pipeline = macro.add_radius_assign(pipeline, 1.5, False, condition_src="Classification==102", + condition_ref=macro.build_condition("Classification", [4,5,6,9,17,64,100]), condition_out="Classification=100") + + # remise à zero du code 102 + pipeline |= pdal.Filter.assign(value="Classification=2", where="Classification==102") + + ## 3 - gestion des ponts + + + + # bouche trou : on élimine les points sol (2) au milieu du pont en les mettant à 102 + pipeline = macro.add_radius_assign(pipeline, 1.5, False, condition_src="Classification==2", condition_ref="Classification==17", condition_out="Classification=102") + pipeline = macro.add_radius_assign(pipeline, 1.5, False, condition_src="Classification==102", + condition_ref=macro.build_condition("Classification", [2,3,4,5]), condition_out="Classification=2") + + # bouche trou : on élimine les points basse végétation (3) au milieu du pont en les mettant à 103 + pipeline = macro.add_radius_assign(pipeline, 1.5, False, condition_src="Classification==3", condition_ref="Classification==17", condition_out="Classification=103") + pipeline = macro.add_radius_assign(pipeline, 1.5, False, condition_src="Classification==103", + condition_ref=macro.build_condition("Classification", [2,3,4,5]), condition_out="Classification=3") + + # bouche trou : on élimine les points moyenne végétation (4) au milieu du pont en les mettant à 104 + pipeline = macro.add_radius_assign(pipeline, 1.5, False, condition_src="Classification==4", condition_ref="Classification==17", condition_out="Classification=104") + pipeline = macro.add_radius_assign(pipeline, 1.5, False, condition_src="Classification==104", + condition_ref=macro.build_condition("Classification", [2,3,4,5]), condition_out="Classification=4") + + # bouche trou : on élimine les points haute végétation (5) au milieu du pont en les mettant à 105 + pipeline = macro.add_radius_assign(pipeline, 1.5, False, condition_src="Classification==5", condition_ref="Classification==17", condition_out="Classification=105") + pipeline = macro.add_radius_assign(pipeline, 1.5, False, condition_src="Classification==105", + condition_ref=macro.build_condition("Classification", [2,3,4,5]), condition_out="Classification=5") + + # bouche trou : on élimine les points eau (9) au milieu du pont en les mettant à 109 + pipeline = macro.add_radius_assign(pipeline, 1.5, False, condition_src="Classification==9", condition_ref="Classification==17", condition_out="Classification=109") + pipeline = macro.add_radius_assign(pipeline, 1.5, False, condition_src="Classification==109", + condition_ref="Classification==9", condition_out="Classification=9") + + # step 15 et supression des points ?? + + # 4 - export du nuage + pipeline |= pdal.Writer.las(extra_dims="all",forward="all",filename=args.output_las) + + # export des DSM/DTM + pipeline |= pdal.Writer.gdal(gdaldriver="GTiff", output_type="max", resolution=2.0, filename=args.output_dtm, where=macro.build_condition("Classification", [2,66])) + pipeline |= pdal.Writer.gdal(gdaldriver="GTiff", output_type="max", resolution=2.0, filename=args.output_dsm, where=macro.build_condition("Classification", [2,3,4,5,17,64])) + + pipeline.execute() + diff --git a/macro/ex_filtering_points_with_add_dimensions.py b/macro/ex_filtering_points_with_add_dimensions.py new file mode 100755 index 0000000..107a263 --- /dev/null +++ b/macro/ex_filtering_points_with_add_dimensions.py @@ -0,0 +1,78 @@ +import argparse +import pdal +import macro + +""" +This tool shows how to use functions of macro in a pdal pipeline +""" + +def parse_args(): + parser = argparse.ArgumentParser("Tool to apply pdal pipelines for DSM and DTM calculation (with add dimensions for the concerned points)") + parser.add_argument("--input", "-i", type=str, required=True, help="Input las file") + parser.add_argument("--output_las", "-o", type=str, required=True, help="Output cloud las file") + parser.add_argument("--output_dsm", "-s", type=str, required=True, help="Output dsm tiff file") + parser.add_argument("--output_dtm", "-t", type=str, required=True, help="Output dtm tiff file") + return parser.parse_args() + + +if __name__ == "__main__": + args = parse_args() + + pipeline = pdal.Reader.las(args.input) + + # 0 - ajout de dimensions temporaires + pipeline |= pdal.Filter.ferry(dimensions=f"=>PT_GRID_DSM, =>PT_VEG_DSM, =>PT_GRID_DTM, =>PT_ON_BRIDGE") + + + ## 1 - recherche des points max de végétation (4,5) sur une grille régulière, avec prise en compte des points sol (2) et basse + ## vegetation (3) proche de la végétation + ## pour le calcul du DSM + + pipeline |= pdal.Filter.assign(value=["PT_VEG_DSM = 1 WHERE " + macro.build_condition("Classification", [4,5])]) + + # bouche trou : assigne les points sol à l'intérieur de la veget (4,5) + pipeline = macro.add_radius_assign(pipeline, 1, False, condition_src="Classification==2", condition_ref=macro.build_condition("Classification", [4,5]), condition_out="PT_VEG_DSM=1") + pipeline = macro.add_radius_assign(pipeline, 1, False, condition_src="PT_VEG_DSM==1 && Classification==2", condition_ref="Classification==2", condition_out="PT_VEG_DSM=0") + + # selection des points de veget basse proche de la veget haute + pipeline = macro.add_radius_assign(pipeline, 1, False, condition_src="Classification==3", condition_ref="Classification==5", condition_out="PT_VEG_DSM=1") + + # max des points de veget (PT_VEG_DSM==1) sur une grille régulière : + pipeline |= pdal.Filter.gridDecimation(resolution=0.75, value="PT_GRID_DSM=1", output_type="max", where="PT_VEG_DSM==1") + + + ## 2 - sélection des points pour DTM et DSM + + # selection de points DTM (max) sur une grille régulière + pipeline |= pdal.Filter.gridDecimation(resolution=0.5, value="PT_GRID_DTM=1", output_type="max", where="Classification==2") + + # selection de points DSM (max) sur une grille régulière + pipeline |= pdal.Filter.gridDecimation(resolution=0.5, value="PT_GRID_DSM=1", output_type="max", + where="(" + macro.build_condition("Classification", [6,9,17,64]) + ") || PT_GRID_DSM==1") + + # assigne des points sol sélectionnés : les points proches de la végétation, des ponts, de l'eau, 64 + pipeline = macro.add_radius_assign(pipeline, 1.5, False, condition_src="PT_GRID_DTM==1", + condition_ref=macro.build_condition("Classification", [4,5,6,9,17,64]), + condition_out="PT_GRID_DSM=1") + + + ## 3 - gestion des ponts + # bouche trou : on filtre les points (2,3,4,5,9) au milieu du pont en les mettant à PT_ON_BRIDGE=1 + + pipeline = macro.add_radius_assign(pipeline, 1.5, False, condition_src=macro.build_condition("Classification", [2,3,4,5,9]), condition_ref="Classification==17", condition_out="PT_ON_BRIDGE=1") + pipeline = macro.add_radius_assign(pipeline, 1.5, False, condition_src="PT_ON_BRIDGE==1", + condition_ref=macro.build_condition("Classification", [2,3,4,5]), condition_out="PT_ON_BRIDGE=0") + pipeline |= pdal.Filter.assign(value=["PT_GRID_DSM=0 WHERE PT_ON_BRIDGE==1"]) + + + ## 4 - point pour DTM servent au DSM également + pipeline |= pdal.Filter.assign(value=["PT_GRID_DSM=1 WHERE PT_GRID_DTM==1"]) + + ## 5 - export du nuage et des DSM + + pipeline |= pdal.Writer.las(extra_dims="all", forward="all", filename=args.output_las) + pipeline |= pdal.Writer.gdal(gdaldriver="GTiff", output_type="max", resolution=2.0, filename=args.output_dtm, where="PT_GRID_DTM==1") + pipeline |= pdal.Writer.gdal(gdaldriver="GTiff", output_type="max", resolution=2.0, filename=args.output_dsm, where="PT_GRID_DSM==1") + + pipeline.execute() + diff --git a/macro/macro.py b/macro/macro.py new file mode 100755 index 0000000..71d3684 --- /dev/null +++ b/macro/macro.py @@ -0,0 +1,63 @@ +import argparse +import pdal + +""" +Some useful filters combinations for complete pdal pipeline +""" + + +def add_radius_assign(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, =>SRC_DOMAIN, =>radius_search") + pipeline |= pdal.Filter.assign(value=["SRS_DOMAIN = 0", f"SRC_DOMAIN = 1 WHERE {condition_src}", + "REF_DOMAIN = 0", f"REF_DOMAIN = 1 WHERE {condition_ref}", + "radius_search = 0"]) + pipeline |= pdal.Filter.radius_assign(radius=radius, src_domain="SRC_DOMAIN",reference_domain="REF_DOMAIN", + output_dimension="radius_search", is3d=search_3d) + pipeline |= pdal.Filter.assign(value=condition_out,where="radius_search==1") + 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 + 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/GridDecimationFilter.cpp similarity index 85% rename from src/filter_grid_decimation/grid_decimationFilter.cpp rename to src/filter_grid_decimation/GridDecimationFilter.cpp index 60f03a9..1366dcb 100755 --- a/src/filter_grid_decimation/grid_decimationFilter.cpp +++ b/src/filter_grid_decimation/GridDecimationFilter.cpp @@ -5,7 +5,7 @@ * ****************************************************************************/ -#include "grid_decimationFilter.hpp" +#include "GridDecimationFilter.hpp" #include #include @@ -18,7 +18,7 @@ namespace pdal static StaticPluginInfo const s_info { - "filters.grid_decimation", + "filters.grid_decimation_deprecated", // better to use the pdal gridDecimation plugIN "keep max or min points in a grid", "", }; @@ -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_dimension", "Name of the added dimension", m_args->m_nameOutDimension, "grid" ); args.add("output_wkt", "Export the grid as wkt", m_args->m_nameWktgrid, "" ); } @@ -61,8 +61,8 @@ void GridDecimationFilter::ready(PointTableRef table) if (m_args->m_methodKeep != "max" && m_args->m_methodKeep != "min") throwError("The output_type must be 'max' or 'min'."); - if (m_args->m_nameAddAttribut.empty()) - throwError("The output_name_attribut must be given."); + if (m_args->m_nameOutDimension.empty()) + throwError("The output_dimension must be given."); if (!m_args->m_nameWktgrid.empty()) std::remove(m_args->m_nameWktgrid.c_str()); @@ -70,8 +70,7 @@ void GridDecimationFilter::ready(PointTableRef table) void GridDecimationFilter::addDimensions(PointLayoutPtr layout) { - m_args->m_dim = layout->registerOrAssignDim(m_args->m_nameAddAttribut, - Dimension::Type::Double); + m_args->m_dim = layout->registerOrAssignDim(m_args->m_nameOutDimension, Dimension::Type::Unsigned8); } void GridDecimationFilter::processOne(BOX2D bounds, PointRef& point, PointViewPtr view) @@ -81,8 +80,8 @@ void GridDecimationFilter::processOne(BOX2D bounds, PointRef& point, PointViewPt double y = point.getFieldAs(Dimension::Id::Y); int id = point.getFieldAs(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::ceil((x - bounds.minx) / m_args->m_edgeLength); + double d_height_pt = std::ceil((y - bounds.miny) / m_args->m_edgeLength); int width = static_cast(d_width_pt); int height = static_cast(d_height_pt); @@ -108,8 +107,8 @@ void GridDecimationFilter::processOne(BOX2D bounds, PointRef& point, PointViewPt void GridDecimationFilter::createGrid(BOX2D bounds) { - double d_width = std::floor((bounds.maxx - bounds.minx) / m_args->m_edgeLength) + 1; - double d_height = std::floor((bounds.maxy - bounds.miny) / m_args->m_edgeLength) + 1; + double d_width = std::ceil((bounds.maxx - bounds.minx) / m_args->m_edgeLength); + double d_height = std::ceil((bounds.maxy - bounds.miny) / m_args->m_edgeLength); if (d_width < 0.0 || d_width > (std::numeric_limits::max)()) throwError("Grid width out of range."); diff --git a/src/filter_grid_decimation/grid_decimationFilter.hpp b/src/filter_grid_decimation/GridDecimationFilter.hpp similarity index 95% rename from src/filter_grid_decimation/grid_decimationFilter.hpp rename to src/filter_grid_decimation/GridDecimationFilter.hpp index 781d8eb..231a6a1 100755 --- a/src/filter_grid_decimation/grid_decimationFilter.hpp +++ b/src/filter_grid_decimation/GridDecimationFilter.hpp @@ -32,7 +32,7 @@ class PDAL_DLL GridDecimationFilter : public Filter { std::string m_methodKeep; // type of output (min, max) double m_edgeLength; // lenght of grid - std::string m_nameAddAttribut; // name of the new attribut + std::string m_nameOutDimension; // name of the new dimension std::string m_nameWktgrid; // export wkt grid Dimension::Id m_dim; }; diff --git a/src/filter_radius_assign/CMakeLists.txt b/src/filter_radius_assign/CMakeLists.txt new file mode 100755 index 0000000..2a49e7c --- /dev/null +++ b/src/filter_radius_assign/CMakeLists.txt @@ -0,0 +1,16 @@ + +file( GLOB_RECURSE GD_SRCS + ${CMAKE_SOURCE_DIR}/src/filter_radius_assign/*.hpp + ${CMAKE_SOURCE_DIR}/src/filter_radius_assign/*.cpp) + +PDAL_CREATE_PLUGIN( + TYPE filter + NAME radius_assign + VERSION 1.0 + SOURCES ${GD_SRCS} +) + +install(TARGETS + pdal_plugin_filter_radius_assign +) + diff --git a/src/filter_radius_assign/RadiusAssignFilter.cpp b/src/filter_radius_assign/RadiusAssignFilter.cpp new file mode 100644 index 0000000..ca54ee0 --- /dev/null +++ b/src/filter_radius_assign/RadiusAssignFilter.cpp @@ -0,0 +1,140 @@ +#include "RadiusAssignFilter.hpp" + +#include +#include +#include + +#include + +#include +#include + +namespace pdal +{ + +static PluginInfo const s_info = PluginInfo( + "filters.radius_assign", + "Re-assign some point dimension based on KNN voting", + "" ); + +CREATE_SHARED_STAGE(RadiusAssignFilter, s_info) + +std::string RadiusAssignFilter::getName() const { return s_info.name; } + +RadiusAssignFilter::RadiusAssignFilter() : +m_args(new RadiusAssignFilter::RadiusAssignArgs) +{} + + +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). 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::prepared(PointTableRef table) +{ + PointLayoutPtr layout(table.layout()); +} + +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) + { + double Zref = point.getFieldAs(Dimension::Id::Z); + if (m_args->m_max2d_below>0 || m_args->m_max2d_above>0) + { + 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; + } + } + + 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); + } + + 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_assign/RadiusAssignFilter.hpp b/src/filter_radius_assign/RadiusAssignFilter.hpp new file mode 100644 index 0000000..d94f6a8 --- /dev/null +++ b/src/filter_radius_assign/RadiusAssignFilter.hpp @@ -0,0 +1,54 @@ +#pragma once + +#include +#include +#include + +extern "C" int32_t RadiusAssignFilter_ExitFunc(); +extern "C" PF_ExitFunc RadiusAssignFilter_InitPlugin(); + +namespace pdal +{ + +class PDAL_DLL RadiusAssignFilter : public Filter +{ +public: + RadiusAssignFilter(); + ~RadiusAssignFilter(); + + static void * create(); + static int32_t destroy(void *); + std::string getName() const; + +private: + + struct RadiusAssignArgs + { + std::string m_referenceDomain; + std::string m_srcDomain; + double m_radius; + PointIdList m_ptsToUpdate; + std::string m_outputDimension; + Dimension::Id m_dim; + bool search3d; + Dimension::Id m_dim_ref, m_dim_src; + double m_max2d_above, m_max2d_below; + }; + 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); + + RadiusAssignFilter& operator=(const RadiusAssignFilter&) = delete; + RadiusAssignFilter(const RadiusAssignFilter&) = delete; +}; + +} // namespace pdal diff --git a/test/test_grid_decimation.py b/test/test_grid_decimation.py index 9491d5a..c5c3d36 100755 --- a/test/test_grid_decimation.py +++ b/test/test_grid_decimation.py @@ -1,24 +1,22 @@ 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 run_filter(type): -def test_grid_decimation(): ini_las = "test/data/4_6.las" resolution = 10 tmp_out_las = tempfile.NamedTemporaryFile(suffix=".las").name tmp_out_wkt = tempfile.NamedTemporaryFile(suffix=".wkt").name - filter = "filters.grid_decimation" + filter = "filters.grid_decimation_deprecated" utils.pdal_has_plugin(filter) bounds = li.las_get_xy_bounds(ini_las) @@ -26,15 +24,14 @@ 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}, { "type": filter, "resolution": resolution, - "output_type": "max", - "output_name_attribut": "grid", + "output_type": type, + "output_dimension": "grid", "output_wkt": tmp_out_wkt, }, { @@ -57,7 +54,7 @@ def test_grid_decimation(): if pt["grid"] > 0: nb_pts_grid += 1 - assert nb_pts_grid == 3196 + assert nb_pts_grid == 3231 assert nb_pts_grid <= nb_dalle data = [] @@ -67,3 +64,13 @@ def test_grid_decimation(): data.append(line[0]) assert len(data) == nb_dalle + + return nb_pts_grid + +def test_grid_decimation_max(): + run_filter("max") + +def test_grid_decimation_max(): + run_filter("min") + + diff --git a/test/test_radius_assign.py b/test/test_radius_assign.py new file mode 100755 index 0000000..9e884de --- /dev/null +++ b/test/test_radius_assign.py @@ -0,0 +1,145 @@ +import json +import tempfile +from test import utils + +from random import * +import numpy as np +import pdal +import random as rand +from math import * + +import pytest + + +def distance2d(pt1, pt2): + return round(sqrt((pt1[0] - pt2[0]) ** 2 + (pt1[1] - pt2[1]) ** 2), 2) + + +def distance3d(pt1, pt2): + return round(sqrt((pt1[0] - pt2[0]) ** 2 + (pt1[1] - pt2[1]) ** 2 + (pt1[2] - pt2[2]) ** 2), 2) + +def run_filter(arrays_las, distance_radius, search_3d, distance_cylinder=0. ): + + filter = "filters.radius_assign" + utils.pdal_has_plugin(filter) + + with tempfile.NamedTemporaryFile(suffix="_las_tmp.las") as las: + pipeline = pdal.Writer.las(filename=las.name).pipeline(arrays_las) + pipeline.execute() + + PIPELINE = [ + {"type": "readers.las", "filename": las.name}, + { + "type": "filters.ferry", + "dimensions": "=>SRC_DOMAIN" + }, + { + "type": "filters.ferry", + "dimensions": "=>REF_DOMAIN" + }, + { + "type": "filters.assign", + "value": [ + "SRC_DOMAIN = 1 WHERE Classification==2", + "SRC_DOMAIN = 0 WHERE Classification!=2", + "REF_DOMAIN = 1 WHERE Classification==1", + "REF_DOMAIN = 0 WHERE Classification!=1", + ], + }, + { + "type": filter, + "radius": distance_radius, + "src_domain": "SRC_DOMAIN", + "reference_domain": "REF_DOMAIN", + "output_dimension": "radius_search", + "is3d": search_3d, + "max2d_above": distance_cylinder, + "max2d_below": distance_cylinder, + } + ] + + pipeline = pdal.Pipeline(json.dumps(PIPELINE)) + pipeline.execute() + arrays = pipeline.arrays + array = arrays[0] + + nb_pts_radius_search = 0 + for pt in array: + if pt["radius_search"] > 0: + nb_pts_radius_search += 1 + + return nb_pts_radius_search + + +def build_random_points_around_one_point(test_function): + + pt_x = 1639825.15 + pt_y = 1454924.63 + pt_z = 7072.17 + pt_ini = (pt_x, pt_y, pt_z, 1) + + dtype = [('X', '