graphic design tool
git clone https://git.lucas.co/cce-designer.git
src/hull.rs (5.4K)
1 //! The convex hull of a point cloud, as a closed triangle mesh — the `hull`
2 //! node, and what hou-control's `developer_embryo` does with its scattered
3 //! points (its `shrinkwrap`).
4 //!
5 //! The incremental algorithm rather than quickhull: a starting tetrahedron
6 //! from the extreme points, then each remaining point in turn either lies
7 //! inside the hull so far or sees some of its faces — those are removed, and
8 //! the ring of edges where visible meets hidden (the horizon) is fanned to
9 //! the new point. Points within a tolerance of a face are treated as inside,
10 //! which is what the HDA's Remove Inline Points does: a hull with a thousand
11 //! coplanar slivers is not a better hull. It is O(points × faces), which for
12 //! a thousand scattered points is nothing; a hull of a million would want
13 //! the conflict lists.
14
15 use crate::detail::Detail;
16 use glam::Vec3;
17
18 /// The convex hull of `points`, as a closed triangle mesh over only the
19 /// points that lie on it — `None` when the points do not span a volume.
20 pub fn convex_hull(points: &[Vec3]) -> Option<Detail> {
21 let faces = hull_faces(points)?;
22 let mut remap = vec![u32::MAX; points.len()];
23 let mut d = Detail::new();
24 for f in &faces {
25 let mut ids = [0u32; 3];
26 for (k, &pi) in f.iter().enumerate() {
27 if remap[pi] == u32::MAX {
28 remap[pi] = d.add_point(points[pi]);
29 }
30 ids[k] = remap[pi];
31 }
32 d.add_prim(&ids);
33 }
34 Some(d)
35 }
36
37 /// The hull's faces as index triples into `points`, wound outward.
38 fn hull_faces(points: &[Vec3]) -> Option<Vec<[usize; 3]>> {
39 if points.len() < 4 {
40 return None;
41 }
42 let (lo, hi) = points
43 .iter()
44 .fold((Vec3::splat(f32::MAX), Vec3::splat(f32::MIN)), |(lo, hi), &p| (lo.min(p), hi.max(p)));
45 let diag = (hi - lo).length();
46 if !(diag > 0.0) {
47 return None;
48 }
49 // "On the face" and "no volume" are both judged against the cloud's own
50 // size, so the hull of a millimetre embryo and of a metre one build the
51 // same way.
52 let eps = diag * 1e-5;
53
54 // The starting tetrahedron: the two points furthest apart along an axis,
55 // the point furthest from that line, the point furthest from that plane.
56 let mut ext = [0usize; 6];
57 for (i, p) in points.iter().enumerate() {
58 for a in 0..3 {
59 if p[a] < points[ext[a]][a] {
60 ext[a] = i;
61 }
62 if p[a] > points[ext[a + 3]][a] {
63 ext[a + 3] = i;
64 }
65 }
66 }
67 let (mut i0, mut i1, mut best) = (0, 0, -1.0);
68 for &a in &ext {
69 for &b in &ext {
70 let d = (points[a] - points[b]).length();
71 if d > best {
72 best = d;
73 i0 = a;
74 i1 = b;
75 }
76 }
77 }
78 if best <= eps {
79 return None;
80 }
81 let dir = (points[i1] - points[i0]).normalize();
82 let (mut i2, mut best) = (0, -1.0);
83 for (i, &p) in points.iter().enumerate() {
84 let off = p - points[i0];
85 let d = (off - dir * off.dot(dir)).length();
86 if d > best {
87 best = d;
88 i2 = i;
89 }
90 }
91 if best <= eps {
92 return None;
93 }
94 let n = (points[i1] - points[i0]).cross(points[i2] - points[i0]).normalize();
95 let (mut i3, mut best) = (0, -1.0);
96 for (i, &p) in points.iter().enumerate() {
97 let d = (p - points[i0]).dot(n).abs();
98 if d > best {
99 best = d;
100 i3 = i;
101 }
102 }
103 if best <= eps {
104 return None;
105 }
106
107 // Wind the four faces so every normal points away from the centroid.
108 let centroid = (points[i0] + points[i1] + points[i2] + points[i3]) / 4.0;
109 let mut faces: Vec<[usize; 3]> = Vec::new();
110 for f in [[i0, i1, i2], [i0, i1, i3], [i0, i2, i3], [i1, i2, i3]] {
111 let n = face_normal(points, f);
112 if (points[f[0]] - centroid).dot(n) < 0.0 {
113 faces.push([f[0], f[2], f[1]]);
114 } else {
115 faces.push(f);
116 }
117 }
118
119 let mut edges: Vec<(usize, usize)> = Vec::new();
120 for (pi, &p) in points.iter().enumerate() {
121 if pi == i0 || pi == i1 || pi == i2 || pi == i3 {
122 continue;
123 }
124 // The faces this point looks at from outside.
125 let visible: Vec<bool> = faces
126 .iter()
127 .map(|f| (p - points[f[0]]).dot(face_normal(points, *f)) > eps)
128 .collect();
129 if !visible.iter().any(|&v| v) {
130 continue;
131 }
132 // The horizon: directed edges of visible faces whose reverse is not
133 // an edge of a visible face. The winding of the visible face gives
134 // the new face's winding for free.
135 edges.clear();
136 for (f, &v) in faces.iter().zip(&visible) {
137 if v {
138 edges.push((f[0], f[1]));
139 edges.push((f[1], f[2]));
140 edges.push((f[2], f[0]));
141 }
142 }
143 let horizon: Vec<(usize, usize)> =
144 edges.iter().copied().filter(|&(a, b)| !edges.contains(&(b, a))).collect();
145 let mut kept: Vec<[usize; 3]> =
146 faces.iter().zip(&visible).filter(|(_, &v)| !v).map(|(f, _)| *f).collect();
147 for (a, b) in horizon {
148 kept.push([a, b, pi]);
149 }
150 faces = kept;
151 }
152 Some(faces)
153 }
154
155 fn face_normal(points: &[Vec3], f: [usize; 3]) -> Vec3 {
156 (points[f[1]] - points[f[0]]).cross(points[f[2]] - points[f[0]]).normalize_or_zero()
157 }