551 lines
29 KiB
Rust
551 lines
29 KiB
Rust
//! Multi-section blends: one face that runs through many cross-sections and still leaves the
|
||
//! first and last of them the way a support surface wants.
|
||
//!
|
||
//! The sections run in order along the face's **u**, each sampled at the same `samples` stations
|
||
//! along the face's **v**. Section `i` sits at `u = i/(n−1)`, and the face passes through it
|
||
//! exactly: its v-curve there is that section's own cross-section curve through its stations.
|
||
//! Along v the stations are interpolated by a B-spline of degree `min(3, samples − 1)` — the same
|
||
//! interpolation [`crate::bridge`] builds its boundary curves with — so every station of every
|
||
//! section lies on the face.
|
||
//!
|
||
//! Along u the face is a B-spline whose knots are clamped at both ends and repeated at every
|
||
//! interior section, so every span is one Bézier; a span's leading control rows are the Hermite
|
||
//! rows of `bridge::section` — position, tangent (`leave`), curvature (`bend`) and change of
|
||
//! curvature (`bend_change`), scaled to that span's parameter interval. The degree is the largest
|
||
//! `c + c′ − 1` over the spans (`c` conditions at the section the span leaves, `c′` at the one it
|
||
//! arrives at), and a span whose two ends carry fewer conditions than that is elevated exactly,
|
||
//! which moves no point of the curve. A blend of two sections is therefore the [`crate::bridge`]
|
||
//! of those two boundaries, only with the two parameters the other way round.
|
||
//!
|
||
//! **The first and last section** take their tangent, curvature and change of curvature from the
|
||
//! caller's frames, exactly as a bridge takes them from a support surface's boundary, and hold
|
||
//! them to the requested [`Continuity`] (G0 up to G3). The `leave` direction of the last section's
|
||
//! frames points back into the blend, towards the interior, as a bridge's far boundary does.
|
||
//!
|
||
//! **Interior sections** take nothing from the caller's frames but their positions: the run
|
||
//! direction's tangent there is the secant of the two neighbouring sections' stations, its
|
||
//! curvature their second difference and its change of curvature their third difference, all
|
||
//! along the run. The two spans that meet at such a section are given exactly those rows, so they
|
||
//! agree in ∂/∂u, ∂²/∂u² and ∂³/∂u³ there and the joint is **C³** in the run parameter (hence G³)
|
||
//! whenever a fourth station is there to take the third difference over — that is, whenever the
|
||
//! blend runs through at least four sections. A blend of exactly three sections has only its one
|
||
//! interior section's two neighbours to read, so no third difference can be taken: the interior
|
||
//! section then keeps position, tangent and curvature only and its joint is **C²** (G²). Dropping
|
||
//! interior conditions, never the surface, is also how this blend stays within
|
||
//! [`MAX_SPLINE_DEGREE`]: a section family asks for at most `c + c′ − 1 ≤ 7` at the four conditions
|
||
//! a section can carry today, and a caller that ever asked for more would lose interior
|
||
//! conditions (change of curvature first, then curvature, then tangent) rather than fail.
|
||
//!
|
||
//! `crown` scales how fast the face leaves the first and the last section (0 the plain blend, up
|
||
//! to 3 fuller, down to −0.9 flatter), exactly as a bridge's crown does. It shapes the middle of
|
||
//! each end span and leaves every joint's conditions alone.
|
||
|
||
use solidcraft_geom::Vec3;
|
||
use solidcraft_geom::bspline::{CurveN, HPoint, MAX_SPLINE_DEGREE, NurbsSurface};
|
||
use truck_modeling as mt;
|
||
use truck_modeling::{ParametricSurface3D, SearchNearestParameter};
|
||
|
||
use crate::body::{Body, from_p3, p3};
|
||
use crate::bridge::{Continuity, Frame};
|
||
use crate::surfaces::surface_body;
|
||
use crate::{KernelError, Result, guard};
|
||
|
||
/// Most cross-sections one blend may run through.
|
||
pub const MAX_BLEND_SECTIONS: usize = 50;
|
||
/// Fewest stations one cross-section may be sampled at (a straight line across it).
|
||
pub const MIN_BLEND_STATIONS: usize = 2;
|
||
/// Most stations one cross-section may be sampled at. The control net (and the interpolation that
|
||
/// builds it) grows with the station count, and no cross-section curve needs more than this.
|
||
pub const MAX_BLEND_STATIONS: usize = 257;
|
||
|
||
/// How many rows of the control net one interpolation solve carries at a time: a row is three
|
||
/// coordinates, so this many rows are solved together.
|
||
const ROWS_PER_SOLVE: usize = 16;
|
||
|
||
/// A blend surface through `sections`, in order, that leaves the first and the last of them with
|
||
/// the continuity their frames ask for.
|
||
///
|
||
/// Each section is one cross-section: `samples` stations (`samples` must match every section's
|
||
/// length, so that the stations correspond from one section to the next). Section `i` is
|
||
/// interpolated at `u = i/(n−1)`, and its stations run along `v`, so the face of the body holds
|
||
/// every station of every section. `orders[0]` and `orders[1]` are what the face holds at the
|
||
/// first and the last section, read from their frames; the interior sections join the face
|
||
/// smoothly by finite differences along the run (see the module comment for the order each joint
|
||
/// holds). `crown` scales how fast the face leaves the first (`crown[0]`) and the last
|
||
/// (`crown[1]`) section and must be within −0.9 and 3.
|
||
///
|
||
/// Refused, without ever panicking: fewer than 2 or more than [`MAX_BLEND_SECTIONS`] sections,
|
||
/// sections of unequal station counts or a `samples` below [`MIN_BLEND_STATIONS`] or above
|
||
/// [`MAX_BLEND_STATIONS`], a frame that is not finite or a `leave` with no direction, a crown
|
||
/// outside its range, and two neighbouring sections that meet at a station (no room for a blend).
|
||
pub fn loft_blend(sections: &[Vec<Frame>], orders: [Continuity; 2], crown: [f64; 2], samples: usize) -> Result<Body> {
|
||
let blend = Blend::new(sections, orders, crown, samples)?;
|
||
let net = blend.net()?;
|
||
guard("loft blend surface", || crate::spline_conv::tensor_face(&net).map(|f| surface_body(vec![f])))
|
||
}
|
||
|
||
/// The frames of one cross-section picked as an edge chain of a body, sampled at `stations`
|
||
/// stations: the chain itself, with the neighbouring face's tangent plane, curvature vector and
|
||
/// change of curvature read along `leave` — the direction the blend runs on from this section —
|
||
/// exactly as [`crate::bridge::frames`] reads a support surface's own boundary.
|
||
///
|
||
/// Hand the result to [`loft_blend`] as one section. The face the section repeats is the body's
|
||
/// face that holds the whole chain and faces the run best; `leave` must have a part in that face's
|
||
/// tangent plane (a direction normal to it has no tangent plane to leave in). The last section's
|
||
/// `leave` runs back into the blend, towards the section before it.
|
||
pub fn section_frames(body: &Body, edges: &[Vec3], leave: Vec3, stations: usize) -> Result<Vec<Frame>> {
|
||
guard("loft blend section", || {
|
||
let chain = crate::bridge::chain(body, edges)?;
|
||
let leave =
|
||
leave.normalized().filter(|l| l.is_finite()).ok_or_else(|| KernelError::Invalid("the leave direction must be a finite vector".into()))?;
|
||
let first = chain.first().ok_or_else(|| KernelError::Failed("loft blend: no chain".into()))?;
|
||
let mid = (from_p3(first.front().point()) + from_p3(first.back().point())) * 0.5;
|
||
// The face the chain lies on that the blend leaves through: the one whose tangent plane
|
||
// holds as much of `leave` as possible.
|
||
let mut best: Option<(f64, mt::Face)> = None;
|
||
for face in body.solid.face_iter() {
|
||
if !chain.iter().all(|e| face.edge_iter().any(|f| f.id() == e.id())) {
|
||
continue;
|
||
}
|
||
let surf = face.surface();
|
||
let Some((u, v)) = surf.search_nearest_parameter(p3(mid), None, 100) else { continue };
|
||
let normal = surf.normal(u, v);
|
||
let n = Vec3::new(normal.x, normal.y, normal.z).normalized().unwrap_or(Vec3::ZERO);
|
||
let score = (leave - n * leave.dot(n)).len();
|
||
if best.as_ref().is_none_or(|(b, _)| score > *b) {
|
||
best = Some((score, face.clone()));
|
||
}
|
||
}
|
||
let (score, face) = best.ok_or_else(|| KernelError::Invalid("no face of the body holds the picked edges".into()))?;
|
||
if score < 0.2 {
|
||
return Err(KernelError::Invalid("the direction towards the next section lies in no face at the picked edges".into()));
|
||
}
|
||
crate::bridge::frames(&face, &chain, leave, stations)
|
||
})
|
||
}
|
||
|
||
/// What the face has to do at one section: where it is, and the run direction's first, second and
|
||
/// third derivatives at each of the section's stations. The third is carried only when the
|
||
/// section holds a change-of-curvature condition.
|
||
struct Shape {
|
||
p: Vec<Vec3>,
|
||
d1: Vec<Vec3>,
|
||
d2: Vec<Vec3>,
|
||
d3: Option<Vec<Vec3>>,
|
||
}
|
||
|
||
/// The blend's inputs, validated, with the conditions and the run direction's degree worked out.
|
||
struct Blend<'a> {
|
||
sections: &'a [Vec<Frame>],
|
||
crown: [f64; 2],
|
||
samples: usize,
|
||
run: usize,
|
||
/// The run parameter interval between two neighbouring sections: `1/(run − 1)`.
|
||
step: f64,
|
||
degree: usize,
|
||
/// Conditions held at each section, one per control row of the spans that end there.
|
||
conditions: Vec<usize>,
|
||
}
|
||
|
||
impl<'a> Blend<'a> {
|
||
fn new(sections: &'a [Vec<Frame>], orders: [Continuity; 2], crown: [f64; 2], samples: usize) -> Result<Blend<'a>> {
|
||
if !(MIN_BLEND_STATIONS..=MAX_BLEND_STATIONS).contains(&samples) {
|
||
return Err(KernelError::Invalid(format!(
|
||
"a blend needs {MIN_BLEND_STATIONS} to {MAX_BLEND_STATIONS} stations per cross-section, not {samples}"
|
||
)));
|
||
}
|
||
let run = sections.len();
|
||
if !(2..=MAX_BLEND_SECTIONS).contains(&run) {
|
||
return Err(KernelError::Invalid(format!("a blend runs through 2 to {MAX_BLEND_SECTIONS} cross-sections, not {run}")));
|
||
}
|
||
for c in crown {
|
||
if !c.is_finite() || !(-0.9..=3.0).contains(&c) {
|
||
return Err(KernelError::Invalid("a blend's crown must be between -0.9 and 3".into()));
|
||
}
|
||
}
|
||
for (i, section) in sections.iter().enumerate() {
|
||
if section.len() != samples {
|
||
return Err(KernelError::Invalid(format!(
|
||
"every cross-section of a blend must be sampled at the same {samples} stations (section {i} has {})",
|
||
section.len()
|
||
)));
|
||
}
|
||
for (s, f) in section.iter().enumerate() {
|
||
if ![f.p, f.leave, f.bend, f.bend_change].into_iter().all(Vec3::is_finite) {
|
||
return Err(KernelError::Invalid(format!("the frame at section {i}, station {s} is not finite")));
|
||
}
|
||
if f.leave.normalized().is_none() {
|
||
return Err(KernelError::Invalid(format!("the leave direction at section {i}, station {s} is zero")));
|
||
}
|
||
}
|
||
}
|
||
// The face runs from one section to the next: where two of them meet at a station there is
|
||
// no room for a blend between them.
|
||
for i in 0..run - 1 {
|
||
for s in 0..samples {
|
||
if !(point(sections, i, s)?.dist(point(sections, i + 1, s)?) > 1e-9) {
|
||
return Err(KernelError::Invalid(format!("sections {} and {} meet at station {s}: there is no room for a blend", i + 1, i + 2)));
|
||
}
|
||
}
|
||
}
|
||
// The ends hold what their orders ask for. An interior section would like all four
|
||
// conditions (position, tangent, curvature, change of curvature); where no fourth station
|
||
// is there to take the third difference over it keeps the three its neighbours can give,
|
||
// and its joint is C² instead of C³.
|
||
let mut conditions = vec![2usize; run];
|
||
if let Some(slot) = conditions.first_mut() {
|
||
*slot = orders[0].conditions();
|
||
}
|
||
if let Some(slot) = conditions.last_mut() {
|
||
*slot = orders[1].conditions();
|
||
}
|
||
for i in 1..run.saturating_sub(1) {
|
||
if let Some(slot) = conditions.get_mut(i) {
|
||
*slot = if third_difference_available(i, run) { 4 } else { 3 };
|
||
}
|
||
}
|
||
// The run direction's degree is the largest `c + c′ − 1` of a span. Four conditions a
|
||
// section is at most what a section can carry, so this fits today; should a caller ever
|
||
// ask for more, the interior joints give up conditions rather than the blend failing.
|
||
let mut degree = span_degree(&conditions);
|
||
while degree > MAX_SPLINE_DEGREE {
|
||
let next = (1..run.saturating_sub(1))
|
||
.filter(|&i| conditions.get(i).is_some_and(|&c| c > 1))
|
||
.max_by_key(|&i| conditions.get(i).copied().unwrap_or(0));
|
||
let Some(i) = next else { break };
|
||
if let Some(c) = conditions.get_mut(i) {
|
||
*c = c.saturating_sub(1);
|
||
}
|
||
degree = span_degree(&conditions);
|
||
}
|
||
if degree == 0 || degree > MAX_SPLINE_DEGREE {
|
||
return Err(KernelError::Failed("loft blend: the requested continuity needs a degree above the kernel's".into()));
|
||
}
|
||
Ok(Blend { sections, crown, samples, run, step: 1.0 / (run - 1) as f64, degree, conditions })
|
||
}
|
||
|
||
/// What the face has to do at section `i`: the run direction's derivatives at every station.
|
||
fn shape(&self, i: usize) -> Result<Shape> {
|
||
let samples = self.samples;
|
||
let p: Vec<Vec3> = (0..samples).map(|s| point(self.sections, i, s)).collect::<Result<Vec<_>>>()?;
|
||
let conditions = self.conditions.get(i).copied().unwrap_or(1);
|
||
let mut d1 = vec![Vec3::ZERO; samples];
|
||
let mut d2 = vec![Vec3::ZERO; samples];
|
||
let mut d3 = (conditions >= 4).then(|| vec![Vec3::ZERO; samples]);
|
||
if i == 0 {
|
||
for s in 0..samples {
|
||
let f = frame(self.sections, 0, s)?;
|
||
let leave = f.leave.normalized().ok_or_else(|| KernelError::Failed("loft blend: the leave direction".into()))?;
|
||
// The speed a bridge gives a boundary: a third of the station's chord to the
|
||
// section it leaves for, the crown, over this span's parameter interval.
|
||
let speed = f.p.dist(point(self.sections, 1, s)?) / 3.0 * (1.0 + self.crown[0]) / self.step;
|
||
if let (Some(v1), Some(v2)) = (d1.get_mut(s), d2.get_mut(s)) {
|
||
*v1 = leave * speed;
|
||
*v2 = f.bend * (speed * speed);
|
||
}
|
||
}
|
||
} else if i + 1 == self.run {
|
||
for s in 0..samples {
|
||
let f = frame(self.sections, i, s)?;
|
||
let leave = f.leave.normalized().ok_or_else(|| KernelError::Failed("loft blend: the leave direction".into()))?;
|
||
let speed = f.p.dist(point(self.sections, i - 1, s)?) / 3.0 * (1.0 + self.crown[1]) / self.step;
|
||
if let (Some(v1), Some(v2)) = (d1.get_mut(s), d2.get_mut(s)) {
|
||
// The run leaves the last section outwards; the caller's frames there point
|
||
// back into it, so the odd derivatives turn round.
|
||
*v1 = leave * -speed;
|
||
*v2 = f.bend * (speed * speed);
|
||
}
|
||
}
|
||
} else {
|
||
let h = self.step;
|
||
for s in 0..samples {
|
||
let before = point(self.sections, i - 1, s)?;
|
||
let here = point(self.sections, i, s)?;
|
||
let after = point(self.sections, i + 1, s)?;
|
||
if let (Some(v1), Some(v2)) = (d1.get_mut(s), d2.get_mut(s)) {
|
||
*v1 = (after - before) * (1.0 / (2.0 * h));
|
||
*v2 = (after - here * 2.0 + before) * (1.0 / (h * h));
|
||
}
|
||
}
|
||
}
|
||
// The change of curvature: the caller's bend_change at an end, a third difference of the
|
||
// neighbouring sections at an interior one.
|
||
if let Some(third) = d3.as_mut() {
|
||
for s in 0..samples {
|
||
let (Some(&a), Some(&b)) = (d1.get(s), d2.get(s)) else {
|
||
return Err(KernelError::Failed("loft blend: a section's derivatives".into()));
|
||
};
|
||
let value = if i == 0 {
|
||
let f = frame(self.sections, 0, s)?;
|
||
third_derivative(a, b, f.bend_change)?
|
||
} else if i + 1 == self.run {
|
||
let f = frame(self.sections, i, s)?;
|
||
-third_derivative(-a, b, f.bend_change)?
|
||
} else {
|
||
self.third_difference(i, s)?
|
||
};
|
||
if let Some(slot) = third.get_mut(s) {
|
||
*slot = value;
|
||
}
|
||
}
|
||
}
|
||
Ok(Shape { p, d1, d2, d3 })
|
||
}
|
||
|
||
/// The third difference of the sections' stations at an interior section: the estimate of how
|
||
/// the curvature vector changes along the run. `i` is an interior section, and the neighbours
|
||
/// reach far enough for a stencil (see [`third_difference_available`]).
|
||
fn third_difference(&self, i: usize, s: usize) -> Result<Vec3> {
|
||
let h = self.step;
|
||
let cube = h * h * h;
|
||
let before = i.checked_sub(2);
|
||
let value = match (before, i + 2 < self.run) {
|
||
// Both sides: the symmetric five-point difference.
|
||
(Some(k), true) => {
|
||
(point(self.sections, i + 2, s)? - point(self.sections, i + 1, s)? * 2.0 + point(self.sections, i - 1, s)? * 2.0
|
||
- point(self.sections, k, s)?)
|
||
* (1.0 / (2.0 * cube))
|
||
}
|
||
// Only ahead: the forward four-point stencil at the second of them.
|
||
(None, true) => {
|
||
(point(self.sections, i + 2, s)? - point(self.sections, i + 1, s)? * 3.0 + point(self.sections, i, s)? * 3.0
|
||
- point(self.sections, i - 1, s)?)
|
||
* (1.0 / cube)
|
||
}
|
||
// Only behind: the same the other way round.
|
||
(Some(k), false) => {
|
||
(point(self.sections, k, s)? - point(self.sections, i - 1, s)? * 3.0 + point(self.sections, i, s)? * 3.0
|
||
- point(self.sections, i + 1, s)?)
|
||
* (1.0 / cube)
|
||
}
|
||
(None, false) => return Err(KernelError::Failed("loft blend: no fourth station for a change of curvature".into())),
|
||
};
|
||
if !value.is_finite() {
|
||
return Err(KernelError::Failed("loft blend: the sections' third difference is not finite".into()));
|
||
}
|
||
Ok(value)
|
||
}
|
||
|
||
/// The Bézier rows of the span from section `j` to section `j + 1`: one control row per
|
||
/// condition at each end, elevated to the run direction's degree so the whole face is one
|
||
/// tensor surface. Row `k` of the span is control row `j·degree + k` of the face.
|
||
fn span(&self, j: usize, left: &Shape, right: &Shape) -> Result<Vec<Vec<Vec3>>> {
|
||
let (Some(&kl), Some(&kr)) = (self.conditions.get(j), self.conditions.get(j + 1)) else {
|
||
return Err(KernelError::Failed("loft blend: a section with no conditions".into()));
|
||
};
|
||
let m = kl.saturating_add(kr).saturating_sub(1);
|
||
if !(1..=MAX_SPLINE_DEGREE).contains(&m) {
|
||
return Err(KernelError::Failed("loft blend: a span's degree".into()));
|
||
}
|
||
let k = m as f64;
|
||
let delta = self.step;
|
||
let one = delta / k;
|
||
let two = if m >= 2 { delta * delta / (k * (k - 1.0)) } else { 0.0 };
|
||
let three = if m >= 3 { delta * delta * delta / (k * (k - 1.0) * (k - 2.0)) } else { 0.0 };
|
||
let mut out = vec![vec![Vec3::ZERO; self.samples]; m + 1];
|
||
let missing = || KernelError::Failed("loft blend: a section's derivatives".into());
|
||
for s in 0..self.samples {
|
||
// The near end: one row per condition, position first.
|
||
let p = *left.p.get(s).ok_or_else(missing)?;
|
||
put(&mut out, 0, s, p)?;
|
||
if kl >= 2 {
|
||
let d1 = *left.d1.get(s).ok_or_else(missing)?;
|
||
let q1 = p + d1 * one;
|
||
put(&mut out, 1, s, q1)?;
|
||
if kl >= 3 {
|
||
let d2 = *left.d2.get(s).ok_or_else(missing)?;
|
||
let q2 = d2 * two + q1 * 2.0 - p;
|
||
put(&mut out, 2, s, q2)?;
|
||
if kl >= 4 {
|
||
let d3 = left.d3.as_ref().and_then(|d| d.get(s)).copied().ok_or_else(missing)?;
|
||
put(&mut out, 3, s, d3 * three + q2 * 3.0 - q1 * 3.0 + p)?;
|
||
}
|
||
}
|
||
}
|
||
// The far end, its own conditions read backwards.
|
||
let e = *right.p.get(s).ok_or_else(missing)?;
|
||
put(&mut out, m, s, e)?;
|
||
if kr >= 2 {
|
||
let d1 = *right.d1.get(s).ok_or_else(missing)?;
|
||
let r1 = e - d1 * one;
|
||
put(&mut out, m - 1, s, r1)?;
|
||
if kr >= 3 {
|
||
let d2 = *right.d2.get(s).ok_or_else(missing)?;
|
||
let r2 = d2 * two + r1 * 2.0 - e;
|
||
put(&mut out, m - 2, s, r2)?;
|
||
if kr >= 4 {
|
||
let d3 = right.d3.as_ref().and_then(|d| d.get(s)).copied().ok_or_else(missing)?;
|
||
put(&mut out, m - 3, s, r2 * 3.0 - r1 * 3.0 + e - d3 * three)?;
|
||
}
|
||
}
|
||
}
|
||
}
|
||
// A span whose two ends carry fewer conditions than the run direction's degree is
|
||
// elevated: the same curve, written in more control rows.
|
||
while out.len().saturating_sub(1) < self.degree {
|
||
out = elevate(&out);
|
||
}
|
||
Ok(out)
|
||
}
|
||
|
||
/// The run direction's knots: clamped at both ends, repeated `degree` times at every interior
|
||
/// section — so the face interpolates each section and every span in between is one Bézier.
|
||
fn run_knots(&self) -> Vec<f64> {
|
||
let degree = self.degree;
|
||
let mut knots = vec![0.0; degree + 1];
|
||
for i in 1..self.run.saturating_sub(1) {
|
||
knots.extend(std::iter::repeat_n(i as f64 / (self.run - 1) as f64, degree));
|
||
}
|
||
knots.extend(std::iter::repeat_n(1.0, degree + 1));
|
||
knots
|
||
}
|
||
|
||
/// The face's control net: u runs with the sections, v along each section's curve.
|
||
fn net(&self) -> Result<NurbsSurface> {
|
||
let shapes: Vec<Shape> = (0..self.run).map(|i| self.shape(i)).collect::<Result<Vec<_>>>()?;
|
||
// One span per pair of neighbouring sections. The first row of a span is the last row of
|
||
// the one before it — the shared section — so only the first span writes it.
|
||
let mut rows: Vec<Vec<Vec3>> = Vec::with_capacity((self.run - 1) * self.degree + 1);
|
||
for j in 0..self.run.saturating_sub(1) {
|
||
let (Some(left), Some(right)) = (shapes.get(j), shapes.get(j + 1)) else {
|
||
return Err(KernelError::Failed("loft blend: a section's shape".into()));
|
||
};
|
||
for row in self.span(j, left, right)?.into_iter().skip(usize::from(j > 0)) {
|
||
rows.push(row);
|
||
}
|
||
}
|
||
let count_u = rows.len();
|
||
// Along each section's curve the rows are the interpolation through that section's
|
||
// stations: the same parameters and degree for every row, so every row has one knot
|
||
// vector. Rows are solved in batches of `ROWS_PER_SOLVE` (three coordinates each), which
|
||
// shares one elimination between them.
|
||
let params: Vec<f64> = (0..self.samples).map(|s| s as f64 / (self.samples - 1) as f64).collect();
|
||
let degree_v = 3.min(self.samples - 1);
|
||
let mut ctrl = vec![HPoint::ZERO; count_u * self.samples];
|
||
let mut knots_v: Option<Vec<f64>> = None;
|
||
for (batch, chunk) in rows.chunks(ROWS_PER_SOLVE).enumerate() {
|
||
const WIDTH: usize = 3 * ROWS_PER_SOLVE;
|
||
let mut points: Vec<[f64; WIDTH]> = Vec::with_capacity(self.samples);
|
||
for s in 0..self.samples {
|
||
let mut q = [0.0f64; WIDTH];
|
||
for (r, row) in chunk.iter().enumerate() {
|
||
if let Some(p) = row.get(s) {
|
||
for (slot, value) in q.iter_mut().skip(3 * r).zip([p.x, p.y, p.z]) {
|
||
*slot = value;
|
||
}
|
||
}
|
||
}
|
||
points.push(q);
|
||
}
|
||
let curve = CurveN::<WIDTH>::interpolate(¶ms, &points, degree_v).map_err(|e| KernelError::Failed(format!("loft blend: {e}")))?;
|
||
if knots_v.is_none() {
|
||
knots_v = Some(curve.knots.clone());
|
||
}
|
||
for (s, c) in curve.ctrl.iter().enumerate() {
|
||
for (r, xyz) in c.chunks(3).enumerate() {
|
||
let i = batch * ROWS_PER_SOLVE + r;
|
||
// Coordinates past the net's last row are this batch's own padding: their
|
||
// controls are not part of the surface.
|
||
if i >= count_u {
|
||
break;
|
||
}
|
||
if let [x, y, z] = *xyz {
|
||
let Some(slot) = s.checked_mul(count_u).and_then(|base| base.checked_add(i)).and_then(|k| ctrl.get_mut(k)) else {
|
||
continue;
|
||
};
|
||
*slot = HPoint::weighted(Vec3::new(x, y, z), 1.0);
|
||
}
|
||
}
|
||
}
|
||
}
|
||
let knots_v = knots_v.ok_or_else(|| KernelError::Failed("loft blend: no row could be interpolated".into()))?;
|
||
NurbsSurface::from_homogeneous([self.degree, degree_v], [self.run_knots(), knots_v], ctrl)
|
||
.map_err(|e| KernelError::Failed(format!("loft blend: {e}")))
|
||
}
|
||
}
|
||
|
||
/// The frame at section `i`, station `s`.
|
||
fn frame(sections: &[Vec<Frame>], i: usize, s: usize) -> Result<&Frame> {
|
||
sections.get(i).and_then(|section| section.get(s)).ok_or_else(|| KernelError::Failed("loft blend: a section is missing a station".into()))
|
||
}
|
||
|
||
/// The point of section `i` at station `s`.
|
||
fn point(sections: &[Vec<Frame>], i: usize, s: usize) -> Result<Vec3> {
|
||
frame(sections, i, s).map(|f| f.p)
|
||
}
|
||
|
||
/// `rows[k][s] = v`, an error when the net has no such row or station.
|
||
fn put(rows: &mut [Vec<Vec3>], k: usize, s: usize, v: Vec3) -> Result<()> {
|
||
match rows.get_mut(k).and_then(|row| row.get_mut(s)) {
|
||
Some(slot) => {
|
||
*slot = v;
|
||
Ok(())
|
||
}
|
||
None => Err(KernelError::Failed("loft blend: the control net is too small".into())),
|
||
}
|
||
}
|
||
|
||
/// One exact degree elevation of a Bézier control net of rows: the same curve, one control row
|
||
/// more.
|
||
fn elevate(rows: &[Vec<Vec3>]) -> Vec<Vec<Vec3>> {
|
||
let degree = rows.len().saturating_sub(1);
|
||
let up = degree as f64 + 1.0;
|
||
let mut out: Vec<Vec<Vec3>> = Vec::with_capacity(rows.len() + 1);
|
||
if let Some(first) = rows.first() {
|
||
out.push(first.clone());
|
||
}
|
||
for i in 1..=degree {
|
||
if let (Some(before), Some(here)) = (rows.get(i - 1), rows.get(i)) {
|
||
let a = i as f64 / up;
|
||
out.push(before.iter().zip(here).map(|(x, y)| *x * a + *y * (1.0 - a)).collect());
|
||
}
|
||
}
|
||
if let Some(last) = rows.last() {
|
||
out.push(last.clone());
|
||
}
|
||
out
|
||
}
|
||
|
||
/// Whether a third difference for the interior section `i` can be taken from the sections around
|
||
/// it: a stencil needs a station on one side of it and two on the other, and with three sections
|
||
/// there is none.
|
||
fn third_difference_available(i: usize, run: usize) -> bool {
|
||
i >= 2 || i + 2 < run
|
||
}
|
||
|
||
/// The degree a family of sections asks for: the largest `c + c′ − 1` over its spans.
|
||
fn span_degree(conditions: &[usize]) -> usize {
|
||
let mut degree = 0;
|
||
for w in conditions.windows(2) {
|
||
if let [a, b] = w {
|
||
degree = degree.max(a.saturating_add(*b).saturating_sub(1));
|
||
}
|
||
}
|
||
degree
|
||
}
|
||
|
||
/// The third derivative a curve needs at a boundary so that the arclength derivative of its
|
||
/// curvature vector is `target`, from its first two derivatives there. Only the normal-plane part
|
||
/// of the third derivative can hold that condition; the rest is a reparameterisation of the
|
||
/// cross-section. This is `bridge`'s own construction, kept here so the bridge core stays
|
||
/// untouched.
|
||
fn third_derivative(d1: Vec3, d2: Vec3, target: Vec3) -> Result<Vec3> {
|
||
let sigma = d1.len();
|
||
if !(sigma > 1e-12) || !sigma.is_finite() {
|
||
return Err(KernelError::Failed("loft blend: no room for the curvature condition".into()));
|
||
}
|
||
let t = d1 * (1.0 / sigma);
|
||
let tp = d2 * (1.0 / sigma) - d1 * (d1.dot(d2) / (sigma * sigma * sigma));
|
||
let nn = d2 - t * d2.dot(t);
|
||
let r = tp * (-d2.dot(t)) - t * d2.dot(tp);
|
||
let x = (target * (sigma.powi(5)) - r * (sigma * sigma) + nn * (2.0 * d1.dot(d2))) * (1.0 / (sigma * sigma));
|
||
if !x.is_finite() {
|
||
return Err(KernelError::Failed("loft blend: the change of curvature is not finite".into()));
|
||
}
|
||
Ok(x - t * x.dot(t))
|
||
}
|
||
|
||
#[cfg(test)]
|
||
#[path = "loft_blend_tests.rs"]
|
||
mod tests;
|