runmat_plot/plots/

surface.rs

//! 3D surface plot implementation
//!
//! High-performance GPU-accelerated 3D surface rendering.

use crate::context::shared_wgpu_context;
use crate::core::{
    BoundingBox, DrawCall, GpuVertexBuffer, ImageData, Material, PipelineType, RenderData, Vertex,
};
use crate::gpu::axis::OwnedAxisData;
use crate::gpu::{util::readback_scalar_buffer_f64, ScalarType};
use glam::{Vec3, Vec4};
use std::fmt::Write;
use std::sync::Arc;

/// High-performance GPU-accelerated 3D surface plot
#[derive(Debug, Clone)]
pub struct SurfacePlot {
    /// Grid data (Z values at X,Y coordinates)
    pub x_data: Vec<f64>,
    pub y_data: Vec<f64>,
    pub z_data: Option<Vec<Vec<f64>>>, // Host data when available
    /// Optional full coordinate grids for parametric surfaces where X/Y are not separable axes.
    pub x_grid: Option<Vec<Vec<f64>>>,
    pub y_grid: Option<Vec<Vec<f64>>>,
    /// Grid resolution for rendering/index generation (kept even for GPU-backed plots).
    x_len: usize,
    y_len: usize,

    /// Surface properties
    pub colormap: ColorMap,
    pub shading_mode: ShadingMode,
    pub wireframe: bool,
    pub alpha: f32,
    /// If true, render Z at 0 (flat), but color-map using Z values
    pub flatten_z: bool,

    /// If true, this flattened surface should behave like a 2D image for camera/UI decisions.
    pub image_mode: bool,

    /// Optional color limits override for mapping Z -> color (caxis)
    pub color_limits: Option<(f64, f64)>,

    /// Optional per-vertex color grid (for RGB images); if set, overrides colormap mapping
    pub color_grid: Option<Vec<Vec<Vec4>>>, // [x_index][y_index] -> RGBA

    /// Lighting and material
    pub lighting_enabled: bool,
    pub ambient_strength: f32,
    pub diffuse_strength: f32,
    pub specular_strength: f32,
    pub shininess: f32,

    /// Metadata
    pub label: Option<String>,
    pub visible: bool,

    /// Generated rendering data (cached)
    vertices: Option<Vec<Vertex>>,
    indices: Option<Vec<u32>>,
    bounds: Option<BoundingBox>,
    dirty: bool,
    gpu_vertices: Option<GpuVertexBuffer>,
    gpu_vertex_count: Option<usize>,
    gpu_bounds: Option<BoundingBox>,
    gpu_source: Option<SurfaceGpuSource>,
    gpu_color_grid_source: Option<SurfaceGpuColorGridSource>,
}

#[derive(Clone, Debug)]
pub struct SurfaceGpuSource {
    pub x_axis: OwnedAxisData,
    pub y_axis: OwnedAxisData,
    pub z_buffer: Arc<wgpu::Buffer>,
    pub x_len: usize,
    pub y_len: usize,
    pub scalar: ScalarType,
}

#[derive(Clone, Debug)]
pub struct SurfaceGpuColorGridSource {
    pub image_buffer: Arc<wgpu::Buffer>,
    pub rows: usize,
    pub cols: usize,
    pub channels: usize,
    pub scalar: ScalarType,
}

/// Color mapping schemes
#[derive(Clone, PartialEq)]
pub enum ColorMap {
    /// MATLAB-compatible colormaps
    Jet,
    Hot,
    Cool,
    Spring,
    Summer,
    Autumn,
    Winter,
    Gray,
    Bone,
    Copper,
    Pink,
    Lines,

    /// Scientific colormaps
    Viridis,
    Plasma,
    Inferno,
    Magma,
    Turbo,

    /// Perceptually uniform
    Parula,

    /// MATLAB colorcube colormap
    ColorCube,

    /// Custom color ranges
    Custom(Vec4, Vec4), // (min_color, max_color)

    /// Explicit RGB lookup table from an m-by-3 MATLAB colormap matrix
    Listed(Arc<[[f32; 3]]>),
}

impl std::fmt::Debug for ColorMap {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self {
            Self::Jet => f.write_str("Jet"),
            Self::Hot => f.write_str("Hot"),
            Self::Cool => f.write_str("Cool"),
            Self::Spring => f.write_str("Spring"),
            Self::Summer => f.write_str("Summer"),
            Self::Autumn => f.write_str("Autumn"),
            Self::Winter => f.write_str("Winter"),
            Self::Gray => f.write_str("Gray"),
            Self::Bone => f.write_str("Bone"),
            Self::Copper => f.write_str("Copper"),
            Self::Pink => f.write_str("Pink"),
            Self::Lines => f.write_str("Lines"),
            Self::Viridis => f.write_str("Viridis"),
            Self::Plasma => f.write_str("Plasma"),
            Self::Inferno => f.write_str("Inferno"),
            Self::Magma => f.write_str("Magma"),
            Self::Turbo => f.write_str("Turbo"),
            Self::Parula => f.write_str("Parula"),
            Self::ColorCube => f.write_str("ColorCube"),
            Self::Custom(min, max) => f.debug_tuple("Custom").field(min).field(max).finish(),
            Self::Listed(colors) => f.debug_tuple("Listed").field(&colors.len()).finish(),
        }
    }
}

impl ColorMap {
    pub const CANONICAL_NAMES: &[&str] = &[
        "parula",
        "colorcube",
        "viridis",
        "plasma",
        "inferno",
        "magma",
        "turbo",
        "jet",
        "hot",
        "cool",
        "spring",
        "summer",
        "autumn",
        "winter",
        "gray",
        "bone",
        "copper",
        "pink",
        "lines",
    ];

    pub const ALIASES: &[&str] = &["grey"];

    pub fn from_name(name: &str) -> Option<Self> {
        match name.trim().to_ascii_lowercase().as_str() {
            "parula" => Some(Self::Parula),
            "colorcube" => Some(Self::ColorCube),
            "viridis" => Some(Self::Viridis),
            "plasma" => Some(Self::Plasma),
            "inferno" => Some(Self::Inferno),
            "magma" => Some(Self::Magma),
            "turbo" => Some(Self::Turbo),
            "jet" => Some(Self::Jet),
            "hot" => Some(Self::Hot),
            "cool" => Some(Self::Cool),
            "spring" => Some(Self::Spring),
            "summer" => Some(Self::Summer),
            "autumn" => Some(Self::Autumn),
            "winter" => Some(Self::Winter),
            "gray" | "grey" => Some(Self::Gray),
            "bone" => Some(Self::Bone),
            "copper" => Some(Self::Copper),
            "pink" => Some(Self::Pink),
            "lines" => Some(Self::Lines),
            _ => None,
        }
    }

    pub fn from_rgb_rows(colors: Vec<[f32; 3]>) -> Option<Self> {
        if colors.is_empty() {
            return None;
        }
        Some(Self::Listed(Arc::from(colors.into_boxed_slice())))
    }

    pub fn to_serialized_token(&self) -> String {
        match self {
            Self::Listed(colors) => {
                let mut token = String::from("listed:");
                for (idx, color) in colors.iter().enumerate() {
                    if idx > 0 {
                        token.push(';');
                    }
                    let _ = write!(token, "{:?},{:?},{:?}", color[0], color[1], color[2]);
                }
                token
            }
            Self::Custom(min, max) => format!(
                "custom:{:?},{:?},{:?},{:?};{:?},{:?},{:?},{:?}",
                min.x, min.y, min.z, min.w, max.x, max.y, max.z, max.w
            ),
            _ => format!("{self:?}"),
        }
    }

    pub fn from_serialized_token(token: &str) -> Option<Self> {
        let trimmed = token.trim();
        if let Some(rest) = trimmed.strip_prefix("listed:") {
            let mut colors = Vec::new();
            for row in rest.split(';') {
                if row.is_empty() {
                    return None;
                }
                let mut parts = row.split(',');
                let r = parts.next()?.parse::<f32>().ok()?;
                let g = parts.next()?.parse::<f32>().ok()?;
                let b = parts.next()?.parse::<f32>().ok()?;
                if parts.next().is_some() {
                    return None;
                }
                if ![r, g, b]
                    .iter()
                    .all(|value| value.is_finite() && (0.0..=1.0).contains(value))
                {
                    return None;
                }
                colors.push([r, g, b]);
            }
            return Self::from_rgb_rows(colors);
        }
        if let Some(rest) = trimmed.strip_prefix("custom:") {
            let mut rows = rest.split(';');
            let min = parse_vec4_token(rows.next()?)?;
            let max = parse_vec4_token(rows.next()?)?;
            if rows.next().is_some() {
                return None;
            }
            return Some(Self::Custom(min, max));
        }
        Self::from_name(trimmed)
    }
}

fn parse_vec4_token(token: &str) -> Option<Vec4> {
    let mut parts = token.split(',');
    let x = parts.next()?.parse::<f32>().ok()?;
    let y = parts.next()?.parse::<f32>().ok()?;
    let z = parts.next()?.parse::<f32>().ok()?;
    let w = parts.next()?.parse::<f32>().ok()?;
    if parts.next().is_some() || ![x, y, z, w].iter().all(|value| value.is_finite()) {
        return None;
    }
    Some(Vec4::new(x, y, z, w))
}

/// Surface shading modes
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum ShadingMode {
    /// Flat shading (per-face normals)
    Flat,
    /// Smooth shading (interpolated normals)
    Smooth,
    /// Faceted (flat with visible edges)
    Faceted,
    /// No shading (just color mapping)
    None,
}

impl Default for ColorMap {
    fn default() -> Self {
        Self::Viridis
    }
}

impl Default for ShadingMode {
    fn default() -> Self {
        Self::Smooth
    }
}

impl SurfacePlot {
    pub async fn export_scene_grid_data(
        &self,
    ) -> Result<(Vec<f64>, Vec<f64>, Vec<Vec<f64>>), String> {
        if let Some(z) = &self.z_data {
            return Ok((self.x_data.clone(), self.y_data.clone(), z.clone()));
        }

        if let Some(source) = &self.gpu_source {
            let context = shared_wgpu_context().ok_or_else(|| {
                "surface plot has GPU source data but no shared WGPU context is installed"
                    .to_string()
            })?;
            let x = source
                .x_axis
                .export_f64(&context.device, &context.queue, source.x_len, source.scalar)
                .await?;
            let y = source
                .y_axis
                .export_f64(&context.device, &context.queue, source.y_len, source.scalar)
                .await?;
            let z_flat = readback_scalar_buffer_f64(
                &context.device,
                &context.queue,
                &source.z_buffer,
                source.x_len * source.y_len,
                source.scalar,
            )
            .await?;
            let mut z = Vec::with_capacity(source.x_len);
            for row in 0..source.x_len {
                let start = row * source.y_len;
                z.push(
                    z_flat
                        .get(start..start + source.y_len)
                        .ok_or_else(|| "surface GPU source grid is out of range".to_string())?
                        .to_vec(),
                );
            }
            return Ok((x, y, z));
        }

        if self.gpu_vertices.is_some() {
            return Err(
                "surface plot has GPU render vertices but no exportable source data".to_string(),
            );
        }

        Ok((Vec::new(), Vec::new(), Vec::new()))
    }

    pub async fn export_scene_color_grid(&self) -> Result<Option<Vec<Vec<Vec4>>>, String> {
        if let Some(grid) = &self.color_grid {
            return Ok(Some(grid.clone()));
        }

        let Some(source) = &self.gpu_color_grid_source else {
            return Ok(None);
        };
        let context = shared_wgpu_context().ok_or_else(|| {
            "surface image has GPU color data but no shared WGPU context is installed".to_string()
        })?;
        let values = readback_scalar_buffer_f64(
            &context.device,
            &context.queue,
            &source.image_buffer,
            source.rows * source.cols * source.channels,
            source.scalar,
        )
        .await?;
        let mut grid = vec![vec![Vec4::ZERO; source.rows]; source.cols];
        let plane = source.rows * source.cols;
        for (col, grid_row) in grid.iter_mut().enumerate() {
            for (row, color) in grid_row.iter_mut().enumerate() {
                let base = row + source.rows * col;
                let r = values.get(base).copied().unwrap_or(0.0) as f32;
                let g = values.get(base + plane).copied().unwrap_or(0.0) as f32;
                let b = values.get(base + (2 * plane)).copied().unwrap_or(0.0) as f32;
                let a = if source.channels == 4 {
                    values.get(base + (3 * plane)).copied().unwrap_or(1.0) as f32
                } else {
                    1.0
                };
                *color = Vec4::new(r, g, b, a);
            }
        }
        Ok(Some(grid))
    }

    /// Create a new surface plot from meshgrid data
    pub fn new(x_data: Vec<f64>, y_data: Vec<f64>, z_data: Vec<Vec<f64>>) -> Result<Self, String> {
        // Validate dimensions
        if z_data.len() != x_data.len() {
            return Err(format!(
                "Z data rows ({}) must match X data length ({})",
                z_data.len(),
                x_data.len()
            ));
        }

        for (i, row) in z_data.iter().enumerate() {
            if row.len() != y_data.len() {
                return Err(format!(
                    "Z data row {} length ({}) must match Y data length ({})",
                    i,
                    row.len(),
                    y_data.len()
                ));
            }
        }

        Ok(Self {
            x_len: x_data.len(),
            y_len: y_data.len(),
            x_data,
            y_data,
            z_data: Some(z_data),
            x_grid: None,
            y_grid: None,
            colormap: ColorMap::default(),
            shading_mode: ShadingMode::default(),
            wireframe: false,
            alpha: 1.0,
            flatten_z: false,
            image_mode: false,
            color_limits: None,
            color_grid: None,
            lighting_enabled: true,
            ambient_strength: 0.2,
            diffuse_strength: 0.8,
            specular_strength: 0.5,
            shininess: 32.0,
            label: None,
            visible: true,
            vertices: None,
            indices: None,
            bounds: None,
            dirty: true,
            gpu_vertices: None,
            gpu_vertex_count: None,
            gpu_bounds: None,
            gpu_source: None,
            gpu_color_grid_source: None,
        })
    }

    /// Create a new surface plot from full X/Y/Z coordinate grids.
    pub fn from_coordinate_grids(
        x_grid: Vec<Vec<f64>>,
        y_grid: Vec<Vec<f64>>,
        z_grid: Vec<Vec<f64>>,
    ) -> Result<Self, String> {
        validate_coordinate_grids(&x_grid, &y_grid, &z_grid)?;
        let x_len = z_grid.len();
        let y_len = z_grid.first().map_or(0, Vec::len);
        Ok(Self {
            x_data: (0..x_len).map(|i| i as f64 + 1.0).collect(),
            y_data: (0..y_len).map(|i| i as f64 + 1.0).collect(),
            z_data: Some(z_grid),
            x_grid: Some(x_grid),
            y_grid: Some(y_grid),
            x_len,
            y_len,
            colormap: ColorMap::default(),
            shading_mode: ShadingMode::default(),
            wireframe: false,
            alpha: 1.0,
            flatten_z: false,
            image_mode: false,
            color_limits: None,
            color_grid: None,
            lighting_enabled: true,
            ambient_strength: 0.2,
            diffuse_strength: 0.8,
            specular_strength: 0.5,
            shininess: 32.0,
            label: None,
            visible: true,
            vertices: None,
            indices: None,
            bounds: None,
            dirty: true,
            gpu_vertices: None,
            gpu_vertex_count: None,
            gpu_bounds: None,
            gpu_source: None,
            gpu_color_grid_source: None,
        })
    }

    /// Create a surface plot backed by a GPU vertex buffer.
    pub fn from_gpu_buffer(
        x_len: usize,
        y_len: usize,
        buffer: GpuVertexBuffer,
        vertex_count: usize,
        bounds: BoundingBox,
    ) -> Self {
        Self {
            x_data: Vec::new(),
            y_data: Vec::new(),
            z_data: None,
            x_grid: None,
            y_grid: None,
            x_len,
            y_len,
            colormap: ColorMap::default(),
            shading_mode: ShadingMode::default(),
            wireframe: false,
            alpha: 1.0,
            flatten_z: false,
            image_mode: false,
            color_limits: None,
            color_grid: None,
            lighting_enabled: true,
            ambient_strength: 0.2,
            diffuse_strength: 0.8,
            specular_strength: 0.5,
            shininess: 32.0,
            label: None,
            visible: true,
            vertices: None,
            indices: None,
            bounds: Some(bounds),
            dirty: false,
            gpu_vertices: Some(buffer),
            gpu_vertex_count: Some(vertex_count),
            gpu_bounds: Some(bounds),
            gpu_source: None,
            gpu_color_grid_source: None,
        }
    }

    pub fn with_gpu_source(mut self, source: SurfaceGpuSource) -> Self {
        self.gpu_source = Some(source);
        self
    }

    pub fn with_gpu_color_grid_source(mut self, source: SurfaceGpuColorGridSource) -> Self {
        self.gpu_color_grid_source = Some(source);
        self
    }

    pub fn update_axis_data(
        &mut self,
        x_data: Vec<f64>,
        y_data: Vec<f64>,
        z_data: Vec<Vec<f64>>,
    ) -> Result<(), String> {
        if z_data.len() != x_data.len() {
            return Err(format!(
                "Z data rows ({}) must match X data length ({})",
                z_data.len(),
                x_data.len()
            ));
        }
        for (idx, row) in z_data.iter().enumerate() {
            if row.len() != y_data.len() {
                return Err(format!(
                    "Z data row {idx} length ({}) must match Y data length ({})",
                    row.len(),
                    y_data.len()
                ));
            }
        }

        self.x_len = x_data.len();
        self.y_len = y_data.len();
        self.x_data = x_data;
        self.y_data = y_data;
        self.z_data = Some(z_data);
        self.x_grid = None;
        self.y_grid = None;
        self.reset_source_data();
        Ok(())
    }

    pub fn update_coordinate_grids(
        &mut self,
        x_grid: Vec<Vec<f64>>,
        y_grid: Vec<Vec<f64>>,
        z_grid: Vec<Vec<f64>>,
    ) -> Result<(), String> {
        validate_coordinate_grids(&x_grid, &y_grid, &z_grid)?;
        let x_len = z_grid.len();
        let y_len = z_grid.first().map_or(0, Vec::len);
        self.x_data = (0..x_len).map(|i| i as f64 + 1.0).collect();
        self.y_data = (0..y_len).map(|i| i as f64 + 1.0).collect();
        self.z_data = Some(z_grid);
        self.x_grid = Some(x_grid);
        self.y_grid = Some(y_grid);
        self.x_len = x_len;
        self.y_len = y_len;
        self.reset_source_data();
        Ok(())
    }

    fn drop_gpu_if_possible(&mut self) {
        if self.gpu_vertices.is_some() && self.z_data.is_some() {
            self.invalidate_gpu_data();
        }
    }

    fn reset_source_data(&mut self) {
        self.vertices = None;
        self.indices = None;
        self.bounds = None;
        self.gpu_color_grid_source = None;
        self.invalidate_gpu_data();
        if self.color_grid.as_ref().is_some_and(|grid| {
            grid.len() != self.x_len || grid.iter().any(|row| row.len() != self.y_len)
        }) {
            self.color_grid = None;
        }
        self.dirty = true;
    }

    /// Create surface from a function
    pub fn from_function<F>(
        x_range: (f64, f64),
        y_range: (f64, f64),
        resolution: (usize, usize),
        func: F,
    ) -> Result<Self, String>
    where
        F: Fn(f64, f64) -> f64,
    {
        let (x_res, y_res) = resolution;
        if x_res < 2 || y_res < 2 {
            return Err("Resolution must be at least 2x2".to_string());
        }

        let x_data: Vec<f64> = (0..x_res)
            .map(|i| x_range.0 + (x_range.1 - x_range.0) * i as f64 / (x_res - 1) as f64)
            .collect();

        let y_data: Vec<f64> = (0..y_res)
            .map(|j| y_range.0 + (y_range.1 - y_range.0) * j as f64 / (y_res - 1) as f64)
            .collect();

        let z_data: Vec<Vec<f64>> = x_data
            .iter()
            .map(|&x| y_data.iter().map(|&y| func(x, y)).collect())
            .collect();

        Self::new(x_data, y_data, z_data)
    }

    fn invalidate_gpu_data(&mut self) {
        self.gpu_vertices = None;
        self.gpu_vertex_count = None;
        self.gpu_bounds = None;
        self.gpu_source = None;
    }

    /// Set color mapping
    pub fn with_colormap(mut self, colormap: ColorMap) -> Self {
        self.colormap = colormap;
        self.dirty = true;
        self.drop_gpu_if_possible();
        self
    }

    /// Set shading mode
    pub fn with_shading(mut self, shading: ShadingMode) -> Self {
        self.shading_mode = shading;
        self.dirty = true;
        self.drop_gpu_if_possible();
        self
    }

    /// Enable/disable wireframe
    pub fn with_wireframe(mut self, enabled: bool) -> Self {
        self.wireframe = enabled;
        self.dirty = true;
        self.drop_gpu_if_possible();
        self
    }

    /// Set transparency
    pub fn with_alpha(mut self, alpha: f32) -> Self {
        self.alpha = alpha.clamp(0.0, 1.0);
        self.dirty = true;
        self.drop_gpu_if_possible();
        self
    }

    /// Render surface flat in Z while mapping colors from Z values (for imagesc/imshow)
    pub fn with_flatten_z(mut self, enabled: bool) -> Self {
        self.flatten_z = enabled;
        self.dirty = true;
        self.drop_gpu_if_possible();
        self
    }

    pub fn with_image_mode(mut self, enabled: bool) -> Self {
        self.image_mode = enabled;
        self.dirty = true;
        self.drop_gpu_if_possible();
        self
    }

    /// Override color mapping limits (caxis)
    pub fn with_color_limits(mut self, limits: Option<(f64, f64)>) -> Self {
        self.color_limits = limits;
        self.dirty = true;
        self.drop_gpu_if_possible();
        self
    }

    /// Mutably set color mapping limits (caxis)
    pub fn set_color_limits(&mut self, limits: Option<(f64, f64)>) {
        self.color_limits = limits;
        self.dirty = true;
        self.drop_gpu_if_possible();
    }

    /// Provide explicit per-vertex colors (RGB[A])
    pub fn with_color_grid(mut self, grid: Vec<Vec<Vec4>>) -> Self {
        self.color_grid = Some(grid);
        self.gpu_color_grid_source = None;
        self.dirty = true;
        self.drop_gpu_if_possible();
        self
    }

    /// Set plot label for legends
    pub fn with_label<S: Into<String>>(mut self, label: S) -> Self {
        self.label = Some(label.into());
        self
    }

    /// Get the number of grid points
    pub fn len(&self) -> usize {
        self.x_len * self.y_len
    }

    /// Check if the surface has no data
    pub fn is_empty(&self) -> bool {
        self.x_len == 0 || self.y_len == 0
    }

    /// Get the bounding box of the surface
    pub fn bounds(&mut self) -> BoundingBox {
        if self.dirty || self.bounds.is_none() {
            self.compute_bounds();
        }
        self.bounds.unwrap()
    }

    /// Compute bounding box
    fn compute_bounds(&mut self) {
        if let Some(bounds) = self.gpu_bounds {
            self.bounds = Some(bounds);
            return;
        }

        let mut min_x = f32::INFINITY;
        let mut max_x = f32::NEG_INFINITY;
        let mut min_y = f32::INFINITY;
        let mut max_y = f32::NEG_INFINITY;
        let mut min_z = f32::INFINITY;
        let mut max_z = f32::NEG_INFINITY;

        if let (Some(x_grid), Some(y_grid)) = (&self.x_grid, &self.y_grid) {
            for row in x_grid {
                for &x in row {
                    min_x = min_x.min(x as f32);
                    max_x = max_x.max(x as f32);
                }
            }
            for row in y_grid {
                for &y in row {
                    min_y = min_y.min(y as f32);
                    max_y = max_y.max(y as f32);
                }
            }
        } else {
            for &x in &self.x_data {
                min_x = min_x.min(x as f32);
                max_x = max_x.max(x as f32);
            }

            for &y in &self.y_data {
                min_y = min_y.min(y as f32);
                max_y = max_y.max(y as f32);
            }
        }

        if let Some(rows) = &self.z_data {
            for row in rows {
                for &z in row {
                    if z.is_finite() {
                        min_z = min_z.min(z as f32);
                        max_z = max_z.max(z as f32);
                    }
                }
            }
        }

        self.bounds = Some(BoundingBox::new(
            Vec3::new(min_x, min_y, min_z),
            Vec3::new(max_x, max_y, max_z),
        ));
    }

    /// Get plot statistics for debugging
    pub fn statistics(&self) -> SurfaceStatistics {
        let grid_size = self.x_len * self.y_len;
        let triangle_count = if self.x_len > 1 && self.y_len > 1 {
            (self.x_len - 1) * (self.y_len - 1) * 2
        } else {
            0
        };

        SurfaceStatistics {
            grid_points: grid_size,
            triangle_count,
            x_resolution: self.x_len,
            y_resolution: self.y_len,
            memory_usage: self.estimated_memory_usage(),
        }
    }

    /// Estimate memory usage in bytes
    pub fn estimated_memory_usage(&self) -> usize {
        let coordinate_grid_size = self
            .x_grid
            .as_ref()
            .zip(self.y_grid.as_ref())
            .map_or(0, |(x, y)| {
                x.iter().map(Vec::len).sum::<usize>() + y.iter().map(Vec::len).sum::<usize>()
            });
        let data_size = std::mem::size_of::<f64>()
            * (coordinate_grid_size
                + if coordinate_grid_size == 0 {
                    self.x_data.len() + self.y_data.len()
                } else {
                    0
                }
                + self
                    .z_data
                    .as_ref()
                    .map_or(0, |z| z.iter().map(Vec::len).sum()));

        let vertices_size = self
            .vertices
            .as_ref()
            .map_or(0, |v| v.len() * std::mem::size_of::<Vertex>());

        let indices_size = self
            .indices
            .as_ref()
            .map_or(0, |i| i.len() * std::mem::size_of::<u32>());

        let gpu_size = self.gpu_vertex_count.unwrap_or(0) * std::mem::size_of::<Vertex>();

        data_size + vertices_size + indices_size + gpu_size
    }

    /// Generate vertices for surface mesh
    fn generate_vertices(&mut self) -> &Vec<Vertex> {
        if self.gpu_vertices.is_some() {
            if self.vertices.is_none() {
                self.vertices = Some(Vec::new());
            }
            return self.vertices.as_ref().unwrap();
        }

        if self.dirty || self.vertices.is_none() {
            log::trace!(
                target: "runmat_plot",
                "surface gen vertices {} x {}",
                self.x_data.len(),
                self.y_data.len()
            );

            let mut vertices = Vec::new();

            // Determine color mapping range
            let z_rows = self
                .z_data
                .as_ref()
                .expect("CPU surface data missing during vertex generation");
            let (min_z, max_z) = if let Some((lo, hi)) = self.color_limits {
                (lo, hi)
            } else {
                let mut min_z = f64::INFINITY;
                let mut max_z = f64::NEG_INFINITY;
                for row in z_rows {
                    for &z in row {
                        if z.is_finite() {
                            min_z = min_z.min(z);
                            max_z = max_z.max(z);
                        }
                    }
                }
                (min_z, max_z)
            };
            let z_range = (max_z - min_z).max(f64::MIN_POSITIVE);

            // Generate vertices for each grid point
            for i in 0..self.x_len {
                for j in 0..self.y_len {
                    let (x, y) = self.coordinate_at(i, j);
                    let z = z_rows[i][j];
                    let z_pos = if self.flatten_z { 0.0 } else { z as f32 };
                    let position = Vec3::new(x as f32, y as f32, z_pos);

                    // Simple normal calculation (can be improved with proper gradients)
                    let normal = Vec3::new(0.0, 0.0, 1.0); // Placeholder

                    // Determine color: explicit grid (RGB) or colormap from Z
                    let color = if let Some(grid) = &self.color_grid {
                        let c = grid[i][j];
                        Vec4::new(c.x, c.y, c.z, c.w)
                    } else {
                        let t = ((z - min_z) / z_range) as f32;
                        let color_rgb = self.colormap.map_value(t.clamp(0.0, 1.0));
                        Vec4::new(color_rgb.x, color_rgb.y, color_rgb.z, self.alpha)
                    };

                    vertices.push(Vertex {
                        position: position.to_array(),
                        normal: normal.to_array(),
                        color: color.to_array(),
                        tex_coords: [
                            i as f32 / (self.x_len - 1).max(1) as f32,
                            j as f32 / (self.y_len - 1).max(1) as f32,
                        ],
                    });
                }
            }

            log::trace!(target: "runmat_plot", "surface vertices={}", vertices.len());
            self.vertices = Some(vertices);
        }
        self.vertices.as_ref().unwrap()
    }

    fn coordinate_at(&self, i: usize, j: usize) -> (f64, f64) {
        match (&self.x_grid, &self.y_grid) {
            (Some(x_grid), Some(y_grid)) => (x_grid[i][j], y_grid[i][j]),
            _ => (self.x_data[i], self.y_data[j]),
        }
    }

    /// Generate indices for surface triangulation
    fn generate_indices(&mut self) -> &Vec<u32> {
        if self.dirty || self.indices.is_none() {
            log::trace!(target: "runmat_plot", "surface generating indices");

            let mut indices = Vec::new();
            let x_res = self.x_len;
            let y_res = self.y_len;

            // Generate triangle indices for surface mesh
            for i in 0..x_res - 1 {
                for j in 0..y_res - 1 {
                    let base = (i * y_res + j) as u32;
                    let next_row = base + y_res as u32;

                    // Two triangles per quad
                    // Triangle 1: (i,j), (i+1,j), (i,j+1)
                    indices.push(base);
                    indices.push(next_row);
                    indices.push(base + 1);

                    // Triangle 2: (i+1,j), (i+1,j+1), (i,j+1)
                    indices.push(next_row);
                    indices.push(next_row + 1);
                    indices.push(base + 1);
                }
            }

            log::trace!(target: "runmat_plot", "surface indices={}", indices.len());
            self.indices = Some(indices);
            self.dirty = false;
        }
        self.indices.as_ref().unwrap()
    }

    fn generate_wireframe_indices(&self) -> Vec<u32> {
        let mut indices = Vec::new();
        if self.x_len < 2 || self.y_len < 2 {
            return indices;
        }

        // Horizontal grid edges (along Y for each X row)
        for i in 0..self.x_len {
            for j in 0..(self.y_len - 1) {
                let a = (i * self.y_len + j) as u32;
                let b = (i * self.y_len + j + 1) as u32;
                indices.push(a);
                indices.push(b);
            }
        }

        // Vertical grid edges (along X for each Y column)
        for i in 0..(self.x_len - 1) {
            for j in 0..self.y_len {
                let a = (i * self.y_len + j) as u32;
                let b = ((i + 1) * self.y_len + j) as u32;
                indices.push(a);
                indices.push(b);
            }
        }

        indices
    }

    /// Generate complete render data for the graphics pipeline
    pub fn render_data(&mut self) -> RenderData {
        log::debug!(
            target: "runmat_plot",
            "surface render_data start: {} x {}",
            self.x_len,
            self.y_len
        );

        if self.image_mode && self.z_data.is_some() && self.gpu_vertices.is_none() {
            return self.image_render_data();
        }

        let using_gpu = self.gpu_vertices.is_some();
        let bounds = self.bounds();
        let vertices = if using_gpu {
            Vec::new()
        } else {
            self.generate_vertices().clone()
        };
        let indices = if self.wireframe {
            self.generate_wireframe_indices()
        } else {
            self.generate_indices().clone()
        };

        let material = Material {
            albedo: Vec4::new(1.0, 1.0, 1.0, self.alpha),
            ..Default::default()
        };

        let vertex_count = if using_gpu {
            self.gpu_vertex_count.unwrap_or(0)
        } else {
            vertices.len()
        };

        log::debug!(
            target: "runmat_plot",
            "surface render_data generated: vertex_count={} (gpu={}), indices={}",
            vertex_count,
            using_gpu,
            indices.len()
        );

        let draw_call = DrawCall {
            vertex_offset: 0,
            vertex_count,
            index_offset: Some(0),
            index_count: Some(indices.len()),
            instance_count: 1,
        };

        log::trace!(target: "runmat_plot", "surface render_data done");

        RenderData {
            pipeline_type: if self.wireframe {
                PipelineType::Lines
            } else {
                PipelineType::Triangles
            },
            vertices,
            indices: Some(indices),

            gpu_vertices: self.gpu_vertices.clone(),
            bounds: Some(bounds),
            material,
            draw_calls: vec![draw_call],
            image: None,
        }
    }

    fn image_render_data(&mut self) -> RenderData {
        let bounds = self.bounds();
        let x_min = bounds.min.x;
        let x_max = bounds.max.x;
        let y_min = bounds.min.y;
        let y_max = bounds.max.y;
        let z_rows = self
            .z_data
            .as_ref()
            .expect("image-mode surfaces require host color data");
        let width = self.x_len.max(1);
        let height = self.y_len.max(1);
        let color_limits = self.color_limits.or_else(|| {
            let mut min_z = f64::INFINITY;
            let mut max_z = f64::NEG_INFINITY;
            for row in z_rows {
                for &z in row {
                    if z.is_finite() {
                        min_z = min_z.min(z);
                        max_z = max_z.max(z);
                    }
                }
            }
            if min_z.is_finite() && max_z.is_finite() {
                Some((min_z, max_z))
            } else {
                None
            }
        });
        let (min_z, max_z) = color_limits.unwrap_or((0.0, 1.0));
        let z_range = (max_z - min_z).max(f64::MIN_POSITIVE);
        let mut data = Vec::with_capacity(width * height * 4);

        for row in 0..height {
            let y_idx = height - 1 - row;
            for x_idx in 0..width {
                let color = if let Some(grid) = &self.color_grid {
                    grid[x_idx][y_idx]
                } else {
                    let z = z_rows[x_idx][y_idx];
                    let t = ((z - min_z) / z_range) as f32;
                    let rgb = self.colormap.map_value(t.clamp(0.0, 1.0));
                    Vec4::new(rgb.x, rgb.y, rgb.z, self.alpha)
                };
                data.push((color.x.clamp(0.0, 1.0) * 255.0).round() as u8);
                data.push((color.y.clamp(0.0, 1.0) * 255.0).round() as u8);
                data.push((color.z.clamp(0.0, 1.0) * 255.0).round() as u8);
                data.push((color.w.clamp(0.0, 1.0) * 255.0).round() as u8);
            }
        }

        let vertices = vec![
            Vertex {
                position: [x_min, y_min, 0.0],
                normal: [0.0, 0.0, 1.0],
                color: [1.0, 1.0, 1.0, self.alpha],
                tex_coords: [0.0, 1.0],
            },
            Vertex {
                position: [x_max, y_min, 0.0],
                normal: [0.0, 0.0, 1.0],
                color: [1.0, 1.0, 1.0, self.alpha],
                tex_coords: [1.0, 1.0],
            },
            Vertex {
                position: [x_max, y_max, 0.0],
                normal: [0.0, 0.0, 1.0],
                color: [1.0, 1.0, 1.0, self.alpha],
                tex_coords: [1.0, 0.0],
            },
            Vertex {
                position: [x_min, y_max, 0.0],
                normal: [0.0, 0.0, 1.0],
                color: [1.0, 1.0, 1.0, self.alpha],
                tex_coords: [0.0, 0.0],
            },
        ];
        let indices = vec![0, 1, 2, 0, 2, 3];

        RenderData {
            pipeline_type: PipelineType::Textured,
            vertices,
            indices: Some(indices.clone()),
            gpu_vertices: None,
            bounds: Some(bounds),
            material: Material {
                albedo: Vec4::new(1.0, 1.0, 1.0, self.alpha),
                ..Default::default()
            },
            draw_calls: vec![DrawCall {
                vertex_offset: 0,
                vertex_count: 4,
                index_offset: Some(0),
                index_count: Some(indices.len()),
                instance_count: 1,
            }],
            image: Some(ImageData::Rgba8 {
                width: width as u32,
                height: height as u32,
                data,
            }),
        }
    }
}

fn validate_coordinate_grids(
    x_grid: &[Vec<f64>],
    y_grid: &[Vec<f64>],
    z_grid: &[Vec<f64>],
) -> Result<(), String> {
    if z_grid.is_empty() {
        return Err("Z coordinate grid must not be empty".to_string());
    }
    if x_grid.len() != z_grid.len() || y_grid.len() != z_grid.len() {
        return Err("X, Y, and Z coordinate grids must have the same row count".to_string());
    }
    let y_len = z_grid[0].len();
    if y_len == 0 {
        return Err("Z coordinate grid rows must not be empty".to_string());
    }
    for (idx, ((x_row, y_row), z_row)) in x_grid
        .iter()
        .zip(y_grid.iter())
        .zip(z_grid.iter())
        .enumerate()
    {
        if x_row.len() != y_len || y_row.len() != y_len || z_row.len() != y_len {
            return Err(format!(
                "X, Y, and Z coordinate grid row {idx} must all have length {y_len}"
            ));
        }
        if x_row
            .iter()
            .chain(y_row.iter())
            .any(|value| !value.is_finite())
        {
            return Err("X and Y coordinate grids must contain finite values".to_string());
        }
    }
    Ok(())
}

/// Surface plot performance and data statistics
#[derive(Debug, Clone)]
pub struct SurfaceStatistics {
    pub grid_points: usize,
    pub triangle_count: usize,
    pub x_resolution: usize,
    pub y_resolution: usize,
    pub memory_usage: usize,
}

impl ColorMap {
    /// Map a normalized value [0,1] to a color
    pub fn map_value(&self, t: f32) -> Vec3 {
        let t = t.clamp(0.0, 1.0);

        match self {
            ColorMap::Jet => self.jet_colormap(t),
            ColorMap::Hot => self.hot_colormap(t),
            ColorMap::Cool => self.cool_colormap(t),
            ColorMap::Spring => self.spring_colormap(t),
            ColorMap::Summer => self.summer_colormap(t),
            ColorMap::Autumn => self.autumn_colormap(t),
            ColorMap::Winter => self.winter_colormap(t),
            ColorMap::Gray => Vec3::splat(t),
            ColorMap::Bone => self.bone_colormap(t),
            ColorMap::Copper => self.copper_colormap(t),
            ColorMap::Pink => self.pink_colormap(t),
            ColorMap::Lines => self.lines_colormap(t),
            ColorMap::Viridis => self.viridis_colormap(t),
            ColorMap::Plasma => self.plasma_colormap(t),
            ColorMap::Inferno => self.inferno_colormap(t),
            ColorMap::Magma => self.magma_colormap(t),
            ColorMap::Turbo => self.turbo_colormap(t),
            ColorMap::Parula => self.parula_colormap(t),
            ColorMap::ColorCube => self.colorcube_colormap(t),
            ColorMap::Custom(min_color, max_color) => {
                min_color.truncate().lerp(max_color.truncate(), t)
            }
            ColorMap::Listed(colors) => {
                let last = colors.len().saturating_sub(1);
                let index = if t >= 1.0 {
                    last
                } else {
                    (t.clamp(0.0, 1.0) * colors.len() as f32).floor() as usize
                };
                let color = colors[index.min(last)];
                Vec3::new(color[0], color[1], color[2])
            }
        }
    }

    /// MATLAB Jet colormap
    fn jet_colormap(&self, t: f32) -> Vec3 {
        let r = (1.5 - 4.0 * (t - 0.75).abs()).clamp(0.0, 1.0);
        let g = (1.5 - 4.0 * (t - 0.5).abs()).clamp(0.0, 1.0);
        let b = (1.5 - 4.0 * (t - 0.25).abs()).clamp(0.0, 1.0);
        Vec3::new(r, g, b)
    }

    /// Hot colormap (black -> red -> yellow -> white)
    fn hot_colormap(&self, t: f32) -> Vec3 {
        if t < 1.0 / 3.0 {
            Vec3::new(3.0 * t, 0.0, 0.0)
        } else if t < 2.0 / 3.0 {
            Vec3::new(1.0, 3.0 * t - 1.0, 0.0)
        } else {
            Vec3::new(1.0, 1.0, 3.0 * t - 2.0)
        }
    }

    /// Cool colormap (cyan -> magenta)
    fn cool_colormap(&self, t: f32) -> Vec3 {
        Vec3::new(t, 1.0 - t, 1.0)
    }

    /// Viridis colormap (perceptually uniform)
    fn viridis_colormap(&self, t: f32) -> Vec3 {
        // Simplified Viridis approximation
        let r = (0.267004 + t * (0.993248 - 0.267004)).clamp(0.0, 1.0);
        let g = (0.004874 + t * (0.906157 - 0.004874)).clamp(0.0, 1.0);
        let b = (0.329415 + t * (0.143936 - 0.329415) + t * t * 0.5).clamp(0.0, 1.0);
        Vec3::new(r, g, b)
    }

    /// Plasma colormap (perceptually uniform)
    fn plasma_colormap(&self, t: f32) -> Vec3 {
        // Simplified Plasma approximation
        let r = (0.050383 + t * (0.940015 - 0.050383)).clamp(0.0, 1.0);
        let g = (0.029803 + t * (0.975158 - 0.029803) * (1.0 - t)).clamp(0.0, 1.0);
        let b = (0.527975 + t * (0.131326 - 0.527975)).clamp(0.0, 1.0);
        Vec3::new(r, g, b)
    }

    /// Spring colormap (magenta -> yellow)
    fn spring_colormap(&self, t: f32) -> Vec3 {
        Vec3::new(1.0, t, 1.0 - t)
    }

    /// Summer colormap (green -> yellow)
    fn summer_colormap(&self, t: f32) -> Vec3 {
        Vec3::new(t, 0.5 + 0.5 * t, 0.4)
    }

    /// Autumn colormap (red -> yellow)
    fn autumn_colormap(&self, t: f32) -> Vec3 {
        Vec3::new(1.0, t, 0.0)
    }

    /// Winter colormap (blue -> green)
    fn winter_colormap(&self, t: f32) -> Vec3 {
        Vec3::new(0.0, t, 1.0 - 0.5 * t)
    }

    /// Bone colormap (black -> white with blue tint)
    fn bone_colormap(&self, t: f32) -> Vec3 {
        if t < 3.0 / 8.0 {
            Vec3::new(7.0 / 8.0 * t, 7.0 / 8.0 * t, 29.0 / 24.0 * t)
        } else {
            Vec3::new(
                (29.0 + 7.0 * t) / 24.0,
                (29.0 + 7.0 * t) / 24.0,
                (29.0 + 7.0 * t) / 24.0,
            )
        }
    }

    /// Copper colormap (black -> copper)
    fn copper_colormap(&self, t: f32) -> Vec3 {
        Vec3::new((1.25 * t).min(1.0), 0.7812 * t, 0.4975 * t)
    }

    /// Pink colormap (black -> pink -> white)
    fn pink_colormap(&self, t: f32) -> Vec3 {
        let sqrt_t = t.sqrt();
        if t < 3.0 / 8.0 {
            Vec3::new(14.0 / 9.0 * sqrt_t, 2.0 / 3.0 * sqrt_t, 2.0 / 3.0 * sqrt_t)
        } else {
            Vec3::new(
                2.0 * sqrt_t - 1.0 / 3.0,
                8.0 / 9.0 * sqrt_t + 1.0 / 3.0,
                8.0 / 9.0 * sqrt_t + 1.0 / 3.0,
            )
        }
    }

    /// Lines colormap (cycling through basic colors)
    fn lines_colormap(&self, t: f32) -> Vec3 {
        let _phase = (t * 7.0) % 1.0; // For future use in color transitions
        let index = (t * 7.0) as usize % 7;
        match index {
            0 => Vec3::new(0.0, 0.0, 1.0),    // Blue
            1 => Vec3::new(0.0, 0.5, 0.0),    // Green
            2 => Vec3::new(1.0, 0.0, 0.0),    // Red
            3 => Vec3::new(0.0, 0.75, 0.75),  // Cyan
            4 => Vec3::new(0.75, 0.0, 0.75),  // Magenta
            5 => Vec3::new(0.75, 0.75, 0.0),  // Yellow
            _ => Vec3::new(0.25, 0.25, 0.25), // Dark gray
        }
    }

    /// Inferno colormap (perceptually uniform)
    fn inferno_colormap(&self, t: f32) -> Vec3 {
        // Simplified Inferno approximation
        let r = (0.001462 + t * (0.988362 - 0.001462)).clamp(0.0, 1.0);
        let g = (0.000466 + t * t * (0.982895 - 0.000466)).clamp(0.0, 1.0);
        let b = (0.013866 + t * (1.0 - t) * (0.416065 - 0.013866)).clamp(0.0, 1.0);
        Vec3::new(r, g, b)
    }

    /// Magma colormap (perceptually uniform)
    fn magma_colormap(&self, t: f32) -> Vec3 {
        // Simplified Magma approximation
        let r = (0.001462 + t * (0.987053 - 0.001462)).clamp(0.0, 1.0);
        let g = (0.000466 + t * t * (0.991438 - 0.000466)).clamp(0.0, 1.0);
        let b = (0.013866 + t * (0.644237 - 0.013866) * (1.0 - t)).clamp(0.0, 1.0);
        Vec3::new(r, g, b)
    }

    /// Turbo colormap (improved rainbow)
    fn turbo_colormap(&self, t: f32) -> Vec3 {
        // Simplified Turbo approximation (Google's improved rainbow)
        let r = if t < 0.5 {
            (0.13 + 0.87 * (2.0 * t).powf(0.25)).clamp(0.0, 1.0)
        } else {
            (0.8685 + 0.1315 * (2.0 * (1.0 - t)).powf(0.25)).clamp(0.0, 1.0)
        };

        let g = if t < 0.25 {
            4.0 * t
        } else if t < 0.75 {
            1.0
        } else {
            1.0 - 4.0 * (t - 0.75)
        }
        .clamp(0.0, 1.0);

        let b = if t < 0.5 {
            (0.8 * (1.0 - 2.0 * t).powf(0.25)).clamp(0.0, 1.0)
        } else {
            (0.1 + 0.9 * (2.0 * t - 1.0).powf(0.25)).clamp(0.0, 1.0)
        };

        Vec3::new(r, g, b)
    }

    /// Parula colormap (MATLAB's default)
    fn parula_colormap(&self, t: f32) -> Vec3 {
        // Simplified Parula approximation
        let r = if t < 0.25 {
            0.2081 * (1.0 - t)
        } else if t < 0.5 {
            t - 0.25
        } else if t < 0.75 {
            1.0
        } else {
            1.0 - 0.5 * (t - 0.75)
        }
        .clamp(0.0, 1.0);

        let g = if t < 0.125 {
            0.1663 * t / 0.125
        } else if t < 0.375 {
            0.1663 + (0.7079 - 0.1663) * (t - 0.125) / 0.25
        } else if t < 0.625 {
            0.7079 + (0.9839 - 0.7079) * (t - 0.375) / 0.25
        } else {
            0.9839 * (1.0 - (t - 0.625) / 0.375)
        }
        .clamp(0.0, 1.0);

        let b = if t < 0.25 {
            0.5 + 0.5 * t / 0.25
        } else if t < 0.5 {
            1.0
        } else {
            1.0 - 2.0 * (t - 0.5)
        }
        .clamp(0.0, 1.0);

        Vec3::new(r, g, b)
    }

    /// Colorcube colormap.
    fn colorcube_colormap(&self, t: f32) -> Vec3 {
        if t >= 1.0 {
            return Vec3::ONE;
        }
        let levels = 6.0_f32;
        let index = (t.clamp(0.0, 1.0) * levels.powi(3)).floor() as u32;
        let r = (index % 6) as f32 / 5.0;
        let g = ((index / 6) % 6) as f32 / 5.0;
        let b = ((index / 36) % 6) as f32 / 5.0;
        Vec3::new(r, g, b)
    }

    /// Default colormap fallback
    #[allow(dead_code)] // Fallback method for colormap errors
    fn default_colormap(&self, t: f32) -> Vec3 {
        // Use a simple RGB transition as fallback
        if t < 0.5 {
            Vec3::new(0.0, 2.0 * t, 1.0 - 2.0 * t)
        } else {
            Vec3::new(2.0 * (t - 0.5), 1.0 - 2.0 * (t - 0.5), 0.0)
        }
    }
}

/// MATLAB-compatible surface plot creation utilities
pub mod matlab_compat {
    use super::*;

    /// Create a surface plot (equivalent to MATLAB's `surf(X, Y, Z)`)
    pub fn surf(x: Vec<f64>, y: Vec<f64>, z: Vec<Vec<f64>>) -> Result<SurfacePlot, String> {
        SurfacePlot::new(x, y, z)
    }

    /// Create a mesh plot (wireframe surface)
    pub fn mesh(x: Vec<f64>, y: Vec<f64>, z: Vec<Vec<f64>>) -> Result<SurfacePlot, String> {
        Ok(SurfacePlot::new(x, y, z)?
            .with_wireframe(true)
            .with_shading(ShadingMode::None))
    }

    /// Create surface from meshgrid
    pub fn meshgrid_surf(
        x_range: (f64, f64),
        y_range: (f64, f64),
        resolution: (usize, usize),
        func: impl Fn(f64, f64) -> f64,
    ) -> Result<SurfacePlot, String> {
        SurfacePlot::from_function(x_range, y_range, resolution, func)
    }

    /// Create surface with specific colormap
    pub fn surf_with_colormap(
        x: Vec<f64>,
        y: Vec<f64>,
        z: Vec<Vec<f64>>,
        colormap: &str,
    ) -> Result<SurfacePlot, String> {
        let cmap =
            ColorMap::from_name(colormap).ok_or_else(|| format!("Unknown colormap: {colormap}"))?;

        Ok(SurfacePlot::new(x, y, z)?.with_colormap(cmap))
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_surface_plot_creation() {
        let x = vec![0.0, 1.0, 2.0];
        let y = vec![0.0, 1.0];
        let z = vec![vec![0.0, 1.0], vec![1.0, 2.0], vec![2.0, 3.0]];

        let surface = SurfacePlot::new(x, y, z).unwrap();

        assert_eq!(surface.x_data.len(), 3);
        assert_eq!(surface.y_data.len(), 2);
        let rows = surface.z_data.as_ref().unwrap();
        assert_eq!(rows.len(), 3);
        assert_eq!(rows[0].len(), 2);
        assert!(surface.visible);
    }

    #[test]
    fn test_surface_from_function() {
        let surface =
            SurfacePlot::from_function((-2.0, 2.0), (-2.0, 2.0), (10, 10), |x, y| x * x + y * y)
                .unwrap();

        assert_eq!(surface.x_data.len(), 10);
        assert_eq!(surface.y_data.len(), 10);
        let rows = surface.z_data.as_ref().unwrap();
        assert_eq!(rows.len(), 10);

        // Check that function is evaluated correctly
        assert_eq!(rows[0][0], 8.0); // (-2)^2 + (-2)^2 = 8
    }

    #[test]
    fn test_surface_validation() {
        let x = vec![0.0, 1.0];
        let y = vec![0.0, 1.0, 2.0];
        let z = vec![
            vec![0.0, 1.0], // Wrong: should have 3 elements to match y
            vec![1.0, 2.0],
        ];

        assert!(SurfacePlot::new(x, y, z).is_err());
    }

    #[test]
    fn test_surface_styling() {
        let x = vec![0.0, 1.0];
        let y = vec![0.0, 1.0];
        let z = vec![vec![0.0, 1.0], vec![1.0, 2.0]];

        let surface = SurfacePlot::new(x, y, z)
            .unwrap()
            .with_colormap(ColorMap::Hot)
            .with_wireframe(true)
            .with_alpha(0.8)
            .with_label("Test Surface");

        assert_eq!(surface.colormap, ColorMap::Hot);
        assert!(surface.wireframe);
        assert_eq!(surface.alpha, 0.8);
        assert_eq!(surface.label, Some("Test Surface".to_string()));
    }

    #[test]
    fn image_mode_surface_uses_textured_render_data() {
        let x = vec![0.0, 1.0];
        let y = vec![10.0, 20.0];
        let z = vec![vec![0.0, 0.25], vec![0.75, 1.0]];

        let mut surface = SurfacePlot::new(x, y, z)
            .unwrap()
            .with_image_mode(true)
            .with_colormap(ColorMap::Gray)
            .with_color_limits(Some((0.0, 1.0)));
        let render_data = surface.render_data();

        assert_eq!(render_data.pipeline_type, PipelineType::Textured);
        assert_eq!(render_data.vertices.len(), 4);
        assert_eq!(
            render_data.indices.as_deref(),
            Some(&[0, 1, 2, 0, 2, 3][..])
        );

        let Some(ImageData::Rgba8 {
            width,
            height,
            data,
        }) = render_data.image
        else {
            panic!("image-mode surfaces should carry an RGBA texture payload");
        };
        assert_eq!((width, height), (2, 2));
        assert_eq!(data.len(), 16);

        // Image data is row-major, top-to-bottom. The top image row corresponds to
        // the highest Y data row.
        assert_eq!(&data[0..4], &[64, 64, 64, 255]);
        assert_eq!(&data[4..8], &[255, 255, 255, 255]);
        assert_eq!(&data[8..12], &[0, 0, 0, 255]);
        assert_eq!(&data[12..16], &[191, 191, 191, 255]);
    }

    #[test]
    fn image_mode_surface_preserves_non_square_rgba_pixel_order() {
        let x = vec![0.0, 1.0, 2.0];
        let y = vec![10.0, 20.0];
        let z = vec![vec![1.0, 2.0], vec![3.0, 4.0], vec![5.0, 6.0]];

        let mut surface = SurfacePlot::new(x, y, z)
            .unwrap()
            .with_image_mode(true)
            .with_colormap(ColorMap::Gray)
            .with_color_limits(Some((1.0, 6.0)));
        let render_data = surface.render_data();

        let Some(ImageData::Rgba8 {
            width,
            height,
            data,
        }) = render_data.image
        else {
            panic!("image-mode surfaces should carry an RGBA texture payload");
        };
        assert_eq!((width, height), (3, 2));
        assert_eq!(
            data,
            vec![
                51, 51, 51, 255, 153, 153, 153, 255, 255, 255, 255, 255, // top
                0, 0, 0, 255, 102, 102, 102, 255, 204, 204, 204, 255, // bottom
            ]
        );
    }

    #[test]
    fn test_colormap_mapping() {
        let jet = ColorMap::Jet;

        // Test boundary values
        let color_0 = jet.map_value(0.0);
        let color_1 = jet.map_value(1.0);

        assert!(color_0.x >= 0.0 && color_0.x <= 1.0);
        assert!(color_1.x >= 0.0 && color_1.x <= 1.0);

        // Test that different values give different colors
        let color_mid = jet.map_value(0.5);
        assert_ne!(color_0, color_mid);
        assert_ne!(color_mid, color_1);
    }

    #[test]
    fn test_surface_statistics() {
        let x = vec![0.0, 1.0, 2.0, 3.0];
        let y = vec![0.0, 1.0, 2.0];
        let z = vec![
            vec![0.0, 1.0, 2.0],
            vec![1.0, 2.0, 3.0],
            vec![2.0, 3.0, 4.0],
            vec![3.0, 4.0, 5.0],
        ];

        let surface = SurfacePlot::new(x, y, z).unwrap();
        let stats = surface.statistics();

        assert_eq!(stats.grid_points, 12); // 4 * 3
        assert_eq!(stats.triangle_count, 12); // (4-1) * (3-1) * 2
        assert_eq!(stats.x_resolution, 4);
        assert_eq!(stats.y_resolution, 3);
        assert!(stats.memory_usage > 0);
    }

    #[test]
    fn test_matlab_compat() {
        use super::matlab_compat::*;

        let x = vec![0.0, 1.0];
        let y = vec![0.0, 1.0];
        let z = vec![vec![0.0, 1.0], vec![1.0, 2.0]];

        let surface = surf(x.clone(), y.clone(), z.clone()).unwrap();
        assert!(!surface.wireframe);

        let mesh_plot = mesh(x.clone(), y.clone(), z.clone()).unwrap();
        assert!(mesh_plot.wireframe);

        let colormap_surface = surf_with_colormap(x, y, z, "viridis").unwrap();
        assert_eq!(colormap_surface.colormap, ColorMap::Viridis);
    }

    #[test]
    fn colormap_from_name_accepts_canonical_names_and_aliases() {
        let cases = [
            ("parula", ColorMap::Parula),
            ("colorcube", ColorMap::ColorCube),
            ("viridis", ColorMap::Viridis),
            ("plasma", ColorMap::Plasma),
            ("inferno", ColorMap::Inferno),
            ("magma", ColorMap::Magma),
            ("turbo", ColorMap::Turbo),
            ("jet", ColorMap::Jet),
            ("hot", ColorMap::Hot),
            ("cool", ColorMap::Cool),
            ("spring", ColorMap::Spring),
            ("summer", ColorMap::Summer),
            ("autumn", ColorMap::Autumn),
            ("winter", ColorMap::Winter),
            ("gray", ColorMap::Gray),
            ("grey", ColorMap::Gray),
            ("bone", ColorMap::Bone),
            ("copper", ColorMap::Copper),
            ("pink", ColorMap::Pink),
            ("lines", ColorMap::Lines),
        ];

        for (name, expected) in cases {
            assert_eq!(ColorMap::from_name(name), Some(expected), "{name}");
        }
        for name in ColorMap::CANONICAL_NAMES
            .iter()
            .chain(ColorMap::ALIASES.iter())
            .copied()
        {
            assert!(
                ColorMap::from_name(name).is_some(),
                "colormap table entry should parse: {name}"
            );
        }
    }

    #[test]
    fn colormap_from_name_normalizes_and_rejects_unknown_names() {
        assert_eq!(ColorMap::from_name(" Turbo "), Some(ColorMap::Turbo));
        assert_eq!(ColorMap::from_name("GREY"), Some(ColorMap::Gray));
        assert_eq!(ColorMap::from_name("hsv"), None);
        assert!(!ColorMap::CANONICAL_NAMES.contains(&"hsv"));
        assert!(!ColorMap::ALIASES.contains(&"hsv"));
        assert_eq!(ColorMap::from_name("not-a-colormap"), None);
    }

    #[test]
    fn listed_colormap_serialization_preserves_rows() {
        let map = ColorMap::from_rgb_rows(vec![[0.0, 0.5, 1.0], [1.0, 0.25, 0.0]])
            .expect("listed colormap");
        let token = map.to_serialized_token();
        assert!(token.starts_with("listed:"));

        let parsed = ColorMap::from_serialized_token(&token).expect("parse listed colormap");
        let ColorMap::Listed(colors) = parsed else {
            panic!("expected listed colormap");
        };
        assert_eq!(colors.as_ref(), &[[0.0, 0.5, 1.0], [1.0, 0.25, 0.0]]);
    }
}