use crate::topology::EdgeRecord;
use crate::{
fit_polyline, intersect_analytic_pair, intersect_surfaces, project_point_to_surface, NurbsCurve,
NurbsSurface, SurfaceIntersectionOptions, Vec3, Vec4,
};
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub(crate) enum RimLane {
Analytic,
Marched,
}
pub(crate) struct RimSection {
pub(crate) curve: NurbsCurve,
pub(crate) polyline: Vec<Vec3>,
pub(crate) closed: bool,
}
pub(crate) struct RimIntersection {
pub(crate) sections: Vec<RimSection>,
pub(crate) lane: RimLane,
pub(crate) residual: f64,
}
impl RimIntersection {
pub(crate) fn curves(&self) -> Vec<NurbsCurve> {
self.sections
.iter()
.map(|section| section.curve.clone())
.collect()
}
pub(crate) fn nearest_section(&self, reference: &[Vec3]) -> Result<&RimSection, String> {
let mut best: Option<(f64, &RimSection)> = None;
for section in &self.sections {
let [t0, t1] = section.curve.domain()?;
let mut samples = Vec::with_capacity(SECTION_MATCH_SAMPLES);
for index in 0..SECTION_MATCH_SAMPLES {
let fraction = index as f64 / (SECTION_MATCH_SAMPLES - 1) as f64;
samples.push(section.curve.evaluate(t0 + (t1 - t0) * fraction)?);
}
let mut worst = 0.0f64;
for point in reference {
let nearest = samples
.iter()
.map(|sample| sample.sub(*point).length())
.fold(f64::INFINITY, f64::min);
worst = worst.max(nearest);
}
if best.map(|(score, _)| worst < score).unwrap_or(true) {
best = Some((worst, section));
}
}
best.map(|(_, section)| section)
.ok_or_else(|| "no section to match against the old boundary".to_string())
}
}
const SECTION_MATCH_SAMPLES: usize = 24;
#[derive(Clone, Debug)]
pub(crate) enum ReintersectRefusal {
Separated,
Residual { drift: f64, limit: f64 },
Failed(String),
}
impl ReintersectRefusal {
pub(crate) fn describe(&self) -> String {
match self {
Self::Separated => "the carriers do not meet".to_string(),
Self::Residual { drift, limit } => format!(
"the marched section drifts {drift:.3e} off a carrier (limit {limit:.3e})"
),
Self::Failed(message) => message.clone(),
}
}
}
pub(crate) struct MarchPolicy {
pub(crate) tolerance: f64,
pub(crate) residual_tolerance: f64,
pub(crate) seeds: Vec<Vec3>,
}
const RESIDUAL_SAMPLES: usize = 33;
const MAXIMUM_FIT_POINTS: usize = 80;
pub(crate) fn reintersect_carriers(
first: &NurbsSurface,
second: &NurbsSurface,
policy: &MarchPolicy,
) -> Result<RimIntersection, ReintersectRefusal> {
match intersect_analytic_pair(first, second, policy.tolerance) {
Some(curves) if !curves.is_empty() => {
return Ok(RimIntersection {
sections: curves
.into_iter()
.map(|curve| RimSection {
curve,
polyline: Vec::new(),
closed: false,
})
.collect(),
lane: RimLane::Analytic,
residual: 0.0,
})
}
Some(_) => return Err(ReintersectRefusal::Separated),
None => {}
}
let branches = intersect_surfaces(
first,
second,
&SurfaceIntersectionOptions {
tolerance: policy.tolerance,
seed_points: policy.seeds.clone(),
..SurfaceIntersectionOptions::default()
},
)
.map_err(ReintersectRefusal::Failed)?;
let mut sections: Vec<RimSection> = Vec::new();
let mut worst_residual = 0.0f64;
for branch in &branches {
if branch.points.len() < 2 {
continue;
}
let length: f64 = branch
.points
.windows(2)
.map(|pair| pair[1].sub(pair[0]).length())
.sum();
if length <= policy.tolerance * 100.0 {
continue;
}
let fit = fit_polyline(
&branch.points,
policy.tolerance.max(1e-7),
MAXIMUM_FIT_POINTS,
false,
)
.map_err(ReintersectRefusal::Failed)?;
let curve = if branch.closed {
close_exactly(&fit.curve, policy.tolerance).map_err(ReintersectRefusal::Failed)?
} else {
fit.curve
};
let residual =
section_residual(&curve, first, second).map_err(ReintersectRefusal::Failed)?;
worst_residual = worst_residual.max(residual);
sections.push(RimSection {
curve,
polyline: branch.points.clone(),
closed: branch.closed,
});
}
if sections.is_empty() {
return Err(ReintersectRefusal::Separated);
}
if worst_residual > policy.residual_tolerance {
return Err(ReintersectRefusal::Residual {
drift: worst_residual,
limit: policy.residual_tolerance,
});
}
Ok(RimIntersection {
sections,
lane: RimLane::Marched,
residual: worst_residual,
})
}
pub(crate) fn section_corner(section: &RimSection, old: Vec3) -> Result<Vec3, String> {
if !section.closed || section.polyline.len() < 4 {
return Err("a section that is not a traced closed loop has no corner samples".into());
}
let ring = §ion.polyline[..section.polyline.len() - 1];
ring.iter()
.copied()
.min_by(|a, b| {
a.sub(old)
.length()
.total_cmp(&b.sub(old).length())
})
.ok_or_else(|| "empty section ring".to_string())
}
pub(crate) fn arc_of_section(
section: &RimSection,
from: Vec3,
to: Vec3,
through: Vec3,
tolerance: f64,
) -> Result<NurbsCurve, String> {
if !section.closed || section.polyline.len() < 4 {
return Err("a section that is not a traced closed loop cannot be cut into arcs".into());
}
let ring = §ion.polyline[..section.polyline.len() - 1];
let nearest_index = |target: Vec3| -> usize {
let mut best = 0usize;
let mut best_distance = f64::INFINITY;
for (index, point) in ring.iter().enumerate() {
let distance = point.sub(target).length();
if distance < best_distance {
best_distance = distance;
best = index;
}
}
best
};
let start = nearest_index(from);
let end = nearest_index(to);
let count = ring.len();
let forward_len = (end + count - start) % count;
if forward_len < 2 || count - forward_len < 2 {
return Err(
"the two ends of an arc land on the same place of its section — refusing".into(),
);
}
let run = |begin: usize, step: isize| -> Vec<Vec3> {
let mut points = Vec::new();
let mut index = begin as isize;
loop {
points.push(ring[index.rem_euclid(count as isize) as usize]);
if index.rem_euclid(count as isize) as usize == end && points.len() > 1 {
break;
}
index += step;
}
points
};
let forward = run(start, 1);
let backward = run(start, -1);
let closeness = |points: &[Vec3]| -> f64 {
points
.iter()
.map(|point| point.sub(through).length())
.fold(f64::INFINITY, f64::min)
};
let mut chosen = if closeness(&forward) <= closeness(&backward) {
forward
} else {
backward
};
let last = chosen.len() - 1;
chosen[0] = from;
chosen[last] = to;
if chosen.len() < 3 {
return Err("a rebuilt arc has too few samples to fit — refusing".into());
}
Ok(fit_polyline(&chosen, tolerance.max(1e-7), MAXIMUM_FIT_POINTS, false)?.curve)
}
fn close_exactly(curve: &NurbsCurve, tolerance: f64) -> Result<NurbsCurve, String> {
let controls = &curve.control_points;
let (Some(first), Some(last)) = (controls.first(), controls.last()) else {
return Ok(curve.clone());
};
let gap = last.point()?.sub(first.point()?).length();
if gap == 0.0 {
return Ok(curve.clone());
}
if gap > tolerance.max(1e-9) * 10.0 {
return Err(format!(
"a marched section reported itself closed but its fit leaves a {gap:.3e} gap"
));
}
let mut controls = controls.clone();
let head = controls[0];
let tail_weight = controls[controls.len() - 1].w;
let head_point = head.point()?;
let index = controls.len() - 1;
controls[index] = Vec4::from_point(head_point, tail_weight);
NurbsCurve::new(curve.degree, curve.knots.clone(), controls)
}
fn section_residual(
curve: &NurbsCurve,
first: &NurbsSurface,
second: &NurbsSurface,
) -> Result<f64, String> {
let [t0, t1] = curve.domain()?;
let mut worst = 0.0f64;
for index in 0..RESIDUAL_SAMPLES {
let fraction = index as f64 / (RESIDUAL_SAMPLES - 1) as f64;
let point = curve.evaluate(t0 + (t1 - t0) * fraction)?;
worst = worst.max(project_point_to_surface(first, point)?.distance);
worst = worst.max(project_point_to_surface(second, point)?.distance);
}
Ok(worst)
}
pub(crate) fn edge_seeds(edge: &EdgeRecord, count: usize) -> Result<Vec<Vec3>, String> {
let count = count.max(2);
(0..count)
.map(|index| {
let fraction = index as f64 / (count - 1) as f64;
edge.curve.evaluate(edge.t0 + (edge.t1 - edge.t0) * fraction)
})
.collect()
}
pub(crate) fn match_marched_rim_direction(
rim: NurbsCurve,
previous: &EdgeRecord,
) -> Result<NurbsCurve, String> {
let [d0, _] = rim.domain()?;
let head = rim.evaluate(d0)?;
let projection = crate::project_point_to_curve(&previous.curve, head)?;
let parameter = projection.u.clamp(previous.t0.min(previous.t1), previous.t1.max(previous.t0));
let incoming = previous.curve.derivatives(parameter, 1)?;
let rebuilt = rim.derivatives(d0, 1)?;
if incoming[1].dot(rebuilt[1]) < 0.0 {
return rim.reversed();
}
Ok(rim)
}