use crate::bezier::{de_casteljau_rational_2d, de_casteljau_rational_3d, interpolate_frame};
use crate::vec3;
pub type VertexId = usize;
pub type EdgeId = usize;
pub type FaceId = usize;
const VERTEX_TOL: f64 = 1e-10;
#[derive(Debug, Clone)]
pub enum Surface {
Plane {
origin: [f64; 3],
normal: [f64; 3],
},
Cylinder {
origin: [f64; 3],
axis: [f64; 3],
radius: f64,
},
Cone {
origin: [f64; 3],
axis: [f64; 3],
half_angle: f64,
},
Sphere {
center: [f64; 3],
radius: f64,
},
Ellipsoid {
center: [f64; 3],
rx: f64,
ry: f64,
rz: f64,
},
Torus {
center: [f64; 3],
axis: [f64; 3],
major_radius: f64,
minor_radius: f64,
},
SurfaceOfRevolution {
center: [f64; 3],
axis: [f64; 3],
frame_u: [f64; 3],
frame_v: [f64; 3],
profile_control_points: Vec<[f64; 2]>,
profile_weights: Vec<f64>,
profile_degree: u32,
n_profile_spans: u32,
theta_start: f64,
theta_range: f64,
},
SurfaceOfSweep {
spine_control_points: Vec<[f64; 3]>,
spine_weights: Vec<f64>,
spine_degree: u32,
profile_control_points: Vec<[f64; 2]>,
profile_weights: Vec<f64>,
profile_degree: u32,
n_profile_spans: u32,
frames: Vec<[[f64; 3]; 3]>,
},
NurbsSurface {
data: Box<neco_nurbs::NurbsSurface3D>,
},
}
impl Surface {
pub fn evaluate(&self, u: f64, v: f64) -> [f64; 3] {
match self {
Surface::Plane { origin, normal } => {
let n = vec3::normalized(*normal);
let up = if n[0].abs() < 0.9 {
[1.0, 0.0, 0.0]
} else {
[0.0, 1.0, 0.0]
};
let u_vec = vec3::normalized(vec3::cross(n, up));
let v_vec = vec3::cross(n, u_vec);
vec3::add(
vec3::add(*origin, vec3::scale(u_vec, u)),
vec3::scale(v_vec, v),
)
}
Surface::Cylinder {
origin,
axis,
radius,
} => {
let axis_len = vec3::length(*axis);
let axis_n = vec3::normalized(*axis);
let (bu, bv) = vec3::orthonormal_basis(*axis);
let h = v.clamp(0.0, axis_len);
vec3::add(
vec3::add(*origin, vec3::scale(axis_n, h)),
vec3::add(
vec3::scale(bu, radius * u.cos()),
vec3::scale(bv, radius * u.sin()),
),
)
}
Surface::Cone {
origin,
axis,
half_angle,
} => {
let axis_len = vec3::length(*axis);
let axis_n = vec3::normalized(*axis);
let (bu, bv) = vec3::orthonormal_basis(*axis);
let t = v.clamp(0.0, axis_len);
let r = t * half_angle.tan();
vec3::add(
vec3::add(*origin, vec3::scale(axis_n, t)),
vec3::add(vec3::scale(bu, r * u.cos()), vec3::scale(bv, r * u.sin())),
)
}
Surface::Sphere { center, radius } => {
let sin_v = v.sin();
let cos_v = v.cos();
[
center[0] + radius * sin_v * u.cos(),
center[1] + radius * sin_v * u.sin(),
center[2] + radius * cos_v,
]
}
Surface::Ellipsoid { center, rx, ry, rz } => {
let sin_v = v.sin();
let cos_v = v.cos();
[
center[0] + rx * sin_v * u.cos(),
center[1] + ry * sin_v * u.sin(),
center[2] + rz * cos_v,
]
}
Surface::Torus {
center,
axis,
major_radius,
minor_radius,
} => {
let axis_n = vec3::normalized(*axis);
let (bu, bv) = vec3::orthonormal_basis(*axis);
let r = major_radius + minor_radius * v.cos();
let dir = vec3::add(vec3::scale(bu, u.cos()), vec3::scale(bv, u.sin()));
vec3::add(
vec3::add(*center, vec3::scale(dir, r)),
vec3::scale(axis_n, minor_radius * v.sin()),
)
}
Surface::SurfaceOfRevolution {
center,
axis,
frame_u,
frame_v,
profile_control_points,
profile_weights,
profile_degree,
n_profile_spans,
..
} => {
let (r, z) = eval_revolution_profile(
profile_control_points,
profile_weights,
*profile_degree,
*n_profile_spans,
v,
);
let axis_n = vec3::normalized(*axis);
let bu = vec3::normalized(*frame_u);
let bv = vec3::normalized(*frame_v);
let theta = u;
vec3::add(
vec3::add(
*center,
vec3::add(
vec3::scale(bu, r * theta.cos()),
vec3::scale(bv, r * theta.sin()),
),
),
vec3::scale(axis_n, z),
)
}
Surface::SurfaceOfSweep {
spine_control_points,
spine_weights,
spine_degree,
profile_control_points,
profile_weights,
profile_degree,
n_profile_spans,
frames,
} => {
let deg = *spine_degree as usize;
let n_spine_spans = if deg > 0 {
(spine_control_points.len() - 1) / deg
} else {
1
};
let (spine_pt, normal, binormal) = if n_spine_spans <= 1 {
let pt = de_casteljau_rational_3d(spine_control_points, spine_weights, u);
let (n, b) = interpolate_frame(frames, u);
(pt, n, b)
} else {
let u_clamped = u.clamp(0.0, 1.0);
let span_f = u_clamped * n_spine_spans as f64;
let span_idx = (span_f as usize).min(n_spine_spans - 1);
let local_u = span_f - span_idx as f64;
let start_cp = span_idx * deg;
let end_cp = start_cp + deg + 1;
let pt = de_casteljau_rational_3d(
&spine_control_points[start_cp..end_cp],
&spine_weights[start_cp..end_cp],
local_u,
);
let (n, b) = interpolate_frame(&frames[span_idx..span_idx + 2], local_u);
(pt, n, b)
};
let n_prof = *n_profile_spans as usize;
let p_deg = *profile_degree as usize;
let cps_per_span = p_deg + 1;
let (px, py) = if n_prof <= 1 {
de_casteljau_rational_2d(profile_control_points, profile_weights, v)
} else {
let v_clamped = v.clamp(0.0, 1.0);
let span_f = v_clamped * n_prof as f64;
let span_idx = (span_f as usize).min(n_prof - 1);
let local_v = span_f - span_idx as f64;
let start = span_idx * cps_per_span;
let end = start + cps_per_span;
de_casteljau_rational_2d(
&profile_control_points[start..end],
&profile_weights[start..end],
local_v,
)
};
[
spine_pt[0] + px * normal[0] + py * binormal[0],
spine_pt[1] + px * normal[1] + py * binormal[1],
spine_pt[2] + px * normal[2] + py * binormal[2],
]
}
Surface::NurbsSurface { data } => data.evaluate(u, v),
}
}
pub fn normal_at(&self, u: f64, v: f64) -> [f64; 3] {
match self {
Surface::Plane { normal, .. } => vec3::normalized(*normal),
Surface::Cylinder { axis, .. } => {
let (bu, bv) = vec3::orthonormal_basis(*axis);
vec3::normalized(vec3::add(
vec3::scale(bu, u.cos()),
vec3::scale(bv, u.sin()),
))
}
Surface::Cone {
axis, half_angle, ..
} => {
let (bu, bv) = vec3::orthonormal_basis(*axis);
let cos_ha = half_angle.cos();
let sin_ha = half_angle.sin();
let a = vec3::normalized(*axis);
let radial = vec3::add(vec3::scale(bu, u.cos()), vec3::scale(bv, u.sin()));
vec3::normalized(vec3::add(
vec3::scale(radial, cos_ha),
vec3::scale(a, -sin_ha),
))
}
Surface::Sphere { center, radius } => {
let p = self.evaluate(u, v);
[
(p[0] - center[0]) / radius,
(p[1] - center[1]) / radius,
(p[2] - center[2]) / radius,
]
}
Surface::Ellipsoid { center, rx, ry, rz } => {
let p = self.evaluate(u, v);
let nx = (p[0] - center[0]) / (rx * rx);
let ny = (p[1] - center[1]) / (ry * ry);
let nz = (p[2] - center[2]) / (rz * rz);
let len = (nx * nx + ny * ny + nz * nz).sqrt();
if len < 1e-30 {
[0.0, 0.0, 1.0]
} else {
[nx / len, ny / len, nz / len]
}
}
Surface::Torus {
center,
axis,
major_radius,
..
} => {
let (bu, bv) = vec3::orthonormal_basis(*axis);
let p = self.evaluate(u, v);
let dir = vec3::add(vec3::scale(bu, u.cos()), vec3::scale(bv, u.sin()));
let ring_center = vec3::add(*center, vec3::scale(dir, *major_radius));
vec3::normalized(vec3::sub(p, ring_center))
}
Surface::SurfaceOfRevolution {
axis,
frame_u,
frame_v,
profile_control_points,
profile_weights,
profile_degree,
n_profile_spans,
..
} => {
let (r, _z) = eval_revolution_profile(
profile_control_points,
profile_weights,
*profile_degree,
*n_profile_spans,
v,
);
let axis_n = vec3::normalized(*axis);
let bu = vec3::normalized(*frame_u);
let bv = vec3::normalized(*frame_v);
let theta = u;
let cos_t = theta.cos();
let sin_t = theta.sin();
let ds_du = vec3::add(vec3::scale(bu, -r * sin_t), vec3::scale(bv, r * cos_t));
let eps = 1e-6;
let v_lo = (v - eps).max(0.0);
let v_hi = (v + eps).min(1.0);
let (r_lo, z_lo) = eval_revolution_profile(
profile_control_points,
profile_weights,
*profile_degree,
*n_profile_spans,
v_lo,
);
let (r_hi, z_hi) = eval_revolution_profile(
profile_control_points,
profile_weights,
*profile_degree,
*n_profile_spans,
v_hi,
);
let dt = v_hi - v_lo;
let dr = (r_hi - r_lo) / dt;
let dz = (z_hi - z_lo) / dt;
let ds_dv = vec3::add(
vec3::add(vec3::scale(bu, dr * cos_t), vec3::scale(bv, dr * sin_t)),
vec3::scale(axis_n, dz),
);
let normal = vec3::cross(ds_du, ds_dv);
let len = vec3::length(normal);
if len < 1e-30 {
axis_n
} else {
vec3::scale(normal, 1.0 / len)
}
}
Surface::SurfaceOfSweep { .. } => {
let eps = 1e-6;
let u_lo = (u - eps).max(0.0);
let u_hi = (u + eps).min(1.0);
let v_lo = (v - eps).max(0.0);
let v_hi = (v + eps).min(1.0);
let du = vec3::sub(self.evaluate(u_hi, v), self.evaluate(u_lo, v));
let dv = vec3::sub(self.evaluate(u, v_hi), self.evaluate(u, v_lo));
let n = vec3::cross(du, dv);
let len = vec3::length(n);
if len < 1e-12 {
[0.0, 0.0, 1.0]
} else {
vec3::scale(n, 1.0 / len)
}
}
Surface::NurbsSurface { data } => data.normal(u, v),
}
}
pub fn to_nurbs_surface(&self) -> Option<neco_nurbs::NurbsSurface3D> {
match self {
Surface::SurfaceOfSweep {
spine_control_points,
spine_weights,
spine_degree,
profile_control_points,
profile_weights,
profile_degree,
n_profile_spans,
frames,
} => {
let s_deg = *spine_degree as usize;
let n_spine_spans = if s_deg > 0 {
(spine_control_points.len() - 1) / s_deg
} else {
1
};
let n_prof = *n_profile_spans as usize;
let (v_degree, prof_cps, prof_ws) = if n_prof <= 1 {
let deg = profile_control_points.len() - 1;
(deg, profile_control_points.clone(), profile_weights.clone())
} else {
let p_deg = *profile_degree as usize;
let cps_per_span = p_deg + 1;
let mut cps = Vec::new();
let mut ws = Vec::new();
for span_j in 0..n_prof {
let base = span_j * cps_per_span;
let start_k = if span_j == 0 { 0 } else { 1 };
for k in start_k..cps_per_span {
cps.push(profile_control_points[base + k]);
ws.push(profile_weights[base + k]);
}
}
(p_deg, cps, ws)
};
let n_spine_pts = n_spine_spans + 1;
let mut spine_pts = Vec::with_capacity(n_spine_pts);
let mut spine_wts = Vec::with_capacity(n_spine_pts);
let mut normals = Vec::with_capacity(n_spine_pts);
let mut binormals = Vec::with_capacity(n_spine_pts);
for si in 0..n_spine_pts {
let pt = if n_spine_spans <= 1 {
let t = si as f64 / n_spine_spans.max(1) as f64;
de_casteljau_rational_3d(spine_control_points, spine_weights, t)
} else if si < n_spine_spans {
spine_control_points[si * s_deg]
} else {
spine_control_points[spine_control_points.len() - 1]
};
let w = if n_spine_spans <= 1 {
if si == 0 {
spine_weights[0]
} else {
spine_weights[spine_weights.len() - 1]
}
} else if si < n_spine_spans {
spine_weights[si * s_deg]
} else {
spine_weights[spine_weights.len() - 1]
};
let n = frames[si][0];
let b = frames[si][1];
spine_pts.push(pt);
spine_wts.push(w);
normals.push(n);
binormals.push(b);
}
let n_v_cps = prof_cps.len();
let mut rows_cp: Vec<Vec<[f64; 3]>> = Vec::with_capacity(n_spine_pts);
let mut rows_w: Vec<Vec<f64>> = Vec::with_capacity(n_spine_pts);
for i in 0..n_spine_pts {
let sp = spine_pts[i];
let sw = spine_wts[i];
let n = &normals[i];
let b = &binormals[i];
let mut row_cp = Vec::with_capacity(n_v_cps);
let mut row_w = Vec::with_capacity(n_v_cps);
for j in 0..n_v_cps {
let [px, py] = prof_cps[j];
let pw = prof_ws[j];
let cp = [
sp[0] + px * n[0] + py * b[0],
sp[1] + px * n[1] + py * b[1],
sp[2] + px * n[2] + py * b[2],
];
row_cp.push(cp);
row_w.push(sw * pw);
}
rows_cp.push(row_cp);
rows_w.push(row_w);
}
let mut knots_u = vec![0.0; 2];
for i in 1..n_spine_spans {
let t = i as f64 / n_spine_spans as f64;
knots_u.push(t);
}
knots_u.extend(vec![1.0; 2]);
let knots_v = if n_prof <= 1 {
let mut kv = vec![0.0; v_degree + 1];
kv.extend(vec![1.0; v_degree + 1]);
kv
} else {
let p_deg = *profile_degree as usize;
let mut kv = vec![0.0; p_deg + 1];
for i in 1..n_prof {
let t = i as f64 / n_prof as f64;
for _ in 0..p_deg {
kv.push(t);
}
}
kv.extend(vec![1.0; p_deg + 1]);
kv
};
Some(neco_nurbs::NurbsSurface3D {
degree_u: 1,
degree_v: v_degree,
control_points: rows_cp,
weights: rows_w,
knots_u,
knots_v,
})
}
Surface::SurfaceOfRevolution {
center,
axis,
frame_u,
frame_v,
profile_control_points,
profile_weights,
profile_degree,
n_profile_spans,
theta_start,
theta_range,
} => {
use std::f64::consts::FRAC_PI_2;
let axis_n = vec3::normalized(*axis);
let bu = vec3::normalized(*frame_u);
let bv = vec3::normalized(*frame_v);
let n_u_spans = (*theta_range / FRAC_PI_2).ceil().max(1.0) as usize;
let span_angle = *theta_range / n_u_spans as f64;
let half_angle = span_angle / 2.0;
let w_mid = half_angle.cos();
let n_u_cps = 2 * n_u_spans + 1;
let mut u_cos = Vec::with_capacity(n_u_cps);
let mut u_sin = Vec::with_capacity(n_u_cps);
let mut u_weights = Vec::with_capacity(n_u_cps);
for span_i in 0..n_u_spans {
let theta0 = *theta_start + span_i as f64 * span_angle;
let theta_mid = theta0 + half_angle;
let theta1 = theta0 + span_angle;
if span_i == 0 {
u_cos.push(theta0.cos());
u_sin.push(theta0.sin());
u_weights.push(1.0);
}
u_cos.push(theta_mid.cos() / w_mid);
u_sin.push(theta_mid.sin() / w_mid);
u_weights.push(w_mid);
u_cos.push(theta1.cos());
u_sin.push(theta1.sin());
u_weights.push(1.0);
}
let mut knots_u = vec![*theta_start; 3];
for i in 1..n_u_spans {
let t = *theta_start + i as f64 * span_angle;
knots_u.push(t);
knots_u.push(t);
}
let theta_end = *theta_start + *theta_range;
knots_u.extend(vec![theta_end; 3]);
let n_prof = *n_profile_spans as usize;
let (v_degree, prof_cps, prof_ws) = if n_prof <= 1 {
let deg = profile_control_points.len() - 1;
(deg, profile_control_points.clone(), profile_weights.clone())
} else {
let p_deg = *profile_degree as usize;
let cps_per_span = p_deg + 1;
let mut cps = Vec::new();
let mut ws = Vec::new();
for span_j in 0..n_prof {
let base = span_j * cps_per_span;
let start_k = if span_j == 0 { 0 } else { 1 };
for k in start_k..cps_per_span {
cps.push(profile_control_points[base + k]);
ws.push(profile_weights[base + k]);
}
}
(p_deg, cps, ws)
};
let knots_v = if n_prof <= 1 {
let mut kv = vec![0.0; v_degree + 1];
kv.extend(vec![1.0; v_degree + 1]);
kv
} else {
let p_deg = *profile_degree as usize;
let mut kv = vec![0.0; p_deg + 1];
for i in 1..n_prof {
let t = i as f64 / n_prof as f64;
for _ in 0..p_deg {
kv.push(t);
}
}
kv.extend(vec![1.0; p_deg + 1]);
kv
};
let n_v_cps = prof_cps.len();
let mut rows_cp: Vec<Vec<[f64; 3]>> = Vec::with_capacity(n_u_cps);
let mut rows_w: Vec<Vec<f64>> = Vec::with_capacity(n_u_cps);
for ui in 0..n_u_cps {
let c = u_cos[ui];
let s = u_sin[ui];
let uw = u_weights[ui];
let mut row_cp = Vec::with_capacity(n_v_cps);
let mut row_w = Vec::with_capacity(n_v_cps);
for vj in 0..n_v_cps {
let [r, z] = prof_cps[vj];
let pw = prof_ws[vj];
let cp = [
center[0] + r * (c * bu[0] + s * bv[0]) + z * axis_n[0],
center[1] + r * (c * bu[1] + s * bv[1]) + z * axis_n[1],
center[2] + r * (c * bu[2] + s * bv[2]) + z * axis_n[2],
];
row_cp.push(cp);
row_w.push(uw * pw);
}
rows_cp.push(row_cp);
rows_w.push(row_w);
}
Some(neco_nurbs::NurbsSurface3D {
degree_u: 2,
degree_v: v_degree,
control_points: rows_cp,
weights: rows_w,
knots_u,
knots_v,
})
}
_ => None,
}
}
pub fn is_analytical(&self) -> bool {
!matches!(self, Surface::NurbsSurface { .. })
}
pub fn param_range(&self) -> (f64, f64, f64, f64) {
use std::f64::consts::{PI, TAU};
match self {
Surface::Plane { .. } => (0.0, 1.0, 0.0, 1.0),
Surface::Cylinder { axis, .. } => (0.0, TAU, 0.0, vec3::length(*axis)),
Surface::Cone { axis, .. } => (0.0, TAU, 0.0, vec3::length(*axis)),
Surface::Sphere { .. } => (0.0, TAU, 0.0, PI),
Surface::Ellipsoid { .. } => (0.0, TAU, 0.0, PI),
Surface::Torus { .. } => (0.0, TAU, 0.0, TAU),
Surface::SurfaceOfRevolution {
theta_start,
theta_range,
..
} => (*theta_start, *theta_start + *theta_range, 0.0, 1.0),
Surface::SurfaceOfSweep { .. } => (0.0, 1.0, 0.0, 1.0),
Surface::NurbsSurface { data } => {
let (u0, u1) = data.u_range();
let (v0, v1) = data.v_range();
(u0, u1, v0, v1)
}
}
}
pub fn inverse_project(&self, p: &[f64; 3]) -> Option<(f64, f64)> {
match self {
Surface::Plane { origin, normal } => {
let n = vec3::normalized(*normal);
let up = if n[0].abs() < 0.9 {
[1.0, 0.0, 0.0]
} else {
[0.0, 1.0, 0.0]
};
let u_vec = vec3::normalized(vec3::cross(n, up));
let v_vec = vec3::cross(n, u_vec);
let d = vec3::sub(*p, *origin);
Some((vec3::dot(d, u_vec), vec3::dot(d, v_vec)))
}
Surface::Cylinder { origin, axis, .. } => {
let axis_len = vec3::length(*axis);
if axis_len < 1e-30 {
return None;
}
let axis_n = vec3::normalized(*axis);
let (bu, bv) = vec3::orthonormal_basis(*axis);
let q = vec3::sub(*p, *origin);
let h = vec3::dot(q, axis_n).clamp(0.0, axis_len);
let q_perp = vec3::sub(q, vec3::scale(axis_n, h));
let u = vec3::dot(q_perp, bv)
.atan2(vec3::dot(q_perp, bu))
.rem_euclid(std::f64::consts::TAU);
Some((u, h))
}
Surface::Cone { origin, axis, .. } => {
let axis_len = vec3::length(*axis);
if axis_len < 1e-30 {
return None;
}
let axis_n = vec3::normalized(*axis);
let (bu, bv) = vec3::orthonormal_basis(*axis);
let q = vec3::sub(*p, *origin);
let h = vec3::dot(q, axis_n).clamp(0.0, axis_len);
let q_perp = vec3::sub(q, vec3::scale(axis_n, h));
let u = vec3::dot(q_perp, bv)
.atan2(vec3::dot(q_perp, bu))
.rem_euclid(std::f64::consts::TAU);
Some((u, h))
}
Surface::Sphere { center, radius } => {
let d = vec3::sub(*p, *center);
let len = vec3::length(d);
if len < 1e-30 || radius.abs() < 1e-30 {
return None;
}
let u = d[1].atan2(d[0]).rem_euclid(std::f64::consts::TAU);
let cos_v = (d[2] / radius).clamp(-1.0, 1.0);
let v = cos_v.acos();
Some((u, v))
}
Surface::Ellipsoid { center, rx, ry, rz } => {
if rx.abs() < 1e-30 || ry.abs() < 1e-30 || rz.abs() < 1e-30 {
return None;
}
let x = (p[0] - center[0]) / rx;
let y = (p[1] - center[1]) / ry;
let z = (p[2] - center[2]) / rz;
let u = y.atan2(x).rem_euclid(std::f64::consts::TAU);
let v = z.clamp(-1.0, 1.0).acos();
Some((u, v))
}
Surface::Torus {
center,
axis,
major_radius,
..
} => {
let axis_n = vec3::normalized(*axis);
let (bu, bv) = vec3::orthonormal_basis(*axis);
let q = vec3::sub(*p, *center);
let axial = vec3::dot(q, axis_n);
let q_perp = vec3::sub(q, vec3::scale(axis_n, axial));
let r_perp = vec3::length(q_perp);
if r_perp < 1e-30 {
return None;
}
let u = vec3::dot(q_perp, bv)
.atan2(vec3::dot(q_perp, bu))
.rem_euclid(std::f64::consts::TAU);
let v = axial
.atan2(r_perp - major_radius)
.rem_euclid(std::f64::consts::TAU);
Some((u, v))
}
Surface::SurfaceOfRevolution {
center,
axis,
frame_u,
frame_v,
profile_control_points,
profile_weights,
profile_degree,
n_profile_spans,
theta_start,
theta_range,
} => {
let axis_n = vec3::normalized(*axis);
let bu = vec3::normalized(*frame_u);
let bv = vec3::normalized(*frame_v);
let q = vec3::sub(*p, *center);
let theta_raw = vec3::dot(q, bv).atan2(vec3::dot(q, bu));
let theta = if *theta_range >= std::f64::consts::TAU - 1e-12 {
let offset = (theta_raw - theta_start).rem_euclid(std::f64::consts::TAU);
theta_start + offset
} else {
let offset = (theta_raw - theta_start).rem_euclid(std::f64::consts::TAU);
if offset <= *theta_range {
theta_start + offset
} else if offset.min(std::f64::consts::TAU - offset)
<= (offset - theta_range).min(std::f64::consts::TAU - offset + theta_range)
{
*theta_start
} else {
*theta_start + *theta_range
}
};
let radial_dir =
vec3::add(vec3::scale(bu, theta.cos()), vec3::scale(bv, theta.sin()));
let rho = vec3::dot(q, radial_dir);
let z = vec3::dot(q, axis_n);
let v = find_closest_v_on_profile(
rho,
z,
profile_control_points,
profile_weights,
*profile_degree,
*n_profile_spans,
);
Some((theta, v))
}
Surface::SurfaceOfSweep { .. } => {
let (u0, u1, v0, v1) = self.param_range();
Some(grid_project_surface(self, p, u0, u1, v0, v1))
}
Surface::NurbsSurface { data } => {
let (u0, u1) = data.u_range();
let (v0, v1) = data.v_range();
Some(grid_project_surface(self, p, u0, u1, v0, v1))
}
}
}
}
fn grid_project_surface(
surface: &Surface,
p: &[f64; 3],
u0: f64,
u1: f64,
v0: f64,
v1: f64,
) -> (f64, f64) {
let nu = 16;
let nv = 16;
let mut best_u = u0;
let mut best_v = v0;
let mut best_dist = f64::INFINITY;
for iu in 0..=nu {
for iv in 0..=nv {
let u = u0 + (u1 - u0) * (iu as f64 / nu as f64);
let v = v0 + (v1 - v0) * (iv as f64 / nv as f64);
let pt = surface.evaluate(u, v);
let d = vec3::distance(pt, *p);
if d < best_dist {
best_dist = d;
best_u = u;
best_v = v;
}
}
}
let eps = 1e-8;
for _ in 0..20 {
let pt = surface.evaluate(best_u, best_v);
let diff = vec3::sub(pt, *p);
if vec3::length(diff) < 1e-12 {
break;
}
let du = vec3::scale(
vec3::sub(
surface.evaluate(best_u + eps, best_v),
surface.evaluate(best_u - eps, best_v),
),
0.5 / eps,
);
let dv = vec3::scale(
vec3::sub(
surface.evaluate(best_u, best_v + eps),
surface.evaluate(best_u, best_v - eps),
),
0.5 / eps,
);
let a11 = vec3::dot(du, du);
let a12 = vec3::dot(du, dv);
let a22 = vec3::dot(dv, dv);
let b1 = -vec3::dot(du, diff);
let b2 = -vec3::dot(dv, diff);
let det = a11 * a22 - a12 * a12;
if det.abs() < 1e-30 {
break;
}
let delta_u = (a22 * b1 - a12 * b2) / det;
let delta_v = (a11 * b2 - a12 * b1) / det;
best_u = (best_u + delta_u).clamp(u0, u1);
best_v = (best_v + delta_v).clamp(v0, v1);
}
(best_u, best_v)
}
pub fn eval_revolution_profile(
profile_control_points: &[[f64; 2]],
profile_weights: &[f64],
profile_degree: u32,
n_profile_spans: u32,
v: f64,
) -> (f64, f64) {
let n_prof = n_profile_spans as usize;
let p_deg = profile_degree as usize;
let cps_per_span = p_deg + 1;
if n_prof <= 1 {
de_casteljau_rational_2d(profile_control_points, profile_weights, v)
} else {
let v_clamped = v.clamp(0.0, 1.0);
let span_f = v_clamped * n_prof as f64;
let span_idx = (span_f as usize).min(n_prof - 1);
let local_v = span_f - span_idx as f64;
let start = span_idx * cps_per_span;
let end = start + cps_per_span;
de_casteljau_rational_2d(
&profile_control_points[start..end],
&profile_weights[start..end],
local_v,
)
}
}
pub fn find_closest_v_on_profile(
rho: f64,
z: f64,
profile_control_points: &[[f64; 2]],
profile_weights: &[f64],
profile_degree: u32,
n_profile_spans: u32,
) -> f64 {
let n_prof = n_profile_spans as usize;
let n_search = if n_prof <= 1 { 32 } else { n_prof * 16 };
let mut best_v = 0.0;
let mut best_dist_sq = f64::INFINITY;
for i in 0..=n_search {
let v = i as f64 / n_search as f64;
let (r, zz) = eval_revolution_profile(
profile_control_points,
profile_weights,
profile_degree,
n_profile_spans,
v,
);
let dr = rho - r;
let dz = z - zz;
let d2 = dr * dr + dz * dz;
if d2 < best_dist_sq {
best_dist_sq = d2;
best_v = v;
}
}
let eps = 1e-8;
for _ in 0..20 {
let (r, zz) = eval_revolution_profile(
profile_control_points,
profile_weights,
profile_degree,
n_profile_spans,
best_v,
);
let (r_p, zz_p) = eval_revolution_profile(
profile_control_points,
profile_weights,
profile_degree,
n_profile_spans,
(best_v + eps).min(1.0),
);
let (r_m, zz_m) = eval_revolution_profile(
profile_control_points,
profile_weights,
profile_degree,
n_profile_spans,
(best_v - eps).max(0.0),
);
let dr_dv = (r_p - r_m) / (2.0 * eps);
let dz_dv = (zz_p - zz_m) / (2.0 * eps);
let residual_r = rho - r;
let residual_z = z - zz;
let f_val = residual_r * (-dr_dv) + residual_z * (-dz_dv);
let df_val = dr_dv * dr_dv + dz_dv * dz_dv;
if df_val.abs() < 1e-20 {
break;
}
let dv = f_val / df_val;
best_v = (best_v - dv).clamp(0.0, 1.0);
if dv.abs() < 1e-12 {
break;
}
}
best_v
}
#[derive(Debug, Clone)]
pub enum Curve3D {
Line {
start: [f64; 3],
end: [f64; 3],
},
Arc {
center: [f64; 3],
axis: [f64; 3],
start: [f64; 3],
end: [f64; 3],
radius: f64,
},
Ellipse {
center: [f64; 3],
axis_u: [f64; 3],
axis_v: [f64; 3],
t_start: f64,
t_end: f64,
},
NurbsCurve3D {
degree: usize,
control_points: Vec<[f64; 3]>,
weights: Vec<f64>,
knots: Vec<f64>,
},
}
impl Curve3D {
pub fn evaluate(&self, t: f64) -> [f64; 3] {
match self {
Curve3D::Line { start, end } => {
vec3::add(*start, vec3::scale(vec3::sub(*end, *start), t))
}
Curve3D::Arc {
center,
axis,
start,
radius,
..
} => {
let r_vec = vec3::sub(*start, *center);
let tangent = vec3::scale(vec3::normalized(vec3::cross(*axis, r_vec)), *radius);
let cos_t = t.cos();
let sin_t = t.sin();
vec3::add(
vec3::add(*center, vec3::scale(r_vec, cos_t)),
vec3::scale(tangent, sin_t),
)
}
Curve3D::Ellipse {
center,
axis_u,
axis_v,
..
} => vec3::add(
vec3::add(*center, vec3::scale(*axis_u, t.cos())),
vec3::scale(*axis_v, t.sin()),
),
Curve3D::NurbsCurve3D {
degree,
control_points,
weights,
knots,
} => {
let n = control_points.len();
let p = *degree;
let t_min = knots[p];
let t_max = knots[n];
let t = t.clamp(t_min, t_max);
let k = {
if t >= knots[n] {
n - 1
} else {
let mut lo = p;
let mut hi = n;
while lo < hi {
let mid = (lo + hi) / 2;
if t < knots[mid] {
hi = mid;
} else {
lo = mid + 1;
}
}
lo - 1
}
};
let mut hx = vec![0.0; p + 1];
let mut hy = vec![0.0; p + 1];
let mut hz = vec![0.0; p + 1];
let mut hw = vec![0.0; p + 1];
for j in 0..=p {
let idx = k - p + j;
let w = weights[idx];
hx[j] = control_points[idx][0] * w;
hy[j] = control_points[idx][1] * w;
hz[j] = control_points[idx][2] * w;
hw[j] = w;
}
for r in 1..=p {
for j in (r..=p).rev() {
let left = k + j - p;
let right = k + 1 + j - r;
let denom = knots[right] - knots[left];
if denom.abs() < 1e-30 {
continue;
}
let alpha = (t - knots[left]) / denom;
hx[j] = (1.0 - alpha) * hx[j - 1] + alpha * hx[j];
hy[j] = (1.0 - alpha) * hy[j - 1] + alpha * hy[j];
hz[j] = (1.0 - alpha) * hz[j - 1] + alpha * hz[j];
hw[j] = (1.0 - alpha) * hw[j - 1] + alpha * hw[j];
}
}
let w = hw[p];
if w.abs() < 1e-30 {
control_points.last().copied().unwrap_or([0.0; 3])
} else {
[hx[p] / w, hy[p] / w, hz[p] / w]
}
}
}
}
pub fn midpoint(&self) -> [f64; 3] {
match self {
Curve3D::Line { start, end } => vec3::scale(vec3::add(*start, *end), 0.5),
Curve3D::Arc { .. } => self.evaluate(0.5),
Curve3D::Ellipse { t_start, t_end, .. } => self.evaluate((t_start + t_end) * 0.5),
Curve3D::NurbsCurve3D { degree, knots, .. } => {
let p = *degree;
let n = knots.len() - p - 1;
let mid = (knots[p] + knots[n]) * 0.5;
self.evaluate(mid)
}
}
}
pub fn param_range(&self) -> (f64, f64) {
match self {
Curve3D::Line { .. } => (0.0, 1.0),
Curve3D::Arc {
start,
end,
center,
axis,
radius,
} => {
let r_start = vec3::sub(*start, *center);
let tangent = vec3::scale(vec3::normalized(vec3::cross(*axis, r_start)), *radius);
let r_end = vec3::sub(*end, *center);
let r_start_len = vec3::length(r_start);
let r_end_len = vec3::length(r_end);
let cos_angle = vec3::dot(r_start, r_end) / (r_start_len * r_end_len);
let sin_angle = vec3::dot(tangent, r_end) / (vec3::length(tangent) * r_end_len);
let angle = sin_angle.atan2(cos_angle);
let angle = if angle.abs() < 1e-10 && vec3::length(vec3::sub(*start, *end)) < 1e-10
{
std::f64::consts::TAU
} else {
angle.max(0.0)
};
(0.0, angle)
}
Curve3D::Ellipse { t_start, t_end, .. } => (*t_start, *t_end),
Curve3D::NurbsCurve3D { degree, knots, .. } => {
let p = *degree;
let n = knots.len() - p - 1;
(knots[p], knots[n])
}
}
}
pub fn to_polyline(&self, chord_tol: f64) -> Vec<[f64; 3]> {
match self {
Curve3D::Line { start, end } => vec![*start, *end],
_ => {
let (t0, t1) = self.param_range();
let mut points = vec![self.evaluate(t0)];
self.adaptive_sample(t0, t1, chord_tol, &mut points);
points
}
}
}
fn adaptive_sample(&self, t0: f64, t1: f64, tol: f64, points: &mut Vec<[f64; 3]>) {
if (t1 - t0).abs() < 1e-10 {
points.push(self.evaluate(t1));
return;
}
let mid_t = (t0 + t1) * 0.5;
let p0 = self.evaluate(t0);
let p1 = self.evaluate(t1);
let mid_curve = self.evaluate(mid_t);
let mid_chord = vec3::scale(vec3::add(p0, p1), 0.5);
let chord_height = vec3::length(vec3::sub(mid_curve, mid_chord));
if chord_height > tol {
self.adaptive_sample(t0, mid_t, tol, points);
self.adaptive_sample(mid_t, t1, tol, points);
} else {
points.push(self.evaluate(t1));
}
}
}
#[derive(Debug, Clone, Copy)]
pub struct EdgeRef {
pub edge_id: EdgeId,
pub forward: bool,
}
#[derive(Debug, Clone)]
pub struct Edge {
pub v_start: VertexId,
pub v_end: VertexId,
pub curve: Curve3D,
}
#[derive(Debug, Clone)]
pub struct Face {
pub loop_edges: Vec<EdgeRef>,
pub surface: Surface,
pub orientation_reversed: bool,
}
#[derive(Clone, Debug)]
pub struct SubFace {
pub surface: Surface,
pub polygon: Vec<[f64; 3]>,
pub candidate_curves: Vec<Curve3D>,
pub flipped: bool,
pub source_shell: usize,
pub source_face: usize,
}
impl SubFace {
pub fn centroid(&self) -> [f64; 3] {
let n = self.polygon.len() as f64;
let (sx, sy, sz) = self.polygon.iter().fold((0.0, 0.0, 0.0), |(x, y, z), p| {
(x + p[0], y + p[1], z + p[2])
});
[sx / n, sy / n, sz / n]
}
}
#[derive(Debug, Clone, Default)]
pub struct Shell {
pub vertices: Vec<[f64; 3]>,
pub edges: Vec<Edge>,
pub faces: Vec<Face>,
}
impl Shell {
pub fn new() -> Self {
Self {
vertices: Vec::new(),
edges: Vec::new(),
faces: Vec::new(),
}
}
pub fn add_vertex(&mut self, p: [f64; 3]) -> VertexId {
for (i, v) in self.vertices.iter().enumerate() {
if vec3::distance(*v, p) < VERTEX_TOL {
return i;
}
}
let id = self.vertices.len();
self.vertices.push(p);
id
}
pub fn add_edge(&mut self, v_start: VertexId, v_end: VertexId, curve: Curve3D) -> EdgeId {
let id = self.edges.len();
self.edges.push(Edge {
v_start,
v_end,
curve,
});
id
}
pub fn bounding_box(&self) -> ([f64; 3], [f64; 3]) {
let mut min = [f64::INFINITY; 3];
let mut max = [f64::NEG_INFINITY; 3];
let mut extend = |p: [f64; 3]| {
for i in 0..3 {
min[i] = min[i].min(p[i]);
max[i] = max[i].max(p[i]);
}
};
for v in &self.vertices {
extend(*v);
}
for edge in &self.edges {
match &edge.curve {
Curve3D::Line { .. } => {}
curve => {
let polyline = curve.to_polyline(0.01);
for p in &polyline {
extend(*p);
}
}
}
}
(min, max)
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn sphere_evaluate_poles() {
let s = Surface::Sphere {
center: [0.0, 0.0, 0.0],
radius: 1.0,
};
let north = s.evaluate(0.0, 0.0);
assert!((north[2] - 1.0).abs() < 1e-10);
let south = s.evaluate(0.0, std::f64::consts::PI);
assert!((south[2] + 1.0).abs() < 1e-10);
}
#[test]
fn sphere_normal_outward() {
let s = Surface::Sphere {
center: [1.0, 2.0, 3.0],
radius: 2.0,
};
let n = s.normal_at(0.0, std::f64::consts::FRAC_PI_2);
let p = s.evaluate(0.0, std::f64::consts::FRAC_PI_2);
let expected_x = (p[0] - 1.0) / 2.0;
assert!((n[0] - expected_x).abs() < 1e-10);
}
#[test]
fn cylinder_evaluate() {
let s = Surface::Cylinder {
origin: [0.0, 0.0, 0.0],
axis: [0.0, 0.0, 2.0],
radius: 1.0,
};
let p = s.evaluate(0.0, 1.0);
assert!((p[2] - 1.0).abs() < 1e-10);
let r = (p[0] * p[0] + p[1] * p[1]).sqrt();
assert!((r - 1.0).abs() < 1e-10);
}
fn make_multispan_revolution() -> Surface {
let profile_control_points = vec![
[1.0, 0.0],
[1.0, 1.0],
[1.0, 2.0],
[2.0, 3.0],
[2.0, 3.0],
[2.0, 4.0],
[2.0, 5.0],
[1.0, 6.0],
];
let profile_weights = vec![1.0; 8];
Surface::SurfaceOfRevolution {
center: [0.0, 0.0, 0.0],
axis: [0.0, 0.0, 1.0],
frame_u: [1.0, 0.0, 0.0],
frame_v: [0.0, 1.0, 0.0],
profile_control_points,
profile_weights,
profile_degree: 3,
n_profile_spans: 2,
theta_start: 0.0,
theta_range: std::f64::consts::TAU,
}
}
#[test]
fn revolution_multispan_evaluate_endpoints() {
let s = make_multispan_revolution();
let p0 = s.evaluate(0.0, 0.0);
assert!((p0[0] - 1.0).abs() < 1e-10);
assert!(p0[2].abs() < 1e-10);
let p1 = s.evaluate(0.0, 1.0);
assert!((p1[0] - 1.0).abs() < 1e-10);
assert!((p1[2] - 6.0).abs() < 1e-10);
}
#[test]
fn revolution_multispan_evaluate_midpoint() {
let s = make_multispan_revolution();
let pm = s.evaluate(0.0, 0.5);
assert!((pm[0] - 2.0).abs() < 1e-10);
assert!((pm[2] - 3.0).abs() < 1e-10);
}
#[test]
fn revolution_multispan_normal_at() {
let s = make_multispan_revolution();
let n = s.normal_at(0.0, 0.5);
let len = vec3::length(n);
assert!((len - 1.0).abs() < 1e-6);
}
#[test]
fn param_range_sphere() {
let s = Surface::Sphere {
center: [0.0, 0.0, 0.0],
radius: 1.0,
};
let (u0, u1, v0, v1) = s.param_range();
assert!(u0.abs() < 1e-10);
assert!((u1 - std::f64::consts::TAU).abs() < 1e-10);
assert!(v0.abs() < 1e-10);
assert!((v1 - std::f64::consts::PI).abs() < 1e-10);
}
#[test]
fn sphere_inverse_project_round_trip() {
let s = Surface::Sphere {
center: [1.0, 2.0, 3.0],
radius: 2.0,
};
let p = s.evaluate(0.7, 1.2);
let (u, v) = s.inverse_project(&p).unwrap();
let q = s.evaluate(u, v);
assert!(vec3::distance(p, q) < 1e-10);
}
#[test]
fn torus_inverse_project_round_trip() {
let s = Surface::Torus {
center: [0.0, 0.0, 0.0],
axis: [0.0, 0.0, 1.0],
major_radius: 3.0,
minor_radius: 1.0,
};
let p = s.evaluate(1.1, 2.2);
let (u, v) = s.inverse_project(&p).unwrap();
let q = s.evaluate(u, v);
assert!(vec3::distance(p, q) < 1e-10);
}
#[test]
fn revolution_inverse_project_round_trip() {
let s = make_multispan_revolution();
let p = s.evaluate(1.3, 0.4);
let (u, v) = s.inverse_project(&p).unwrap();
let q = s.evaluate(u, v);
assert!(vec3::distance(p, q) < 1e-6);
}
}