From b73f01cad624a9164558e667a249953d5b3dd1ee Mon Sep 17 00:00:00 2001 From: Florian Fontan Date: Wed, 22 Jul 2026 07:38:02 +0200 Subject: [PATCH 1/2] Make rasterization() return a dense grid instead of a sparse cell list RasterizedGrid replaces the sparse std::vector: a dense, column-major std::vector (Empty/Full/Border) plus a column/row offset and an at(column, row) accessor. rasterization()'s sweep now writes into grid slots instead of push_back-ing sparse entries, so untouched cells default to Empty rather than simply being absent. cells_to_shapes(const RasterizedGrid&, ...) replaces the vector overload. This is prep for storing each cell's clipped shape part and coverage ratio alongside its state. --- include/shape/rasterization.hpp | 50 ++++++++++--- src/rasterization.cpp | 58 ++++++++------- test/rasterization_test.cpp | 120 +++++++++++++++++--------------- 3 files changed, 138 insertions(+), 90 deletions(-) diff --git a/include/shape/rasterization.hpp b/include/shape/rasterization.hpp index 1263541..30bed6b 100644 --- a/include/shape/rasterization.hpp +++ b/include/shape/rasterization.hpp @@ -14,19 +14,53 @@ struct Cell RowId row = 0; }; -struct IntersectedCell +enum class CellState { - Cell cell; - bool full = true; + Empty, + Full, + Border, +}; + +/** + * Dense rasterization grid. + * + * Cell (column, row) is stored at + * cells[(column - column_offset) * number_of_rows + (row - row_offset)], + * i.e. column-major, matching the order in which rasterization() computes + * it. + */ +struct RasterizedGrid +{ + ColumnId column_offset = 0; + RowId row_offset = 0; + ColumnId number_of_columns = 0; + RowId number_of_rows = 0; + + std::vector cells; + + CellState& at( + ColumnId column, + RowId row) + { + return cells[(column - column_offset) * number_of_rows + (row - row_offset)]; + } + + CellState at( + ColumnId column, + RowId row) const + { + return cells[(column - column_offset) * number_of_rows + (row - row_offset)]; + } }; /** * Shape rasterization. * - * Return all the cells fully inside or intersecting the border of the given - * shape. + * Return a dense grid covering the bounding box of the given shape, with + * each cell marked Empty, Full (entirely inside the shape) or Border + * (intersecting the shape's boundary). */ -std::vector rasterization( +RasterizedGrid rasterization( const ShapeWithHoles& shape, LengthDbl cell_width, LengthDbl cell_height); @@ -54,10 +88,10 @@ MultiShapeWithHoles cells_to_shapes( LengthDbl cell_height); /** - * Convert a list of intersected cells into shapes with holes. + * Convert a rasterized grid into shapes with holes. */ MultiShapeWithHoles cells_to_shapes( - const std::vector& cells, + const RasterizedGrid& grid, LengthDbl cell_width, LengthDbl cell_height, bool only_full = false); diff --git a/src/rasterization.cpp b/src/rasterization.cpp index d2aa3dd..1b56e01 100644 --- a/src/rasterization.cpp +++ b/src/rasterization.cpp @@ -170,13 +170,11 @@ void fill_columns_intersections( } } -std::vector shape::rasterization( +RasterizedGrid shape::rasterization( const ShapeWithHoles& shape, LengthDbl cell_width, LengthDbl cell_height) { - std::vector cells; - AxisAlignedBoundingBox aabb = shape.compute_min_max(); ColumnId column_min = find_cell(cell_width, cell_height, {aabb.x_min, aabb.y_min}).column; ColumnId column_max = find_cell(cell_width, cell_height, {aabb.x_max, aabb.y_max}).column; @@ -187,6 +185,15 @@ std::vector shape::rasterization( std::cout << "row min " << row_min << " max " << row_max << std::endl; #endif + RasterizedGrid grid; + grid.column_offset = column_min; + grid.row_offset = row_min; + grid.number_of_columns = column_max - column_min + 1; + grid.number_of_rows = row_max - row_min + 1; + grid.cells = std::vector( + grid.number_of_columns * grid.number_of_rows, + CellState::Empty); + // General case: column_inters is guaranteed non-empty for each column. std::vector> column_intersections(column_max - column_min + 1); fill_columns_intersections( @@ -217,16 +224,13 @@ std::vector shape::rasterization( if (number_of_intersections == 0) { for (ColumnId column = column_min; column <= column_max; ++column) { for (RowId row = row_min; row <= row_max; ++row) { - IntersectedCell cell; - cell.cell = {column, row}; - cell.full = false; - cells.push_back(cell); + grid.at(column, row) = CellState::Border; } } #ifdef RASTERIZATION_ENABLE_DEBUG - std::cout << "cells.size() " << cells.size() << std::endl; + std::cout << "grid.cells.size() " << grid.cells.size() << std::endl; #endif - return cells; + return grid; } for (ColumnId column = column_min; column <= column_max; ++column) { @@ -272,19 +276,11 @@ std::vector shape::rasterization( << std::endl; #endif if (is_above_inside_prev) { - for (RowId row = row_hi_prev + 1; row < row_lo; ++row) { - IntersectedCell cell; - cell.cell = {column, row}; - cell.full = true; - cells.push_back(cell); - } - } - for (RowId row = (std::max)(row_lo, row_hi_prev + 1); row <= row_hi; ++row) { - IntersectedCell cell; - cell.cell = {column, row}; - cell.full = false; - cells.push_back(cell); + for (RowId row = row_hi_prev + 1; row < row_lo; ++row) + grid.at(column, row) = CellState::Full; } + for (RowId row = (std::max)(row_lo, row_hi_prev + 1); row <= row_hi; ++row) + grid.at(column, row) = CellState::Border; row_hi_prev = row_hi; if (column_intersection.y_max > y_max_prev) { is_above_inside_prev = column_intersection.is_above_inside; @@ -295,7 +291,7 @@ std::vector shape::rasterization( } } - return cells; + return grid; } Shape shape::cell_to_shape( @@ -320,15 +316,25 @@ MultiShapeWithHoles shape::cells_to_shapes( } MultiShapeWithHoles shape::cells_to_shapes( - const std::vector& cells, + const RasterizedGrid& grid, LengthDbl cell_width, LengthDbl cell_height, bool only_full) { std::vector union_input; - for (const IntersectedCell& cell: cells) - if (cell.full || !only_full) - union_input.push_back({cell_to_shape(cell.cell, cell_width, cell_height)}); + for (ColumnId column = grid.column_offset; + column < grid.column_offset + grid.number_of_columns; + ++column) { + for (RowId row = grid.row_offset; + row < grid.row_offset + grid.number_of_rows; + ++row) { + CellState state = grid.at(column, row); + if (state == CellState::Full + || (state == CellState::Border && !only_full)) { + union_input.push_back({cell_to_shape({column, row}, cell_width, cell_height)}); + } + } + } return compute_union(union_input); } diff --git a/test/rasterization_test.cpp b/test/rasterization_test.cpp index 186f4cc..144c0b2 100644 --- a/test/rasterization_test.cpp +++ b/test/rasterization_test.cpp @@ -55,11 +55,6 @@ void PrintTo(const RasterizationTestParams& params, std::ostream* os) << " cell_height " << params.cell_height << "\n"; } -bool operator==(const Cell& a, const Cell& b) -{ - return a.column == b.column && a.row == b.row; -} - class RasterizationTest: public testing::TestWithParam { }; TEST_P(RasterizationTest, Rasterization) @@ -71,68 +66,81 @@ TEST_P(RasterizationTest, Rasterization) Writer writer; writer.add_shape_with_holes(test_params.shape).write_json("rasterization_input.json"); #endif - std::vector cells = rasterization( + RasterizedGrid grid = rasterization( test_params.shape, test_params.cell_width, test_params.cell_height); #ifdef RASTERIZATION_TEST_DEBUG - for (const IntersectedCell& cell: cells) - writer.add_shape(cell_to_shape(cell.cell, test_params.cell_width, test_params.cell_height)); + for (ColumnId column = grid.column_offset; + column < grid.column_offset + grid.number_of_columns; + ++column) { + for (RowId row = grid.row_offset; + row < grid.row_offset + grid.number_of_rows; + ++row) { + if (grid.at(column, row) != CellState::Empty) + writer.add_shape(cell_to_shape({column, row}, test_params.cell_width, test_params.cell_height)); + } + } writer.write_json("rasterization_output.json"); #endif - std::cout << "cells (" << cells.size() << ")" << std::endl; - for (const IntersectedCell& ic: cells) { - std::cout << " col=" << ic.cell.column - << " row=" << ic.cell.row - << " full=" << ic.full << std::endl; - } - - // Property 1: no duplicate cells. - std::vector cell_keys; - for (const IntersectedCell& ic: cells) - cell_keys.push_back(ic.cell); - std::sort( - cell_keys.begin(), - cell_keys.end(), - [](const Cell& a, const Cell& b) { - if (a.column != b.column) - return a.column < b.column; - return a.row < b.row; - }); - for (size_t i = 1; i < cell_keys.size(); ++i) { - EXPECT_FALSE(cell_keys[i] == cell_keys[i - 1]) - << "duplicate cell col=" << cell_keys[i].column - << " row=" << cell_keys[i].row; - } - - // Property 2: all full=true cells have their center strictly inside the shape. - for (const IntersectedCell& ic: cells) { - if (!ic.full) - continue; - Point center = { - (ic.cell.column + 0.5) * test_params.cell_width, - (ic.cell.row + 0.5) * test_params.cell_height}; - EXPECT_TRUE(test_params.shape.contains(center, true)) - << "full cell col=" << ic.cell.column - << " row=" << ic.cell.row - << " center (" << center.x << "," << center.y << ")" - << " is not strictly inside the shape"; + std::cout << "grid " << grid.number_of_columns << "x" << grid.number_of_rows + << " offset (" << grid.column_offset << "," << grid.row_offset << ")" << std::endl; + for (ColumnId column = grid.column_offset; + column < grid.column_offset + grid.number_of_columns; + ++column) { + for (RowId row = grid.row_offset; + row < grid.row_offset + grid.number_of_rows; + ++row) { + CellState state = grid.at(column, row); + if (state == CellState::Empty) + continue; + std::cout << " col=" << column + << " row=" << row + << " state=" << (state == CellState::Full? "Full": "Border") << std::endl; + } } - // Property 3: all returned cells intersect the original shape. - for (const IntersectedCell& cell: cells) { - Shape cell_shape = cell_to_shape({cell.cell}, test_params.cell_width, test_params.cell_height); - EXPECT_TRUE(intersect(test_params.shape, cell_shape, true)) - << "cell col=" << cell.cell.column - << " row=" << cell.cell.row - << " does not intersect the shape"; + for (ColumnId column = grid.column_offset; + column < grid.column_offset + grid.number_of_columns; + ++column) { + for (RowId row = grid.row_offset; + row < grid.row_offset + grid.number_of_rows; + ++row) { + CellState state = grid.at(column, row); + + // Property 1: Full cells have their center strictly inside the shape. + if (state == CellState::Full) { + Point center = { + (column + 0.5) * test_params.cell_width, + (row + 0.5) * test_params.cell_height}; + EXPECT_TRUE(test_params.shape.contains(center, true)) + << "full cell col=" << column + << " row=" << row + << " center (" << center.x << "," << center.y << ")" + << " is not strictly inside the shape"; + } + + Shape cell_shape = cell_to_shape({column, row}, test_params.cell_width, test_params.cell_height); + bool cell_intersects = intersect(test_params.shape, cell_shape, true); + + // Property 2: Full and Border cells intersect the original shape. + // Property 3: Empty cells do not intersect the original shape. + if (state != CellState::Empty) { + EXPECT_TRUE(cell_intersects) + << "cell col=" << column + << " row=" << row + << " does not intersect the shape"; + } else { + EXPECT_FALSE(cell_intersects) + << "empty cell col=" << column + << " row=" << row + << " intersects the shape"; + } + } } - std::vector all_cells; - for (const IntersectedCell& cell: cells) - all_cells.push_back(cell.cell); - ShapeWithHoles cells_union = cells_to_shapes(all_cells, test_params.cell_width, test_params.cell_height).shapes_with_holes.front(); + ShapeWithHoles cells_union = cells_to_shapes(grid, test_params.cell_width, test_params.cell_height).shapes_with_holes.front(); std::vector union_output = compute_union({cells_union, test_params.shape}).shapes_with_holes; EXPECT_EQ(union_output.size(), 1); EXPECT_TRUE(equal(union_output.front(), cells_union)); From e93553bbeb26b518b7f75f1a10ac6115acf4b45f Mon Sep 17 00:00:00 2001 From: Florian Fontan Date: Wed, 22 Jul 2026 08:15:37 +0200 Subject: [PATCH 2/2] Store each cell's exact shape part and coverage ratio in RasterizedGrid Splits the former Cell struct in two: CellId (column/row, used by find_cell/cell_to_shape/the vector cells_to_shapes overload) and a new Cell (coverage in [0.0, 1.0] plus a MultiShapeWithHoles holding the part of the shape within that cell), stored densely in RasterizedGrid. Full cells get coverage = 1.0 with no stored shape, since it is trivially the cell's own rectangle. Border cells are batched and clipped in a single pass: their rectangles are grouped into one MultiShapeWithHoles and intersected against the original shape via the new compute_intersection_faces (a boolean_operations addition that returns the individual intersection faces instead of merging them with compute_union, so each face's cell can still be recovered via find_point_strictly_inside + find_cell -- safe because cell edges are part of the arrangement, so no face straddles two cells). Each border cell's coverage is then its accumulated face area over the cell area. cells_to_shapes(const RasterizedGrid&, ...) takes a CellsToShapesMode (Outer/Exact/Inner) instead of a bool, controlling how partially covered cells are handled; Exact (the new per-cell precise shape) is the default. Strengthened the rasterization test's final check from "union doesn't grow when combined with the original shape" to an exact equality between the reconstructed union and the original shape, and added coverage-range and coverage/stored-area consistency checks. --- include/shape/boolean_operations.hpp | 11 ++ include/shape/rasterization.hpp | 58 +++++--- src/boolean_operations.cpp | 11 +- src/rasterization.cpp | 195 +++++++++++++++++---------- test/rasterization_test.cpp | 51 +++++-- 5 files changed, 220 insertions(+), 106 deletions(-) diff --git a/include/shape/boolean_operations.hpp b/include/shape/boolean_operations.hpp index e4f65ab..4a52c2d 100644 --- a/include/shape/boolean_operations.hpp +++ b/include/shape/boolean_operations.hpp @@ -35,6 +35,17 @@ void compute_intersection_export_inputs( MultiShapeWithHoles compute_intersection( const std::vector& multi_shapes); +/** + * Same as compute_intersection, but return the individual faces of the + * intersection instead of merging them into a single MultiShapeWithHoles. + * + * Faces of touching/adjacent regions are not unioned together, so each + * face's provenance (e.g. which grid cell it came from) can still be + * recovered afterwards, typically via find_point_strictly_inside(). + */ +std::vector compute_intersection_faces( + const std::vector& multi_shapes); + /** * Compute the difference between two multi-shapes. */ diff --git a/include/shape/rasterization.hpp b/include/shape/rasterization.hpp index 30bed6b..b29289c 100644 --- a/include/shape/rasterization.hpp +++ b/include/shape/rasterization.hpp @@ -8,17 +8,25 @@ namespace shape using RowId = int64_t; using ColumnId = int64_t; -struct Cell +struct CellId { ColumnId column = 0; RowId row = 0; }; -enum class CellState +struct Cell { - Empty, - Full, - Border, + /** Fraction of the cell's area covered by the shape, in [0.0, 1.0]. */ + double coverage = 0.0; + + /** + * The part of the shape that lies within this cell. + * + * Only populated for cells not entirely covered (coverage < 1.0); a + * fully covered cell's shape is simply its whole rectangle, so there is + * no need to store it (see cell_to_shape). + */ + MultiShapeWithHoles shape; }; /** @@ -36,16 +44,16 @@ struct RasterizedGrid ColumnId number_of_columns = 0; RowId number_of_rows = 0; - std::vector cells; + std::vector cells; - CellState& at( + Cell& at( ColumnId column, RowId row) { return cells[(column - column_offset) * number_of_rows + (row - row_offset)]; } - CellState at( + const Cell& at( ColumnId column, RowId row) const { @@ -56,9 +64,10 @@ struct RasterizedGrid /** * Shape rasterization. * - * Return a dense grid covering the bounding box of the given shape, with - * each cell marked Empty, Full (entirely inside the shape) or Border - * (intersecting the shape's boundary). + * Return a dense grid covering the bounding box of the given shape. Each + * cell's coverage is 0.0 (outside), 1.0 (entirely inside) or in between + * (intersecting the shape's boundary), and its shape attribute holds the + * part of the original shape that lies within that cell. */ RasterizedGrid rasterization( const ShapeWithHoles& shape, @@ -72,28 +81,45 @@ void rasterization_export_inputs( LengthDbl cell_height); /** - * Convert a cell to a shape. + * Convert a cell position to a shape. */ Shape cell_to_shape( - const Cell& cell, + const CellId& cell, LengthDbl cell_width, LengthDbl cell_height); /** - * Convert a list of cells into shapes with holes. + * Convert a list of cell positions into shapes with holes. */ MultiShapeWithHoles cells_to_shapes( - const std::vector& cells, + const std::vector& cells, LengthDbl cell_width, LengthDbl cell_height); +/** + * How cells_to_shapes(const RasterizedGrid&, ...) should handle partially + * covered (0.0 < coverage < 1.0) cells. + */ +enum class CellsToShapesMode +{ + /** Include a partial cell as its whole rectangle (over-approximation). */ + Outer, + /** Include a partial cell as its exact clipped shape. */ + Exact, + /** Exclude partial cells entirely (under-approximation). */ + Inner, +}; + /** * Convert a rasterized grid into shapes with holes. + * + * Fully covered cells are always included as their whole rectangle; mode + * controls how partially covered cells are handled. */ MultiShapeWithHoles cells_to_shapes( const RasterizedGrid& grid, LengthDbl cell_width, LengthDbl cell_height, - bool only_full = false); + CellsToShapesMode mode = CellsToShapesMode::Exact); } diff --git a/src/boolean_operations.cpp b/src/boolean_operations.cpp index c67a080..236a7a7 100644 --- a/src/boolean_operations.cpp +++ b/src/boolean_operations.cpp @@ -1544,7 +1544,7 @@ void shape::compute_intersection_export_inputs( file << std::setw(4) << json << std::endl; } -MultiShapeWithHoles shape::compute_intersection( +std::vector shape::compute_intersection_faces( const std::vector& multi_shapes) { // An empty multi-shape represents an empty region: intersecting with it @@ -1564,12 +1564,17 @@ MultiShapeWithHoles shape::compute_intersection( } } - std::vector faces = compute_boolean_operation( + return compute_boolean_operation( shapes, BooleanOperation::Intersection, 1, group_ids); - return compute_union(faces); +} + +MultiShapeWithHoles shape::compute_intersection( + const std::vector& multi_shapes) +{ + return compute_union(compute_intersection_faces(multi_shapes)); } MultiShapeWithHoles shape::compute_difference( diff --git a/src/rasterization.cpp b/src/rasterization.cpp index 1b56e01..ac13fa4 100644 --- a/src/rasterization.cpp +++ b/src/rasterization.cpp @@ -17,12 +17,12 @@ using namespace shape; -Cell find_cell( +CellId find_cell( LengthDbl cell_width, LengthDbl cell_height, const Point& point) { - Cell cell; + CellId cell; cell.column = (ColumnId)std::floor(point.x / cell_width); cell.row = (RowId)std::floor(point.y / cell_height); return cell; @@ -31,7 +31,7 @@ Cell find_cell( bool contains( LengthDbl cell_width, LengthDbl cell_height, - const Cell& cell, + const CellId& cell, const Point& point) { if (strictly_lesser(point.x, cell.column * cell_width)) @@ -48,7 +48,7 @@ bool contains( bool strictly_contains( LengthDbl cell_width, LengthDbl cell_height, - const Cell& cell, + const CellId& cell, const Point& point) { if (!strictly_greater(point.x, cell.column * cell_width)) @@ -126,7 +126,7 @@ void fill_columns_intersections( ++point_pos) { const ShapePoint& point = intersection_points[point_pos]; ShapePoint point_between = shape.find_point_between(point_prev, point); - Cell cell = find_cell(cell_width, cell_height, point_between.point); + CellId cell = find_cell(cell_width, cell_height, point_between.point); if (!equal(point_between.point.x, cell.column * cell_width) && !equal(point_between.point.x, (cell.column + 1) * cell_width)) { AxisAlignedBoundingBox aabb = shape.compute_min_max(point_prev, point); @@ -190,9 +190,20 @@ RasterizedGrid shape::rasterization( grid.row_offset = row_min; grid.number_of_columns = column_max - column_min + 1; grid.number_of_rows = row_max - row_min + 1; - grid.cells = std::vector( - grid.number_of_columns * grid.number_of_rows, - CellState::Empty); + grid.cells = std::vector(grid.number_of_columns * grid.number_of_rows); + + // Cells fully inside the shape: known immediately, no need for the + // intersection computation below; their shape is simply their whole + // rectangle, so there is no need to store it (see cell_to_shape). + auto mark_full = [&](ColumnId column, RowId row) { + grid.at(column, row).coverage = 1.0; + }; + // Cells intersecting the shape's boundary: their exact shape part and + // coverage are computed once all of them are known, below. + std::vector border_cells; + auto mark_border = [&](ColumnId column, RowId row) { + border_cells.push_back({column, row}); + }; // General case: column_inters is guaranteed non-empty for each column. std::vector> column_intersections(column_max - column_min + 1); @@ -222,80 +233,106 @@ RasterizedGrid shape::rasterization( // Single-column or single-row case: fill_columns_intersections produces no // records because the shape crosses no horizontal grid lines. if (number_of_intersections == 0) { + for (ColumnId column = column_min; column <= column_max; ++column) + for (RowId row = row_min; row <= row_max; ++row) + mark_border(column, row); + } else { for (ColumnId column = column_min; column <= column_max; ++column) { - for (RowId row = row_min; row <= row_max; ++row) { - grid.at(column, row) = CellState::Border; - } - } -#ifdef RASTERIZATION_ENABLE_DEBUG - std::cout << "grid.cells.size() " << grid.cells.size() << std::endl; -#endif - return grid; - } - - for (ColumnId column = column_min; column <= column_max; ++column) { #ifdef RASTERIZATION_ENABLE_DEBUG - std::cout << "column " << column << std::endl; + std::cout << "column " << column << std::endl; #endif - const std::vector& column_inters = - column_intersections[column - column_min]; + const std::vector& column_inters = + column_intersections[column - column_min]; - // Sort intersections by ascending min(y_in, y_out). - std::vector sorted_intersections = column_inters; - std::sort( - sorted_intersections.begin(), - sorted_intersections.end(), - []( - const ColumnIntersection& column_intersection_1, - const ColumnIntersection& column_intersection_2) - { - return column_intersection_1.y_min < column_intersection_2.y_min; - }); + // Sort intersections by ascending min(y_in, y_out). + std::vector sorted_intersections = column_inters; + std::sort( + sorted_intersections.begin(), + sorted_intersections.end(), + []( + const ColumnIntersection& column_intersection_1, + const ColumnIntersection& column_intersection_2) + { + return column_intersection_1.y_min < column_intersection_2.y_min; + }); - LengthDbl y_max_prev = -std::numeric_limits::infinity(); - RowId row_hi_prev = row_min - 1; - bool is_above_inside_prev = false; - for (ElementPos intersection_pos = 0; - intersection_pos < sorted_intersections.size(); - ++intersection_pos) { - const ColumnIntersection& column_intersection = sorted_intersections[intersection_pos]; - RowId row_lo = (RowId)std::floor(column_intersection.y_min / cell_height); - if (equal((row_lo + 1) * cell_height, column_intersection.y_min)) - row_lo++; - RowId row_hi = (RowId)std::floor(column_intersection.y_max / cell_height); - if (equal(row_hi * cell_height, column_intersection.y_max)) - row_hi--; + LengthDbl y_max_prev = -std::numeric_limits::infinity(); + RowId row_hi_prev = row_min - 1; + bool is_above_inside_prev = false; + for (ElementPos intersection_pos = 0; + intersection_pos < sorted_intersections.size(); + ++intersection_pos) { + const ColumnIntersection& column_intersection = sorted_intersections[intersection_pos]; + RowId row_lo = (RowId)std::floor(column_intersection.y_min / cell_height); + if (equal((row_lo + 1) * cell_height, column_intersection.y_min)) + row_lo++; + RowId row_hi = (RowId)std::floor(column_intersection.y_max / cell_height); + if (equal(row_hi * cell_height, column_intersection.y_max)) + row_hi--; #ifdef RASTERIZATION_ENABLE_DEBUG - std::cout << "i " << intersection_pos - << " is_above_inside_prev " << is_above_inside_prev - << " row_hi_prev " << row_hi_prev - << " y_min " << column_intersection.y_min - << " y_max " << column_intersection.y_max - << " row_lo " << row_lo - << " row_hi " << row_hi - << std::endl; + std::cout << "i " << intersection_pos + << " is_above_inside_prev " << is_above_inside_prev + << " row_hi_prev " << row_hi_prev + << " y_min " << column_intersection.y_min + << " y_max " << column_intersection.y_max + << " row_lo " << row_lo + << " row_hi " << row_hi + << std::endl; #endif - if (is_above_inside_prev) { - for (RowId row = row_hi_prev + 1; row < row_lo; ++row) - grid.at(column, row) = CellState::Full; - } - for (RowId row = (std::max)(row_lo, row_hi_prev + 1); row <= row_hi; ++row) - grid.at(column, row) = CellState::Border; - row_hi_prev = row_hi; - if (column_intersection.y_max > y_max_prev) { - is_above_inside_prev = column_intersection.is_above_inside; - y_max_prev = column_intersection.y_max; - } else { - is_above_inside_prev = is_above_inside_prev & column_intersection.is_above_inside; + if (is_above_inside_prev) { + for (RowId row = row_hi_prev + 1; row < row_lo; ++row) + mark_full(column, row); + } + for (RowId row = (std::max)(row_lo, row_hi_prev + 1); row <= row_hi; ++row) + mark_border(column, row); + row_hi_prev = row_hi; + if (column_intersection.y_max > y_max_prev) { + is_above_inside_prev = column_intersection.is_above_inside; + y_max_prev = column_intersection.y_max; + } else { + is_above_inside_prev = is_above_inside_prev & column_intersection.is_above_inside; + } } } } + // Compute the exact shape part and coverage of each border cell by + // intersecting the original shape with the union of the border cells' + // rectangles, without merging the resulting faces together: since cell + // boundaries are part of the arrangement, no face can straddle two + // cells, so each face's interior point identifies exactly one cell. + if (!border_cells.empty()) { + MultiShapeWithHoles shape_group; + shape_group.shapes_with_holes.push_back(shape); + + MultiShapeWithHoles border_cells_group; + for (const CellId& cell: border_cells) { + border_cells_group.shapes_with_holes.push_back( + {cell_to_shape(cell, cell_width, cell_height)}); + } + + std::vector faces = compute_intersection_faces( + {shape_group, border_cells_group}); + for (const ShapeWithHoles& face: faces) { + Point point = face.find_point_strictly_inside(); + CellId cell = find_cell(cell_width, cell_height, point); + grid.at(cell.column, cell.row).shape.shapes_with_holes.push_back(face); + } + + for (const CellId& cell_id: border_cells) { + Cell& cell = grid.at(cell_id.column, cell_id.row); + AreaDbl area = 0.0; + for (const ShapeWithHoles& face: cell.shape.shapes_with_holes) + area += face.compute_area(); + cell.coverage = area / (cell_width * cell_height); + } + } + return grid; } Shape shape::cell_to_shape( - const Cell& cell, + const CellId& cell, LengthDbl cell_width, LengthDbl cell_height) { @@ -305,12 +342,12 @@ Shape shape::cell_to_shape( } MultiShapeWithHoles shape::cells_to_shapes( - const std::vector& cells, + const std::vector& cells, LengthDbl cell_width, LengthDbl cell_height) { std::vector union_input; - for (const Cell& cell: cells) + for (const CellId& cell: cells) union_input.push_back({cell_to_shape(cell, cell_width, cell_height)}); return compute_union(union_input); } @@ -319,7 +356,7 @@ MultiShapeWithHoles shape::cells_to_shapes( const RasterizedGrid& grid, LengthDbl cell_width, LengthDbl cell_height, - bool only_full) + CellsToShapesMode mode) { std::vector union_input; for (ColumnId column = grid.column_offset; @@ -328,10 +365,22 @@ MultiShapeWithHoles shape::cells_to_shapes( for (RowId row = grid.row_offset; row < grid.row_offset + grid.number_of_rows; ++row) { - CellState state = grid.at(column, row); - if (state == CellState::Full - || (state == CellState::Border && !only_full)) { + const Cell& cell = grid.at(column, row); + if (equal(cell.coverage, 1.0)) { union_input.push_back({cell_to_shape({column, row}, cell_width, cell_height)}); + } else if (cell.coverage > 0.0) { + switch (mode) { + case CellsToShapesMode::Outer: { + union_input.push_back({cell_to_shape({column, row}, cell_width, cell_height)}); + break; + } case CellsToShapesMode::Exact: { + for (const ShapeWithHoles& part: cell.shape.shapes_with_holes) + union_input.push_back(part); + break; + } case CellsToShapesMode::Inner: { + break; + } + } } } } diff --git a/test/rasterization_test.cpp b/test/rasterization_test.cpp index 144c0b2..fea953c 100644 --- a/test/rasterization_test.cpp +++ b/test/rasterization_test.cpp @@ -77,7 +77,7 @@ TEST_P(RasterizationTest, Rasterization) for (RowId row = grid.row_offset; row < grid.row_offset + grid.number_of_rows; ++row) { - if (grid.at(column, row) != CellState::Empty) + if (grid.at(column, row).coverage > 0.0) writer.add_shape(cell_to_shape({column, row}, test_params.cell_width, test_params.cell_height)); } } @@ -92,25 +92,33 @@ TEST_P(RasterizationTest, Rasterization) for (RowId row = grid.row_offset; row < grid.row_offset + grid.number_of_rows; ++row) { - CellState state = grid.at(column, row); - if (state == CellState::Empty) + const Cell& cell = grid.at(column, row); + if (cell.coverage <= 0.0) continue; std::cout << " col=" << column << " row=" << row - << " state=" << (state == CellState::Full? "Full": "Border") << std::endl; + << " coverage=" << cell.coverage << std::endl; } } + LengthDbl cell_area = test_params.cell_width * test_params.cell_height; for (ColumnId column = grid.column_offset; column < grid.column_offset + grid.number_of_columns; ++column) { for (RowId row = grid.row_offset; row < grid.row_offset + grid.number_of_rows; ++row) { - CellState state = grid.at(column, row); + const Cell& cell = grid.at(column, row); - // Property 1: Full cells have their center strictly inside the shape. - if (state == CellState::Full) { + // Property 1: coverage stays within [0, 1]. + EXPECT_GE(cell.coverage, 0.0) + << "col=" << column << " row=" << row; + EXPECT_LE(cell.coverage, 1.0 + 1e-9) + << "col=" << column << " row=" << row; + + // Property 2: fully covered cells have their center strictly + // inside the shape. + if (equal(cell.coverage, 1.0)) { Point center = { (column + 0.5) * test_params.cell_width, (row + 0.5) * test_params.cell_height}; @@ -124,9 +132,9 @@ TEST_P(RasterizationTest, Rasterization) Shape cell_shape = cell_to_shape({column, row}, test_params.cell_width, test_params.cell_height); bool cell_intersects = intersect(test_params.shape, cell_shape, true); - // Property 2: Full and Border cells intersect the original shape. - // Property 3: Empty cells do not intersect the original shape. - if (state != CellState::Empty) { + // Property 3: covered cells (coverage > 0) intersect the shape; + // uncovered cells (coverage == 0) do not. + if (cell.coverage > 0.0) { EXPECT_TRUE(cell_intersects) << "cell col=" << column << " row=" << row @@ -137,13 +145,28 @@ TEST_P(RasterizationTest, Rasterization) << " row=" << row << " intersects the shape"; } + + // Property 4: a partially covered cell's stored shape area + // matches its coverage ratio. + if (cell.coverage > 0.0 && !equal(cell.coverage, 1.0)) { + AreaDbl shape_area = 0.0; + for (const ShapeWithHoles& part: cell.shape.shapes_with_holes) + shape_area += part.compute_area(); + EXPECT_TRUE(equal(shape_area / cell_area, cell.coverage)) + << "cell col=" << column + << " row=" << row + << " stored shape area / cell area (" << (shape_area / cell_area) << ")" + << " does not match coverage (" << cell.coverage << ")"; + } } } - ShapeWithHoles cells_union = cells_to_shapes(grid, test_params.cell_width, test_params.cell_height).shapes_with_holes.front(); - std::vector union_output = compute_union({cells_union, test_params.shape}).shapes_with_holes; - EXPECT_EQ(union_output.size(), 1); - EXPECT_TRUE(equal(union_output.front(), cells_union)); + // The union of all cells' exact shape parts must reconstruct the + // original shape. + MultiShapeWithHoles cells_union = cells_to_shapes( + grid, test_params.cell_width, test_params.cell_height); + ASSERT_EQ(cells_union.shapes_with_holes.size(), 1); + EXPECT_TRUE(equal(cells_union.shapes_with_holes.front(), test_params.shape)); } INSTANTIATE_TEST_SUITE_P(