fenris_geometry/
sdf.rs

1use fenris_traits::Real;
2use nalgebra::{OPoint, Point2, Scalar, Vector2, U2};
3
4use crate::{AxisAlignedBoundingBox2d, BoundedGeometry};
5use numeric_literals::replace_float_literals;
6
7pub trait SignedDistanceFunction2d<T>
8where
9    T: Scalar,
10{
11    fn eval(&self, x: &Point2<T>) -> T;
12    fn gradient(&self, x: &Point2<T>) -> Option<Vector2<T>>;
13
14    fn union<Other>(self, other: Other) -> SdfUnion<Self, Other>
15    where
16        Self: Sized,
17        Other: Sized + SignedDistanceFunction2d<T>,
18    {
19        SdfUnion {
20            left: self,
21            right: other,
22        }
23    }
24}
25
26pub trait BoundedSdf<T>: SignedDistanceFunction2d<T> + BoundedGeometry<T, Dimension = U2>
27where
28    T: Scalar,
29{
30}
31
32impl<X, T> BoundedSdf<T> for X
33where
34    T: Scalar,
35    X: SignedDistanceFunction2d<T> + BoundedGeometry<T, Dimension = U2>,
36{
37}
38
39#[derive(Copy, Clone, Debug)]
40pub struct SdfCircle<T>
41where
42    T: Scalar,
43{
44    pub radius: T,
45    pub center: Vector2<T>,
46}
47
48#[derive(Copy, Clone, Debug)]
49pub struct SdfUnion<Left, Right> {
50    pub left: Left,
51    pub right: Right,
52}
53
54#[derive(Copy, Clone, Debug)]
55pub struct SdfAxisAlignedBox<T>
56where
57    T: Scalar,
58{
59    pub aabb: AxisAlignedBoundingBox2d<T>,
60}
61
62impl<T> BoundedGeometry<T> for SdfCircle<T>
63where
64    T: Real,
65{
66    type Dimension = U2;
67
68    fn bounding_box(&self) -> AxisAlignedBoundingBox2d<T> {
69        let eps = self.radius * T::from_f64(0.01).unwrap();
70        AxisAlignedBoundingBox2d::new(
71            OPoint::from(self.center - Vector2::repeat(T::one()) * (self.radius + eps)),
72            OPoint::from(self.center + Vector2::repeat(T::one()) * (self.radius + eps)),
73        )
74    }
75}
76
77impl<T> SignedDistanceFunction2d<T> for SdfCircle<T>
78where
79    T: Real,
80{
81    fn eval(&self, x: &Point2<T>) -> T {
82        let y = x - self.center;
83        y.coords.norm() - self.radius
84    }
85
86    fn gradient(&self, x: &Point2<T>) -> Option<Vector2<T>> {
87        let y = x - self.center;
88        let y_norm = y.coords.norm();
89
90        if y_norm == T::zero() {
91            None
92        } else {
93            Some(y.coords / y_norm)
94        }
95    }
96}
97
98impl<T, Left, Right> BoundedGeometry<T> for SdfUnion<Left, Right>
99where
100    T: Real,
101    Left: BoundedGeometry<T, Dimension = U2>,
102    Right: BoundedGeometry<T, Dimension = U2>,
103{
104    type Dimension = U2;
105
106    fn bounding_box(&self) -> AxisAlignedBoundingBox2d<T> {
107        self.left.bounding_box().enclose(&self.right.bounding_box())
108    }
109}
110
111impl<T, Left, Right> SignedDistanceFunction2d<T> for SdfUnion<Left, Right>
112where
113    T: Real,
114    Left: SignedDistanceFunction2d<T>,
115    Right: SignedDistanceFunction2d<T>,
116{
117    fn eval(&self, x: &Point2<T>) -> T {
118        self.left.eval(x).min(self.right.eval(x))
119    }
120
121    fn gradient(&self, x: &Point2<T>) -> Option<Vector2<T>> {
122        // TODO: Is this actually correct? It might give funky results if exactly
123        // at points where either SDF is non-differentiable
124        if self.left.eval(x) < self.right.eval(x) {
125            self.left.gradient(x)
126        } else {
127            self.right.gradient(x)
128        }
129    }
130}
131
132impl<T> BoundedGeometry<T> for SdfAxisAlignedBox<T>
133where
134    T: Real,
135{
136    type Dimension = U2;
137
138    fn bounding_box(&self) -> AxisAlignedBoundingBox2d<T> {
139        self.aabb
140    }
141}
142
143impl<T> SignedDistanceFunction2d<T> for SdfAxisAlignedBox<T>
144where
145    T: Real,
146{
147    #[replace_float_literals(T::from_f64(literal).unwrap())]
148    fn eval(&self, x: &Point2<T>) -> T {
149        let b = self.aabb.extents() / 2.0;
150        let p = x - self.aabb.center();
151        let d = p.abs() - b;
152
153        // TODO: Use d.max() when fixed. See https://github.com/rustsim/nalgebra/issues/620
154        d.sup(&Vector2::zeros()).norm() + T::min(T::zero(), d[d.imax()])
155    }
156
157    #[replace_float_literals(T::from_f64(literal).expect("Literal must fit in T"))]
158    fn gradient(&self, x: &Point2<T>) -> Option<Vector2<T>> {
159        // TODO: Replace finite differences with "proper" gradient
160        // Note: arbitrary "step"/resolution h
161        let h = 1e-4;
162        let mut gradient = Vector2::zeros();
163        for i in 0..2 {
164            let mut dx = Vector2::zeros();
165            dx[i] = h;
166            gradient[i] = (self.eval(&(x + dx)) - self.eval(&(x - dx))) / (2.0 * h)
167        }
168        gradient.normalize_mut();
169        Some(gradient)
170    }
171}