diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 91a22ec..d1b8b12 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -1,10 +1,10 @@ repos: - repo: https://github.com/ambv/black - rev: 24.4.2 + rev: 24.8.0 hooks: - id: black - repo: https://github.com/pycqa/flake8 - rev: 7.0.0 + rev: 7.1.1 hooks: - id: flake8 - repo: https://github.com/pycqa/isort diff --git a/RELEASE_NOTES.md b/RELEASE_NOTES.md index 1ff0699..93ba4d8 100644 --- a/RELEASE_NOTES.md +++ b/RELEASE_NOTES.md @@ -1,3 +1,11 @@ +### 0.3.0 + +Update algorithm of [mark_points_to_use_for_digital_models_with_new_dimension](pdal_ign_macro/mark_points_to_use_for_digital_models_with_new_dimension.py) +In details : +- manage water and virtuals points, +- update building consideration +- 2 levels of water + ### 0.2.1 Fix (and test) arguments parsing in [mark_points_to_use_for_digital_models_with_new_dimension](pdal_ign_macro/mark_points_to_use_for_digital_models_with_new_dimension.py) diff --git a/pdal_ign_macro/mark_points_to_use_for_digital_models_with_new_dimension.py b/pdal_ign_macro/mark_points_to_use_for_digital_models_with_new_dimension.py index 0fe1e19..e1750a9 100755 --- a/pdal_ign_macro/mark_points_to_use_for_digital_models_with_new_dimension.py +++ b/pdal_ign_macro/mark_points_to_use_for_digital_models_with_new_dimension.py @@ -92,21 +92,26 @@ def define_marking_pipeline(input_las, output_las, dsm_dimension, dtm_dimension) "PT_ON_BUILDING", "PT_ON_VEGET", "PT_ON_SOL", + "PT_ON_VIRT", ] added_dimensions = [dtm_dimension, dsm_dimension] + temporary_dimensions pipeline |= pdal.Filter.ferry(dimensions="=>" + ", =>".join(added_dimensions)) - # 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) proches de la végétation - # pour le calcul du DSM + ################################################################################################################### + # 1 - Gestion de la végétation pour le calcul du DSM: + ################################################################################################################### + # Recherche des points max de végétation (4,5) sur une grille régulière = PT_VEG_DSM=1, + # Prise en compte des points sol (2) proche de la végétation = PT_VEG_DSM=1 + # Isolement des classes (hors sol) sous la végetation comme les bâtis, l'eau, les ponts et les divers-bâtis + # Prise en compte des points basse vegetation (3) proche de la végétation = PT_VEG_DSM=1 + # Passage des premiers points "veget" PT_VEG_DSM=1 dans le taguage pour le MNS dsm_dimension=1 + # 1.1 Point veget max 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) + # 1.2 bouche trou : assigne les points sol à l'intérieur de la veget (4,5) pipeline = macro.add_radius_assign( pipeline, 1, @@ -115,16 +120,6 @@ def define_marking_pipeline(input_las, output_las, dsm_dimension, dtm_dimension) condition_ref=macro.build_condition("Classification", [4, 5]), condition_out="PT_VEG_DSM=1", ) - pipeline = macro.add_radius_assign( - pipeline, - 1, - False, - condition_src=macro.build_condition("Classification", [6, 9, 17, 67]), - condition_ref=macro.build_condition("Classification", [4, 5]), - condition_out="PT_ON_VEGET=1", - max2d_above=0, # ne pas prendre les points qui sont au dessus des points pont (condition_ref) - max2d_below=-1, # prendre tous les points qui sont en dessous des points pont (condition_ref) - ) pipeline = macro.add_radius_assign( pipeline, 1, @@ -133,38 +128,32 @@ def define_marking_pipeline(input_las, output_las, dsm_dimension, dtm_dimension) condition_ref="Classification==2 && PT_VEG_DSM==0", condition_out="PT_VEG_DSM=0", ) + # 1.3 Isolement en PT_ON_VEGET=1 des éléments sous la végétation (hors sol) pipeline = macro.add_radius_assign( pipeline, 1, False, - condition_src="PT_ON_VEGET==1 && " + macro.build_condition("Classification", [6, 67]), - condition_ref="PT_ON_VEGET==0 &&" + macro.build_condition("Classification", [6, 67]), - condition_out="PT_ON_VEGET=0", - max2d_above=0.5, # ne pas prendre les points qui sont au dessus des points pont (condition_ref) - max2d_below=0.5, # prendre tous les points qui sont en dessous des points pont (condition_ref) - ) - pipeline = macro.add_radius_assign( - pipeline, - 1, - False, - condition_src="PT_ON_VEGET==1 && Classification==17", - condition_ref="Classification==17 && PT_ON_VEGET==0", - condition_out="PT_ON_VEGET=0", - max2d_above=0.5, # ne pas prendre les points qui sont au dessus des points pont (condition_ref) - max2d_below=0.5, # prendre tous les points qui sont en dessous des points pont (condition_ref) + condition_src=macro.build_condition("Classification", [6, 9, 17, 67]), + condition_ref=macro.build_condition("Classification", [4, 5]), + condition_out="PT_ON_VEGET=1", + max2d_above=0, + max2d_below=-1, ) pipeline = macro.add_radius_assign( pipeline, 1, False, - condition_src="PT_ON_VEGET==1 && Classification==9", - condition_ref="Classification==9 && PT_ON_VEGET==0", + condition_src="PT_ON_VEGET==1 && ( " + + macro.build_condition("Classification", [6, 9, 17, 67]) + + " )", + condition_ref="PT_ON_VEGET==0 && ( " + + macro.build_condition("Classification", [6, 9, 17, 67]) + + " )", condition_out="PT_ON_VEGET=0", - max2d_above=0.5, # ne pas prendre les points qui sont au dessus des points pont (condition_ref) - max2d_below=0.5, # prendre tous les points qui sont en dessous des points pont (condition_ref) + max2d_above=0.5, + max2d_below=0.5, ) - - # selection des points de veget basse proche de la veget haute + # 1.4 selection des points de veget basse proche de la veget haute pipeline = macro.add_radius_assign( pipeline, 1, @@ -173,17 +162,19 @@ def define_marking_pipeline(input_las, output_las, dsm_dimension, dtm_dimension) condition_ref="Classification==5", condition_out="PT_VEG_DSM=1", ) - + # 1.5 Premiers points tagués pour le MNS # max des points de veget (PT_VEG_DSM==1) sur une grille régulière : # TODO: remplacer par GridDecimation une fois le correctif mergé dans PDAL pipeline |= pdal.Filter.grid_decimation_deprecated( resolution=0.75, output_dimension=dsm_dimension, output_type="max", where="PT_VEG_DSM==1" ) + ################################################################################################################### + # 2 - Gestion de l'eau + ################################################################################################################### + # L'eau sous la roche en dévers (typique Corse) + # Gestion de l'eau sur les masques hydro qui ne doivent pas se supperposer aux points virtuels 66 - # 2 - sélection des points pour DTM et DSM - - ########################################################################################################################################### - # Voir pour l'eau sous le sol dans les dévers corse + # 2.1 L'eau sous la roche pipeline = macro.add_radius_assign( pipeline, 1.25, @@ -204,29 +195,67 @@ def define_marking_pipeline(input_las, output_las, dsm_dimension, dtm_dimension) max2d_above=0.5, max2d_below=0.5, ) - ########################################################################################################################################### + # 2.2 Gestion de l'eau sur les masques hydro + pipeline = macro.add_radius_assign( + pipeline, + 1, + False, + condition_src="Classification==9", + condition_ref="Classification==66", + condition_out="PT_ON_VIRT=1", + ) + pipeline = macro.add_radius_assign( + pipeline, + 1, + False, + condition_src="PT_ON_VIRT==1 && Classification==9", + condition_ref="PT_ON_VIRT==0 && Classification==9", + condition_out="PT_ON_VIRT=0", + ) + ################################################################################################################### + # 3 - sélection des premiers points pour MNT et MNS + ################################################################################################################### + # Initialisation pour le MNT + # Initialisation pour le MNS - # selection de points DTM (max) sur une grille régulière + # 3.1 Pour le MNT (le point sol max sur uen grille de 50cm) # TODO: remplacer par GridDecimation une fois le correctif mergé dans PDAL pipeline |= pdal.Filter.grid_decimation_deprecated( resolution=0.5, output_dimension=dtm_dimension, output_type="max", - where="(Classification==2 || PT_ON_SOL==0 && Classification==9)", + where="(Classification==2)", ) - - # selection de points DSM (max) sur une grille régulière + # 3.2 Pour les MNS (Pour le moment: Les bâtis, ponts, veget) # TODO: remplacer par GridDecimation une fois le correctif mergé dans PDAL pipeline |= pdal.Filter.grid_decimation_deprecated( resolution=0.5, output_dimension=dsm_dimension, output_type="max", - where="(PT_ON_VEGET==0 && " - + macro.build_condition("Classification", [6, 9, 17, 64]) - + f" || {dsm_dimension}==1)", + where="(PT_ON_VEGET==0 && (" + + macro.build_condition("Classification", [6, 17, 67]) + + f") || {dsm_dimension}==1)", + ) + # 3.3 Pour les points "eau" on prendra le point le plus bas de la grille de 50cm et qui ne sont ni sous la roche ni près de pts virtuels + pipeline |= pdal.Filter.grid_decimation_deprecated( + resolution=0.5, + output_dimension=dtm_dimension, + output_type="min", + where="(PT_ON_SOL==0 && PT_ON_VIRT==0 && Classification==9)", ) + pipeline |= pdal.Filter.grid_decimation_deprecated( + resolution=0.5, + output_dimension=dsm_dimension, + output_type="min", + where="(PT_ON_VEGET==0 && PT_ON_SOL==0 && PT_ON_VIRT==0 && Classification==9)", + ) + ################################################################################################################### + # 4 - Gestion des points sol sous la veget,bâtis et ponts pour le MNS + ################################################################################################################### + # On enlève les points sol sous la véget, le bati et les ponts du taguage pour les MNS + # Particularité de reprise des points sol au plus près des bâtiments - # assigne des points sol sélectionnés : les points proches de la végétation, des ponts, de l'eau, 67 + # 4.1 Isolement des points sols sous la véget, le bati et les ponts pipeline = macro.add_radius_assign( pipeline, 1.5, @@ -235,16 +264,15 @@ def define_marking_pipeline(input_las, output_las, dsm_dimension, dtm_dimension) condition_ref=macro.build_condition("Classification", [4, 5, 6, 17, 67]), condition_out=f"{dsm_dimension}=0", ) - # Test proximité batiment + # 4.2 Particularité de reprise des points sol au plus près des bâtiments pipeline = macro.add_radius_assign( pipeline, 1.25, False, condition_src="Classification==2 && PT_VEG_DSM==0", - condition_ref="Classification==6", + condition_ref=macro.build_condition("Classification", [6, 67]), condition_out="PT_ON_BUILDING=1", ) - # BUFFER INVERSE Se mettre pipeline = macro.add_radius_assign( pipeline, 1, @@ -253,8 +281,10 @@ def define_marking_pipeline(input_las, output_las, dsm_dimension, dtm_dimension) condition_ref="Classification==2 && PT_ON_BUILDING==0 && PT_VEG_DSM==0", condition_out=f"{dsm_dimension}=1", ) - # 3 - gestion des ponts - # bouche trou : on filtre les points au milieu du pont en les mettant à PT_ON_BRIDGE=1 + ################################################################################################################### + # 5 - Gestion des classes sous les ponts pour être détaguées pour le MNS dsm_dimension=0 + ################################################################################################################### + pipeline = macro.add_radius_assign( pipeline, 1.5, @@ -262,7 +292,7 @@ def define_marking_pipeline(input_las, output_las, dsm_dimension, dtm_dimension) condition_src=macro.build_condition("Classification", [2, 3, 4, 5, 6, 9, 67]), condition_ref="Classification==17", condition_out="PT_ON_BRIDGE=1", - max2d_above=0, # ne pas prendre les points qui sont au dessus des points pont (condition_ref) + max2d_above=0, # ne pas prendre les points qui sont au dessus des points pont (condition_ref) max2d_below=-1, # prendre tous les points qui sont en dessous des points pont (condition_ref) ) pipeline = macro.add_radius_assign( @@ -277,20 +307,63 @@ def define_marking_pipeline(input_las, output_las, dsm_dimension, dtm_dimension) max2d_below=0.5, # prendre tous les points qui sont en dessous des points pont (condition_ref) ) pipeline |= pdal.Filter.assign(value=[f"{dsm_dimension}=0 WHERE PT_ON_BRIDGE==1"]) + ################################################################################################################### + # 6 - Ajout des point pour MNT (sol) qui servent au MNS également + ################################################################################################################### - # 4 - point pour DTM servent au DSM également - # HOMOGENEISER L UTILISATION DE PT_VEG_DSM POUR LES POINT SOL SOUS VEGET AVEC PT_ON_VEGET pipeline |= pdal.Filter.assign( value=[ f"{dsm_dimension}=1 WHERE ({dtm_dimension}==1 && PT_VEG_DSM==0 && PT_ON_BRIDGE==0 && PT_ON_BUILDING==0 && PT_ON_VEGET==0)" ] ) - # 5 - Ajout de la classe 66 pts virtuels dans DTM et DSM - pipeline |= pdal.Filter.assign(value=[f"{dtm_dimension}=1 WHERE (Classification==66)"]) + ################################################################################################################### + # 7 - Gestion des pts virtuels 66 pont et eau / MNT et MNS + ################################################################################################################### + # Taguage pour les MNT de tous les points virtuels + # Gestion des pts virtuels qui sont sous la végétation ou autres pour le MNS + # Taguage pour les MNS des points vituels "eau" seulement - # 6 - export du nuage et des DSM - pipeline |= pdal.Writer.las(extra_dims="all", forward="all", filename=output_las) + # 7.1 Taguade pour les MNT des points virtuels ponts et eau + pipeline |= pdal.Filter.assign(value=[f"{dtm_dimension}=1 WHERE Classification==66"]) + # 7.2 gestion des pts 66 "eau" sous le sursol + pipeline = macro.add_radius_assign( + pipeline, + 0.5, + False, + condition_src="Classification==66", + condition_ref=macro.build_condition("Classification", [4, 5, 6, 17, 67]), + condition_out="PT_ON_VEGET=1", + ) + pipeline = macro.add_radius_assign( + pipeline, + 0.5, + False, + condition_src="PT_ON_VEGET==1 && Classification==66", + condition_ref="PT_ON_VEGET==0 && Classification==66", + condition_out="PT_ON_VEGET=0", + ) + # 7.3 Taguage pour les MNS des points virtuels eau seulement + pipeline = macro.add_radius_assign( + pipeline, + 0.5, + False, + condition_src="Classification==66", + condition_ref="Classification==17", + condition_out="PT_ON_BRIDGE=1", + ) + pipeline |= pdal.Filter.assign( + value=[ + f"{dsm_dimension}=1 WHERE (Classification==66 && PT_ON_VEGET==0 && PT_ON_BRIDGE==0)" + ] + ) + + ################################################################################################################## + # 8 - export du nuage et des DSM + ################################################################################################################### + pipeline |= pdal.Writer.las( + extra_dims="all", forward="all", filename=output_las, minor_version="4" + ) return pipeline, temporary_dimensions diff --git a/pdal_ign_macro/version.py b/pdal_ign_macro/version.py index 19f10f4..74bd34b 100755 --- a/pdal_ign_macro/version.py +++ b/pdal_ign_macro/version.py @@ -1,4 +1,4 @@ -__version__ = "0.2.1" +__version__ = "0.3.0" if __name__ == "__main__": diff --git a/src/filter_radius_assign/RadiusAssignFilter.cpp b/src/filter_radius_assign/RadiusAssignFilter.cpp index e670956..4cbb28a 100644 --- a/src/filter_radius_assign/RadiusAssignFilter.cpp +++ b/src/filter_radius_assign/RadiusAssignFilter.cpp @@ -83,12 +83,12 @@ void RadiusAssignFilter::doOneNoDomain(PointRef &point) double Zref = point.getFieldAs(Dimension::Id::Z); if (m_args->m_max2d_below>=0 || m_args->m_max2d_above>=0) { - bool take (false); + bool take (true); 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 (m_args->m_max2d_below>=0 && (Zref-Zpt)>m_args->m_max2d_below) {take=false; break;} + if (m_args->m_max2d_above>=0 && (Zpt-Zref)>m_args->m_max2d_above) {take=false; break;} } if (!take) return; } diff --git a/test/data/mnx/input/crop_1.laz b/test/data/mnx/input/crop_1.laz new file mode 100644 index 0000000..3478412 Binary files /dev/null and b/test/data/mnx/input/crop_1.laz differ diff --git a/test/data/mnx/input/crop_2.laz b/test/data/mnx/input/crop_2.laz new file mode 100644 index 0000000..c070343 Binary files /dev/null and b/test/data/mnx/input/crop_2.laz differ diff --git a/test/data/mnx/input/crop_3.laz b/test/data/mnx/input/crop_3.laz new file mode 100644 index 0000000..d7348b7 Binary files /dev/null and b/test/data/mnx/input/crop_3.laz differ diff --git a/test/data/mnx/reference/crop_1.laz b/test/data/mnx/reference/crop_1.laz new file mode 100644 index 0000000..38048ad Binary files /dev/null and b/test/data/mnx/reference/crop_1.laz differ diff --git a/test/data/mnx/reference/crop_2.laz b/test/data/mnx/reference/crop_2.laz new file mode 100644 index 0000000..44d5fef Binary files /dev/null and b/test/data/mnx/reference/crop_2.laz differ diff --git a/test/data/mnx/reference/crop_3.laz b/test/data/mnx/reference/crop_3.laz new file mode 100644 index 0000000..fa6a21a Binary files /dev/null and b/test/data/mnx/reference/crop_3.laz differ diff --git a/test/pdal_ign_macro/test_mark_points_to_use_for_digital_models_with_new_dimension.py b/test/pdal_ign_macro/test_mark_points_to_use_for_digital_models_with_new_dimension.py index 07ff379..89954e8 100644 --- a/test/pdal_ign_macro/test_mark_points_to_use_for_digital_models_with_new_dimension.py +++ b/test/pdal_ign_macro/test_mark_points_to_use_for_digital_models_with_new_dimension.py @@ -3,6 +3,7 @@ import numpy as np import pdal +import pytest from pdal_ign_macro.mark_points_to_use_for_digital_models_with_new_dimension import ( main, @@ -131,3 +132,61 @@ def test_parse_args(): parsed_args_keys = args.__dict__.keys() main_parameters = inspect.signature(main).parameters.keys() assert parsed_args_keys == main_parameters + + +@pytest.mark.parametrize( + "crop", + [ + "crop_1.laz", + # "crop_2.laz", ToDo : rebuild the reference for crop_2 which is false + "crop_3.laz", + ], +) +def test_algo_mark_points_for_dm_with_reference(crop): + ini_las = "test/data/mnx/input/" + crop + dsm_dimension = "dsm_marker" + dtm_dimension = "dtm_marker" + with tempfile.NamedTemporaryFile(suffix="_after.las") as las_output: + + main( + ini_las, + las_output.name, + dsm_dimension, + dtm_dimension, + "", + "", + keep_temporary_dims=False, + skip_buffer=True, + ) + + def sort_points(points): + sorted_index = np.lexsort((points["Y"], points["X"], points["Z"], points["GpsTime"])) + np.set_printoptions(precision=30) + sorted_points = points[sorted_index] + return sorted_points + + pipeline_calc = pdal.Pipeline() | pdal.Reader.las(filename=las_output.name) + pipeline_calc.execute() + arr_result = pipeline_calc.arrays[0] + + # read reference : + ref_las = "test/data/mnx/reference/" + crop + + pipeline_ref = pdal.Pipeline() | pdal.Reader.las(filename=ref_las) + pipeline_ref.execute() + arr_reference = pipeline_ref.arrays[0] + + assert len(arr_result) == len(arr_reference) + + arr_result = sort_points(arr_result) + arr_reference = sort_points(arr_reference) + + arr_result_dimensions = list(arr_result.dtype.fields.keys()) + arr_ref_dimensions = list(arr_reference.dtype.fields.keys()) + + assert arr_result_dimensions == arr_ref_dimensions + + for dim in arr_result_dimensions: + diff_mask = np.where(arr_result[dim] == arr_reference[dim], 0, 1) + nb_pts_incorrect = np.count_nonzero(diff_mask) + assert nb_pts_incorrect == 0 diff --git a/test/test_radius_assign.py b/test/test_radius_assign.py index 7f35205..0a7c1bd 100755 --- a/test/test_radius_assign.py +++ b/test/test_radius_assign.py @@ -14,6 +14,7 @@ pt_ini = (pt_x, pt_y, pt_z, 1) numeric_precision = 4 +distance_radius = 1 def distance2d(pt1, pt2): @@ -74,7 +75,7 @@ def run_filter(arrays_las, distance_radius, search_3d, limit_z_above=-1, limit_w return nb_pts_radius_search -def build_random_points_around_one_point(test_function, distance_radius): +def build_random_points_around_one_point(test_function, points=[]): dtype = [("X", "= 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 + if (limit_z_above >= 0) and ((pt_ini[2] - pt[2]) > limit_z_above): + return 0 + if (limit_z_below >= 0) and ((pt[2] - pt_ini[2]) > limit_z_below): + return 0 + return 1 + else: + return 0 + + arrays_las, nb_points_take_2d = build_random_points_around_one_point(func_test) + + assert len(arrays_las) > 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, limit_z_above, limit_z_below ) + assert nb_pts_radius_2d_cylinder == nb_points_take_2d