feat(terrain): build BAS basement foundations

3 files changed, 2953 insertions, 24 deletions

diff --git a/crates/fparkan-terrain/src/contour.rs b/crates/fparkan-terrain/src/contour.rs
new file mode 100644
index 0000000..de5077a
--- /dev/null
+++ b/crates/fparkan-terrain/src/contour.rs
@@ -0,0 +1,809 @@
+/// One triangle fragment retained after subtracting a building contour.
+/// `source_triangle` indexes the triangle array passed to
+/// [`subtract_contour`].
+#[derive(Clone, Copy, Debug, PartialEq)]
+pub(crate) struct TerrainTriangleFragment {
+    /// Source triangle index, for preserving its material and vertex channels.
+    pub source_triangle: usize,
+    /// Fragment triangle vertices in the input coordinate space.
+    pub positions: [[f32; 3]; 3],
+}
+
+/// Error returned when a terrain contour cannot be subtracted safely.
+#[derive(Clone, Debug, Eq, PartialEq)]
+pub(crate) enum TerrainContourError {
+    /// The contour is non-finite, degenerate, or self-intersecting.
+    UnsupportedContour,
+    /// A source triangle has non-finite coordinates or a degenerate XY plane.
+    InvalidSourceTriangle {
+        /// Index of the source triangle.
+        triangle: usize,
+    },
+}
+
+impl std::fmt::Display for TerrainContourError {
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        match self {
+            Self::UnsupportedContour => write!(
+                f,
+                "terrain contour must be finite, simple, and non-degenerate"
+            ),
+            Self::InvalidSourceTriangle { triangle } => {
+                write!(
+                    f,
+                    "terrain source triangle {triangle} is not a finite XY surface"
+                )
+            }
+        }
+    }
+}
+
+impl std::error::Error for TerrainContourError {}
+
+/// Subtracts a finite simple contour using the native per-source-face point set:
+/// the three source vertices plus the actual contour segments clipped to that
+/// source triangle. It does not extend cuts across the triangle or triangulate
+/// the contour into unrelated internal edges.
+///
+/// # Errors
+///
+/// Returns [`TerrainContourError::UnsupportedContour`] for malformed or
+/// self-intersecting contours, and
+/// [`TerrainContourError::InvalidSourceTriangle`] when a source face cannot be
+/// triangulated or lifted onto its original plane.
+pub(crate) fn subtract_contour(
+    triangles: &[[[f32; 3]; 3]],
+    contour: &[[f32; 2]],
+) -> Result<Vec<TerrainTriangleFragment>, TerrainContourError> {
+    let ring = normalize_ring(contour)?;
+    let contour_min: [f64; 2] = std::array::from_fn(|axis| {
+        ring.iter()
+            .map(|point| point[axis])
+            .fold(f64::INFINITY, f64::min)
+    });
+    let contour_max: [f64; 2] = std::array::from_fn(|axis| {
+        ring.iter()
+            .map(|point| point[axis])
+            .fold(f64::NEG_INFINITY, f64::max)
+    });
+    let mut output = Vec::new();
+
+    for (source_triangle, source) in triangles.iter().copied().enumerate() {
+        if !source.iter().flatten().all(|value| value.is_finite()) {
+            return Err(TerrainContourError::InvalidSourceTriangle {
+                triangle: source_triangle,
+            });
+        }
+        let source_xy = source.map(|point| [f64::from(point[0]), f64::from(point[1])]);
+        let source_area = cross2(source_xy[0], source_xy[1], source_xy[2]);
+        if source_area.abs() <= 1.0e-10 {
+            return Err(TerrainContourError::InvalidSourceTriangle {
+                triangle: source_triangle,
+            });
+        }
+        let source_min: [f64; 2] = std::array::from_fn(|axis| {
+            source_xy
+                .iter()
+                .map(|point| point[axis])
+                .fold(f64::INFINITY, f64::min)
+        });
+        let source_max: [f64; 2] = std::array::from_fn(|axis| {
+            source_xy
+                .iter()
+                .map(|point| point[axis])
+                .fold(f64::NEG_INFINITY, f64::max)
+        });
+        if (0..2).any(|axis| {
+            source_max[axis] < contour_min[axis] || contour_max[axis] < source_min[axis]
+        }) {
+            output.push(TerrainTriangleFragment {
+                source_triangle,
+                positions: source,
+            });
+            continue;
+        }
+
+        let mut vertices = source_xy.to_vec();
+        let initial = if source_area > 0.0 {
+            [0, 1, 2]
+        } else {
+            [0, 2, 1]
+        };
+        let mut triangulation = vec![initial];
+        let mut clipped_segments = Vec::new();
+        for index in 0..ring.len() {
+            let first = ring[index];
+            let second = ring[(index + 1) % ring.len()];
+            if let Some((start, end)) = clip_segment_to_triangle(first, second, source_xy) {
+                let start_id = insert_unique_point(&mut vertices, start);
+                let end_id = insert_unique_point(&mut vertices, end);
+                if start_id != end_id {
+                    clipped_segments.push((start_id, end_id));
+                }
+            }
+        }
+
+        if clipped_segments.is_empty() {
+            let center = [
+                (source_xy[0][0] + source_xy[1][0] + source_xy[2][0]) / 3.0,
+                (source_xy[0][1] + source_xy[1][1] + source_xy[2][1]) / 3.0,
+            ];
+            if !point_in_ring(center, &ring) {
+                output.push(TerrainTriangleFragment {
+                    source_triangle,
+                    positions: source,
+                });
+            }
+            continue;
+        }
+
+        for point_id in 3..vertices.len() {
+            insert_triangulation_point(point_id, &vertices, &mut triangulation).map_err(|()| {
+                TerrainContourError::InvalidSourceTriangle {
+                    triangle: source_triangle,
+                }
+            })?;
+        }
+
+        let mut constraints = source_boundary_constraints(&vertices, source_xy);
+        for (first, second) in clipped_segments {
+            let mut on_segment = vertices
+                .iter()
+                .copied()
+                .enumerate()
+                .filter_map(|(index, point)| {
+                    point_on_segment(point, vertices[first], vertices[second]).then_some((
+                        segment_parameter(point, vertices[first], vertices[second]),
+                        index,
+                    ))
+                })
+                .collect::<Vec<_>>();
+            on_segment.sort_by(|left, right| left.0.total_cmp(&right.0));
+            on_segment.dedup_by_key(|(_, index)| *index);
+            for pair in on_segment.windows(2) {
+                let edge = canonical_edge(pair[0].1, pair[1].1);
+                if edge.0 == edge.1 {
+                    continue;
+                }
+                recover_constraint(edge, &mut triangulation, &vertices, &constraints).map_err(
+                    |()| TerrainContourError::InvalidSourceTriangle {
+                        triangle: source_triangle,
+                    },
+                )?;
+                constraints.insert(edge);
+            }
+        }
+        legalize_constrained_delaunay(&mut triangulation, &vertices, &constraints).map_err(
+            |()| TerrainContourError::InvalidSourceTriangle {
+                triangle: source_triangle,
+            },
+        )?;
+
+        for triangle in triangulation {
+            let center = triangle.map(|index| vertices[index]);
+            let center = [
+                (center[0][0] + center[1][0] + center[2][0]) / 3.0,
+                (center[0][1] + center[1][1] + center[2][1]) / 3.0,
+            ];
+            if point_in_ring(center, &ring) {
+                continue;
+            }
+            let mut positions: [[f32; 3]; 3] = triangle
+                .map(|index| lift_to_source_plane(vertices[index], source))
+                .into_iter()
+                .collect::<Option<Vec<_>>>()
+                .and_then(|points| points.try_into().ok())
+                .ok_or(TerrainContourError::InvalidSourceTriangle {
+                    triangle: source_triangle,
+                })?;
+            if source_area < 0.0 {
+                positions.swap(1, 2);
+            }
+            output.push(TerrainTriangleFragment {
+                source_triangle,
+                positions,
+            });
+        }
+    }
+    Ok(output)
+}
+
+fn normalize_ring(contour: &[[f32; 2]]) -> Result<Vec<[f64; 2]>, TerrainContourError> {
+    let mut ring = contour
+        .iter()
+        .map(|point| [f64::from(point[0]), f64::from(point[1])])
+        .collect::<Vec<_>>();
+    if ring.iter().flatten().any(|value| !value.is_finite()) {
+        return Err(TerrainContourError::UnsupportedContour);
+    }
+    if ring.len() > 1 && distance_squared(ring[0], *ring.last().expect("nonempty")) <= 1.0e-16 {
+        ring.pop();
+    }
+    let mut normalized = Vec::with_capacity(ring.len());
+    for point in ring {
+        if normalized
+            .last()
+            .is_none_or(|previous| distance_squared(*previous, point) > 1.0e-16)
+        {
+            normalized.push(point);
+        }
+    }
+    if normalized.len() > 1
+        && distance_squared(normalized[0], *normalized.last().expect("nonempty")) <= 1.0e-16
+    {
+        normalized.pop();
+    }
+    if normalized.len() < 3 {
+        return Err(TerrainContourError::UnsupportedContour);
+    }
+    let signed_area = polygon_area(&normalized);
+    if !signed_area.is_finite() || signed_area.abs() <= 1.0e-10 {
+        return Err(TerrainContourError::UnsupportedContour);
+    }
+    for first in 0..normalized.len() {
+        let first_next = (first + 1) % normalized.len();
+        for second in (first + 1)..normalized.len() {
+            let second_next = (second + 1) % normalized.len();
+            if first_next == second || second_next == first {
+                continue;
+            }
+            if segments_intersect(
+                normalized[first],
+                normalized[first_next],
+                normalized[second],
+                normalized[second_next],
+            ) {
+                return Err(TerrainContourError::UnsupportedContour);
+            }
+        }
+    }
+    if signed_area < 0.0 {
+        normalized.reverse();
+    }
+    Ok(normalized)
+}
+
+fn clip_segment_to_triangle(
+    first: [f64; 2],
+    second: [f64; 2],
+    triangle: [[f64; 2]; 3],
+) -> Option<([f64; 2], [f64; 2])> {
+    let mut triangle = triangle;
+    if cross2(triangle[0], triangle[1], triangle[2]) < 0.0 {
+        triangle.swap(1, 2);
+    }
+    let mut start = 0.0_f64;
+    let mut end = 1.0_f64;
+    for edge in 0..3 {
+        let edge_first = triangle[edge];
+        let edge_second = triangle[(edge + 1) % 3];
+        let at_start = cross2(edge_first, edge_second, first);
+        let at_end = cross2(edge_first, edge_second, second);
+        let change = at_end - at_start;
+        if change.abs() <= 1.0e-14 {
+            if at_start < -1.0e-10 {
+                return None;
+            }
+            continue;
+        }
+        let crossing = -at_start / change;
+        if change > 0.0 {
+            start = start.max(crossing);
+        } else {
+            end = end.min(crossing);
+        }
+        if start > end + 1.0e-12 {
+            return None;
+        }
+    }
+    start = start.clamp(0.0, 1.0);
+    end = end.clamp(0.0, 1.0);
+    if end - start <= 1.0e-12 {
+        return None;
+    }
+    Some((
+        interpolate(first, second, start),
+        interpolate(first, second, end),
+    ))
+}
+
+fn insert_unique_point(vertices: &mut Vec<[f64; 2]>, point: [f64; 2]) -> usize {
+    // Terrain positions are f32 throughout the native data and renderer. Make
+    // clipped crossings use that same representable coordinate before building
+    // the local CDT, so distinct f64 intersections cannot collapse to one XY
+    // vertex only after fragment emission.
+    let point = [f64::from(point[0] as f32), f64::from(point[1] as f32)];
+    if let Some(index) = vertices
+        .iter()
+        .position(|candidate| distance_squared(*candidate, point) <= 1.0e-18)
+    {
+        return index;
+    }
+    let index = vertices.len();
+    vertices.push(point);
+    index
+}
+
+fn insert_triangulation_point(
+    point_id: usize,
+    vertices: &[[f64; 2]],
+    triangles: &mut Vec<[usize; 3]>,
+) -> Result<(), ()> {
+    let point = vertices[point_id];
+    let edge = triangulation_edges(triangles)
+        .keys()
+        .copied()
+        .find(|&(first, second)| point_on_open_segment(point, vertices[first], vertices[second]));
+    if let Some(edge) = edge {
+        let mut split = Vec::new();
+        triangles.retain(|triangle| {
+            if triangle_has_edge(*triangle, edge) {
+                let opposite = triangle
+                    .iter()
+                    .copied()
+                    .find(|index| *index != edge.0 && *index != edge.1)
+                    .expect("triangle has opposite vertex");
+                split.push(opposite);
+                false
+            } else {
+                true
+            }
+        });
+        if split.is_empty() {
+            return Err(());
+        }
+        for opposite in split {
+            let first = make_ccw([edge.0, point_id, opposite], vertices).ok_or(())?;
+            let second = make_ccw([point_id, edge.1, opposite], vertices).ok_or(())?;
+            triangles.extend([first, second]);
+        }
+        return Ok(());
+    }
+
+    let containing = triangles.iter().position(|triangle| {
+        point_in_triangle(
+            point,
+            vertices[triangle[0]],
+            vertices[triangle[1]],
+            vertices[triangle[2]],
+        )
+    });
+    let Some(containing) = containing else {
+        return Err(());
+    };
+    let triangle = triangles.swap_remove(containing);
+    for edge in 0..3 {
+        let split = make_ccw(
+            [triangle[edge], triangle[(edge + 1) % 3], point_id],
+            vertices,
+        )
+        .ok_or(())?;
+        triangles.push(split);
+    }
+    Ok(())
+}
+
+fn source_boundary_constraints(
+    vertices: &[[f64; 2]],
+    triangle: [[f64; 2]; 3],
+) -> std::collections::BTreeSet<(usize, usize)> {
+    let mut constraints = std::collections::BTreeSet::new();
+    for edge in 0..3 {
+        let first = triangle[edge];
+        let second = triangle[(edge + 1) % 3];
+        let mut points = vertices
+            .iter()
+            .copied()
+            .enumerate()
+            .filter_map(|(index, point)| {
+                point_on_segment(point, first, second)
+                    .then_some((segment_parameter(point, first, second), index))
+            })
+            .collect::<Vec<_>>();
+        points.sort_by(|left, right| left.0.total_cmp(&right.0));
+        points.dedup_by_key(|(_, index)| *index);
+        for pair in points.windows(2) {
+            constraints.insert(canonical_edge(pair[0].1, pair[1].1));
+        }
+    }
+    constraints
+}
+
+fn recover_constraint(
+    target: (usize, usize),
+    triangles: &mut [[usize; 3]],
+    vertices: &[[f64; 2]],
+    constraints: &std::collections::BTreeSet<(usize, usize)>,
+) -> Result<(), ()> {
+    let limit = vertices
+        .len()
+        .saturating_mul(vertices.len())
+        .saturating_mul(8)
+        .max(1);
+    for _ in 0..limit {
+        if edge_exists(triangles, target) {
+            return Ok(());
+        }
+        let mut crossed = triangulation_edges(triangles)
+            .into_iter()
+            .filter_map(|(edge, adjacent)| {
+                if adjacent.len() != 2
+                    || edge.0 == target.0
+                    || edge.0 == target.1
+                    || edge.1 == target.0
+                    || edge.1 == target.1
+                {
+                    return None;
+                }
+                proper_intersection_parameter(
+                    vertices[target.0],
+                    vertices[target.1],
+                    vertices[edge.0],
+                    vertices[edge.1],
+                )
+                .map(|parameter| (parameter, edge, adjacent))
+            })
+            .collect::<Vec<_>>();
+        crossed.sort_by(|left, right| left.0.total_cmp(&right.0));
+        let mut flipped = false;
+        for (_, edge, adjacent) in crossed {
+            if constraints.contains(&edge) {
+                return Err(());
+            }
+            let first_index = adjacent[0].0;
+            let second_index = adjacent[1].0;
+            let first_opposite = adjacent[0].1;
+            let second_opposite = adjacent[1].1;
+            if !convex_quadrilateral(
+                vertices[edge.0],
+                vertices[edge.1],
+                vertices[first_opposite],
+                vertices[second_opposite],
+            ) {
+                continue;
+            }
+            let replacement = canonical_edge(first_opposite, second_opposite);
+            if constraints.contains(&replacement)
+                || proper_intersection_parameter(
+                    vertices[target.0],
+                    vertices[target.1],
+                    vertices[replacement.0],
+                    vertices[replacement.1],
+                )
+                .is_some()
+            {
+                continue;
+            }
+            let (a, b) = directed_edge(triangles[first_index], edge).ok_or(())?;
+            triangles[first_index] =
+                make_ccw([first_opposite, second_opposite, b], vertices).ok_or(())?;
+            triangles[second_index] =
+                make_ccw([second_opposite, first_opposite, a], vertices).ok_or(())?;
+            flipped = true;
+            break;
+        }
+        if !flipped {
+            return Err(());
+        }
+    }
+    Err(())
+}
+
+fn legalize_constrained_delaunay(
+    triangles: &mut [[usize; 3]],
+    vertices: &[[f64; 2]],
+    constraints: &std::collections::BTreeSet<(usize, usize)>,
+) -> Result<(), ()> {
+    let limit = triangles
+        .len()
+        .saturating_mul(triangles.len())
+        .saturating_mul(8)
+        .max(1);
+    for _ in 0..limit {
+        let edges = triangulation_edges(triangles);
+        let mut flipped = false;
+        for (edge, adjacent) in edges {
+            if adjacent.len() != 2 || constraints.contains(&edge) {
+                continue;
+            }
+            let first_index = adjacent[0].0;
+            let second_index = adjacent[1].0;
+            let c = adjacent[0].1;
+            let d = adjacent[1].1;
+            if !convex_quadrilateral(vertices[edge.0], vertices[edge.1], vertices[c], vertices[d]) {
+                continue;
+            }
+            let (a, b) = directed_edge(triangles[first_index], edge).ok_or(())?;
+            if in_circle(vertices[a], vertices[b], vertices[c], vertices[d]) <= 1.0e-10 {
+                continue;
+            }
+            let replacement = canonical_edge(c, d);
+            if constraints.iter().any(|constraint| {
+                proper_intersection_parameter(
+                    vertices[replacement.0],
+                    vertices[replacement.1],
+                    vertices[constraint.0],
+                    vertices[constraint.1],
+                )
+                .is_some()
+            }) {
+                continue;
+            }
+            triangles[first_index] = make_ccw([c, d, b], vertices).ok_or(())?;
+            triangles[second_index] = make_ccw([d, c, a], vertices).ok_or(())?;
+            flipped = true;
+            break;
+        }
+        if !flipped {
+            return Ok(());
+        }
+    }
+    Err(())
+}
+
+fn triangulation_edges(
+    triangles: &[[usize; 3]],
+) -> std::collections::BTreeMap<(usize, usize), Vec<(usize, usize)>> {
+    let mut edges = std::collections::BTreeMap::new();
+    for (triangle_index, triangle) in triangles.iter().copied().enumerate() {
+        for edge in 0..3 {
+            let first = triangle[edge];
+            let second = triangle[(edge + 1) % 3];
+            let opposite = triangle[(edge + 2) % 3];
+            edges
+                .entry(canonical_edge(first, second))
+                .or_insert_with(Vec::new)
+                .push((triangle_index, opposite));
+        }
+    }
+    edges
+}
+
+fn triangle_has_edge(triangle: [usize; 3], edge: (usize, usize)) -> bool {
+    triangle.contains(&edge.0) && triangle.contains(&edge.1)
+}
+
+fn edge_exists(triangles: &[[usize; 3]], edge: (usize, usize)) -> bool {
+    triangles
+        .iter()
+        .copied()
+        .any(|triangle| triangle_has_edge(triangle, edge))
+}
+
+fn directed_edge(triangle: [usize; 3], edge: (usize, usize)) -> Option<(usize, usize)> {
+    (0..3).find_map(|index| {
+        let first = triangle[index];
+        let second = triangle[(index + 1) % 3];
+        (canonical_edge(first, second) == edge).then_some((first, second))
+    })
+}
+
+fn canonical_edge(first: usize, second: usize) -> (usize, usize) {
+    if first < second {
+        (first, second)
+    } else {
+        (second, first)
+    }
+}
+
+fn make_ccw(triangle: [usize; 3], vertices: &[[f64; 2]]) -> Option<[usize; 3]> {
+    let area = cross2(
+        vertices[triangle[0]],
+        vertices[triangle[1]],
+        vertices[triangle[2]],
+    );
+    if !area.is_finite() || area.abs() <= 1.0e-12 {
+        return None;
+    }
+    Some(if area > 0.0 {
+        triangle
+    } else {
+        [triangle[0], triangle[2], triangle[1]]
+    })
+}
+
+fn point_in_triangle(point: [f64; 2], a: [f64; 2], b: [f64; 2], c: [f64; 2]) -> bool {
+    let edges = [(a, b), (b, c), (c, a)];
+    edges.into_iter().all(|(first, second)| {
+        let scale = distance_squared(first, second).sqrt();
+        cross2(first, second, point) >= -linear_tolerance(first, second, point) * scale
+    })
+}
+
+fn point_on_open_segment(point: [f64; 2], first: [f64; 2], second: [f64; 2]) -> bool {
+    if !point_on_segment(point, first, second) {
+        return false;
+    }
+    let t = segment_parameter(point, first, second);
+    t > 1.0e-9 && t < 1.0 - 1.0e-9
+}
+
+fn point_on_segment(point: [f64; 2], first: [f64; 2], second: [f64; 2]) -> bool {
+    let length_squared = distance_squared(first, second);
+    if length_squared <= 1.0e-20 {
+        return distance_squared(point, first) <= 1.0e-18;
+    }
+    let cross = cross2(first, second, point).abs();
+    let dot = (point[0] - first[0]) * (second[0] - first[0])
+        + (point[1] - first[1]) * (second[1] - first[1]);
+    let tolerance = linear_tolerance(first, second, point);
+    let length = length_squared.sqrt();
+    cross <= tolerance * length
+        && dot >= -tolerance * length
+        && dot <= length_squared + tolerance * length
+}
+
+fn linear_tolerance(first: [f64; 2], second: [f64; 2], point: [f64; 2]) -> f64 {
+    let coordinate_scale = first
+        .into_iter()
+        .chain(second)
+        .chain(point)
+        .map(f64::abs)
+        .fold(1.0_f64, f64::max);
+    2.0 * f64::from(f32::EPSILON) * coordinate_scale
+}
+
+fn segment_parameter(point: [f64; 2], first: [f64; 2], second: [f64; 2]) -> f64 {
+    let delta = [second[0] - first[0], second[1] - first[1]];
+    let length_squared = delta[0] * delta[0] + delta[1] * delta[1];
+    if length_squared <= 1.0e-20 {
+        0.0
+    } else {
+        ((point[0] - first[0]) * delta[0] + (point[1] - first[1]) * delta[1]) / length_squared
+    }
+}
+
+fn proper_intersection_parameter(
+    first: [f64; 2],
+    second: [f64; 2],
+    edge_first: [f64; 2],
+    edge_second: [f64; 2],
+) -> Option<f64> {
+    let first_delta = [second[0] - first[0], second[1] - first[1]];
+    let edge_delta = [
+        edge_second[0] - edge_first[0],
+        edge_second[1] - edge_first[1],
+    ];
+    let denominator = first_delta[0] * edge_delta[1] - first_delta[1] * edge_delta[0];
+    if denominator.abs() <= 1.0e-12 {
+        return None;
+    }
+    let offset = [edge_first[0] - first[0], edge_first[1] - first[1]];
+    let t = (offset[0] * edge_delta[1] - offset[1] * edge_delta[0]) / denominator;
+    let u = (offset[0] * first_delta[1] - offset[1] * first_delta[0]) / denominator;
+    (t > 1.0e-9 && t < 1.0 - 1.0e-9 && u > 1.0e-9 && u < 1.0 - 1.0e-9).then_some(t)
+}
+
+fn convex_quadrilateral(a: [f64; 2], b: [f64; 2], c: [f64; 2], d: [f64; 2]) -> bool {
+    cross2(a, b, c) * cross2(a, b, d) < -1.0e-12 && cross2(c, d, a) * cross2(c, d, b) < -1.0e-12
+}
+
+fn in_circle(a: [f64; 2], b: [f64; 2], c: [f64; 2], d: [f64; 2]) -> f64 {
+    let a = [a[0] - d[0], a[1] - d[1]];
+    let b = [b[0] - d[0], b[1] - d[1]];
+    let c = [c[0] - d[0], c[1] - d[1]];
+    let a2 = a[0] * a[0] + a[1] * a[1];
+    let b2 = b[0] * b[0] + b[1] * b[1];
+    let c2 = c[0] * c[0] + c[1] * c[1];
+    a2 * (b[0] * c[1] - b[1] * c[0]) - b2 * (a[0] * c[1] - a[1] * c[0])
+        + c2 * (a[0] * b[1] - a[1] * b[0])
+}
+
+fn segments_intersect(a: [f64; 2], b: [f64; 2], c: [f64; 2], d: [f64; 2]) -> bool {
+    let ab_c = cross2(a, b, c);
+    let ab_d = cross2(a, b, d);
+    let cd_a = cross2(c, d, a);
+    let cd_b = cross2(c, d, b);
+    let eps = 1.0e-8;
+    (ab_c.abs() <= eps && point_on_segment(c, a, b))
+        || (ab_d.abs() <= eps && point_on_segment(d, a, b))
+        || (cd_a.abs() <= eps && point_on_segment(a, c, d))
+        || (cd_b.abs() <= eps && point_on_segment(b, c, d))
+        || (ab_c * ab_d < 0.0 && cd_a * cd_b < 0.0)
+}
+
+fn point_in_ring(point: [f64; 2], ring: &[[f64; 2]]) -> bool {
+    let mut winding = 0_i32;
+    for index in 0..ring.len() {
+        let first = ring[index];
+        let second = ring[(index + 1) % ring.len()];
+        let side = cross2(first, second, point);
+        if first[1] <= point[1] {
+            if second[1] > point[1] && side > 0.0 {
+                winding += 1;
+            }
+        } else if second[1] <= point[1] && side < 0.0 {
+            winding -= 1;
+        }
+    }
+    winding != 0
+}
+
+fn polygon_area(points: &[[f64; 2]]) -> f64 {
+    (0..points.len())
+        .map(|index| {
+            let first = points[index];
+            let second = points[(index + 1) % points.len()];
+            first[0] * second[1] - second[0] * first[1]
+        })
+        .sum::<f64>()
+        * 0.5
+}
+
+fn lift_to_source_plane(point: [f64; 2], source: [[f32; 3]; 3]) -> Option<[f32; 3]> {
+    for vertex in source {
+        if distance_squared(point, [f64::from(vertex[0]), f64::from(vertex[1])]) <= 1.0e-18 {
+            return Some(vertex);
+        }
+    }
+    let a = source[0].map(f64::from);
+    let b = source[1].map(f64::from);
+    let c = source[2].map(f64::from);
+    let denominator = cross2([a[0], a[1]], [b[0], b[1]], [c[0], c[1]]);
+    if denominator.abs() <= 1.0e-10 {
+        return None;
+    }
+    let second_weight = cross2([a[0], a[1]], point, [c[0], c[1]]) / denominator;
+    let third_weight = cross2([a[0], a[1]], [b[0], b[1]], point) / denominator;
+    let first_weight = 1.0 - second_weight - third_weight;
+    let z = first_weight * a[2] + second_weight * b[2] + third_weight * c[2];
+    let position = [point[0] as f32, point[1] as f32, z as f32];
+    position
+        .iter()
+        .all(|value| value.is_finite())
+        .then_some(position)
+}
+
+fn cross2(first: [f64; 2], second: [f64; 2], third: [f64; 2]) -> f64 {
+    (second[0] - first[0]) * (third[1] - first[1]) - (second[1] - first[1]) * (third[0] - first[0])
+}
+
+fn interpolate(first: [f64; 2], second: [f64; 2], amount: f64) -> [f64; 2] {
+    [
+        first[0] + (second[0] - first[0]) * amount,
+        first[1] + (second[1] - first[1]) * amount,
+    ]
+}
+
+fn distance_squared(first: [f64; 2], second: [f64; 2]) -> f64 {
+    let delta_x = first[0] - second[0];
+    let delta_y = first[1] - second[1];
+    delta_x * delta_x + delta_y * delta_y
+}
+
+#[cfg(test)]
+mod tests {
+    use super::subtract_contour;
+
+    #[test]
+    fn translated_thin_outside_strip_keeps_its_area() {
+        let base = 10_000.0_f32;
+        let source = [[
+            [base, base, 0.0],
+            [base + 10.0, base, 0.0],
+            [base, base + 10.0, 0.0],
+        ]];
+        let outer = [
+            [base + 0.003, base + 0.003],
+            [base + 1.0, base + 0.003],
+            [base + 1.0, base + 0.997],
+            [base + 0.003, base + 0.997],
+        ];
+        let output = subtract_contour(&source, &outer).expect("simple contour");
+        let retained_area = output
+            .iter()
+            .map(|fragment| {
+                let [a, b, c] = fragment.positions;
+                let ab = [f64::from(b[0] - a[0]), f64::from(b[1] - a[1])];
+                let ac = [f64::from(c[0] - a[0]), f64::from(c[1] - a[1])];
+                (ab[0] * ac[1] - ab[1] * ac[0]).abs() * 0.5
+            })
+            .sum::<f64>();
+        let expected_area =
+            50.0 - (f64::from(outer[1][0] - outer[0][0]) * f64::from(outer[2][1] - outer[1][1]));
+        assert!(
+            (retained_area - expected_area).abs() < 0.001,
+            "retained area {retained_area}, expected {expected_area}"
+        );
+    }
+}
diff --git a/crates/fparkan-terrain/src/foundation.rs b/crates/fparkan-terrain/src/foundation.rs
new file mode 100644
index 0000000..594bbfe
--- /dev/null
+++ b/crates/fparkan-terrain/src/foundation.rs
@@ -0,0 +1,1806 @@
+//! Native-shaped building basement foundation geometry.
+
+use crate::contour::subtract_contour;
+use fparkan_terrain_format::{FullSurfaceMask, LandMeshDocument, TerrainLod};
+
+const LINE_EPSILON: f64 = 1.0e-5;
+const POINT_EPSILON_SQUARED: f64 = 1.0e-12;
+
+/// A terrain or generated basement triangle shared by preview consumers.
+#[derive(Clone, Copy, Debug, PartialEq)]
+pub struct BasementFoundationFragment {
+    /// Original Land.msh face for terrain fragments; `None` for generated band faces.
+    pub source_face: Option<usize>,
+    /// Raw terrain or generated basement face flags.
+    pub surface_mask: FullSurfaceMask,
+    /// Triangle coordinates shared by rendering, shadows, and terrain queries.
+    pub positions: [[f32; 3]; 3],
+}
+
+/// Canonical output of the native-shaped terrain basement preparation.
+#[derive(Clone, Debug, PartialEq)]
+pub struct BasementFoundation {
+    pub fragments: Vec<BasementFoundationFragment>,
+    /// Native placement-height inputs in footprint order.
+    pub placements: Vec<BasementPlacement>,
+}
+
+/// One building footprint in native preparation order.
+#[derive(Clone, Copy, Debug)]
+pub struct BasementFootprint<'a> {
+    /// Authored outside ring in world XY coordinates.
+    pub outer: &'a [[f32; 2]],
+    /// Authored inner ring, retaining its source Z for placement calculations.
+    pub inner: &'a [[f32; 3]],
+    /// Whether native getter mode flattens the generated inner boundary to mean outer Z.
+    pub flatten_inner_to_mean_outer_z: bool,
+}
+
+/// Native placement-height inputs captured for one building.
+#[derive(Clone, Copy, Debug, PartialEq)]
+pub struct BasementPlacement {
+    /// Native transform adjustment: first authored inner Z minus mean sampled outer terrain Z.
+    pub placement_delta: f32,
+}
+
+/// Error returned while preparing a basement foundation.
+#[derive(Clone, Debug, PartialEq)]
+pub enum BasementFoundationError {
+    /// Terrain Node38/Slot68 data cannot resolve LOD0 faces.
+    InvalidLodMapping(fparkan_terrain_format::TerrainLodError),
+    /// A selected LOD0 face refers to a missing vertex.
+    InvalidTerrainFace { face: usize, vertex: u16 },
+    /// A contour, nested-ring relationship, or generated annulus is invalid.
+    InvalidContour,
+    /// The LOD0 mesh has no usable terrain faces.
+    EmptyLod0,
+    /// The outer contour did not cross any LOD0 face.
+    NoOuterCrossings,
+    /// A point of the authored outer contour lies outside the LOD0 terrain.
+    OuterPointOutsideLod0 { point: usize },
+    /// An authored outer segment leaves the current terrain surface.
+    UncoveredOuterSegment { segment: usize },
+    /// Footprint outer contours overlap; native overlap semantics remain unverified.
+    OverlappingFootprints,
+}
+
+impl std::fmt::Display for BasementFoundationError {
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        match self {
+            Self::InvalidLodMapping(error) => write!(f, "invalid terrain LOD mapping: {error}"),
+            Self::InvalidTerrainFace { face, vertex } => {
+                write!(f, "terrain face {face} references missing vertex {vertex}")
+            }
+            Self::InvalidContour => write!(f, "invalid basement contour geometry"),
+            Self::EmptyLod0 => write!(f, "terrain has no LOD0 faces"),
+            Self::NoOuterCrossings => write!(f, "outer contour does not cross LOD0 terrain"),
+            Self::OuterPointOutsideLod0 { point } => {
+                write!(f, "outer contour point {point} is outside LOD0 terrain")
+            }
+            Self::UncoveredOuterSegment { segment } => {
+                write!(
+                    f,
+                    "outer contour segment {segment} leaves the terrain surface"
+                )
+            }
+            Self::OverlappingFootprints => {
+                write!(f, "overlapping basement footprints are not supported")
+            }
+        }
+    }
+}
+
+impl std::error::Error for BasementFoundationError {
+    fn source(&self) -> Option<&(dyn std::error::Error + 'static)> {
+        match self {
+            Self::InvalidLodMapping(error) => Some(error),
+            Self::InvalidTerrainFace { .. }
+            | Self::InvalidContour
+            | Self::EmptyLod0
+            | Self::NoOuterCrossings => None,
+            Self::OuterPointOutsideLod0 { .. } | Self::UncoveredOuterSegment { .. } => None,
+            Self::OverlappingFootprints => None,
+        }
+    }
+}
+
+/// Строит основания ландшафта в порядке нативной подготовки размещений.
+///
+/// Вершины внешнего контура нормализуются по топологии LOD0: сначала проверяются
+/// вершины граней и соседних граней, затем применяется проекция, используемая
+/// нативным `FUN_10006240`. Это косинусная формула с отдельной проверкой
+/// расстояния, а не обычная ортогональная проекция. При обходе сегмента
+/// вставляются только вершины LOD0, прошедшие тот же нативный предикат и
+/// находящиеся не ближе 0.5 мировых единиц к обоим концам; геометрическое
+/// совпадение с отрезком само по себе недостаточно.
+///
+/// Для каждого контура функция упорядочивает пересечения с плоскостями текущих
+/// канонических треугольников, вычисляет среднюю высоту внешней границы и строит
+/// ограниченную триангуляцию кольца. Явный флаг выбирает нативную ветвь высоты
+/// внутренней границы. Непересекающиеся контуры последовательно вырезаются из
+/// текущей поверхности; пересекающиеся контуры отклоняются, поскольку их
+/// взаимодействие в нативном коде пока не установлено.
+///
+/// Возвращённый набор фрагментов является общей геометрией для рендеринга,
+/// теней и запросов поверхности. Фрагменты исходной поверхности сохраняют
+/// номера и маски исходных граней; у сгенерированных граней маска соответствует
+/// записи `AddBasementFaceProc`.
+///
+/// The returned fragment vector is the canonical geometry for rendering,
+/// shadows, and terrain queries. Source faces preserve their original face IDs
+/// and masks. Generated faces use the raw mask written by AddBasementFaceProc.
+pub fn build_basement_foundations(
+    mesh: &LandMeshDocument,
+    footprints: &[BasementFootprint<'_>],
+) -> Result<BasementFoundation, BasementFoundationError> {
+    let assignments = mesh
+        .face_lod_assignments()
+        .map_err(BasementFoundationError::InvalidLodMapping)?;
+    let mut lod0_face_indices = Vec::new();
+    let mut terrain = Vec::new();
+    for assignment in assignments
+        .iter()
+        .filter(|assignment| assignment.lod == TerrainLod::Lod0)
+    {
+        let face = &mesh.faces[assignment.face_index];
+        let positions = terrain_face_positions(mesh, assignment.face_index)?;
+        terrain.push(BasementFoundationFragment {
+            source_face: Some(assignment.face_index),
+            surface_mask: face.flags,
+            positions,
+        });
+        lod0_face_indices.push(assignment.face_index);
+    }
+    if terrain.is_empty() {
+        return Err(BasementFoundationError::EmptyLod0);
+    }
+
+    let mut bands = Vec::new();
+    let mut placements = Vec::with_capacity(footprints.len());
+    let mut previous_outer_rings: Vec<Vec<[f64; 2]>> = Vec::with_capacity(footprints.len());
+    for footprint in footprints {
+        let (outer_points, outer_closed) = normalize_ring(footprint.outer)?;
+        let snapped_outer = normalize_outer_points(&outer_points, &lod0_face_indices, mesh)?;
+        if !simple_ring(&snapped_outer) {
+            return Err(BasementFoundationError::InvalidContour);
+        }
+        if previous_outer_rings
+            .iter()
+            .any(|previous| rings_overlap(previous, &snapped_outer))
+        {
+            return Err(BasementFoundationError::OverlappingFootprints);
+        }
+        previous_outer_rings.push(snapped_outer.clone());
+        let mut outer_ring = snapped_outer
+            .iter()
+            .map(|point| [point[0] as f32, point[1] as f32])
+            .collect::<Vec<_>>();
+        if outer_closed {
+            outer_ring.push(outer_ring[0]);
+        }
+        let samples = terrain
+            .iter()
+            .map(|fragment| {
+                (
+                    fragment
+                        .source_face
+                        .expect("terrain fragments have source faces"),
+                    fragment.surface_mask,
+                    fragment.positions,
+                )
+            })
+            .collect::<Vec<_>>();
+        let mut outer_crossings = contour_crossings(&snapped_outer, &samples)?;
+        if outer_crossings.is_empty() {
+            return Err(BasementFoundationError::NoOuterCrossings);
+        }
+        // Crossings are open here; restore the authored duplicate only afterward.
+        let mean_outer_z = outer_crossings
+            .iter()
+            .map(|crossing| f64::from(crossing[2]))
+            .sum::<f64>()
+            / outer_crossings.len() as f64;
+        if !mean_outer_z.is_finite() {
+            return Err(BasementFoundationError::InvalidContour);
+        }
+        let mean_outer_z = mean_outer_z as f32;
+        if outer_closed {
+            let closing = *outer_crossings
+                .first()
+                .ok_or(BasementFoundationError::NoOuterCrossings)?;
+            outer_crossings.push(closing);
+        }
+
+        let inner_xy = footprint
+            .inner
+            .iter()
+            .map(|point| [point[0], point[1]])
+            .collect::<Vec<_>>();
+        let (inner_points, inner_closed) = normalize_ring(&inner_xy)?;
+        if footprint
+            .inner
+            .iter()
+            .flatten()
+            .any(|value| !value.is_finite())
+            || !nested_rings(&snapped_outer, &inner_points)
+        {
+            return Err(BasementFoundationError::InvalidContour);
+        }
+        let authored_inner_z = footprint
+            .inner
+            .first()
+            .ok_or(BasementFoundationError::InvalidContour)?[2];
+        let placement_delta = authored_inner_z - mean_outer_z;
+        if !placement_delta.is_finite() {
+            return Err(BasementFoundationError::InvalidContour);
+        }
+        placements.push(BasementPlacement { placement_delta });
+        let inner_vertex_count = if inner_closed {
+            footprint.inner.len() - 1
+        } else {
+            footprint.inner.len()
+        };
+        let mut inner_ring = footprint.inner[..inner_vertex_count]
+            .iter()
+            .map(|point| {
+                [
+                    point[0],
+                    point[1],
+                    if footprint.flatten_inner_to_mean_outer_z {
+                        mean_outer_z
+                    } else {
+                        point[2]
+                    },
+                ]
+            })
+            .collect::<Vec<_>>();
+        if inner_closed {
+            inner_ring.push(inner_ring[0]);
+        }
+        let outer_surface = crossings_as_ring(&outer_crossings, outer_closed)?;
+        let inner_surface = inner_ring_without_closure(&inner_ring);
+        for positions in triangulate_annulus(&outer_surface, inner_surface)? {
+            bands.push(BasementFoundationFragment {
+                source_face: None,
+                surface_mask: basement_surface_mask(positions),
+                positions,
+            });
+        }
+
+        let source_positions = terrain
+            .iter()
+            .map(|fragment| fragment.positions)
+            .collect::<Vec<_>>();
+        let outside = subtract_contour(&source_positions, &outer_ring)
+            .map_err(|_| BasementFoundationError::InvalidContour)?;
+        let clipped = outside
+            .into_iter()
+            .map(|fragment| {
+                let source = terrain[fragment.source_triangle];
+                BasementFoundationFragment {
+                    source_face: source.source_face,
+                    surface_mask: source.surface_mask,
+                    positions: fragment.positions,
+                }
+            })
+            .collect();
+        terrain = clipped;
+    }
+    terrain.extend(bands);
+    Ok(BasementFoundation {
+        fragments: terrain,
+        placements,
+    })
+}
+
+fn terrain_face_positions(
+    mesh: &LandMeshDocument,
+    face_index: usize,
+) -> Result<[[f32; 3]; 3], BasementFoundationError> {
+    let face = &mesh.faces[face_index];
+    let resolve = |vertex: u16| {
+        mesh.positions.get(usize::from(vertex)).copied().ok_or(
+            BasementFoundationError::InvalidTerrainFace {
+                face: face_index,
+                vertex,
+            },
+        )
+    };
+    Ok([
+        resolve(face.vertices[0])?,
+        resolve(face.vertices[1])?,
+        resolve(face.vertices[2])?,
+    ])
+}
+
+fn normalize_ring(ring: &[[f32; 2]]) -> Result<(Vec<[f64; 2]>, bool), BasementFoundationError> {
+    if ring.iter().flatten().any(|value| !value.is_finite()) {
+        return Err(BasementFoundationError::InvalidContour);
+    }
+    let mut points = ring
+        .iter()
+        .map(|point| [f64::from(point[0]), f64::from(point[1])])
+        .collect::<Vec<_>>();
+    let closed = points.len() > 1
+        && distance_squared(points[0], *points.last().expect("nonempty")) <= POINT_EPSILON_SQUARED;
+    if closed {
+        points.pop();
+    }
+    if points.len() < 3 || !simple_ring(&points) {
+        return Err(BasementFoundationError::InvalidContour);
+    }
+    Ok((points, closed))
+}
+
+fn normalize_outer_points(
+    points: &[[f64; 2]],
+    lod0_faces: &[usize],
+    mesh: &LandMeshDocument,
+) -> Result<Vec<[f64; 2]>, BasementFoundationError> {
+    // shortcut: first-face lookup and neighbor ordering still use source LOD0
+    // IDs for every placement. The crossing sampler and cutter use the current
+    // canonical fragment list, but dynamic neighbor-ID reconstruction is not
+    // modeled yet, so later placements are not claimed as native-identical.
+    let mut output = points
+        .iter()
+        .enumerate()
+        .map(|(point_index, point)| {
+            native_normalize_outer_point(
+                mesh,
+                lod0_faces,
+                [point[0] as f32, point[1] as f32],
+                point_index,
+            )
+        })
+        .collect::<Result<Vec<_>, _>>()?
+        .into_iter()
+        .map(|point| [f64::from(point[0]), f64::from(point[1])])
+        .collect::<Vec<_>>();
+
+    // The native contour pass asks for terrain faces overlapping the current
+    // segment's AABB, then tries their three vertices in face order. On one
+    // accepted vertex it inserts the full terrain point and resumes from the
+    // shortened segment at the same index.
+    let mut segment_index = 0;
+    while segment_index < output.len() {
+        let first = output[segment_index].map(|value| value as f32);
+        let second = output[(segment_index + 1) % output.len()].map(|value| value as f32);
+        let min = [first[0].min(second[0]), first[1].min(second[1])];
+        let max = [first[0].max(second[0]), first[1].max(second[1])];
+        let mut inserted = false;
+        'faces: for &face_index in lod0_faces {
+            let Some(face) = mesh.faces.get(face_index) else {
+                continue;
+            };
+            let triangle = face
+                .vertices
+                .map(|vertex| mesh.positions.get(usize::from(vertex)).copied());
+            let [Some(a), Some(b), Some(c)] = triangle else {
+                continue;
+            };
+            let face_min = [a[0].min(b[0]).min(c[0]), a[1].min(b[1]).min(c[1])];
+            let face_max = [a[0].max(b[0]).max(c[0]), a[1].max(b[1]).max(c[1])];
+            if (0..2).any(|axis| face_max[axis] < min[axis] || max[axis] < face_min[axis]) {
+                continue;
+            }
+            for vertex in [a, b, c] {
+                if !native_segment_inserts_vertex(first, second, vertex) {
+                    continue;
+                }
+                let vertex_xy = [f64::from(vertex[0]), f64::from(vertex[1])];
+                if output
+                    .iter()
+                    .any(|existing| distance_squared(*existing, vertex_xy) <= POINT_EPSILON_SQUARED)
+                {
+                    continue;
+                }
+                output.insert(segment_index + 1, vertex_xy);
+                inserted = true;
+                break 'faces;
+            }
+        }
+        if !inserted {
+            segment_index += 1;
+        }
+    }
+    Ok(output)
+}
+
+fn native_normalize_outer_point(
+    mesh: &LandMeshDocument,
+    lod0_faces: &[usize],
+    point: [f32; 2],
+    point_index: usize,
+) -> Result<[f32; 2], BasementFoundationError> {
+    let containing_face = lod0_faces.iter().copied().find(|&face_index| {
+        let Some(face) = mesh.faces.get(face_index) else {
+            return false;
+        };
+        let triangle = face
+            .vertices
+            .map(|vertex| mesh.positions[usize::from(vertex)]);
+        point_in_triangle_xy(point, triangle)
+    });
+    let containing_face = containing_face
+        .ok_or(BasementFoundationError::OuterPointOutsideLod0 { point: point_index })?;
+    let mut candidates = vec![containing_face];
+    for neighbor in mesh.faces[containing_face].neighbors.iter().flatten() {
+        let neighbor = usize::from(*neighbor);
+        if lod0_faces.contains(&neighbor) && !candidates.contains(&neighbor) {
+            candidates.push(neighbor);
+        }
+    }
+    let normalized = point;
+    for face_index in candidates {
+        let face = &mesh.faces[face_index];
+        let vertices = face
+            .vertices
+            .map(|vertex| mesh.positions.get(usize::from(vertex)).copied());
+        let [Some(a), Some(b), Some(c)] = vertices else {
+            return Err(BasementFoundationError::InvalidTerrainFace {
+                face: face_index,
+                vertex: face.vertices[0],
+            });
+        };
+        let vertices = [a, b, c];
+        for (first, second) in [(0, 1), (1, 2), (2, 0)] {
+            if let Some(candidate) =
+                native_normalize_edge(normalized, vertices[first], vertices[second])?
+            {
+                return Ok(candidate);
+            }
+        }
+    }
+    Ok(normalized)
+}
+
+fn point_in_triangle_xy(point: [f32; 2], triangle: [[f32; 3]; 3]) -> bool {
+    let sign = |a: [f32; 3], b: [f32; 3]| {
+        (point[0] - b[0]) * (a[1] - b[1]) - (a[0] - b[0]) * (point[1] - b[1])
+    };
+    let d1 = sign(triangle[0], triangle[1]);
+    let d2 = sign(triangle[1], triangle[2]);
+    let d3 = sign(triangle[2], triangle[0]);
+    let has_negative = d1 < -1.0e-4 || d2 < -1.0e-4 || d3 < -1.0e-4;
+    let has_positive = d1 > 1.0e-4 || d2 > 1.0e-4 || d3 > 1.0e-4;
+    !(has_negative && has_positive)
+}
+
+fn native_normalize_edge(
+    point: [f32; 2],
+    first: [f32; 3],
+    second: [f32; 3],
+) -> Result<Option<[f32; 2]>, BasementFoundationError> {
+    if distance_xy(point, [first[0], first[1]]) < 0.5 {
+        let point = [first[0], first[1]];
+        return Ok(Some(point));
+    }
+    if distance_xy(point, [second[0], second[1]]) < 0.5 {
+        let point = [second[0], second[1]];
+        return Ok(Some(point));
+    }
+    let Some((offset, _)) = native_edge_projection(point, first, second) else {
+        return Ok(None);
+    };
+    let mut projected = [point[0] + offset[0], point[1] + offset[1]];
+    let Some((check_offset, _)) = native_edge_projection(projected, first, second) else {
+        return Err(BasementFoundationError::InvalidContour);
+    };
+    if check_offset[0].hypot(check_offset[1]) >= 0.0001 {
+        return Err(BasementFoundationError::InvalidContour);
+    }
+    if distance_xy(projected, [first[0], first[1]]) < 0.5 {
+        projected = [first[0], first[1]];
+    }
+    if distance_xy(projected, [second[0], second[1]]) < 0.5 {
+        projected = [second[0], second[1]];
+    }
+    Ok(Some(projected))
+}
+
+fn native_segment_inserts_vertex(first: [f32; 2], second: [f32; 2], point: [f32; 3]) -> bool {
+    if distance_xy([point[0], point[1]], first) < 0.5
+        || distance_xy([point[0], point[1]], second) < 0.5
+    {
+        return false;
+    }
+    native_edge_projection(
+        [point[0], point[1]],
+        [first[0], first[1], 0.0],
+        [second[0], second[1], 0.0],
+    )
+    .is_some()
+}
+
+fn native_edge_projection(
+    point: [f32; 2],
+    first: [f32; 3],
+    second: [f32; 3],
+) -> Option<([f32; 2], f32)> {
+    let ux = point[0] - first[0];
+    let uy = point[1] - first[1];
+    let ex = second[0] - first[0];
+    let ey = second[1] - first[1];
+    let edge_length = (ex * ex + ey * ey).sqrt();
+    let point_length = (ux * ux + uy * uy).sqrt();
+    if edge_length == 0.0 || point_length == 0.0 {
+        return None;
+    }
+    let t = ((ex * ux + ey * uy) / edge_length) / point_length;
+    if !(0.0..=1.0).contains(&t) {
+        return None;
+    }
+    let offset = [
+        (t * ex + first[0]) - point[0],
+        (t * ey + first[1]) - point[1],
+    ];
+    (offset[0].hypot(offset[1]) < 0.5).then_some((offset, t))
+}
+
+fn distance_xy(first: [f32; 2], second: [f32; 2]) -> f32 {
+    ((first[0] - second[0]).powi(2) + (first[1] - second[1]).powi(2)).sqrt()
+}
+
+fn contour_crossings(
+    ring: &[[f64; 2]],
+    terrain_faces: &[(usize, FullSurfaceMask, [[f32; 3]; 3])],
+) -> Result<Vec<[f32; 3]>, BasementFoundationError> {
+    let mut crossings = Vec::new();
+    for segment_index in 0..ring.len() {
+        let first = ring[segment_index];
+        let second = ring[(segment_index + 1) % ring.len()];
+        let mut segment_crossings = Vec::new();
+        let mut coverage = Vec::new();
+        for (source_face, _, triangle) in terrain_faces {
+            let Some((start, end)) = segment_triangle_interval(first, second, *triangle) else {
+                continue;
+            };
+            coverage.push((start, end));
+            for t in [start, end] {
+                let point = [
+                    first[0] + (second[0] - first[0]) * t,
+                    first[1] + (second[1] - first[1]) * t,
+                ];
+                let Some(position) = lift_to_face(point, *triangle) else {
+                    continue;
+                };
+                segment_crossings.push((t, *source_face, position));
+            }
+        }
+        segment_crossings.sort_by(|left, right| {
+            left.0
+                .total_cmp(&right.0)
+                .then_with(|| left.1.cmp(&right.1))
+        });
+        coverage.sort_by(|left, right| {
+            left.0
+                .total_cmp(&right.0)
+                .then_with(|| left.1.total_cmp(&right.1))
+        });
+        let mut covered_until = 0.0_f64;
+        for (start, end) in coverage {
+            if start > covered_until + 1.0e-8 {
+                return Err(BasementFoundationError::UncoveredOuterSegment {
+                    segment: segment_index,
+                });
+            }
+            covered_until = covered_until.max(end);
+        }
+        if covered_until < 1.0 - 1.0e-8 {
+            return Err(BasementFoundationError::UncoveredOuterSegment {
+                segment: segment_index,
+            });
+        }
+        for (_, _, position) in segment_crossings {
+            let duplicate_of_previous = crossings.last().is_some_and(|previous: &[f32; 3]| {
+                distance_squared(
+                    [f64::from(previous[0]), f64::from(previous[1])],
+                    [f64::from(position[0]), f64::from(position[1])],
+                ) <= POINT_EPSILON_SQUARED
+            });
+            if duplicate_of_previous {
+                continue;
+            }
+            crossings.push(position);
+        }
+    }
+    if crossings.len() > 1
+        && distance_squared(
+            [f64::from(crossings[0][0]), f64::from(crossings[0][1])],
+            [
+                f64::from(crossings.last().expect("nonempty")[0]),
+                f64::from(crossings.last().expect("nonempty")[1]),
+            ],
+        ) <= POINT_EPSILON_SQUARED
+    {
+        crossings.pop();
+    }
+    Ok(crossings)
+}
+
+fn segment_triangle_interval(
+    first: [f64; 2],
+    second: [f64; 2],
+    triangle: [[f32; 3]; 3],
+) -> Option<(f64, f64)> {
+    let vertices = triangle.map(|point| [f64::from(point[0]), f64::from(point[1])]);
+    let area = cross2(vertices[0], vertices[1], vertices[2]);
+    if !area.is_finite() || area.abs() <= 1.0e-12 {
+        return None;
+    }
+    let orientation = area.signum();
+    let delta = [second[0] - first[0], second[1] - first[1]];
+    let mut lower: f64 = 0.0;
+    let mut upper: f64 = 1.0;
+    for edge in 0..3 {
+        let edge_first = vertices[edge];
+        let edge_second = vertices[(edge + 1) % 3];
+        let start_distance = orientation * cross2(edge_first, edge_second, first);
+        let end = [first[0] + delta[0], first[1] + delta[1]];
+        let end_distance = orientation * cross2(edge_first, edge_second, end);
+        let slope = end_distance - start_distance;
+        if slope.abs() <= 1.0e-15 {
+            if start_distance < -LINE_EPSILON {
+                return None;
+            }
+            continue;
+        }
+        // Tolerance decides whether a nearly outside interval is accepted;
+        // the actual boundary sample remains on the geometric edge (distance 0).
+        let crossing = -start_distance / slope;
+        if slope > 0.0 {
+            lower = lower.max(crossing);
+        } else {
+            upper = upper.min(crossing);
+        }
+        if lower > upper {
+            return None;
+        }
+    }
+    (upper >= 0.0 && lower <= 1.0).then_some((lower.clamp(0.0, 1.0), upper.clamp(0.0, 1.0)))
+}
+
+fn lift_to_face(point: [f64; 2], triangle: [[f32; 3]; 3]) -> Option<[f32; 3]> {
+    let vertices = triangle.map(|vertex| vertex.map(f64::from));
+    let denominator = cross2(
+        [vertices[0][0], vertices[0][1]],
+        [vertices[1][0], vertices[1][1]],
+        [vertices[2][0], vertices[2][1]],
+    );
+    if denominator.abs() <= 1.0e-12 {
+        return None;
+    }
+    let weight_b = cross2(
+        [vertices[0][0], vertices[0][1]],
+        point,
+        [vertices[2][0], vertices[2][1]],
+    ) / denominator;
+    let weight_c = cross2(
+        [vertices[0][0], vertices[0][1]],
+        [vertices[1][0], vertices[1][1]],
+        point,
+    ) / denominator;
+    let weight_a = 1.0 - weight_b - weight_c;
+    let z = weight_a * vertices[0][2] + weight_b * vertices[1][2] + weight_c * vertices[2][2];
+    let result = [point[0] as f32, point[1] as f32, z as f32];
+    result
+        .iter()
+        .all(|value| value.is_finite())
+        .then_some(result)
+}
+
+fn crossings_as_ring(
+    crossings: &[[f32; 3]],
+    closed: bool,
+) -> Result<Vec<[f32; 3]>, BasementFoundationError> {
+    let length = crossings.len().saturating_sub(usize::from(closed));
+    let ring = crossings.iter().take(length).copied().collect::<Vec<_>>();
+    (ring.len() >= 3)
+        .then_some(ring)
+        .ok_or(BasementFoundationError::InvalidContour)
+}
+
+fn inner_ring_without_closure(ring: &[[f32; 3]]) -> Vec<[f32; 3]> {
+    let mut result = ring.to_vec();
+    if result.len() > 1 && same_xy(result[0], *result.last().expect("nonempty")) {
+        result.pop();
+    }
+    result
+}
+
+fn triangulate_annulus(
+    outer: &[[f32; 3]],
+    inner: Vec<[f32; 3]>,
+) -> Result<Vec<[[f32; 3]; 3]>, BasementFoundationError> {
+    let mut outer = remove_consecutive_duplicates(outer);
+    let mut inner = remove_consecutive_duplicates(&inner);
+    if outer.len() < 3 || inner.len() < 3 || !nested_rings_3d(&outer, &inner) {
+        return Err(BasementFoundationError::InvalidContour);
+    }
+    if signed_area_3d(&outer) < 0.0 {
+        outer.reverse();
+    }
+    if signed_area_3d(&inner) > 0.0 {
+        inner.reverse();
+    }
+    let polygon = bridge_hole(&outer, &inner)?;
+    let output =
+        legalize_constrained_edges(ear_clip_weakly_simple_polygon(&polygon)?, &outer, &inner)?;
+    if output.iter().any(|triangle| {
+        !triangle.iter().flatten().all(|value| value.is_finite())
+            || cross2(xy(triangle[0]), xy(triangle[1]), xy(triangle[2])).abs() <= 1.0e-9
+    }) {
+        return Err(BasementFoundationError::InvalidContour);
+    }
+    Ok(output)
+}
+
+fn bridge_hole(
+    outer: &[[f32; 3]],
+    inner: &[[f32; 3]],
+) -> Result<Vec<[f32; 3]>, BasementFoundationError> {
+    let hole_index = (0..inner.len())
+        .max_by(|left, right| {
+            inner[*left][0]
+                .total_cmp(&inner[*right][0])
+                .then_with(|| inner[*right][1].total_cmp(&inner[*left][1]))
+        })
+        .ok_or(BasementFoundationError::InvalidContour)?;
+    let hole = inner[hole_index];
+    let mut candidates = (0..outer.len())
+        .map(|index| (distance_squared(xy(hole), xy(outer[index])), index))
+        .collect::<Vec<_>>();
+    candidates.sort_by(|left, right| {
+        left.0
+            .total_cmp(&right.0)
+            .then_with(|| outer[right.1][0].total_cmp(&outer[left.1][0]))
+    });
+    let mut bridged_outer = None;
+    for (_, candidate_index) in candidates {
+        let candidate = outer[candidate_index];
+        if visible_bridge(xy(hole), xy(candidate), outer, inner) {
+            bridged_outer = Some((outer.to_vec(), candidate_index));
+            break;
+        }
+    }
+    let (outer, outer_index) = bridged_outer.ok_or(BasementFoundationError::InvalidContour)?;
+    let mut polygon = Vec::with_capacity(outer.len() + inner.len() + 2);
+    for step in 0..outer.len() {
+        let index = (outer_index + step) % outer.len();
+        polygon.push(outer[index]);
+        if step == 0 {
+            polygon.push(inner[hole_index]);
+            for hole_step in 1..inner.len() {
+                polygon.push(inner[(hole_index + hole_step) % inner.len()]);
+            }
+            polygon.push(inner[hole_index]);
+            polygon.push(outer[outer_index]);
+        }
+    }
+    Ok(polygon)
+}
+
+fn visible_bridge(
+    hole: [f64; 2],
+    outer: [f64; 2],
+    outer_ring: &[[f32; 3]],
+    inner_ring: &[[f32; 3]],
+) -> bool {
+    let midpoint = [(hole[0] + outer[0]) * 0.5, (hole[1] + outer[1]) * 0.5];
+    let outer_xy = outer_ring.iter().copied().map(xy).collect::<Vec<_>>();
+    let inner_xy = inner_ring.iter().copied().map(xy).collect::<Vec<_>>();
+    if !point_in_ring(midpoint, &outer_xy) || point_in_ring(midpoint, &inner_xy) {
+        return false;
+    }
+    if (0..outer_xy.len()).any(|index| {
+        segment_crosses_open_edge(
+            hole,
+            outer,
+            outer_xy[index],
+            outer_xy[(index + 1) % outer_xy.len()],
+        )
+    }) {
+        return false;
+    }
+    !(0..inner_xy.len()).any(|index| {
+        segment_crosses_open_edge(
+            hole,
+            outer,
+            inner_xy[index],
+            inner_xy[(index + 1) % inner_xy.len()],
+        )
+    })
+}
+
+fn segment_crosses_open_edge(
+    first: [f64; 2],
+    second: [f64; 2],
+    edge_first: [f64; 2],
+    edge_second: [f64; 2],
+) -> bool {
+    if same_point(first, edge_first)
+        || same_point(first, edge_second)
+        || same_point(second, edge_first)
+        || same_point(second, edge_second)
+    {
+        return false;
+    }
+    segments_intersect(first, second, edge_first, edge_second)
+}
+
+fn same_point(first: [f64; 2], second: [f64; 2]) -> bool {
+    distance_squared(first, second) <= POINT_EPSILON_SQUARED
+}
+
+fn ear_clip_weakly_simple_polygon(
+    polygon: &[[f32; 3]],
+) -> Result<Vec<[[f32; 3]; 3]>, BasementFoundationError> {
+    let mut remaining = (0..polygon.len()).collect::<Vec<_>>();
+    let mut triangles = Vec::with_capacity(polygon.len().saturating_sub(2));
+    while remaining.len() > 3 {
+        let mut ear = None;
+        for slot in 0..remaining.len() {
+            let previous = remaining[(slot + remaining.len() - 1) % remaining.len()];
+            let current = remaining[slot];
+            let next = remaining[(slot + 1) % remaining.len()];
+            let a = xy(polygon[previous]);
+            let b = xy(polygon[current]);
+            let c = xy(polygon[next]);
+            if cross2(a, b, c) <= 1.0e-9 || !diagonal_is_clear(previous, next, &remaining, polygon)
+            {
+                continue;
+            }
+            let contains_point = remaining.iter().copied().any(|candidate| {
+                candidate != previous
+                    && candidate != current
+                    && candidate != next
+                    && !same_xy(polygon[candidate], polygon[previous])
+                    && !same_xy(polygon[candidate], polygon[current])
+                    && !same_xy(polygon[candidate], polygon[next])
+                    && point_in_or_on_triangle(xy(polygon[candidate]), a, b, c)
+            });
+            if !contains_point {
+                ear = Some((slot, [polygon[previous], polygon[current], polygon[next]]));
+                break;
+            }
+        }
+        if let Some((slot, triangle)) = ear {
+            triangles.push(triangle);
+            remaining.remove(slot);
+        } else {
+            return Err(BasementFoundationError::InvalidContour);
+        }
+    }
+    if remaining.len() == 3 {
+        let triangle = [
+            polygon[remaining[0]],
+            polygon[remaining[1]],
+            polygon[remaining[2]],
+        ];
+        if cross2(xy(triangle[0]), xy(triangle[1]), xy(triangle[2])) > 1.0e-9 {
+            triangles.push(triangle);
+        }
+    }
+    (!triangles.is_empty())
+        .then_some(triangles)
+        .ok_or(BasementFoundationError::InvalidContour)
+}
+
+fn diagonal_is_clear(
+    first_index: usize,
+    second_index: usize,
+    polygon_indices: &[usize],
+    polygon: &[[f32; 3]],
+) -> bool {
+    let first = xy(polygon[first_index]);
+    let second = xy(polygon[second_index]);
+    for edge in 0..polygon_indices.len() {
+        let edge_first_index = polygon_indices[edge];
+        let edge_second_index = polygon_indices[(edge + 1) % polygon_indices.len()];
+        if edge_first_index == first_index
+            || edge_second_index == first_index
+            || edge_first_index == second_index
+            || edge_second_index == second_index
+        {
+            continue;
+        }
+        if segment_crosses_open_edge(
+            first,
+            second,
+            xy(polygon[edge_first_index]),
+            xy(polygon[edge_second_index]),
+        ) {
+            return false;
+        }
+    }
+    true
+}
+
+fn point_in_or_on_triangle(point: [f64; 2], a: [f64; 2], b: [f64; 2], c: [f64; 2]) -> bool {
+    let epsilon = 1.0e-9;
+    cross2(a, b, point) >= -epsilon
+        && cross2(b, c, point) >= -epsilon
+        && cross2(c, a, point) >= -epsilon
+}
+
+fn nested_rings(outer: &[[f64; 2]], inner: &[[f64; 2]]) -> bool {
+    if !simple_ring(outer) || !simple_ring(inner) {
+        return false;
+    }
+    if rings_intersect(outer, inner) {
+        return false;
+    }
+    point_in_ring(inner[0], outer)
+}
+
+fn nested_rings_3d(outer: &[[f32; 3]], inner: &[[f32; 3]]) -> bool {
+    let outer_xy = outer.iter().copied().map(xy).collect::<Vec<_>>();
+    let inner_xy = inner.iter().copied().map(xy).collect::<Vec<_>>();
+    nested_rings(&outer_xy, &inner_xy)
+}
+
+fn remove_consecutive_duplicates(ring: &[[f32; 3]]) -> Vec<[f32; 3]> {
+    let mut result = Vec::with_capacity(ring.len());
+    for point in ring {
+        if result
+            .last()
+            .is_none_or(|previous| !same_xy(*previous, *point))
+        {
+            result.push(*point);
+        }
+    }
+    if result.len() > 1 && same_xy(result[0], *result.last().expect("nonempty")) {
+        result.pop();
+    }
+    result
+}
+
+fn simple_ring(ring: &[[f64; 2]]) -> bool {
+    if ring.len() < 3 || ring.iter().flatten().any(|value| !value.is_finite()) {
+        return false;
+    }
+    let area = ring_area(ring);
+    if !area.is_finite() || area.abs() <= 1.0e-10 {
+        return false;
+    }
+    for first in 0..ring.len() {
+        let first_next = (first + 1) % ring.len();
+        if distance_squared(ring[first], ring[first_next]) <= POINT_EPSILON_SQUARED {
+            return false;
+        }
+        for second in (first + 1)..ring.len() {
+            let second_next = (second + 1) % ring.len();
+            if first == second || first_next == second || second_next == first {
+                continue;
+            }
+            if segments_intersect(
+                ring[first],
+                ring[first_next],
+                ring[second],
+                ring[second_next],
+            ) {
+                return false;
+            }
+        }
+    }
+    true
+}
+
+fn rings_intersect(first: &[[f64; 2]], second: &[[f64; 2]]) -> bool {
+    (0..first.len()).any(|first_index| {
+        (0..second.len()).any(|second_index| {
+            segments_intersect(
+                first[first_index],
+                first[(first_index + 1) % first.len()],
+                second[second_index],
+                second[(second_index + 1) % second.len()],
+            )
+        })
+    })
+}
+
+fn rings_overlap(first: &[[f64; 2]], second: &[[f64; 2]]) -> bool {
+    rings_intersect(first, second)
+        || point_in_ring(first[0], second)
+        || point_in_ring(second[0], first)
+}
+
+fn point_in_ring(point: [f64; 2], ring: &[[f64; 2]]) -> bool {
+    let mut inside = false;
+    for index in 0..ring.len() {
+        let first = ring[index];
+        let second = ring[(index + 1) % ring.len()];
+        if point_on_segment(point, first, second) {
+            return true;
+        }
+        if (first[1] > point[1]) != (second[1] > point[1]) {
+            let intersection_x =
+                first[0] + (point[1] - first[1]) * (second[0] - first[0]) / (second[1] - first[1]);
+            if point[0] < intersection_x {
+                inside = !inside;
+            }
+        }
+    }
+    inside
+}
+
+fn point_on_segment(point: [f64; 2], first: [f64; 2], second: [f64; 2]) -> bool {
+    let dx = second[0] - first[0];
+    let dy = second[1] - first[1];
+    let length = (dx * dx + dy * dy).sqrt();
+    cross2(first, second, point).abs() <= LINE_EPSILON * length
+        && point[0] >= first[0].min(second[0]) - LINE_EPSILON
+        && point[0] <= first[0].max(second[0]) + LINE_EPSILON
+        && point[1] >= first[1].min(second[1]) - LINE_EPSILON
+        && point[1] <= first[1].max(second[1]) + LINE_EPSILON
+}
+
+fn segments_intersect(a: [f64; 2], b: [f64; 2], c: [f64; 2], d: [f64; 2]) -> bool {
+    let ab_c = cross2(a, b, c);
+    let ab_d = cross2(a, b, d);
+    let cd_a = cross2(c, d, a);
+    let cd_b = cross2(c, d, b);
+    let epsilon = 1.0e-8;
+    (ab_c.abs() <= epsilon && on_segment(c, a, b, epsilon))
+        || (ab_d.abs() <= epsilon && on_segment(d, a, b, epsilon))
+        || (cd_a.abs() <= epsilon && on_segment(a, c, d, epsilon))
+        || (cd_b.abs() <= epsilon && on_segment(b, c, d, epsilon))
+        || (ab_c * ab_d < 0.0 && cd_a * cd_b < 0.0)
+}
+
+fn on_segment(point: [f64; 2], first: [f64; 2], second: [f64; 2], epsilon: f64) -> bool {
+    point[0] >= first[0].min(second[0]) - epsilon
+        && point[0] <= first[0].max(second[0]) + epsilon
+        && point[1] >= first[1].min(second[1]) - epsilon
+        && point[1] <= first[1].max(second[1]) + epsilon
+}
+
+fn ring_area(ring: &[[f64; 2]]) -> f64 {
+    (0..ring.len())
+        .map(|index| {
+            let first = ring[index];
+            let second = ring[(index + 1) % ring.len()];
+            first[0] * second[1] - second[0] * first[1]
+        })
+        .sum::<f64>()
+        * 0.5
+}
+
+fn signed_area_3d(ring: &[[f32; 3]]) -> f64 {
+    ring_area(&ring.iter().copied().map(xy).collect::<Vec<_>>())
+}
+
+fn xy(position: [f32; 3]) -> [f64; 2] {
+    [f64::from(position[0]), f64::from(position[1])]
+}
+
+fn same_xy(first: [f32; 3], second: [f32; 3]) -> bool {
+    distance_squared(xy(first), xy(second)) <= POINT_EPSILON_SQUARED
+}
+
+fn cross2(first: [f64; 2], second: [f64; 2], third: [f64; 2]) -> f64 {
+    (second[0] - first[0]) * (third[1] - first[1]) - (second[1] - first[1]) * (third[0] - first[0])
+}
+
+fn distance_squared(first: [f64; 2], second: [f64; 2]) -> f64 {
+    let dx = first[0] - second[0];
+    let dy = first[1] - second[1];
+    dx * dx + dy * dy
+}
+
+pub(crate) fn basement_surface_mask(positions: [[f32; 3]; 3]) -> FullSurfaceMask {
+    let first = positions[1];
+    let origin = positions[0];
+    let second = positions[2];
+    let normal_z = (first[0] - origin[0]) * (second[1] - origin[1])
+        - (first[1] - origin[1]) * (second[0] - origin[0]);
+    FullSurfaceMask(if normal_z < 0.0 { 0x0300 } else { 0x8300 })
+}
+
+fn legalize_constrained_edges(
+    mut triangles: Vec<[[f32; 3]; 3]>,
+    outer: &[[f32; 3]],
+    inner: &[[f32; 3]],
+) -> Result<Vec<[[f32; 3]; 3]>, BasementFoundationError> {
+    let max_flips = triangles
+        .len()
+        .saturating_mul(triangles.len())
+        .saturating_mul(4)
+        .max(1);
+    for _ in 0..max_flips {
+        let mut flipped = false;
+        'pairs: for first in 0..triangles.len() {
+            for second in (first + 1)..triangles.len() {
+                let Some((a, b, c, d)) = shared_edge(triangles[first], triangles[second]) else {
+                    continue;
+                };
+                if ring_has_edge(outer, a, b) || ring_has_edge(inner, a, b) {
+                    continue;
+                }
+                if !convex_quadrilateral(a, b, c, d) || in_circle(a, b, c, d) <= 1.0e-8 {
+                    continue;
+                }
+                let first_triangle = orient_ccw([c, d, b]);
+                let second_triangle = orient_ccw([d, c, a]);
+                if cross2(
+                    xy(first_triangle[0]),
+                    xy(first_triangle[1]),
+                    xy(first_triangle[2]),
+                ) <= 1.0e-9
+                    || cross2(
+                        xy(second_triangle[0]),
+                        xy(second_triangle[1]),
+                        xy(second_triangle[2]),
+                    ) <= 1.0e-9
+                {
+                    return Err(BasementFoundationError::InvalidContour);
+                }
+                triangles[first] = first_triangle;
+                triangles[second] = second_triangle;
+                flipped = true;
+                break 'pairs;
+            }
+        }
+        if !flipped {
+            return Ok(triangles);
+        }
+    }
+    Err(BasementFoundationError::InvalidContour)
+}
+
+fn shared_edge(
+    first: [[f32; 3]; 3],
+    second: [[f32; 3]; 3],
+) -> Option<([f32; 3], [f32; 3], [f32; 3], [f32; 3])> {
+    for edge in 0..3 {
+        let a = first[edge];
+        let b = first[(edge + 1) % 3];
+        if let Some(opposite) = second
+            .iter()
+            .copied()
+            .find(|point| !same_xy(*point, a) && !same_xy(*point, b))
+        {
+            let has_a = second.iter().any(|point| same_xy(*point, a));
+            let has_b = second.iter().any(|point| same_xy(*point, b));
+            if has_a && has_b {
+                let first_opposite = first[(edge + 2) % 3];
+                return Some((a, b, first_opposite, opposite));
+            }
+        }
+    }
+    None
+}
+
+fn ring_has_edge(ring: &[[f32; 3]], a: [f32; 3], b: [f32; 3]) -> bool {
+    (0..ring.len()).any(|index| {
+        let first = ring[index];
+        let second = ring[(index + 1) % ring.len()];
+        (same_xy(first, a) && same_xy(second, b)) || (same_xy(first, b) && same_xy(second, a))
+    })
+}
+
+fn convex_quadrilateral(a: [f32; 3], b: [f32; 3], c: [f32; 3], d: [f32; 3]) -> bool {
+    let epsilon = 1.0e-9;
+    let c_side = cross2(xy(a), xy(b), xy(c));
+    let d_side = cross2(xy(a), xy(b), xy(d));
+    let a_side = cross2(xy(c), xy(d), xy(a));
+    let b_side = cross2(xy(c), xy(d), xy(b));
+    c_side * d_side < -epsilon && a_side * b_side < -epsilon
+}
+
+fn in_circle(a: [f32; 3], b: [f32; 3], c: [f32; 3], d: [f32; 3]) -> f64 {
+    let a = [
+        f64::from(a[0]) - f64::from(d[0]),
+        f64::from(a[1]) - f64::from(d[1]),
+    ];
+    let b = [
+        f64::from(b[0]) - f64::from(d[0]),
+        f64::from(b[1]) - f64::from(d[1]),
+    ];
+    let c = [
+        f64::from(c[0]) - f64::from(d[0]),
+        f64::from(c[1]) - f64::from(d[1]),
+    ];
+    (a[0] * a[0] + a[1] * a[1]) * (b[0] * c[1] - b[1] * c[0])
+        - (b[0] * b[0] + b[1] * b[1]) * (a[0] * c[1] - a[1] * c[0])
+        + (c[0] * c[0] + c[1] * c[1]) * (a[0] * b[1] - a[1] * b[0])
+}
+
+fn orient_ccw(mut triangle: [[f32; 3]; 3]) -> [[f32; 3]; 3] {
+    if cross2(xy(triangle[0]), xy(triangle[1]), xy(triangle[2])) < 0.0 {
+        triangle.swap(1, 2);
+    }
+    triangle
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use fparkan_terrain_format::{TerrainFace28, TerrainSlotTable};
+
+    #[test]
+    fn foundation_samples_lod0_planes_and_flattens_or_preserves_authored_inner_z() {
+        let mesh = two_lod_mesh();
+        let outer = [[0.4, 0.0], [4.0, 0.0], [4.0, 4.0], [0.0, 4.0], [0.4, 0.0]];
+        let inner = [
+            [1.0, 1.0, 10.0],
+            [3.0, 1.0, 11.0],
+            [3.0, 3.0, 12.0],
+            [1.0, 3.0, 13.0],
+            [1.0, 1.0, 10.0],
+        ];
+        let (open_outer, _) = normalize_ring(&outer).expect("outer ring");
+        let snapped = normalize_outer_points(&open_outer, &[0, 1, 2], &mesh)
+            .expect("outer points lie within LOD0");
+        assert_eq!(snapped.first(), Some(&[0.0, 0.0]));
+        assert!(!snapped.contains(&[2.0, 0.0]));
+        assert!(native_segment_inserts_vertex(
+            [0.0, 0.0],
+            [4.0, 0.0],
+            [3.96875, 0.4991, 12.0]
+        ));
+        assert!(!native_segment_inserts_vertex(
+            [0.0, 0.0],
+            [4.0, 0.0],
+            [2.0, 0.0, 12.0]
+        ));
+        let terrain = [
+            (
+                0,
+                mesh.faces[0].flags,
+                terrain_face_positions(&mesh, 0).unwrap(),
+            ),
+            (
+                1,
+                mesh.faces[1].flags,
+                terrain_face_positions(&mesh, 1).unwrap(),
+            ),
+            (
+                2,
+                mesh.faces[2].flags,
+                terrain_face_positions(&mesh, 2).unwrap(),
+            ),
+        ];
+        let crossings = contour_crossings(&snapped, &terrain).expect("LOD0 crossings");
+        assert!(crossings.len() > 4);
+        assert_eq!(crossings.first().unwrap()[2], 0.0);
+        assert!(crossings
+            .iter()
+            .any(|crossing| same_xy(*crossing, [2.0, 0.0, 2.0])));
+        let expected_mean =
+            crossings.iter().map(|crossing| crossing[2]).sum::<f32>() / crossings.len() as f32;
+
+        let flat = build_basement_foundations(
+            &mesh,
+            &[BasementFootprint {
+                outer: &outer,
+                inner: &inner,
+                flatten_inner_to_mean_outer_z: true,
+            }],
+        )
+        .expect("flat foundation");
+        let placement = flat.placements[0];
+        assert!((placement.placement_delta - (10.0 - expected_mean)).abs() < 1.0e-6);
+        let flat_band = flat
+            .fragments
+            .iter()
+            .filter(|fragment| fragment.source_face.is_none())
+            .collect::<Vec<_>>();
+        assert!(flat
+            .fragments
+            .iter()
+            .filter_map(|fragment| fragment.source_face)
+            .all(|face| face < 3));
+        assert!(!flat_band.is_empty());
+        assert!(flat_band
+            .iter()
+            .all(|fragment| fragment.surface_mask == FullSurfaceMask(0x8300)));
+        for point in &inner[..4] {
+            assert!(flat_band
+                .iter()
+                .flat_map(|fragment| fragment.positions)
+                .any(|vertex| {
+                    same_xy(vertex, *point) && (vertex[2] - expected_mean).abs() < 1.0e-6
+                }));
+        }
+        let flat_area = flat_band
+            .iter()
+            .map(|fragment| {
+                let [a, b, c] = fragment.positions;
+                cross2(xy(a), xy(b), xy(c)).abs() * 0.5
+            })
+            .sum::<f64>();
+        assert!((flat_area - 12.0).abs() < 1.0e-5);
+
+        let authored = build_basement_foundations(
+            &mesh,
+            &[BasementFootprint {
+                outer: &outer,
+                inner: &inner,
+                flatten_inner_to_mean_outer_z: false,
+            }],
+        )
+        .expect("authored-height foundation");
+        let authored_band = authored
+            .fragments
+            .iter()
+            .filter(|fragment| fragment.source_face.is_none())
+            .collect::<Vec<_>>();
+        for point in &inner[..4] {
+            assert!(authored_band
+                .iter()
+                .flat_map(|fragment| fragment.positions)
+                .any(|vertex| {
+                    same_xy(vertex, *point) && (vertex[2] - point[2]).abs() < 1.0e-6
+                }));
+        }
+    }
+
+    #[test]
+    fn foundation_rejects_outer_contours_outside_or_across_uncovered_terrain() {
+        let mesh = two_lod_mesh();
+        let outside = [[5.0, 1.0], [5.5, 1.0], [5.0, 1.5]];
+        let inner = [[5.1, 1.1, 0.0], [5.3, 1.1, 0.0], [5.1, 1.3, 0.0]];
+        assert_eq!(
+            build_basement_foundations(
+                &mesh,
+                &[BasementFootprint {
+                    outer: &outside,
+                    inner: &inner,
+                    flatten_inner_to_mean_outer_z: true,
+                }],
+            ),
+            Err(BasementFoundationError::OuterPointOutsideLod0 { point: 0 })
+        );
+
+        let mut islands = two_lod_mesh();
+        islands.positions = vec![
+            [0.0, 0.0, 0.0],
+            [1.0, 0.0, 0.0],
+            [0.0, 1.0, 0.0],
+            [2.0, 2.0, 0.0],
+            [3.0, 2.0, 0.0],
+            [2.0, 3.0, 0.0],
+        ];
+        islands.faces.truncate(2);
+        islands.faces[0].vertices = [0, 1, 2];
+        islands.faces[1].vertices = [3, 4, 5];
+        islands.slots.slots_raw[0][2..4].copy_from_slice(&2_u16.to_le_bytes());
+        islands.slots.slots_raw[1][0..2].copy_from_slice(&2_u16.to_le_bytes());
+        islands.slots.slots_raw[1][2..4].copy_from_slice(&0_u16.to_le_bytes());
+        let outer = [[0.1, 0.1], [0.7, 0.1], [2.1, 2.1]];
+        let inner = [[0.2, 0.2, 0.0], [0.4, 0.2, 0.0], [0.2, 0.4, 0.0]];
+        assert_eq!(
+            build_basement_foundations(
+                &islands,
+                &[BasementFootprint {
+                    outer: &outer,
+                    inner: &inner,
+                    flatten_inner_to_mean_outer_z: true,
+                }],
+            ),
+            Err(BasementFoundationError::UncoveredOuterSegment { segment: 1 })
+        );
+    }
+
+    #[test]
+    fn shared_diagonal_crossing_is_exact_and_counted_once() {
+        let first = [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [1.0, 1.0, 1.0]];
+        let second = [[0.0, 0.0, 0.0], [1.0, 1.0, 1.0], [0.0, 1.0, 2.0]];
+        let edge_start = [0.8, 0.5];
+        let edge_end = [0.5, 0.8];
+        let first_interval =
+            segment_triangle_interval(edge_start, edge_end, first).expect("first face interval");
+        let second_interval =
+            segment_triangle_interval(edge_start, edge_end, second).expect("second face interval");
+        assert!((first_interval.1 - 0.5).abs() < 1.0e-9);
+        assert!((second_interval.0 - 0.5).abs() < 1.0e-9);
+
+        let ring = [[0.5, 0.2], [0.8, 0.5], [0.5, 0.8], [0.2, 0.5]];
+        let crossings = contour_crossings(
+            &ring,
+            &[
+                (0, FullSurfaceMask(1), first),
+                (1, FullSurfaceMask(2), second),
+            ],
+        )
+        .expect("shared-edge crossings");
+        let shared = crossings
+            .iter()
+            .filter(|crossing| {
+                distance_squared(
+                    [f64::from(crossing[0]), f64::from(crossing[1])],
+                    [0.65, 0.65],
+                ) <= POINT_EPSILON_SQUARED
+            })
+            .collect::<Vec<_>>();
+        assert_eq!(
+            shared.len(),
+            1,
+            "shared-edge crossing duplicated: {shared:?}"
+        );
+        assert!((shared[0][2] - 0.65).abs() < 1.0e-6);
+    }
+
+    #[test]
+    fn annulus_preserves_collinear_outer_height_breakpoints_and_is_constrained_delaunay() {
+        let outer = [
+            [0.0, 0.0, 0.0],
+            [2.0, 0.0, 3.0],
+            [4.0, 0.0, 10.0],
+            [4.0, 4.0, 12.0],
+            [0.0, 4.0, 0.0],
+        ];
+        let inner = vec![
+            [1.0, 1.0, 5.0],
+            [3.0, 1.0, 5.0],
+            [3.0, 3.0, 5.0],
+            [1.0, 3.0, 5.0],
+        ];
+        let triangles = triangulate_annulus(&outer, inner.clone()).expect("annulus");
+        for edge in outer.windows(2) {
+            assert!(
+                triangles
+                    .iter()
+                    .any(|triangle| has_xy_edge(*triangle, edge[0], edge[1])),
+                "missing outer boundary segment {:?} -> {:?}",
+                edge[0],
+                edge[1]
+            );
+        }
+        for point in outer {
+            assert!(
+                triangles
+                    .iter()
+                    .flat_map(|triangle| triangle.iter())
+                    .any(|vertex| {
+                        same_xy(*vertex, point) && (vertex[2] - point[2]).abs() < 1.0e-6
+                    }),
+                "lost outer breakpoint {point:?}"
+            );
+        }
+        assert_locally_delaunay(&triangles, &outer, &inner);
+    }
+
+    #[test]
+    fn constrained_annulus_handles_concave_rings_with_different_vertex_counts() {
+        let outer = [
+            [0.0, 0.0, 1.0],
+            [10.0, 0.0, 2.0],
+            [10.0, 3.0, 3.0],
+            [4.0, 3.0, 4.0],
+            [4.0, 10.0, 5.0],
+            [0.0, 10.0, 6.0],
+        ];
+        let inner = vec![
+            [1.0, 1.0, 3.0],
+            [3.0, 1.0, 3.0],
+            [3.0, 1.8, 3.0],
+            [2.2, 1.8, 3.0],
+            [2.2, 2.4, 3.0],
+            [3.0, 2.4, 3.0],
+            [3.0, 3.0, 3.0],
+            [1.0, 3.0, 3.0],
+        ];
+        let band = triangulate_annulus(&outer, inner.clone()).expect("constrained annulus");
+        let area = band
+            .iter()
+            .map(|triangle| cross2(xy(triangle[0]), xy(triangle[1]), xy(triangle[2])).abs() * 0.5)
+            .sum::<f64>();
+        let expected = ring_area(&outer.iter().copied().map(xy).collect::<Vec<_>>()).abs()
+            - ring_area(&inner.iter().copied().map(xy).collect::<Vec<_>>()).abs();
+        assert!(
+            (area - expected).abs() < 1.0e-5,
+            "area {area}, expected {expected}"
+        );
+
+        let outer_xy = outer.iter().copied().map(xy).collect::<Vec<_>>();
+        let inner_xy = inner.iter().copied().map(xy).collect::<Vec<_>>();
+        for triangle in &band {
+            let centroid = triangle_centroid(*triangle);
+            assert!(point_in_ring(centroid, &outer_xy));
+            assert!(!point_in_ring(centroid, &inner_xy));
+        }
+        for first in 0..band.len() {
+            for second in (first + 1)..band.len() {
+                assert!(!positive_area_overlap(band[first], band[second]));
+            }
+        }
+        assert_locally_delaunay(&band, &outer, &inner);
+    }
+
+    #[test]
+    fn multiple_disjoint_footprints_clip_canonical_outside_once() {
+        let mesh = two_lod_mesh();
+        let first_outer = [[0.6, 0.6], [1.7, 0.6], [1.7, 1.7], [0.6, 1.7]];
+        let first_inner = [
+            [0.9, 0.9, 5.0],
+            [1.4, 0.9, 5.0],
+            [1.4, 1.4, 5.0],
+            [0.9, 1.4, 5.0],
+        ];
+        let second_outer = [[2.3, 2.3], [3.4, 2.3], [3.4, 3.4], [2.3, 3.4]];
+        let second_inner = [
+            [2.6, 2.6, 7.0],
+            [3.1, 2.6, 7.0],
+            [3.1, 3.1, 7.0],
+            [2.6, 3.1, 7.0],
+        ];
+        let output = build_basement_foundations(
+            &mesh,
+            &[
+                BasementFootprint {
+                    outer: &first_outer,
+                    inner: &first_inner,
+                    flatten_inner_to_mean_outer_z: true,
+                },
+                BasementFootprint {
+                    outer: &second_outer,
+                    inner: &second_inner,
+                    flatten_inner_to_mean_outer_z: true,
+                },
+            ],
+        )
+        .expect("two disjoint foundations");
+        assert_eq!(output.placements.len(), 2);
+        let terrain = output
+            .fragments
+            .iter()
+            .filter(|fragment| fragment.source_face.is_some())
+            .collect::<Vec<_>>();
+        let bands = output
+            .fragments
+            .iter()
+            .filter(|fragment| fragment.source_face.is_none())
+            .collect::<Vec<_>>();
+        assert!(!terrain.is_empty());
+        assert!(!bands.is_empty());
+        for fragment in &terrain {
+            let centroid = triangle_centroid(fragment.positions);
+            assert!(!point_strictly_in_ring(centroid, &ring_f64(&first_outer)));
+            assert!(!point_strictly_in_ring(centroid, &ring_f64(&second_outer)));
+        }
+        for first in 0..terrain.len() {
+            for second in (first + 1)..terrain.len() {
+                assert!(!positive_area_overlap(
+                    terrain[first].positions,
+                    terrain[second].positions
+                ));
+            }
+        }
+        for band in bands {
+            let centroid = triangle_centroid(band.positions);
+            assert!(
+                point_in_ring(centroid, &ring_f64(&first_outer))
+                    || point_in_ring(centroid, &ring_f64(&second_outer))
+            );
+            assert!(!point_in_ring(
+                centroid,
+                &ring_f64(&first_inner.iter().map(|p| [p[0], p[1]]).collect::<Vec<_>>())
+            ));
+            assert!(!point_in_ring(
+                centroid,
+                &ring_f64(
+                    &second_inner
+                        .iter()
+                        .map(|p| [p[0], p[1]])
+                        .collect::<Vec<_>>()
+                )
+            ));
+        }
+    }
+
+    #[test]
+    fn overlapping_footprints_are_explicitly_rejected() {
+        let mesh = two_lod_mesh();
+        let outer_a = [[0.5, 0.5], [2.5, 0.5], [2.5, 2.5], [0.5, 2.5]];
+        let inner_a = [
+            [1.0, 1.0, 3.0],
+            [2.0, 1.0, 3.0],
+            [2.0, 2.0, 3.0],
+            [1.0, 2.0, 3.0],
+        ];
+        let outer_b = [[1.5, 1.5], [3.5, 1.5], [3.5, 3.5], [1.5, 3.5]];
+        let inner_b = [
+            [2.0, 2.0, 3.0],
+            [3.0, 2.0, 3.0],
+            [3.0, 3.0, 3.0],
+            [2.0, 3.0, 3.0],
+        ];
+        assert_eq!(
+            build_basement_foundations(
+                &mesh,
+                &[
+                    BasementFootprint {
+                        outer: &outer_a,
+                        inner: &inner_a,
+                        flatten_inner_to_mean_outer_z: true
+                    },
+                    BasementFootprint {
+                        outer: &outer_b,
+                        inner: &inner_b,
+                        flatten_inner_to_mean_outer_z: true
+                    },
+                ],
+            ),
+            Err(BasementFoundationError::OverlappingFootprints)
+        );
+    }
+
+    #[test]
+    fn generated_flags_follow_native_normal_z_branch() {
+        let positive = [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]];
+        let negative = [positive[0], positive[2], positive[1]];
+        assert_eq!(basement_surface_mask(positive), FullSurfaceMask(0x8300));
+        assert_eq!(basement_surface_mask(negative), FullSurfaceMask(0x0300));
+    }
+
+    #[test]
+    fn foundation_rejects_inner_ring_that_is_not_nested() {
+        let mesh = two_lod_mesh();
+        let outer = [[0.0, 0.0], [4.0, 0.0], [4.0, 4.0], [0.0, 4.0]];
+        let inner = [
+            [3.0, 1.0, 0.0],
+            [5.0, 1.0, 0.0],
+            [5.0, 3.0, 0.0],
+            [3.0, 3.0, 0.0],
+        ];
+        assert_eq!(
+            build_basement_foundations(
+                &mesh,
+                &[BasementFootprint {
+                    outer: &outer,
+                    inner: &inner,
+                    flatten_inner_to_mean_outer_z: true
+                }],
+            ),
+            Err(BasementFoundationError::InvalidContour)
+        );
+    }
+
+    fn assert_locally_delaunay(
+        triangles: &[[[f32; 3]; 3]],
+        outer: &[[f32; 3]],
+        inner: &[[f32; 3]],
+    ) {
+        for first in 0..triangles.len() {
+            for second in (first + 1)..triangles.len() {
+                let Some((a, b, c, d)) = shared_edge(triangles[first], triangles[second]) else {
+                    continue;
+                };
+                if ring_has_edge(outer, a, b)
+                    || ring_has_edge(inner, a, b)
+                    || !convex_quadrilateral(a, b, c, d)
+                {
+                    continue;
+                }
+                assert!(
+                    in_circle(a, b, c, d) <= 1.0e-8,
+                    "illegal diagonal {a:?}->{b:?}"
+                );
+            }
+        }
+    }
+
+    fn has_xy_edge(triangle: [[f32; 3]; 3], first: [f32; 3], second: [f32; 3]) -> bool {
+        (0..3).any(|index| {
+            (same_xy(triangle[index], first) && same_xy(triangle[(index + 1) % 3], second))
+                || (same_xy(triangle[index], second) && same_xy(triangle[(index + 1) % 3], first))
+        })
+    }
+
+    fn triangle_centroid(triangle: [[f32; 3]; 3]) -> [f64; 2] {
+        [
+            triangle
+                .iter()
+                .map(|point| f64::from(point[0]))
+                .sum::<f64>()
+                / 3.0,
+            triangle
+                .iter()
+                .map(|point| f64::from(point[1]))
+                .sum::<f64>()
+                / 3.0,
+        ]
+    }
+
+    fn point_strictly_in_ring(point: [f64; 2], ring: &[[f64; 2]]) -> bool {
+        !((0..ring.len())
+            .any(|index| point_on_segment(point, ring[index], ring[(index + 1) % ring.len()])))
+            && point_in_ring(point, ring)
+    }
+
+    fn ring_f64(ring: &[[f32; 2]]) -> Vec<[f64; 2]> {
+        ring.iter()
+            .map(|point| [f64::from(point[0]), f64::from(point[1])])
+            .collect()
+    }
+
+    fn two_lod_mesh() -> LandMeshDocument {
+        let positions = vec![
+            [0.0, 0.0, 0.0],
+            [2.0, 0.0, 2.0],
+            [4.0, 0.0, 4.0],
+            [0.0, 4.0, 4.0],
+            [4.0, 4.0, 8.0],
+            [0.0, 0.0, 100.0],
+            [4.0, 0.0, 100.0],
+            [0.0, 4.0, 100.0],
+        ];
+        let face = |flags, vertices| TerrainFace28 {
+            flags: FullSurfaceMask(flags),
+            material_tag: 0,
+            aux_tag: 0,
+            vertices,
+            neighbors: [None; 3],
+            tail_raw: [0; 8],
+            raw: [0; 28],
+        };
+        let faces = vec![
+            face(1, [0, 1, 3]),
+            face(2, [1, 4, 3]),
+            face(4, [1, 2, 4]),
+            face(8, [5, 6, 7]),
+        ];
+        let mut nodes_raw = vec![u8::MAX; 38];
+        nodes_raw[8..10].copy_from_slice(&0_u16.to_le_bytes());
+        nodes_raw[10..12].copy_from_slice(&1_u16.to_le_bytes());
+        let mut slots_raw = vec![[0_u8; 68]; 2];
+        slots_raw[0][2..4].copy_from_slice(&3_u16.to_le_bytes());
+        slots_raw[1][0..2].copy_from_slice(&3_u16.to_le_bytes());
+        slots_raw[1][2..4].copy_from_slice(&1_u16.to_le_bytes());
+        LandMeshDocument {
+            streams: Vec::new(),
+            nodes_raw,
+            slots: TerrainSlotTable {
+                header_raw: Vec::new(),
+                slots_raw,
+            },
+            positions,
+            normals: Vec::new(),
+            uv0: Vec::new(),
+            accelerator: Vec::new(),
+            aux14: Vec::new(),
+            aux18: Vec::new(),
+            faces,
+        }
+    }
+
+    fn positive_area_overlap(first: [[f32; 3]; 3], second: [[f32; 3]; 3]) -> bool {
+        for triangle in [first, second] {
+            for edge in 0..3 {
+                let a = xy(triangle[edge]);
+                let b = xy(triangle[(edge + 1) % 3]);
+                let axis = [b[1] - a[1], a[0] - b[0]];
+                let project = |vertices: [[f32; 3]; 3]| {
+                    vertices.map(|vertex| {
+                        let point = xy(vertex);
+                        point[0] * axis[0] + point[1] * axis[1]
+                    })
+                };
+                let first_projection = project(first);
+                let second_projection = project(second);
+                let first_min = first_projection
+                    .iter()
+                    .copied()
+                    .fold(f64::INFINITY, f64::min);
+                let first_max = first_projection
+                    .iter()
+                    .copied()
+                    .fold(f64::NEG_INFINITY, f64::max);
+                let second_min = second_projection
+                    .iter()
+                    .copied()
+                    .fold(f64::INFINITY, f64::min);
+                let second_max = second_projection
+                    .iter()
+                    .copied()
+                    .fold(f64::NEG_INFINITY, f64::max);
+                if first_max <= second_min + 1.0e-8 || second_max <= first_min + 1.0e-8 {
+                    return false;
+                }
+            }
+        }
+        true
+    }
+}
diff --git a/crates/fparkan-terrain/src/lib.rs b/crates/fparkan-terrain/src/lib.rs
index a014b61..db55506 100644
--- a/crates/fparkan-terrain/src/lib.rs
+++ b/crates/fparkan-terrain/src/lib.rs
@@ -4,6 +4,15 @@
 use fparkan_terrain_format::{LandMapDocument, LandMeshDocument};
 use std::collections::VecDeque;
 
+mod contour;
+mod foundation;
+
+use foundation::basement_surface_mask;
+pub use foundation::{
+    build_basement_foundations, BasementFootprint, BasementFoundation, BasementFoundationError,
+    BasementFoundationFragment, BasementPlacement,
+};
+
 /// Terrain material-selector contract preserved from the decoded map mesh.
 ///
 /// Applications access this semantic terrain API through this terrain crate
@@ -21,6 +30,7 @@ pub struct TerrainWorld {
     surfaces: Vec<RuntimeTriangle>,
     surface_bvh: SurfaceBvh,
     source_mesh: Option<LandMeshDocument>,
+    basement_foundation: Option<BasementFoundation>,
 }
 
 /// Surface hit.
@@ -91,6 +101,13 @@ pub enum TerrainError {
         /// Position count.
         position_count: usize,
     },
+    /// A basement fragment does not preserve a valid LOD0 source face and mask.
+    InvalidFoundationFragment {
+        /// Fragment index.
+        fragment: usize,
+    },
+    /// Foundation fragments require a retained source Land.msh document.
+    FoundationRequiresMesh,
     /// Terrain Node38/Slot68 records do not assign every face to one LOD.
     InvalidLodMapping(TerrainLodError),
 }
@@ -119,6 +136,12 @@ impl std::fmt::Display for TerrainError {
                 f,
                 "terrain face {face} references vertex {vertex} outside {position_count} positions"
             ),
+            Self::InvalidFoundationFragment { fragment } => {
+                write!(f, "invalid basement foundation fragment {fragment}")
+            }
+            Self::FoundationRequiresMesh => {
+                write!(f, "basement foundation requires a source Land.msh")
+            }
             Self::InvalidLodMapping(error) => write!(f, "invalid terrain LOD mapping: {error}"),
         }
     }
@@ -132,7 +155,9 @@ impl std::error::Error for TerrainError {
             | Self::TooManyAreals { .. }
             | Self::InvalidGridReference { .. }
             | Self::InvalidArealReference { .. }
-            | Self::InvalidSurfaceVertex { .. } => None,
+            | Self::InvalidSurfaceVertex { .. }
+            | Self::InvalidFoundationFragment { .. }
+            | Self::FoundationRequiresMesh => None,
         }
     }
 }
@@ -254,6 +279,7 @@ impl TerrainWorld {
             surfaces: Vec::new(),
             surface_bvh: SurfaceBvh::default(),
             source_mesh: None,
+            basement_foundation: None,
         })
     }
 
@@ -290,6 +316,87 @@ impl TerrainWorld {
         Ok(world)
     }
 
+    /// Returns a new world whose surface queries use exactly the canonical
+    /// source and generated triangles from a basement foundation.
+    ///
+    /// # Errors
+    ///
+    /// Returns [`TerrainError`] if this world has no source mesh or any
+    /// fragment does not preserve valid LOD0 identity, flags, and geometry.
+    pub fn with_basement_foundation(
+        mut self,
+        foundation: BasementFoundation,
+    ) -> Result<Self, TerrainError> {
+        let mesh = self
+            .source_mesh
+            .as_ref()
+            .ok_or(TerrainError::FoundationRequiresMesh)?;
+        let assignments = mesh
+            .face_lod_assignments()
+            .map_err(TerrainError::InvalidLodMapping)?;
+        let is_lod0 = assignments
+            .iter()
+            .map(|assignment| assignment.lod == TerrainLod::Lod0)
+            .collect::<Vec<_>>();
+        let mut surfaces = Vec::with_capacity(foundation.fragments.len());
+        let mut generated = 0_usize;
+        for (fragment_index, fragment) in foundation.fragments.iter().enumerate() {
+            if !valid_foundation_triangle(fragment.positions) {
+                return Err(TerrainError::InvalidFoundationFragment {
+                    fragment: fragment_index,
+                });
+            }
+            let face = if let Some(source_face) = fragment.source_face {
+                let Some(source) = mesh.faces.get(source_face) else {
+                    return Err(TerrainError::InvalidFoundationFragment {
+                        fragment: fragment_index,
+                    });
+                };
+                if !is_lod0.get(source_face).copied().unwrap_or(false)
+                    || source.flags != fragment.surface_mask
+                {
+                    return Err(TerrainError::InvalidFoundationFragment {
+                        fragment: fragment_index,
+                    });
+                }
+                source_face
+            } else {
+                if basement_surface_mask(fragment.positions) != fragment.surface_mask {
+                    return Err(TerrainError::InvalidFoundationFragment {
+                        fragment: fragment_index,
+                    });
+                }
+                let face = mesh.faces.len().checked_add(generated).ok_or(
+                    TerrainError::InvalidFoundationFragment {
+                        fragment: fragment_index,
+                    },
+                )?;
+                generated =
+                    generated
+                        .checked_add(1)
+                        .ok_or(TerrainError::InvalidFoundationFragment {
+                            fragment: fragment_index,
+                        })?;
+                face
+            };
+            surfaces.push(RuntimeTriangle {
+                face,
+                mask: fragment.surface_mask,
+                vertices: fragment.positions,
+            });
+        }
+        self.surfaces = surfaces;
+        self.surface_bvh = SurfaceBvh::from_surfaces(&self.surfaces);
+        self.basement_foundation = Some(foundation);
+        Ok(self)
+    }
+
+    /// Returns the canonical basement geometry used by the surface query index.
+    #[must_use]
+    pub fn basement_foundation(&self) -> Option<&BasementFoundation> {
+        self.basement_foundation.as_ref()
+    }
+
     /// Returns the number of navigation areas.
     #[must_use]
     pub fn areal_count(&self) -> usize {
@@ -740,26 +847,31 @@ impl RuntimeTriangle {
     }
 
     fn raycast(&self, origin: [f32; 3], direction: [f32; 3]) -> Option<f32> {
-        let edge1 = sub3(self.vertices[1], self.vertices[0]);
-        let edge2 = sub3(self.vertices[2], self.vertices[0]);
-        let pvec = cross3(direction, edge2);
-        let det = dot3(edge1, pvec);
-        if det.abs() <= 1.0e-6 {
+        let edge1 = sub3_f64(self.vertices[1], self.vertices[0]);
+        let edge2 = sub3_f64(self.vertices[2], self.vertices[0]);
+        let direction = direction.map(f64::from);
+        let pvec = cross3_f64(direction, edge2);
+        let det = dot3_f64(edge1, pvec);
+        if det == 0.0 || !det.is_finite() {
             return None;
         }
         let inv_det = 1.0 / det;
-        let tvec = sub3(origin, self.vertices[0]);
-        let u = dot3(tvec, pvec) * inv_det;
+        let tvec = sub3_f64(origin, self.vertices[0]);
+        let u = dot3_f64(tvec, pvec) * inv_det;
         if !(-1.0e-5..=1.00001).contains(&u) {
             return None;
         }
-        let qvec = cross3(tvec, edge1);
-        let v = dot3(direction, qvec) * inv_det;
+        let qvec = cross3_f64(tvec, edge1);
+        let v = dot3_f64(direction, qvec) * inv_det;
         if v < -1.0e-5 || u + v > 1.00001 {
             return None;
         }
-        let distance = dot3(edge2, qvec) * inv_det;
-        (distance >= 0.0).then_some(distance)
+        let distance = dot3_f64(edge2, qvec) * inv_det;
+        if distance >= 0.0 && distance.is_finite() && distance <= f64::from(f32::MAX) {
+            Some(distance as f32)
+        } else {
+            None
+        }
     }
 }
 
@@ -919,6 +1031,21 @@ fn build_surfaces(mesh: &LandMeshDocument) -> Result<Vec<RuntimeTriangle>, Terra
     Ok(triangles)
 }
 
+fn valid_foundation_triangle(vertices: [[f32; 3]; 3]) -> bool {
+    if !vertices.iter().flatten().all(|value| value.is_finite()) {
+        return false;
+    }
+    let first = [
+        f64::from(vertices[1][0]) - f64::from(vertices[0][0]),
+        f64::from(vertices[1][1]) - f64::from(vertices[0][1]),
+    ];
+    let second = [
+        f64::from(vertices[2][0]) - f64::from(vertices[0][0]),
+        f64::from(vertices[2][1]) - f64::from(vertices[0][1]),
+    ];
+    first[0] * second[1] - first[1] * second[0] != 0.0
+}
+
 fn surface_vertex(
     mesh: &LandMeshDocument,
     face: usize,
@@ -940,25 +1067,41 @@ fn barycentric_2d(
     second: [f32; 2],
     third: [f32; 2],
 ) -> Option<[f32; 3]> {
-    let edge_second = [second[0] - first[0], second[1] - first[1]];
-    let edge_third = [third[0] - first[0], third[1] - first[1]];
-    let point_delta = [point[0] - first[0], point[1] - first[1]];
+    let edge_second = [
+        f64::from(second[0]) - f64::from(first[0]),
+        f64::from(second[1]) - f64::from(first[1]),
+    ];
+    let edge_third = [
+        f64::from(third[0]) - f64::from(first[0]),
+        f64::from(third[1]) - f64::from(first[1]),
+    ];
+    let point_delta = [
+        f64::from(point[0]) - f64::from(first[0]),
+        f64::from(point[1]) - f64::from(first[1]),
+    ];
     let denom = edge_second[0] * edge_third[1] - edge_third[0] * edge_second[1];
-    if denom.abs() <= 1.0e-6 {
+    if denom == 0.0 {
         return None;
     }
-    let inv = 1.0 / denom;
-    let second_weight = (point_delta[0] * edge_third[1] - edge_third[0] * point_delta[1]) * inv;
-    let third_weight = (edge_second[0] * point_delta[1] - point_delta[0] * edge_second[1]) * inv;
+    let second_weight = (point_delta[0] * edge_third[1] - edge_third[0] * point_delta[1]) / denom;
+    let third_weight = (edge_second[0] * point_delta[1] - point_delta[0] * edge_second[1]) / denom;
     let first_weight = 1.0 - second_weight - third_weight;
-    Some([first_weight, second_weight, third_weight])
+    let weights = [first_weight, second_weight, third_weight];
+    if weights.iter().any(|weight| !weight.is_finite()) {
+        return None;
+    }
+    Some(weights.map(|weight| weight as f32))
 }
 
-fn sub3(a: [f32; 3], b: [f32; 3]) -> [f32; 3] {
-    [a[0] - b[0], a[1] - b[1], a[2] - b[2]]
+fn sub3_f64(a: [f32; 3], b: [f32; 3]) -> [f64; 3] {
+    [
+        f64::from(a[0]) - f64::from(b[0]),
+        f64::from(a[1]) - f64::from(b[1]),
+        f64::from(a[2]) - f64::from(b[2]),
+    ]
 }
 
-fn cross3(a: [f32; 3], b: [f32; 3]) -> [f32; 3] {
+fn cross3_f64(a: [f64; 3], b: [f64; 3]) -> [f64; 3] {
     [
         a[1] * b[2] - a[2] * b[1],
         a[2] * b[0] - a[0] * b[2],
@@ -966,7 +1109,7 @@ fn cross3(a: [f32; 3], b: [f32; 3]) -> [f32; 3] {
     ]
 }
 
-fn dot3(a: [f32; 3], b: [f32; 3]) -> f32 {
+fn dot3_f64(a: [f64; 3], b: [f64; 3]) -> f64 {
     a[0] * b[0] + a[1] * b[1] + a[2] * b[2]
 }
 
@@ -1127,6 +1270,177 @@ mod tests {
     }
 
     #[test]
+    fn foundation_world_queries_use_the_canonical_cut_and_generated_band() {
+        use crate::{build_basement_foundations, BasementFootprint};
+
+        let mesh = synthetic_land_mesh();
+        let outer = [[0.1, 0.1], [0.9, 0.1], [0.9, 0.9], [0.1, 0.9]];
+        let inner = [
+            [0.4, 0.4, 1.0],
+            [0.6, 0.4, 1.0],
+            [0.6, 0.6, 1.0],
+            [0.4, 0.6, 1.0],
+        ];
+        let foundation = build_basement_foundations(
+            &mesh,
+            &[BasementFootprint {
+                outer: &outer,
+                inner: &inner,
+                flatten_inner_to_mean_outer_z: true,
+            }],
+        )
+        .expect("foundation");
+        let expected_fragments = foundation.fragments.clone();
+        let expected_fragment_count = foundation.fragments.len();
+        let world = TerrainWorld::from_land_msh(&mesh)
+            .expect("source world")
+            .with_basement_foundation(foundation)
+            .expect("foundation world");
+
+        assert_eq!(world.surface_count(), expected_fragment_count);
+        assert_eq!(
+            world
+                .basement_foundation()
+                .map(|value| value.fragments.as_slice()),
+            Some(expected_fragments.as_slice())
+        );
+        assert!(world
+            .height_at([0.05, 0.05])
+            .expect("band height")
+            .is_some());
+        assert!(world.height_at([0.3, 0.3]).expect("band height").is_some());
+        assert_eq!(world.height_at([0.5, 0.5]).expect("hole height"), None);
+        assert_eq!(
+            world
+                .raycast([0.5, 0.5, 4.0], [0.0, 0.0, -1.0], FullSurfaceMask(0),)
+                .expect("hole ray"),
+            None
+        );
+
+        let mut generated_face = mesh.faces.len();
+        for (runtime, fragment) in world.surfaces.iter().zip(&expected_fragments) {
+            assert_eq!(runtime.mask, fragment.surface_mask);
+            assert_eq!(runtime.vertices, fragment.positions);
+            if let Some(source_face) = fragment.source_face {
+                assert_eq!(runtime.face, source_face);
+                assert_eq!(runtime.mask, mesh.faces[source_face].flags);
+            } else {
+                assert_eq!(runtime.face, generated_face);
+                generated_face += 1;
+            }
+        }
+        assert!(generated_face > mesh.faces.len());
+    }
+
+    #[test]
+    fn foundation_install_accepts_thin_nonzero_fragments_for_surface_queries() {
+        let mesh = synthetic_land_mesh();
+        let fragment = BasementFoundationFragment {
+            source_face: Some(0),
+            surface_mask: mesh.faces[0].flags,
+            positions: [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 0.000_000_5, 0.0]],
+        };
+        let foundation = BasementFoundation {
+            fragments: vec![fragment],
+            placements: Vec::new(),
+        };
+        let world = TerrainWorld::from_land_msh(&mesh)
+            .expect("source world")
+            .with_basement_foundation(foundation)
+            .expect("thin but non-degenerate triangle");
+
+        assert_eq!(world.surface_count(), 1);
+        assert_eq!(
+            world.height_at([0.25, 0.000_000_125]).expect("height"),
+            Some(0.0)
+        );
+        assert!(world
+            .raycast(
+                [0.25, 0.000_000_125, 1.0],
+                [0.0, 0.0, -1.0],
+                FullSurfaceMask(0),
+            )
+            .expect("raycast")
+            .is_some());
+    }
+
+    #[test]
+    fn foundation_install_rejects_invalid_source_identity_masks_and_generated_flags() {
+        let mesh = synthetic_land_mesh();
+        let source_triangle = [[0.0, 0.0, 0.0], [1.0, 0.0, 1.0], [0.0, 1.0, 1.0]];
+        let generated_triangle = [[0.0, 0.0, 0.0], [0.25, 0.0, 0.0], [0.0, 0.25, 0.0]];
+        let invalid = [
+            BasementFoundationFragment {
+                source_face: Some(mesh.faces.len()),
+                surface_mask: mesh.faces[0].flags,
+                positions: source_triangle,
+            },
+            BasementFoundationFragment {
+                source_face: Some(0),
+                surface_mask: FullSurfaceMask(mesh.faces[0].flags.0 ^ 1),
+                positions: source_triangle,
+            },
+            BasementFoundationFragment {
+                source_face: None,
+                surface_mask: FullSurfaceMask(0x0300),
+                positions: generated_triangle,
+            },
+        ];
+
+        for (index, fragment) in invalid.into_iter().enumerate() {
+            let foundation = BasementFoundation {
+                fragments: vec![fragment],
+                placements: Vec::new(),
+            };
+            assert!(
+                matches!(
+                    TerrainWorld::from_land_msh(&mesh)
+                        .expect("source world")
+                        .with_basement_foundation(foundation),
+                    Err(TerrainError::InvalidFoundationFragment { fragment: 0 })
+                ),
+                "invalid fragment case {index}"
+            );
+        }
+    }
+
+    #[test]
+    fn foundation_install_rejects_lod1_source_faces() {
+        use fparkan_terrain_format::TerrainLod;
+
+        let mut mesh = synthetic_land_mesh();
+        mesh.positions.push([0.0, 0.0, 2.0]);
+        mesh.positions.push([1.0, 0.0, 2.0]);
+        mesh.positions.push([0.0, 1.0, 2.0]);
+        mesh.faces
+            .push(terrain_face(FullSurfaceMask(0x0000_0003), [4, 5, 6]));
+        mesh.nodes_raw[10..12].copy_from_slice(&1_u16.to_le_bytes());
+        let mut lod1_slot = [0_u8; 68];
+        lod1_slot[0..2].copy_from_slice(&2_u16.to_le_bytes());
+        lod1_slot[2..4].copy_from_slice(&1_u16.to_le_bytes());
+        mesh.slots.slots_raw.push(lod1_slot);
+        assert_eq!(
+            mesh.face_lod_assignments().expect("LOD assignment")[2].lod,
+            TerrainLod::Lod1
+        );
+
+        let foundation = BasementFoundation {
+            fragments: vec![BasementFoundationFragment {
+                source_face: Some(2),
+                surface_mask: mesh.faces[2].flags,
+                positions: [[0.0, 0.0, 2.0], [1.0, 0.0, 2.0], [0.0, 1.0, 2.0]],
+            }],
+            placements: Vec::new(),
+        };
+        assert!(matches!(
+            TerrainWorld::from_land_msh(&mesh)
+                .expect("source world")
+                .with_basement_foundation(foundation),
+            Err(TerrainError::InvalidFoundationFragment { fragment: 0 })
+        ));
+    }
+
+    #[test]
     fn surface_bvh_reduces_candidates_without_changing_surface_math() {
         let surfaces: Vec<_> = (0_u16..16)
             .map(|face| {