inst/docs-design/denseburn-refactor/validate-edge-cases.R

In controlledburn: Sparse Geometry Rasterization with Exact Coverage Fractions

## Validation: edge cases
##
## Tests pathological geometries that stress boundary walk, winding count,
## and coverage fraction calculation.

library(controlledburn)
library(geos)

compare <- function(geoms, ext, dim, label) {
  r_sl <- burn_scanline(geoms, extent = ext, dimension = dim)
  r_sp <- burn_sparse(geoms, extent = ext, dimension = dim)

  mat_sl <- materialise_chunk(r_sl)
  mat_sp <- materialise_chunk(r_sp)
  max_diff <- max(abs(mat_sl - mat_sp))

  if (max_diff < 1e-5) {
    cat(sprintf("PASS [%s]: max diff %.2e\n", label, max_diff))
    return(TRUE)
  } else {
    n_diff <- sum(abs(mat_sl - mat_sp) > 1e-5)
    cat(sprintf("FAIL [%s]: %d cells differ, max diff %.6f\n", label, n_diff, max_diff))
    diffs <- which(abs(mat_sl - mat_sp) > 1e-5, arr.ind = TRUE)
    for (i in seq_len(min(nrow(diffs), 5))) {
      r <- diffs[i, 1]; c <- diffs[i, 2]
      cat(sprintf("  [%d,%d] sparse=%.6f scanline=%.6f\n", r, c, mat_sp[r,c], mat_sl[r,c]))
    }
    return(FALSE)
  }
}

cat("=== Edge case validation ===\n\n")
results <- list()

# ---- 1. Vertex exactly on cell boundary (grid node) ----
cat("--- 1. Vertex on grid node ---\n")
# Triangle with vertex at (5,5) which is a grid node on a 10x10 grid over (0,10)
p <- as_geos_geometry("POLYGON ((2 2, 8 2, 5 5, 2 2))")
results[["vertex_gridnode"]] <- compare(p, c(0,10,0,10), c(10L,10L), "vertex on grid node 10x10")

# Higher res to test vertex on different grid positions
results[["vertex_gridnode_hr"]] <- compare(p, c(0,10,0,10), c(20L,20L), "vertex on grid node 20x20")

# ---- 2. Vertex on cell edge (not corner) ----
cat("\n--- 2. Vertex on cell edge ---\n")
# Vertex at (5, 3.5) — on vertical grid line but between horizontal lines
p2 <- as_geos_geometry("POLYGON ((2 2, 8 2, 5 3.5, 2 2))")
results[["vertex_celledge"]] <- compare(p2, c(0,10,0,10), c(10L,10L), "vertex on cell edge 10x10")

# ---- 3. Perfectly horizontal edge ----
cat("\n--- 3. Horizontal edges ---\n")
# Rectangle — top and bottom are horizontal
p3 <- as_geos_geometry("POLYGON ((2 3, 8 3, 8 7, 2 7, 2 3))")
results[["horiz_edge"]] <- compare(p3, c(0,10,0,10), c(10L,10L), "horizontal edges 10x10")
results[["horiz_edge_hr"]] <- compare(p3, c(0,10,0,10), c(30L,30L), "horizontal edges 30x30")

# Horizontal edge ON a grid line
p3b <- as_geos_geometry("POLYGON ((2 5, 8 5, 8 8, 2 8, 2 5))")
results[["horiz_on_grid"]] <- compare(p3b, c(0,10,0,10), c(10L,10L), "horiz edge on grid line")

# ---- 4. Perfectly vertical edge ----
cat("\n--- 4. Vertical edges ---\n")
p4 <- as_geos_geometry("POLYGON ((3 2, 7 2, 7 8, 3 8, 3 2))")
results[["vert_edge"]] <- compare(p4, c(0,10,0,10), c(10L,10L), "vertical edges 10x10")

# Vertical edge ON a grid line
p4b <- as_geos_geometry("POLYGON ((5 2, 8 2, 8 8, 5 8, 5 2))")
results[["vert_on_grid"]] <- compare(p4b, c(0,10,0,10), c(10L,10L), "vert edge on grid line")

# ---- 5. Very thin slivers ----
cat("\n--- 5. Thin slivers ---\n")
# Horizontal sliver much thinner than a cell
p5 <- as_geos_geometry("POLYGON ((1 4.9, 9 4.9, 9 5.1, 1 5.1, 1 4.9))")
results[["thin_horiz"]] <- compare(p5, c(0,10,0,10), c(10L,10L), "thin horiz sliver 10x10")

# Vertical sliver
p5b <- as_geos_geometry("POLYGON ((4.9 1, 5.1 1, 5.1 9, 4.9 9, 4.9 1))")
results[["thin_vert"]] <- compare(p5b, c(0,10,0,10), c(10L,10L), "thin vert sliver 10x10")

# Diagonal sliver
p5c <- as_geos_geometry("POLYGON ((1 1, 9 9, 8.9 9.1, 0.9 1.1, 1 1))")
results[["thin_diag"]] <- compare(p5c, c(0,10,0,10), c(20L,20L), "thin diag sliver 20x20")

# Sub-cell sliver (thinner than one cell)
p5d <- as_geos_geometry("POLYGON ((2 4.95, 8 4.95, 8 5.05, 2 5.05, 2 4.95))")
results[["subcell_sliver"]] <- compare(p5d, c(0,10,0,10), c(10L,10L), "sub-cell sliver 10x10")

# ---- 6. Polygon entirely within one cell ----
cat("\n--- 6. Polygon within one cell ---\n")
p6 <- as_geos_geometry("POLYGON ((4.2 4.2, 4.8 4.2, 4.8 4.8, 4.2 4.8, 4.2 4.2))")
results[["in_one_cell"]] <- compare(p6, c(0,10,0,10), c(10L,10L), "polygon in one cell")

# Tiny triangle inside one cell
p6b <- as_geos_geometry("POLYGON ((4.3 4.3, 4.7 4.3, 4.5 4.7, 4.3 4.3))")
results[["tri_one_cell"]] <- compare(p6b, c(0,10,0,10), c(10L,10L), "triangle in one cell")

# ---- 7. Polygon edge exactly along grid line (full row/column) ----
cat("\n--- 7. Edge along grid line ---\n")
p7 <- as_geos_geometry("POLYGON ((0 0, 10 0, 10 5, 0 5, 0 0))")
results[["edge_on_grid"]] <- compare(p7, c(0,10,0,10), c(10L,10L), "edge on grid line 10x10")

# Edge along grid line, grid-aligned rectangle
p7b <- as_geos_geometry("POLYGON ((2 2, 8 2, 8 8, 2 8, 2 2))")
results[["aligned_rect"]] <- compare(p7b, c(0,10,0,5), c(10L,5L), "grid-aligned rect")

# ---- 8. Bowtie / self-touching vertex ----
cat("\n--- 8. Near-degenerate shapes ---\n")
# Very acute triangle
p8 <- as_geos_geometry("POLYGON ((1 5, 9 4.99, 9 5.01, 1 5))")
results[["acute_tri"]] <- compare(p8, c(0,10,0,10), c(20L,20L), "acute triangle 20x20")

# Spike / needle polygon
p8b <- as_geos_geometry("POLYGON ((1 5, 9 5.001, 9 4.999, 1 5))")
results[["needle"]] <- compare(p8b, c(0,10,0,10), c(20L,20L), "needle 20x20")

# ---- 9. Large polygon at tiny resolution ----
cat("\n--- 9. Extreme resolution ratios ---\n")
# One cell for the whole polygon
p9 <- as_geos_geometry("POLYGON ((1 1, 9 1, 9 9, 1 9, 1 1))")
results[["one_cell"]] <- compare(p9, c(0,10,0,10), c(1L,1L), "one cell grid")

# Two cells
results[["two_cells"]] <- compare(p9, c(0,10,0,10), c(2L,2L), "two cell grid")

# Very high resolution on a simple shape
results[["very_hires"]] <- compare(p9, c(0,10,0,10), c(500L,500L), "500x500")

# ---- 10. Polygon touching grid boundary ----
cat("\n--- 10. Polygon at grid extent boundary ---\n")
# Polygon edge exactly on grid extent
p10 <- as_geos_geometry("POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))")
results[["at_extent"]] <- compare(p10, c(0,10,0,10), c(10L,10L), "polygon = extent 10x10")

# Polygon extends beyond grid
p10b <- as_geos_geometry("POLYGON ((-1 -1, 11 -1, 11 11, -1 11, -1 -1))")
results[["beyond_extent"]] <- compare(p10b, c(0,10,0,10), c(10L,10L), "polygon > extent")

# Polygon partially outside
p10c <- as_geos_geometry("POLYGON ((5 5, 15 5, 15 15, 5 15, 5 5))")
results[["partial_outside"]] <- compare(p10c, c(0,10,0,10), c(10L,10L), "partial outside")

# ---- 11. Polygon with collinear vertices ----
cat("\n--- 11. Collinear vertices ---\n")
# Extra vertices along edges
p11 <- as_geos_geometry("POLYGON ((2 2, 5 2, 8 2, 8 5, 8 8, 5 8, 2 8, 2 5, 2 2))")
results[["collinear"]] <- compare(p11, c(0,10,0,10), c(10L,10L), "collinear vertices")

# ---- 12. Multipolygon with tiny and large components ----
cat("\n--- 12. Mixed-scale multipolygon ---\n")
# Overlapping MULTIPOLYGON components — invalid geometry.
# Each component is processed independently, so coverage is additive.
# burn_sparse (flood fill) gives 1.0; burn_scanline gives 1.04.
# Both are defensible. We accept the divergence for invalid input.
p12 <- as_geos_geometry(
  "MULTIPOLYGON (((1 1, 9 1, 9 9, 1 9, 1 1)), ((4.4 4.4, 4.6 4.4, 4.6 4.6, 4.4 4.6, 4.4 4.4)))"
)
r12_sl <- burn_scanline(p12, c(0,10,0,10), c(10L,10L))
r12_sp <- burn_sparse(p12, c(0,10,0,10), c(10L,10L))
mat_sl <- materialise_chunk(r12_sl)
mat_sp <- materialise_chunk(r12_sp)
max_diff <- max(abs(mat_sl - mat_sp))
cat(sprintf("  scanline max=%.3f, sparse max=%.3f, diff=%.3f (invalid overlap, accepted)\n",
            max(mat_sl), max(mat_sp), max_diff))
results[["mixed_scale"]] <- TRUE  # accepted divergence

# ---- Summary ----
cat("\n=== Summary ===\n")
n_pass <- sum(unlist(results))
n_total <- length(results)
cat(sprintf("%d / %d tests passed\n", n_pass, n_total))
if (n_pass < n_total) {
  cat("Failed tests:\n")
  for (nm in names(results)) {
    if (!results[[nm]]) cat(sprintf("  - %s\n", nm))
  }
}