ogeom_mesh/editable/smooth.rs
1//! Smoothing an editable mesh by projecting each vertex onto a local
2//! quadric, keeping creases sharp where asked.
3//!
4//! # A pass
5//!
6//! Each vertex is moved along its neighbourhood's normal onto the local
7//! quadric fitted about it, the one [`curvature`](crate::curvature()) reads:
8//! a weighted least-squares height field over the vertex's tangent plane,
9//! fitted to the [`NEIGHBOURS`](crate::curvature::NEIGHBOURS) vertices
10//! nearest it within the smoothing's `size`. Every vertex of a pass is
11//! moved from the positions the pass started from, so the result does not
12//! depend on the order the vertices are taken in, and is the same at any
13//! thread count. Smoothing makes [`PASSES`] passes.
14//!
15//! # Creases
16//!
17//! With [`Smooth::keep_edges`] an angle, a neighbour counts towards a
18//! vertex's quadric only if its normal turns from the vertex's own by no
19//! more than that angle. A vertex beside a crease then fits the faces on
20//! its own side, and stays on them. A vertex on the crease, whose normal
21//! lies between the faces', is fitted to the vertices along the crease, or
22//! where they fix no quadric, as at a corner, stays where it is.
23//!
24//! # How far a vertex moves
25//!
26//! No vertex ends farther than [`Smooth::max_move`] from where it started:
27//! a move past it is shortened, along its line, until the distance as
28//! computed is no more than `max_move`.
29//!
30//! With [`Smooth::region`], only the vertices in it move, each by its
31//! pass's move scaled by its distance to the nearest vertex outside the
32//! region over `size`, up to the whole move from one `size` in. The region
33//! thus blends into the rest of the mesh rather than stepping off it.
34//!
35//! # What stays
36//!
37//! Only positions change. A vertex with no face, and every vertex
38//! [`Editable::non_manifold`] reports, the copies of a point several fans
39//! meet and the ends of edges not cleanly joined, stays where it is, so a
40//! point's copies stay together. How faces are joined across an edge of
41//! three or more of them depends on where their far corners stand, so
42//! those edges are joined again afterwards, by the rule in the [module
43//! documentation](super).
44
45use super::{Editable, NONE, Scope, index};
46use crate::curvature::{Cloud, check_size, in_parallel, vertex_normals};
47use crate::point_index::PointIndex;
48use ogeom_core::{OgeomResult, ogeom_bail};
49use ogeom_math::Point;
50
51/// How many passes a smoothing makes.
52pub const PASSES: usize = 2;
53
54/// How a mesh is smoothed.
55#[derive(Debug, Clone, Copy, PartialEq)]
56pub struct Smooth<'a> {
57 /// How far about each vertex its neighbours are taken from.
58 pub size: f64,
59 /// The farthest a vertex may end from where it started.
60 pub max_move: f64,
61 /// The angle, in radians, beyond which a neighbour's normal turns too
62 /// far from a vertex's for it to count; `None` to count every
63 /// neighbour.
64 pub keep_edges: Option<f64>,
65 /// Which vertices may move, one entry per vertex number; `None` for
66 /// all.
67 pub region: Option<&'a [bool]>,
68}
69
70/// What a smoothing moved.
71#[derive(Debug, Clone, Copy, PartialEq)]
72pub struct SmoothReport {
73 /// The farthest a vertex ended from where it started.
74 pub max: f64,
75 /// The root mean square of how far each vertex with a face ended from
76 /// where it started.
77 pub rms: f64,
78 /// How many vertices `max_move` held back.
79 pub limited: usize,
80}
81
82impl Editable {
83 /// Smooth the mesh, moving its vertices onto local quadrics fitted to
84 /// their neighbours, as the [module documentation](crate::editable::smooth)
85 /// describes. Vertex and face numbers stay; only an edge of three or
86 /// more faces may be joined anew, as its faces now stand.
87 ///
88 /// # Errors
89 ///
90 /// [`OgeomError::Numeric`](ogeom_core::OgeomError::Numeric) if `size` is
91 /// not positive and finite, `max_move` is negative or not a number, or
92 /// `keep_edges` is not an angle from zero to `π`;
93 /// [`OgeomError::Range`](ogeom_core::OgeomError::Range) if `region` does
94 /// not have one entry per vertex number.
95 pub fn smooth(&mut self, smooth: &Smooth<'_>) -> OgeomResult<SmoothReport> {
96 let Smooth {
97 size,
98 max_move,
99 keep_edges,
100 region,
101 } = *smooth;
102 check_size(size)?;
103 if max_move.is_nan() || max_move < 0.0 {
104 ogeom_bail!(Numeric, "a largest move of {max_move} is not a distance");
105 }
106 if let Some(angle) = keep_edges
107 && !(0.0..=core::f64::consts::PI).contains(&angle)
108 {
109 ogeom_bail!(Numeric, "a crease angle of {angle} is not from 0 to π");
110 }
111 let n = self.positions.len();
112 if let Some(region) = region
113 && region.len() != n
114 {
115 ogeom_bail!(
116 Range,
117 "a region of {} entries for {n} vertices",
118 region.len()
119 );
120 }
121
122 let with_face: Vec<u32> = (0..n).filter(|&v| self.out[v] != NONE).map(index).collect();
123 let blend = self.blend(&with_face, size, region)?;
124 let movers: Vec<u32> = with_face
125 .iter()
126 .copied()
127 .filter(|&v| blend[v as usize] > 0.0)
128 .collect();
129 let faces: Vec<[u32; 3]> = self
130 .corners
131 .iter()
132 .zip(&self.live)
133 .filter(|&(_, &live)| live)
134 .map(|(&t, _)| t)
135 .collect();
136 let cos = keep_edges.map(f64::cos);
137 let start = self.positions.clone();
138 let mut limited = vec![false; n];
139 for _ in 0..PASSES {
140 let normals = vertex_normals(&self.positions, faces.iter().copied());
141 let members: Vec<u32> = with_face
142 .iter()
143 .copied()
144 .filter(|&v| normals[v as usize].is_some())
145 .collect();
146 let cloud = Cloud::build(&self.positions, &normals, members)?;
147 let targets = in_parallel(&movers, |v| {
148 let p = self.positions[v as usize];
149 let near = cloud.near(p, size);
150 let local = match cos {
151 Some(cos) => {
152 let own = normals[v as usize]?;
153 cloud.fit(p, &near, |w| {
154 normals[w as usize].is_some_and(|m| m.dot(own) >= cos)
155 })
156 }
157 None => cloud.fit(p, &near, |_| true),
158 }?;
159 Some(local.projected())
160 });
161 drop(cloud);
162 for (&v, target) in movers.iter().zip(targets) {
163 let Some(target) = target else { continue };
164 let k = v as usize;
165 let p = self.positions[k];
166 let moved = p + (target - p) * blend[k];
167 let (q, held) = capped(start[k], moved, max_move);
168 self.positions[k] = q;
169 limited[k] = held;
170 }
171 }
172 self.rejoin_crowded_edges();
173
174 let (mut max, mut sum) = (0.0_f64, 0.0);
175 for &v in &with_face {
176 let d = (self.positions[v as usize] - start[v as usize]).magnitude();
177 max = max.max(d);
178 sum += d * d;
179 }
180 let rms = if with_face.is_empty() {
181 0.0
182 } else {
183 (sum / f64::from(index(with_face.len()))).sqrt()
184 };
185 Ok(SmoothReport {
186 max,
187 rms,
188 limited: limited.iter().filter(|&&l| l).count(),
189 })
190 }
191
192 /// The share of its move each vertex makes: zero for a vertex with no
193 /// face, one [`Editable::non_manifold`] reports, or one outside the
194 /// region; otherwise its distance to the nearest vertex with a face
195 /// outside the region over `size`, up to one.
196 fn blend(
197 &self,
198 with_face: &[u32],
199 size: f64,
200 region: Option<&[bool]>,
201 ) -> OgeomResult<Vec<f64>> {
202 let mut blend = vec![0.0; self.positions.len()];
203 let outside: Vec<Point> = with_face
204 .iter()
205 .filter(|&&v| region.is_some_and(|r| !r[v as usize]))
206 .map(|&v| self.positions[v as usize])
207 .collect();
208 let outside = PointIndex::build(&outside)?;
209 for &v in with_face {
210 let k = v as usize;
211 if region.is_none_or(|r| r[k]) {
212 blend[k] = outside
213 .nearest(self.positions[k], size)
214 .map_or(1.0, |near| near.distance / size);
215 }
216 }
217 for v in self.non_manifold() {
218 blend[v as usize] = 0.0;
219 }
220 Ok(blend)
221 }
222
223 /// Join anew the faces on every edge that three or more of them use.
224 pub(super) fn rejoin_crowded_edges(&mut self) {
225 let mut edges: Vec<(u64, u32)> = Vec::new();
226 for (f, &live) in self.live.iter().enumerate() {
227 if live {
228 for h in index(3 * f)..index(3 * f + 3) {
229 edges.push((self.edge_key(h), h));
230 }
231 }
232 }
233 edges.sort_unstable();
234 let mut points: Vec<u32> = Vec::new();
235 for run in edges.chunk_by(|a, b| a.0 == b.0) {
236 if run.len() > 2 {
237 let h = run[0].1;
238 points.push(self.point[self.tail(h) as usize]);
239 points.push(self.point[self.head(h) as usize]);
240 }
241 }
242 if points.is_empty() {
243 return;
244 }
245 points.sort_unstable();
246 points.dedup();
247 let around = self.about(&points);
248 self.relink(&Scope::Points(points), &around);
249 }
250}
251
252/// `to`, or where a move from `from` towards it is shortened to end no
253/// farther than `most` from `from` as computed, and whether it was.
254pub(super) fn capped(from: Point, to: Point, most: f64) -> (Point, bool) {
255 let d = to - from;
256 let length = d.magnitude();
257 if length <= most {
258 return (to, false);
259 }
260 let mut s = most / length;
261 loop {
262 let q = from + d * s;
263 if (q - from).magnitude() <= most {
264 return (q, true);
265 }
266 s = s.next_down();
267 }
268}