/// 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}"
        );
    }
}