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    /// Generate complete render data for the graphics pipeline
544    pub fn render_data(&mut self) -> RenderData {
545        log::debug!(
546            target: "runmat_plot",
547            "surface render_data start: {} x {}",
548            self.x_len,
549            self.y_len
550        );
551
552        let using_gpu = self.gpu_vertices.is_some();
553        let bounds = self.bounds();
554        let vertices = if using_gpu {
555            Vec::new()
556        } else {
557            self.generate_vertices().clone()
558        };
559        let indices = self.generate_indices().clone();
560
561        let material = Material {
562            albedo: Vec4::new(1.0, 1.0, 1.0, self.alpha),
563            ..Default::default()
564        };
565
566        let vertex_count = if using_gpu {
567            self.gpu_vertex_count.unwrap_or(0)
568        } else {
569            vertices.len()
570        };
571
572        log::debug!(
573            target: "runmat_plot",
574            "surface render_data generated: vertex_count={} (gpu={}), indices={}",
575            vertex_count,
576            using_gpu,
577            indices.len()
578        );
579
580        let draw_call = DrawCall {
581            vertex_offset: 0,
582            vertex_count,
583            index_offset: Some(0),
584            index_count: Some(indices.len()),
585            instance_count: 1,
586        };
587
588        log::trace!(target: "runmat_plot", "surface render_data done");
589
590        RenderData {
591            pipeline_type: if self.wireframe {
592                PipelineType::Lines
593            } else {
594                PipelineType::Triangles
595            },
596            vertices,
597            indices: Some(indices),
598
599            gpu_vertices: self.gpu_vertices.clone(),
600            bounds: Some(bounds),
601            material,
602            draw_calls: vec![draw_call],
603            image: None,
604        }
605    }
606}
607
608/// Surface plot performance and data statistics
609#[derive(Debug, Clone)]
610pub struct SurfaceStatistics {
611    pub grid_points: usize,
612    pub triangle_count: usize,
613    pub x_resolution: usize,
614    pub y_resolution: usize,
615    pub memory_usage: usize,
616}
617
618impl ColorMap {
619    /// Map a normalized value [0,1] to a color
620    pub fn map_value(&self, t: f32) -> Vec3 {
621        let t = t.clamp(0.0, 1.0);
622
623        match self {
624            ColorMap::Jet => self.jet_colormap(t),
625            ColorMap::Hot => self.hot_colormap(t),
626            ColorMap::Cool => self.cool_colormap(t),
627            ColorMap::Spring => self.spring_colormap(t),
628            ColorMap::Summer => self.summer_colormap(t),
629            ColorMap::Autumn => self.autumn_colormap(t),
630            ColorMap::Winter => self.winter_colormap(t),
631            ColorMap::Gray => Vec3::splat(t),
632            ColorMap::Bone => self.bone_colormap(t),
633            ColorMap::Copper => self.copper_colormap(t),
634            ColorMap::Pink => self.pink_colormap(t),
635            ColorMap::Lines => self.lines_colormap(t),
636            ColorMap::Viridis => self.viridis_colormap(t),
637            ColorMap::Plasma => self.plasma_colormap(t),
638            ColorMap::Inferno => self.inferno_colormap(t),
639            ColorMap::Magma => self.magma_colormap(t),
640            ColorMap::Turbo => self.turbo_colormap(t),
641            ColorMap::Parula => self.parula_colormap(t),
642            ColorMap::Custom(min_color, max_color) => {
643                min_color.truncate().lerp(max_color.truncate(), t)
644            }
645        }
646    }
647
648    /// MATLAB Jet colormap
649    fn jet_colormap(&self, t: f32) -> Vec3 {
650        let r = (1.5 - 4.0 * (t - 0.75).abs()).clamp(0.0, 1.0);
651        let g = (1.5 - 4.0 * (t - 0.5).abs()).clamp(0.0, 1.0);
652        let b = (1.5 - 4.0 * (t - 0.25).abs()).clamp(0.0, 1.0);
653        Vec3::new(r, g, b)
654    }
655
656    /// Hot colormap (black -> red -> yellow -> white)
657    fn hot_colormap(&self, t: f32) -> Vec3 {
658        if t < 1.0 / 3.0 {
659            Vec3::new(3.0 * t, 0.0, 0.0)
660        } else if t < 2.0 / 3.0 {
661            Vec3::new(1.0, 3.0 * t - 1.0, 0.0)
662        } else {
663            Vec3::new(1.0, 1.0, 3.0 * t - 2.0)
664        }
665    }
666
667    /// Cool colormap (cyan -> magenta)
668    fn cool_colormap(&self, t: f32) -> Vec3 {
669        Vec3::new(t, 1.0 - t, 1.0)
670    }
671
672    /// Viridis colormap (perceptually uniform)
673    fn viridis_colormap(&self, t: f32) -> Vec3 {
674        // Simplified Viridis approximation
675        let r = (0.267004 + t * (0.993248 - 0.267004)).clamp(0.0, 1.0);
676        let g = (0.004874 + t * (0.906157 - 0.004874)).clamp(0.0, 1.0);
677        let b = (0.329415 + t * (0.143936 - 0.329415) + t * t * 0.5).clamp(0.0, 1.0);
678        Vec3::new(r, g, b)
679    }
680
681    /// Plasma colormap (perceptually uniform)
682    fn plasma_colormap(&self, t: f32) -> Vec3 {
683        // Simplified Plasma approximation
684        let r = (0.050383 + t * (0.940015 - 0.050383)).clamp(0.0, 1.0);
685        let g = (0.029803 + t * (0.975158 - 0.029803) * (1.0 - t)).clamp(0.0, 1.0);
686        let b = (0.527975 + t * (0.131326 - 0.527975)).clamp(0.0, 1.0);
687        Vec3::new(r, g, b)
688    }
689
690    /// Spring colormap (magenta -> yellow)
691    fn spring_colormap(&self, t: f32) -> Vec3 {
692        Vec3::new(1.0, t, 1.0 - t)
693    }
694
695    /// Summer colormap (green -> yellow)
696    fn summer_colormap(&self, t: f32) -> Vec3 {
697        Vec3::new(t, 0.5 + 0.5 * t, 0.4)
698    }
699
700    /// Autumn colormap (red -> yellow)
701    fn autumn_colormap(&self, t: f32) -> Vec3 {
702        Vec3::new(1.0, t, 0.0)
703    }
704
705    /// Winter colormap (blue -> green)
706    fn winter_colormap(&self, t: f32) -> Vec3 {
707        Vec3::new(0.0, t, 1.0 - 0.5 * t)
708    }
709
710    /// Bone colormap (black -> white with blue tint)
711    fn bone_colormap(&self, t: f32) -> Vec3 {
712        if t < 3.0 / 8.0 {
713            Vec3::new(7.0 / 8.0 * t, 7.0 / 8.0 * t, 29.0 / 24.0 * t)
714        } else {
715            Vec3::new(
716                (29.0 + 7.0 * t) / 24.0,
717                (29.0 + 7.0 * t) / 24.0,
718                (29.0 + 7.0 * t) / 24.0,
719            )
720        }
721    }
722
723    /// Copper colormap (black -> copper)
724    fn copper_colormap(&self, t: f32) -> Vec3 {
725        Vec3::new((1.25 * t).min(1.0), 0.7812 * t, 0.4975 * t)
726    }
727
728    /// Pink colormap (black -> pink -> white)
729    fn pink_colormap(&self, t: f32) -> Vec3 {
730        let sqrt_t = t.sqrt();
731        if t < 3.0 / 8.0 {
732            Vec3::new(14.0 / 9.0 * sqrt_t, 2.0 / 3.0 * sqrt_t, 2.0 / 3.0 * sqrt_t)
733        } else {
734            Vec3::new(
735                2.0 * sqrt_t - 1.0 / 3.0,
736                8.0 / 9.0 * sqrt_t + 1.0 / 3.0,
737                8.0 / 9.0 * sqrt_t + 1.0 / 3.0,
738            )
739        }
740    }
741
742    /// Lines colormap (cycling through basic colors)
743    fn lines_colormap(&self, t: f32) -> Vec3 {
744        let _phase = (t * 7.0) % 1.0; // For future use in color transitions
745        let index = (t * 7.0) as usize % 7;
746        match index {
747            0 => Vec3::new(0.0, 0.0, 1.0),    // Blue
748            1 => Vec3::new(0.0, 0.5, 0.0),    // Green
749            2 => Vec3::new(1.0, 0.0, 0.0),    // Red
750            3 => Vec3::new(0.0, 0.75, 0.75),  // Cyan
751            4 => Vec3::new(0.75, 0.0, 0.75),  // Magenta
752            5 => Vec3::new(0.75, 0.75, 0.0),  // Yellow
753            _ => Vec3::new(0.25, 0.25, 0.25), // Dark gray
754        }
755    }
756
757    /// Inferno colormap (perceptually uniform)
758    fn inferno_colormap(&self, t: f32) -> Vec3 {
759        // Simplified Inferno approximation
760        let r = (0.001462 + t * (0.988362 - 0.001462)).clamp(0.0, 1.0);
761        let g = (0.000466 + t * t * (0.982895 - 0.000466)).clamp(0.0, 1.0);
762        let b = (0.013866 + t * (1.0 - t) * (0.416065 - 0.013866)).clamp(0.0, 1.0);
763        Vec3::new(r, g, b)
764    }
765
766    /// Magma colormap (perceptually uniform)
767    fn magma_colormap(&self, t: f32) -> Vec3 {
768        // Simplified Magma approximation
769        let r = (0.001462 + t * (0.987053 - 0.001462)).clamp(0.0, 1.0);
770        let g = (0.000466 + t * t * (0.991438 - 0.000466)).clamp(0.0, 1.0);
771        let b = (0.013866 + t * (0.644237 - 0.013866) * (1.0 - t)).clamp(0.0, 1.0);
772        Vec3::new(r, g, b)
773    }
774
775    /// Turbo colormap (improved rainbow)
776    fn turbo_colormap(&self, t: f32) -> Vec3 {
777        // Simplified Turbo approximation (Google's improved rainbow)
778        let r = if t < 0.5 {
779            (0.13 + 0.87 * (2.0 * t).powf(0.25)).clamp(0.0, 1.0)
780        } else {
781            (0.8685 + 0.1315 * (2.0 * (1.0 - t)).powf(0.25)).clamp(0.0, 1.0)
782        };
783
784        let g = if t < 0.25 {
785            4.0 * t
786        } else if t < 0.75 {
787            1.0
788        } else {
789            1.0 - 4.0 * (t - 0.75)
790        }
791        .clamp(0.0, 1.0);
792
793        let b = if t < 0.5 {
794            (0.8 * (1.0 - 2.0 * t).powf(0.25)).clamp(0.0, 1.0)
795        } else {
796            (0.1 + 0.9 * (2.0 * t - 1.0).powf(0.25)).clamp(0.0, 1.0)
797        };
798
799        Vec3::new(r, g, b)
800    }
801
802    /// Parula colormap (MATLAB's default)
803    fn parula_colormap(&self, t: f32) -> Vec3 {
804        // Simplified Parula approximation
805        let r = if t < 0.25 {
806            0.2081 * (1.0 - t)
807        } else if t < 0.5 {
808            t - 0.25
809        } else if t < 0.75 {
810            1.0
811        } else {
812            1.0 - 0.5 * (t - 0.75)
813        }
814        .clamp(0.0, 1.0);
815
816        let g = if t < 0.125 {
817            0.1663 * t / 0.125
818        } else if t < 0.375 {
819            0.1663 + (0.7079 - 0.1663) * (t - 0.125) / 0.25
820        } else if t < 0.625 {
821            0.7079 + (0.9839 - 0.7079) * (t - 0.375) / 0.25
822        } else {
823            0.9839 * (1.0 - (t - 0.625) / 0.375)
824        }
825        .clamp(0.0, 1.0);
826
827        let b = if t < 0.25 {
828            0.5 + 0.5 * t / 0.25
829        } else if t < 0.5 {
830            1.0
831        } else {
832            1.0 - 2.0 * (t - 0.5)
833        }
834        .clamp(0.0, 1.0);
835
836        Vec3::new(r, g, b)
837    }
838
839    /// Default colormap fallback
840    #[allow(dead_code)] // Fallback method for colormap errors
841    fn default_colormap(&self, t: f32) -> Vec3 {
842        // Use a simple RGB transition as fallback
843        if t < 0.5 {
844            Vec3::new(0.0, 2.0 * t, 1.0 - 2.0 * t)
845        } else {
846            Vec3::new(2.0 * (t - 0.5), 1.0 - 2.0 * (t - 0.5), 0.0)
847        }
848    }
849}
850
851/// MATLAB-compatible surface plot creation utilities
852pub mod matlab_compat {
853    use super::*;
854
855    /// Create a surface plot (equivalent to MATLAB's `surf(X, Y, Z)`)
856    pub fn surf(x: Vec<f64>, y: Vec<f64>, z: Vec<Vec<f64>>) -> Result<SurfacePlot, String> {
857        SurfacePlot::new(x, y, z)
858    }
859
860    /// Create a mesh plot (wireframe surface)
861    pub fn mesh(x: Vec<f64>, y: Vec<f64>, z: Vec<Vec<f64>>) -> Result<SurfacePlot, String> {
862        Ok(SurfacePlot::new(x, y, z)?
863            .with_wireframe(true)
864            .with_shading(ShadingMode::None))
865    }
866
867    /// Create surface from meshgrid
868    pub fn meshgrid_surf(
869        x_range: (f64, f64),
870        y_range: (f64, f64),
871        resolution: (usize, usize),
872        func: impl Fn(f64, f64) -> f64,
873    ) -> Result<SurfacePlot, String> {
874        SurfacePlot::from_function(x_range, y_range, resolution, func)
875    }
876
877    /// Create surface with specific colormap
878    pub fn surf_with_colormap(
879        x: Vec<f64>,
880        y: Vec<f64>,
881        z: Vec<Vec<f64>>,
882        colormap: &str,
883    ) -> Result<SurfacePlot, String> {
884        let cmap = match colormap {
885            "jet" => ColorMap::Jet,
886            "hot" => ColorMap::Hot,
887            "cool" => ColorMap::Cool,
888            "viridis" => ColorMap::Viridis,
889            "plasma" => ColorMap::Plasma,
890            "gray" | "grey" => ColorMap::Gray,
891            _ => return Err(format!("Unknown colormap: {colormap}")),
892        };
893
894        Ok(SurfacePlot::new(x, y, z)?.with_colormap(cmap))
895    }
896}
897
898#[cfg(test)]
899mod tests {
900    use super::*;
901
902    #[test]
903    fn test_surface_plot_creation() {
904        let x = vec![0.0, 1.0, 2.0];
905        let y = vec![0.0, 1.0];
906        let z = vec![vec![0.0, 1.0], vec![1.0, 2.0], vec![2.0, 3.0]];
907
908        let surface = SurfacePlot::new(x, y, z).unwrap();
909
910        assert_eq!(surface.x_data.len(), 3);
911        assert_eq!(surface.y_data.len(), 2);
912        let rows = surface.z_data.as_ref().unwrap();
913        assert_eq!(rows.len(), 3);
914        assert_eq!(rows[0].len(), 2);
915        assert!(surface.visible);
916    }
917
918    #[test]
919    fn test_surface_from_function() {
920        let surface =
921            SurfacePlot::from_function((-2.0, 2.0), (-2.0, 2.0), (10, 10), |x, y| x * x + y * y)
922                .unwrap();
923
924        assert_eq!(surface.x_data.len(), 10);
925        assert_eq!(surface.y_data.len(), 10);
926        let rows = surface.z_data.as_ref().unwrap();
927        assert_eq!(rows.len(), 10);
928
929        // Check that function is evaluated correctly
930        assert_eq!(rows[0][0], 8.0); // (-2)^2 + (-2)^2 = 8
931    }
932
933    #[test]
934    fn test_surface_validation() {
935        let x = vec![0.0, 1.0];
936        let y = vec![0.0, 1.0, 2.0];
937        let z = vec![
938            vec![0.0, 1.0], // Wrong: should have 3 elements to match y
939            vec![1.0, 2.0],
940        ];
941
942        assert!(SurfacePlot::new(x, y, z).is_err());
943    }
944
945    #[test]
946    fn test_surface_styling() {
947        let x = vec![0.0, 1.0];
948        let y = vec![0.0, 1.0];
949        let z = vec![vec![0.0, 1.0], vec![1.0, 2.0]];
950
951        let surface = SurfacePlot::new(x, y, z)
952            .unwrap()
953            .with_colormap(ColorMap::Hot)
954            .with_wireframe(true)
955            .with_alpha(0.8)
956            .with_label("Test Surface");
957
958        assert_eq!(surface.colormap, ColorMap::Hot);
959        assert!(surface.wireframe);
960        assert_eq!(surface.alpha, 0.8);
961        assert_eq!(surface.label, Some("Test Surface".to_string()));
962    }
963
964    #[test]
965    fn test_colormap_mapping() {
966        let jet = ColorMap::Jet;
967
968        // Test boundary values
969        let color_0 = jet.map_value(0.0);
970        let color_1 = jet.map_value(1.0);
971
972        assert!(color_0.x >= 0.0 && color_0.x <= 1.0);
973        assert!(color_1.x >= 0.0 && color_1.x <= 1.0);
974
975        // Test that different values give different colors
976        let color_mid = jet.map_value(0.5);
977        assert_ne!(color_0, color_mid);
978        assert_ne!(color_mid, color_1);
979    }
980
981    #[test]
982    fn test_surface_statistics() {
983        let x = vec![0.0, 1.0, 2.0, 3.0];
984        let y = vec![0.0, 1.0, 2.0];
985        let z = vec![
986            vec![0.0, 1.0, 2.0],
987            vec![1.0, 2.0, 3.0],
988            vec![2.0, 3.0, 4.0],
989            vec![3.0, 4.0, 5.0],
990        ];
991
992        let surface = SurfacePlot::new(x, y, z).unwrap();
993        let stats = surface.statistics();
994
995        assert_eq!(stats.grid_points, 12); // 4 * 3
996        assert_eq!(stats.triangle_count, 12); // (4-1) * (3-1) * 2
997        assert_eq!(stats.x_resolution, 4);
998        assert_eq!(stats.y_resolution, 3);
999        assert!(stats.memory_usage > 0);
1000    }
1001
1002    #[test]
1003    fn test_matlab_compat() {
1004        use super::matlab_compat::*;
1005
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 = surf(x.clone(), y.clone(), z.clone()).unwrap();
1011        assert!(!surface.wireframe);
1012
1013        let mesh_plot = mesh(x.clone(), y.clone(), z.clone()).unwrap();
1014        assert!(mesh_plot.wireframe);
1015
1016        let colormap_surface = surf_with_colormap(x, y, z, "viridis").unwrap();
1017        assert_eq!(colormap_surface.colormap, ColorMap::Viridis);
1018    }
1019}
runmat_plot/plots/surface.rs

runmat_plot/plots/
surface.rs
