#include "vde/brep/brep_boolean.h"
#include "vde/brep/brep_face_split.h"
#include "vde/core/plane.h"
#include "vde/core/line.h"
#include "vde/core/aabb.h"
#include "vde/mesh/halfedge_mesh.h"
#include "vde/curves/surface_intersection.h"
#include <algorithm>
#include <set>
#include <map>
#include <unordered_map>
#include <cmath>
#include <limits>
#include <mutex>

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

namespace vde::brep {

using core::Point3D;
using core::Vector3D;
using core::Plane3D;
using core::Line3Dd;
using core::Ray3Dd;
using core::AABB3D;

namespace {

// ── Ray-casting classification ──────────────────────────
enum ClassResult { IN = -1, ON = 0, OUT = 1 };

// ── Mesh caching ─────────────────────────────────────────
/// Cache to_mesh() results keyed by BrepModel pointer.
/// Avoids re-tessellating the same body for every face classification.
static std::unordered_map<const BrepModel*, mesh::HalfedgeMesh> mesh_cache;

/// Get the tessellated mesh for a body, reusing cached result.
static const mesh::HalfedgeMesh& get_mesh(const BrepModel& body) {
    auto it = mesh_cache.find(&body);
    if (it != mesh_cache.end()) return it->second;
    auto [inserted_it, _] = mesh_cache.emplace(&body, body.to_mesh(0.05));
    return inserted_it->second;
}

/// Minimum signed distance from point p to the mesh.
/// Also returns the ray-cast in/out classification.
ClassResult classify_point_mesh(const BrepModel& body, const Point3D& p,
                                  double tol = 1e-6) {
    const mesh::HalfedgeMesh& mesh = get_mesh(body);
    if (mesh.num_faces() == 0) return OUT;

    // ── First: check if point lies ON the mesh surface ──
    double min_dist_sq = std::numeric_limits<double>::max();
    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);

        // Compute closest point on triangle to p
        Vector3D e0 = v1 - v0;
        Vector3D e1 = v2 - v0;
        Vector3D dv = v0 - p;

        double a = e0.dot(e0);
        double b = e0.dot(e1);
        double c = e1.dot(e1);
        double d = e0.dot(dv);
        double e = e1.dot(dv);

        double det = a * c - b * b;
        double s = b * e - c * d;
        double t = b * d - a * e;

        if (s + t <= det) {
            if (s < 0) {
                if (t < 0) {
                    // Region 4: closest to v0
                    double dist = dv.squaredNorm();
                    min_dist_sq = std::min(min_dist_sq, dist);
                } else {
                    // Region 3: closest to edge v0-v2
                    double dist = (v0 + (t/det) * e1 - p).squaredNorm();
                    min_dist_sq = std::min(min_dist_sq, dist);
                }
            } else if (t < 0) {
                // Region 5: closest to edge v0-v1
                double dist = (v0 + (s/det) * e0 - p).squaredNorm();
                min_dist_sq = std::min(min_dist_sq, dist);
            } else {
                // Region 0: inside triangle
                double dist = (v0 + e0 * (s/det) + e1 * (t/det) - p).squaredNorm();
                min_dist_sq = std::min(min_dist_sq, dist);
            }
        } else {
            if (s < 0) {
                // Region 2: closest to v2
                double dist = (v2 - p).squaredNorm();
                min_dist_sq = std::min(min_dist_sq, dist);
            } else if (t < 0) {
                // Region 6: closest to v1
                double dist = (v1 - p).squaredNorm();
                min_dist_sq = std::min(min_dist_sq, dist);
            } else {
                // Region 1: closest to edge v1-v2
                double numer = a + d - b - e;
                double denom = a - 2*b + c;
                double w = (denom > 1e-12) ? std::clamp(numer / denom, 0.0, 1.0) : 0.0;
                double dist = (v1 + w * (v2 - v1) - p).squaredNorm();
                min_dist_sq = std::min(min_dist_sq, dist);
            }
        }
    }

    // If very close to surface, classify as ON
    double on_threshold = 1e-4;  // 0.1mm
    if (min_dist_sq < on_threshold * on_threshold) {
        return ON;
    }

    // ── Ray casting for IN/OUT ──
    // Cast ray in +X direction and count intersections
    Ray3Dd ray(p, Vector3D::UnitX());
    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);

        // Möller–Trumbore ray-triangle intersection
        Vector3D e1 = v1 - v0;
        Vector3D e2 = v2 - v0;
        Vector3D h = ray.direction().cross(e2);
        double a = e1.dot(h);

        if (std::abs(a) < tol) continue;
        double f = 1.0 / a;
        Vector3D s = ray.origin() - v0;
        double u = f * s.dot(h);

        if (u < 0.0 || u > 1.0) continue;

        Vector3D q = s.cross(e1);
        double v = f * ray.direction().dot(q);

        if (v < 0.0 || u + v > 1.0) continue;

        double t = f * e2.dot(q);
        if (t > tol) hits++;
    }

    return (hits % 2 == 1) ? IN : OUT;
}

/// Classify a face fragment relative to another body.
/// Samples the face at its centroid and uses ray casting.
ClassResult classify_face_fragment(const BrepModel& face_body,
                                    const BrepModel& other_body) {
    // Quick AABB rejection: if fragment doesn't intersect other body, it's outside
    AABB3D bb;
    for (size_t vi = 0; vi < face_body.num_vertices(); ++vi)
        bb.expand(face_body.vertex(static_cast<int>(vi)).point);

    auto other_bbox = other_body.bounds();
    if (!bb.intersects(other_bbox)) return OUT;

    Point3D centroid = bb.center();
    return classify_point_mesh(other_body, centroid);
}

// ── Plane-plane intersection ────────────────────────────
/// Compute intersection line between two planes.
/// Returns true if planes intersect (non-parallel).
bool intersect_plane_plane(const Plane3D& p1, const Plane3D& p2,
                           Line3Dd& out_line) {
    Vector3D n1 = p1.normal();
    Vector3D n2 = p2.normal();
    Vector3D dir = n1.cross(n2);
    double len_sq = dir.squaredNorm();

    if (len_sq < 1e-12) return false; // parallel

    dir.normalize();

    // Find a point on the intersection line by solving:
    // n1·p = -d1,  n2·p = -d2,  dir·p = 0
    // Use Cramer's rule on the 3x3 system
    double d1 = p1.d();
    double d2 = p2.d();

    Eigen::Matrix3d M;
    M.row(0) = n1.transpose();
    M.row(1) = n2.transpose();
    M.row(2) = dir.transpose();

    Vector3D rhs(-d1, -d2, 0.0);
    Vector3D p = M.inverse() * rhs;

    out_line = Line3Dd(p, dir);
    return true;
}

// ── Get face plane ──────────────────────────────────────
Plane3D face_plane(const BrepModel& body, int face_id) {
    const auto& face = body.face(face_id);
    const auto& surf = body.surface(face.surface_id);
    Point3D p = surf.evaluate(0.5, 0.5);
    Vector3D n = surf.normal(0.5, 0.5);
    if (n.norm() < 1e-12) n = Vector3D::UnitZ();
    return Plane3D(p, n);
}

// ── Face normal from 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;
}

// ── Compute face centroid ───────────────────────────────
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;
}

// ── Face-face intersection test ─────────────────────────
bool faces_intersect(const BrepModel& body_a, int face_a,
                      const BrepModel& body_b, int face_b) {
    // Quick AABB test
    AABB3D bb_a, bb_b;
    auto ea = body_a.face_edges(face_a);
    for (int ei : ea) {
        auto& e = body_a.edge(ei);
        bb_a.expand(body_a.vertex(e.v_start).point);
        bb_a.expand(body_a.vertex(e.v_end).point);
    }
    auto eb = body_b.face_edges(face_b);
    for (int ei : eb) {
        auto& e = body_b.edge(ei);
        bb_b.expand(body_b.vertex(e.v_start).point);
        bb_b.expand(body_b.vertex(e.v_end).point);
    }
    return bb_a.intersects(bb_b);
}

// ── Find or create edge between two vertices ────────────
/// Returns the edge ID for an edge between v0 and v1.
/// Reuses an existing edge if one already exists, otherwise creates new.
int find_or_add_edge(BrepModel& result, int v0, int v1,
                     std::map<std::pair<int,int>, int>& edge_cache) {
    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;
}

// ── Sew faces into BrepModel ────────────────────────────
BrepModel sew_faces(const std::vector<BrepModel>& face_bodies) {
    BrepModel result;
    std::vector<int> all_face_ids;
    // Cache: (min_vertex, max_vertex) → edge_id for edge sharing
    std::map<std::pair<int,int>, int> edge_cache;

    for (auto& fb : face_bodies) {
        // Merge vertices (deduplicate by proximity)
        std::map<int, int> old_to_new;  // fragment vertex ID → result vertex index
        for (size_t vi = 0; vi < fb.num_vertices(); ++vi) {
            auto& v = fb.vertex(static_cast<int>(vi));
            int best_idx = -1;
            double best_dist = 1e-6;
            for (size_t ri = 0; ri < result.num_vertices(); ++ri) {
                double d = (result.vertex(static_cast<int>(ri)).point - v.point).norm();
                if (d < best_dist) { best_dist = d; best_idx = static_cast<int>(ri); }
            }
            if (best_idx >= 0) {
                old_to_new[static_cast<int>(vi)] = best_idx;
            } else {
                old_to_new[static_cast<int>(vi)] = result.add_vertex(v.point);
            }
        }

        // 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
        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;
            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;
                // Use shared edges when vertices match
                new_edges.push_back(find_or_add_edge(result, vs, ve, edge_cache));
            }
            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;
}

// ── Extract a face as a separate 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);

    // Map original vertex indices to new vertex indices in fragment
    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);
        }
    }

    // Add the surface
    int sid = frag.add_surface(body.surface(f.surface_id));

    // Add edges and loop
    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);

    std::string name = "face_" + std::to_string(face_id);
    int nf = frag.add_face(sid, {loop});
    int sh = frag.add_shell({nf}, false);
    frag.add_body({sh}, name);
    return frag;
}

// ── Split and classify faces by face-plane cutting ─────
/// Split each face of body_a by the planes of body_b's faces,
/// then classify each resulting fragment against body_b.
/// Returns vector of (fragment, classification) pairs.
std::vector<std::pair<BrepModel, ClassResult>>
split_and_classify_faces(const BrepModel& a, const BrepModel& b) {
    std::vector<std::pair<BrepModel, ClassResult>> result;

    for (size_t fi = 0; fi < a.num_faces(); ++fi) {
        std::vector<BrepModel> current;
        current.push_back(extract_face_fragment(a, static_cast<int>(fi)));

        // Split by each face plane of B
        for (size_t bfi = 0; bfi < b.num_faces(); ++bfi) {
            Point3D plane_pt = face_centroid(b, static_cast<int>(bfi));
            Vector3D plane_n = face_normal(b, static_cast<int>(bfi));

            std::vector<BrepModel> next;
            for (auto& frag : current) {
                auto splits = split_face_by_plane(frag, 0, plane_pt, plane_n);
                for (auto& s : splits) next.push_back(std::move(s));
            }
            current = std::move(next);
            if (current.empty()) break;  // all fragments eliminated
        }

        // Classify each remaining fragment
        for (auto& frag : current) {
            auto cls = classify_face_fragment(frag, b);
            result.emplace_back(std::move(frag), cls);
        }
    }
    return result;
}

/// Create a deep copy of a BrepModel, preserving all topology exactly.
BrepModel copy_brep(const BrepModel& src) {
    BrepModel result;
    // Copy vertices (key by array index, since e.v_start/e.v_end store indices)
    std::map<int, int> old_vid_to_new_idx;
    for (size_t vi = 0; vi < src.num_vertices(); ++vi) {
        auto& v = src.vertex(static_cast<int>(vi));
        old_vid_to_new_idx[static_cast<int>(vi)] = result.add_vertex(v.point);
    }
    // Copy surfaces
    std::map<int, int> old_sid_to_new_idx;
    for (size_t si = 0; si < src.num_surfaces(); ++si) {
        old_sid_to_new_idx[static_cast<int>(si)] = result.add_surface(src.surface(static_cast<int>(si)));
    }
    // Copy faces
    std::vector<int> all_face_ids;
    for (size_t fi = 0; fi < src.num_faces(); ++fi) {
        auto& f = src.face(static_cast<int>(fi));
        auto es = src.face_edges(static_cast<int>(fi));
        std::vector<int> new_edges;
        for (int ei : es) {
            auto& e = src.edge(ei);
            int vs = old_vid_to_new_idx.at(e.v_start);
            int ve = old_vid_to_new_idx.at(e.v_end);
            new_edges.push_back(result.add_edge(vs, ve));
        }
        int loop = result.add_loop(new_edges, true);
        int sid = old_sid_to_new_idx.at(f.surface_id);
        all_face_ids.push_back(result.add_face(sid, {loop}));
    }
    // Copy body
    if (!all_face_ids.empty()) {
        int sh = result.add_shell(all_face_ids, true);
        result.add_body({sh}, "copy");
    }
    return result;
}

// ── Face bounding box (from edges/vertices) ────────────
core::AABB3D face_bounds(const BrepModel& body, int face_id) {
    core::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;
}

// ── SSI-based face splitting ────────────────────────────
/// Split a face fragment by its surface-surface intersection with another face.
/// Uses intersect_surfaces() to detect and locate the intersection,
/// then creates a cutting plane from the SSI points for clean face splitting.
/// Falls back to returning the fragment unchanged when SSI fails or no intersection exists.
///
/// @param frag_a Fragment BrepModel (assumed to have exactly 1 face at index 0)
/// @param body_b The other body whose face we're intersecting against
/// @param face_b_id Face index within body_b
std::vector<BrepModel> split_face_by_face_ssi(
    const BrepModel& frag_a,
    const BrepModel& body_b, int face_b_id) {

    if (frag_a.num_faces() == 0) return {};

    const auto& surf_a = frag_a.surface(frag_a.face(0).surface_id);
    const auto& surf_b = body_b.surface(body_b.face(face_b_id).surface_id);

    // ── Guard 1: Quick AABB reject ──
    AABB3D bb_a;
    for (size_t vi = 0; vi < frag_a.num_vertices(); ++vi)
        bb_a.expand(frag_a.vertex(static_cast<int>(vi)).point);
    auto bb_b = face_bounds(body_b, face_b_id);
    if (!bb_a.intersects(bb_b)) {
        return {frag_a};
    }

    // ── Guard 2: Compute SSI ──
    auto isect_curves = curves::intersect_surfaces(surf_a, surf_b, 4, 1e-3);

    // Find the first valid intersection curve with enough points
    const curves::IntersectionCurve* valid_curve = nullptr;
    for (const auto& curve : isect_curves) {
        if (curve.points.size() >= 4) {
            valid_curve = &curve;
            break;
        }
    }

    if (!valid_curve) {
        // Guard 3: No SSI intersection → return fragment unchanged
        // Centroid classification will determine IN/OUT later
        return {frag_a};
    }

    // ── SSI confirmed intersection: fit cutting plane from SSI points ──
    Point3D centroid(0, 0, 0);
    for (const auto& p : valid_curve->points) centroid += p;
    centroid /= static_cast<double>(valid_curve->points.size());

    // Compute plane normal: use cross product of two spread directions, fall back to B face normal
    Vector3D plane_n(0,0,0);
    if (valid_curve->points.size() >= 3) {
        Vector3D d1 = valid_curve->points[1] - valid_curve->points[0];
        Vector3D d2 = valid_curve->points.back() - valid_curve->points[0];
        plane_n = d1.cross(d2);
    }
    if (plane_n.norm() < 1e-12) {
        plane_n = face_normal(body_b, face_b_id);
    }

    // Split the fragment by the SSI-derived plane
    auto split_result = split_face_by_plane(frag_a, 0, centroid, plane_n);

    if (split_result.empty()) {
        return {frag_a};
    }

    return split_result;
}

}  // anonymous namespace

// ── Public API ──────────────────────────────────────────

BrepModel brep_union(const BrepModel& a, const BrepModel& b) {
    // Empty body handling
    if (a.num_faces() == 0) return b;
    if (b.num_faces() == 0) return a;

    // Quick AABB check for disjoint bodies
    if (!a.bounds().intersects(b.bounds())) {
        // Disjoint: merge both bodies
        return sew_faces({a, b});
    }

    // ── v3.6 SSI-based approach ──
    // For each face of A: split against each face of B using SSI,
    // then classify fragments. Keep OUT and ON (from A).
    // For each face of B: same against A, keep only OUT.
    std::vector<BrepModel> keep;
    std::mutex keep_mutex;

    // A faces: SSI-split against B, classify, keep OUT/ON
#ifdef VDE_USE_OPENMP
#pragma omp parallel for schedule(dynamic)
#endif
    for (size_t fia = 0; fia < a.num_faces(); ++fia) {
        std::vector<BrepModel> fragments;
        fragments.push_back(extract_face_fragment(a, static_cast<int>(fia)));

        // Split by SSI against each face of B
        for (size_t fib = 0; fib < b.num_faces(); ++fib) {
            std::vector<BrepModel> new_fragments;
            for (auto& frag : fragments) {
                auto splits = split_face_by_face_ssi(
                    frag, b, static_cast<int>(fib));
                if (splits.size() > 1) {
                    // SSI confirmed intersection → use split fragments
                    for (auto& s : splits) new_fragments.push_back(std::move(s));
                } else {
                    // No SSI intersection or single fragment → keep as-is
                    new_fragments.push_back(std::move(frag));
                }
            }
            fragments = std::move(new_fragments);
            if (fragments.empty()) break;
        }

        // Classify each fragment against B
        for (auto& frag : fragments) {
            auto cls = classify_face_fragment(frag, b);
            if (cls == OUT || cls == ON) {
                std::lock_guard<std::mutex> lock(keep_mutex);
                keep.push_back(std::move(frag));
            }
        }
    }

    // B faces: SSI-split against A, classify, keep OUT only (ON from A already)
#ifdef VDE_USE_OPENMP
#pragma omp parallel for schedule(dynamic)
#endif
    for (size_t fib = 0; fib < b.num_faces(); ++fib) {
        std::vector<BrepModel> fragments;
        fragments.push_back(extract_face_fragment(b, static_cast<int>(fib)));

        for (size_t fia = 0; fia < a.num_faces(); ++fia) {
            std::vector<BrepModel> new_fragments;
            for (auto& frag : fragments) {
                auto splits = split_face_by_face_ssi(
                    frag, a, static_cast<int>(fia));
                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;
        }

        for (auto& frag : fragments) {
            auto cls = classify_face_fragment(frag, a);
            if (cls == OUT) {
                std::lock_guard<std::mutex> lock(keep_mutex);
                keep.push_back(std::move(frag));
            }
        }
    }

    if (keep.empty()) return BrepModel();
    return sew_faces(keep);
}

BrepModel brep_intersection(const BrepModel& a, const BrepModel& b) {
    // Empty body handling
    if (a.num_faces() == 0 || b.num_faces() == 0) return BrepModel();

    // Quick AABB check for disjoint bodies
    if (!a.bounds().intersects(b.bounds())) return BrepModel();

    // ── v3.6 SSI-based approach ──
    // For intersection: keep fragments IN the other body,
    // and fragments ON the boundary (only from A to avoid duplicates).
    std::vector<BrepModel> keep;
    std::mutex keep_mutex;

    // A faces: SSI-split against B, classify, keep IN/ON
#ifdef VDE_USE_OPENMP
#pragma omp parallel for schedule(dynamic)
#endif
    for (size_t fia = 0; fia < a.num_faces(); ++fia) {
        std::vector<BrepModel> fragments;
        fragments.push_back(extract_face_fragment(a, static_cast<int>(fia)));

        for (size_t fib = 0; fib < b.num_faces(); ++fib) {
            std::vector<BrepModel> new_fragments;
            for (auto& frag : fragments) {
                auto splits = split_face_by_face_ssi(
                    frag, b, static_cast<int>(fib));
                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;
        }

        for (auto& frag : fragments) {
            auto cls = classify_face_fragment(frag, b);
            if (cls == IN || cls == ON) {
                std::lock_guard<std::mutex> lock(keep_mutex);
                keep.push_back(std::move(frag));
            }
        }
    }

    // B faces: SSI-split against A, classify, keep IN only
#ifdef VDE_USE_OPENMP
#pragma omp parallel for schedule(dynamic)
#endif
    for (size_t fib = 0; fib < b.num_faces(); ++fib) {
        std::vector<BrepModel> fragments;
        fragments.push_back(extract_face_fragment(b, static_cast<int>(fib)));

        for (size_t fia = 0; fia < a.num_faces(); ++fia) {
            std::vector<BrepModel> new_fragments;
            for (auto& frag : fragments) {
                auto splits = split_face_by_face_ssi(
                    frag, a, static_cast<int>(fia));
                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;
        }

        for (auto& frag : fragments) {
            auto cls = classify_face_fragment(frag, a);
            if (cls == IN) {
                std::lock_guard<std::mutex> lock(keep_mutex);
                keep.push_back(std::move(frag));
            }
        }
    }

    if (keep.empty()) return BrepModel();
    return sew_faces(keep);
}

BrepModel brep_difference(const BrepModel& a, const BrepModel& b) {
    // Empty body handling
    if (a.num_faces() == 0) return BrepModel();
    if (b.num_faces() == 0) return copy_brep(a);

    // Quick AABB check for disjoint bodies
    if (!a.bounds().intersects(b.bounds())) return copy_brep(a);

    // ── v3.6 SSI-based approach ──
    // For difference A \ B:
    // - Keep A fragments that are OUT of B (discard IN and ON)
    // - Keep B fragments that are IN A (forms the inner cut surface)
    std::vector<BrepModel> keep;
    std::mutex keep_mutex;

    // A faces: SSI-split against B, classify, keep OUT
#ifdef VDE_USE_OPENMP
#pragma omp parallel for schedule(dynamic)
#endif
    for (size_t fia = 0; fia < a.num_faces(); ++fia) {
        std::vector<BrepModel> fragments;
        fragments.push_back(extract_face_fragment(a, static_cast<int>(fia)));

        for (size_t fib = 0; fib < b.num_faces(); ++fib) {
            std::vector<BrepModel> new_fragments;
            for (auto& frag : fragments) {
                auto splits = split_face_by_face_ssi(
                    frag, b, static_cast<int>(fib));
                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;
        }

        for (auto& frag : fragments) {
            auto cls = classify_face_fragment(frag, b);
            if (cls == OUT) {
                std::lock_guard<std::mutex> lock(keep_mutex);
                keep.push_back(std::move(frag));
            }
        }
    }

    // B faces: SSI-split against A, classify, keep IN (forms inner cut surface)
#ifdef VDE_USE_OPENMP
#pragma omp parallel for schedule(dynamic)
#endif
    for (size_t fib = 0; fib < b.num_faces(); ++fib) {
        std::vector<BrepModel> fragments;
        fragments.push_back(extract_face_fragment(b, static_cast<int>(fib)));

        for (size_t fia = 0; fia < a.num_faces(); ++fia) {
            std::vector<BrepModel> new_fragments;
            for (auto& frag : fragments) {
                auto splits = split_face_by_face_ssi(
                    frag, a, static_cast<int>(fia));
                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;
        }

        for (auto& frag : fragments) {
            auto cls = classify_face_fragment(frag, a);
            if (cls == IN) {
                std::lock_guard<std::mutex> lock(keep_mutex);
                keep.push_back(std::move(frag));
            }
        }
    }

    if (keep.empty()) return BrepModel();
    return sew_faces(keep);
}

}  // namespace vde::brep