axiolid_overlay/
visibility.rs

1//! The visibility polygon of a point inside a region with holes (#184).
2//!
3//! # A sweep with exact decisions
4//!
5//! Around the viewpoint `v` every boundary vertex has a direction. Between
6//! two consecutive directions -- a wedge -- no vertex lies, and boundary
7//! edges do not cross, so one edge is nearest throughout the wedge and
8//! bounds what is seen there. The visibility polygon is the fan of those
9//! pieces.
10//!
11//! Every choice is exact: the order of directions (a half-plane and an
12//! orientation), which edges a wedge's middle ray meets, and which of them
13//! is nearest along it. The middle ray points along `(w1 - v) + (w2 - v)`,
14//! the sum of its two bounding directions, so it is exact too; a wedge is
15//! narrower than a half-turn whenever the viewpoint is strictly inside, so
16//! that sum points into it. Only the output points where a wedge's
17//! boundary ray meets its edge -- the ends of shadows -- are rounded, once,
18//! and the ring is presented like every other region operation.
19
20use axiolid_core::{Point2, Tolerance};
21use axiolid_exact::{certify, Arith, SignExpr};
22use axiolid_guarantees::Sign;
23
24use crate::arc::ArcRing;
25use crate::arrangement::ArcArrangement;
26use crate::region::Region;
27use crate::OverlayError;
28
29/// Why no visibility polygon was built.
30#[derive(Debug, Clone, PartialEq, Eq)]
31#[non_exhaustive]
32pub enum VisibilityError {
33    /// The viewpoint is on the region's boundary or outside it.
34    NotInside,
35    /// The result failed the overlay's own checks.
36    Overlay(OverlayError),
37}
38
39impl From<OverlayError> for VisibilityError {
40    fn from(error: OverlayError) -> Self {
41        Self::Overlay(error)
42    }
43}
44
45impl Region {
46    /// The part of the region in sight of `viewpoint`: every point the
47    /// straight segment from the viewpoint reaches without leaving the
48    /// region. Walls and holes cast shadows.
49    ///
50    /// # Errors
51    ///
52    /// [`VisibilityError::NotInside`] for a viewpoint on the boundary or
53    /// outside.
54    pub fn visibility_polygon(
55        &self,
56        viewpoint: Point2,
57        tolerance: Tolerance,
58    ) -> Result<Self, VisibilityError> {
59        if !viewpoint.is_finite() || !strictly_inside(self, viewpoint) {
60            return Err(VisibilityError::NotInside);
61        }
62        let edges: Vec<(Point2, Point2)> = self
63            .boundary_rings()
64            .iter()
65            .flat_map(|r| {
66                let n = r.points.len();
67                (0..n).map(move |i| (r.points[i], r.points[(i + 1) % n]))
68            })
69            .filter(|(a, b)| a != b)
70            .collect();
71        let v = viewpoint;
72        // One vertex per direction, counter-clockwise from +x.
73        let mut around: Vec<Point2> = Vec::new();
74        for &(a, _) in &edges {
75            let at = around.partition_point(|&w| before(v, w, a));
76            if at < around.len() && same_direction(v, around[at], a) {
77                continue;
78            }
79            around.insert(at, a);
80        }
81        let k = around.len();
82        let mut ring: Vec<Point2> = Vec::with_capacity(2 * k);
83        for i in 0..k {
84            let (w1, w2) = (around[i], around[(i + 1) % k]);
85            let Some(edge) = nearest(v, w1, w2, &edges) else {
86                // Every wedge meets the boundary when the viewpoint is inside.
87                return Err(VisibilityError::NotInside);
88            };
89            for w in [w1, w2] {
90                let p = on_ray(v, w, edge);
91                if ring.last() != Some(&p) {
92                    ring.push(p);
93                }
94            }
95        }
96        while ring.len() > 1 && ring.first() == ring.last() {
97            ring.pop();
98        }
99        if ring.len() < 3 {
100            return Ok(Self::empty());
101        }
102        let arrangement = ArcArrangement::new(&[ArcRing::from_points(&ring)], tolerance)?;
103        Ok(crate::minkowski::region_of(
104            &arrangement,
105            |f| f[0],
106            true,
107            tolerance,
108        )?)
109    }
110}
111
112/// Where the ray from `v` towards `w` meets `edge`: an endpoint exactly
113/// when it lies on the ray, else rounded once.
114fn on_ray(v: Point2, w: Point2, (a, b): (Point2, Point2)) -> Point2 {
115    for end in [a, b] {
116        if sign(&Orient { a: v, b: w, c: end }) == Sign::Zero {
117            return end;
118        }
119    }
120    let (d, e) = (w - v, b - a);
121    let t = (a - v).perp_dot(e) / d.perp_dot(e);
122    v + d * t
123}
124
125/// The nearest edge along the wedge's middle ray from `v`, between the
126/// directions to `w1` and `w2`.
127fn nearest(
128    v: Point2,
129    w1: Point2,
130    w2: Point2,
131    edges: &[(Point2, Point2)],
132) -> Option<(Point2, Point2)> {
133    let mut best: Option<(Point2, Point2)> = None;
134    for &edge in edges {
135        if sign(&Hits { v, w1, w2, edge }) != Sign::Positive {
136            continue;
137        }
138        if best.is_none_or(|b| {
139            sign(&Nearer {
140                v,
141                w1,
142                w2,
143                near: edge,
144                far: b,
145            }) == Sign::Positive
146        }) {
147            best = Some(edge);
148        }
149    }
150    best
151}
152
153/// Whether the direction from `v` to `a` comes strictly before that to
154/// `b`, counter-clockwise from `+x`.
155fn before(v: Point2, a: Point2, b: Point2) -> bool {
156    let half = |p: Point2| u8::from(!(p.y > v.y || (p.y == v.y && p.x > v.x)));
157    let (ha, hb) = (half(a), half(b));
158    if ha != hb {
159        return ha < hb;
160    }
161    sign(&Orient { a: v, b: a, c: b }) == Sign::Positive
162}
163
164fn same_direction(v: Point2, a: Point2, b: Point2) -> bool {
165    !before(v, a, b) && !before(v, b, a)
166}
167
168/// Strictly inside the region: off every boundary ring and with an odd
169/// number of rings round it, exactly.
170fn strictly_inside(region: &Region, p: Point2) -> bool {
171    let mut winding = 0usize;
172    for ring in region.boundary_rings() {
173        let n = ring.points.len();
174        let mut crossings = 0usize;
175        for i in 0..n {
176            let (a, b) = (ring.points[i], ring.points[(i + 1) % n]);
177            let s = sign(&Orient { a, b, c: p });
178            let between = p.x >= a.x.min(b.x)
179                && p.x <= a.x.max(b.x)
180                && p.y >= a.y.min(b.y)
181                && p.y <= a.y.max(b.y);
182            if s == Sign::Zero && between {
183                return false;
184            }
185            // Crossings of the rightward ray, half-open in y.
186            let upward = a.y <= p.y && b.y > p.y;
187            let downward = b.y <= p.y && a.y > p.y;
188            if (upward && s == Sign::Positive) || (downward && s == Sign::Negative) {
189                crossings += 1;
190            }
191        }
192        winding += crossings % 2;
193    }
194    winding % 2 == 1
195}
196
197/// The inputs are finite, so the exact tier always decides.
198fn sign<E: SignExpr>(e: &E) -> Sign {
199    certify(e).unwrap_or(Sign::Zero)
200}
201
202struct Orient {
203    a: Point2,
204    b: Point2,
205    c: Point2,
206}
207
208impl SignExpr for Orient {
209    fn sign_in<T: Arith>(&self) -> Option<Sign> {
210        let f = T::from_f64;
211        let (ux, uy) = (f(self.b.x).sub(&f(self.a.x)), f(self.b.y).sub(&f(self.a.y)));
212        let (vx, vy) = (f(self.c.x).sub(&f(self.a.x)), f(self.c.y).sub(&f(self.a.y)));
213        ux.mul(&vy).sub(&uy.mul(&vx)).sign()
214    }
215}
216
217/// `(p - q)` in `T`.
218fn diff<T: Arith>(p: Point2, q: Point2) -> [T; 2] {
219    [
220        T::from_f64(p.x).sub(&T::from_f64(q.x)),
221        T::from_f64(p.y).sub(&T::from_f64(q.y)),
222    ]
223}
224
225fn cross<T: Arith>(a: &[T; 2], b: &[T; 2]) -> T {
226    a[0].mul(&b[1]).sub(&a[1].mul(&b[0]))
227}
228
229/// The middle direction of a wedge, `(w1 - v) + (w2 - v)`.
230fn middle<T: Arith>(v: Point2, w1: Point2, w2: Point2) -> [T; 2] {
231    let (p, q) = (diff::<T>(w1, v), diff::<T>(w2, v));
232    [p[0].add(&q[0]), p[1].add(&q[1])]
233}
234
235/// The ray from `v` along `u` meets the line of `(a, b)` at parameter
236/// `N / D`: `N = (a - v) x (b - a)`, `D = u x (b - a)`.
237fn parameter<T: Arith>(v: Point2, u: &[T; 2], (a, b): (Point2, Point2)) -> (T, T) {
238    let e = diff::<T>(b, a);
239    (cross(&diff::<T>(a, v), &e), cross(u, &e))
240}
241
242/// Positive when the wedge's middle ray crosses the edge's interior ahead
243/// of the viewpoint.
244struct Hits {
245    v: Point2,
246    w1: Point2,
247    w2: Point2,
248    edge: (Point2, Point2),
249}
250
251impl SignExpr for Hits {
252    fn sign_in<T: Arith>(&self) -> Option<Sign> {
253        let u = middle::<T>(self.v, self.w1, self.w2);
254        let (a, b) = self.edge;
255        let sa = cross(&u, &diff::<T>(a, self.v)).sign()?;
256        let sb = cross(&u, &diff::<T>(b, self.v)).sign()?;
257        let straddles = matches!(
258            (sa, sb),
259            (Sign::Positive, Sign::Negative) | (Sign::Negative, Sign::Positive)
260        );
261        if !straddles {
262            return Some(Sign::Negative);
263        }
264        let (n, d) = parameter::<T>(self.v, &u, self.edge);
265        let (sn, sd) = (n.sign()?, d.sign()?);
266        Some(if sn != Sign::Zero && sn == sd {
267            Sign::Positive
268        } else {
269            Sign::Negative
270        })
271    }
272}
273
274/// Positive when `near` is strictly nearer than `far` along the wedge's
275/// middle ray; both are known to cross it ahead.
276struct Nearer {
277    v: Point2,
278    w1: Point2,
279    w2: Point2,
280    near: (Point2, Point2),
281    far: (Point2, Point2),
282}
283
284impl SignExpr for Nearer {
285    fn sign_in<T: Arith>(&self) -> Option<Sign> {
286        let u = middle::<T>(self.v, self.w1, self.w2);
287        let (n1, d1) = parameter::<T>(self.v, &u, self.near);
288        let (n2, d2) = parameter::<T>(self.v, &u, self.far);
289        // t1 < t2  <=>  (n2 d1 - n1 d2) has the sign of d1 d2.
290        let gap = n2.mul(&d1).sub(&n1.mul(&d2)).sign()?;
291        let same = d1.sign()? == d2.sign()?;
292        Some(match (gap, same) {
293            (Sign::Zero, _) => Sign::Zero,
294            (s, true) => s,
295            (Sign::Positive, false) => Sign::Negative,
296            (_, false) => Sign::Positive,
297        })
298    }
299}
axiolid_overlay/visibility.rs

axiolid_overlay/
visibility.rs
