1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
|
// Copyright 2022 the Kurbo Authors
// SPDX-License-Identifier: Apache-2.0 OR MIT
//! Test case that demonstrates quadratic/quadratic intersection using
//! new quartic solver.
use arrayvec::ArrayVec;
use kurbo::{common::solve_quartic, ParamCurve, Point, QuadBez, Shape};
use rand::{thread_rng, Rng};
fn rand_point() -> Point {
let mut rng = thread_rng();
Point::new(rng.gen_range(0.0..500.0), rng.gen_range(0.0..500.0))
}
fn rand_quad() -> QuadBez {
QuadBez::new(rand_point(), rand_point(), rand_point())
}
struct ImplicitQuad {
x2: f64,
xy: f64,
y2: f64,
x: f64,
y: f64,
c: f64,
}
impl ImplicitQuad {
fn from_quad(q: QuadBez) -> Self {
let Point { x: ax, y: ay } = q.p0;
let Point { x: bx, y: by } = q.p1;
let Point { x: cx, y: cy } = q.p2;
let (u0, u1, u2) = (by - cy, cx - bx, bx * cy - by * cx);
let (v0, v1, v2) = (cy - ay, ax - cx, cx * ay - cy * ax);
let (w0, w1, w2) = (ay - by, bx - ax, ax * by - ay * bx);
ImplicitQuad {
x2: 4. * u0 * w0 - v0 * v0,
xy: 4. * (u0 * w1 + u1 * w0) - 2. * v0 * v1,
y2: 4. * u1 * w1 - v1 * v1,
x: 4. * (u0 * w2 + u2 * w0) - 2. * v0 * v2,
y: 4. * (u1 * w2 + u2 * w1) - 2. * v1 * v2,
c: 4. * u2 * w2 - v2 * v2,
}
}
fn eval(&self, x: f64, y: f64) -> f64 {
self.x2 * x * x + self.xy * x * y + self.y2 * y * y + self.x * x + self.y * y + self.c
}
}
fn intersect_quads(q0: QuadBez, q1: QuadBez) -> ArrayVec<f64, 4> {
let iq = ImplicitQuad::from_quad(q0);
let c = q1.p0.to_vec2();
let b = (q1.p1 - q1.p0) * 2.0;
let a = c - q1.p1.to_vec2() * 2.0 + q1.p2.to_vec2();
let c0 = iq.eval(c.x, c.y);
let c1 = iq.x * b.x
+ iq.y * b.y
+ 2. * iq.x2 * (b.x * c.x)
+ 2. * iq.y2 * (b.y * c.y)
+ iq.xy * (b.x * c.y + b.y * c.x);
let c2 = iq.x * a.x
+ iq.y * a.y
+ iq.x2 * (2. * a.x * c.x + b.x * b.x)
+ iq.xy * (a.x * c.y + b.x * b.y + a.y * c.x)
+ iq.y2 * (2. * a.y * c.y + b.y * b.y);
let c3 = iq.x2 * 2. * a.x * b.x + iq.xy * (a.x * b.y + b.x * a.y) + iq.y2 * 2. * a.y * b.y;
let c4 = iq.x2 * a.x * a.x + iq.xy * a.x * a.y + iq.y2 * a.y * a.y;
let ts = solve_quartic(c0, c1, c2, c3, c4);
println!("ts: {:?}", ts);
ts
}
fn main() {
let q0 = rand_quad();
let q1 = rand_quad();
println!(
" <path d='{}' stroke='#000' fill='none' />",
q0.to_path(1e-9).to_svg()
);
println!(
" <path d='{}' stroke='#008' fill='none' />",
q1.to_path(1e-9).to_svg()
);
for t in intersect_quads(q0, q1) {
if (0.0..=1.0).contains(&t) {
let p = q1.eval(t);
println!(" <circle cx='{}' cy='{}' r='3' />", p.x, p.y);
}
}
}
|
