proof_engine/surfaces/
parametric.rs

1//! # Parametric Surface Rendering
2//!
3//! Mathematical surfaces evaluated over (u, v) parameter space, tessellated into
4//! triangle meshes suitable for glyph-grid rendering.
5//!
6//! Every surface implements the [`Surface`] trait, which provides `evaluate(u, v) -> Vec3`
7//! and `normal(u, v) -> Vec3`. The [`SurfaceMesh`] struct tessellates any surface into
8//! a triangle mesh with computed normals, tangents, and UV coordinates.
9//!
10//! ## Included surfaces
11//!
12//! - Sphere, Torus, MobiusStrip, KleinBottle
13//! - BoySurface, RomanSurface, CrossCap
14//! - TrefoilKnot, FigureEight
15//! - Catenoid, Helicoid, EnneperSurface, DiniSurface
16//! - FunctionSurface (user-defined closure)
17
18use glam::Vec3;
19use std::f32::consts::{PI, TAU};
20
21// ─────────────────────────────────────────────────────────────────────────────
22// Surface trait
23// ─────────────────────────────────────────────────────────────────────────────
24
25/// A parametric surface defined over (u, v) in [0, 1] x [0, 1].
26pub trait Surface {
27    /// Evaluate the surface position at parameter (u, v).
28    fn evaluate(&self, u: f32, v: f32) -> Vec3;
29
30    /// Compute the surface normal at parameter (u, v).
31    /// Default implementation uses central differences.
32    fn normal(&self, u: f32, v: f32) -> Vec3 {
33        let eps = 1e-4;
34        let du = self.evaluate((u + eps).min(1.0), v) - self.evaluate((u - eps).max(0.0), v);
35        let dv = self.evaluate(u, (v + eps).min(1.0)) - self.evaluate(u, (v - eps).max(0.0));
36        du.cross(dv).normalize_or_zero()
37    }
38
39    /// Compute the partial derivative with respect to u at (u, v).
40    fn tangent_u(&self, u: f32, v: f32) -> Vec3 {
41        let eps = 1e-4;
42        let forward = self.evaluate((u + eps).min(1.0), v);
43        let backward = self.evaluate((u - eps).max(0.0), v);
44        (forward - backward).normalize_or_zero()
45    }
46
47    /// Compute the partial derivative with respect to v at (u, v).
48    fn tangent_v(&self, u: f32, v: f32) -> Vec3 {
49        let eps = 1e-4;
50        let forward = self.evaluate(u, (v + eps).min(1.0));
51        let backward = self.evaluate(u, (v - eps).max(0.0));
52        (forward - backward).normalize_or_zero()
53    }
54
55    /// Parameter domain for u: (min, max). Defaults to (0, 1).
56    fn u_range(&self) -> (f32, f32) { (0.0, 1.0) }
57
58    /// Parameter domain for v: (min, max). Defaults to (0, 1).
59    fn v_range(&self) -> (f32, f32) { (0.0, 1.0) }
60
61    /// Whether the surface wraps in u. Used for seamless tessellation.
62    fn wraps_u(&self) -> bool { false }
63
64    /// Whether the surface wraps in v. Used for seamless tessellation.
65    fn wraps_v(&self) -> bool { false }
66
67    /// Human-readable name for debugging.
68    fn name(&self) -> &str { "Surface" }
69}
70
71// ─────────────────────────────────────────────────────────────────────────────
72// Sphere
73// ─────────────────────────────────────────────────────────────────────────────
74
75/// A unit sphere centered at the origin.
76#[derive(Debug, Clone)]
77pub struct Sphere {
78    pub radius: f32,
79    pub center: Vec3,
80}
81
82impl Sphere {
83    pub fn new(radius: f32) -> Self {
84        Self { radius, center: Vec3::ZERO }
85    }
86
87    pub fn with_center(radius: f32, center: Vec3) -> Self {
88        Self { radius, center }
89    }
90}
91
92impl Default for Sphere {
93    fn default() -> Self { Self::new(1.0) }
94}
95
96impl Surface for Sphere {
97    fn evaluate(&self, u: f32, v: f32) -> Vec3 {
98        let theta = u * TAU;
99        let phi = v * PI;
100        let sp = phi.sin();
101        self.center + Vec3::new(
102            self.radius * sp * theta.cos(),
103            self.radius * phi.cos(),
104            self.radius * sp * theta.sin(),
105        )
106    }
107
108    fn normal(&self, u: f32, v: f32) -> Vec3 {
109        let theta = u * TAU;
110        let phi = v * PI;
111        let sp = phi.sin();
112        Vec3::new(sp * theta.cos(), phi.cos(), sp * theta.sin()).normalize_or_zero()
113    }
114
115    fn wraps_u(&self) -> bool { true }
116    fn name(&self) -> &str { "Sphere" }
117}
118
119// ─────────────────────────────────────────────────────────────────────────────
120// Torus
121// ─────────────────────────────────────────────────────────────────────────────
122
123/// A torus with major radius R and minor radius r.
124#[derive(Debug, Clone)]
125pub struct Torus {
126    pub major_radius: f32,
127    pub minor_radius: f32,
128    pub center: Vec3,
129}
130
131impl Torus {
132    pub fn new(major_radius: f32, minor_radius: f32) -> Self {
133        Self { major_radius, minor_radius, center: Vec3::ZERO }
134    }
135}
136
137impl Default for Torus {
138    fn default() -> Self { Self::new(1.0, 0.3) }
139}
140
141impl Surface for Torus {
142    fn evaluate(&self, u: f32, v: f32) -> Vec3 {
143        let theta = u * TAU;
144        let phi = v * TAU;
145        let r = self.major_radius + self.minor_radius * phi.cos();
146        self.center + Vec3::new(
147            r * theta.cos(),
148            self.minor_radius * phi.sin(),
149            r * theta.sin(),
150        )
151    }
152
153    fn normal(&self, u: f32, v: f32) -> Vec3 {
154        let theta = u * TAU;
155        let phi = v * TAU;
156        Vec3::new(
157            phi.cos() * theta.cos(),
158            phi.sin(),
159            phi.cos() * theta.sin(),
160        ).normalize_or_zero()
161    }
162
163    fn wraps_u(&self) -> bool { true }
164    fn wraps_v(&self) -> bool { true }
165    fn name(&self) -> &str { "Torus" }
166}
167
168// ─────────────────────────────────────────────────────────────────────────────
169// Mobius Strip
170// ─────────────────────────────────────────────────────────────────────────────
171
172/// A Mobius strip — a surface with only one side.
173#[derive(Debug, Clone)]
174pub struct MobiusStrip {
175    pub radius: f32,
176    pub width: f32,
177}
178
179impl MobiusStrip {
180    pub fn new(radius: f32, width: f32) -> Self {
181        Self { radius, width }
182    }
183}
184
185impl Default for MobiusStrip {
186    fn default() -> Self { Self::new(1.0, 0.4) }
187}
188
189impl Surface for MobiusStrip {
190    fn evaluate(&self, u: f32, v: f32) -> Vec3 {
191        let theta = u * TAU;
192        let s = (v - 0.5) * self.width;
193        let half_theta = theta * 0.5;
194        let r = self.radius + s * half_theta.cos();
195        Vec3::new(
196            r * theta.cos(),
197            s * half_theta.sin(),
198            r * theta.sin(),
199        )
200    }
201
202    fn wraps_u(&self) -> bool { true }
203    fn name(&self) -> &str { "MobiusStrip" }
204}
205
206// ─────────────────────────────────────────────────────────────────────────────
207// Klein Bottle
208// ─────────────────────────────────────────────────────────────────────────────
209
210/// The Klein bottle — a non-orientable closed surface.
211/// Uses the "figure-8" immersion in R^3.
212#[derive(Debug, Clone)]
213pub struct KleinBottle {
214    pub scale: f32,
215}
216
217impl KleinBottle {
218    pub fn new(scale: f32) -> Self { Self { scale } }
219}
220
221impl Default for KleinBottle {
222    fn default() -> Self { Self::new(1.0) }
223}
224
225impl Surface for KleinBottle {
226    fn evaluate(&self, u: f32, v: f32) -> Vec3 {
227        let theta = u * TAU;
228        let phi = v * TAU;
229        let s = self.scale;
230        // Figure-8 Klein bottle immersion
231        let r = 4.0 * (1.0 - 0.5 * theta.sin());
232        let x = s * (6.0 * theta.cos() * (1.0 + theta.sin())
233                     + r * theta.cos() * phi.cos());
234        let y = s * (16.0 * theta.sin()
235                     + r * theta.sin() * phi.cos());
236        let z = s * r * phi.sin();
237        Vec3::new(x * 0.05, y * 0.05, z * 0.05)
238    }
239
240    fn wraps_u(&self) -> bool { true }
241    fn wraps_v(&self) -> bool { true }
242    fn name(&self) -> &str { "KleinBottle" }
243}
244
245// ─────────────────────────────────────────────────────────────────────────────
246// Boy's Surface
247// ─────────────────────────────────────────────────────────────────────────────
248
249/// Boy's surface — an immersion of the real projective plane in R^3.
250#[derive(Debug, Clone)]
251pub struct BoySurface {
252    pub scale: f32,
253}
254
255impl BoySurface {
256    pub fn new(scale: f32) -> Self { Self { scale } }
257}
258
259impl Default for BoySurface {
260    fn default() -> Self { Self::new(1.0) }
261}
262
263impl Surface for BoySurface {
264    fn evaluate(&self, u: f32, v: f32) -> Vec3 {
265        // Bryant-Kusner parametrization via spherical coords
266        let theta = u * PI;
267        let phi = v * TAU;
268        let st = theta.sin();
269        let ct = theta.cos();
270        let sp = phi.sin();
271        let cp = phi.cos();
272
273        let denom = 2.0_f32.sqrt();
274        let x1 = st * cp;
275        let x2 = st * sp;
276        let x3 = ct;
277
278        // Boy surface parametrization (Apery)
279        let g0 = x1 * x1 - x2 * x2;
280        let g1 = x1 * x2;
281        let g2 = x2 * x3;
282        let g3 = x3 * x1;
283
284        let px = denom * (2.0 * g3 + g0) / (x1 * x1 + x2 * x2 + x3 * x3 + 1.0);
285        let py = denom * (2.0 * g2 + g1) / (x1 * x1 + x2 * x2 + x3 * x3 + 1.0);
286        let pz = (3.0 * x3 * x3 - 1.0) / (2.0 * (x1 * x1 + x2 * x2 + x3 * x3 + 1.0));
287
288        Vec3::new(px, py, pz) * self.scale
289    }
290
291    fn name(&self) -> &str { "BoySurface" }
292}
293
294// ─────────────────────────────────────────────────────────────────────────────
295// Roman Surface (Steiner surface)
296// ─────────────────────────────────────────────────────────────────────────────
297
298/// Steiner's Roman surface — a self-intersecting immersion of the real projective plane.
299#[derive(Debug, Clone)]
300pub struct RomanSurface {
301    pub scale: f32,
302}
303
304impl RomanSurface {
305    pub fn new(scale: f32) -> Self { Self { scale } }
306}
307
308impl Default for RomanSurface {
309    fn default() -> Self { Self::new(1.0) }
310}
311
312impl Surface for RomanSurface {
313    fn evaluate(&self, u: f32, v: f32) -> Vec3 {
314        let theta = u * PI;
315        let phi = v * TAU;
316        let st = theta.sin();
317        let ct = theta.cos();
318        let sp = phi.sin();
319        let cp = phi.cos();
320        let s = self.scale;
321        Vec3::new(
322            s * st * st * sp * cp,
323            s * st * cp * ct,
324            s * st * sp * ct,
325        )
326    }
327
328    fn name(&self) -> &str { "RomanSurface" }
329}
330
331// ─────────────────────────────────────────────────────────────────────────────
332// Cross-Cap
333// ─────────────────────────────────────────────────────────────────────────────
334
335/// Cross-cap — another immersion of the real projective plane in R^3.
336#[derive(Debug, Clone)]
337pub struct CrossCap {
338    pub scale: f32,
339}
340
341impl CrossCap {
342    pub fn new(scale: f32) -> Self { Self { scale } }
343}
344
345impl Default for CrossCap {
346    fn default() -> Self { Self::new(1.0) }
347}
348
349impl Surface for CrossCap {
350    fn evaluate(&self, u: f32, v: f32) -> Vec3 {
351        let theta = u * PI;
352        let phi = v * TAU;
353        let st = theta.sin();
354        let ct = theta.cos();
355        let sp = phi.sin();
356        let cp = phi.cos();
357        let s = self.scale;
358        Vec3::new(
359            s * 0.5 * st * (2.0 * phi).sin(),
360            s * st * cp,
361            s * ct,
362        )
363    }
364
365    fn name(&self) -> &str { "CrossCap" }
366}
367
368// ─────────────────────────────────────────────────────────────────────────────
369// Trefoil Knot (tube surface around a trefoil curve)
370// ─────────────────────────────────────────────────────────────────────────────
371
372/// Trefoil knot — a tube surface of given radius wrapped around the trefoil curve.
373#[derive(Debug, Clone)]
374pub struct TrefoilKnot {
375    pub tube_radius: f32,
376    pub scale: f32,
377}
378
379impl TrefoilKnot {
380    pub fn new(tube_radius: f32, scale: f32) -> Self {
381        Self { tube_radius, scale }
382    }
383}
384
385impl Default for TrefoilKnot {
386    fn default() -> Self { Self::new(0.15, 1.0) }
387}
388
389impl TrefoilKnot {
390    /// Evaluate the trefoil centerline at parameter t in [0, TAU].
391    fn curve(&self, t: f32) -> Vec3 {
392        let s = self.scale;
393        Vec3::new(
394            s * (t.sin() + 2.0 * (2.0 * t).sin()),
395            s * (t.cos() - 2.0 * (2.0 * t).cos()),
396            s * (-(3.0 * t).sin()),
397        )
398    }
399
400    /// Approximate tangent by central differences.
401    fn curve_tangent(&self, t: f32) -> Vec3 {
402        let eps = 1e-4;
403        (self.curve(t + eps) - self.curve(t - eps)).normalize_or_zero()
404    }
405}
406
407impl Surface for TrefoilKnot {
408    fn evaluate(&self, u: f32, v: f32) -> Vec3 {
409        let t = u * TAU;
410        let phi = v * TAU;
411        let center = self.curve(t);
412        let tangent = self.curve_tangent(t);
413
414        // Build a Frenet frame
415        let up = if tangent.y.abs() < 0.99 { Vec3::Y } else { Vec3::X };
416        let normal = tangent.cross(up).normalize_or_zero();
417        let binormal = tangent.cross(normal).normalize_or_zero();
418
419        center + self.tube_radius * (phi.cos() * normal + phi.sin() * binormal)
420    }
421
422    fn wraps_u(&self) -> bool { true }
423    fn wraps_v(&self) -> bool { true }
424    fn name(&self) -> &str { "TrefoilKnot" }
425}
426
427// ─────────────────────────────────────────────────────────────────────────────
428// Figure-Eight Knot (tube surface)
429// ─────────────────────────────────────────────────────────────────────────────
430
431/// Figure-eight knot — a tube surface around the figure-eight curve.
432#[derive(Debug, Clone)]
433pub struct FigureEight {
434    pub tube_radius: f32,
435    pub scale: f32,
436}
437
438impl FigureEight {
439    pub fn new(tube_radius: f32, scale: f32) -> Self {
440        Self { tube_radius, scale }
441    }
442}
443
444impl Default for FigureEight {
445    fn default() -> Self { Self::new(0.12, 1.0) }
446}
447
448impl FigureEight {
449    fn curve(&self, t: f32) -> Vec3 {
450        let s = self.scale;
451        let ct = t.cos();
452        let st = t.sin();
453        let s2t = (2.0 * t).sin();
454        Vec3::new(
455            s * (2.0 + ct) * (3.0 * t).cos(),
456            s * (2.0 + ct) * (3.0 * t).sin(),
457            s * st * 3.0,
458        )
459    }
460
461    fn curve_tangent(&self, t: f32) -> Vec3 {
462        let eps = 1e-4;
463        (self.curve(t + eps) - self.curve(t - eps)).normalize_or_zero()
464    }
465}
466
467impl Surface for FigureEight {
468    fn evaluate(&self, u: f32, v: f32) -> Vec3 {
469        let t = u * TAU;
470        let phi = v * TAU;
471        let center = self.curve(t);
472        let tangent = self.curve_tangent(t);
473
474        let up = if tangent.y.abs() < 0.99 { Vec3::Y } else { Vec3::X };
475        let normal = tangent.cross(up).normalize_or_zero();
476        let binormal = tangent.cross(normal).normalize_or_zero();
477
478        center + self.tube_radius * (phi.cos() * normal + phi.sin() * binormal)
479    }
480
481    fn wraps_u(&self) -> bool { true }
482    fn wraps_v(&self) -> bool { true }
483    fn name(&self) -> &str { "FigureEight" }
484}
485
486// ─────────────────────────────────────────────────────────────────────────────
487// Catenoid
488// ─────────────────────────────────────────────────────────────────────────────
489
490/// A catenoid — a minimal surface of revolution formed by a catenary.
491#[derive(Debug, Clone)]
492pub struct Catenoid {
493    pub c: f32,
494    pub height_range: f32,
495}
496
497impl Catenoid {
498    pub fn new(c: f32, height_range: f32) -> Self {
499        Self { c, height_range }
500    }
501}
502
503impl Default for Catenoid {
504    fn default() -> Self { Self::new(1.0, 2.0) }
505}
506
507impl Surface for Catenoid {
508    fn evaluate(&self, u: f32, v: f32) -> Vec3 {
509        let theta = u * TAU;
510        let t = (v - 0.5) * self.height_range;
511        let ch = (t / self.c).cosh();
512        Vec3::new(
513            self.c * ch * theta.cos(),
514            t,
515            self.c * ch * theta.sin(),
516        )
517    }
518
519    fn normal(&self, u: f32, v: f32) -> Vec3 {
520        let theta = u * TAU;
521        let t = (v - 0.5) * self.height_range;
522        let sh = (t / self.c).sinh();
523        let ch = (t / self.c).cosh();
524        Vec3::new(
525            -theta.cos() / ch,
526            sh / ch,
527            -theta.sin() / ch,
528        ).normalize_or_zero()
529    }
530
531    fn wraps_u(&self) -> bool { true }
532    fn name(&self) -> &str { "Catenoid" }
533}
534
535// ─────────────────────────────────────────────────────────────────────────────
536// Helicoid
537// ─────────────────────────────────────────────────────────────────────────────
538
539/// A helicoid — a minimal surface swept by a line rotating about and translating along an axis.
540#[derive(Debug, Clone)]
541pub struct Helicoid {
542    pub pitch: f32,
543    pub radius: f32,
544    pub height_range: f32,
545}
546
547impl Helicoid {
548    pub fn new(pitch: f32, radius: f32, height_range: f32) -> Self {
549        Self { pitch, radius, height_range }
550    }
551}
552
553impl Default for Helicoid {
554    fn default() -> Self { Self::new(1.0, 1.0, 4.0) }
555}
556
557impl Surface for Helicoid {
558    fn evaluate(&self, u: f32, v: f32) -> Vec3 {
559        let theta = u * TAU * 2.0; // two full turns
560        let r = (v - 0.5) * 2.0 * self.radius;
561        Vec3::new(
562            r * theta.cos(),
563            self.pitch * theta / TAU * self.height_range * 0.5,
564            r * theta.sin(),
565        )
566    }
567
568    fn normal(&self, u: f32, v: f32) -> Vec3 {
569        let theta = u * TAU * 2.0;
570        let denom = (self.pitch * self.pitch + ((v - 0.5) * 2.0 * self.radius).powi(2)).sqrt();
571        if denom < 1e-8 {
572            return Vec3::Y;
573        }
574        Vec3::new(
575            self.pitch * theta.sin() / denom,
576            -((v - 0.5) * 2.0 * self.radius) / denom,
577            -self.pitch * theta.cos() / denom,
578        ).normalize_or_zero()
579    }
580
581    fn name(&self) -> &str { "Helicoid" }
582}
583
584// ─────────────────────────────────────────────────────────────────────────────
585// Enneper Surface
586// ─────────────────────────────────────────────────────────────────────────────
587
588/// Enneper's minimal surface — self-intersects at higher ranges.
589#[derive(Debug, Clone)]
590pub struct EnneperSurface {
591    pub scale: f32,
592    pub range: f32,
593}
594
595impl EnneperSurface {
596    pub fn new(scale: f32, range: f32) -> Self {
597        Self { scale, range }
598    }
599}
600
601impl Default for EnneperSurface {
602    fn default() -> Self { Self::new(0.3, 2.0) }
603}
604
605impl Surface for EnneperSurface {
606    fn evaluate(&self, u: f32, v: f32) -> Vec3 {
607        let uu = (u - 0.5) * 2.0 * self.range;
608        let vv = (v - 0.5) * 2.0 * self.range;
609        let s = self.scale;
610        Vec3::new(
611            s * (uu - uu.powi(3) / 3.0 + uu * vv * vv),
612            s * (vv - vv.powi(3) / 3.0 + vv * uu * uu),
613            s * (uu * uu - vv * vv),
614        )
615    }
616
617    fn name(&self) -> &str { "EnneperSurface" }
618}
619
620// ─────────────────────────────────────────────────────────────────────────────
621// Dini's Surface
622// ─────────────────────────────────────────────────────────────────────────────
623
624/// Dini's surface — a twisted pseudosphere with constant negative curvature.
625#[derive(Debug, Clone)]
626pub struct DiniSurface {
627    pub a: f32,
628    pub b: f32,
629    pub turns: f32,
630}
631
632impl DiniSurface {
633    pub fn new(a: f32, b: f32, turns: f32) -> Self {
634        Self { a, b, turns }
635    }
636}
637
638impl Default for DiniSurface {
639    fn default() -> Self { Self::new(1.0, 0.2, 4.0) }
640}
641
642impl Surface for DiniSurface {
643    fn evaluate(&self, u: f32, v: f32) -> Vec3 {
644        let theta = u * TAU * self.turns;
645        // v in (0, PI) but avoid singularity at 0
646        let phi = v * (PI - 0.02) + 0.01;
647        let a = self.a;
648        let b = self.b;
649        Vec3::new(
650            a * theta.cos() * phi.sin(),
651            a * theta.sin() * phi.sin(),
652            a * (phi.cos() + (phi / 2.0).tan().max(-100.0).min(100.0).ln()) + b * theta,
653        )
654    }
655
656    fn name(&self) -> &str { "DiniSurface" }
657}
658
659// ─────────────────────────────────────────────────────────────────────────────
660// FunctionSurface (user-defined)
661// ─────────────────────────────────────────────────────────────────────────────
662
663/// A surface defined by a user-provided closure.
664pub struct FunctionSurface {
665    pub func: Box<dyn Fn(f32, f32) -> Vec3 + Send + Sync>,
666    label: String,
667    wrap_u: bool,
668    wrap_v: bool,
669}
670
671impl FunctionSurface {
672    pub fn new<F>(func: F) -> Self
673    where
674        F: Fn(f32, f32) -> Vec3 + Send + Sync + 'static,
675    {
676        Self {
677            func: Box::new(func),
678            label: "FunctionSurface".to_string(),
679            wrap_u: false,
680            wrap_v: false,
681        }
682    }
683
684    pub fn with_name(mut self, name: &str) -> Self {
685        self.label = name.to_string();
686        self
687    }
688
689    pub fn with_wrapping(mut self, wrap_u: bool, wrap_v: bool) -> Self {
690        self.wrap_u = wrap_u;
691        self.wrap_v = wrap_v;
692        self
693    }
694}
695
696impl Surface for FunctionSurface {
697    fn evaluate(&self, u: f32, v: f32) -> Vec3 {
698        (self.func)(u, v)
699    }
700
701    fn wraps_u(&self) -> bool { self.wrap_u }
702    fn wraps_v(&self) -> bool { self.wrap_v }
703    fn name(&self) -> &str { &self.label }
704}
705
706// ─────────────────────────────────────────────────────────────────────────────
707// Vertex
708// ─────────────────────────────────────────────────────────────────────────────
709
710/// A vertex in a surface mesh.
711#[derive(Debug, Clone, Copy)]
712pub struct SurfaceVertex {
713    pub position: Vec3,
714    pub normal: Vec3,
715    pub tangent: Vec3,
716    pub bitangent: Vec3,
717    pub uv: [f32; 2],
718    /// Parameter-space coordinates that produced this vertex.
719    pub param_u: f32,
720    pub param_v: f32,
721}
722
723impl Default for SurfaceVertex {
724    fn default() -> Self {
725        Self {
726            position: Vec3::ZERO,
727            normal: Vec3::Y,
728            tangent: Vec3::X,
729            bitangent: Vec3::Z,
730            uv: [0.0, 0.0],
731            param_u: 0.0,
732            param_v: 0.0,
733        }
734    }
735}
736
737/// A triangle index triple.
738#[derive(Debug, Clone, Copy)]
739pub struct TriangleIndex {
740    pub a: u32,
741    pub b: u32,
742    pub c: u32,
743}
744
745impl TriangleIndex {
746    pub fn new(a: u32, b: u32, c: u32) -> Self { Self { a, b, c } }
747}
748
749// ─────────────────────────────────────────────────────────────────────────────
750// SurfaceMesh
751// ─────────────────────────────────────────────────────────────────────────────
752
753/// A tessellated mesh of a parametric surface.
754#[derive(Debug, Clone)]
755pub struct SurfaceMesh {
756    pub vertices: Vec<SurfaceVertex>,
757    pub indices: Vec<TriangleIndex>,
758    pub resolution_u: usize,
759    pub resolution_v: usize,
760    /// Axis-aligned bounding box: (min, max).
761    pub aabb_min: Vec3,
762    pub aabb_max: Vec3,
763}
764
765impl SurfaceMesh {
766    /// Tessellate a surface at the given resolution.
767    ///
768    /// `resolution_u` and `resolution_v` define the number of subdivisions along
769    /// each parameter axis. The resulting mesh has `(res_u + 1) * (res_v + 1)` vertices
770    /// and `2 * res_u * res_v` triangles.
771    pub fn tessellate(surface: &dyn Surface, resolution_u: usize, resolution_v: usize) -> Self {
772        let res_u = resolution_u.max(2);
773        let res_v = resolution_v.max(2);
774        let num_verts = (res_u + 1) * (res_v + 1);
775        let mut vertices = Vec::with_capacity(num_verts);
776        let mut aabb_min = Vec3::splat(f32::MAX);
777        let mut aabb_max = Vec3::splat(f32::MIN);
778
779        // Generate vertices
780        for j in 0..=res_v {
781            let v = j as f32 / res_v as f32;
782            for i in 0..=res_u {
783                let u = i as f32 / res_u as f32;
784                let pos = surface.evaluate(u, v);
785                let norm = surface.normal(u, v);
786                let tan_u = surface.tangent_u(u, v);
787                let bitan = norm.cross(tan_u).normalize_or_zero();
788
789                aabb_min = aabb_min.min(pos);
790                aabb_max = aabb_max.max(pos);
791
792                vertices.push(SurfaceVertex {
793                    position: pos,
794                    normal: norm,
795                    tangent: tan_u,
796                    bitangent: bitan,
797                    uv: [u, v],
798                    param_u: u,
799                    param_v: v,
800                });
801            }
802        }
803
804        // Generate triangle indices
805        let num_triangles = res_u * res_v * 2;
806        let mut indices = Vec::with_capacity(num_triangles);
807        for j in 0..res_v {
808            for i in 0..res_u {
809                let tl = (j * (res_u + 1) + i) as u32;
810                let tr = tl + 1;
811                let bl = tl + (res_u + 1) as u32;
812                let br = bl + 1;
813                indices.push(TriangleIndex::new(tl, bl, tr));
814                indices.push(TriangleIndex::new(tr, bl, br));
815            }
816        }
817
818        let mut mesh = Self {
819            vertices,
820            indices,
821            resolution_u: res_u,
822            resolution_v: res_v,
823            aabb_min,
824            aabb_max,
825        };
826
827        // Recompute smooth normals by averaging face normals at shared vertices
828        mesh.recompute_smooth_normals();
829        mesh
830    }
831
832    /// Tessellate with adaptive resolution that places more samples where curvature is higher.
833    pub fn tessellate_adaptive(
834        surface: &dyn Surface,
835        base_resolution_u: usize,
836        base_resolution_v: usize,
837        curvature_threshold: f32,
838    ) -> Self {
839        // Sample curvature to find high-curvature regions
840        let sample_res = 16.min(base_resolution_u).max(4);
841        let mut max_curvature = 0.0_f32;
842        let mut curvature_map: Vec<f32> = Vec::with_capacity(sample_res * sample_res);
843
844        for j in 0..sample_res {
845            let v = j as f32 / sample_res as f32;
846            for i in 0..sample_res {
847                let u = i as f32 / sample_res as f32;
848                let c = Self::estimate_curvature(surface, u, v);
849                max_curvature = max_curvature.max(c);
850                curvature_map.push(c);
851            }
852        }
853
854        // For now, if curvature is uniformly low, use base resolution.
855        // Otherwise, double it.
856        let factor = if max_curvature > curvature_threshold { 2 } else { 1 };
857        Self::tessellate(surface, base_resolution_u * factor, base_resolution_v * factor)
858    }
859
860    /// Estimate Gaussian curvature at a point using finite differences.
861    fn estimate_curvature(surface: &dyn Surface, u: f32, v: f32) -> f32 {
862        let eps = 1e-3;
863        let n0 = surface.normal(u, v);
864        let nu = surface.normal((u + eps).min(1.0), v);
865        let nv = surface.normal(u, (v + eps).min(1.0));
866        let du = (nu - n0).length() / eps;
867        let dv = (nv - n0).length() / eps;
868        (du * du + dv * dv).sqrt()
869    }
870
871    /// Recompute vertex normals by averaging adjacent face normals.
872    pub fn recompute_smooth_normals(&mut self) {
873        // Zero all normals
874        for v in &mut self.vertices {
875            v.normal = Vec3::ZERO;
876        }
877
878        // Accumulate face normals
879        for tri in &self.indices {
880            let a = tri.a as usize;
881            let b = tri.b as usize;
882            let c = tri.c as usize;
883            let pa = self.vertices[a].position;
884            let pb = self.vertices[b].position;
885            let pc = self.vertices[c].position;
886            let face_normal = (pb - pa).cross(pc - pa);
887            // Weight by face area (non-normalized cross product = 2 * area)
888            self.vertices[a].normal += face_normal;
889            self.vertices[b].normal += face_normal;
890            self.vertices[c].normal += face_normal;
891        }
892
893        // Normalize and recompute tangent frames
894        for v in &mut self.vertices {
895            v.normal = v.normal.normalize_or_zero();
896            // Recompute tangent frame from normal
897            let up = if v.normal.y.abs() < 0.99 { Vec3::Y } else { Vec3::X };
898            v.tangent = v.normal.cross(up).normalize_or_zero();
899            v.bitangent = v.normal.cross(v.tangent).normalize_or_zero();
900        }
901    }
902
903    /// Recompute tangent vectors using UV-aligned tangent calculation (MikkTSpace-style).
904    pub fn recompute_tangents(&mut self) {
905        // Zero tangents and bitangents
906        for v in &mut self.vertices {
907            v.tangent = Vec3::ZERO;
908            v.bitangent = Vec3::ZERO;
909        }
910
911        for tri in &self.indices {
912            let i0 = tri.a as usize;
913            let i1 = tri.b as usize;
914            let i2 = tri.c as usize;
915
916            let p0 = self.vertices[i0].position;
917            let p1 = self.vertices[i1].position;
918            let p2 = self.vertices[i2].position;
919
920            let uv0 = self.vertices[i0].uv;
921            let uv1 = self.vertices[i1].uv;
922            let uv2 = self.vertices[i2].uv;
923
924            let dp1 = p1 - p0;
925            let dp2 = p2 - p0;
926            let duv1 = [uv1[0] - uv0[0], uv1[1] - uv0[1]];
927            let duv2 = [uv2[0] - uv0[0], uv2[1] - uv0[1]];
928
929            let det = duv1[0] * duv2[1] - duv1[1] * duv2[0];
930            if det.abs() < 1e-8 {
931                continue;
932            }
933            let inv_det = 1.0 / det;
934
935            let tangent = (dp1 * duv2[1] - dp2 * duv1[1]) * inv_det;
936            let bitangent = (dp2 * duv1[0] - dp1 * duv2[0]) * inv_det;
937
938            self.vertices[i0].tangent += tangent;
939            self.vertices[i1].tangent += tangent;
940            self.vertices[i2].tangent += tangent;
941
942            self.vertices[i0].bitangent += bitangent;
943            self.vertices[i1].bitangent += bitangent;
944            self.vertices[i2].bitangent += bitangent;
945        }
946
947        // Gram-Schmidt orthonormalize
948        for v in &mut self.vertices {
949            let n = v.normal;
950            let t = v.tangent;
951            let b = v.bitangent;
952
953            // Orthonormalize tangent
954            let t_orth = (t - n * n.dot(t)).normalize_or_zero();
955            // Compute bitangent with correct handedness
956            let b_sign = if n.cross(t).dot(b) < 0.0 { -1.0 } else { 1.0 };
957            let b_orth = n.cross(t_orth) * b_sign;
958
959            v.tangent = t_orth;
960            v.bitangent = b_orth;
961        }
962    }
963
964    /// Transform all vertices by a 4x4-like operation (scale + translate).
965    pub fn transform(&mut self, scale: f32, offset: Vec3) {
966        self.aabb_min = Vec3::splat(f32::MAX);
967        self.aabb_max = Vec3::splat(f32::MIN);
968        for v in &mut self.vertices {
969            v.position = v.position * scale + offset;
970            self.aabb_min = self.aabb_min.min(v.position);
971            self.aabb_max = self.aabb_max.max(v.position);
972        }
973    }
974
975    /// Return the center of the bounding box.
976    pub fn center(&self) -> Vec3 {
977        (self.aabb_min + self.aabb_max) * 0.5
978    }
979
980    /// Return the extent (half-size) of the bounding box.
981    pub fn extent(&self) -> Vec3 {
982        (self.aabb_max - self.aabb_min) * 0.5
983    }
984
985    /// Total number of triangles in the mesh.
986    pub fn triangle_count(&self) -> usize {
987        self.indices.len()
988    }
989
990    /// Total number of vertices in the mesh.
991    pub fn vertex_count(&self) -> usize {
992        self.vertices.len()
993    }
994
995    /// Compute the surface area of the mesh by summing triangle areas.
996    pub fn surface_area(&self) -> f32 {
997        let mut area = 0.0_f32;
998        for tri in &self.indices {
999            let a = self.vertices[tri.a as usize].position;
1000            let b = self.vertices[tri.b as usize].position;
1001            let c = self.vertices[tri.c as usize].position;
1002            area += (b - a).cross(c - a).length() * 0.5;
1003        }
1004        area
1005    }
1006
1007    /// Generate flat normals (each face gets its own normal, vertices are duplicated).
1008    pub fn make_flat_shaded(&mut self) {
1009        let mut new_verts = Vec::with_capacity(self.indices.len() * 3);
1010        let mut new_indices = Vec::with_capacity(self.indices.len());
1011
1012        for tri in &self.indices {
1013            let va = self.vertices[tri.a as usize];
1014            let vb = self.vertices[tri.b as usize];
1015            let vc = self.vertices[tri.c as usize];
1016
1017            let face_normal = (vb.position - va.position)
1018                .cross(vc.position - va.position)
1019                .normalize_or_zero();
1020
1021            let base = new_verts.len() as u32;
1022
1023            let mut a = va;
1024            let mut b = vb;
1025            let mut c = vc;
1026            a.normal = face_normal;
1027            b.normal = face_normal;
1028            c.normal = face_normal;
1029
1030            new_verts.push(a);
1031            new_verts.push(b);
1032            new_verts.push(c);
1033            new_indices.push(TriangleIndex::new(base, base + 1, base + 2));
1034        }
1035
1036        self.vertices = new_verts;
1037        self.indices = new_indices;
1038    }
1039
1040    /// Subdivide each triangle into 4 triangles (Loop-style without smoothing).
1041    pub fn subdivide(&mut self) {
1042        let mut edge_midpoints: std::collections::HashMap<(u32, u32), u32> =
1043            std::collections::HashMap::new();
1044
1045        let mut new_indices = Vec::with_capacity(self.indices.len() * 4);
1046
1047        for tri in &self.indices {
1048            let mid_ab = Self::get_or_create_midpoint(
1049                &mut self.vertices, &mut edge_midpoints, tri.a, tri.b,
1050            );
1051            let mid_bc = Self::get_or_create_midpoint(
1052                &mut self.vertices, &mut edge_midpoints, tri.b, tri.c,
1053            );
1054            let mid_ca = Self::get_or_create_midpoint(
1055                &mut self.vertices, &mut edge_midpoints, tri.c, tri.a,
1056            );
1057
1058            new_indices.push(TriangleIndex::new(tri.a, mid_ab, mid_ca));
1059            new_indices.push(TriangleIndex::new(mid_ab, tri.b, mid_bc));
1060            new_indices.push(TriangleIndex::new(mid_ca, mid_bc, tri.c));
1061            new_indices.push(TriangleIndex::new(mid_ab, mid_bc, mid_ca));
1062        }
1063
1064        self.indices = new_indices;
1065        self.recompute_smooth_normals();
1066    }
1067
1068    fn get_or_create_midpoint(
1069        vertices: &mut Vec<SurfaceVertex>,
1070        edge_map: &mut std::collections::HashMap<(u32, u32), u32>,
1071        a: u32,
1072        b: u32,
1073    ) -> u32 {
1074        let key = if a < b { (a, b) } else { (b, a) };
1075        if let Some(&idx) = edge_map.get(&key) {
1076            return idx;
1077        }
1078        let va = vertices[a as usize];
1079        let vb = vertices[b as usize];
1080        let mid = SurfaceVertex {
1081            position: (va.position + vb.position) * 0.5,
1082            normal: (va.normal + vb.normal).normalize_or_zero(),
1083            tangent: (va.tangent + vb.tangent).normalize_or_zero(),
1084            bitangent: (va.bitangent + vb.bitangent).normalize_or_zero(),
1085            uv: [(va.uv[0] + vb.uv[0]) * 0.5, (va.uv[1] + vb.uv[1]) * 0.5],
1086            param_u: (va.param_u + vb.param_u) * 0.5,
1087            param_v: (va.param_v + vb.param_v) * 0.5,
1088        };
1089        let idx = vertices.len() as u32;
1090        vertices.push(mid);
1091        edge_map.insert(key, idx);
1092        idx
1093    }
1094
1095    /// Extract positions as a flat array of f32 triples (for GPU upload).
1096    pub fn positions_flat(&self) -> Vec<f32> {
1097        let mut out = Vec::with_capacity(self.vertices.len() * 3);
1098        for v in &self.vertices {
1099            out.push(v.position.x);
1100            out.push(v.position.y);
1101            out.push(v.position.z);
1102        }
1103        out
1104    }
1105
1106    /// Extract normals as a flat array of f32 triples.
1107    pub fn normals_flat(&self) -> Vec<f32> {
1108        let mut out = Vec::with_capacity(self.vertices.len() * 3);
1109        for v in &self.vertices {
1110            out.push(v.normal.x);
1111            out.push(v.normal.y);
1112            out.push(v.normal.z);
1113        }
1114        out
1115    }
1116
1117    /// Extract UVs as a flat array of f32 pairs.
1118    pub fn uvs_flat(&self) -> Vec<f32> {
1119        let mut out = Vec::with_capacity(self.vertices.len() * 2);
1120        for v in &self.vertices {
1121            out.push(v.uv[0]);
1122            out.push(v.uv[1]);
1123        }
1124        out
1125    }
1126
1127    /// Extract indices as a flat array of u32.
1128    pub fn indices_flat(&self) -> Vec<u32> {
1129        let mut out = Vec::with_capacity(self.indices.len() * 3);
1130        for tri in &self.indices {
1131            out.push(tri.a);
1132            out.push(tri.b);
1133            out.push(tri.c);
1134        }
1135        out
1136    }
1137
1138    /// Generate wireframe edges (deduplicated edge list).
1139    pub fn wireframe_edges(&self) -> Vec<(u32, u32)> {
1140        let mut edges = std::collections::HashSet::new();
1141        for tri in &self.indices {
1142            let mut add_edge = |a: u32, b: u32| {
1143                let key = if a < b { (a, b) } else { (b, a) };
1144                edges.insert(key);
1145            };
1146            add_edge(tri.a, tri.b);
1147            add_edge(tri.b, tri.c);
1148            add_edge(tri.c, tri.a);
1149        }
1150        edges.into_iter().collect()
1151    }
1152
1153    /// Sample the mesh at a given (u, v) by bilinear interpolation of grid vertices.
1154    pub fn sample(&self, u: f32, v: f32) -> Vec3 {
1155        let fu = u.clamp(0.0, 1.0) * self.resolution_u as f32;
1156        let fv = v.clamp(0.0, 1.0) * self.resolution_v as f32;
1157        let iu = (fu as usize).min(self.resolution_u - 1);
1158        let iv = (fv as usize).min(self.resolution_v - 1);
1159        let su = fu - iu as f32;
1160        let sv = fv - iv as f32;
1161
1162        let stride = self.resolution_u + 1;
1163        let i00 = iv * stride + iu;
1164        let i10 = i00 + 1;
1165        let i01 = i00 + stride;
1166        let i11 = i01 + 1;
1167
1168        let p00 = self.vertices[i00].position;
1169        let p10 = self.vertices[i10.min(self.vertices.len() - 1)].position;
1170        let p01 = self.vertices[i01.min(self.vertices.len() - 1)].position;
1171        let p11 = self.vertices[i11.min(self.vertices.len() - 1)].position;
1172
1173        let top = p00 * (1.0 - su) + p10 * su;
1174        let bottom = p01 * (1.0 - su) + p11 * su;
1175        top * (1.0 - sv) + bottom * sv
1176    }
1177
1178    /// Convert the mesh into a grid of characters for text-mode rendering.
1179    /// Maps the surface onto a width x height character grid using depth buffering.
1180    pub fn to_glyph_grid(&self, width: usize, height: usize) -> Vec<Vec<char>> {
1181        let mut grid = vec![vec![' '; width]; height];
1182        let mut depth = vec![vec![f32::MAX; width]; height];
1183
1184        let shading_chars = ['.', ':', '-', '=', '+', '*', '#', '%', '@'];
1185
1186        let center = self.center();
1187        let ext = self.extent();
1188        let max_ext = ext.x.max(ext.y).max(ext.z).max(0.001);
1189
1190        for tri in &self.indices {
1191            let va = &self.vertices[tri.a as usize];
1192            let vb = &self.vertices[tri.b as usize];
1193            let vc = &self.vertices[tri.c as usize];
1194
1195            // Project each vertex to screen space
1196            for vert in [va, vb, vc] {
1197                let rel = vert.position - center;
1198                let sx = ((rel.x / max_ext + 1.0) * 0.5 * (width - 1) as f32) as i32;
1199                let sy = ((1.0 - (rel.y / max_ext + 1.0) * 0.5) * (height - 1) as f32) as i32;
1200                let sz = rel.z / max_ext;
1201
1202                if sx >= 0 && sx < width as i32 && sy >= 0 && sy < height as i32 {
1203                    let ux = sx as usize;
1204                    let uy = sy as usize;
1205                    if sz < depth[uy][ux] {
1206                        depth[uy][ux] = sz;
1207                        // Use normal dot product with light direction for shading
1208                        let light = Vec3::new(0.577, 0.577, 0.577);
1209                        let intensity = vert.normal.dot(light).abs();
1210                        let idx = (intensity * (shading_chars.len() - 1) as f32) as usize;
1211                        let idx = idx.min(shading_chars.len() - 1);
1212                        grid[uy][ux] = shading_chars[idx];
1213                    }
1214                }
1215            }
1216        }
1217
1218        grid
1219    }
1220}
1221
1222// ─────────────────────────────────────────────────────────────────────────────
1223// Surface builder utilities
1224// ─────────────────────────────────────────────────────────────────────────────
1225
1226/// Create a surface of revolution from a 2D profile curve.
1227///
1228/// The profile is a function `f(v) -> (radius, height)` where v is in [0, 1].
1229/// The surface revolves around the Y axis.
1230pub struct RevolutionSurface {
1231    pub profile: Box<dyn Fn(f32) -> (f32, f32) + Send + Sync>,
1232    label: String,
1233}
1234
1235impl RevolutionSurface {
1236    pub fn new<F>(profile: F) -> Self
1237    where
1238        F: Fn(f32) -> (f32, f32) + Send + Sync + 'static,
1239    {
1240        Self {
1241            profile: Box::new(profile),
1242            label: "RevolutionSurface".to_string(),
1243        }
1244    }
1245
1246    pub fn with_name(mut self, name: &str) -> Self {
1247        self.label = name.to_string();
1248        self
1249    }
1250}
1251
1252impl Surface for RevolutionSurface {
1253    fn evaluate(&self, u: f32, v: f32) -> Vec3 {
1254        let theta = u * TAU;
1255        let (r, h) = (self.profile)(v);
1256        Vec3::new(r * theta.cos(), h, r * theta.sin())
1257    }
1258
1259    fn wraps_u(&self) -> bool { true }
1260    fn name(&self) -> &str { &self.label }
1261}
1262
1263/// Create a ruled surface between two space curves.
1264///
1265/// Given curves `a(t)` and `b(t)` for t in [0,1], the ruled surface is
1266/// `S(u, v) = (1-v) * a(u) + v * b(u)`.
1267pub struct RuledSurface {
1268    pub curve_a: Box<dyn Fn(f32) -> Vec3 + Send + Sync>,
1269    pub curve_b: Box<dyn Fn(f32) -> Vec3 + Send + Sync>,
1270    label: String,
1271}
1272
1273impl RuledSurface {
1274    pub fn new<A, B>(curve_a: A, curve_b: B) -> Self
1275    where
1276        A: Fn(f32) -> Vec3 + Send + Sync + 'static,
1277        B: Fn(f32) -> Vec3 + Send + Sync + 'static,
1278    {
1279        Self {
1280            curve_a: Box::new(curve_a),
1281            curve_b: Box::new(curve_b),
1282            label: "RuledSurface".to_string(),
1283        }
1284    }
1285
1286    pub fn with_name(mut self, name: &str) -> Self {
1287        self.label = name.to_string();
1288        self
1289    }
1290}
1291
1292impl Surface for RuledSurface {
1293    fn evaluate(&self, u: f32, v: f32) -> Vec3 {
1294        let a = (self.curve_a)(u);
1295        let b = (self.curve_b)(u);
1296        a * (1.0 - v) + b * v
1297    }
1298
1299    fn name(&self) -> &str { &self.label }
1300}
1301
1302/// Compose two surfaces by blending them with a weight function.
1303pub struct BlendedSurface<A: Surface, B: Surface> {
1304    pub surface_a: A,
1305    pub surface_b: B,
1306    pub blend_fn: Box<dyn Fn(f32, f32) -> f32 + Send + Sync>,
1307}
1308
1309impl<A: Surface, B: Surface> BlendedSurface<A, B> {
1310    pub fn new<F>(a: A, b: B, blend: F) -> Self
1311    where
1312        F: Fn(f32, f32) -> f32 + Send + Sync + 'static,
1313    {
1314        Self {
1315            surface_a: a,
1316            surface_b: b,
1317            blend_fn: Box::new(blend),
1318        }
1319    }
1320
1321    /// Uniform blend (constant t).
1322    pub fn uniform(a: A, b: B, t: f32) -> Self {
1323        Self::new(a, b, move |_u, _v| t)
1324    }
1325}
1326
1327impl<A: Surface, B: Surface> Surface for BlendedSurface<A, B> {
1328    fn evaluate(&self, u: f32, v: f32) -> Vec3 {
1329        let t = (self.blend_fn)(u, v).clamp(0.0, 1.0);
1330        let pa = self.surface_a.evaluate(u, v);
1331        let pb = self.surface_b.evaluate(u, v);
1332        pa * (1.0 - t) + pb * t
1333    }
1334
1335    fn name(&self) -> &str { "BlendedSurface" }
1336}
1337
1338/// Displace a surface along its normals by a scalar field.
1339pub struct DisplacedSurface<S: Surface> {
1340    pub base: S,
1341    pub displacement: Box<dyn Fn(f32, f32) -> f32 + Send + Sync>,
1342}
1343
1344impl<S: Surface> DisplacedSurface<S> {
1345    pub fn new<F>(base: S, displacement: F) -> Self
1346    where
1347        F: Fn(f32, f32) -> f32 + Send + Sync + 'static,
1348    {
1349        Self {
1350            base,
1351            displacement: Box::new(displacement),
1352        }
1353    }
1354}
1355
1356impl<S: Surface> Surface for DisplacedSurface<S> {
1357    fn evaluate(&self, u: f32, v: f32) -> Vec3 {
1358        let pos = self.base.evaluate(u, v);
1359        let n = self.base.normal(u, v);
1360        let d = (self.displacement)(u, v);
1361        pos + n * d
1362    }
1363
1364    fn name(&self) -> &str { "DisplacedSurface" }
1365}
1366
1367// ─────────────────────────────────────────────────────────────────────────────
1368// Surface catalog / registry
1369// ─────────────────────────────────────────────────────────────────────────────
1370
1371/// Enumeration of built-in surface types for catalog purposes.
1372#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
1373pub enum SurfaceKind {
1374    Sphere,
1375    Torus,
1376    MobiusStrip,
1377    KleinBottle,
1378    BoySurface,
1379    RomanSurface,
1380    CrossCap,
1381    TrefoilKnot,
1382    FigureEight,
1383    Catenoid,
1384    Helicoid,
1385    Enneper,
1386    Dini,
1387}
1388
1389impl SurfaceKind {
1390    /// All built-in surface types.
1391    pub fn all() -> &'static [SurfaceKind] {
1392        &[
1393            SurfaceKind::Sphere,
1394            SurfaceKind::Torus,
1395            SurfaceKind::MobiusStrip,
1396            SurfaceKind::KleinBottle,
1397            SurfaceKind::BoySurface,
1398            SurfaceKind::RomanSurface,
1399            SurfaceKind::CrossCap,
1400            SurfaceKind::TrefoilKnot,
1401            SurfaceKind::FigureEight,
1402            SurfaceKind::Catenoid,
1403            SurfaceKind::Helicoid,
1404            SurfaceKind::Enneper,
1405            SurfaceKind::Dini,
1406        ]
1407    }
1408
1409    /// Human-readable name.
1410    pub fn name(self) -> &'static str {
1411        match self {
1412            SurfaceKind::Sphere => "Sphere",
1413            SurfaceKind::Torus => "Torus",
1414            SurfaceKind::MobiusStrip => "Mobius Strip",
1415            SurfaceKind::KleinBottle => "Klein Bottle",
1416            SurfaceKind::BoySurface => "Boy's Surface",
1417            SurfaceKind::RomanSurface => "Roman Surface",
1418            SurfaceKind::CrossCap => "Cross-Cap",
1419            SurfaceKind::TrefoilKnot => "Trefoil Knot",
1420            SurfaceKind::FigureEight => "Figure-Eight",
1421            SurfaceKind::Catenoid => "Catenoid",
1422            SurfaceKind::Helicoid => "Helicoid",
1423            SurfaceKind::Enneper => "Enneper Surface",
1424            SurfaceKind::Dini => "Dini's Surface",
1425        }
1426    }
1427
1428    /// Create a default-parameterized mesh for this surface at given resolution.
1429    pub fn tessellate(self, res_u: usize, res_v: usize) -> SurfaceMesh {
1430        match self {
1431            SurfaceKind::Sphere => SurfaceMesh::tessellate(&Sphere::default(), res_u, res_v),
1432            SurfaceKind::Torus => SurfaceMesh::tessellate(&Torus::default(), res_u, res_v),
1433            SurfaceKind::MobiusStrip => SurfaceMesh::tessellate(&MobiusStrip::default(), res_u, res_v),
1434            SurfaceKind::KleinBottle => SurfaceMesh::tessellate(&KleinBottle::default(), res_u, res_v),
1435            SurfaceKind::BoySurface => SurfaceMesh::tessellate(&BoySurface::default(), res_u, res_v),
1436            SurfaceKind::RomanSurface => SurfaceMesh::tessellate(&RomanSurface::default(), res_u, res_v),
1437            SurfaceKind::CrossCap => SurfaceMesh::tessellate(&CrossCap::default(), res_u, res_v),
1438            SurfaceKind::TrefoilKnot => SurfaceMesh::tessellate(&TrefoilKnot::default(), res_u, res_v),
1439            SurfaceKind::FigureEight => SurfaceMesh::tessellate(&FigureEight::default(), res_u, res_v),
1440            SurfaceKind::Catenoid => SurfaceMesh::tessellate(&Catenoid::default(), res_u, res_v),
1441            SurfaceKind::Helicoid => SurfaceMesh::tessellate(&Helicoid::default(), res_u, res_v),
1442            SurfaceKind::Enneper => SurfaceMesh::tessellate(&EnneperSurface::default(), res_u, res_v),
1443            SurfaceKind::Dini => SurfaceMesh::tessellate(&DiniSurface::default(), res_u, res_v),
1444        }
1445    }
1446}
1447
1448// ─────────────────────────────────────────────────────────────────────────────
1449// Iso-curve extraction
1450// ─────────────────────────────────────────────────────────────────────────────
1451
1452/// Extract iso-parameter curves from a surface.
1453pub struct IsoCurves;
1454
1455impl IsoCurves {
1456    /// Extract a curve at constant u, varying v from 0 to 1.
1457    pub fn at_constant_u(surface: &dyn Surface, u: f32, samples: usize) -> Vec<Vec3> {
1458        let n = samples.max(2);
1459        (0..n).map(|i| {
1460            let v = i as f32 / (n - 1) as f32;
1461            surface.evaluate(u, v)
1462        }).collect()
1463    }
1464
1465    /// Extract a curve at constant v, varying u from 0 to 1.
1466    pub fn at_constant_v(surface: &dyn Surface, v: f32, samples: usize) -> Vec<Vec3> {
1467        let n = samples.max(2);
1468        (0..n).map(|i| {
1469            let u = i as f32 / (n - 1) as f32;
1470            surface.evaluate(u, v)
1471        }).collect()
1472    }
1473
1474    /// Extract a grid of iso-curves in both directions.
1475    pub fn grid(
1476        surface: &dyn Surface,
1477        u_curves: usize,
1478        v_curves: usize,
1479        samples: usize,
1480    ) -> (Vec<Vec<Vec3>>, Vec<Vec<Vec3>>) {
1481        let u_lines: Vec<Vec<Vec3>> = (0..u_curves)
1482            .map(|i| {
1483                let u = i as f32 / (u_curves - 1).max(1) as f32;
1484                Self::at_constant_u(surface, u, samples)
1485            })
1486            .collect();
1487
1488        let v_lines: Vec<Vec<Vec3>> = (0..v_curves)
1489            .map(|i| {
1490                let v = i as f32 / (v_curves - 1).max(1) as f32;
1491                Self::at_constant_v(surface, v, samples)
1492            })
1493            .collect();
1494
1495        (u_lines, v_lines)
1496    }
1497}
1498
1499// ─────────────────────────────────────────────────────────────────────────────
1500// Ray-surface intersection
1501// ─────────────────────────────────────────────────────────────────────────────
1502
1503/// Result of a ray-mesh intersection test.
1504#[derive(Debug, Clone, Copy)]
1505pub struct RayHit {
1506    pub distance: f32,
1507    pub position: Vec3,
1508    pub normal: Vec3,
1509    pub uv: [f32; 2],
1510    pub triangle_index: usize,
1511}
1512
1513impl SurfaceMesh {
1514    /// Ray-mesh intersection using Moller-Trumbore algorithm.
1515    /// Returns the closest hit, if any.
1516    pub fn ray_intersect(&self, ray_origin: Vec3, ray_dir: Vec3) -> Option<RayHit> {
1517        let mut closest: Option<RayHit> = None;
1518        let eps = 1e-7;
1519
1520        for (tri_idx, tri) in self.indices.iter().enumerate() {
1521            let v0 = self.vertices[tri.a as usize].position;
1522            let v1 = self.vertices[tri.b as usize].position;
1523            let v2 = self.vertices[tri.c as usize].position;
1524
1525            let e1 = v1 - v0;
1526            let e2 = v2 - v0;
1527            let h = ray_dir.cross(e2);
1528            let a = e1.dot(h);
1529
1530            if a.abs() < eps {
1531                continue;
1532            }
1533
1534            let f = 1.0 / a;
1535            let s = ray_origin - v0;
1536            let u_bary = f * s.dot(h);
1537            if !(0.0..=1.0).contains(&u_bary) {
1538                continue;
1539            }
1540
1541            let q = s.cross(e1);
1542            let v_bary = f * ray_dir.dot(q);
1543            if v_bary < 0.0 || u_bary + v_bary > 1.0 {
1544                continue;
1545            }
1546
1547            let t = f * e2.dot(q);
1548            if t < eps {
1549                continue;
1550            }
1551
1552            if closest.as_ref().map_or(true, |c| t < c.distance) {
1553                let w = 1.0 - u_bary - v_bary;
1554                let va = &self.vertices[tri.a as usize];
1555                let vb = &self.vertices[tri.b as usize];
1556                let vc = &self.vertices[tri.c as usize];
1557
1558                let normal = (va.normal * w + vb.normal * u_bary + vc.normal * v_bary)
1559                    .normalize_or_zero();
1560                let uv = [
1561                    va.uv[0] * w + vb.uv[0] * u_bary + vc.uv[0] * v_bary,
1562                    va.uv[1] * w + vb.uv[1] * u_bary + vc.uv[1] * v_bary,
1563                ];
1564
1565                closest = Some(RayHit {
1566                    distance: t,
1567                    position: ray_origin + ray_dir * t,
1568                    normal,
1569                    uv,
1570                    triangle_index: tri_idx,
1571                });
1572            }
1573        }
1574
1575        closest
1576    }
1577}
1578
1579// ─────────────────────────────────────────────────────────────────────────────
1580// Tests
1581// ─────────────────────────────────────────────────────────────────────────────
1582
1583#[cfg(test)]
1584mod tests {
1585    use super::*;
1586
1587    #[test]
1588    fn sphere_poles() {
1589        let s = Sphere::default();
1590        let north = s.evaluate(0.0, 0.0);
1591        let south = s.evaluate(0.0, 1.0);
1592        assert!((north.y - 1.0).abs() < 1e-5);
1593        assert!((south.y + 1.0).abs() < 1e-5);
1594    }
1595
1596    #[test]
1597    fn torus_evaluate() {
1598        let t = Torus::default();
1599        let p = t.evaluate(0.0, 0.0);
1600        assert!((p.x - 1.3).abs() < 0.01); // major + minor
1601    }
1602
1603    #[test]
1604    fn tessellate_sphere() {
1605        let s = Sphere::default();
1606        let mesh = SurfaceMesh::tessellate(&s, 16, 8);
1607        assert_eq!(mesh.vertex_count(), 17 * 9);
1608        assert_eq!(mesh.triangle_count(), 16 * 8 * 2);
1609    }
1610
1611    #[test]
1612    fn surface_mesh_area() {
1613        let s = Sphere::new(1.0);
1614        let mesh = SurfaceMesh::tessellate(&s, 64, 32);
1615        let area = mesh.surface_area();
1616        // 4 * pi * r^2 ~ 12.566
1617        assert!((area - 4.0 * PI).abs() < 0.5);
1618    }
1619
1620    #[test]
1621    fn ray_intersect_sphere() {
1622        let s = Sphere::new(1.0);
1623        let mesh = SurfaceMesh::tessellate(&s, 32, 16);
1624        let hit = mesh.ray_intersect(Vec3::new(0.0, 0.0, -5.0), Vec3::Z);
1625        assert!(hit.is_some());
1626        let h = hit.unwrap();
1627        assert!((h.distance - 4.0).abs() < 0.2);
1628    }
1629
1630    #[test]
1631    fn surface_kind_catalog() {
1632        assert_eq!(SurfaceKind::all().len(), 13);
1633        for kind in SurfaceKind::all() {
1634            let mesh = kind.tessellate(8, 8);
1635            assert!(mesh.vertex_count() > 0);
1636            assert!(mesh.triangle_count() > 0);
1637        }
1638    }
1639
1640    #[test]
1641    fn function_surface() {
1642        let fs = FunctionSurface::new(|u, v| Vec3::new(u, v, 0.0));
1643        let p = fs.evaluate(0.5, 0.5);
1644        assert!((p.x - 0.5).abs() < 1e-5);
1645    }
1646
1647    #[test]
1648    fn wireframe_edges() {
1649        let s = Sphere::default();
1650        let mesh = SurfaceMesh::tessellate(&s, 4, 4);
1651        let edges = mesh.wireframe_edges();
1652        assert!(!edges.is_empty());
1653    }
1654
1655    #[test]
1656    fn subdivide_mesh() {
1657        let s = Sphere::default();
1658        let mut mesh = SurfaceMesh::tessellate(&s, 4, 4);
1659        let original_tris = mesh.triangle_count();
1660        mesh.subdivide();
1661        assert_eq!(mesh.triangle_count(), original_tris * 4);
1662    }
1663
1664    #[test]
1665    fn flat_shaded_mesh() {
1666        let s = Torus::default();
1667        let mut mesh = SurfaceMesh::tessellate(&s, 8, 8);
1668        let tri_count = mesh.triangle_count();
1669        mesh.make_flat_shaded();
1670        assert_eq!(mesh.vertex_count(), tri_count * 3);
1671    }
1672
1673    #[test]
1674    fn iso_curves_extraction() {
1675        let s = Sphere::default();
1676        let curve = IsoCurves::at_constant_u(&s, 0.0, 32);
1677        assert_eq!(curve.len(), 32);
1678    }
1679}
proof_engine/surfaces/parametric.rs

proof_engine/surfaces/
parametric.rs
