diff --git a/Cargo.lock b/Cargo.lock index f8c06aa04a..8f09808c75 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -863,6 +863,20 @@ dependencies = [ "unicode-segmentation", ] +[[package]] +name = "convex-hull-nodes" +version = "0.1.0" +dependencies = [ + "core-types", + "dyn-any", + "glam", + "graphic-types", + "kurbo", + "node-macro", + "poly-cool", + "vector-types", +] + [[package]] name = "cookie" version = "0.18.1" @@ -1223,7 +1237,7 @@ dependencies = [ "libc", "option-ext", "redox_users", - "windows-sys 0.59.0", + "windows-sys 0.61.2", ] [[package]] @@ -1474,7 +1488,7 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "778e2ac28f6c47af28e4907f13ffd1e1ddbd400980a9abd7c8df189bf578a5ad" dependencies = [ "libc", - "windows-sys 0.59.0", + "windows-sys 0.60.2", ] [[package]] @@ -2115,6 +2129,7 @@ dependencies = [ "blending-nodes", "brush-nodes", "bytemuck", + "convex-hull-nodes", "core-types", "dyn-any", "glam", @@ -3950,7 +3965,7 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "7d8fae84b431384b68627d0f9b3b1245fcf9f46f6c0e3dc902e9dce64edd1967" dependencies = [ "libc", - "windows-sys 0.45.0", + "windows-sys 0.61.2", ] [[package]] @@ -4306,6 +4321,15 @@ version = "0.4.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "2f3a9f18d041e6d0e102a0a46750538147e5e8992d3b4873aaafee2520b00ce3" +[[package]] +name = "poly-cool" +version = "0.3.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cb79f376772fbca123c950f4c2a74558bae8f8fd7d35037c093d5d53d45b46c2" +dependencies = [ + "arrayvec", +] + [[package]] name = "polycool" version = "0.4.0" @@ -4510,7 +4534,7 @@ dependencies = [ "once_cell", "socket2", "tracing", - "windows-sys 0.59.0", + "windows-sys 0.60.2", ] [[package]] @@ -4986,7 +5010,7 @@ dependencies = [ "errno", "libc", "linux-raw-sys", - "windows-sys 0.59.0", + "windows-sys 0.60.2", ] [[package]] @@ -5054,7 +5078,7 @@ dependencies = [ "security-framework", "security-framework-sys", "webpki-root-certs", - "windows-sys 0.59.0", + "windows-sys 0.61.2", ] [[package]] @@ -5787,7 +5811,7 @@ dependencies = [ "getrandom 0.3.3", "once_cell", "rustix", - "windows-sys 0.59.0", + "windows-sys 0.61.2", ] [[package]] @@ -7059,7 +7083,7 @@ version = "0.1.10" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "0978bf7171b3d90bac376700cb56d606feb40f251a475a5d6634613564460b22" dependencies = [ - "windows-sys 0.59.0", + "windows-sys 0.60.2", ] [[package]] diff --git a/Cargo.toml b/Cargo.toml index 3cc49a3614..3decdb5abc 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -164,6 +164,7 @@ url = "2.5" tokio = { version = "1.29", features = ["fs", "macros", "io-std", "rt", "rt-multi-thread"] } # Linebender ecosystem (BEGIN) kurbo = { version = "0.13", features = ["serde"] } +poly-cool = "0.3" vello = "0.9" vello_encoding = "0.9" resvg = "0.47" diff --git a/node-graph/nodes/convex_hull/Cargo.toml b/node-graph/nodes/convex_hull/Cargo.toml index 9658e54644..c571597e6a 100644 --- a/node-graph/nodes/convex_hull/Cargo.toml +++ b/node-graph/nodes/convex_hull/Cargo.toml @@ -12,11 +12,7 @@ dyn-any = { workspace = true } core-types = { workspace = true } graphic-types = { workspace = true } node-macro = { workspace = true } -glam = { workspace = true } -specta = { workspace = true } -log = { workspace = true } -path-bool = { workspace = true } -serde = { workspace = true } vector-types = { workspace = true } +glam = { workspace = true } kurbo = { workspace = true } -convex_hull = { path = "../../../../convex_hull" } +poly-cool = { workspace = true } diff --git a/node-graph/nodes/convex_hull/src/hull.rs b/node-graph/nodes/convex_hull/src/hull.rs new file mode 100644 index 0000000000..abf443c5df --- /dev/null +++ b/node-graph/nodes/convex_hull/src/hull.rs @@ -0,0 +1,1192 @@ +// Exact convex hull of a closed Bézier spline. +// +// Architecture: support-function envelope. Every boundary feature (convex +// arc, corner, line segment) contributes a support function h(θ) = ⟨p(θ), n(θ)⟩ +// over an interval of outward-normal angles θ. The hull is the upper envelope +// of these functions over [0, 2π); envelope pieces are hull arcs and envelope +// transitions are bitangent line segments. +// +// Each curved piece contributes *two* candidates (one per normal side), so no +// global orientation or convexity classification is needed — wrong-side +// candidates simply never win the envelope. This makes the algorithm robust +// to input orientation and self-intersecting curves. + +use kurbo::{CubicBez, ParamCurve, ParamCurveDeriv, Point, Vec2}; +use poly_cool::PolyDyn; +use std::f64::consts::PI; + +const TOL: f64 = 1e-10; +/// Angular tolerance (radians) for interval membership and breakpoint dedup. +const ANG_EPS: f64 = 1e-9; +/// Angular tolerance for matching normals of a bitangent solution. +const ANG_MATCH: f64 = 1e-6; + +const TAU: f64 = 2.0 * PI; + +// ─── Public Types ─── + +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum PieceKind { + /// Monotone-curvature curved arc. + Arc, + /// Straight line segment. + Line, + /// Zero-length corner at a tangent-discontinuous junction. + Corner, +} + +/// A boundary piece produced by decomposition: monotone arc, line, or corner. +#[derive(Clone, Debug)] +pub struct MonotoneArc { + pub bezier: CubicBez, + pub kind: PieceKind, + /// Tangent angle at start/end, unwrapped along the piece + /// (theta_end - theta_start = signed tangent turn). + pub theta_start: f64, + pub theta_end: f64, + /// Index of the original spline segment this piece came from. + pub original_segment: usize, + /// Parameter range within the original segment. + pub original_t_end: f64, +} + +#[derive(Clone, Debug)] +pub enum HullSegment { + /// A portion of piece `arc_index` on the hull. `t_start` may exceed + /// `t_end` when the hull traverses the piece against its parameterization. + Arc { + arc_index: usize, + t_start: f64, + t_end: f64, + }, + /// A bitangent line segment bridging between two hull contacts. + Line { start: Point, end: Point }, +} + +// ─── Basic Geometry Helpers ─── + +fn is_collinear(cb: &CubicBez) -> bool { + let d1 = cb.p1 - cb.p0; + let d2 = cb.p2 - cb.p0; + let d3 = cb.p3 - cb.p0; + let len = d1.hypot().max(d2.hypot()).max(d3.hypot()).max(1e-15); + (d1.x * d2.y - d1.y * d2.x).abs() < 1e-8 * len + && (d1.x * d3.y - d1.y * d3.x).abs() < 1e-8 * len +} + +fn is_degenerate_point(cb: &CubicBez) -> bool { + let tol = 1e-12; + (cb.p0 - cb.p3).hypot() < tol + && (cb.p0 - cb.p1).hypot() < tol + && (cb.p0 - cb.p2).hypot() < tol +} + +/// Tangent direction at the start of a cubic, robust to degenerate +/// parameterizations (p0 == p1 etc.): first nonzero control-point difference. +fn start_tangent_dir(cb: &CubicBez) -> Vec2 { + for v in [cb.p1 - cb.p0, cb.p2 - cb.p0, cb.p3 - cb.p0] { + if v.hypot() > 1e-12 { + return v; + } + } + Vec2::new(1.0, 0.0) +} + +/// Tangent direction at the end of a cubic, robust to degenerate +/// parameterizations. +fn end_tangent_dir(cb: &CubicBez) -> Vec2 { + for v in [cb.p3 - cb.p2, cb.p3 - cb.p1, cb.p3 - cb.p0] { + if v.hypot() > 1e-12 { + return v; + } + } + Vec2::new(1.0, 0.0) +} + +/// Tangent direction at parameter t, falling back to control-point geometry +/// near degenerate endpoints. +fn tangent_dir_at(cb: &CubicBez, t: f64) -> Vec2 { + let d = cb.deriv().eval(t).to_vec2(); + if d.hypot() > 1e-9 * chord_scale(cb) { + return d; + } + if t < 0.5 { + start_tangent_dir(cb) + } else { + end_tangent_dir(cb) + } +} + +fn chord_scale(cb: &CubicBez) -> f64 { + (cb.p1 - cb.p0) + .hypot() + .max((cb.p2 - cb.p0).hypot()) + .max((cb.p3 - cb.p0).hypot()) + .max(1e-15) +} + +fn unwrap_angle(angle: f64, reference: f64) -> f64 { + let mut a = angle; + while a - reference > PI { + a -= TAU; + } + while a - reference < -PI { + a += TAU; + } + a +} + +fn normalize_angle(a: f64) -> f64 { + let mut x = a % TAU; + if x < 0.0 { + x += TAU; + } + x +} + +/// Circular forward distance from `from` to `to` in [0, 2π). +fn ang_forward(from: f64, to: f64) -> f64 { + normalize_angle(to - from) +} + +// ─── Decomposition ─── + +/// Split multiple closed splines (loops) into pieces. Corner insertion wraps +/// within each loop; `original_segment` indexes into the flattened segment +/// list across all loops. +pub fn split_loops_into_arcs(loops: &[&[CubicBez]]) -> Vec { + let mut all = Vec::new(); + let mut seg_offset = 0; + for segments in loops { + decompose_loop(segments, seg_offset, &mut all); + seg_offset += segments.len(); + } + all +} + +fn decompose_loop(segments: &[CubicBez], seg_offset: usize, out: &mut Vec) { + let mut pieces: Vec = Vec::new(); + + for (rel_idx, cb) in segments.iter().enumerate() { + let seg_idx = seg_offset + rel_idx; + if is_degenerate_point(cb) { + continue; + } + if is_collinear(cb) { + let dir = (cb.p3 - cb.p0).atan2(); + pieces.push(MonotoneArc { + bezier: *cb, + kind: PieceKind::Line, + theta_start: dir, + theta_end: dir, + original_segment: seg_idx, + original_t_end: 1.0, + }); + continue; + } + + let mut cuts: Vec = vec![0.0]; + let mut infl: Vec = cb + .inflections() + .iter() + .copied() + .filter(|&t| t > 1e-6 && t < 1.0 - 1e-6) + .collect(); + infl.sort_by(|a, b| a.partial_cmp(b).unwrap()); + cuts.extend(infl); + cuts.push(1.0); + + for w in cuts.windows(2) { + let (t0, t1) = (w[0], w[1]); + if t1 - t0 < TOL { + continue; + } + let sub = cb.subsegment(t0..t1); + if is_degenerate_point(&sub) { + continue; + } + if is_collinear(&sub) { + let dir = (sub.p3 - sub.p0).atan2(); + pieces.push(MonotoneArc { + bezier: sub, + kind: PieceKind::Line, + theta_start: dir, + theta_end: dir, + original_segment: seg_idx, + original_t_end: t1, + }); + continue; + } + + // Unwrap the tangent angle along the piece via interior samples so + // theta_end - theta_start is the true signed turn. + let mut theta = start_tangent_dir(&sub).atan2(); + let theta_s = theta; + for k in 1..=8 { + let t = k as f64 / 8.0; + let d = tangent_dir_at(&sub, t); + theta = unwrap_angle(d.atan2(), theta); + } + let theta_e = unwrap_angle(end_tangent_dir(&sub).atan2(), theta); + + pieces.push(MonotoneArc { + bezier: sub, + kind: PieceKind::Arc, + theta_start: theta_s, + theta_end: theta_e, + original_segment: seg_idx, + original_t_end: t1, + }); + } + } + + // Insert corner pieces at tangent-discontinuous junctions (wrapping + // within this loop). + let n = pieces.len(); + for i in 0..n { + out.push(pieces[i].clone()); + let next = &pieces[(i + 1) % n]; + let theta_out = end_tangent_dir(&pieces[i].bezier).atan2(); + let theta_in = start_tangent_dir(&next.bezier).atan2(); + let mut gap = theta_in - theta_out; + while gap > PI { + gap -= TAU; + } + while gap <= -PI { + gap += TAU; + } + if gap.abs() < 1e-8 { + continue; + } + let vertex = pieces[i].bezier.p3; + out.push(MonotoneArc { + bezier: CubicBez::new(vertex, vertex, vertex, vertex), + kind: PieceKind::Corner, + theta_start: theta_out, + theta_end: theta_out + gap, + original_segment: pieces[i].original_segment, + original_t_end: pieces[i].original_t_end, + }); + } +} + +// ─── Support Candidates ─── + +/// One support-function candidate: a piece viewed with one choice of outward +/// normal side. +#[derive(Clone, Debug)] +struct Candidate { + piece: usize, + kind: PieceKind, + /// +1: normal = tangent rotated -90°; -1: normal = tangent rotated +90°. + side: f64, + /// Interval of normal angles covered, as (start ∈ [0,2π), length ≥ 0). + /// Corners cover the full circle (len = 2π). Lines have len = 0. + ang_start: f64, + ang_len: f64, + /// Curve parameter at interval start / end (arcs only). + t_at_start: f64, + t_at_end: f64, +} + +impl Candidate { + fn contains(&self, theta: f64, eps: f64) -> bool { + if self.kind == PieceKind::Corner { + return true; + } + let off = ang_forward(self.ang_start, theta); + off <= self.ang_len + eps || off >= TAU - eps + } +} + +fn normal_angle(tangent: Vec2, side: f64) -> f64 { + Vec2::new(side * tangent.y, -side * tangent.x).atan2() +} + +fn build_candidates(pieces: &[MonotoneArc]) -> Vec { + let mut cands = Vec::new(); + for (i, p) in pieces.iter().enumerate() { + match p.kind { + PieceKind::Corner => { + cands.push(Candidate { + piece: i, + kind: PieceKind::Corner, + side: 1.0, + ang_start: 0.0, + ang_len: TAU, + t_at_start: 0.0, + t_at_end: 0.0, + }); + } + PieceKind::Line => { + for side in [1.0, -1.0] { + let dir = p.bezier.p3 - p.bezier.p0; + cands.push(Candidate { + piece: i, + kind: PieceKind::Line, + side, + ang_start: normalize_angle(normal_angle(dir, side)), + ang_len: 0.0, + t_at_start: 0.0, + t_at_end: 1.0, + }); + } + } + PieceKind::Arc => { + let turn = p.theta_end - p.theta_start; + for side in [1.0, -1.0] { + // Normal angle at t=0 / t=1 for this side. + let n0 = p.theta_start - side * PI / 2.0; + let n1 = p.theta_end - side * PI / 2.0; + let (start, len, t_s, t_e) = if turn >= 0.0 { + (n0, turn, 0.0, 1.0) + } else { + (n1, -turn, 1.0, 0.0) + }; + cands.push(Candidate { + piece: i, + kind: PieceKind::Arc, + side, + ang_start: normalize_angle(start), + ang_len: len.min(TAU), + t_at_start: t_s, + t_at_end: t_e, + }); + } + } + } + } + cands +} + +/// Do two candidates' angular intervals overlap on the circle (open overlap)? +fn intervals_overlap(a: &Candidate, b: &Candidate) -> bool { + if a.kind == PieceKind::Corner || b.kind == PieceKind::Corner { + return true; + } + if a.ang_len >= TAU - ANG_EPS || b.ang_len >= TAU - ANG_EPS { + return true; + } + let off = ang_forward(a.ang_start, b.ang_start); + off < a.ang_len + ANG_EPS || TAU - off < b.ang_len + ANG_EPS +} + +// ─── Support Evaluation ─── + +/// Solve a quadratic a·t² + b·t + c = 0, returning real roots. +fn solve_quadratic(a: f64, b: f64, c: f64, scale: f64) -> Vec { + if a.abs() < 1e-14 * scale { + if b.abs() < 1e-14 * scale { + return vec![]; + } + return vec![-c / b]; + } + let disc = b * b - 4.0 * a * c; + if disc < 0.0 { + return vec![]; + } + let sq = disc.sqrt(); + // Numerically stable form + let q = -0.5 * (b + b.signum() * sq); + let mut roots = vec![q / a]; + if q.abs() > 1e-300 { + roots.push(c / q); + } else { + roots.push(if a != 0.0 { -b / a - roots[0] } else { roots[0] }); + } + roots +} + +/// Power basis coefficients for γ'(t) = c0 + c1·t + c2·t². +fn deriv_power_basis(cb: &CubicBez) -> [Vec2; 3] { + let p0 = cb.p0.to_vec2(); + let p1 = cb.p1.to_vec2(); + let p2 = cb.p2.to_vec2(); + let p3 = cb.p3.to_vec2(); + let q0 = 3.0 * (p1 - p0); + let q1 = 3.0 * (p2 - p1); + let q2 = 3.0 * (p3 - p2); + [q0, 2.0 * (q1 - q0), q0 - 2.0 * q1 + q2] +} + +/// Power basis coefficients for γ(t) = d0 + d1·t + d2·t² + d3·t³. +fn curve_power_basis(cb: &CubicBez) -> [Vec2; 4] { + let p0 = cb.p0.to_vec2(); + let p1 = cb.p1.to_vec2(); + let p2 = cb.p2.to_vec2(); + let p3 = cb.p3.to_vec2(); + [ + p0, + 3.0 * (p1 - p0), + 3.0 * (p0 - 2.0 * p1 + p2), + -p0 + 3.0 * p1 - 3.0 * p2 + p3, + ] +} + +/// Invert the Gauss map of an arc candidate: parameter t whose outward normal +/// (for this candidate's side) equals theta. +fn arc_param_at_normal(cand: &Candidate, piece: &MonotoneArc, theta: f64) -> f64 { + let off = ang_forward(cand.ang_start, theta); + let off = if off > TAU - ANG_EPS { 0.0 } else { off }; + if off < ANG_EPS { + return cand.t_at_start; + } + if off > cand.ang_len - ANG_EPS { + return cand.t_at_end; + } + + let cb = &piece.bezier; + // Tangent direction required at the solution. + let u = Vec2::new(-theta.sin() * cand.side, theta.cos() * cand.side); + let d = deriv_power_basis(cb); + let scale = chord_scale(cb); + let roots = solve_quadratic(d[2].cross(u), d[1].cross(u), d[0].cross(u), scale); + + let expected = cand.t_at_start + (off / cand.ang_len) * (cand.t_at_end - cand.t_at_start); + let mut best: Option = None; + for r in roots { + if !(-1e-9..=1.0 + 1e-9).contains(&r) { + continue; + } + let rc = r.clamp(0.0, 1.0); + let dv = d[0] + rc * d[1] + rc * rc * d[2]; + if dv.dot(u) <= 0.0 { + continue; + } + match best { + Some(b) if (b - expected).abs() <= (rc - expected).abs() => {} + _ => best = Some(rc), + } + } + best.unwrap_or(expected.clamp(0.0, 1.0)) +} + +/// Support value and contact point of a candidate at normal angle theta. +/// Assumes `cand.contains(theta)`. +fn support_at(cand: &Candidate, piece: &MonotoneArc, theta: f64) -> (f64, Point, f64) { + let n = Vec2::new(theta.cos(), theta.sin()); + match cand.kind { + PieceKind::Corner => { + let p = piece.bezier.p0; + (p.to_vec2().dot(n), p, 0.0) + } + PieceKind::Line => { + let p = piece.bezier.p0; + (p.to_vec2().dot(n), p, 0.0) + } + PieceKind::Arc => { + let t = arc_param_at_normal(cand, piece, theta); + let p = piece.bezier.eval(t); + (p.to_vec2().dot(n), p, t) + } + } +} + +// ─── Polynomial Utilities (for crossings) ─── + +fn poly_mul(a: &[f64], b: &[f64]) -> Vec { + if a.is_empty() || b.is_empty() { + return vec![]; + } + let mut result = vec![0.0; a.len() + b.len() - 1]; + for (i, &ai) in a.iter().enumerate() { + for (j, &bj) in b.iter().enumerate() { + result[i + j] += ai * bj; + } + } + result +} + +fn poly_add(a: &[f64], b: &[f64]) -> Vec { + let len = a.len().max(b.len()); + let mut result = vec![0.0; len]; + for (i, &v) in a.iter().enumerate() { + result[i] += v; + } + for (i, &v) in b.iter().enumerate() { + result[i] += v; + } + result +} + +fn poly_sub(a: &[f64], b: &[f64]) -> Vec { + let len = a.len().max(b.len()); + let mut result = vec![0.0; len]; + for (i, &v) in a.iter().enumerate() { + result[i] += v; + } + for (i, &v) in b.iter().enumerate() { + result[i] -= v; + } + result +} + +fn poly_scale(a: &[f64], s: f64) -> Vec { + a.iter().map(|&c| c * s).collect() +} + +fn is_zero_poly(p: &[f64]) -> bool { + p.is_empty() || p.iter().all(|&c| c.abs() < 1e-20) +} + +fn trim_poly(p: &[f64]) -> &[f64] { + let mut len = p.len(); + while len > 1 && p[len - 1].abs() < 1e-20 { + len -= 1; + } + &p[..len] +} + +fn poly_div_exact(num: &[f64], den: &[f64]) -> Vec { + if is_zero_poly(num) { + return vec![0.0]; + } + let den_trimmed = trim_poly(den); + let num_trimmed = trim_poly(num); + if den_trimmed.len() == 1 { + return poly_scale(num_trimmed, 1.0 / den_trimmed[0]); + } + if num_trimmed.len() < den_trimmed.len() { + return vec![0.0]; + } + let mut remainder = num_trimmed.to_vec(); + let mut quotient = vec![0.0; remainder.len() - den_trimmed.len() + 1]; + let lead_den = *den_trimmed.last().unwrap(); + for i in (0..quotient.len()).rev() { + let idx = i + den_trimmed.len() - 1; + let coeff = remainder[idx] / lead_den; + quotient[i] = coeff; + for (j, &d) in den_trimmed.iter().enumerate() { + remainder[i + j] -= coeff * d; + } + } + trim_poly("ient).to_vec() +} + +fn eval_poly(coeffs: &[f64], t: f64) -> f64 { + if coeffs.is_empty() { + return 0.0; + } + let mut result = coeffs[coeffs.len() - 1]; + for i in (0..coeffs.len() - 1).rev() { + result = result * t + coeffs[i]; + } + result +} + +fn eval_bivariate(eq: &[Vec], t1: f64, t2: f64) -> f64 { + let mut result = 0.0; + let mut t2_pow = 1.0; + for coeffs_t1 in eq { + result += eval_poly(coeffs_t1, t1) * t2_pow; + t2_pow *= t2; + } + result +} + +/// Bivariate system for the bitangent conditions between two cubics: +/// eq1: γ'_i(t1) × γ'_j(t2) = 0 (tangents parallel) +/// eq2: γ'_i(t1) × (γ_j(t2) - γ_i(t1)) = 0 (chord aligned with tangent) +/// Each equation is a vector of polynomials in t1, indexed by power of t2. +fn build_bitangent_system(arc_i: &CubicBez, arc_j: &CubicBez) -> (Vec>, Vec>) { + let di = deriv_power_basis(arc_i); + let dj = deriv_power_basis(arc_j); + let gi = curve_power_basis(arc_i); + let gj = curve_power_basis(arc_j); + + let dix: Vec = di.iter().map(|v| v.x).collect(); + let diy: Vec = di.iter().map(|v| v.y).collect(); + let djx: Vec = dj.iter().map(|v| v.x).collect(); + let djy: Vec = dj.iter().map(|v| v.y).collect(); + let gix: Vec = gi.iter().map(|v| v.x).collect(); + let giy: Vec = gi.iter().map(|v| v.y).collect(); + let gjx: Vec = gj.iter().map(|v| v.x).collect(); + let gjy: Vec = gj.iter().map(|v| v.y).collect(); + + let max_t2_eq1 = djy.len().max(djx.len()); + let mut eq1: Vec> = vec![vec![]; max_t2_eq1]; + for k in 0..djy.len() { + eq1[k] = poly_scale(&dix, djy[k]); + } + for k in 0..djx.len() { + let term = poly_scale(&diy, djx[k]); + eq1[k] = if eq1[k].is_empty() { + poly_scale(&term, -1.0) + } else { + poly_sub(&eq1[k], &term) + }; + } + + let max_t2_eq2 = gjy.len().max(gjx.len()); + let mut eq2: Vec> = vec![vec![]; max_t2_eq2]; + for k in 0..gjy.len() { + eq2[k] = poly_add(&eq2[k], &poly_scale(&dix, gjy[k])); + } + for k in 0..gjx.len() { + eq2[k] = poly_sub(&eq2[k], &poly_scale(&diy, gjx[k])); + } + let pure_t1 = poly_sub(&poly_mul(&diy, &gix), &poly_mul(&dix, &giy)); + eq2[0] = poly_add(&eq2[0], &pure_t1); + + (eq1, eq2) +} + +/// Sylvester resultant of two bivariate polynomials w.r.t. t2 (result: poly in t1). +fn sylvester_resultant(f: &[Vec], g: &[Vec]) -> Vec { + let m = f.len() - 1; + let n = g.len() - 1; + let size = m + n; + let mut matrix: Vec>> = vec![vec![vec![]; size]; size]; + for i in 0..n { + for k in 0..=m { + let col = i + k; + if col < size { + matrix[i][col] = f[m - k].clone(); + } + } + } + for i in 0..m { + for k in 0..=n { + let col = i + k; + if col < size { + matrix[n + i][col] = g[n - k].clone(); + } + } + } + poly_matrix_determinant(&mut matrix, size) +} + +/// Fraction-free (Bareiss) elimination for a matrix of polynomials. +fn poly_matrix_determinant(matrix: &mut Vec>>, n: usize) -> Vec { + let mut prev_pivot = vec![1.0]; + for col in 0..n { + let mut pivot_row = None; + for row in col..n { + if !is_zero_poly(&matrix[row][col]) { + pivot_row = Some(row); + break; + } + } + let pivot_row = match pivot_row { + Some(r) => r, + None => return vec![0.0], + }; + if pivot_row != col { + matrix.swap(pivot_row, col); + } + let pivot = matrix[col][col].clone(); + for row in (col + 1)..n { + for j in (col + 1)..n { + let term1 = poly_mul(&pivot, &matrix[row][j]); + let term2 = poly_mul(&matrix[row][col], &matrix[col][j]); + matrix[row][j] = poly_div_exact(&poly_sub(&term1, &term2), &prev_pivot); + } + matrix[row][col] = vec![0.0]; + } + prev_pivot = pivot; + } + matrix[n - 1][n - 1].clone() +} + +fn back_substitute_t2(eq1: &[Vec], t1: f64) -> Vec { + let coeffs: Vec = eq1.iter().map(|c| eval_poly(c, t1)).collect(); + let trimmed = trim_poly(&coeffs); + if is_zero_poly(trimmed) { + return vec![]; + } + PolyDyn::new(trimmed.iter().copied()).roots_between(0.0, 1.0, TOL) +} + +// ─── Crossings ─── + +/// A point where two candidates' support functions meet: a common supporting +/// direction with equal support value (a bitangency). +#[derive(Clone, Debug)] +struct Crossing { + cand_a: usize, + cand_b: usize, + theta: f64, + point_a: Point, + point_b: Point, + t_a: f64, + t_b: f64, +} + +/// All bitangent (t1, t2) solutions between two cubics. +fn solve_arc_arc(cb_i: &CubicBez, cb_j: &CubicBez) -> Vec<(f64, f64)> { + let (eq1, eq2) = build_bitangent_system(cb_i, cb_j); + let resultant = sylvester_resultant(&eq1, &eq2); + let trimmed = trim_poly(&resultant); + if is_zero_poly(trimmed) { + return vec![]; + } + let t1_roots = PolyDyn::new(trimmed.iter().copied()).roots_between(0.0, 1.0, TOL); + let mut out = Vec::new(); + for &t1 in &t1_roots { + for t2 in back_substitute_t2(&eq1, t1) { + if !(-TOL..=1.0 + TOL).contains(&t2) { + continue; + } + let t1c = t1.clamp(0.0, 1.0); + let t2c = t2.clamp(0.0, 1.0); + // Residual check with relative tolerance (eq2 scales as O(L²)). + let residual = eval_bivariate(&eq2, t1c, t2c); + let d = cb_i.deriv().eval(t1c).to_vec2(); + let disp = cb_j.eval(t2c) - cb_i.eval(t1c); + let scale = d.hypot() * disp.hypot(); + if residual.abs() < 1e-6 * scale.max(1.0) { + out.push((t1c, t2c)); + } + } + } + out +} + +/// All parameters t where the tangent line of `cb` at t passes through P. +fn solve_point_tangency(cb: &CubicBez, p: Point) -> Vec { + let g = curve_power_basis(cb); + let d = deriv_power_basis(cb); + let dpx = [g[0].x - p.x, g[1].x, g[2].x, g[3].x]; + let dpy = [g[0].y - p.y, g[1].y, g[2].y, g[3].y]; + let dx: Vec = d.iter().map(|v| v.x).collect(); + let dy: Vec = d.iter().map(|v| v.y).collect(); + let eq = poly_sub(&poly_mul(&dx, &dpy), &poly_mul(&dy, &dpx)); + let trimmed = trim_poly(&eq); + if is_zero_poly(trimmed) { + return vec![]; + } + PolyDyn::new(trimmed.iter().copied()).roots_between(0.0, 1.0, TOL) +} + +fn find_crossings(cands: &[Candidate], pieces: &[MonotoneArc], scale: f64) -> Vec { + let geom_eps = 1e-9 * scale; + let mut crossings = Vec::new(); + + for a in 0..cands.len() { + for b in (a + 1)..cands.len() { + let (ca, cb) = (&cands[a], &cands[b]); + if ca.piece == cb.piece { + continue; + } + if !intervals_overlap(ca, cb) { + continue; + } + let (pa, pb) = (&pieces[ca.piece], &pieces[cb.piece]); + + match (ca.kind, cb.kind) { + (PieceKind::Arc, PieceKind::Arc) => { + // Dedup: the same geometric pair is solved once per (a,b) + // candidate pair, but the algebraic solve doesn't depend on + // sides. Cache would help; correctness first. + for (t1, t2) in solve_arc_arc(&pa.bezier, &pb.bezier) { + let q1 = pa.bezier.eval(t1); + let q2 = pb.bezier.eval(t2); + if (q2 - q1).hypot() < geom_eps { + continue; + } + let d1 = tangent_dir_at(&pa.bezier, t1); + let d2 = tangent_dir_at(&pb.bezier, t2); + let th_a = normal_angle(d1, ca.side); + let th_b = normal_angle(d2, cb.side); + if ang_diff(th_a, th_b) > ANG_MATCH { + continue; + } + let theta = normalize_angle(th_a); + if !ca.contains(theta, ANG_EPS) || !cb.contains(theta, ANG_EPS) { + continue; + } + crossings.push(Crossing { + cand_a: a, + cand_b: b, + theta, + point_a: q1, + point_b: q2, + t_a: t1, + t_b: t2, + }); + } + } + (PieceKind::Arc, PieceKind::Corner) | (PieceKind::Arc, PieceKind::Line) + | (PieceKind::Corner, PieceKind::Arc) | (PieceKind::Line, PieceKind::Arc) => { + // Point-vs-arc tangency. For lines, both endpoints act as + // points; but line-endpoint junction corners already cover + // hull-relevant transitions, and a line interior can never + // be tangent from outside. Only corners need solving here. + let (arc_idx, pt_idx, arc_cand, pt_cand) = if ca.kind == PieceKind::Arc { + (ca.piece, cb.piece, a, b) + } else { + (cb.piece, ca.piece, b, a) + }; + if pieces[pt_idx].kind == PieceKind::Line { + continue; + } + let arc_piece = &pieces[arc_idx]; + let p = pieces[pt_idx].bezier.p0; + for t in solve_point_tangency(&arc_piece.bezier, p) { + let q = arc_piece.bezier.eval(t); + if (q - p).hypot() < geom_eps { + continue; + } + // Spurious roots where the derivative vanishes. + let draw = arc_piece.bezier.deriv().eval(t).to_vec2(); + let dir = tangent_dir_at(&arc_piece.bezier, t); + let chord = q - p; + if draw.hypot() < 1e-9 * chord_scale(&arc_piece.bezier) + && dir.cross(chord).abs() > 1e-6 * dir.hypot() * chord.hypot() + { + continue; + } + let arc_side = cands[arc_cand].side; + let theta = normalize_angle(normal_angle(dir, arc_side)); + if !cands[arc_cand].contains(theta, ANG_EPS) { + continue; + } + let (ta, tb, qa, qb) = if arc_cand == a { + (t, 0.0, q, p) + } else { + (0.0, t, p, q) + }; + let _ = pt_cand; + crossings.push(Crossing { + cand_a: a, + cand_b: b, + theta, + point_a: qa, + point_b: qb, + t_a: ta, + t_b: tb, + }); + } + } + (PieceKind::Corner, PieceKind::Corner) => { + let pi = pa.bezier.p0; + let pj = pb.bezier.p0; + let v = pj - pi; + if v.hypot() < geom_eps { + continue; + } + for theta in [ + normalize_angle(v.atan2() + PI / 2.0), + normalize_angle(v.atan2() - PI / 2.0), + ] { + crossings.push(Crossing { + cand_a: a, + cand_b: b, + theta, + point_a: pi, + point_b: pj, + t_a: 0.0, + t_b: 0.0, + }); + } + } + _ => {} // line-line, line-corner: transitions occur at + // breakpoints already contributed by the line angles. + } + } + } + crossings +} + +fn ang_diff(a: f64, b: f64) -> f64 { + let d = normalize_angle(a - b); + d.min(TAU - d) +} + +// ─── Envelope Sweep & Assembly ─── + +struct Envelope { + /// Breakpoint angles, sorted, covering the circle. + breaks: Vec, + /// Winner candidate index for each interval (breaks[i], breaks[i+1]). + winners: Vec, +} + +fn sweep_envelope(cands: &[Candidate], pieces: &[MonotoneArc], crossings: &[Crossing], scale: f64) -> Envelope { + let mut breaks: Vec = Vec::new(); + for c in cands { + if c.kind == PieceKind::Corner { + continue; + } + breaks.push(normalize_angle(c.ang_start)); + breaks.push(normalize_angle(c.ang_start + c.ang_len)); + } + for x in crossings { + breaks.push(x.theta); + } + breaks.sort_by(|a, b| a.partial_cmp(b).unwrap()); + breaks.dedup_by(|a, b| (*a - *b).abs() < ANG_EPS); + if breaks.is_empty() { + breaks.push(0.0); + } + // Circular dedup of first/last. + if breaks.len() > 1 && (TAU - breaks[breaks.len() - 1] + breaks[0]).abs() < ANG_EPS { + breaks.pop(); + } + + let eps_h = 1e-9 * scale; + let m = breaks.len(); + let mut winners = Vec::with_capacity(m); + for i in 0..m { + let a = breaks[i]; + let b = if i + 1 < m { breaks[i + 1] } else { breaks[0] + TAU }; + let mid = normalize_angle(a + ang_forward(a, normalize_angle(b)) / 2.0); + let mut best: Option<(f64, usize)> = None; + for (ci, c) in cands.iter().enumerate() { + if c.kind == PieceKind::Line { + continue; + } + if !c.contains(mid, ANG_EPS) { + continue; + } + let (h, _, _) = support_at(c, &pieces[c.piece], mid); + match best { + Some((bh, _)) if h <= bh + eps_h => { + // Tie: prefer arcs over corners (a corner coincident with + // an arc endpoint should yield to the arc). + if h >= bh - eps_h + && c.kind == PieceKind::Arc + && cands[best.unwrap().1].kind == PieceKind::Corner + { + best = Some((h, ci)); + } + } + _ => best = Some((h, ci)), + } + } + let w = best.expect("no active candidate — coverage gap").1; + winners.push(w); + } + Envelope { breaks, winners } +} + +/// One entry of a transition tie-set: a contact point on the common support line. +struct TieEntry { + cand: usize, + proj: f64, + point: Point, + t: f64, +} + +fn assemble(cands: &[Candidate], pieces: &[MonotoneArc], crossings: &[Crossing], env: &Envelope, scale: f64) -> Vec { + let m = env.breaks.len(); + let geom_eps = 1e-7 * scale; + let eps_h = 1e-7 * scale; + + // Merge consecutive intervals with the same winner into runs. + // runs: (winner, theta_from, theta_to) with theta_to lifted ≥ theta_from. + let mut run_bounds: Vec = Vec::new(); // indices into breaks where winner changes + for i in 0..m { + let prev = env.winners[(i + m - 1) % m]; + if env.winners[i] != prev { + run_bounds.push(i); + } + } + + if run_bounds.is_empty() { + // Single winner covers everything: hull is that single closed piece + // (or one arc traversed fully) — emit all its pieces. + let w = env.winners[0]; + let c = &cands[w]; + if c.kind == PieceKind::Arc { + // A closed convex curve decomposed into one arc — unusual but emit fully. + return vec![HullSegment::Arc { + arc_index: c.piece, + t_start: c.t_at_start, + t_end: c.t_at_end, + }]; + } + return vec![]; + } + + // For crossing lookup at transitions. + let find_crossing = |x: usize, y: usize, theta: f64| -> Option<&Crossing> { + crossings + .iter() + .filter(|c| { + ((c.cand_a == x && c.cand_b == y) || (c.cand_a == y && c.cand_b == x)) + && ang_diff(c.theta, theta) < 1e-7 + }) + .min_by(|p, q| { + ang_diff(p.theta, theta) + .partial_cmp(&ang_diff(q.theta, theta)) + .unwrap() + }) + }; + + let mut segments: Vec = Vec::new(); + let nb = run_bounds.len(); + + for ri in 0..nb { + // Run: winner w from break run_bounds[ri] to run_bounds[(ri+1) % nb]. + let start_bi = run_bounds[ri]; + let end_bi = run_bounds[(ri + 1) % nb]; + let w = env.winners[start_bi]; + let c = &cands[w]; + let theta_in = env.breaks[start_bi]; + let theta_out = env.breaks[end_bi]; + + // Entry/exit contact parameters for arc winners. + if c.kind == PieceKind::Arc { + let prev_w = env.winners[(start_bi + m - 1) % m]; + let next_w = env.winners[end_bi]; + let t_in = find_crossing(prev_w, w, theta_in) + .map(|x| if x.cand_a == w { x.t_a } else { x.t_b }) + .unwrap_or_else(|| arc_param_at_normal(c, &pieces[c.piece], theta_in)); + let t_out = find_crossing(w, next_w, theta_out) + .map(|x| if x.cand_a == w { x.t_a } else { x.t_b }) + .unwrap_or_else(|| arc_param_at_normal(c, &pieces[c.piece], theta_out)); + if (t_out - t_in).abs() > 1e-12 { + segments.push(HullSegment::Arc { + arc_index: c.piece, + t_start: t_in, + t_end: t_out, + }); + } + } + + // Transition at theta_out between w and the next run's winner. + let next_w = env.winners[end_bi]; + let theta_c = theta_out; + let n = Vec2::new(theta_c.cos(), theta_c.sin()); + let tau = Vec2::new(-theta_c.sin(), theta_c.cos()); + + // Gather tie set: every candidate achieving max support at theta_c. + // Line pieces at this angle are tracked separately so consecutive + // nodes connected by an input line emit that line, not a bitangent. + let mut max_h = f64::NEG_INFINITY; + let mut entries: Vec = Vec::new(); + let mut line_ties: Vec = Vec::new(); // candidate indices + for (ci, cc) in cands.iter().enumerate() { + if cc.kind == PieceKind::Line { + let d = ang_diff(normalize_angle(cc.ang_start), theta_c); + if d > 1e-7 { + continue; + } + // Both endpoints of the line are contacts. + let (q0, q3) = (pieces[cc.piece].bezier.p0, pieces[cc.piece].bezier.p3); + let h = q0.to_vec2().dot(n); + max_h = max_h.max(h); + line_ties.push(ci); + entries.push(TieEntry { + cand: ci, + proj: q0.to_vec2().dot(tau), + point: q0, + t: 0.0, + }); + entries.push(TieEntry { + cand: ci, + proj: q3.to_vec2().dot(tau), + point: q3, + t: 1.0, + }); + continue; + } + if !cc.contains(theta_c, 1e-7) { + continue; + } + let (h, p, t) = support_at(cc, &pieces[cc.piece], theta_c); + max_h = max_h.max(h); + entries.push(TieEntry { + cand: ci, + proj: p.to_vec2().dot(tau), + point: p, + t, + }); + } + // The exiting/entering winners define this transition: use their exact + // algebraic crossing contacts when available, and always keep them + // (by continuity of the envelope they are at the max). + if let Some(x) = find_crossing(w, next_w, theta_c) { + for e in entries.iter_mut() { + let (p, t) = if x.cand_a == e.cand { + (x.point_a, x.t_a) + } else if x.cand_b == e.cand { + (x.point_b, x.t_b) + } else { + continue; + }; + e.point = p; + e.t = t; + e.proj = p.to_vec2().dot(tau); + } + } + // Keep only entries at the max support level (winners always kept). + let mut tie: Vec = entries + .into_iter() + .filter(|e| { + e.cand == w || e.cand == next_w || e.point.to_vec2().dot(n) >= max_h - eps_h + }) + .collect(); + tie.sort_by(|a, b| a.proj.partial_cmp(&b.proj).unwrap()); + + // Deduplicate coincident contact points, preferring the exiting winner + // first and the entering winner last. + let mut nodes: Vec = Vec::new(); + for e in tie { + if let Some(last) = nodes.last() { + if (e.point - last.point).hypot() < geom_eps { + // Same geometric node: keep the more relevant candidate. + let keep_new = e.cand == w || e.cand == next_w; + let keep_old = last.cand == w || last.cand == next_w; + if keep_new && !keep_old { + nodes.pop(); + nodes.push(e); + } + continue; + } + } + nodes.push(e); + } + + // Thread through the nodes in traversal order. A gap between + // consecutive nodes is an input line piece if one spans it, otherwise + // a bitangent line. + for k in 0..nodes.len().saturating_sub(1) { + let (e1, e2) = (&nodes[k], &nodes[k + 1]); + let spanning_line = line_ties.iter().copied().find(|&ci| { + let lb = &pieces[cands[ci].piece].bezier; + ((lb.p0 - e1.point).hypot() < geom_eps && (lb.p3 - e2.point).hypot() < geom_eps) + || ((lb.p3 - e1.point).hypot() < geom_eps + && (lb.p0 - e2.point).hypot() < geom_eps) + }); + if let Some(ci) = spanning_line { + let lb = &pieces[cands[ci].piece].bezier; + let forward = (lb.p0 - e1.point).hypot() < geom_eps; + segments.push(HullSegment::Arc { + arc_index: cands[ci].piece, + t_start: if forward { 0.0 } else { 1.0 }, + t_end: if forward { 1.0 } else { 0.0 }, + }); + } else { + segments.push(HullSegment::Line { start: e1.point, end: e2.point }); + } + } + } + + segments +} + +// ─── Public API ─── + +/// Convex hull of multiple closed splines. Hull `arc_index` values reference +/// the pieces returned by `split_loops_into_arcs` for the same loops. +pub fn convex_hull_loops(loops: &[&[CubicBez]]) -> Vec { + let pieces = split_loops_into_arcs(loops); + if pieces.is_empty() { + return vec![]; + } + + let mut scale: f64 = 0.0; + for p in &pieces { + for cp in [p.bezier.p0, p.bezier.p1, p.bezier.p2, p.bezier.p3] { + scale = scale.max(cp.x.abs()).max(cp.y.abs()); + } + } + let scale = scale.max(1e-9); + + let cands = build_candidates(&pieces); + let crossings = find_crossings(&cands, &pieces, scale); + let env = sweep_envelope(&cands, &pieces, &crossings, scale); + assemble(&cands, &pieces, &crossings, &env, scale) +} diff --git a/node-graph/nodes/convex_hull/src/lib.rs b/node-graph/nodes/convex_hull/src/lib.rs index 4d4f6b4967..8b5e50fddd 100644 --- a/node-graph/nodes/convex_hull/src/lib.rs +++ b/node-graph/nodes/convex_hull/src/lib.rs @@ -1,94 +1,14 @@ use core_types::Ctx; -use core_types::table::{Table, TableRow, TableRowRef}; -use glam::{DAffine2, DVec2}; -use graphic_types::Vector; -use graphic_types::vector_types::subpath::{ManipulatorGroup, PathSegPoints, Subpath, pathseg_points}; +use core_types::list::{ATTR_EDITOR_MERGED_LAYERS, ATTR_TRANSFORM, Item, List}; +use glam::DAffine2; +use graphic_types::vector_types::subpath::Subpath; use graphic_types::vector_types::vector::PointId; use graphic_types::vector_types::vector::algorithms::merge_by_distance::MergeByDistanceExt; -pub use path_bool as path_bool_lib; -use path_bool::{FillRule, PathBooleanOperation}; -use std::ops::Mul; +use graphic_types::{Graphic, IntoGraphicList, Vector}; +use vector_types::kurbo::{Affine, CubicBez, Line as KurboLine, ParamCurve, PathSeg as KurboPathSeg}; -use ::convex_hull::{HullSegment, MonotoneArc, convex_hull as compute_convex_hull, split_into_arcs}; -use kurbo::{CubicBez, Line as KurboLine, ParamCurve, PathSeg as KurboPathSeg, Point as KurboPoint}; - -// ─── Graham's Scan Convex Hull ─── - -/// Compute the convex hull of a set of 2D points using Graham's scan. -/// Returns points in counter-clockwise order. -fn graham_scan_hull(points: &[DVec2]) -> Vec { - if points.len() <= 2 { - return points.to_vec(); - } - - // Find the lowest-y point (leftmost if tied) - let mut pivot_idx = 0; - for (i, p) in points.iter().enumerate() { - if p.y < points[pivot_idx].y || (p.y == points[pivot_idx].y && p.x < points[pivot_idx].x) { - pivot_idx = i; - } - } - let pivot = points[pivot_idx]; - - // Sort remaining points by polar angle from pivot - let mut indexed: Vec<(usize, DVec2)> = points.iter().copied().enumerate().filter(|&(i, _)| i != pivot_idx).collect(); - indexed.sort_by(|&(_, a), &(_, b)| { - let da = a - pivot; - let db = b - pivot; - let angle_a = da.y.atan2(da.x); - let angle_b = db.y.atan2(db.x); - angle_a.partial_cmp(&angle_b).unwrap().then_with(|| { - // If same angle, closer point first - da.length_squared().partial_cmp(&db.length_squared()).unwrap() - }) - }); - - // Build hull using cross-product left-turn test - let mut hull = vec![pivot]; - for (_, p) in indexed { - while hull.len() >= 2 { - let a = hull[hull.len() - 2]; - let b = hull[hull.len() - 1]; - let cross = (b - a).perp_dot(p - b); - if cross <= 0.0 { - hull.pop(); - } else { - break; - } - } - hull.push(p); - } - - hull -} - -// ─── Kurbo PathSeg → CubicBez Conversion ─── - -/// Convert any `kurbo::PathSeg` to a `CubicBez`. -fn pathseg_to_cubicbez(seg: KurboPathSeg) -> CubicBez { - match seg { - KurboPathSeg::Cubic(cb) => cb, - KurboPathSeg::Quad(qb) => { - // Degree elevation: quadratic → cubic - let p0 = qb.p0; - let p3 = qb.p2; - let q1 = qb.p1; - let p1 = KurboPoint::new(p0.x + 2.0 / 3.0 * (q1.x - p0.x), p0.y + 2.0 / 3.0 * (q1.y - p0.y)); - let p2 = KurboPoint::new(p3.x + 2.0 / 3.0 * (q1.x - p3.x), p3.y + 2.0 / 3.0 * (q1.y - p3.y)); - CubicBez::new(p0, p1, p2, p3) - } - KurboPathSeg::Line(l) => { - // Place control points at 1/3 and 2/3 along the line - let p0 = l.p0; - let p3 = l.p1; - let p1 = KurboPoint::new(p0.x + (p3.x - p0.x) / 3.0, p0.y + (p3.y - p0.y) / 3.0); - let p2 = KurboPoint::new(p0.x + 2.0 * (p3.x - p0.x) / 3.0, p0.y + 2.0 * (p3.y - p0.y) / 3.0); - CubicBez::new(p0, p1, p2, p3) - } - } -} - -// ─── Subpath → Vec Conversion ─── +mod hull; +use hull::{HullSegment, MonotoneArc, convex_hull_loops, split_loops_into_arcs}; /// Check if a CubicBez is degenerate (all control points at essentially the same location). fn is_degenerate_cubic(cb: &CubicBez) -> bool { @@ -99,49 +19,13 @@ fn is_degenerate_cubic(cb: &CubicBez) -> bool { (d03.x * d03.x + d03.y * d03.y) < EPS_SQ && (d01.x * d01.x + d01.y * d01.y) < EPS_SQ && (d02.x * d02.x + d02.y * d02.y) < EPS_SQ } -/// Convert a `Subpath` into a `Vec` for the convex hull library. -/// Filters out degenerate zero-length segments. -fn subpath_to_cubicbez_vec(subpath: &Subpath) -> Vec { - subpath.iter().map(pathseg_to_cubicbez).filter(|cb| !is_degenerate_cubic(cb)).collect() +/// Convert a `Subpath` into a closed loop of `CubicBez` segments in world space. +/// Open subpaths are closed with a line segment; degenerate segments are dropped. +fn subpath_to_loop(subpath: &Subpath, transform: DAffine2) -> Vec { + let affine = Affine::new(transform.to_cols_array()); + subpath.iter_closed().map(|seg| (affine * seg).to_cubic()).filter(|cb| !is_degenerate_cubic(cb)).collect() } -// ─── Winding Direction ─── - -/// Compute the signed area of a closed cubic bezier path by sampling. -/// Positive = CCW in standard math coords, Negative = CW. -fn signed_area_of_cubic_path(segments: &[CubicBez]) -> f64 { - let mut area = 0.0; - let n = 16; - for seg in segments { - for i in 0..n { - let t0 = i as f64 / n as f64; - let t1 = (i + 1) as f64 / n as f64; - let p0 = seg.eval(t0); - let p1 = seg.eval(t1); - area += p0.x * p1.y - p1.x * p0.y; - } - } - area / 2.0 -} - -/// Reverse a cubic bezier path (reverse segment order + swap endpoints within each segment). -fn reverse_cubic_path(segments: &[CubicBez]) -> Vec { - segments.iter().rev().map(|cb| CubicBez::new(cb.p3, cb.p2, cb.p1, cb.p0)).collect() -} - -// ─── Select Outer Subpath ─── - -/// Select the outermost subpath from a Vector by choosing the one with the largest absolute area. -fn select_outer_subpath(vector: &Vector) -> Option> { - vector.stroke_bezier_paths().max_by(|a, b| { - let area_a = a.area_centroid_and_area(None, None).map(|(_, area)| area.abs()).unwrap_or(0.0); - let area_b = b.area_centroid_and_area(None, None).map(|(_, area)| area.abs()).unwrap_or(0.0); - area_a.partial_cmp(&area_b).unwrap_or(std::cmp::Ordering::Equal) - }) -} - -// ─── Hull Segments → Subpath ─── - /// Convert hull segments back into a `Subpath`. fn hull_segments_to_subpath(segments: &[HullSegment], arcs: &[MonotoneArc]) -> Option> { let mut kurbo_segs: Vec = segments @@ -162,296 +46,137 @@ fn hull_segments_to_subpath(segments: &[HullSegment], arcs: &[MonotoneArc]) -> O // Subpath::from_beziers requires at least 2 segments for a closed path. // If we have only 1, split it at the midpoint. if kurbo_segs.len() == 1 { - let seg = kurbo_segs[0]; - match seg { - KurboPathSeg::Cubic(cb) => { - let first_half = cb.subsegment(0.0..0.5); - let second_half = cb.subsegment(0.5..1.0); - kurbo_segs = vec![KurboPathSeg::Cubic(first_half), KurboPathSeg::Cubic(second_half)]; - } - KurboPathSeg::Line(l) => { - let mid = KurboPoint::new((l.p0.x + l.p1.x) / 2.0, (l.p0.y + l.p1.y) / 2.0); - kurbo_segs = vec![KurboPathSeg::Line(KurboLine::new(l.p0, mid)), KurboPathSeg::Line(KurboLine::new(mid, l.p1))]; - } - KurboPathSeg::Quad(qb) => { - let cb = pathseg_to_cubicbez(KurboPathSeg::Quad(qb)); - let first_half = cb.subsegment(0.0..0.5); - let second_half = cb.subsegment(0.5..1.0); - kurbo_segs = vec![KurboPathSeg::Cubic(first_half), KurboPathSeg::Cubic(second_half)]; - } - } + let cb = kurbo_segs[0].to_cubic(); + kurbo_segs = vec![KurboPathSeg::Cubic(cb.subsegment(0.0..0.5)), KurboPathSeg::Cubic(cb.subsegment(0.5..1.0))]; } Some(Subpath::from_beziers(&kurbo_segs, true)) } -// ─── Main Node ─── +/// Exact convex hull of a set of closed loops, as a subpath. +fn compute_hull_subpath(loops: &[&[CubicBez]]) -> Option> { + if loops.is_empty() { + return None; + } + let arcs = split_loops_into_arcs(loops); + let hull_segments = convex_hull_loops(loops); + hull_segments_to_subpath(&hull_segments, &arcs) +} #[node_macro::node(category("Vector: Modifier"), path(core_types::vector))] -async fn convex_hull(_: impl Ctx, content: Table) -> Table { - // Handle empty input - if content.is_empty() { - return Table::default(); +async fn convex_hull(_: impl Ctx, #[implementations(List, List)] content: I) -> List { + let content = content.into_graphic_list(); + let flattened: List = content.clone().into_flattened_list(); + if flattened.is_empty() { + return List::default(); } - // Step 1: Collect one representative point per subpath (in world space) - let mut hull_points: Vec = Vec::new(); - for row in content.iter() { - let transform = *row.transform; - for subpath in row.element.stroke_bezier_paths() { - if let Some(first) = subpath.manipulator_groups().first() { - hull_points.push(transform.transform_point2(first.anchor)); + // Collect every subpath of every element as a closed loop in world space. + // The hull library handles multiple loops, occlusion, nesting, arbitrary + // winding, and self-intersections directly, so no union or winding + // normalization is needed. + let mut loops: Vec> = Vec::new(); + for index in 0..flattened.len() { + let Some(element) = flattened.element(index) else { continue }; + let transform: DAffine2 = flattened.attribute_cloned_or_default(ATTR_TRANSFORM, index); + for subpath in element.stroke_bezier_paths() { + let segs = subpath_to_loop(&subpath, transform); + if !segs.is_empty() { + loops.push(segs); } } } - // Step 2: Union all input shapes - let mut result_vector_table = union(content.iter()); - - // Step 3: Flatten union result to world space (apply transform, set to IDENTITY) - let style; - { - let Some(result_row) = result_vector_table.iter_mut().next() else { - return Table::default(); - }; - let transform = *result_row.transform; - *result_row.transform = DAffine2::IDENTITY; - Vector::transform(result_row.element, transform); - result_row.element.style.set_stroke_transform(DAffine2::IDENTITY); - - // Step 4: Save style - style = result_row.element.style.clone(); - } - - // Step 5: If the union has multiple disjoint subpaths AND we have ≥3 hull points, - // build a polyline convex hull and boolean-union it with the result to connect everything. - let subpath_count = result_vector_table.iter().next().map(|r| r.element.stroke_bezier_paths().count()).unwrap_or(0); - log::debug!("subpath_count: {}", subpath_count); - - if subpath_count > 1 && hull_points.len() >= 3 { - let poly_points = graham_scan_hull(&hull_points); - if poly_points.len() >= 3 { - // Build a polyline subpath from the hull points - let poly_subpath = Subpath::::from_anchors(poly_points.into_iter(), true); - let poly_vector = Vector::from_subpath(poly_subpath); - - // Boolean union the current result with the polyline - let current_vector = &result_vector_table.iter().next().unwrap().element; - let upper_path = to_path(current_vector, DAffine2::IDENTITY); - let lower_path = to_path(&poly_vector, DAffine2::IDENTITY); - - #[allow(unused_unsafe)] - let union_result_paths = unsafe { boolean_union(upper_path, lower_path) }; - let union_result = from_path(&union_result_paths); - - // Replace the result vector's geometry - let result_row = result_vector_table.iter_mut().next().unwrap(); - result_row.element.colinear_manipulators = union_result.colinear_manipulators; - result_row.element.point_domain = union_result.point_domain; - result_row.element.segment_domain = union_result.segment_domain; - result_row.element.region_domain = union_result.region_domain; - } - } - - // Step 6: Select the outer boundary subpath (largest by area) - let outer_subpath = { - let result_row = result_vector_table.iter().next().unwrap(); - select_outer_subpath(result_row.element) + let loop_refs: Vec<&[CubicBez]> = loops.iter().map(|l| l.as_slice()).collect(); + let Some(hull_subpath) = compute_hull_subpath(&loop_refs) else { + // Degenerate input (e.g. only zero-length segments): pass through. + return flattened; }; - let Some(outer_subpath) = outer_subpath else { - return result_vector_table; - }; + // Carry the first input item's paint attributes onto the hull. The hull + // geometry is already in world space, so the cloned transform attribute + // must be reset to identity or it would be applied a second time. + let paint_attributes = flattened.clone_item_attributes(0); + let hull_vector = Vector::from_subpath(hull_subpath); + let mut result = List::new_from_item(Item::from_parts(hull_vector, paint_attributes)); + result.set_attribute(ATTR_TRANSFORM, 0, DAffine2::IDENTITY); - // Step 7: Convert to Vec - let cubic_segments = subpath_to_cubicbez_vec(&outer_subpath); - if cubic_segments.is_empty() { - return result_vector_table; - } + // Snapshot the input layers as `editor:merged_layers` so the renderer can + // recurse into them and keep the original layers' overlays and click + // targets in place (same pattern as Boolean Operation and Flatten Path). + // No transform pre-compensation is needed since item 0's transform is + // identity. + result.set_attribute(ATTR_EDITOR_MERGED_LAYERS, 0, content); - // The hull library expects CCW winding. Graphite paths are typically CW in screen coords - // (Y-down), so we reverse if the signed area is negative (CW in math coords). - let cubic_segments = if signed_area_of_cubic_path(&cubic_segments) < 0.0 { - reverse_cubic_path(&cubic_segments) - } else { - cubic_segments - }; - - log::debug!("path: {:?}", cubic_segments); - - // Step 8: Run the curved convex hull algorithm - let arcs = split_into_arcs(&cubic_segments); - log::debug!("arcs: {:?}", arcs); - let hull_segments = compute_convex_hull(&cubic_segments); - log::debug!("segments: {:?}", hull_segments); - - if hull_segments.is_empty() { - // Fallback: return the union result as-is - return result_vector_table; - } - - // Step 9: Reconstruct hull as Subpath - let Some(hull_subpath) = hull_segments_to_subpath(&hull_segments, &arcs) else { - return result_vector_table; - }; - log::debug!("hull_subpath: {:?}", hull_subpath); - - // Step 10: Create Vector from hull subpath, apply saved style - let mut hull_vector = Vector::from_subpath(hull_subpath); - hull_vector.style = style; - - // Step 11: Build result table - let mut result: Table = Table::new_from_element(hull_vector); - if let Some(row) = result.iter_mut().next() { - // Step 11: Clean up with merge_by_distance_spatial - row.element.merge_by_distance_spatial(*row.transform, 0.0001); + if let Some(element) = result.element_mut(0) { + element.merge_by_distance_spatial(DAffine2::IDENTITY, 0.0001); } result } -// ─── Boolean Operations (shared helpers) ─── +#[cfg(test)] +mod tests { + use super::*; + use vector_types::kurbo::Point as KurboPoint; -fn union<'a>(vector: impl DoubleEndedIterator>) -> Table { - // Reverse the vector table rows so that the result style is the style of the first vector row - let mut vector_reversed = vector.rev(); - - let mut result_vector_table = Table::new_from_row(vector_reversed.next().map(|x| x.into_cloned()).unwrap_or_default()); - let mut first_row = result_vector_table.iter_mut().next().expect("Expected the one row we just pushed"); - - // Loop over all vector table rows and union it with the result - let default = TableRow::default(); - let mut second_vector = Some(vector_reversed.next().unwrap_or(default.as_ref())); - while let Some(lower_vector) = second_vector { - let transform_of_lower_into_space_of_upper = first_row.transform.inverse() * *lower_vector.transform; - - let result = &mut first_row.element; - - let upper_path_string = to_path(result, DAffine2::IDENTITY); - let lower_path_string = to_path(lower_vector.element, transform_of_lower_into_space_of_upper); - - #[allow(unused_unsafe)] - let boolean_operation_string = unsafe { boolean_union(upper_path_string, lower_path_string) }; - let boolean_operation_result = from_path(&boolean_operation_string); - - result.colinear_manipulators = boolean_operation_result.colinear_manipulators; - result.point_domain = boolean_operation_result.point_domain; - result.segment_domain = boolean_operation_result.segment_domain; - result.region_domain = boolean_operation_result.region_domain; - - second_vector = vector_reversed.next(); + fn circle_at(cx: f64, cy: f64, r: f64) -> Vec { + let k = 0.5522847498 * r; + vec![ + CubicBez::new((cx + r, cy), (cx + r, cy + k), (cx + k, cy + r), (cx, cy + r)), + CubicBez::new((cx, cy + r), (cx - k, cy + r), (cx - r, cy + k), (cx - r, cy)), + CubicBez::new((cx - r, cy), (cx - r, cy - k), (cx - k, cy - r), (cx, cy - r)), + CubicBez::new((cx, cy - r), (cx + k, cy - r), (cx + r, cy - k), (cx + r, cy)), + ] } - result_vector_table -} - -fn to_path(vector: &Vector, transform: DAffine2) -> Vec { - let mut path = Vec::new(); - for subpath in vector.stroke_bezier_paths() { - to_path_segments(&mut path, &subpath, transform); - } - path -} - -fn to_path_segments(path: &mut Vec, subpath: &Subpath, transform: DAffine2) { - use path_bool::PathSegment; - let mut global_start = None; - let mut global_end = DVec2::ZERO; - - for bezier in subpath.iter() { - const EPS: f64 = 1e-8; - let transform_point = |pos: DVec2| transform.transform_point2(pos).mul(EPS.recip()).round().mul(EPS); - - let PathSegPoints { p0, p1, p2, p3 } = pathseg_points(bezier); - - let p0 = transform_point(p0); - let p1 = p1.map(transform_point); - let p2 = p2.map(transform_point); - let p3 = transform_point(p3); - - if global_start.is_none() { - global_start = Some(p0); - } - global_end = p3; - - let segment = match (p1, p2) { - (None, None) => PathSegment::Line(p0, p3), - (None, Some(p2)) | (Some(p2), None) => PathSegment::Quadratic(p0, p2, p3), - (Some(p1), Some(p2)) => PathSegment::Cubic(p0, p1, p2, p3), - }; - - path.push(segment); - } - if let Some(start) = global_start { - path.push(PathSegment::Line(global_end, start)); - } -} - -fn from_path(path_data: &[Path]) -> Vector { - const EPSILON: f64 = 1e-5; - - fn is_close(a: DVec2, b: DVec2) -> bool { - (a - b).length_squared() < EPSILON * EPSILON - } - - let mut all_subpaths = Vec::new(); - - for path in path_data.iter().filter(|path| !path.is_empty()) { - let cubics: Vec<[DVec2; 4]> = path.iter().map(|segment| segment.to_cubic()).collect(); - let mut manipulators_list = Vec::new(); - let mut current_start = None; - - for (index, cubic) in cubics.iter().enumerate() { - let [start, handle1, handle2, end] = *cubic; - - if current_start.is_none() || !is_close(start, current_start.unwrap()) { - // Start a new subpath - if !manipulators_list.is_empty() { - all_subpaths.push(Subpath::new(std::mem::take(&mut manipulators_list), true)); - } - // Use the correct in-handle (None) and out-handle for the start point - manipulators_list.push(ManipulatorGroup::new(start, None, Some(handle1))); - } else { - // Update the out-handle of the previous point - if let Some(last) = manipulators_list.last_mut() { - last.out_handle = Some(handle1); + fn assert_closed_and_contains(subpath: &Subpath, points: &[(f64, f64)]) { + assert!(subpath.closed()); + assert!(subpath.manipulator_groups().len() >= 2); + // Sampled containment check against the hull polygon. + let poly: Vec = subpath.iter_closed().flat_map(|seg| (0..32).map(move |k| seg.eval(k as f64 / 32.0))).collect(); + for &(px, py) in points { + // Point-in-polygon via ray casting. + let mut inside = false; + for i in 0..poly.len() { + let a = poly[i]; + let b = poly[(i + 1) % poly.len()]; + if (a.y > py) != (b.y > py) && px < a.x + (b.x - a.x) * (py - a.y) / (b.y - a.y) { + inside = !inside; } } - - // Add the end point with the correct in-handle and out-handle (None) - manipulators_list.push(ManipulatorGroup::new(end, Some(handle2), None)); - - current_start = Some(end); - - // Check if this is the last segment - if index == cubics.len() - 1 { - all_subpaths.push(Subpath::new(manipulators_list, true)); - manipulators_list = Vec::new(); // Reset manipulators for the next path - } + assert!(inside, "point ({px}, {py}) not inside hull"); } } - Vector::from_subpaths(all_subpaths, false) -} + #[test] + fn hull_of_two_disjoint_circles_spans_both() { + let a = circle_at(0.0, 0.0, 1.0); + let b = circle_at(5.0, 0.0, 1.0); + let hull = compute_hull_subpath(&[&a, &b]).unwrap(); + // Previously (union + Graham bridging) two disjoint shapes fell through + // the >= 3 anchor guard and the hull covered only one of them. + assert_closed_and_contains(&hull, &[(0.0, 0.0), (5.0, 0.0), (2.5, 0.5)]); + } -type Path = Vec; + #[test] + fn hull_of_nested_circles_is_outer() { + let outer = circle_at(0.0, 0.0, 2.0); + let inner = circle_at(0.3, 0.1, 1.0); + let hull = compute_hull_subpath(&[&outer, &inner]).unwrap(); + assert_closed_and_contains(&hull, &[(0.0, 0.0), (1.9, 0.0), (0.0, -1.9)]); + } -fn boolean_union(a: Path, b: Path) -> Vec { - path_bool(a, b, PathBooleanOperation::Union) -} - -fn path_bool(a: Path, b: Path, op: PathBooleanOperation) -> Vec { - match path_bool::path_boolean(&a, FillRule::NonZero, &b, FillRule::NonZero, op) { - Ok(results) => results, - Err(e) => { - let a_path = path_bool::path_to_path_data(&a, 0.001); - let b_path = path_bool::path_to_path_data(&b, 0.001); - log::error!("Boolean error {e:?} encountered while processing {a_path}\n {op:?}\n {b_path}"); - Vec::new() - } + #[test] + fn hull_handles_degenerate_line_parameterization() { + // Rectangle with zero-derivative endpoints, as produced by real paths. + let rect = vec![ + CubicBez::new((0.0, 0.0), (0.0, 0.0), (4.0, 0.0), (4.0, 0.0)), + CubicBez::new((4.0, 0.0), (4.0, 0.0), (4.0, 2.0), (4.0, 2.0)), + CubicBez::new((4.0, 2.0), (4.0, 2.0), (0.0, 2.0), (0.0, 2.0)), + CubicBez::new((0.0, 2.0), (0.0, 2.0), (0.0, 0.0), (0.0, 0.0)), + ]; + let hull = compute_hull_subpath(&[&rect]).unwrap(); + assert_closed_and_contains(&hull, &[(2.0, 1.0), (0.1, 0.1), (3.9, 1.9)]); } } - -pub fn boolean_intersect(a: Path, b: Path) -> Vec { - path_bool(a, b, PathBooleanOperation::Intersection) -}