step_ops.py

"""
step_ops.py

STEP file parsing, triangulation, and face unfolding CLI module for Pathstitch.
"""

import sys
import json
import argparse
import os
import math
from typing import Dict, List, Any, Tuple, Optional
import ezdxf

from OCC.Core.STEPControl import STEPControl_Reader
from OCC.Core.IFSelect import IFSelect_RetDone
from OCC.Core.TopExp import TopExp_Explorer
from OCC.Core.TopAbs import TopAbs_SOLID, TopAbs_SHELL, TopAbs_FACE, TopAbs_EDGE
from OCC.Core.BRepMesh import BRepMesh_IncrementalMesh
from OCC.Core.BRep import BRep_Tool
from OCC.Core.TopLoc import TopLoc_Location
from OCC.Core.GProp import GProp_GProps
from OCC.Core.BRepGProp import brepgprop

from pathstitch_core.surface_unfold import get_surface_type, unfold_face_geometry, save_polylines_to_dxf, triangulate_face, parameterize_mesh

def load_step_shape(file_path: str):
    """Loads a STEP, STL, or OBJ file and returns its consolidated shape."""
    if not os.path.exists(file_path):
        raise FileNotFoundError(f"3D file not found: {file_path}")
        
    ext = os.path.splitext(file_path)[1].lower()
    
    if ext == ".stl":
        from OCC.Core.StlAPI import StlAPI_Reader
        from OCC.Core.TopoDS import TopoDS_Shape
        reader = StlAPI_Reader()
        shape = TopoDS_Shape()
        success = reader.Read(shape, file_path)
        if not success or shape.IsNull():
            raise ValueError(f"Failed to read STL file: {file_path}")
        return shape
        
    elif ext == ".obj":
        from OCC.Core.gp import gp_Pnt
        from OCC.Core.BRepBuilderAPI import BRepBuilderAPI_MakePolygon, BRepBuilderAPI_MakeFace
        from OCC.Core.TopoDS import TopoDS_Compound
        from OCC.Core.BRep import BRep_Builder
        
        vertices = []
        builder = BRep_Builder()
        compound = TopoDS_Compound()
        builder.MakeCompound(compound)
        
        has_faces = False
        with open(file_path, "r", encoding="utf-8", errors="ignore") as f:
            for line in f:
                line = line.strip()
                if not line or line.startswith("#"):
                    continue
                parts = line.split()
                if not parts:
                    continue
                cmd = parts[0].lower()
                if cmd == "v":
                    if len(parts) >= 4:
                        try:
                            x = float(parts[1])
                            y = float(parts[2])
                            z = float(parts[3])
                            vertices.append(gp_Pnt(x, y, z))
                        except ValueError:
                            pass
                elif cmd == "f":
                    face_vertices = []
                    for p in parts[1:]:
                        v_idx_str = p.split("/")[0]
                        try:
                            v_idx = int(v_idx_str)
                            if v_idx > 0:
                                v_idx = v_idx - 1
                            elif v_idx < 0:
                                v_idx = len(vertices) + v_idx
                            if 0 <= v_idx < len(vertices):
                                face_vertices.append(vertices[v_idx])
                        except ValueError:
                            pass
                    if len(face_vertices) >= 3:
                        try:
                            poly = BRepBuilderAPI_MakePolygon()
                            for v in face_vertices:
                                poly.Add(v)
                            poly.Close()
                            face_maker = BRepBuilderAPI_MakeFace(poly.Wire())
                            if not face_maker.IsDone():
                                continue
                            face = face_maker.Face()
                            if not face.IsNull():
                                builder.Add(compound, face)
                                has_faces = True
                        except Exception:
                            pass
        if not has_faces:
            raise ValueError(f"No valid faces could be parsed from OBJ file: {file_path}")
        return compound
        
    else:
        # Default to STEP
        reader = STEPControl_Reader()
        status = reader.ReadFile(file_path)
        if status != IFSelect_RetDone:
            raise ValueError(f"STEP control reader failed to read. Status code: {status}")
            
        reader.TransferRoots()
        return reader.OneShape()

def get_solid_bodies(shape) -> List[Any]:
    """Isolates and returns all solid bodies (or shells as fallback)."""
    bodies = []
    
    # 1. Search for Solids
    exp = TopExp_Explorer(shape, TopAbs_SOLID)
    while exp.More():
        bodies.append(exp.Current())
        exp.Next()
        
    # 2. Search for Shells if no Solids found
    if not bodies:
        exp = TopExp_Explorer(shape, TopAbs_SHELL)
        while exp.More():
            bodies.append(exp.Current())
            exp.Next()
            
    # 3. Fallback: treat the entire shape as a single body if it contains any faces
    if not bodies:
        exp = TopExp_Explorer(shape, TopAbs_FACE)
        if exp.More():
            bodies.append(shape)
            
    return bodies

def op_list_bodies(args: Dict[str, Any]) -> Dict[str, Any]:
    """
    Triangulates the STEP file bodies and returns their faces, types,
    and 3D coordinates for rendering in Three.js.
    """
    input_path = args.get("input")
    if not input_path:
        return {"status": "error", "message": "Input path must be specified."}
        
    try:
        shape = load_step_shape(input_path)
        bodies = get_solid_bodies(shape)
        
        bodies_data = []
        global_min = [float("inf"), float("inf"), float("inf")]
        global_max = [float("-inf"), float("-inf"), float("-inf")]
        
        for b_idx, body in enumerate(bodies):
            # Run mesh triangulation (0.3mm linear deflection deflection tolerance is a good speed/detail trade-off)
            BRepMesh_IncrementalMesh(body, 0.05)
            
            # Map all edges of this body
            from OCC.Core.TopTools import TopTools_IndexedMapOfShape
            from OCC.Core.TopoDS import topods
            emap = TopTools_IndexedMapOfShape()
            edge_exp = TopExp_Explorer(body, TopAbs_EDGE)
            while edge_exp.More():
                emap.Add(topods.Edge(edge_exp.Current()))
                edge_exp.Next()
                
            # Build edge to faces adjacency map
            edge_to_faces = {}
            faces_list = []
            face_exp = TopExp_Explorer(body, TopAbs_FACE)
            f_idx = 0
            
            while face_exp.More():
                face = face_exp.Current()
                face_exp.Next()
                
                # Traverse edges of this face to record adjacency
                e_exp = TopExp_Explorer(face, TopAbs_EDGE)
                while e_exp.More():
                    edge = topods.Edge(e_exp.Current())
                    eid = emap.Add(edge)
                    edge_to_faces.setdefault(eid, []).append(f_idx)
                    e_exp.Next()
                
                stype = get_surface_type(face)
                
                # Area calculation
                gprops = GProp_GProps()
                brepgprop.SurfaceProperties(face, gprops)
                area = gprops.Mass()
                
                # Retrieve triangulation
                loc = TopLoc_Location()
                tri = BRep_Tool.Triangulation(face, loc)
                
                vertices = []
                indices = []
                
                if tri:
                    trans = loc.Transformation()
                    for i in range(1, tri.NbNodes() + 1):
                        pnt = tri.Node(i).Transformed(trans)
                        px, py, pz = float(pnt.X()), float(pnt.Y()), float(pnt.Z())
                        vertices.extend([px, py, pz])
                        
                        # Update bounding box
                        global_min[0] = min(global_min[0], px)
                        global_min[1] = min(global_min[1], py)
                        global_min[2] = min(global_min[2], pz)
                        
                        global_max[0] = max(global_max[0], px)
                        global_max[1] = max(global_max[1], py)
                        global_max[2] = max(global_max[2], pz)
                        
                    for i in range(1, tri.NbTriangles() + 1):
                        t = tri.Triangle(i)
                        idx1, idx2, idx3 = t.Get()
                        indices.extend([idx1 - 1, idx2 - 1, idx3 - 1])
                
                faces_list.append({
                    "face_index": f_idx,
                    "type": stype,
                    "area": float(area),
                    "vertices": vertices,
                    "indices": indices
                })
                f_idx += 1
                
            # Now build the edges list for this body
            edges_list = []
            from OCC.Core.BRepAdaptor import BRepAdaptor_Curve
            for e_idx in range(1, emap.Size() + 1):
                edge = topods.Edge(emap.FindKey(e_idx))
                try:
                    adaptor = BRepAdaptor_Curve(edge)
                    u_start = adaptor.FirstParameter()
                    u_end = adaptor.LastParameter()
                    try:
                        from OCC.Core.GeomAbs import GeomAbs_Line
                        is_line = (adaptor.GetType() == GeomAbs_Line)
                    except Exception:
                        is_line = False
                    
                    samples = 1 if is_line else 24
                    pts = []
                    for i in range(samples + 1):
                        t = u_start + (u_end - u_start) * (i / samples)
                        p = adaptor.Value(t)
                        pts.extend([float(p.X()), float(p.Y()), float(p.Z())])
                except Exception:
                    pts = []
                
                edges_list.append({
                    "edge_index": e_idx,
                    "vertices": pts,
                    "faces": edge_to_faces.get(e_idx, [])
                })
                
            bodies_data.append({
                "body_index": b_idx,
                "name": f"Body {b_idx + 1}",
                "faces": faces_list,
                "edges": edges_list
            })
            
        # If no vertices were found, reset bounding box
        if global_min[0] == float("inf"):
            global_min = [0.0, 0.0, 0.0]
            global_max = [0.0, 0.0, 0.0]
            
        bbox = {
            "min": global_min,
            "max": global_max,
            "center": [
                (global_min[0] + global_max[0]) / 2.0,
                (global_min[1] + global_max[1]) / 2.0,
                (global_min[2] + global_max[2]) / 2.0
            ]
        }
        
        return {
            "status": "ok",
            "data": {
                "bodies": bodies_data,
                "bbox": bbox
            }
        }
        
    except Exception as e:
        import traceback
        return {"status": "error", "message": f"Failed to list solid bodies: {str(e)}\n{traceback.format_exc()}"}

def get_dxf_bounds(msp) -> Optional[Tuple[float, float, float, float]]:
    """Calculates the 2D bounding box of all renderable geometries in the modelspace."""
    from ezdxf.path import make_path
    min_x, min_y = float('inf'), float('inf')
    max_x, max_y = float('-inf'), float('-inf')
    found = False
    
    for ent in msp:
        if ent.dxftype() == "LINE":
            min_x = min(min_x, ent.dxf.start.x, ent.dxf.end.x)
            max_x = max(max_x, ent.dxf.start.x, ent.dxf.end.x)
            min_y = min(min_y, ent.dxf.start.y, ent.dxf.end.y)
            max_y = max(max_y, ent.dxf.start.y, ent.dxf.end.y)
            found = True
        elif ent.dxftype() in ("CIRCLE", "ARC"):
            cx, cy = ent.dxf.center.x, ent.dxf.center.y
            r = ent.dxf.radius
            min_x = min(min_x, cx - r)
            max_x = max(max_x, cx + r)
            min_y = min(min_y, cy - r)
            max_y = max(max_y, cy + r)
            found = True
        elif ent.dxftype() in ("LWPOLYLINE", "POLYLINE"):
            for p in ent.get_points() if hasattr(ent, 'get_points') else ent.points:
                min_x = min(min_x, p[0])
                max_x = max(max_x, p[0])
                min_y = min(min_y, p[1])
                max_y = max(max_y, p[1])
            found = True
        elif ent.dxftype() in ("SPLINE", "ELLIPSE"):
            try:
                path = make_path(ent)
                for p in path.flattening(distance=0.1):
                    min_x = min(min_x, p.x)
                    max_x = max(max_x, p.x)
                    min_y = min(min_y, p.y)
                    max_y = max(max_y, p.y)
                found = True
            except Exception:
                pass
                
    if not found:
        return None
    return min_x, min_y, max_x, max_y

def op_unfold_face(args: Dict[str, Any]) -> Dict[str, Any]:
    """
    Unfolds a specific face of a specific body to a DXF flat pattern.
    Appends and arranges to the right of any existing DXF geometry if existing_dxf is provided.
    """
    input_path = args.get("input")
    output_path = args.get("output")
    body_idx = args.get("body_index")
    face_idx = args.get("face_index")
    existing_dxf = args.get("existing_dxf")
    
    if not input_path or not output_path or body_idx is None or face_idx is None:
        return {"status": "error", "message": "Missing required arguments: input, output, body_index, face_index."}
        
    try:
        shape = load_step_shape(input_path)
        bodies = get_solid_bodies(shape)
        
        if body_idx < 0 or body_idx >= len(bodies):
            return {"status": "error", "message": f"Body index {body_idx} out of range. Total bodies: {len(bodies)}."}
            
        body = bodies[body_idx]
        
        # Traverse faces to find the target face index
        face_exp = TopExp_Explorer(body, TopAbs_FACE)
        f_idx = 0
        target_face = None
        
        while face_exp.More():
            face = face_exp.Current()
            face_exp.Next()
            if f_idx == face_idx:
                target_face = face
                break
            f_idx += 1
            
        if target_face is None:
            return {"status": "error", "message": f"Face index {face_idx} out of range for body {body_idx}."}
            
        # Call unfolding engine with distortion mode
        distortion_mode = args.get("distortion_mode", "conformal")
        polylines = unfold_face_geometry(target_face, mode=distortion_mode)
        if not polylines:
            return {"status": "error", "message": "No geometry returned from unfolding."}
            
        # Calculate 2D bounds of unfolded shape
        min_x = min(pt[0] for poly in polylines for pt in poly)
        max_x = max(pt[0] for poly in polylines for pt in poly)
        min_y = min(pt[1] for poly in polylines for pt in poly)
        
        # Load or create DXF
        if existing_dxf and os.path.exists(existing_dxf):
            doc = ezdxf.readfile(existing_dxf)
            msp = doc.modelspace()
            bounds = get_dxf_bounds(msp)
            if bounds:
                start_x = bounds[2] + 10.0  # 10mm gap
                start_y = bounds[1]         # Align bottom Y
            else:
                start_x = 0.0
                start_y = 0.0
        else:
            doc = ezdxf.new(dxfversion="R2010")
            msp = doc.modelspace()
            start_x = 0.0
            start_y = 0.0
            
        if "UNFOLDED_3D" not in doc.layers:
            doc.layers.new("UNFOLDED_3D", dxfattribs={"color": 6})

        # Translate to correct position and add to layout
        for poly in polylines:
            translated = []
            for pt in poly:
                tx = pt[0] - min_x + start_x
                ty = pt[1] - min_y + start_y
                translated.append((tx, ty))
            if len(translated) >= 2:
                msp.add_lwpolyline(translated, dxfattribs={"layer": "UNFOLDED_3D"})
                
        doc.saveas(output_path)
        
        return {
            "status": "ok",
            "data": {
                "body_index": body_idx,
                "face_index": face_idx,
                "output": output_path,
                "polylines_count": len(polylines)
            }
        }
        
    except Exception as e:
        return {"status": "error", "message": f"Failed to unfold face: {str(e)}"}

def op_unfold_faces(args: Dict[str, Any]) -> Dict[str, Any]:
    """
    Unfolds multiple faces side-by-side and saves to a combined DXF.
    Appends and arranges to the right of any existing DXF geometry if existing_dxf is provided.
    """
    input_path = args.get("input")
    output_path = args.get("output")
    faces_to_unfold = args.get("faces")
    existing_dxf = args.get("existing_dxf")
    
    if not input_path or not output_path or not faces_to_unfold:
        return {"status": "error", "message": "Missing required arguments: input, output, faces."}
        
    try:
        shape = load_step_shape(input_path)
        bodies = get_solid_bodies(shape)
        
        # Load or create DXF
        if existing_dxf and os.path.exists(existing_dxf):
            doc = ezdxf.readfile(existing_dxf)
            msp = doc.modelspace()
            bounds = get_dxf_bounds(msp)
            if bounds:
                current_x_offset = bounds[2] + 10.0  # 10mm gap
                current_y_offset = bounds[1]
            else:
                current_x_offset = 0.0
                current_y_offset = 0.0
        else:
            doc = ezdxf.new(dxfversion="R2010")
            msp = doc.modelspace()
            current_x_offset = 0.0
            current_y_offset = 0.0
            
        if "UNFOLDED_3D" not in doc.layers:
            doc.layers.new("UNFOLDED_3D", dxfattribs={"color": 6})
            
        gap = 10.0 # 10mm gap between unfolded layouts
        unfolded_count = 0
        
        for item in faces_to_unfold:
            body_idx = item.get("body_index")
            face_idx = item.get("face_index")
            if body_idx is None or face_idx is None:
                continue
                
            if body_idx < 0 or body_idx >= len(bodies):
                continue
            body = bodies[body_idx]
            
            # Find face
            face_exp = TopExp_Explorer(body, TopAbs_FACE)
            f_idx = 0
            target_face = None
            while face_exp.More():
                face = face_exp.Current()
                face_exp.Next()
                if f_idx == face_idx:
                    target_face = face
                    break
                f_idx += 1
                
            if target_face is None:
                continue
                
            # Unfold face geometry with distortion mode
            distortion_mode = args.get("distortion_mode", "conformal")
            polylines = unfold_face_geometry(target_face, mode=distortion_mode)
            if not polylines:
                continue
                
            # Compute 2D bounding box
            min_x = min(pt[0] for poly in polylines for pt in poly)
            max_x = max(pt[0] for poly in polylines for pt in poly)
            min_y = min(pt[1] for poly in polylines for pt in poly)
            
            # Translate to current horizontal offset and align base Y
            translated_polylines = []
            for poly in polylines:
                translated_poly = []
                for pt in poly:
                    tx = pt[0] - min_x + current_x_offset
                    ty = pt[1] - min_y + current_y_offset
                    translated_poly.append((tx, ty))
                translated_polylines.append(translated_poly)
                
            for pts in translated_polylines:
                if len(pts) >= 2:
                    msp.add_lwpolyline(pts, dxfattribs={"layer": "UNFOLDED_3D"})
                    
            width = max_x - min_x
            current_x_offset += width + gap
            unfolded_count += 1
            
        doc.saveas(output_path)
        return {
            "status": "ok",
            "data": {
                "output": output_path,
                "unfolded_count": unfolded_count
            }
        }
    except Exception as e:
        return {"status": "error", "message": f"Failed to unfold multiple faces: {str(e)}"}

def discretize_edge(edge, num_points: int = 30) -> List[Tuple[float, float, float]]:
    from OCC.Core.BRepAdaptor import BRepAdaptor_Curve
    from OCC.Core.BRep import BRep_Tool
    adaptor = BRepAdaptor_Curve(edge)
    first = adaptor.FirstParameter()
    last = adaptor.LastParameter()
    
    if math.isinf(first) or math.isinf(last):
        from OCC.Core.TopExp import TopExp
        v1 = TopExp.FirstVertex(edge)
        v2 = TopExp.LastVertex(edge)
        p1 = BRep_Tool.Pnt(v1)
        p2 = BRep_Tool.Pnt(v2)
        return [(p1.X(), p1.Y(), p1.Z()), (p2.X(), p2.Y(), p2.Z())]
        
    pts = []
    for i in range(num_points):
        t = first + (last - first) * i / (num_points - 1)
        p = adaptor.Value(t)
        pts.append((p.X(), p.Y(), p.Z()))
    return pts

def get_face_plane_basis(face):
    from OCC.Core.BRepAdaptor import BRepAdaptor_Surface
    from OCC.Core.gp import gp_Pnt, gp_Vec
    
    surf = BRepAdaptor_Surface(face)
    u_mid = (surf.FirstUParameter() + surf.LastUParameter()) / 2.0
    v_mid = (surf.FirstVParameter() + surf.LastVParameter()) / 2.0
    
    pnt = gp_Pnt()
    u_vec = gp_Vec()
    v_vec = gp_Vec()
    surf.D1(u_mid, v_mid, pnt, u_vec, v_vec)
    
    normal_vec = u_vec.Crossed(v_vec)
    if normal_vec.Magnitude() < 1e-6:
        normal_vec = gp_Vec(0, 0, 1)
    else:
        normal_vec.Normalize()
        
    u_vec.Normalize()
    v_axis = normal_vec.Crossed(u_vec)
    v_axis.Normalize()
    
    origin = (pnt.X(), pnt.Y(), pnt.Z())
    normal = (normal_vec.X(), normal_vec.Y(), normal_vec.Z())
    u_axis = (u_vec.X(), u_vec.Y(), u_vec.Z())
    v_axis_coords = (v_axis.X(), v_axis.Y(), v_axis.Z())
    
    return origin, normal, u_axis, v_axis_coords

def _translate_shape(shape, dx: float, dy: float, dz: float):
    """Returns a rigidly translated copy of `shape`. Identity (returns the input)
    when the offset is zero, to avoid needless OCC transform churn."""
    if dx == 0.0 and dy == 0.0 and dz == 0.0:
        return shape
    from OCC.Core.gp import gp_Trsf, gp_Vec
    from OCC.Core.BRepBuilderAPI import BRepBuilderAPI_Transform
    trsf = gp_Trsf()
    trsf.SetTranslation(gp_Vec(dx, dy, dz))
    return BRepBuilderAPI_Transform(shape, trsf, True).Shape()


def op_project_edges(args: Dict[str, Any]) -> Dict[str, Any]:
    input_path = args.get("input")
    output_path = args.get("output")
    body_idx = args.get("body_index", 0)
    plane_type = args.get("plane_type", "XY") # "XY", "XZ", "YZ", or "face"
    face_idx = args.get("face_index")
    face_body_idx = args.get("face_body_index", body_idx)
    existing_dxf = args.get("existing_dxf")
    offset = float(args.get("offset", 0.0))
    visible_bodies_indices = args.get("visible_bodies")
    # Per-body manual move offsets from the 3D gizmo (MAS-140), keyed by body
    # index: {"0": [x, y, z], ...}. Applied so the projection matches exactly how
    # the model is arranged in the 3D view — moved bodies project where they sit.
    body_offsets = args.get("body_offsets") or {}

    def _offset_for(idx: int):
        o = body_offsets.get(str(idx))
        if o is None:
            o = body_offsets.get(idx)
        if not o:
            return (0.0, 0.0, 0.0)
        return (float(o[0]), float(o[1]), float(o[2]))

    if not input_path or not output_path:
        return {"status": "error", "message": "Missing input or output path."}

    try:
        shape = load_step_shape(input_path)
        bodies = get_solid_bodies(shape)
        # Apply each body's manual move so all downstream geometry (the face plane
        # basis, the section, and the silhouette fallback) reflects the moved pose.
        bodies = [_translate_shape(b, *_offset_for(i)) for i, b in enumerate(bodies)]

        target_bodies = []
        if visible_bodies_indices is not None:
            for idx in visible_bodies_indices:
                if 0 <= idx < len(bodies):
                    target_bodies.append(bodies[idx])
        else:
            if 0 <= body_idx < len(bodies):
                target_bodies.append(bodies[body_idx])
            else:
                target_bodies = bodies
                
        if not target_bodies:
            return {"status": "error", "message": "No solid bodies found to project."}
            
        origin = (0.0, 0.0, 0.0)
        normal = (0.0, 0.0, 1.0)
        u_axis = (1.0, 0.0, 0.0)
        v_axis = (0.0, 1.0, 0.0)
        
        if plane_type == "XY":
            pass
        elif plane_type == "XZ":
            normal = (0.0, 1.0, 0.0)
            v_axis = (0.0, 0.0, 1.0)
        elif plane_type == "YZ":
            normal = (1.0, 0.0, 0.0)
            u_axis = (0.0, 1.0, 0.0)
            v_axis = (0.0, 0.0, 1.0)
        elif plane_type == "face" and face_idx is not None:
            if 0 <= face_body_idx < len(bodies):
                body = bodies[face_body_idx]
                face_exp = TopExp_Explorer(body, TopAbs_FACE)
                f_idx = 0
                target_face = None
                while face_exp.More():
                    face = face_exp.Current()
                    face_exp.Next()
                    if f_idx == face_idx:
                        target_face = face
                        break
                    f_idx += 1
                if target_face is None:
                    return {"status": "error", "message": f"Face index {face_idx} not found on body {face_body_idx}."}
                basis = get_face_plane_basis(target_face)
                if basis:
                    origin, normal, u_axis, v_axis = basis
            else:
                return {"status": "error", "message": f"Body index {face_body_idx} out of range for face selection."}
                
        # Shift origin along normal by offset
        origin = (
            origin[0] + normal[0] * offset,
            origin[1] + normal[1] * offset,
            origin[2] + normal[2] * offset
        )
        
        from OCC.Core.gp import gp_Ax2, gp_Pnt, gp_Dir, gp_Pln
        from OCC.Core.BRepAlgoAPI import BRepAlgoAPI_Section
        
        proj_origin = gp_Pnt(origin[0], origin[1], origin[2])
        proj_normal = gp_Dir(normal[0], normal[1], normal[2])
        proj_u = gp_Dir(u_axis[0], u_axis[1], u_axis[2])
        
        # Build projection plane Pln for intersection section
        proj_pln = gp_Pln(proj_origin, proj_normal)
        
        section_edges = []
        for body in target_bodies:
            try:
                sec = BRepAlgoAPI_Section(body, proj_pln)
                sec.Build()
                if sec.IsDone():
                    sec_shape = sec.Shape()
                    edge_exp = TopExp_Explorer(sec_shape, TopAbs_EDGE)
                    while edge_exp.More():
                        edge = edge_exp.Current()
                        section_edges.append(edge)
                        edge_exp.Next()
            except Exception:
                pass
                
        polylines = []
        seen_projections = set()
        
        if section_edges:
            # Intersection Mode: project 3D section curves onto local u, v axes
            for edge in section_edges:
                try:
                    pts3d = discretize_edge(edge)
                    pts2d = []
                    for pt in pts3d:
                        dx = pt[0] - origin[0]
                        dy = pt[1] - origin[1]
                        dz = pt[2] - origin[2]
                        u = dx * u_axis[0] + dy * u_axis[1] + dz * u_axis[2]
                        v = dx * v_axis[0] + dy * v_axis[1] + dz * v_axis[2]
                        pts2d.append((u, v))
                        
                    if len(pts2d) < 2:
                        continue
                    xs = [p[0] for p in pts2d]
                    ys = [p[1] for p in pts2d]
                    if (max(xs) - min(xs)) < 1e-6 and (max(ys) - min(ys)) < 1e-6:
                        continue
                    key_fwd = tuple((round(p[0], 4), round(p[1], 4)) for p in pts2d)
                    key = min(key_fwd, key_fwd[::-1])
                    if key in seen_projections:
                        continue
                    seen_projections.add(key)
                    polylines.append(pts2d)
                except Exception:
                    pass
                    
        # MAS-126: the plane imports ONLY the geometry it actually intersects. The
        # full-silhouette projection of all visible bodies is a fallback used
        # exclusively when the plane intersects nothing at all — never when a
        # section was found (even if those section curves degenerated away).
        any_intersection = len(section_edges) > 0
        if not polylines and not any_intersection:
            # Silhouette/HLR Mode: Project visible outlines/boundaries
            from OCC.Core.HLRAlgo import HLRAlgo_Projector
            from OCC.Core.HLRBRep import HLRBRep_Algo, HLRBRep_HLRToShape
            
            proj_axes = gp_Ax2(proj_origin, proj_normal, proj_u)
            projector = HLRAlgo_Projector(proj_axes)
            
            hlr = HLRBRep_Algo()
            for body in target_bodies:
                hlr.Add(body)
            hlr.Projector(projector)
            hlr.Update()
            hlr.Hide()
            
            hlr_to_shape = HLRBRep_HLRToShape(hlr)
            
            compounds = []
            v_comp = hlr_to_shape.VCompound()
            out_comp = hlr_to_shape.OutLineVCompound()
            
            if v_comp is not None:
                compounds.append(v_comp)
            if out_comp is not None:
                compounds.append(out_comp)
                
            for comp in compounds:
                edge_exp = TopExp_Explorer(comp, TopAbs_EDGE)
                while edge_exp.More():
                    edge = edge_exp.Current()
                    edge_exp.Next()
                    
                    try:
                        pts3d = discretize_edge(edge)
                        pts2d = []
                        for pt in pts3d:
                            pts2d.append((pt[0], pt[1]))
                        if len(pts2d) < 2:
                            continue
                        xs = [p[0] for p in pts2d]
                        ys = [p[1] for p in pts2d]
                        if (max(xs) - min(xs)) < 1e-6 and (max(ys) - min(ys)) < 1e-6:
                            continue
                        key_fwd = tuple((round(p[0], 4), round(p[1], 4)) for p in pts2d)
                        key = min(key_fwd, key_fwd[::-1])
                        if key in seen_projections:
                            continue
                        seen_projections.add(key)
                        polylines.append(pts2d)
                    except Exception:
                        pass
                        
        if not polylines:
            return {"status": "error", "message": "No projectable edges found."}
            
        min_x = min(pt[0] for poly in polylines for pt in poly)
        min_y = min(pt[1] for poly in polylines for pt in poly)
        
        if existing_dxf and os.path.exists(existing_dxf):
            doc = ezdxf.readfile(existing_dxf)
            msp = doc.modelspace()
            bounds = get_dxf_bounds(msp)
            if bounds:
                start_x = bounds[2] + 10.0
                start_y = bounds[1]
            else:
                start_x = 0.0
                start_y = 0.0
        else:
            doc = ezdxf.new(dxfversion="R2010")
            msp = doc.modelspace()
            start_x = 0.0
            start_y = 0.0
            
        if "PROJECTED_SKETCH" not in doc.layers:
            doc.layers.new("PROJECTED_SKETCH", dxfattribs={"color": 5})
            
        for poly in polylines:
            translated = []
            for pt in poly:
                tx = pt[0] - min_x + start_x
                ty = pt[1] - min_y + start_y
                translated.append((tx, ty))
            msp.add_lwpolyline(translated, dxfattribs={"layer": "PROJECTED_SKETCH"})
            
        doc.saveas(output_path)
        return {
            "status": "ok",
            "data": {
                "output": output_path,
                "polylines_count": len(polylines)
            }
        }
    except Exception as e:
        import traceback
        return {"status": "error", "message": f"Projection failed: {str(e)}\n{traceback.format_exc()}"}

def op_combine_steps(args: Dict[str, Any]) -> Dict[str, Any]:
    """Merges two STEP files into one (MAS-125): every solid/shell from `input`
    and `incoming` is packed into a single compound and written to `output`, so a
    model dragged into an already-loaded 3D workspace is appended rather than
    replacing it. The viewport's loader handles the side-by-side distribution."""
    input_path = args.get("input")
    incoming_path = args.get("incoming")
    output_path = args.get("output")
    if not input_path or not incoming_path or not output_path:
        return {"status": "error", "message": "Missing input, incoming, or output path."}
    if not os.path.exists(input_path):
        return {"status": "error", "message": f"Existing STEP not found: {input_path}"}
    if not os.path.exists(incoming_path):
        return {"status": "error", "message": f"Incoming STEP not found: {incoming_path}"}

    try:
        from OCC.Core.TopoDS import TopoDS_Compound
        from OCC.Core.BRep import BRep_Builder
        from OCC.Core.STEPControl import STEPControl_Writer, STEPControl_AsIs
        from OCC.Core.Bnd import Bnd_Box
        from OCC.Core.BRepBndLib import brepbndlib

        def _bbox(shape):
            box = Bnd_Box()
            brepbndlib.Add(shape, box)
            if box.IsVoid():
                return None
            return box.Get()  # (xmin, ymin, zmin, xmax, ymax, zmax)

        existing_bodies = get_solid_bodies(load_step_shape(input_path))
        incoming_bodies = get_solid_bodies(load_step_shape(incoming_path))

        # Place the incoming model beside the existing geometry, every time, so
        # repeated imports lay out as a tidy, non-overlapping row instead of
        # stacking at the same origin (the "spacing gets disturbed" report).
        # The shift is a single rigid translation of the whole incoming group, so
        # a multi-body file keeps its internal arrangement; the EXISTING bodies are
        # never moved, so already-placed models stay exactly where they are.
        GAP = 20.0
        ex_boxes = [b for b in (_bbox(s) for s in existing_bodies) if b]
        in_boxes = [b for b in (_bbox(s) for s in incoming_bodies) if b]
        if ex_boxes and in_boxes:
            ex_xmax = max(b[3] for b in ex_boxes)
            in_xmin = min(b[0] for b in in_boxes)
            dx = (ex_xmax + GAP) - in_xmin
            incoming_bodies = [_translate_shape(b, dx, 0.0, 0.0) for b in incoming_bodies]

        builder = BRep_Builder()
        compound = TopoDS_Compound()
        builder.MakeCompound(compound)

        total = 0
        for body in existing_bodies:
            builder.Add(compound, body)
            total += 1
        for body in incoming_bodies:
            builder.Add(compound, body)
            total += 1

        writer = STEPControl_Writer()
        writer.Transfer(compound, STEPControl_AsIs)
        writer.Write(output_path)
        return {"status": "ok", "data": {"output": output_path, "body_count": total}}
    except Exception as e:
        import traceback
        return {"status": "error", "message": f"Failed to combine STEP files: {str(e)}\n{traceback.format_exc()}"}


def op_face_distortion(args: Dict[str, Any]) -> Dict[str, Any]:
    """Computes the 2D parameterization of a face under a given distortion mode,
    and returns the per-vertex distortion values (area ratios) to color the 3D mesh."""
    import numpy as np
    
    input_path = args.get("input")
    body_idx = args.get("body_index")
    face_idx = args.get("face_index")
    mode = args.get("distortion_mode", "conformal")
    
    if not input_path or body_idx is None or face_idx is None:
        return {"status": "error", "message": "Missing required arguments: input, body_index, face_index."}
        
    try:
        shape = load_step_shape(input_path)
        bodies = get_solid_bodies(shape)
        if body_idx < 0 or body_idx >= len(bodies):
            return {"status": "error", "message": f"Body index {body_idx} out of range."}
        body = bodies[body_idx]
        
        # Traverse faces to find the target face index
        face_exp = TopExp_Explorer(body, TopAbs_FACE)
        f_idx = 0
        target_face = None
        while face_exp.More():
            face = face_exp.Current()
            face_exp.Next()
            if f_idx == face_idx:
                target_face = face
                break
            f_idx += 1
            
        if target_face is None:
            return {"status": "error", "message": f"Face index {face_idx} out of range."}
            
        BRepMesh_IncrementalMesh(target_face, 0.05)
        verts3d, tris = triangulate_face(target_face)
        if not tris:
            return {"status": "ok", "distortion": []}
            
        n = len(verts3d)
        
        stype = get_surface_type(target_face)
        if stype in ("Plane", "Cylinder", "Cone") and mode == "conformal":
            return {"status": "ok", "distortion": [0.0] * n}
            
        uv = parameterize_mesh(verts3d, tris, mode)
        
        from collections import defaultdict
        vert_distortion = defaultdict(list)
        
        V = np.array(verts3d)
        for t_idx, (a, b, c) in enumerate(tris):
            p0, p1, p2 = V[a], V[b], V[c]
            a3d = 0.5 * np.linalg.norm(np.cross(p1 - p0, p2 - p0))
            a3d = max(a3d, 1e-12)
            
            u0, u1, u2 = uv[a], uv[b], uv[c]
            cross = (u1[0] - u0[0]) * (u2[1] - u0[1]) - (u2[0] - u0[0]) * (u1[1] - u0[1])
            a2d = 0.5 * abs(cross)
            
            # Symmetric area distortion: max(a2d/a3d, a3d/a2d) - 1.0
            dist = max(a2d / a3d, a3d / a2d) - 1.0
            
            vert_distortion[a].append(dist)
            vert_distortion[b].append(dist)
            vert_distortion[c].append(dist)
            
        distortion = []
        for i in range(n):
            d_list = vert_distortion.get(i, [0.0])
            distortion.append(float(np.mean(d_list)))
            
        return {"status": "ok", "distortion": distortion}
    except Exception as e:
        import traceback
        return {"status": "error", "message": f"Failed to compute distortion: {str(e)}\n{traceback.format_exc()}"}


# Imported at the bottom: net_unfold imports step_ops helpers lazily inside its
# op, so this ordering avoids a circular import.
from pathstitch_core.net_unfold import op_unfold_connected

# Module-level dispatch table, shared by the CLI (`main`) and the persistent
# worker (`pathstitch_core.worker`).
OPERATIONS = {
    "list_bodies": op_list_bodies,
    "unfold_face": op_unfold_face,
    "unfold_faces": op_unfold_faces,
    "unfold_connected": op_unfold_connected,
    "project_edges": op_project_edges,
    "combine_steps": op_combine_steps,
    "face_distortion": op_face_distortion,
}


def main() -> None:
    """CLI entry point for JSON subprocess interactions."""
    parser = argparse.ArgumentParser(description="Pathstitch STEP operations CLI tool.")
    parser.add_argument("--json", type=str, help="JSON execution configuration.")
    args = parser.parse_args()

    config_str = ""
    if args.json:
        config_str = args.json
    else:
        config_str = sys.stdin.read()

    try:
        config = json.loads(config_str)
    except Exception as e:
        print(json.dumps({"status": "error", "message": f"Failed to parse input JSON: {str(e)}"}))
        sys.exit(1)

    op = config.get("op")
    op_args = config.get("args", {})

    if op not in OPERATIONS:
        print(json.dumps({"status": "error", "message": f"Unknown operation: {op}"}))
        sys.exit(1)

    try:
        result = OPERATIONS[op](op_args)
        print(json.dumps(result))
    except Exception as e:
        print(json.dumps({"status": "error", "message": f"Operation failed: {str(e)}"}))
        sys.exit(1)

if __name__ == "__main__":
    main()