proof_engine/geometry/
geodesic.rs

1//! Geodesic distance computation on surfaces — shortest path along a curved surface.
2
3use glam::Vec3;
4use std::collections::{BinaryHeap, HashMap};
5use std::cmp::Ordering;
6use super::GeoMesh;
7
8/// Geodesic distance solver using Dijkstra on mesh edges.
9pub struct GeodesicSolver;
10
11#[derive(Debug, Clone)]
12struct DijkNode { vertex: u32, dist: f32 }
13impl PartialEq for DijkNode { fn eq(&self, other: &Self) -> bool { self.dist == other.dist } }
14impl Eq for DijkNode {}
15impl PartialOrd for DijkNode { fn partial_cmp(&self, other: &Self) -> Option<Ordering> { other.dist.partial_cmp(&self.dist) } }
16impl Ord for DijkNode { fn cmp(&self, other: &Self) -> Ordering { self.partial_cmp(other).unwrap_or(Ordering::Equal) } }
17
18impl GeodesicSolver {
19    /// Compute geodesic distances from a source vertex to all other vertices.
20    pub fn distances(mesh: &GeoMesh, source: u32) -> Vec<f32> {
21        let n = mesh.vertices.len();
22        let adj = Self::build_adjacency(mesh);
23        let mut dist = vec![f32::MAX; n];
24        let mut heap = BinaryHeap::new();
25
26        dist[source as usize] = 0.0;
27        heap.push(DijkNode { vertex: source, dist: 0.0 });
28
29        while let Some(DijkNode { vertex, dist: d }) = heap.pop() {
30            if d > dist[vertex as usize] { continue; }
31            if let Some(neighbors) = adj.get(&vertex) {
32                for &(nb, edge_len) in neighbors {
33                    let new_dist = d + edge_len;
34                    if new_dist < dist[nb as usize] {
35                        dist[nb as usize] = new_dist;
36                        heap.push(DijkNode { vertex: nb, dist: new_dist });
37                    }
38                }
39            }
40        }
41        dist
42    }
43
44    /// Compute shortest path between two vertices (sequence of vertex indices).
45    pub fn shortest_path(mesh: &GeoMesh, from: u32, to: u32) -> Vec<u32> {
46        let n = mesh.vertices.len();
47        let adj = Self::build_adjacency(mesh);
48        let mut dist = vec![f32::MAX; n];
49        let mut prev = vec![u32::MAX; n];
50        let mut heap = BinaryHeap::new();
51
52        dist[from as usize] = 0.0;
53        heap.push(DijkNode { vertex: from, dist: 0.0 });
54
55        while let Some(DijkNode { vertex, dist: d }) = heap.pop() {
56            if vertex == to { break; }
57            if d > dist[vertex as usize] { continue; }
58            if let Some(neighbors) = adj.get(&vertex) {
59                for &(nb, edge_len) in neighbors {
60                    let new_dist = d + edge_len;
61                    if new_dist < dist[nb as usize] {
62                        dist[nb as usize] = new_dist;
63                        prev[nb as usize] = vertex;
64                        heap.push(DijkNode { vertex: nb, dist: new_dist });
65                    }
66                }
67            }
68        }
69
70        // Reconstruct path
71        let mut path = Vec::new();
72        let mut current = to;
73        while current != u32::MAX && current != from {
74            path.push(current);
75            current = prev[current as usize];
76        }
77        if current == from { path.push(from); }
78        path.reverse();
79        path
80    }
81
82    fn build_adjacency(mesh: &GeoMesh) -> HashMap<u32, Vec<(u32, f32)>> {
83        let mut adj: HashMap<u32, Vec<(u32, f32)>> = HashMap::new();
84        for tri in &mesh.triangles {
85            let edges = [(tri.a, tri.b), (tri.b, tri.c), (tri.c, tri.a)];
86            for (a, b) in edges {
87                let d = (mesh.vertices[a as usize] - mesh.vertices[b as usize]).length();
88                adj.entry(a).or_default().push((b, d));
89                adj.entry(b).or_default().push((a, d));
90            }
91        }
92        adj
93    }
94}
95
96#[cfg(test)]
97mod tests {
98    use super::*;
99    use glam::Vec2;
100
101    fn simple_quad_mesh() -> GeoMesh {
102        let mut mesh = GeoMesh::new();
103        mesh.add_vertex(Vec3::new(0.0, 0.0, 0.0), Vec3::Y, Vec2::ZERO);
104        mesh.add_vertex(Vec3::new(1.0, 0.0, 0.0), Vec3::Y, Vec2::ZERO);
105        mesh.add_vertex(Vec3::new(1.0, 1.0, 0.0), Vec3::Y, Vec2::ZERO);
106        mesh.add_vertex(Vec3::new(0.0, 1.0, 0.0), Vec3::Y, Vec2::ZERO);
107        mesh.add_triangle(0, 1, 2);
108        mesh.add_triangle(0, 2, 3);
109        mesh
110    }
111
112    #[test]
113    fn distance_to_self_is_zero() {
114        let mesh = simple_quad_mesh();
115        let dists = GeodesicSolver::distances(&mesh, 0);
116        assert_eq!(dists[0], 0.0);
117    }
118
119    #[test]
120    fn distance_to_neighbor() {
121        let mesh = simple_quad_mesh();
122        let dists = GeodesicSolver::distances(&mesh, 0);
123        assert!((dists[1] - 1.0).abs() < 0.01); // direct edge
124    }
125
126    #[test]
127    fn shortest_path_basic() {
128        let mesh = simple_quad_mesh();
129        let path = GeodesicSolver::shortest_path(&mesh, 0, 2);
130        assert!(!path.is_empty());
131        assert_eq!(*path.first().unwrap(), 0);
132        assert_eq!(*path.last().unwrap(), 2);
133    }
134}
proof_engine/geometry/geodesic.rs

proof_engine/geometry/
geodesic.rs
