aboutsummaryrefslogtreecommitdiff
path: root/crates/fparkan-terrain
diff options
context:
space:
mode:
Diffstat (limited to 'crates/fparkan-terrain')
-rw-r--r--crates/fparkan-terrain/src/contour.rs809
-rw-r--r--crates/fparkan-terrain/src/foundation.rs1806
-rw-r--r--crates/fparkan-terrain/src/lib.rs362
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| {
Copyright © 2026 Valentin A. Popov. All Rights Reverse Engineered.