"""
OpenCASCADE-based geometry kernel for Fluency CAD.

This module provides a concrete implementation of the geometry kernel
using OCP (OpenCASCADE Python bindings).
"""

import logging
from typing import List, Tuple, Optional, Any, Dict
import numpy as np

from fluency.geometry.base import (
    GeometryKernel,
    GeometryObject,
    Point2D,
    Point3D,
)

logger = logging.getLogger(__name__)


class OCCGeometryObject(GeometryObject):
    """Geometry object wrapper for OpenCASCADE shapes."""

    def __init__(self, shape: Any = None, metadata: Optional[Dict] = None):
        super().__init__(shape, metadata)


class OCGeometryKernel(GeometryKernel):
    """
    OpenCASCADE-based geometry kernel implementation.

    This kernel uses OCP (OpenCASCADE Python bindings) for all geometry
    operations.
    """

    def __init__(self) -> None:
        self._tolerance: float = 0.001
        self._mesh_tolerance: float = 0.1

    def _get_shape(self, obj: GeometryObject) -> Any:
        """Extract the underlying OCC shape from a GeometryObject.

        Returns *None* if the object carries no shape (e.g. an empty sketch) —
        callers should check for None before using the result.
        """
        if isinstance(obj, OCCGeometryObject):
            if obj.shape is not None and hasattr(obj.shape, "wrapped"):
                return obj.shape.wrapped
            return obj.shape
        # Non-OCCGeometryObject: return its shape if present, else None.
        # (Use explicit identity/truth checks — some OCP TopoDS objects have a
        # falsy __bool__, so ``obj.shape if obj.shape`` is unsafe.)
        shape = getattr(obj, "shape", None)
        return shape if shape is not None else None

    def create_point(self, x: float, y: float) -> GeometryObject:
        """Create a 2D point."""
        from OCP.gp import gp_Pnt

        return OCCGeometryObject(gp_Pnt(x, y, 0))

    def create_line(self, start: Point2D, end: Point2D) -> GeometryObject:
        """Create a 2D line segment."""
        from OCP.gp import gp_Pnt
        from OCP.BRepBuilderAPI import BRepBuilderAPI_MakeEdge

        edge = BRepBuilderAPI_MakeEdge(gp_Pnt(start.x, start.y, 0), gp_Pnt(end.x, end.y, 0)).Edge()
        return OCCGeometryObject(edge, {"type": "line"})

    def create_circle(self, center: Point2D, radius: float) -> GeometryObject:
        """Create a 2D circle."""
        from OCP.gp import gp_Pnt, gp_Dir, gp_Ax2, gp_Circ
        from OCP.BRepBuilderAPI import BRepBuilderAPI_MakeEdge

        circ = gp_Circ(gp_Ax2(gp_Pnt(center.x, center.y, 0), gp_Dir(0, 0, 1)), radius)
        edge = BRepBuilderAPI_MakeEdge(circ).Edge()
        return OCCGeometryObject(edge, {"type": "circle"})

    def create_arc(
        self, center: Point2D, radius: float , start_angle: float, end_angle: float
    ) -> GeometryObject:
        """Create a 2D arc."""
        import math
        from OCP.gp import gp_Pnt, gp_Dir, gp_Ax2, gp_Circ
        from OCP.BRepBuilderAPI import BRepBuilderAPI_MakeEdge

        start_rad = math.radians(start_angle)
        end_rad = math.radians(end_angle)

        circ = gp_Circ(gp_Ax2(gp_Pnt(center.x, center.y, 0), gp_Dir(0, 0, 1)), radius)
        edge = BRepBuilderAPI_MakeEdge(circ, start_rad, end_rad).Edge()
        return OCCGeometryObject(edge, {"type": "arc"})

    def create_polygon(self, points: List[Point2D]) -> GeometryObject:
        """Create a closed polygon from points."""
        from OCP.gp import gp_Pnt
        from OCP.BRepBuilderAPI import BRepBuilderAPI_MakePolygon

        if len(points) < 3:
            raise ValueError("Polygon requires at least 3 points")

        mp = BRepBuilderAPI_MakePolygon()
        for pt in points:
            mp.Add(gp_Pnt(pt.x, pt.y, 0))
        mp.Close()
        return OCCGeometryObject(mp.Wire(), {"type": "polygon"})

    def create_rectangle(
        self, width: float, height: float, center: Optional[Point2D] = None
    ) -> GeometryObject:
        """Create a rectangle."""
        from OCP.gp import gp_Pnt
        from OCP.BRepBuilderAPI import BRepBuilderAPI_MakePolygon

        cx = center.x if center else 0
        cy = center.y if center else 0

        hw = width / 2.0
        hh = height / 2.0

        mp = BRepBuilderAPI_MakePolygon()
        mp.Add(gp_Pnt(cx - hw, cy - hh, 0))
        mp.Add(gp_Pnt(cx + hw, cy - hh, 0))
        mp.Add(gp_Pnt(cx + hw, cy + hh, 0))
        mp.Add(gp_Pnt(cx - hw, cy + hh, 0))
        mp.Close()
        return OCCGeometryObject(mp.Wire(), {"type": "rectangle"})

    def extrude(
        self,
        sketch: GeometryObject,
        height: float,
        direction: Tuple[float, float, float] = (0, 0, 1),
        symmetric: bool = False,
    ) -> GeometryObject:
        """Extrude a sketch face into a 3D solid along the sketch plane normal.

        The sketch's plane normal is read from ``sketch.metadata["normal"]``
        (set by ``OCCSketch.build_face_geometry``); it defaults to +Z for
        legacy objects that don't carry one.  *direction* is accepted for API
        compatibility but ignored — the plane normal is authoritative.  A
        negative *height* extrudes against the normal.
        """
        from OCP.gp import gp_Vec
        from OCP.BRepPrimAPI import BRepPrimAPI_MakePrism
        from OCP.BRepAlgoAPI import BRepAlgoAPI_Fuse
        from OCP.TopoDS import TopoDS_Shape

        # Defensive: figure out the actual shape from whatever the caller
        # hands us, and surface a clear error if we can't get one.
        if isinstance(sketch, OCCGeometryObject):
            face = self._get_shape(sketch)
        elif isinstance(sketch, TopoDS_Shape):
            face = sketch
        else:
            face = self._get_shape(sketch)
        if face is None:
            raise ValueError(
                "Cannot extrude: sketch has no geometry. Draw a closed profile before extruding."
            )
        # If the wrapper class itself leaked through somehow, surface a
        # clear error instead of letting BRepPrimAPI_MakePrism raise an
        # opaque TypeError.
        if isinstance(face, OCCGeometryObject):
            raise ValueError(
                "Cannot extrude: sketch geometry is a wrapper, not a shape. "
                "This is a bug — please report it."
            )
        # ``face`` may be a TopoDS_Face (new path) or a compound/wire.
        # If it's not already a face, build one.
        face = self._ensure_face(face)
        if face is None:
            raise ValueError(
                "Cannot extrude: sketch geometry is not a valid face. "
                "Ensure the profile is closed (no open ends)."
            )

        normal = self._sketch_normal(sketch)
        nx, ny, nz = normal

        def _prism(h: float):
            vec = gp_Vec(nx * h, ny * h, nz * h)
            maker = BRepPrimAPI_MakePrism(face, vec, False, True)
            maker.Build()
            return maker.Shape()

        if symmetric:
            half = height / 2.0
            pos = _prism(half)
            neg = _prism(-half)
            fuse = BRepAlgoAPI_Fuse(pos, neg)
            fuse.Build()
            solid = fuse.Shape()
        else:
            solid = _prism(height)

        return OCCGeometryObject(solid, {"type": "extrusion", "normal": normal})

    @staticmethod
    def _sketch_normal(obj: GeometryObject) -> Tuple[float, float, float]:
        """Return the normal stored on a sketch-derived geometry object, else +Z."""
        import numpy as np

        meta = getattr(obj, "metadata", None) or {}
        n = meta.get("normal")
        if n is None:
            return (0.0, 0.0, 1.0)
        arr = np.asarray(n, dtype=float)
        norm = float(np.linalg.norm(arr))
        if norm < 1e-12:
            return (0.0, 0.0, 1.0)
        arr = arr / norm
        return (float(arr[0]), float(arr[1]), float(arr[2]))

    @staticmethod
    def _ensure_face(shape: Any) -> Any:
        """Return a ``TopoDS_Face`` from *shape*, or *None* if impossible.

        If *shape* is already a face, return it unchanged; otherwise try to
        build a planar face from it (wire/edge/compound).  Returns *None* for
        empty/invalid input so callers can surface a clear error instead of
        feeding a non-face to ``BRepPrimAPI_MakePrism``.
        """
        from OCP.TopoDS import TopoDS_Face
        from OCP.BRepBuilderAPI import BRepBuilderAPI_MakeFace

        if shape is None:
            return None
        if isinstance(shape, TopoDS_Face):
            return shape
        try:
            maker = BRepBuilderAPI_MakeFace(shape, True)
            maker.Build()
            if maker.IsDone():
                return maker.Face()
        except Exception:
            pass
        return None

    @staticmethod
    def find_coplanar_face(
        shape: Any,
        origin: Tuple[float, float, float],
        normal: Tuple[float, float, float],
        ref_center: Optional[Tuple[float, float, float]] = None,
        angle_tol_deg: float = 5.0,
        dist_tol: float = 1e-3,
    ) -> Optional[Tuple[Any, Tuple[float, float, float]]]:
        """Find a planar face on *shape* coplanar with the given plane.

        Iterates the faces of *shape* and returns the first planar face whose
        plane normal is parallel to *normal* (within *angle_tol_deg* degrees)
        and whose plane passes through *origin* (within *dist_tol* distance).
        When several faces match, the one whose surface centre is closest to
        *ref_center* (if provided) is preferred.

        Returns ``(face, center)`` where *center* is the surface centroid as a
        3-tuple, or *None* if no matching face is found.
        """
        import math

        from OCP.TopExp import TopExp_Explorer
        from OCP.TopAbs import TopAbs_FACE
        from OCP.TopoDS import TopoDS
        from OCP.BRepAdaptor import BRepAdaptor_Surface
        from OCP.GeomAbs import GeomAbs_Plane
        from OCP.BRepGProp import BRepGProp
        from OCP.GProp import GProp_GProps

        import numpy as np

        if shape is None:
            return None
        n = np.asarray(normal, dtype=float)
        n = n / (np.linalg.norm(n) + 1e-30)
        ox, oy, oz = origin

        cos_tol = math.cos(math.radians(angle_tol_deg))
        candidates: list = []

        explorer = TopExp_Explorer(shape, TopAbs_FACE)
        while explorer.More():
            face = TopoDS.Face_s(explorer.Current())
            try:
                surf = BRepAdaptor_Surface(face)
                if surf.GetType() != GeomAbs_Plane:
                    explorer.Next()
                    continue
                plane = surf.Plane()
                pn = np.array(
                    [
                        plane.Axis().Direction().X(),
                        plane.Axis().Direction().Y(),
                        plane.Axis().Direction().Z(),
                    ],
                    dtype=float,
                )
                # Check normals parallel (same or opposite direction)
                cos_angle = abs(float(np.dot(n, pn)))
                if cos_angle < cos_tol:
                    explorer.Next()
                    continue
                # Check distance from plane to origin
                pp = plane.Location()
                d = abs(float(np.dot(n, np.array([pp.X() - ox, pp.Y() - oy, pp.Z() - oz]))))
                if d > dist_tol:
                    explorer.Next()
                    continue
                # Surface centroid via GProp (SurfaceProperties for faces)
                props = GProp_GProps()
                BRepGProp.SurfaceProperties_s(face, props)
                c = props.CentreOfMass()
                center = (float(c.X()), float(c.Y()), float(c.Z()))
                candidates.append((face, center))
            except Exception:
                pass
            explorer.Next()

        if not candidates:
            return None
        if ref_center is not None and len(candidates) > 1:
            rc = np.asarray(ref_center, dtype=float)
            best = min(candidates, key=lambda fc: float(np.linalg.norm(np.asarray(fc[1]) - rc)))
            return best
        return candidates[0]

    def revolve(
        self,
        sketch: GeometryObject,
        angle: float = 360.0,
        axis: Tuple[float, float, float] = (0, 0, 1),
        origin: Tuple[float, float, float] = (0, 0, 0),
    ) -> GeometryObject:
        """Revolve a sketch face around an axis."""
        import math

        # Get the OCC shape directly (a TopoDS_Face for new sketch geometry).
        shape = self._get_shape(sketch)
        face = self._ensure_face(shape)

        from OCP.gp import gp_Ax1, gp_Pnt, gp_Dir
        from OCP.BRepPrimAPI import BRepPrimAPI_MakeRevol

        # Revolve the face around the axis
        revolve_axis = gp_Ax1(gp_Pnt(*origin), gp_Dir(*axis))
        angle_rad = math.radians(angle)
        revolver = BRepPrimAPI_MakeRevol(face, revolve_axis, angle_rad)
        revolver.Build()
        solid_shape = revolver.Shape()

        return OCCGeometryObject(solid_shape, {"type": "revolution"})

    def loft(self, profiles: List[GeometryObject], ruled: bool = False) -> GeometryObject:
        """Create a loft between multiple profiles."""
        from OCP.BRepOffsetAPI import BRepOffsetAPI_ThruSections
        from OCP.TopExp import TopExp_Explorer
        from OCP.TopAbs import TopAbs_WIRE
        from OCP.TopoDS import TopoDS

        if len(profiles) < 2:
            raise ValueError("Loft requires at least 2 profiles")

        loft_maker = BRepOffsetAPI_ThruSections(True, ruled)
        for profile in profiles:
            shape = self._get_shape(profile)
            explorer = TopExp_Explorer(shape, TopAbs_WIRE)
            while explorer.More():
                wire = TopoDS.Wire_s(explorer.Current())
                loft_maker.AddWire(wire)
                explorer.Next()

        loft_maker.Build()
        solid = loft_maker.Shape()
        return OCCGeometryObject(solid, {"type": "loft"})

    def sweep(
        self, profile: GeometryObject, path: GeometryObject, is_frenet: bool = False
    ) -> GeometryObject:
        """Sweep a profile along a path."""
        from OCP.BRepOffsetAPI import BRepOffsetAPI_MakePipeShell
        from OCP.TopExp import TopExp_Explorer
        from OCP.TopAbs import TopAbs_WIRE
        from OCP.TopoDS import TopoDS

        profile_shape = self._get_shape(profile)
        path_shape = self._get_shape(path)

        def _first_wire(shape):
            exp = TopExp_Explorer(shape, TopAbs_WIRE)
            if exp.More():
                return TopoDS.Wire_s(exp.Current())
            raise ValueError("No wire found in shape for sweep")

        profile_wire = _first_wire(profile_shape)
        path_wire = _first_wire(path_shape)

        pipe = BRepOffsetAPI_MakePipeShell(path_wire)
        pipe.Add(profile_wire, False, False)
        if is_frenet:
            pipe.SetMode(True)
        pipe.Build()
        solid = pipe.Shape()
        return OCCGeometryObject(solid, {"type": "sweep"})

    def boolean_union(self, *bodies: GeometryObject) -> GeometryObject:
        """Union multiple bodies."""
        if len(bodies) < 2:
            return bodies[0] if bodies else OCCGeometryObject(None)

        result = self._get_shape(bodies[0])
        for body in bodies[1:]:
            shape = self._get_shape(body)
            from OCP.BRepAlgoAPI import BRepAlgoAPI_Fuse

            fuse = BRepAlgoAPI_Fuse(result, shape)
            fuse.Build()
            result = fuse.Shape()

        return OCCGeometryObject(result, {"type": "union"})

    def boolean_difference(self, base: GeometryObject, tool: GeometryObject) -> GeometryObject:
        """Subtract tool from base."""
        base_shape = self._get_shape(base)
        tool_shape = self._get_shape(tool)

        from OCP.BRepAlgoAPI import BRepAlgoAPI_Cut

        cut = BRepAlgoAPI_Cut(base_shape, tool_shape)
        cut.Build()

        return OCCGeometryObject(cut.Shape(), {"type": "difference"})

    def boolean_intersection(self, body1: GeometryObject, body2: GeometryObject) -> GeometryObject:
        """Intersect two bodies."""
        shape1 = self._get_shape(body1)
        shape2 = self._get_shape(body2)

        from OCP.BRepAlgoAPI import BRepAlgoAPI_Common

        common = BRepAlgoAPI_Common(shape1, shape2)
        common.Build()

        return OCCGeometryObject(common.Shape(), {"type": "intersection"})

    def fillet(
        self, body: GeometryObject, radius: float, edges: Optional[List[Any]] = None
    ) -> GeometryObject:
        """Apply fillet to edges."""
        shape = self._get_shape(body)
        from OCP.BRepFilletAPI import BRepFilletAPI_MakeFillet

        fillet = BRepFilletAPI_MakeFillet(shape)

        if edges:
            for edge in edges:
                fillet.Add(radius, edge)
        else:
            from OCP.TopExp import TopExp_Explorer
            from OCP.TopAbs import TopAbs_EDGE
            from OCP.TopoDS import TopoDS

            explorer = TopExp_Explorer(shape, TopAbs_EDGE)
            while explorer.More():
                fillet.Add(radius, TopoDS.Edge_s(explorer.Current()))
                explorer.Next()

        fillet.Build()
        return OCCGeometryObject(fillet.Shape(), {"type": "fillet"})

    def chamfer(
        self, body: GeometryObject, size: float, edges: Optional[List[Any]] = None
    ) -> GeometryObject:
        """Apply chamfer to edges."""
        shape = self._get_shape(body)
        from OCP.BRepFilletAPI import BRepFilletAPI_MakeChamfer

        chamfer = BRepFilletAPI_MakeChamfer(shape)

        if edges:
            for edge in edges:
                chamfer.Add(size, edge)
        else:
            from OCP.TopExp import TopExp_Explorer
            from OCP.TopAbs import TopAbs_EDGE
            from OCP.TopoDS import TopoDS

            explorer = TopExp_Explorer(shape, TopAbs_EDGE)
            while explorer.More():
                chamfer.Add(size, TopoDS.Edge_s(explorer.Current()))
                explorer.Next()

        chamfer.Build()
        return OCCGeometryObject(chamfer.Shape(), {"type": "chamfer"})

    def shell(
        self, body: GeometryObject, thickness: float, faces_to_remove: Optional[List[Any]] = None
    ) -> GeometryObject:
        """Create a shell (hollow body)."""
        shape = self._get_shape(body)
        from OCP.BRepOffsetAPI import BRepOffsetAPI_MakeThickSolid
        from OCP.TopTools import TopTools_ListOfShape

        faces_list = TopTools_ListOfShape()
        if faces_to_remove:
            for face in faces_to_remove:
                faces_list.Append(face)

        shell_maker = BRepOffsetAPI_MakeThickSolid()
        shell_maker.MakeThickSolidByJoin(shape, faces_list, thickness, 0.001)
        shell_maker.Build()
        return OCCGeometryObject(shell_maker.Shape(), {"type": "shell"})

    def offset(self, face: GeometryObject, distance: float) -> GeometryObject:
        """Offset a face or surface."""
        shape = self._get_shape(face)
        from OCP.BRepOffsetAPI import BRepOffsetAPI_MakeOffset

        offset_maker = BRepOffsetAPI_MakeOffset(shape, False)
        offset_maker.Perform(distance)

        return OCCGeometryObject(offset_maker.Shape(), {"type": "offset"})

    def translate(self, body: GeometryObject, vector: Tuple[float, float, float]) -> GeometryObject:
        """Translate a body."""
        shape = self._get_shape(body)
        from OCP.BRepBuilderAPI import BRepBuilderAPI_Transform
        from OCP.gp import gp_Trsf, gp_Vec

        transform = gp_Trsf()
        transform.SetTranslation(gp_Vec(*vector))
        transformer = BRepBuilderAPI_Transform(shape, transform)
        return OCCGeometryObject(transformer.Shape(), {"type": "translated"})

    def rotate(
        self,
        body: GeometryObject,
        axis: Tuple[float, float, float],
        angle: float,
        origin: Tuple[float, float, float] = (0, 0, 0),
    ) -> GeometryObject:
        """Rotate a body around an axis."""
        shape = self._get_shape(body)
        from OCP.BRepBuilderAPI import BRepBuilderAPI_Transform
        from OCP.gp import gp_Trsf, gp_Ax1, gp_Pnt, gp_Dir

        ax1 = gp_Ax1(gp_Pnt(*origin), gp_Dir(*axis))
        transform = gp_Trsf()
        transform.SetRotation(ax1, angle)
        transformer = BRepBuilderAPI_Transform(shape, transform)
        return OCCGeometryObject(transformer.Shape(), {"type": "rotated"})

    def scale(self, body: GeometryObject, factor: float) -> GeometryObject:
        """Scale a body uniformly."""
        shape = self._get_shape(body)
        from OCP.BRepBuilderAPI import BRepBuilderAPI_Transform
        from OCP.gp import gp_Trsf

        transform = gp_Trsf()
        transform.SetScale(factor)
        transformer = BRepBuilderAPI_Transform(shape, transform)

        return OCCGeometryObject(transformer.Shape(), {"type": "scaled"})

    def mirror(
        self,
        body: GeometryObject,
        plane_normal: Tuple[float, float, float],
        plane_origin: Tuple[float, float, float] = (0, 0, 0),
    ) -> GeometryObject:
        """Mirror a body across a plane."""
        shape = self._get_shape(body)
        from OCP.BRepBuilderAPI import BRepBuilderAPI_Transform
        from OCP.gp import gp_Trsf, gp_Ax2, gp_Pnt, gp_Dir

        ax2 = gp_Ax2(gp_Pnt(*plane_origin), gp_Dir(*plane_normal))
        transform = gp_Trsf()
        transform.SetMirror(ax2)
        transformer = BRepBuilderAPI_Transform(shape, transform)

        return OCCGeometryObject(transformer.Shape(), {"type": "mirrored"})

    def pattern(
        self,
        body: GeometryObject,
        pattern_type: str = "linear",
        count: int = 2,
        direction: Tuple[float, float, float] = (1, 0, 0),
        spacing: float = 10.0,
        axis: Tuple[float, float, float] = (0, 0, 1),
        origin: Tuple[float, float, float] = (0.0, 0.0, 0.0),
        angle: float = 360.0,
    ) -> GeometryObject:
        """Repeat *body* in a linear or circular array (pattern).

        Linear: *count* copies spaced *spacing* mm apart along
        *direction* (a negative spacing reverses the direction).

        Circular: *count* copies rotated evenly around *axis* passing
        through *origin*, distributed over a total angular span of
        *angle* degrees (step = angle / count).  ``angle=360`` gives the
        classic evenly-spaced full-circle bolt pattern.

        Returns the union (compound when the copies don't touch) of the
        original solid and all its copies — disjoint copies keep their
        separate volumes inside one result object, touching copies fuse.
        """
        count = max(1, int(count))
        if count <= 1:
            return body

        import math as _math

        instances: list = [body]
        if pattern_type == "circular":
            # Normalize the rotation axis.
            ax = float(axis[0]), float(axis[1]), float(axis[2])
            norm = _math.sqrt(ax[0] * ax[0] + ax[1] * ax[1] + ax[2] * ax[2])
            if norm < 1e-12:
                ax = (0.0, 0.0, 1.0)
            else:
                ax = (ax[0] / norm, ax[1] / norm, ax[2] / norm)
            step = _math.radians(float(angle)) / count
            for i in range(1, count):
                instances.append(self.rotate(body, ax, step * i, origin))
        else:
            d = float(direction[0]), float(direction[1]), float(direction[2])
            norm = _math.sqrt(d[0] * d[0] + d[1] * d[1] + d[2] * d[2])
            if norm < 1e-12:
                d = (1.0, 0.0, 0.0)
            else:
                d = (d[0] / norm, d[1] / norm, d[2] / norm)
            step = float(spacing)
            for i in range(1, count):
                instances.append(
                    self.translate(
                        body,
                        (d[0] * step * i, d[1] * step * i, d[2] * step * i),
                    )
                )

        return self.boolean_union(*instances)

    def export_step(self, body: GeometryObject, filepath: str, schema: str = "AP214") -> bool:
        """Export to STEP format."""
        try:
            shape = self._get_shape(body)
            from OCP.STEPControl import STEPControl_Writer, STEPControl_AsIs
            from OCP.Interface import Interface_Static

            writer = STEPControl_Writer()
            if schema == "AP214":
                Interface_Static.SetCVal_s("write.step.schema", "AP214")
            elif schema == "AP203":
                Interface_Static.SetCVal_s("write.step.schema", "AP203")

            writer.Transfer(shape, STEPControl_AsIs)
            return writer.Write(filepath)
        except Exception as e:
            print(f"STEP export error: {e}")
            return False

    def export_iges(self, body: GeometryObject, filepath: str) -> bool:
        """Export to IGES format."""
        try:
            shape = self._get_shape(body)
            from OCP.IGESControl import IGESControl_Writer
            from OCP.Interface import Interface_Static

            Interface_Static.SetCVal_s("write.iges.schema", "5.3")
            writer = IGESControl_Writer()
            writer.AddShape(shape)
            return writer.Write(filepath)
        except Exception as e:
            print(f"IGES export error: {e}")
            return False

    def export_stl(
        self, body: GeometryObject, filepath: str, tolerance: float = 0.1, ascii_mode: bool = False
    ) -> bool:
        """Export to STL format."""
        try:
            shape = self._get_shape(body)
            from OCP.StlAPI import StlAPI_Writer
            from OCP.BRepMesh import BRepMesh_IncrementalMesh

            mesh = BRepMesh_IncrementalMesh(shape, tolerance)
            mesh.Perform()

            writer = StlAPI_Writer()
            writer.ASCIIMode = ascii_mode
            return writer.Write(shape, filepath)
        except Exception as e:
            print(f"STL export error: {e}")
            return False

    def import_step(self, filepath: str) -> GeometryObject:
        """Import from STEP format."""
        from OCP.STEPControl import STEPControl_Reader
        from OCP.IFSelect import IFSelect_RetDone

        reader = STEPControl_Reader()
        status = reader.ReadFile(filepath)

        if status != IFSelect_RetDone:
            raise ValueError(f"Failed to read STEP file: {filepath}")

        reader.TransferRoots()
        shape = reader.OneShape()
        return OCCGeometryObject(shape, {"type": "imported_step"})

    def import_step_components(self, filepath: str) -> list:
        """Import a STEP file and return each solid as a separate ``(name, shape)`` pair.

        The STEP reader transfers the entire root shape, then we iterate
        over individual ``TopAbs_SOLID`` entities so that each solid gets
        its own ``OCCGeometryObject``.  If the file contains only a single
        solid the list will have one entry.

        Returns a list of ``(name, OCCGeometryObject)`` tuples.
        """
        from OCP.STEPControl import STEPControl_Reader
        from OCP.IFSelect import IFSelect_RetDone
        from OCP.TopExp import TopExp_Explorer
        from OCP.TopAbs import TopAbs_SOLID
        from OCP.TopoDS import TopoDS

        reader = STEPControl_Reader()
        status = reader.ReadFile(filepath)

        if status != IFSelect_RetDone:
            raise ValueError(f"Failed to read STEP file: {filepath}")

        reader.TransferRoots()
        shape = reader.OneShape()

        # Extract individual solids
        solids: list = []
        explorer = TopExp_Explorer(shape, TopAbs_SOLID)
        idx = 0
        while explorer.More():
            solid = TopoDS.Solid_s(explorer.Current())
            idx += 1
            solids.append(
                (
                    f"Part {idx}",
                    OCCGeometryObject(solid, {"type": "imported_step"}),
                )
            )
            explorer.Next()

        # Fallback: no individual solids found — return the whole shape
        if not solids:
            solids = [("Imported", OCCGeometryObject(shape, {"type": "imported_step"}))]

        return solids

    def import_iges(self, filepath: str) -> GeometryObject:
        """Import from IGES format."""
        from OCP.IGESControl import IGESControl_Reader
        from OCP.IFSelect import IFSelect_RetDone

        reader = IGESControl_Reader()
        status = reader.ReadFile(filepath)

        if status != IFSelect_RetDone:
            raise ValueError(f"Failed to read IGES file: {filepath}")

        reader.TransferRoots()
        shape = reader.OneShape()

        return OCCGeometryObject(shape, {"type": "imported_iges"})

    def get_mesh(
        self, body: GeometryObject, tolerance: float = 0.1
    ) -> Tuple[np.ndarray, np.ndarray]:
        """Get triangulated mesh for rendering."""
        shape = self._get_shape(body)

        from OCP.BRepMesh import BRepMesh_IncrementalMesh
        from OCP.TopExp import TopExp_Explorer
        from OCP.TopAbs import TopAbs_FACE
        from OCP.BRep import BRep_Tool
        from OCP.TopLoc import TopLoc_Location

        # Use finer angular deflection (0.15 rad ≈ 24 segments/circle) so
        # curved surfaces like cylinders render smoothly instead of faceted.
        mesh = BRepMesh_IncrementalMesh(shape, tolerance, False, 0.15, True)
        mesh.Perform()

        vertices_list: List[List[float]] = []
        faces_list: List[List[int]] = []
        vertex_offset = 0

        from OCP.TopoDS import TopoDS
        from OCP.TopAbs import TopAbs_Orientation

        explorer = TopExp_Explorer(shape, TopAbs_FACE)
        while explorer.More():
            face = TopoDS.Face_s(explorer.Current())
            location = TopLoc_Location()
            triangulation = BRep_Tool.Triangulation_s(face, location)

            if triangulation is not None:
                n_vertices = triangulation.NbNodes()
                for i in range(1, n_vertices + 1):
                    p = triangulation.Node(i)
                    vertices_list.append([p.X(), p.Y(), p.Z()])

                n_triangles = triangulation.NbTriangles()
                # REVERSED faces store triangle winding in the natural (surface)
                # orientation — we must flip it so the computed normals point
                # outward (away from solid interior).  TopAbs_REVERSED = 1.
                reverse_winding = face.Orientation() == TopAbs_Orientation.TopAbs_REVERSED
                for i in range(1, n_triangles + 1):
                    tri = triangulation.Triangle(i)
                    v0, v1, v2 = (
                        tri.Value(1) - 1 + vertex_offset,
                        tri.Value(2) - 1 + vertex_offset,
                        tri.Value(3) - 1 + vertex_offset,
                    )
                    if reverse_winding:
                        # Swap last two vertices to flip winding direction.
                        v1, v2 = v2, v1
                    faces_list.append([v0, v1, v2])

                vertex_offset += n_vertices

            explorer.Next()

        return np.array(vertices_list, dtype=np.float32), np.array(faces_list, dtype=np.int32)

    def get_edges(self, body: GeometryObject) -> Tuple[np.ndarray, np.ndarray]:
        """Get edge wireframe for rendering."""
        shape = self._get_shape(body)

        from OCP.TopExp import TopExp_Explorer
        from OCP.TopAbs import TopAbs_EDGE
        from OCP.BRepAdaptor import BRepAdaptor_Curve
        from OCP.GeomAbs import GeomAbs_Line

        vertices_list: List[List[float]] = []
        edges_list: List[List[int]] = []
        vertex_offset = 0

        def discretize_edge(edge: Any, num_points: int = 20) -> List[List[float]]:
            curve = BRepAdaptor_Curve(edge)
            curve_type = curve.GetType()

            points = []

            if curve_type == GeomAbs_Line:
                first = curve.FirstParameter()
                last = curve.LastParameter()
                p1 = curve.Value(first)
                p2 = curve.Value(last)
                points = [[p1.X(), p1.Y(), p1.Z()], [p2.X(), p2.Y(), p2.Z()]]
            else:
                first = curve.FirstParameter()
                last = curve.LastParameter()

                for i in range(num_points + 1):
                    t = first + (last - first) * i / num_points
                    p = curve.Value(t)
                    points.append([p.X(), p.Y(), p.Z()])

            return points

        explorer = TopExp_Explorer(shape, TopAbs_EDGE)
        while explorer.More():
            from OCP.TopoDS import TopoDS

            edge = TopoDS.Edge_s(explorer.Current())
            edge_points = discretize_edge(edge)

            for i, pt in enumerate(edge_points):
                vertices_list.append(pt)
                if i < len(edge_points) - 1:
                    edges_list.append([vertex_offset + i, vertex_offset + i + 1])

            vertex_offset += len(edge_points)
            explorer.Next()

        return np.array(vertices_list, dtype=np.float32), np.array(edges_list, dtype=np.int32)

    def get_bounding_box(self, body: GeometryObject) -> Tuple[Point3D, Point3D]:
        """Get the bounding box of a body."""
        shape = self._get_shape(body)

        from OCP.Bnd import Bnd_Box
        from OCP.BRepBndLib import BRepBndLib

        bbox = Bnd_Box()
        BRepBndLib.AddClose_s(shape, bbox)

        xmin, ymin, zmin, xmax, ymax, zmax = bbox.Get()

        return Point3D(xmin, ymin, zmin), Point3D(xmax, ymax, zmax)

    def get_volume(self, body: GeometryObject) -> float:
        """Calculate the volume of a solid body."""
        shape = self._get_shape(body)

        from OCP.GProp import GProp_GProps
        from OCP.BRepGProp import BRepGProp

        props = GProp_GProps()
        BRepGProp.VolumeProperties_s(shape, props)

        return props.Mass()

    def get_surface_area(self, body: GeometryObject) -> float:
        """Calculate the surface area of a body."""
        shape = self._get_shape(body)

        from OCP.GProp import GProp_GProps
        from OCP.BRepGProp import BRepGProp

        props = GProp_GProps()
        BRepGProp.SurfaceProperties_s(shape, props)

        return props.Mass()

    def get_center_of_mass(self, body: GeometryObject) -> Point3D:
        """Calculate the center of mass of a solid body."""
        shape = self._get_shape(body)

        from OCP.GProp import GProp_GProps
        from OCP.BRepGProp import BRepGProp

        props = GProp_GProps()
        BRepGProp.VolumeProperties_s(shape, props)

        cg = props.CentreOfMass()
        return Point3D(cg.X(), cg.Y(), cg.Z())

    def create_thread(
        self,
        body: GeometryObject,
        cylindrical_face: Any,
        nominal_diameter: float,
        pitch: float,
        thread_length: Optional[float] = None,
        internal: bool = False,
    ) -> Optional[GeometryObject]:
        """Cut (or add) an ISO metric thread on the cylindrical face of *body*.

        The geometry is driven by the PICKED face's actual radius and axis
        (``nominal_diameter`` is only metadata used for the feature record).
        External threads cut the ISO groove trapezoid (7P/8 at the surface,
        P/4 at the root, 5H/8 deep) out of the shaft; internal threads fuse
        the ISO ridge trapezoid (3P/4 at the wall, P/8 crest) into the hole.
        """
        import math

        from OCP.BRepAdaptor import BRepAdaptor_Surface
        from OCP.GeomAbs import GeomAbs_Cylinder
        from OCP.TopoDS import TopoDS
        from OCP.gp import gp_Pnt, gp_Pnt2d, gp_Dir2d
        from OCP.BRepBuilderAPI import (
            BRepBuilderAPI_MakeEdge,
            BRepBuilderAPI_MakeWire,
        )

        # ── 1. Cylinder parameters from the picked face ─────────────────
        try:
            surf = BRepAdaptor_Surface(cylindrical_face)
        except Exception:
            try:
                surf = BRepAdaptor_Surface(TopoDS.Face_s(cylindrical_face))
            except Exception as exc:
                logger.warning(f"create_thread: cannot adapt face: {exc}")
                return None

        if surf.GetType() != GeomAbs_Cylinder:
            logger.warning("create_thread: face is not cylindrical")
            return None

        cyl = surf.Cylinder()                       # gp_Cylinder
        radius = cyl.Radius()                       # ACTUAL picked radius
        ax3 = cyl.Position()                        # gp_Ax3 (location, Z, X)
        loc = ax3.Location()
        zdir = ax3.Direction()
        xdir = ax3.XDirection()
        axis_origin = np.array([loc.X(), loc.Y(), loc.Z()])
        axis_dir = np.array([zdir.X(), zdir.Y(), zdir.Z()])
        axis_dir = axis_dir / np.linalg.norm(axis_dir)
        axis_x = np.array([xdir.X(), xdir.Y(), xdir.Z()])
        axis_x = axis_x / np.linalg.norm(axis_x)
        axis_y = np.cross(axis_dir, axis_x)

        u_start = surf.FirstUParameter()            # angular start of face
        v1, v2 = surf.FirstVParameter(), surf.LastVParameter()
        v_lo, v_hi = min(v1, v2), max(v1, v2)
        face_height = v_hi - v_lo

        if not thread_length or thread_length <= 0:
            thread_length = face_height
        thread_length = min(thread_length, face_height)

        num_turns = thread_length / pitch
        if num_turns < 0.05:
            logger.warning("create_thread: thread too short for one turn")
            return None

        # ── 2. ISO metric profile dimensions ────────────────────────────
        # Basic profile (H = P·√3/2, thread engagement depth 5H/8):
        #  • external shaft: groove cut is a trapezoid 7P/8 wide at the
        #    surface narrowing to P/4 at the root.
        #  • internal hole:  ridge fused onto the wall is a trapezoid 3P/4
        #    wide at the wall narrowing to P/8 at the inner crest, leaving
        #    the 7P/8-wide groove open at the bore.
        H = pitch * math.sqrt(3.0) / 2.0
        depth = (5.0 / 8.0) * H
        overcut = max(0.1 * depth, 0.02)            # overhang past the surface
        if internal:
            w_surf = 3.0 * pitch / 4.0
            w_deep = pitch / 8.0
        else:
            w_surf = 7.0 * pitch / 8.0
            w_deep = pitch / 4.0

        # ── 3. Helix spine ON the picked cylinder's surface ─────────────
        # The swept profile sits in the helix's normal plane, tilted by the
        # lead angle; its end caps therefore stick out past the spine ends
        # by roughly half the profile width along the axis.  For a CUT that
        # is harmless (the groove simply runs to the part edge), but a FUSE
        # would leave the protruding cap as floating material outside the
        # part, so inset the internal helix by exactly that amount.
        lead = math.atan2(pitch, 2.0 * math.pi * radius)
        cap_axial = (w_surf / 2.0) * math.cos(lead)  # cap half-extent along axis

        if internal:
            v_start = v_lo + cap_axial
            v_end = min(v_lo + thread_length, v_hi) - cap_axial
        else:
            # extend one pitch past each face end so the groove runs off
            # the part edges cleanly
            v_start = v_lo - pitch
            v_end = min(v_lo + thread_length + pitch, v_hi + pitch)
        thread_span = v_end - v_start
        if thread_span < 0.5 * pitch:
            logger.warning("create_thread: part too short for a thread")
            return None
        turns_ext = thread_span / pitch

        spine_wire = None

        # 3a. TRUE helix: a 2D straight line on the cylinder surface.
        #
        # NOTE 1: gp_Dir2d NORMALIZES its argument, so the 2D line
        # parameter t advances the point by t·|(2π, pitch)| in (u, v)
        # space — scale the trim range so t = n turns covers exactly
        # n revolutions plus n·pitch of axial travel.
        # NOTE 2: the edge from a pcurve has no 3D curve; the pipe sweep
        # needs one, so force it with BRepLib.BuildCurves3d.
        spine_wire = None
        try:
            from OCP.Geom import Geom_CylindricalSurface
            from OCP.Geom2d import Geom2d_Line, Geom2d_TrimmedCurve
            from OCP.BRepLib import BRepLib

            dir_len = math.hypot(2.0 * math.pi, pitch)
            cyl_surf = Geom_CylindricalSurface(cyl)
            line2d = Geom2d_Line(
                gp_Pnt2d(u_start, v_start), gp_Dir2d(2.0 * math.pi, pitch)
            )
            seg = Geom2d_TrimmedCurve(line2d, 0.0, turns_ext * dir_len)
            helix_edge = BRepBuilderAPI_MakeEdge(seg, cyl_surf).Edge()
            BRepLib.BuildCurves3d_s(helix_edge)
            spine_wire = BRepBuilderAPI_MakeWire(helix_edge).Wire()
            logger.info("create_thread: using exact helix spine")
        except Exception as exc:
            logger.info(f"create_thread: exact helix failed ({exc})")

        # 3b. Fallback: smooth BSpline through sampled helix points
        # (only if the exact construction is unavailable).
        if spine_wire is None:
            try:
                from OCP.GeomAPI import GeomAPI_PointsToBSpline
                from OCP.TColgp import TColgp_Array1OfPnt
                from OCP.GeomAbs import GeomAbs_C2

                pts_per_turn = 96
                n_total = max(int(turns_ext * pts_per_turn) + 1, 2)
                arr = TColgp_Array1OfPnt(1, n_total)
                for i in range(1, n_total + 1):
                    u = u_start + ((i - 1) / pts_per_turn) * 2.0 * math.pi
                    v = v_start + ((i - 1) / pts_per_turn) * pitch
                    p = (
                        axis_origin
                        + radius * (math.cos(u) * axis_x + math.sin(u) * axis_y)
                        + v * axis_dir
                    )
                    arr.SetValue(i, gp_Pnt(float(p[0]), float(p[1]), float(p[2])))
                bspline = GeomAPI_PointsToBSpline(arr, 3, 8, GeomAbs_C2, 1e-5)
                bs_edge = BRepBuilderAPI_MakeEdge(bspline.Curve()).Edge()
                spine_wire = BRepBuilderAPI_MakeWire(bs_edge).Wire()
                logger.info("create_thread: using BSpline helix fallback")
            except Exception as exc:
                logger.warning(f"create_thread: BSpline helix failed ({exc})")

        if spine_wire is None:
            logger.warning("create_thread: no usable helix spine")
            return None

        # Start frame (same for both spine types — computed analytically).
        def _cyl_pt(u: float, v: float) -> np.ndarray:
            return (
                axis_origin
                + radius * (math.cos(u) * axis_x + math.sin(u) * axis_y)
                + v * axis_dir
            )

        start_S = _cyl_pt(u_start, v_start)
        start_T = (
            2.0 * math.pi * radius
            * (-math.sin(u_start) * axis_x + math.cos(u_start) * axis_y)
            + pitch * axis_dir
        )
        start_T = start_T / np.linalg.norm(start_T)
        start_R = math.cos(u_start) * axis_x + math.sin(u_start) * axis_y  # outward

        # Profile width direction: perpendicular to tangent in the surface
        # plane (≈ axial direction).  Trapezoid is symmetric so sign is fine.
        binormal = np.cross(start_T, start_R)
        binormal = binormal / np.linalg.norm(binormal)

        # ── 4. Trapezoidal profile at the spine start ───────────────────
        # Built directly in world coords: base sits *overcut* OUTSIDE the
        # surface so the boolean fuses/cuts cleanly across it; the working
        # end reaches *depth* INSIDE the surface.
        def _mk(b: float, r: float) -> gp_Pnt:
            p = start_S + b * binormal + r * start_R
            return gp_Pnt(float(p[0]), float(p[1]), float(p[2]))

        p0 = _mk(-w_surf / 2.0, overcut)
        p1 = _mk(-w_deep / 2.0, -depth)
        p2 = _mk(+w_deep / 2.0, -depth)
        p3 = _mk(+w_surf / 2.0, overcut)

        prof_wb = BRepBuilderAPI_MakeWire()
        for a, b in ((p0, p1), (p1, p2), (p2, p3), (p3, p0)):
            prof_wb.Add(BRepBuilderAPI_MakeEdge(a, b).Edge())
        profile_wire = prof_wb.Wire()

        # ── 5. Sweep the profile along the helix ────────────────────────
        from OCP.BRepOffsetAPI import BRepOffsetAPI_MakePipeShell

        try:
            pipe = BRepOffsetAPI_MakePipeShell(spine_wire)
            pipe.SetMode(True)                       # Frenet frame
            pipe.Add(profile_wire, False, False)
            pipe.Build()
            if not pipe.IsDone():
                logger.warning("create_thread: pipe sweep failed")
                return None
            solid_ok = False
            try:
                solid_ok = bool(pipe.MakeSolid())    # cap the tube ends
            except Exception as exc:
                logger.info(f"create_thread: MakeSolid unavailable ({exc})")
            tool_shape = pipe.Shape()
            if not solid_ok:
                logger.warning("create_thread: sweep is not a solid")
        except Exception as exc:
            logger.warning(f"create_thread: sweep failed: {exc}")
            return None

        # ── 6. Boolean cut (shaft) or fuse (hole) ───────────────────────
        body_shape = self._get_shape(body)
        if body_shape is None:
            logger.warning("create_thread: body has no shape")
            return None

        tool = OCCGeometryObject(tool_shape)
        vol_before = self.get_volume(body)

        if internal:
            result = self.boolean_union(body, tool)
        else:
            result = self.boolean_difference(body, tool)

        if result is None or self._get_shape(result) is None:
            logger.warning("create_thread: boolean op produced no shape")
            return None

        try:
            vol_after = self.get_volume(result)
        except Exception:
            vol_after = -1.0

        if internal and vol_after <= vol_before:
            logger.warning(
                f"create_thread: fuse did not add volume "
                f"({vol_before:.4f} → {vol_after:.4f}) — tool missed the body?"
            )
            return None
        if not internal and vol_after >= vol_before:
            logger.warning(
                f"create_thread: cut did not remove volume "
                f"({vol_before:.4f} → {vol_after:.4f}) — tool missed the body?"
            )
            return None

        logger.info(
            f"create_thread: {'internal' if internal else 'external'} thread OK, "
            f"volume {vol_before:.4f} → {vol_after:.4f}"
        )
        return result

    def detect_cylindrical_face(
        self,
        face: Any,
    ) -> Optional[Dict[str, Any]]:
        """Check if *face* is cylindrical and return its parameters.

        The *face* can be a ``TopoDS_Face`` (from the picker) or a
        ``TopoDS_Shape`` that contains a face.  We try several paths to
        extract the underlying cylindrical surface.

        Returns a dict with keys ``radius``, ``axis_origin``, ``axis_dir``,
        ``height``, or *None* if the face isn't cylindrical.
        """
        import logging
        import numpy as np
        from OCP.BRepAdaptor import BRepAdaptor_Surface
        from OCP.GeomAbs import GeomAbs_Cylinder
        from OCP.TopoDS import TopoDS

        _log = logging.getLogger(__name__)

        # ── Resolve the actual face from whatever the caller handed us ──
        actual_face: Any = None

        # Try direct BRepAdaptor_Surface first — the picker already returns
        # a valid TopoDS_Face, and calling TopoDS.Face_s() again on an
        # already-downcast face can fail in some OCP versions.
        try:
            surf = BRepAdaptor_Surface(face)
            surf_type_test = surf.GetType()
            actual_face = face
        except Exception:
            pass

        if actual_face is None:
            # Fallback: try the explicit TopoDS.Face_s downcast path.
            try:
                candidate = TopoDS.Face_s(face)
                _ = BRepAdaptor_Surface(candidate)
                actual_face = candidate
            except Exception:
                pass

        if actual_face is None:
            _log.warning("detect_cylindrical_face: could not resolve face from pick result")
            return None

        # ── Probe the surface type ──
        try:
            surf = BRepAdaptor_Surface(actual_face)
            surf_type = surf.GetType()
            if surf_type != GeomAbs_Cylinder:
                type_names = {
                    0: "Plane", 1: "Cylinder", 2: "Cone", 3: "Sphere",
                    4: "Torus", 5: "Bezier", 6: "BSpline", 7: "Revolution",
                    8: "Extrusion", 9: "Offset", 10: "Other",
                }
                type_name = type_names.get(int(surf_type), f"Unknown({int(surf_type)})")
                _log.warning(
                    f"detect_cylindrical_face: face is {type_name}, not a Cylinder"
                )
                return None

            cyl = surf.Cylinder()
            radius = cyl.Radius()
            axis = cyl.Axis()
            origin = axis.Location()
            direction = axis.Direction()
            # BRepAdaptor_Surface uses FirstUParameter/LastUParameter etc.
            u1 = surf.FirstUParameter()
            u2 = surf.LastUParameter()
            v1 = surf.FirstVParameter()
            v2 = surf.LastVParameter()
            height = abs(v2 - v1)
            return {
                "radius": radius,
                "diameter": 2.0 * radius,
                "axis_origin": (origin.X(), origin.Y(), origin.Z()),
                "axis_dir": (direction.X(), direction.Y(), direction.Z()),
                "height": height,
            }
        except Exception as exc:
            _log.warning(f"detect_cylindrical_face: surface probe failed: {exc}")
            return None