# TODO: oxigdal-algorithms
> **Purpose:** High-performance SIMD-optimized raster and vector algorithms for OxiGDAL — Pure Rust geospatial processing.
> **Status (2026-05-16):** 63,448 Rust LoC · 1,478 tests · 3 real stubs
> **Roadmap:** v0.1.5 → v0.2.0 → v1.0.0
## High Priority (next slice — verified gaps)
- [x] Replace planar-graph polygon split stub with real face-finding implementation (completed 2026-05-17)
- **Done:** Real DCEL-based polygon split at `src/vector/topology/split.rs`. `HalfEdge`/`Dcel` structs, `split_half_edge` (inserts intersection vertex), `splice_diagonal` (DCEL face-split), `collect_face_cycles` (left-turn walk), `signed_area_ring`, `point_in_ring`. The outer (most-negative-area) face is discarded. `SplitOptions::min_area`, `preserve_orientation`, `keep_holes` applied at the end.
- **Tests added (5):** test_split_square_by_vertical_returns_two_halves, test_split_square_by_diagonal_returns_two_triangles, test_split_line_misses_polygon_returns_original, test_split_concave_polygon_three_pieces, test_split_respects_min_area_filter.
- [x] Wire `CrsTransformer` to `oxigdal-proj` (completed 2026-05-17)
- **Done:** `CrsTransformer` in `src/vector/transform.rs` now holds `Option<oxigdal_proj::Transformer>` behind `#[cfg(feature = "crs-transform")]`. `build_proj_transformer()` tries a real proj transform and falls back gracefully for unknown EPSG codes. `transform_coordinate` delegates to `transform_via_proj()` when the feature is enabled; the hardcoded WGS84↔WebMerc fallback still operates without the feature. `Cargo.toml` gains `oxigdal-proj = { workspace = true, optional = true }` and `crs-transform = ["dep:oxigdal-proj", "std"]` feature.
- **Tests added (4):** test_crs_transformer_4326_to_3857_point, test_crs_transformer_wgs84_identity_passthrough, test_crs_transformer_polygon_preserves_vertex_count, test_crs_transformer_unknown_epsg_falls_back_gracefully. Total: 1448 tests pass.
- [x] Implement real Common Subexpression Elimination (CSE) for raster algebra DSL (completed 2026-05-17)
- **Done:** `src/raster/calculator/ast.rs` — added `Expr::CacheRef(u32)` variant; manual `impl Hash for Expr` hashing `f64` via `to_bits()`; manual `impl PartialEq + Eq for Expr`. `src/raster/calculator/optimizer.rs` — `eliminate_common_subexpressions` now implements hash-consing: `collect_subexpr_counts` (recursive post-order walk, ignores leaf nodes Number/Band/CacheRef), rewrites repeated subtrees to `CacheRef(slot_id)` (pre-order, frequency-sorted for determinism). `Optimizer::optimize` returns `(Expr, Vec<Expr>)` tuple. `src/raster/calculator/evaluator.rs` — `Evaluator` holds `cache_slots: &[Expr]`; `eval_pixel` takes `pixel_cache: &mut Vec<Option<f64>>` and lazily memoizes CacheRef slots. `ops.rs` — both `evaluate` and `evaluate_parallel` destructure the tuple, allocate fresh `pixel_cache` per pixel.
- **Tests added (5):** test_cse_eliminates_repeated_slope, test_cse_no_change_on_unique_subexprs, test_cse_preserves_semantics_via_eval_equality, test_cse_cache_ref_evaluation_caches_once, test_cse_handles_nested_repeated_subexpressions. Total: 1453 tests pass.
- [x] Add parallel raster processing via `rayon` feature (hillshade, slope, focal stats) (completed 2026-05-17)
- **Done:** `src/parallel/raster.rs` — `hillshade_parallel`, `slope_parallel`, `focal_parallel` (with `FocalOp` enum: Mean/Median/Min/Max/Sum/Range/StdDev), all gated on `#[cfg(feature = "parallel")]`. Strip-based decomposition with halo rows; rayon parallel iteration over output strips; each strip runs existing scalar kernel on extracted sub-buffer; results match scalar within FP rounding. Re-exported from `src/parallel/mod.rs`.
- **Tests added (9):** test_hillshade_parallel_matches_scalar, test_slope_parallel_matches_scalar, test_focal_parallel_mean_matches_scalar, test_focal_parallel_large_window_halo_correctness, test_parallel_with_single_thread_matches_serial, test_parallel_handles_small_raster_below_chunk_size, plus 3 additional correctness tests. Total: 1500 tests pass.
## Medium Priority (planned — design sketched)
- [x] Implement Snap Rounding for robust geometric operations with floating point (completed 2026-05-17)
- **Done:** New `src/vector/snap_rounding.rs` (~340 LoC). `SnapRoundingOptions { precision: f64 }`, `SnappedSegment { a, b }`, `SnapRoundingResult { segments }`. `snap_coordinate(p, precision)` — rounds to grid. `snap_linestring(coords, precision)` — snaps all vertices + deduplicates consecutive duplicates. `snap_round(segs, options) -> SnapRoundingResult` — iterative Hobby snap rounding: snap all vertices to grid, brute-force O(n²) pairwise `intersect_segment_segment` search, split segments at snapped intersection points, repeat until no new splits (convergence). Handles collinear and T-junction cases. Re-exported from `src/vector/mod.rs` and `src/lib.rs`.
- **Tests added (10):** test_snap_round_single_line_no_crossing, test_snap_round_two_crossing_lines_produces_intersection_vertex, test_snap_round_t_junction_single_split, test_snap_round_no_intersections_returns_original_segments, test_snap_round_multiple_lines_no_infinite_loop, test_snap_round_collinear_segments_no_false_intersection, test_snap_round_intersection_vertex_lies_on_grid, test_snap_linestring_removes_consecutive_duplicates, test_snap_round_invalid_precision (inline), plus additional inline unit tests. Total: 1505 tests pass.
- [x] Add weighted Voronoi diagrams (power diagrams) — Aurenhammer 1987 (completed 2026-05-17)
- **Done:** `src/vector/voronoi.rs` extended (~322 LoC added, 570 total). `WeightedPoint { point, weight }` with `new()`, `unweighted()`. `PowerCell { site_index, polygon, is_empty }`. `PowerDiagram { cells }`. `PowerDiagramOptions { bounding_box: Option<(f64,f64,f64,f64)> }`. `power_diagram(weighted_points, options) -> Result<PowerDiagram>` — half-plane intersection O(N²·V): for each site, start from 4-sided bbox polygon, clip against radical axis of every other site via Sutherland-Hodgman; empty cell detection when polygon degenerates to 0 vertices. `weighted_bisector(xi,yi,wi,xj,yj,wj) -> (a,b,c)` — radical axis coefficients for `a·x+b·y ≤ c` keeping site i's power-distance region. Private helpers: `half_plane_clip`, `intersect_segment_with_halfplane`, `compute_power_bbox`, `bbox_to_polygon`. Re-exported from `vector/mod.rs` and `lib.rs`.
- **Tests added (7):** test_power_diagram_empty_input_returns_empty_diagram, test_power_diagram_single_site_returns_full_bbox, test_power_diagram_two_equal_weights_bisector_is_midpoint, test_power_diagram_unequal_weights_shifts_boundary, test_power_diagram_heavy_site_dominates_to_empty_cell, test_power_diagram_bounding_box_clips_cells, test_weighted_bisector_perpendicular_for_equal_weights, test_power_diagram_cells_partition_bbox. Total: 1513 tests pass.
- [x] Implement constrained Delaunay triangulation (CDT) — Sloan 1987/1993. (completed 2026-05-19)
- **Done:** Extended `src/vector/delaunay.rs` (+160 LoC, total ~470). Replaced no-op `violates_constraints` stub with real segment-triangle intersection check. New pub helpers: `segment_segment_intersect_exclusive()` (strict open-interval t,u ∈ (0,1)), `cross_sign()`, `point_in_triangle_strict()` (barycentric sign test), `triangle_has_edge()`. `constrained_delaunay_with_recovery()` — Sloan-style: deduplicates constraints, skips already-present edges, collects intersecting triangles, Lawson edge-flip loop bounded by `4 * |intersecting| + 4` (Sloan §3.3); returns `AlgorithmError::InvalidInput` on convergence failure. `constrained_delaunay()` delegates to `constrained_delaunay_with_recovery()`. `constrained_delaunay_with_recovery` re-exported from `vector/mod.rs`.
- **Tests added (14 in `tests/constrained_delaunay_test.rs`):** segment intersection cases, point-in-triangle, edge membership, CDT integration (no-op, square diagonal, 5-point recovery, empty constraints, termination bound).
- [x] Add morphological operations: opening, closing, top-hat, black-hat for raster (completed 2026-05-18)
- **Done:** New `src/raster/morphology_compound.rs` (~310 LoC) — slice-based `&[f32]` compound operators alongside the existing buffer-based `RasterBuffer` API. `BorderMode` enum (`Replicate` default, `Reflect`, `Constant(f32)`). `MorphologyOptions { kernel_size, iterations, border_mode }` with `new`/`with_iterations`/`with_border_mode` builders. Slice primitives `erode_slice`/`dilate_slice` with border-mode support. Compound operators: `opening(raster, w, h, kernel_size) = erode→dilate`, `closing = dilate→erode`, `top_hat = raster - opening`, `black_hat = closing - raster`. Option-bag variants `opening_with`/`closing_with`/`top_hat_with`/`black_hat_with` apply N iterations. Existing buffer-based `top_hat`/`black_hat` re-exported as `top_hat_buffer`/`black_hat_buffer` (grep confirmed no external consumers used the simple names) to free up the slice-based API. Re-exported from `src/raster/mod.rs` and `src/lib.rs`.
- **Tests added (14):** test_opening_removes_small_bright_features, test_opening_preserves_large_bright_features, test_closing_fills_small_gaps, test_closing_preserves_large_gaps, test_top_hat_extracts_thin_bright_features, test_black_hat_extracts_thin_dark_features, test_opening_idempotent_on_second_application, test_closing_idempotent_on_second_application, test_opening_then_closing_equals_alternating_filter, test_morphology_options_iterations_compose, plus 4 additional. Total: 1532 tests pass.
- [x] Implement line offset (parallel curves) for cartographic styling (completed 2026-05-17).
- **Done:** New `src/vector/offset.rs` (~360 LoC). `JoinStyle { Miter, Bevel, Round }`, `OffsetOptions { miter_limit: f64, join_style, simplify_tolerance: Option<f64> }`, `OffsetResult { coords, was_simplified }`. `offset_linestring(coords, distance, options) -> Result<OffsetResult>` — per-segment left-normal computation; endpoint special-casing (start/end use first/last valid normal); interior vertices use bisector/miter join (bevel fallback when miter_length/distance > miter_limit), bevel cut (two points at `p+d*n_{i-1}` and `p+d*n_i`), or round arc (8 segs/π); zero-length segments skipped; optional DP simplification. `offset_polygon_rings(rings, distance, options) -> Result<Vec<...>>` — closed-ring variant where all vertices are interior; CCW detection via shoelace negates distance so positive = outward. Re-exported from `src/vector/mod.rs` and `src/lib.rs`.
- **Tests added (8+):** test_offset_horizontal_line_left_positive, test_offset_horizontal_line_right_negative, test_offset_zero_distance_returns_input, test_offset_right_angle_miter_within_limit, test_offset_right_angle_bevel_style, test_offset_miter_limit_clamps_sharp_angle, test_offset_insufficient_vertices_errors, test_offset_polygon_ring_square. Total: 1490 tests pass.
- [x] Add Frechet distance and Hausdorff distance between geometries — Alt & Godau 1995 (completed 2026-05-17).
- **Done:** `src/vector/distance.rs` extended (+111 LoC, ~722 total). `point_dist(a,b)` and `point_to_segment_dist(p,a,b)` primitives. `discrete_frechet_distance(a,b) -> f64` — Eiter-Mannila (1994) O(n·m) DP table: `ca[i][j] = max(d(a[i],b[j]), min(ca[i-1][j], ca[i][j-1], ca[i-1][j-1]))`. `directed_hausdorff(a,b) -> f64` — max over A of min distance to B (vertex variant). `hausdorff_distance(a,b) -> f64` — `max(directed(A,B), directed(B,A))`. `directed_hausdorff_to_segs` and `hausdorff_distance_to_segments` — variant using point-to-segment distance for tighter bounds. All four functions re-exported from `src/vector/mod.rs` and `src/lib.rs`.
- **Tests added (8):** test_frechet_identical_paths_zero, test_frechet_translated_paths_constant_offset, test_frechet_different_length_paths, test_frechet_single_point_paths, test_hausdorff_symmetric_for_two_squares, test_hausdorff_distance_to_segments_le_vertex_variant, test_hausdorff_degenerate_single_point, test_frechet_vs_hausdorff_ordering. Total: 1471 tests pass.
- [x] Implement minimum bounding geometry (rotated rectangle, smallest enclosing circle) (completed 2026-05-17)
- **Done:** New `src/vector/min_bounds.rs`. `RotatedRect { center, width, height, angle_rad }` with `corners() -> [(f64,f64);4]` and `area()`. `Circle { center, radius }` with `contains()` (1e-10 epsilon). `aabb(points) -> Option<(f64,f64,f64,f64)>`. `min_area_rotated_rect(points) -> Option<RotatedRect>` — rotating calipers on convex hull: for each hull edge, project all hull vertices onto `(u,v)` frame, compute AABB extent, track minimum-area orientation, reconstruct world-space center as `c_u * u_hat + c_v * v_hat`. `smallest_enclosing_circle(points) -> Circle` — recursive Welzl with `lcg_shuffle` (Knuth MMIX LCG, deterministic seed from XOR of coordinate bits, no `rand` crate). Base cases: 0=zero, 1=point, 2=diameter, 3=circumcircle (falls back to diameter on collinear). Re-exported from `src/vector/mod.rs` and `src/lib.rs`.
- **Tests added (7):** test_min_rect_axis_aligned_square_returns_square, test_min_rect_rotated_45deg_recovers_orientation, test_min_rect_random_hull_area_no_greater_than_aabb, test_welzl_circle_three_collinear_points_is_diameter, test_welzl_circle_unit_square_radius_sqrt_half, test_welzl_circle_large_random_point_set_contains_all_points, test_welzl_empty_input_returns_zero_circle. Total: 1463 tests pass.
- [ ] Add viewshed analysis with Earth curvature and refraction correction
- **Files:** `src/raster/viewshed.rs` (file exists but basic line-of-sight only)
- **Why deferred:** Atmospheric refraction model selection requires user input; document IACA-like default.
## Low Priority / Future (speculative — one-liners only)
- [x] Implement DSL expression compiler to bytecode stack VM (completed 2026-05-19)
- **Done:** New `src/raster/calculator/bytecode.rs` (~1000 LoC). `OpCode` enum (48 variants). `CompiledProgram { ops, constants, required_bands, cache_slot_count, estimated_stack_depth }`. `compile_expr` (post-order AST traversal), `compile_program`, `estimate_stack_depth` (static simulation), `eval_bytecode` (stack machine; stack underflow → error; div-by-zero → NaN). `RasterCalculator::evaluate_bytecode` full pipeline. Re-exported from mod.rs and lib.rs.
- **Tests added (23 in `tests/bytecode_test.rs`):** compile-level tests, stack machine evaluation, NDVI parity with tree-walk evaluator, edge cases (div-by-zero, mismatched dimensions).
- [x] Add TIN interpolation (natural neighbor, IDW from triangulation) (completed 2026-05-18).
- **Done:** New `src/raster/tin_interp.rs` (~430 LoC). `TinPoint { x, y, z }` + `Tin { points, triangles }` with `triangle_count()`, `point_count()`. `build_tin(&[TinPoint]) -> Result<Tin>` — adapts to existing `delaunay_triangulation(&[Point], &DelaunayOptions) -> Result<DelaunayTriangulation>`, extracts `vertices: [usize; 3]` from each `Triangle`. `find_triangle(tin, qx, qy)` — orientation-agnostic barycentric containment with `BARY_EPS = 1e-9`, skips degenerate triangles. `interpolate_idw_tin(tin, qx, qy, power) -> Option<f64>` — Shepard IDW from enclosing triangle's 3 vertices with `COINCIDENT_EPS = 1e-12` vertex-coincidence early return; defensively clamps non-positive/non-finite power to 1.0. `interpolate_natural_neighbor(tin, qx, qy) -> Option<f64>` — barycentric interpolation (Sibson area-stealing documented as v0.2.0 enhancement). `TinInterpMethod` enum (Idw { power }, NaturalNeighbor). `rasterize_tin(tin, min_x, min_y, max_x, max_y, width, height, method) -> Vec<f32>` — row-major output, row 0 at top (max_y), NaN outside hull, all-NaN buffer on degenerate dims. Re-exported from `raster/mod.rs` and `lib.rs`.
- **Tests added (21):** 9 inline unit tests + 12 integration tests covering three-points-one-triangle, too-few-points, IDW at vertex, IDW inside triangle, IDW outside hull, IDW power=2 quadratic falloff, natural-neighbor at vertex and inside, rasterize dim correctness for both methods, method dispatch, planar input returns constant z. Total: 1553 tests pass.
- [x] Implement point cloud thinning algorithms (grid, random, Poisson disk) (completed 2026-05-18).
- **Done:** New `src/raster/point_cloud_thin.rs` (~514 LoC). `ThinPoint3 { x, y, z }` with `new(x, y, z)` ctor. `ThinningMethod` enum (Grid { cell_size }, Random { target_count, seed }, PoissonDisk { min_distance, seed }). `ThinningStats { input_count, kept_count, reduction_ratio }` with `new(input, kept)` computing `1 - kept/input`. `thin_grid(&[ThinPoint3], cell_size) -> Vec<ThinPoint3>` — voxel-bucket via `HashMap<(i64,i64,i64), usize>`, first-encountered wins per voxel. `thin_random(&[ThinPoint3], target, seed) -> Vec<ThinPoint3>` — Fisher-Yates shuffle driven by Knuth MMIX 64-bit LCG (a=6364136223846793005, c=1442695040888963407); no `rand` crate. `thin_poisson_disk(&[ThinPoint3], min_distance, seed) -> Vec<ThinPoint3>` — Bridson-style greedy acceptance with cubic spatial hash (side = min_distance), 27-neighbour scan, LCG-shuffled candidate order. `thin_with_stats(points, method) -> (Vec<ThinPoint3>, ThinningStats)` dispatcher. NaN-finite guards on inputs. Re-exported from `raster/mod.rs` and `lib.rs`.
- **Tests added (23):** 6 inline unit + 17 integration tests covering empty-returns-empty, one-per-voxel, large-cell-keeps-one, distinct-voxels, target-count-respected, deterministic-same-seed, different-seed-different-result, geq-target-keeps-all, all-kept-pairs-min-distance, zero-min-distance-keeps-all, dense-input-thins, stats-correctness. Total: 1576 tests pass.
- [x] Add map generalization operators (collapse, exaggerate, displace). (completed 2026-05-19)
- **Done:** New `src/vector/generalization.rs` (~500 LoC). **Collapse** (ICA): `CollapseOptions { min_polygon_area, min_linestring_length }` with no-op default; `should_collapse_polygon`, `collapse_polygon_to_point` (area-weighted centroid), `collapse_linestring_to_point` (arc-length midpoint). **Exaggerate** (ICA): `ExaggerateAnchor { Centroid, BoundingBoxCenter, Custom(Coordinate) }`; `ExaggerateOptions { scale_factor, anchor }`; `exaggerate_coord` (anchor + scale*(point−anchor)), `exaggerate_linestring`, `exaggerate_polygon` (preserves ring orientation). **Displace** (ICA): `DisplaceOptions { min_distance, max_iterations, damping, convergence_tol }`; `DisplaceStats { iterations, final_max_displacement, converged }`; `displace_points` (O(n²) simultaneous-delta damped repulsion); `displace_polygons_by_centroid` (translates polygons by centroid deltas). Re-exported via `vector/mod.rs` and `lib.rs`.
- **Tests added (14):** collapse small/large polygon, threshold-zero keeps all, should_collapse predicate, linestring short/long, exaggerate coord 2x/1x, polygon ring count + centroid-anchor doubles extent, displace two overlapping points pushed apart, already-separated unchanged, converges within max_iter, polygons translated by centroid delta. Total: 14 pass.
- [x] Implement network routing (Dijkstra, A* on vector networks). Already implemented; verified 2026-05-20 (`src/vector/network/shortest_path.rs`: `dijkstra_search`, `astar_search`, `dijkstra_single_source` + turn-restricted variants).
- [x] Add cost-distance analysis with anisotropic friction surfaces. Already implemented; verified 2026-05-20 (`src/raster/cost_distance.rs`: `cost_distance`, `cost_distance_anisotropic`, `least_cost_path`).
- [x] Implement geometric median and other robust location estimators (completed 2026-05-19).
- **Done:** New `src/vector/robust_location.rs` (~390 LoC). `RobustLocationOptions { max_iter: 200, tol: 1e-10, coincidence_eps: 1e-12 }` with builder methods. `geometric_median(points) -> Option<(f64,f64)>` and `geometric_median_with_options` — Weiszfeld iterative algorithm: initial estimate = spatial mean; iteration `x_{k+1} = (Σ p_i/d_i) / (Σ 1/d_i)` where `d_i = max(||p_i-x_k||, coincidence_eps)`; convergence `||step|| < tol`; Vardi-Zhang coincidence handling (perturb by eps along x-axis when coincident). `weighted_geometric_median(points, weights, opts)` — same with per-point weights, None on length mismatch. `geometric_median_3d` — 3D variant. `l1_median` — coordinate-wise sort-and-median (O(n log n)). `spatial_mean` — arithmetic centroid. Re-exported from `vector/mod.rs` and `lib.rs`.
- **Tests added (14 in `tests/robust_location_test.rs` + 6 inline):** test_geometric_median_empty/single/two_points, test_equilateral_triangle_fermat_point, test_unit_square_at_centroid, test_robust_to_single_outlier, test_weighted_geometric_median_pulls_toward_heavier/uniform, test_mismatched_lengths_returns_none, test_l1_median_odd/even_count, test_spatial_mean_matches_centroid, test_3d_octahedron_at_origin, test_max_iter_zero_returns_initial_guess. Total: 1609 tests pass.
- [x] Implement RANSAC robust line and plane fitting (completed 2026-05-20).
- **Done:** New `src/vector/ransac.rs` (~883 LoC). Complements the Weiszfeld `geometric_median` above for bimodal / outlier-heavy data. `RansacOptions { max_iterations: 1000, inlier_threshold: 1e-3, min_inlier_count: 2, seed: 0xDEADBEEF, confidence: 0.99 }`; `RansacResult<M> { model, inliers, iterations, converged }`. `RansacLineModel { a, b, c }` (`a²+b²=1`) with `from_two_points`/`distance_to`/`direction`/`point_on_line`; `RansacPlaneModel { a, b, c, d }` (`a²+b²+c²=1`) with `from_three_points`/`distance_to`/`normal`. `ransac_fit_line`/`ransac_fit_plane` run a Knuth-MMIX-LCG sample-and-score loop (no `rand`; `state*6364136223846793005 + 1442695040888963407`, index from upper bits) with adaptive `N = log(1−p)/log(1−w^k)` early stopping (NaN-guarded: ratio≥1→1, ratio≤0/non-finite→max) and least-squares refit via covariance eigendecomposition (closed-form 2-D; deflated power iteration in 3-D), keeping the refit only when it retains ≥ the sampled inlier count and falling back on degenerate (all-identical) input. Reuses `Point`/`AlgorithmError::InvalidInput`. Re-exported from `vector/mod.rs`.
- **Tests added (21):** 14 in `tests/ransac_test.rs` (two-point line, perfect / 30%-outlier / single-outlier recovery, too-few-points InvalidInput, collinear, plane analogues, options default, seed reproducibility + divergence on ambiguous input, adaptive early termination) + 7 inline. Total: 1702 tests pass.
- [x] Add streaming/chunked raster processing for datasets exceeding RAM (completed 2026-05-19).
- **Done:** New `src/raster/streaming/` submodule (3 files, ~615 LoC total). `RasterError { OutOfBounds, InvalidChunkSize, IoError, AlgorithmError }`. `RasterSource` trait (`width`, `height`, `read_window`). `InMemoryRasterSource` — `Vec<f32>`-backed impl with zero-padding for out-of-bounds halo reads. `Chunk { x, y, width, height, halo, data }` — data is `(width+2·halo)×(height+2·halo)` with halo-aware scatter-write for edge chunks. `ChunkIterator<R>` over `ChunkedRaster<R>` — row-major, computes `left_pad`/`top_pad` for correct halo placement at boundaries. `process_streaming<R,F>(raster, kernel)` — applies any closure over each chunk and stitches inner results. `extract_inner` helper strips halo from result buffer. `streaming_hillshade`, `streaming_slope`, `streaming_focal_mean` wrappers — run existing kernels on halo-augmented chunks, strip inner region. Re-exported from `raster/mod.rs` and `lib.rs`.
- **Tests added (13):** test_inmemory_source_read_window_full, test_read_window_partial, test_out_of_bounds_returns_error, test_iter_chunk_count, test_chunk_dimensions_with_halo, test_edge_chunks_zero_halo_padding, test_process_streaming_identity, test_process_streaming_constant_kernel, test_streaming_hillshade_within_tolerance, test_streaming_slope_within_tolerance, test_streaming_focal_mean_radius_2, test_total_pixels_equals_raster_size, test_non_multiple_chunk_size. Total: 1589 tests pass.
## Cross-crate dependencies
- **Blocks:** `oxigdal-qc` (TopologyChecker uses `intersect_linestrings_sweep`), `oxigdal-terrain` (consumes hillshade/slope), `oxigdal-services` (vector ops)
- **Blocked by:** `oxigdal-core` (geometry types), `oxigdal-index` (R-tree for spatial joins)
## Recently completed (verbatim)
- [x] Implement Weiler-Atherton polygon clipping for robust polygon intersection (completed 2026-04-19)
- **Goal:** Production-grade polygon clipping engine supporting Intersection, Union, Difference, XOR. Handles holes, degeneracies, self-intersecting polygons. difference.rs refactored from 1859 to 1650 lines (production code 996 lines).
- **Design:** ClipOperation enum, ClipVertex with parametric alpha + cross-links, Weiler-Atherton 5-phase algorithm (intersection detection, entry/exit labeling, polygon tracing, hole handling, degeneracy handling). Area validation fallback for edge cases.
- **Files:** CREATED clipping.rs (1621 lines), MODIFIED difference.rs/intersection.rs/union_ops.rs/mod.rs/lib.rs
- **Tests:** 19 new tests: overlapping squares, holes, containment, L-shaped concave, self-intersecting bowtie, degenerate shared edge, XOR identical/disjoint, area invariants, clip_multi, validation errors
- [x] Add geodetic area calculation using Karney's method (completed 2026-04-19)
- **Goal:** Karney's ellipsoidal geodesic area algorithm, sub-millimeter accuracy, antimeridian/polar support. New AreaMethod::KarneyGeodesic variant.
- **Design:** Geodesic struct wrapping geographiclib-rs (Pure Rust), PolygonArea with signed/unsigned modes, InverseResult with S12 area element.
- **Files:** CREATE geodesic.rs (981 lines), MODIFY area.rs/mod.rs/Cargo.toml
- **Tests:** 25 tests: equator unit square (rel err <1e-6 vs GeographicLib), ellipsoid area, high-lat, antimeridian, polar, CW/CCW, diamond, perimeter, signed winding, degenerate collinear, open ring, hole subtraction, AreaMethod integration
- [x] Implement topology-preserving simplification (shared edges between adjacent polygons) (completed 2026-04-19)
- **Goal:** Shared edges simplified once, both polygons get identical result. Junction vertices preserved.
- **Design:** Edge decomposition + canonical EdgeKey (quantized 1e6), chain building per ring, Douglas-Peucker on shared chains with shared_cache (canonical poly-pair ordering), non-shared chains simplified independently, ring reassembly with junction dedup, self-intersection validation with retry at reduced tolerance.
- **Files:** CREATE src/vector/topology_simplify.rs (1419 lines), MODIFY src/vector/mod.rs (module declared + re-exports)
- **Tests:** 24 tests passing: adjacent squares basic, 2x2 grid, non-adjacent polygons, jagged shared edge, polygon with hole, junction preservation, self-intersection prevention, bowtie rejection, tolerance=0 no-op, negative tolerance error, shared-edge consistency between reverse-winding polygons, large (100-pt) circle, options/quantize/edge-key/DP helpers
- [x] Add Overlaps and Crosses spatial predicates with full DE-9IM matrix (completed 2026-04-19)
- **Goal:** Complete DE-9IM matrix computation. Fix CrossesPredicate for Polygon (OGC compliance). All named predicates: equals, disjoint, intersects, touches, crosses, within, contains, overlaps, covers, coveredBy.
- **Design:** De9im struct ([Dimension; 9]), relate() function, pattern matching, named predicate methods. Fix: Polygon/Polygon crosses() returns false.
- **Files:** CREATE de9im.rs, MODIFY contains.rs/mod.rs
- **Tests:** 37 tests (35 unit + 1 doctest + 1 crosses-fix): disjoint/overlapping/contained/touching squares, line crossing polygon, point-polygon, equals, pattern round-trip, trait impls
- [x] Implement raster contour generation using marching squares algorithm (completed 2026-04-19)
- **Goal:** Iso-contour extraction from raster grids as LineString geometries. Multi-level, linear interpolation, saddle disambiguation, NoData, GeoTransform.
- **Design:** 4-bit case index, EDGE_TABLE lookup, saddle disambiguation via center value, greedy line assembly with quantized endpoint hash map, ContourOptions/ContourLevel/ContourResult structs.
- **Files:** CREATE raster/contour.rs, MODIFY raster/mod.rs
- **Tests:** 48 tests (47 unit + 1 doctest): flat raster, 2x2 gradient, 3x3 peak, saddle point, multi-level, NoData, GeoTransform, edge cases, line assembly, simplification, 100x100 radial stress test
- [x] Add SIMD-optimized bilinear/bicubic resampling kernels (AVX2 + NEON) (completed 2026-04-19)
- **Goal:** AVX2 (8 pixels/iter) + NEON (4 pixels/iter) intrinsics for bilinear and bicubic. 3-6x throughput target.
- **Design:** Gather-lerp-lerp bilinear with FMA, separable 4x4 bicubic with Catmull-Rom and vaddvq_f32 horizontal sum. Three sub-modules: avx2_impl, neon_impl, scalar_impl with runtime dispatch.
- **Files:** MODIFIED src/simd/resampling.rs (1402 lines)
- **Tests:** 29 tests: identity, up/downsample, odd sizes, asymmetric, gradient, monotonicity, accuracy vs scalar within 1e-4 (FMA rounding tolerance)
- [x] Implement raster polygonization (connected component labeling) (completed 2026-04-19)
- **Goal:** Full raster-to-polygon via CCL + boundary tracing + polygon assembly. Equivalent to GDAL GDALPolygonize().
- **Design:** Two-pass union-find CCL (path compression + union by rank), 4/8-connectivity, dual boundary strategies (pixel-edge rectilinear matching GDAL; Moore-Neighbor pixel-center), shoelace ring classification, hole assignment via smallest-containing-exterior PIP, optional Douglas-Peucker simplification, min-area filter, GeoTransform pixel-to-world.
- **Files:** CREATED src/raster/polygonize/{mod.rs (1310L), union_find.rs (226L), boundary.rs (1093L)}; MODIFIED src/raster/mod.rs (+2 lines for re-exports)
- **Tests:** 80 tests (12 union_find + 25 boundary + 43 polygonize): single pixel, uniform, checkerboard (4-conn=9 components, 8-conn=2), donut with hole (verifies interior ring), multi-component, NoData (NaN + custom), GeoTransform, 500x500 stress test, min_area filter, Douglas-Peucker simplify, 1xN/Nx1 grids, RasterBuffer integration.
---
*Last audited: 2026-05-16*