kasane_logic/geometry/shapes/polygon/
geometry_relation.rs

1use crate::{Coordinate, Ecef, IntoCoordinates, IntoLines, IntoTriangles, Line, Polygon, Triangle};
2
3impl IntoCoordinates for Polygon {
4    fn into_coordinates(self) -> impl Iterator<Item = Coordinate> {
5        self.vertices.into_iter()
6    }
7
8    fn iter_coordinates(&self) -> impl Iterator<Item = Coordinate> {
9        self.vertices.clone().into_iter()
10    }
11}
12
13impl IntoLines for Polygon {
14    fn into_lines(self) -> impl Iterator<Item = Line> {
15        let triangles: Vec<Triangle> = self.into_triangles().collect();
16        triangles
17            .into_iter()
18            .flat_map(|triangle| triangle.into_lines())
19    }
20
21    fn iter_lines(&self) -> impl Iterator<Item = Line> {
22        let triangles: Vec<Triangle> = self.iter_triangles().collect();
23        triangles
24            .into_iter()
25            .flat_map(|triangle| triangle.into_lines())
26    }
27}
28
29impl IntoTriangles for Polygon {
30    fn into_triangles(self) -> impl Iterator<Item = Triangle> {
31        self.iter_triangles().collect::<Vec<_>>().into_iter()
32    }
33
34    fn iter_triangles(&self) -> impl Iterator<Item = Triangle> {
35        let n = self.vertices.len();
36        if n < 3 {
37            return Vec::<Triangle>::new().into_iter();
38        }
39
40        if n == 3 {
41            return vec![Triangle::new([
42                self.vertices[0],
43                self.vertices[1],
44                self.vertices[2],
45            ])]
46            .into_iter();
47        }
48
49        // 計算用に全て ECEF に変換
50        let ecef_points: Vec<Ecef> = self.vertices.iter().map(|&c| c.into()).collect();
51
52        // 投影軸の決定
53        let (u_axis, v_axis) = get_projection_axes(&ecef_points);
54
55        // 2D投影
56        let points_2d: Vec<(f64, f64)> = ecef_points
57            .iter()
58            .map(|p| p.project_2d(u_axis, v_axis))
59            .collect();
60
61        // 回転方向の検知
62        let area = signed_area(&points_2d);
63        let area_sign = if area > 0.0 { 1.0 } else { -1.0 };
64
65        // 耳切りループ
66        let mut indices: Vec<usize> = (0..n).collect();
67        let mut result = Vec::with_capacity(n - 2);
68        let mut count = 0;
69        let max_iters = n * n;
70
71        while indices.len() > 3 && count < max_iters {
72            let mut ear_found = false;
73            let len = indices.len();
74
75            for i in 0..len {
76                let prev_idx = indices[(i + len - 1) % len];
77                let curr_idx = indices[i];
78                let next_idx = indices[(i + 1) % len];
79
80                if is_ear(
81                    prev_idx, curr_idx, next_idx, &indices, &points_2d, area_sign,
82                ) {
83                    result.push(Triangle::new([
84                        self.vertices[prev_idx],
85                        self.vertices[curr_idx],
86                        self.vertices[next_idx],
87                    ]));
88                    indices.remove(i);
89                    ear_found = true;
90                    break;
91                }
92            }
93            if !ear_found {
94                break;
95            }
96            count += 1;
97        }
98
99        if indices.len() == 3 {
100            result.push(Triangle::new([
101                self.vertices[indices[0]],
102                self.vertices[indices[1]],
103                self.vertices[indices[2]],
104            ]));
105        }
106
107        result.into_iter()
108    }
109}
110
111/// Newell's Method による法線概算と投影軸の選択
112fn get_projection_axes(pts: &[Ecef]) -> (usize, usize) {
113    let mut nx = 0.0;
114    let mut ny = 0.0;
115    let mut nz = 0.0;
116
117    let len = pts.len();
118
119    for i in 0..len {
120        let curr = pts[i];
121        let next = pts[(i + 1) % len];
122
123        nx += (curr.y() - next.y()) * (curr.z() + next.z());
124        ny += (curr.z() - next.z()) * (curr.x() + next.x());
125        nz += (curr.x() - next.x()) * (curr.y() + next.y());
126    }
127
128    let ax = nx.abs();
129    let ay = ny.abs();
130    let az = nz.abs();
131
132    if ax >= ay && ax >= az {
133        (1, 2)
134    } else if ay >= ax && ay >= az {
135        (0, 2)
136    } else {
137        (0, 1)
138    }
139}
140
141fn signed_area(pts: &[(f64, f64)]) -> f64 {
142    let mut area = 0.0;
143    for i in 0..pts.len() {
144        let j = (i + 1) % pts.len();
145        area += pts[i].0 * pts[j].1;
146        area -= pts[j].0 * pts[i].1;
147    }
148    area / 2.0
149}
150
151fn is_ear(
152    p: usize,
153    c: usize,
154    n: usize,
155    indices: &[usize],
156    pts: &[(f64, f64)],
157    area_sign: f64,
158) -> bool {
159    let a = pts[p];
160    let b = pts[c];
161    let c_pt = pts[n];
162
163    // 2D外積 (z成分)
164    let cross = (b.0 - a.0) * (c_pt.1 - b.1) - (b.1 - a.1) * (c_pt.0 - b.0);
165
166    // 凹判定
167    if cross * area_sign <= -1e-10 {
168        return false;
169    }
170
171    // 包含判定
172    for &idx in indices {
173        if idx == p || idx == c || idx == n {
174            continue;
175        }
176        if is_point_in_triangle(pts[idx], a, b, c_pt) {
177            return false;
178        }
179    }
180    true
181}
182
183fn is_point_in_triangle(p: (f64, f64), a: (f64, f64), b: (f64, f64), c: (f64, f64)) -> bool {
184    let area2 = 0.5 * (-b.1 * c.0 + a.1 * (-b.0 + c.0) + a.0 * (b.1 - c.1) + b.0 * c.1);
185    if area2.abs() < 1e-12 {
186        return false;
187    }
188    let s = 1.0 / (2.0 * area2) * (a.1 * c.0 - a.0 * c.1 + (c.1 - a.1) * p.0 + (a.0 - c.0) * p.1);
189    let t = 1.0 / (2.0 * area2) * (a.0 * b.1 - a.1 * b.0 + (a.1 - b.1) * p.0 + (b.0 - a.0) * p.1);
190    s > 0.0 && t > 0.0 && (1.0 - s - t) > 0.0
191}
kasane_logic/geometry/shapes/polygon/geometry_relation.rs

kasane_logic/geometry/shapes/polygon/
geometry_relation.rs
