add_sleeves.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Add Sleeves to Holes - Visual Hole Selection and Sleeve Installation
====================================================================
Combined tool for detecting holes and adding sleeves to reduce hole diameter.
Follows the proven workflow: detect tessellated -> replace with analytical -> add sleeves.

Usage:
    python add_sleeves.py model.step
    python add_sleeves.py model.step --min-diameter 3.0 --max-diameter 5.0
    python add_sleeves.py model.step --sleeve-target-diameter 3.6
"""

import sys
import math
import argparse
import numpy as np
from pathlib import Path
from typing import List, Dict, Tuple, Optional
from collections import defaultdict

try:
    from build123d import *
    from OCP.BRep import BRep_Builder
    from OCP.TopoDS import TopoDS_Compound
    from OCP.BRepAdaptor import BRepAdaptor_Surface, BRepAdaptor_Curve
    from OCP.GeomAbs import GeomAbs_Cylinder
    try:
        from ocp_vscode import show
        OCP_VSCODE_AVAILABLE = True
    except ImportError:
        OCP_VSCODE_AVAILABLE = False
        print("Note: ocp_vscode not available - visual display disabled")
except ImportError as e:
    print(f"Error: Required packages not found: {e}")
    print("Please install: pip install build123d ocp_vscode")
    sys.exit(1)


# === STEP 1: TESSELLATED HOLE DETECTION (from detect_holes.py) ===

def get_edge_endpoints(edge_wrapped) -> Tuple[np.ndarray, np.ndarray]:
    """Extract endpoints from edge using OCP."""
    try:
        curve_adaptor = BRepAdaptor_Curve(edge_wrapped)
        first_param = curve_adaptor.FirstParameter()
        last_param = curve_adaptor.LastParameter()
        
        p1 = curve_adaptor.Value(first_param)
        p2 = curve_adaptor.Value(last_param)
        
        return np.array([p1.X(), p1.Y(), p1.Z()]), np.array([p2.X(), p2.Y(), p2.Z()])
    except Exception:
        return None, None


def fit_circle_to_points(points: List[np.ndarray]) -> Tuple[Optional[np.ndarray], float, np.ndarray, float]:
    """Fit a circle to 3D points using least-squares."""
    try:
        points = np.array(points)
        centroid = points.mean(axis=0)
        centered = points - centroid
        
        # Find best-fit plane using SVD
        _, _, Vt = np.linalg.svd(centered)
        normal = Vt[-1]
        
        # Create local 2D coordinate system on the plane
        if abs(normal[2]) < 0.9:
            u = np.cross(normal, [0, 0, 1])
        else:
            u = np.cross(normal, [1, 0, 0])
        u = u / np.linalg.norm(u)
        v = np.cross(normal, u)
        
        # Project to 2D
        points_2d = np.column_stack([np.dot(centered, u), np.dot(centered, v)])
        
        # Fit circle in 2D using least squares
        A = np.column_stack([2*points_2d[:, 0], 2*points_2d[:, 1], np.ones(len(points_2d))])
        b = points_2d[:, 0]**2 + points_2d[:, 1]**2
        
        result, _, _, _ = np.linalg.lstsq(A, b, rcond=None)
        cx_2d, cy_2d, c = result
        r = math.sqrt(max(0, c + cx_2d**2 + cy_2d**2))
        
        # Convert back to 3D
        center_3d = centroid + cx_2d * u + cy_2d * v
        
        # Calculate fitting error
        distances = [np.linalg.norm(p - center_3d) for p in points]
        error = np.std(distances) / r if r > 0 else float('inf')
        
        return center_3d, r, normal, error
        
    except Exception:
        return None, 0.0, np.array([0, 0, 1]), float('inf')


def detect_tessellated_holes(model: Solid, min_diameter: float = 4.0, max_diameter: float = 4.5, max_error: float = 0.02, min_segments: int = 6) -> List[Dict]:
    """Detect cylindrical holes in tessellated geometry."""
    # Extract all edges
    edge_data = [get_edge_endpoints(e.wrapped) for e in model.edges()]
    edge_data = [(p1, p2) for p1, p2 in edge_data if p1 is not None and p2 is not None]
    
    # Build vertex-to-edge adjacency
    def point_key(p, precision=0.001):
        return tuple(np.round(p / precision) * precision)
    
    vertex_to_edges = defaultdict(list)
    for i, (p1, p2) in enumerate(edge_data):
        vertex_to_edges[point_key(p1)].append((i, p1, p2))
        vertex_to_edges[point_key(p2)].append((i, p2, p1))
    
    def find_loop_from_edge(start_idx, used_globally, max_size=100):
        """Find a closed edge loop starting from given edge."""
        p1, p2 = edge_data[start_idx]
        start_key = point_key(p1)
        visited = {start_idx}
        loop_points = [p1]
        current_point = p2
        current_key = point_key(p2)
        
        while len(loop_points) < max_size:
            if current_key == start_key:
                # Closed loop found
                edges_in_loop = list(visited)
                return edges_in_loop, loop_points
            
            loop_points.append(current_point)
            
            # Find next edge
            candidates = [
                (edge_idx, next_point) for edge_idx, _, next_point in vertex_to_edges[current_key]
                if edge_idx not in visited and edge_idx not in used_globally
            ]
            
            if not candidates:
                break
            
            next_edge_idx, next_point = candidates[0]
            visited.add(next_edge_idx)
            current_point = next_point
            current_key = point_key(current_point)
        
        return [], []
    
    all_circles = []
    used_edges = set()
    
    for start_idx in range(len(edge_data)):
        if start_idx in used_edges:
            continue
            
        edges_in_loop, points = find_loop_from_edge(start_idx, used_edges)
        
        if edges_in_loop and len(points) >= min_segments:
            center, radius, axis, error = fit_circle_to_points(points)
            diameter = radius * 2
            
            if (center is not None and 
                error < max_error and 
                min_diameter <= diameter <= max_diameter):
                
                all_circles.append({
                    'center': center,
                    'diameter': diameter,
                    'radius': radius,
                    'axis': axis,
                    'error': error,
                    'num_segments': len(points)
                })
                used_edges.update(edges_in_loop)
    
    return all_circles


# === STEP 2: HOLE REPLACEMENT (from replace_holes.py) ===

def create_cylinder_at_position(x: float, y: float, z: float, axis: List[float], radius: float, length: float) -> Solid:
    """Create a cylinder centered at specified position with given axis direction and length.

    The cylinder is centered at (x, y, z), extending length/2 in both directions along the axis.
    """
    try:
        # Normalize axis vector
        axis_vec = np.array(axis)
        axis_vec = axis_vec / np.linalg.norm(axis_vec)

        # Create cylinder along Z axis, centered at origin (extends from -length/2 to +length/2)
        with BuildPart() as cyl:
            with BuildSketch(Plane.XY.offset(-length/2)):
                Circle(radius)
            extrude(amount=length)

        # Calculate rotation to align with axis direction
        z_axis = np.array([0, 0, 1])

        # If axis is already aligned with Z, no rotation needed
        if np.allclose(axis_vec, z_axis) or np.allclose(axis_vec, -z_axis):
            if np.dot(axis_vec, z_axis) < 0:
                cylinder = cyl.part.rotate(Axis.X, 180)
            else:
                cylinder = cyl.part
        else:
            # Calculate rotation axis (cross product)
            rotation_axis = np.cross(z_axis, axis_vec)
            rotation_axis = rotation_axis / np.linalg.norm(rotation_axis)

            # Calculate rotation angle
            cos_angle = np.dot(z_axis, axis_vec)
            angle = math.degrees(math.acos(np.clip(cos_angle, -1.0, 1.0)))

            # Apply rotation
            rotation_axis_build123d = Vector(*rotation_axis)
            cylinder = cyl.part.rotate(Axis(origin=Vector(0, 0, 0), direction=rotation_axis_build123d), angle)

        # Translate to final position (cylinder center is now at x, y, z)
        cylinder = cylinder.translate(Vector(x, y, z))

        return cylinder

    except Exception as e:
        print(f"Error creating cylinder: {e}")
        return None


def create_compound_from_cylinders(cylinders: List[Solid]) -> Solid:
    """Create a compound solid from multiple cylinders."""
    try:
        if not cylinders:
            return None
        
        if len(cylinders) == 1:
            return cylinders[0]
        
        # Use build123d's compound creation
        compound = cylinders[0]
        for cylinder in cylinders[1:]:
            compound = compound + cylinder
        
        return compound
        
    except Exception as e:
        print(f"Error creating compound: {e}")
        return None


def replace_tessellated_with_analytical(model: Solid, tessellated_holes: List[Dict], cylinder_length: float = 20.0) -> Solid:
    """Replace tessellated holes with analytical cylinders."""
    if not tessellated_holes:
        return model
    
    print(f"Replacing {len(tessellated_holes)} tessellated holes with analytical cylinders...")
    
    # Create cylinders for tessellated holes
    cylinders = []
    for i, hole in enumerate(tessellated_holes, 1):
        print(f"  Cylinder {i}/{len(tessellated_holes)} @ ({hole['center'][0]:.2f}, {hole['center'][1]:.2f}, {hole['center'][2]:.2f})")
        
        cylinder = create_cylinder_at_position(
            hole['center'][0], hole['center'][1], hole['center'][2],
            hole['axis'], hole['radius'], cylinder_length
        )
        
        if cylinder is not None:
            cylinders.append(cylinder)
        else:
            print(f"  WARNING: Failed to create cylinder {i}")
    
    if not cylinders:
        print("Error: No cylinders created")
        return model
    
    print(f"Successfully created {len(cylinders)} cylinders")
    
    # Perform boolean cut operation
    print("Performing boolean cut operation...")
    try:
        original_volume = model.volume
        
        # Create compound of all cylinders
        combined_tool = create_compound_from_cylinders(cylinders)
        
        if combined_tool is None:
            print("Error: Failed to create compound tool")
            return model
        
        # Single boolean operation
        result = model - combined_tool
        
        result_volume = result.volume
        volume_diff = original_volume - result_volume
        
        print(f"  Original volume: {original_volume:.2f} mm³")
        print(f"  Result volume: {result_volume:.2f} mm³")
        print(f"  Volume removed: {volume_diff:.2f} mm³")
        
        return result
        
    except Exception as e:
        print(f"Error during boolean operation: {e}")
        return model


# === STEP 3: ANALYTICAL CYLINDER DETECTION (from holeshrink.py) ===

def get_cylinder_info(face):
    """Extract cylinder axis, center point, and radius from a cylindrical face."""
    try:
        # Check if face is cylindrical
        if face.geom_type != GeomType.CYLINDER:
            return None

        # Get radius using face.radius property
        radius = face.radius
        if radius is None:
            return None

        # Get the underlying surface for axis information using OCP
        surf_adaptor = BRepAdaptor_Surface(face.wrapped, True)

        if surf_adaptor.GetType() == GeomAbs_Cylinder:
            cyl = surf_adaptor.Cylinder()
            axis = cyl.Axis()

            # Get axis direction
            direction = axis.Direction()
            axis_vector = Vector(direction.X(), direction.Y(), direction.Z())

            # Get axis location (point on axis)
            location = axis.Location()
            axis_point = Vector(location.X(), location.Y(), location.Z())

            # Get face center
            center = face.center()

            return {
                'axis_vector': axis_vector,
                'axis_point': axis_point,
                'radius': radius,
                'face_center': center,
                'face': face
            }
    except Exception as e:
        pass

    return None


def find_analytical_cylinders(model: Solid, target_diameter: float, tolerance: float = 0.15) -> List[Dict]:
    """Find all cylindrical faces matching the target diameter."""
    target_radius = target_diameter / 2
    hole_faces = []

    print(f"   Searching for analytical cylinders with radius {target_radius:.2f}mm +/- {tolerance:.2f}mm")

    for face in model.faces():
        cyl_info = get_cylinder_info(face)
        if cyl_info:
            radius = cyl_info['radius']
            if abs(radius - target_radius) < tolerance:
                hole_faces.append(cyl_info)

    print(f"   Found {len(hole_faces)} cylindrical faces")
    return hole_faces


def group_coaxial_faces(faces_info: List[Dict], axis_tolerance: float = 0.5, proximity_tolerance: float = 2.0, debug: bool = False) -> List[List[Dict]]:
    """Group cylindrical faces that share the same axis and are spatially close.

    Uses face centers to determine if faces belong to the same hole.
    Two faces belong to the same hole if:
    1. Their axes are parallel
    2. Their face centers lie on the same axis line (within tolerance)
    3. They are spatially close to each other
    """
    holes = []
    used = set()

    if debug:
        print(f"   DEBUG: Grouping {len(faces_info)} faces with axis_tol={axis_tolerance}, prox_tol={proximity_tolerance}")

    for i, face1 in enumerate(faces_info):
        if i in used:
            continue

        # Start a new hole group
        hole_group = [face1]
        used.add(i)

        axis1 = face1['axis_vector']
        center1 = face1['face_center']

        if debug:
            print(f"   DEBUG: Starting group from face {i}: center=({center1.X:.2f}, {center1.Y:.2f}, {center1.Z:.2f})")

        # Find all faces coaxial with this one AND spatially close
        for j, face2 in enumerate(faces_info):
            if j in used or j <= i:
                continue

            axis2 = face2['axis_vector']
            center2 = face2['face_center']

            # Check if axes are parallel (or anti-parallel)
            axis_parallel = abs(abs(axis1.dot(axis2)) - 1.0) < 0.01

            if axis_parallel:
                # Check if face centers lie on the same axis line
                # Vector from center1 to center2
                center_diff = center2 - center1
                center_distance = center_diff.length

                if center_distance < 0.01:
                    # Centers are essentially the same point - same hole
                    hole_group.append(face2)
                    used.add(j)
                    if debug:
                        print(f"   DEBUG:   Added face {j} (same center)")
                else:
                    # Check if center_diff is parallel to axis (centers on same line)
                    # Cross product of axis and center_diff should be zero if parallel
                    cross = axis1.cross(center_diff)
                    perpendicular_distance = cross.length

                    if debug and perpendicular_distance < axis_tolerance:
                        print(f"   DEBUG:   Face {j}: perp_dist={perpendicular_distance:.3f}, center_dist={center_distance:.2f}")

                    # The perpendicular distance from center2 to the axis through center1
                    # If this is small, centers are on the same axis line
                    if perpendicular_distance < axis_tolerance:
                        # Centers are on the same axis line, now check proximity
                        if center_distance < proximity_tolerance:
                            hole_group.append(face2)
                            used.add(j)
                            if debug:
                                print(f"   DEBUG:   Added face {j} (on axis, close)")
                        elif debug:
                            print(f"   DEBUG:   Rejected face {j} (on axis but too far: {center_distance:.2f} > {proximity_tolerance})")

        if debug:
            print(f"   DEBUG: Group has {len(hole_group)} faces")
        holes.append(hole_group)

    return holes


def analyze_analytical_hole(hole_faces_info: List[Dict]) -> Dict:
    """Analyze a group of coaxial faces to determine hole parameters."""
    # Use first face for axis information
    first_face = hole_faces_info[0]
    axis_vector = first_face['axis_vector'].normalized()
    axis_point = first_face['axis_point']
    radius = first_face['radius']

    # Find the extent of the hole along its axis
    # Project all face centers onto the axis
    projections = []

    for face_info in hole_faces_info:
        center = face_info['face_center']
        # Vector from axis point to face center
        to_center = center - axis_point
        # Project onto axis
        projection_length = to_center.dot(axis_vector)
        projections.append(projection_length)

    # Get min and max projections
    min_proj = min(projections)
    max_proj = max(projections)

    # Hole length
    length = abs(max_proj - min_proj)

    # If length is very small, this might be a single face - use face bounding box
    if length < 0.5:
        # Get the bounding box of all faces
        all_points = []
        for face_info in hole_faces_info:
            bbox = face_info['face'].bounding_box()
            all_points.extend([
                Vector(bbox.min.X, bbox.min.Y, bbox.min.Z),
                Vector(bbox.max.X, bbox.max.Y, bbox.max.Z)
            ])

        # Project bounding box points onto axis
        projections = []
        for point in all_points:
            to_point = point - axis_point
            projection_length = to_point.dot(axis_vector)
            projections.append(projection_length)

        min_proj = min(projections)
        max_proj = max(projections)
        length = abs(max_proj - min_proj)

    # Add some margin to ensure sleeve covers the hole
    length = length + 0.2

    # Center point of the hole
    mid_proj = (min_proj + max_proj) / 2
    center_point = axis_point + axis_vector * mid_proj

    return {
        'center': center_point,
        'axis': axis_vector,
        'length': length,
        'radius': radius,
        'diameter': radius * 2,
        'num_faces': len(hole_faces_info)
    }


# === STEP 4: SLEEVE CREATION (from holeshrink.py) ===

def create_sleeve_for_hole(hole_info: Dict, new_radius: float) -> Solid:
    """Create a sleeve insert for a hole."""
    center = hole_info['center']
    axis = hole_info['axis']
    length = hole_info['length']
    outer_radius = hole_info['radius']

    # Build the sleeve
    with BuildPart() as sleeve_part:
        # Create outer cylinder
        with BuildSketch(Plane(origin=center, z_dir=axis)) as sk:
            Circle(outer_radius)
        extrude(amount=length/2, both=True)

        # Subtract inner cylinder
        with BuildSketch(Plane(origin=center, z_dir=axis)) as sk:
            Circle(new_radius)
        extrude(amount=length/2 + 0.1, both=True, mode=Mode.SUBTRACT)

    return sleeve_part.part


# === STEP 5: VISUAL DISPLAY ===

def display_holes_with_model(model: Solid, holes: List[Dict]) -> bool:
    """Display model with hole markers using ocp_vscode."""
    if not OCP_VSCODE_AVAILABLE:
        return False
    
    try:
        print(f"\nDisplaying model with {len(holes)} hole markers...")
        print("Hole markers are colored spheres at hole centers.")
        
        # Create hole markers as small spheres
        markers = []
        for i, hole in enumerate(holes):
            center = hole['center']
            if hasattr(center, 'to_tuple'):
                center = center.to_tuple()
            
            # Create small sphere at hole center
            with BuildPart() as marker_part:
                with BuildSketch(Plane.XY):
                    Circle(0.5)  # Small marker sphere
                extrude(amount=1.0)
                
            # Move to hole center
            marker = marker_part.part.translate(Vector(*center))
            markers.append(marker)
            
            # Print hole info with color
            colors = ["red", "green", "blue", "yellow", "cyan", "magenta", "orange", "purple"]
            color = colors[i % len(colors)]
            hole_type = hole.get('type', 'unknown')
            print(f"  #{i+1}: {color} marker at ({center[0]:.2f}, {center[1]:.2f}, {center[2]:.2f}) - Ø{hole['diameter']:.2f} mm ({hole_type})")
        
        print("Using port 3939")
        print("++++")
        
        # Create marker names
        marker_names = [f"Hole {i+1}" for i in range(len(markers))]
        
        # Define distinct colors for each hole (RGBA format)
        marker_colors = []
        for i in range(len(markers)):
            # Use predefined colors cycling
            base_colors = [
                (1.0, 0.3, 0.2, 0.8),  # Red
                (0.3, 1.0, 0.2, 0.8),  # Green
                (0.2, 0.3, 1.0, 0.8),  # Blue
                (1.0, 1.0, 0.2, 0.8),  # Yellow
                (0.2, 1.0, 1.0, 0.8),  # Cyan
                (1.0, 0.2, 1.0, 0.8),  # Magenta
                (1.0, 0.7, 0.2, 0.8),  # Orange
                (0.7, 0.2, 1.0, 0.8),  # Purple
            ]
            marker_colors.append(base_colors[i % len(base_colors)])
        
        # Show with proper names and colors
        show(
            model, *markers,
            names=["Model"] + marker_names,
            colors=[(0.7, 0.7, 0.75, 0.4)] + marker_colors  # Semi-transparent model + colored markers
        )
        
        return True
        
    except Exception as e:
        print(f"Error displaying model: {e}")
        return False


def display_model_with_sleeves(model: Solid, sleeves: List[Solid]) -> bool:
    """Display model with sleeves in different color."""
    if not OCP_VSCODE_AVAILABLE:
        return False
    
    try:
        print(f"\nDisplaying model with {len(sleeves)} sleeves...")
        print("Sleeves are shown in red color.")
        print("Using port 3939")
        print("++++")
        
        # Create sleeve names
        sleeve_names = [f"Sleeve {i+1}" for i in range(len(sleeves))]
        
        # Define colors (semi-transparent model, red sleeves)
        model_color = (0.7, 0.7, 0.75, 0.4)  # Semi-transparent gray
        sleeve_colors = [(0.9, 0.3, 0.2, 0.8)] * len(sleeves)  # Red sleeves
        
        # Show with proper names and colors
        show(
            model, *sleeves,
            names=["Model"] + sleeve_names,
            colors=[model_color] + sleeve_colors
        )
        
        return True
        
    except Exception as e:
        print(f"Visual display not available: {e}")
        return False


# === USER INTERACTION ===

def get_user_hole_selection(holes: List[Dict], auto_select_all: bool = False) -> List[int]:
    """Get user selection of which holes to process."""
    if not holes:
        return []
    
    # Check if we're in interactive mode
    auto_select_all_final = auto_select_all
    try:
        # Windows PowerShell/CMD check - stdin.isatty() is more reliable
        if hasattr(sys.stdin, 'isatty') and sys.stdin.isatty():
            auto_select_all_final = False  # We're in a real terminal
        else:
            auto_select_all_final = True
    except:
        auto_select_all_final = True
    
    if auto_select_all_final:
        print("Auto-selecting all holes (non-interactive mode)")
        return list(range(len(holes)))
    
    print()
    print("  #   Diameter (mm)     X        Y        Z      Type        Segments")
    print("  " + "-" * 70)
    
    for i, hole in enumerate(holes, 1):
        if hasattr(hole['center'], 'to_tuple'):
            c = hole['center'].to_tuple()
        else:
            c = hole['center']
        # Handle both analytical and tessellated holes
        segments = hole.get('num_segments', hole.get('num_faces', 0))
        hole_type = hole.get('type', 'unknown')
        print(f"{i:3d}     {hole['diameter']:7.2f}    {c[0]:7.2f}  {c[1]:7.2f}  {c[2]:7.2f}  {hole_type:10s}   {segments:3d}")
    
    print()
    print("Select which holes to install sleeves for:")
    print("  - Enter hole numbers separated by commas (e.g., 1,2,3)")
    print("  - Enter 'all' to select all holes")
    print("  - Enter 'none' or empty to skip sleeve installation")
    
    while True:
        try:
            selection = input("\nYour selection: ").strip().lower()
            
            if selection == "" or selection == "none":
                return []
            
            if selection == "all":
                return list(range(len(holes)))
            
            # Parse comma-separated numbers
            indices = []
            for part in selection.split(','):
                part = part.strip()
                if part:
                    try:
                        num = int(part)
                        if 1 <= num <= len(holes):
                            indices.append(num - 1)  # Convert to 0-based index
                        else:
                            print(f"Error: '{num}' is not a valid hole number (1-{len(holes)})")
                    except ValueError:
                        print(f"Error: '{part}' is not a valid hole number")
                        continue
            
            if indices:
                return sorted(list(set(indices)))  # Remove duplicates and sort
            else:
                print("No valid hole numbers entered.")
                
        except (KeyboardInterrupt, EOFError):
            print("\nAuto-selecting all holes (input not available)")
            return list(range(len(holes)))
        except Exception as e:
            print(f"Error parsing selection: {e}")
            # Fallback to all holes
            return list(range(len(holes)))


# === MAIN WORKFLOW ===

def main():
    parser = argparse.ArgumentParser(
        description="Add sleeves to holes in STEP files to reduce hole diameter",
        formatter_class=argparse.RawDescriptionHelpFormatter,
        epilog="""
Examples:
  %(prog)s model.step                                    # Default M4 holes (4.0-4.5 mm)
  %(prog)s model.step --min-diameter 3.0 --max-diameter 5.0
  %(prog)s model.step --sleeve-target-diameter 3.6 --sleeve-thickness 0.25
        """
    )
    
    parser.add_argument('step_file', help='Input STEP file')
    
    # Hole detection parameters
    parser.add_argument('--min-diameter', type=float, default=4.0,
                        help='Minimum hole diameter in mm (default: 4.0)')
    parser.add_argument('--max-diameter', type=float, default=4.5,
                        help='Maximum hole diameter in mm (default: 4.5)')
    parser.add_argument('--max-error', type=float, default=0.02,
                        help='Maximum circle fitting error (default: 0.02)')
    
    # Sleeve parameters
    parser.add_argument('--sleeve-target-diameter', type=float, default=3.6,
                        help='Target inner diameter for sleeves in mm (default: 3.6)')
    parser.add_argument('--sleeve-thickness', type=float, default=0.5,
                        help='Minimum sleeve wall thickness in mm (default: 0.5)')
    
    # Replacement parameters
    parser.add_argument('--cylinder-length', type=float, default=20.0,
                        help='Cylinder length for hole replacement in mm (default: 20.0)')
    
    # Output options
    parser.add_argument('--output', '-o', type=str,
                        help='Output STEP file (default: input_with_sleeves.step)')
    parser.add_argument('--no-visual', action='store_true',
                        help='Skip visual display (for non-interactive use)')
    parser.add_argument('--auto-all', action='store_true',
                        help='Automatically select all detected holes')
    
    args = parser.parse_args()
    
    # Validate input file
    step_path = Path(args.step_file)
    if not step_path.exists():
        print(f"Error: STEP file not found: {step_path}")
        sys.exit(1)
    
    # Determine output filename
    if args.output:
        output_path = Path(args.output)
    else:
        output_path = step_path.with_stem(step_path.stem + "_with_sleeves")
    
    print(f"Loading STEP file: {step_path}")
    
    # Load model
    try:
        model = import_step(str(step_path))
        print(f"Model loaded successfully")
        print(f"Model volume: {model.volume:.2f} mm³")
    except Exception as e:
        print(f"Error loading STEP file: {e}")
        sys.exit(1)
    
    print()
    print(f"Detecting holes (Ø {args.min_diameter:.1f}-{args.max_diameter:.1f} mm)...")
    
    # STEP 1: Find analytical cylinders
    print("Looking for analytical cylinders...")
    analytical_holes = []
    try:
        # Calculate average target diameter for analytical search
        avg_target = (args.min_diameter + args.max_diameter) / 2
        analytical_faces = find_analytical_cylinders(model, avg_target, tolerance=0.25)
        
        if analytical_faces:
            print(f"Found {len(analytical_faces)} analytical cylindrical faces")
            
            # Group coaxial faces into holes
            grouped_holes = group_coaxial_faces(analytical_faces)
            print(f"Grouped into {len(grouped_holes)} unique analytical holes")
            
            # Analyze each hole group
            for hole_group in grouped_holes:
                hole_info = analyze_analytical_hole(hole_group)
                hole_info['center'] = hole_info['center'].to_tuple()
                hole_info['axis'] = hole_info['axis'].to_tuple()
                hole_info['type'] = 'analytical'
                analytical_holes.append(hole_info)
        
        # STEP 2: Find tessellated holes (don't replace yet)
        print("Looking for tessellated holes...")
        tessellated_holes = detect_tessellated_holes(
            model,
            min_diameter=args.min_diameter,
            max_diameter=args.max_diameter,
            max_error=args.max_error
        )
        
        # Mark tessellated holes
        for hole in tessellated_holes:
            hole['type'] = 'tessellated'
        
        print(f"Found {len(tessellated_holes)} tessellated holes")
        
        # Combine both types for display and selection
        all_holes = analytical_holes + tessellated_holes
        
    except Exception as e:
        print(f"Error detecting holes: {e}")
        sys.exit(1)
    
    if not all_holes:
        print("No holes detected matching the criteria.")
        return 0
    
    print(f"Total holes found: {len(all_holes)} ({len(analytical_holes)} analytical, {len(tessellated_holes)} tessellated)")
    
    # STEP 3: Display visual interface and get user selection
    if not args.no_visual:
        visual_displayed = display_holes_with_model(model, all_holes)
        if not visual_displayed:
            print("Visual display not available. Using text interface only.")
    else:
        print("Visual display disabled by --no-visual flag.")
    
    # Get user selection for sleeve installation
    selected_indices = get_user_hole_selection(all_holes, auto_select_all=args.auto_all)
    
    if not selected_indices:
        print("No holes selected for sleeve installation. Exiting.")
        return 0
    
    # Filter selected holes
    selected_holes = [all_holes[i] for i in selected_indices]
    
    # STEP 4: Replace tessellated holes with analytical cylinders for selected holes only
    tessellated_selected = [hole for hole in selected_holes if hole.get('type') == 'tessellated']
    analytical_selected = [hole for hole in selected_holes if hole.get('type') == 'analytical']
    
    if tessellated_selected:
        print(f"\nReplacing {len(tessellated_selected)} selected tessellated holes with analytical cylinders...")
        model = replace_tessellated_with_analytical(model, tessellated_selected, args.cylinder_length)
        
        # Re-detect analytical cylinders for the replaced holes
        print("Re-detecting analytical cylinders in updated model...")
        analytical_faces = find_analytical_cylinders(model, avg_target, tolerance=0.25)
        
        if analytical_faces:
            grouped_holes = group_coaxial_faces(analytical_faces)
            
            # Find which analytical holes correspond to our selected tessellated holes
            # Only keep ONE analytical hole for each tessellated hole
            new_analytical_holes = []
            used_tessellated_indices = set()
            
            for hole_group in grouped_holes:
                hole_info = analyze_analytical_hole(hole_group)
                hole_center = hole_info['center'].to_tuple()
                
                # Find the closest tessellated hole that hasn't been used yet
                best_distance = float('inf')
                best_tessellated_idx = None
                
                for i, tessellated_hole in enumerate(tessellated_selected):
                    if i in used_tessellated_indices:
                        continue  # This tessellated hole already matched
                        
                    tess_center = tessellated_hole['center']
                    if hasattr(tess_center, 'to_tuple'):
                        tess_center = tess_center.to_tuple()
                    
                    # Calculate distance between centers
                    distance = ((hole_center[0] - tess_center[0])**2 + 
                              (hole_center[1] - tess_center[1])**2 + 
                              (hole_center[2] - tess_center[2])**2)**0.5
                    
                    # If close and closer than previous candidates
                    if distance < 5.0 and distance < best_distance:
                        best_distance = distance
                        best_tessellated_idx = i
                
                # If we found a close match, add this analytical hole
                if best_tessellated_idx is not None:
                    tess_center = tessellated_selected[best_tessellated_idx]['center']
                    if hasattr(tess_center, 'to_tuple'):
                        tess_center = tess_center.to_tuple()
                    
                    
                    hole_info['center'] = hole_info['center'].to_tuple()
                    hole_info['axis'] = hole_info['axis'].to_tuple()
                    hole_info['type'] = 'analytical'
                    new_analytical_holes.append(hole_info)
                    used_tessellated_indices.add(best_tessellated_idx)
            
            print(f"Found {len(new_analytical_holes)} analytical holes matching replaced tessellated holes")
            all_selected_analytical = analytical_selected + new_analytical_holes
        else:
            print("Warning: Could not find analytical cylinders after replacement!")
            all_selected_analytical = analytical_selected
    else:
        all_selected_analytical = analytical_selected
    
    print(f"\nProcessing {len(all_selected_analytical)} holes for sleeve installation...")
    
    # STEP 5: Create sleeves for selected holes
    sleeve_target_radius = args.sleeve_target_diameter / 2
    
    print(f"Creating sleeves to reduce hole diameter to Ø{args.sleeve_target_diameter:.1f}mm")
    print()
    
    sleeves = []
    for i, hole in enumerate(all_selected_analytical, 1):
        print(f"  Sleeve {i}/{len(all_selected_analytical)}: {hole['num_faces']} faces, "
              f"length={hole['length']:.2f}mm, "
              f"axis=({hole['axis'][0]:.2f}, {hole['axis'][1]:.2f}, {hole['axis'][2]:.2f})")
        
        # Convert back to Vector objects for sleeve creation
        hole_info_vectors = hole.copy()
        hole_info_vectors['center'] = Vector(*hole['center'])
        hole_info_vectors['axis'] = Vector(*hole['axis'])
        
        sleeve = create_sleeve_for_hole(hole_info_vectors, sleeve_target_radius)
        sleeves.append(sleeve)
    
    print(f"Successfully created {len(sleeves)} sleeves")
    print()
    
    # Combine all sleeves with model
    print("Combining model with sleeves...")
    try:
        original_volume = model.volume
        
        if sleeves:
            all_sleeves = sleeves[0]
            for sleeve in sleeves[1:]:
                all_sleeves = all_sleeves + sleeve
            
            # Combine with original model
            result = model + all_sleeves
            
            result_volume = result.volume
            volume_added = result_volume - original_volume
            
            print(f"  Original volume: {original_volume:.2f} mm³")
            print(f"  Result volume: {result_volume:.2f} mm³")
            print(f"  Volume added: {volume_added:.2f} mm³")
        else:
            result = model
            
    except Exception as e:
        print(f"Error combining model with sleeves: {e}")
        return 1
    
    # Display result with sleeves if visual is enabled
    if not args.no_visual and sleeves:
        display_model_with_sleeves(model, sleeves)
    
    # Export result
    try:
        export_step(result, str(output_path))
        print(f"\nSaved: {output_path}")
        print("SUCCESS: Sleeve installation completed successfully")
        print(f"Holes processed: {len(all_selected_analytical)}")
        print(f"Sleeves installed: {len(sleeves)}")
        print(f"New hole diameter: Ø{args.sleeve_target_diameter:.1f}mm")
        return 0
            
    except Exception as e:
        print(f"Error saving result: {e}")
        return 1


if __name__ == "__main__":
    sys.exit(main())