axiolid_reference/
primitive.rs

1//! Tessellation of CSG primitives into closed solids.
2//!
3//! Every variant is analytic, so each is emitted directly. A block has no
4//! curvature to sample; a sphere's is exactly known. The chord budget still
5//! decides the radial segment count so a primitive and a swept face of the
6//! same radius agree on what a tolerance means.
7
8use axiolid_contracts::{GeomError, GeomResult};
9use axiolid_core::{Point3, Scalar, Tolerance};
10use axiolid_mesh::TriMesh;
11use axiolid_primitive::Primitive;
12
13/// Radial segments needed for `radius` under a chord budget.
14///
15/// Same sagitta rule the curve flattener uses: r(1 - cos(pi/n)) <= tol.
16/// Clamped low so a degenerate tolerance cannot ask for an unbounded mesh.
17fn segments(radius: Scalar, tolerance: Scalar) -> usize {
18    if !(radius.is_finite() && radius > 0.0 && tolerance.is_finite() && tolerance > 0.0) {
19        return 3;
20    }
21    let ratio = 1.0 - (tolerance / radius).min(1.0);
22    let n = (core::f64::consts::PI / ratio.acos().max(1e-9)).ceil();
23    (n as usize).clamp(3, 4096)
24}
25
26/// Tessellate one CSG primitive into a closed, outward-wound solid.
27///
28/// Outward winding is not decoration: `axiolid-mesh-boolean-boolmesh` and the clash
29/// containment test both read signed volume, and an inverted primitive
30/// silently produces negative volume and wrong verdicts.
31pub fn tessellate_primitive(primitive: &Primitive, tolerance: Tolerance) -> GeomResult<TriMesh> {
32    let tol = tolerance.linear();
33    match primitive {
34        Primitive::Block { x, y, z } => block(*x, *y, *z),
35        Primitive::Sphere { radius } => sphere(*radius, tol),
36        Primitive::Cylinder { radius, height } => cylinder(*radius, *height, tol),
37        Primitive::Cone { radius, height } => cone(*radius, *height, tol),
38        Primitive::Pyramid { x, y, height } => pyramid(*x, *y, *height),
39        Primitive::Torus {
40            major_radius,
41            minor_radius,
42        } => torus(*major_radius, *minor_radius, tol),
43        Primitive::Wedge {
44            x,
45            y,
46            height,
47            top_x_min,
48            top_x_max,
49            top_y_min,
50            top_y_max,
51        } => wedge(
52            [*x, *y, *height],
53            [*top_x_min, *top_x_max],
54            [*top_y_min, *top_y_max],
55        ),
56        // The enum is non_exhaustive: a new family is unsupported, never
57        // silently approximated by the nearest one.
58        _ => Err(GeomError::Unsupported {
59            backend: crate::ScalarBoolean::ID,
60            operation: axiolid_contracts::Operation::Tessellation,
61        }),
62    }
63}
64
65/// Validate a positive finite extent, naming the offender.
66fn positive(value: Scalar, what: &str) -> GeomResult<Scalar> {
67    if !value.is_finite() || value <= 0.0 {
68        return Err(GeomError::InvalidInput(format!(
69            "{what} must be positive and finite, got {value}"
70        )));
71    }
72    Ok(value)
73}
74
75/// Axis-aligned block centred on the local origin.
76fn block(x: Scalar, y: Scalar, z: Scalar) -> GeomResult<TriMesh> {
77    let (hx, hy, hz) = (
78        positive(x, "block x")? / 2.0,
79        positive(y, "block y")? / 2.0,
80        positive(z, "block z")? / 2.0,
81    );
82    let p = vec![
83        Point3::new(-hx, -hy, -hz),
84        Point3::new(hx, -hy, -hz),
85        Point3::new(hx, hy, -hz),
86        Point3::new(-hx, hy, -hz),
87        Point3::new(-hx, -hy, hz),
88        Point3::new(hx, -hy, hz),
89        Point3::new(hx, hy, hz),
90        Point3::new(-hx, hy, hz),
91    ];
92    // Outward winding, verified by the positive-volume test rather than by
93    // reading the index list.
94    let i = vec![
95        0, 2, 1, 0, 3, 2, 4, 5, 6, 4, 6, 7, 0, 1, 5, 0, 5, 4, 1, 2, 6, 1, 6, 5, 2, 3, 7, 2, 7, 6,
96        3, 0, 4, 3, 4, 7,
97    ];
98    Ok(TriMesh::new(p, i))
99}
100
101/// Rectangular pyramid: base on z = 0, apex on +z.
102fn pyramid(x: Scalar, y: Scalar, height: Scalar) -> GeomResult<TriMesh> {
103    let (hx, hy) = (
104        positive(x, "pyramid x")? / 2.0,
105        positive(y, "pyramid y")? / 2.0,
106    );
107    let h = positive(height, "pyramid height")?;
108    let p = vec![
109        Point3::new(-hx, -hy, 0.0),
110        Point3::new(hx, -hy, 0.0),
111        Point3::new(hx, hy, 0.0),
112        Point3::new(-hx, hy, 0.0),
113        Point3::new(0.0, 0.0, h),
114    ];
115    let i = vec![0, 2, 1, 0, 3, 2, 0, 1, 4, 1, 2, 4, 2, 3, 4, 3, 0, 4];
116    Ok(TriMesh::new(p, i))
117}
118
119/// Wedge: base `[0, x] x [0, y]` at z = 0, top `[x0, x1] x [y0, y1]` at
120/// z = `height`.
121///
122/// Every face is planar: the x sides join edges parallel to y, the y sides
123/// edges parallel to x. A top collapsed to a segment or a point shares
124/// vertices, so the faces that lose area are dropped and the rest meet at
125/// the shared vertices rather than along zero-length edges.
126fn wedge(
127    [x, y, height]: [Scalar; 3],
128    [x0, x1]: [Scalar; 2],
129    [y0, y1]: [Scalar; 2],
130) -> GeomResult<TriMesh> {
131    let (x, y, h) = (
132        positive(x, "wedge x")?,
133        positive(y, "wedge y")?,
134        positive(height, "wedge height")?,
135    );
136    for (value, what) in [
137        (x0, "wedge top x min"),
138        (x1, "wedge top x max"),
139        (y0, "wedge top y min"),
140        (y1, "wedge top y max"),
141    ] {
142        if !value.is_finite() {
143            return Err(GeomError::InvalidInput(format!(
144                "{what} must be finite, got {value}"
145            )));
146        }
147    }
148    if x0 > x1 || y0 > y1 {
149        return Err(GeomError::InvalidInput(format!(
150            "wedge top must have min <= max, got x {x0}..{x1}, y {y0}..{y1}"
151        )));
152    }
153    let corners = [
154        Point3::new(0.0, 0.0, 0.0),
155        Point3::new(x, 0.0, 0.0),
156        Point3::new(x, y, 0.0),
157        Point3::new(0.0, y, 0.0),
158        Point3::new(x0, y0, h),
159        Point3::new(x1, y0, h),
160        Point3::new(x1, y1, h),
161        Point3::new(x0, y1, h),
162    ];
163    // Outward faces over the corners above: bottom, top, then the sides at
164    // y = 0, x = x, y = y and x = 0.
165    const FACES: [[usize; 4]; 6] = [
166        [0, 3, 2, 1],
167        [4, 5, 6, 7],
168        [0, 1, 5, 4],
169        [1, 2, 6, 5],
170        [2, 3, 7, 6],
171        [3, 0, 4, 7],
172    ];
173    let mut p: Vec<Point3> = Vec::with_capacity(8);
174    let index: Vec<u32> = corners
175        .iter()
176        .map(|&c| match p.iter().position(|&q| q == c) {
177            Some(k) => k as u32,
178            None => {
179                p.push(c);
180                (p.len() - 1) as u32
181            }
182        })
183        .collect();
184    let mut i = Vec::with_capacity(36);
185    for face in FACES {
186        let mut ring: Vec<u32> = Vec::with_capacity(4);
187        for corner in face {
188            let v = index[corner];
189            if ring.last() != Some(&v) {
190                ring.push(v);
191            }
192        }
193        if ring.len() > 1 && ring.first() == ring.last() {
194            ring.pop();
195        }
196        if ring.len() < 3 {
197            continue;
198        }
199        // Each face is convex, so a fan triangulates it.
200        for k in 1..ring.len() - 1 {
201            i.extend([ring[0], ring[k], ring[k + 1]]);
202        }
203    }
204    Ok(TriMesh::new(p, i))
205}
206
207/// Ring torus about +z, tube centre circle of radius `major` in z = 0.
208///
209/// A grid of `n` steps round the axis by `m` round the tube, each sized by
210/// the same chord rule as the other curved primitives: `n` for the outer
211/// equator, the largest circle round the axis. Each grid cell is a planar
212/// trapezoid (its two edges round the axis are parallel chords), split in
213/// two.
214fn torus(major: Scalar, minor: Scalar, tol: Scalar) -> GeomResult<TriMesh> {
215    let big = positive(major, "torus major radius")?;
216    let r = positive(minor, "torus minor radius")?;
217    if r >= big {
218        let kind = if r == big {
219            "a horn torus (minor radius equal to major)"
220        } else {
221            "a spindle torus (minor radius above major)"
222        };
223        return Err(GeomError::InvalidInput(format!(
224            "torus minor radius {r} must be below major radius {big}: {kind} \
225             does not bound a two-manifold solid"
226        )));
227    }
228    let n = segments(big + r, tol);
229    let m = segments(r, tol);
230    let mut p = Vec::with_capacity(n * m);
231    for i in 0..n {
232        let theta = core::f64::consts::TAU * (i as Scalar) / (n as Scalar);
233        for j in 0..m {
234            let phi = core::f64::consts::TAU * (j as Scalar) / (m as Scalar);
235            let rho = big + r * phi.cos();
236            p.push(Point3::new(
237                rho * theta.cos(),
238                rho * theta.sin(),
239                r * phi.sin(),
240            ));
241        }
242    }
243    let at = |i: usize, j: usize| ((i % n) * m + (j % m)) as u32;
244    let mut idx = Vec::with_capacity(n * m * 6);
245    for i in 0..n {
246        for j in 0..m {
247            // Round the axis, then round the tube: outward, since the
248            // axis tangent crossed with the tube tangent points away from
249            // the tube's centre.
250            let (a, b, c, d) = (at(i, j), at(i + 1, j), at(i + 1, j + 1), at(i, j + 1));
251            idx.extend([a, b, c]);
252            idx.extend([a, c, d]);
253        }
254    }
255    Ok(TriMesh::new(p, idx))
256}
257
258/// A ring of `n` points at `radius`, height `z`.
259fn ring(radius: Scalar, z: Scalar, n: usize) -> Vec<Point3> {
260    (0..n)
261        .map(|k| {
262            let a = core::f64::consts::TAU * (k as Scalar) / (n as Scalar);
263            Point3::new(radius * a.cos(), radius * a.sin(), z)
264        })
265        .collect()
266}
267
268/// Cylinder along +z, base on z = 0.
269fn cylinder(radius: Scalar, height: Scalar, tol: Scalar) -> GeomResult<TriMesh> {
270    let r = positive(radius, "cylinder radius")?;
271    let h = positive(height, "cylinder height")?;
272    let n = segments(r, tol);
273    let mut p = ring(r, 0.0, n);
274    p.extend(ring(r, h, n));
275    p.push(Point3::new(0.0, 0.0, 0.0));
276    p.push(Point3::new(0.0, 0.0, h));
277    let (bc, tc) = (2 * n, 2 * n + 1);
278    let mut i = Vec::with_capacity(n * 12);
279    for k in 0..n {
280        let (a, b) = (k, (k + 1) % n);
281        // Side quad, then the two caps. The base fan is wound opposite to
282        // the top so both face away from the enclosed volume.
283        i.extend([a as u32, b as u32, (b + n) as u32]);
284        i.extend([a as u32, (b + n) as u32, (a + n) as u32]);
285        i.extend([bc as u32, b as u32, a as u32]);
286        i.extend([tc as u32, (a + n) as u32, (b + n) as u32]);
287    }
288    Ok(TriMesh::new(p, i))
289}
290
291/// Cone along +z: base ring on z = 0, apex at height.
292fn cone(radius: Scalar, height: Scalar, tol: Scalar) -> GeomResult<TriMesh> {
293    let r = positive(radius, "cone radius")?;
294    let h = positive(height, "cone height")?;
295    let n = segments(r, tol);
296    let mut p = ring(r, 0.0, n);
297    p.push(Point3::new(0.0, 0.0, 0.0));
298    p.push(Point3::new(0.0, 0.0, h));
299    let (base, apex) = (n, n + 1);
300    let mut i = Vec::with_capacity(n * 6);
301    for k in 0..n {
302        let (a, b) = (k as u32, ((k + 1) % n) as u32);
303        i.extend([base as u32, b, a]);
304        i.extend([a, b, apex as u32]);
305    }
306    Ok(TriMesh::new(p, i))
307}
308
309/// Sphere centred on the local origin, as a UV mesh.
310fn sphere(radius: Scalar, tol: Scalar) -> GeomResult<TriMesh> {
311    let r = positive(radius, "sphere radius")?;
312    let n = segments(r, tol);
313    // Half as many stacks as segments: the polar direction spans PI, not TAU,
314    // so equal counts would oversample it by 2x for the same chord error.
315    let stacks = (n / 2).max(2);
316    let mut p = Vec::with_capacity((stacks + 1) * n);
317    for i in 0..=stacks {
318        let v = core::f64::consts::PI * (i as Scalar) / (stacks as Scalar);
319        for j in 0..n {
320            let u = core::f64::consts::TAU * (j as Scalar) / (n as Scalar);
321            p.push(Point3::new(
322                r * v.sin() * u.cos(),
323                r * v.sin() * u.sin(),
324                r * v.cos(),
325            ));
326        }
327    }
328    let mut idx = Vec::with_capacity(stacks * n * 6);
329    for i in 0..stacks {
330        for j in 0..n {
331            let jn = (j + 1) % n;
332            // A pole row's n entries are the same point, so they must map
333            // to ONE index or no triangle around the pole shares an edge.
334            let row = |r: usize, k: usize| -> u32 {
335                if r == 0 || r == stacks {
336                    (r * n) as u32
337                } else {
338                    (r * n + k) as u32
339                }
340            };
341            let a = row(i, j);
342            let b = row(i, jn);
343            let c = row(i + 1, j);
344            let d = row(i + 1, jn);
345            // Pole rows collapse to a point; skip the degenerate half.
346            if i > 0 {
347                idx.extend([a, c, b]);
348            }
349            if i + 1 < stacks {
350                idx.extend([b, c, d]);
351            }
352        }
353    }
354    Ok(TriMesh::new(p, idx))
355}
axiolid_reference/primitive.rs

axiolid_reference/
primitive.rs
