poincare-lib 0.7.0

use std::collections::HashMap;
use std::f64::consts::PI;

use glam::Vec3;
use viewport_lib::{AttributeData, BuiltinColourmap, MeshData};

use crate::{
    ArrowAnnotation, ColormapSource, ColourMode, CurveInterpolation, CurveInterpolationKind,
    Diagnostic, DiagnosticKind, DiagnosticLocation, PlotMetadata, PlotSpec, PlotStyle,
    PointAnnotation, SliceAxis, TableDataSet, default_slice_position, eval_curve_point,
    eval_with_vars, parse_curve_expr, parse_expr_with_vars, sample_curve_points,
};

/// The type of analysis to perform on a plot.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum AnalysisKind {
    PointCloudStatistics,
    DataQualityChecks,
    InterpolateCurve,
    FitCurve,
    DifferentiateCurve,
    AxisDerivativeCurve,
    IntegralCurve,
    ArcLengthCurve,
    CurvatureCurve,
    TangentField,
    NormalField,
    BinormalField,
    FrenetFrame,
    BishopFrame,
    DarbouxFrame,
    SurfaceAlignedFrame,
    ExtractPoints,
    ScalarSlice,
    VectorSlice,
    GradientField,
    DivergenceField,
    CurlField,
    SurfaceNormals,
    SurfaceCurvature,
    SurfaceArea,
    SurfaceMeshQuality,
    CurveSurfaceMeasurement,
    SurfaceIntersection,
}

/// The form of data returned by an analysis.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum AnalysisOutputKind {
    /// Produces new derived plots (e.g. a normal field or fitted curve).
    PlotSpec,
    /// Produces a flat key/value numeric summary.
    NumericReport,
    /// Produces a labelled data table.
    Table,
    /// Produces a mix of plots, reports, and tables.
    Composite,
}

/// What kind of data an analysis operates on.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum AnalysisTargetKind {
    /// Works from the symbolic plot definition (e.g. expression parsing).
    Definition,
    /// Works from sampled point data (e.g. CSV imports, point clouds).
    SampledData,
    /// Works from tessellated mesh geometry.
    Geometry,
    /// Requires two plots as input (e.g. intersection analysis).
    PlotPair,
}

/// Describes one analysis that is available for a given plot, as returned by [`available_analyses`].
#[derive(Clone, Debug, PartialEq)]
pub struct AnalysisCapability {
    pub kind: AnalysisKind,
    pub target_kind: AnalysisTargetKind,
    pub output_kind: AnalysisOutputKind,
    /// Named parameters this analysis accepts (passed via [`AnalysisRequest::parameters`]).
    pub parameters: Vec<&'static str>,
}

/// Specifies which plot(s) an analysis should run on.
#[derive(Clone, Debug, PartialEq)]
pub enum AnalysisTarget {
    Plot { index: usize, name: Option<String> },
    PlotPair { first: usize, second: usize },
}

/// Input to [`run_analysis`].
#[derive(Clone, Debug, PartialEq)]
pub struct AnalysisRequest {
    pub kind: AnalysisKind,
    pub target: AnalysisTarget,
    /// Key/value string parameters; valid keys depend on the analysis kind.
    pub parameters: Vec<(String, String)>,
}

/// Metadata attached to every [`AnalysisOutput`] describing what produced it.
#[derive(Clone, Debug, PartialEq)]
pub struct AnalysisProvenance {
    pub kind: AnalysisKind,
    pub source_plots: Vec<String>,
    pub parameters: Vec<(String, String)>,
    pub notes: Vec<String>,
}

/// A flat key/value numeric summary returned by some analyses.
#[derive(Clone, Debug, PartialEq)]
pub struct AnalysisReport {
    pub title: String,
    pub values: Vec<(String, String)>,
}

/// A labelled data table returned by some analyses.
#[derive(Clone, Debug, PartialEq)]
pub struct AnalysisTable {
    pub title: String,
    pub columns: Vec<String>,
    pub rows: Vec<Vec<String>>,
}

/// One sample from a moving frame (Frenet, Bishop, or Darboux) along a curve.
#[derive(Clone, Debug, PartialEq)]
pub struct FrameSample {
    /// Curve parameter value at this sample.
    pub parameter: f32,
    pub position: [f32; 3],
    pub tangent: [f32; 3],
    pub normal: [f32; 3],
    pub binormal: [f32; 3],
}

/// A sampled moving frame field along a curve, used to drive animated frame overlays.
#[derive(Clone, Debug, PartialEq)]
pub struct FrameField {
    pub title: String,
    pub source_plot: String,
    pub frame_kind: AnalysisKind,
    pub samples: Vec<FrameSample>,
}

/// Result returned by [`run_analysis`].
#[derive(Clone, Debug)]
pub enum AnalysisOutput {
    DerivedPlots {
        plots: Vec<PlotSpec>,
        provenance: AnalysisProvenance,
    },
    Report {
        report: AnalysisReport,
        provenance: AnalysisProvenance,
    },
    Table {
        table: AnalysisTable,
        provenance: AnalysisProvenance,
    },
    Composite {
        plots: Vec<PlotSpec>,
        reports: Vec<AnalysisReport>,
        tables: Vec<AnalysisTable>,
        diagnostics: Vec<Diagnostic>,
        frame_fields: Vec<FrameField>,
        provenance: AnalysisProvenance,
    },
}

/// Controls what kind of point groups [`sample_groups`] extracts from a plot.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum SampleGroupsKind {
    /// Evaluate the plot as a parametric curve.
    Curve,
    /// Return the raw polyline segments from the plot.
    Polyline,
    /// Return the source control points for interpolation.
    InterpolationSource,
    /// Return imported sample data as-is.
    SampleData,
}

/// Error returned by analysis functions when an operation is unsupported for the
/// given plot type, or when input data fails validation.
#[derive(Clone, Debug)]
pub struct AnalysisError {
    pub diagnostic: Diagnostic,
}

impl AnalysisError {
    pub fn unsupported(message: impl Into<String>) -> Self {
        Self {
            diagnostic: Diagnostic::error(DiagnosticKind::Build, message),
        }
    }

    pub fn invalid(message: impl Into<String>) -> Self {
        Self {
            diagnostic: Diagnostic::error(DiagnosticKind::Validation, message),
        }
    }
}

impl std::fmt::Display for AnalysisError {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        write!(f, "{}", self.diagnostic)
    }
}

impl std::error::Error for AnalysisError {}

/// Returns the analyses available for a given plot based on its definition type and metadata.
pub fn available_analyses(plot: &PlotSpec) -> Vec<AnalysisCapability> {
    let metadata = plot.metadata();
    let mut capabilities = capabilities_for_metadata(&metadata);
    if supports_point_cloud_statistics(plot) {
        capabilities.push(AnalysisCapability {
            kind: AnalysisKind::PointCloudStatistics,
            target_kind: AnalysisTargetKind::SampledData,
            output_kind: AnalysisOutputKind::Composite,
            parameters: vec![],
        });
    }
    if supports_data_quality_analysis(plot) {
        capabilities.push(AnalysisCapability {
            kind: AnalysisKind::DataQualityChecks,
            target_kind: AnalysisTargetKind::SampledData,
            output_kind: AnalysisOutputKind::Composite,
            parameters: vec![],
        });
    }
    capabilities
}

/// Extract sampled point groups from a plot for use in curve analysis.
pub fn sample_groups(
    plot: &PlotSpec,
    kind: SampleGroupsKind,
) -> Result<Vec<Vec<[f32; 3]>>, AnalysisError> {
    match kind {
        SampleGroupsKind::Curve => curve_sample_groups(plot),
        SampleGroupsKind::Polyline => polyline_sample_groups(plot),
        SampleGroupsKind::InterpolationSource => interpolation_source_groups(plot),
        SampleGroupsKind::SampleData => sample_data_groups(plot),
    }
}

/// Run an analysis on a plot and return the result.
///
/// Call [`available_analyses`] first to check which analyses are valid for the plot.
pub fn run_analysis(
    plot: &PlotSpec,
    request: &AnalysisRequest,
) -> Result<AnalysisOutput, AnalysisError> {
    let params = parameter_map(&request.parameters);
    let provenance = AnalysisProvenance {
        kind: request.kind,
        source_plots: vec![plot.name.clone()],
        parameters: request.parameters.clone(),
        notes: Vec::new(),
    };

    let plots = match request.kind {
        AnalysisKind::PointCloudStatistics => return make_point_cloud_statistics_output(plot),
        AnalysisKind::DataQualityChecks => return make_data_quality_output(plot),
        AnalysisKind::ScalarSlice => vec![make_scalar_slice_plot(
            plot,
            parse_axis(params.get("axis").map(String::as_str)).unwrap_or(SliceAxis::Z),
            params
                .get("position")
                .and_then(|value| value.parse::<f64>().ok()),
            params
                .get("contours")
                .and_then(|value| value.parse::<usize>().ok()),
        )?],
        AnalysisKind::VectorSlice => vec![make_vector_slice_plot(
            plot,
            parse_axis(params.get("axis").map(String::as_str)).unwrap_or(SliceAxis::Z),
            params
                .get("position")
                .and_then(|value| value.parse::<f64>().ok()),
        )?],
        AnalysisKind::GradientField => vec![make_gradient_plot(plot)?],
        AnalysisKind::DivergenceField => vec![make_divergence_plot(plot)?],
        AnalysisKind::CurlField => vec![make_curl_plot(plot)?],
        AnalysisKind::DifferentiateCurve => vec![make_curve_derivative_plot(plot)?],
        AnalysisKind::FitCurve => {
            return make_curve_fit_output(plot, build_curve_fit_options(&params));
        }
        AnalysisKind::AxisDerivativeCurve => vec![make_axis_derivative_plot(
            plot,
            parse_axis_index(params.get("numerator_axis").map(String::as_str)).unwrap_or(1),
            parse_axis_index(params.get("denominator_axis").map(String::as_str)).unwrap_or(0),
            params.get("output_name").cloned(),
        )?],
        AnalysisKind::IntegralCurve => vec![make_curve_integral_plot(
            plot,
            params
                .get("normalize_integral")
                .is_none_or(|value| matches!(value.as_str(), "1" | "true" | "yes")),
        )?],
        AnalysisKind::ArcLengthCurve => vec![make_curve_arc_length_plot(plot)?],
        AnalysisKind::CurvatureCurve => vec![make_curve_curvature_plot(plot)?],
        AnalysisKind::TangentField => vec![make_curve_tangent_plot(plot)?],
        AnalysisKind::NormalField => vec![make_curve_normal_plot(plot)?],
        AnalysisKind::BinormalField => vec![make_curve_binormal_plot(plot)?],
        AnalysisKind::FrenetFrame => return make_curve_frame_output(plot, request.kind, &params),
        AnalysisKind::BishopFrame => return make_curve_frame_output(plot, request.kind, &params),
        AnalysisKind::DarbouxFrame
        | AnalysisKind::SurfaceAlignedFrame
        | AnalysisKind::CurveSurfaceMeasurement => {
            return Err(AnalysisError::unsupported(
                "Curve-surface analyses require both curve and surface context.",
            ));
        }
        AnalysisKind::ExtractPoints => vec![make_extracted_points_plot(plot)?],
        AnalysisKind::InterpolateCurve => make_interpolated_plots(
            plot,
            build_interpolation(&params),
            params.get("output_name").cloned(),
        )?,
        AnalysisKind::SurfaceNormals
        | AnalysisKind::SurfaceCurvature
        | AnalysisKind::SurfaceArea
        | AnalysisKind::SurfaceMeshQuality => {
            return Err(AnalysisError::unsupported(
                "Surface geometry analyses operate on cached sampled surface meshes.",
            ));
        }
        AnalysisKind::SurfaceIntersection => {
            return Err(AnalysisError::unsupported(
                "Surface intersection remains a geometry-level app workflow.",
            ));
        }
    };

    Ok(AnalysisOutput::DerivedPlots { plots, provenance })
}

fn supports_point_cloud_statistics(plot: &PlotSpec) -> bool {
    sample_data_groups(plot)
        .map(|groups| groups.iter().map(Vec::len).sum::<usize>() >= 1)
        .unwrap_or(false)
}

fn supports_data_quality_analysis(plot: &PlotSpec) -> bool {
    sample_data_groups(plot)
        .map(|groups| groups.iter().map(Vec::len).sum::<usize>() >= 2)
        .unwrap_or(false)
}

pub fn run_surface_mesh_analysis(
    source: &PlotSpec,
    kind: AnalysisKind,
    meshes: &[&MeshData],
    parameters: &[(String, String)],
) -> Result<AnalysisOutput, AnalysisError> {
    let combined = CombinedSurfaceMesh::from_meshes(meshes)?;
    let params = parameter_map(parameters);
    match kind {
        AnalysisKind::SurfaceNormals => make_surface_normals_output(
            source,
            &combined,
            params
                .get("max_samples")
                .and_then(|value| value.parse::<usize>().ok())
                .unwrap_or(256),
            params
                .get("vector_scale")
                .and_then(|value| value.parse::<f32>().ok())
                .unwrap_or(1.0),
        ),
        AnalysisKind::SurfaceCurvature => make_surface_curvature_output(source, &combined, &params),
        AnalysisKind::SurfaceArea => make_surface_area_output(source, &combined),
        AnalysisKind::SurfaceMeshQuality => make_surface_mesh_quality_output(source, &combined),
        _ => Err(AnalysisError::unsupported(
            "Not a surface-mesh analysis kind.",
        )),
    }
}

pub fn run_curve_surface_frame_analysis(
    curve_source: &PlotSpec,
    surface_source_name: &str,
    kind: AnalysisKind,
    meshes: &[&MeshData],
    parameters: &[(String, String)],
) -> Result<AnalysisOutput, AnalysisError> {
    let combined = CombinedSurfaceMesh::from_meshes(meshes)?;
    let params = parameter_map(parameters);
    let groups = curve_sample_groups(curve_source)?;
    let frame_fields = groups
        .iter()
        .enumerate()
        .filter_map(|(index, group)| {
            build_curve_surface_frame_field(
                group,
                curve_source,
                surface_source_name,
                kind,
                &combined,
                index,
                groups.len(),
            )
            .ok()
        })
        .collect::<Vec<_>>();
    if frame_fields.is_empty() {
        return Err(AnalysisError::invalid(
            "Surface-coupled frame analysis requires a sampled curve and a usable sampled surface mesh.",
        ));
    }

    let _max_samples = params
        .get("max_samples")
        .and_then(|value| value.parse::<usize>().ok())
        .unwrap_or(128);
    let _vector_scale = params
        .get("vector_scale")
        .and_then(|value| value.parse::<f32>().ok())
        .unwrap_or(1.0)
        .max(0.01);

    let reports = frame_fields
        .iter()
        .map(|field| frame_field_report(field))
        .collect::<Vec<_>>();
    let tables = frame_fields
        .iter()
        .map(|field| frame_field_table(field))
        .collect::<Vec<_>>();

    Ok(AnalysisOutput::Composite {
        plots: Vec::new(),
        reports,
        tables,
        diagnostics: Vec::new(),
        frame_fields,
        provenance: AnalysisProvenance {
            kind,
            source_plots: vec![curve_source.name.clone(), surface_source_name.to_string()],
            parameters: parameters.to_vec(),
            notes: vec![
                "Surface-coupled frame normals are sampled from the nearest cached surface vertices.".to_string(),
            ],
        },
    })
}

pub fn run_curve_surface_measurement_analysis(
    curve_source: &PlotSpec,
    surface_source_name: &str,
    meshes: &[&MeshData],
    parameters: &[(String, String)],
) -> Result<AnalysisOutput, AnalysisError> {
    let combined = CombinedSurfaceMesh::from_meshes(meshes)?;
    let params = parameter_map(parameters);
    let max_samples = params
        .get("max_samples")
        .and_then(|value| value.parse::<usize>().ok())
        .unwrap_or(512)
        .max(2);
    let vector_scale = params
        .get("vector_scale")
        .and_then(|value| value.parse::<f32>().ok())
        .unwrap_or(1.0)
        .max(0.05);
    make_curve_surface_measurement_output(
        curve_source,
        surface_source_name,
        &combined,
        max_samples,
        vector_scale,
    )
}

fn capabilities_for_metadata(metadata: &PlotMetadata) -> Vec<AnalysisCapability> {
    let mut capabilities = Vec::new();

    if metadata.style_caps.line {
        capabilities.extend([
            AnalysisCapability {
                kind: AnalysisKind::InterpolateCurve,
                target_kind: AnalysisTargetKind::SampledData,
                output_kind: AnalysisOutputKind::PlotSpec,
                parameters: vec![
                    "output_name",
                    "interpolation_kind",
                    "samples_per_segment",
                    "closed",
                    "smoothing_window",
                ],
            },
            AnalysisCapability {
                kind: AnalysisKind::FitCurve,
                target_kind: AnalysisTargetKind::SampledData,
                output_kind: AnalysisOutputKind::Composite,
                parameters: vec![
                    "fit_method",
                    "output_name",
                    "degree",
                    "harmonics",
                    "smoothing_window",
                    "samples_per_segment",
                    "show_control_points",
                    "show_residual_plot",
                ],
            },
            AnalysisCapability {
                kind: AnalysisKind::DifferentiateCurve,
                target_kind: AnalysisTargetKind::Definition,
                output_kind: AnalysisOutputKind::PlotSpec,
                parameters: vec![],
            },
            AnalysisCapability {
                kind: AnalysisKind::AxisDerivativeCurve,
                target_kind: AnalysisTargetKind::SampledData,
                output_kind: AnalysisOutputKind::PlotSpec,
                parameters: vec!["numerator_axis", "denominator_axis", "output_name"],
            },
            AnalysisCapability {
                kind: AnalysisKind::IntegralCurve,
                target_kind: AnalysisTargetKind::Definition,
                output_kind: AnalysisOutputKind::PlotSpec,
                parameters: vec!["normalize_integral"],
            },
            AnalysisCapability {
                kind: AnalysisKind::ArcLengthCurve,
                target_kind: AnalysisTargetKind::SampledData,
                output_kind: AnalysisOutputKind::PlotSpec,
                parameters: vec![],
            },
            AnalysisCapability {
                kind: AnalysisKind::CurvatureCurve,
                target_kind: AnalysisTargetKind::SampledData,
                output_kind: AnalysisOutputKind::PlotSpec,
                parameters: vec![],
            },
            AnalysisCapability {
                kind: AnalysisKind::TangentField,
                target_kind: AnalysisTargetKind::SampledData,
                output_kind: AnalysisOutputKind::PlotSpec,
                parameters: vec![],
            },
            AnalysisCapability {
                kind: AnalysisKind::NormalField,
                target_kind: AnalysisTargetKind::SampledData,
                output_kind: AnalysisOutputKind::PlotSpec,
                parameters: vec![],
            },
            AnalysisCapability {
                kind: AnalysisKind::BinormalField,
                target_kind: AnalysisTargetKind::SampledData,
                output_kind: AnalysisOutputKind::PlotSpec,
                parameters: vec![],
            },
            AnalysisCapability {
                kind: AnalysisKind::ExtractPoints,
                target_kind: AnalysisTargetKind::SampledData,
                output_kind: AnalysisOutputKind::PlotSpec,
                parameters: vec![],
            },
            AnalysisCapability {
                kind: AnalysisKind::FrenetFrame,
                target_kind: AnalysisTargetKind::SampledData,
                output_kind: AnalysisOutputKind::Composite,
                parameters: vec!["max_samples", "vector_scale"],
            },
            AnalysisCapability {
                kind: AnalysisKind::BishopFrame,
                target_kind: AnalysisTargetKind::SampledData,
                output_kind: AnalysisOutputKind::Composite,
                parameters: vec!["max_samples", "vector_scale"],
            },
            AnalysisCapability {
                kind: AnalysisKind::DarbouxFrame,
                target_kind: AnalysisTargetKind::PlotPair,
                output_kind: AnalysisOutputKind::Composite,
                parameters: vec!["max_samples", "vector_scale"],
            },
            AnalysisCapability {
                kind: AnalysisKind::SurfaceAlignedFrame,
                target_kind: AnalysisTargetKind::PlotPair,
                output_kind: AnalysisOutputKind::Composite,
                parameters: vec!["max_samples", "vector_scale"],
            },
            AnalysisCapability {
                kind: AnalysisKind::CurveSurfaceMeasurement,
                target_kind: AnalysisTargetKind::PlotPair,
                output_kind: AnalysisOutputKind::Composite,
                parameters: vec!["max_samples", "vector_scale"],
            },
        ]);
    }

    if metadata.coordinate_semantics == crate::CoordinateSemantics::CartesianVolume {
        capabilities.extend([
            AnalysisCapability {
                kind: AnalysisKind::ScalarSlice,
                target_kind: AnalysisTargetKind::Definition,
                output_kind: AnalysisOutputKind::PlotSpec,
                parameters: vec!["axis", "position", "contours"],
            },
            AnalysisCapability {
                kind: AnalysisKind::VectorSlice,
                target_kind: AnalysisTargetKind::Definition,
                output_kind: AnalysisOutputKind::PlotSpec,
                parameters: vec!["axis", "position"],
            },
        ]);
    }

    if metadata.required_variables == ["x".to_string(), "y".to_string(), "z".to_string()] {
        capabilities.extend([
            AnalysisCapability {
                kind: AnalysisKind::GradientField,
                target_kind: AnalysisTargetKind::Definition,
                output_kind: AnalysisOutputKind::PlotSpec,
                parameters: vec![],
            },
            AnalysisCapability {
                kind: AnalysisKind::DivergenceField,
                target_kind: AnalysisTargetKind::Definition,
                output_kind: AnalysisOutputKind::PlotSpec,
                parameters: vec![],
            },
            AnalysisCapability {
                kind: AnalysisKind::CurlField,
                target_kind: AnalysisTargetKind::Definition,
                output_kind: AnalysisOutputKind::PlotSpec,
                parameters: vec![],
            },
        ]);
    }

    if metadata.supports_surface_intersection {
        capabilities.extend([
            AnalysisCapability {
                kind: AnalysisKind::SurfaceNormals,
                target_kind: AnalysisTargetKind::Geometry,
                output_kind: AnalysisOutputKind::PlotSpec,
                parameters: vec![],
            },
            AnalysisCapability {
                kind: AnalysisKind::SurfaceCurvature,
                target_kind: AnalysisTargetKind::Geometry,
                output_kind: AnalysisOutputKind::Composite,
                parameters: vec![],
            },
            AnalysisCapability {
                kind: AnalysisKind::SurfaceArea,
                target_kind: AnalysisTargetKind::Geometry,
                output_kind: AnalysisOutputKind::NumericReport,
                parameters: vec![],
            },
            AnalysisCapability {
                kind: AnalysisKind::SurfaceMeshQuality,
                target_kind: AnalysisTargetKind::Geometry,
                output_kind: AnalysisOutputKind::Composite,
                parameters: vec![],
            },
            AnalysisCapability {
                kind: AnalysisKind::SurfaceIntersection,
                target_kind: AnalysisTargetKind::PlotPair,
                output_kind: AnalysisOutputKind::PlotSpec,
                parameters: vec!["samples", "tolerance"],
            },
        ]);
    }

    capabilities
}

fn make_scalar_slice_plot(
    source: &PlotSpec,
    axis: SliceAxis,
    position: Option<f64>,
    contour_count: Option<usize>,
) -> Result<PlotSpec, AnalysisError> {
    let (expression, parameters) = match &source.definition {
        crate::PlotDefinition::ExprVolume {
            expression,
            parameters,
            ..
        }
        | crate::PlotDefinition::ExprIsosurface {
            expression,
            parameters,
            ..
        } => (expression.clone(), parameters.clone()),
        _ => {
            return Err(AnalysisError::unsupported(
                "Scalar slices require a scalar volume or isosurface source.",
            ));
        }
    };
    Ok(PlotSpec {
        name: format!("{} Slice {}", axis.label(), source.name),
        visible: true,
        domain: source.domain.clone(),
        resolution: source.resolution,
        style: PlotStyle {
            colour_mode: ColourMode::ByAttribute {
                name: "value".to_string(),
                kind: viewport_lib::AttributeKind::Vertex,
            },
            two_sided: true,
            ..source.style.clone()
        },
        definition: crate::PlotDefinition::ScalarSlice {
            expression,
            parameters,
            axis,
            position: position.unwrap_or_else(|| default_slice_position(&source.domain, axis)),
            contour_values: evenly_spaced_isovalues(contour_count.unwrap_or(8)),
            contour_style: PlotStyle {
                colour_mode: ColourMode::Solid([1.0, 0.95, 0.35, 1.0]),
                line_width: 2.0,
                ..PlotStyle::default()
            },
        },
    })
}

fn make_vector_slice_plot(
    source: &PlotSpec,
    axis: SliceAxis,
    position: Option<f64>,
) -> Result<PlotSpec, AnalysisError> {
    let (expression, parameters) = match &source.definition {
        crate::PlotDefinition::ExprVectorField {
            expression,
            parameters,
        } => (expression.clone(), parameters.clone()),
        _ => {
            return Err(AnalysisError::unsupported(
                "Vector slices require a vector field source.",
            ));
        }
    };
    Ok(PlotSpec {
        name: format!("{} Slice {}", axis.label(), source.name),
        visible: true,
        domain: source.domain.clone(),
        resolution: source.resolution,
        style: source.style.clone(),
        definition: crate::PlotDefinition::VectorSlice {
            expression,
            parameters,
            axis,
            position: position.unwrap_or_else(|| default_slice_position(&source.domain, axis)),
        },
    })
}

fn make_gradient_plot(source: &PlotSpec) -> Result<PlotSpec, AnalysisError> {
    let (expression, parameters) = match &source.definition {
        crate::PlotDefinition::ExprVolume {
            expression,
            parameters,
            ..
        }
        | crate::PlotDefinition::ExprIsosurface {
            expression,
            parameters,
            ..
        } => (expression.clone(), parameters.clone()),
        _ => {
            return Err(AnalysisError::unsupported(
                "Gradient plots require a scalar volume or isosurface source.",
            ));
        }
    };
    Ok(PlotSpec {
        name: format!("Gradient {}", source.name),
        visible: true,
        domain: source.domain.clone(),
        resolution: source.resolution,
        style: PlotStyle {
            colour_mode: ColourMode::ByAttribute {
                name: "magnitude".to_string(),
                kind: viewport_lib::AttributeKind::Vertex,
            },
            glyph_scale: 0.8,
            shading: crate::ShadingMode::Unlit,
            ..PlotStyle::default()
        },
        definition: crate::PlotDefinition::GradientField {
            expression,
            parameters,
        },
    })
}

fn make_divergence_plot(source: &PlotSpec) -> Result<PlotSpec, AnalysisError> {
    let (expression, parameters) = match &source.definition {
        crate::PlotDefinition::ExprVectorField {
            expression,
            parameters,
        } => (expression.clone(), parameters.clone()),
        _ => {
            return Err(AnalysisError::unsupported(
                "Divergence plots require a vector field source.",
            ));
        }
    };
    Ok(PlotSpec {
        name: format!("Divergence {}", source.name),
        visible: true,
        domain: source.domain.clone(),
        resolution: source.resolution,
        style: PlotStyle {
            opacity: 0.3,
            transfer_function: Some(crate::TransferFunction {
                opacity_scale: 0.4,
                threshold: None,
            }),
            ..PlotStyle::default()
        },
        definition: crate::PlotDefinition::DivergenceField {
            expression,
            parameters,
            vol_resolution: [64, 64, 64],
        },
    })
}

fn make_curl_plot(source: &PlotSpec) -> Result<PlotSpec, AnalysisError> {
    let (expression, parameters) = match &source.definition {
        crate::PlotDefinition::ExprVectorField {
            expression,
            parameters,
        } => (expression.clone(), parameters.clone()),
        _ => {
            return Err(AnalysisError::unsupported(
                "Curl plots require a vector field source.",
            ));
        }
    };
    Ok(PlotSpec {
        name: format!("Curl {}", source.name),
        visible: true,
        domain: source.domain.clone(),
        resolution: source.resolution,
        style: PlotStyle {
            colour_mode: ColourMode::ByAttribute {
                name: "magnitude".to_string(),
                kind: viewport_lib::AttributeKind::Vertex,
            },
            glyph_scale: 0.8,
            shading: crate::ShadingMode::Unlit,
            ..PlotStyle::default()
        },
        definition: crate::PlotDefinition::CurlField {
            expression,
            parameters,
        },
    })
}

fn make_curve_derivative_plot(source: &PlotSpec) -> Result<PlotSpec, AnalysisError> {
    let groups = curve_sample_groups(source)?;
    let derived_groups = match &source.definition {
        crate::PlotDefinition::ExprCartesianLine {
            dep_var, ind_var, ..
        } => groups
            .iter()
            .map(|group| derivative_cartesian_line_group(group, dep_var.as_str(), ind_var.as_str()))
            .filter(|group| group.len() >= 2)
            .collect(),
        _ => groups
            .iter()
            .map(|group| derivative_curve_group(group))
            .filter(|group| group.len() >= 2)
            .collect(),
    };
    derived_polyline_plot(
        source,
        format!("Derivative {}", source.name),
        [1.0, 0.55, 0.25, 1.0],
        derived_groups,
    )
}

fn make_curve_tangent_plot(source: &PlotSpec) -> Result<PlotSpec, AnalysisError> {
    let groups = curve_sample_groups(source)?;
    let derived_groups = groups
        .iter()
        .map(|group| tangent_curve_group(group))
        .filter(|group| group.len() >= 2)
        .collect();
    derived_polyline_plot(
        source,
        format!("Tangent {}", source.name),
        [0.25, 0.85, 0.45, 1.0],
        derived_groups,
    )
}

fn make_curve_integral_plot(
    source: &PlotSpec,
    normalize_integral: bool,
) -> Result<PlotSpec, AnalysisError> {
    let groups = curve_sample_groups(source)?;
    let derived_groups = match &source.definition {
        crate::PlotDefinition::ExprCartesianLine {
            dep_var, ind_var, ..
        } => groups
            .iter()
            .map(|group| {
                integral_cartesian_line_group(
                    group,
                    dep_var.as_str(),
                    ind_var.as_str(),
                    normalize_integral,
                )
            })
            .filter(|group| group.len() >= 2)
            .collect(),
        _ => groups
            .iter()
            .map(|group| integral_curve_group(group, normalize_integral))
            .filter(|group| group.len() >= 2)
            .collect(),
    };
    derived_polyline_plot(
        source,
        format!("Integral {}", source.name),
        [0.45, 0.7, 1.0, 1.0],
        derived_groups,
    )
}

fn make_curve_arc_length_plot(source: &PlotSpec) -> Result<PlotSpec, AnalysisError> {
    let groups = curve_sample_groups(source)?;
    let derived_groups = match &source.definition {
        crate::PlotDefinition::ExprCartesianLine {
            dep_var, ind_var, ..
        } => groups
            .iter()
            .map(|group| {
                scalar_plot_cartesian_line_group(
                    group,
                    dep_var.as_str(),
                    ind_var.as_str(),
                    &cumulative_arc_lengths(group),
                )
            })
            .filter(|group| group.len() >= 2)
            .collect(),
        _ => groups
            .iter()
            .map(|group| scalar_curve_group(group, &cumulative_arc_lengths(group)))
            .filter(|group| group.len() >= 2)
            .collect(),
    };
    derived_polyline_plot(
        source,
        format!("Arc Length {}", source.name),
        [0.95, 0.85, 0.3, 1.0],
        derived_groups,
    )
}

fn make_curve_curvature_plot(source: &PlotSpec) -> Result<PlotSpec, AnalysisError> {
    let groups = curve_sample_groups(source)?;
    let derived_groups = match &source.definition {
        crate::PlotDefinition::ExprCartesianLine {
            dep_var, ind_var, ..
        } => groups
            .iter()
            .map(|group| {
                scalar_plot_cartesian_line_group(
                    group,
                    dep_var.as_str(),
                    ind_var.as_str(),
                    &curvature_values(group),
                )
            })
            .filter(|group| group.len() >= 2)
            .collect(),
        _ => groups
            .iter()
            .map(|group| scalar_curve_group(group, &curvature_values(group)))
            .filter(|group| group.len() >= 2)
            .collect(),
    };
    derived_polyline_plot(
        source,
        format!("Curvature {}", source.name),
        [0.8, 0.45, 1.0, 1.0],
        derived_groups,
    )
}

fn make_curve_normal_plot(source: &PlotSpec) -> Result<PlotSpec, AnalysisError> {
    let groups = curve_sample_groups(source)?;
    let derived_groups = groups
        .iter()
        .map(|group| normal_curve_group(group))
        .filter(|group| group.len() >= 2)
        .collect();
    derived_polyline_plot(
        source,
        format!("Normal {}", source.name),
        [0.3, 0.8, 1.0, 1.0],
        derived_groups,
    )
}

fn make_curve_binormal_plot(source: &PlotSpec) -> Result<PlotSpec, AnalysisError> {
    let groups = curve_sample_groups(source)?;
    let derived_groups = groups
        .iter()
        .map(|group| binormal_curve_group(group))
        .filter(|group| group.len() >= 2)
        .collect();
    derived_polyline_plot(
        source,
        format!("Binormal {}", source.name),
        [1.0, 0.45, 0.65, 1.0],
        derived_groups,
    )
}

fn make_curve_frame_output(
    source: &PlotSpec,
    frame_kind: AnalysisKind,
    params: &HashMap<String, String>,
) -> Result<AnalysisOutput, AnalysisError> {
    let groups = curve_sample_groups(source)?;
    let frame_fields = groups
        .iter()
        .enumerate()
        .filter_map(|(index, group)| {
            build_curve_frame_field(group, source, frame_kind, index, groups.len()).ok()
        })
        .collect::<Vec<_>>();
    if frame_fields.is_empty() {
        return Err(AnalysisError::invalid(
            "Moving-frame analysis requires a sampled curve with at least two distinct points.",
        ));
    }

    let max_samples = params
        .get("max_samples")
        .and_then(|value| value.parse::<usize>().ok())
        .unwrap_or(128);
    let vector_scale = params
        .get("vector_scale")
        .and_then(|value| value.parse::<f32>().ok())
        .unwrap_or(1.0)
        .max(0.01);

    let reports = frame_fields
        .iter()
        .map(|field| frame_field_report(field))
        .collect::<Vec<_>>();
    let tables = frame_fields
        .iter()
        .map(|field| frame_field_table(field))
        .collect::<Vec<_>>();

    Ok(AnalysisOutput::Composite {
        plots: Vec::new(),
        reports,
        tables,
        diagnostics: Vec::new(),
        frame_fields,
        provenance: AnalysisProvenance {
            kind: frame_kind,
            source_plots: vec![source.name.clone()],
            parameters: vec![
                ("max_samples".to_string(), max_samples.to_string()),
                ("vector_scale".to_string(), format!("{vector_scale:.4}")),
            ],
            notes: vec![
                "Frame parameters are cumulative sampled arc length values.".to_string(),
                "Bishop frames use rotation-minimizing transport across sampled tangents."
                    .to_string(),
            ],
        },
    })
}

fn make_axis_derivative_plot(
    source: &PlotSpec,
    numerator_axis: usize,
    denominator_axis: usize,
    output_name: Option<String>,
) -> Result<PlotSpec, AnalysisError> {
    if numerator_axis == denominator_axis {
        return Err(AnalysisError::invalid(
            "Numerator and denominator axes must be different.",
        ));
    }
    let groups = curve_sample_groups(source)?;
    let derived_groups: Vec<Vec<[f32; 3]>> = groups
        .iter()
        .map(|group| axis_derivative_group(group, numerator_axis, denominator_axis))
        .filter(|group| group.len() >= 2)
        .collect();
    if derived_groups.is_empty() {
        return Err(AnalysisError::invalid(
            "Could not compute an axis derivative from the selected curve.",
        ));
    }
    derived_polyline_plot(
        source,
        output_name.unwrap_or_else(|| {
            format!(
                "d{}/d{} {}",
                axis_name(numerator_axis),
                axis_name(denominator_axis),
                source.name
            )
        }),
        [1.0, 0.7, 0.25, 1.0],
        derived_groups,
    )
}

fn make_extracted_points_plot(source: &PlotSpec) -> Result<PlotSpec, AnalysisError> {
    let positions = polyline_sample_groups(source)?
        .into_iter()
        .flatten()
        .collect::<Vec<_>>();
    Ok(PlotSpec {
        name: format!("Points from {}", source.name),
        visible: true,
        domain: source.domain.clone(),
        resolution: source.resolution,
        style: PlotStyle {
            colour_mode: ColourMode::Solid([0.35, 0.85, 1.0, 1.0]),
            point_size: 8.0,
            ..PlotStyle::default()
        },
        definition: crate::PlotDefinition::PointAnnotations {
            points: make_point_annotations(&positions, "Point"),
            show_labels: false,
        },
    })
}

fn make_interpolated_plots(
    source: &PlotSpec,
    interpolation: CurveInterpolation,
    output_name: Option<String>,
) -> Result<Vec<PlotSpec>, AnalysisError> {
    let groups = interpolation_source_groups(source)?;
    if groups.is_empty() || groups.iter().all(Vec::is_empty) {
        return Err(AnalysisError::invalid(
            "The selected plot does not have usable point samples.",
        ));
    }
    if groups.iter().all(|group| group.len() < 2) {
        return Err(AnalysisError::invalid(
            "At least two points are required to interpolate a curve.",
        ));
    }

    let base_name = output_name.unwrap_or_else(|| format!("Interpolated {}", source.name));
    let style = PlotStyle {
        colour_mode: ColourMode::Solid([0.95, 0.7, 0.2, 1.0]),
        line_width: 2.5,
        ..PlotStyle::default()
    };

    if groups.len() == 1 {
        let points = groups.into_iter().next().unwrap_or_default();
        return Ok(vec![PlotSpec {
            name: base_name,
            visible: true,
            domain: source.domain.clone(),
            resolution: source.resolution,
            style,
            definition: crate::PlotDefinition::InterpolatedCurve {
                points,
                interpolation,
            },
        }]);
    }

    Ok(groups
        .into_iter()
        .enumerate()
        .filter(|(_, group)| group.len() >= 2)
        .map(|(index, group)| PlotSpec {
            name: format!("{base_name} {}", index + 1),
            visible: true,
            domain: source.domain.clone(),
            resolution: source.resolution,
            style: style.clone(),
            definition: crate::PlotDefinition::InterpolatedCurve {
                points: group,
                interpolation,
            },
        })
        .collect())
}

fn interpolation_source_groups(plot: &PlotSpec) -> Result<Vec<Vec<[f32; 3]>>, AnalysisError> {
    match &plot.definition {
        crate::PlotDefinition::PointAnnotations { points, .. } => {
            Ok(vec![points.iter().map(|point| point.position).collect()])
        }
        crate::PlotDefinition::ExprCurve { .. }
        | crate::PlotDefinition::ExprCartesianLine { .. }
        | crate::PlotDefinition::HelixCurve => curve_sample_groups(plot),
        crate::PlotDefinition::ImportedTable { definition } => match definition.validate() {
            Ok(TableDataSet::Curve { groups, .. }) => Ok(groups
                .iter()
                .map(|group| group.iter().map(|point| point.to_array()).collect())
                .collect()),
            Ok(TableDataSet::Scatter { points, .. }) => {
                Ok(vec![points.iter().map(|point| point.to_array()).collect()])
            }
            Ok(_) => Err(AnalysisError::unsupported(
                "Interpolation is not available for this imported table target.",
            )),
            Err(errors) => Err(table_errors(errors)),
        },
        crate::PlotDefinition::DerivedPolylineGroups { groups } => Ok(groups.clone()),
        crate::PlotDefinition::InterpolatedCurve { points, .. } => Ok(vec![points.clone()]),
        _ => Err(AnalysisError::unsupported(
            "Interpolation is available for point and ordered sample plots.",
        )),
    }
}

fn polyline_sample_groups(plot: &PlotSpec) -> Result<Vec<Vec<[f32; 3]>>, AnalysisError> {
    match &plot.definition {
        crate::PlotDefinition::ExprCurve { .. }
        | crate::PlotDefinition::ExprCartesianLine { .. }
        | crate::PlotDefinition::HelixCurve => curve_sample_groups(plot),
        crate::PlotDefinition::ImportedTable { definition } => match definition.validate() {
            Ok(TableDataSet::Curve { groups, .. }) => Ok(groups
                .iter()
                .map(|group| group.iter().map(|point| point.to_array()).collect())
                .collect()),
            Ok(_) => Err(AnalysisError::unsupported(
                "Point extraction is only available for imported curve tables.",
            )),
            Err(errors) => Err(table_errors(errors)),
        },
        crate::PlotDefinition::DerivedPolylineGroups { groups } => Ok(groups.clone()),
        crate::PlotDefinition::InterpolatedCurve {
            points,
            interpolation,
        } => {
            let sampled = sample_curve_points(
                &points
                    .iter()
                    .map(|point| Vec3::from_array(*point))
                    .collect::<Vec<_>>(),
                *interpolation,
            );
            Ok(vec![
                sampled.into_iter().map(|point| point.to_array()).collect(),
            ])
        }
        _ => Err(AnalysisError::unsupported(
            "Point extraction is available for polyline and interpolated curve plots.",
        )),
    }
}

fn sample_data_groups(plot: &PlotSpec) -> Result<Vec<Vec<[f32; 3]>>, AnalysisError> {
    match &plot.definition {
        crate::PlotDefinition::ScatterCloud => Ok(vec![scatter_cloud_points()]),
        crate::PlotDefinition::PointAnnotations { points, .. } => {
            Ok(vec![points.iter().map(|point| point.position).collect()])
        }
        crate::PlotDefinition::ArrowAnnotations { arrows, .. } => {
            Ok(vec![arrows.iter().map(|arrow| arrow.origin).collect()])
        }
        crate::PlotDefinition::ExprCurve { .. }
        | crate::PlotDefinition::ExprCartesianLine { .. }
        | crate::PlotDefinition::HelixCurve => curve_sample_groups(plot),
        crate::PlotDefinition::ImportedTable { definition } => match definition.validate() {
            Ok(TableDataSet::Curve { groups, .. }) => Ok(groups
                .iter()
                .map(|group| group.iter().map(|point| point.to_array()).collect())
                .collect()),
            Ok(TableDataSet::Scatter { points, .. }) => {
                Ok(vec![points.iter().map(|point| point.to_array()).collect()])
            }
            Ok(TableDataSet::VectorField { samples, .. }) => Ok(vec![
                samples
                    .iter()
                    .map(|sample| sample.position.to_array())
                    .collect(),
            ]),
            Ok(TableDataSet::SurfaceGrid { xs, ys, zs }) => {
                if xs.is_empty() || ys.is_empty() || zs.is_empty() {
                    return Err(AnalysisError::invalid(
                        "Imported surface grid does not contain enough samples.",
                    ));
                }
                let ny = ys.len();
                let zs_ref = &zs;
                let points = xs
                    .iter()
                    .enumerate()
                    .flat_map(|(ix, x)| {
                        ys.iter().enumerate().filter_map(move |(iy, y)| {
                            let index = ix.checked_mul(ny)?.checked_add(iy)?;
                            let z = *zs_ref.get(index)?;
                            Some([*x as f32, *y as f32, z as f32])
                        })
                    })
                    .collect::<Vec<_>>();
                Ok(vec![points])
            }
            Err(errors) => Err(table_errors(errors)),
        },
        crate::PlotDefinition::DerivedPolylineGroups { groups } => Ok(groups.clone()),
        crate::PlotDefinition::InterpolatedCurve {
            points,
            interpolation,
        } => {
            let sampled = sample_curve_points(
                &points
                    .iter()
                    .map(|point| Vec3::from_array(*point))
                    .collect::<Vec<_>>(),
                *interpolation,
            );
            Ok(vec![
                sampled.into_iter().map(|point| point.to_array()).collect(),
            ])
        }
        _ => Err(AnalysisError::unsupported(
            "Sample statistics require point-like or ordered sampled data.",
        )),
    }
}

fn curve_sample_groups(plot: &PlotSpec) -> Result<Vec<Vec<[f32; 3]>>, AnalysisError> {
    match &plot.definition {
        crate::PlotDefinition::HelixCurve => {
            let steps = plot.resolution.u.max(2) as usize;
            let points = (0..steps)
                .map(|i| {
                    let t = 20.0 * PI * i as f64 / (steps - 1) as f64;
                    Vec3::new(
                        (t.cos() * 3.0) as f32,
                        (t.sin() * 3.0) as f32,
                        (t * 0.15) as f32,
                    )
                    .to_array()
                })
                .collect();
            Ok(vec![points])
        }
        crate::PlotDefinition::ExprCurve {
            expression,
            parameters,
            t_range,
        } => {
            let parsed = parse_curve_expr(expression).map_err(parse_error)?;
            let steps = plot.resolution.u.max(2) as usize;
            let (t0, t1) = *t_range;
            let points = (0..steps)
                .map(|i| {
                    let t = t0 + (i as f64 / (steps - 1) as f64) * (t1 - t0);
                    let p = eval_curve_point(&parsed, t, parameters);
                    Vec3::new(p.x as f32, p.y as f32, p.z as f32).to_array()
                })
                .collect();
            Ok(vec![points])
        }
        crate::PlotDefinition::ExprCartesianLine {
            dep_var,
            ind_var,
            expression,
            parameters,
        } => {
            let parsed =
                parse_expr_with_vars(expression, &[ind_var.as_str()]).map_err(parse_error)?;
            let steps = plot.resolution.u.max(2) as usize;
            let (t0, t1) = (*plot.domain.x.start(), *plot.domain.x.end());
            let dep = dep_var.clone();
            let ind = ind_var.clone();
            let points = (0..steps)
                .map(|i| {
                    let t = t0 + (i as f64 / (steps - 1) as f64) * (t1 - t0);
                    let vars: Vec<(&str, f64)> = parameters
                        .iter()
                        .map(|(n, v)| (n.as_str(), *v))
                        .chain(std::iter::once((ind.as_str(), t)))
                        .collect();
                    let val = eval_with_vars(&parsed, &vars);
                    cartesian_line_point(dep.as_str(), ind.as_str(), t as f32, val as f32)
                        .to_array()
                })
                .collect();
            Ok(vec![points])
        }
        crate::PlotDefinition::ImportedTable { definition } => match definition.validate() {
            Ok(TableDataSet::Curve { groups, .. }) => Ok(groups
                .iter()
                .map(|group| group.iter().map(|point| point.to_array()).collect())
                .collect()),
            Ok(_) => Err(AnalysisError::unsupported(
                "Curve calculus tools require curve-like sample data.",
            )),
            Err(errors) => Err(table_errors(errors)),
        },
        crate::PlotDefinition::DerivedPolylineGroups { groups } => Ok(groups.clone()),
        crate::PlotDefinition::InterpolatedCurve {
            points,
            interpolation,
        } => {
            let sampled = sample_curve_points(
                &points
                    .iter()
                    .map(|point| Vec3::from_array(*point))
                    .collect::<Vec<_>>(),
                *interpolation,
            );
            Ok(vec![
                sampled.into_iter().map(|point| point.to_array()).collect(),
            ])
        }
        _ => Err(AnalysisError::unsupported(
            "Curve calculus tools are available for curve and polyline plots.",
        )),
    }
}

fn build_curve_frame_field(
    group: &[[f32; 3]],
    source: &PlotSpec,
    frame_kind: AnalysisKind,
    group_index: usize,
    group_count: usize,
) -> Result<FrameField, AnalysisError> {
    let points = group
        .iter()
        .map(|point| Vec3::from_array(*point))
        .collect::<Vec<_>>();
    if points.len() < 2 {
        return Err(AnalysisError::invalid(
            "Moving-frame analysis requires at least two sampled points.",
        ));
    }

    let tangents = sampled_tangents(&points);
    if tangents
        .iter()
        .all(|tangent| tangent.length_squared() <= 1.0e-8)
    {
        return Err(AnalysisError::invalid(
            "Moving-frame analysis requires a curve with non-zero tangent variation.",
        ));
    }
    let (normals, binormals) = match frame_kind {
        AnalysisKind::FrenetFrame => frenet_frames(&points, &tangents),
        AnalysisKind::BishopFrame => bishop_frames(&tangents),
        _ => {
            return Err(AnalysisError::unsupported(
                "Unsupported moving-frame analysis kind.",
            ));
        }
    };
    let parameters = cumulative_arc_lengths(group);
    let title = if group_count > 1 {
        format!(
            "{} {} {}",
            frame_kind_label(frame_kind),
            source.name,
            group_index + 1
        )
    } else {
        format!("{} {}", frame_kind_label(frame_kind), source.name)
    };

    Ok(FrameField {
        title,
        source_plot: source.name.clone(),
        frame_kind,
        samples: points
            .iter()
            .enumerate()
            .map(|(index, point)| FrameSample {
                parameter: parameters.get(index).copied().unwrap_or_default(),
                position: point.to_array(),
                tangent: tangents[index].to_array(),
                normal: normals[index].to_array(),
                binormal: binormals[index].to_array(),
            })
            .collect(),
    })
}

fn build_curve_surface_frame_field(
    group: &[[f32; 3]],
    source: &PlotSpec,
    surface_name: &str,
    frame_kind: AnalysisKind,
    mesh: &CombinedSurfaceMesh,
    group_index: usize,
    group_count: usize,
) -> Result<FrameField, AnalysisError> {
    let points = group
        .iter()
        .map(|point| Vec3::from_array(*point))
        .collect::<Vec<_>>();
    if points.len() < 2 {
        return Err(AnalysisError::invalid(
            "Surface-coupled frame analysis requires at least two sampled points.",
        ));
    }
    let tangents = sampled_tangents(&points);
    let surface_normals = points
        .iter()
        .map(|point| nearest_surface_normal(mesh, *point))
        .collect::<Vec<_>>();
    let (normals, binormals) = match frame_kind {
        AnalysisKind::SurfaceAlignedFrame => surface_aligned_frames(&tangents, &surface_normals),
        AnalysisKind::DarbouxFrame => darboux_frames(&tangents, &surface_normals),
        _ => {
            return Err(AnalysisError::unsupported(
                "Unsupported surface-coupled frame analysis kind.",
            ));
        }
    };
    let parameters = cumulative_arc_lengths(group);
    let title = if group_count > 1 {
        format!(
            "{} {} on {} {}",
            frame_kind_label(frame_kind),
            source.name,
            surface_name,
            group_index + 1
        )
    } else {
        format!(
            "{} {} on {}",
            frame_kind_label(frame_kind),
            source.name,
            surface_name
        )
    };

    Ok(FrameField {
        title,
        source_plot: source.name.clone(),
        frame_kind,
        samples: points
            .iter()
            .enumerate()
            .map(|(index, point)| FrameSample {
                parameter: parameters.get(index).copied().unwrap_or_default(),
                position: point.to_array(),
                tangent: tangents[index].to_array(),
                normal: normals[index].to_array(),
                binormal: binormals[index].to_array(),
            })
            .collect(),
    })
}

fn frame_field_report(field: &FrameField) -> AnalysisReport {
    let total_length = field
        .samples
        .last()
        .map(|sample| sample.parameter)
        .unwrap_or(0.0);
    AnalysisReport {
        title: field.title.clone(),
        values: vec![
            (
                "Frame Kind".to_string(),
                frame_kind_label(field.frame_kind).to_string(),
            ),
            ("Sample Count".to_string(), field.samples.len().to_string()),
            ("Arc Length".to_string(), format_float(total_length)),
        ],
    }
}

fn frame_field_table(field: &FrameField) -> AnalysisTable {
    AnalysisTable {
        title: format!("{} Samples", field.title),
        columns: vec![
            "row".to_string(),
            "s".to_string(),
            "x".to_string(),
            "y".to_string(),
            "z".to_string(),
            "tx".to_string(),
            "ty".to_string(),
            "tz".to_string(),
            "nx".to_string(),
            "ny".to_string(),
            "nz".to_string(),
            "bx".to_string(),
            "by".to_string(),
            "bz".to_string(),
        ],
        rows: field
            .samples
            .iter()
            .enumerate()
            .map(|(index, sample)| {
                vec![
                    (index + 1).to_string(),
                    format_float(sample.parameter),
                    format_float(sample.position[0]),
                    format_float(sample.position[1]),
                    format_float(sample.position[2]),
                    format_float(sample.tangent[0]),
                    format_float(sample.tangent[1]),
                    format_float(sample.tangent[2]),
                    format_float(sample.normal[0]),
                    format_float(sample.normal[1]),
                    format_float(sample.normal[2]),
                    format_float(sample.binormal[0]),
                    format_float(sample.binormal[1]),
                    format_float(sample.binormal[2]),
                ]
            })
            .collect(),
    }
}

fn frame_kind_label(kind: AnalysisKind) -> &'static str {
    match kind {
        AnalysisKind::FrenetFrame => "Frenet Frame",
        AnalysisKind::BishopFrame => "Bishop Frame",
        AnalysisKind::DarbouxFrame => "Darboux Frame",
        AnalysisKind::SurfaceAlignedFrame => "Surface-Aligned Frame",
        _ => "Frame",
    }
}

fn sampled_tangents(points: &[Vec3]) -> Vec<Vec3> {
    (0..points.len())
        .map(|index| {
            let prev = if index == 0 {
                points[index]
            } else {
                points[index - 1]
            };
            let next = if index + 1 >= points.len() {
                points[index]
            } else {
                points[index + 1]
            };
            let tangent = if index == 0 {
                next - points[index]
            } else if index + 1 >= points.len() {
                points[index] - prev
            } else {
                next - prev
            };
            tangent.normalize_or_zero()
        })
        .collect()
}

fn frenet_frames(points: &[Vec3], tangents: &[Vec3]) -> (Vec<Vec3>, Vec<Vec3>) {
    let mut normals = Vec::with_capacity(points.len());
    let mut binormals = Vec::with_capacity(points.len());
    let mut previous_normal: Option<Vec3> = None;

    for index in 0..points.len() {
        let tangent = tangents[index];
        let prev_tangent = tangents
            .get(index.saturating_sub(1))
            .copied()
            .unwrap_or(tangent);
        let next_tangent = tangents
            .get((index + 1).min(tangents.len() - 1))
            .copied()
            .unwrap_or(tangent);
        let tangent_delta = if index == 0 {
            next_tangent - tangent
        } else if index + 1 >= tangents.len() {
            tangent - prev_tangent
        } else {
            next_tangent - prev_tangent
        };
        let mut normal = orthogonalized(tangent_delta, tangent);
        if normal.length_squared() <= 1.0e-8 {
            normal = previous_normal
                .map(|prev| orthogonalized(prev, tangent))
                .filter(|candidate| candidate.length_squared() > 1.0e-8)
                .unwrap_or_else(|| arbitrary_perpendicular(tangent));
        }
        normal = normal.normalize_or_zero();
        let mut binormal = tangent.cross(normal);
        if binormal.length_squared() <= 1.0e-8 {
            normal = arbitrary_perpendicular(tangent);
            binormal = tangent.cross(normal);
        }
        binormal = binormal.normalize_or_zero();
        normal = binormal.cross(tangent).normalize_or_zero();
        previous_normal = Some(normal);
        normals.push(normal);
        binormals.push(binormal);
    }

    (normals, binormals)
}

fn bishop_frames(tangents: &[Vec3]) -> (Vec<Vec3>, Vec<Vec3>) {
    let first_tangent = tangents
        .iter()
        .copied()
        .find(|tangent| tangent.length_squared() > 1.0e-8)
        .unwrap_or(Vec3::X);
    let mut normals = Vec::with_capacity(tangents.len());
    let mut binormals = Vec::with_capacity(tangents.len());
    let mut current_normal = arbitrary_perpendicular(first_tangent);

    for (index, tangent) in tangents.iter().copied().enumerate() {
        if index > 0 {
            let previous_tangent = tangents[index - 1];
            current_normal = rotate_minimizing(previous_tangent, tangent, current_normal);
            current_normal = orthogonalized(current_normal, tangent).normalize_or_zero();
            if current_normal.length_squared() <= 1.0e-8 {
                current_normal = arbitrary_perpendicular(tangent);
            }
        } else {
            current_normal = orthogonalized(current_normal, tangent).normalize_or_zero();
        }
        let binormal = tangent.cross(current_normal).normalize_or_zero();
        current_normal = binormal.cross(tangent).normalize_or_zero();
        normals.push(current_normal);
        binormals.push(binormal);
    }

    (normals, binormals)
}

fn surface_aligned_frames(tangents: &[Vec3], surface_normals: &[Vec3]) -> (Vec<Vec3>, Vec<Vec3>) {
    let mut normals = Vec::with_capacity(tangents.len());
    let mut binormals = Vec::with_capacity(tangents.len());
    for (tangent, surface_normal) in tangents
        .iter()
        .copied()
        .zip(surface_normals.iter().copied())
    {
        let mut normal = orthogonalized(surface_normal, tangent).normalize_or_zero();
        if normal.length_squared() <= 1.0e-8 {
            normal = arbitrary_perpendicular(tangent);
        }
        let binormal = tangent.cross(normal).normalize_or_zero();
        normals.push(normal);
        binormals.push(binormal);
    }
    (normals, binormals)
}

fn darboux_frames(tangents: &[Vec3], surface_normals: &[Vec3]) -> (Vec<Vec3>, Vec<Vec3>) {
    let mut normals = Vec::with_capacity(tangents.len());
    let mut binormals = Vec::with_capacity(tangents.len());
    for (tangent, surface_normal) in tangents
        .iter()
        .copied()
        .zip(surface_normals.iter().copied())
    {
        let mut surface_normal = orthogonalized(surface_normal, tangent).normalize_or_zero();
        if surface_normal.length_squared() <= 1.0e-8 {
            surface_normal = arbitrary_perpendicular(tangent);
        }
        let mut geodesic = surface_normal.cross(tangent).normalize_or_zero();
        if geodesic.length_squared() <= 1.0e-8 {
            geodesic = arbitrary_perpendicular(tangent);
        }
        normals.push(geodesic);
        binormals.push(surface_normal);
    }
    (normals, binormals)
}

fn nearest_surface_normal(mesh: &CombinedSurfaceMesh, point: Vec3) -> Vec3 {
    let mut best_distance = f32::INFINITY;
    let mut best_normal = Vec3::Z;
    for triangle in mesh.indices.chunks_exact(3) {
        let a = mesh.positions[triangle[0] as usize];
        let b = mesh.positions[triangle[1] as usize];
        let c = mesh.positions[triangle[2] as usize];
        let normal = (b - a).cross(c - a).normalize_or_zero();
        if normal.length_squared() <= 1.0e-8 {
            continue;
        }
        let centroid = (a + b + c) / 3.0;
        let distance = centroid.distance_squared(point);
        if distance < best_distance {
            best_distance = distance;
            best_normal = normal;
        }
    }
    best_normal
}

fn rotate_minimizing(previous_tangent: Vec3, tangent: Vec3, normal: Vec3) -> Vec3 {
    if previous_tangent.length_squared() <= 1.0e-8 || tangent.length_squared() <= 1.0e-8 {
        return normal;
    }
    let dot = previous_tangent.dot(tangent).clamp(-1.0, 1.0);
    if dot > 0.9999 {
        return normal;
    }
    if dot < -0.9999 {
        let axis = arbitrary_perpendicular(previous_tangent);
        return glam::Quat::from_axis_angle(axis.normalize_or_zero(), std::f32::consts::PI)
            * normal;
    }
    glam::Quat::from_rotation_arc(previous_tangent, tangent) * normal
}

fn orthogonalized(vector: Vec3, tangent: Vec3) -> Vec3 {
    vector - tangent * vector.dot(tangent)
}

fn arbitrary_perpendicular(tangent: Vec3) -> Vec3 {
    let reference = if tangent.z.abs() < 0.9 {
        Vec3::Z
    } else {
        Vec3::Y
    };
    tangent.cross(reference).normalize_or_zero()
}

fn make_point_cloud_statistics_output(plot: &PlotSpec) -> Result<AnalysisOutput, AnalysisError> {
    let groups = sample_data_groups(plot)?;
    let samples = flatten_sample_groups(&groups);
    if samples.is_empty() {
        return Err(AnalysisError::invalid(
            "Point-cloud statistics require at least one sample.",
        ));
    }

    let centroid = samples.iter().copied().sum::<Vec3>() / samples.len() as f32;
    let (bbox_min, bbox_max) = bounds_for_points(&samples);
    let extent = bbox_max - bbox_min;
    let covariance = covariance_matrix(&samples, centroid);
    let variance = Vec3::new(
        covariance.x_axis.x,
        covariance.y_axis.y,
        covariance.z_axis.z,
    );
    let principal_components = principal_components(covariance);

    let reports = vec![AnalysisReport {
        title: format!("Point Statistics {}", plot.name),
        values: vec![
            ("Sample Count".to_string(), samples.len().to_string()),
            ("Sequence Count".to_string(), groups.len().to_string()),
            ("Centroid".to_string(), format_vec3(centroid)),
            ("Bounds Min".to_string(), format_vec3(bbox_min)),
            ("Bounds Max".to_string(), format_vec3(bbox_max)),
            ("Extent".to_string(), format_vec3(extent)),
            (
                "Variance".to_string(),
                format!("{:.5}, {:.5}, {:.5}", variance.x, variance.y, variance.z),
            ),
        ],
    }];
    let tables = vec![
        sample_positions_table(&groups),
        AnalysisTable {
            title: "Covariance Matrix".to_string(),
            columns: vec![
                "axis".to_string(),
                "x".to_string(),
                "y".to_string(),
                "z".to_string(),
            ],
            rows: vec![
                vec![
                    "covariance row 1".to_string(),
                    format_float(covariance.x_axis.x),
                    format_float(covariance.y_axis.x),
                    format_float(covariance.z_axis.x),
                ],
                vec![
                    "covariance row 2".to_string(),
                    format_float(covariance.x_axis.y),
                    format_float(covariance.y_axis.y),
                    format_float(covariance.z_axis.y),
                ],
                vec![
                    "covariance row 3".to_string(),
                    format_float(covariance.x_axis.z),
                    format_float(covariance.y_axis.z),
                    format_float(covariance.z_axis.z),
                ],
            ],
        },
        AnalysisTable {
            title: "Principal Components".to_string(),
            columns: vec![
                "component".to_string(),
                "eigenvalue".to_string(),
                "direction".to_string(),
            ],
            rows: principal_components
                .iter()
                .enumerate()
                .map(|(index, (eigenvalue, direction))| {
                    vec![
                        format!("PC{}", index + 1),
                        format_float(*eigenvalue),
                        format_vec3(*direction),
                    ]
                })
                .collect(),
        },
    ];
    let plots = vec![
        PlotSpec {
            name: format!("Centroid {}", plot.name),
            visible: true,
            domain: plot.domain.clone(),
            resolution: plot.resolution,
            style: PlotStyle {
                colour_mode: ColourMode::Solid([1.0, 0.86, 0.3, 1.0]),
                point_size: 9.0,
                ..PlotStyle::default()
            },
            definition: crate::PlotDefinition::PointAnnotations {
                points: vec![PointAnnotation {
                    position: centroid.to_array(),
                    label: "Centroid".to_string(),
                }],
                show_labels: true,
            },
        },
        PlotSpec {
            name: format!("PCA Axes {}", plot.name),
            visible: true,
            domain: plot.domain.clone(),
            resolution: plot.resolution,
            style: PlotStyle {
                colour_mode: ColourMode::Colormap {
                    colormap: ColormapSource::Builtin(BuiltinColourmap::RdBu),
                    scalar_range: Some((-1.0, 1.0)),
                },
                glyph_scale: 1.0,
                shading: crate::ShadingMode::Unlit,
                ..PlotStyle::default()
            },
            definition: crate::PlotDefinition::ArrowAnnotations {
                arrows: principal_components
                    .iter()
                    .enumerate()
                    .map(|(index, (eigenvalue, direction))| ArrowAnnotation {
                        origin: centroid.to_array(),
                        vector: (*direction * (eigenvalue.max(0.0).sqrt() * 2.5)).to_array(),
                        label: format!("PC{}", index + 1),
                    })
                    .collect(),
                show_labels: true,
            },
        },
    ];
    Ok(AnalysisOutput::Composite {
        plots,
        reports,
        tables,
        diagnostics: Vec::new(),
        frame_fields: Vec::new(),
        provenance: AnalysisProvenance {
            kind: AnalysisKind::PointCloudStatistics,
            source_plots: vec![plot.name.clone()],
            parameters: Vec::new(),
            notes: vec!["Computed from sampled point positions.".to_string()],
        },
    })
}

fn scatter_cloud_points() -> Vec<[f32; 3]> {
    (0..200)
        .map(|i| {
            glam::Vec3::new(
                (i as f32 * 0.37).sin() * 5.0,
                (i as f32 * 0.73).cos() * 5.0,
                (i as f32 * 0.11).sin() * 5.0,
            )
            .to_array()
        })
        .collect()
}

fn make_data_quality_output(plot: &PlotSpec) -> Result<AnalysisOutput, AnalysisError> {
    let groups = sample_data_groups(plot)?;
    let indexed = flatten_indexed_sample_groups(&groups);
    if indexed.len() < 2 {
        return Err(AnalysisError::invalid(
            "Data quality checks require at least two samples.",
        ));
    }
    let samples = indexed
        .iter()
        .map(|sample| sample.position)
        .collect::<Vec<_>>();
    let centroid = samples.iter().copied().sum::<Vec3>() / samples.len() as f32;
    let (bbox_min, bbox_max) = bounds_for_points(&samples);
    let diag = bbox_max.distance(bbox_min).max(1.0e-4);
    let duplicate_rows = exact_duplicate_rows(&indexed);

    let heavy_checks = samples.len() <= 4_000;
    let nearest = if heavy_checks {
        nearest_neighbor_stats(&indexed)
    } else {
        None
    };
    let near_duplicate_epsilon = diag * 1.0e-3;
    let near_duplicate_samples = nearest
        .as_ref()
        .map(|stats| {
            stats
                .nearest_distances
                .iter()
                .enumerate()
                .filter(|(_, distance)| **distance > 0.0 && **distance <= near_duplicate_epsilon)
                .map(|(index, distance)| (indexed[index].clone(), *distance))
                .collect::<Vec<_>>()
        })
        .unwrap_or_default();
    let (distance_mean, distance_std) = mean_std(
        &samples
            .iter()
            .map(|point| point.distance(centroid))
            .collect::<Vec<_>>(),
    );
    let outlier_threshold = distance_mean + distance_std * 2.5;
    let outliers = indexed
        .iter()
        .filter(|sample| sample.position.distance(centroid) > outlier_threshold)
        .cloned()
        .collect::<Vec<_>>();
    let sparse_samples = nearest
        .as_ref()
        .map(|stats| {
            stats
                .nearest_distances
                .iter()
                .enumerate()
                .filter(|(_, distance)| **distance > stats.mean + stats.std * 2.0)
                .map(|(index, _)| indexed[index].clone())
                .collect::<Vec<_>>()
        })
        .unwrap_or_default();
    let monotonicity_rows = groups
        .iter()
        .enumerate()
        .filter(|(_, group)| group.len() >= 2)
        .map(|(index, group)| {
            vec![
                format!("sequence {}", index + 1),
                monotonicity_label(group, 0),
                monotonicity_label(group, 1),
                monotonicity_label(group, 2),
            ]
        })
        .collect::<Vec<_>>();

    let mut diagnostics = Vec::new();
    if !duplicate_rows.is_empty() {
        diagnostics.push(Diagnostic::warning(
            DiagnosticKind::Validation,
            format!(
                "Detected {} exact duplicate sample(s).",
                duplicate_rows.len()
            ),
        ));
    }
    if !near_duplicate_samples.is_empty() {
        diagnostics.push(Diagnostic::warning(
            DiagnosticKind::Validation,
            format!(
                "Detected {} near-duplicate sample(s) within {:.4}.",
                near_duplicate_samples.len(),
                near_duplicate_epsilon
            ),
        ));
    }
    if !outliers.is_empty() {
        diagnostics.push(Diagnostic::warning(
            DiagnosticKind::Validation,
            format!("Detected {} positional outlier(s).", outliers.len()),
        ));
    }
    if !sparse_samples.is_empty() {
        diagnostics.push(Diagnostic::warning(
            DiagnosticKind::Validation,
            format!("Detected {} sparse-region sample(s).", sparse_samples.len()),
        ));
    }
    if !heavy_checks {
        diagnostics.push(Diagnostic::warning(
            DiagnosticKind::Build,
            "Skipped nearest-neighbor heavy checks for datasets above 4000 samples.",
        ));
    }
    diagnostics.extend(
        duplicate_rows
            .iter()
            .take(12)
            .map(|sample| sample_warning(sample, "Exact duplicate sample.".to_string())),
    );
    diagnostics.extend(
        near_duplicate_samples
            .iter()
            .take(12)
            .map(|(sample, distance)| {
                sample_warning(
                    sample,
                    format!(
                        "Near-duplicate sample with nearest-neighbour distance {:.4}.",
                        distance
                    ),
                )
            }),
    );
    diagnostics.extend(sparse_samples.iter().take(12).map(|sample| {
        let nearest_distance = nearest
            .as_ref()
            .and_then(|stats| stats.nearest_distances.get(sample.sample_index))
            .copied()
            .unwrap_or_default();
        sample_warning(
            sample,
            format!(
                "Sparse-region sample with nearest-neighbour distance {:.4}.",
                nearest_distance
            ),
        )
    }));
    diagnostics.extend(outliers.iter().take(12).map(|sample| {
        sample_warning(
            sample,
            format!(
                "Outlier distance {:.4} exceeds threshold {:.4}.",
                sample.position.distance(centroid),
                outlier_threshold
            ),
        )
    }));

    let mut reports = vec![AnalysisReport {
        title: format!("Data Quality {}", plot.name),
        values: vec![
            ("Sample Count".to_string(), samples.len().to_string()),
            (
                "Exact Duplicates".to_string(),
                duplicate_rows.len().to_string(),
            ),
            (
                "Near Duplicates".to_string(),
                near_duplicate_samples.len().to_string(),
            ),
            ("Outliers".to_string(), outliers.len().to_string()),
            (
                "Sparse Samples".to_string(),
                sparse_samples.len().to_string(),
            ),
        ],
    }];
    if let Some(stats) = &nearest {
        reports.push(AnalysisReport {
            title: "Spacing Diagnostics".to_string(),
            values: vec![
                ("Nearest Min".to_string(), format_float(stats.min)),
                ("Nearest Mean".to_string(), format_float(stats.mean)),
                ("Nearest Median".to_string(), format_float(stats.median)),
                ("Nearest Max".to_string(), format_float(stats.max)),
                ("Nearest Std".to_string(), format_float(stats.std)),
            ],
        });
    }

    let mut tables = vec![sample_positions_table(&groups)];
    if !duplicate_rows.is_empty() {
        tables.push(flagged_samples_table(
            "Exact Duplicates",
            "duplicate_count",
            duplicate_rows
                .iter()
                .map(|sample| {
                    (
                        sample,
                        sample_position_occurrences(&indexed, sample.position),
                    )
                })
                .take(250),
        ));
    }
    if !near_duplicate_samples.is_empty() {
        tables.push(flagged_samples_table(
            "Near Duplicates",
            "nearest_distance",
            near_duplicate_samples
                .iter()
                .map(|(sample, distance)| (sample, *distance))
                .take(250),
        ));
    }
    if !sparse_samples.is_empty() {
        tables.push(flagged_samples_table(
            "Sparse Samples",
            "nearest_distance",
            sparse_samples
                .iter()
                .map(|sample| {
                    let nearest_distance = nearest
                        .as_ref()
                        .and_then(|stats| stats.nearest_distances.get(sample.sample_index))
                        .copied()
                        .unwrap_or_default();
                    (sample, nearest_distance)
                })
                .take(250),
        ));
    }
    if !outliers.is_empty() {
        tables.push(flagged_samples_table(
            "Outliers",
            "distance_from_centroid",
            outliers
                .iter()
                .map(|sample| (sample, sample.position.distance(centroid)))
                .take(250),
        ));
    }
    if !monotonicity_rows.is_empty() {
        tables.push(AnalysisTable {
            title: "Monotonicity by Sequence".to_string(),
            columns: vec![
                "sequence".to_string(),
                "x monotonic".to_string(),
                "y monotonic".to_string(),
                "z monotonic".to_string(),
            ],
            rows: monotonicity_rows,
        });
    }

    let mut plots = Vec::new();
    if !outliers.is_empty() {
        plots.push(PlotSpec {
            name: format!("Outliers {}", plot.name),
            visible: true,
            domain: plot.domain.clone(),
            resolution: plot.resolution,
            style: PlotStyle {
                colour_mode: ColourMode::Solid([1.0, 0.35, 0.35, 1.0]),
                point_size: 8.0,
                ..PlotStyle::default()
            },
            definition: crate::PlotDefinition::PointAnnotations {
                points: outliers
                    .iter()
                    .map(|sample| PointAnnotation {
                        position: sample.position.to_array(),
                        label: String::new(),
                    })
                    .collect(),
                show_labels: false,
            },
        });
    }

    Ok(AnalysisOutput::Composite {
        plots,
        reports,
        tables,
        diagnostics,
        frame_fields: Vec::new(),
        provenance: AnalysisProvenance {
            kind: AnalysisKind::DataQualityChecks,
            source_plots: vec![plot.name.clone()],
            parameters: Vec::new(),
            notes: vec!["Ordered-sequence checks operate per sampled sequence.".to_string()],
        },
    })
}

struct CombinedSurfaceMesh {
    positions: Vec<Vec3>,
    indices: Vec<u32>,
    angle_distortion: Vec<f32>,
    area_distortion: Vec<f32>,
}

impl CombinedSurfaceMesh {
    fn from_meshes(meshes: &[&MeshData]) -> Result<Self, AnalysisError> {
        if meshes.is_empty() {
            return Err(AnalysisError::invalid(
                "Surface analysis requires at least one sampled surface mesh.",
            ));
        }

        let mut positions = Vec::new();
        let mut indices = Vec::new();
        let mut angle_distortion = Vec::new();
        let mut area_distortion = Vec::new();

        for mesh in meshes {
            if mesh.positions.is_empty() || mesh.indices.len() < 3 {
                continue;
            }
            let base = positions.len() as u32;
            positions.extend(mesh.positions.iter().copied().map(Vec3::from_array));
            indices.extend(mesh.indices.iter().map(|index| base + *index));

            if let Some(AttributeData::Face(values)) = mesh.attributes.get("angle_distortion") {
                angle_distortion.extend(values.iter().copied());
            }
            if let Some(AttributeData::Face(values)) = mesh.attributes.get("area_distortion") {
                area_distortion.extend(values.iter().copied());
            }
        }

        if positions.is_empty() || indices.len() < 3 {
            return Err(AnalysisError::invalid(
                "Surface analysis requires cached triangle mesh geometry.",
            ));
        }

        Ok(Self {
            positions,
            indices,
            angle_distortion,
            area_distortion,
        })
    }
}

struct SurfaceCurvatureSummary {
    boundary_vertices: Vec<bool>,
    mean_curvature: Vec<f32>,
    gaussian_curvature: Vec<f32>,
    principal_max: Vec<f32>,
    principal_min: Vec<f32>,
    neighbors: Vec<Vec<usize>>,
}

#[derive(Clone)]
struct SurfaceExtremum {
    vertex_index: usize,
    position: Vec3,
    value: f32,
    label: String,
}

#[derive(Clone)]
struct MeshQualityFaceRow {
    face_index: usize,
    centroid: Vec3,
    area: f32,
    min_angle_deg: f32,
    max_angle_deg: f32,
    aspect_ratio: f32,
    angle_distortion: Option<f32>,
    area_distortion: Option<f32>,
    score: f32,
}

fn make_surface_normals_output(
    source: &PlotSpec,
    mesh: &CombinedSurfaceMesh,
    max_samples: usize,
    vector_scale: f32,
) -> Result<AnalysisOutput, AnalysisError> {
    let normals = compute_vertex_normals(&mesh.positions, &mesh.indices);
    let (bbox_min, bbox_max) = bounds_for_points(&mesh.positions);
    let diagonal = bbox_max.distance(bbox_min);
    let scale = (diagonal * 0.06 * vector_scale.max(0.05)).max(0.01);
    let sampled = sampled_vertex_indices(mesh.positions.len(), max_samples.max(1));
    let plot = PlotSpec {
        name: format!("Surface Normals {}", source.name),
        visible: true,
        domain: source.domain.clone(),
        resolution: source.resolution,
        style: PlotStyle {
            colour_mode: ColourMode::Colormap {
                colormap: ColormapSource::Builtin(BuiltinColourmap::RdBu),
                scalar_range: Some((-1.0, 1.0)),
            },
            glyph_scale: 1.0,
            shading: crate::ShadingMode::Unlit,
            ..PlotStyle::default()
        },
        definition: crate::PlotDefinition::ArrowAnnotations {
            arrows: sampled
                .into_iter()
                .map(|index| ArrowAnnotation {
                    origin: mesh.positions[index].to_array(),
                    vector: (normals[index] * scale).to_array(),
                    label: String::new(),
                })
                .collect(),
            show_labels: false,
        },
    };
    Ok(AnalysisOutput::DerivedPlots {
        plots: vec![plot],
        provenance: AnalysisProvenance {
            kind: AnalysisKind::SurfaceNormals,
            source_plots: vec![source.name.clone()],
            parameters: vec![
                ("max_samples".to_string(), max_samples.to_string()),
                ("vector_scale".to_string(), format_float(vector_scale)),
            ],
            notes: vec!["Normals sampled from cached surface vertices.".to_string()],
        },
    })
}

fn make_surface_area_output(
    source: &PlotSpec,
    mesh: &CombinedSurfaceMesh,
) -> Result<AnalysisOutput, AnalysisError> {
    let face_areas = triangle_areas(&mesh.positions, &mesh.indices);
    let total_area = face_areas.iter().sum::<f32>();
    let (min_area, mean_area, max_area) = min_mean_max(&face_areas);
    Ok(AnalysisOutput::Report {
        report: AnalysisReport {
            title: format!("Surface Area {}", source.name),
            values: vec![
                ("Vertex Count".to_string(), mesh.positions.len().to_string()),
                (
                    "Triangle Count".to_string(),
                    (mesh.indices.len() / 3).to_string(),
                ),
                ("Surface Area".to_string(), format_float(total_area)),
                ("Triangle Area Min".to_string(), format_float(min_area)),
                ("Triangle Area Mean".to_string(), format_float(mean_area)),
                ("Triangle Area Max".to_string(), format_float(max_area)),
            ],
        },
        provenance: AnalysisProvenance {
            kind: AnalysisKind::SurfaceArea,
            source_plots: vec![source.name.clone()],
            parameters: Vec::new(),
            notes: vec!["Surface area is estimated from cached triangle geometry.".to_string()],
        },
    })
}

fn make_surface_curvature_output(
    source: &PlotSpec,
    mesh: &CombinedSurfaceMesh,
    params: &HashMap<String, String>,
) -> Result<AnalysisOutput, AnalysisError> {
    let summary = estimate_surface_curvatures(mesh);
    let quantity = params
        .get("quantity")
        .map(String::as_str)
        .unwrap_or("mean_curvature");
    let (quantity_label, values) = match quantity {
        "gaussian_curvature" => ("Gaussian Curvature", summary.gaussian_curvature.clone()),
        "k_max" => ("Principal Max Curvature", summary.principal_max.clone()),
        "k_min" => ("Principal Min Curvature", summary.principal_min.clone()),
        _ => ("Mean Curvature", summary.mean_curvature.clone()),
    };
    let show_extrema = params
        .get("show_extrema")
        .is_none_or(|value| matches!(value.as_str(), "1" | "true" | "yes"));
    let ridge_extrema = local_principal_maxima(mesh, &summary, 8);
    let valley_extrema = local_principal_minima(mesh, &summary, 8);
    let gaussian_peak =
        top_vertex_extrema(mesh, &summary.gaussian_curvature, 2, "Gaussian Peak", true);
    let gaussian_pit =
        top_vertex_extrema(mesh, &summary.gaussian_curvature, 2, "Gaussian Pit", false);

    let mut marker_points = Vec::new();
    marker_points.extend(ridge_extrema.iter().cloned());
    marker_points.extend(valley_extrema.iter().cloned());
    marker_points.extend(gaussian_peak.iter().cloned());
    marker_points.extend(gaussian_pit.iter().cloned());

    let mean_abs = summary
        .mean_curvature
        .iter()
        .map(|value| value.abs())
        .collect::<Vec<_>>();
    let gaussian_abs = summary
        .gaussian_curvature
        .iter()
        .map(|value| value.abs())
        .collect::<Vec<_>>();
    let principal_abs = summary
        .principal_max
        .iter()
        .zip(summary.principal_min.iter())
        .map(|(k1, k2)| k1.abs().max(k2.abs()))
        .collect::<Vec<_>>();

    let (range_min, _, range_max) = min_mean_max(&values);
    let mut plots = vec![PlotSpec {
        name: format!("{quantity_label} {}", source.name),
        visible: true,
        domain: source.domain.clone(),
        resolution: source.resolution,
        style: PlotStyle {
            colour_mode: ColourMode::Colormap {
                colormap: ColormapSource::Builtin(BuiltinColourmap::Coolwarm),
                scalar_range: Some((range_min, range_max)),
            },
            opacity: 0.92,
            two_sided: true,
            ..PlotStyle::default()
        },
        definition: crate::PlotDefinition::DerivedSurfaceMesh {
            positions: mesh
                .positions
                .iter()
                .map(|position| position.to_array())
                .collect(),
            indices: mesh.indices.clone(),
            values: values.clone(),
            value_name: quantity.to_string(),
        },
    }];

    if show_extrema && !marker_points.is_empty() {
        plots.push(PlotSpec {
            name: format!("Curvature Markers {}", source.name),
            visible: true,
            domain: source.domain.clone(),
            resolution: source.resolution,
            style: PlotStyle {
                colour_mode: ColourMode::Solid([1.0, 0.6, 0.18, 1.0]),
                point_size: 8.0,
                ..PlotStyle::default()
            },
            definition: crate::PlotDefinition::PointAnnotations {
                points: marker_points
                    .iter()
                    .map(|marker| PointAnnotation {
                        position: marker.position.to_array(),
                        label: String::new(),
                    })
                    .collect(),
                show_labels: false,
            },
        });
    }

    let reports = vec![AnalysisReport {
        title: format!("Surface Curvature {}", source.name),
        values: vec![
            ("Vertex Count".to_string(), mesh.positions.len().to_string()),
            (
                "Triangle Count".to_string(),
                (mesh.indices.len() / 3).to_string(),
            ),
            (
                "Boundary Vertices".to_string(),
                summary
                    .boundary_vertices
                    .iter()
                    .filter(|value| **value)
                    .count()
                    .to_string(),
            ),
            (
                "Mean |H|".to_string(),
                format_float(mean_abs.iter().sum::<f32>() / mean_abs.len().max(1) as f32),
            ),
            (
                "Max |H|".to_string(),
                format_float(mean_abs.iter().copied().fold(0.0, f32::max)),
            ),
            (
                "Mean |K|".to_string(),
                format_float(gaussian_abs.iter().sum::<f32>() / gaussian_abs.len().max(1) as f32),
            ),
            (
                "Max |K|".to_string(),
                format_float(gaussian_abs.iter().copied().fold(0.0, f32::max)),
            ),
            (
                "Max |Principal|".to_string(),
                format_float(principal_abs.iter().copied().fold(0.0, f32::max)),
            ),
        ],
    }];

    let mut tables = vec![AnalysisTable {
        title: "Vertex Curvature Samples".to_string(),
        columns: vec![
            "vertex".to_string(),
            "x".to_string(),
            "y".to_string(),
            "z".to_string(),
            "mean_curvature".to_string(),
            "gaussian_curvature".to_string(),
            "k_max".to_string(),
            "k_min".to_string(),
            "boundary".to_string(),
        ],
        rows: mesh
            .positions
            .iter()
            .enumerate()
            .map(|(index, position)| {
                vec![
                    (index + 1).to_string(),
                    format_float(position.x),
                    format_float(position.y),
                    format_float(position.z),
                    format_float(summary.mean_curvature[index]),
                    format_float(summary.gaussian_curvature[index]),
                    format_float(summary.principal_max[index]),
                    format_float(summary.principal_min[index]),
                    summary.boundary_vertices[index].to_string(),
                ]
            })
            .collect(),
    }];

    let mut extrema_rows = Vec::new();
    extrema_rows.extend(ridge_extrema.iter().cloned());
    extrema_rows.extend(valley_extrema.iter().cloned());
    extrema_rows.extend(gaussian_peak.iter().cloned());
    extrema_rows.extend(gaussian_pit.iter().cloned());
    if !extrema_rows.is_empty() {
        tables.push(AnalysisTable {
            title: "Curvature Extrema".to_string(),
            columns: vec![
                "vertex".to_string(),
                "label".to_string(),
                "value".to_string(),
                "x".to_string(),
                "y".to_string(),
                "z".to_string(),
            ],
            rows: extrema_rows
                .iter()
                .map(|row| {
                    vec![
                        (row.vertex_index + 1).to_string(),
                        row.label.clone(),
                        format_float(row.value),
                        format_float(row.position.x),
                        format_float(row.position.y),
                        format_float(row.position.z),
                    ]
                })
                .collect(),
        });
    }

    Ok(AnalysisOutput::Composite {
        plots,
        reports,
        tables,
        diagnostics: Vec::new(),
        frame_fields: Vec::new(),
        provenance: AnalysisProvenance {
            kind: AnalysisKind::SurfaceCurvature,
            source_plots: vec![source.name.clone()],
            parameters: vec![
                ("quantity".to_string(), quantity.to_string()),
                ("show_extrema".to_string(), show_extrema.to_string()),
            ],
            notes: vec![
                "Curvature estimates use discrete triangle-mesh angle-defect and cotangent-weight formulas."
                    .to_string(),
            ],
        },
    })
}

struct CurveSurfaceSample {
    group_index: usize,
    sample_index: usize,
    position: Vec3,
    projected: Vec3,
    distance: f32,
}

fn make_curve_surface_measurement_output(
    curve_source: &PlotSpec,
    surface_source_name: &str,
    mesh: &CombinedSurfaceMesh,
    max_samples: usize,
    vector_scale: f32,
) -> Result<AnalysisOutput, AnalysisError> {
    let groups = curve_sample_groups(curve_source)?;
    let sampled_groups = groups
        .iter()
        .map(|group| sample_curve_group_for_measurement(group, max_samples))
        .filter(|group| group.len() >= 2)
        .collect::<Vec<_>>();
    if sampled_groups.is_empty() {
        return Err(AnalysisError::invalid(
            "Curve-surface measurement requires at least two sampled curve points.",
        ));
    }

    let mut rows = Vec::new();
    let mut projected_groups = Vec::new();
    let mut deviation_arrows = Vec::new();
    let mut original_length = 0.0_f32;
    let mut projected_length = 0.0_f32;

    let (bbox_min, bbox_max) = bounds_for_points(&mesh.positions);
    let diagonal = bbox_max.distance(bbox_min).max(1.0);
    let arrow_stride = (sampled_groups.iter().map(Vec::len).sum::<usize>() / 64).max(1);

    for (group_index, group) in sampled_groups.iter().enumerate() {
        let mut projected_group = Vec::with_capacity(group.len());
        let mut previous_curve = None;
        let mut previous_projected = None;
        for (sample_index, point) in group.iter().copied().enumerate() {
            let projected = closest_point_on_surface(mesh, point);
            let distance = point.distance(projected);
            if let Some(previous) = previous_curve {
                original_length += point.distance(previous);
            }
            if let Some(previous) = previous_projected {
                projected_length += projected.distance(previous);
            }
            previous_curve = Some(point);
            previous_projected = Some(projected);
            projected_group.push(projected.to_array());
            rows.push(CurveSurfaceSample {
                group_index,
                sample_index,
                position: point,
                projected,
                distance,
            });
            if sample_index % arrow_stride == 0 && distance > diagonal * 1.0e-5 {
                deviation_arrows.push(ArrowAnnotation {
                    origin: projected.to_array(),
                    vector: ((point - projected) * vector_scale).to_array(),
                    label: String::new(),
                });
            }
        }
        projected_groups.push(projected_group);
    }

    let distances = rows.iter().map(|row| row.distance).collect::<Vec<_>>();
    let (min_distance, mean_distance, max_distance) = min_mean_max(&distances);
    let rms_distance =
        (distances.iter().map(|value| value * value).sum::<f32>() / distances.len() as f32).sqrt();

    let mut plots = vec![PlotSpec {
        name: format!("Projected {} on {}", curve_source.name, surface_source_name),
        visible: true,
        domain: curve_source.domain.clone(),
        resolution: curve_source.resolution,
        style: PlotStyle {
            colour_mode: ColourMode::Solid([0.2, 0.9, 0.95, 1.0]),
            line_width: 3.0,
            ..PlotStyle::default()
        },
        definition: crate::PlotDefinition::DerivedPolylineGroups {
            groups: projected_groups,
        },
    }];
    if !deviation_arrows.is_empty() {
        plots.push(PlotSpec {
            name: format!("Deviation {} to {}", curve_source.name, surface_source_name),
            visible: true,
            domain: curve_source.domain.clone(),
            resolution: curve_source.resolution,
            style: PlotStyle {
                colour_mode: ColourMode::Solid([1.0, 0.35, 0.25, 1.0]),
                glyph_scale: 1.0,
                shading: crate::ShadingMode::Unlit,
                ..PlotStyle::default()
            },
            definition: crate::PlotDefinition::ArrowAnnotations {
                arrows: deviation_arrows,
                show_labels: false,
            },
        });
    }

    Ok(AnalysisOutput::Composite {
        plots,
        reports: vec![AnalysisReport {
            title: format!("Curve-Surface Measurement {}", curve_source.name),
            values: vec![
                ("Target Surface".to_string(), surface_source_name.to_string()),
                ("Sample Count".to_string(), rows.len().to_string()),
                ("Curve Length".to_string(), format_float(original_length)),
                ("Projected Length".to_string(), format_float(projected_length)),
                (
                    "Length Difference".to_string(),
                    format_float(projected_length - original_length),
                ),
                ("Deviation Min".to_string(), format_float(min_distance)),
                ("Deviation Mean".to_string(), format_float(mean_distance)),
                ("Deviation RMS".to_string(), format_float(rms_distance)),
                ("Deviation Max".to_string(), format_float(max_distance)),
            ],
        }],
        tables: vec![AnalysisTable {
            title: "Curve Surface Samples".to_string(),
            columns: vec![
                "group".to_string(),
                "sample".to_string(),
                "x".to_string(),
                "y".to_string(),
                "z".to_string(),
                "projected_x".to_string(),
                "projected_y".to_string(),
                "projected_z".to_string(),
                "deviation".to_string(),
            ],
            rows: rows
                .iter()
                .map(|row| {
                    vec![
                        (row.group_index + 1).to_string(),
                        (row.sample_index + 1).to_string(),
                        format_float(row.position.x),
                        format_float(row.position.y),
                        format_float(row.position.z),
                        format_float(row.projected.x),
                        format_float(row.projected.y),
                        format_float(row.projected.z),
                        format_float(row.distance),
                    ]
                })
                .collect(),
        }],
        diagnostics: Vec::new(),
        frame_fields: Vec::new(),
        provenance: AnalysisProvenance {
            kind: AnalysisKind::CurveSurfaceMeasurement,
            source_plots: vec![curve_source.name.clone(), surface_source_name.to_string()],
            parameters: vec![
                ("max_samples".to_string(), max_samples.to_string()),
                ("vector_scale".to_string(), format_float(vector_scale)),
            ],
            notes: vec![
                "Curve points are projected to the nearest point on cached target-surface triangles."
                    .to_string(),
            ],
        },
    })
}

fn sample_curve_group_for_measurement(group: &[[f32; 3]], max_samples: usize) -> Vec<Vec3> {
    if group.len() <= max_samples {
        return group.iter().copied().map(Vec3::from_array).collect();
    }
    sampled_vertex_indices(group.len(), max_samples)
        .into_iter()
        .map(|index| Vec3::from_array(group[index]))
        .collect()
}

fn closest_point_on_surface(mesh: &CombinedSurfaceMesh, point: Vec3) -> Vec3 {
    let mut closest = mesh.positions.first().copied().unwrap_or(Vec3::ZERO);
    let mut best_distance_sq = point.distance_squared(closest);
    for tri in mesh.indices.chunks_exact(3) {
        let a = mesh.positions[tri[0] as usize];
        let b = mesh.positions[tri[1] as usize];
        let c = mesh.positions[tri[2] as usize];
        let candidate = closest_point_on_triangle(point, a, b, c);
        let distance_sq = point.distance_squared(candidate);
        if distance_sq < best_distance_sq {
            best_distance_sq = distance_sq;
            closest = candidate;
        }
    }
    closest
}

fn closest_point_on_triangle(point: Vec3, a: Vec3, b: Vec3, c: Vec3) -> Vec3 {
    let ab = b - a;
    let ac = c - a;
    let ap = point - a;
    let d1 = ab.dot(ap);
    let d2 = ac.dot(ap);
    if d1 <= 0.0 && d2 <= 0.0 {
        return a;
    }

    let bp = point - b;
    let d3 = ab.dot(bp);
    let d4 = ac.dot(bp);
    if d3 >= 0.0 && d4 <= d3 {
        return b;
    }

    let vc = d1 * d4 - d3 * d2;
    if vc <= 0.0 && d1 >= 0.0 && d3 <= 0.0 {
        let v = d1 / (d1 - d3);
        return a + ab * v;
    }

    let cp = point - c;
    let d5 = ab.dot(cp);
    let d6 = ac.dot(cp);
    if d6 >= 0.0 && d5 <= d6 {
        return c;
    }

    let vb = d5 * d2 - d1 * d6;
    if vb <= 0.0 && d2 >= 0.0 && d6 <= 0.0 {
        let w = d2 / (d2 - d6);
        return a + ac * w;
    }

    let va = d3 * d6 - d5 * d4;
    if va <= 0.0 && (d4 - d3) >= 0.0 && (d5 - d6) >= 0.0 {
        let w = (d4 - d3) / ((d4 - d3) + (d5 - d6));
        return b + (c - b) * w;
    }

    let denom = 1.0 / (va + vb + vc);
    let v = vb * denom;
    let w = vc * denom;
    a + ab * v + ac * w
}

fn make_surface_mesh_quality_output(
    source: &PlotSpec,
    mesh: &CombinedSurfaceMesh,
) -> Result<AnalysisOutput, AnalysisError> {
    let faces = mesh_quality_rows(mesh);
    let poor_faces = faces
        .iter()
        .filter(|row| row.score > 0.0)
        .cloned()
        .collect::<Vec<_>>();
    let area_values = faces.iter().map(|row| row.area).collect::<Vec<_>>();
    let aspect_values = faces.iter().map(|row| row.aspect_ratio).collect::<Vec<_>>();
    let min_angles = faces
        .iter()
        .map(|row| row.min_angle_deg)
        .collect::<Vec<_>>();
    let angle_distortion_values = faces
        .iter()
        .filter_map(|row| row.angle_distortion)
        .collect::<Vec<_>>();
    let area_distortion_values = faces
        .iter()
        .filter_map(|row| row.area_distortion)
        .collect::<Vec<_>>();

    let mut diagnostics = Vec::new();
    if !poor_faces.is_empty() {
        diagnostics.push(Diagnostic::warning(
            DiagnosticKind::Validation,
            format!(
                "Detected {} mesh-quality warning face(s).",
                poor_faces.len()
            ),
        ));
        diagnostics.extend(poor_faces.iter().take(12).map(|row| {
            Diagnostic::warning(
                DiagnosticKind::Validation,
                format!(
                    "Face quality warning: aspect {:.3}, min angle {:.2}deg.",
                    row.aspect_ratio, row.min_angle_deg
                ),
            )
            .with_location(
                DiagnosticLocation::new()
                    .with_component("face")
                    .with_row(row.face_index + 1),
            )
        }));
    }

    let reports = vec![AnalysisReport {
        title: format!("Surface Mesh Quality {}", source.name),
        values: vec![
            ("Triangle Count".to_string(), faces.len().to_string()),
            ("Warning Faces".to_string(), poor_faces.len().to_string()),
            (
                "Triangle Area Min".to_string(),
                format_float(area_values.iter().copied().fold(f32::INFINITY, f32::min)),
            ),
            (
                "Triangle Area Mean".to_string(),
                format_float(area_values.iter().sum::<f32>() / area_values.len().max(1) as f32),
            ),
            (
                "Aspect Ratio Mean".to_string(),
                format_float(aspect_values.iter().sum::<f32>() / aspect_values.len().max(1) as f32),
            ),
            (
                "Aspect Ratio Max".to_string(),
                format_float(aspect_values.iter().copied().fold(0.0, f32::max)),
            ),
            (
                "Min Angle".to_string(),
                format_float(min_angles.iter().copied().fold(f32::INFINITY, f32::min)),
            ),
            (
                "Mean Angle Distortion".to_string(),
                format_optional_mean(&angle_distortion_values),
            ),
            (
                "Mean Area Distortion".to_string(),
                format_optional_mean(&area_distortion_values),
            ),
        ],
    }];

    let mut tables = vec![AnalysisTable {
        title: "Mesh Face Statistics".to_string(),
        columns: vec![
            "face".to_string(),
            "cx".to_string(),
            "cy".to_string(),
            "cz".to_string(),
            "area".to_string(),
            "min_angle_deg".to_string(),
            "max_angle_deg".to_string(),
            "aspect_ratio".to_string(),
            "angle_distortion".to_string(),
            "area_distortion".to_string(),
            "warning_score".to_string(),
        ],
        rows: faces
            .iter()
            .map(|row| {
                vec![
                    (row.face_index + 1).to_string(),
                    format_float(row.centroid.x),
                    format_float(row.centroid.y),
                    format_float(row.centroid.z),
                    format_float(row.area),
                    format_float(row.min_angle_deg),
                    format_float(row.max_angle_deg),
                    format_float(row.aspect_ratio),
                    row.angle_distortion
                        .map(format_float)
                        .unwrap_or_else(|| "n/a".to_string()),
                    row.area_distortion
                        .map(format_float)
                        .unwrap_or_else(|| "n/a".to_string()),
                    format_float(row.score),
                ]
            })
            .collect(),
    }];

    if !poor_faces.is_empty() {
        tables.push(AnalysisTable {
            title: "Mesh Quality Alerts".to_string(),
            columns: vec![
                "face".to_string(),
                "area".to_string(),
                "min_angle_deg".to_string(),
                "aspect_ratio".to_string(),
                "warning_score".to_string(),
            ],
            rows: poor_faces
                .iter()
                .take(250)
                .map(|row| {
                    vec![
                        (row.face_index + 1).to_string(),
                        format_float(row.area),
                        format_float(row.min_angle_deg),
                        format_float(row.aspect_ratio),
                        format_float(row.score),
                    ]
                })
                .collect(),
        });
    }

    let plots = if poor_faces.is_empty() {
        Vec::new()
    } else {
        vec![PlotSpec {
            name: format!("Mesh Quality Alerts {}", source.name),
            visible: true,
            domain: source.domain.clone(),
            resolution: source.resolution,
            style: PlotStyle {
                colour_mode: ColourMode::Solid([1.0, 0.35, 0.3, 1.0]),
                point_size: 7.5,
                ..PlotStyle::default()
            },
            definition: crate::PlotDefinition::PointAnnotations {
                points: poor_faces
                    .iter()
                    .take(32)
                    .map(|row| PointAnnotation {
                        position: row.centroid.to_array(),
                        label: String::new(),
                    })
                    .collect(),
                show_labels: false,
            },
        }]
    };

    Ok(AnalysisOutput::Composite {
        plots,
        reports,
        tables,
        diagnostics,
        frame_fields: Vec::new(),
        provenance: AnalysisProvenance {
            kind: AnalysisKind::SurfaceMeshQuality,
            source_plots: vec![source.name.clone()],
            parameters: Vec::new(),
            notes: vec![
                "Mesh quality uses per-face area, angle, aspect-ratio, and distortion summaries."
                    .to_string(),
            ],
        },
    })
}

fn estimate_surface_curvatures(mesh: &CombinedSurfaceMesh) -> SurfaceCurvatureSummary {
    let vertex_count = mesh.positions.len();
    let mut vertex_normals_accum = vec![Vec3::ZERO; vertex_count];
    let mut vertex_area = vec![0.0_f32; vertex_count];
    let mut angle_sum = vec![0.0_f32; vertex_count];
    let mut edge_weights: HashMap<(usize, usize), f32> = HashMap::new();
    let mut edge_counts: HashMap<(usize, usize), u32> = HashMap::new();

    for tri in mesh.indices.chunks_exact(3) {
        let i0 = tri[0] as usize;
        let i1 = tri[1] as usize;
        let i2 = tri[2] as usize;
        let p0 = mesh.positions[i0];
        let p1 = mesh.positions[i1];
        let p2 = mesh.positions[i2];

        let face_cross = (p1 - p0).cross(p2 - p0);
        let area = 0.5 * face_cross.length();
        let face_normal = face_cross.normalize_or_zero();
        vertex_normals_accum[i0] += face_normal * area;
        vertex_normals_accum[i1] += face_normal * area;
        vertex_normals_accum[i2] += face_normal * area;
        vertex_area[i0] += area / 3.0;
        vertex_area[i1] += area / 3.0;
        vertex_area[i2] += area / 3.0;

        let angle0 = angle_between_3d_local(p1 - p0, p2 - p0);
        let angle1 = angle_between_3d_local(p2 - p1, p0 - p1);
        let angle2 = angle_between_3d_local(p0 - p2, p1 - p2);
        angle_sum[i0] += angle0;
        angle_sum[i1] += angle1;
        angle_sum[i2] += angle2;

        accumulate_edge_weight(
            &mut edge_weights,
            &mut edge_counts,
            i1,
            i2,
            cotangent(angle0),
        );
        accumulate_edge_weight(
            &mut edge_weights,
            &mut edge_counts,
            i0,
            i2,
            cotangent(angle1),
        );
        accumulate_edge_weight(
            &mut edge_weights,
            &mut edge_counts,
            i0,
            i1,
            cotangent(angle2),
        );
    }

    let vertex_normals = vertex_normals_accum
        .into_iter()
        .map(|normal| normal.normalize_or_zero())
        .collect::<Vec<_>>();
    let mut boundary_vertices = vec![false; vertex_count];
    let mut neighbors = vec![Vec::new(); vertex_count];
    let mut mean_curvature_normal_sum = vec![Vec3::ZERO; vertex_count];

    for (&(a, b), &weight) in &edge_weights {
        neighbors[a].push(b);
        neighbors[b].push(a);
        mean_curvature_normal_sum[a] += weight * (mesh.positions[a] - mesh.positions[b]);
        mean_curvature_normal_sum[b] += weight * (mesh.positions[b] - mesh.positions[a]);
        if edge_counts.get(&(a, b)).copied().unwrap_or(0) == 1 {
            boundary_vertices[a] = true;
            boundary_vertices[b] = true;
        }
    }

    let mut mean_curvature = vec![0.0_f32; vertex_count];
    let mut gaussian_curvature = vec![0.0_f32; vertex_count];
    let mut principal_max = vec![0.0_f32; vertex_count];
    let mut principal_min = vec![0.0_f32; vertex_count];

    for index in 0..vertex_count {
        let area = vertex_area[index].max(1.0e-6);
        let mean_normal = mean_curvature_normal_sum[index] / (2.0 * area);
        let signed_mean = 0.5 * mean_normal.dot(vertex_normals[index]);
        let gaussian = ((if boundary_vertices[index] {
            std::f32::consts::PI
        } else {
            std::f32::consts::TAU
        }) - angle_sum[index])
            / area;
        let discriminant = (signed_mean * signed_mean - gaussian).max(0.0).sqrt();
        mean_curvature[index] = signed_mean;
        gaussian_curvature[index] = gaussian;
        principal_max[index] = signed_mean + discriminant;
        principal_min[index] = signed_mean - discriminant;
    }

    SurfaceCurvatureSummary {
        boundary_vertices,
        mean_curvature,
        gaussian_curvature,
        principal_max,
        principal_min,
        neighbors,
    }
}

fn mesh_quality_rows(mesh: &CombinedSurfaceMesh) -> Vec<MeshQualityFaceRow> {
    mesh.indices
        .chunks_exact(3)
        .enumerate()
        .map(|(face_index, tri)| {
            let p0 = mesh.positions[tri[0] as usize];
            let p1 = mesh.positions[tri[1] as usize];
            let p2 = mesh.positions[tri[2] as usize];
            let edge_lengths = [(p1 - p0).length(), (p2 - p1).length(), (p0 - p2).length()];
            let longest = edge_lengths.iter().copied().fold(0.0, f32::max);
            let shortest = edge_lengths
                .iter()
                .copied()
                .fold(f32::INFINITY, f32::min)
                .max(1.0e-6);
            let aspect_ratio = longest / shortest;
            let angles = triangle_angles_3d_local([p0, p1, p2]).map(f32::to_degrees);
            let min_angle_deg = angles.iter().copied().fold(f32::INFINITY, f32::min);
            let max_angle_deg = angles.iter().copied().fold(0.0, f32::max);
            let area = 0.5 * (p1 - p0).cross(p2 - p0).length();
            let angle_distortion = mesh.angle_distortion.get(face_index).copied();
            let area_distortion = mesh.area_distortion.get(face_index).copied();
            let score = mesh_quality_score(area, min_angle_deg, max_angle_deg, aspect_ratio);
            MeshQualityFaceRow {
                face_index,
                centroid: (p0 + p1 + p2) / 3.0,
                area,
                min_angle_deg,
                max_angle_deg,
                aspect_ratio,
                angle_distortion,
                area_distortion,
                score,
            }
        })
        .collect()
}

fn local_principal_maxima(
    mesh: &CombinedSurfaceMesh,
    summary: &SurfaceCurvatureSummary,
    limit: usize,
) -> Vec<SurfaceExtremum> {
    local_extrema(
        mesh,
        &summary.principal_max,
        &summary.neighbors,
        limit,
        true,
    )
    .into_iter()
    .filter(|extremum| extremum.value > 0.0)
    .map(|mut extremum| {
        extremum.label = "Ridge".to_string();
        extremum
    })
    .collect()
}

fn local_principal_minima(
    mesh: &CombinedSurfaceMesh,
    summary: &SurfaceCurvatureSummary,
    limit: usize,
) -> Vec<SurfaceExtremum> {
    local_extrema(
        mesh,
        &summary.principal_min,
        &summary.neighbors,
        limit,
        false,
    )
    .into_iter()
    .filter(|extremum| extremum.value < 0.0)
    .map(|mut extremum| {
        extremum.label = "Valley".to_string();
        extremum
    })
    .collect()
}

fn local_extrema(
    mesh: &CombinedSurfaceMesh,
    values: &[f32],
    neighbors: &[Vec<usize>],
    limit: usize,
    maxima: bool,
) -> Vec<SurfaceExtremum> {
    let mut extrema = values
        .iter()
        .enumerate()
        .filter_map(|(index, value)| {
            let neighbour_values = neighbors[index]
                .iter()
                .filter_map(|neighbor| values.get(*neighbor))
                .copied()
                .collect::<Vec<_>>();
            if neighbour_values.is_empty() {
                return None;
            }
            let is_extremum = if maxima {
                neighbour_values.iter().all(|other| *value >= *other)
                    && neighbour_values.iter().any(|other| *value > *other)
            } else {
                neighbour_values.iter().all(|other| *value <= *other)
                    && neighbour_values.iter().any(|other| *value < *other)
            };
            is_extremum.then_some(SurfaceExtremum {
                vertex_index: index,
                position: mesh.positions[index],
                value: *value,
                label: String::new(),
            })
        })
        .collect::<Vec<_>>();

    extrema.sort_by(|a, b| {
        b.value
            .abs()
            .partial_cmp(&a.value.abs())
            .unwrap_or(std::cmp::Ordering::Equal)
    });
    extrema.truncate(limit);
    extrema
}

fn top_vertex_extrema(
    mesh: &CombinedSurfaceMesh,
    values: &[f32],
    limit: usize,
    label: &str,
    descending: bool,
) -> Vec<SurfaceExtremum> {
    let mut rows = values
        .iter()
        .enumerate()
        .map(|(index, value)| SurfaceExtremum {
            vertex_index: index,
            position: mesh.positions[index],
            value: *value,
            label: label.to_string(),
        })
        .collect::<Vec<_>>();
    rows.sort_by(|a, b| {
        if descending {
            b.value
                .partial_cmp(&a.value)
                .unwrap_or(std::cmp::Ordering::Equal)
        } else {
            a.value
                .partial_cmp(&b.value)
                .unwrap_or(std::cmp::Ordering::Equal)
        }
    });
    rows.truncate(limit);
    rows
}

fn compute_vertex_normals(positions: &[Vec3], indices: &[u32]) -> Vec<Vec3> {
    let mut normals = vec![Vec3::ZERO; positions.len()];
    for tri in indices.chunks_exact(3) {
        let i0 = tri[0] as usize;
        let i1 = tri[1] as usize;
        let i2 = tri[2] as usize;
        let cross = (positions[i1] - positions[i0]).cross(positions[i2] - positions[i0]);
        normals[i0] += cross;
        normals[i1] += cross;
        normals[i2] += cross;
    }
    normals
        .into_iter()
        .map(|normal| normal.normalize_or_zero())
        .collect()
}

fn triangle_areas(positions: &[Vec3], indices: &[u32]) -> Vec<f32> {
    indices
        .chunks_exact(3)
        .map(|tri| {
            let p0 = positions[tri[0] as usize];
            let p1 = positions[tri[1] as usize];
            let p2 = positions[tri[2] as usize];
            0.5 * (p1 - p0).cross(p2 - p0).length()
        })
        .collect()
}

fn sampled_vertex_indices(vertex_count: usize, max_samples: usize) -> Vec<usize> {
    if vertex_count <= max_samples {
        return (0..vertex_count).collect();
    }
    let step = ((vertex_count as f32 / max_samples as f32).ceil() as usize).max(1);
    (0..vertex_count).step_by(step).collect()
}

fn accumulate_edge_weight(
    edge_weights: &mut HashMap<(usize, usize), f32>,
    edge_counts: &mut HashMap<(usize, usize), u32>,
    a: usize,
    b: usize,
    weight: f32,
) {
    let key = if a < b { (a, b) } else { (b, a) };
    *edge_weights.entry(key).or_insert(0.0) += weight;
    *edge_counts.entry(key).or_insert(0) += 1;
}

fn cotangent(angle: f32) -> f32 {
    let sin = angle.sin().abs().max(1.0e-6);
    angle.cos() / sin
}

fn triangle_angles_3d_local(points: [Vec3; 3]) -> [f32; 3] {
    [
        angle_between_3d_local(points[1] - points[0], points[2] - points[0]),
        angle_between_3d_local(points[0] - points[1], points[2] - points[1]),
        angle_between_3d_local(points[0] - points[2], points[1] - points[2]),
    ]
}

fn angle_between_3d_local(a: Vec3, b: Vec3) -> f32 {
    let denom = a.length() * b.length();
    if denom <= 1.0e-8 {
        0.0
    } else {
        (a.dot(b) / denom).clamp(-1.0, 1.0).acos()
    }
}

fn mesh_quality_score(area: f32, min_angle_deg: f32, max_angle_deg: f32, aspect_ratio: f32) -> f32 {
    let mut score = 0.0;
    if area <= 1.0e-8 {
        score += 10.0;
    }
    if min_angle_deg < 15.0 {
        score += (15.0 - min_angle_deg) / 15.0;
    }
    if max_angle_deg > 150.0 {
        score += (max_angle_deg - 150.0) / 30.0;
    }
    if aspect_ratio > 4.0 {
        score += (aspect_ratio - 4.0) / 2.0;
    }
    score
}

fn min_mean_max(values: &[f32]) -> (f32, f32, f32) {
    if values.is_empty() {
        return (0.0, 0.0, 0.0);
    }
    (
        values.iter().copied().fold(f32::INFINITY, f32::min),
        values.iter().sum::<f32>() / values.len() as f32,
        values.iter().copied().fold(f32::NEG_INFINITY, f32::max),
    )
}

fn format_optional_mean(values: &[f32]) -> String {
    if values.is_empty() {
        "n/a".to_string()
    } else {
        format_float(values.iter().sum::<f32>() / values.len() as f32)
    }
}

#[derive(Clone)]
struct IndexedSample {
    sequence_index: usize,
    sequence_row: usize,
    sample_index: usize,
    position: Vec3,
}

struct NearestNeighborSummary {
    nearest_distances: Vec<f32>,
    min: f32,
    mean: f32,
    median: f32,
    max: f32,
    std: f32,
}

fn flatten_sample_groups(groups: &[Vec<[f32; 3]>]) -> Vec<Vec3> {
    groups
        .iter()
        .flat_map(|group| group.iter().map(|point| Vec3::from_array(*point)))
        .collect()
}

fn flatten_indexed_sample_groups(groups: &[Vec<[f32; 3]>]) -> Vec<IndexedSample> {
    groups
        .iter()
        .enumerate()
        .flat_map(|(sequence_index, group)| {
            group
                .iter()
                .enumerate()
                .map(move |(row, point)| IndexedSample {
                    sequence_index: sequence_index + 1,
                    sequence_row: row + 1,
                    sample_index: 0,
                    position: Vec3::from_array(*point),
                })
        })
        .enumerate()
        .map(|(sample_index, mut sample)| {
            sample.sample_index = sample_index;
            sample
        })
        .collect()
}

fn bounds_for_points(points: &[Vec3]) -> (Vec3, Vec3) {
    let mut min = Vec3::splat(f32::INFINITY);
    let mut max = Vec3::splat(f32::NEG_INFINITY);
    for point in points {
        min = min.min(*point);
        max = max.max(*point);
    }
    (min, max)
}

fn covariance_matrix(points: &[Vec3], centroid: Vec3) -> glam::Mat3 {
    let mut xx = 0.0_f32;
    let mut xy = 0.0_f32;
    let mut xz = 0.0_f32;
    let mut yy = 0.0_f32;
    let mut yz = 0.0_f32;
    let mut zz = 0.0_f32;
    for point in points {
        let d = *point - centroid;
        xx += d.x * d.x;
        xy += d.x * d.y;
        xz += d.x * d.z;
        yy += d.y * d.y;
        yz += d.y * d.z;
        zz += d.z * d.z;
    }
    let scale = 1.0 / points.len().max(1) as f32;
    glam::Mat3::from_cols_array(&[
        xx * scale,
        xy * scale,
        xz * scale,
        xy * scale,
        yy * scale,
        yz * scale,
        xz * scale,
        yz * scale,
        zz * scale,
    ])
}

fn principal_components(covariance: glam::Mat3) -> [(f32, Vec3); 3] {
    let mut components = [(0.0_f32, Vec3::X), (0.0_f32, Vec3::Y), (0.0_f32, Vec3::Z)];
    let mut matrix = covariance;
    for component in &mut components {
        let direction = power_iteration(matrix);
        let eigenvalue = direction.dot(matrix * direction).max(0.0);
        *component = (eigenvalue, direction);
        matrix -= eigenvalue * outer_product(direction);
    }
    components.sort_by(|a, b| b.0.partial_cmp(&a.0).unwrap_or(std::cmp::Ordering::Equal));
    if components[2].1.length_squared() < 1.0e-6 {
        components[2].1 = components[0].1.cross(components[1].1).normalize_or_zero();
    }
    components
}

fn power_iteration(matrix: glam::Mat3) -> Vec3 {
    let mut v = Vec3::new(0.71, 0.53, 0.47).normalize();
    for _ in 0..16 {
        let next = matrix * v;
        if next.length_squared() <= 1.0e-8 {
            break;
        }
        v = next.normalize();
    }
    v.normalize_or_zero()
}

fn outer_product(v: Vec3) -> glam::Mat3 {
    glam::Mat3::from_cols(
        Vec3::new(v.x * v.x, v.x * v.y, v.x * v.z),
        Vec3::new(v.y * v.x, v.y * v.y, v.y * v.z),
        Vec3::new(v.z * v.x, v.z * v.y, v.z * v.z),
    )
}

fn exact_duplicate_rows(samples: &[IndexedSample]) -> Vec<IndexedSample> {
    let mut seen = std::collections::HashMap::new();
    let mut duplicates = Vec::new();
    for sample in samples {
        let key = (
            sample.position.x.to_bits(),
            sample.position.y.to_bits(),
            sample.position.z.to_bits(),
        );
        if seen.insert(key, sample.sequence_row).is_some() {
            duplicates.push(sample.clone());
        }
    }
    duplicates
}

fn nearest_neighbor_stats(samples: &[IndexedSample]) -> Option<NearestNeighborSummary> {
    if samples.len() < 2 {
        return None;
    }
    let mut nearest_distances = vec![f32::INFINITY; samples.len()];
    for i in 0..samples.len() {
        for j in (i + 1)..samples.len() {
            let distance = samples[i].position.distance(samples[j].position);
            nearest_distances[i] = nearest_distances[i].min(distance);
            nearest_distances[j] = nearest_distances[j].min(distance);
        }
    }
    let filtered = nearest_distances
        .iter()
        .copied()
        .filter(|distance| distance.is_finite())
        .collect::<Vec<_>>();
    if filtered.is_empty() {
        return None;
    }
    let mut sorted = filtered.clone();
    sorted.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
    let (mean, std) = mean_std(&filtered);
    Some(NearestNeighborSummary {
        min: *sorted.first().unwrap_or(&0.0),
        mean,
        median: sorted[sorted.len() / 2],
        max: *sorted.last().unwrap_or(&0.0),
        std,
        nearest_distances,
    })
}

fn mean_std(values: &[f32]) -> (f32, f32) {
    if values.is_empty() {
        return (0.0, 0.0);
    }
    let mean = values.iter().sum::<f32>() / values.len() as f32;
    let variance = values
        .iter()
        .map(|value| {
            let d = *value - mean;
            d * d
        })
        .sum::<f32>()
        / values.len() as f32;
    (mean, variance.sqrt())
}

fn monotonicity_label(group: &[[f32; 3]], axis: usize) -> String {
    if axis_is_monotonic(group, axis) {
        let first = axis_value(Vec3::from_array(group[0]), axis);
        let last = axis_value(Vec3::from_array(*group.last().unwrap_or(&group[0])), axis);
        if last >= first {
            "increasing".to_string()
        } else {
            "decreasing".to_string()
        }
    } else {
        "mixed".to_string()
    }
}

fn format_float(value: f32) -> String {
    format!("{value:.5}")
}

fn format_vec3(value: Vec3) -> String {
    format!("{:.5}, {:.5}, {:.5}", value.x, value.y, value.z)
}

fn sample_positions_table(groups: &[Vec<[f32; 3]>]) -> AnalysisTable {
    AnalysisTable {
        title: "Sample Positions".to_string(),
        columns: vec![
            "sequence".to_string(),
            "row".to_string(),
            "x".to_string(),
            "y".to_string(),
            "z".to_string(),
        ],
        rows: groups
            .iter()
            .enumerate()
            .flat_map(|(sequence_idx, group)| {
                group.iter().enumerate().map(move |(row_idx, point)| {
                    vec![
                        (sequence_idx + 1).to_string(),
                        (row_idx + 1).to_string(),
                        format_float(point[0]),
                        format_float(point[1]),
                        format_float(point[2]),
                    ]
                })
            })
            .collect(),
    }
}

fn flagged_samples_table<'a>(
    title: impl Into<String>,
    metric_label: &str,
    rows: impl IntoIterator<Item = (&'a IndexedSample, f32)>,
) -> AnalysisTable {
    AnalysisTable {
        title: title.into(),
        columns: vec![
            "sequence".to_string(),
            "row".to_string(),
            "sample".to_string(),
            "x".to_string(),
            "y".to_string(),
            "z".to_string(),
            metric_label.to_string(),
        ],
        rows: rows
            .into_iter()
            .map(|(sample, metric)| {
                vec![
                    sample.sequence_index.to_string(),
                    sample.sequence_row.to_string(),
                    (sample.sample_index + 1).to_string(),
                    format_float(sample.position.x),
                    format_float(sample.position.y),
                    format_float(sample.position.z),
                    format_float(metric),
                ]
            })
            .collect(),
    }
}

fn sample_position_occurrences(samples: &[IndexedSample], position: Vec3) -> f32 {
    samples
        .iter()
        .filter(|sample| sample.position == position)
        .count() as f32
}

fn sample_warning(sample: &IndexedSample, message: String) -> Diagnostic {
    Diagnostic::warning(DiagnosticKind::Validation, message).with_location(
        DiagnosticLocation::new()
            .with_component(format!("sequence {} sample", sample.sequence_index))
            .with_row(sample.sequence_row),
    )
}

fn fit_curve_group(
    source: &PlotSpec,
    group: &[[f32; 3]],
    options: &CurveFitOptions,
) -> Result<CurveFitResult, AnalysisError> {
    match &source.definition {
        crate::PlotDefinition::ExprCartesianLine {
            dep_var, ind_var, ..
        } => fit_cartesian_line_group(group, dep_var.as_str(), ind_var.as_str(), options),
        _ => fit_parametric_curve_group(group, options),
    }
}

fn fit_cartesian_line_group(
    group: &[[f32; 3]],
    dep_var: &str,
    ind_var: &str,
    options: &CurveFitOptions,
) -> Result<CurveFitResult, AnalysisError> {
    if group.len() < 2 {
        return Err(AnalysisError::invalid(
            "Need at least two points to fit a curve.",
        ));
    }
    let xs: Vec<f64> = group
        .iter()
        .filter_map(|point| {
            cartesian_axis_value(Vec3::from_array(*point), ind_var).map(|v| v as f64)
        })
        .collect();
    let ys: Vec<f64> = group
        .iter()
        .filter_map(|point| {
            cartesian_axis_value(Vec3::from_array(*point), dep_var).map(|v| v as f64)
        })
        .collect();
    if xs.len() != group.len() || ys.len() != group.len() {
        return Err(AnalysisError::invalid(
            "Could not read curve axes for fitting.",
        ));
    }
    let evaluation_count = resampled_point_count(group.len(), options.samples_per_segment);
    let x_eval = evenly_spaced_range(*xs.first().unwrap(), *xs.last().unwrap(), evaluation_count);

    match options.method {
        CurveFitMethod::Polynomial => {
            let fit = polynomial_fit(&xs, &ys, options.degree, None)?;
            let fitted_group = x_eval
                .iter()
                .map(|x| {
                    cartesian_line_point(dep_var, ind_var, *x as f32, fit.evaluate(*x) as f32)
                        .to_array()
                })
                .collect();
            let residual_values = xs
                .iter()
                .zip(&ys)
                .map(|(x, y)| (fit.evaluate(*x) - *y) as f32)
                .collect();
            Ok(CurveFitResult {
                fitted_group,
                residual_values,
            })
        }
        CurveFitMethod::RobustPolynomial => {
            let fit = robust_polynomial_fit_scalar(&xs, &ys, options.degree)?;
            let fitted_group = x_eval
                .iter()
                .map(|x| {
                    cartesian_line_point(dep_var, ind_var, *x as f32, fit.evaluate(*x) as f32)
                        .to_array()
                })
                .collect();
            let residual_values = xs
                .iter()
                .zip(&ys)
                .map(|(x, y)| (fit.evaluate(*x) - *y) as f32)
                .collect();
            Ok(CurveFitResult {
                fitted_group,
                residual_values,
            })
        }
        CurveFitMethod::Fourier => {
            let fit = fourier_fit(&xs, &ys, options.harmonics)?;
            let fitted_group = x_eval
                .iter()
                .map(|x| {
                    cartesian_line_point(dep_var, ind_var, *x as f32, fit.evaluate(*x) as f32)
                        .to_array()
                })
                .collect();
            let residual_values = xs
                .iter()
                .zip(&ys)
                .map(|(x, y)| (fit.evaluate(*x) - *y) as f32)
                .collect();
            Ok(CurveFitResult {
                fitted_group,
                residual_values,
            })
        }
        CurveFitMethod::Spline => {
            let smoothed = smooth_scalar_values(&ys, options.smoothing_window);
            let control_points: Vec<[f32; 3]> = xs
                .iter()
                .zip(smoothed.iter())
                .map(|(x, y)| {
                    cartesian_line_point(dep_var, ind_var, *x as f32, *y as f32).to_array()
                })
                .collect();
            let fitted_group = sample_curve_points(
                &control_points
                    .iter()
                    .map(|point| Vec3::from_array(*point))
                    .collect::<Vec<_>>(),
                CurveInterpolation {
                    kind: CurveInterpolationKind::CentripetalCatmullRom,
                    samples_per_segment: options.samples_per_segment,
                    closed: false,
                    smoothing_window: options.smoothing_window,
                },
            )
            .into_iter()
            .map(|point| point.to_array())
            .collect();
            let residual_values = ys
                .iter()
                .zip(smoothed.iter())
                .map(|(observed, fitted)| (*fitted - *observed) as f32)
                .collect();
            Ok(CurveFitResult {
                fitted_group,
                residual_values,
            })
        }
    }
}

fn fit_parametric_curve_group(
    group: &[[f32; 3]],
    options: &CurveFitOptions,
) -> Result<CurveFitResult, AnalysisError> {
    if group.len() < 2 {
        return Err(AnalysisError::invalid(
            "Need at least two points to fit a curve.",
        ));
    }
    let ts = evenly_spaced_parameter_values(group.len());
    let evaluation_ts = evenly_spaced_parameter_values(resampled_point_count(
        group.len(),
        options.samples_per_segment,
    ));
    let positions: Vec<Vec3> = group.iter().map(|point| Vec3::from_array(*point)).collect();

    match options.method {
        CurveFitMethod::Polynomial => {
            let fit = polynomial_fit_vector(&ts, &positions, options.degree, None)?;
            let fitted_group = evaluation_ts
                .iter()
                .map(|t| fit.evaluate(*t).to_array())
                .collect();
            let residual_values = ts
                .iter()
                .zip(positions.iter())
                .map(|(t, observed)| fit.evaluate(*t).distance(*observed))
                .collect();
            Ok(CurveFitResult {
                fitted_group,
                residual_values,
            })
        }
        CurveFitMethod::RobustPolynomial => {
            let fit = robust_polynomial_fit_vector(&ts, &positions, options.degree)?;
            let fitted_group = evaluation_ts
                .iter()
                .map(|t| fit.evaluate(*t).to_array())
                .collect();
            let residual_values = ts
                .iter()
                .zip(positions.iter())
                .map(|(t, observed)| fit.evaluate(*t).distance(*observed))
                .collect();
            Ok(CurveFitResult {
                fitted_group,
                residual_values,
            })
        }
        CurveFitMethod::Fourier => {
            let fit = fourier_fit_vector(&ts, &positions, options.harmonics)?;
            let fitted_group = evaluation_ts
                .iter()
                .map(|t| fit.evaluate(*t).to_array())
                .collect();
            let residual_values = ts
                .iter()
                .zip(positions.iter())
                .map(|(t, observed)| fit.evaluate(*t).distance(*observed))
                .collect();
            Ok(CurveFitResult {
                fitted_group,
                residual_values,
            })
        }
        CurveFitMethod::Spline => {
            let smoothed = smooth_vec3_values(&positions, options.smoothing_window);
            let fitted_group = sample_curve_points(
                &smoothed,
                CurveInterpolation {
                    kind: CurveInterpolationKind::CentripetalCatmullRom,
                    samples_per_segment: options.samples_per_segment,
                    closed: false,
                    smoothing_window: options.smoothing_window,
                },
            )
            .into_iter()
            .map(|point| point.to_array())
            .collect();
            let residual_values = positions
                .iter()
                .zip(smoothed.iter())
                .map(|(observed, fitted)| fitted.distance(*observed))
                .collect();
            Ok(CurveFitResult {
                fitted_group,
                residual_values,
            })
        }
    }
}

fn derived_polyline_plot(
    source: &PlotSpec,
    name: String,
    color: [f32; 4],
    groups: Vec<Vec<[f32; 3]>>,
) -> Result<PlotSpec, AnalysisError> {
    if groups.is_empty() {
        return Err(AnalysisError::invalid(
            "The selected plot did not produce enough samples for this analysis.",
        ));
    }
    Ok(PlotSpec {
        name,
        visible: true,
        domain: source.domain.clone(),
        resolution: source.resolution,
        style: PlotStyle {
            colour_mode: ColourMode::Solid(color),
            line_width: 2.25,
            ..PlotStyle::default()
        },
        definition: crate::PlotDefinition::DerivedPolylineGroups { groups },
    })
}

fn make_curve_fit_output(
    source: &PlotSpec,
    options: CurveFitOptions,
) -> Result<AnalysisOutput, AnalysisError> {
    let groups = curve_sample_groups(source)?;
    if groups.is_empty() || groups.iter().all(|group| group.len() < 2) {
        return Err(AnalysisError::invalid(
            "Curve fitting requires at least one sampled curve with two or more points.",
        ));
    }

    let mut fitted_groups = Vec::new();
    let mut residual_groups = Vec::new();
    let mut control_points = Vec::new();
    let mut residual_samples = Vec::new();
    let mut total_points = 0usize;

    for group in &groups {
        if group.len() < 2 {
            continue;
        }
        total_points += group.len();
        control_points.extend(group.iter().copied());
        let fit = fit_curve_group(source, group, &options)?;
        if fit.fitted_group.len() >= 2 {
            fitted_groups.push(fit.fitted_group);
        }
        if fit.residual_values.len() >= 2 {
            let residual_group = match &source.definition {
                crate::PlotDefinition::ExprCartesianLine {
                    dep_var, ind_var, ..
                } => scalar_plot_cartesian_line_group(
                    group,
                    dep_var.as_str(),
                    ind_var.as_str(),
                    &fit.residual_values,
                ),
                _ => scalar_curve_group(group, &fit.residual_values),
            };
            if residual_group.len() >= 2 {
                residual_groups.push(residual_group);
            }
        }
        residual_samples.extend(fit.residual_values.iter().map(|value| value.abs()));
    }

    if fitted_groups.is_empty() {
        return Err(AnalysisError::invalid(
            "Curve fitting could not produce a fitted curve from the selected source.",
        ));
    }

    let mut plots = vec![PlotSpec {
        name: options.output_name.clone().unwrap_or_else(|| {
            format!("{} {}", curve_fit_method_label(options.method), source.name)
        }),
        visible: true,
        domain: source.domain.clone(),
        resolution: source.resolution,
        style: PlotStyle {
            colour_mode: ColourMode::Solid([1.0, 0.72, 0.25, 1.0]),
            line_width: 2.5,
            ..PlotStyle::default()
        },
        definition: crate::PlotDefinition::DerivedPolylineGroups {
            groups: fitted_groups,
        },
    }];

    if options.show_control_points && !control_points.is_empty() {
        plots.push(PlotSpec {
            name: format!("Control Points {}", source.name),
            visible: true,
            domain: source.domain.clone(),
            resolution: source.resolution,
            style: PlotStyle {
                colour_mode: ColourMode::Solid([0.35, 0.85, 1.0, 1.0]),
                point_size: 7.0,
                ..PlotStyle::default()
            },
            definition: crate::PlotDefinition::PointAnnotations {
                points: control_points
                    .iter()
                    .enumerate()
                    .map(|(index, position)| PointAnnotation {
                        position: *position,
                        label: format!("Control {}", index + 1),
                    })
                    .collect(),
                show_labels: false,
            },
        });
    }

    if options.show_residual_plot && !residual_groups.is_empty() {
        plots.push(PlotSpec {
            name: format!("Residuals {}", source.name),
            visible: true,
            domain: source.domain.clone(),
            resolution: source.resolution,
            style: PlotStyle {
                colour_mode: ColourMode::Solid([1.0, 0.35, 0.35, 1.0]),
                line_width: 2.0,
                ..PlotStyle::default()
            },
            definition: crate::PlotDefinition::DerivedPolylineGroups {
                groups: residual_groups,
            },
        });
    }

    let rms = if residual_samples.is_empty() {
        0.0
    } else {
        (residual_samples
            .iter()
            .map(|value| value * value)
            .sum::<f32>()
            / residual_samples.len() as f32)
            .sqrt()
    };
    let max_residual = residual_samples.iter().copied().fold(0.0_f32, f32::max);

    let mut report_values = vec![
        (
            "Method".to_string(),
            curve_fit_method_label(options.method).to_string(),
        ),
        ("Samples".to_string(), total_points.to_string()),
        ("RMS Residual".to_string(), format!("{rms:.6}")),
        ("Max Residual".to_string(), format!("{max_residual:.6}")),
    ];
    match options.method {
        CurveFitMethod::Polynomial | CurveFitMethod::RobustPolynomial => {
            report_values.push(("Degree".to_string(), options.degree.to_string()));
        }
        CurveFitMethod::Fourier => {
            report_values.push(("Harmonics".to_string(), options.harmonics.to_string()));
        }
        CurveFitMethod::Spline => {
            report_values.push((
                "Smoothing Window".to_string(),
                options.smoothing_window.to_string(),
            ));
        }
    }

    Ok(AnalysisOutput::Composite {
        plots,
        reports: vec![AnalysisReport {
            title: format!("Curve Fit {}", source.name),
            values: report_values,
        }],
        tables: Vec::new(),
        diagnostics: Vec::new(),
        frame_fields: Vec::new(),
        provenance: AnalysisProvenance {
            kind: AnalysisKind::FitCurve,
            source_plots: vec![source.name.clone()],
            parameters: vec![
                ("fit_method".to_string(), curve_fit_method_key(options.method).to_string()),
                ("degree".to_string(), options.degree.to_string()),
                ("harmonics".to_string(), options.harmonics.to_string()),
                ("smoothing_window".to_string(), options.smoothing_window.to_string()),
                (
                    "samples_per_segment".to_string(),
                    options.samples_per_segment.to_string(),
                ),
            ],
            notes: vec![
                "Derived output includes a fitted curve plus optional control-point preview and residual plot.".to_string(),
            ],
        },
    })
}

#[derive(Clone, Copy, Debug, PartialEq, Eq)]
enum CurveFitMethod {
    Polynomial,
    RobustPolynomial,
    Spline,
    Fourier,
}

#[derive(Clone, Debug)]
struct CurveFitOptions {
    method: CurveFitMethod,
    output_name: Option<String>,
    degree: usize,
    harmonics: usize,
    smoothing_window: u32,
    samples_per_segment: u32,
    show_control_points: bool,
    show_residual_plot: bool,
}

#[derive(Clone, Debug)]
struct CurveFitResult {
    fitted_group: Vec<[f32; 3]>,
    residual_values: Vec<f32>,
}

fn parameter_map(parameters: &[(String, String)]) -> HashMap<String, String> {
    parameters.iter().cloned().collect()
}

fn build_curve_fit_options(params: &HashMap<String, String>) -> CurveFitOptions {
    CurveFitOptions {
        method: parse_curve_fit_method(params.get("fit_method").map(String::as_str))
            .unwrap_or(CurveFitMethod::Polynomial),
        output_name: params.get("output_name").cloned(),
        degree: params
            .get("degree")
            .and_then(|value| value.parse::<usize>().ok())
            .unwrap_or(5)
            .clamp(1, 12),
        harmonics: params
            .get("harmonics")
            .and_then(|value| value.parse::<usize>().ok())
            .unwrap_or(3)
            .clamp(1, 16),
        smoothing_window: normalized_window_value(
            params
                .get("smoothing_window")
                .and_then(|value| value.parse::<u32>().ok())
                .unwrap_or(7),
        ),
        samples_per_segment: params
            .get("samples_per_segment")
            .and_then(|value| value.parse::<u32>().ok())
            .unwrap_or(8)
            .clamp(1, 64),
        show_control_points: params
            .get("show_control_points")
            .is_none_or(|value| matches!(value.as_str(), "1" | "true" | "yes")),
        show_residual_plot: params
            .get("show_residual_plot")
            .is_none_or(|value| matches!(value.as_str(), "1" | "true" | "yes")),
    }
}

fn build_interpolation(params: &HashMap<String, String>) -> CurveInterpolation {
    CurveInterpolation {
        kind: parse_interpolation_kind(params.get("interpolation_kind").map(String::as_str))
            .unwrap_or(CurveInterpolationKind::Linear),
        samples_per_segment: params
            .get("samples_per_segment")
            .and_then(|value| value.parse::<u32>().ok())
            .unwrap_or(1),
        closed: params
            .get("closed")
            .is_some_and(|value| matches!(value.as_str(), "1" | "true" | "yes")),
        smoothing_window: normalized_window_value(
            params
                .get("smoothing_window")
                .and_then(|value| value.parse::<u32>().ok())
                .unwrap_or(5),
        ),
    }
}

fn parse_curve_fit_method(value: Option<&str>) -> Option<CurveFitMethod> {
    match value? {
        "polynomial" => Some(CurveFitMethod::Polynomial),
        "robust_polynomial" => Some(CurveFitMethod::RobustPolynomial),
        "spline" => Some(CurveFitMethod::Spline),
        "fourier" => Some(CurveFitMethod::Fourier),
        _ => None,
    }
}

fn curve_fit_method_label(method: CurveFitMethod) -> &'static str {
    match method {
        CurveFitMethod::Polynomial => "Fit (Polynomial)",
        CurveFitMethod::RobustPolynomial => "Fit (Robust Polynomial)",
        CurveFitMethod::Spline => "Fit (Spline / Smoothed Catmull-Rom)",
        CurveFitMethod::Fourier => "Fit (Fourier Series)",
    }
}

fn curve_fit_method_key(method: CurveFitMethod) -> &'static str {
    match method {
        CurveFitMethod::Polynomial => "polynomial",
        CurveFitMethod::RobustPolynomial => "robust_polynomial",
        CurveFitMethod::Spline => "spline",
        CurveFitMethod::Fourier => "fourier",
    }
}

#[derive(Clone, Debug)]
struct PolynomialFit {
    coeffs: Vec<f64>,
    x_min: f64,
    x_max: f64,
}

impl PolynomialFit {
    fn evaluate(&self, x: f64) -> f64 {
        evaluate_polynomial(
            &self.coeffs,
            normalize_domain_value(x, self.x_min, self.x_max),
        )
    }
}

#[derive(Clone, Debug)]
struct PolynomialVectorFit {
    x: PolynomialFit,
    y: PolynomialFit,
    z: PolynomialFit,
}

impl PolynomialVectorFit {
    fn evaluate(&self, t: f64) -> Vec3 {
        Vec3::new(
            self.x.evaluate(t) as f32,
            self.y.evaluate(t) as f32,
            self.z.evaluate(t) as f32,
        )
    }
}

#[derive(Clone, Debug)]
struct FourierFit {
    constant: f64,
    cos_coeffs: Vec<f64>,
    sin_coeffs: Vec<f64>,
    x_min: f64,
    x_max: f64,
}

impl FourierFit {
    fn evaluate(&self, x: f64) -> f64 {
        let theta = normalized_angle(x, self.x_min, self.x_max);
        let mut value = self.constant;
        for (index, (cos_coeff, sin_coeff)) in self
            .cos_coeffs
            .iter()
            .zip(self.sin_coeffs.iter())
            .enumerate()
        {
            let k = index as f64 + 1.0;
            value += *cos_coeff * (k * theta).cos() + *sin_coeff * (k * theta).sin();
        }
        value
    }
}

#[derive(Clone, Debug)]
struct FourierVectorFit {
    x: FourierFit,
    y: FourierFit,
    z: FourierFit,
}

impl FourierVectorFit {
    fn evaluate(&self, t: f64) -> Vec3 {
        Vec3::new(
            self.x.evaluate(t) as f32,
            self.y.evaluate(t) as f32,
            self.z.evaluate(t) as f32,
        )
    }
}

fn polynomial_fit(
    xs: &[f64],
    ys: &[f64],
    degree: usize,
    weights: Option<&[f64]>,
) -> Result<PolynomialFit, AnalysisError> {
    let x_min = *xs
        .first()
        .ok_or_else(|| AnalysisError::invalid("No samples for fit."))?;
    let x_max = *xs
        .last()
        .ok_or_else(|| AnalysisError::invalid("No samples for fit."))?;
    let normalized: Vec<f64> = xs
        .iter()
        .map(|x| normalize_domain_value(*x, x_min, x_max))
        .collect();
    let coeffs = weighted_polynomial_coefficients(&normalized, ys, degree, weights)?;
    Ok(PolynomialFit {
        coeffs,
        x_min,
        x_max,
    })
}

fn robust_polynomial_fit_scalar(
    xs: &[f64],
    ys: &[f64],
    degree: usize,
) -> Result<PolynomialFit, AnalysisError> {
    let mut weights = vec![1.0_f64; xs.len()];
    let mut fit = polynomial_fit(xs, ys, degree, Some(&weights))?;
    for _ in 0..5 {
        let residuals: Vec<f64> = xs
            .iter()
            .zip(ys.iter())
            .map(|(x, y)| (fit.evaluate(*x) - *y).abs())
            .collect();
        weights = huber_weights(&residuals);
        fit = polynomial_fit(xs, ys, degree, Some(&weights))?;
    }
    Ok(fit)
}

fn polynomial_fit_vector(
    ts: &[f64],
    positions: &[Vec3],
    degree: usize,
    weights: Option<&[f64]>,
) -> Result<PolynomialVectorFit, AnalysisError> {
    let xs: Vec<f64> = positions.iter().map(|position| position.x as f64).collect();
    let ys: Vec<f64> = positions.iter().map(|position| position.y as f64).collect();
    let zs: Vec<f64> = positions.iter().map(|position| position.z as f64).collect();
    Ok(PolynomialVectorFit {
        x: polynomial_fit(ts, &xs, degree, weights)?,
        y: polynomial_fit(ts, &ys, degree, weights)?,
        z: polynomial_fit(ts, &zs, degree, weights)?,
    })
}

fn robust_polynomial_fit_vector(
    ts: &[f64],
    positions: &[Vec3],
    degree: usize,
) -> Result<PolynomialVectorFit, AnalysisError> {
    let mut weights = vec![1.0_f64; ts.len()];
    let mut fit = polynomial_fit_vector(ts, positions, degree, Some(&weights))?;
    for _ in 0..5 {
        let residuals: Vec<f64> = ts
            .iter()
            .zip(positions.iter())
            .map(|(t, position)| fit.evaluate(*t).distance(*position) as f64)
            .collect();
        weights = huber_weights(&residuals);
        fit = polynomial_fit_vector(ts, positions, degree, Some(&weights))?;
    }
    Ok(fit)
}

fn fourier_fit(xs: &[f64], ys: &[f64], harmonics: usize) -> Result<FourierFit, AnalysisError> {
    let x_min = *xs
        .first()
        .ok_or_else(|| AnalysisError::invalid("No samples for fit."))?;
    let x_max = *xs
        .last()
        .ok_or_else(|| AnalysisError::invalid("No samples for fit."))?;
    let design = xs
        .iter()
        .map(|x| fourier_basis(normalized_angle(*x, x_min, x_max), harmonics))
        .collect::<Vec<_>>();
    let coeffs = linear_least_squares(&design, ys, None)?;
    Ok(FourierFit {
        constant: coeffs[0],
        cos_coeffs: (0..harmonics).map(|index| coeffs[1 + index * 2]).collect(),
        sin_coeffs: (0..harmonics).map(|index| coeffs[2 + index * 2]).collect(),
        x_min,
        x_max,
    })
}

fn fourier_fit_vector(
    ts: &[f64],
    positions: &[Vec3],
    harmonics: usize,
) -> Result<FourierVectorFit, AnalysisError> {
    let xs: Vec<f64> = positions.iter().map(|position| position.x as f64).collect();
    let ys: Vec<f64> = positions.iter().map(|position| position.y as f64).collect();
    let zs: Vec<f64> = positions.iter().map(|position| position.z as f64).collect();
    Ok(FourierVectorFit {
        x: fourier_fit(ts, &xs, harmonics)?,
        y: fourier_fit(ts, &ys, harmonics)?,
        z: fourier_fit(ts, &zs, harmonics)?,
    })
}

fn weighted_polynomial_coefficients(
    xs: &[f64],
    ys: &[f64],
    degree: usize,
    weights: Option<&[f64]>,
) -> Result<Vec<f64>, AnalysisError> {
    let design = xs
        .iter()
        .map(|x| polynomial_basis(*x, degree))
        .collect::<Vec<_>>();
    linear_least_squares(&design, ys, weights)
}

fn linear_least_squares(
    design: &[Vec<f64>],
    ys: &[f64],
    weights: Option<&[f64]>,
) -> Result<Vec<f64>, AnalysisError> {
    if design.is_empty() || ys.is_empty() || design.len() != ys.len() {
        return Err(AnalysisError::invalid("Invalid least-squares system."));
    }
    let cols = design[0].len();
    let mut ata = vec![vec![0.0_f64; cols]; cols];
    let mut atb = vec![0.0_f64; cols];
    for (row_index, row) in design.iter().enumerate() {
        let weight = weights
            .and_then(|values| values.get(row_index).copied())
            .unwrap_or(1.0)
            .max(1.0e-8);
        for i in 0..cols {
            atb[i] += weight * row[i] * ys[row_index];
            for j in 0..cols {
                ata[i][j] += weight * row[i] * row[j];
            }
        }
    }
    solve_linear_system(ata, atb)
}

fn solve_linear_system(
    mut matrix: Vec<Vec<f64>>,
    mut rhs: Vec<f64>,
) -> Result<Vec<f64>, AnalysisError> {
    let size = rhs.len();
    for pivot in 0..size {
        let mut best_row = pivot;
        let mut best_value = matrix[pivot][pivot].abs();
        for row in (pivot + 1)..size {
            let value = matrix[row][pivot].abs();
            if value > best_value {
                best_value = value;
                best_row = row;
            }
        }
        if best_value <= 1.0e-10 {
            return Err(AnalysisError::invalid(
                "Curve fit is ill-conditioned for the selected method and parameters.",
            ));
        }
        if best_row != pivot {
            matrix.swap(pivot, best_row);
            rhs.swap(pivot, best_row);
        }
        let scale = matrix[pivot][pivot];
        for col in pivot..size {
            matrix[pivot][col] /= scale;
        }
        rhs[pivot] /= scale;
        for row in 0..size {
            if row == pivot {
                continue;
            }
            let factor = matrix[row][pivot];
            if factor.abs() <= 1.0e-12 {
                continue;
            }
            for col in pivot..size {
                matrix[row][col] -= factor * matrix[pivot][col];
            }
            rhs[row] -= factor * rhs[pivot];
        }
    }
    Ok(rhs)
}

fn polynomial_basis(x: f64, degree: usize) -> Vec<f64> {
    let mut basis = vec![1.0_f64; degree + 1];
    for index in 1..=degree {
        basis[index] = basis[index - 1] * x;
    }
    basis
}

fn fourier_basis(theta: f64, harmonics: usize) -> Vec<f64> {
    let mut basis = Vec::with_capacity(1 + harmonics * 2);
    basis.push(1.0);
    for harmonic in 1..=harmonics {
        let angle = theta * harmonic as f64;
        basis.push(angle.cos());
        basis.push(angle.sin());
    }
    basis
}

fn normalize_domain_value(x: f64, x_min: f64, x_max: f64) -> f64 {
    let range = (x_max - x_min).abs();
    if range <= f64::EPSILON {
        0.0
    } else {
        ((x - x_min) / range) * 2.0 - 1.0
    }
}

fn normalized_angle(x: f64, x_min: f64, x_max: f64) -> f64 {
    let range = (x_max - x_min).abs();
    if range <= f64::EPSILON {
        0.0
    } else {
        ((x - x_min) / range) * 2.0 * PI
    }
}

fn evaluate_polynomial(coeffs: &[f64], x: f64) -> f64 {
    coeffs.iter().rev().fold(0.0, |acc, coeff| acc * x + coeff)
}

fn huber_weights(residuals: &[f64]) -> Vec<f64> {
    let scale = residual_scale(residuals).max(1.0e-6);
    residuals
        .iter()
        .map(|residual| {
            let normalized = residual.abs() / (1.5 * scale);
            if normalized <= 1.0 {
                1.0
            } else {
                1.0 / normalized
            }
        })
        .collect()
}

fn residual_scale(residuals: &[f64]) -> f64 {
    if residuals.is_empty() {
        return 1.0;
    }
    let mut sorted = residuals
        .iter()
        .map(|value| value.abs())
        .collect::<Vec<_>>();
    sorted.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
    sorted[sorted.len() / 2]
}

fn smooth_scalar_values(values: &[f64], smoothing_window: u32) -> Vec<f64> {
    let radius = normalized_window_value(smoothing_window) as isize / 2;
    (0..values.len())
        .map(|index| {
            let mut total = 0.0_f64;
            let mut count = 0.0_f64;
            for offset in -radius..=radius {
                let sample_index =
                    (index as isize + offset).clamp(0, values.len() as isize - 1) as usize;
                total += values[sample_index];
                count += 1.0;
            }
            total / count
        })
        .collect()
}

fn smooth_vec3_values(values: &[Vec3], smoothing_window: u32) -> Vec<Vec3> {
    let radius = normalized_window_value(smoothing_window) as isize / 2;
    (0..values.len())
        .map(|index| {
            let mut total = Vec3::ZERO;
            let mut count = 0.0_f32;
            for offset in -radius..=radius {
                let sample_index =
                    (index as isize + offset).clamp(0, values.len() as isize - 1) as usize;
                total += values[sample_index];
                count += 1.0;
            }
            total / count
        })
        .collect()
}

fn resampled_point_count(control_point_count: usize, samples_per_segment: u32) -> usize {
    let segments = control_point_count.saturating_sub(1);
    (segments * samples_per_segment.max(1) as usize + 1).max(control_point_count)
}

fn evenly_spaced_range(start: f64, end: f64, count: usize) -> Vec<f64> {
    if count <= 1 {
        return vec![start];
    }
    (0..count)
        .map(|index| start + (end - start) * index as f64 / (count - 1) as f64)
        .collect()
}

fn evenly_spaced_parameter_values(count: usize) -> Vec<f64> {
    evenly_spaced_range(0.0, 1.0, count)
}

fn parse_interpolation_kind(value: Option<&str>) -> Option<CurveInterpolationKind> {
    match value? {
        "linear" => Some(CurveInterpolationKind::Linear),
        "catmull_rom" => Some(CurveInterpolationKind::CatmullRom),
        "centripetal_catmull_rom" => Some(CurveInterpolationKind::CentripetalCatmullRom),
        "moving_average" => Some(CurveInterpolationKind::MovingAverage),
        "savitzky_golay" => Some(CurveInterpolationKind::SavitzkyGolay),
        _ => None,
    }
}

fn parse_axis(value: Option<&str>) -> Option<SliceAxis> {
    match value?.to_ascii_lowercase().as_str() {
        "x" => Some(SliceAxis::X),
        "y" => Some(SliceAxis::Y),
        "z" => Some(SliceAxis::Z),
        _ => None,
    }
}

fn parse_axis_index(value: Option<&str>) -> Option<usize> {
    match value?.to_ascii_lowercase().as_str() {
        "x" | "0" => Some(0),
        "y" | "1" => Some(1),
        "z" | "2" => Some(2),
        _ => None,
    }
}

fn evenly_spaced_isovalues(count: usize) -> Vec<f32> {
    let count = count.max(1);
    if count == 1 {
        return vec![0.0];
    }
    (0..count)
        .map(|i| -0.9 + 1.8 * i as f32 / (count - 1) as f32)
        .collect()
}

fn parse_error(error: impl ToString) -> AnalysisError {
    AnalysisError {
        diagnostic: Diagnostic::error(DiagnosticKind::Parse, error.to_string()),
    }
}

fn table_errors(errors: Vec<crate::TableValidationError>) -> AnalysisError {
    let summary = errors
        .into_iter()
        .map(|error| error.display())
        .collect::<Vec<_>>()
        .join("; ");
    AnalysisError::invalid(summary)
}

fn make_point_annotations(points: &[[f32; 3]], _prefix: &str) -> Vec<PointAnnotation> {
    points
        .iter()
        .map(|position| PointAnnotation {
            position: *position,
            label: String::new(),
        })
        .collect()
}

fn axis_name(axis: usize) -> &'static str {
    match axis {
        0 => "x",
        1 => "y",
        _ => "z",
    }
}

fn normalized_window_value(window: u32) -> u32 {
    let mut normalized = window.max(3);
    if normalized % 2 == 0 {
        normalized += 1;
    }
    normalized
}

fn derivative_curve_group(group: &[[f32; 3]]) -> Vec<[f32; 3]> {
    if group.len() < 2 {
        return Vec::new();
    }
    if let Some(layout) = detect_planar_graph_layout(group) {
        return derivative_planar_graph_group(group, layout);
    }
    (0..group.len())
        .map(|index| finite_difference(group, index).to_array())
        .collect()
}

fn tangent_curve_group(group: &[[f32; 3]]) -> Vec<[f32; 3]> {
    if group.len() < 2 {
        return Vec::new();
    }
    (0..group.len())
        .map(|index| {
            finite_difference(group, index)
                .normalize_or_zero()
                .to_array()
        })
        .collect()
}

fn finite_difference(group: &[[f32; 3]], index: usize) -> Vec3 {
    let current = Vec3::from_array(group[index]);
    if index == 0 {
        let next = Vec3::from_array(group[1]);
        return next - current;
    }
    if index + 1 == group.len() {
        let prev = Vec3::from_array(group[index - 1]);
        return current - prev;
    }
    let prev = Vec3::from_array(group[index - 1]);
    let next = Vec3::from_array(group[index + 1]);
    (next - prev) * 0.5
}

fn integral_curve_group(group: &[[f32; 3]], normalize_integral: bool) -> Vec<[f32; 3]> {
    if group.len() < 2 {
        return Vec::new();
    }
    if let Some(layout) = detect_planar_graph_layout(group) {
        return integral_planar_graph_group(group, layout, normalize_integral);
    }
    let mut out = Vec::with_capacity(group.len());
    let mut accum = Vec3::ZERO;
    out.push(accum.to_array());
    let dt = 1.0_f32 / (group.len() - 1) as f32;
    for pair in group.windows(2) {
        let a = Vec3::from_array(pair[0]);
        let b = Vec3::from_array(pair[1]);
        accum += (a + b) * 0.5 * dt;
        out.push(accum.to_array());
    }
    out
}

fn cumulative_arc_lengths(group: &[[f32; 3]]) -> Vec<f32> {
    if group.is_empty() {
        return Vec::new();
    }
    let mut out = Vec::with_capacity(group.len());
    let mut total = 0.0_f32;
    out.push(total);
    for pair in group.windows(2) {
        let a = Vec3::from_array(pair[0]);
        let b = Vec3::from_array(pair[1]);
        total += b.distance(a);
        out.push(total);
    }
    out
}

fn curvature_values(group: &[[f32; 3]]) -> Vec<f32> {
    if group.len() < 3 {
        return vec![0.0; group.len()];
    }
    let mut values = vec![0.0_f32; group.len()];
    for index in 1..(group.len() - 1) {
        let a = Vec3::from_array(group[index - 1]);
        let b = Vec3::from_array(group[index]);
        let c = Vec3::from_array(group[index + 1]);
        let ab = b - a;
        let bc = c - b;
        let ac = c - a;
        let denom = ab.length() * bc.length() * ac.length();
        if denom > 1.0e-6 {
            values[index] = 2.0 * ab.cross(ac).length() / denom;
        }
    }
    values[0] = values[1];
    values[group.len() - 1] = values[group.len() - 2];
    values
}

fn scalar_curve_group(group: &[[f32; 3]], values: &[f32]) -> Vec<[f32; 3]> {
    if group.len() != values.len() || group.len() < 2 {
        return Vec::new();
    }
    let denom = (group.len() - 1) as f32;
    values
        .iter()
        .enumerate()
        .map(|(index, value)| Vec3::new(index as f32 / denom, *value, 0.0).to_array())
        .collect()
}

fn scalar_plot_cartesian_line_group(
    group: &[[f32; 3]],
    dep_var: &str,
    ind_var: &str,
    values: &[f32],
) -> Vec<[f32; 3]> {
    if group.len() != values.len() || group.len() < 2 {
        return Vec::new();
    }
    group
        .iter()
        .zip(values.iter())
        .filter_map(|(point, value)| {
            let independent = cartesian_axis_value(Vec3::from_array(*point), ind_var)?;
            Some(cartesian_line_point(dep_var, ind_var, independent, *value).to_array())
        })
        .collect()
}

fn axis_derivative_group(
    group: &[[f32; 3]],
    numerator_axis: usize,
    denominator_axis: usize,
) -> Vec<[f32; 3]> {
    if group.len() < 2 {
        return Vec::new();
    }
    (0..group.len())
        .filter_map(|index| {
            let point = Vec3::from_array(group[index]);
            let denominator = axis_value(point, denominator_axis);
            let derivative = axis_derivative_value(group, index, numerator_axis, denominator_axis)?;
            let constant_axis =
                (0..3).find(|axis| *axis != numerator_axis && *axis != denominator_axis);
            let constant_value = constant_axis.map(|axis| average_axis_value(group, axis));
            Some(
                point_on_axes(
                    denominator_axis,
                    numerator_axis,
                    denominator,
                    derivative,
                    constant_axis.zip(constant_value),
                )
                .to_array(),
            )
        })
        .collect()
}

fn derivative_cartesian_line_group(
    group: &[[f32; 3]],
    dep_var: &str,
    ind_var: &str,
) -> Vec<[f32; 3]> {
    if group.len() < 2 {
        return Vec::new();
    }
    (0..group.len())
        .filter_map(|index| {
            let current = Vec3::from_array(group[index]);
            let current_ind = cartesian_axis_value(current, ind_var)?;
            let derivative = scalar_derivative(group, index, dep_var, ind_var)?;
            Some(cartesian_line_point(dep_var, ind_var, current_ind, derivative).to_array())
        })
        .collect()
}

fn scalar_derivative(
    group: &[[f32; 3]],
    index: usize,
    dep_var: &str,
    ind_var: &str,
) -> Option<f32> {
    if group.len() < 2 {
        return None;
    }
    let (a, b) = if index == 0 {
        (0, 1)
    } else if index + 1 == group.len() {
        (group.len() - 2, group.len() - 1)
    } else {
        (index - 1, index + 1)
    };
    let pa = Vec3::from_array(group[a]);
    let pb = Vec3::from_array(group[b]);
    let da = cartesian_axis_value(pa, dep_var)?;
    let db = cartesian_axis_value(pb, dep_var)?;
    let ia = cartesian_axis_value(pa, ind_var)?;
    let ib = cartesian_axis_value(pb, ind_var)?;
    let denom = ib - ia;
    if denom.abs() <= 1.0e-6 {
        return Some(0.0);
    }
    Some((db - da) / denom)
}

fn axis_derivative_value(
    group: &[[f32; 3]],
    index: usize,
    numerator_axis: usize,
    denominator_axis: usize,
) -> Option<f32> {
    if group.len() < 2 {
        return None;
    }
    let (a, b) = if index == 0 {
        (0, 1)
    } else if index + 1 == group.len() {
        (group.len() - 2, group.len() - 1)
    } else {
        (index - 1, index + 1)
    };
    let pa = Vec3::from_array(group[a]);
    let pb = Vec3::from_array(group[b]);
    let na = axis_value(pa, numerator_axis);
    let nb = axis_value(pb, numerator_axis);
    let da = axis_value(pa, denominator_axis);
    let db = axis_value(pb, denominator_axis);
    let denom = db - da;
    if denom.abs() <= 1.0e-6 {
        return Some(0.0);
    }
    Some((nb - na) / denom)
}

fn cartesian_axis_value(point: Vec3, axis: &str) -> Option<f32> {
    match axis {
        "x" => Some(point.x),
        "y" => Some(point.y),
        "z" => Some(point.z),
        _ => None,
    }
}

fn axis_value(point: Vec3, axis: usize) -> f32 {
    match axis {
        0 => point.x,
        1 => point.y,
        _ => point.z,
    }
}

fn set_axis_value(point: &mut Vec3, axis: usize, value: f32) {
    match axis {
        0 => point.x = value,
        1 => point.y = value,
        _ => point.z = value,
    }
}

fn point_on_axes(
    independent_axis: usize,
    dependent_axis: usize,
    independent: f32,
    dependent: f32,
    constant_axis: Option<(usize, f32)>,
) -> Vec3 {
    let mut point = Vec3::ZERO;
    set_axis_value(&mut point, independent_axis, independent);
    set_axis_value(&mut point, dependent_axis, dependent);
    if let Some((axis, value)) = constant_axis {
        set_axis_value(&mut point, axis, value);
    }
    point
}

#[derive(Clone, Copy, Debug)]
struct PlanarGraphLayout {
    independent_axis: usize,
    dependent_axis: usize,
    constant_axis: usize,
    constant_value: f32,
}

fn detect_planar_graph_layout(group: &[[f32; 3]]) -> Option<PlanarGraphLayout> {
    if group.len() < 2 {
        return None;
    }
    let ranges = (0..3)
        .map(|axis| axis_range(group, axis))
        .collect::<Vec<_>>();
    let constant_axis = ranges
        .iter()
        .enumerate()
        .min_by(|a, b| a.1.partial_cmp(b.1).unwrap_or(std::cmp::Ordering::Equal))
        .map(|(axis, _)| axis)?;
    if ranges[constant_axis] > 1.0e-3 {
        return None;
    }
    let varying_axes = (0..3)
        .filter(|axis| *axis != constant_axis)
        .collect::<Vec<_>>();
    let monotonic_axes = varying_axes
        .iter()
        .copied()
        .filter(|axis| axis_is_monotonic(group, *axis))
        .collect::<Vec<_>>();
    let independent_axis = if monotonic_axes.len() == 1 {
        monotonic_axes[0]
    } else {
        *varying_axes.iter().max_by(|a, b| {
            ranges[**a]
                .partial_cmp(&ranges[**b])
                .unwrap_or(std::cmp::Ordering::Equal)
        })?
    };
    let dependent_axis = varying_axes
        .into_iter()
        .find(|axis| *axis != independent_axis)?;
    Some(PlanarGraphLayout {
        independent_axis,
        dependent_axis,
        constant_axis,
        constant_value: average_axis_value(group, constant_axis),
    })
}

fn derivative_planar_graph_group(group: &[[f32; 3]], layout: PlanarGraphLayout) -> Vec<[f32; 3]> {
    (0..group.len())
        .filter_map(|index| {
            let point = Vec3::from_array(group[index]);
            let independent = axis_value(point, layout.independent_axis);
            let derivative = axis_derivative_value(
                group,
                index,
                layout.dependent_axis,
                layout.independent_axis,
            )?;
            Some(
                point_on_axes(
                    layout.independent_axis,
                    layout.dependent_axis,
                    independent,
                    derivative,
                    Some((layout.constant_axis, layout.constant_value)),
                )
                .to_array(),
            )
        })
        .collect()
}

fn integral_planar_graph_group(
    group: &[[f32; 3]],
    layout: PlanarGraphLayout,
    normalize_integral: bool,
) -> Vec<[f32; 3]> {
    let mut out = Vec::with_capacity(group.len());
    let start_independent = axis_value(Vec3::from_array(group[0]), layout.independent_axis);
    let mut accum = 0.0_f32;
    out.push(
        point_on_axes(
            layout.independent_axis,
            layout.dependent_axis,
            start_independent,
            accum,
            Some((layout.constant_axis, layout.constant_value)),
        )
        .to_array(),
    );
    for pair in group.windows(2) {
        let a = Vec3::from_array(pair[0]);
        let b = Vec3::from_array(pair[1]);
        let ia = axis_value(a, layout.independent_axis);
        let ib = axis_value(b, layout.independent_axis);
        let da = axis_value(a, layout.dependent_axis);
        let db = axis_value(b, layout.dependent_axis);
        accum += (da + db) * 0.5 * (ib - ia);
        out.push(
            point_on_axes(
                layout.independent_axis,
                layout.dependent_axis,
                ib,
                accum,
                Some((layout.constant_axis, layout.constant_value)),
            )
            .to_array(),
        );
    }
    if normalize_integral {
        normalize_planar_graph_group(&mut out, layout.dependent_axis);
    }
    out
}

fn normalize_planar_graph_group(group: &mut [[f32; 3]], dependent_axis: usize) {
    if group.is_empty() {
        return;
    }
    let mean = group
        .iter()
        .map(|point| axis_value(Vec3::from_array(*point), dependent_axis))
        .sum::<f32>()
        / group.len() as f32;
    for point in group {
        let mut vec = Vec3::from_array(*point);
        let value = axis_value(vec, dependent_axis);
        set_axis_value(&mut vec, dependent_axis, value - mean);
        *point = vec.to_array();
    }
}

fn normalize_cartesian_line_group(group: &mut [[f32; 3]], dep_var: &str) {
    let Some(dep_axis) = parse_axis_name(dep_var) else {
        return;
    };
    normalize_planar_graph_group(group, dep_axis);
}

fn parse_axis_name(axis: &str) -> Option<usize> {
    match axis {
        "x" => Some(0),
        "y" => Some(1),
        "z" => Some(2),
        _ => None,
    }
}

fn axis_range(group: &[[f32; 3]], axis: usize) -> f32 {
    let mut min = f32::INFINITY;
    let mut max = f32::NEG_INFINITY;
    for point in group {
        let value = axis_value(Vec3::from_array(*point), axis);
        min = min.min(value);
        max = max.max(value);
    }
    max - min
}

fn average_axis_value(group: &[[f32; 3]], axis: usize) -> f32 {
    group
        .iter()
        .map(|point| axis_value(Vec3::from_array(*point), axis))
        .sum::<f32>()
        / group.len() as f32
}

fn axis_is_monotonic(group: &[[f32; 3]], axis: usize) -> bool {
    let epsilon = 1.0e-6_f32;
    let mut nondecreasing = true;
    let mut nonincreasing = true;
    for pair in group.windows(2) {
        let a = axis_value(Vec3::from_array(pair[0]), axis);
        let b = axis_value(Vec3::from_array(pair[1]), axis);
        if b + epsilon < a {
            nondecreasing = false;
        }
        if b > a + epsilon {
            nonincreasing = false;
        }
    }
    nondecreasing || nonincreasing
}

fn cartesian_line_point(dep_var: &str, ind_var: &str, independent: f32, dependent: f32) -> Vec3 {
    match (dep_var, ind_var) {
        ("y", "x") => Vec3::new(independent, dependent, 0.0),
        ("z", "x") => Vec3::new(independent, 0.0, dependent),
        ("z", "y") => Vec3::new(0.0, independent, dependent),
        ("x", "y") => Vec3::new(dependent, independent, 0.0),
        ("x", "z") => Vec3::new(dependent, 0.0, independent),
        ("y", "z") => Vec3::new(0.0, dependent, independent),
        _ => Vec3::new(independent, dependent, 0.0),
    }
}

fn integral_cartesian_line_group(
    group: &[[f32; 3]],
    dep_var: &str,
    ind_var: &str,
    normalize_integral: bool,
) -> Vec<[f32; 3]> {
    if group.len() < 2 {
        return Vec::new();
    }
    let mut out = Vec::with_capacity(group.len());
    let start_independent =
        cartesian_axis_value(Vec3::from_array(group[0]), ind_var).unwrap_or(0.0);
    let mut accum = 0.0_f32;
    out.push(cartesian_line_point(dep_var, ind_var, start_independent, accum).to_array());
    for pair in group.windows(2) {
        let a = Vec3::from_array(pair[0]);
        let b = Vec3::from_array(pair[1]);
        let ia = match cartesian_axis_value(a, ind_var) {
            Some(value) => value,
            None => continue,
        };
        let ib = match cartesian_axis_value(b, ind_var) {
            Some(value) => value,
            None => continue,
        };
        let da = match cartesian_axis_value(a, dep_var) {
            Some(value) => value,
            None => continue,
        };
        let db = match cartesian_axis_value(b, dep_var) {
            Some(value) => value,
            None => continue,
        };
        accum += (da + db) * 0.5 * (ib - ia);
        out.push(cartesian_line_point(dep_var, ind_var, ib, accum).to_array());
    }
    if normalize_integral {
        normalize_cartesian_line_group(&mut out, dep_var);
    }
    out
}

fn normal_curve_group(group: &[[f32; 3]]) -> Vec<[f32; 3]> {
    if group.len() < 3 {
        return Vec::new();
    }
    let tangents = tangent_vectors(group);
    tangents
        .iter()
        .enumerate()
        .map(|(index, _)| {
            let dt = if index == 0 {
                tangents[1] - tangents[0]
            } else if index + 1 == tangents.len() {
                tangents[index] - tangents[index - 1]
            } else {
                (tangents[index + 1] - tangents[index - 1]) * 0.5
            };
            dt.normalize_or_zero().to_array()
        })
        .collect()
}

fn binormal_curve_group(group: &[[f32; 3]]) -> Vec<[f32; 3]> {
    if group.len() < 3 {
        return Vec::new();
    }
    let tangents = tangent_vectors(group);
    let normals = normal_curve_group(group);
    tangents
        .iter()
        .zip(normals.iter())
        .map(|(tangent, normal)| {
            tangent
                .cross(Vec3::from_array(*normal))
                .normalize_or_zero()
                .to_array()
        })
        .collect()
}

fn tangent_vectors(group: &[[f32; 3]]) -> Vec<Vec3> {
    if group.len() < 2 {
        return Vec::new();
    }
    (0..group.len())
        .map(|index| finite_difference(group, index).normalize_or_zero())
        .collect()
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::{Domain, PlotDefinition, PlotSpec, PlotStyle, Resolution};
    use viewport_lib::MeshData;

    fn approx_eq(a: f32, b: f32) -> bool {
        (a - b).abs() <= 1.0e-3
    }

    fn point_annotations_plot(points: &[[f32; 3]]) -> PlotSpec {
        PlotSpec {
            name: "Points".to_string(),
            visible: true,
            domain: Domain::default(),
            resolution: Resolution::default(),
            style: PlotStyle::default(),
            definition: PlotDefinition::PointAnnotations {
                points: points
                    .iter()
                    .enumerate()
                    .map(|(index, position)| PointAnnotation {
                        position: *position,
                        label: format!("P{}", index + 1),
                    })
                    .collect(),
                show_labels: false,
            },
        }
    }

    fn square_surface_plot() -> PlotSpec {
        PlotSpec {
            name: "Surface".to_string(),
            visible: true,
            domain: Domain::default(),
            resolution: Resolution::default(),
            style: PlotStyle::default(),
            definition: PlotDefinition::GridSurface,
        }
    }

    fn square_surface_mesh() -> MeshData {
        let mut mesh = MeshData::default();
        mesh.positions = vec![
            [0.0, 0.0, 0.0],
            [1.0, 0.0, 0.0],
            [1.0, 1.0, 0.0],
            [0.0, 1.0, 0.0],
        ];
        mesh.indices = vec![0, 1, 2, 0, 2, 3];
        mesh
    }

    #[test]
    fn derivative_of_planar_xz_curve_stays_in_xz_plane() {
        let group = vec![[0.0, 4.0, 0.0], [1.0, 4.0, 1.0], [2.0, 4.0, 0.0]];

        let derived = derivative_curve_group(&group);

        assert_eq!(derived.len(), 3);
        assert!(approx_eq(derived[1][0], 1.0));
        assert!(approx_eq(derived[1][1], 4.0));
        assert!(approx_eq(derived[1][2], 0.0));
    }

    #[test]
    fn axis_derivative_preserves_constant_plane_axis() {
        let group = vec![[0.0, 4.0, 0.0], [1.0, 4.0, 1.0], [2.0, 4.0, 0.0]];

        let derived = axis_derivative_group(&group, 2, 0);

        assert_eq!(derived.len(), 3);
        assert!(approx_eq(derived[1][0], 1.0));
        assert!(approx_eq(derived[1][1], 4.0));
        assert!(approx_eq(derived[1][2], 0.0));
    }

    #[test]
    fn integral_of_planar_xz_curve_stays_in_xz_plane() {
        let group = vec![[0.0, 4.0, 0.0], [1.0, 4.0, 1.0], [2.0, 4.0, 0.0]];

        let derived = integral_curve_group(&group, false);

        assert_eq!(derived.len(), 3);
        assert!(approx_eq(derived[0][1], 4.0));
        assert!(approx_eq(derived[1][1], 4.0));
        assert!(approx_eq(derived[2][1], 4.0));
        assert!(approx_eq(derived[2][2], 1.0));
    }

    #[test]
    fn repeated_integral_of_sine_stays_centered_when_normalized() {
        let group = (-80..=80)
            .map(|i| {
                let x = i as f32 * 0.1;
                [x, 4.0, x.sin()]
            })
            .collect::<Vec<_>>();

        let first = integral_curve_group(&group, true);
        let second = integral_curve_group(&first, true);

        assert_eq!(second.len(), group.len());
        let mean = second.iter().map(|point| point[2]).sum::<f32>() / second.len() as f32;
        let left = second.first().map(|point| point[2]).unwrap_or_default();
        let right = second.last().map(|point| point[2]).unwrap_or_default();

        assert!(mean.abs() < 1.0e-3);
        assert!((left - right).abs() < 0.2);
    }

    #[test]
    fn phase_three_quality_output_titles_and_links_flagged_samples() {
        let plot = point_annotations_plot(&[
            [0.0, 0.0, 0.0],
            [0.0, 0.0, 0.0],
            [0.0001, 0.0, 0.0],
            [2.0, 0.0, 0.0],
            [20.0, 0.0, 0.0],
        ]);

        let output = make_data_quality_output(&plot).expect("quality output");
        let AnalysisOutput::Composite {
            reports: _,
            tables,
            diagnostics,
            ..
        } = output
        else {
            panic!("expected composite output");
        };

        let table_titles = tables
            .iter()
            .map(|table| table.title.as_str())
            .collect::<Vec<_>>();
        assert!(table_titles.contains(&"Sample Positions"));
        assert!(table_titles.contains(&"Exact Duplicates"));
        assert!(table_titles.contains(&"Near Duplicates"));

        let linked_messages = diagnostics
            .iter()
            .filter_map(|diagnostic| {
                diagnostic.location.as_ref().map(|location| {
                    (
                        diagnostic.message.as_str(),
                        location.component.as_deref(),
                        location.row,
                    )
                })
            })
            .collect::<Vec<_>>();

        assert!(linked_messages.iter().any(|(message, component, row)| {
            message.contains("Exact duplicate")
                && *component == Some("sequence 1 sample")
                && *row == Some(2)
        }));
        assert!(linked_messages.iter().any(|(message, component, row)| {
            message.contains("Near-duplicate")
                && *component == Some("sequence 1 sample")
                && *row == Some(3)
        }));
    }

    #[test]
    fn surface_area_and_curvature_on_plane_are_stable() {
        let plot = square_surface_plot();
        let mesh = square_surface_mesh();

        let area = run_surface_mesh_analysis(&plot, AnalysisKind::SurfaceArea, &[&mesh], &[])
            .expect("surface area");
        let AnalysisOutput::Report { report, .. } = area else {
            panic!("expected report output");
        };
        assert!(
            report
                .values
                .iter()
                .any(|(label, value)| { label == "Surface Area" && value == "1.00000" })
        );

        let curvature =
            run_surface_mesh_analysis(&plot, AnalysisKind::SurfaceCurvature, &[&mesh], &[])
                .expect("surface curvature");
        let AnalysisOutput::Composite { tables, .. } = curvature else {
            panic!("expected composite output");
        };
        let sample_table = tables
            .iter()
            .find(|table| table.title == "Vertex Curvature Samples")
            .expect("curvature sample table");
        assert!(
            sample_table
                .rows
                .iter()
                .all(|row| { row[4].parse::<f32>().unwrap_or(1.0).abs() < 1.0e-3 })
        );
    }

    #[test]
    fn frenet_frame_output_emits_tables_and_axis_plots() {
        let plot = PlotSpec {
            name: "Curve".to_string(),
            visible: true,
            domain: Domain::default(),
            resolution: Resolution { u: 32, v: 1 },
            style: PlotStyle::default(),
            definition: PlotDefinition::HelixCurve,
        };

        let request = AnalysisRequest {
            kind: AnalysisKind::FrenetFrame,
            target: AnalysisTarget::Plot {
                index: 0,
                name: Some(plot.name.clone()),
            },
            parameters: vec![
                ("max_samples".to_string(), "24".to_string()),
                ("vector_scale".to_string(), "0.5".to_string()),
            ],
        };
        let output = run_analysis(&plot, &request).expect("frenet frame output");
        let AnalysisOutput::Composite {
            plots,
            tables,
            frame_fields,
            ..
        } = output
        else {
            panic!("expected composite output");
        };

        assert!(plots.is_empty());
        assert_eq!(frame_fields.len(), 1);
        assert_eq!(frame_fields[0].frame_kind, AnalysisKind::FrenetFrame);
        assert!(tables.iter().any(|table| table.title.contains("Samples")));
    }

    #[test]
    fn bishop_frames_stay_orthonormal_on_planar_curve() {
        let points = (0..8)
            .map(|index| {
                let t = index as f32 * 0.3;
                [t.cos(), t.sin(), 0.0]
            })
            .collect::<Vec<_>>();

        let field = build_curve_frame_field(
            &points,
            &point_annotations_plot(&points),
            AnalysisKind::BishopFrame,
            0,
            1,
        )
        .expect("bishop frame field");

        for sample in field.samples {
            let t = Vec3::from_array(sample.tangent);
            let n = Vec3::from_array(sample.normal);
            let b = Vec3::from_array(sample.binormal);
            assert!(approx_eq(t.length(), 1.0));
            assert!(approx_eq(n.length(), 1.0));
            assert!(approx_eq(b.length(), 1.0));
            assert!(t.dot(n).abs() < 1.0e-3);
            assert!(t.dot(b).abs() < 1.0e-3);
            assert!(n.dot(b).abs() < 1.0e-3);
        }
    }

    #[test]
    fn surface_aligned_frame_uses_surface_context() {
        let curve = PlotSpec {
            name: "Curve".to_string(),
            visible: true,
            domain: Domain::default(),
            resolution: Resolution::default(),
            style: PlotStyle::default(),
            definition: PlotDefinition::DerivedPolylineGroups {
                groups: vec![vec![[0.0, 0.0, 0.0], [0.5, 0.0, 0.0], [1.0, 0.0, 0.0]]],
            },
        };
        let surface = square_surface_mesh();
        let output = run_curve_surface_frame_analysis(
            &curve,
            "Surface",
            AnalysisKind::SurfaceAlignedFrame,
            &[&surface],
            &[("max_samples".to_string(), "8".to_string())],
        )
        .expect("surface aligned frame");
        let AnalysisOutput::Composite { frame_fields, .. } = output else {
            panic!("expected composite output");
        };
        let field = frame_fields.first().expect("frame field");
        assert_eq!(field.frame_kind, AnalysisKind::SurfaceAlignedFrame);
        let first = field.samples.first().expect("sample");
        assert!(first.normal[2].abs() > 0.9);
    }
}