Merge pull request #9 from DPIRD-DMA/perf/median-combine

perf: median-consensus variant merge + arc-length start rotation

Author: wrignj08

CHANGELOG.md

@@ -4,6 +4,13 @@ All notable changes to this project will be documented in this file.
 
 ## [Unreleased]
 
+### Changed
+- Start-point variants are now merged by a per-point median consensus instead of a geometric union. Douglas-Peucker is start-vertex dependent, so the four rotated variants disagree slightly about where a feature's vertices land; the union superimposed those disagreements, whereas the median resamples the variants to a common point count, brings them into phase via an FFT cross-correlation, and takes the coordinate-wise median — a start-invariant consensus that outvotes a single bad variant instead of averaging toward it. This is also faster than the union it replaces (~1.3x end-to-end at the default 3 iterations, ~1.4x at 5, on `examples/Water.gpkg`). The median sits at the consensus interior rather than the outer envelope, so the merged ring is marginally smaller before area preservation restores it (with `preserve_area=True`); with `preserve_area=False` the small reduction is left as-is.
+- Faster smoothing with no meaningful change to output. Start-point variants now begin at evenly spaced arc-length positions, computed by interpolating a single new start vertex, rather than rotating to a vertex index on a fully densified ring. The old path segmentized the whole ring only so that index-based rotation would land on evenly spaced points — but every vertex it added was collinear and immediately stripped by the following Douglas-Peucker. Feeding DP just the original vertices plus one cuts its input roughly 3.4x and the simplify step ~3.5x on `examples/Water.gpkg` (total single-core CPU work ~4.4s → ~2.7s, about 1.25x faster end-to-end at the default 3 iterations). The per-point median phase-alignment and Chaikin corner cutting also replace `numpy.roll` with a slice-based cyclic shift (bit-identical, fewer allocations). Output is unchanged within tolerance (per-polygon area drift ≤ 0.01%, worst concave turn unchanged).
+
+### Fixed
+- Eliminated the residual "dimple" fold that appeared where the start-point variants disagreed about a shallow feature near the simplify tolerance. Their union superimposed the disagreement as a double peak with a valley that the area-preservation buffer then sharpened into a fold slipping under the concave-turn seal threshold. The per-point median consensus (see Changed) resolves the variants to one boundary instead of overlaying them, removing the artifact (raster-noise fuzz harness: 4 → 0 residual folds; `examples/Water.gpkg`: 8 → 0).
+
 ## [0.3.2] - 2026-06-30
 
 ### Changed

README.md

@@ -196,22 +196,21 @@ smoothed = smoothify(
 
 ## How It Works
 
-Smoothify uses an advanced multi-step smoothing pipeline:
+Smoothify uses an advanced multi-step smoothing pipeline. The numbered steps below correspond to the panels in the figure:
 
 <p align="left">
   <img src="https://raw.githubusercontent.com/DPIRD-DMA/Smoothify/main/images/pipeline_steps.png" alt="Smoothify pipeline steps" width="800">
 </p>
 
 
-1. Joins touching holes (for Polygons, when `merge_holes=True`) so they smooth as one opening
-2. Adds intermediate vertices along line segments (segmentize)
-3. Generates multiple rotated variants (for Polygons) to avoid artifacts
-4. Simplifies each variant to remove noise
-5. Applies Chaikin corner cutting to smooth
-6. Merges all variants via union to eliminate start-point artifacts
-7. Applies final smoothing pass
-8. Optionally restores original area via buffering (for Polygons)
-9. Detects and repairs any sharp folds left by features near the smoothing scale (e.g. one-pixel-wide arms), using a small morphological opening/closing bounded at `segment_length / 4`
+1. **Pixelated input** — a polygon straight from raster-to-vector conversion, with a stair-stepped boundary
+2. **Multiple variants** (for Polygons) that start at evenly spaced arc-length positions, each simplified to strip staircase noise, so no artifact is tied to a fixed start vertex
+3. **Chaikin corner cutting** applied to each variant
+4. **Per-point median merge** — a start-invariant consensus that resolves the variants' disagreements
+5. **Final smoothing pass** on the merged result
+6. **Restore original area** via buffering (for Polygons, when `preserve_area=True`)
+
+Two steps are not shown in the figure: before step 2, touching holes are joined (for Polygons, when `merge_holes=True`) so they smooth as one opening; and after step 6, the result is checked for sharp concave folds left by features near the smoothing scale (e.g. one-pixel-wide arms) and repaired with a small morphological opening/closing bounded at `segment_length / 4`.
 
 ## Invalid Geometries
 

examples/Water_Smoothed.gpkg

@@ -20,7 +20,9 @@
 from smoothify.smoothify_core import (
     _CHAIKIN_SEGMENT_FACTOR,
     _chaikin_corner_cutting,
+    _combine_variants,
     _generate_starting_point_variants,
+    _polygonal_only,
     _preserve_area_with_buffer,
 )
 
@@ -51,7 +53,12 @@ def pixel_blob() -> Polygon:
     ]
     merged = unary_union(squares)
     assert isinstance(merged, Polygon)
-    return merged
+    # unary_union leaves a vertex at every pixel boundary along straight runs;
+    # drop the collinear ones (tolerance 0) so only the true staircase corners
+    # remain, matching what a polygonized raster's meaningful vertices are.
+    corners = merged.simplify(0)
+    assert isinstance(corners, Polygon)
+    return corners
 
 
 def draw(ax, geom, color="black", lw=1.8, dots=False, alpha=1.0, style="-"):
@@ -78,34 +85,35 @@ def main() -> None:
     original = pixel_blob()
 
     # --- the simplified pipeline, capturing intermediates -------------------
-    # 1. densify so the start-point rotation has vertices to rotate to
-    densified = original.segmentize(SEGMENT_LENGTH / 2)
-
-    # 2. rotated start-point variants, each simplified (noise removal) and
-    #    re-segmentized so corner rounding stays capped at segment_length
+    # 1. variants that start at evenly spaced arc-length positions, each
+    #    simplified (noise removal) and re-segmentized so corner rounding stays
+    #    capped at segment_length. No up-front densify is needed: the start
+    #    points are interpolated directly at their arc-length positions.
     variants = []
-    for variant in _generate_starting_point_variants(densified, n_starting_points=4):
+    for variant in _generate_starting_point_variants(original, n_starting_points=4):
         variant = variant.simplify(
             tolerance=SEGMENT_LENGTH, preserve_topology=True
         ).segmentize(SEGMENT_LENGTH * _CHAIKIN_SEGMENT_FACTOR)
         variants.append(variant)
 
-    # 3. Chaikin corner cutting per variant
+    # 2. Chaikin corner cutting per variant
     smoothed_variants = [
-        make_valid(_chaikin_corner_cutting(v, num_iterations=2)) for v in variants
+        _polygonal_only(make_valid(_chaikin_corner_cutting(v, num_iterations=2)))
+        for v in variants
     ]
 
-    # 4. union of the variants removes each one's start-point artifact
-    merged = make_valid(unary_union(smoothed_variants)).simplify(
+    # 3. per-point median merge: a start-invariant consensus of the variants
+    merged = _combine_variants(smoothed_variants).simplify(
         tolerance=SEGMENT_LENGTH / 5, preserve_topology=True
     )
+    assert isinstance(merged, Polygon)
 
-    # 5. final smoothing pass
+    # 4. final smoothing pass
     final_smooth = _chaikin_corner_cutting(
         merged.segmentize(SEGMENT_LENGTH * _CHAIKIN_SEGMENT_FACTOR), num_iterations=3
     )
 
-    # 6. restore the original area by buffering
+    # 5. restore the original area by buffering
     final = _preserve_area_with_buffer(
         final_smooth, target_area=original.area, tolerance=original.area * 1e-4
     )
@@ -121,32 +129,34 @@ def main() -> None:
 
     ax = panels[1]
     reference(ax, original)
-    draw(ax, densified, color="black", lw=1.0, dots=True)
-    ax.set_title("2. Densify\n(segmentize at segment_length / 2)")
+    for v, c, s in zip(variants, VARIANT_COLORS, VARIANT_STYLES, strict=True):
+        draw(ax, v, color=c, lw=1.6, alpha=0.9, style=s)
+    ax.set_title("2. Start at 4 arc-length\npositions, simplify each")
 
     ax = panels[2]
     reference(ax, original)
-    for v, c, s in zip(variants, VARIANT_COLORS, VARIANT_STYLES, strict=True):
+    for v, c, s in zip(smoothed_variants, VARIANT_COLORS, VARIANT_STYLES, strict=True):
         draw(ax, v, color=c, lw=1.6, alpha=0.9, style=s)
-    ax.set_title("3. Rotate start point 4 ways,\nsimplify each variant")
+    ax.set_title("3. Chaikin corner cutting\nper variant")
 
     ax = panels[3]
     reference(ax, original)
     for v, c, s in zip(smoothed_variants, VARIANT_COLORS, VARIANT_STYLES, strict=True):
-        draw(ax, v, color=c, lw=1.6, alpha=0.9, style=s)
-    ax.set_title("4. Chaikin corner cutting\nper variant")
+        draw(ax, v, color=c, lw=1.0, alpha=0.35, style=s)
+    draw(ax, merged, color="black", lw=2.0)
+    ax.set_title("4. Per-point median merge\n(start-invariant consensus)")
 
     ax = panels[4]
     reference(ax, original)
-    draw(ax, merged, color="black", lw=1.8)
-    ax.set_title("5. Union of variants\n(removes start-point artifacts)")
+    draw(ax, final_smooth, color="black", lw=1.8)
+    ax.set_title("5. Final smoothing pass")
 
     ax = panels[5]
     reference(ax, original)
     gpd.GeoSeries([final]).plot(ax=ax, color="#d3eed3", alpha=0.8)
     draw(ax, final, color="#1a7a1a", lw=2.0)
     err = abs(final.area - original.area) / original.area
-    ax.set_title(f"6. Final smooth + restore area\n(area error {err:.4%})")
+    ax.set_title(f"6. Restore original area\n(area error {err:.4%})")
 
     minx, miny, maxx, maxy = original.buffer(1.4).bounds
     for ax in panels: