From d2fcb2f1ab93ef555d9e07206282e13d50970edd Mon Sep 17 00:00:00 2001 From: Antoine Lavenant Date: Fri, 7 Jun 2024 15:08:19 +0200 Subject: [PATCH] fix code + test --- .../GridDecimationFilter.cpp | 21 +-- test/test_grid_decimation.py | 139 ++++++++---------- 2 files changed, 72 insertions(+), 88 deletions(-) diff --git a/src/filter_grid_decimation/GridDecimationFilter.cpp b/src/filter_grid_decimation/GridDecimationFilter.cpp index f01fab6..366caab 100755 --- a/src/filter_grid_decimation/GridDecimationFilter.cpp +++ b/src/filter_grid_decimation/GridDecimationFilter.cpp @@ -80,13 +80,16 @@ void GridDecimationFilter::processOne(BOX2D bounds, PointRef& point, PointViewPt double y = point.getFieldAs(Dimension::Id::Y); int id = point.getFieldAs(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 = (x - bounds.minx) / m_args->m_edgeLength; - double d_height_pt = (y - bounds.miny) / m_args->m_edgeLength; + // 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); - int width = static_cast(d_width_pt); - int height = static_cast(d_height_pt); + // 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() ); @@ -111,9 +114,9 @@ void GridDecimationFilter::processOne(BOX2D bounds, PointRef& point, PointViewPt void GridDecimationFilter::createGrid(BOX2D bounds) { - 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); - + 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)()) diff --git a/test/test_grid_decimation.py b/test/test_grid_decimation.py index f2f76a3..2e01d08 100755 --- a/test/test_grid_decimation.py +++ b/test/test_grid_decimation.py @@ -1,88 +1,30 @@ import csv import json import math -import re import tempfile from test import utils +import pytest import pdal import pdaltools.las_info as li +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 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): +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) bounds = li.las_get_xy_bounds(ini_las) - d_width = math.ceil((bounds[0][1] - bounds[0][0]) / resolution) - d_height = math.ceil((bounds[1][1] - bounds[1][0]) / resolution) + d_width = math.floor((bounds[0][1] - bounds[0][0]) / resolution) + 1 + d_height = math.floor((bounds[1][1] - bounds[1][0]) / resolution) + 1 nb_dalle = d_width * d_height PIPELINE = [ @@ -94,11 +36,6 @@ def run_filter(type): "output_dimension": "grid", "output_wkt": tmp_out_wkt, }, - { - "type": "writers.las", - "extra_dims": "all", - "filename": tmp_out_las, - }, ] pipeline = pdal.Pipeline(json.dumps(PIPELINE)) @@ -115,19 +52,63 @@ def run_filter(type): assert nb_pts_grid == nb_dalle + for l in range(d_height): + for c in range(d_width): + + cell = [bounds[0][0] + c*resolution, bounds[0][0] + (c+1)*resolution, + bounds[1][0] + l*resolution, bounds[1][0] + (l+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]) - read_las_crop(tmp_out_las, line[0], type) assert len(data) == nb_dalle - -def test_grid_decimation_max(): - run_filter("max") - - -def test_grid_decimation_min(): - run_filter("min") +@pytest.mark.parametrize( + "resolution", + [ + (10.1), + (10.), + (9.8) + ], +) +def test_grid_decimation_max(resolution): + run_filter("max", resolution) + +@pytest.mark.parametrize( + "resolution", + [ + (10.3), + (10.), + (9.9) + ], +) +def test_grid_decimation_min(resolution): + run_filter("min", resolution)