Proximity and Terrain Toolkit

Proximity and Terrain Toolkit

A field team has collected elevation samples across a flood-prone valley. Each record contains a horizontal position and a measured height. Analysts need quick answers to three kinds of question: What is the estimated elevation here? Which sensor is nearest? Where are the gaps, crowded clusters, and unreliable extrapolations? Your capstone is to build one coherent toolkit that answers all three from the same geometric model.

The project succeeds when every output becomes a useful input, cross-check, or index for another stage. The planar convex hull defines the trusted interpolation domain. The closest pair detects suspicious crowding. Delaunay triangles supply interpolation; their Voronoi dual supplies nearest-site regions and coverage diagnostics. A quadtree limits interactive-query candidates. Part 19's lifting viewpoint explains why the structures agree.

Problem brief and deliverables

Build a Python library plus a small command-line driver. It ingests scattered 2D samples (s_i=(x_i,y_i,z_i)), builds an immutable terrain snapshot, and supports nearest-sample, interpolated-height, window, and diagnostic queries. The implementation may use a trusted triangulation library for production, but it must expose the structures and invariants described here; a black-box interpolate(points) call does not satisfy the brief.

Your final submission should contain:

1. a loader and normalization report;
2. convex hull, closest-pair, Delaunay, Voronoi, and quadtree modules;
3. a query facade that never leaks internal indices without their source sample IDs;
4. deterministic fixtures and randomized property tests;
5. benchmark results for build time, memory, and query latency;
6. one visual inspection export, such as SVG or GeoJSON, showing samples, hull, triangles, Voronoi edges, and flagged queries.

Input and output contracts

Keep the file contract deliberately boring. CSV rows contain id,x,y,z; IDs are nonempty strings, coordinates are finite numbers, and all rows use one documented projected coordinate system. Longitude and latitude must be projected before Euclidean distance or triangulation. Optional metadata can travel beside a sample, but cannot influence predicates unless the specification says so.

Serialize status separately from the numeric result. Suggested height statuses are INTERPOLATED, ON_EDGE, ON_SAMPLE, OUTSIDE_HULL, and DEGENERATE. Never encode an outside-domain result as height zero, and never silently replace it with the nearest sample.

Architecture: one snapshot, several views

Use a central array of canonical samples and refer to positions by integer index. Hull vertices, triangle corners, Voronoi sites, and quadtree payloads then share identity. Build a new snapshot off-thread when data changes and swap it atomically.

The build pipeline is:

1. Parse, reject nonfinite values, project if necessary, and establish a scale-aware predicate policy.
2. Resolve duplicate (x,y) positions before any combinatorial algorithm.
3. Compute the convex hull and closest pair over canonical positions.
4. Construct Delaunay triangles with consistently counterclockwise vertices.
5. Derive finite Voronoi edges from adjacency, clipping rays and cells to a display rectangle when needed.
6. Index samples and triangle bounding boxes in separate quadtree payload layers.
7. Validate cross-structure identities, then publish the snapshot.

Phase 1: sanitize and establish scale

Exact duplicate positions with different heights are not a triangulation problem; they are a data-governance problem. Choose one explicit rule: reject the group, retain the newest trusted measurement, or aggregate heights and record the spread. Aggregation is often useful, but a large spread must raise a sensor-conflict warning. Do not merge merely close points using an arbitrary epsilon: that can change the hull and Delaunay topology.

Run closest pair after exact duplicates are resolved. Its distance (d_{min}) is a diagnostic, not a universal tolerance. Report pairs below an application threshold, and use the median nearest-neighbor distance as a more stable characteristic scale. Then compute the monotone-chain hull. Fewer than three noncollinear positions produce a DEGENERATE snapshot that can still answer nearest and window queries but cannot interpolate a surface.

Phase 2: triangulation, dual, and indices

Construct the Delaunay triangulation using the approach from Part 17. Store each triangle as three counterclockwise sample indices plus three neighbor triangle indices, using None across hull edges. The neighbor table enables a walking point-location query: start from the previous query's triangle, test edge orientations, and cross an edge on whose exterior side the query lies. For cold queries, ask the quadtree for triangles whose axis-aligned bounding boxes contain the point, then perform exact containment tests.

Derive Voronoi vertices as triangle circumcenters. Adjacent triangles contribute a finite dual edge between circumcenters; a hull edge contributes an outward ray. For visualization and coverage scoring, clip unbounded cells to the project rectangle. Keep this clipping boundary conceptually separate from the interpolation hull: a clipped Voronoi region outside the sample hull does not authorize terrain interpolation there.

Use two quadtree collections. The sample index accelerates nearest and window queries. The triangle index stores a triangle in every leaf intersected by its bounding box, or at the smallest node that fully contains the box. The latter avoids uncontrolled duplication, but queries must examine ancestor payloads while descending. Record the choice because it changes traversal.

Phase 3: height interpolation

For a query (q=(x,y)) inside triangle (abc), find barycentric weights satisfying

\[ q=\lambda_a a+\lambda_b b+\lambda_c c, \qquad \lambda_a+\lambda_b+\lambda_c=1. \]

Oriented areas make the computation direct. Let (D=\operatorname{orient}(a,b,c)), positive because triangles are stored counterclockwise:

\[ \lambda_a=\frac{\operatorname{orient}(b,c,q)}{D},\quad \lambda_b=\frac{\operatorname{orient}(c,a,q)}{D},\quad \lambda_c=\frac{\operatorname{orient}(a,b,q)}{D}. \]

Then (z(q)=\lambda_a z_a+\lambda_b z_b+\lambda_c z_c). Inside points have nonnegative weights; an edge point has one zero weight. This piecewise-linear terrain is continuous across shared edges because both neighboring triangles reduce to the same interpolation along that edge, although slope generally changes there.

from dataclasses import dataclass

@dataclass(frozen=True)
class HeightResult:
    status: str
    height: float | None
    triangle: tuple[str, str, str] | None
    weights: tuple[float, float, float] | None

def orient(a, b, p):
    return (b[0] - a[0]) * (p[1] - a[1]) - \
           (b[1] - a[1]) * (p[0] - a[0])

def interpolate(q, samples, ids, tol=1e-12):
    """samples are three (x, y, z) vertices in CCW order."""
    a, b, c = samples
    den = orient(a, b, c)
    if den <= 0:
        return HeightResult("DEGENERATE", None, None, None)
    w = (orient(b, c, q) / den,
         orient(c, a, q) / den,
         orient(a, b, q) / den)
    if min(w) < -tol:
        return HeightResult("OUTSIDE_TRIANGLE", None, None, w)
    # Clamp only after the robust containment decision.
    w = tuple(0.0 if abs(t) <= tol else t for t in w)
    total = sum(w)
    w = tuple(t / total for t in w)
    z = sum(t * p[2] for t, p in zip(w, samples))
    status = "ON_EDGE" if 0.0 in w else "INTERPOLATED"
    return HeightResult(status, z, tuple(ids), w)

The example tolerance is pedagogical. In the real toolkit, make triangle containment decisions with the same robust orientation predicate used to build topology. Floating-point weights may then be clamped only to remove harmless signed-zero-scale noise, never to turn a truly exterior point into an interior one.

The default outside-hull policy is refusal: return OUTSIDE_HULL, no height, and distance to the hull boundary. Offer nearest-sample or nearest-edge extrapolation only as an opt-in strategy with status EXTRAPOLATED and a maximum distance. This prevents an unbounded plane from being inferred from one skinny boundary triangle. On a hull edge, interpolate along that edge and return ON_EDGE; at a sample, return its recorded height exactly.

Precision and non-uniqueness policies

Robustness is a system property. Reuse one orientation implementation for hull turns, triangle orientation, walking, and barycentric signs, and one in-circle implementation for Delaunay decisions. A filtered predicate—fast floating point followed by exact or extended precision when the sign is uncertain—preserves performance without allowing modules to disagree.

Four cases deserve explicit tests:

Cocircular points may produce different valid triangulations, so tests should not demand one diagonal unless your symbolic tie-break specifies it. Voronoi vertices can have degree greater than three in the mathematical diagram. A triangulation that chooses a diagonal still yields the same interpolated values at vertices, but not necessarily the same interior surface when heights are arbitrary; reproducibility therefore requires a stable ordering or simulation-of-simplicity rule.

Part 19's 3D hull and lifting view supplies the conceptual audit: lift ((x,y)) to ((x,y,x^2+y^2)), and the lower convex hull projects to Delaunay triangles. If an in-circle result and the lifted orientation disagree, the implementation has mixed sign conventions or precision models.

Validation and property tests

Start with tiny fixtures whose geometry can be checked by hand: one triangle, a square with a chosen cocircular tie-break, a square plus center, collinear samples, repeated positions, and queries on every vertex and edge. Then add seeded random tests.

High-value properties include:

• the hull is counterclockwise and every sample lies on or inside every hull half-plane;
• closest-pair distance matches an (O(n^2)) oracle for small inputs;
• each Delaunay interior edge has two incident triangles and each boundary edge has one;
• no sample lies strictly inside a triangle's circumcircle, according to the chosen tie policy;
• triangle adjacency and Voronoi dual edges are reciprocal;
• barycentric weights sum to one and reconstruct both (x,y) and interpolated (z);
• heights at samples are exact, and two triangles agree along a shared edge;
• quadtree nearest and window answers match brute force;
• permutation of input rows preserves results, except where a documented stable ID tie-break applies.

Use metamorphic tests too: translation must not change distances or topology; uniform scaling multiplies distances and areas by the expected factors; adding a constant to every height adds it to every interpolated height without changing planar structures.

Benchmark the whole query path

Measure uniform, clustered, grid-like, and ridge-shaped point sets at (10^3), (10^4), and (10^5) samples. Report median and 95th-percentile build time, snapshot size, nearest-query latency, height-query latency, quadtree nodes visited, and triangle candidates tested. Warm and cold point-location workloads should be separate: coherent cursor movement favors triangle walking, while random queries favor spatial seeding.

Compare against brute force for nearest queries and a linear triangle scan for height queries. Also vary quadtree bucket capacities such as 4, 16, and 64. A benchmark result without the data distribution, predicate mode, hardware, and random seed is not reproducible. Complexity provides expectations—(O(n\log n)) construction for hull, closest pair, and Delaunay, with output-linear Voronoi extraction—but measured constants determine whether the index pays for itself.

Suggested implementation milestones

Extensions after the baseline

Useful next steps include contour extraction, slope and aspect from each triangle's fitted plane, natural-neighbor interpolation, constrained Delaunay edges for ridges and rivers, and sensor-error propagation. Streaming data invites snapshot rebuilds or dynamic quadtrees; very large domains invite tiled builds with reconciled boundaries. Three-dimensional tetrahedralization and volumetric interpolation demand stronger predicates and much more memory.

Before starting, revisit Part 1: predicates and robustness, Part 2: intersections and distances, Part 3: planar convex hulls, Part 4: closest pair, Part 14: quadtrees, Part 16: Voronoi diagrams, Part 17: Delaunay triangulations, and Part 19: 3D hulls and half-spaces. Together they turn this capstone from a collection of formulas into an auditable geometry system.