Skip to content

Skeleton

Wolren edited this page Apr 30, 2026 · 3 revisions

Skeleton Algorithm

Skeleton uses medial-axis skeleton decomposition for candidate generation, combined with SDF-guided boundary expansion and SDF-based containment certification. Fully BCRS-free — replaces the superseded contained family.

Algorithm Flow

flowchart TD
    A[Input polygon] --> B[Fast-path?<br/>Rectangle / simple convex]
    B -->|yes| C[Return instantly]
    B -->|no| D[Rasterize polygon to grid]
    D --> E[Compute distance transform]
    E --> F[Detect ridge points]
    F --> G[Cluster into regions]
    G --> H[Filter by clearance threshold]
    H --> I[PCA on each region → skeleton seeds]
    I --> J[Augment with nearest edge angles]
    J --> K[Grid solve per seed (40×40)]
    K --> L[SDF-guided binary expansion]
    L --> M[Delta test ±2°]
    M --> N[SDF containment certification]
    N --> O[Select best certified]
Loading

Stage-by-Stage

Stage 0: Fast-Path Detection

Before any skeleton work, check if the polygon is a trivial case:

  • Rectangle (4 vertices, right angles, no holes): return the polygon itself, O(1).
  • Simple convex (≤ 8 vertices, < 0.5% concavity, no holes): grid-solve at each unique edge angle (60×60), SDF expand, certify, return best — no skeleton extraction needed.

These fast-paths skip the entire SDF expansion pipeline for ~2-3% of real-world polygons.

Stage 1: Geometry Preparation

Validate polygon geometry and extract basic properties.

  1. Validate polygon: check is_valid and is_empty
  2. If multipart: use largest component by area
  3. Extract bounds: minx, miny, maxx, maxy
  4. Compute aspect ratio: h / w (long axis / short axis)

Stage 2: Medial-Axis Skeleton Extraction

The medial axis (skeleton) of a polygon is the set of points equidistant from the boundary.

Step 1: Rasterize polygon

Create a dense grid covering the polygon bounds. Long axis gets 150 cells, short axis scales proportionally (minimum 3 cells).

If polygon height exceeds width, swap dimensions to ensure consistent processing.

Step 2: Compute distance transform

For each grid cell inside the polygon, compute the Euclidean distance to the nearest boundary point via scipy.ndimage.distance_transform_edt. This creates an interior distance field where higher values = more clearance.

Step 3: Detect ridge points

Ridge points are local maxima of the distance field — they form the skeleton branches. For each cell with positive distance, check if it exceeds all neighbours in a 3×3 neighbourhood. Mark qualifying cells as ridge points.

Step 4: Cluster into regions

Apply scipy.ndimage.label for connected components labelling. Each cluster is one skeleton branch or chamber. Discard regions with < 5 pixels as noise.

Step 5: Filter by clearance

Keep only skeleton regions above the 60th percentile of clearance values. This prunes narrow corridors and noise spurs.

Step 6: Principal Component Analysis

For each surviving region, compute PCA on the real-world pixel coordinates:

  • Centre the coordinates
  • Build the 2×2 covariance matrix
  • Extract eigenvector for the largest eigenvalue → principal direction
  • Convert to angle in [0°, 90°)

Stage 3: Hybrid Seed Generation

Build candidate orientations from skeleton regions and augment with edge-angle fillers.

Step 1: Build skeleton seeds

For each skeleton region (up to top_k), store: centre, angle, maximum clearance, and score = clearance² × total clearance. Sort by score descending. Apply spatial diversity — discard seeds within 0.5× clearance of a higher-ranked seed.

Step 2: Augment with edge angles

Extract all unique edge orientations from the polygon exterior (atan2 of edge vectors, modulo 90°). Include 0° and 45° as fillers.

For each skeleton seed, find the nearest two edge angles within ±10° and add them. This fills the precision gap — skeleton PCA gives structural orientation, edge angles give exact boundary alignment.

Step 3: Fallback

If fewer than 2 skeleton seeds survive, fall through to pure edge-angle candidates (same as bcrs_worker's heuristic).

Stage 4: Coarse Grid Solve

For each candidate angle θ:

  1. Rotate polygon by −θ around polygon centroid
  2. Run uniform grid solve at 40×40 resolution:
    • Build binary occupancy mask
    • Apply stack-based largest-rectangle-in-histogram (LRH)
    • Return the best axis-aligned rectangle
  3. SDF-guided binary expansion pushes the rectangle to the boundary
  4. ±2° delta test catches angular optima missed by the seed angle

Stage 5: SDF-Guided Boundary Expansion

For each side (left, right, bottom, top):

  1. Query SDF at the side's midpoint → distance to nearest boundary
  2. Upper bound = minimum of (distance to bbox edge, abs(SDF value))
  3. Binary search over expansion distance (10 steps) to find maximum safe expansion
  4. Repeat for all four sides for 3 coordinate-ascent iterations

Uses prep.covers(box(...)) for the inside test — the single most efficient GEOS operation for whole-box containment.

Stage 6: SDF Containment Certification

Verify containment via SDF at 4 corners + 4 edge midpoints (8 sample points):

  • If max SDF ≤ 1e-7: certified contained
  • Otherwise: symmetric shrink by (max_sdf + 1e-7)
  • If shrink > 20% of shorter side, or shrunk rectangle fails: reject

Stage 7: Selection

Among all seeds, select the highest-area certified rectangle. Record: rank (0 = best), ratio, and diagnostics.


Algorithm Parameters

Parameter Default Description
TOP_K 3 Maximum candidates (3 = accuracy, 1 = speed)
MAX_RATIO 0 Max aspect ratio (0 = unlimited)
ALWAYS_RETURN True Enable shrink fallback
ANGLE_STEP 5 Fallback sweep step (degrees)

GRID_COARSE and GRID_FINE are accepted but ignored — skeleton uses its own fixed 40×40 grid.


Performance (QGIS, N_WORKERS=1)

Variant Time @290 Fill rate (mean)
skeleton (top_k=3) 24.64s 54.96%
skeleton (top_k=1) 8.41s ~53.5%

Dropping from top_k=3 to top_k=1 gives ~3× speed with ~1.5% fill-rate loss.


Output Fields

Field Description
area Rectangle area in CRS units
angle_deg Rotation angle (0-90°)
ratio Aspect ratio (long:short)
cand_rank Seed rank (0 = best)
s2_gain Area gain over baseline
best_effort 1 if shrink fallback used

Comparison with BCRS

Aspect Skeleton BCRS
Seed generation Skeleton PCA + edge angles Edge histogram
Grid solve 40×40 uniform 40×40 + 120×120 + vertex-coordinate
Angle refinement ±2° delta test Brent bounded optimisation
Expansion SDF-guided BCRS + SDF
Certification SDF SDF
Fill rate (mean) 54.96% 55.74%
Time @290 (top_k=3) 24.64s 26.33s
Time @290 (top_k=1) 8.41s 10.82s

Skeleton is BCRS‑free — no vertex-coordinate raster, no Brent polish, no grid histogram. It trades ~1% fill rate for a different seed generation path that avoids the BCRS overhead entirely. At top_k=1, the skeleton is the fastest certified‑contained solver in the pack.


References

[1] Blum, H. (1967). A transformation for extracting new descriptors of shape. In W. Wathon (Ed.), Models for the Perception of Speech and Visual Form (pp. 362-380). MIT Press.

[2] Osher, S., & Fedkiw, R. (2003). Level Set Methods and Dynamic Implicit Surfaces. Springer.


Navigation

Clone this wiki locally