Improve rasterization

include/shape/boolean_operations.hpp

@@ -35,6 +35,17 @@ void compute_intersection_export_inputs(
 MultiShapeWithHoles compute_intersection(
         const std::vector<MultiShapeWithHoles>& 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<ShapeWithHoles> compute_intersection_faces(
+        const std::vector<MultiShapeWithHoles>& multi_shapes);
+
 /**
  * Compute the difference between two multi-shapes.
  */

include/shape/rasterization.hpp

@@ -8,25 +8,68 @@ namespace shape
 using RowId = int64_t;
 using ColumnId = int64_t;
 
-struct Cell
+struct CellId
 {
     ColumnId column = 0;
     RowId row = 0;
 };
 
-struct IntersectedCell
+struct Cell
+{
+    /** 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;
+};
+
+/**
+ * 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
 {
-    Cell cell;
-    bool full = true;
+    ColumnId column_offset = 0;
+    RowId row_offset = 0;
+    ColumnId number_of_columns = 0;
+    RowId number_of_rows = 0;
+
+    std::vector<Cell> cells;
+
+    Cell& at(
+            ColumnId column,
+            RowId row)
+    {
+        return cells[(column - column_offset) * number_of_rows + (row - row_offset)];
+    }
+
+    const Cell& 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. 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.
  */
-std::vector<IntersectedCell> rasterization(
+RasterizedGrid rasterization(
         const ShapeWithHoles& shape,
         LengthDbl cell_width,
         LengthDbl cell_height);
@@ -38,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<Cell>& cells,
+        const std::vector<CellId>& cells,
         LengthDbl cell_width,
         LengthDbl cell_height);
 
 /**
- * Convert a list of intersected cells into shapes with holes.
+ * 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 std::vector<IntersectedCell>& cells,
+        const RasterizedGrid& grid,
         LengthDbl cell_width,
         LengthDbl cell_height,
-        bool only_full = false);
+        CellsToShapesMode mode = CellsToShapesMode::Exact);
 
 }

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<ShapeWithHoles> shape::compute_intersection_faces(
         const std::vector<MultiShapeWithHoles>& multi_shapes)
 {
     // An empty multi-shape represents an empty region: intersecting with it
@@ -1564,12 +1564,17 @@ MultiShapeWithHoles shape::compute_intersection(
         }
     }
 
-    std::vector<ShapeWithHoles> 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<MultiShapeWithHoles>& multi_shapes)
+{
+    return compute_union(compute_intersection_faces(multi_shapes));
 }
 
 MultiShapeWithHoles shape::compute_difference(