algebraeon_geometry/
oriented_simplex.rs

1use crate::{
2    ambient_space::AffineSpace,
3    simplex::Simplex,
4    vector::{DotProduct, Vector},
5};
6use algebraeon_rings::{
7    matrix::{Matrix, MatrixStructure},
8    structure::{FieldSignature, OrderedRingSignature},
9};
10
11#[derive(Clone)]
12pub struct OrientedSimplex<'f, FS: OrderedRingSignature + FieldSignature> {
13    simplex: Simplex<'f, FS>,
14    orientation: Option<OrientedSimplexOrientation<'f, FS>>, //None iff simplex is null
15}
16
17impl<'f, FS: OrderedRingSignature + FieldSignature> std::fmt::Debug for OrientedSimplex<'f, FS> {
18    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
19        f.debug_struct("OrientedSimplex")
20            .field("simplex", &self.simplex)
21            .field("orientation", &self.orientation)
22            .finish()
23    }
24}
25
26impl<'f, FS: OrderedRingSignature + FieldSignature> OrientedSimplex<'f, FS> {
27    pub fn new_with_positive_point(
28        ambient_space: AffineSpace<'f, FS>,
29        points: Vec<Vector<'f, FS>>,
30        ref_point: &Vector<'f, FS>,
31    ) -> Result<Self, &'static str> {
32        assert_eq!(ref_point.ambient_space(), ambient_space);
33        if points.len() != ambient_space.linear_dimension().unwrap() {
34            return Err("Oriented simplex must have dimension one less than the ambient space");
35        }
36        let n = points.len();
37        if n == 0 {
38            Ok(Self {
39                simplex: ambient_space.simplex(points)?,
40                orientation: None,
41            })
42        } else {
43            let root = &points[0];
44            let positive_normal = {
45                let mat = Matrix::construct(n - 1, n, |r, c| {
46                    ambient_space
47                        .field()
48                        .sub(points[r + 1].coordinate(c), points[0].coordinate(c))
49                });
50                let kernel = MatrixStructure::<FS, &FS>::new(ambient_space.field())
51                    .col_kernel(mat)
52                    .basis();
53                if kernel.len() != 1 {
54                    return Err("points are not affine independent");
55                }
56                ambient_space.vector(kernel.into_iter().next().unwrap())
57            };
58            let flip = match ambient_space.field().cmp(
59                &(ref_point - root).dot(&positive_normal),
60                &ambient_space.field().zero(),
61            ) {
62                std::cmp::Ordering::Less => true,
63                std::cmp::Ordering::Equal => {
64                    return Err("ref_point lines inside the hyperplane");
65                }
66                std::cmp::Ordering::Greater => false,
67            };
68            let plane_point = root.clone();
69
70            Ok(Self {
71                simplex: ambient_space.simplex(points)?,
72                orientation: Some(OrientedSimplexOrientation {
73                    flip,
74                    plane_point,
75                    positive_normal,
76                }),
77            })
78        }
79    }
80
81    pub fn new_with_negative_point(
82        ambient_space: AffineSpace<'f, FS>,
83        points: Vec<Vector<'f, FS>>,
84        ref_point: &Vector<'f, FS>,
85    ) -> Result<Self, &'static str> {
86        let mut ans = Self::new_with_positive_point(ambient_space, points, ref_point)?;
87        ans.flip();
88        Ok(ans)
89    }
90
91    pub fn positive_point(&self) -> Option<Vector<'f, FS>> {
92        if let Some(orientation) = self.orientation.as_ref() {
93            match orientation.flip {
94                false => Some(&orientation.plane_point + &orientation.positive_normal),
95                true => Some(&orientation.plane_point - &orientation.positive_normal),
96            }
97        } else {
98            None
99        }
100    }
101
102    pub fn negative_point(&self) -> Option<Vector<'f, FS>> {
103        let positive_point = self.positive_point()?;
104        Some({
105            let pt: &Vector<'f, FS> = &self.simplex.points()[0];
106            &(pt + pt) - &positive_point
107        })
108    }
109
110    pub fn ambient_space(&self) -> AffineSpace<'f, FS> {
111        self.simplex().ambient_space()
112    }
113
114    pub fn simplex(&self) -> &Simplex<'f, FS> {
115        &self.simplex
116    }
117
118    pub fn into_simplex(self) -> Simplex<'f, FS> {
119        self.simplex
120    }
121
122    pub fn into_oriented_hyperplane(self) -> OrientedHyperplane<'f, FS> {
123        OrientedHyperplane {
124            oriented_simplex: self,
125        }
126    }
127
128    pub fn flip(&mut self) {
129        if let Some(OrientedSimplexOrientation { flip, .. }) = &mut self.orientation {
130            *flip = !*flip;
131        }
132        if let Some(OrientedSimplexOrientation {
133            flip,
134            plane_point,
135            positive_normal,
136        }) = &self.orientation
137        {
138            debug_assert_eq!(
139                self.classify_point(&(plane_point + positive_normal)),
140                match flip {
141                    false => OrientationSide::Positive,
142                    true => OrientationSide::Negative,
143                }
144            );
145        }
146    }
147
148    fn classify_point_quantitatively(&self, point: &Vector<'f, FS>) -> FS::Set {
149        match &self.orientation {
150            Some(OrientedSimplexOrientation {
151                flip,
152                plane_point,
153                positive_normal,
154            }) => {
155                let mut value = positive_normal.dot(&(point - plane_point));
156                if *flip {
157                    value = self.ambient_space().field().neg(&value);
158                }
159                value
160            }
161            None => self.ambient_space().field().zero(),
162        }
163    }
164
165    pub fn classify_point(&self, point: &Vector<'f, FS>) -> OrientationSide {
166        let space = self.ambient_space();
167        let field = space.field();
168        let value = self.classify_point_quantitatively(point);
169        match field.cmp(&value, &field.zero()) {
170            std::cmp::Ordering::Less => OrientationSide::Negative,
171            std::cmp::Ordering::Equal => OrientationSide::Neutral,
172            std::cmp::Ordering::Greater => OrientationSide::Positive,
173        }
174    }
175}
176
177#[derive(Debug, Clone, Copy, PartialEq, Eq)]
178pub enum OrientationSide {
179    Positive,
180    Neutral,
181    Negative,
182}
183
184/// a simplex spanning a hyperplane with a positive side and a negative side
185/// in 3d space it is a triangle
186/// in 2d space it is a line
187/// in 1d space it is a point
188/// in 0d space it is a null simplex. The orientation looses meaning but it is nice to still count this case.
189#[derive(Debug, Clone)]
190struct OrientedSimplexOrientation<'f, FS: OrderedRingSignature + FieldSignature> {
191    flip: bool,                      // flip the orientation if necessary
192    plane_point: Vector<'f, FS>,     // A point on the hyperplane
193    positive_normal: Vector<'f, FS>, // normal vector pointing out of the positive side
194}
195
196#[derive(Clone)]
197pub struct OrientedHyperplane<'f, FS: OrderedRingSignature + FieldSignature> {
198    oriented_simplex: OrientedSimplex<'f, FS>,
199}
200
201impl<'f, FS: OrderedRingSignature + FieldSignature> std::fmt::Debug for OrientedHyperplane<'f, FS> {
202    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
203        f.debug_struct("OrientedHyperplane")
204            .field("oriented_simplex", &self.oriented_simplex)
205            .finish()
206    }
207}
208
209// impl<FS: OrderedRingStructure + FieldStructure, AffineSpace<'f, FS>: Borrow<AffineSpace<FS>> + Clone>
210//     From<OrientedSimplex<FS>> for OrientedHyperplane<FS>
211// {
212//     fn from(oriented_simplex: OrientedSimplex<FS>) -> Self {
213//         Self { oriented_simplex }
214//     }
215// }
216
217pub enum OrientedHyperplaneIntersectLineSegmentResult<'f, FS: OrderedRingSignature + FieldSignature>
218{
219    PositivePositive,
220    NegativeNegative,
221    PositiveNeutral,
222    NegativeNeutral,
223    NeutralPositive,
224    NeutralNegative,
225    NeutralNeutral,
226    PositiveNegative { intersection_point: Vector<'f, FS> },
227    NegativePositive { intersection_point: Vector<'f, FS> },
228}
229
230impl<'f, FS: OrderedRingSignature + FieldSignature> OrientedHyperplane<'f, FS> {
231    pub fn ambient_space(&self) -> AffineSpace<'f, FS> {
232        self.oriented_simplex.ambient_space()
233    }
234
235    pub fn classify_point(&self, point: &Vector<'f, FS>) -> OrientationSide {
236        self.oriented_simplex.classify_point(point)
237    }
238
239    pub fn intersect_line(
240        &self,
241        a: &Vector<'f, FS>,
242        b: &Vector<'f, FS>,
243    ) -> OrientedHyperplaneIntersectLineSegmentResult<'f, FS> {
244        let space = self.ambient_space();
245        let field = space.field();
246
247        let a_val = self.oriented_simplex.classify_point_quantitatively(a);
248        let b_val = self.oriented_simplex.classify_point_quantitatively(b);
249
250        match (
251            field.cmp(&a_val, &field.zero()),
252            field.cmp(&b_val, &field.zero()),
253        ) {
254            (std::cmp::Ordering::Less, std::cmp::Ordering::Greater) => {
255                let t = field
256                    .try_divide(&a_val, &field.add(&a_val, &field.neg(&b_val)))
257                    .unwrap();
258                {
259                    debug_assert_eq!(
260                        self.oriented_simplex.classify_point(a),
261                        OrientationSide::Negative
262                    );
263                    debug_assert_eq!(
264                        self.oriented_simplex.classify_point(b),
265                        OrientationSide::Positive
266                    );
267                    OrientedHyperplaneIntersectLineSegmentResult::NegativePositive {
268                        intersection_point: a + &(b - a).scalar_mul(&t),
269                    }
270                }
271            }
272            (std::cmp::Ordering::Greater, std::cmp::Ordering::Less) => {
273                let t = field
274                    .try_divide(&a_val, &field.add(&a_val, &field.neg(&b_val)))
275                    .unwrap();
276                {
277                    debug_assert_eq!(
278                        self.oriented_simplex.classify_point(a),
279                        OrientationSide::Positive
280                    );
281                    debug_assert_eq!(
282                        self.oriented_simplex.classify_point(b),
283                        OrientationSide::Negative
284                    );
285                    OrientedHyperplaneIntersectLineSegmentResult::PositiveNegative {
286                        intersection_point: a + &(b - a).scalar_mul(&t),
287                    }
288                }
289            }
290            (std::cmp::Ordering::Less, std::cmp::Ordering::Less) => {
291                OrientedHyperplaneIntersectLineSegmentResult::NegativeNegative
292            }
293            (std::cmp::Ordering::Less, std::cmp::Ordering::Equal) => {
294                OrientedHyperplaneIntersectLineSegmentResult::NegativeNeutral
295            }
296            (std::cmp::Ordering::Equal, std::cmp::Ordering::Less) => {
297                OrientedHyperplaneIntersectLineSegmentResult::NeutralNegative
298            }
299            (std::cmp::Ordering::Equal, std::cmp::Ordering::Equal) => {
300                OrientedHyperplaneIntersectLineSegmentResult::NeutralNeutral
301            }
302            (std::cmp::Ordering::Equal, std::cmp::Ordering::Greater) => {
303                OrientedHyperplaneIntersectLineSegmentResult::NeutralPositive
304            }
305            (std::cmp::Ordering::Greater, std::cmp::Ordering::Equal) => {
306                OrientedHyperplaneIntersectLineSegmentResult::PositiveNeutral
307            }
308            (std::cmp::Ordering::Greater, std::cmp::Ordering::Greater) => {
309                OrientedHyperplaneIntersectLineSegmentResult::PositivePositive
310            }
311        }
312    }
313
314    pub fn flip(&mut self) {
315        self.oriented_simplex.flip();
316    }
317}
318
319#[cfg(test)]
320mod tests {
321    use algebraeon_nzq::Rational;
322
323    use super::*;
324
325    #[test]
326    fn make_oriented_simplex() {
327        let space = AffineSpace::new_linear(Rational::structure_ref(), 2);
328        let v1 = space.vector([1, 0]);
329        let v2 = space.vector([0, 1]);
330        let v3 = space.vector([2, 3]);
331        let s_pos =
332            OrientedSimplex::new_with_positive_point(space, vec![v1.clone(), v2.clone()], &v3)
333                .unwrap();
334        let s_neg = OrientedSimplex::new_with_negative_point(space, vec![v1, v2], &v3).unwrap();
335
336        assert_ne!(
337            s_pos.orientation.unwrap().flip,
338            s_neg.orientation.unwrap().flip
339        );
340    }
341}
algebraeon_geometry/oriented_simplex.rs

algebraeon_geometry/
oriented_simplex.rs
