Update grid decimation

src/filter_grid_decimation/GridDecimationFilter.cpp

@@ -79,15 +79,22 @@ 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);
-
-    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<int>(d_width_pt);
-    int height = static_cast<int>(d_height_pt);
-
-    auto ptRefid = this->grid[ std::make_pair(width,height) ];
-
+
+    // if x==(xmax of the cell), we assume the point are in the upper cell
+    // if y==(ymax of the cell), we assume the point are in the right cell
+    int width = static_cast<int>((x - bounds.minx) / m_args->m_edgeLength);
+    int height = static_cast<int>((y - bounds.miny) / m_args->m_edgeLength);
+
+    // to avoid numeric pb with the division (append if the point is on the grid)
+    if (x < bounds.minx+width*m_args->m_edgeLength) width--;
+    if (y < bounds.miny+height*m_args->m_edgeLength) height--;
+    if (x >= bounds.minx+(width+1)*m_args->m_edgeLength) width++;
+    if (y >= bounds.miny+(height+1)*m_args->m_edgeLength) height++;
+
+    auto mptRefid = this->grid.find( std::make_pair(width,height) );
+    assert( mptRefid != this->grid.end() );
+    auto ptRefid = mptRefid->second;
+
     if (ptRefid==-1)
     {
         this->grid[ std::make_pair(width,height) ] = point.pointId();
@@ -107,9 +114,9 @@ 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);
-    
+    size_t d_width = std::floor((bounds.maxx - bounds.minx) / m_args->m_edgeLength) + 1;
+    size_t d_height = std::floor((bounds.maxy - bounds.miny) / m_args->m_edgeLength) + 1;
+
     if (d_width < 0.0 || d_width > (std::numeric_limits<int>::max)())
         throwError("Grid width out of range.");
     if (d_height < 0.0 || d_height > (std::numeric_limits<int>::max)())
@@ -123,7 +130,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 +150,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

@@ -6,15 +6,19 @@
 
 import pdal
 import pdaltools.las_info as li
+import pytest
 
 
-def run_filter(type):
+def contains(bounds, x, y):
+    # to be coherent with the grid decimation algorithm
+    return bounds[0] <= x and x < bounds[1] and bounds[2] <= y and y < bounds[3]
+
+
+def run_filter(type, resolution):
 
     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
+    tmp_out_wkt = tempfile.NamedTemporaryFile(suffix=f"_{resolution}.wkt").name
 
     filter = "filters.grid_decimation_deprecated"
     utils.pdal_has_plugin(filter)
@@ -34,12 +38,6 @@ def run_filter(type):
             "output_dimension": "grid",
             "output_wkt": tmp_out_wkt,
         },
-        {
-            "type": "writers.las",
-            "extra_dims": "all",
-            "filename": tmp_out_las,
-            "where": "grid==1",
-        },
     ]
 
     pipeline = pdal.Pipeline(json.dumps(PIPELINE))
@@ -54,26 +52,63 @@ 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
+
+    for lig in range(d_height):
+        for col in range(d_width):
+
+            cell = [
+                bounds[0][0] + col * resolution,
+                bounds[0][0] + (col + 1) * resolution,
+                bounds[1][0] + lig * resolution,
+                bounds[1][0] + (lig + 1) * resolution,
+            ]
+
+            nbThreadPtsCrop = 0
+            ZRef = 0
+            ZRefGrid = 0
+
+            for pt in array:
+                x = pt["X"]
+                y = pt["Y"]
+                if not contains(cell, x, y):
+                    continue
+
+                z = pt["Z"]
+                if type == "max":
+                    if ZRef == 0 or z > ZRef:
+                        ZRef = z
+                elif type == "min":
+                    if ZRef == 0 or z < ZRef:
+                        ZRef = z
+
+                if pt["grid"] > 0:
+                    nbThreadPtsCrop += 1
+                    ZRefGrid = z
+
+            assert nbThreadPtsCrop == 1
+            assert ZRef == ZRefGrid
 
     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])
 
-    # 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")
+@pytest.mark.parametrize(
+    "resolution",
+    [(10.1), (10.0), (9.8)],
+)
+def test_grid_decimation_max(resolution):
+    run_filter("max", resolution)
 
 
-# fix dans une autre branche en cours
-# def test_grid_decimation_min():
-#    run_filter("min")
+@pytest.mark.parametrize(
+    "resolution",
+    [(10.3), (10.0), (9.9)],
+)
+def test_grid_decimation_min(resolution):
+    run_filter("min", resolution)