use super::*;
pub fn sweep_profile_along_path(
profile: &[NurbsCurve],
path: &NurbsCurve,
name: Option<&str>,
) -> Result<BrepSolid, String> {
sweep_profile_along_path_stations(
profile,
path,
name,
32,
0.0,
None,
SectionPlacement::Transplant,
)
}
pub fn sweep_profile_twisted(
profile: &[NurbsCurve],
path: &NurbsCurve,
twist_angle: f64,
name: Option<&str>,
) -> Result<BrepSolid, String> {
use std::f64::consts::{FRAC_PI_2, TAU};
if !twist_angle.is_finite() {
return Err("sweep_profile_twisted: twist angle must be finite".into());
}
const MAX_TURNS: f64 = 16.0;
if twist_angle.abs() > MAX_TURNS * TAU {
return Err(format!(
"sweep_profile_twisted: twist of {:.3} turns exceeds the {MAX_TURNS}-turn \
limit the 1024-station cap can resolve at 16 stations per quarter turn; \
split the sweep or reduce the twist",
twist_angle.abs() / TAU
));
}
let quarter_turns = (twist_angle.abs() / FRAC_PI_2).ceil() as usize;
let stations = (quarter_turns * 16).clamp(32, 1024);
sweep_profile_along_path_stations(
profile,
path,
name,
stations,
twist_angle,
None,
SectionPlacement::Transplant,
)
.map_err(|error| format!("sweep_profile_twisted: {error}"))
}
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum SectionPlacement {
Transplant,
Rigid,
}
#[derive(Debug, Clone, Copy)]
pub struct ProfileAnchor {
pub origin: Vec3,
pub normal: Vec3,
pub pu: Vec3,
pub pv: Vec3,
}
pub fn sweep_profile_along_path_anchored(
profile: &[NurbsCurve],
path: &NurbsCurve,
name: Option<&str>,
anchor: ProfileAnchor,
) -> Result<BrepSolid, String> {
sweep_profile_along_path_stations(
profile,
path,
name,
32,
0.0,
Some(anchor),
SectionPlacement::Transplant,
)
}
pub fn sweep_profile_twisted_anchored(
profile: &[NurbsCurve],
path: &NurbsCurve,
twist_angle: f64,
name: Option<&str>,
anchor: ProfileAnchor,
) -> Result<BrepSolid, String> {
use std::f64::consts::{FRAC_PI_2, TAU};
if !twist_angle.is_finite() {
return Err("sweep_profile_twisted: twist angle must be finite".into());
}
const MAX_TURNS: f64 = 16.0;
if twist_angle.abs() > MAX_TURNS * TAU {
return Err(format!(
"sweep_profile_twisted: twist of {:.3} turns exceeds the {MAX_TURNS}-turn \
limit the 1024-station cap can resolve at 16 stations per quarter turn; \
split the sweep or reduce the twist",
twist_angle.abs() / TAU
));
}
let quarter_turns = (twist_angle.abs() / FRAC_PI_2).ceil() as usize;
let stations = (quarter_turns * 16).clamp(32, 1024);
sweep_profile_along_path_stations(
profile,
path,
name,
stations,
twist_angle,
Some(anchor),
SectionPlacement::Transplant,
)
.map_err(|error| format!("sweep_profile_twisted: {error}"))
}
pub fn sweep_profile_along_chain(
profile: &[NurbsCurve],
chain: &[NurbsCurve],
segment_names: &[String],
twist_angle: f64,
name: Option<&str>,
anchor: Option<ProfileAnchor>,
placement_mode: SectionPlacement,
corner_advice: &str,
) -> Result<BrepSolid, String> {
use std::f64::consts::{FRAC_PI_2, TAU};
if chain.is_empty() {
return Err("sweepSolid: path chain is empty".into());
}
if let ([single], SectionPlacement::Transplant) = (chain, placement_mode) {
return match (anchor, twist_angle == 0.0) {
(None, true) => sweep_profile_along_path(profile, single, name),
(Some(anchor), true) => {
sweep_profile_along_path_anchored(profile, single, name, anchor)
}
(None, false) => sweep_profile_twisted(profile, single, twist_angle, name),
(Some(anchor), false) => {
sweep_profile_twisted_anchored(profile, single, twist_angle, name, anchor)
}
};
}
if !twist_angle.is_finite() {
return Err("sweep_profile_along_chain: twist angle must be finite".into());
}
const MAX_TURNS: f64 = 16.0;
if twist_angle.abs() > MAX_TURNS * TAU {
return Err(format!(
"sweep_profile_along_chain: twist of {:.3} turns exceeds the {MAX_TURNS}-turn \
limit the 1024-station cap can resolve at 16 stations per quarter turn; \
split the sweep or reduce the twist",
twist_angle.abs() / TAU
));
}
let stations = if twist_angle == 0.0 {
32
} else {
let quarter_turns = (twist_angle.abs() / FRAC_PI_2).ceil() as usize;
(quarter_turns * 16).clamp(32, 1024)
};
let _ = name;
sweep_profile_along_chain_stations(
profile,
chain,
segment_names,
stations,
twist_angle,
anchor,
placement_mode,
corner_advice,
)
}
pub fn sweep_profile_along_chain_with_stations(
profile: &[NurbsCurve],
chain: &[NurbsCurve],
segment_names: &[String],
stations: usize,
corner_advice: &str,
) -> Result<BrepSolid, String> {
sweep_profile_along_chain_stations(
profile,
chain,
segment_names,
stations.clamp(2, 1024),
0.0,
None,
SectionPlacement::Transplant,
corner_advice,
)
}
pub fn profile_anchor(profile: &[NurbsCurve]) -> Result<ProfileAnchor, String> {
let tolerance = 1e-6;
if profile.len() < 2 {
return Err("sweepSolid: profile needs at least 2 curves forming a closed loop".into());
}
let mut samples = Vec::new();
for (index, curve) in profile.iter().enumerate() {
let [start, end] = curve.domain()?;
let next = &profile[(index + 1) % profile.len()];
let next_start = next.domain()?[0];
if curve
.evaluate(end)?
.sub(next.evaluate(next_start)?)
.length()
> tolerance
{
return Err(format!(
"sweepSolid: profile is not closed at curve {index}"
));
}
for sample in 0..16 {
samples.push(curve.evaluate(start + (end - start) * sample as f64 / 16.0)?);
}
}
let normal = crate::polygon::newell_normal(&samples);
let centroid = samples.iter().fold(Vec3::default(), |sum, &point| sum.add(point));
let np = normal
.normalized()
.map_err(|_| "sweepSolid: profile is degenerate (zero enclosed area)".to_string())?;
let origin = centroid.scale(1.0 / samples.len() as f64);
if samples
.iter()
.any(|point| point.sub(origin).dot(np).abs() > tolerance * 100.0)
{
return Err("sweepSolid: profile is not planar".into());
}
let pu = np.perpendicular()?;
let pv = np.cross(pu).normalized()?;
Ok(ProfileAnchor {
origin,
normal: np,
pu,
pv,
})
}
fn sweep_profile_along_path_stations(
profile: &[NurbsCurve],
path: &NurbsCurve,
name: Option<&str>,
stations: usize,
twist_angle: f64,
anchor: Option<ProfileAnchor>,
placement_mode: SectionPlacement,
) -> Result<BrepSolid, String> {
let _ = name;
let tolerance = 1e-6;
if stations < 2 {
return Err("sweepSolid: need at least 2 stations".into());
}
let placement = resolve_placement(profile, anchor, placement_mode)?;
let [t0, t1] = path.domain()?;
if (t1 - t0).abs() <= tolerance {
return Err("sweepSolid: path domain is degenerate".into());
}
let mut points = Vec::with_capacity(stations);
let mut tangents = Vec::with_capacity(stations);
for index in 0..stations {
let t = t0 + (t1 - t0) * index as f64 / (stations - 1) as f64;
let derivatives = path.derivatives(t, 1)?;
let tangent = derivatives[1]
.normalized()
.map_err(|_| format!("sweepSolid: path tangent is degenerate at station {index}"))?;
points.push(derivatives[0]);
tangents.push(tangent);
}
sweep_sections_through_samples(
profile,
&points,
&tangents,
twist_angle,
placement,
placement_mode,
)
}
fn resolve_placement(
profile: &[NurbsCurve],
anchor: Option<ProfileAnchor>,
placement_mode: SectionPlacement,
) -> Result<ProfileAnchor, String> {
let computed = profile_anchor(profile)?;
Ok(match placement_mode {
SectionPlacement::Rigid => computed,
SectionPlacement::Transplant => anchor.unwrap_or(computed),
})
}
const MAX_JOINT_TANGENT_BREAK: f64 = 0.02;
fn sweep_profile_along_chain_stations(
profile: &[NurbsCurve],
chain: &[NurbsCurve],
segment_names: &[String],
stations: usize,
twist_angle: f64,
anchor: Option<ProfileAnchor>,
placement_mode: SectionPlacement,
corner_advice: &str,
) -> Result<BrepSolid, String> {
let tolerance = 1e-6;
if chain.is_empty() {
return Err("sweepSolid: path chain is empty".into());
}
if stations < 2 {
return Err("sweepSolid: need at least 2 stations".into());
}
let placement = resolve_placement(profile, anchor, placement_mode)?;
let mut domains = Vec::with_capacity(chain.len());
let mut lengths = Vec::with_capacity(chain.len());
for (index, curve) in chain.iter().enumerate() {
let [t0, t1] = curve.domain()?;
if (t1 - t0).abs() <= tolerance {
return Err(format!(
"sweepSolid: path segment {index} has a degenerate domain"
));
}
let mut length = 0.0;
let mut previous = curve.evaluate(t0)?;
for sample in 1..=16 {
let point = curve.evaluate(t0 + (t1 - t0) * sample as f64 / 16.0)?;
length += point.sub(previous).length();
previous = point;
}
domains.push([t0, t1]);
lengths.push(length);
}
let total_length: f64 = lengths.iter().sum();
if total_length <= tolerance {
return Err("sweepSolid: path chain has zero length".into());
}
let mut budget: Vec<usize> = lengths
.iter()
.map(|length| {
let share = (stations as f64 * length / total_length).round() as usize;
share.max(2)
})
.collect();
let joints = chain.len() - 1;
let emitted: usize = budget.iter().sum::<usize>() - joints;
const MAX_STATIONS: usize = 1024;
if emitted > MAX_STATIONS {
return Err(format!(
"sweepSolid: a {}-segment path needs at least {emitted} stations, past the \
{MAX_STATIONS}-station cap the loft solve can carry; sweep it in fewer, \
longer pieces",
chain.len()
));
}
if emitted < stations {
let longest = lengths
.iter()
.enumerate()
.max_by(|a, b| a.1.total_cmp(b.1))
.map(|(index, _)| index)
.unwrap_or(0);
budget[longest] += stations - emitted;
}
let mut points: Vec<Vec3> = Vec::new();
let mut tangents: Vec<Vec3> = Vec::new();
for (index, curve) in chain.iter().enumerate() {
let [t0, t1] = domains[index];
let count = budget[index];
if index > 0 {
let departing = curve.derivatives(t0, 1)?[1].normalized().map_err(|_| {
format!("sweepSolid: path tangent is degenerate at the start of segment {index}")
})?;
let arriving = *tangents
.last()
.expect("a previous segment emitted stations");
let break_angle = arriving.dot(departing).clamp(-1.0, 1.0).acos();
if break_angle > MAX_JOINT_TANGENT_BREAK {
let previous_name = segment_names
.get(index - 1)
.map(String::as_str)
.unwrap_or("<unnamed>");
let name = segment_names
.get(index)
.map(String::as_str)
.unwrap_or("<unnamed>");
return Err(format!(
"the path must be tangent-continuous: segments '{previous_name}' and \
'{name}' meet at a {:.1}° corner. The profile is skinned between sampled \
stations along the path, so a corner would be ROUNDED OFF rather than \
mitred. {corner_advice}",
break_angle.to_degrees()
));
}
}
let first = usize::from(index > 0);
for sample in first..count {
let t = t0 + (t1 - t0) * sample as f64 / (count - 1) as f64;
let derivatives = curve.derivatives(t, 1)?;
let tangent = derivatives[1].normalized().map_err(|_| {
format!("sweepSolid: path tangent is degenerate on segment {index}")
})?;
points.push(derivatives[0]);
tangents.push(tangent);
}
}
sweep_sections_through_samples(
profile,
&points,
&tangents,
twist_angle,
placement,
placement_mode,
)
}
fn sweep_sections_through_samples(
profile: &[NurbsCurve],
points: &[Vec3],
tangents: &[Vec3],
twist_angle: f64,
placement: ProfileAnchor,
placement_mode: SectionPlacement,
) -> Result<BrepSolid, String> {
let tolerance = 1e-6;
let stations = points.len();
if stations < 2 {
return Err("sweepSolid: need at least 2 stations".into());
}
let ProfileAnchor { origin, pu, pv, .. } = placement;
let twist_fractions: Option<Vec<f64>> = if twist_angle != 0.0 {
let mut cumulative = vec![0.0; stations];
let mut total = 0.0;
for index in 1..stations {
total += points[index].sub(points[index - 1]).length();
cumulative[index] = total;
}
if total <= tolerance {
return Err("sweepSolid: path has zero length; cannot distribute the twist".into());
}
for length in &mut cumulative {
*length /= total;
}
Some(cumulative)
} else {
None
};
let mut r_axes = Vec::with_capacity(stations);
let mut s_axes = Vec::with_capacity(stations);
let r0 = tangents[0].perpendicular()?; s_axes.push(tangents[0].cross(r0).normalized()?);
r_axes.push(r0);
for index in 0..stations - 1 {
let t_next = tangents[index + 1];
let v1 = points[index + 1].sub(points[index]);
let c1 = v1.dot(v1);
let r_candidate = if c1 <= 1e-18 {
r_axes[index]
} else {
let reflected_r = r_axes[index].sub(v1.scale(2.0 / c1 * v1.dot(r_axes[index])));
let reflected_t = tangents[index].sub(v1.scale(2.0 / c1 * v1.dot(tangents[index])));
let v2 = t_next.sub(reflected_t);
let c2 = v2.dot(v2);
if c2 <= 1e-18 {
reflected_r
} else {
reflected_r.sub(v2.scale(2.0 / c2 * v2.dot(reflected_r)))
}
};
let r_next = r_candidate
.sub(t_next.scale(r_candidate.dot(t_next)))
.normalized()
.map_err(|_| format!("sweepSolid: frame degenerated at station {index}"))?;
s_axes.push(t_next.cross(r_next).normalized()?);
r_axes.push(r_next);
}
{
let travel: f64 = points
.windows(2)
.map(|pair| pair[1].sub(pair[0]).length())
.sum();
let closure = points[stations - 1].sub(points[0]).length();
if travel > tolerance && closure <= 1e-3 * travel {
return Err(format!(
"sweepSolid: the path returns to where it started ({closure:.3e} apart after \
travelling {travel:.3e}), so the sweep's two end caps would land on top of each \
other. A closed path swept this way IS a revolution — build it with Revolve \
about the same axis, or sweep the run in two halves and union them"
));
}
}
let path_start = points[0];
let (t0, r0, s0) = (tangents[0], r_axes[0], s_axes[0]);
let rigid = placement_mode == SectionPlacement::Rigid;
let mut sections: Vec<Vec<NurbsCurve>> = Vec::with_capacity(stations);
let mut station_points: Vec<Vec<Vec3>> = Vec::with_capacity(if rigid { stations } else { 0 });
let mut section_normals: Vec<Vec3> = Vec::with_capacity(if rigid { stations } else { 0 });
let mut anchor_track: Vec<Vec3> = Vec::with_capacity(if rigid { stations } else { 0 });
for station in 0..stations {
let station_origin = points[station];
let (ri, si) = match &twist_fractions {
Some(fractions) => {
let phi = twist_angle * fractions[station];
let (sin_phi, cos_phi) = phi.sin_cos();
let r = r_axes[station];
let s = s_axes[station];
(
r.scale(cos_phi).add(s.scale(sin_phi)),
s.scale(cos_phi).sub(r.scale(sin_phi)),
)
}
None => (r_axes[station], s_axes[station]),
};
let tk = tangents[station];
let place = |x: Vec3| -> Vec3 {
if rigid {
let local = x.sub(path_start);
station_origin
.add(tk.scale(local.dot(t0)))
.add(ri.scale(local.dot(r0)))
.add(si.scale(local.dot(s0)))
} else {
let local = x.sub(origin);
station_origin
.add(ri.scale(local.dot(pu)))
.add(si.scale(local.dot(pv)))
}
};
if rigid {
let n0 = placement.normal;
section_normals.push(
tk.scale(n0.dot(t0))
.add(ri.scale(n0.dot(r0)))
.add(si.scale(n0.dot(s0))),
);
anchor_track.push(place(origin));
station_points.push(Vec::with_capacity(
profile.iter().map(|c| c.control_points.len()).sum(),
));
}
let mut section = Vec::with_capacity(profile.len());
for curve in profile {
let control_points = curve
.control_points
.iter()
.map(|point| {
let weight = point.w;
let euclidean = Vec3::new(point.x / weight, point.y / weight, point.z / weight);
let world = place(euclidean);
if rigid {
station_points[station].push(world);
}
Vec4 {
x: world.x * weight,
y: world.y * weight,
z: world.z * weight,
w: weight,
}
})
.collect();
section.push(NurbsCurve::new(
curve.degree,
curve.knots.clone(),
control_points,
)?);
}
sections.push(section);
}
if rigid {
let mut most_positive = 0.0f64;
let mut most_negative = 0.0f64;
let mut longest_step = 0.0f64;
for station in 0..stations - 1 {
let normal = section_normals[station];
for (before, after) in station_points[station]
.iter()
.zip(&station_points[station + 1])
{
let step = after.sub(*before);
longest_step = longest_step.max(step.length());
let advance = step.dot(normal);
most_positive = most_positive.max(advance);
most_negative = most_negative.min(advance);
}
}
let noise = 1e-9 * longest_step.max(tolerance);
if most_positive > noise && most_negative < -noise {
return Err(
"sweepSolid: the profile sweeps back through itself — part of it advances along \
the path while part of it retreats. A turning path carries the profile around \
the turn's own axis, so a profile that STRADDLES that axis folds into itself; \
move the profile clear of the axis, or sweep the run in pieces"
.into(),
);
}
if most_positive <= noise && most_negative >= -noise {
return Err(
"sweepSolid: the path does not advance through the profile — it runs inside the \
profile plane, so the sweep encloses no volume"
.into(),
);
}
for station in 0..stations - 1 {
let step = anchor_track[station + 1].sub(anchor_track[station]);
let length = step.length();
if length <= noise {
continue;
}
if step.dot(section_normals[station]).abs() / length < 0.1 {
return Err(
"sweepSolid: the path is nearly parallel to the profile plane, which sweeps a \
sliver rather than a solid"
.into(),
);
}
}
}
loft_profile_brep(§ions)
.map_err(|error| format!("sweepSolid: loft through swept sections failed: {error}"))
}
#[allow(clippy::too_many_arguments)]
pub fn helix_sample_points(
axis_origin: Vec3,
axis_direction: Vec3,
reference: Option<Vec3>,
start_radius: f64,
end_radius: f64,
pitch: f64,
turns: f64,
start_angle: f64,
left_handed: bool,
) -> Result<(Vec<Vec3>, Vec<f64>), String> {
use std::f64::consts::TAU;
let w = axis_direction
.normalized()
.map_err(|_| "helix: axis direction is degenerate".to_string())?;
for (name, value) in [
("radius", start_radius),
("end radius", end_radius),
("pitch", pitch),
("turns", turns),
("start angle", start_angle),
] {
if !value.is_finite() {
return Err(format!("helix: {name} must be a finite number"));
}
}
if start_radius < 0.0 || end_radius < 0.0 {
return Err("helix: radius must not be negative".into());
}
if pitch < 0.0 {
return Err("helix: pitch must not be negative".into());
}
if turns <= 0.0 {
return Err("helix: turns must be positive".into());
}
if turns > 256.0 {
return Err("helix: turns must be at most 256".into());
}
if start_radius == 0.0 && end_radius == 0.0 && pitch == 0.0 {
return Err("helix: zero radius and zero pitch describe a single point".into());
}
let u = match reference {
Some(reference) => {
let radial = reference.sub(w.scale(reference.dot(w)));
radial.normalized().or_else(|_| w.perpendicular())?
}
None => w.perpendicular()?,
};
let v = w.cross(u).normalized()?;
let total_angle = turns * TAU;
let height = pitch * turns;
let sign = if left_handed { -1.0 } else { 1.0 };
let count = ((turns * 64.0).ceil() as usize + 1).clamp(9, 4097);
let mut points = Vec::with_capacity(count);
let mut parameters = Vec::with_capacity(count);
for index in 0..count {
let s = index as f64 / (count - 1) as f64;
let theta = start_angle + sign * total_angle * s;
let radius = start_radius + (end_radius - start_radius) * s;
points.push(
axis_origin
.add(u.scale(radius * theta.cos()))
.add(v.scale(radius * theta.sin()))
.add(w.scale(height * s)),
);
parameters.push(s);
}
Ok((points, parameters))
}
#[allow(clippy::too_many_arguments)]
pub fn fit_helix_curve(
axis_origin: Vec3,
axis_direction: Vec3,
reference: Option<Vec3>,
start_radius: f64,
end_radius: f64,
pitch: f64,
turns: f64,
start_angle: f64,
left_handed: bool,
) -> Result<NurbsCurve, String> {
let (points, parameters) = helix_sample_points(
axis_origin,
axis_direction,
reference,
start_radius,
end_radius,
pitch,
turns,
start_angle,
left_handed,
)?;
interpolate_curve(&points, 3, ¶meters)
}
pub fn sweep_profile_helix(
profile: &[NurbsCurve],
axis_origin: Vec3,
axis_direction: Vec3,
helix_radius: f64,
pitch: f64,
turns: f64,
name: Option<&str>,
) -> Result<BrepSolid, String> {
use std::f64::consts::TAU;
let w = axis_direction
.normalized()
.map_err(|_| "sweep_profile_helix: axis direction is degenerate".to_string())?;
if !(helix_radius.is_finite() && helix_radius > 0.0) {
return Err("sweep_profile_helix: helix radius must be positive".into());
}
if !(pitch.is_finite() && pitch > 0.0) {
return Err("sweep_profile_helix: pitch must be positive".into());
}
if !(turns.is_finite() && turns > 0.0) {
return Err("sweep_profile_helix: turns must be positive".into());
}
if turns > 64.0 {
return Err("sweep_profile_helix: turns must be at most 64".into());
}
let mut samples = Vec::new();
for curve in profile {
let [start, end] = curve.domain()?;
for sample in 0..16 {
samples.push(curve.evaluate(start + (end - start) * sample as f64 / 16.0)?);
}
}
if !samples.is_empty() {
let mut centroid = Vec3::default();
for point in &samples {
centroid = centroid.add(*point);
}
let origin = centroid.scale(1.0 / samples.len() as f64);
let extent = samples
.iter()
.map(|point| point.sub(origin).length())
.fold(0.0, f64::max);
let c = pitch / TAU; let curvature_radius = (helix_radius * helix_radius + c * c) / helix_radius;
if extent >= curvature_radius {
return Err(format!(
"sweep_profile_helix: profile extent {extent:.6} reaches the helix \
curvature radius {curvature_radius:.6}; the tube would fold through \
itself — increase the helix radius or pitch, or shrink the profile"
));
}
if turns >= 1.0 {
let circumference = TAU * helix_radius;
let turn_length = (circumference * circumference + pitch * pitch).sqrt();
let coil_gap = pitch * circumference / turn_length;
if coil_gap <= 2.0 * extent {
return Err(format!(
"sweep_profile_helix: consecutive turns would self-intersect — \
coil gap {coil_gap:.6} does not clear the profile diameter {:.6}; \
increase the pitch or shrink the profile",
2.0 * extent
));
}
}
}
let path = fit_helix_curve(
axis_origin,
w,
None,
helix_radius,
helix_radius,
pitch,
turns,
0.0,
false,
)
.map_err(|error| format!("sweep_profile_helix: helix path fit failed: {error}"))?;
let stations = ((turns * 32.0).ceil() as usize).clamp(32, 1024);
sweep_profile_along_path_stations(
profile,
&path,
name,
stations,
0.0,
None,
SectionPlacement::Transplant,
)
.map_err(|error| format!("sweep_profile_helix: {error}"))
}