#include "vde/brep/ssi_boolean.h"
#include "vde/brep/brep_face_split.h"
#include "vde/brep/brep_heal.h"
#include "vde/brep/tolerant_modeling.h"
#include "vde/curves/surface_intersection.h"
#include "vde/core/exact_predicates.h"
#include "vde/core/aabb.h"
#include "vde/core/point.h"
#include "vde/mesh/halfedge_mesh.h"
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <iterator>
#include <limits>
#include <map>
#include <mutex>
#include <random>
#include <set>
#include <unordered_map>
#include <unordered_set>
#include <vector>

#ifdef VDE_USE_OPENMP
#include <omp.h>
#endif

using vde::core::AABB3D;
using vde::core::Point2D;
using vde::core::Point3D;
using vde::core::Vector2D;
using vde::core::Vector3D;
using vde::core::Transform3D;
using vde::core::Plane3D;
namespace vde::brep {

using core::Point3D;
using core::Vector3D;
using core::AABB3D;

// Forward declarations (defined after anonymous namespace)
AABB3D face_bounds(const BrepModel& body, int face_id);
Vector3D face_normal(const BrepModel& body, int face_id);
struct SSICurveData;
using ClassResult;
using SSICurveData;

/// SSI curve data for a face pair
struct SSICurveData {
    int face_a;
    int face_b;
    std::vector<Point3D> points;
    Vector3D cutting_normal;
    Point3D cutting_center;
    bool valid = false;
};

namespace {

// ═══════════════════════════════════════════════════════════
//  Mesh caching
// ═══════════════════════════════════════════════════════════

/// ClassResult is now defined in ssi_boolean.h (public API)
// ═══════════════════════════════════════════════════════════
//  Constants & types
BrepModel extract_face_fragment(const BrepModel& body, int face_id);
ClassResult classify_point_adaptive(
    const BrepModel& body, const Point3D& p, int& num_exact_upgrades);

// ═══════════════════════════════════════════════════════════
//  Mesh caching (shared with brep_boolean.cpp pattern)
// ═══════════════════════════════════════════════════════════

static std::unordered_map<const BrepModel*, mesh::HalfedgeMesh> mesh_cache;
static std::mutex mesh_cache_mutex;

static const mesh::HalfedgeMesh& get_mesh(const BrepModel& body) {
    {
        std::lock_guard<std::mutex> lock(mesh_cache_mutex);
        auto it = mesh_cache.find(&body);
        if (it != mesh_cache.end()) return it->second;
    }
    auto mesh = body.to_mesh(0.05);
    std::lock_guard<std::mutex> lock(mesh_cache_mutex);
    auto [it, _] = mesh_cache.emplace(&body, std::move(mesh));
    return it->second;
}

// ═══════════════════════════════════════════════════════════
//  PrecisionTracker implementation (public API)
// ═══════════════════════════════════════════════════════════

}  // anonymous namespace — public functions follow

void PrecisionTracker::record(const std::string& operation, double precision_loss) {
    accumulated_ += std::abs(precision_loss);
    count_++;
    last_op_ = operation;
}

bool PrecisionTracker::exceeds_threshold(double threshold) const {
    return accumulated_ > threshold;
}

void PrecisionTracker::reset() {
    accumulated_ = 0.0;
    count_ = 0;
    last_op_.clear();
}

// ═══════════════════════════════════════════════════════════
//  Degeneracy detection: 7 categories
// ═══════════════════════════════════════════════════════════

/// Compute face-plane coefficients (ax+by+cz+d=0) from face normal and a point on the face
struct FacePlane {
    double a, b, c, d;  // ax + by + cz + d = 0
};

FacePlane face_plane(const BrepModel& body, int face_id) {
    Vector3D n = face_normal(body, face_id);
    n.normalize();
    // Use first edge's start vertex as a point on the face
    auto es = body.face_edges(face_id);
    Point3D p;
    if (!es.empty()) {
        p = body.vertex(body.edge(es[0]).v_start).point;
    }
    double d = -(n.x() * p.x() + n.y() * p.y() + n.z() * p.z());
    return {n.x(), n.y(), n.z(), d};
}

// ── Category 1: Surface overlap (coplanar) ──

SurfaceOverlapResult detect_surface_overlap(
    const BrepModel& a, const BrepModel& b,
    int face_a, int face_b, double tol)
{
    SurfaceOverlapResult result;

    auto fp_a = face_plane(a, face_a);
    auto fp_b = face_plane(b, face_b);

    // Compare normals: |dot| ≈ 1 means parallel (coplanar or opposite)
    double dot = std::abs(fp_a.a * fp_b.a + fp_a.b * fp_b.b + fp_a.c * fp_b.c);
    if (dot < 0.9999) return result;  // not parallel enough

    // Compare plane distances
    Vector3D na(fp_a.a, fp_a.b, fp_a.c);
    Vector3D nb(fp_b.a, fp_b.b, fp_b.c);

    // Get a point on face B and compute signed distance to plane A
    auto es = b.face_edges(face_b);
    if (es.empty()) return result;
    Point3D p = b.vertex(b.edge(es[0]).v_start).point;
    double dist = std::abs(fp_a.a * p.x() + fp_a.b * p.y() + fp_a.c * p.z() + fp_a.d);

    if (dist < tol * 10.0) {
        // Check bounding box overlap
        AABB3D bb_a = face_bounds(a, face_a);
        AABB3D bb_b = face_bounds(b, face_b);
        if (bb_a.intersects(bb_b)) {
            result.overlapping = true;
        }
    }

    return result;
}

// ── Category 2: Curve overlap (collinear) ──

CurveOverlapResult detect_curve_overlap(
    const std::vector<Point3D>& curve_a,
    const std::vector<Point3D>& curve_b,
    double tol)
{
    CurveOverlapResult result;
    if (curve_a.size() < 2 || curve_b.size() < 2) return result;

    // Compute direction vectors from endpoints
    Vector3D dir_a = (curve_a.back() - curve_a.front()).normalized();
    Vector3D dir_b = (curve_b.back() - curve_b.front()).normalized();

    double dot = std::abs(dir_a.dot(dir_b));
    if (dot < 0.9999) return result;

    // Check if they lie on the same line: project each point
    Vector3D perp = dir_a.cross(Vector3D(1, 0, 0));
    if (perp.norm() < 1e-6) perp = dir_a.cross(Vector3D(0, 1, 0));
    perp.normalize();

    double max_dev = 0.0;
    for (const auto& p : curve_a) {
        double dev = std::abs((p - curve_a.front()).dot(perp));
        max_dev = std::max(max_dev, dev);
    }
    for (const auto& p : curve_b) {
        double dev = std::abs((p - curve_a.front()).dot(perp));
        max_dev = std::max(max_dev, dev);
    }
    if (max_dev > tol) return result;

    // Both curves are collinear; compute parameter ranges along dir_a
    auto proj_a0 = (curve_a.front() - curve_a.front()).dot(dir_a);  // 0
    auto proj_a1 = (curve_a.back() - curve_a.front()).dot(dir_a);
    auto proj_b0 = (curve_b.front() - curve_a.front()).dot(dir_a);
    auto proj_b1 = (curve_b.back() - curve_a.front()).dot(dir_a);

    double lo_a = std::min(proj_a0, proj_a1), hi_a = std::max(proj_a0, proj_a1);
    double lo_b = std::min(proj_b0, proj_b1), hi_b = std::max(proj_b0, proj_b1);

    // Check overlap
    double overlap_lo = std::max(lo_a, lo_b);
    double overlap_hi = std::min(hi_a, hi_b);
    if (overlap_hi > overlap_lo + tol) {
        result.overlapping = true;
        result.t_start_a = (overlap_lo - lo_a) / (hi_a - lo_a + 1e-16);
        result.t_end_a   = (overlap_hi - lo_a) / (hi_a - lo_a + 1e-16);
        result.t_start_b = (overlap_lo - lo_b) / (hi_b - lo_b + 1e-16);
        result.t_end_b   = (overlap_hi - lo_b) / (hi_b - lo_b + 1e-16);
    }

    return result;
}

// ── Category 3: High-order contact (G2+) ──

HighOrderContactResult detect_high_order_contact(
    const BrepModel& a, const BrepModel& b,
    int face_a, int face_b, double tol)
{
    HighOrderContactResult result;

    const auto& sa = a.surface(a.face(face_a).surface_id);
    const auto& sb = b.surface(b.face(face_b).surface_id);

    // Sample surface normals at multiple points
    double ua_mid = (sa.knots_u().front() + sa.knots_u().back()) * 0.5;
    double va_mid = (sa.knots_v().front() + sa.knots_v().back()) * 0.5;

    // Check multiple sample points for near-tangent contact
    Point3D pa = sa.evaluate(ua_mid, va_mid);
    Vector3D na = sa.normal(ua_mid, va_mid);
    if (na.norm() < 1e-12) return result;
    na.normalize();

    // Find the closest point on surface B to pa
    // Simple approach: sample B's surface and find closest
    double best_dist = std::numeric_limits<double>::max();
    Point3D best_pb;
    Vector3D best_nb;
    for (int su = 0; su <= 4; ++su) {
        for (int sv = 0; sv <= 4; ++sv) {
            double ub = sb.knots_u().front() + (sb.knots_u().back() - sb.knots_u().front()) * su / 4.0;
            double vb = sb.knots_v().front() + (sb.knots_v().back() - sb.knots_v().front()) * sv / 4.0;
            Point3D pb = sb.evaluate(ub, vb);
            double d = (pa - pb).norm();
            if (d < best_dist) {
                best_dist = d;
                best_pb = pb;
                best_nb = sb.normal(ub, vb);
            }
        }
    }

    if (best_nb.norm() < 1e-12) return result;
    best_nb.normalize();

    // G0: point on surface → G1: tangent → G2: curvature match
    bool g0_contact = (best_dist < tol * 10.0);
    bool g1_contact = (std::abs(na.dot(best_nb)) > 0.9999);

    if (g0_contact && g1_contact) {
        // High-order contact detected
        result.has_contact = true;
        result.contact_point = pa;
        result.offset_direction = na;
        result.offset_magnitude = tol * 100.0;  // micro-offset
    }

    return result;
}

// ── Category 4: Boundary contact ──

BoundaryContactResult detect_boundary_contact(
    const std::vector<Point3D>& intersection_points,
    const BrepModel& a, const BrepModel& b,
    int face_a, int face_b, double tol)
{
    BoundaryContactResult result;
    if (intersection_points.size() < 2) return result;

    auto es_a = a.face_edges(face_a);
    auto es_b = b.face_edges(face_b);

    // Build boundary segments for both faces
    auto build_segments = [&](const BrepModel& body, const std::vector<int>& edges) {
        std::vector<std::pair<Point3D, Point3D>> segs;
        for (int ei : edges) {
            const auto& e = body.edge(ei);
            segs.emplace_back(body.vertex(e.v_start).point,
                              body.vertex(e.v_end).point);
        }
        return segs;
    };

    auto segs_a = build_segments(a, es_a);
    auto segs_b = build_segments(b, es_b);

    // Check if intersection endpoints are near any boundary segment
    auto point_near_segment = [&](const Point3D& p, const std::vector<std::pair<Point3D, Point3D>>& segs) -> bool {
        for (const auto& [s, e] : segs) {
            Vector3D se = e - s;
            double len_sq = se.squaredNorm();
            if (len_sq < 1e-16) {
                if ((p - s).norm() < tol) return true;
                continue;
            }
            double t = (p - s).dot(se) / len_sq;
            t = std::max(0.0, std::min(1.0, t));
            Point3D closest = s + se * t;
            if ((p - closest).norm() < tol) return true;
        }
        return false;
    };

    bool start_near_boundary = point_near_segment(intersection_points.front(), segs_a) ||
                                point_near_segment(intersection_points.front(), segs_b);
    bool end_near_boundary   = point_near_segment(intersection_points.back(), segs_a) ||
                                point_near_segment(intersection_points.back(), segs_b);

    if (start_near_boundary || end_near_boundary) {
        result.has_boundary_contact = true;

        // Extend endpoints along the tangent direction
        Vector3D tangent = (intersection_points.back() - intersection_points.front()).normalized();
        double extend_len = tol * 100.0;

        result.extended_start = intersection_points.front() - tangent * extend_len;
        result.extended_end   = intersection_points.back() + tangent * extend_len;
    }

    return result;
}

// ── Category 5: Vertex contact (multiple faces converging at a point) ──

VertexContactResult detect_vertex_contact(
    const std::vector<Point3D>& intersection_points,
    const BrepModel& a, const BrepModel& b,
    int face_a, int face_b, double tol)
{
    VertexContactResult result;

    // Count how many distinct vertices from both models converge near each intersection point
    std::vector<Point3D> all_verts;
    for (size_t vi = 0; vi < a.num_vertices(); ++vi)
        all_verts.push_back(a.vertex(static_cast<int>(vi)).point);
    for (size_t vi = 0; vi < b.num_vertices(); ++vi)
        all_verts.push_back(b.vertex(static_cast<int>(vi)).point);

    // For each intersection point, find nearby vertices
    for (const auto& ip : intersection_points) {
        std::vector<int> nearby_faces_a, nearby_faces_b;
        int verts_nearby = 0;

        for (size_t vi = 0; vi < a.num_vertices(); ++vi) {
            if ((a.vertex(static_cast<int>(vi)).point - ip).norm() < tol * 10.0) {
                verts_nearby++;
                // Find faces incident to this vertex
                auto ve = a.vertex_edges(static_cast<int>(vi));
                std::set<int> faces;
                for (int ei : ve) {
                    auto ef = a.edge_faces(ei);
                    faces.insert(ef.begin(), ef.end());
                }
                for (int fid : faces) {
                    if (std::find(nearby_faces_a.begin(), nearby_faces_a.end(), fid) == nearby_faces_a.end())
                        nearby_faces_a.push_back(fid);
                }
            }
        }
        for (size_t vi = 0; vi < b.num_vertices(); ++vi) {
            if ((b.vertex(static_cast<int>(vi)).point - ip).norm() < tol * 10.0) {
                verts_nearby++;
                auto ve = b.vertex_edges(static_cast<int>(vi));
                std::set<int> faces;
                for (int ei : ve) {
                    auto ef = b.edge_faces(ei);
                    faces.insert(ef.begin(), ef.end());
                }
                for (int fid : faces) {
                    if (std::find(nearby_faces_b.begin(), nearby_faces_b.end(), fid) == nearby_faces_b.end())
                        nearby_faces_b.push_back(fid);
                }
            }
        }

        // If 3+ faces converge at one point (counting both bodies)
        if (nearby_faces_a.size() + nearby_faces_b.size() >= 3) {
            result.has_vertex_contact = true;
            result.contact_point = ip;
            result.involved_faces_a = nearby_faces_a;
            result.involved_faces_b = nearby_faces_b;
            result.sphere_radius = tol * 100.0;
            break;
        }
    }

    return result;
}

// ── Category 6: Non-manifold contact ──

NonManifoldContactResult detect_non_manifold_contact(
    const BrepModel& body, int edge_id)
{
    NonManifoldContactResult result;
    result.edge_id = edge_id;

    auto faces = body.edge_faces(edge_id);
    result.incident_faces = faces;
    result.is_non_manifold = (faces.size() >= 3);

    return result;
}

// ═══════════════════════════════════════════════════════════
//  Face utilities (shared — used by both public and internal code)
// ═══════════════════════════════════════════════════════════

}  // anonymous namespace — shared face utilities follow

/// Get bounding box of a face from its edges/vertices
AABB3D face_bounds(const BrepModel& body, int face_id) {
    AABB3D bb;
    auto es = body.face_edges(face_id);
    for (int ei : es) {
        const auto& e = body.edge(ei);
        bb.expand(body.vertex(e.v_start).point);
        bb.expand(body.vertex(e.v_end).point);
    }
    return bb;
}

/// Compute face centroid from edge vertices
Point3D face_centroid(const BrepModel& body, int face_id) {
    auto edges = body.face_edges(face_id);
    Point3D c(0, 0, 0);
    int count = 0;
    for (int ei : edges) {
        const auto& e = body.edge(ei);
        c += body.vertex(e.v_start).point; count++;
        c += body.vertex(e.v_end).point; count++;
    }
    if (count > 0) c /= static_cast<double>(count);
    return c;
}

/// Compute face normal from its surface
Vector3D face_normal(const BrepModel& body, int face_id) {
    const auto& face = body.face(face_id);
    const auto& surf = body.surface(face.surface_id);
    Vector3D n = surf.normal(0.5, 0.5);
    if (n.norm() < 1e-12) n = Vector3D::UnitZ();
    return n;
}

// ═══════════════════════════════════════════════════════════
//  Extract a face as a standalone BrepModel fragment
// ═══════════════════════════════════════════════════════════

BrepModel extract_face_fragment(const BrepModel& body, int face_id) {
    BrepModel frag;
    auto& f = body.face(face_id);
    auto es = body.face_edges(face_id);

    std::map<int, int> old_to_new;
    for (int ei : es) {
        auto& e = body.edge(ei);
        if (old_to_new.find(e.v_start) == old_to_new.end()) {
            old_to_new[e.v_start] = frag.add_vertex(body.vertex(e.v_start).point);
        }
        if (old_to_new.find(e.v_end) == old_to_new.end()) {
            old_to_new[e.v_end] = frag.add_vertex(body.vertex(e.v_end).point);
        }
    }

    int sid = frag.add_surface(body.surface(f.surface_id));
    std::vector<int> new_edges;
    for (int ei : es) {
        auto& e = body.edge(ei);
        new_edges.push_back(frag.add_edge(old_to_new[e.v_start], old_to_new[e.v_end]));
    }
    int loop = frag.add_loop(new_edges, true);
    int nf = frag.add_face(sid, {loop});
    int sh = frag.add_shell({nf}, false);
    frag.add_body({sh}, "face_" + std::to_string(face_id));
    return frag;
}

// ═══════════════════════════════════════════════════════════
//  Step 1: Parallel SSI computation for all face pairs
// ═══════════════════════════════════════════════════════════

/// Compute all SSI curves between body A and body B.
/// Returns a list of SSICurveData, one per intersecting face pair.
/// Uses OpenMP parallelization over face pairs.
std::vector<SSICurveData> compute_all_ssi(
    const BrepModel& a, const BrepModel& b,
    int& num_exact_upgrades)
{
    std::vector<SSICurveData> all_isects;
    std::mutex isect_mutex;

    // Build face-AABB cache for quick rejection
    std::vector<AABB3D> bb_a(a.num_faces());
    for (size_t i = 0; i < a.num_faces(); ++i) {
        bb_a[i] = face_bounds(a, static_cast<int>(i));
    }
    std::vector<AABB3D> bb_b(b.num_faces());
    for (size_t j = 0; j < b.num_faces(); ++j) {
        bb_b[j] = face_bounds(b, static_cast<int>(j));
    }

    size_t total_pairs = a.num_faces() * b.num_faces();
    // Flatten face pairs for parallel iteration
    std::vector<std::pair<int,int>> pairs;
    pairs.reserve(total_pairs);
    for (size_t i = 0; i < a.num_faces(); ++i) {
        for (size_t j = 0; j < b.num_faces(); ++j) {
            pairs.emplace_back(static_cast<int>(i), static_cast<int>(j));
        }
    }

#ifdef VDE_USE_OPENMP
#pragma omp parallel
    {
        std::vector<SSICurveData> local;
#pragma omp for schedule(dynamic) nowait
        for (size_t pi = 0; pi < pairs.size(); ++pi) {
#else
        std::vector<SSICurveData> local;
        for (size_t pi = 0; pi < pairs.size(); ++pi) {
#endif
            int ia = pairs[pi].first;
            int ib = pairs[pi].second;

            // Quick AABB rejection
            if (!bb_a[ia].intersects(bb_b[ib])) continue;

            const auto& surf_a = a.surface(a.face(ia).surface_id);
            const auto& surf_b = b.surface(b.face(ib).surface_id);

            auto curves = curves::intersect_surfaces(surf_a, surf_b, 4, 1e-3);

            for (const auto& curve : curves) {
                if (curve.points.size() < 4) continue;

                SSICurveData d;
                d.face_a = ia;
                d.face_b = ib;
                d.points = curve.points;

                // Compute centroid
                d.cutting_center = Point3D(0, 0, 0);
                for (const auto& p : d.points) d.cutting_center += p;
                d.cutting_center /= static_cast<double>(d.points.size());

                // Compute cutting plane normal from SSI points
                Vector3D n(0, 0, 0);
                if (d.points.size() >= 3) {
                    Vector3D d1 = d.points[1] - d.points[0];
                    Vector3D d2 = d.points.back() - d.points[0];
                    n = d1.cross(d2);
                }
                if (n.norm() < 1e-12) {
                    n = face_normal(b, ib);
                }
                d.cutting_normal = n;
                d.valid = true;
                local.push_back(std::move(d));
            }
        }

#ifdef VDE_USE_OPENMP
#pragma omp critical
        {
            all_isects.insert(all_isects.end(),
                std::make_move_iterator(local.begin()),
                std::make_move_iterator(local.end()));
        }
    }
#else
        all_isects = std::move(local);
#endif

    return all_isects;
}

// ═══════════════════════════════════════════════════════════
//  Step 2: Split faces using SSI curves
// ═══════════════════════════════════════════════════════════

/// Build a per-face list of cutting planes from SSI curves
struct CuttingPlane {
    Point3D point;
    Vector3D normal;
};

std::unordered_map<int, std::vector<CuttingPlane>>
build_face_cutting_planes(const std::vector<SSICurveData>& ssi_curves,
                          bool for_body_a) {
    std::unordered_map<int, std::vector<CuttingPlane>> result;
    for (const auto& c : ssi_curves) {
        if (!c.valid) continue;
        int face_id = for_body_a ? c.face_a : c.face_b;
        result[face_id].push_back({c.cutting_center, c.cutting_normal});
    }
    // Deduplicate: merge close cutting planes within each face
    for (auto& [fid, planes] : result) {
        std::vector<CuttingPlane> deduped;
        for (auto& p : planes) {
            bool duplicate = false;
            for (auto& existing : deduped) {
                double dist = (p.point - existing.point).norm();
                double dot = std::abs(p.normal.normalized().dot(existing.normal.normalized()));
                if (dist < 1e-4 && dot > 0.99) {
                    duplicate = true;
                    break;
                }
            }
            if (!duplicate) deduped.push_back(p);
        }
        planes = std::move(deduped);
    }
    return result;
}

/// Split all faces of a body using SSI-derived cutting planes,
/// then classify each resulting fragment against the other body.
/// Returns a vector of (fragment, classification) pairs.
std::vector<std::pair<BrepModel, ClassResult>>
split_and_classify_faces_ssi(
    const BrepModel& body,
    const std::unordered_map<int, std::vector<CuttingPlane>>& cutting_planes,
    const BrepModel& other_body,
    int& num_exact_upgrades)
{
    std::vector<std::pair<BrepModel, ClassResult>> result;
    std::mutex result_mutex;

#ifdef VDE_USE_OPENMP
#pragma omp parallel for schedule(dynamic)
#endif
    for (size_t fi = 0; fi < body.num_faces(); ++fi) {
        int fid = static_cast<int>(fi);
        std::vector<BrepModel> fragments;
        fragments.push_back(extract_face_fragment(body, fid));

        // Split by SSI-derived cutting planes for this face
        auto it = cutting_planes.find(fid);
        if (it != cutting_planes.end()) {
            // Static max fragment limit to prevent explosion
            static constexpr size_t MAX_FRAGMENTS = 30;
            int no_progress = 0;
            for (const auto& plane : it->second) {
                size_t before = fragments.size();
                std::vector<BrepModel> new_fragments;
                for (auto& frag : fragments) {
                    auto splits = split_face_by_plane(frag, 0, plane.point, plane.normal);
                    if (splits.size() > 1) {
                        for (auto& s : splits) new_fragments.push_back(std::move(s));
                    } else {
                        new_fragments.push_back(std::move(frag));
                    }
                }
                fragments = std::move(new_fragments);
                if (fragments.empty()) break;
                if (fragments.size() > MAX_FRAGMENTS) break;
                if (fragments.size() == before && ++no_progress >= 10) break;
                else if (fragments.size() != before) no_progress = 0;
            }
        }

        // Classify each fragment using multi-point sampling
        std::vector<std::pair<BrepModel, ClassResult>> local;
        for (auto& frag : fragments) {
            // Use the new multi-point classify_face_fragment
            // Extract face_id from the fragment (fragment has exactly one face at index 0)
            auto cls = classify_face_fragment(body, other_body, fid, num_exact_upgrades);
            local.emplace_back(std::move(frag), cls);
        }

#ifdef VDE_USE_OPENMP
        std::lock_guard<std::mutex> lock(result_mutex);
#endif
        for (auto& p : local) result.push_back(std::move(p));
    }

    return result;
}

// ═══════════════════════════════════════════════════════════
//  Step 3: Classification with adaptive precision
// ═══════════════════════════════════════════════════════════

/// Fast classification: sample surface points + AABB sphere test.
/// Returns true if ALL samples consistently classify as uniform_result.
/// This avoids expensive mesh construction for clearly IN/OUT faces.
bool classify_face_uniform(const BrepModel& face_body,
                            const BrepModel& other_body,
                            ClassResult& uniform_result) {
    if (face_body.num_faces() == 0) return false;
    const auto& surf = face_body.surface(face_body.face(0).surface_id);
    double u0 = surf.knots_u().front(), u1 = surf.knots_u().back();
    double v0 = surf.knots_v().front(), v1 = surf.knots_v().back();

    double us[] = { (u0+u1)*0.5, u0, u1, u0, u1 };
    double vs[] = { (v0+v1)*0.5, v0, v0, v1, v1 };

    auto bb = other_body.bounds();
    Point3D center = bb.center();
    double radius = (bb.max() - bb.min()).norm() * 0.5;
    double radius_sq = radius * radius;

    ClassResult first;
    {
        Point3D p = surf.evaluate(us[0], vs[0]);
        if (!bb.contains(p)) { first = OUT; }
        else if ((p - center).squaredNorm() <= radius_sq * 0.5) { first = IN; }
        else return false;
    }
    for (int i = 1; i < 5; ++i) {
        Point3D p = surf.evaluate(us[i], vs[i]);
        ClassResult cls;
        if (!bb.contains(p)) cls = OUT;
        else if ((p - center).squaredNorm() <= radius_sq * 0.5) cls = IN;
        else return false;
        if (cls != first) return false;
    }
    uniform_result = first;
    return true;
}

/// Double-precision ray cast classification for a point against a mesh.
/// Returns IN/OUT (never ON — ON is handled separately).
/// Uses exact_orient3d for robust coplanarity tests on boundary triangles.
ClassResult classify_point_double(const BrepModel& body, const Point3D& p) {
    const mesh::HalfedgeMesh& mesh = get_mesh(body);
    if (mesh.num_faces() == 0) return OUT;

    // Ray in +X direction, count intersections
    Vector3D ray_dir(1, 0, 0);
    int hits = 0;
    double tol = 1e-8;

    for (size_t fi = 0; fi < mesh.num_faces(); ++fi) {
        auto& face = mesh.face(fi);
        int h0 = face.halfedge_index;
        int h1 = mesh.halfedge(h0).next_index;
        int h2 = mesh.halfedge(h1).next_index;

        Point3D v0 = mesh.vertex(mesh.halfedge(h0).vertex_index);
        Point3D v1 = mesh.vertex(mesh.halfedge(h1).vertex_index);
        Point3D v2 = mesh.vertex(mesh.halfedge(h2).vertex_index);

        // Möller–Trumbore
        Vector3D e1 = v1 - v0;
        Vector3D e2 = v2 - v0;
        Vector3D h = ray_dir.cross(e2);
        double a = e1.dot(h);
        if (std::abs(a) < tol) continue;
        double f = 1.0 / a;
        Vector3D s = p - v0;
        double u = f * s.dot(h);
        if (u < 0.0 || u > 1.0) continue;
        Vector3D q = s.cross(e1);
        double v = f * ray_dir.dot(q);
        if (v < 0.0 || u + v > 1.0) continue;
        double t = f * e2.dot(q);

        // Tolerance check: use exact_orient3d for near-boundary cases
        if (t > -tol && t < tol) {
            // Near-edge or near-vertex: use exact_orient3d to resolve
            double orient = core::exact_orient3d(
                v0.x(), v0.y(), v0.z(),
                v1.x(), v1.y(), v1.z(),
                v2.x(), v2.y(), v2.z(),
                p.x(), p.y(), p.z());
            if (std::abs(orient) < 1e-12) {
                continue; // coplanar, skip this triangle
            }
        }

        if (t > tol) hits++;
    }
    return (hits % 2 == 1) ? IN : OUT;
}

/// Check if a point is ON the surface of a mesh (within tolerance).
bool point_on_surface_double(const BrepModel& body, const Point3D& p,
                              double tol = 1e-4) {
    const mesh::HalfedgeMesh& mesh = get_mesh(body);
    double tol_sq = tol * tol;

    for (size_t fi = 0; fi < mesh.num_faces(); ++fi) {
        auto& face = mesh.face(fi);
        int h0 = face.halfedge_index;
        int h1 = mesh.halfedge(h0).next_index;
        int h2 = mesh.halfedge(h1).next_index;

        Point3D v0 = mesh.vertex(mesh.halfedge(h0).vertex_index);
        Point3D v1 = mesh.vertex(mesh.halfedge(h1).vertex_index);
        Point3D v2 = mesh.vertex(mesh.halfedge(h2).vertex_index);

        Vector3D e0 = v1 - v0, e1 = v2 - v0, dv = v0 - p;
        double a = e0.dot(e0), b = e0.dot(e1), c = e1.dot(e1);
        double d = e0.dot(dv), ee = e1.dot(dv);
        double det = a * c - b * b;

        double s = b * ee - c * d;
        double t = b * d - a * ee;

        if (s + t <= det) {
            if (s < 0 && t < 0) { if (dv.squaredNorm() < tol_sq) return true; }
            else if (s < 0) { if ((v0 + (t/det)*e1 - p).squaredNorm() < tol_sq) return true; }
            else if (t < 0) { if ((v0 + (s/det)*e0 - p).squaredNorm() < tol_sq) return true; }
            else { if ((v0 + e0*(s/det) + e1*(t/det) - p).squaredNorm() < tol_sq) return true; }
        } else {
            if (s < 0) { if ((v2 - p).squaredNorm() < tol_sq) return true; }
            else if (t < 0) { if ((v1 - p).squaredNorm() < tol_sq) return true; }
            else { /* closest to edge v1-v2 — skip for ON check */ }
        }
    }
    return false;
}

/// Adaptive classification: try double first, upgrade to GMP exact if uncertain.
ClassResult classify_point_adaptive(
    const BrepModel& body, const Point3D& p, int& num_exact_upgrades)
{
    // Step 1: Check ON surface (double precision is adequate here)
    if (point_on_surface_double(body, p)) {
        return ON;
    }

    // Step 2: Try double-precision ray cast
    ClassResult result = classify_point_double(body, p);

    // Step 3: check if the result is uncertain (near boundary)
    // Re-cast with a slightly offset origin to check consistency
    Point3D p2(p.x() + 1e-8, p.y(), p.z());
    ClassResult verify = classify_point_double(body, p2);

    if (result != verify) {
        // Inconsistent: upgrade to exact
        num_exact_upgrades++;
#ifdef VDE_USE_GMP
        // Build polyhedron data for exact_point_in_polyhedron
        const mesh::HalfedgeMesh& mesh = get_mesh(body);
        std::vector<Point3D> verts;
        std::vector<std::vector<int>> faces;

        // Map vertex indices to compact array for exact predicates
        std::map<int, int> vmap;
        for (size_t fi = 0; fi < mesh.num_faces(); ++fi) {
            auto& face = mesh.face(fi);
            int h0 = face.halfedge_index;
            int h1 = mesh.halfedge(h0).next_index;
            int h2 = mesh.halfedge(h1).next_index;
            int v0 = mesh.halfedge(h0).vertex_index;
            int v1 = mesh.halfedge(h1).vertex_index;
            int v2 = mesh.halfedge(h2).vertex_index;

            std::vector<int> fv;
            for (int vi : {v0, v1, v2}) {
                auto it = vmap.find(vi);
                int idx;
                if (it == vmap.end()) {
                    idx = static_cast<int>(verts.size());
                    vmap[vi] = idx;
                    verts.push_back(mesh.vertex(vi));
                } else {
                    idx = it->second;
                }
                fv.push_back(idx);
            }
            faces.push_back(std::move(fv));
        }

        bool inside = core::exact_point_in_polyhedron(p, verts, faces);
        return inside ? IN : OUT;
#else
        // GMP not available: use more samples for verification
        // Cast rays in 4 directions and majority-vote
        Vector3D dirs[] = {
            Vector3D(1,0,0), Vector3D(-1,0,0),
            Vector3D(0,1,0), Vector3D(0,-1,0)
        };
        const mesh::HalfedgeMesh& mesh = get_mesh(body);
        int in_votes = 0;
        for (const auto& d : dirs) {
            int hits = 0;
            for (size_t fi = 0; fi < mesh.num_faces(); ++fi) {
                auto& face = mesh.face(fi);
                int h0 = face.halfedge_index;
                int h1 = mesh.halfedge(h0).next_index;
                int h2 = mesh.halfedge(h1).next_index;
                Point3D v0 = mesh.vertex(mesh.halfedge(h0).vertex_index);
                Point3D v1 = mesh.vertex(mesh.halfedge(h1).vertex_index);
                Point3D v2 = mesh.vertex(mesh.halfedge(h2).vertex_index);
                Vector3D e1 = v1 - v0, e2 = v2 - v0;
                Vector3D h = d.cross(e2);
                double a = e1.dot(h);
                if (std::abs(a) < 1e-8) continue;
                double f = 1.0 / a;
                Vector3D s = p - v0;
                double u = f * s.dot(h);
                if (u < 0 || u > 1) continue;
                Vector3D q = s.cross(e1);
                double v = f * d.dot(q);
                if (v < 0 || u + v > 1) continue;
                double t = f * e2.dot(q);
                if (t > 1e-8) hits++;
            }
            if (hits % 2 == 1) in_votes++;
        }
        return (in_votes >= 3) ? IN : OUT;
#endif
    }

    return result;
}

/// Multi-point sampling classification for a face fragment (public API).
/// Replaces simple centroid classification.
///
/// Algorithm:
/// 1. Compute face bounding box and approximate area
/// 2. Uniformly sample N = sqrt(area / tol^2) points (minimum 4)
/// 3. Each point independently classified via classify_point_adaptive()
/// 4. Vote: if >75% agree, return result
/// 5. If vote is inconclusive, sample 2N more points and re-vote
/// 6. Boundary points (distance < tol) use exact_orient3d for robust determination

}  // anonymous namespace — classify_face_fragment is public

ClassResult classify_face_fragment(
    const BrepModel& body, const BrepModel& other_body,
    int face_id, int& num_exact_upgrades)
{
    // Quick AABB rejection
    if (body.num_faces() == 0) return OUT;
    if (other_body.num_faces() == 0) return OUT;

    AABB3D bb = face_bounds(body, face_id);
    auto other_bbox = other_body.bounds();
    if (!bb.intersects(other_bbox)) return OUT;

    // Compute approximate face area from bounding box
    Vector3D diag = bb.max() - bb.min();
    double approx_area = diag.x() * diag.y() + diag.x() * diag.z() + diag.y() * diag.z();
    if (approx_area < 1e-16) {
        // Degenerate face: fall back to centroid
        Point3D c = bb.center();
        return classify_point_adaptive(other_body, c, num_exact_upgrades);
    }

    double tol = 1e-6;
    int N = static_cast<int>(std::sqrt(approx_area / (tol * tol)));
    if (N < 4) N = 4;
    if (N > 256) N = 256;  // cap to prevent excessive sampling

    const auto& surf = body.surface(body.face(face_id).surface_id);
    double u0 = surf.knots_u().front(), u1 = surf.knots_u().back();
    double v0 = surf.knots_v().front(), v1 = surf.knots_v().back();
    double du = u1 - u0, dv = v1 - v0;
    if (du < 1e-12 || dv < 1e-12) {
        Point3D c = bb.center();
        return classify_point_adaptive(other_body, c, num_exact_upgrades);
    }

    // Generator for uniform (u,v) samples
    std::mt19937 rng(42);  // fixed seed for reproducibility
    std::uniform_real_distribution<double> dist(0.0, 1.0);

    auto sample_and_vote = [&](int n_samples) -> ClassResult {
        int in_count = 0, out_count = 0, on_count = 0;

        for (int i = 0; i < n_samples; ++i) {
            // Uniform sample in parameter space
            double u = u0 + dist(rng) * du;
            double v = v0 + dist(rng) * dv;
            Point3D p = surf.evaluate(u, v);

            // Check if point is near boundary of the other body
            const mesh::HalfedgeMesh& mesh = get_mesh(other_body);
            bool near_boundary = false;
            for (size_t fi = 0; fi < mesh.num_faces() && !near_boundary; ++fi) {
                auto& face = mesh.face(fi);
                int h0 = face.halfedge_index;
                int h1 = mesh.halfedge(h0).next_index;
                int h2 = mesh.halfedge(h1).next_index;
                Point3D v0_m = mesh.vertex(mesh.halfedge(h0).vertex_index);
                Point3D v1_m = mesh.vertex(mesh.halfedge(h1).vertex_index);
                Point3D v2_m = mesh.vertex(mesh.halfedge(h2).vertex_index);
                // Fast distance-to-triangle check using barycentric projection
                Vector3D e0 = v1_m - v0_m, e1 = v2_m - v0_m, dv_vec = v0_m - p;
                double a = e0.dot(e0), b = e0.dot(e1), c = e1.dot(e1);
                double d = e0.dot(dv_vec), e = e1.dot(dv_vec);
                double det = a * c - b * b;
                if (std::abs(det) < 1e-12) continue;
                double s = (b * e - c * d) / det;
                double t = (b * d - a * e) / det;
                if (s >= 0 && t >= 0 && s + t <= 1) {
                    Point3D closest = v0_m + e0 * s + e1 * t;
                    if ((p - closest).norm() < tol * 10.0) {
                        near_boundary = true;
                    }
                }
            }

            if (near_boundary) {
                // Use exact_orient3d for robust boundary classification
                // Cast 3 rays and majority-vote with exact predicates
                Vector3D dirs[] = {
                    Vector3D(1, 0, 0), Vector3D(0, 1, 0), Vector3D(0, 0, 1)
                };
                int in_votes = 0;
                for (const auto& ray_dir : dirs) {
                    int hits = 0;
                    for (size_t fi = 0; fi < mesh.num_faces(); ++fi) {
                        auto& face = mesh.face(fi);
                        int h0 = face.halfedge_index;
                        int h1 = mesh.halfedge(h0).next_index;
                        int h2 = mesh.halfedge(h1).next_index;
                        Point3D v0_m2 = mesh.vertex(mesh.halfedge(h0).vertex_index);
                        Point3D v1_m2 = mesh.vertex(mesh.halfedge(h1).vertex_index);
                        Point3D v2_m2 = mesh.vertex(mesh.halfedge(h2).vertex_index);
                        Vector3D e1 = v1_m2 - v0_m2, e2 = v2_m2 - v0_m2;
                        Vector3D h = ray_dir.cross(e2);
                        double a2 = e1.dot(h);
                        if (std::abs(a2) < 1e-12) continue;
                        double f = 1.0 / a2;
                        Vector3D s2 = p - v0_m2;
                        double u2 = f * s2.dot(h);
                        if (u2 < 0 || u2 > 1) continue;
                        Vector3D q = s2.cross(e1);
                        double v2 = f * ray_dir.dot(q);
                        if (v2 < 0 || u2 + v2 > 1) continue;
                        double t2 = f * e2.dot(q);
                        if (t2 > 1e-12) {
                            // Exact orient3d check for near-coplanar triangles
                            double orient = core::exact_orient3d(
                                v0_m2.x(), v0_m2.y(), v0_m2.z(),
                                v1_m2.x(), v1_m2.y(), v1_m2.z(),
                                v2_m2.x(), v2_m2.y(), v2_m2.z(),
                                p.x(), p.y(), p.z());
                            // If coplanar, classify as boundary intersection
                            if (std::abs(orient) < 1e-12) {
                                hits++;
                            }
                        }
                        if (t2 > tol) hits++;
                    }
                    if (hits % 2 == 1) in_votes++;
                }
                if (in_votes >= 2) in_count++;
                else out_count++;
                continue;
            }

            // Normal classification
            auto cls = classify_point_adaptive(other_body, p, num_exact_upgrades);
            if (cls == IN) in_count++;
            else if (cls == OUT) out_count++;
            else on_count++;
        }

        int total = in_count + out_count + on_count;
        if (total == 0) return ON;

        double in_ratio = static_cast<double>(in_count) / total;
        double out_ratio = static_cast<double>(out_count) / total;
        double on_ratio = static_cast<double>(on_count) / total;

        // ON-dominated: treat as ON
        if (on_ratio > 0.75) return ON;
        // IN-dominated
        if (in_ratio > 0.75) return IN;
        // OUT-dominated
        if (out_ratio > 0.75) return OUT;

        // Inconclusive
        if (in_count > out_count) return IN;
        if (out_count > in_count) return OUT;
        return ON;
    };

    // First pass: N samples
    ClassResult result = sample_and_vote(N);

    // If result is ON or the vote was close, re-sample with 2N points
    if (result == ON) {
        result = sample_and_vote(2 * N);
    }

    return result;
}

// ═══════════════════════════════════════════════════════════
//  Step 4: Sew kept face fragments (internal)
// ═══════════════════════════════════════════════════════════

namespace {  // back to anonymous namespace

/// Sew a collection of face fragments into a unified BrepModel.
/// Deduplicates vertices by proximity, reuses shared edges.
/// Handles singular points: edges with (v_start == v_end) from
/// tangent contacts are recorded but skipped during sewing.
BrepModel sew_face_fragments(const std::vector<BrepModel>& fragments) {
    if (fragments.empty()) return BrepModel();

    BrepModel result;
    std::map<std::pair<int,int>, int> edge_cache;
    std::vector<int> all_face_ids;

    // Spatial hash for vertex dedup
    static constexpr double Q = 1e5;
    auto quantize = [](double v) -> int64_t { return static_cast<int64_t>(v * Q); };
    std::unordered_map<uint64_t, int> pos_to_idx;

    auto hash_pos = [&](const Point3D& p) -> uint64_t {
        uint64_t h = static_cast<uint64_t>(quantize(p.x()));
        h = h * 0x9e3779b9 + static_cast<uint64_t>(quantize(p.y()));
        h = h * 0x9e3779b9 + static_cast<uint64_t>(quantize(p.z()));
        return h;
    };

    // Exact predicate cache: deduplicate with exact_orient3d for near-coincident points
    auto exact_same_point = [&](const Point3D& a, const Point3D& b) -> bool {
        double dist_sq = (a - b).squaredNorm();
        if (dist_sq < 1e-12) return true;
        if (dist_sq > 1e-8) return false;
        // Near threshold: use exact_orient3d with a reference point
        double orient = core::exact_orient3d(
            a.x(), a.y(), a.z(),
            b.x(), b.y(), b.z(),
            a.x()+1.0, a.y(), a.z(),
            a.x(), a.y()+1.0, a.z());
        return std::abs(orient) < 1e-12;
    };

    auto find_or_add_edge = [&](int v0, int v1) -> int {
        // Degenerate edge check: tangent contact produces v0==v1
        if (v0 == v1) return -1; // skip degenerate edges
        auto key = std::make_pair(std::min(v0,v1), std::max(v0,v1));
        auto it = edge_cache.find(key);
        if (it != edge_cache.end()) return it->second;
        int ei = result.add_edge(v0, v1);
        edge_cache[key] = ei;
        return ei;
    };

    for (const auto& fb : fragments) {
        // Merge vertices
        std::map<int, int> old_to_new;
        const double tol_sq = 1e-12;
        for (size_t vi = 0; vi < fb.num_vertices(); ++vi) {
            auto& v = fb.vertex(static_cast<int>(vi));
            uint64_t h = hash_pos(v.point);
            auto it = pos_to_idx.find(h);
            if (it != pos_to_idx.end() &&
                (result.vertex(it->second).point - v.point).squaredNorm() < tol_sq) {
                old_to_new[static_cast<int>(vi)] = it->second;
            } else {
                int new_idx = result.add_vertex(v.point);
                old_to_new[static_cast<int>(vi)] = new_idx;
                pos_to_idx[h] = new_idx;
            }
        }

        // Merge surfaces
        std::map<int, int> surf_map;
        for (size_t fi = 0; fi < fb.num_faces(); ++fi) {
            auto& f = fb.face(static_cast<int>(fi));
            if (surf_map.find(f.surface_id) == surf_map.end()) {
                auto& s = fb.surface(f.surface_id);
                surf_map[f.surface_id] = result.add_surface(s);
            }
        }

        // Copy faces with shared edges (skip degenerate edges from tangent contacts)
        for (size_t fi = 0; fi < fb.num_faces(); ++fi) {
            auto& f = fb.face(static_cast<int>(fi));
            auto fe = fb.face_edges(static_cast<int>(fi));
            std::vector<int> new_edges;
            bool has_degenerate_edge = false;
            for (int ei : fe) {
                auto& e = fb.edge(ei);
                int vs = old_to_new.count(e.v_start) ? old_to_new[e.v_start] : 0;
                int ve = old_to_new.count(e.v_end) ? old_to_new[e.v_end] : 0;
                int new_ei = find_or_add_edge(vs, ve);
                if (new_ei < 0) {
                    has_degenerate_edge = true;
                    continue;
                }
                new_edges.push_back(new_ei);
            }
            // Skip faces that cannot form valid loops due to degenerate edges
            if (new_edges.size() < 2 && has_degenerate_edge) continue;
            if (new_edges.empty()) continue;
            int loop = result.add_loop(new_edges, true);
            int sid = surf_map.count(f.surface_id) ? surf_map[f.surface_id] : 0;
            all_face_ids.push_back(result.add_face(sid, {loop}));
        }
    }

    if (!all_face_ids.empty()) {
        int sh = result.add_shell(all_face_ids, true);
        result.add_body({sh}, "result");
    }
    return result;
}

// ═══════════════════════════════════════════════════════════
//  Core SSI boolean pipeline
// ═══════════════════════════════════════════════════════════

/// Run the full SSI boolean pipeline for a given operation.
/// @param op_type: 0=union, 1=intersection, 2=difference
/// @param keep_predicate: fn(ClassResult, bool is_from_a) → bool to keep
SSIBooleanResult ssi_boolean_core(
    const BrepModel& a, const BrepModel& b, int op_type)
{
    SSIBooleanResult diag;

    // ── Empty body handling ──
    if (a.num_faces() == 0 && b.num_faces() == 0) {
        return diag;  // empty result
    }
    if (a.num_faces() == 0) {
        if (op_type == 0) { diag.result = b; return diag; }  // union
        return diag;  // intersection/difference: empty
    }
    if (b.num_faces() == 0) {
        if (op_type == 0) { diag.result = a; return diag; }  // union
        if (op_type == 2) { diag.result = a; return diag; }  // difference A\∅=A
        return diag;  // intersection: empty
    }

    // ── Quick AABB disjoint check ──
    if (!a.bounds().intersects(b.bounds())) {
        if (op_type == 0) {
            // Union of disjoint bodies: sew both
            diag.result = sew_face_fragments({a, b});
        } else if (op_type == 2) {
            // Difference: A stays
            diag.result = a;
        }
        // Intersection of disjoint: empty (already default)
        return diag;
    }

    // ═══════════════════════════════════════════════
    //  Degeneracy detection: 7-category check (P0+P1 robustness)
    //  Each detected degeneracy is recorded in the precision tracker
    // ═══════════════════════════════════════════════
    bool has_degeneracy = false;
    {
        // Category 1 & 3: Surface overlap + high-order contact detection
        for (size_t i = 0; i < a.num_faces() && !has_degeneracy; ++i) {
            for (size_t j = 0; j < b.num_faces() && !has_degeneracy; ++j) {
                int fi = static_cast<int>(i), fj = static_cast<int>(j);

                auto overlap = detect_surface_overlap(a, b, fi, fj);
                if (overlap.overlapping) {
                    has_degeneracy = true;
                    diag.precision_tracker.record("surface_overlap", 1e-8);
                    break;
                }

                auto hoc = detect_high_order_contact(a, b, fi, fj);
                if (hoc.has_contact) {
                    has_degeneracy = true;
                    diag.precision_tracker.record("high_order_contact", 1e-7);
                    break;
                }
            }
        }

        // Category 6: Non-manifold edge detection
        for (size_t ei = 0; ei < a.num_edges() && !has_degeneracy; ++ei) {
            auto nm = detect_non_manifold_contact(a, static_cast<int>(ei));
            if (nm.is_non_manifold) {
                has_degeneracy = true;
                diag.precision_tracker.record("non_manifold_contact", 1e-8);
            }
        }
        for (size_t ei = 0; ei < b.num_edges() && !has_degeneracy; ++ei) {
            auto nm = detect_non_manifold_contact(b, static_cast<int>(ei));
            if (nm.is_non_manifold) {
                has_degeneracy = true;
                diag.precision_tracker.record("non_manifold_contact", 1e-8);
            }
        }
    }
    // Record overall degeneracy presence
    if (has_degeneracy) {
        diag.num_exact_upgrades += 1;
    }

    // Category 4 & 5: Will be checked per-face-pair during SSI computation

    // ═══════════════════════════════════════════════
    //  Step 1: SSI computation
    // ═══════════════════════════════════════════════
    auto ssi_t0 = std::chrono::steady_clock::now();
    auto ssi_curves = compute_all_ssi(a, b, diag.num_exact_upgrades);
    auto ssi_t1 = std::chrono::steady_clock::now();
    diag.ssi_time_ms = std::chrono::duration<double, std::milli>(ssi_t1 - ssi_t0).count();
    diag.num_ssi_curves = static_cast<int>(ssi_curves.size());
    diag.precision_tracker.record("ssi_compute", 1e-9 * ssi_curves.size());

    // Category 2 & 4 & 5: Check SSI curves for overlap/boundary/vertex contact
    for (const auto& c : ssi_curves) {
        if (!c.valid) continue;

        // Category 4: Boundary contact detection
        auto bc = detect_boundary_contact(c.points, a, b, c.face_a, c.face_b);
        if (bc.has_boundary_contact) {
            diag.precision_tracker.record("boundary_contact", 1e-8);
        }

        // Category 5: Vertex contact detection
        auto vc = detect_vertex_contact(c.points, a, b, c.face_a, c.face_b);
        if (vc.has_vertex_contact) {
            diag.precision_tracker.record("vertex_contact", 1e-8);
        }
    }

    // Category 2: Check for overlapping SSI curves
    for (size_t ci = 0; ci < ssi_curves.size(); ++ci) {
        for (size_t cj = ci + 1; cj < ssi_curves.size(); ++cj) {
            auto cov = detect_curve_overlap(ssi_curves[ci].points, ssi_curves[cj].points);
            if (cov.overlapping) {
                diag.precision_tracker.record("curve_overlap", 1e-8);
            }
        }
    }

    // ═══════════════════════════════════════════════
    //  Step 2: Split faces using SSI cutting planes
    // ═══════════════════════════════════════════════
    auto split_t0 = std::chrono::steady_clock::now();

    auto a_cut_planes = build_face_cutting_planes(ssi_curves, true);
    auto b_cut_planes = build_face_cutting_planes(ssi_curves, false);

    auto a_fragments = split_and_classify_faces_ssi(a, a_cut_planes, b, diag.num_exact_upgrades);
    auto b_fragments = split_and_classify_faces_ssi(b, b_cut_planes, a, diag.num_exact_upgrades);

    auto split_t1 = std::chrono::steady_clock::now();
    diag.split_time_ms = std::chrono::duration<double, std::milli>(split_t1 - split_t0).count();
    diag.num_fragments = static_cast<int>(a_fragments.size() + b_fragments.size());
    diag.precision_tracker.record("face_split", 1e-8 * diag.num_fragments);

    // ═══════════════════════════════════════════════
    //  Step 3: Classification statistics
    // ═══════════════════════════════════════════════
    auto classify_t0 = std::chrono::steady_clock::now();

    // Classification already happened inside split_and_classify_faces_ssi.
    // Now apply boolean operation filter.
    std::vector<BrepModel> keep_fragments;

    // Predicate: which fragments to keep based on operation type
    auto should_keep = [op_type](ClassResult cls, bool from_a) -> bool {
        switch (op_type) {
            case 0: // union
                // Keep A fragments OUT/ON of B; keep B fragments OUT of A
                if (from_a) return (cls == OUT || cls == ON);
                else        return (cls == OUT);
            case 1: // intersection
                // Keep A fragments IN/ON of B; keep B fragments IN of A
                if (from_a) return (cls == IN || cls == ON);
                else        return (cls == IN);
            case 2: // difference A \ B
                // Keep A fragments OUT of B; keep B fragments IN of A (inner surface)
                if (from_a) return (cls == OUT);
                else        return (cls == IN);
            default: return false;
        }
    };

    for (auto& [frag, cls] : a_fragments) {
        if (cls == IN) diag.num_classified_in++;
        else if (cls == OUT) diag.num_classified_out++;
        else diag.num_classified_on++;

        if (should_keep(cls, true)) {
            keep_fragments.push_back(std::move(frag));
            diag.num_kept++;
        }
    }
    for (auto& [frag, cls] : b_fragments) {
        if (cls == IN) diag.num_classified_in++;
        else if (cls == OUT) diag.num_classified_out++;
        else diag.num_classified_on++;

        if (should_keep(cls, false)) {
            keep_fragments.push_back(std::move(frag));
            diag.num_kept++;
        }
    }

    auto classify_t1 = std::chrono::steady_clock::now();
    diag.classify_time_ms = std::chrono::duration<double, std::milli>(classify_t1 - classify_t0).count();
    diag.precision_tracker.record("classify", 1e-9 * diag.num_fragments);

    // ═══════════════════════════════════════════════
    //  Step 4: Sew kept fragments + heal orientation
    // ═══════════════════════════════════════════════
    auto sew_t0 = std::chrono::steady_clock::now();

    if (!keep_fragments.empty()) {
        diag.result = sew_face_fragments(keep_fragments);

        // Heal orientation for consistent face normals
        if (diag.result.num_faces() > 0) {
            heal_orientation(diag.result);
        }
    }

    auto sew_t1 = std::chrono::steady_clock::now();
    diag.sew_time_ms = std::chrono::duration<double, std::milli>(sew_t1 - sew_t0).count();
    diag.precision_tracker.record("sew", 1e-9 * diag.num_kept);

    // Check accumulated precision loss against threshold
    // If excessive, record warning in diagnostic (caller can check)
    if (diag.precision_tracker.exceeds_threshold(1e-4)) {
        diag.num_exact_upgrades += 1;  // flag excessive precision loss
    }

    return diag;
}

}  // anonymous namespace

// ═══════════════════════════════════════════════════════════
//  Public API
// ═══════════════════════════════════════════════════════════

SSIBooleanResult ssi_boolean_union(const BrepModel& a, const BrepModel& b) {
    return ssi_boolean_core(a, b, 0);
}

SSIBooleanResult ssi_boolean_intersection(const BrepModel& a, const BrepModel& b) {
    return ssi_boolean_core(a, b, 1);
}

SSIBooleanResult ssi_boolean_difference(const BrepModel& a, const BrepModel& b) {
    return ssi_boolean_core(a, b, 2);
}

}  // namespace vde::brep