use super::*;
const MITER_TRANSVERSALITY_FLOOR: f64 = 0.1;
fn fold_angle() -> f64 {
2.0 * MITER_TRANSVERSALITY_FLOOR.acos()
}
pub(super) fn sweep_section_mitered(
section: &[NurbsCurve],
section_normal: Vec3,
path: &SweepPath,
twist: f64,
) -> Result<BrepSolid, String> {
let plan = plan_mitre(section, section_normal, path, twist)?;
let MitrePlan {
closed,
segments,
count,
ring_count,
rings,
reversed_winding,
start_frame: (x_axis, y_axis),
end_frame: (end_x, end_y),
bands,
closure,
..
} = plan;
let shift = closure.map_or(0, |closure| closure.curve_shift);
let landing = |segment: usize, corner: usize| {
let next_ring = (segment + 1) % ring_count;
if next_ring == 0 {
(0, (corner + shift) % count)
} else {
(next_ring, corner)
}
};
let vertex_id = |ring: usize, corner: usize| (ring * count + corner) as u64 + 1;
let ring_edge_id = |ring: usize, corner: usize| 1000 + (ring * count + corner) as u64;
let side_edge_id = |segment: usize, corner: usize| {
1000 + (ring_count * count + segment * count + corner) as u64
};
let mut vertices = Vec::with_capacity(ring_count * count);
let mut ring_points: Vec<Vec<Vec3>> = Vec::with_capacity(ring_count);
for (index, curves) in rings.iter().take(ring_count).enumerate() {
let mut points = Vec::with_capacity(count);
for curve in curves {
let [t0, _] = curve.domain()?;
points.push(curve.evaluate(t0)?);
}
for (corner, point) in points.iter().enumerate() {
vertices.push(VertexRecord {
id: vertex_id(index, corner),
point: *point,
});
}
ring_points.push(points);
}
let mut edges = Vec::with_capacity((ring_count + segments) * count);
for (index, curves) in rings.iter().take(ring_count).enumerate() {
for (corner, curve) in curves.iter().enumerate() {
let [t0, t1] = curve.domain()?;
edges.push(EdgeRecord {
id: ring_edge_id(index, corner),
curve: curve.clone(),
t0,
t1,
start_vertex_id: vertex_id(index, corner),
end_vertex_id: vertex_id(index, (corner + 1) % count),
degenerate: false,
name: None,
});
}
}
for segment in 0..segments {
for corner in 0..count {
let (next_ring, arriving) = landing(segment, corner);
let curve = match &bands {
None => make_line(
ring_points[segment][corner],
ring_points[next_ring][arriving],
)?,
Some(bands) => twisted_side(
&rings[segment][corner],
&rings[next_ring][arriving],
&bands[segment],
)?,
};
edges.push(EdgeRecord {
id: side_edge_id(segment, corner),
curve,
t0: 0.0,
t1: 1.0,
start_vertex_id: vertex_id(segment, corner),
end_vertex_id: vertex_id(next_ring, arriving),
degenerate: false,
name: None,
});
}
}
let mut next_id = 1_000_000_u64;
let mut faces: Vec<FaceRecord> = Vec::with_capacity(segments * count + 2);
for segment in 0..segments {
let mut walls = Vec::with_capacity(count);
for corner in 0..count {
let (next_ring, arriving) = landing(segment, corner);
let bottom = &rings[segment][corner];
let top = &rings[next_ring][arriving];
let [t0, t1] = bottom.domain()?;
let surface = match &bands {
None => ruled_between(bottom, top)?,
Some(bands) => twisted_wall(bottom, top, &bands[segment])?,
};
let coedges = vec![
CoedgeRecord {
id: next_id,
edge_id: ring_edge_id(segment, corner),
forward: true,
pcurve: parameter_line(t0, 0.0, t1, 0.0)?,
},
CoedgeRecord {
id: next_id + 1,
edge_id: side_edge_id(segment, (corner + 1) % count),
forward: true,
pcurve: parameter_line(t1, 0.0, t1, 1.0)?,
},
CoedgeRecord {
id: next_id + 2,
edge_id: ring_edge_id(next_ring, arriving),
forward: false,
pcurve: parameter_line(t1, 1.0, t0, 1.0)?,
},
CoedgeRecord {
id: next_id + 3,
edge_id: side_edge_id(segment, corner),
forward: false,
pcurve: parameter_line(t0, 1.0, t0, 0.0)?,
},
];
next_id += 4;
walls.push(FaceRecord {
id: next_id + 1,
surface,
same_sense: true,
loops: vec![LoopRecord {
id: next_id,
coedges,
}],
name: None,
});
next_id += 2;
}
if reversed_winding {
walls.reverse();
}
faces.extend(walls);
}
let cap = |index: usize,
ex: Vec3,
ey: Vec3,
forward: bool,
next_id: &mut u64|
-> Result<FaceRecord, String> {
let curves = &rings[index];
let anchor = ring_points[index][0];
let mut min_x = f64::INFINITY;
let mut min_y = f64::INFINITY;
let mut max_x = f64::NEG_INFINITY;
let mut max_y = f64::NEG_INFINITY;
for curve in curves {
let [t0, t1] = curve.domain()?;
for sample in 0..=16 {
let point = curve.evaluate(t0 + (t1 - t0) * sample as f64 / 16.0)?;
let delta = point.sub(anchor);
min_x = min_x.min(delta.dot(ex));
max_x = max_x.max(delta.dot(ex));
min_y = min_y.min(delta.dot(ey));
max_y = max_y.max(delta.dot(ey));
}
}
let padding = (max_x - min_x).max(max_y - min_y) * 0.05 + 1e-6;
let origin = anchor
.add(ex.scale(min_x - padding))
.add(ey.scale(min_y - padding));
let surface = make_plane(
origin,
ex,
ey,
max_x - min_x + 2.0 * padding,
max_y - min_y + 2.0 * padding,
)?;
let mut coedges = Vec::with_capacity(count);
if forward {
for (corner, curve) in curves.iter().enumerate() {
coedges.push(CoedgeRecord {
id: *next_id,
edge_id: ring_edge_id(index, corner),
forward: true,
pcurve: curve_to_plane_parameters(curve, origin, ex, ey)?,
});
*next_id += 1;
}
} else {
for corner in (0..count).rev() {
coedges.push(CoedgeRecord {
id: *next_id,
edge_id: ring_edge_id(index, corner),
forward: false,
pcurve: curve_to_plane_parameters(&curves[corner], origin, ex, ey)?.reversed()?,
});
*next_id += 1;
}
}
let loop_id = *next_id;
*next_id += 2;
Ok(FaceRecord {
id: loop_id + 1,
surface,
same_sense: forward,
loops: vec![LoopRecord {
id: loop_id,
coedges,
}],
name: None,
})
};
if !closed {
faces.push(cap(0, x_axis, y_axis, false, &mut next_id)?);
faces.push(cap(segments, end_x, end_y, true, &mut next_id)?);
}
let solid = BrepSolid {
id: next_id + 1,
vertices,
edges,
shells: vec![ShellRecord {
id: next_id,
faces,
}],
genus: i64::from(closed),
};
let issues = solid.validate();
if !issues.is_empty() {
return Err(format!(
"sweepSolid: the mitred {} produced invalid topology: {issues:?}",
if closed { "frame" } else { "sweep" }
));
}
crate::accept_sound(solid, "sweepSolid")
}
fn plan_mitre(
section: &[NurbsCurve],
section_normal: Vec3,
path: &SweepPath,
requested: f64,
) -> Result<MitrePlan, String> {
let tolerance = 1e-6;
let count = section.len();
if count < 2 {
return Err("sweepSolid: a mitred sweep needs a section of at least 2 curves".into());
}
if requested != 0.0 && !path.closed {
return Err(
"sweepSolid: a twist on a mitred run is rolled along the sides of a CLOSED frame; \
this path is open"
.into(),
);
}
if !path.cornered_polyline() {
return Err(
"sweepSolid: the mitre lane takes a path of STRAIGHT segments with at least one \
corner"
.into(),
);
}
let segments = path.len();
let closed = path.closed;
if closed && segments < 3 {
return Err(format!(
"sweepSolid: a CLOSED cornered path needs at least 3 segments to be a frame; this \
one has {segments}. Two straight segments that close on each other are the same \
line twice"
));
}
let ring_count = if closed { segments } else { segments + 1 };
let mut joints: Vec<Vec3> = Vec::with_capacity(segments + 1);
let mut dirs: Vec<Vec3> = Vec::with_capacity(segments);
let mut lengths: Vec<f64> = Vec::with_capacity(segments);
let [first_t0, _] = path.curves[0].domain()?;
joints.push(path.curves[0].evaluate(first_t0)?);
for index in 0..segments {
let curve = &path.curves[index];
let [t0, t1] = curve.domain()?;
let chord = curve.evaluate(t1)?.sub(curve.evaluate(t0)?);
let length = chord.length();
if length <= tolerance {
return Err(format!(
"sweepSolid: path segment '{}' has zero length",
path.name(index)
));
}
dirs.push(chord.scale(1.0 / length));
lengths.push(length);
joints.push(joints[index].add(chord));
}
let travel: f64 = lengths.iter().sum();
let mut plane_point: Vec<Vec3> = vec![Vec3::default(); segments + 1];
let mut plane_normal: Vec<Vec3> = vec![Vec3::default(); segments + 1];
let mut turns: Vec<f64> = vec![0.0; segments + 1];
for index in usize::from(!closed)..segments {
let before = (index + segments - 1) % segments;
let arriving = dirs[before];
let departing = dirs[index];
let bisector = arriving.add(departing);
let cos_half = bisector.length() * 0.5;
let turn = arriving.dot(departing).clamp(-1.0, 1.0).acos();
if cos_half < MITER_TRANSVERSALITY_FLOOR {
return Err(format!(
"sweepSolid: the joint between path segments '{}' and '{}' is a COLLAPSED bend — \
it turns {:.3}°, at or past the {:.3}° fold angle a mitre can carry. The mitre \
plane bisects the two directions, so at {:.3}° it makes only {:.2}° with the \
path and the section's footprint on it is stretched {:.1}× across the bend; at \
180° the path doubles back on itself and there is no bisector plane at all. \
Open the bend, or sweep the two legs separately",
path.name(before),
path.name(index),
turn.to_degrees(),
fold_angle().to_degrees(),
turn.to_degrees(),
(90.0 - turn.to_degrees() / 2.0),
1.0 / cos_half.max(f64::MIN_POSITIVE),
));
}
plane_normal[index] = bisector.scale(1.0 / (2.0 * cos_half));
plane_point[index] = joints[index];
turns[index] = turn;
}
if closed {
plane_point[segments] = plane_point[0];
plane_normal[segments] = plane_normal[0];
turns[segments] = turns[0];
}
if closed {
for first in 0..segments {
for second in (first + 2)..segments {
if first == 0 && second == segments - 1 {
continue;
}
let gap = segment_gap(
joints[first],
joints[first + 1],
joints[second],
joints[second + 1],
);
if gap <= tolerance * travel {
return Err(format!(
"sweepSolid: this closed path CROSSES ITSELF — segments '{}' and '{}' \
pass within {gap:.6e} of each other, and they do not share a joint. \
Every joint of it still has a mitre, so the frame would be built and \
would pass through itself where those two sides meet. Redraw the loop \
so its sides only meet at its corners",
path.name(first),
path.name(second),
));
}
}
}
}
let mut normal = section_normal
.normalized()
.map_err(|_| "sweepSolid: the section's plane normal is degenerate".to_string())?;
if normal.dot(dirs[0]) < 0.0 {
normal = normal.scale(-1.0);
}
if normal.dot(dirs[0]) < MITER_TRANSVERSALITY_FLOOR {
return Err(
"sweepSolid: the path is nearly parallel to the section plane, which sweeps a sliver \
rather than a solid"
.into(),
);
}
let x_axis = normal.perpendicular()?;
let y_axis = normal.cross(x_axis).normalized()?;
let section_origin = {
let [t0, _] = section[0].domain()?;
section[0].evaluate(t0)?
};
let mut ring: Vec<NurbsCurve> = section.to_vec();
let reversed_winding = profile_area(&ring, section_origin, x_axis, y_axis)? < 0.0;
if reversed_winding {
ring = ring
.iter()
.rev()
.map(NurbsCurve::reversed)
.collect::<Result<_, _>>()?;
}
if !closed {
plane_point[0] = section_origin;
plane_normal[0] = normal;
}
let mut end_x = x_axis;
let mut end_y = y_axis;
if !closed {
let mut end_normal = normal;
let mut anchor_offset = section_origin.sub(joints[0]);
for index in 1..segments {
let (from, to) = (dirs[index - 1], dirs[index]);
end_normal = rotate_minimal(end_normal, from, to);
end_x = rotate_minimal(end_x, from, to);
end_y = rotate_minimal(end_y, from, to);
anchor_offset = rotate_minimal(anchor_offset, from, to);
}
plane_normal[segments] = end_normal;
plane_point[segments] = joints[segments].add(anchor_offset);
}
let mut rings: Vec<Vec<NurbsCurve>> = Vec::with_capacity(segments + 1);
if closed {
let mut first = Vec::with_capacity(count);
for curve in &ring {
first.push(project_curve_along(
curve,
dirs[0],
plane_point[0],
plane_normal[0],
)?);
}
rings.push(first);
} else {
rings.push(ring);
}
for index in 0..segments {
let mut next = Vec::with_capacity(count);
for curve in &rings[index] {
next.push(project_curve_along(
curve,
dirs[index],
plane_point[index + 1],
plane_normal[index + 1],
)?);
}
rings.push(next);
}
let room_floor = tolerance * travel;
let mut closure = None;
let mut bands: Option<Vec<TwistBand>> = None;
if closed {
let holonomy = frame_holonomy(&dirs)?;
let round_off = 64.0 * f64::EPSILON * segments as f64;
let requested = if requested == 0.0 {
ClosingTwist {
twist: 0.0,
shift: 0,
}
} else {
let slice = rings[0]
.iter()
.map(|curve| project_curve_along(curve, dirs[0], joints[0], dirs[0]))
.collect::<Result<Vec<_>, String>>()?;
closing_twist(&slice, joints[0], dirs[0], requested, room_floor)?
};
if holonomy.abs() <= round_off && requested.twist == 0.0 {
let residual = ring_residual(&rings[segments], &rings[0])?;
if residual > room_floor {
return Err(format!(
"sweepSolid: this closed path's frame does NOT close. Carried round the loop \
and projected back onto the closing joint's own bisector plane, the section \
returns {residual:.6e} away from the ring it started as, past the \
{room_floor:.6e} floor — so the last mitre would not share one loop with the \
first and the frame would carry a seam. Its holonomy reads {:.3e}°, which is \
the joint rotations' own rounding ({round_off:.1e} rad), so the loop should \
have closed as a planar one does",
holonomy.to_degrees(),
));
}
closure = Some(SweepClosure {
holonomy,
requested_twist: 0.0,
applied_twist: 0.0,
curve_shift: 0,
});
} else {
let counter = if holonomy.abs() <= round_off { 0.0 } else { -holonomy };
let twist = if requested.twist == 0.0 {
counter
} else {
requested.twist + counter
};
let shares: Vec<f64> = lengths.iter().map(|length| twist * length / travel).collect();
let mut twisted: Vec<Vec<NurbsCurve>> = Vec::with_capacity(segments + 1);
twisted.push(rings[0].clone());
for index in 0..segments {
let mut next = Vec::with_capacity(count);
for curve in &twisted[index] {
let rolled = roll_curve_about(curve, joints[index], dirs[index], shares[index])?;
next.push(project_curve_along(
&rolled,
dirs[index],
plane_point[index + 1],
plane_normal[index + 1],
)?);
}
twisted.push(next);
}
let landed: Vec<NurbsCurve> = (0..count)
.map(|index| twisted[0][(index + requested.shift) % count].clone())
.collect();
let residual = ring_residual(&twisted[segments], &landed)?;
if residual > room_floor {
return Err(format!(
"sweepSolid: this closed path's frame does NOT close under its twist. Its \
holonomy is {:.6}°, a twist of {:.6}° (requested {:.6}°) was shared over its \
segments, yet the section returns {residual:.6e} from the ring it started \
as, past the {room_floor:.6e} floor",
holonomy.to_degrees(),
twist.to_degrees(),
requested.twist.to_degrees(),
));
}
let max_tilt = (1.0 / (MITER_TRANSVERSALITY_FLOOR * MITER_TRANSVERSALITY_FLOOR) - 1.0).sqrt();
let mut planned = Vec::with_capacity(segments);
for index in 0..segments {
let (origin, direction, length) = (joints[index], dirs[index], lengths[index]);
let (start, end) = (length / 4.0, 3.0 * length / 4.0);
let share = shares[index];
let mut behind = f64::NEG_INFINITY;
let mut reach = 0.0_f64;
for curve in &twisted[index] {
for control in &curve.control_points {
let offset = control.point()?.sub(origin);
let along = offset.dot(direction);
behind = behind.max(along);
reach = reach.max(offset.sub(direction.scale(along)).length());
}
}
let mut ahead = f64::INFINITY;
for curve in &twisted[index + 1] {
for control in &curve.control_points {
ahead = ahead.min(control.point()?.sub(origin).dot(direction));
}
}
let pieces = (share.abs() / std::f64::consts::FRAC_PI_2).ceil().max(1.0);
let rate = 3.0 * (share.abs() / pieces / 2.0).tan() / ((end - start) / pieces);
if reach * rate > max_tilt {
return Err(format!(
"{MITRED_TWIST_ROOM_REFUSAL}: path segment '{}' is too SHORT for its share of \
the twist. This loop's holonomy is {:.6}° and the twist requested on top of \
its counter-twist is {:.6}°, and the segment's {:.6}° share rolls over the \
middle half of its {length:.6} length in {pieces} piece(s), a peak rate of \
{rate:.6} rad per unit length (bounded as 3·tan(|piece|/2) over the piece's \
run). At the section's farthest reach from the path, {reach:.6}, that tilts \
the wall to ρ·ω = {:.6}, past the {max_tilt:.6} the {MITER_TRANSVERSALITY_FLOOR} \
transversality floor allows — the floor the mitre planes and the prism \
builder refuse a grazing sweep at. Lengthen the segment, narrow the section, \
or twist less",
path.name(index),
holonomy.to_degrees(),
requested.twist.to_degrees(),
share.to_degrees(),
reach * rate,
));
}
if behind > start - room_floor || ahead < end + room_floor {
return Err(format!(
"{MITRED_TWIST_ROOM_REFUSAL}: path segment '{}' has NO UNTRIMMED MIDDLE for \
its share of the twist. This loop comes back rolled by its holonomy, \
{:.6}°, with {:.6}° of twist requested on top, and the frame closes by \
rolling each segment through a share of the counter-twist plus that twist \
in proportion to its length — {:.6}° here — over the middle half of the \
segment, from {start:.6} to {end:.6} along its {length:.6}. The mitre ring \
behind reaches {behind:.6} along it and the ring ahead starts at {ahead:.6}, \
so a trim reaches into that run, and a roll inside a mitre's reach would \
pull its shared loop apart. The share is not moved onto a neighbour. \
Lengthen the segment, or narrow the section across its bends",
path.name(index),
holonomy.to_degrees(),
requested.twist.to_degrees(),
share.to_degrees(),
));
}
planned.push(TwistBand {
origin,
direction,
start,
end,
roll: share,
pieces: pieces as usize,
});
}
rings = twisted;
bands = Some(planned);
closure = Some(SweepClosure {
holonomy,
requested_twist: requested.twist,
applied_twist: twist,
curve_shift: requested.shift,
});
}
}
if bands.is_none() {
for index in 0..segments {
let denominator = dirs[index].dot(plane_normal[index + 1]);
let mut room = f64::INFINITY;
for curve in &rings[index] {
for control in &curve.control_points {
let point = control.point()?;
room = room.min(plane_point[index + 1].sub(point).dot(plane_normal[index + 1]) / denominator);
}
}
if room <= room_floor {
let mut bounds: Vec<String> = Vec::with_capacity(2);
if closed || index > 0 {
bounds.push(format!(
"{:.3}° where it meets '{}'",
turns[index].to_degrees(),
path.name((index + segments - 1) % segments)
));
}
if closed || index + 1 < segments {
bounds.push(format!(
"{:.3}° where it meets '{}'",
turns[index + 1].to_degrees(),
path.name((index + 1) % segments)
));
}
return Err(format!(
"sweepSolid: the bend is too TIGHT for the section at path segment '{}'. The \
segment is {:.6} long and {} the section {:.6} back into it ({}), \
leaving {:.6e} of wall at the innermost point of the section: the bisector plane \
would cut the section a SECOND time, at the segment's far end, and the inner \
side of the bend would intersect itself.{} Lengthen the segment, open the bend, \
or narrow the section across the bend",
path.name(index),
lengths[index],
if closed { "its two mitres fold" } else { "its mitre folds" },
lengths[index] - room,
bounds.join(" and "),
room,
if closed {
" On a FRAME every segment is trimmed at BOTH ends, so the two mitres OVERLAP \
here: a section wider across the bend than a side is long leaves that side \
no wall at all."
} else {
""
},
));
}
}
}
Ok(MitrePlan {
closed,
segments,
count,
ring_count,
rings,
reversed_winding,
start_frame: (x_axis, y_axis),
end_frame: (end_x, end_y),
bands,
closure,
})
}
pub(crate) const MITRED_TWIST_ROOM_REFUSAL: &str =
"sweepSolid: a spatial frame's counter-twist has no room";
#[derive(Debug, Clone, Copy)]
struct TwistBand {
origin: Vec3,
direction: Vec3,
start: f64,
end: f64,
roll: f64,
pieces: usize,
}
struct MitrePlan {
closed: bool,
segments: usize,
count: usize,
ring_count: usize,
rings: Vec<Vec<NurbsCurve>>,
reversed_winding: bool,
start_frame: (Vec3, Vec3),
end_frame: (Vec3, Vec3),
bands: Option<Vec<TwistBand>>,
closure: Option<SweepClosure>,
}
pub(super) fn mitred_frame_closure(
section: &[NurbsCurve],
section_normal: Vec3,
path: &SweepPath,
twist: f64,
) -> Result<SweepClosure, String> {
plan_mitre(section, section_normal, path, twist)?
.closure
.ok_or_else(|| "sweepSolid: an OPEN mitred run has no lap to close".to_string())
}
fn ring_residual(returned: &[NurbsCurve], start: &[NurbsCurve]) -> Result<f64, String> {
let mut residual = 0.0_f64;
for (a_curve, b_curve) in returned.iter().zip(start) {
for (a, b) in a_curve.control_points.iter().zip(&b_curve.control_points) {
residual = residual.max(a.point()?.sub(b.point()?).length());
}
}
Ok(residual)
}
fn roll_curve_about(
curve: &NurbsCurve,
origin: Vec3,
direction: Vec3,
angle: f64,
) -> Result<NurbsCurve, String> {
let (sin, cos) = angle.sin_cos();
let control_points = curve
.control_points
.iter()
.map(|control| {
let offset = control.point()?.sub(origin);
let along = direction.scale(offset.dot(direction));
let across = offset.sub(along);
let rolled = along
.add(across.scale(cos))
.add(direction.cross(across).scale(sin));
Ok(Vec4::from_point(origin.add(rolled), control.w))
})
.collect::<Result<Vec<_>, String>>()?;
NurbsCurve::new(curve.degree, curve.knots.clone(), control_points)
}
const TWISTED_WALL_DEGREE: usize = 7;
fn twisted_wall_knots(pieces: usize) -> Vec<f64> {
let degree = TWISTED_WALL_DEGREE;
let mut knots = Vec::with_capacity(2 * (degree + 1) + (pieces + 1) * degree);
knots.extend(std::iter::repeat(0.0).take(degree + 1));
knots.extend(std::iter::repeat(0.25).take(degree));
for piece in 1..pieces {
knots.extend(std::iter::repeat(0.25 + 0.5 * piece as f64 / pieces as f64).take(degree));
}
knots.extend(std::iter::repeat(0.75).take(degree));
knots.extend(std::iter::repeat(1.0).take(degree + 1));
knots
}
fn binomial(n: usize, k: usize) -> f64 {
(0..k).fold(1.0, |value, index| value * (n - index) as f64 / (index + 1) as f64)
}
fn bernstein_product(a: &[f64], b: &[f64]) -> Vec<f64> {
let (n, m) = (a.len() - 1, b.len() - 1);
(0..=n + m)
.map(|k| {
let low = k.saturating_sub(m);
let high = k.min(n);
(low..=high)
.map(|i| binomial(n, i) * binomial(m, k - i) * a[i] * b[k - i])
.sum::<f64>()
/ binomial(n + m, k)
})
.collect()
}
fn bernstein_elevate(a: &[f64], by: usize) -> Vec<f64> {
bernstein_product(a, &vec![1.0; by + 1])
}
fn twisted_column(
bottom: Vec3,
weight: f64,
top: Vec3,
band: &TwistBand,
) -> Vec<Vec4> {
let degree = TWISTED_WALL_DEGREE;
let TwistBand {
origin,
direction,
start,
end,
roll,
pieces,
} = *band;
let piece_roll = roll / pieces as f64;
let piece_run = (end - start) / pieces as f64;
let u = {
let half = (piece_roll / 2.0).tan();
[0.0, 0.0, half, half]
};
let squared = bernstein_product(&u, &u);
let w6: Vec<f64> = squared.iter().map(|value| 1.0 + value).collect();
let c6: Vec<f64> = squared.iter().map(|value| 1.0 - value).collect();
let s6 = bernstein_elevate(&u.map(|value| 2.0 * value), 3);
let w = bernstein_elevate(&w6, 1);
let c = bernstein_elevate(&c6, 1);
let s = bernstein_elevate(&s6, 1);
let l = bernstein_product(&w6, &[0.0, piece_run]);
let first = bottom.add(direction.scale(start - bottom.sub(origin).dot(direction)));
let mut column = Vec::with_capacity((pieces + 2) * degree + 1);
for index in 0..=degree {
let t = index as f64 / degree as f64;
column.push(Vec4::from_point(bottom.add(first.sub(bottom).scale(t)), weight));
}
let mut last = first;
let mut carried = weight;
for piece in 0..pieces {
let hub = origin.add(direction.scale(start + piece_run * piece as f64));
let offset = if piece == 0 {
last.sub(hub)
} else {
let offset = last.sub(hub);
offset.sub(direction.scale(offset.dot(direction)))
};
let across = direction.cross(offset);
let rolled = |index: usize| {
hub.add(
offset
.scale(c[index])
.add(across.scale(s[index]))
.add(direction.scale(l[index]))
.scale(1.0 / w[index]),
)
};
for index in 1..=degree {
column.push(Vec4::from_point(rolled(index), carried * w[index]));
}
last = rolled(degree);
carried *= w[degree];
}
for index in 1..=degree {
let t = index as f64 / degree as f64;
column.push(Vec4::from_point(last.add(top.sub(last).scale(t)), carried));
}
column
}
fn twisted_wall(
bottom: &NurbsCurve,
top: &NurbsCurve,
band: &TwistBand,
) -> Result<NurbsSurface, String> {
if bottom.degree != top.degree
|| bottom.control_points.len() != top.control_points.len()
|| bottom.knots.len() != top.knots.len()
{
return Err("twisted_wall: rows are not representation-compatible".into());
}
let grid = bottom
.control_points
.iter()
.zip(&top.control_points)
.map(|(a, b)| Ok(twisted_column(a.point()?, a.w, b.point()?, band)))
.collect::<Result<Vec<_>, String>>()?;
NurbsSurface::new(
bottom.degree,
TWISTED_WALL_DEGREE,
bottom.knots.clone(),
twisted_wall_knots(band.pieces),
grid,
)
}
fn twisted_side(
bottom: &NurbsCurve,
top: &NurbsCurve,
band: &TwistBand,
) -> Result<NurbsCurve, String> {
let (first, last) = match (bottom.control_points.first(), top.control_points.first()) {
(Some(first), Some(last)) => (first, last),
_ => return Err("twisted_side: a ring curve has no control points".into()),
};
NurbsCurve::new(
TWISTED_WALL_DEGREE,
twisted_wall_knots(band.pieces),
twisted_column(first.point()?, first.w, last.point()?, band),
)
}
fn rotate_minimal(x: Vec3, from: Vec3, to: Vec3) -> Vec3 {
let axis = from.cross(to);
let sin = axis.length();
if sin <= f64::EPSILON {
return x;
}
let cos = from.dot(to);
let axis = axis.scale(1.0 / sin);
x.scale(cos)
.add(axis.cross(x).scale(sin))
.add(axis.scale(axis.dot(x) * (1.0 - cos)))
}
fn frame_holonomy(dirs: &[Vec3]) -> Result<f64, String> {
let reference = dirs[0]
.perpendicular()
.map_err(|_| "sweepSolid: the first path segment's direction is degenerate".to_string())?;
let across = dirs[0].cross(reference);
let mut carried = reference;
for index in 1..dirs.len() {
carried = rotate_minimal(carried, dirs[index - 1], dirs[index]);
}
carried = rotate_minimal(carried, dirs[dirs.len() - 1], dirs[0]);
Ok(carried.dot(across).atan2(carried.dot(reference)))
}
fn segment_gap(p0: Vec3, p1: Vec3, q0: Vec3, q1: Vec3) -> f64 {
let u = p1.sub(p0);
let v = q1.sub(q0);
let w = p0.sub(q0);
let (a, b, c) = (u.dot(u), u.dot(v), v.dot(v));
let (d, e) = (u.dot(w), v.dot(w));
let determinant = a * c - b * b;
let mut gap = f64::INFINITY;
if determinant > 0.0 {
let s = (b * e - c * d) / determinant;
let t = (a * e - b * d) / determinant;
if (0.0..=1.0).contains(&s) && (0.0..=1.0).contains(&t) {
gap = p0.add(u.scale(s)).sub(q0.add(v.scale(t))).length();
}
}
for (point, start, direction, squared) in
[(p0, q0, v, c), (p1, q0, v, c), (q0, p0, u, a), (q1, p0, u, a)]
{
let parameter = if squared > 0.0 {
(point.sub(start).dot(direction) / squared).clamp(0.0, 1.0)
} else {
0.0
};
gap = gap.min(point.sub(start.add(direction.scale(parameter))).length());
}
gap
}
fn project_curve_along(
curve: &NurbsCurve,
direction: Vec3,
plane_point: Vec3,
plane_normal: Vec3,
) -> Result<NurbsCurve, String> {
let denominator = direction.dot(plane_normal);
if denominator.abs() < MITER_TRANSVERSALITY_FLOOR {
return Err(
"sweepSolid: a mitre plane grazes the segment it trims (the transversality floor is \
checked before this point; reaching it means the planes and the directions \
disagree)"
.into(),
);
}
let control_points = curve
.control_points
.iter()
.map(|control| {
let point = control.point()?;
let advance = plane_point.sub(point).dot(plane_normal) / denominator;
Ok(Vec4::from_point(point.add(direction.scale(advance)), control.w))
})
.collect::<Result<Vec<_>, String>>()?;
NurbsCurve::new(curve.degree, curve.knots.clone(), control_points)
}