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

runmat_plot/plots/
surface.rs
