opensubdiv-petite 0.3.1

//! Integration with the monstertruck CAD kernel for B-rep surface generation.
//!
//! This module provides `From` trait implementations to convert OpenSubdiv
//! patches to monstertruck's surface representations, enabling high-order
//! surface export to STEP format.

use crate::bfr::SurfaceFactory as BfrSurfaceFactory;
use crate::far::{PatchEvalResult, PatchTable, PatchType, TopologyRefiner};
use monstertruck::geometry::prelude::{BsplineCurve, BsplineSurface, KnotVector, ParametricCurve};
#[cfg(feature = "monstertruck_export_boundary")]
use monstertruck::modeling::{Curve, Edge, Vertex, Wire};
use monstertruck::modeling::{
    EuclideanSpace, Face, InnerSpace, MetricSpace, Point3, Shell, Surface, Vector3,
};
#[cfg(feature = "rayon")]
use rayon::prelude::*;
use std::{convert::TryFrom, panic};
use thiserror::Error;

/// Type alias for results in this module.
pub type Result<T> = std::result::Result<T, MonstertruckError>;

/// Error type for monstertruck integration.
#[derive(Debug, Clone, Error)]
pub enum MonstertruckError {
    /// Unsupported patch type
    #[error("Unsupported patch type: {0:?}")]
    UnsupportedPatchType(PatchType),

    /// BFR surface conversion failed
    #[error("BFR conversion failed: {0}")]
    BfrConversionFailed(String),

    /// Invalid control point configuration
    #[error("Invalid control points configuration")]
    InvalidControlPoints,

    /// Patch evaluation failed
    #[error("Patch evaluation failed")]
    EvaluationFailed,

    /// Invalid knot vector
    #[error("Invalid knot vector")]
    InvalidKnotVector,
}

/// How to handle extraordinary vertices (valence != 4) during STEP export.
///
/// Extraordinary vertices require special treatment because the Catmull-Clark
/// limit surface at these points cannot be exactly represented as a bicubic
/// B-spline patch.
#[derive(Debug, Clone, Copy, Default, PartialEq, Eq)]
pub enum GregoryAccuracy {
    /// Use BSplineBasis end caps (~0.1% deviation, guaranteed compatibility).
    ///
    /// This is the recommended default. OpenSubdiv generates B-spline patches
    /// at extraordinary vertices with small approximation error. Results in
    /// patches that are guaranteed to be valid B-splines for CAD export.
    #[default]
    BSplineEndCaps,

    /// Evaluate Gregory patches at 8x8 grid and fit B-spline (higher accuracy).
    ///
    /// When higher accuracy is needed at extraordinary vertices, evaluate the
    /// true Gregory patch at a denser grid and fit a B-spline surface to the
    /// sampled points using least-squares. Results in smaller geometric error
    /// but requires more computation.
    HighPrecision,
}

/// Direct shell conversion mode for `OpenSubdiv -> monstertruck`.
#[derive(Debug, Clone, Copy, Default, PartialEq, Eq)]
pub enum MonstertruckShellFlavor {
    /// Merge adjacent regular patches into the coarsest shell possible and
    /// append non-regular patch fallback surfaces.
    #[default]
    CoarseMerged,

    /// Export disconnected per-patch faces.
    DisconnectedPatches,

    /// Build stitched shared-edge topology from the raw patch layout.
    ///
    /// This path does not yet honor [`GregoryAccuracy`].
    StitchedPatches,
}

/// Options for direct `OpenSubdiv -> monstertruck` shell conversion.
#[derive(Debug, Clone, Copy)]
pub struct MonstertruckShellOptions {
    /// How to handle extraordinary vertices.
    pub gregory_accuracy: GregoryAccuracy,

    /// Tolerance for coarse merged regular-patch joining.
    pub stitch_tolerance: f64,

    /// Which shell construction path to use.
    pub flavor: MonstertruckShellFlavor,
}

impl Default for MonstertruckShellOptions {
    fn default() -> Self {
        Self {
            gregory_accuracy: GregoryAccuracy::BSplineEndCaps,
            stitch_tolerance: 1e-6,
            flavor: MonstertruckShellFlavor::CoarseMerged,
        }
    }
}

/// Compatibility options for STEP-oriented callers.
///
/// Controls how OpenSubdiv patches are converted to monstertruck B-spline
/// surfaces for STEP file output.
#[derive(Debug, Clone)]
pub struct StepExportOptions {
    /// How to handle extraordinary vertices (default: BSplineEndCaps).
    ///
    /// See [`GregoryAccuracy`] for details on each option.
    pub gregory_accuracy: GregoryAccuracy,

    /// Tolerance for vertex/edge stitching (default: 1e-6).
    ///
    /// When `stitch_edges` is true, vertices closer than this tolerance
    /// are considered the same vertex for B-rep construction.
    pub stitch_tolerance: f64,

    /// Create shared edges between adjacent patches (default: false).
    ///
    /// When true, builds a proper B-rep with shared vertices and edges
    /// between adjacent patches. When false, creates disconnected faces
    /// (simpler but less suitable for CAD operations).
    pub stitch_edges: bool,

    /// Use superpatch merging to combine adjacent regular patches.
    ///
    /// When true (default), adjacent regular quad patches are merged into
    /// larger B-spline surfaces with shared control vertices. This is
    /// critical for creased models where subdivision creates many small
    /// patches that can be efficiently combined.
    pub use_superpatches: bool,
}

impl Default for StepExportOptions {
    fn default() -> Self {
        Self {
            gregory_accuracy: GregoryAccuracy::BSplineEndCaps,
            stitch_tolerance: 1e-6,
            stitch_edges: false,
            use_superpatches: true,
        }
    }
}

impl From<StepExportOptions> for MonstertruckShellOptions {
    fn from(options: StepExportOptions) -> Self {
        let flavor = if options.use_superpatches {
            MonstertruckShellFlavor::CoarseMerged
        } else if options.stitch_edges {
            MonstertruckShellFlavor::StitchedPatches
        } else {
            MonstertruckShellFlavor::DisconnectedPatches
        };

        Self {
            gregory_accuracy: options.gregory_accuracy,
            stitch_tolerance: options.stitch_tolerance,
            flavor,
        }
    }
}

/// Apply OpenSubdiv regular-patch boundary adjustments to control points.
///
/// OpenSubdiv evaluates regular patches by first computing uniform bicubic
/// B-spline weights, then applying `adjustBSplineBoundaryWeights` based on the
/// boundary mask (see far/patchBasis.cpp). To export an equivalent surface
/// using the plain uniform basis, we transpose that weight transform and apply
/// it to the control points instead.
fn adjust_regular_control_points(
    control_matrix: Vec<Vec<Point3>>,
    boundary_mask: i32,
) -> Vec<Vec<Point3>> {
    if boundary_mask == 0 {
        return control_matrix;
    }

    // Flatten 4×4 control matrix in the same order as OpenSubdiv weights:
    // w[4*i + j] = sWeights[j] * tWeights[i], where i is v-row, j is u-col.
    let cps: Vec<_> = control_matrix
        .iter()
        .flat_map(|row| row.iter().copied())
        .collect();

    // Build transformation matrix (16×16) starting as identity, then run the
    // same boundary adjustments on its rows as adjustBSplineBoundaryWeights.
    let mut trans = [[0.0f64; 16]; 16];
    for (i, row) in trans.iter_mut().enumerate() {
        row[i] = 1.0;
    }

    // bit 0: v-min (bottom)
    if (boundary_mask & 0b0001) != 0 {
        for i in 0..4 {
            let row0 = trans[i];
            for k in 0..16 {
                trans[i + 8][k] -= row0[k];
                trans[i + 4][k] += row0[k] * 2.0;
            }
            trans[i] = [0.0; 16];
        }
    }

    // bit 1: u-max (right)
    if (boundary_mask & 0b0010) != 0 {
        for i in (0..16).step_by(4) {
            let row3 = trans[i + 3];
            for k in 0..16 {
                trans[i + 1][k] -= row3[k];
                trans[i + 2][k] += row3[k] * 2.0;
            }
            trans[i + 3] = [0.0; 16];
        }
    }

    // bit 2: v-max (top)
    if (boundary_mask & 0b0100) != 0 {
        for i in 0..4 {
            let row12 = trans[i + 12];
            for k in 0..16 {
                trans[i + 4][k] -= row12[k];
                trans[i + 8][k] += row12[k] * 2.0;
            }
            trans[i + 12] = [0.0; 16];
        }
    }

    // bit 3: u-min (left)
    if (boundary_mask & 0b1000) != 0 {
        for i in (0..16).step_by(4) {
            let row0 = trans[i];
            for k in 0..16 {
                trans[i + 2][k] -= row0[k];
                trans[i + 1][k] += row0[k] * 2.0;
            }
            trans[i] = [0.0; 16];
        }
    }

    // Apply transpose(trans) to control points: P' = trans^T * P.
    let mut new_cps = vec![Point3::origin(); 16];
    for new_idx in 0..16 {
        let mut acc = Point3::origin();
        for old_idx in 0..16 {
            let w = trans[old_idx][new_idx];
            if w != 0.0 {
                let p = cps[old_idx];
                acc.x += p.x * w;
                acc.y += p.y * w;
                acc.z += p.z * w;
            }
        }
        new_cps[new_idx] = acc;
    }

    // Unflatten back to 4×4.
    let mut adjusted = vec![vec![Point3::origin(); 4]; 4];
    for row in 0..4 {
        for col in 0..4 {
            adjusted[row][col] = new_cps[row * 4 + col];
        }
    }
    adjusted
}

/// A wrapper around a single patch with its associated data
pub struct PatchRef<'a> {
    pub patch_table: &'a PatchTable,
    pub patch_index: usize,
    pub control_points: &'a [[f32; 3]],
}

/// A wrapper around patches with control points for conversion
pub struct PatchTableWithControlPointsRef<'a> {
    pub patch_table: &'a PatchTable,
    pub control_points: &'a [[f32; 3]],
}

impl<'a> PatchRef<'a> {
    /// Create a new patch reference
    pub fn new(
        patch_table: &'a PatchTable,
        patch_index: usize,
        control_points: &'a [[f32; 3]],
    ) -> Self {
        Self {
            patch_table,
            patch_index,
            control_points,
        }
    }

    /// Get the patch type for this patch.
    pub fn patch_type(&self) -> std::result::Result<PatchType, MonstertruckError> {
        let (_, _, patch_type) = self.patch_info()?;
        Ok(patch_type)
    }

    /// Check if this patch is a Gregory patch (at an extraordinary vertex).
    pub fn is_gregory(&self) -> bool {
        matches!(
            self.patch_type(),
            Ok(PatchType::GregoryBasis) | Ok(PatchType::GregoryTriangle)
        )
    }

    /// Check if this patch is a regular B-spline patch.
    pub fn is_regular(&self) -> bool {
        matches!(self.patch_type(), Ok(PatchType::Regular))
    }

    /// Get the boundary mask for this patch.
    ///
    /// Returns a bitmask indicating which edges are boundaries (including
    /// infinite creases):
    /// - bit 0 (1): v-min edge (bottom)
    /// - bit 1 (2): u-max edge (right)
    /// - bit 2 (4): v-max edge (top)
    /// - bit 3 (8): u-min edge (left)
    pub fn boundary_mask(&self) -> i32 {
        if let Ok((array_index, local_index, _)) = self.patch_info() {
            self.patch_table
                .patch_param(array_index, local_index)
                .map(|p| p.boundary())
                .unwrap_or(0)
        } else {
            0
        }
    }

    /// Get patch array information
    fn patch_info(&self) -> std::result::Result<(usize, usize, PatchType), MonstertruckError> {
        let mut current_index = self.patch_index;

        for array_idx in 0..self.patch_table.patch_array_count() {
            let array_size = self.patch_table.patch_array_patch_count(array_idx);
            if current_index < array_size {
                let desc = self
                    .patch_table
                    .patch_array_descriptor(array_idx)
                    .ok_or(MonstertruckError::InvalidControlPoints)?;
                return Ok((array_idx, current_index, desc.patch_type()));
            }
            current_index -= array_size;
        }

        eprintln!("Failed to find patch {} in patch table", self.patch_index);
        Err(MonstertruckError::InvalidControlPoints)
    }

    /// Extract control points for this patch
    fn control_points(&self) -> std::result::Result<Vec<Vec<Point3>>, MonstertruckError> {
        let (array_index, local_index, patch_type) = self.patch_info()?;
        let boundary_mask = self.boundary_mask();

        // AIDEV-NOTE: Gregory patch support
        // Currently we only support Regular B-spline patches and Gregory patches.
        // Gregory patches are used at extraordinary vertices (valence != 4).
        // For now, we approximate Gregory patches as B-spline patches.
        match patch_type {
            PatchType::Regular => {
                self.extract_regular_patch_control_points(array_index, local_index, boundary_mask)
            }
            PatchType::GregoryBasis => {
                self.extract_gregory_basis_patch_control_points(array_index, local_index)
            }
            PatchType::GregoryTriangle => {
                self.extract_gregory_triangle_patch_control_points(array_index, local_index)
            }
            _ => Err(MonstertruckError::UnsupportedPatchType(patch_type)),
        }
    }

    /// Extract control points for a regular B-spline patch
    fn extract_regular_patch_control_points(
        &self,
        array_index: usize,
        local_index: usize,
        boundary_mask: i32,
    ) -> std::result::Result<Vec<Vec<Point3>>, MonstertruckError> {
        const REGULAR_PATCH_SIZE: usize = 4;
        let desc = self
            .patch_table
            .patch_array_descriptor(array_index)
            .ok_or(MonstertruckError::InvalidControlPoints)?;

        if desc.control_vertex_count() != REGULAR_PATCH_SIZE * REGULAR_PATCH_SIZE {
            return Err(MonstertruckError::InvalidControlPoints);
        }

        let cv_indices = self
            .patch_table
            .patch_array_vertices(array_index)
            .ok_or(MonstertruckError::InvalidControlPoints)?;

        let start = local_index * desc.control_vertex_count();
        if start + desc.control_vertex_count() > cv_indices.len() {
            eprintln!(
                "Patch {} (array {}, local {}): start={}, cvs_needed={}, available={}",
                self.patch_index,
                array_index,
                local_index,
                start,
                desc.control_vertex_count(),
                cv_indices.len()
            );
            return Err(MonstertruckError::InvalidControlPoints);
        }
        let patch_cvs = &cv_indices[start..start + desc.control_vertex_count()];

        let mut control_matrix =
            vec![vec![Point3::origin(); REGULAR_PATCH_SIZE]; REGULAR_PATCH_SIZE];

        for (i, &cv_idx) in patch_cvs.iter().enumerate() {
            let row = i / REGULAR_PATCH_SIZE;
            let col = i % REGULAR_PATCH_SIZE;

            let idx: usize = cv_idx.into();
            if idx >= self.control_points.len() {
                eprintln!(
                    "Patch {}: Control vertex index {} is out of bounds (max {})",
                    self.patch_index,
                    idx,
                    self.control_points.len() - 1
                );
                return Err(MonstertruckError::InvalidControlPoints);
            }

            let cp = &self.control_points[idx];
            control_matrix[row][col] = Point3::new(cp[0] as f64, cp[1] as f64, cp[2] as f64);
        }

        // AIDEV-NOTE: Check if this patch is adjacent to an extraordinary vertex
        // If so, we may need to adjust the control points to ensure proper continuity
        // This is a workaround for when OpenSubdiv generates Regular patches instead of
        // Gregory patches

        Ok(adjust_regular_control_points(control_matrix, boundary_mask))
    }

    /// Extract control points for a Gregory basis patch (20 control points)
    fn extract_gregory_basis_patch_control_points(
        &self,
        _array_index: usize,
        _local_index: usize,
    ) -> std::result::Result<Vec<Vec<Point3>>, MonstertruckError> {
        // AIDEV-NOTE: Gregory basis patch approximation
        // Gregory basis patches have 20 control points arranged in a special pattern.
        // For now, we evaluate the patch at a 4x4 grid to create an approximation.
        // This is not ideal but allows us to export something at extraordinary
        // vertices.

        // Evaluate the patch at 16 points to create a 4x4 control point grid
        let mut control_matrix = vec![vec![Point3::origin(); 4]; 4];

        #[allow(clippy::needless_range_loop)]
        for i in 0..4 {
            for j in 0..4 {
                // Map to parameter space [0,1]
                let u = i as f32 / 3.0;
                let v = j as f32 / 3.0;

                // Evaluate the patch at this parameter location
                if let Some(result) =
                    self.patch_table
                        .evaluate_point(self.patch_index, u, v, self.control_points)
                {
                    control_matrix[i][j] = Point3::new(
                        result.point[0] as f64,
                        result.point[1] as f64,
                        result.point[2] as f64,
                    );
                } else {
                    return Err(MonstertruckError::EvaluationFailed);
                }
            }
        }

        Ok(control_matrix)
    }

    /// Extract control points for a Gregory triangle patch (18 control points)
    fn extract_gregory_triangle_patch_control_points(
        &self,
        _array_index: usize,
        _local_index: usize,
    ) -> std::result::Result<Vec<Vec<Point3>>, MonstertruckError> {
        // AIDEV-NOTE: Gregory triangle patch approximation
        // Gregory triangle patches have 18 control points for triangular domains.
        // For now, we evaluate the patch at a 4x4 grid to create a quad approximation.
        // This converts the triangular patch to a degenerate quad patch.

        // Evaluate the patch at 16 points to create a 4x4 control point grid
        let mut control_matrix = vec![vec![Point3::origin(); 4]; 4];

        #[allow(clippy::needless_range_loop)]
        for i in 0..4 {
            for j in 0..4 {
                // For triangular patches, we need to ensure u + v <= 1
                let u = i as f32 / 3.0;
                let v = j as f32 / 3.0;

                // If we're outside the triangular domain, collapse to the edge
                let (u_eval, v_eval) = if u + v > 1.0 {
                    // Project back onto the triangle edge
                    let sum = u + v;
                    (u / sum, v / sum)
                } else {
                    (u, v)
                };

                // Evaluate the patch at this parameter location
                if let Some(result) = self.patch_table.evaluate_point(
                    self.patch_index,
                    u_eval,
                    v_eval,
                    self.control_points,
                ) {
                    control_matrix[i][j] = Point3::new(
                        result.point[0] as f64,
                        result.point[1] as f64,
                        result.point[2] as f64,
                    );
                } else {
                    return Err(MonstertruckError::EvaluationFailed);
                }
            }
        }

        Ok(control_matrix)
    }

    /// Convert a Gregory patch to a B-spline surface using high-precision
    /// sampling.
    ///
    /// This evaluates the Gregory patch at an 8×8 grid and creates a B-spline
    /// surface that approximates the original patch more accurately than the
    /// standard 4×4 sampling.
    ///
    /// # How it works
    ///
    /// 1. Evaluates the Gregory patch at 64 points (8×8 grid)
    /// 2. Creates a B-spline surface with 8×8 control points
    /// 3. Uses uniform knot vectors appropriate for the larger control grid
    ///
    /// Note: This is still an approximation since Gregory patches cannot be
    /// exactly represented as B-splines. However, the denser sampling captures
    /// more of the surface curvature at extraordinary vertices.
    pub fn to_bspline_high_precision(&self) -> Result<BsplineSurface<Point3>> {
        const GRID_SIZE: usize = 8;

        // Evaluate the Gregory patch at an 8×8 grid.
        let samples: Vec<Vec<Point3>> = (0..GRID_SIZE)
            .map(|i| {
                let u = i as f32 / (GRID_SIZE - 1) as f32;
                (0..GRID_SIZE)
                    .map(|j| {
                        let v = j as f32 / (GRID_SIZE - 1) as f32;
                        self.patch_table
                            .evaluate_point(self.patch_index, u, v, self.control_points)
                            .map(|result| {
                                Point3::new(
                                    result.point[0] as f64,
                                    result.point[1] as f64,
                                    result.point[2] as f64,
                                )
                            })
                            .ok_or(MonstertruckError::EvaluationFailed)
                    })
                    .collect::<Result<Vec<_>>>()
            })
            .collect::<Result<Vec<_>>>()?;

        // Create knot vectors for an 8×8 control point grid with degree 3.
        // For n control points and degree p, we need n + p + 1 knots.
        // With 8 control points and degree 3, we need 12 knots.
        // Use uniform spacing for smooth C² continuity.
        let knots = KnotVector::from(vec![
            -3.0, -2.0, -1.0, 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0,
        ]);

        Ok(BsplineSurface::new((knots.clone(), knots), samples))
    }
}

/// Convert a regular B-spline patch to a monstertruck BsplineSurface
impl<'a> TryFrom<PatchRef<'a>> for BsplineSurface<Point3> {
    type Error = MonstertruckError;

    fn try_from(patch: PatchRef<'a>) -> std::result::Result<Self, Self::Error> {
        let control_matrix = patch.control_points()?;

        // AIDEV-NOTE: OpenSubdiv B-spline patch knot vectors
        // Export patches exactly as OpenSubdiv defines them: uniform knots with
        // phantom rows/columns. Boundary masks affect evaluation inside OSD, but
        // the control net is authored for this uniform basis, so we keep it.
        let u_knots = KnotVector::from(vec![-3.0, -2.0, -1.0, 0.0, 1.0, 2.0, 3.0, 4.0]);
        let v_knots = KnotVector::from(vec![-3.0, -2.0, -1.0, 0.0, 1.0, 2.0, 3.0, 4.0]);

        Ok(BsplineSurface::new((u_knots, v_knots), control_matrix))
    }
}

/// Create a Face with explicit B-spline boundary curves from a 4×4 control
/// matrix.
///
/// This helper function extracts the boundary control points from a bicubic
/// B-spline patch and creates proper boundary curves for the face. This is
/// required for some STEP viewers that need explicit boundary definitions.
#[cfg(feature = "monstertruck_export_boundary")]
fn create_face_with_boundary(
    control_matrix: &[Vec<Point3>],
    surface: BsplineSurface<Point3>,
) -> Face {
    use monstertruck::geometry::prelude::BsplineCurve;

    // AIDEV-NOTE: Boundary control point extraction.
    // For OpenSubdiv B-spline patches with uniform knot vectors,
    // we extract the boundary control points to define B-spline boundary curves.

    // Bottom edge (row 0): (0,0), (0,1), (0,2), (0,3).
    let bottom_cps = vec![
        control_matrix[0][0],
        control_matrix[0][1],
        control_matrix[0][2],
        control_matrix[0][3],
    ];

    // Right edge (column 3): (0,3), (1,3), (2,3), (3,3).
    let right_cps = vec![
        control_matrix[0][3],
        control_matrix[1][3],
        control_matrix[2][3],
        control_matrix[3][3],
    ];

    // Top edge (row 3, reversed): (3,3), (3,2), (3,1), (3,0).
    let top_cps = vec![
        control_matrix[3][3],
        control_matrix[3][2],
        control_matrix[3][1],
        control_matrix[3][0],
    ];

    // Left edge (column 0, reversed): (3,0), (2,0), (1,0), (0,0).
    let left_cps = vec![
        control_matrix[3][0],
        control_matrix[2][0],
        control_matrix[1][0],
        control_matrix[0][0],
    ];

    // Create the same knot vector as the surface.
    let edge_knots = KnotVector::from(vec![-3.0, -2.0, -1.0, 0.0, 1.0, 2.0, 3.0, 4.0]);

    // Create B-spline curves for edges.
    let bottom_curve = BsplineCurve::new(edge_knots.clone(), bottom_cps);
    let right_curve = BsplineCurve::new(edge_knots.clone(), right_cps);
    let top_curve = BsplineCurve::new(edge_knots.clone(), top_cps);
    let left_curve = BsplineCurve::new(edge_knots, left_cps);

    // Create vertices at the corner positions.
    let v00 = Vertex::new(control_matrix[0][0]); // Bottom-left.
    let v10 = Vertex::new(control_matrix[0][3]); // Bottom-right.
    let v11 = Vertex::new(control_matrix[3][3]); // Top-right.
    let v01 = Vertex::new(control_matrix[3][0]); // Top-left.

    // Create edges with B-spline curves.
    let e0 = Edge::new(&v00, &v10, Curve::BsplineCurve(bottom_curve));
    let e1 = Edge::new(&v10, &v11, Curve::BsplineCurve(right_curve));
    let e2 = Edge::new(&v11, &v01, Curve::BsplineCurve(top_curve));
    let e3 = Edge::new(&v01, &v00, Curve::BsplineCurve(left_curve));

    // Create wire and face.
    let wire = Wire::from(vec![e0, e1, e2, e3]);
    Face::new(vec![wire], Surface::BsplineSurface(surface))
}

/// Convert all regular patches to B-spline surfaces
impl<'a> TryFrom<PatchTableWithControlPointsRef<'a>> for Vec<BsplineSurface<Point3>> {
    type Error = MonstertruckError;

    fn try_from(
        patches: PatchTableWithControlPointsRef<'a>,
    ) -> std::result::Result<Self, Self::Error> {
        let mut surfaces = Vec::new();
        let mut patch_index = 0;

        for array_idx in 0..patches.patch_table.patch_array_count() {
            if let Some(desc) = patches.patch_table.patch_array_descriptor(array_idx) {
                let patch_type = desc.patch_type();
                // Handle Regular, GregoryBasis, and GregoryTriangle patches
                if matches!(
                    patch_type,
                    PatchType::Regular | PatchType::GregoryBasis | PatchType::GregoryTriangle
                ) {
                    for _ in 0..patches.patch_table.patch_array_patch_count(array_idx) {
                        let patch =
                            PatchRef::new(patches.patch_table, patch_index, patches.control_points);
                        match BsplineSurface::try_from(patch) {
                            Ok(surface) => surfaces.push(surface),
                            Err(e) => eprintln!(
                                "Failed to convert patch {} (type {:?}): {:?}",
                                patch_index, patch_type, e
                            ),
                        }
                        patch_index += 1;
                    }
                } else {
                    eprintln!(
                        "Skipping patch array {} with type {:?} ({} patches)",
                        array_idx,
                        patch_type,
                        patches.patch_table.patch_array_patch_count(array_idx)
                    );
                    patch_index += patches.patch_table.patch_array_patch_count(array_idx);
                }
            }
        }

        if surfaces.is_empty() {
            // No non-regular patches to convert; return an empty list so callers
            // can rely on BFR-regular surfaces alone.
            Ok(surfaces)
        } else {
            Ok(surfaces)
        }
    }
}

/// Convert only non-regular patches to B-spline surfaces (skip regular to allow
/// BFR substitution).
pub fn patch_table_surfaces_non_regular(
    patch_table: &PatchTable,
    control_points: &[[f32; 3]],
) -> Result<Vec<BsplineSurface<Point3>>> {
    patch_table_surfaces_non_regular_with_options(
        patch_table,
        control_points,
        GregoryAccuracy::BSplineEndCaps,
    )
}

/// Convert only non-regular patches to B-spline surfaces with configurable
/// Gregory handling.
pub fn patch_table_surfaces_non_regular_with_options(
    patch_table: &PatchTable,
    control_points: &[[f32; 3]],
    gregory_accuracy: GregoryAccuracy,
) -> Result<Vec<BsplineSurface<Point3>>> {
    let mut surfaces = Vec::new();
    let mut patch_index = 0;

    for array_idx in 0..patch_table.patch_array_count() {
        if let Some(desc) = patch_table.patch_array_descriptor(array_idx) {
            let patch_type = desc.patch_type();
            let num_patches = patch_table.patch_array_patch_count(array_idx);

            if patch_type == PatchType::Regular {
                patch_index += num_patches;
                continue;
            }

            if matches!(
                patch_type,
                PatchType::GregoryBasis
                    | PatchType::GregoryTriangle
                    | PatchType::Gregory
                    | PatchType::GregoryBoundary
                    | PatchType::GregoryCorner
            ) {
                for _ in 0..num_patches {
                    let patch = PatchRef::new(patch_table, patch_index, control_points);
                    let surface = if patch.is_gregory()
                        && gregory_accuracy == GregoryAccuracy::HighPrecision
                    {
                        patch.to_bspline_high_precision()
                    } else {
                        BsplineSurface::try_from(patch)
                    };
                    match surface {
                        Ok(surface) => surfaces.push(surface),
                        Err(e) => eprintln!(
                            "Failed to convert non-regular patch {} (type {:?}): {:?}",
                            patch_index, patch_type, e
                        ),
                    }
                    patch_index += 1;
                }
            } else {
                // Unsupported types skipped
                patch_index += num_patches;
            }
        }
    }

    Ok(surfaces)
}

/// Merge adjacent regular patches into larger bicubic superpatches to reduce
/// patch counts while keeping curvature.
pub fn superpatch_surfaces(
    patch_table: &PatchTable,
    control_points: &[[f32; 3]],
    tol: f64,
) -> Result<Vec<BsplineSurface<Point3>>> {
    const DEGREE: usize = 3;

    #[derive(Clone)]
    struct RegPatch {
        _index: usize,
        control: Vec<Vec<Point3>>, // 4x4, row-major (v-major)
        boundary_mask: i32,        // Boundary flags from PatchParam
    }

    #[derive(Default, Clone)]
    struct Adjacency {
        right: Option<usize>,
        bottom: Option<usize>,
    }

    #[derive(Clone)]
    struct Superpatch {
        control: Vec<Vec<Point3>>, // u-major
        width_cells: usize,
        height_cells: usize,
        origin_x: i32,
        origin_y: i32,
        component: usize,
        boundary_mask: i32, // Boundary flags for single-patch superpatches
    }

    /// Edge data from a superpatch: (left, right, bottom, top).
    type SuperpatchEdges = (Vec<Point3>, Vec<Point3>, Vec<Point3>, Vec<Point3>);

    fn edge_row(control: &[Vec<Point3>], edge: &str) -> [Point3; 4] {
        match edge {
            "bottom" => [control[3][0], control[3][1], control[3][2], control[3][3]],
            "top" => [control[0][0], control[0][1], control[0][2], control[0][3]],
            "left" => [control[0][0], control[1][0], control[2][0], control[3][0]],
            "right" => [control[0][3], control[1][3], control[2][3], control[3][3]],
            _ => unreachable!(),
        }
    }

    fn rows_match(a: &[Point3; 4], b: &[Point3; 4], tol: f64) -> bool {
        a.iter()
            .zip(b.iter())
            .all(|(p, q)| (p - q).magnitude2() <= tol * tol)
    }

    fn superpatch_edges(sp: &Superpatch) -> SuperpatchEdges {
        let u_max = sp.control.len().saturating_sub(1);
        let v_max = sp
            .control
            .first()
            .map(|c| c.len().saturating_sub(1))
            .unwrap_or(0);

        let left = sp.control.first().cloned().unwrap_or_default();
        let right = sp.control.get(u_max).cloned().unwrap_or_default();

        let (top, bottom): (Vec<_>, Vec<_>) = sp
            .control
            .iter()
            .map(|col| {
                (
                    *col.first().unwrap_or(&Point3::origin()),
                    *col.get(v_max).unwrap_or(&Point3::origin()),
                )
            })
            .unzip();

        (left, right, bottom, top)
    }

    fn edges_match(a: &[Point3], b: &[Point3], tol: f64) -> bool {
        if a.len() != b.len() {
            return false;
        }
        a.iter()
            .zip(b.iter())
            .all(|(p, q)| (*p - *q).magnitude2() <= tol * tol)
    }

    fn merge_horizontal(a: &Superpatch, b: &Superpatch) -> Superpatch {
        // Concatenate control columns: all of a, then b (skipping first shared column).
        let control: Vec<Vec<_>> = a
            .control
            .iter()
            .chain(b.control.iter().skip(1))
            .cloned()
            .collect();

        // Combine boundary masks: keep a's left, top, bottom; keep b's right.
        // bit 0 (v-min/bottom), bit 1 (u-max/right), bit 2 (v-max/top), bit 3
        // (u-min/left)
        let boundary_mask = (a.boundary_mask & 0b1101) | (b.boundary_mask & 0b0010);

        Superpatch {
            control,
            width_cells: a.width_cells + b.width_cells,
            height_cells: a.height_cells,
            origin_x: a.origin_x.min(b.origin_x),
            origin_y: a.origin_y,
            component: a.component,
            boundary_mask,
        }
    }

    fn merge_vertical(top: &Superpatch, bottom: &Superpatch) -> Superpatch {
        let control: Vec<Vec<_>> = top
            .control
            .iter()
            .zip(bottom.control.iter())
            .map(|(top_col, bottom_col)| {
                top_col
                    .iter()
                    .chain(bottom_col.iter().skip(1))
                    .copied()
                    .collect()
            })
            .collect();

        // Combine boundary masks: keep top's left, right, top; keep bottom's bottom.
        // bit 0 (v-min/bottom), bit 1 (u-max/right), bit 2 (v-max/top), bit 3
        // (u-min/left)
        let boundary_mask = (top.boundary_mask & 0b1110) | (bottom.boundary_mask & 0b0001);

        Superpatch {
            control,
            width_cells: top.width_cells,
            height_cells: top.height_cells + bottom.height_cells,
            origin_x: top.origin_x,
            origin_y: top.origin_y,
            component: top.component,
            boundary_mask,
        }
    }

    // Collect regular patches and build adjacency.
    let mut regular = Vec::<RegPatch>::new();
    let mut patch_index = 0usize;

    for array_idx in 0..patch_table.patch_array_count() {
        if let Some(desc) = patch_table.patch_array_descriptor(array_idx) {
            let patch_type = desc.patch_type();
            let num_patches = patch_table.patch_array_patch_count(array_idx);
            if patch_type == PatchType::Regular {
                for _ in 0..num_patches {
                    let patch = PatchRef::new(patch_table, patch_index, control_points);
                    if let Ok(cm) = patch.control_points() {
                        regular.push(RegPatch {
                            _index: patch_index,
                            control: cm,
                            boundary_mask: patch.boundary_mask(),
                        });
                    }
                    patch_index += 1;
                }
            } else {
                patch_index += num_patches;
            }
        }
    }

    let mut adjacency = vec![Adjacency::default(); regular.len()];

    // AIDEV-NOTE: Boundary edge constants for superpatch merging.
    // Patches with boundary edges (including infinite creases) have clamped knot
    // vectors and should NOT be merged across those edges.
    const BOUNDARY_V_MIN: i32 = 0b0001; // bottom edge (v=0)
    const BOUNDARY_U_MAX: i32 = 0b0010; // right edge (u=1)
    const BOUNDARY_V_MAX: i32 = 0b0100; // top edge (v=1)
    const BOUNDARY_U_MIN: i32 = 0b1000; // left edge (u=0)

    for i in 0..regular.len() {
        let r_i = &regular[i];
        let bottom_i = edge_row(&r_i.control, "bottom");
        let right_i = edge_row(&r_i.control, "right");

        for (j, r_j) in regular.iter().enumerate() {
            if i == j {
                continue;
            }
            let top_j = edge_row(&r_j.control, "top");
            let left_j = edge_row(&r_j.control, "left");

            // Bottom adjacency: i's bottom connects to j's top.
            // Skip if either edge is a boundary (clamped knots incompatible).
            let i_bottom_boundary = r_i.boundary_mask & BOUNDARY_V_MIN != 0;
            let j_top_boundary = r_j.boundary_mask & BOUNDARY_V_MAX != 0;
            if adjacency[i].bottom.is_none()
                && !i_bottom_boundary
                && !j_top_boundary
                && rows_match(&bottom_i, &top_j, tol)
            {
                adjacency[i].bottom = Some(j);
            }

            // Right adjacency: i's right connects to j's left.
            // Skip if either edge is a boundary (clamped knots incompatible).
            let i_right_boundary = r_i.boundary_mask & BOUNDARY_U_MAX != 0;
            let j_left_boundary = r_j.boundary_mask & BOUNDARY_U_MIN != 0;
            if adjacency[i].right.is_none()
                && !i_right_boundary
                && !j_left_boundary
                && rows_match(&right_i, &left_j, tol)
            {
                adjacency[i].right = Some(j);
            }
        }
    }

    // Reverse adjacency to walk both directions.
    let mut left_of = vec![None; regular.len()];
    let mut top_of = vec![None; regular.len()];
    for (i, adj) in adjacency.iter().enumerate() {
        if let Some(r) = adj.right {
            left_of[r] = Some(i);
        }
        if let Some(b) = adj.bottom {
            top_of[b] = Some(i);
        }
    }

    // Assign grid coordinates to patches using adjacency.
    let mut coords = vec![None; regular.len()];
    let mut components = vec![None; regular.len()];
    let mut component_id = 0usize;

    for start in 0..regular.len() {
        if coords[start].is_some() {
            continue;
        }
        coords[start] = Some((0i32, 0i32));
        components[start] = Some(component_id);
        let mut queue = std::collections::VecDeque::new();
        queue.push_back(start);
        while let Some(i) = queue.pop_front() {
            let (x, y) = coords[i].unwrap();
            if let Some(r) = adjacency[i].right {
                if coords[r].is_none() {
                    coords[r] = Some((x + 1, y));
                    components[r] = Some(component_id);
                    queue.push_back(r);
                }
            }
            if let Some(b) = adjacency[i].bottom {
                if coords[b].is_none() {
                    coords[b] = Some((x, y + 1));
                    components[b] = Some(component_id);
                    queue.push_back(b);
                }
            }
            if let Some(l) = left_of[i] {
                if coords[l].is_none() {
                    coords[l] = Some((x - 1, y));
                    components[l] = Some(component_id);
                    queue.push_back(l);
                }
            }
            if let Some(t) = top_of[i] {
                if coords[t].is_none() {
                    coords[t] = Some((x, y - 1));
                    components[t] = Some(component_id);
                    queue.push_back(t);
                }
            }
        }
        component_id += 1;
    }

    // Build coordinate maps per component.
    let mut comp_maps: std::collections::HashMap<
        usize,
        std::collections::HashMap<(i32, i32), usize>,
    > = std::collections::HashMap::new();
    for (idx, coord) in coords.iter().enumerate() {
        if let (Some((x, y)), Some(comp)) = (coord, components[idx]) {
            comp_maps.entry(comp).or_default().insert((*x, *y), idx);
        }
    }

    // Build superpatch rectangles from grid coordinates.
    let mut superpatches = Vec::<Superpatch>::new();
    let mut visited = vec![false; regular.len()];

    for idx in 0..regular.len() {
        if visited[idx] {
            continue;
        }
        let (x0, y0) = coords[idx].unwrap_or((0, 0));
        let comp = components[idx].unwrap_or(usize::MAX);
        let Some(coord_map) = comp_maps.get(&comp) else {
            continue;
        };

        // Determine rectangle width.
        let mut width = 0usize;
        while coord_map.contains_key(&(x0 + width as i32, y0)) {
            width += 1;
        }
        // Determine rectangle height.
        let mut height = 0usize;
        'rows: loop {
            let y = y0 + height as i32;
            for u in 0..width {
                if !coord_map.contains_key(&(x0 + u as i32, y)) {
                    break 'rows;
                }
            }
            height += 1;
        }

        if width == 0 || height == 0 {
            continue;
        }

        let ctrl_u = width * DEGREE + 1;
        let ctrl_v = height * DEGREE + 1;
        let mut grid: Vec<Vec<Option<Point3>>> = vec![vec![None; ctrl_u]; ctrl_v];
        let mut valid = true;

        for v_off in 0..height {
            for u_off in 0..width {
                let x = x0 + u_off as i32;
                let y = y0 + v_off as i32;
                if let Some(&p_idx) = coord_map.get(&(x, y)) {
                    let patch = &regular[p_idx];
                    let off_u = u_off * DEGREE;
                    let off_v = v_off * DEGREE;
                    for i in 0..4 {
                        for j in 0..4 {
                            let dest = &mut grid[off_v + i][off_u + j];
                            let value = patch.control[i][j];
                            if let Some(existing) = dest {
                                if (*existing - value).magnitude2() > tol * tol {
                                    eprintln!(
                                        "Superpatch mismatch at patch {}, slot ({},{}): existing vs new differ",
                                        p_idx,
                                        off_v + i,
                                        off_u + j
                                    );
                                    valid = false;
                                }
                            } else {
                                *dest = Some(value);
                            }
                        }
                    }
                }
            }
        }

        if !valid {
            // Fall back to individual 1x1 patches for this rectangle to avoid corrupt
            // merges.
            for v_off in 0..height {
                for u_off in 0..width {
                    if let Some(&p_idx) = coord_map.get(&(x0 + u_off as i32, y0 + v_off as i32)) {
                        let patch = &regular[p_idx];
                        // Transpose from v-major to u-major.
                        let control: Vec<Vec<Point3>> = (0..4)
                            .map(|u| (0..4).map(|v| patch.control[v][u]).collect())
                            .collect();
                        superpatches.push(Superpatch {
                            control,
                            width_cells: 1,
                            height_cells: 1,
                            origin_x: x0 + u_off as i32,
                            origin_y: y0 + v_off as i32,
                            component: comp,
                            boundary_mask: patch.boundary_mask,
                        });
                        visited[p_idx] = true;
                    }
                }
            }
            continue;
        }

        // Transpose to u-major for monstertruck surfaces.
        let mut control_matrix = vec![vec![Point3::origin(); ctrl_v]; ctrl_u];
        for v in 0..ctrl_v {
            for u in 0..ctrl_u {
                control_matrix[u][v] = grid[v][u].unwrap_or_else(Point3::origin);
            }
        }

        // Compute combined boundary mask from outer edge patches.
        // bit 0 (v-min/bottom): from bottom row (v_off = height-1)
        // bit 1 (u-max/right): from right column (u_off = width-1)
        // bit 2 (v-max/top): from top row (v_off = 0)
        // bit 3 (u-min/left): from left column (u_off = 0)
        let mut combined_boundary = 0i32;
        for v_off in 0..height {
            for u_off in 0..width {
                if let Some(&p_idx) = coord_map.get(&(x0 + u_off as i32, y0 + v_off as i32)) {
                    visited[p_idx] = true;
                    let mask = regular[p_idx].boundary_mask;
                    // Left column contributes left boundary.
                    if u_off == 0 {
                        combined_boundary |= mask & 0b1000;
                    }
                    // Right column contributes right boundary.
                    if u_off == width - 1 {
                        combined_boundary |= mask & 0b0010;
                    }
                    // Top row contributes top boundary.
                    if v_off == 0 {
                        combined_boundary |= mask & 0b0100;
                    }
                    // Bottom row contributes bottom boundary.
                    if v_off == height - 1 {
                        combined_boundary |= mask & 0b0001;
                    }
                }
            }
        }

        superpatches.push(Superpatch {
            control: control_matrix,
            width_cells: width,
            height_cells: height,
            origin_x: x0,
            origin_y: y0,
            component: comp,
            boundary_mask: combined_boundary,
        });
    }

    // Hierarchical merging: always try to merge the smallest superpatches into
    // larger ones so power-of-two layouts can coalesce progressively.
    loop {
        // Sort ascending by area to prioritize smaller patches.
        superpatches.sort_by_key(|sp| sp.width_cells * sp.height_cells);

        let mut used = vec![false; superpatches.len()];
        let mut next = Vec::<Superpatch>::new();
        let mut merged_any = false;

        for i in 0..superpatches.len() {
            if used[i] {
                continue;
            }
            let sp_i = &superpatches[i];
            let (_left_i, right_i, bottom_i, _top_i) = superpatch_edges(sp_i);
            let mut merged = false;

            for j in (i + 1)..superpatches.len() {
                if used[j] {
                    continue;
                }
                let sp_j = &superpatches[j];
                let (left_j, _right_j, _bottom_j, top_j) = superpatch_edges(sp_j);

                // Horizontal merge: same component and height, touching in grid.
                if sp_i.component == sp_j.component
                    && sp_i.height_cells == sp_j.height_cells
                    && sp_i.origin_y == sp_j.origin_y
                    && sp_i.origin_x + sp_i.width_cells as i32 == sp_j.origin_x
                    && edges_match(&right_i, &left_j, tol)
                {
                    let combined = merge_horizontal(sp_i, sp_j);
                    used[i] = true;
                    used[j] = true;
                    merged = true;
                    merged_any = true;
                    next.push(combined);
                    break;
                }

                // Vertical merge: same component and width, touching in grid.
                if sp_i.component == sp_j.component
                    && sp_i.width_cells == sp_j.width_cells
                    && sp_i.origin_x == sp_j.origin_x
                    && sp_i.origin_y + sp_i.height_cells as i32 == sp_j.origin_y
                    && edges_match(&bottom_i, &top_j, tol)
                {
                    let combined = merge_vertical(sp_i, sp_j);
                    used[i] = true;
                    used[j] = true;
                    merged = true;
                    merged_any = true;
                    next.push(combined);
                    break;
                }
            }

            if !merged {
                used[i] = true;
                next.push(sp_i.clone());
            }
        }

        if !merged_any {
            superpatches = next;
            break;
        }
        superpatches = next;
    }

    let surfaces = superpatches
        .into_iter()
        .map(|sp| {
            let ctrl_u = sp.control.len();
            let ctrl_v = sp.control.first().map(|c| c.len()).unwrap_or(0);
            // Use uniform knot vectors - OpenSubdiv control points are computed for this.
            let knot_u: Vec<f64> = (0..(ctrl_u + DEGREE + 1))
                .map(|k| k as f64 - DEGREE as f64)
                .collect();
            let knot_v: Vec<f64> = (0..(ctrl_v + DEGREE + 1))
                .map(|k| k as f64 - DEGREE as f64)
                .collect();
            BsplineSurface::new(
                (KnotVector::from(knot_u), KnotVector::from(knot_v)),
                sp.control,
            )
        })
        .collect();

    Ok(surfaces)
}

// AIDEV-NOTE: Commented out full B-rep Shell implementation with shared edges
// This implementation creates a proper B-rep with shared vertices and edges,
// but for debugging we're using a simpler disconnected patch approach below.
/*
/// Convert patches to a complete Shell with shared topology
impl<'a> TryFrom<PatchTableWithControlPointsRef<'a>> for Shell {
    type Error = MonstertruckError;

    fn try_from(patches: PatchTableWithControlPoints<'a>) -> std::result::Result<Self, Self::Error> {
        let surfaces: Vec<BsplineSurface<Point3>> = patches.try_into()?;

        use std::collections::HashMap;
        use monstertruck::geometry::prelude::BsplineCurve;

        // AIDEV-NOTE: Create proper B-rep with shared vertices and edges
        // Following the pattern from monstertruck-topology's cube example, we need to:
        // 1. Create all vertices first
        // 2. Create all edges between vertices
        // 3. Build faces using these edges with proper orientation

        // Tolerance for position comparison
        const TOLERANCE: f64 = 1e-10;

        // First pass: collect all unique corner points and create vertices
        let mut vertex_map: HashMap<[i64; 3], Vertex> = HashMap::new();
        let mut surface_corners = Vec::new();

        for surface in &surfaces {
            // Get the four corner points
            let p00 = surface.evaluate(0.0, 0.0);
            let p10 = surface.evaluate(1.0, 0.0);
            let p11 = surface.evaluate(1.0, 1.0);
            let p01 = surface.evaluate(0.0, 1.0);

            // Get or create vertices
            let mut get_or_create_vertex = |point: Point3| -> Vertex {
                let key = [
                    (point.x / TOLERANCE).round() as i64,
                    (point.y / TOLERANCE).round() as i64,
                    (point.z / TOLERANCE).round() as i64,
                ];

                vertex_map.entry(key)
                    .or_insert_with(|| Vertex::new(point))
                    .clone()
            };

            let v00 = get_or_create_vertex(p00);
            let v10 = get_or_create_vertex(p10);
            let v11 = get_or_create_vertex(p11);
            let v01 = get_or_create_vertex(p01);

            surface_corners.push((v00, v10, v11, v01, p00, p10, p11, p01));
        }

        // Second pass: create all unique edges
        type EdgeKey = ([i64; 3], [i64; 3]);
        let mut edge_map: HashMap<EdgeKey, Edge> = HashMap::new();

        let make_edge_key = |p0: Point3, p1: Point3| -> EdgeKey {
            let k0 = [
                (p0.x / TOLERANCE).round() as i64,
                (p0.y / TOLERANCE).round() as i64,
                (p0.z / TOLERANCE).round() as i64,
            ];
            let k1 = [
                (p1.x / TOLERANCE).round() as i64,
                (p1.y / TOLERANCE).round() as i64,
                (p1.z / TOLERANCE).round() as i64,
            ];
            // Always order vertices consistently for the key
            if k0 <= k1 { (k0, k1) } else { (k1, k0) }
        };

        // Collect all edges from all patches
        for (v00, v10, v11, v01, p00, p10, p11, p01) in &surface_corners {
            // Helper to create or get edge
            let mut get_or_create_edge = |v0: &Vertex, v1: &Vertex, p0: Point3, p1: Point3| {
                let key = make_edge_key(p0, p1);
                edge_map.entry(key)
                    .or_insert_with(|| {
                        // Always create edge in consistent direction based on key
                        if make_edge_key(p0, p1) == (key.0, key.1) {
                            Edge::new(v0, v1, Curve::BsplineCurve(
                                BsplineCurve::new(KnotVector::bezier_knot(1), vec![p0, p1])
                            ))
                        } else {
                            Edge::new(v1, v0, Curve::BsplineCurve(
                                BsplineCurve::new(KnotVector::bezier_knot(1), vec![p1, p0])
                            ))
                        }
                    });
            };

            // Create all four edges for this patch
            get_or_create_edge(v00, v10, *p00, *p10);
            get_or_create_edge(v10, v11, *p10, *p11);
            get_or_create_edge(v11, v01, *p11, *p01);
            get_or_create_edge(v01, v00, *p01, *p00);
        }

        // Third pass: create faces using the shared edges
        let mut faces = Vec::new();

        for (i, (surface, (v00, v10, v11, v01, p00, p10, p11, p01))) in surfaces.into_iter().zip(surface_corners).enumerate() {
            // Get the edges for this face
            let bottom_edge = edge_map.get(&make_edge_key(p00, p10))
                .ok_or(MonstertruckError::InvalidControlPoints)?;
            let right_edge = edge_map.get(&make_edge_key(p10, p11))
                .ok_or(MonstertruckError::InvalidControlPoints)?;
            let top_edge = edge_map.get(&make_edge_key(p11, p01))
                .ok_or(MonstertruckError::InvalidControlPoints)?;
            let left_edge = edge_map.get(&make_edge_key(p01, p00))
                .ok_or(MonstertruckError::InvalidControlPoints)?;

            // Calculate face normal at the center to determine proper orientation
            let center_u = 0.5;
            let center_v = 0.5;
            let _center_pt = surface.evaluate(center_u, center_v);
            let du = surface.derivative_u(center_u, center_v);
            let dv = surface.derivative_v(center_u, center_v);
            let normal = du.cross(dv);

            // Compute the expected outward normal based on corner points
            // Using (p10-p00) x (p01-p00) which should point outward for CCW winding
            let edge1 = Vector3::new(p10.x - p00.x, p10.y - p00.y, p10.z - p00.z);
            let edge2 = Vector3::new(p01.x - p00.x, p01.y - p00.y, p01.z - p00.z);
            let expected_normal = edge1.cross(edge2);

            // Check if surface normal matches expected normal
            let dot = normal.dot(expected_normal);
            let needs_inversion = dot < 0.0;

            if needs_inversion {
                eprintln!("Warning: Face {} has inverted normal (dot = {})", i, dot);
            }

            // Determine proper orientation for each edge
            let bottom = if bottom_edge.front() == &v00 {
                bottom_edge.clone()
            } else {
                bottom_edge.inverse()
            };

            let right = if right_edge.front() == &v10 {
                right_edge.clone()
            } else {
                right_edge.inverse()
            };

            let top = if top_edge.front() == &v11 {
                top_edge.clone()
            } else {
                top_edge.inverse()
            };

            let left = if left_edge.front() == &v01 {
                left_edge.clone()
            } else {
                left_edge.inverse()
            };

            // Create wire and face
            let wire = if needs_inversion {
                // Reverse the edge order to flip the face normal
                Wire::from(vec![bottom.inverse(), left.inverse(), top.inverse(), right.inverse()])
            } else {
                Wire::from(vec![bottom, right, top, left])
            };

            let mut face = Face::new(vec![wire], Surface::BsplineSurface(surface));
            if needs_inversion {
                face.invert();
            }
            faces.push(face);
        }

        let shell = Shell::from(faces);

        Ok(shell)
    }
}
*/

/// Convert patches to a simple Shell with disconnected faces
impl<'a> TryFrom<PatchTableWithControlPointsRef<'a>> for Shell {
    type Error = MonstertruckError;

    fn try_from(
        patches: PatchTableWithControlPointsRef<'a>,
    ) -> std::result::Result<Self, Self::Error> {
        // Create faces directly from patches to keep access to control points
        let mut faces = Vec::new();
        let mut patch_index = 0;

        for array_idx in 0..patches.patch_table.patch_array_count() {
            if let Some(desc) = patches.patch_table.patch_array_descriptor(array_idx) {
                let patch_type = desc.patch_type();
                let num_patches = patches.patch_table.patch_array_patch_count(array_idx);

                // Handle Regular, GregoryBasis, and GregoryTriangle patches.
                if matches!(
                    patch_type,
                    PatchType::Regular | PatchType::GregoryBasis | PatchType::GregoryTriangle
                ) {
                    for _local_idx in 0..num_patches {
                        let patch =
                            PatchRef::new(patches.patch_table, patch_index, patches.control_points);

                        // Get control matrix before try_into() consumes the patch.
                        #[cfg(feature = "monstertruck_export_boundary")]
                        let control_matrix = match patch.control_points() {
                            Ok(cp) => cp,
                            Err(_e) => {
                                patch_index += 1;
                                continue;
                            }
                        };

                        // Convert to monstertruck surface.
                        let surface: BsplineSurface<Point3> = match patch.try_into() {
                            Ok(s) => s,
                            Err(_e) => {
                                patch_index += 1;
                                continue;
                            }
                        };

                        #[cfg(feature = "monstertruck_export_boundary")]
                        {
                            let face = create_face_with_boundary(&control_matrix, surface);
                            faces.push(face);
                        }

                        #[cfg(not(feature = "monstertruck_export_boundary"))]
                        {
                            // Create face without explicit boundary - let monstertruck determine
                            // it.
                            let face = Face::new(vec![], Surface::BsplineSurface(surface));
                            faces.push(face);
                        }

                        patch_index += 1;
                    }
                } else {
                    patch_index += patches.patch_table.patch_array_patch_count(array_idx);
                }
            }
        }
        Ok(Shell::from(faces))
    }
}

/// Convert patches to a vector of individual Shells (one face per shell)
impl<'a> TryFrom<PatchTableWithControlPointsRef<'a>> for Vec<Shell> {
    type Error = MonstertruckError;

    fn try_from(
        patches: PatchTableWithControlPointsRef<'a>,
    ) -> std::result::Result<Self, Self::Error> {
        // Create one shell per surface for disconnected export
        let mut shells = Vec::new();
        let mut patch_index = 0;

        for array_idx in 0..patches.patch_table.patch_array_count() {
            if let Some(desc) = patches.patch_table.patch_array_descriptor(array_idx) {
                let patch_type = desc.patch_type();
                // Handle Regular, GregoryBasis, and GregoryTriangle patches
                if matches!(
                    patch_type,
                    PatchType::Regular | PatchType::GregoryBasis | PatchType::GregoryTriangle
                ) {
                    for _ in 0..patches.patch_table.patch_array_patch_count(array_idx) {
                        let patch =
                            PatchRef::new(patches.patch_table, patch_index, patches.control_points);

                        // Get control matrix before try_into() consumes the patch.
                        #[cfg(feature = "monstertruck_export_boundary")]
                        let control_matrix = patch.control_points()?;

                        // Convert to monstertruck surface
                        let surface: BsplineSurface<Point3> = patch.try_into()?;

                        #[cfg(feature = "monstertruck_export_boundary")]
                        {
                            let face = create_face_with_boundary(&control_matrix, surface);
                            shells.push(Shell::from(vec![face]));
                        }

                        #[cfg(not(feature = "monstertruck_export_boundary"))]
                        {
                            let face = Face::new(vec![], Surface::BsplineSurface(surface));
                            shells.push(Shell::from(vec![face]));
                        }

                        patch_index += 1;
                    }
                } else {
                    patch_index += patches.patch_table.patch_array_patch_count(array_idx);
                }
            }
        }

        Ok(shells)
    }
}

/// Convert patch evaluation result to Point3
impl From<PatchEvalResult> for Point3 {
    fn from(result: PatchEvalResult) -> Self {
        Point3::new(
            result.point[0] as f64,
            result.point[1] as f64,
            result.point[2] as f64,
        )
    }
}

/// Helper function to convert array to Vector3
pub fn array_to_vector3(v: &[f32; 3]) -> Vector3 {
    Vector3::new(v[0] as f64, v[1] as f64, v[2] as f64)
}

fn base_quad_shell_from_inf_sharp_edges(
    refiner: &crate::far::TopologyRefiner,
    control_points: &[[f32; 3]],
) -> Option<Shell> {
    const INF_SHARP_THRESHOLD: f32 = 10.0;
    let base = refiner.level(0)?;

    if base
        .face_vertices_iter()
        .any(|face_vertices| face_vertices.len() != 4)
    {
        return None;
    }

    if (0..base.edge_count()).any(|edge_index| {
        base.edge_sharpness(crate::Index::from(edge_index as u32)) < INF_SHARP_THRESHOLD
    }) {
        return None;
    }

    if control_points.len() < base.vertex_count() {
        return None;
    }

    let knot_vector = KnotVector::from(vec![-3.0, -2.0, -1.0, 0.0, 1.0, 2.0, 3.0, 4.0]);
    let params = [0.0, 1.0 / 3.0, 2.0 / 3.0, 1.0];
    let base_points: Vec<Point3> = control_points[..base.vertex_count()]
        .iter()
        .map(|point| Point3::new(point[0] as f64, point[1] as f64, point[2] as f64))
        .collect();

    let faces: Vec<_> = base
        .face_vertices_iter()
        .map(|face_vertices| {
            let corners: Vec<Point3> = face_vertices
                .iter()
                .map(|&index| base_points[usize::from(index)])
                .collect();
            let p00 = corners[0];
            let p10 = corners[1];
            let p11 = corners[2];
            let p01 = corners[3];
            let control_matrix: Vec<Vec<_>> = params
                .iter()
                .map(|&v| {
                    params
                        .iter()
                        .map(|&u| {
                            let u0 = 1.0 - u;
                            let v0 = 1.0 - v;
                            Point3::from_vec(
                                p00.to_vec() * (u0 * v0)
                                    + p10.to_vec() * (u * v0)
                                    + p11.to_vec() * (u * v)
                                    + p01.to_vec() * (u0 * v),
                            )
                        })
                        .collect()
                })
                .collect();

            Face::new(
                vec![],
                Surface::BsplineSurface(BsplineSurface::new(
                    (knot_vector.clone(), knot_vector.clone()),
                    control_matrix,
                )),
            )
        })
        .collect();

    Some(Shell::from(faces))
}

/// Build B-spline surfaces for regular faces using BFR (avoids
/// over-refinement). Irregular faces are skipped.
pub fn bfr_regular_surfaces(
    refiner: &crate::far::TopologyRefiner,
    control_points: &[[f32; 3]],
    approx_smooth: i32,
    approx_sharp: i32,
) -> Result<Vec<BsplineSurface<Point3>>> {
    let factory = BfrSurfaceFactory::new(refiner, approx_smooth, approx_sharp).map_err(|e| {
        MonstertruckError::BfrConversionFailed(format!("factory creation failed: {:?}", e))
    })?;

    factory
        .build_regular_surfaces(refiner, control_points)
        .map_err(|e| MonstertruckError::BfrConversionFailed(format!("{:?}", e)))
}

/// Create a triangular patch as a degenerate quad B-spline surface
/// This is used to fill gaps near extraordinary vertices
pub fn create_triangular_patch(
    p0: Point3,
    p1: Point3,
    p2: Point3,
    center: Point3,
) -> BsplineSurface<Point3> {
    // AIDEV-NOTE: Triangular patch creation for extraordinary vertices
    // When OpenSubdiv doesn't generate Gregory patches, we need to create
    // triangular patches to fill the gaps. We do this by creating a degenerate
    // quad patch where one edge collapses to a point.

    // Create control points for a degenerate quad patch
    // The patch will have one corner at the extraordinary vertex (center)
    // and form a triangle with p0, p1, p2

    // Compute intermediate control points using cubic interpolation
    let c01 = Point3::from_vec((p0.to_vec() * 2.0 + p1.to_vec()) / 3.0);
    let c10 = Point3::from_vec((p1.to_vec() * 2.0 + p0.to_vec()) / 3.0);
    let c02 = Point3::from_vec((p0.to_vec() * 2.0 + p2.to_vec()) / 3.0);
    let _c20 = Point3::from_vec((p2.to_vec() * 2.0 + p0.to_vec()) / 3.0);
    let c12 = Point3::from_vec((p1.to_vec() * 2.0 + p2.to_vec()) / 3.0);
    let _c21 = Point3::from_vec((p2.to_vec() * 2.0 + p1.to_vec()) / 3.0);

    // Center control points
    let cc = Point3::from_vec((p0.to_vec() + p1.to_vec() + p2.to_vec() + center.to_vec()) / 4.0);
    let cc0 = Point3::from_vec((center.to_vec() * 2.0 + p0.to_vec()) / 3.0);
    let cc1 = Point3::from_vec((center.to_vec() * 2.0 + p1.to_vec()) / 3.0);
    let cc2 = Point3::from_vec((center.to_vec() * 2.0 + p2.to_vec()) / 3.0);

    // Build 4x4 control point matrix
    // Row 0: degenerate to center point
    let control_matrix = vec![
        vec![center, center, center, center],
        vec![cc0, cc, cc1, cc2],
        vec![
            c02,
            Point3::from_vec((c02.to_vec() + c10.to_vec()) / 2.0),
            c10,
            c12,
        ],
        vec![p0, c01, p1, p2],
    ];

    // Use the same knot vectors as regular patches
    let u_knots = KnotVector::from(vec![-3.0, -2.0, -1.0, 0.0, 1.0, 2.0, 3.0, 4.0]);
    let v_knots = KnotVector::from(vec![-3.0, -2.0, -1.0, 0.0, 1.0, 2.0, 3.0, 4.0]);

    BsplineSurface::new((u_knots, v_knots), control_matrix)
}

/// Extension trait for PatchTable to provide conversion methods
pub trait PatchTableExt {
    /// Create a wrapper for conversion to monstertruck surfaces
    fn with_control_points<'a>(
        &'a self,
        control_points: &'a [[f32; 3]],
    ) -> PatchTableWithControlPointsRef<'a>;

    /// Get a specific patch for conversion
    fn patch<'a>(&'a self, index: usize, control_points: &'a [[f32; 3]]) -> PatchRef<'a>;

    /// Convert patches to a monstertruck shell with the given control points
    fn to_monstertruck_shell(&self, control_points: &[[f32; 3]]) -> Result<Shell>;

    /// Convert patches to monstertruck surfaces with the given control points
    fn to_monstertruck_surfaces(
        &self,
        control_points: &[[f32; 3]],
    ) -> Result<Vec<BsplineSurface<Point3>>>;

    /// Convert patches to monstertruck surfaces with configurable Gregory
    /// accuracy.
    ///
    /// This method allows specifying how Gregory patches (at extraordinary
    /// vertices) should be converted:
    /// - `BSplineEndCaps`: Standard 4×4 sampling (faster, slight approximation)
    /// - `HighPrecision`: 8×8 sampling for better accuracy at extraordinary
    ///   vertices
    fn to_monstertruck_surfaces_with_options(
        &self,
        control_points: &[[f32; 3]],
        gregory_accuracy: GregoryAccuracy,
    ) -> Result<Vec<BsplineSurface<Point3>>>;

    /// Prefer BFR for regular faces and fall back to PatchTable for non-regular
    /// patches.
    fn to_monstertruck_surfaces_bfr_mixed(
        &self,
        refiner: &TopologyRefiner,
        control_points: &[[f32; 3]],
        approx_smooth: i32,
        approx_sharp: i32,
    ) -> Result<Vec<BsplineSurface<Point3>>>;

    /// Build the coarsest direct `monstertruck` shell practical for
    /// interactive consumers.
    ///
    /// This uses BFR for regular regions and patch-table fallback for
    /// extraordinary regions with zero additional approximation refinement.
    fn to_monstertruck_shell_coarse(
        &self,
        refiner: &TopologyRefiner,
        control_points: &[[f32; 3]],
    ) -> Result<Shell>;

    /// Build a shell using BFR for regular faces and PatchTable for irregular
    /// patches.
    fn to_monstertruck_shell_bfr_mixed(
        &self,
        refiner: &TopologyRefiner,
        control_points: &[[f32; 3]],
        approx_smooth: i32,
        approx_sharp: i32,
    ) -> Result<Shell>;

    /// Build a stitched B-rep shell with shared vertices/edges between patches.
    fn to_monstertruck_shell_stitched(&self, control_points: &[[f32; 3]]) -> Result<Shell>;

    /// Build a `monstertruck` shell with explicit conversion options.
    ///
    /// This is the direct conversion entry point for consumers such as Akatela.
    /// STEP export should use this path and serialize the returned shell.
    fn to_monstertruck_shell_with_options(
        &self,
        control_points: &[[f32; 3]],
        options: MonstertruckShellOptions,
    ) -> Result<Shell>;

    /// Convert patches to individual shells (one per patch) for disconnected
    /// export
    fn to_monstertruck_shells(&self, control_points: &[[f32; 3]]) -> Result<Vec<Shell>>;

    /// Convert patches to a shell with gap filling for extraordinary vertices
    fn to_monstertruck_shell_with_gap_filling(&self, control_points: &[[f32; 3]]) -> Result<Shell>;

    /// Export patches to a `monstertruck` shell for STEP output.
    ///
    /// This is a compatibility wrapper over
    /// [`PatchTableExt::to_monstertruck_shell_with_options`].
    ///
    /// # Example
    ///
    /// ```ignore
    /// use opensubdiv_petite::monstertruck::{PatchTableExt, StepExportOptions, GregoryAccuracy};
    ///
    /// // Default export (BSpline end caps, no stitching, superpatch merging)
    /// let shell = patch_table.to_step_shell(&vertices, Default::default())?;
    ///
    /// // High precision Gregory fitting
    /// let shell = patch_table.to_step_shell(&vertices, StepExportOptions {
    ///     gregory_accuracy: GregoryAccuracy::HighPrecision,
    ///     ..Default::default()
    /// })?;
    ///
    /// // Full B-rep with stitched edges
    /// let shell = patch_table.to_step_shell(&vertices, StepExportOptions {
    ///     stitch_edges: true,
    ///     ..Default::default()
    /// })?;
    /// ```
    fn to_step_shell(
        &self,
        control_points: &[[f32; 3]],
        options: StepExportOptions,
    ) -> Result<Shell>;

    /// Internal fallback for `to_step_shell`.
    #[doc(hidden)]
    fn to_step_shell_fallback(
        &self,
        control_points: &[[f32; 3]],
        options: &StepExportOptions,
    ) -> Result<Shell>;
}

impl PatchTableExt for PatchTable {
    fn with_control_points<'a>(
        &'a self,
        control_points: &'a [[f32; 3]],
    ) -> PatchTableWithControlPointsRef<'a> {
        PatchTableWithControlPointsRef {
            patch_table: self,
            control_points,
        }
    }

    fn patch<'a>(&'a self, index: usize, control_points: &'a [[f32; 3]]) -> PatchRef<'a> {
        PatchRef::new(self, index, control_points)
    }

    fn to_monstertruck_shell(&self, control_points: &[[f32; 3]]) -> Result<Shell> {
        let wrapper = self.with_control_points(control_points);
        Shell::try_from(wrapper)
    }

    fn to_monstertruck_surfaces(
        &self,
        control_points: &[[f32; 3]],
    ) -> Result<Vec<BsplineSurface<Point3>>> {
        let wrapper = self.with_control_points(control_points);
        Vec::<BsplineSurface<Point3>>::try_from(wrapper)
    }

    fn to_monstertruck_surfaces_with_options(
        &self,
        control_points: &[[f32; 3]],
        gregory_accuracy: GregoryAccuracy,
    ) -> Result<Vec<BsplineSurface<Point3>>> {
        let convert_patch = |patch_index: usize| {
            let patch_ref = self.patch(patch_index, control_points);
            if patch_ref.is_gregory() && gregory_accuracy == GregoryAccuracy::HighPrecision {
                // Use high-precision 8×8 sampling for Gregory patches.
                patch_ref.to_bspline_high_precision()
            } else {
                // Use standard conversion for regular patches or when using
                // BSplineEndCaps accuracy (the patch table should already have
                // B-spline end caps if configured that way).
                BsplineSurface::try_from(patch_ref)
            }
        };

        #[cfg(feature = "rayon")]
        {
            (0..self.patch_count())
                .into_par_iter()
                .map(convert_patch)
                .collect()
        }
        #[cfg(not(feature = "rayon"))]
        {
            (0..self.patch_count()).map(convert_patch).collect()
        }
    }

    /// Prefer BFR for regular faces and fall back to PatchTable for non-regular
    /// patches. BFR approximation levels control how far sharp/smooth
    /// features refine; use 0/0 to keep base quads coarse.
    fn to_monstertruck_surfaces_bfr_mixed(
        &self,
        refiner: &TopologyRefiner,
        control_points: &[[f32; 3]],
        approx_smooth: i32,
        approx_sharp: i32,
    ) -> Result<Vec<BsplineSurface<Point3>>> {
        let mut surfaces =
            match bfr_regular_surfaces(refiner, control_points, approx_smooth, approx_sharp) {
                Ok(s) => s,
                Err(e) => {
                    eprintln!("BFR regular surface build failed: {:?}", e);
                    Vec::new()
                }
            };

        let mut fallback = match patch_table_surfaces_non_regular(self, control_points) {
            Ok(f) => f,
            Err(e) => {
                eprintln!("Non-regular PatchTable conversion failed: {:?}", e);
                Vec::new()
            }
        };
        surfaces.append(&mut fallback);

        if surfaces.is_empty() {
            // No BFR or non-regular outputs; fall back to full PatchTable
            // conversion so we still return curved 4x4 patches.
            self.to_monstertruck_surfaces(control_points)
        } else {
            Ok(surfaces)
        }
    }

    fn to_monstertruck_shell_coarse(
        &self,
        refiner: &TopologyRefiner,
        control_points: &[[f32; 3]],
    ) -> Result<Shell> {
        if let Some(shell) = base_quad_shell_from_inf_sharp_edges(refiner, control_points) {
            Ok(shell)
        } else {
            self.to_monstertruck_shell_bfr_mixed(refiner, control_points, 0, 0)
        }
    }

    fn to_monstertruck_shells(&self, control_points: &[[f32; 3]]) -> Result<Vec<Shell>> {
        let wrapper = self.with_control_points(control_points);
        Vec::<Shell>::try_from(wrapper)
    }

    fn to_monstertruck_shell_with_gap_filling(&self, control_points: &[[f32; 3]]) -> Result<Shell> {
        // AIDEV-NOTE: Gap filling for extraordinary vertices
        // This method detects gaps in the patch layout and fills them with
        // triangular patches. This is a workaround for when OpenSubdiv
        // doesn't generate Gregory patches at extraordinary vertices.

        // First, convert regular patches.
        let wrapper = self.with_control_points(control_points);
        let shell = Shell::try_from(wrapper)?;

        Ok(shell)
    }

    fn to_monstertruck_shell_bfr_mixed(
        &self,
        refiner: &TopologyRefiner,
        control_points: &[[f32; 3]],
        approx_smooth: i32,
        approx_sharp: i32,
    ) -> Result<Shell> {
        let surfaces = self.to_monstertruck_surfaces_bfr_mixed(
            refiner,
            control_points,
            approx_smooth,
            approx_sharp,
        )?;
        let faces: Vec<Face> = surfaces
            .into_iter()
            .map(|s| Face::new(vec![], Surface::BsplineSurface(s)))
            .collect();
        Ok(Shell::from(faces))
    }

    fn to_monstertruck_shell_stitched(&self, control_points: &[[f32; 3]]) -> Result<Shell> {
        // AIDEV-NOTE: Tolerant weld mode stitches edges using curve geometry samples
        // instead of control cages so slight mismatches still share topology.
        const STITCH_TOL: f64 = 1.0e-6;
        const STITCH_SAMPLES: usize = 7;

        #[derive(Clone)]
        struct EdgeSamples {
            start: Point3,
            end: Point3,
            samples: Vec<Point3>,
        }

        #[derive(Clone)]
        struct EdgeEntry {
            samples: EdgeSamples,
            edge: monstertruck::modeling::Edge,
        }

        enum OrientationMatch {
            Aligned,
            Reversed,
        }

        fn points_within(a: &Point3, b: &Point3, tol_sq: f64) -> bool {
            a.distance2(*b) <= tol_sq
        }

        fn sample_curve_points(curve: &BsplineCurve<Point3>, sample_count: usize) -> Vec<Point3> {
            let start = curve.knot(0);
            let end = curve.knot(curve.knot_vector().len() - 1);

            if sample_count <= 1 {
                return vec![curve.evaluate((start + end) * 0.5)];
            }

            let step = (end - start) / (sample_count as f64 - 1.0);
            (0..sample_count)
                .map(|i| curve.evaluate(start + step * i as f64))
                .collect()
        }

        fn samples_match(
            reference: &EdgeSamples,
            candidate: &EdgeSamples,
            tol_sq: f64,
            reversed: bool,
        ) -> bool {
            let (start_ref, end_ref) = if reversed {
                (&reference.end, &reference.start)
            } else {
                (&reference.start, &reference.end)
            };

            if !points_within(start_ref, &candidate.start, tol_sq)
                || !points_within(end_ref, &candidate.end, tol_sq)
            {
                return false;
            }

            let mut iter: Box<dyn Iterator<Item = (&Point3, &Point3)>> = if reversed {
                Box::new(reference.samples.iter().zip(candidate.samples.iter().rev()))
            } else {
                Box::new(reference.samples.iter().zip(candidate.samples.iter()))
            };

            iter.all(|(a, b)| points_within(a, b, tol_sq))
        }

        fn find_matching_edge(
            entries: &[EdgeEntry],
            candidate: &EdgeSamples,
            tol_sq: f64,
        ) -> Option<(monstertruck::modeling::Edge, OrientationMatch)> {
            entries.iter().find_map(|entry| {
                if samples_match(&entry.samples, candidate, tol_sq, false) {
                    Some((entry.edge.clone(), OrientationMatch::Aligned))
                } else if samples_match(&entry.samples, candidate, tol_sq, true) {
                    Some((entry.edge.clone(), OrientationMatch::Reversed))
                } else {
                    None
                }
            })
        }

        let mut faces = Vec::new();
        let mut patch_index = 0usize;
        let mut vertex_pool: Vec<(Point3, monstertruck::modeling::Vertex)> = Vec::new();
        let mut edge_entries: Vec<EdgeEntry> = Vec::new();
        let tol_sq = STITCH_TOL * STITCH_TOL;
        let mut invalid_count = 0usize;

        let wrapper = self.with_control_points(control_points);

        for array_idx in 0..wrapper.patch_table.patch_array_count() {
            if let Some(desc) = wrapper.patch_table.patch_array_descriptor(array_idx) {
                let patch_type = desc.patch_type();
                if !matches!(
                    patch_type,
                    PatchType::Regular | PatchType::GregoryBasis | PatchType::GregoryTriangle
                ) {
                    patch_index += wrapper.patch_table.patch_array_patch_count(array_idx);
                    continue;
                }

                let num_patches = wrapper.patch_table.patch_array_patch_count(array_idx);
                for _ in 0..num_patches {
                    let patch =
                        PatchRef::new(wrapper.patch_table, patch_index, wrapper.control_points);
                    let control_matrix = match patch.control_points() {
                        Ok(cp) => cp,
                        Err(_) => {
                            patch_index += 1;
                            continue;
                        }
                    };

                    let surface: BsplineSurface<Point3> = match patch.try_into() {
                        Ok(s) => s,
                        Err(_) => {
                            patch_index += 1;
                            continue;
                        }
                    };

                    let find_vertex = |pool: &mut Vec<(Point3, monstertruck::modeling::Vertex)>,
                                       p: Point3| {
                        if let Some((_, v)) =
                            pool.iter().find(|(pos, _)| points_within(pos, &p, tol_sq))
                        {
                            v.clone()
                        } else {
                            let v = monstertruck::modeling::Vertex::new(p);
                            pool.push((p, v.clone()));
                            v
                        }
                    };

                    let make_edge_curve = |cps: Vec<Point3>| -> BsplineCurve<Point3> {
                        let edge_knots =
                            KnotVector::from(vec![-3.0, -2.0, -1.0, 0.0, 1.0, 2.0, 3.0, 4.0]);
                        BsplineCurve::new(edge_knots, cps)
                    };

                    let bottom_curve = make_edge_curve(vec![
                        control_matrix[0][0],
                        control_matrix[0][1],
                        control_matrix[0][2],
                        control_matrix[0][3],
                    ]);
                    let right_curve = make_edge_curve(vec![
                        control_matrix[0][3],
                        control_matrix[1][3],
                        control_matrix[2][3],
                        control_matrix[3][3],
                    ]);
                    let top_curve = make_edge_curve(vec![
                        control_matrix[3][3],
                        control_matrix[3][2],
                        control_matrix[3][1],
                        control_matrix[3][0],
                    ]);
                    let left_curve = make_edge_curve(vec![
                        control_matrix[3][0],
                        control_matrix[2][0],
                        control_matrix[1][0],
                        control_matrix[0][0],
                    ]);

                    let bottom_samples = sample_curve_points(&bottom_curve, STITCH_SAMPLES);
                    let right_samples = sample_curve_points(&right_curve, STITCH_SAMPLES);
                    let top_samples = sample_curve_points(&top_curve, STITCH_SAMPLES);
                    let left_samples = sample_curve_points(&left_curve, STITCH_SAMPLES);

                    let distinct_vertex =
                        |pool: &mut Vec<(Point3, monstertruck::modeling::Vertex)>,
                         p: Point3,
                         avoid: &monstertruck::modeling::Vertex| {
                            if points_within(&avoid.point(), &p, tol_sq) {
                                let v = monstertruck::modeling::Vertex::new(p);
                                pool.push((p, v.clone()));
                                v
                            } else {
                                find_vertex(pool, p)
                            }
                        };

                    let v00 = find_vertex(
                        &mut vertex_pool,
                        *bottom_samples.first().unwrap_or(&control_matrix[0][0]),
                    );
                    let v10 = distinct_vertex(
                        &mut vertex_pool,
                        *bottom_samples.last().unwrap_or(&control_matrix[0][3]),
                        &v00,
                    );
                    let v11 = distinct_vertex(
                        &mut vertex_pool,
                        *right_samples.last().unwrap_or(&control_matrix[3][3]),
                        &v10,
                    );
                    let v01 = distinct_vertex(
                        &mut vertex_pool,
                        *left_samples.first().unwrap_or(&control_matrix[3][0]),
                        &v11,
                    );

                    let edges = [
                        (
                            v00.clone(),
                            v10.clone(),
                            bottom_curve,
                            EdgeSamples {
                                start: *bottom_samples.first().unwrap_or(&v00.point()),
                                end: *bottom_samples.last().unwrap_or(&v10.point()),
                                samples: bottom_samples,
                            },
                        ),
                        (
                            v10.clone(),
                            v11.clone(),
                            right_curve,
                            EdgeSamples {
                                start: *right_samples.first().unwrap_or(&v10.point()),
                                end: *right_samples.last().unwrap_or(&v11.point()),
                                samples: right_samples,
                            },
                        ),
                        (
                            v11.clone(),
                            v01.clone(),
                            top_curve,
                            EdgeSamples {
                                start: *top_samples.first().unwrap_or(&v11.point()),
                                end: *top_samples.last().unwrap_or(&v01.point()),
                                samples: top_samples,
                            },
                        ),
                        (
                            v01.clone(),
                            v00.clone(),
                            left_curve,
                            EdgeSamples {
                                start: *left_samples.first().unwrap_or(&v01.point()),
                                end: *left_samples.last().unwrap_or(&v00.point()),
                                samples: left_samples,
                            },
                        ),
                    ];

                    let edge_count = edges.len();
                    let edges = edges.into_iter().collect::<Vec<_>>();
                    let mut wire_edges = Vec::with_capacity(4);
                    let mut first_front: Option<monstertruck::modeling::Vertex> = None;
                    let mut prev_end: Option<monstertruck::modeling::Vertex> = None;

                    let linear_curve = |v0: &monstertruck::modeling::Vertex,
                                        v1: &monstertruck::modeling::Vertex|
                     -> monstertruck::modeling::Curve {
                        monstertruck::modeling::Curve::BsplineCurve(BsplineCurve::new(
                            KnotVector::bezier_knot(1),
                            vec![v0.point(), v1.point()],
                        ))
                    };

                    let build_edge = |v0: &monstertruck::modeling::Vertex,
                                      v1: &monstertruck::modeling::Vertex,
                                      curve: monstertruck::modeling::Curve|
                     -> monstertruck::modeling::Edge {
                        if v0.id() == v1.id() {
                            let dup = monstertruck::modeling::Vertex::new(v1.point());
                            monstertruck::modeling::Edge::new(v0, &dup, curve)
                        } else {
                            monstertruck::modeling::Edge::new(v0, v1, curve)
                        }
                    };

                    for (edge_idx, (start_vertex, end_vertex, _curve, samples)) in
                        edges.into_iter().enumerate()
                    {
                        let candidate = samples.clone();
                        let start_vertex = if let Some(prev) = &prev_end {
                            if points_within(&prev.point(), &candidate.start, tol_sq) {
                                prev.clone()
                            } else {
                                start_vertex
                            }
                        } else {
                            start_vertex
                        };

                        let end_vertex =
                            if points_within(&start_vertex.point(), &candidate.end, tol_sq) {
                                monstertruck::modeling::Vertex::new(candidate.end)
                            } else {
                                end_vertex
                            };

                        let mut edge = if let Some((existing_edge, orientation)) =
                            find_matching_edge(&edge_entries, &candidate, tol_sq)
                        {
                            match orientation {
                                OrientationMatch::Aligned => existing_edge,
                                OrientationMatch::Reversed => existing_edge.inverse(),
                            }
                        } else {
                            build_edge(
                                &start_vertex,
                                &end_vertex,
                                linear_curve(&start_vertex, &end_vertex),
                            )
                        };

                        if let Some(prev) = &prev_end {
                            if !points_within(&prev.point(), &edge.front().point(), tol_sq)
                                && points_within(&prev.point(), &edge.back().point(), tol_sq)
                            {
                                edge = edge.inverse();
                            }
                        }

                        if edge_idx == edge_count - 1 {
                            if let Some(first) = &first_front {
                                if points_within(&edge.back().point(), &first.point(), tol_sq)
                                    && edge.back().id() != first.id()
                                {
                                    let curve = edge.curve().clone();
                                    edge = build_edge(&edge.front().clone(), first, curve);
                                } else if points_within(
                                    &edge.front().point(),
                                    &first.point(),
                                    tol_sq,
                                ) && edge.front().id() != first.id()
                                {
                                    let curve = edge.curve().clone();
                                    edge = build_edge(first, &edge.back().clone(), curve);
                                }
                            }
                        }

                        if let Some(prev) = &prev_end {
                            if edge.front().id() != prev.id()
                                && points_within(&edge.front().point(), &prev.point(), tol_sq)
                            {
                                let curve = edge.curve().clone();
                                edge = build_edge(prev, &edge.back().clone(), curve);
                            }
                        }

                        // If we created a brand new edge, track it for reuse.
                        if !edge_entries
                            .iter()
                            .any(|entry| entry.edge.id() == edge.id())
                        {
                            edge_entries.push(EdgeEntry {
                                samples: candidate,
                                edge: edge.clone(),
                            });
                        }

                        if first_front.is_none() {
                            first_front = Some(edge.front().clone());
                        }
                        prev_end = Some(edge.back().clone());
                        wire_edges.push(edge);
                    }

                    let surface_clone = surface.clone();
                    let enforce_continuity = |edges: Vec<monstertruck::modeling::Edge>| {
                        if edges.is_empty() {
                            return edges;
                        }

                        let mut fixed = Vec::with_capacity(edges.len());
                        let first_front = edges[0].front().clone();
                        let mut prev_back = edges[0].back().clone();
                        fixed.push(edges[0].clone());

                        for edge in edges.into_iter().skip(1) {
                            let mut edge = edge;
                            if !points_within(&edge.front().point(), &prev_back.point(), tol_sq)
                                && points_within(&edge.back().point(), &prev_back.point(), tol_sq)
                            {
                                edge = edge.inverse();
                            }

                            if edge.front().id() != prev_back.id() {
                                let curve = edge.curve().clone();
                                edge = build_edge(&prev_back, &edge.back().clone(), curve);
                            }

                            prev_back = edge.back().clone();
                            fixed.push(edge);
                        }

                        let last_idx = fixed.len().saturating_sub(1);
                        if let Some(last) = fixed.get(last_idx) {
                            if !points_within(&last.back().point(), &first_front.point(), tol_sq)
                                && points_within(
                                    &last.front().point(),
                                    &first_front.point(),
                                    tol_sq,
                                )
                            {
                                let mut last_edge = last.clone().inverse();
                                if last_edge.front().id() != first_front.id() {
                                    let curve = last_edge.curve().clone();
                                    last_edge =
                                        build_edge(&first_front, &last_edge.back().clone(), curve);
                                }
                                fixed[last_idx] = last_edge;
                            } else if last.back().id() != first_front.id() {
                                let curve = last.curve().clone();
                                fixed[last_idx] =
                                    build_edge(&last.front().clone(), &first_front, curve);
                            }
                        }

                        fixed
                    };

                    let validate_wire =
                        |edges: &[monstertruck::modeling::Edge],
                         tol_sq: f64|
                         -> std::result::Result<(), &'static str> {
                            if edges.len() != 4 {
                                return Err("edge count mismatch");
                            }

                            let mut edge_ids = std::collections::HashSet::new();
                            if !edges.iter().all(|e| edge_ids.insert(e.id())) {
                                return Err("duplicate edges");
                            }

                            for i in 0..edges.len() {
                                let a_back = edges[i].back().point();
                                let b_front = edges[(i + 1) % edges.len()].front().point();
                                if !points_within(&a_back, &b_front, tol_sq) {
                                    return Err("edge adjacency gap");
                                }
                            }

                            Ok(())
                        };

                    let mut wire_edges = enforce_continuity(wire_edges);
                    let closed = wire_edges.iter().enumerate().all(|(i, e)| {
                        e.back().id() == wire_edges[(i + 1) % wire_edges.len()].front().id()
                    });

                    if !closed {
                        eprintln!("Stitched wire not closed for patch {}", patch_index);
                        if wire_edges.len() == 4 {
                            let v0 = wire_edges[0].front().clone();
                            let v1 = wire_edges[0].back().clone();
                            let v2 = wire_edges[1].back().clone();
                            let v3 = wire_edges[2].back().clone();

                            let mut corners = [v0.clone(), v1.clone(), v2.clone(), v3.clone()];
                            for i in 0..corners.len() {
                                let next = (i + 1) % corners.len();
                                if corners[i].id() == corners[next].id() {
                                    corners[next] =
                                        monstertruck::modeling::Vertex::new(corners[next].point());
                                }
                            }

                            wire_edges = vec![
                                build_edge(
                                    &corners[0],
                                    &corners[1],
                                    linear_curve(&corners[0], &corners[1]),
                                ),
                                build_edge(
                                    &corners[1],
                                    &corners[2],
                                    linear_curve(&corners[1], &corners[2]),
                                ),
                                build_edge(
                                    &corners[2],
                                    &corners[3],
                                    linear_curve(&corners[2], &corners[3]),
                                ),
                                build_edge(
                                    &corners[3],
                                    &corners[0],
                                    linear_curve(&corners[3], &corners[0]),
                                ),
                            ];

                            let closed = wire_edges.iter().enumerate().all(|(i, e)| {
                                e.back().id() == wire_edges[(i + 1) % wire_edges.len()].front().id()
                            });
                            if !closed {
                                invalid_count += 1;
                                patch_index += 1;
                                continue;
                            }
                        } else {
                            invalid_count += 1;
                            patch_index += 1;
                            continue;
                        }
                    }

                    match validate_wire(&wire_edges, tol_sq) {
                        Ok(()) => match panic::catch_unwind(panic::AssertUnwindSafe(|| {
                            Face::new(
                                vec![monstertruck::modeling::Wire::from(wire_edges)],
                                Surface::BsplineSurface(surface_clone),
                            )
                        })) {
                            Ok(face) => faces.push(face),
                            Err(_) => {
                                eprintln!(
                                    "Stitched wire failed to build face for patch {}",
                                    patch_index
                                );
                                invalid_count += 1;
                            }
                        },
                        Err(reason) => {
                            eprintln!(
                                "Stitched wire invalid for patch {}: {}",
                                patch_index, reason
                            );
                            for (i, e) in wire_edges.iter().enumerate() {
                                let f = e.front().point();
                                let b = e.back().point();
                                eprintln!(
                                    "  edge {}: front id {:?} ({:.6},{:.6},{:.6}) -> back id {:?} ({:.6},{:.6},{:.6})",
                                    i,
                                    e.front().id(),
                                    f.x,
                                    f.y,
                                    f.z,
                                    e.back().id(),
                                    b.x,
                                    b.y,
                                    b.z
                                );
                            }
                            invalid_count += 1;
                        }
                    }
                    patch_index += 1;
                }
            }
        }

        if invalid_count > 0 {
            eprintln!(
                "Stitched export skipped welding on {} patches with invalid wires.",
                invalid_count
            );
        }

        if faces.is_empty() {
            Err(MonstertruckError::InvalidControlPoints)
        } else {
            Ok(Shell::from(faces))
        }
    }

    fn to_monstertruck_shell_with_options(
        &self,
        control_points: &[[f32; 3]],
        options: MonstertruckShellOptions,
    ) -> Result<Shell> {
        match options.flavor {
            MonstertruckShellFlavor::CoarseMerged => {
                let mut surfaces =
                    superpatch_surfaces(self, control_points, options.stitch_tolerance)
                        .unwrap_or_default();
                let mut fallback = patch_table_surfaces_non_regular_with_options(
                    self,
                    control_points,
                    options.gregory_accuracy,
                )
                .unwrap_or_default();
                surfaces.append(&mut fallback);

                let surfaces = if surfaces.is_empty() {
                    self.to_monstertruck_surfaces_with_options(
                        control_points,
                        options.gregory_accuracy,
                    )?
                } else {
                    surfaces
                };

                let faces: Vec<_> = surfaces
                    .into_iter()
                    .map(|surface| Face::new(vec![], Surface::BsplineSurface(surface)))
                    .collect();
                Ok(Shell::from(faces))
            }
            MonstertruckShellFlavor::DisconnectedPatches => {
                let surfaces = self.to_monstertruck_surfaces_with_options(
                    control_points,
                    options.gregory_accuracy,
                )?;
                let faces: Vec<_> = surfaces
                    .into_iter()
                    .map(|surface| Face::new(vec![], Surface::BsplineSurface(surface)))
                    .collect();
                Ok(Shell::from(faces))
            }
            MonstertruckShellFlavor::StitchedPatches => {
                self.to_monstertruck_shell_stitched(control_points)
            }
        }
    }

    fn to_step_shell(
        &self,
        control_points: &[[f32; 3]],
        options: StepExportOptions,
    ) -> Result<Shell> {
        self.to_monstertruck_shell_with_options(control_points, options.into())
    }

    fn to_step_shell_fallback(
        &self,
        control_points: &[[f32; 3]],
        options: &StepExportOptions,
    ) -> Result<Shell> {
        self.to_monstertruck_shell_with_options(control_points, options.clone().into())
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::far::{
        AdaptiveRefinementOptions, EndCapType, PatchTable, PatchTableOptions, PrimvarRefiner,
        TopologyDescriptor, TopologyRefiner, TopologyRefinerOptions,
    };

    fn creased_cube_patch_table() -> (TopologyRefiner, PatchTable, Vec<[f32; 3]>) {
        let vertex_positions = vec![
            [-0.5, -0.5, 0.5],
            [0.5, -0.5, 0.5],
            [-0.5, 0.5, 0.5],
            [0.5, 0.5, 0.5],
            [-0.5, 0.5, -0.5],
            [0.5, 0.5, -0.5],
            [-0.5, -0.5, -0.5],
            [0.5, -0.5, -0.5],
        ];

        let face_vertex_counts = vec![4, 4, 4, 4, 4, 4];
        let face_vertex_indices = vec![
            0, 1, 3, 2, //
            2, 3, 5, 4, //
            4, 5, 7, 6, //
            6, 7, 1, 0, //
            0, 2, 4, 6, //
            1, 7, 5, 3, //
        ];
        let crease_indices = vec![
            0, 1, //
            1, 3, //
            3, 2, //
            2, 0, //
            2, 4, //
            4, 6, //
            6, 0, //
            1, 7, //
            7, 5, //
            5, 3, //
            4, 5, //
            6, 7, //
        ];
        let crease_weights = vec![10.0; crease_indices.len() / 2];

        let mut descriptor = TopologyDescriptor::new(
            vertex_positions.len(),
            &face_vertex_counts,
            &face_vertex_indices,
        )
        .expect("Failed to create topology descriptor.");
        descriptor.creases(&crease_indices, &crease_weights);

        let mut refiner = TopologyRefiner::new(descriptor, TopologyRefinerOptions::default())
            .expect("Failed to create topology refiner.");
        refiner.refine_adaptive(
            AdaptiveRefinementOptions {
                isolation_level: 3,
                ..Default::default()
            },
            &[],
        );

        let patch_table = PatchTable::new(
            &refiner,
            Some(
                PatchTableOptions::new()
                    .end_cap_type(EndCapType::GregoryBasis)
                    .use_inf_sharp_patch(true),
            ),
        )
        .expect("Failed to create patch table.");

        let primvar_refiner =
            PrimvarRefiner::new(&refiner).expect("Failed to create primvar refiner.");
        let mut all_vertices = vertex_positions.clone();
        let mut level_start = 0usize;

        (1..refiner.refinement_levels()).for_each(|level| {
            let prev_level_count = refiner
                .level(level - 1)
                .map(|topology_level| topology_level.vertex_count())
                .unwrap_or(0);
            let src_data: Vec<f32> = all_vertices[level_start..level_start + prev_level_count]
                .iter()
                .flat_map(|vertex| vertex.iter().copied())
                .collect();

            if let Some(refined) = primvar_refiner.interpolate(level, 3, &src_data) {
                let level_vertices: Vec<[f32; 3]> = refined
                    .chunks_exact(3)
                    .map(|chunk| [chunk[0], chunk[1], chunk[2]])
                    .collect();
                all_vertices.extend_from_slice(&level_vertices);
            }

            level_start += prev_level_count;
        });

        if let Some(stencil_table) = patch_table.local_point_stencil_table() {
            let mut local_points = Vec::with_capacity(patch_table.local_point_count());

            (0..3).for_each(|dim| {
                let src_dim: Vec<f32> = all_vertices.iter().map(|vertex| vertex[dim]).collect();
                let dst_dim = stencil_table.update_values(&src_dim, None, None);

                dst_dim.iter().enumerate().for_each(|(index, &value)| {
                    if dim == 0 {
                        local_points.push([value, 0.0, 0.0]);
                    } else {
                        local_points[index][dim] = value;
                    }
                });
            });

            all_vertices.extend_from_slice(&local_points);
        }

        drop(primvar_refiner);
        (refiner, patch_table, all_vertices)
    }

    #[test]
    fn test_from_traits() {
        // This would require a proper patch table setup
        // Example usage:
        // let patch_table: PatchTable = ...;
        // let control_points: Vec<[f32; 3]> = ...;
        //
        // // Convert a single patch
        // let surface: BsplineSurface<Point3> = patch_table.patch(0,
        // &control_points).try_into()?;
        //
        // // Convert all patches to surfaces
        // let surfaces: Vec<BsplineSurface<Point3>> =
        // patch_table.with_control_points(&control_points).try_into()?;
        //
        // // Convert to shell
        // let shell: Shell =
        // patch_table.with_control_points(&control_points).try_into()?;
    }

    #[test]
    fn coarse_monstertruck_shell_keeps_cube_as_six_faces() -> Result<()> {
        let (refiner, patch_table, all_vertices) = creased_cube_patch_table();
        let shell = patch_table.to_monstertruck_shell_coarse(&refiner, &all_vertices)?;

        assert_eq!(shell.face_iter().count(), 6);

        Ok(())
    }
}