mesh-geometry 0.4.0

High-performance geometry utilities for coastal/ocean FVM
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58

//! Point-in-polygon queries for mesh-geometry.

use crate::{Float, Point2};

/// Returns true if `pt` is inside the planar polygon `poly` (winding-number).
/// Assumes `poly` is closed (first != last) and simple (non-self-intersecting).
pub fn point_in_polygon<T: Float>(pt: Point2<T>, poly: &[Point2<T>]) -> bool {
    let mut winding: i32 = 0;
    let n = poly.len();
    for i in 0..n {
        let p1 = poly[i];
        let p2 = poly[(i + 1) % n];
        if p1.y <= pt.y {
            if p2.y > pt.y && is_left(p1, p2, pt) > T::zero() {
                winding += 1;
            }
        } else {
            if p2.y <= pt.y && is_left(p1, p2, pt) < T::zero() {
                winding -= 1;
            }
        }
    }
    winding != 0
}

#[inline]
fn is_left<T: Float>(a: Point2<T>, b: Point2<T>, c: Point2<T>) -> T {
    (b.x - a.x) * (c.y - a.y) - (c.x - a.x) * (b.y - a.y)
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::Point2;

    #[test]
    fn inside_triangle() {
        let tri = [
            Point2::new(0.0, 0.0),
            Point2::new(5.0, 0.0),
            Point2::new(0.0, 5.0),
        ];
        assert!(point_in_polygon(Point2::new(1.0,1.0), &tri));
        assert!(!point_in_polygon(Point2::new(5.0,5.0), &tri));
    }

    #[test]
    fn inside_square_and_hole() {
        let square = [
            Point2::new(-1.0, -1.0),
            Point2::new( 1.0, -1.0),
            Point2::new( 1.0,  1.0),
            Point2::new(-1.0,  1.0),
        ];
        assert!(point_in_polygon(Point2::new(0.0,0.0), &square));
        assert!(!point_in_polygon(Point2::new(2.0,0.0), &square));
    }
}