diff --git a/src/filter_grid_decimation/GridDecimationFilter.cpp b/src/filter_grid_decimation/GridDecimationFilter.cpp index 1366dcb..e82acd5 100755 --- a/src/filter_grid_decimation/GridDecimationFilter.cpp +++ b/src/filter_grid_decimation/GridDecimationFilter.cpp @@ -79,15 +79,19 @@ 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); + + // 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(d_width_pt); int height = static_cast(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::max)()) throwError("Grid width out of range."); @@ -123,7 +127,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 +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); 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..f2f76a3 100755 --- a/test/test_grid_decimation.py +++ b/test/test_grid_decimation.py @@ -1,6 +1,7 @@ import csv import json import math +import re import tempfile from test import utils @@ -8,6 +9,65 @@ 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")