BREP_kernel 0.3.1

use super::*;

pub(super) fn nearly_equal_points(a: &NurbsSurface, b: &NurbsSurface, scale: f64) -> bool {
    if a.degree_u != b.degree_u
        || a.degree_v != b.degree_v
        || a.knots_u.len() != b.knots_u.len()
        || a.knots_v.len() != b.knots_v.len()
        || a.control_points.len() != b.control_points.len()
    {
        return false;
    }
    // Numerical knot dedup for exact-reconstruction recognition (see curve.rs).
    let knot_tolerance = crate::curve::KNOT_DEDUP_EPS;
    if a.knots_u
        .iter()
        .zip(&b.knots_u)
        .any(|(x, y)| (x - y).abs() > knot_tolerance)
        || a.knots_v
            .iter()
            .zip(&b.knots_v)
            .any(|(x, y)| (x - y).abs() > knot_tolerance)
    {
        return false;
    }
    let tolerance = RECOGNITION_TOLERANCE * scale.max(1.0);
    a.control_points
        .iter()
        .zip(&b.control_points)
        .all(|(row_a, row_b)| {
            row_a.len() == row_b.len()
                && row_a.iter().zip(row_b).all(|(p, q)| {
                    (p.x - q.x).abs() <= tolerance
                        && (p.y - q.y).abs() <= tolerance
                        && (p.z - q.z).abs() <= tolerance
                        && (p.w - q.w).abs() <= RECOGNITION_TOLERANCE
                })
        })
}

pub(super) fn surface_scale(surface: &NurbsSurface) -> f64 {
    surface
        .control_points
        .iter()
        .flatten()
        .map(|p| (p.x.abs() / p.w).max(p.y.abs() / p.w).max(p.z.abs() / p.w))
        .fold(0.0, f64::max)
}

pub(super) fn full_circle_knots(spans: usize) -> Vec<f64> {
    let mut knots = vec![0.0, 0.0, 0.0];
    for index in 1..spans {
        let knot = index as f64 / spans as f64;
        knots.extend([knot, knot]);
    }
    knots.extend([1.0, 1.0, 1.0]);
    knots
}

fn recognize_plane(surface: &NurbsSurface) -> Option<AnalyticSurface> {
    if !surface.is_affine().unwrap_or(false) {
        return None;
    }
    let u_domain = [surface.knots_u[0], *surface.knots_u.last()?];
    let v_domain = [surface.knots_v[0], *surface.knots_v.last()?];
    let extent_u = u_domain[1] - u_domain[0];
    let extent_v = v_domain[1] - v_domain[0];
    if extent_u <= 0.0 || extent_v <= 0.0 {
        return None;
    }
    let p00 = surface.control_points[0][0].point().ok()?;
    let p10 = surface.control_points[1][0].point().ok()?;
    let p01 = surface.control_points[0][1].point().ok()?;
    Some(AnalyticSurface::Plane {
        origin: p00,
        u_dir: p10.sub(p00).scale(1.0 / extent_u),
        v_dir: p01.sub(p00).scale(1.0 / extent_v),
        u_domain,
        v_domain,
    })
}

fn recognize_revolution(surface: &NurbsSurface) -> Option<AnalyticSurface> {
    if surface.degree_u != 2 || surface.control_points.len() != 9 {
        return None;
    }
    let expected_knots = full_circle_knots(4);
    if surface.knots_u.len() != expected_knots.len()
        || surface
            .knots_u
            .iter()
            .zip(&expected_knots)
            .any(|(a, b)| (a - b).abs() > 1e-12)
    {
        return None;
    }
    let scale = surface_scale(surface);
    let rows = &surface.control_points;
    // Antipodal full-circle control points straddle the axis; their midpoint
    // lies on it for every generatrix column.
    let mut frame: Option<RevolutionFrame> = None;
    for column in 0..rows[0].len() {
        let p0 = rows[0][column].point().ok()?;
        let p180 = rows[4][column].point().ok()?;
        let p90 = rows[2][column].point().ok()?;
        let center = p0.add(p180).scale(0.5);
        let radial = p0.sub(center);
        let radius = radial.length();
        if radius <= RECOGNITION_TOLERANCE * scale.max(1.0) {
            continue;
        }
        let toward_90 = p90.sub(center);
        let axis = radial.cross(toward_90).normalized().ok()?;
        let x_axis = radial.scale(1.0 / radius);
        let y_axis = axis.cross(x_axis);
        frame = Some(RevolutionFrame {
            origin: center,
            axis,
            x_axis,
            y_axis,
        });
        break;
    }
    let frame = frame?;
    let generatrix =
        NurbsCurve::new(surface.degree_v, surface.knots_v.clone(), rows[0].clone()).ok()?;
    let rebuilt =
        make_revolution(frame.origin, frame.axis, &generatrix, std::f64::consts::TAU).ok()?;
    if !nearly_equal_points(surface, &rebuilt, scale) {
        return None;
    }
    classify_generatrix(surface, frame, &generatrix, scale)
}

fn classify_generatrix(
    surface: &NurbsSurface,
    frame: RevolutionFrame,
    generatrix: &NurbsCurve,
    scale: f64,
) -> Option<AnalyticSurface> {
    let tolerance = RECOGNITION_TOLERANCE * scale.max(1.0);
    let controls = &generatrix.control_points;
    if generatrix.degree == 1
        && controls.len() == 2
        && (controls[0].w - 1.0).abs() <= RECOGNITION_TOLERANCE
        && (controls[1].w - 1.0).abs() <= RECOGNITION_TOLERANCE
        && (surface.knots_v[0], *surface.knots_v.last()?) == (0.0, 1.0)
        && generatrix_is_meridional_half_ray(generatrix, &frame)
    {
        let p0 = controls[0].point().ok()?;
        let p1 = controls[1].point().ok()?;
        let (_, rho0, z0) = frame.cylindrical(p0);
        let (_, rho1, z1) = frame.cylindrical(p1);
        let height = z1 - z0;
        if height.abs() <= tolerance {
            return None;
        }
        // Rebase the frame origin to the generatrix start's axial position
        // so v = axial / height exactly.
        let origin = frame.origin.add(frame.axis.scale(z0));
        return Some(AnalyticSurface::RuledRevolution {
            frame: RevolutionFrame { origin, ..frame },
            rho0,
            rho1,
            height,
        });
    }
    if generatrix.degree == 2 && controls.len() == 5 {
        // Polar meridian: poles at both ends, equator at the middle column.
        let south = controls[0].point().ok()?;
        let north = controls[4].point().ok()?;
        let (_, rho_south, z_south) = frame.cylindrical(south);
        let (_, rho_north, z_north) = frame.cylindrical(north);
        if rho_south <= tolerance && rho_north <= tolerance {
            let radius = (z_north - z_south) * 0.5;
            if radius <= tolerance {
                return None;
            }
            let center = frame.origin.add(frame.axis.scale(z_south + radius));
            let meridian = make_arc(
                center,
                frame.x_axis,
                frame.axis,
                radius,
                -std::f64::consts::FRAC_PI_2,
                std::f64::consts::FRAC_PI_2,
            )
            .ok()?;
            if curves_nearly_equal(generatrix, &meridian, scale) {
                return Some(AnalyticSurface::Sphere {
                    frame: RevolutionFrame {
                        origin: center,
                        ..frame
                    },
                    radius,
                });
            }
        }
        return None;
    }
    if generatrix.degree == 2 && controls.len() == 9 {
        let tube_start = controls[0].point().ok()?;
        let tube_opposite = controls[4].point().ok()?;
        let tube_center = tube_start.add(tube_opposite).scale(0.5);
        let minor = tube_start.sub(tube_center).length();
        let (_, major, axial) = frame.cylindrical(tube_center);
        if minor <= tolerance || major <= minor || axial.abs() > tolerance {
            return None;
        }
        let tube = make_arc(
            tube_center,
            frame.x_axis,
            frame.axis,
            minor,
            0.0,
            std::f64::consts::TAU,
        )
        .ok()?;
        if curves_nearly_equal(generatrix, &tube, scale) {
            return Some(AnalyticSurface::Torus {
                frame,
                major_radius: major,
                minor_radius: minor,
            });
        }
        return None;
    }
    None
}

fn curves_nearly_equal(a: &NurbsCurve, b: &NurbsCurve, scale: f64) -> bool {
    if a.degree != b.degree
        || a.knots.len() != b.knots.len()
        || a.control_points.len() != b.control_points.len()
    {
        return false;
    }
    if a.knots
        .iter()
        .zip(&b.knots)
        .any(|(x, y)| (x - y).abs() > 1e-12)
    {
        return false;
    }
    let tolerance = RECOGNITION_TOLERANCE * scale.max(1.0);
    a.control_points
        .iter()
        .zip(&b.control_points)
        .all(|(p, q)| {
            (p.x - q.x).abs() <= tolerance
                && (p.y - q.y).abs() <= tolerance
                && (p.z - q.z).abs() <= tolerance
                && (p.w - q.w).abs() <= RECOGNITION_TOLERANCE
        })
}

pub fn recognize(surface: &NurbsSurface) -> Option<AnalyticSurface> {
    recognize_plane(surface)
        .or_else(|| recognize_revolution(surface))
        .or_else(|| recognize_general_revolution(surface))
}

/// Fallback recognizer: any `make_revolution` product the specialized
/// quadric recognizers above rejected (arbitrary generatrix, partial
/// sweep), validated by exact reconstruction in `revolution_structure`.
fn recognize_general_revolution(surface: &NurbsSurface) -> Option<AnalyticSurface> {
    let structure = revolution_structure(surface)?;
    let spans = (surface.control_points.len() - 1) / 2;
    Some(AnalyticSurface::Revolution {
        frame: structure.frame,
        spans,
        sweep: structure.sweep,
        generatrix: structure.generatrix,
    })
}

/// Whether the generatrix lies in the axis frame's starting meridian
/// half-plane.  Only then can a surface projection be reduced to projection
/// onto the unrotated 3D generatrix.
pub(super) fn generatrix_is_meridional_half_ray(
    generatrix: &NurbsCurve,
    frame: &RevolutionFrame,
) -> bool {
    let mut scale = 0.0_f64;
    for control in &generatrix.control_points {
        let Ok(point) = control.point() else {
            return false;
        };
        let offset = point.sub(frame.origin);
        let radial = offset.sub(frame.axis.scale(offset.dot(frame.axis)));
        scale = scale.max(radial.length());
    }
    let tolerance = RECOGNITION_TOLERANCE * scale.max(1.0);
    generatrix.control_points.iter().all(|control| {
        let Ok(point) = control.point() else {
            return false;
        };
        let offset = point.sub(frame.origin);
        let radial = offset.sub(frame.axis.scale(offset.dot(frame.axis)));
        radial.dot(frame.y_axis).abs() <= tolerance
            && radial.dot(frame.x_axis) >= -tolerance
    })
}