runmat_plot/plots/
surface.rs

1//! 3D surface plot implementation
2//!
3//! High-performance GPU-accelerated 3D surface rendering.
4
5use crate::core::{
6    BoundingBox, DrawCall, GpuVertexBuffer, Material, PipelineType, RenderData, Vertex,
7};
8use glam::{Vec3, Vec4};
9
10/// High-performance GPU-accelerated 3D surface plot
11#[derive(Debug, Clone)]
12pub struct SurfacePlot {
13    /// Grid data (Z values at X,Y coordinates)
14    pub x_data: Vec<f64>,
15    pub y_data: Vec<f64>,
16    pub z_data: Option<Vec<Vec<f64>>>, // Host data when available
17    /// Grid resolution for rendering/index generation (kept even for GPU-backed plots).
18    x_len: usize,
19    y_len: usize,
20
21    /// Surface properties
22    pub colormap: ColorMap,
23    pub shading_mode: ShadingMode,
24    pub wireframe: bool,
25    pub alpha: f32,
26    /// If true, render Z at 0 (flat), but color-map using Z values
27    pub flatten_z: bool,
28
29    /// If true, this flattened surface should behave like a 2D image for camera/UI decisions.
30    pub image_mode: bool,
31
32    /// Optional color limits override for mapping Z -> color (caxis)
33    pub color_limits: Option<(f64, f64)>,
34
35    /// Optional per-vertex color grid (for RGB images); if set, overrides colormap mapping
36    pub color_grid: Option<Vec<Vec<Vec4>>>, // [x_index][y_index] -> RGBA
37
38    /// Lighting and material
39    pub lighting_enabled: bool,
40    pub ambient_strength: f32,
41    pub diffuse_strength: f32,
42    pub specular_strength: f32,
43    pub shininess: f32,
44
45    /// Metadata
46    pub label: Option<String>,
47    pub visible: bool,
48
49    /// Generated rendering data (cached)
50    vertices: Option<Vec<Vertex>>,
51    indices: Option<Vec<u32>>,
52    bounds: Option<BoundingBox>,
53    dirty: bool,
54    gpu_vertices: Option<GpuVertexBuffer>,
55    gpu_vertex_count: Option<usize>,
56    gpu_bounds: Option<BoundingBox>,
57}
58
59/// Color mapping schemes
60#[derive(Debug, Clone, Copy, PartialEq)]
61pub enum ColorMap {
62    /// MATLAB-compatible colormaps
63    Jet,
64    Hot,
65    Cool,
66    Spring,
67    Summer,
68    Autumn,
69    Winter,
70    Gray,
71    Bone,
72    Copper,
73    Pink,
74    Lines,
75
76    /// Scientific colormaps
77    Viridis,
78    Plasma,
79    Inferno,
80    Magma,
81    Turbo,
82
83    /// Perceptually uniform
84    Parula,
85
86    /// Custom color ranges
87    Custom(Vec4, Vec4), // (min_color, max_color)
88}
89
90impl ColorMap {
91    pub const CANONICAL_NAMES: &[&str] = &[
92        "parula", "viridis", "plasma", "inferno", "magma", "turbo", "jet", "hot", "cool", "spring",
93        "summer", "autumn", "winter", "gray", "bone", "copper", "pink", "lines",
94    ];
95
96    pub const ALIASES: &[&str] = &["grey"];
97
98    pub fn from_name(name: &str) -> Option<Self> {
99        match name.trim().to_ascii_lowercase().as_str() {
100            "parula" => Some(Self::Parula),
101            "viridis" => Some(Self::Viridis),
102            "plasma" => Some(Self::Plasma),
103            "inferno" => Some(Self::Inferno),
104            "magma" => Some(Self::Magma),
105            "turbo" => Some(Self::Turbo),
106            "jet" => Some(Self::Jet),
107            "hot" => Some(Self::Hot),
108            "cool" => Some(Self::Cool),
109            "spring" => Some(Self::Spring),
110            "summer" => Some(Self::Summer),
111            "autumn" => Some(Self::Autumn),
112            "winter" => Some(Self::Winter),
113            "gray" | "grey" => Some(Self::Gray),
114            "bone" => Some(Self::Bone),
115            "copper" => Some(Self::Copper),
116            "pink" => Some(Self::Pink),
117            "lines" => Some(Self::Lines),
118            _ => None,
119        }
120    }
121}
122
123/// Surface shading modes
124#[derive(Debug, Clone, Copy, PartialEq)]
125pub enum ShadingMode {
126    /// Flat shading (per-face normals)
127    Flat,
128    /// Smooth shading (interpolated normals)
129    Smooth,
130    /// Faceted (flat with visible edges)
131    Faceted,
132    /// No shading (just color mapping)
133    None,
134}
135
136impl Default for ColorMap {
137    fn default() -> Self {
138        Self::Viridis
139    }
140}
141
142impl Default for ShadingMode {
143    fn default() -> Self {
144        Self::Smooth
145    }
146}
147
148impl SurfacePlot {
149    /// Create a new surface plot from meshgrid data
150    pub fn new(x_data: Vec<f64>, y_data: Vec<f64>, z_data: Vec<Vec<f64>>) -> Result<Self, String> {
151        // Validate dimensions
152        if z_data.len() != x_data.len() {
153            return Err(format!(
154                "Z data rows ({}) must match X data length ({})",
155                z_data.len(),
156                x_data.len()
157            ));
158        }
159
160        for (i, row) in z_data.iter().enumerate() {
161            if row.len() != y_data.len() {
162                return Err(format!(
163                    "Z data row {} length ({}) must match Y data length ({})",
164                    i,
165                    row.len(),
166                    y_data.len()
167                ));
168            }
169        }
170
171        Ok(Self {
172            x_len: x_data.len(),
173            y_len: y_data.len(),
174            x_data,
175            y_data,
176            z_data: Some(z_data),
177            colormap: ColorMap::default(),
178            shading_mode: ShadingMode::default(),
179            wireframe: false,
180            alpha: 1.0,
181            flatten_z: false,
182            image_mode: false,
183            color_limits: None,
184            color_grid: None,
185            lighting_enabled: true,
186            ambient_strength: 0.2,
187            diffuse_strength: 0.8,
188            specular_strength: 0.5,
189            shininess: 32.0,
190            label: None,
191            visible: true,
192            vertices: None,
193            indices: None,
194            bounds: None,
195            dirty: true,
196            gpu_vertices: None,
197            gpu_vertex_count: None,
198            gpu_bounds: None,
199        })
200    }
201
202    /// Create a surface plot backed by a GPU vertex buffer.
203    pub fn from_gpu_buffer(
204        x_len: usize,
205        y_len: usize,
206        buffer: GpuVertexBuffer,
207        vertex_count: usize,
208        bounds: BoundingBox,
209    ) -> Self {
210        Self {
211            x_data: Vec::new(),
212            y_data: Vec::new(),
213            z_data: None,
214            x_len,
215            y_len,
216            colormap: ColorMap::default(),
217            shading_mode: ShadingMode::default(),
218            wireframe: false,
219            alpha: 1.0,
220            flatten_z: false,
221            image_mode: false,
222            color_limits: None,
223            color_grid: None,
224            lighting_enabled: true,
225            ambient_strength: 0.2,
226            diffuse_strength: 0.8,
227            specular_strength: 0.5,
228            shininess: 32.0,
229            label: None,
230            visible: true,
231            vertices: None,
232            indices: None,
233            bounds: Some(bounds),
234            dirty: false,
235            gpu_vertices: Some(buffer),
236            gpu_vertex_count: Some(vertex_count),
237            gpu_bounds: Some(bounds),
238        }
239    }
240
241    fn drop_gpu_if_possible(&mut self) {
242        if self.gpu_vertices.is_some() && self.z_data.is_some() {
243            self.invalidate_gpu_data();
244        }
245    }
246
247    /// Create surface from a function
248    pub fn from_function<F>(
249        x_range: (f64, f64),
250        y_range: (f64, f64),
251        resolution: (usize, usize),
252        func: F,
253    ) -> Result<Self, String>
254    where
255        F: Fn(f64, f64) -> f64,
256    {
257        let (x_res, y_res) = resolution;
258        if x_res < 2 || y_res < 2 {
259            return Err("Resolution must be at least 2x2".to_string());
260        }
261
262        let x_data: Vec<f64> = (0..x_res)
263            .map(|i| x_range.0 + (x_range.1 - x_range.0) * i as f64 / (x_res - 1) as f64)
264            .collect();
265
266        let y_data: Vec<f64> = (0..y_res)
267            .map(|j| y_range.0 + (y_range.1 - y_range.0) * j as f64 / (y_res - 1) as f64)
268            .collect();
269
270        let z_data: Vec<Vec<f64>> = x_data
271            .iter()
272            .map(|&x| y_data.iter().map(|&y| func(x, y)).collect())
273            .collect();
274
275        Self::new(x_data, y_data, z_data)
276    }
277
278    fn invalidate_gpu_data(&mut self) {
279        self.gpu_vertices = None;
280        self.gpu_vertex_count = None;
281        self.gpu_bounds = None;
282    }
283
284    /// Set color mapping
285    pub fn with_colormap(mut self, colormap: ColorMap) -> Self {
286        self.colormap = colormap;
287        self.dirty = true;
288        self.drop_gpu_if_possible();
289        self
290    }
291
292    /// Set shading mode
293    pub fn with_shading(mut self, shading: ShadingMode) -> Self {
294        self.shading_mode = shading;
295        self.dirty = true;
296        self.drop_gpu_if_possible();
297        self
298    }
299
300    /// Enable/disable wireframe
301    pub fn with_wireframe(mut self, enabled: bool) -> Self {
302        self.wireframe = enabled;
303        self.dirty = true;
304        self.drop_gpu_if_possible();
305        self
306    }
307
308    /// Set transparency
309    pub fn with_alpha(mut self, alpha: f32) -> Self {
310        self.alpha = alpha.clamp(0.0, 1.0);
311        self.dirty = true;
312        self.drop_gpu_if_possible();
313        self
314    }
315
316    /// Render surface flat in Z while mapping colors from Z values (for imagesc/imshow)
317    pub fn with_flatten_z(mut self, enabled: bool) -> Self {
318        self.flatten_z = enabled;
319        self.dirty = true;
320        self.drop_gpu_if_possible();
321        self
322    }
323
324    pub fn with_image_mode(mut self, enabled: bool) -> Self {
325        self.image_mode = enabled;
326        self.dirty = true;
327        self.drop_gpu_if_possible();
328        self
329    }
330
331    /// Override color mapping limits (caxis)
332    pub fn with_color_limits(mut self, limits: Option<(f64, f64)>) -> Self {
333        self.color_limits = limits;
334        self.dirty = true;
335        self.drop_gpu_if_possible();
336        self
337    }
338
339    /// Mutably set color mapping limits (caxis)
340    pub fn set_color_limits(&mut self, limits: Option<(f64, f64)>) {
341        self.color_limits = limits;
342        self.dirty = true;
343        self.drop_gpu_if_possible();
344    }
345
346    /// Provide explicit per-vertex colors (RGB[A])
347    pub fn with_color_grid(mut self, grid: Vec<Vec<Vec4>>) -> Self {
348        self.color_grid = Some(grid);
349        self.dirty = true;
350        self.drop_gpu_if_possible();
351        self
352    }
353
354    /// Set plot label for legends
355    pub fn with_label<S: Into<String>>(mut self, label: S) -> Self {
356        self.label = Some(label.into());
357        self
358    }
359
360    /// Get the number of grid points
361    pub fn len(&self) -> usize {
362        self.x_len * self.y_len
363    }
364
365    /// Check if the surface has no data
366    pub fn is_empty(&self) -> bool {
367        self.x_len == 0 || self.y_len == 0
368    }
369
370    /// Get the bounding box of the surface
371    pub fn bounds(&mut self) -> BoundingBox {
372        if self.dirty || self.bounds.is_none() {
373            self.compute_bounds();
374        }
375        self.bounds.unwrap()
376    }
377
378    /// Compute bounding box
379    fn compute_bounds(&mut self) {
380        if let Some(bounds) = self.gpu_bounds {
381            self.bounds = Some(bounds);
382            return;
383        }
384
385        let mut min_x = f32::INFINITY;
386        let mut max_x = f32::NEG_INFINITY;
387        let mut min_y = f32::INFINITY;
388        let mut max_y = f32::NEG_INFINITY;
389        let mut min_z = f32::INFINITY;
390        let mut max_z = f32::NEG_INFINITY;
391
392        for &x in &self.x_data {
393            min_x = min_x.min(x as f32);
394            max_x = max_x.max(x as f32);
395        }
396
397        for &y in &self.y_data {
398            min_y = min_y.min(y as f32);
399            max_y = max_y.max(y as f32);
400        }
401
402        if let Some(rows) = &self.z_data {
403            for row in rows {
404                for &z in row {
405                    if z.is_finite() {
406                        min_z = min_z.min(z as f32);
407                        max_z = max_z.max(z as f32);
408                    }
409                }
410            }
411        }
412
413        self.bounds = Some(BoundingBox::new(
414            Vec3::new(min_x, min_y, min_z),
415            Vec3::new(max_x, max_y, max_z),
416        ));
417    }
418
419    /// Get plot statistics for debugging
420    pub fn statistics(&self) -> SurfaceStatistics {
421        let grid_size = self.x_len * self.y_len;
422        let triangle_count = if self.x_len > 1 && self.y_len > 1 {
423            (self.x_len - 1) * (self.y_len - 1) * 2
424        } else {
425            0
426        };
427
428        SurfaceStatistics {
429            grid_points: grid_size,
430            triangle_count,
431            x_resolution: self.x_len,
432            y_resolution: self.y_len,
433            memory_usage: self.estimated_memory_usage(),
434        }
435    }
436
437    /// Estimate memory usage in bytes
438    pub fn estimated_memory_usage(&self) -> usize {
439        let data_size = std::mem::size_of::<f64>()
440            * (self.x_data.len()
441                + self.y_data.len()
442                + self
443                    .z_data
444                    .as_ref()
445                    .map_or(0, |z| z.len() * self.y_data.len()));
446
447        let vertices_size = self
448            .vertices
449            .as_ref()
450            .map_or(0, |v| v.len() * std::mem::size_of::<Vertex>());
451
452        let indices_size = self
453            .indices
454            .as_ref()
455            .map_or(0, |i| i.len() * std::mem::size_of::<u32>());
456
457        let gpu_size = self.gpu_vertex_count.unwrap_or(0) * std::mem::size_of::<Vertex>();
458
459        data_size + vertices_size + indices_size + gpu_size
460    }
461
462    /// Generate vertices for surface mesh
463    fn generate_vertices(&mut self) -> &Vec<Vertex> {
464        if self.gpu_vertices.is_some() {
465            if self.vertices.is_none() {
466                self.vertices = Some(Vec::new());
467            }
468            return self.vertices.as_ref().unwrap();
469        }
470
471        if self.dirty || self.vertices.is_none() {
472            log::trace!(
473                target: "runmat_plot",
474                "surface gen vertices {} x {}",
475                self.x_data.len(),
476                self.y_data.len()
477            );
478
479            let mut vertices = Vec::new();
480
481            // Determine color mapping range
482            let z_rows = self
483                .z_data
484                .as_ref()
485                .expect("CPU surface data missing during vertex generation");
486            let (min_z, max_z) = if let Some((lo, hi)) = self.color_limits {
487                (lo, hi)
488            } else {
489                let mut min_z = f64::INFINITY;
490                let mut max_z = f64::NEG_INFINITY;
491                for row in z_rows {
492                    for &z in row {
493                        if z.is_finite() {
494                            min_z = min_z.min(z);
495                            max_z = max_z.max(z);
496                        }
497                    }
498                }
499                (min_z, max_z)
500            };
501            let z_range = (max_z - min_z).max(f64::MIN_POSITIVE);
502
503            // Generate vertices for each grid point
504            for (i, &x) in self.x_data.iter().enumerate() {
505                for (j, &y) in self.y_data.iter().enumerate() {
506                    let z = z_rows[i][j];
507                    let z_pos = if self.flatten_z { 0.0 } else { z as f32 };
508                    let position = Vec3::new(x as f32, y as f32, z_pos);
509
510                    // Simple normal calculation (can be improved with proper gradients)
511                    let normal = Vec3::new(0.0, 0.0, 1.0); // Placeholder
512
513                    // Determine color: explicit grid (RGB) or colormap from Z
514                    let color = if let Some(grid) = &self.color_grid {
515                        let c = grid[i][j];
516                        Vec4::new(c.x, c.y, c.z, c.w)
517                    } else {
518                        let t = ((z - min_z) / z_range) as f32;
519                        let color_rgb = self.colormap.map_value(t.clamp(0.0, 1.0));
520                        Vec4::new(color_rgb.x, color_rgb.y, color_rgb.z, self.alpha)
521                    };
522
523                    vertices.push(Vertex {
524                        position: position.to_array(),
525                        normal: normal.to_array(),
526                        color: color.to_array(),
527                        tex_coords: [
528                            i as f32 / (self.x_data.len() - 1).max(1) as f32,
529                            j as f32 / (self.y_data.len() - 1).max(1) as f32,
530                        ],
531                    });
532                }
533            }
534
535            log::trace!(target: "runmat_plot", "surface vertices={}", vertices.len());
536            self.vertices = Some(vertices);
537        }
538        self.vertices.as_ref().unwrap()
539    }
540
541    /// Generate indices for surface triangulation
542    fn generate_indices(&mut self) -> &Vec<u32> {
543        if self.dirty || self.indices.is_none() {
544            log::trace!(target: "runmat_plot", "surface generating indices");
545
546            let mut indices = Vec::new();
547            let x_res = self.x_len;
548            let y_res = self.y_len;
549
550            // Generate triangle indices for surface mesh
551            for i in 0..x_res - 1 {
552                for j in 0..y_res - 1 {
553                    let base = (i * y_res + j) as u32;
554                    let next_row = base + y_res as u32;
555
556                    // Two triangles per quad
557                    // Triangle 1: (i,j), (i+1,j), (i,j+1)
558                    indices.push(base);
559                    indices.push(next_row);
560                    indices.push(base + 1);
561
562                    // Triangle 2: (i+1,j), (i+1,j+1), (i,j+1)
563                    indices.push(next_row);
564                    indices.push(next_row + 1);
565                    indices.push(base + 1);
566                }
567            }
568
569            log::trace!(target: "runmat_plot", "surface indices={}", indices.len());
570            self.indices = Some(indices);
571            self.dirty = false;
572        }
573        self.indices.as_ref().unwrap()
574    }
575
576    fn generate_wireframe_indices(&self) -> Vec<u32> {
577        let mut indices = Vec::new();
578        if self.x_len < 2 || self.y_len < 2 {
579            return indices;
580        }
581
582        // Horizontal grid edges (along Y for each X row)
583        for i in 0..self.x_len {
584            for j in 0..(self.y_len - 1) {
585                let a = (i * self.y_len + j) as u32;
586                let b = (i * self.y_len + j + 1) as u32;
587                indices.push(a);
588                indices.push(b);
589            }
590        }
591
592        // Vertical grid edges (along X for each Y column)
593        for i in 0..(self.x_len - 1) {
594            for j in 0..self.y_len {
595                let a = (i * self.y_len + j) as u32;
596                let b = ((i + 1) * self.y_len + j) as u32;
597                indices.push(a);
598                indices.push(b);
599            }
600        }
601
602        indices
603    }
604
605    /// Generate complete render data for the graphics pipeline
606    pub fn render_data(&mut self) -> RenderData {
607        log::debug!(
608            target: "runmat_plot",
609            "surface render_data start: {} x {}",
610            self.x_len,
611            self.y_len
612        );
613
614        let using_gpu = self.gpu_vertices.is_some();
615        let bounds = self.bounds();
616        let vertices = if using_gpu {
617            Vec::new()
618        } else {
619            self.generate_vertices().clone()
620        };
621        let indices = if self.wireframe {
622            self.generate_wireframe_indices()
623        } else {
624            self.generate_indices().clone()
625        };
626
627        let material = Material {
628            albedo: Vec4::new(1.0, 1.0, 1.0, self.alpha),
629            ..Default::default()
630        };
631
632        let vertex_count = if using_gpu {
633            self.gpu_vertex_count.unwrap_or(0)
634        } else {
635            vertices.len()
636        };
637
638        log::debug!(
639            target: "runmat_plot",
640            "surface render_data generated: vertex_count={} (gpu={}), indices={}",
641            vertex_count,
642            using_gpu,
643            indices.len()
644        );
645
646        let draw_call = DrawCall {
647            vertex_offset: 0,
648            vertex_count,
649            index_offset: Some(0),
650            index_count: Some(indices.len()),
651            instance_count: 1,
652        };
653
654        log::trace!(target: "runmat_plot", "surface render_data done");
655
656        RenderData {
657            pipeline_type: if self.wireframe {
658                PipelineType::Lines
659            } else {
660                PipelineType::Triangles
661            },
662            vertices,
663            indices: Some(indices),
664
665            gpu_vertices: self.gpu_vertices.clone(),
666            bounds: Some(bounds),
667            material,
668            draw_calls: vec![draw_call],
669            image: None,
670        }
671    }
672}
673
674/// Surface plot performance and data statistics
675#[derive(Debug, Clone)]
676pub struct SurfaceStatistics {
677    pub grid_points: usize,
678    pub triangle_count: usize,
679    pub x_resolution: usize,
680    pub y_resolution: usize,
681    pub memory_usage: usize,
682}
683
684impl ColorMap {
685    /// Map a normalized value [0,1] to a color
686    pub fn map_value(&self, t: f32) -> Vec3 {
687        let t = t.clamp(0.0, 1.0);
688
689        match self {
690            ColorMap::Jet => self.jet_colormap(t),
691            ColorMap::Hot => self.hot_colormap(t),
692            ColorMap::Cool => self.cool_colormap(t),
693            ColorMap::Spring => self.spring_colormap(t),
694            ColorMap::Summer => self.summer_colormap(t),
695            ColorMap::Autumn => self.autumn_colormap(t),
696            ColorMap::Winter => self.winter_colormap(t),
697            ColorMap::Gray => Vec3::splat(t),
698            ColorMap::Bone => self.bone_colormap(t),
699            ColorMap::Copper => self.copper_colormap(t),
700            ColorMap::Pink => self.pink_colormap(t),
701            ColorMap::Lines => self.lines_colormap(t),
702            ColorMap::Viridis => self.viridis_colormap(t),
703            ColorMap::Plasma => self.plasma_colormap(t),
704            ColorMap::Inferno => self.inferno_colormap(t),
705            ColorMap::Magma => self.magma_colormap(t),
706            ColorMap::Turbo => self.turbo_colormap(t),
707            ColorMap::Parula => self.parula_colormap(t),
708            ColorMap::Custom(min_color, max_color) => {
709                min_color.truncate().lerp(max_color.truncate(), t)
710            }
711        }
712    }
713
714    /// MATLAB Jet colormap
715    fn jet_colormap(&self, t: f32) -> Vec3 {
716        let r = (1.5 - 4.0 * (t - 0.75).abs()).clamp(0.0, 1.0);
717        let g = (1.5 - 4.0 * (t - 0.5).abs()).clamp(0.0, 1.0);
718        let b = (1.5 - 4.0 * (t - 0.25).abs()).clamp(0.0, 1.0);
719        Vec3::new(r, g, b)
720    }
721
722    /// Hot colormap (black -> red -> yellow -> white)
723    fn hot_colormap(&self, t: f32) -> Vec3 {
724        if t < 1.0 / 3.0 {
725            Vec3::new(3.0 * t, 0.0, 0.0)
726        } else if t < 2.0 / 3.0 {
727            Vec3::new(1.0, 3.0 * t - 1.0, 0.0)
728        } else {
729            Vec3::new(1.0, 1.0, 3.0 * t - 2.0)
730        }
731    }
732
733    /// Cool colormap (cyan -> magenta)
734    fn cool_colormap(&self, t: f32) -> Vec3 {
735        Vec3::new(t, 1.0 - t, 1.0)
736    }
737
738    /// Viridis colormap (perceptually uniform)
739    fn viridis_colormap(&self, t: f32) -> Vec3 {
740        // Simplified Viridis approximation
741        let r = (0.267004 + t * (0.993248 - 0.267004)).clamp(0.0, 1.0);
742        let g = (0.004874 + t * (0.906157 - 0.004874)).clamp(0.0, 1.0);
743        let b = (0.329415 + t * (0.143936 - 0.329415) + t * t * 0.5).clamp(0.0, 1.0);
744        Vec3::new(r, g, b)
745    }
746
747    /// Plasma colormap (perceptually uniform)
748    fn plasma_colormap(&self, t: f32) -> Vec3 {
749        // Simplified Plasma approximation
750        let r = (0.050383 + t * (0.940015 - 0.050383)).clamp(0.0, 1.0);
751        let g = (0.029803 + t * (0.975158 - 0.029803) * (1.0 - t)).clamp(0.0, 1.0);
752        let b = (0.527975 + t * (0.131326 - 0.527975)).clamp(0.0, 1.0);
753        Vec3::new(r, g, b)
754    }
755
756    /// Spring colormap (magenta -> yellow)
757    fn spring_colormap(&self, t: f32) -> Vec3 {
758        Vec3::new(1.0, t, 1.0 - t)
759    }
760
761    /// Summer colormap (green -> yellow)
762    fn summer_colormap(&self, t: f32) -> Vec3 {
763        Vec3::new(t, 0.5 + 0.5 * t, 0.4)
764    }
765
766    /// Autumn colormap (red -> yellow)
767    fn autumn_colormap(&self, t: f32) -> Vec3 {
768        Vec3::new(1.0, t, 0.0)
769    }
770
771    /// Winter colormap (blue -> green)
772    fn winter_colormap(&self, t: f32) -> Vec3 {
773        Vec3::new(0.0, t, 1.0 - 0.5 * t)
774    }
775
776    /// Bone colormap (black -> white with blue tint)
777    fn bone_colormap(&self, t: f32) -> Vec3 {
778        if t < 3.0 / 8.0 {
779            Vec3::new(7.0 / 8.0 * t, 7.0 / 8.0 * t, 29.0 / 24.0 * t)
780        } else {
781            Vec3::new(
782                (29.0 + 7.0 * t) / 24.0,
783                (29.0 + 7.0 * t) / 24.0,
784                (29.0 + 7.0 * t) / 24.0,
785            )
786        }
787    }
788
789    /// Copper colormap (black -> copper)
790    fn copper_colormap(&self, t: f32) -> Vec3 {
791        Vec3::new((1.25 * t).min(1.0), 0.7812 * t, 0.4975 * t)
792    }
793
794    /// Pink colormap (black -> pink -> white)
795    fn pink_colormap(&self, t: f32) -> Vec3 {
796        let sqrt_t = t.sqrt();
797        if t < 3.0 / 8.0 {
798            Vec3::new(14.0 / 9.0 * sqrt_t, 2.0 / 3.0 * sqrt_t, 2.0 / 3.0 * sqrt_t)
799        } else {
800            Vec3::new(
801                2.0 * sqrt_t - 1.0 / 3.0,
802                8.0 / 9.0 * sqrt_t + 1.0 / 3.0,
803                8.0 / 9.0 * sqrt_t + 1.0 / 3.0,
804            )
805        }
806    }
807
808    /// Lines colormap (cycling through basic colors)
809    fn lines_colormap(&self, t: f32) -> Vec3 {
810        let _phase = (t * 7.0) % 1.0; // For future use in color transitions
811        let index = (t * 7.0) as usize % 7;
812        match index {
813            0 => Vec3::new(0.0, 0.0, 1.0),    // Blue
814            1 => Vec3::new(0.0, 0.5, 0.0),    // Green
815            2 => Vec3::new(1.0, 0.0, 0.0),    // Red
816            3 => Vec3::new(0.0, 0.75, 0.75),  // Cyan
817            4 => Vec3::new(0.75, 0.0, 0.75),  // Magenta
818            5 => Vec3::new(0.75, 0.75, 0.0),  // Yellow
819            _ => Vec3::new(0.25, 0.25, 0.25), // Dark gray
820        }
821    }
822
823    /// Inferno colormap (perceptually uniform)
824    fn inferno_colormap(&self, t: f32) -> Vec3 {
825        // Simplified Inferno approximation
826        let r = (0.001462 + t * (0.988362 - 0.001462)).clamp(0.0, 1.0);
827        let g = (0.000466 + t * t * (0.982895 - 0.000466)).clamp(0.0, 1.0);
828        let b = (0.013866 + t * (1.0 - t) * (0.416065 - 0.013866)).clamp(0.0, 1.0);
829        Vec3::new(r, g, b)
830    }
831
832    /// Magma colormap (perceptually uniform)
833    fn magma_colormap(&self, t: f32) -> Vec3 {
834        // Simplified Magma approximation
835        let r = (0.001462 + t * (0.987053 - 0.001462)).clamp(0.0, 1.0);
836        let g = (0.000466 + t * t * (0.991438 - 0.000466)).clamp(0.0, 1.0);
837        let b = (0.013866 + t * (0.644237 - 0.013866) * (1.0 - t)).clamp(0.0, 1.0);
838        Vec3::new(r, g, b)
839    }
840
841    /// Turbo colormap (improved rainbow)
842    fn turbo_colormap(&self, t: f32) -> Vec3 {
843        // Simplified Turbo approximation (Google's improved rainbow)
844        let r = if t < 0.5 {
845            (0.13 + 0.87 * (2.0 * t).powf(0.25)).clamp(0.0, 1.0)
846        } else {
847            (0.8685 + 0.1315 * (2.0 * (1.0 - t)).powf(0.25)).clamp(0.0, 1.0)
848        };
849
850        let g = if t < 0.25 {
851            4.0 * t
852        } else if t < 0.75 {
853            1.0
854        } else {
855            1.0 - 4.0 * (t - 0.75)
856        }
857        .clamp(0.0, 1.0);
858
859        let b = if t < 0.5 {
860            (0.8 * (1.0 - 2.0 * t).powf(0.25)).clamp(0.0, 1.0)
861        } else {
862            (0.1 + 0.9 * (2.0 * t - 1.0).powf(0.25)).clamp(0.0, 1.0)
863        };
864
865        Vec3::new(r, g, b)
866    }
867
868    /// Parula colormap (MATLAB's default)
869    fn parula_colormap(&self, t: f32) -> Vec3 {
870        // Simplified Parula approximation
871        let r = if t < 0.25 {
872            0.2081 * (1.0 - t)
873        } else if t < 0.5 {
874            t - 0.25
875        } else if t < 0.75 {
876            1.0
877        } else {
878            1.0 - 0.5 * (t - 0.75)
879        }
880        .clamp(0.0, 1.0);
881
882        let g = if t < 0.125 {
883            0.1663 * t / 0.125
884        } else if t < 0.375 {
885            0.1663 + (0.7079 - 0.1663) * (t - 0.125) / 0.25
886        } else if t < 0.625 {
887            0.7079 + (0.9839 - 0.7079) * (t - 0.375) / 0.25
888        } else {
889            0.9839 * (1.0 - (t - 0.625) / 0.375)
890        }
891        .clamp(0.0, 1.0);
892
893        let b = if t < 0.25 {
894            0.5 + 0.5 * t / 0.25
895        } else if t < 0.5 {
896            1.0
897        } else {
898            1.0 - 2.0 * (t - 0.5)
899        }
900        .clamp(0.0, 1.0);
901
902        Vec3::new(r, g, b)
903    }
904
905    /// Default colormap fallback
906    #[allow(dead_code)] // Fallback method for colormap errors
907    fn default_colormap(&self, t: f32) -> Vec3 {
908        // Use a simple RGB transition as fallback
909        if t < 0.5 {
910            Vec3::new(0.0, 2.0 * t, 1.0 - 2.0 * t)
911        } else {
912            Vec3::new(2.0 * (t - 0.5), 1.0 - 2.0 * (t - 0.5), 0.0)
913        }
914    }
915}
916
917/// MATLAB-compatible surface plot creation utilities
918pub mod matlab_compat {
919    use super::*;
920
921    /// Create a surface plot (equivalent to MATLAB's `surf(X, Y, Z)`)
922    pub fn surf(x: Vec<f64>, y: Vec<f64>, z: Vec<Vec<f64>>) -> Result<SurfacePlot, String> {
923        SurfacePlot::new(x, y, z)
924    }
925
926    /// Create a mesh plot (wireframe surface)
927    pub fn mesh(x: Vec<f64>, y: Vec<f64>, z: Vec<Vec<f64>>) -> Result<SurfacePlot, String> {
928        Ok(SurfacePlot::new(x, y, z)?
929            .with_wireframe(true)
930            .with_shading(ShadingMode::None))
931    }
932
933    /// Create surface from meshgrid
934    pub fn meshgrid_surf(
935        x_range: (f64, f64),
936        y_range: (f64, f64),
937        resolution: (usize, usize),
938        func: impl Fn(f64, f64) -> f64,
939    ) -> Result<SurfacePlot, String> {
940        SurfacePlot::from_function(x_range, y_range, resolution, func)
941    }
942
943    /// Create surface with specific colormap
944    pub fn surf_with_colormap(
945        x: Vec<f64>,
946        y: Vec<f64>,
947        z: Vec<Vec<f64>>,
948        colormap: &str,
949    ) -> Result<SurfacePlot, String> {
950        let cmap =
951            ColorMap::from_name(colormap).ok_or_else(|| format!("Unknown colormap: {colormap}"))?;
952
953        Ok(SurfacePlot::new(x, y, z)?.with_colormap(cmap))
954    }
955}
956
957#[cfg(test)]
958mod tests {
959    use super::*;
960
961    #[test]
962    fn test_surface_plot_creation() {
963        let x = vec![0.0, 1.0, 2.0];
964        let y = vec![0.0, 1.0];
965        let z = vec![vec![0.0, 1.0], vec![1.0, 2.0], vec![2.0, 3.0]];
966
967        let surface = SurfacePlot::new(x, y, z).unwrap();
968
969        assert_eq!(surface.x_data.len(), 3);
970        assert_eq!(surface.y_data.len(), 2);
971        let rows = surface.z_data.as_ref().unwrap();
972        assert_eq!(rows.len(), 3);
973        assert_eq!(rows[0].len(), 2);
974        assert!(surface.visible);
975    }
976
977    #[test]
978    fn test_surface_from_function() {
979        let surface =
980            SurfacePlot::from_function((-2.0, 2.0), (-2.0, 2.0), (10, 10), |x, y| x * x + y * y)
981                .unwrap();
982
983        assert_eq!(surface.x_data.len(), 10);
984        assert_eq!(surface.y_data.len(), 10);
985        let rows = surface.z_data.as_ref().unwrap();
986        assert_eq!(rows.len(), 10);
987
988        // Check that function is evaluated correctly
989        assert_eq!(rows[0][0], 8.0); // (-2)^2 + (-2)^2 = 8
990    }
991
992    #[test]
993    fn test_surface_validation() {
994        let x = vec![0.0, 1.0];
995        let y = vec![0.0, 1.0, 2.0];
996        let z = vec![
997            vec![0.0, 1.0], // Wrong: should have 3 elements to match y
998            vec![1.0, 2.0],
999        ];
1000
1001        assert!(SurfacePlot::new(x, y, z).is_err());
1002    }
1003
1004    #[test]
1005    fn test_surface_styling() {
1006        let x = vec![0.0, 1.0];
1007        let y = vec![0.0, 1.0];
1008        let z = vec![vec![0.0, 1.0], vec![1.0, 2.0]];
1009
1010        let surface = SurfacePlot::new(x, y, z)
1011            .unwrap()
1012            .with_colormap(ColorMap::Hot)
1013            .with_wireframe(true)
1014            .with_alpha(0.8)
1015            .with_label("Test Surface");
1016
1017        assert_eq!(surface.colormap, ColorMap::Hot);
1018        assert!(surface.wireframe);
1019        assert_eq!(surface.alpha, 0.8);
1020        assert_eq!(surface.label, Some("Test Surface".to_string()));
1021    }
1022
1023    #[test]
1024    fn test_colormap_mapping() {
1025        let jet = ColorMap::Jet;
1026
1027        // Test boundary values
1028        let color_0 = jet.map_value(0.0);
1029        let color_1 = jet.map_value(1.0);
1030
1031        assert!(color_0.x >= 0.0 && color_0.x <= 1.0);
1032        assert!(color_1.x >= 0.0 && color_1.x <= 1.0);
1033
1034        // Test that different values give different colors
1035        let color_mid = jet.map_value(0.5);
1036        assert_ne!(color_0, color_mid);
1037        assert_ne!(color_mid, color_1);
1038    }
1039
1040    #[test]
1041    fn test_surface_statistics() {
1042        let x = vec![0.0, 1.0, 2.0, 3.0];
1043        let y = vec![0.0, 1.0, 2.0];
1044        let z = vec![
1045            vec![0.0, 1.0, 2.0],
1046            vec![1.0, 2.0, 3.0],
1047            vec![2.0, 3.0, 4.0],
1048            vec![3.0, 4.0, 5.0],
1049        ];
1050
1051        let surface = SurfacePlot::new(x, y, z).unwrap();
1052        let stats = surface.statistics();
1053
1054        assert_eq!(stats.grid_points, 12); // 4 * 3
1055        assert_eq!(stats.triangle_count, 12); // (4-1) * (3-1) * 2
1056        assert_eq!(stats.x_resolution, 4);
1057        assert_eq!(stats.y_resolution, 3);
1058        assert!(stats.memory_usage > 0);
1059    }
1060
1061    #[test]
1062    fn test_matlab_compat() {
1063        use super::matlab_compat::*;
1064
1065        let x = vec![0.0, 1.0];
1066        let y = vec![0.0, 1.0];
1067        let z = vec![vec![0.0, 1.0], vec![1.0, 2.0]];
1068
1069        let surface = surf(x.clone(), y.clone(), z.clone()).unwrap();
1070        assert!(!surface.wireframe);
1071
1072        let mesh_plot = mesh(x.clone(), y.clone(), z.clone()).unwrap();
1073        assert!(mesh_plot.wireframe);
1074
1075        let colormap_surface = surf_with_colormap(x, y, z, "viridis").unwrap();
1076        assert_eq!(colormap_surface.colormap, ColorMap::Viridis);
1077    }
1078
1079    #[test]
1080    fn colormap_from_name_accepts_canonical_names_and_aliases() {
1081        let cases = [
1082            ("parula", ColorMap::Parula),
1083            ("viridis", ColorMap::Viridis),
1084            ("plasma", ColorMap::Plasma),
1085            ("inferno", ColorMap::Inferno),
1086            ("magma", ColorMap::Magma),
1087            ("turbo", ColorMap::Turbo),
1088            ("jet", ColorMap::Jet),
1089            ("hot", ColorMap::Hot),
1090            ("cool", ColorMap::Cool),
1091            ("spring", ColorMap::Spring),
1092            ("summer", ColorMap::Summer),
1093            ("autumn", ColorMap::Autumn),
1094            ("winter", ColorMap::Winter),
1095            ("gray", ColorMap::Gray),
1096            ("grey", ColorMap::Gray),
1097            ("bone", ColorMap::Bone),
1098            ("copper", ColorMap::Copper),
1099            ("pink", ColorMap::Pink),
1100            ("lines", ColorMap::Lines),
1101        ];
1102
1103        for (name, expected) in cases {
1104            assert_eq!(ColorMap::from_name(name), Some(expected), "{name}");
1105        }
1106        for name in ColorMap::CANONICAL_NAMES
1107            .iter()
1108            .chain(ColorMap::ALIASES.iter())
1109            .copied()
1110        {
1111            assert!(
1112                ColorMap::from_name(name).is_some(),
1113                "colormap table entry should parse: {name}"
1114            );
1115        }
1116    }
1117
1118    #[test]
1119    fn colormap_from_name_normalizes_and_rejects_unknown_names() {
1120        assert_eq!(ColorMap::from_name(" Turbo "), Some(ColorMap::Turbo));
1121        assert_eq!(ColorMap::from_name("GREY"), Some(ColorMap::Gray));
1122        assert_eq!(ColorMap::from_name("hsv"), None);
1123        assert!(!ColorMap::CANONICAL_NAMES.contains(&"hsv"));
1124        assert!(!ColorMap::ALIASES.contains(&"hsv"));
1125        assert_eq!(ColorMap::from_name("not-a-colormap"), None);
1126    }
1127}
runmat_plot/plots/surface.rs

runmat_plot/plots/
surface.rs
