1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
//! Hermite data: where the surface crosses each cell edge, **and which way it
//! faces there**.
//!
//! Marching Cubes and Surface Nets need only the crossing positions. Dual
//! contouring needs the normals too, and that difference is the entire reason it
//! can reproduce a sharp corner: a position says the surface passes through
//! here, a position *and* a normal says it passes through here *going that way*,
//! and three of those intersect at a point that a corner actually is.
//!
//! # Source
//!
//! Ju, Losasso, Schaefer & Warren, *"Dual Contouring of Hermite Data"*, SIGGRAPH
//! 2002, `10.1145/566570.566586`. §2.2 defines the input as "the pairs
//! `pᵢ, nᵢ` which "correspond to the intersections (and unit normals) of the
//! contour with the edges of the cube" — which is exactly
//! [`HermiteCrossing`].
//!
//! # Why the position comes from the same interpolation as everywhere else
//!
//! The crossing is placed by the same linear interpolation the crate's shared
//! `cube::edge_crossing` gives Marching Cubes and Surface Nets, rather than by a
//! bisection or Newton refinement. That is deliberate and it is
//! about measurement, not cost: if dual contouring used *better* crossings than
//! Surface Nets, then E-104's side-by-side would be comparing two changes at
//! once — different crossings *and* different vertex placement — and could not
//! attribute the corner it recovers to either. Same crossings in, so the only
//! variable left is what the algorithm does with them.
//!
//! A refined root is a legitimate accuracy improvement later, but it belongs to
//! every extractor at once or to none.
use crate::cube::{EDGE_CORNERS, EDGE_COUNT, corner_offset, edge_crossing, is_inside};
use crate::equivariant::sum_equivariant;
use crate::vec3;
use crate::{Real, Sdf};
/// One edge crossing: where the surface is, and which way it faces.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct HermiteCrossing<R: Real> {
/// World position of the crossing.
pub position: [R; 3],
/// Unit surface normal there, pointing away from the solid.
///
/// Taken from the field's own gradient rather than estimated from geometry.
/// For an exact signed distance field this is the true normal; for a field
/// that is not a distance function it is still the direction of steepest
/// increase, which is what the vertex solve wants.
pub normal: [R; 3],
}
/// The crossings on one cell's twelve edges.
///
/// Up to twelve, in the crate's own edge numbering, with a bitmask saying which
/// are present. Fixed-size rather than a `Vec`, because a dual contourer builds
/// one of these per surface cell, and allocating per cell is exactly the cost
/// the reusable-buffer rule exists to avoid.
#[derive(Clone, Copy, Debug)]
pub struct HermiteCell<R: Real> {
crossings: [HermiteCrossing<R>; EDGE_COUNT],
/// Bit `e` set when edge `e` is cut.
mask: u16,
}
impl<R: Real> HermiteCell<R> {
/// Sample the crossings of one cell.
///
/// `corner_values` are the eight corner samples in this crate's corner
/// order, `cell_origin` is the world position of corner 0, and `cell_size`
/// is the edge length. The caller supplies the corner values because it has
/// already sampled them — re-sampling here would double the field
/// evaluations, which dominate extraction.
///
/// Costs one [`Sdf::gradient`] per crossing and nothing else.
pub fn from_corners<S: Sdf<Scalar = R>>(
sdf: &S,
corner_values: &[R; 8],
cell_origin: [R; 3],
cell_size: R,
) -> Self {
let empty = HermiteCrossing {
position: [R::ZERO; 3],
normal: [R::ZERO; 3],
};
let mut crossings = [empty; EDGE_COUNT];
let mut mask = 0u16;
for (edge, [lo, hi]) in EDGE_CORNERS.iter().copied().enumerate() {
let (a, b) = (corner_values[lo as usize], corner_values[hi as usize]);
if is_inside(a) == is_inside(b) {
continue;
}
let t = edge_crossing(a, b);
let (lo_offset, hi_offset) = (corner_offset(lo), corner_offset(hi));
let mut position = [R::ZERO; 3];
for (axis, slot) in position.iter_mut().enumerate() {
let from = R::from_f64(f64::from(lo_offset[axis]));
let to = R::from_f64(f64::from(hi_offset[axis]));
*slot = cell_origin[axis] + cell_size * (from + (to - from) * t);
}
let gradient = sdf.gradient(position);
let length = vec3::length(gradient);
// A zero gradient at a crossing means the field is degenerate there.
// It cannot happen for an exact distance field, where |grad| is 1.
debug_assert!(length > R::ZERO, "zero gradient at an edge crossing");
let normal = vec3::scale(gradient, length.recip());
crossings[edge] = HermiteCrossing { position, normal };
mask |= 1 << edge;
}
Self { crossings, mask }
}
/// Build a cell from explicit `(edge, crossing)` pairs.
///
/// Internal, and it exists for one reason: the dual-contouring vertex rule
/// is a piece of linear algebra whose interesting cases — three orthogonal
/// planes, a rank-1 flat region, a lattice rotation applied to both input
/// and expected output — are stated as *planes*, not as a field and a grid.
/// Reaching them through [`from_corners`](Self::from_corners) would mean
/// hand-designing an SDF for each, and the test would then be exercising the
/// SDF as much as the solve.
///
/// Edges outside `0..EDGE_COUNT` are ignored rather than rejected: this
/// takes a fixed test fixture, not user input.
#[cfg(test)]
pub(crate) fn from_crossings(pairs: &[(u8, HermiteCrossing<R>)]) -> Self {
let mut cell = Self {
crossings: [HermiteCrossing {
position: [R::ZERO; 3],
normal: [R::ZERO; 3],
}; EDGE_COUNT],
mask: 0,
};
for &(edge, crossing) in pairs {
if (edge as usize) < EDGE_COUNT {
cell.crossings[edge as usize] = crossing;
cell.mask |= 1 << edge;
}
}
cell
}
/// The same crossings, restricted to the edges in `edges`.
///
/// Bit `e` of `edges` keeps edge `e`; edges this cell does not cut stay
/// uncut. Everything derived from the cell — [`len`](Self::len),
/// [`iter`](Self::iter), [`centroid`](Self::centroid), and the vertex solve —
/// then sees only that subset.
///
/// This is what turns one cell into several vertices. A cell the surface
/// passes through twice has its cut edges partitioned into **surface
/// components**, and solving each component's edges separately puts a vertex
/// on each sheet instead of one vertex shared between them. See
/// [`ManifoldDualContouring`](crate::manifold_dual_contouring::ManifoldDualContouring),
/// which is the only caller that needs it and the reason it exists.
#[must_use]
pub fn restricted(&self, edges: u16) -> Self {
Self {
crossings: self.crossings,
mask: self.mask & edges,
}
}
/// How many edges the surface crosses.
#[must_use]
pub fn len(&self) -> usize {
self.mask.count_ones() as usize
}
/// `true` when the surface misses this cell entirely.
#[must_use]
pub fn is_empty(&self) -> bool {
self.mask == 0
}
/// Whether edge `edge` is cut.
#[must_use]
pub fn contains(&self, edge: u8) -> bool {
edge < EDGE_COUNT as u8 && self.mask & (1 << edge) != 0
}
/// The crossing on one edge, if there is one.
#[must_use]
pub fn get(&self, edge: u8) -> Option<&HermiteCrossing<R>> {
if self.contains(edge) {
Some(&self.crossings[edge as usize])
} else {
None
}
}
/// Every crossing, in edge order.
///
/// Edge order rather than insertion order, so the sequence is a function of
/// the cell's corner signs alone — which is what keeps a vertex solve built
/// on it deterministic.
pub fn iter(&self) -> impl Iterator<Item = &HermiteCrossing<R>> {
(0..EDGE_COUNT as u8).filter_map(|edge| self.get(edge))
}
/// The average of the crossing positions.
///
/// This is precisely the vertex Surface Nets would place, and dual
/// contouring's starting point: the QEF is minimised subject to staying near
/// it, and it is the answer when the normals carry no directional
/// information at all. Returns `None` for an empty cell.
///
/// # Why the sum is not taken in edge order
///
/// A lattice rotation permutes the edge labels, and IEEE addition is not
/// associative, so a running sum over [`iter`](Self::iter) would make this
/// point — and every vertex the QEF places relative to it — a function of
/// how the cell happens to sit on the grid rather than of the cell. Each
/// component is reduced over a slot per edge instead, smallest magnitude
/// first (A-016; the mechanism and the measurement are in
/// the crate-internal `equivariant` module).
#[must_use]
pub fn centroid(&self) -> Option<[R; 3]> {
if self.is_empty() {
return None;
}
let mut axes = [[R::ZERO; EDGE_COUNT]; 3];
for edge in 0..EDGE_COUNT {
let Some(crossing) = self.get(edge as u8) else {
continue;
};
for (slot, value) in axes.iter_mut().zip(crossing.position) {
slot[edge] = value;
}
}
let inverse = R::from_f64(self.len() as f64).recip();
let sum = axes.map(sum_equivariant);
Some([sum[0] * inverse, sum[1] * inverse, sum[2] * inverse])
}
}
#[cfg(test)]
mod tests;