Update radius (#25)

* update ccode radius assign - cylinder * update script MNx * update tests MNx * update tests * update algo radius, update some variable name. * Update changelog --------- Co-authored-by: Lea Vauchier <lea.vauchier@ign.fr>
IGNF · Oct 16, 2024 · 59c36bd · 59c36bd
1 parent ea63dbf
commit 59c36bd
Show file tree

Hide file tree

Showing 15 changed files with 83 additions and 73 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
diff --git a/RELEASE_NOTES.md b/RELEASE_NOTES.md
@@ -1,7 +1,11 @@
-### 0.3.0
+- improve code readability in the radius_assign filter (Z limits part).
+- fix the script for MNx pre-processing (Z limits were inverted)
+- improve tests and test data
+
+## 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 : 
+In details :
 - manage water and virtuals points,
 - update building consideration
 - 2 levels of water

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
@@ -88,9 +88,9 @@ def define_marking_pipeline(input_las, output_las, dsm_dimension, dtm_dimension)
     # 0 - ajout de dimensions temporaires et de sortie
     temporary_dimensions = [
         "PT_VEG_DSM",
-        "PT_ON_BRIDGE",
-        "PT_ON_BUILDING",
-        "PT_ON_VEGET",
+        "PT_UNDER_BRIDGE",
+        "PT_CLOSED_BUILDING",
+        "PT_UNDER_VEGET",
         "PT_ON_SOL",
         "PT_ON_VIRT",
     ]
@@ -128,28 +128,28 @@ 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)
+    # 1.3 Isolement en PT_UNDER_VEGET=1 des éléments sous la végétation (hors sol)
     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,
-        max2d_below=-1,
+        condition_out="PT_UNDER_VEGET=1",
+        max2d_above=-1,
+        max2d_below=0,
     )
     pipeline = macro.add_radius_assign(
         pipeline,
         1,
         False,
-        condition_src="PT_ON_VEGET==1 && ( "
+        condition_src="PT_UNDER_VEGET==1 && ( "
         + macro.build_condition("Classification", [6, 9, 17, 67])
         + " )",
-        condition_ref="PT_ON_VEGET==0 && ( "
+        condition_ref="PT_UNDER_VEGET==0 && ( "
         + macro.build_condition("Classification", [6, 9, 17, 67])
         + " )",
-        condition_out="PT_ON_VEGET=0",
+        condition_out="PT_UNDER_VEGET=0",
         max2d_above=0.5,
         max2d_below=0.5,
     )
@@ -182,8 +182,8 @@ def define_marking_pipeline(input_las, output_las, dsm_dimension, dtm_dimension)
         condition_src="Classification==9",
         condition_ref="Classification==2",
         condition_out="PT_ON_SOL=1",
-        max2d_above=0,
-        max2d_below=-1,
+        max2d_above=-1,
+        max2d_below=0,
     )
     pipeline = macro.add_radius_assign(
         pipeline,
@@ -232,7 +232,7 @@ def define_marking_pipeline(input_las, output_las, dsm_dimension, dtm_dimension)
         resolution=0.5,
         output_dimension=dsm_dimension,
         output_type="max",
-        where="(PT_ON_VEGET==0 && ("
+        where="(PT_UNDER_VEGET==0 && ("
         + macro.build_condition("Classification", [6, 17, 67])
         + f") || {dsm_dimension}==1)",
     )
@@ -247,7 +247,7 @@ def define_marking_pipeline(input_las, output_las, dsm_dimension, dtm_dimension)
         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)",
+        where="(PT_UNDER_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
@@ -271,14 +271,14 @@ def define_marking_pipeline(input_las, output_las, dsm_dimension, dtm_dimension)
         False,
         condition_src="Classification==2 && PT_VEG_DSM==0",
         condition_ref=macro.build_condition("Classification", [6, 67]),
-        condition_out="PT_ON_BUILDING=1",
+        condition_out="PT_CLOSED_BUILDING=1",
     )
     pipeline = macro.add_radius_assign(
         pipeline,
         1,
         False,
-        condition_src=f"Classification==2 && {dsm_dimension}==0 && PT_ON_BUILDING==1 && {dtm_dimension}==1",
-        condition_ref="Classification==2 && PT_ON_BUILDING==0 && PT_VEG_DSM==0",
+        condition_src=f"Classification==2 && {dsm_dimension}==0 && PT_CLOSED_BUILDING==1 && {dtm_dimension}==1",
+        condition_ref="Classification==2 && PT_CLOSED_BUILDING==0 && PT_VEG_DSM==0",
         condition_out=f"{dsm_dimension}=1",
     )
     ###################################################################################################################
@@ -291,29 +291,29 @@ def define_marking_pipeline(input_las, output_las, dsm_dimension, dtm_dimension)
         False,
         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_below=-1,  # prendre tous les points qui sont en dessous des points pont (condition_ref)
+        condition_out="PT_UNDER_BRIDGE=1",
+        max2d_above=-1,  # prendre les points (condition_src) qui on des points ponts au dessus d'eux (condition_ref)
+        max2d_below=0,
     )
     pipeline = macro.add_radius_assign(
         pipeline,
         1.25,
         False,
-        condition_src="PT_ON_BRIDGE==1",
-        condition_ref="PT_ON_BRIDGE==0 && "
+        condition_src="PT_UNDER_BRIDGE==1",
+        condition_ref="PT_UNDER_BRIDGE==0 && "
         + macro.build_condition("Classification", [2, 3, 4, 5, 6, 9, 67]),
-        condition_out="PT_ON_BRIDGE=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_out="PT_UNDER_BRIDGE=0",
+        max2d_above=0.5,
+        max2d_below=0.5,
     )
-    pipeline |= pdal.Filter.assign(value=[f"{dsm_dimension}=0 WHERE PT_ON_BRIDGE==1"])
+    pipeline |= pdal.Filter.assign(value=[f"{dsm_dimension}=0 WHERE PT_UNDER_BRIDGE==1"])
     ###################################################################################################################
     # 6 - Ajout des point pour MNT (sol) qui servent au MNS également
     ###################################################################################################################
 
     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)"
+            f"{dsm_dimension}=1 WHERE ({dtm_dimension}==1 && PT_VEG_DSM==0 && PT_UNDER_BRIDGE==0 && PT_CLOSED_BUILDING==0 && PT_UNDER_VEGET==0)"
         ]
     )
 
@@ -333,15 +333,15 @@ def define_marking_pipeline(input_las, output_las, dsm_dimension, dtm_dimension)
         False,
         condition_src="Classification==66",
         condition_ref=macro.build_condition("Classification", [4, 5, 6, 17, 67]),
-        condition_out="PT_ON_VEGET=1",
+        condition_out="PT_UNDER_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",
+        condition_src="PT_UNDER_VEGET==1 && Classification==66",
+        condition_ref="PT_UNDER_VEGET==0 && Classification==66",
+        condition_out="PT_UNDER_VEGET=0",
     )
     # 7.3 Taguage pour les MNS des points virtuels eau seulement
     pipeline = macro.add_radius_assign(
@@ -350,11 +350,11 @@ def define_marking_pipeline(input_las, output_las, dsm_dimension, dtm_dimension)
         False,
         condition_src="Classification==66",
         condition_ref="Classification==17",
-        condition_out="PT_ON_BRIDGE=1",
+        condition_out="PT_UNDER_BRIDGE=1",
     )
     pipeline |= pdal.Filter.assign(
         value=[
-            f"{dsm_dimension}=1 WHERE (Classification==66 && PT_ON_VEGET==0 && PT_ON_BRIDGE==0)"
+            f"{dsm_dimension}=1 WHERE (Classification==66 && PT_UNDER_VEGET==0 && PT_UNDER_BRIDGE==0)"
         ]
     )
 

diff --git a/src/filter_radius_assign/RadiusAssignFilter.cpp b/src/filter_radius_assign/RadiusAssignFilter.cpp
@@ -68,33 +68,34 @@ void RadiusAssignFilter::ready(PointTableRef)
     m_args->m_ptsToUpdate.clear();
 }
 
-void RadiusAssignFilter::doOneNoDomain(PointRef &point)
+void RadiusAssignFilter::doOneNoDomain(PointRef &pointSrc)
 {
     // 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 (m_args->search3d) iNeighbors = refView->build3dIndex().radius(pointSrc, m_args->m_radius);
+    else iNeighbors = refView->build2dIndex().radius(pointSrc, m_args->m_radius);
+    
     if (iNeighbors.size() == 0)
         return;
 
-    if (!m_args->search3d)
+    if (!m_args->search3d && (m_args->m_max2d_below>=0 || m_args->m_max2d_above>=0))
     {
-        double Zref = point.getFieldAs<double>(Dimension::Id::Z);
-        if (m_args->m_max2d_below>=0 || m_args->m_max2d_above>=0)
+        double Zsrc = pointSrc.getFieldAs<double>(Dimension::Id::Z);
+
+        bool take (false);
+        for (PointId ptId : iNeighbors)
         {
-            bool take (true);
-            for (PointId ptId : iNeighbors)
-            {
-                double Zpt = refView->point(ptId).getFieldAs<double>(Dimension::Id::Z);
-                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;
+            double Zref = refView->point(ptId).getFieldAs<double>(Dimension::Id::Z);
+            if (m_args->m_max2d_above>=0 && Zref>Zsrc && (Zref-Zsrc)>m_args->m_max2d_above) continue;
+            if (m_args->m_max2d_below>=0 && Zsrc>Zref && (Zsrc-Zref)>m_args->m_max2d_below) continue;
+            take = true;
+            break;
         }
+
+        if (!take) return;
     }
 
-    m_args->m_ptsToUpdate.push_back(point.pointId());
+    m_args->m_ptsToUpdate.push_back(pointSrc.pointId());
 }
 
 bool RadiusAssignFilter::doOne(PointRef& point)

diff --git a/test/data/mnx/input/bat.laz b/test/data/mnx/input/bat.laz
diff --git a/test/data/mnx/input/corse.laz b/test/data/mnx/input/corse.laz
diff --git a/test/data/mnx/input/pont.laz b/test/data/mnx/input/pont.laz
diff --git a/test/data/mnx/reference/bat.laz b/test/data/mnx/reference/bat.laz
diff --git a/test/data/mnx/reference/corse.laz b/test/data/mnx/reference/corse.laz
diff --git a/test/data/mnx/reference/crop_1.laz b/test/data/mnx/reference/crop_1.laz
diff --git a/test/data/mnx/reference/crop_2.laz b/test/data/mnx/reference/crop_2.laz
diff --git a/test/data/mnx/reference/crop_3.laz b/test/data/mnx/reference/crop_3.laz
diff --git a/test/data/mnx/reference/pont.laz b/test/data/mnx/reference/pont.laz
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
@@ -57,7 +57,12 @@ def test_mark_points_to_use_for_digital_models_with_new_dimension_keep_dimension
         assert all(
             [
                 dim in output_dimensions
-                for dim in ["PT_VEG_DSM", "PT_ON_BRIDGE", "PT_ON_BUILDING", "PT_ON_VEGET"]
+                for dim in [
+                    "PT_VEG_DSM",
+                    "PT_UNDER_BRIDGE",
+                    "PT_CLOSED_BUILDING",
+                    "PT_UNDER_VEGET",
+                ]
             ]
         )
 
@@ -140,6 +145,9 @@ def test_parse_args():
         "crop_1.laz",
         # "crop_2.laz", ToDo : rebuild the reference for crop_2 which is false
         "crop_3.laz",
+        "bat.laz",
+        "pont.laz",
+        "corse.laz",
     ],
 )
 def test_algo_mark_points_for_dm_with_reference(crop):

diff --git 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
+numeric_precision_z = 2
 distance_radius = 1
 
 
@@ -28,6 +29,10 @@ def distance3d(pt1, pt2):
     )
 
 
+def distanceZ(pt1, pt2):
+    return round(pt1[2] - pt2[2], numeric_precision_z)
+
+
 def run_filter(arrays_las, distance_radius, search_3d, limit_z_above=-1, limit_w_below=-1):
 
     filter = "filters.radius_assign"
@@ -154,7 +159,7 @@ def test_radius_assign_2d_cylinder_below():
 
     def func_test(pt):
         distance_i = distance2d(pt_ini, pt)
-        distance_z = pt[2] - pt_ini[2]
+        distance_z = distanceZ(pt, pt_ini)
         if distance_i < distance_radius and distance_z < limit_z_below:
             return 1
         return 0
@@ -177,7 +182,7 @@ def test_radius_assign_2d_cylinder_above():
 
     def func_test(pt):
         distance_i = distance2d(pt_ini, pt)
-        distance_z = pt_ini[2] - pt[2]
+        distance_z = distanceZ(pt_ini, pt)
         if distance_i < distance_radius and distance_z < limit_z_above:
             return 1
         return 0
@@ -197,7 +202,7 @@ def test_radius_assign_2d_cylinder_above_below_null():
 
     def func_test(pt):
         distance_i = distance2d(pt_ini, pt)
-        distance_z = pt_ini[2] - pt[2]
+        distance_z = distanceZ(pt_ini, pt)
         if distance_i < distance_radius and distance_z == 0:
             return 1
         return 0
@@ -217,7 +222,7 @@ def test_radius_assign_2d_cylinder_above_null_bellow_all():
 
     def func_test(pt):
         distance_i = distance2d(pt_ini, pt)
-        distance_z = pt[2] - pt_ini[2]
+        distance_z = distanceZ(pt, pt_ini)
         if distance_i < distance_radius and distance_z <= 0:
             return 1
         return 0
@@ -261,13 +266,12 @@ def test_radius_assign_2d_cylinder_above_and_bellow(execution_number):
 
     def func_test(pt):
         distance_i = distance2d(pt_ini, pt)
-        distance_z = pt_ini[2] - pt[2]
-        if (
-            distance_i < distance_radius
-            and distance_z <= limit_z_above
-            and (-distance_z) <= limit_z_below
-        ):
-            return 1
+        distance_z = distanceZ(pt, pt_ini)  # src - ref
+        if distance_i < distance_radius:
+            if distance_z <= 0 and (-distance_z) <= limit_z_above:  # src est sur ref
+                return 1
+            if distance_z >= 0 and distance_z <= limit_z_below:  # src est sous ref
+                return 1
         return 0
 
     arrays_las, nb_points_take_2d = build_random_points_around_one_point(func_test, points)
@@ -298,9 +302,10 @@ def test_radius_assign_2d_cylinder(limit_z_above, limit_z_below):
     def func_test(pt):
         distance_i = distance2d(pt_ini, pt)
         if distance_i < distance_radius:
-            if (limit_z_above >= 0) and ((pt_ini[2] - pt[2]) > limit_z_above):
+            distance_z = distanceZ(pt, pt_ini)  # src - ref
+            if limit_z_above >= 0 and distance_z <= 0 and (-distance_z) > limit_z_above:
                 return 0
-            if (limit_z_below >= 0) and ((pt[2] - pt_ini[2]) > limit_z_below):
+            if limit_z_below >= 0 and distance_z >= 0 and distance_z > limit_z_below:
                 return 0
             return 1
         else: