oxiphysics_geometry/
implicit_surfaces.rs

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

oxiphysics_geometry/
implicit_surfaces.rs
