Fix over-rounding, 3.5x perf, merge_holes option, fold repairs, quality test suites

.gitignore

@@ -184,3 +184,6 @@ smoothify/_version.py
 
 # uv lockfile (not tracked — library, resolve fresh)
 uv.lock
+
+# Benchmarks
+benchmarks/baseline_output.gpkg

.vscode/settings.json

@@ -12,6 +12,7 @@
         "geodataframe",
         "geopandas",
         "gpkg",
+        "hausdorffs",
         "ipykernel",
         "joblib",
         "lightgreen",

CHANGELOG.md

@@ -4,6 +4,25 @@ All notable changes to this project will be documented in this file.
 
 ## [Unreleased]
 
+### Changed
+- Single-core smoothing is ~3.5x faster on typical raster-derived data (benchmarked on `examples/Water.gpkg`: 6.8s → 1.9s, including the cost of the new fold-repair check below). Output differences are sub-pixel (bounded by the algorithm's own start-point noise floor); area preservation accuracy is unchanged or slightly better. The main changes:
+  - Removed the reversed-direction smoothing variants for polygons: Chaikin corner cutting is direction-invariant on closed rings, so they duplicated the forward variants bit-for-bit and only inflated the variant union (output unchanged).
+  - Holes are now subtracted in a single `difference` call against their union instead of one `intersection` + `difference` pair per hole (output unchanged).
+  - The tiny merge/dissolve buffer now uses mitre joins, which keep corners as single vertices instead of adding ~8 arc vertices per corner (boundary differences at the millimetre scale of the buffer itself).
+  - Pre-union Chaikin smoothing of the start-point variants is capped at 2 iterations; detail beyond that was erased by the post-union simplify anyway, while doubling the vertex count entering the expensive union. The final smoothing pass still runs the full `smooth_iterations`.
+  - The area-preservation root finder brackets the root from one side using its linear estimate, caches evaluations, and uses a step tolerance derived from the area tolerance via the perimeter, roughly halving the number of buffer operations.
+  - Congruent geometries (translated copies of the same shape, common in raster-derived data) are now smoothed once and the result translated to each occurrence; this also reduces work dispatched to parallel workers.
+
+### Added
+- New `merge_holes` option (default `True`): holes that touch or nearly touch (e.g. diagonally adjacent raster cells) are joined before smoothing, so they smooth into one coherent opening instead of separate overlapping shapes leaving a fake land bridge. Pass `merge_holes=False` for the previous per-hole behaviour. Mirrors what `merge_multipolygons` does for shells.
+- `examples/merge_holes_examples.ipynb`: worked examples of `merge_holes` and its interplay with `merge_collection` (touching holes, overlapping donuts, holes split across features).
+- `benchmarks/bench_water.py`: single-core timing/profiling benchmark with baseline output comparison.
+
+### Fixed
+- Fixed sharp concave folds in smoothed output on shapes with features about one `segment_length` wide (e.g. a hole with a one-pixel-wide arm, as produced by `merge_holes` joining a small hole to a larger one). The start-point smoothing variants can disagree about such features, leaving a forked slit in their union that the area-preservation shrink sharpens into a cusp. The final result is now checked for sharp concave turns and, when one is found, recomputed from the variant union sealed with a small closing (dilate-erode at `segment_length / 4`).
+- Fixed sharp cusps where a smoothed hole crosses the independently smoothed exterior and gets clipped by the hole subtraction (previously up to a 180-degree fold at the tangential crossing points). When detected, the result is repaired with a small opening + closing (at `segment_length / 4`), which removes hair-thin material needles and seals thin slits without visibly moving the boundary.
+- Fixed shapes with long straight edges being massively over-rounded (e.g. a large square with a small `segment_length` collapsed into a circle). The simplify steps strip all collinear vertices from straight edges, and Chaikin's corner cuts scale with segment length, so corner rounding grew with edge length instead of `segment_length`. Geometries are now re-segmentized after each simplify step, capping corner rounding at roughly `segment_length` while still smoothing raster staircase artifacts into curves. Applies to all geometry types. Note: outputs are somewhat denser (more vertices) and smoothing is ~20% slower on typical raster-derived data.
+
 ## [0.2.3] - 2026-06-02
 
 ### Added

README.md

@@ -28,7 +28,7 @@ Polygons and lines derived from classified raster data (e.g., ML model predictio
 Smoothify applies an optimized implementation of Chaikin's corner-cutting algorithm along with other geometric processing to create smooth, natural-looking features while:
 
 - Preserving the general shape and area of polygons
-- Supporting all shapley geometry types
+- Supporting all shapely geometry types
 - Handling shapes with interior holes
 - Efficiently processing large datasets with multiprocessing
 
@@ -85,6 +85,7 @@ Example notebooks:
 - [Usage examples](https://github.com/DPIRD-DMA/Smoothify/blob/main/examples/usage_examples.ipynb)
 - [Smoothify vs. Shapely comparison](https://github.com/DPIRD-DMA/Smoothify/blob/main/examples/smoothify_vs_shapely_comparison.ipynb)
 - [Real-world water example](https://github.com/DPIRD-DMA/Smoothify/blob/main/examples/real_world_water_example.ipynb)
+- [Merging holes](https://github.com/DPIRD-DMA/Smoothify/blob/main/examples/merge_holes_examples.ipynb)
 
 
 ### Basic Polygon Smoothing
@@ -189,21 +190,28 @@ smoothed = smoothify(
 | `merge_collection` | bool | True | Whether to merge/dissolve adjacent geometries in collections before smoothing |
 | `merge_field` | str | None | **GeoDataFrame only**: Column name to use for dissolving geometries. Only valid when `merge_collection=True`. If None, dissolves all geometries together. If specified, dissolves geometries grouped by the column values |
 | `merge_multipolygons` | bool | True | Whether to merge adjacent polygons within MultiPolygons before smoothing |
+| `merge_holes` | bool | True | Whether to join holes that touch or nearly touch (e.g. diagonally adjacent raster cells) before smoothing, so they smooth into one coherent opening instead of separate overlapping shapes |
 | `preserve_area` | bool | True | Whether to restore original area after smoothing via buffering (applies to Polygons only) |
 | `area_tolerance` | float | 0.01 | Percentage of original area allowed as error (e.g., 0.01 = 0.01% error = 99.99% preservation). Only affects Polygons when preserve_area=True |
 
 ## How It Works
 
 Smoothify uses an advanced multi-step smoothing pipeline:
 
+<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. Adds intermediate vertices along line segments (segmentize)
-2. Generates multiple rotated variants (for Polygons) to avoid artifacts
-3. Simplifies each variant to remove noise
-4. Applies Chaikin corner cutting to smooth
-5. Merges all variants via union to eliminate start-point artifacts
-6. Applies final smoothing pass
-7. Optionally restores original area via buffering (for Polygons)
+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`
 
 ## Invalid Geometries
 
@@ -224,9 +232,10 @@ smoothed = smoothify(make_valid(polygon), segment_length=1.0)
 ## Performance Considerations
 
 - **Parallel Processing**: For large GeoDataFrames or collections, use `num_cores` = 0 to enable parallel processing
+- **Duplicate Shapes**: Geometries that are translated copies of the same shape (common in raster-derived data, e.g. single-pixel polygons) are automatically smoothed once and the result reused
 - **Smoothing Iterations**: Values of 3-5 typically provide good results. Higher values create smoother output but increase processing time and vertex count
 - **Memory Usage**: Scales with geometry complexity. The algorithm creates multiple variants during smoothing
-- **Optimal segment_length**: Should match the original raster cell size (pixel size) or be slightly larger for best results
+- **Optimal segment_length**: Anything from about half the original raster pixel size and up should produce reasonable output — larger values produce more rounded output, smaller values stay more faithful to the original geometry
 
 ## Running the Tests
 

benchmarks/bench_water.py

@@ -0,0 +1,93 @@
+"""Benchmark smoothify on examples/Water.gpkg (single core).
+
+Usage:
+    python benchmarks/bench_water.py                  # time it
+    python benchmarks/bench_water.py --profile        # cProfile, print top hotspots
+    python benchmarks/bench_water.py --save-baseline  # save output for comparison
+    python benchmarks/bench_water.py --compare        # compare output to baseline
+"""
+
+import argparse
+import cProfile
+import pstats
+import time
+from pathlib import Path
+
+import geopandas as gpd
+import numpy as np
+import shapely
+
+from smoothify import smoothify
+
+ROOT = Path(__file__).resolve().parent.parent
+DATA = ROOT / "examples" / "Water.gpkg"
+BASELINE = ROOT / "benchmarks" / "baseline_output.gpkg"
+
+
+def run(gdf: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
+    return smoothify(gdf, num_cores=1)
+
+
+def quality_report(result: gpd.GeoDataFrame, baseline: gpd.GeoDataFrame) -> None:
+    """Compare result to baseline: per-feature symmetric difference / hausdorff."""
+    if len(result) != len(baseline):
+        print(f"FEATURE COUNT DIFFERS: {len(result)} vs baseline {len(baseline)}")
+        return
+    sym_ratios = []
+    hausdorffs = []
+    for g_new, g_old in zip(result.geometry, baseline.geometry, strict=True):
+        if g_old.is_empty and g_new.is_empty:
+            continue
+        denom = g_old.area if g_old.area > 0 else 1.0
+        sym_ratios.append(g_new.symmetric_difference(g_old).area / denom)
+        hausdorffs.append(g_new.hausdorff_distance(g_old))
+    sym = np.array(sym_ratios)
+    hd = np.array(hausdorffs)
+    print(f"identical features:        {(sym == 0).sum()}/{len(sym)}")
+    print(f"sym-diff area ratio  mean: {sym.mean():.2e}  max: {sym.max():.2e}")
+    print(f"hausdorff dist (m)   mean: {hd.mean():.3f}  max: {hd.max():.3f}")
+
+
+def main() -> None:
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--profile", action="store_true")
+    parser.add_argument("--save-baseline", action="store_true")
+    parser.add_argument("--compare", action="store_true")
+    parser.add_argument("--repeat", type=int, default=1)
+    parser.add_argument("--limit", type=int, default=0, help="only first N features")
+    args = parser.parse_args()
+
+    gdf = gpd.read_file(DATA)
+    if args.limit:
+        gdf = gdf.head(args.limit).copy()
+    nverts = shapely.get_num_coordinates(gdf.geometry.values).sum()
+    print(f"{len(gdf)} features, {nverts} vertices")
+
+    if args.profile:
+        profiler = cProfile.Profile()
+        profiler.enable()
+        result = run(gdf)
+        profiler.disable()
+        stats = pstats.Stats(profiler)
+        stats.sort_stats("cumulative").print_stats(30)
+        stats.sort_stats("tottime").print_stats(30)
+    else:
+        times = []
+        result = None
+        for _ in range(args.repeat):
+            t0 = time.perf_counter()
+            result = run(gdf)
+            times.append(time.perf_counter() - t0)
+        print(f"time: best {min(times):.2f}s  " + " ".join(f"{t:.2f}" for t in times))
+
+    assert result is not None
+    if args.save_baseline:
+        result.to_file(BASELINE, driver="GPKG")
+        print(f"baseline saved to {BASELINE}")
+    elif args.compare:
+        baseline = gpd.read_file(BASELINE)
+        quality_report(result, baseline)
+
+
+if __name__ == "__main__":
+    main()

images/generate_pipeline_graphic.py

@@ -0,0 +1,168 @@
+"""Generate the README pipeline graphic: how a polygon moves through smoothify.
+
+Runs a small pixelated polygon through a simplified version of the real
+smoothing pipeline (using the library's own building blocks, minus edge-case
+handling) and renders every intermediate product as one multi-panel figure.
+
+Usage:
+    python images/generate_pipeline_graphic.py
+"""
+
+from pathlib import Path
+
+import geopandas as gpd
+import matplotlib.pyplot as plt
+import numpy as np
+from shapely import make_valid
+from shapely.geometry import Polygon
+from shapely.ops import unary_union
+
+from smoothify.smoothify_core import (
+    _CHAIKIN_SEGMENT_FACTOR,
+    _chaikin_corner_cutting,
+    _generate_starting_point_variants,
+    _preserve_area_with_buffer,
+)
+
+OUT = Path(__file__).parent / "pipeline_steps.png"
+SEGMENT_LENGTH = 1.0  # "pixel size" of the demo shape
+VARIANT_COLORS = ["#d62728", "#1f77b4", "#2ca02c", "#9467bd"]
+# distinct linestyles so variants stay individually visible even where
+# their outlines coincide
+VARIANT_STYLES = ["-", "--", "-.", ":"]
+
+
+def pixel_blob() -> Polygon:
+    """A pond-like raster blob: a smooth shape rasterized at pixel size 1.
+
+    Rasterizing a real curve gives an authentic staircase boundary at a
+    realistic feature-to-pixel ratio (~20 pixels across). The lobe placement
+    is chosen so all four start-point variants simplify to visibly different
+    polygons (rotated anchors can otherwise land where two variants come out
+    identical and overprint in the figure)."""
+    xs, ys = np.meshgrid(np.arange(26), np.arange(20))
+    cx, cy = xs + 0.5, ys + 0.5
+    lobe_a = ((cx - 8.5) / 8.0) ** 2 + ((cy - 8) / 6.0) ** 2 < 1
+    lobe_b = ((cx - 16) / 6.0) ** 2 + ((cy - 10.5) / 4.5) ** 2 < 1
+    grid = lobe_a | lobe_b
+    squares = [
+        Polygon([(x, y), (x + 1, y), (x + 1, y + 1), (x, y + 1)])
+        for y, x in zip(*np.nonzero(grid), strict=True)
+    ]
+    merged = unary_union(squares)
+    assert isinstance(merged, Polygon)
+    return merged
+
+
+def draw(ax, geom, color="black", lw=1.8, dots=False, alpha=1.0, style="-"):
+    xs, ys = geom.exterior.xy
+    ax.plot(
+        xs,
+        ys,
+        style,
+        color=color,
+        linewidth=lw,
+        alpha=alpha,
+        solid_capstyle="round",
+    )
+    if dots:
+        ax.plot(xs, ys, "o", color=color, markersize=3.5, alpha=alpha)
+
+
+def reference(ax, geom):
+    gpd.GeoSeries([geom]).plot(ax=ax, color="0.92")
+    draw(ax, geom, color="0.75", lw=1.0)
+
+
+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
+    variants = []
+    for variant in _generate_starting_point_variants(densified, 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
+    smoothed_variants = [
+        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(
+        tolerance=SEGMENT_LENGTH / 5, preserve_topology=True
+    )
+
+    # 5. 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
+    final = _preserve_area_with_buffer(
+        final_smooth, target_area=original.area, tolerance=original.area * 1e-4
+    )
+
+    # --- render --------------------------------------------------------------
+    fig, axes = plt.subplots(2, 3, figsize=(13.5, 8.2))
+    panels = axes.ravel()
+
+    ax = panels[0]
+    gpd.GeoSeries([original]).plot(ax=ax, color="#cfe3f5")
+    draw(ax, original, dots=True)
+    ax.set_title("1. Pixelated input\n(polygonized raster)")
+
+    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)")
+
+    ax = panels[2]
+    reference(ax, original)
+    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("3. Rotate start point 4 ways,\nsimplify each 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")
+
+    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)")
+
+    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%})")
+
+    minx, miny, maxx, maxy = original.buffer(1.4).bounds
+    for ax in panels:
+        ax.set_aspect("equal")
+        ax.set_xticks([])
+        ax.set_yticks([])
+        ax.set_xlim(minx, maxx)
+        ax.set_ylim(miny, maxy)
+        for spine in ax.spines.values():
+            spine.set_color("0.8")
+
+    fig.suptitle("How smoothify works", fontsize=15, y=0.99)
+    plt.tight_layout()
+    plt.savefig(OUT, dpi=130, bbox_inches="tight")
+    print(f"saved {OUT}")
+
+
+if __name__ == "__main__":
+    main()