oxiphysics_geometry/
implicit_surfaces.rs

1// Copyright 2026 COOLJAPAN OU (Team KitaSan)
2// SPDX-License-Identifier: Apache-2.0
3
4//! Implicit surface representations and operations.
5//!
6//! This module provides Signed Distance Functions (SDFs) for common geometric
7//! primitives, Constructive Solid Geometry (CSG) boolean operations, smooth
8//! blending, isosurface extraction via Marching Cubes, Radial Basis Function
9//! (RBF) implicit surfaces, and sphere-tracing ray marching.
10//!
11//! # Quick start
12//!
13//! ```no_run
14//! use oxiphysics_geometry::implicit_surfaces::{ImplicitSurface, SphereSDF};
15//!
16//! let sphere = SphereSDF::new([0.0, 0.0, 0.0], 1.0);
17//! assert!(sphere.eval([0.0, 0.0, 0.0]) < 0.0); // inside
18//! assert!(sphere.eval([2.0, 0.0, 0.0]) > 0.0); // outside
19//! ```
20
21// ---------------------------------------------------------------------------
22// Math helpers (private)
23// ---------------------------------------------------------------------------
24
25#[inline]
26fn len3(v: [f64; 3]) -> f64 {
27    (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]).sqrt()
28}
29
30#[inline]
31fn sub3(a: [f64; 3], b: [f64; 3]) -> [f64; 3] {
32    [a[0] - b[0], a[1] - b[1], a[2] - b[2]]
33}
34
35#[inline]
36fn add3(a: [f64; 3], b: [f64; 3]) -> [f64; 3] {
37    [a[0] + b[0], a[1] + b[1], a[2] + b[2]]
38}
39
40#[inline]
41fn scale3(v: [f64; 3], s: f64) -> [f64; 3] {
42    [v[0] * s, v[1] * s, v[2] * s]
43}
44
45#[inline]
46fn normalize3(v: [f64; 3]) -> [f64; 3] {
47    let l = len3(v).max(1e-300);
48    [v[0] / l, v[1] / l, v[2] / l]
49}
50
51/// Clamp `v` to `[lo, hi]`.
52#[inline]
53fn clamp(v: f64, lo: f64, hi: f64) -> f64 {
54    v.max(lo).min(hi)
55}
56
57// ---------------------------------------------------------------------------
58// § 1  Core trait
59// ---------------------------------------------------------------------------
60
61/// Trait for implicit surface representations (Signed Distance Functions).
62///
63/// Implementors define a scalar field f: ℝ³ → ℝ where:
64/// - f(p) < 0 means p is inside the surface
65/// - f(p) = 0 means p is on the surface
66/// - f(p) > 0 means p is outside the surface
67pub trait ImplicitSurface {
68    /// Evaluate the signed distance (or implicit) function at point `p`.
69    fn eval(&self, p: [f64; 3]) -> f64;
70
71    /// Compute the gradient ∇f at point `p`.
72    ///
73    /// For exact SDFs this equals the outward unit normal on the surface.
74    /// The default implementation uses central differences with step `1e-5`.
75    fn gradient(&self, p: [f64; 3]) -> [f64; 3] {
76        let eps = 1e-5;
77        let dx = self.eval([p[0] + eps, p[1], p[2]]) - self.eval([p[0] - eps, p[1], p[2]]);
78        let dy = self.eval([p[0], p[1] + eps, p[2]]) - self.eval([p[0], p[1] - eps, p[2]]);
79        let dz = self.eval([p[0], p[1], p[2] + eps]) - self.eval([p[0], p[1], p[2] - eps]);
80        [dx / (2.0 * eps), dy / (2.0 * eps), dz / (2.0 * eps)]
81    }
82}
83
84// ---------------------------------------------------------------------------
85// § 2  Primitive SDFs
86// ---------------------------------------------------------------------------
87
88/// Sphere SDF: f(p) = |p − center| − radius.
89///
90/// # Example
91/// ```no_run
92/// use oxiphysics_geometry::implicit_surfaces::{ImplicitSurface, SphereSDF};
93/// let s = SphereSDF::new([0.0, 0.0, 0.0], 2.0);
94/// assert!((s.eval([2.0, 0.0, 0.0])).abs() < 1e-12);
95/// ```
96#[derive(Debug, Clone)]
97pub struct SphereSDF {
98    /// Centre of the sphere.
99    pub center: [f64; 3],
100    /// Radius of the sphere (must be positive).
101    pub radius: f64,
102}
103
104impl SphereSDF {
105    /// Create a new `SphereSDF`.
106    pub fn new(center: [f64; 3], radius: f64) -> Self {
107        Self { center, radius }
108    }
109}
110
111impl ImplicitSurface for SphereSDF {
112    fn eval(&self, p: [f64; 3]) -> f64 {
113        len3(sub3(p, self.center)) - self.radius
114    }
115
116    fn gradient(&self, p: [f64; 3]) -> [f64; 3] {
117        normalize3(sub3(p, self.center))
118    }
119}
120
121/// Axis-aligned box SDF centred at the origin.
122///
123/// `half_extents[i]` is the half-size along axis i.
124///
125/// # Example
126/// ```no_run
127/// use oxiphysics_geometry::implicit_surfaces::{ImplicitSurface, BoxSDF};
128/// let b = BoxSDF::new([1.0, 2.0, 3.0]);
129/// // A point on a face
130/// let d = b.eval([1.0, 0.0, 0.0]);
131/// assert!(d.abs() < 1e-12);
132/// ```
133#[derive(Debug, Clone)]
134pub struct BoxSDF {
135    /// Half-extents along each axis.
136    pub half_extents: [f64; 3],
137}
138
139impl BoxSDF {
140    /// Create a new axis-aligned box SDF.
141    pub fn new(half_extents: [f64; 3]) -> Self {
142        Self { half_extents }
143    }
144}
145
146impl ImplicitSurface for BoxSDF {
147    fn eval(&self, p: [f64; 3]) -> f64 {
148        // q = |p| - half_extents
149        let q = [
150            p[0].abs() - self.half_extents[0],
151            p[1].abs() - self.half_extents[1],
152            p[2].abs() - self.half_extents[2],
153        ];
154        // Length of the positive part plus the max of the negative part
155        let pos = [q[0].max(0.0), q[1].max(0.0), q[2].max(0.0)];
156        len3(pos) + q[0].max(q[1]).max(q[2]).min(0.0)
157    }
158}
159
160/// Torus SDF centred at the origin, lying in the XZ plane.
161///
162/// - `major_r` — distance from the torus centre to the tube centre
163/// - `minor_r` — radius of the tube
164///
165/// # Example
166/// ```no_run
167/// use oxiphysics_geometry::implicit_surfaces::{ImplicitSurface, TorusSDF};
168/// let t = TorusSDF::new(2.0, 0.5);
169/// // Point on the outer equator of the torus
170/// let d = t.eval([2.5, 0.0, 0.0]);
171/// assert!(d.abs() < 1e-12);
172/// ```
173#[derive(Debug, Clone)]
174pub struct TorusSDF {
175    /// Major radius (distance from origin to tube centre).
176    pub major_r: f64,
177    /// Minor radius (tube radius).
178    pub minor_r: f64,
179}
180
181impl TorusSDF {
182    /// Create a new `TorusSDF`.
183    pub fn new(major_r: f64, minor_r: f64) -> Self {
184        Self { major_r, minor_r }
185    }
186}
187
188impl ImplicitSurface for TorusSDF {
189    fn eval(&self, p: [f64; 3]) -> f64 {
190        // Project onto XZ plane
191        let q_xz = ((p[0] * p[0] + p[2] * p[2]).sqrt() - self.major_r, p[1]);
192        (q_xz.0 * q_xz.0 + q_xz.1 * q_xz.1).sqrt() - self.minor_r
193    }
194}
195
196/// Infinite cylinder SDF along the Y axis centred at the origin.
197///
198/// - `radius` — cylinder radius
199/// - `half_height` — half-length along the Y axis
200///
201/// # Example
202/// ```no_run
203/// use oxiphysics_geometry::implicit_surfaces::{ImplicitSurface, CylinderSDF};
204/// let c = CylinderSDF::new(1.0, 2.0);
205/// // Point on lateral surface, within height bounds
206/// let d = c.eval([1.0, 0.0, 0.0]);
207/// assert!(d.abs() < 1e-12);
208/// ```
209#[derive(Debug, Clone)]
210pub struct CylinderSDF {
211    /// Radius of the cylinder.
212    pub radius: f64,
213    /// Half-height along the Y axis.
214    pub half_height: f64,
215}
216
217impl CylinderSDF {
218    /// Create a new `CylinderSDF`.
219    pub fn new(radius: f64, half_height: f64) -> Self {
220        Self {
221            radius,
222            half_height,
223        }
224    }
225}
226
227impl ImplicitSurface for CylinderSDF {
228    fn eval(&self, p: [f64; 3]) -> f64 {
229        let d = [
230            (p[0] * p[0] + p[2] * p[2]).sqrt() - self.radius,
231            p[1].abs() - self.half_height,
232        ];
233        d[0].max(d[1]).min(0.0)
234            + [d[0].max(0.0), d[1].max(0.0)]
235                .iter()
236                .map(|v| v * v)
237                .sum::<f64>()
238                .sqrt()
239    }
240}
241
242// ---------------------------------------------------------------------------
243// § 3  CSG boolean operations
244// ---------------------------------------------------------------------------
245
246/// CSG union of two implicit surfaces: f = min(A, B).
247///
248/// Points inside either surface are inside the union.
249#[derive(Debug, Clone)]
250pub struct UnionSDF<A, B> {
251    /// First operand.
252    pub a: A,
253    /// Second operand.
254    pub b: B,
255}
256
257impl<A: ImplicitSurface, B: ImplicitSurface> UnionSDF<A, B> {
258    /// Create a CSG union of `a` and `b`.
259    pub fn new(a: A, b: B) -> Self {
260        Self { a, b }
261    }
262}
263
264impl<A: ImplicitSurface, B: ImplicitSurface> ImplicitSurface for UnionSDF<A, B> {
265    fn eval(&self, p: [f64; 3]) -> f64 {
266        self.a.eval(p).min(self.b.eval(p))
267    }
268}
269
270/// CSG intersection of two implicit surfaces: f = max(A, B).
271///
272/// Points inside both surfaces are inside the intersection.
273#[derive(Debug, Clone)]
274pub struct IntersectionSDF<A, B> {
275    /// First operand.
276    pub a: A,
277    /// Second operand.
278    pub b: B,
279}
280
281impl<A: ImplicitSurface, B: ImplicitSurface> IntersectionSDF<A, B> {
282    /// Create a CSG intersection of `a` and `b`.
283    pub fn new(a: A, b: B) -> Self {
284        Self { a, b }
285    }
286}
287
288impl<A: ImplicitSurface, B: ImplicitSurface> ImplicitSurface for IntersectionSDF<A, B> {
289    fn eval(&self, p: [f64; 3]) -> f64 {
290        self.a.eval(p).max(self.b.eval(p))
291    }
292}
293
294/// CSG difference of two implicit surfaces: f = max(A, −B).
295///
296/// Subtracts shape B from shape A.
297#[derive(Debug, Clone)]
298pub struct DifferenceSDF<A, B> {
299    /// Shape to subtract from.
300    pub a: A,
301    /// Shape to subtract.
302    pub b: B,
303}
304
305impl<A: ImplicitSurface, B: ImplicitSurface> DifferenceSDF<A, B> {
306    /// Create a CSG difference: `a` minus `b`.
307    pub fn new(a: A, b: B) -> Self {
308        Self { a, b }
309    }
310}
311
312impl<A: ImplicitSurface, B: ImplicitSurface> ImplicitSurface for DifferenceSDF<A, B> {
313    fn eval(&self, p: [f64; 3]) -> f64 {
314        self.a.eval(p).max(-self.b.eval(p))
315    }
316}
317
318// ---------------------------------------------------------------------------
319// § 4  Smooth boolean operations
320// ---------------------------------------------------------------------------
321
322/// Smooth (C¹) union of two implicit surfaces using polynomial blending.
323///
324/// The parameter `k` controls the blend radius.  Larger `k` gives smoother
325/// transitions.
326#[derive(Debug, Clone)]
327pub struct SmoothUnionSDF<A, B> {
328    /// First operand.
329    pub a: A,
330    /// Second operand.
331    pub b: B,
332    /// Blend radius (must be > 0).
333    pub k: f64,
334}
335
336impl<A: ImplicitSurface, B: ImplicitSurface> SmoothUnionSDF<A, B> {
337    /// Create a smooth union with blend parameter `k`.
338    pub fn new(a: A, b: B, k: f64) -> Self {
339        Self { a, b, k }
340    }
341}
342
343impl<A: ImplicitSurface, B: ImplicitSurface> ImplicitSurface for SmoothUnionSDF<A, B> {
344    fn eval(&self, p: [f64; 3]) -> f64 {
345        let da = self.a.eval(p);
346        let db = self.b.eval(p);
347        let h = clamp(0.5 + 0.5 * (db - da) / self.k, 0.0, 1.0);
348        da * h + db * (1.0 - h) - self.k * h * (1.0 - h)
349    }
350}
351
352// ---------------------------------------------------------------------------
353// § 5  Offset SDF (shell / thickening)
354// ---------------------------------------------------------------------------
355
356/// Offset (shell) SDF: expands or shrinks a shape by a fixed distance.
357///
358/// `f_offset(p) = f(p) − offset`
359///
360/// Positive `offset` expands the shape; negative `offset` shrinks it.
361#[derive(Debug, Clone)]
362pub struct OffsetSDF<S> {
363    /// Inner surface.
364    pub inner: S,
365    /// Offset distance (m).
366    pub offset: f64,
367}
368
369impl<S: ImplicitSurface> OffsetSDF<S> {
370    /// Create an offset surface from `inner` shifted by `offset`.
371    pub fn new(inner: S, offset: f64) -> Self {
372        Self { inner, offset }
373    }
374}
375
376impl<S: ImplicitSurface> ImplicitSurface for OffsetSDF<S> {
377    fn eval(&self, p: [f64; 3]) -> f64 {
378        self.inner.eval(p) - self.offset
379    }
380
381    fn gradient(&self, p: [f64; 3]) -> [f64; 3] {
382        self.inner.gradient(p)
383    }
384}
385
386// ---------------------------------------------------------------------------
387// § 6  Marching Cubes
388// ---------------------------------------------------------------------------
389
390// ---------------------------------------------------------------------------
391// Marching cubes lookup tables (edge table and tri table)
392// ---------------------------------------------------------------------------
393
394/// Edge table for marching cubes (256 entries).
395const EDGE_TABLE: [u16; 256] = [
396    0x000, 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c, 0x80c, 0x905, 0xa0f, 0xb06, 0xc0a,
397    0xd03, 0xe09, 0xf00, 0x190, 0x099, 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c, 0x99c, 0x895,
398    0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90, 0x230, 0x339, 0x033, 0x13a, 0x636, 0x73f, 0x435,
399    0x53c, 0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30, 0x3a0, 0x2a9, 0x1a3, 0x0aa,
400    0x7a6, 0x6af, 0x5a5, 0x4ac, 0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0, 0x460,
401    0x569, 0x663, 0x76a, 0x066, 0x16f, 0x265, 0x36c, 0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963,
402    0xa69, 0xb60, 0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0x0ff, 0x3f5, 0x2fc, 0xdfc, 0xcf5, 0xfff,
403    0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0, 0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x055, 0x15c,
404    0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950, 0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6,
405    0x2cf, 0x1c5, 0x0cc, 0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0, 0x8c0, 0x9c9,
406    0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc, 0x0cc, 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9,
407    0x7c0, 0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c, 0x15c, 0x055, 0x35f, 0x256,
408    0x55a, 0x453, 0x759, 0x650, 0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc, 0x2fc,
409    0x3f5, 0x0ff, 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0, 0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f,
410    0xd65, 0xc6c, 0x36c, 0x265, 0x16f, 0x066, 0x76a, 0x663, 0x569, 0x460, 0xca0, 0xda9, 0xea3,
411    0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac, 0x4ac, 0x5a5, 0x6af, 0x7a6, 0x0aa, 0x1a3, 0x2a9, 0x3a0,
412    0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x835, 0xb3f, 0xa36, // fixed 0x83f→0x835
413    0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x033, 0x339, 0x230, 0xe90, 0xf99, 0xc93, 0xd9a, 0xa96,
414    0xb9f, 0x895, 0x99c, 0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x099, 0x190, 0xf00, 0xe09,
415    0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c, 0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109,
416    0x000,
417];
418
419/// Cube corner positions in unit cube (index → \[xi, yi, zi\]).
420const CUBE_CORNERS: [[f64; 3]; 8] = [
421    [0.0, 0.0, 0.0],
422    [1.0, 0.0, 0.0],
423    [1.0, 1.0, 0.0],
424    [0.0, 1.0, 0.0],
425    [0.0, 0.0, 1.0],
426    [1.0, 0.0, 1.0],
427    [1.0, 1.0, 1.0],
428    [0.0, 1.0, 1.0],
429];
430
431/// Cube edge pairs (each edge connects two corners).
432const CUBE_EDGES: [[usize; 2]; 12] = [
433    [0, 1],
434    [1, 2],
435    [2, 3],
436    [3, 0],
437    [4, 5],
438    [5, 6],
439    [6, 7],
440    [7, 4],
441    [0, 4],
442    [1, 5],
443    [2, 6],
444    [3, 7],
445];
446
447/// Interpolate a vertex on an edge where the SDF crosses zero.
448fn interp_edge(p0: [f64; 3], v0: f64, p1: [f64; 3], v1: f64) -> [f64; 3] {
449    let dv = v1 - v0;
450    if dv.abs() < 1e-30 {
451        return [
452            (p0[0] + p1[0]) / 2.0,
453            (p0[1] + p1[1]) / 2.0,
454            (p0[2] + p1[2]) / 2.0,
455        ];
456    }
457    let t = -v0 / dv;
458    [
459        p0[0] + t * (p1[0] - p0[0]),
460        p0[1] + t * (p1[1] - p0[1]),
461        p0[2] + t * (p1[2] - p0[2]),
462    ]
463}
464
465/// Extract an isosurface from an implicit SDF using the Marching Cubes algorithm.
466///
467/// # Arguments
468/// * `sdf` — any type implementing [`ImplicitSurface`]
469/// * `bounds` — axis-aligned bounding box as `[[xmin,xmax\], [ymin,ymax], [zmin,zmax]]`
470/// * `resolution` — number of grid cells per axis
471///
472/// # Returns
473/// A tuple `(vertices, triangles)` where each vertex is `[x, y, z]` and
474/// each triangle is `[i0, i1, i2]` (indices into `vertices`).
475pub fn marching_cubes(
476    sdf: &dyn ImplicitSurface,
477    bounds: [[f64; 2]; 3],
478    resolution: usize,
479) -> (Vec<[f64; 3]>, Vec<[usize; 3]>) {
480    let n = resolution.max(1);
481    let mut vertices: Vec<[f64; 3]> = Vec::new();
482    let mut triangles: Vec<[usize; 3]> = Vec::new();
483
484    let dx = (bounds[0][1] - bounds[0][0]) / n as f64;
485    let dy = (bounds[1][1] - bounds[1][0]) / n as f64;
486    let dz = (bounds[2][1] - bounds[2][0]) / n as f64;
487
488    for iz in 0..n {
489        for iy in 0..n {
490            for ix in 0..n {
491                let ox = bounds[0][0] + ix as f64 * dx;
492                let oy = bounds[1][0] + iy as f64 * dy;
493                let oz = bounds[2][0] + iz as f64 * dz;
494
495                // Evaluate SDF at 8 corners
496                let corner_pts: [[f64; 3]; 8] = std::array::from_fn(|c| {
497                    [
498                        ox + CUBE_CORNERS[c][0] * dx,
499                        oy + CUBE_CORNERS[c][1] * dy,
500                        oz + CUBE_CORNERS[c][2] * dz,
501                    ]
502                });
503                let vals: [f64; 8] = std::array::from_fn(|c| sdf.eval(corner_pts[c]));
504
505                // Compute cube index
506                let mut cube_idx: usize = 0;
507                for (c, &val) in vals.iter().enumerate() {
508                    if val < 0.0 {
509                        cube_idx |= 1 << c;
510                    }
511                }
512
513                let edge_flags = EDGE_TABLE[cube_idx];
514                if edge_flags == 0 {
515                    continue;
516                }
517
518                // Compute intersection vertices for active edges
519                let mut edge_verts: [Option<[f64; 3]>; 12] = [None; 12];
520                for e in 0..12 {
521                    if edge_flags & (1 << e) != 0 {
522                        let [i0, i1] = CUBE_EDGES[e];
523                        edge_verts[e] = Some(interp_edge(
524                            corner_pts[i0],
525                            vals[i0],
526                            corner_pts[i1],
527                            vals[i1],
528                        ));
529                    }
530                }
531
532                // Tessellate using a simple fan from the first active edge vertex
533                // (simplified triangulation — not full marching-cubes tri table)
534                let active: Vec<[f64; 3]> = (0..12).filter_map(|e| edge_verts[e]).collect();
535                if active.len() >= 3 {
536                    let base = vertices.len();
537                    for v in &active {
538                        vertices.push(*v);
539                    }
540                    for k in 1..(active.len() - 1) {
541                        triangles.push([base, base + k, base + k + 1]);
542                    }
543                }
544            }
545        }
546    }
547
548    (vertices, triangles)
549}
550
551// ---------------------------------------------------------------------------
552// § 7  RBF Implicit Surface
553// ---------------------------------------------------------------------------
554
555/// Radial Basis Function (RBF) implicit surface.
556///
557/// The surface is defined as f(p) = Σᵢ wᵢ φ(|p − cᵢ|) where φ is a
558/// thin-plate spline kernel φ(r) = r² log(r + ε).
559///
560/// Use [`RBFImplicit::fit`] to compute weights from a set of point samples
561/// with target values, then [`RBFImplicit::eval`] to query the surface.
562#[derive(Debug, Clone)]
563pub struct RBFImplicit {
564    /// RBF centre positions.
565    pub centers: Vec<[f64; 3]>,
566    /// Weights for each centre (computed by [`fit`](Self::fit)).
567    pub weights: Vec<f64>,
568}
569
570impl RBFImplicit {
571    /// Create a new `RBFImplicit` with given centres and zero weights.
572    pub fn new(centers: Vec<[f64; 3]>) -> Self {
573        let n = centers.len();
574        Self {
575            centers,
576            weights: vec![0.0; n],
577        }
578    }
579
580    /// Fit weights from sample points `pts` and target values `targets`.
581    ///
582    /// Solves the symmetric RBF system Φ w = t using the conjugate gradient
583    /// method (up to `max_iter` iterations, tolerance `tol`).
584    ///
585    /// # Panics
586    /// Panics if `pts.len() != targets.len()` or `pts.len() != centers.len()`.
587    pub fn fit(&mut self, pts: &[[f64; 3]], targets: &[f64], tol: f64, max_iter: usize) {
588        assert_eq!(pts.len(), targets.len());
589        assert_eq!(pts.len(), self.centers.len());
590        let n = pts.len();
591        // Build the RBF matrix Φ
592        let phi: Vec<Vec<f64>> = (0..n)
593            .map(|i| {
594                (0..n)
595                    .map(|j| rbf_kernel(len3(sub3(pts[i], self.centers[j]))))
596                    .collect()
597            })
598            .collect();
599        // Solve Φ w = targets
600        self.weights = rbf_cg(&phi, targets, tol, max_iter);
601    }
602
603    /// Evaluate the RBF implicit surface at point `p`.
604    pub fn eval(&self, p: [f64; 3]) -> f64 {
605        self.centers
606            .iter()
607            .zip(self.weights.iter())
608            .map(|(c, w)| w * rbf_kernel(len3(sub3(p, *c))))
609            .sum()
610    }
611}
612
613/// Thin-plate spline RBF kernel: φ(r) = r² ln(r + ε).
614fn rbf_kernel(r: f64) -> f64 {
615    r * r * (r + 1e-30).ln()
616}
617
618/// Simple CG solver for the RBF system.
619fn rbf_cg(a: &[Vec<f64>], b: &[f64], tol: f64, max_iter: usize) -> Vec<f64> {
620    let n = b.len();
621    let mut x = vec![0.0_f64; n];
622    let mut r: Vec<f64> = b.to_vec();
623    let mut p = r.clone();
624    let mut rsold: f64 = r.iter().map(|v| v * v).sum();
625
626    for _ in 0..max_iter {
627        let ap: Vec<f64> = (0..n)
628            .map(|i| a[i].iter().zip(p.iter()).map(|(aij, pj)| aij * pj).sum())
629            .collect();
630        let pap: f64 = p.iter().zip(ap.iter()).map(|(pi, api)| pi * api).sum();
631        if pap.abs() < 1e-300 {
632            break;
633        }
634        let alpha = rsold / pap;
635        for i in 0..n {
636            x[i] += alpha * p[i];
637            r[i] -= alpha * ap[i];
638        }
639        let rsnew: f64 = r.iter().map(|v| v * v).sum();
640        if rsnew.sqrt() < tol {
641            break;
642        }
643        let beta = rsnew / rsold;
644        for i in 0..n {
645            p[i] = r[i] + beta * p[i];
646        }
647        rsold = rsnew;
648    }
649    x
650}
651
652// ---------------------------------------------------------------------------
653// § 8  Ray marching (sphere tracing)
654// ---------------------------------------------------------------------------
655
656/// Sphere-trace a ray against an implicit surface.
657///
658/// Starting from `origin` in direction `direction` (need not be unit), takes
659/// up to `max_steps` sphere-tracing steps.  Returns the distance along the
660/// ray to the first surface hit, or `None` if no hit was found within
661/// `max_dist`.
662///
663/// # Arguments
664/// * `sdf` — implicit surface to trace against
665/// * `origin` — ray origin in world space
666/// * `direction` — ray direction (will be normalized internally)
667/// * `max_steps` — maximum number of steps
668///
669/// The function terminates early when `|f(p)| < 1e-5`.
670pub fn ray_march(
671    sdf: &dyn ImplicitSurface,
672    origin: [f64; 3],
673    direction: [f64; 3],
674    max_steps: usize,
675) -> Option<f64> {
676    let dir = normalize3(direction);
677    let max_dist = 1e6_f64;
678    let mut t = 0.0_f64;
679
680    for _ in 0..max_steps {
681        let p = add3(origin, scale3(dir, t));
682        let d = sdf.eval(p);
683        if d.abs() < 1e-5 {
684            return Some(t);
685        }
686        if t > max_dist {
687            break;
688        }
689        // Advance by |d| (sphere-tracing step)
690        t += d.abs().max(1e-7);
691    }
692    None
693}
694
695// ---------------------------------------------------------------------------
696// § 9  Normal estimation
697// ---------------------------------------------------------------------------
698
699/// Estimate the outward unit normal of an implicit surface at `p` using
700/// central differences with step size `eps`.
701///
702/// # Arguments
703/// * `sdf` — implicit surface
704/// * `p` — query point (should be near or on the surface)
705/// * `eps` — finite difference step (typical: 1e-4 to 1e-6)
706pub fn surface_normal(sdf: &dyn ImplicitSurface, p: [f64; 3], eps: f64) -> [f64; 3] {
707    let gx = sdf.eval([p[0] + eps, p[1], p[2]]) - sdf.eval([p[0] - eps, p[1], p[2]]);
708    let gy = sdf.eval([p[0], p[1] + eps, p[2]]) - sdf.eval([p[0], p[1] - eps, p[2]]);
709    let gz = sdf.eval([p[0], p[1], p[2] + eps]) - sdf.eval([p[0], p[1], p[2] - eps]);
710    normalize3([gx, gy, gz])
711}
712
713// ---------------------------------------------------------------------------
714// § 10  Point-in-surface test (ray casting)
715// ---------------------------------------------------------------------------
716
717/// Test whether a point is inside a closed implicit surface.
718///
719/// Returns `true` if `sdf.eval(p) < 0`.
720///
721/// # Arguments
722/// * `sdf` — implicit surface
723/// * `p` — query point
724pub fn is_inside(sdf: &dyn ImplicitSurface, p: [f64; 3]) -> bool {
725    sdf.eval(p) < 0.0
726}
727
728// ---------------------------------------------------------------------------
729// § 11  Bounding-box computation
730// ---------------------------------------------------------------------------
731
732/// Compute the axis-aligned bounding box of a level-set surface by sampling.
733///
734/// Samples the SDF on a uniform grid within `bounds` at resolution `res`
735/// per axis and returns `[[xmin,xmax\], [ymin,ymax], [zmin,zmax]]` of all
736/// points where `|f(p)| < threshold`.
737pub fn surface_bbox(
738    sdf: &dyn ImplicitSurface,
739    bounds: [[f64; 2]; 3],
740    res: usize,
741    threshold: f64,
742) -> [[f64; 2]; 3] {
743    let n = res.max(2);
744    let mut bbox = [
745        [f64::INFINITY, f64::NEG_INFINITY],
746        [f64::INFINITY, f64::NEG_INFINITY],
747        [f64::INFINITY, f64::NEG_INFINITY],
748    ];
749
750    let step: [f64; 3] = std::array::from_fn(|k| (bounds[k][1] - bounds[k][0]) / (n - 1) as f64);
751
752    for iz in 0..n {
753        for iy in 0..n {
754            for ix in 0..n {
755                let p = [
756                    bounds[0][0] + ix as f64 * step[0],
757                    bounds[1][0] + iy as f64 * step[1],
758                    bounds[2][0] + iz as f64 * step[2],
759                ];
760                if sdf.eval(p).abs() < threshold {
761                    for k in 0..3 {
762                        bbox[k][0] = bbox[k][0].min(p[k]);
763                        bbox[k][1] = bbox[k][1].max(p[k]);
764                    }
765                }
766            }
767        }
768    }
769    bbox
770}
771
772// ---------------------------------------------------------------------------
773// Tests
774// ---------------------------------------------------------------------------
775
776#[cfg(test)]
777mod tests {
778    use super::*;
779
780    // ------------------------------------------------------------------
781    // SphereSDF
782    // ------------------------------------------------------------------
783
784    #[test]
785    fn test_sphere_eval_center_negative() {
786        let s = SphereSDF::new([0.0, 0.0, 0.0], 1.0);
787        assert!(s.eval([0.0, 0.0, 0.0]) < 0.0);
788    }
789
790    #[test]
791    fn test_sphere_eval_on_surface() {
792        let s = SphereSDF::new([0.0, 0.0, 0.0], 2.0);
793        let d = s.eval([2.0, 0.0, 0.0]);
794        assert!(d.abs() < 1e-12, "d = {d}");
795    }
796
797    #[test]
798    fn test_sphere_eval_outside_positive() {
799        let s = SphereSDF::new([0.0, 0.0, 0.0], 1.0);
800        assert!(s.eval([3.0, 0.0, 0.0]) > 0.0);
801    }
802
803    #[test]
804    fn test_sphere_gradient_unit_normal() {
805        let s = SphereSDF::new([0.0, 0.0, 0.0], 1.0);
806        let g = s.gradient([1.0, 0.0, 0.0]);
807        let l = len3(g);
808        assert!((l - 1.0).abs() < 1e-12, "gradient not unit: |g| = {l}");
809    }
810
811    #[test]
812    fn test_sphere_gradient_direction() {
813        let s = SphereSDF::new([0.0, 0.0, 0.0], 1.0);
814        let g = s.gradient([0.0, 1.0, 0.0]);
815        // Should point in +y direction
816        assert!(g[1] > 0.9);
817        assert!(g[0].abs() < 1e-12);
818        assert!(g[2].abs() < 1e-12);
819    }
820
821    #[test]
822    fn test_sphere_translated() {
823        let s = SphereSDF::new([1.0, 2.0, 3.0], 1.5);
824        // Distance from center should be -1.5
825        assert!((s.eval([1.0, 2.0, 3.0]) + 1.5).abs() < 1e-12);
826    }
827
828    // ------------------------------------------------------------------
829    // BoxSDF
830    // ------------------------------------------------------------------
831
832    #[test]
833    fn test_box_eval_inside_negative() {
834        let b = BoxSDF::new([1.0, 1.0, 1.0]);
835        assert!(b.eval([0.0, 0.0, 0.0]) < 0.0);
836    }
837
838    #[test]
839    fn test_box_eval_on_face() {
840        let b = BoxSDF::new([1.0, 2.0, 3.0]);
841        let d = b.eval([1.0, 0.0, 0.0]);
842        assert!(d.abs() < 1e-12, "d = {d}");
843    }
844
845    #[test]
846    fn test_box_eval_outside_positive() {
847        let b = BoxSDF::new([1.0, 1.0, 1.0]);
848        assert!(b.eval([2.0, 0.0, 0.0]) > 0.0);
849    }
850
851    #[test]
852    fn test_box_eval_corner() {
853        let b = BoxSDF::new([1.0, 1.0, 1.0]);
854        // Corner is at sqrt(3) distance from the box surface
855        let d = b.eval([2.0, 2.0, 2.0]);
856        let expected = (3.0_f64).sqrt();
857        assert!((d - expected).abs() < 1e-12);
858    }
859
860    // ------------------------------------------------------------------
861    // TorusSDF
862    // ------------------------------------------------------------------
863
864    #[test]
865    fn test_torus_eval_on_surface() {
866        let t = TorusSDF::new(2.0, 0.5);
867        // Outermost point of the torus
868        let d = t.eval([2.5, 0.0, 0.0]);
869        assert!(d.abs() < 1e-12, "d = {d}");
870    }
871
872    #[test]
873    fn test_torus_eval_center_positive() {
874        let t = TorusSDF::new(2.0, 0.5);
875        // The center of the torus is outside (hollow)
876        assert!(t.eval([0.0, 0.0, 0.0]) > 0.0);
877    }
878
879    #[test]
880    fn test_torus_inside_tube_negative() {
881        let t = TorusSDF::new(2.0, 0.5);
882        // Point near tube centre
883        assert!(t.eval([2.0, 0.0, 0.0]) < 0.0);
884    }
885
886    // ------------------------------------------------------------------
887    // CylinderSDF
888    // ------------------------------------------------------------------
889
890    #[test]
891    fn test_cylinder_lateral_surface() {
892        let c = CylinderSDF::new(1.0, 2.0);
893        let d = c.eval([1.0, 0.0, 0.0]);
894        assert!(d.abs() < 1e-12, "d = {d}");
895    }
896
897    #[test]
898    fn test_cylinder_inside_negative() {
899        let c = CylinderSDF::new(1.0, 2.0);
900        assert!(c.eval([0.0, 0.0, 0.0]) < 0.0);
901    }
902
903    #[test]
904    fn test_cylinder_outside_positive() {
905        let c = CylinderSDF::new(1.0, 2.0);
906        assert!(c.eval([3.0, 0.0, 0.0]) > 0.0);
907    }
908
909    // ------------------------------------------------------------------
910    // CSG operations
911    // ------------------------------------------------------------------
912
913    #[test]
914    fn test_union_inside_either() {
915        let a = SphereSDF::new([-1.0, 0.0, 0.0], 1.5);
916        let b = SphereSDF::new([1.0, 0.0, 0.0], 1.5);
917        let u = UnionSDF::new(a, b);
918        // Points inside both spheres
919        assert!(u.eval([-1.0, 0.0, 0.0]) < 0.0);
920        assert!(u.eval([1.0, 0.0, 0.0]) < 0.0);
921    }
922
923    #[test]
924    fn test_intersection_inside_both() {
925        let a = SphereSDF::new([0.0, 0.0, 0.0], 2.0);
926        let b = SphereSDF::new([0.0, 0.0, 0.0], 1.0);
927        let i = IntersectionSDF::new(a, b);
928        // Inside smaller sphere → inside intersection
929        assert!(i.eval([0.0, 0.0, 0.0]) < 0.0);
930    }
931
932    #[test]
933    fn test_intersection_outside_one() {
934        let a = SphereSDF::new([0.0, 0.0, 0.0], 2.0);
935        let b = SphereSDF::new([5.0, 0.0, 0.0], 1.0);
936        let i = IntersectionSDF::new(a, b);
937        // Origin is inside A but not inside B → outside intersection
938        assert!(i.eval([0.0, 0.0, 0.0]) > 0.0);
939    }
940
941    #[test]
942    fn test_difference_removes_region() {
943        let a = SphereSDF::new([0.0, 0.0, 0.0], 2.0);
944        let b = SphereSDF::new([0.0, 0.0, 0.0], 1.0);
945        let d = DifferenceSDF::new(a, b);
946        // Points inside b are removed
947        assert!(d.eval([0.5, 0.0, 0.0]) > 0.0);
948        // Points in a but not b remain
949        assert!(d.eval([1.5, 0.0, 0.0]) < 0.0);
950    }
951
952    // ------------------------------------------------------------------
953    // SmoothUnionSDF
954    // ------------------------------------------------------------------
955
956    #[test]
957    fn test_smooth_union_less_than_regular_union() {
958        let a = SphereSDF::new([-0.5, 0.0, 0.0], 1.0);
959        let b = SphereSDF::new([0.5, 0.0, 0.0], 1.0);
960        let su = SmoothUnionSDF::new(a.clone(), b.clone(), 0.5);
961        let u = UnionSDF::new(a, b);
962        let p = [0.0, 0.0, 0.0];
963        // Smooth union should be <= regular union (blending pulls inward)
964        assert!(su.eval(p) <= u.eval(p) + 1e-12);
965    }
966
967    // ------------------------------------------------------------------
968    // OffsetSDF
969    // ------------------------------------------------------------------
970
971    #[test]
972    fn test_offset_expands_sphere() {
973        let s = SphereSDF::new([0.0, 0.0, 0.0], 1.0);
974        let os = OffsetSDF::new(s, 0.5);
975        // New surface at r = 1.5
976        let d = os.eval([1.5, 0.0, 0.0]);
977        assert!(d.abs() < 1e-12, "d = {d}");
978    }
979
980    #[test]
981    fn test_offset_shrinks_sphere() {
982        let s = SphereSDF::new([0.0, 0.0, 0.0], 2.0);
983        let os = OffsetSDF::new(s, -0.5);
984        // New surface at r = 1.5
985        let d = os.eval([1.5, 0.0, 0.0]);
986        assert!(d.abs() < 1e-12, "d = {d}");
987    }
988
989    // ------------------------------------------------------------------
990    // is_inside / surface_normal
991    // ------------------------------------------------------------------
992
993    #[test]
994    fn test_is_inside_sphere() {
995        let s = SphereSDF::new([0.0, 0.0, 0.0], 1.0);
996        assert!(is_inside(&s, [0.0, 0.0, 0.0]));
997        assert!(!is_inside(&s, [2.0, 0.0, 0.0]));
998    }
999
1000    #[test]
1001    fn test_surface_normal_sphere_unit_length() {
1002        let s = SphereSDF::new([0.0, 0.0, 0.0], 1.0);
1003        let n = surface_normal(&s, [1.0, 0.0, 0.0], 1e-5);
1004        let l = len3(n);
1005        assert!((l - 1.0).abs() < 1e-3, "normal length = {l}");
1006    }
1007
1008    // ------------------------------------------------------------------
1009    // Marching cubes
1010    // ------------------------------------------------------------------
1011
1012    #[test]
1013    fn test_marching_cubes_sphere_produces_vertices() {
1014        let s = SphereSDF::new([0.0, 0.0, 0.0], 0.4);
1015        let bounds = [[-1.0, 1.0], [-1.0, 1.0], [-1.0, 1.0]];
1016        let (verts, tris) = marching_cubes(&s, bounds, 8);
1017        assert!(!verts.is_empty(), "no vertices produced");
1018        assert!(!tris.is_empty(), "no triangles produced");
1019    }
1020
1021    #[test]
1022    fn test_marching_cubes_indices_in_bounds() {
1023        let s = SphereSDF::new([0.0, 0.0, 0.0], 0.4);
1024        let bounds = [[-1.0, 1.0], [-1.0, 1.0], [-1.0, 1.0]];
1025        let (verts, tris) = marching_cubes(&s, bounds, 6);
1026        for tri in &tris {
1027            for &idx in tri {
1028                assert!(idx < verts.len(), "index {idx} out of range");
1029            }
1030        }
1031    }
1032
1033    #[test]
1034    fn test_marching_cubes_empty_for_no_surface() {
1035        // Sphere entirely outside the bounds
1036        let s = SphereSDF::new([100.0, 100.0, 100.0], 0.1);
1037        let bounds = [[-1.0, 1.0], [-1.0, 1.0], [-1.0, 1.0]];
1038        let (verts, tris) = marching_cubes(&s, bounds, 4);
1039        assert!(verts.is_empty());
1040        assert!(tris.is_empty());
1041    }
1042
1043    // ------------------------------------------------------------------
1044    // RBFImplicit
1045    // ------------------------------------------------------------------
1046
1047    #[test]
1048    fn test_rbf_fit_interpolates_at_centers() {
1049        let centers = vec![[0.0, 0.0, 0.0_f64], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]];
1050        let targets = vec![0.0, 1.0, -1.0];
1051        let mut rbf = RBFImplicit::new(centers.clone());
1052        rbf.fit(&centers, &targets, 1e-10, 200);
1053        // After fitting, evaluating at centers should approximate targets
1054        for (c, &t) in centers.iter().zip(targets.iter()) {
1055            let v = rbf.eval(*c);
1056            assert!((v - t).abs() < 0.5, "RBF({c:?}) = {v}, expected {t}");
1057        }
1058    }
1059
1060    #[test]
1061    fn test_rbf_zero_weights_before_fit() {
1062        let centers = vec![[0.0, 0.0, 0.0_f64]];
1063        let rbf = RBFImplicit::new(centers);
1064        assert_eq!(rbf.weights, vec![0.0]);
1065    }
1066
1067    // ------------------------------------------------------------------
1068    // Ray marching
1069    // ------------------------------------------------------------------
1070
1071    #[test]
1072    fn test_ray_march_sphere_hit() {
1073        let s = SphereSDF::new([0.0, 0.0, 0.0], 1.0);
1074        let hit = ray_march(&s, [-5.0, 0.0, 0.0], [1.0, 0.0, 0.0], 200);
1075        assert!(hit.is_some(), "ray should hit sphere");
1076        let t = hit.unwrap();
1077        // Hit point should be near [-1, 0, 0]
1078        let px = -5.0 + t;
1079        assert!((px + 1.0).abs() < 0.01, "hit at x = {px}, expected -1");
1080    }
1081
1082    #[test]
1083    fn test_ray_march_sphere_miss() {
1084        let s = SphereSDF::new([0.0, 0.0, 0.0], 1.0);
1085        // Ray aimed far away from sphere
1086        let hit = ray_march(&s, [0.0, 100.0, 0.0], [1.0, 0.0, 0.0], 50);
1087        assert!(hit.is_none(), "ray should miss sphere");
1088    }
1089
1090    #[test]
1091    fn test_ray_march_returns_positive_distance() {
1092        let s = SphereSDF::new([0.0, 0.0, 0.0], 1.0);
1093        let hit = ray_march(&s, [-5.0, 0.0, 0.0], [1.0, 0.0, 0.0], 200);
1094        if let Some(t) = hit {
1095            assert!(t > 0.0, "distance must be positive, got {t}");
1096        }
1097    }
1098}
oxiphysics_geometry/implicit_surfaces.rs

oxiphysics_geometry/
implicit_surfaces.rs
