maj du grid_decimation

src/filter_grid_decimation/GridDecimationFilter.cpp

@@ -79,15 +79,19 @@ void GridDecimationFilter::processOne(BOX2D bounds, PointRef& point, PointViewPt
     double x = point.getFieldAs<double>(Dimension::Id::X);
     double y = point.getFieldAs<double>(Dimension::Id::Y);
     int id = point.getFieldAs<double>(Dimension::Id::PointId);
+
+    // if x==xmax of the cell, the point are in the bottom cell
+    // if y==ymax of the cell, the point are in the right cell
 
-    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);
+    double d_width_pt = (x - bounds.minx) / m_args->m_edgeLength;
+    double d_height_pt = (y - bounds.miny) / m_args->m_edgeLength;
 
     int width = static_cast<int>(d_width_pt);
     int height = static_cast<int>(d_height_pt);
 
-    auto ptRefid = this->grid[ std::make_pair(width,height) ];
-
+    auto mptRefid = this->grid.find( std::make_pair(width,height) );
+    auto ptRefid = mptRefid->second;
+
     if (ptRefid==-1)
     {
         this->grid[ std::make_pair(width,height) ] = point.pointId();
@@ -107,8 +111,8 @@ void GridDecimationFilter::processOne(BOX2D bounds, PointRef& point, PointViewPt
 
 void GridDecimationFilter::createGrid(BOX2D bounds)
 {
-    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);
+    double d_width = std::round((bounds.maxx - bounds.minx) / m_args->m_edgeLength);
+    double d_height = std::round((bounds.maxy - bounds.miny) / m_args->m_edgeLength);
 
     if (d_width < 0.0 || d_width > (std::numeric_limits<int>::max)())
         throwError("Grid width out of range.");
@@ -123,7 +127,7 @@ void GridDecimationFilter::createGrid(BOX2D bounds)
     for (size_t l(0); l<height; l++)
         for (size_t c(0); c<width; c++)
         {
-            BOX2D bounds_dalle ( bounds.minx + c*m_args->m_edgeLength, bounds.miny + l*m_args->m_edgeLength,
+            BOX2D bounds_dalle (bounds.minx + c*m_args->m_edgeLength, bounds.miny + l*m_args->m_edgeLength,
                                 bounds.minx + (c+1)*m_args->m_edgeLength, bounds.miny + (l+1)*m_args->m_edgeLength );
             vgrid.push_back(Polygon(bounds_dalle));
             this->grid.insert( std::make_pair( std::make_pair(c,l), -1)  );
@@ -143,7 +147,7 @@ PointViewSet GridDecimationFilter::run(PointViewPtr view)
     BOX2D bounds;
     view->calculateBounds(bounds);
     createGrid(bounds);
-    
+
     for (PointId i = 0; i < view->size(); ++i)
     {
         PointRef point = view->point(i);

test/test_grid_decimation.py

@@ -1,13 +1,73 @@
 import csv
 import json
 import math
+import re
 import tempfile
 from test import utils
 
 import pdal
 import pdaltools.las_info as li
 
 
+def parse_geometry(geometry):
+    regex = r"[0-9-\.]+"
+    parsed_geom = re.findall(regex, geometry)
+    parsed_geom = [float(i) for i in parsed_geom]
+    return (
+        max(parsed_geom[::2]),
+        min(parsed_geom[::2]),
+        max(parsed_geom[1::2]),
+        min(parsed_geom[1::2]),
+    )
+
+
+def read_las_crop(las, line, type):
+
+    tmp_out_las = tempfile.NamedTemporaryFile(suffix=".las").name
+
+    PIPELINE = [
+        {"type": "readers.las", "filename": las, "extra_dims": "grid=uint8"},
+        {
+            "type": "filters.crop",
+            "polygon": line,
+        },
+        {
+            "type": "writers.las",
+            "extra_dims": "all",
+            "filename": tmp_out_las,
+        },
+    ]
+
+    pipeline = pdal.Pipeline(json.dumps(PIPELINE))
+
+    # execute the pipeline
+    pipeline.execute()
+    arrays = pipeline.arrays
+    array = arrays[0]
+
+    max_x, min_x, max_y, min_y = parse_geometry(line)
+
+    ZGrid = 0
+    if type == "max":
+        Z = -9999
+    else:
+        Z = 9999
+    for pt in array:
+        if type == "max" and pt["Z"] > Z:
+            Z = pt["Z"]
+        if type == "min" and pt["Z"] < Z:
+            Z = pt["Z"]
+        if pt["grid"] > 0:
+            # if the point is exactly on the bbox at xmax or ymax, it's one of another cell
+            if pt["X"] == max_x:
+                continue
+            if pt["Y"] == max_y:
+                continue
+            ZGrid = pt["Z"]
+
+    assert ZGrid == Z
+
+
 def run_filter(type):
 
     ini_las = "test/data/4_6.las"
@@ -21,8 +81,8 @@ def run_filter(type):
 
     bounds = li.las_get_xy_bounds(ini_las)
 
-    d_width = math.floor((bounds[0][1] - bounds[0][0]) / resolution) + 1
-    d_height = math.floor((bounds[1][1] - bounds[1][0]) / resolution) + 1
+    d_width = math.ceil((bounds[0][1] - bounds[0][0]) / resolution)
+    d_height = math.ceil((bounds[1][1] - bounds[1][0]) / resolution)
     nb_dalle = d_width * d_height
 
     PIPELINE = [
@@ -38,7 +98,6 @@ def run_filter(type):
             "type": "writers.las",
             "extra_dims": "all",
             "filename": tmp_out_las,
-            "where": "grid==1",
         },
     ]
 
@@ -54,26 +113,21 @@ def run_filter(type):
         if pt["grid"] > 0:
             nb_pts_grid += 1
 
-    # fix dans une autre branche en cours
-    # assert nb_pts_grid == 3231
-    assert nb_pts_grid <= nb_dalle
+    assert nb_pts_grid == nb_dalle
 
     data = []
     with open(tmp_out_wkt, "r") as f:
         reader = csv.reader(f, delimiter="\t")
         for i, line in enumerate(reader):
             data.append(line[0])
+            read_las_crop(tmp_out_las, line[0], type)
 
-    # fix dans une autre branche en cours
-    # assert len(data) == nb_dalle
-
-    return nb_pts_grid
+    assert len(data) == nb_dalle
 
 
 def test_grid_decimation_max():
     run_filter("max")
 
 
-# fix dans une autre branche en cours
-# def test_grid_decimation_min():
-#    run_filter("min")
+def test_grid_decimation_min():
+    run_filter("min")