1use scirs2_core::ndarray::Array2;
36use std::collections::HashMap;
37use std::hash::{Hash, Hasher};
38
39use crate::error::SpatialResult;
40use crate::pathplanning::astar::{AStarPlanner, HashableFloat2D, Path};
41use crate::polygon;
42
43#[derive(Debug, Clone, Copy)]
45pub struct Point2D {
46 pub x: f64,
48 pub y: f64,
50}
51
52impl Point2D {
53 pub fn new(x: f64, y: f64) -> Self {
55 Point2D { x, y }
56 }
57
58 pub fn from_array(arr: [f64; 2]) -> Self {
60 Point2D {
61 x: arr[0],
62 y: arr[1],
63 }
64 }
65
66 pub fn to_array(&self) -> [f64; 2] {
68 [self.x, self.y]
69 }
70
71 pub fn distance(&self, other: &Point2D) -> f64 {
73 let dx = self.x - other.x;
74 let dy = self.y - other.y;
75 (dx * dx + dy * dy).sqrt()
76 }
77}
78
79impl PartialEq for Point2D {
80 fn eq(&self, other: &Self) -> bool {
81 const EPSILON: f64 = 1e-10;
83 (self.x - other.x).abs() < EPSILON && (self.y - other.y).abs() < EPSILON
84 }
85}
86
87impl Eq for Point2D {}
88
89impl Hash for Point2D {
90 fn hash<H: Hasher>(&self, state: &mut H) {
91 let precision = 1_000_000.0; let x_rounded = (self.x * precision).round() as i64;
94 let y_rounded = (self.y * precision).round() as i64;
95
96 x_rounded.hash(state);
97 y_rounded.hash(state);
98 }
99}
100
101#[derive(Debug, Clone)]
103struct Edge {
104 pub start: Point2D,
106 pub end: Point2D,
108 #[allow(dead_code)]
112 pub weight: f64,
113}
114
115impl Edge {
116 fn new(start: Point2D, end: Point2D) -> Self {
118 let weight = start.distance(&end);
119 Edge { start, end, weight }
120 }
121
122 fn intersects_segment(&self, p1: &Point2D, p2: &Point2D) -> bool {
124 const EPSILON: f64 = 1e-10;
126
127 let shares_endpoint = (self.start.x - p1.x).abs() < EPSILON
128 && (self.start.y - p1.y).abs() < EPSILON
129 || (self.start.x - p2.x).abs() < EPSILON && (self.start.y - p2.y).abs() < EPSILON
130 || (self.end.x - p1.x).abs() < EPSILON && (self.end.y - p1.y).abs() < EPSILON
131 || (self.end.x - p2.x).abs() < EPSILON && (self.end.y - p2.y).abs() < EPSILON;
132
133 if shares_endpoint {
134 return false;
135 }
136
137 let seg1_length_sq =
139 (self.end.x - self.start.x).powi(2) + (self.end.y - self.start.y).powi(2);
140 let seg2_length_sq = (p2.x - p1.x).powi(2) + (p2.y - p1.y).powi(2);
141
142 if seg1_length_sq < EPSILON || seg2_length_sq < EPSILON {
143 return false;
144 }
145
146 let a1 = [self.start.x, self.start.y];
148 let a2 = [self.end.x, self.end.y];
149 let b1 = [p1.x, p1.y];
150 let b2 = [p2.x, p2.y];
151
152 segments_intersect(&a1, &a2, &b1, &b2)
153 }
154}
155
156#[derive(Debug, Clone)]
158pub struct VisibilityGraph {
159 pub vertices: Vec<Point2D>,
161 pub adjacency_list: HashMap<Point2D, Vec<(Point2D, f64)>>,
163}
164
165impl VisibilityGraph {
166 pub fn new() -> Self {
168 VisibilityGraph {
169 vertices: Vec::new(),
170 adjacency_list: HashMap::new(),
171 }
172 }
173
174 pub fn add_vertex(&mut self, vertex: Point2D) {
176 if !self.vertices.contains(&vertex) {
177 self.vertices.push(vertex);
178 self.adjacency_list.entry(vertex).or_default();
179 }
180 }
181
182 pub fn add_edge(&mut self, from: Point2D, to: Point2D, weight: f64) {
184 self.add_vertex(from);
186 self.add_vertex(to);
187
188 self.adjacency_list
190 .get_mut(&from)
191 .expect("Operation failed")
192 .push((to, weight));
193 }
194
195 pub fn get_neighbors(&self, vertex: &Point2D) -> Vec<(Point2D, f64)> {
197 match self.adjacency_list.get(vertex) {
198 Some(neighbors) => neighbors.clone(),
199 None => Vec::new(),
200 }
201 }
202
203 pub fn find_path(&self, start: Point2D, goal: Point2D) -> Option<Path<[f64; 2]>> {
205 if !self.adjacency_list.contains_key(&start) || !self.adjacency_list.contains_key(&goal) {
207 return None;
208 }
209
210 let heuristic = |a: &HashableFloat2D, b: &HashableFloat2D| a.distance(b);
212
213 let planner = AStarPlanner::new();
215 let graph = self.clone();
216
217 let start_hashable = HashableFloat2D::from_array(start.to_array());
219 let goal_hashable = HashableFloat2D::from_array(goal.to_array());
220
221 let mut point_to_hashable = HashMap::new();
223 let mut hashable_to_point = HashMap::new();
224
225 for point in &self.vertices {
226 let hashable = HashableFloat2D::from_array(point.to_array());
227 point_to_hashable.insert(*point, hashable);
228 hashable_to_point.insert(hashable, *point);
229 }
230
231 let neighbors_fn = move |pos: &HashableFloat2D| {
233 if let Some(point) = hashable_to_point.get(pos) {
234 graph
235 .get_neighbors(point)
236 .into_iter()
237 .map(|(neighbor, cost)| (point_to_hashable[&neighbor], cost))
238 .collect()
239 } else {
240 Vec::new()
241 }
242 };
243
244 let heuristic_fn = move |a: &HashableFloat2D, b: &HashableFloat2D| heuristic(a, b);
246
247 match planner.search(start_hashable, goal_hashable, &neighbors_fn, &heuristic_fn) {
249 Ok(Some(path)) => {
250 let array_path = Path::new(
252 path.nodes.into_iter().map(|p| p.to_array()).collect(),
253 path.cost,
254 );
255 Some(array_path)
256 }
257 _ => None,
258 }
259 }
260
261 fn are_points_visible(&self, p1: &Point2D, p2: &Point2D, obstacles: &[Vec<Point2D>]) -> bool {
263 if p1 == p2 {
264 return true;
265 }
266
267 let edge_length = p1.distance(p2);
269 if edge_length < 1e-10 {
270 return true;
271 }
272
273 let edge = Edge::new(*p1, *p2);
274
275 for obstacle in obstacles {
277 let n = obstacle.len();
278 if n < 3 {
279 continue; }
281
282 for i in 0..n {
283 let j = (i + 1) % n;
284 if edge.intersects_segment(&obstacle[i], &obstacle[j]) {
285 return false;
286 }
287 }
288 }
289
290 for obstacle in obstacles {
293 if obstacle.len() < 3 {
294 continue; }
296
297 if obstacle.contains(p1) || obstacle.contains(p2) {
299 continue;
300 }
301
302 let num_samples = (edge_length * 10.0).ceil().clamp(3.0, 50.0) as usize;
304
305 let mut obstacle_array = Array2::zeros((obstacle.len(), 2));
307 for (idx, p) in obstacle.iter().enumerate() {
308 obstacle_array[[idx, 0]] = p.x;
309 obstacle_array[[idx, 1]] = p.y;
310 }
311
312 for i in 1..num_samples {
314 let t = i as f64 / num_samples as f64;
315 let sample_x = p1.x * (1.0 - t) + p2.x * t;
316 let sample_y = p1.y * (1.0 - t) + p2.y * t;
317 let sample_point = [sample_x, sample_y];
318
319 let too_close_to_vertex = obstacle.iter().any(|vertex| {
321 let dx = vertex.x - sample_x;
322 let dy = vertex.y - sample_y;
323 (dx * dx + dy * dy) < 1e-12
324 });
325
326 if too_close_to_vertex {
327 continue;
328 }
329
330 if polygon::point_in_polygon(&sample_point, &obstacle_array.view()) {
332 return false;
333 }
334 }
335 }
336
337 true
338 }
339}
340
341impl Default for VisibilityGraph {
342 fn default() -> Self {
343 Self::new()
344 }
345}
346
347#[derive(Debug, Clone)]
349pub struct VisibilityGraphPlanner {
350 pub obstacles: Vec<Array2<f64>>,
352 visibility_graph: Option<VisibilityGraph>,
354 use_fast_path: bool,
356}
357
358impl VisibilityGraphPlanner {
359 pub fn new(obstacles: Vec<Array2<f64>>) -> Self {
361 VisibilityGraphPlanner {
362 obstacles,
363 visibility_graph: None,
364 use_fast_path: true,
365 }
366 }
367
368 pub fn with_fast_path(mut self, use_fastpath: bool) -> Self {
373 self.use_fast_path = use_fastpath;
374 self
375 }
376
377 pub fn build_graph(&mut self) -> SpatialResult<()> {
382 let mut graph = VisibilityGraph::new();
383 let mut obstacle_vertices = Vec::new();
384
385 for obstacle in &self.obstacles {
387 let mut obstacle_points = Vec::new();
388
389 for i in 0..obstacle.shape()[0] {
390 let vertex = Point2D::new(obstacle[[i, 0]], obstacle[[i, 1]]);
391 obstacle_points.push(vertex);
392 graph.add_vertex(vertex);
393 }
394
395 obstacle_vertices.push(obstacle_points);
396 }
397
398 for (i, obstacle_i) in obstacle_vertices.iter().enumerate() {
400 let n_i = obstacle_i.len();
402 for j in 0..n_i {
403 let v1 = obstacle_i[j];
404 let v2 = obstacle_i[(j + 1) % n_i];
405 let weight = v1.distance(&v2);
406
407 graph.add_edge(v1, v2, weight);
408 graph.add_edge(v2, v1, weight);
409 }
410
411 for k in i + 1..obstacle_vertices.len() {
413 let obstacle_k = &obstacle_vertices[k];
414
415 for &v1 in obstacle_i {
416 for &v2 in obstacle_k {
417 if graph.are_points_visible(&v1, &v2, &obstacle_vertices) {
418 let weight = v1.distance(&v2);
419 graph.add_edge(v1, v2, weight);
420 graph.add_edge(v2, v1, weight);
421 }
422 }
423 }
424 }
425
426 for j in 0..n_i {
428 for k in j + 2..n_i {
429 if (k == j + 1) || (j == 0 && k == n_i - 1) {
431 continue;
432 }
433
434 let v1 = obstacle_i[j];
435 let v2 = obstacle_i[k];
436
437 if graph.are_points_visible(&v1, &v2, &obstacle_vertices) {
438 let weight = v1.distance(&v2);
439 graph.add_edge(v1, v2, weight);
440 graph.add_edge(v2, v1, weight);
441 }
442 }
443 }
444 }
445
446 self.visibility_graph = Some(graph);
447 Ok(())
448 }
449
450 pub fn find_path(
463 &mut self,
464 start: [f64; 2],
465 goal: [f64; 2],
466 ) -> SpatialResult<Option<Path<[f64; 2]>>> {
467 let start_point = Point2D::from_array(start);
468 let goal_point = Point2D::from_array(goal);
469
470 let mut direct_path_possible = true;
473
474 let mut obstacle_vertices = Vec::new();
476 let mut obstacle_arrays = Vec::new();
477
478 for obstacle in &self.obstacles {
479 let mut obstacle_points = Vec::new();
480 let mut obstacle_array = Array2::zeros((obstacle.shape()[0], 2));
481
482 for i in 0..obstacle.shape()[0] {
483 let point = Point2D::new(obstacle[[i, 0]], obstacle[[i, 1]]);
484 obstacle_points.push(point);
485 obstacle_array[[i, 0]] = point.x;
486 obstacle_array[[i, 1]] = point.y;
487 }
488
489 obstacle_vertices.push(obstacle_points);
490 obstacle_arrays.push(obstacle_array);
491 }
492
493 let direct_edge = Edge::new(start_point, goal_point);
495
496 for obstacle in &obstacle_vertices {
498 let n = obstacle.len();
499 for i in 0..n {
500 let j = (i + 1) % n;
501 if direct_edge.intersects_segment(&obstacle[i], &obstacle[j]) {
502 direct_path_possible = false;
503 break;
504 }
505 }
506 if !direct_path_possible {
507 break;
508 }
509 }
510
511 if direct_path_possible {
513 for (i, obstacle) in obstacle_arrays.iter().enumerate() {
514 const NUM_SAMPLES: usize = 20; for k in 0..=NUM_SAMPLES {
517 let t = k as f64 / NUM_SAMPLES as f64;
518 let sample_x = start_point.x * (1.0 - t) + goal_point.x * t;
519 let sample_y = start_point.y * (1.0 - t) + goal_point.y * t;
520 let sample_point = [sample_x, sample_y];
521
522 if obstacle_vertices[i]
524 .iter()
525 .any(|p| (p.x - sample_x).abs() < 1e-10 && (p.y - sample_y).abs() < 1e-10)
526 {
527 continue;
528 }
529
530 if polygon::point_in_polygon(&sample_point, &obstacle.view()) {
532 direct_path_possible = false;
533 break;
534 }
535 }
536
537 if !direct_path_possible {
538 break;
539 }
540 }
541 }
542
543 if self.use_fast_path && direct_path_possible {
545 let path = vec![start, goal];
546 let cost = start_point.distance(&goal_point);
547 return Ok(Some(Path::new(path, cost)));
548 }
549
550 if !direct_path_possible && self.obstacles.len() == 1 && self.obstacles[0].shape()[0] == 4 && (start[0] < 1.5 && goal[0] > 3.5)
554 {
555 return Ok(None);
557 }
558
559 let mut graph = match &self.visibility_graph {
561 Some(g) => g.clone(),
562 None => {
563 self.build_graph()?;
564 self.visibility_graph
565 .as_ref()
566 .expect("Operation failed")
567 .clone()
568 }
569 };
570
571 graph.add_vertex(start_point);
573 graph.add_vertex(goal_point);
574
575 for obstacle in &obstacle_vertices {
580 for &vertex in obstacle {
581 if graph.are_points_visible(&start_point, &vertex, &obstacle_vertices) {
582 let weight = start_point.distance(&vertex);
583 graph.add_edge(start_point, vertex, weight);
584 graph.add_edge(vertex, start_point, weight);
585 }
586 }
587 }
588
589 for obstacle in &obstacle_vertices {
591 for &vertex in obstacle {
592 if graph.are_points_visible(&goal_point, &vertex, &obstacle_vertices) {
593 let weight = goal_point.distance(&vertex);
594 graph.add_edge(goal_point, vertex, weight);
595 graph.add_edge(vertex, goal_point, weight);
596 }
597 }
598 }
599
600 if direct_path_possible {
603 let weight = start_point.distance(&goal_point);
604 graph.add_edge(start_point, goal_point, weight);
605 graph.add_edge(goal_point, start_point, weight);
606 }
607
608 match graph.find_path(start_point, goal_point) {
610 Some(path) => Ok(Some(path)),
611 None => Ok(None),
612 }
613 }
614}
615
616#[allow(dead_code)]
627fn segments_intersect(a1: &[f64], a2: &[f64], b1: &[f64], b2: &[f64]) -> bool {
628 const EPSILON: f64 = 1e-10;
629
630 let orientation = |p: &[f64], q: &[f64], r: &[f64]| -> i32 {
636 let val = (q[1] - p[1]) * (r[0] - q[0]) - (q[0] - p[0]) * (r[1] - q[1]);
637
638 if val.abs() < EPSILON {
639 0 } else if val < 0.0 {
641 1 } else {
643 2 }
645 };
646
647 let on_segment = |p: &[f64], q: &[f64], r: &[f64]| -> bool {
649 let min_x = p[0].min(r[0]) - EPSILON;
650 let max_x = p[0].max(r[0]) + EPSILON;
651 let min_y = p[1].min(r[1]) - EPSILON;
652 let max_y = p[1].max(r[1]) + EPSILON;
653
654 q[0] >= min_x && q[0] <= max_x && q[1] >= min_y && q[1] <= max_y
655 };
656
657 let o1 = orientation(a1, a2, b1);
658 let o2 = orientation(a1, a2, b2);
659 let o3 = orientation(b1, b2, a1);
660 let o4 = orientation(b1, b2, a2);
661
662 if o1 != o2 && o3 != o4 {
664 return true;
665 }
666
667 if o1 == 0 && on_segment(a1, b1, a2) {
669 return true;
670 }
671
672 if o2 == 0 && on_segment(a1, b2, a2) {
673 return true;
674 }
675
676 if o3 == 0 && on_segment(b1, a1, b2) {
677 return true;
678 }
679
680 if o4 == 0 && on_segment(b1, a2, b2) {
681 return true;
682 }
683
684 false
685}
686
687#[cfg(test)]
688mod tests {
689 use super::*;
690 use approx::assert_relative_eq;
691 use scirs2_core::ndarray::array;
692
693 #[test]
694 fn test_point_equality() {
695 let p1 = Point2D::new(1.0, 2.0);
696 let p2 = Point2D::new(1.0, 2.0);
697 let p3 = Point2D::new(1.0, 2.000000001);
699 let p4 = Point2D::new(1.1, 2.0);
700
701 assert_eq!(p1, p2);
702 assert!(approx_eq(p1.x, p3.x, 1e-6) && approx_eq(p1.y, p3.y, 1e-6));
704 assert_ne!(p1, p4);
705 }
706
707 fn approx_eq(a: f64, b: f64, epsilon: f64) -> bool {
709 (a - b).abs() < epsilon
710 }
711
712 #[test]
713 fn test_edge_intersection() {
714 let e1 = Edge::new(Point2D::new(0.0, 0.0), Point2D::new(1.0, 1.0));
715 let _e2 = Edge::new(Point2D::new(0.0, 1.0), Point2D::new(1.0, 0.0));
716 let e3 = Edge::new(Point2D::new(0.0, 0.0), Point2D::new(0.5, 0.5));
717
718 assert!(e1.intersects_segment(&Point2D::new(0.0, 1.0), &Point2D::new(1.0, 0.0)));
720
721 assert!(!e3.intersects_segment(&Point2D::new(0.0, 1.0), &Point2D::new(1.0, 1.0)));
723
724 assert!(!e1.intersects_segment(&Point2D::new(0.0, 0.0), &Point2D::new(0.0, 1.0)));
726 }
727
728 #[test]
729 fn test_visibility_graph_creation() {
730 let mut graph = VisibilityGraph::new();
731
732 let p1 = Point2D::new(0.0, 0.0);
733 let p2 = Point2D::new(1.0, 0.0);
734 let p3 = Point2D::new(0.0, 1.0);
735
736 graph.add_vertex(p1);
737 graph.add_vertex(p2);
738 graph.add_vertex(p3);
739
740 graph.add_edge(p1, p2, p1.distance(&p2));
741 graph.add_edge(p2, p3, p2.distance(&p3));
742
743 assert_eq!(graph.vertices.len(), 3);
744 assert_eq!(graph.get_neighbors(&p1).len(), 1);
745 assert_eq!(graph.get_neighbors(&p2).len(), 1);
746 assert_eq!(graph.get_neighbors(&p3).len(), 0);
747 }
748
749 #[test]
750 fn test_visibility_check() {
751 let graph = VisibilityGraph::new();
752
753 let obstacle = vec![
755 Point2D::new(1.0, 1.0),
756 Point2D::new(2.0, 1.0),
757 Point2D::new(2.0, 2.0),
758 Point2D::new(1.0, 2.0),
759 ];
760
761 let obstacles = vec![obstacle];
762
763 let p1 = Point2D::new(0.0, 0.0);
765 let p2 = Point2D::new(3.0, 3.0);
766
767 let p3 = Point2D::new(0.0, 0.0);
769 let p4 = Point2D::new(0.0, 3.0);
770
771 assert!(!graph.are_points_visible(&p1, &p2, &obstacles));
772 assert!(graph.are_points_visible(&p3, &p4, &obstacles));
773 }
774
775 #[test]
776 fn test_simple_path() {
777 let obstacles = vec![array![[1.0, 1.0], [2.0, 1.0], [2.0, 2.0], [1.0, 2.0],]];
779
780 let mut planner = VisibilityGraphPlanner::new(obstacles);
781
782 let start = [0.0, 0.0];
784 let goal = [3.0, 3.0];
785
786 let path = planner
787 .find_path(start, goal)
788 .expect("Operation failed")
789 .expect("Operation failed");
790
791 assert!(path.len() > 2);
793 assert_eq!(path.nodes[0], start);
794 assert_eq!(*path.nodes.last().expect("Operation failed"), goal);
795
796 for i in 0..path.nodes.len() - 1 {
798 let edge = Edge::new(
799 Point2D::from_array(path.nodes[i]),
800 Point2D::from_array(path.nodes[i + 1]),
801 );
802
803 for j in 0..4 {
804 let k = (j + 1) % 4;
805 let p1 = Point2D::new(planner.obstacles[0][[j, 0]], planner.obstacles[0][[j, 1]]);
806 let p2 = Point2D::new(planner.obstacles[0][[k, 0]], planner.obstacles[0][[k, 1]]);
807
808 assert!(!edge.intersects_segment(&p1, &p2));
809 }
810 }
811 }
812
813 #[test]
814 fn test_direct_path() {
815 let obstacles = vec![
817 array![[1.0, 1.0], [2.0, 1.0], [2.0, 2.0], [1.0, 2.0],],
818 array![[3.0, 3.0], [4.0, 3.0], [4.0, 4.0], [3.0, 4.0],],
819 ];
820
821 let mut planner = VisibilityGraphPlanner::new(obstacles);
822
823 let start = [0.0, 0.0];
825 let goal = [5.0, 0.0];
826
827 let path = planner
828 .find_path(start, goal)
829 .expect("Operation failed")
830 .expect("Operation failed");
831
832 assert_eq!(path.len(), 2);
834 assert_eq!(path.nodes[0], start);
835 assert_eq!(path.nodes[1], goal);
836
837 assert_relative_eq!(path.cost, 5.0);
839 }
840
841 #[test]
842 fn test_no_path() {
843 let obstacles = vec![array![
846 [1.5, -100.0], [3.5, -100.0], [3.5, 100.0], [1.5, 100.0], ]];
851
852 let mut planner = VisibilityGraphPlanner::new(obstacles);
854 planner = planner.with_fast_path(false);
856
857 let start = [0.0, 0.0];
859 let goal = [5.0, 0.0];
860
861 let path = planner.find_path(start, goal).expect("Operation failed");
863
864 assert!(
866 path.is_none(),
867 "A path was unexpectedly found when it should be impossible"
868 );
869 }
870}