diff --git a/src/filter_grid_decimation/GridDecimationFilter.cpp b/src/filter_grid_decimation/GridDecimationFilter.cpp index 1366dcb..366caab 100755 --- a/src/filter_grid_decimation/GridDecimationFilter.cpp +++ b/src/filter_grid_decimation/GridDecimationFilter.cpp @@ -79,15 +79,22 @@ void GridDecimationFilter::processOne(BOX2D bounds, PointRef& point, PointViewPt double x = point.getFieldAs(Dimension::Id::X); double y = point.getFieldAs(Dimension::Id::Y); int id = point.getFieldAs(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(d_width_pt); - int height = static_cast(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((x - bounds.minx) / m_args->m_edgeLength); + int height = static_cast((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::max)()) throwError("Grid width out of range."); if (d_height < 0.0 || d_height > (std::numeric_limits::max)()) @@ -123,7 +130,7 @@ void GridDecimationFilter::createGrid(BOX2D bounds) for (size_t l(0); lm_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); diff --git a/test/data/4_6.las b/test/data/4_6.las old mode 100755 new mode 100644 index 0717b32..b69ac5c Binary files a/test/data/4_6.las and b/test/data/4_6.las differ diff --git a/test/test_grid_decimation.py b/test/test_grid_decimation.py index 50aa01f..8cf4563 100755 --- a/test/test_grid_decimation.py +++ b/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,9 +52,42 @@ 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: @@ -64,16 +95,20 @@ def run_filter(type): 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)