Bézier curves/Intersections
You are encouraged to solve this task according to the task description, using any language you may know.
You are given two planar quadratic Bézier curves, having control points and , respectively. They are parabolas intersecting at four points, as shown in the following diagram:
The task is to write a program that finds all four intersection points and prints their coordinates. You may use any algorithm you know of or can think of, including any of those that others have used.
See also
ATS
This program flattens one of the curves (that is, converts it to a piecewise linear approximation) and finds intersections between the line segments and the other curve. This requires solving many quadratic equations, but that can be done by the quadratic formula.
(I do the flattening in part by using a representation of Bézier curves that probably is not widely known. For quadratic splines, the representation amounts to a line plus a quadratic term. When the quadratic term gets small enough, I simply remove it. Here, though, is an interesting blog post that talks about several other methods: https://raphlinus.github.io/graphics/curves/2019/12/23/flatten-quadbez.html )
(* In this program, one of the two curves is "flattened" (converted to
a piecewise linear approximation). Then the problem is reduced to
finding intersections of the other curve with line segments.
I have never seen this method published in the literature, but
somewhere saw it hinted at.
Mainly to increase awareness of the representation, I flatten the
one curve using the symmetric power polynomial basis. See
J. Sánchez-Reyes, ‘The symmetric analogue of the polynomial power
basis’, ACM Transactions on Graphics, vol 16 no 3, July 1997,
page 319.
J. Sánchez-Reyes, ‘Applications of the polynomial s-power basis
in geometry processing’, ACM Transactions on Graphics, vol 19
no 1, January 2000, page 35. *)
#include "share/atspre_staload.hats"
%{^
#include <math.h>
%}
(* One simple way to make a foreign function call. I want to use only
the ATS prelude, but the prelude does not include support for the C
math library. (The bundled libats/libc does, and separately
available ats2-xprelude does.) *)
extern fn sqrt : double -<> double = "mac#sqrt"
macdef huge_val = $extval (double, "HUGE_VAL")
#define NIL list_nil ()
#define :: list_cons
fun eval_bernstein_degree2
(@(q0 : double,
q1 : double,
q2 : double),
t : double)
: double =
let
(* The de Casteljau algorithm. (The Schumaker-Volk algorithm also
is good BTW and is faster. In this program it should make no
noticeable difference, however.) *)
val s = 1.0 - t
val q01 = (s * q0) + (t * q1)
val q12 = (s * q1) + (t * q2)
val q012 = (s * q01) + (t * q12)
in
q012
end
(* @(...) means an unboxed tuple. Also often can be written without
the @, but then might be mistaken for argument parentheses. *)
fun
bernstein2spower_degree2
(@(c0 : double, c1 : double, c2 : double))
: @(double, double, double) =
(* Convert from Bernstein coefficients (control points) to symmetric
power coefficients. *)
@(c0, c1 + c1 - c0 - c2, c2)
fun
spower_portion_degree2
(@(c0 : double, c1 : double, c2 : double),
@(t0 : double, t1 : double))
: @(double, double, double) =
(* Compose spower(c0, c1, c2) with spower(t0, t1). This will map the
portion [t0,t1] onto [0,1]. (I got these expressions with
Maxima, a while back.) *)
let
val t0_t0 = t0 * t0
and t0_t1 = t0 * t1
and t1_t1 = t1 * t1
and c2p1m0 = c2 + c1 - c0
val d0 = c0 + (c2p1m0 * t0) - (c1 * t0_t0)
and d1 = (c1 * t1_t1) - ((c1 + c1) * t0_t1) + (c1 * t0_t0)
and d2 = c0 + (c2p1m0 * t1) - (c1 * t1_t1)
in
@(d0, d1, d2)
end
fun
solve_linear_quadratic
(@(px0 : double, px1 : double),
@(py0 : double, py1 : double),
@(qx0 : double, qx1 : double, qx2 : double),
@(qy0 : double, qy1 : double, qy2 : double))
(* Returns the two real roots, or any numbers outside [0,1], if
there are no real roots. *)
: @(double, double) =
let
(* The coefficients of the quadratic equation can be found by the
following Maxima commands, which implicitize the line segment
and plug in the parametric equations of the parabola:
/* The line. */
xp(t) := px0*(1-t) + px1*t$
yp(t) := py0*(1-t) + py1*t$
/* The quadratic (Bernstein basis). */
xq(t) := qx0*(1-t)**2 + 2*qx1*t*(1-t) + qx2*t**2$
yq(t) := qy0*(1-t)**2 + 2*qy1*t*(1-t) + qy2*t**2$
/* Implicitize and plug in. */
impl(t) := resultant(xq(t)-xp(tau), yq(t)-yp(tau), tau)$
impl(t);
expand(impl(t));
Consequently you get a quadratic equation in t, which can be
solved by the quadratic formula.
Sometimes people solve this problem by projecting the line
segment onto the x- or y-axis, and similarly projecting the
parabola. However, the following is simpler to write, if you
have Maxima to derive it for you. Whether it is better to use
the expanded expression (as here) or not to, I do not know. *)
val px0py1 = px0 * py1
and px1py0 = px1 * py0
and px0qy0 = px0 * qy0
and px0qy1 = px0 * qy1
and px0qy2 = px0 * qy2
and px1qy0 = px1 * qy0
and px1qy1 = px1 * qy1
and px1qy2 = px1 * qy2
and py0qx0 = py0 * qx0
and py0qx1 = py0 * qx1
and py0qx2 = py0 * qx2
and py1qx0 = py1 * qx0
and py1qx1 = py1 * qx1
and py1qx2 = py1 * qx2
val A = ~px1qy2 + px0qy2 - px1qy0 + py1qx0
+ px0qy0 + py1qx2 - py0qx2 - py0qx0
+ 2.0 * (px1qy1 - px0qy1 - py1qx1 + py0qx1)
and B = 2.0 * (~px1qy1 + px0qy1 + px1qy0 - px0qy0
+ py1qx1 - py0qx1 - py1qx0 + py0qx0)
and C = ~px1qy0 + px0qy0 + py1qx0 - py0qx0 - px0py1 + px1py0
val discriminant = (B * B) - (4.0 * A * C)
in
if discriminant < g0i2f 0 then
@(huge_val, huge_val) (* No real solutions. *)
else
let
val sqrt_discr = sqrt (discriminant)
val t1 = (~B - sqrt_discr) / (A + A)
and t2 = (~B + sqrt_discr) / (A + A)
fn
check_t (t : double) : double =
(* The parameter must lie in [0,1], and the intersection
point must lie between (px0,py0) and (px1,py1). We will
check only the x coordinate. *)
if t < 0.0 || 1.0 < t then
huge_val
else
let
val x = eval_bernstein_degree2 (@(qx0, qx1, qx2), t)
in
if x < px0 || px1 < x then
huge_val
else
t
end
in
@(check_t t1, check_t t2)
end
end
fun
flat_enough (@(px0 : double,
px1 : double,
px2 : double),
@(py0 : double,
py1 : double,
py2 : double),
tol : double)
: bool =
(* The quadratic must be given in s-power coefficients. Its px1 and
py1 terms are to be removed. Compare an error estimate to the
segment length. *)
let
(*
The symmetric power polynomials of degree 2 are
1-t
t(1-t)
t
Conversion from quadratic to linear is effected by removal of
the center term, with absolute error bounded by the value of the
center coefficient, divided by 4 (because t(1-t) reaches a
maximum of 1/4, at t=1/2).
*)
val error_squared = 0.125 * ((px1 * px1) + (py1 * py1))
and length_squared = (px2 - px0)**2 + (py2 - py0)**2
in
error_squared / tol <= length_squared * tol
end
fun
find_intersection_parameters
(px : @(double, double, double),
py : @(double, double, double),
qx : @(double, double, double),
qy : @(double, double, double),
tol : double)
: [m : nat | m <= 4] list (double, m) =
let
val px = bernstein2spower_degree2 px
and py = bernstein2spower_degree2 py
fun
loop {n : nat | n <= 4}
(params : list (double, n),
n : int n,
workload : List0 (@(double, double)))
: [m : nat | m <= 4] list (double, m) =
if n = 4 then
params
else
case+ workload of
| NIL => params
| hd :: tl =>
let
val portionx = spower_portion_degree2 (px, hd)
and portiony = spower_portion_degree2 (py, hd)
in
if flat_enough (portionx, portiony, tol) then
let
val @(portionx0, _, portionx2) = portionx
and @(portiony0, _, portiony2) = portiony
val @(root0, root1) =
solve_linear_quadratic (@(portionx0, portionx2),
@(portiony0, portiony2),
qx, qy)
in
if 0.0 <= root0 && root0 <= 1.0 then
begin
if (n < 3) * (0.0 <= root1 && root1 <= 1.0) then
loop (root0 :: root1 :: params, n + 2, tl)
else
loop (root0 :: params, n + 1, tl)
end
else if 0.0 <= root1 && root1 <= 1.0 then
loop (root1 :: params, n + 1, tl)
else
loop (params, n, tl)
end
else
let
val @(t0, t1) = hd
val tmiddle = (0.5 * t0) + (0.5 * t1)
val job1 = @(t0, tmiddle)
and job2 = @(tmiddle, t1)
in
loop (params, n, job1 :: job2 :: tl)
end
end
in
loop (NIL, 0, @(0.0, 1.0) :: NIL)
end
implement
main0 () =
let
val px = @(~1.0, 0.0, 1.0)
val py = @(0.0, 10.0, 0.0)
val qx = @(2.0, ~8.0, 2.0)
val qy = @(1.0, 2.0, 3.0)
val tol = 0.001 (* "Flatness ratio" *)
val t_list = find_intersection_parameters (px, py, qx, qy, tol)
fun
loop {n : nat} .<n>.
(t_list : list (double, n))
: void =
case+ t_list of
| NIL => ()
| t :: tl =>
begin
println! ("(", eval_bernstein_degree2 (qx, t), ", ",
eval_bernstein_degree2 (qy, t), ")");
loop tl
end
in
loop t_list
end- Output:
(0.881023, 1.118975) (0.654983, 2.854983) (-0.681024, 2.681025) (-0.854982, 1.345016)
D
This program subdivides both curves by de Casteljau's algorithm, until only very small subdivisions with overlapping control polygons remain.
- Update: I have added a crude check against accidentally detecting the same intersection twice, similar to the check in the Modula-2 program. I also changed the value of
tol, so that the check sometimes comes out positive.
- A "crude" check seems to me appropriate for a floating-point algorithm such as this. Even so, in a practical application one might not wish to stop after four detections, since the algorithm also might detect "near-intersections".
// The control points of a planar quadratic Bézier curve form a
// triangle--called the "control polygon"--that completely contains
// the curve. Furthermore, the rectangle formed by the minimum and
// maximum x and y values of the control polygon completely contain
// the polygon, and therefore also the curve.
//
// Thus a simple method for narrowing down where intersections might
// be is: subdivide both curves until you find "small enough" regions
// where these rectangles overlap.
import std.algorithm;
import std.container.slist;
import std.math;
import std.range;
import std.stdio;
struct point
{
double x, y;
}
struct quadratic_spline // Non-parametric spline.
{
double c0, c1, c2;
}
struct quadratic_curve // Planar parametric spline.
{
quadratic_spline x, y;
}
void
subdivide_quadratic_spline (quadratic_spline q, double t,
ref quadratic_spline u,
ref quadratic_spline v)
{
// Subdivision by de Casteljau's algorithm.
immutable s = 1 - t;
immutable c0 = q.c0;
immutable c1 = q.c1;
immutable c2 = q.c2;
u.c0 = c0;
v.c2 = c2;
u.c1 = (s * c0) + (t * c1);
v.c1 = (s * c1) + (t * c2);
u.c2 = (s * u.c1) + (t * v.c1);
v.c0 = u.c2;
}
void
subdivide_quadratic_curve (quadratic_curve q, double t,
ref quadratic_curve u,
ref quadratic_curve v)
{
subdivide_quadratic_spline (q.x, t, u.x, v.x);
subdivide_quadratic_spline (q.y, t, u.y, v.y);
}
bool
rectangles_overlap (double xa0, double ya0, double xa1, double ya1,
double xb0, double yb0, double xb1, double yb1)
{
// It is assumed that xa0<=xa1, ya0<=ya1, xb0<=xb1, and yb0<=yb1.
return (xb0 <= xa1 && xa0 <= xb1 && yb0 <= ya1 && ya0 <= yb1);
}
void
test_intersection (quadratic_curve p, quadratic_curve q, double tol,
ref bool exclude, ref bool accept,
ref point intersection)
{
// I will not do a lot of checking for intersections, as one might
// wish to do in a particular application. If the boxes are small
// enough, I will accept the point as an intersection.
immutable pxmin = min (p.x.c0, p.x.c1, p.x.c2);
immutable pymin = min (p.y.c0, p.y.c1, p.y.c2);
immutable pxmax = max (p.x.c0, p.x.c1, p.x.c2);
immutable pymax = max (p.y.c0, p.y.c1, p.y.c2);
immutable qxmin = min (q.x.c0, q.x.c1, q.x.c2);
immutable qymin = min (q.y.c0, q.y.c1, q.y.c2);
immutable qxmax = max (q.x.c0, q.x.c1, q.x.c2);
immutable qymax = max (q.y.c0, q.y.c1, q.y.c2);
exclude = true;
accept = false;
if (rectangles_overlap (pxmin, pymin, pxmax, pymax,
qxmin, qymin, qxmax, qymax))
{
exclude = false;
immutable xmin = max (pxmin, qxmin);
immutable xmax = min (pxmax, qxmax);
assert (xmax >= xmin);
if (xmax - xmin <= tol)
{
immutable ymin = max (pymin, qymin);
immutable ymax = min (pymax, qymax);
assert (ymax >= ymin);
if (ymax - ymin <= tol)
{
accept = true;
intersection = point ((0.5 * xmin) + (0.5 * xmax),
(0.5 * ymin) + (0.5 * ymax));
}
}
}
}
bool
seems_to_be_a_duplicate (point[] intersections, point xy,
double spacing)
{
bool seems_to_be_dup = false;
int i = 0;
while (!seems_to_be_dup && i != intersections.length)
{
immutable pt = intersections[i];
seems_to_be_dup =
fabs (pt.x - xy.x) < spacing && fabs (pt.y - xy.y) < spacing;
i += 1;
}
return seems_to_be_dup;
}
point[]
find_intersections (quadratic_curve p, quadratic_curve q,
double tol, double spacing)
{
point[] intersections;
int num_intersections = 0;
struct workset
{
quadratic_curve p, q;
}
SList!workset workload;
// Initial workload.
workload.insertFront(workset (p, q));
// Quit looking /* after having found four intersections */ or
// emptied the workload.
while (/*num_intersections != 4 &&*/ !workload.empty)
{
auto work = workload.front;
workload.removeFront();
bool exclude;
bool accept;
point intersection;
test_intersection (work.p, work.q, tol, exclude, accept,
intersection);
if (accept)
{
// This is a crude way to avoid detecting the same
// intersection twice: require some space between
// intersections. For, say, font design work, this method
// should be adequate.
if (!seems_to_be_a_duplicate (intersections,
intersection, spacing))
{
intersections.length = num_intersections + 1;
intersections[num_intersections] = intersection;
num_intersections += 1;
}
}
else if (!exclude)
{
quadratic_curve p0, p1, q0, q1;
subdivide_quadratic_curve (work.p, 0.5, p0, p1);
subdivide_quadratic_curve (work.q, 0.5, q0, q1);
workload.insertFront(workset (p0, q0));
workload.insertFront(workset (p0, q1));
workload.insertFront(workset (p1, q0));
workload.insertFront(workset (p1, q1));
}
}
return intersections;
}
int
main ()
{
quadratic_curve p, q;
p.x.c0 = -1.0; p.x.c1 = 0.0; p.x.c2 = 1.0;
p.y.c0 = 0.0; p.y.c1 = 10.0; p.y.c2 = 0.0;
q.x.c0 = 2.0; q.x.c1 = -8.0; q.x.c2 = 2.0;
q.y.c0 = 1.0; q.y.c1 = 2.0; q.y.c2 = 3.0;
immutable tol = 0.0000001;
immutable spacing = 10 * tol;
auto intersections = find_intersections (p, q, tol, spacing);
for (int i = 0; i != intersections.length; i += 1)
printf("(%f, %f)\n", intersections[i].x, intersections[i].y);
return 0;
}
- Output:
(0.654983, 2.854983) (0.881025, 1.118975) (-0.681025, 2.681025) (-0.854983, 1.345017)
Go
/* The control points of a planar quadratic Bézier curve form a
triangle--called the "control polygon"--that completely contains
the curve. Furthermore, the rectangle formed by the minimum and
maximum x and y values of the control polygon completely contain
the polygon, and therefore also the curve.
Thus a simple method for narrowing down where intersections might
be is: subdivide both curves until you find "small enough" regions
where these rectangles overlap.
*/
package main
import (
"fmt"
"log"
"math"
)
type point struct {
x, y float64
}
type quadSpline struct { // Non-parametric spline.
c0, c1, c2 float64
}
type quadCurve struct { // Planar parametric spline.
x, y quadSpline
}
// Subdivision by de Casteljau's algorithm.
func subdivideQuadSpline(q quadSpline, t float64, u, v *quadSpline) {
s := 1.0 - t
c0 := q.c0
c1 := q.c1
c2 := q.c2
u.c0 = c0
v.c2 = c2
u.c1 = s*c0 + t*c1
v.c1 = s*c1 + t*c2
u.c2 = s*u.c1 + t*v.c1
v.c0 = u.c2
}
func subdivideQuadCurve(q quadCurve, t float64, u, v *quadCurve) {
subdivideQuadSpline(q.x, t, &u.x, &v.x)
subdivideQuadSpline(q.y, t, &u.y, &v.y)
}
// It is assumed that xa0 <= xa1, ya0 <= ya1, xb0 <= xb1, and yb0 <= yb1.
func rectsOverlap(xa0, ya0, xa1, ya1, xb0, yb0, xb1, yb1 float64) bool {
return (xb0 <= xa1 && xa0 <= xb1 && yb0 <= ya1 && ya0 <= yb1)
}
func max3(x, y, z float64) float64 { return math.Max(math.Max(x, y), z) }
func min3(x, y, z float64) float64 { return math.Min(math.Min(x, y), z) }
// This accepts the point as an intersection if the boxes are small enough.
func testIntersect(p, q quadCurve, tol float64, exclude, accept *bool, intersect *point) {
pxmin := min3(p.x.c0, p.x.c1, p.x.c2)
pymin := min3(p.y.c0, p.y.c1, p.y.c2)
pxmax := max3(p.x.c0, p.x.c1, p.x.c2)
pymax := max3(p.y.c0, p.y.c1, p.y.c2)
qxmin := min3(q.x.c0, q.x.c1, q.x.c2)
qymin := min3(q.y.c0, q.y.c1, q.y.c2)
qxmax := max3(q.x.c0, q.x.c1, q.x.c2)
qymax := max3(q.y.c0, q.y.c1, q.y.c2)
*exclude = true
*accept = false
if rectsOverlap(pxmin, pymin, pxmax, pymax, qxmin, qymin, qxmax, qymax) {
*exclude = false
xmin := math.Max(pxmin, qxmin)
xmax := math.Min(pxmax, pxmax)
if xmax < xmin {
log.Fatalf("Assertion failure: %f < %f\n", xmax, xmin)
}
if xmax-xmin <= tol {
ymin := math.Max(pymin, qymin)
ymax := math.Min(pymax, qymax)
if ymax < ymin {
log.Fatalf("Assertion failure: %f < %f\n", ymax, ymin)
}
if ymax-ymin <= tol {
*accept = true
intersect.x = 0.5*xmin + 0.5*xmax
intersect.y = 0.5*ymin + 0.5*ymax
}
}
}
}
func findIntersects(p, q quadCurve, tol float64) []point {
var intersects []point
type workset struct {
p, q quadCurve
}
workload := []workset{workset{p, q}}
// Quit looking after having found four intersections or emptied the workload.
for len(intersects) != 4 && len(workload) > 0 {
l := len(workload)
work := workload[l-1]
workload = workload[0 : l-1]
var exclude, accept bool
intersect := point{0, 0}
testIntersect(work.p, work.q, tol, &exclude, &accept, &intersect)
if accept {
intersects = append(intersects, intersect)
} else if !exclude {
var p0, p1, q0, q1 quadCurve
subdivideQuadCurve(work.p, 0.5, &p0, &p1)
subdivideQuadCurve(work.q, 0.5, &q0, &q1)
workload = append(workload, workset{p0, q0})
workload = append(workload, workset{p0, q1})
workload = append(workload, workset{p1, q0})
workload = append(workload, workset{p1, q1})
}
}
return intersects
}
func main() {
var p, q quadCurve
p.x = quadSpline{-1.0, 0.0, 1.0}
p.y = quadSpline{0.0, 10.0, 0.0}
q.x = quadSpline{2.0, -8.0, 2.0}
q.y = quadSpline{1.0, 2.0, 3.0}
intersects := findIntersects(p, q, 0.000001)
for _, intersect := range intersects {
fmt.Printf("(% f, %f)\n", intersect.x, intersect.y)
}
}
- Output:
( 0.654983, 2.854984) ( 0.881025, 1.118975) (-0.681025, 2.681025) (-0.854984, 1.345017)
Maxima
This Maxima batch script finds an implicit equation for one of the curves, plugs the parametric equations of the other curve into the implicit equation, and then solves the resulting quartic equation.
/*
The method of implicitization:
1. Find an implicit equation for one of the curves, in x and y.
2. Plug the parametric equations of the other curve into the implicit
equation.
3. Find the roots of the resulting polynomial equation in t.
4. Plug those roots into the parametric equations of step (2).
*/
/* The Bernstein basis of degree 2. See
https://en.wikipedia.org/w/index.php?title=Bernstein_polynomial&oldid=1144565695
*/
b02(t) := 1 - 2*t + t**2$
b12(t) := 2*t - 2*t**2$
b22(t) := t**2$
/* The convex-up parabola, with its control points as coefficients of
the Bernstein basis. */
xu(t) := -b02(t) + b22(t)$
yu(t) := 10*b12(t)$
/* The convex-left parabola, with its control points as coefficients
of the Bernstein basis. */
xl(t) := 2*b02(t) - 8*b12(t) + 2*b22(t)$
yl(t) := b02(t) + 2*b12(t) + 3*b22(t)$
/* One can implicitize the convex-up Bézier curve by computing the
resultant of x - xu and y - yu.
The method is mentioned at
https://en.wikipedia.org/w/index.php?title=Gr%C3%B6bner_basis&oldid=1152603392#Implicitization_of_a_rational_curve
although they are describing a more general method that I do not
know how to do.
Here I combine forming the resultant with plugging in xl(t) and
yl(t). */
quartic_poly: resultant (xl(t) - xu(tau), yl(t) - yu(tau), tau)$
/* Find all the roots of the quartic equation that lie in [0,1]. */
roots: ev (realroots (quartic_poly = 0), float)$
roots: sublist(roots, lambda([item], 0 <= rhs(t) and rhs(t) <= 1))$
/* Plug them into xl(t) and yl(t). */
for i: 1 thru length(roots) do
block (
display(expand(xl(roots[i]))),
display(expand(yl(roots[i])))
)$
/* As an afterword, I would like to mention some drawbacks of
implicitization.
* It cannot find self-intersections. This is a major problem for
curves of degree 3 or greater.
* It gives you the t-parameter values for only one of the two
curves. If you just need t-parameter values for both curves
(such as to break them up at intersection points), then you
could perform implicitization both ways. But, if you need to
know which t corresponds to which, you need more than just
implicitization. (A method for finding t from given (x,y), for
instance.)
* It requires first constructing a polynomial of degree 4, 9, 16,
etc., and then finding its roots in [0,1]. There are serious
difficulties associated with both of those tasks. */
- Output:
(%i2) b02(t):=1-2*t+t^2
(%i3) b12(t):=2*t-2*t^2
(%i4) b22(t):=t^2
(%i5) xu(t):=-b02(t)+b22(t)
(%i6) yu(t):=10*b12(t)
(%i7) xl(t):=2*b02(t)-8*b12(t)+2*b22(t)
(%i8) yl(t):=b02(t)+2*b12(t)+3*b22(t)
(%i9) quartic_poly:resultant(xl(t)-xu(tau),yl(t)-yu(tau),tau)
(%i10) roots:ev(realroots(quartic_poly = 0),float)
(%i11) roots:sublist(roots,lambda([item],0 <= rhs(t) and rhs(t) <= 1))
(%i12) for i thru length(roots) do
block(display(expand(xl(roots[i]))),display(expand(yl(roots[i]))))
2
20 t - 20 t + 2 = 0.8810253968010677
2 t + 1 = 1.1189749836921692
2
20 t - 20 t + 2 = - 0.8549833297165428
2 t + 1 = 1.3450165390968323
2
20 t - 20 t + 2 = - 0.681024960552616
2 t + 1 = 2.681024968624115
2
20 t - 20 t + 2 = 0.6549829805579925
2 t + 1 = 2.854983389377594
Modula-2
Compile with the "-fiso" flag.
This program is similar to the D but, instead of using control points to form rectangles, uses values of the curves to form the rectangles. Instead of subdivision, there is function evaluation.
(* This program does not do any subdivision, but instead takes
advantage of monotonicity.
It is possible for points accidentally to be counted twice, for
instance if they lie right on an interval boundary. We will avoid
that by the VERY CRUDE (but perhaps often satisfactory) mechanism
of requiring a minimum max norm between intersections. *)
MODULE bezierIntersectionsInModula2;
(* ISO Modula-2 libraries. *)
FROM Storage IMPORT ALLOCATE, DEALLOCATE;
FROM SYSTEM IMPORT TSIZE;
IMPORT SLongIO;
IMPORT STextIO;
(* GNU Modula-2 gm2-libs *)
FROM Assertion IMPORT Assert;
(* Schumaker's and Volk's algorithm for evaluating a Bézier spline in
Bernstein basis. This is faster than de Casteljau, though not quite
as numerical stable. *)
PROCEDURE SchumakerVolk (c0, c1, c2, t : LONGREAL) : LONGREAL;
VAR s, u, v : LONGREAL;
BEGIN
s := 1.0 - t;
IF t <= 0.5 THEN
(* Horner form in the variable u = t/s, taking into account the
binomial coefficients = 1,2,1. *)
u := t / s;
v := c0 + (u * (c1 + c1 + (u * c2)));
(* Multiply by s raised to the degree of the spline. *)
v := v * s * s;
ELSE
(* Horner form in the variable u = s/t, taking into account the
binomial coefficients = 1,2,1. *)
u := s / t;
v := c2 + (u * (c1 + c1 + (u * c0)));
(* Multiply by t raised to the degree of the spline. *)
v := v * t * t;
END;
RETURN v;
END SchumakerVolk;
PROCEDURE FindExtremePoint (c0, c1, c2 : LONGREAL;
VAR LiesInside01 : BOOLEAN;
VAR ExtremePoint : LONGREAL);
VAR numer, denom : LONGREAL;
BEGIN
(* If the spline has c0=c2 but not c0=c1=c2, then treat it as having
an extreme point at 0.5. *)
IF (c0 = c2) AND (c0 <> c1) THEN
LiesInside01 := TRUE;
ExtremePoint := 0.5
ELSE
(* Find the root of the derivative of the spline. *)
LiesInside01 := FALSE;
numer := c0 - c1;
denom := c2 - c1 - c1 + c0;
IF (denom <> 0.0) AND (numer * denom >= 0.0)
AND (numer <= denom) THEN
LiesInside01 := TRUE;
ExtremePoint := numer / denom
END
END
END FindExtremePoint;
TYPE StartIntervCount = [2 .. 4];
StartIntervArray = ARRAY [1 .. 4] OF LONGREAL;
PROCEDURE PossiblyInsertExtremePoint
(c0, c1, c2 : LONGREAL;
VAR numStartInterv : StartIntervCount;
VAR startInterv : StartIntervArray);
VAR liesInside01 : BOOLEAN;
extremePt : LONGREAL;
BEGIN
FindExtremePoint (c0, c1, c2, liesInside01, extremePt);
IF liesInside01 AND (0.0 < extremePt) AND (extremePt < 1.0) THEN
IF numStartInterv = 2 THEN
startInterv[3] := 1.0;
startInterv[2] := extremePt;
numStartInterv := 3
ELSIF extremePt < startInterv[2] THEN
startInterv[4] := 1.0;
startInterv[3] := startInterv[2];
startInterv[2] := extremePt;
numStartInterv := 4
ELSIF extremePt > startInterv[2] THEN
startInterv[4] := 1.0;
startInterv[3] := extremePt;
numStartInterv := 4
END
END
END PossiblyInsertExtremePoint;
PROCEDURE minimum2 (x, y : LONGREAL) : LONGREAL;
VAR w : LONGREAL;
BEGIN
IF x <= y THEN
w := x
ELSE
w := y
END;
RETURN w;
END minimum2;
PROCEDURE maximum2 (x, y : LONGREAL) : LONGREAL;
VAR w : LONGREAL;
BEGIN
IF x >= y THEN
w := x
ELSE
w := y
END;
RETURN w;
END maximum2;
PROCEDURE RectanglesOverlap (xa0, ya0, xa1, ya1 : LONGREAL;
xb0, yb0, xb1, yb1 : LONGREAL) : BOOLEAN;
BEGIN
(* It is assumed that xa0<=xa1, ya0<=ya1, xb0<=xb1, and yb0<=yb1. *)
RETURN ((xb0 <= xa1) AND (xa0 <= xb1)
AND (yb0 <= ya1) AND (ya0 <= yb1))
END RectanglesOverlap;
PROCEDURE TestIntersection (xp0, xp1 : LONGREAL;
yp0, yp1 : LONGREAL;
xq0, xq1 : LONGREAL;
yq0, yq1 : LONGREAL;
tol : LONGREAL;
VAR exclude, accept : BOOLEAN;
VAR x, y : LONGREAL);
VAR xpmin, ypmin, xpmax, ypmax : LONGREAL;
xqmin, yqmin, xqmax, yqmax : LONGREAL;
xmin, xmax, ymin, ymax : LONGREAL;
BEGIN
xpmin := minimum2 (xp0, xp1);
ypmin := minimum2 (yp0, yp1);
xpmax := maximum2 (xp0, xp1);
ypmax := maximum2 (yp0, yp1);
xqmin := minimum2 (xq0, xq1);
yqmin := minimum2 (yq0, yq1);
xqmax := maximum2 (xq0, xq1);
yqmax := maximum2 (yq0, yq1);
exclude := TRUE;
accept := FALSE;
IF RectanglesOverlap (xpmin, ypmin, xpmax, ypmax,
xqmin, yqmin, xqmax, yqmax) THEN
exclude := FALSE;
xmin := maximum2 (xpmin, xqmin);
xmax := minimum2 (xpmax, xqmax);
Assert (xmax >= xmin);
IF xmax - xmin <= tol THEN
ymin := maximum2 (ypmin, yqmin);
ymax := minimum2 (ypmax, yqmax);
Assert (ymax >= ymin);
IF ymax - ymin <= tol THEN
accept := TRUE;
x := (0.5 * xmin) + (0.5 * xmax);
y := (0.5 * ymin) + (0.5 * ymax);
END;
END;
END;
END TestIntersection;
TYPE WorkPile = POINTER TO WorkTask;
WorkTask =
RECORD
tp0, tp1 : LONGREAL;
tq0, tq1 : LONGREAL;
next : WorkPile
END;
PROCEDURE WorkIsDone (workload : WorkPile) : BOOLEAN;
BEGIN
RETURN workload = NIL
END WorkIsDone;
PROCEDURE DeferWork (VAR workload : WorkPile;
tp0, tp1 : LONGREAL;
tq0, tq1 : LONGREAL);
VAR work : WorkPile;
BEGIN
ALLOCATE (work, TSIZE (WorkTask));
work^.tp0 := tp0;
work^.tp1 := tp1;
work^.tq0 := tq0;
work^.tq1 := tq1;
work^.next := workload;
workload := work
END DeferWork;
PROCEDURE DoSomeWork (VAR workload : WorkPile;
VAR tp0, tp1 : LONGREAL;
VAR tq0, tq1 : LONGREAL);
VAR work : WorkPile;
BEGIN
Assert (NOT WorkIsDone (workload));
work := workload;
tp0 := work^.tp0;
tp1 := work^.tp1;
tq0 := work^.tq0;
tq1 := work^.tq1;
workload := work^.next;
DEALLOCATE (work, TSIZE (WorkTask));
END DoSomeWork;
CONST px0 = -1.0; px1 = 0.0; px2 = 1.0;
py0 = 0.0; py1 = 10.0; py2 = 0.0;
qx0 = 2.0; qx1 = -8.0; qx2 = 2.0;
qy0 = 1.0; qy1 = 2.0; qy2 = 3.0;
tol = 0.0000001;
spacing = 100.0 * tol;
TYPE IntersectionCount = [0 .. 4];
IntersectionRange = [1 .. 4];
VAR pxHasExtremePt, pyHasExtremePt : BOOLEAN;
qxHasExtremePt, qyHasExtremePt : BOOLEAN;
pxExtremePt, pyExtremePt : LONGREAL;
qxExtremePt, qyExtremePt : LONGREAL;
pNumStartInterv, qNumStartInterv : StartIntervCount;
pStartInterv, qStartInterv : StartIntervArray;
workload : WorkPile;
i, j : StartIntervCount;
numIntersections, k : IntersectionCount;
intersectionsX : ARRAY IntersectionRange OF LONGREAL;
intersectionsY : ARRAY IntersectionRange OF LONGREAL;
tp0, tp1, tq0, tq1 : LONGREAL;
xp0, xp1, xq0, xq1 : LONGREAL;
yp0, yp1, yq0, yq1 : LONGREAL;
exclude, accept : BOOLEAN;
x, y : LONGREAL;
tpMiddle, tqMiddle : LONGREAL;
PROCEDURE MaybeAddIntersection (x, y : LONGREAL;
spacing : LONGREAL);
VAR i : IntersectionRange;
VAR TooClose : BOOLEAN;
BEGIN
IF numIntersections = 0 THEN
intersectionsX[1] := x;
intersectionsY[1] := y;
numIntersections := 1;
ELSE
TooClose := FALSE;
FOR i := 1 TO numIntersections DO
IF (ABS (x - intersectionsX[i]) < spacing)
AND (ABS (y - intersectionsY[i]) < spacing) THEN
TooClose := TRUE
END
END;
IF NOT TooClose THEN
numIntersections := numIntersections + 1;
intersectionsX[numIntersections] := x;
intersectionsY[numIntersections] := y
END
END
END MaybeAddIntersection;
BEGIN
(* Find monotonic sections of the curves, and use those as the
starting jobs. *)
pNumStartInterv := 2;
pStartInterv[1] := 0.0; pStartInterv[2] := 1.0;
PossiblyInsertExtremePoint (px0, px1, px2,
pNumStartInterv, pStartInterv);
PossiblyInsertExtremePoint (py0, py1, py2,
pNumStartInterv, pStartInterv);
qNumStartInterv := 2;
qStartInterv[1] := 0.0; qStartInterv[2] := 1.0;
PossiblyInsertExtremePoint (qx0, qx1, qx2,
qNumStartInterv, qStartInterv);
PossiblyInsertExtremePoint (qy0, qy1, qy2,
qNumStartInterv, qStartInterv);
workload := NIL;
FOR i := 2 TO pNumStartInterv DO
FOR j := 2 TO qNumStartInterv DO
DeferWork (workload, pStartInterv[i - 1], pStartInterv[i],
qStartInterv[j - 1], qStartInterv[j])
END;
END;
(* Go through the workload, deferring work as necessary. *)
numIntersections := 0;
WHILE (numIntersections <> 4) AND (NOT WorkIsDone (workload)) DO
(* The following code recomputes values of the splines
sometimes. You may wish to store such values in the work pile,
to avoid recomputing them. *)
DoSomeWork (workload, tp0, tp1, tq0, tq1);
xp0 := SchumakerVolk (px0, px1, px2, tp0);
yp0 := SchumakerVolk (py0, py1, py2, tp0);
xp1 := SchumakerVolk (px0, px1, px2, tp1);
yp1 := SchumakerVolk (py0, py1, py2, tp1);
xq0 := SchumakerVolk (qx0, qx1, qx2, tq0);
yq0 := SchumakerVolk (qy0, qy1, qy2, tq0);
xq1 := SchumakerVolk (qx0, qx1, qx2, tq1);
yq1 := SchumakerVolk (qy0, qy1, qy2, tq1);
TestIntersection (xp0, xp1, yp0, yp1,
xq0, xq1, yq0, yq1, tol,
exclude, accept, x, y);
IF accept THEN
MaybeAddIntersection (x, y, spacing)
ELSIF NOT exclude THEN
tpMiddle := (0.5 * tp0) + (0.5 * tp1);
tqMiddle := (0.5 * tq0) + (0.5 * tq1);
DeferWork (workload, tp0, tpMiddle, tq0, tqMiddle);
DeferWork (workload, tp0, tpMiddle, tqMiddle, tq1);
DeferWork (workload, tpMiddle, tp1, tq0, tqMiddle);
DeferWork (workload, tpMiddle, tp1, tqMiddle, tq1);
END
END;
IF numIntersections = 0 THEN
STextIO.WriteString ("no intersections");
STextIO.WriteLn;
ELSE
FOR k := 1 TO numIntersections DO
STextIO.WriteString ("(");
SLongIO.WriteReal (intersectionsX[k], 10);
STextIO.WriteString (", ");
SLongIO.WriteReal (intersectionsY[k], 10);
STextIO.WriteString (")");
STextIO.WriteLn;
END
END
END bezierIntersectionsInModula2.
- Output:
(0.65498343, 2.85498342) (0.88102499, 1.11897501) (-0.6810249, 2.68102500) (-0.8549834, 1.34501657)
Phix
enum X,Y type quadratic_spline(sequence c) -- return apply(c,sq_atom)={1,1,1} -- oops, requires 1.0.3... return true -- this will do instead for 1.0.2 and earlier end type type quadratic_curve(sequence c) -- return apply(c,sq_atom)={{1,1,1},{1,1,1}} -- ditto return true end type function subdivide_quadratic_spline(quadratic_spline q, atom t) // Subdivision by de Casteljau's algorithm. atom {c0,c1,c2} = q, s = 1 - t, u1 = (s * c0) + (t * c1), v1 = (s * c1) + (t * c2), m = (s * u1) + (t * v1) return {{c0,u1,m},{m,v1,c2}} end function function subdivide_quadratic_curve(quadratic_curve q, atom t) sequence {px,qx} = subdivide_quadratic_spline(q[X], t), {py,qy} = subdivide_quadratic_spline(q[Y], t) return {{px,py},{qx,qy}} end function function rectangles_overlap(atom xa0, ya0, xa1, ya1, xb0, yb0, xb1, yb1) // It is assumed that xa0<=xa1, ya0<=ya1, xb0<=xb1, and yb0<=yb1. return xb0 <= xa1 and xa0 <= xb1 and yb0 <= ya1 and ya0 <= yb1 end function function test_intersection(quadratic_curve p, q, atom tolerance) atom pxmin = min(p[X]), pymin = min(p[Y]), pxmax = max(p[X]), pymax = max(p[Y]), qxmin = min(q[X]), qymin = min(q[Y]), qxmax = max(q[X]), qymax = max(q[Y]) if rectangles_overlap(pxmin, pymin, pxmax, pymax, qxmin, qymin, qxmax, qymax) then atom xmin = max(pxmin, qxmin), xmax = min(pxmax, qxmax); assert(xmax >= xmin); if xmax-xmin <= tolerance then atom ymin = max(pymin, qymin), ymax = min (pymax, qymax); assert(ymax >= ymin); if ymax-ymin <= tolerance then -- we found a suitable intersection! return {(xmin+xmax)/2,(ymin+ymax)/2} end if end if return true -- accept/further subdivide end if return false -- exclude end function function find_intersections(quadratic_curve p, q, atom tolerance) sequence intersections = {}, todo = {{p,q}} -- while length(intersections)!=4 and length(todo) do while length(todo) do -- (we may as well be thorough) {{p,q},todo} = {todo[1],todo[2..$]} object insect = test_intersection(p, q, tolerance) if sequence(insect) then intersections &= {insect} elsif insect then sequence {p1,p2} = subdivide_quadratic_curve(p,0.5), {q1,q2} = subdivide_quadratic_curve(q,0.5) todo &= {{p2,q2},{p2,q1},{p1,q2},{p1,q1}} end if end while return intersections end function quadratic_curve p = {{-1,0,1},{0,10,0}}, q = {{2,-8,2},{1,2,3}} sequence intersections = find_intersections(p, q, 0.000001) pp(intersections,{pp_Nest,1})
- Output:
{{0.6549830437,2.854983807},
{0.8810248375,1.118975163},
{-0.6810250282,2.681025028},
{-0.8549838066,1.345016956}}
Wren
/* The control points of a planar quadratic Bézier curve form a
triangle--called the "control polygon"--that completely contains
the curve. Furthermore, the rectangle formed by the minimum and
maximum x and y values of the control polygon completely contain
the polygon, and therefore also the curve.
Thus a simple method for narrowing down where intersections might
be is: subdivide both curves until you find "small enough" regions
where these rectangles overlap.
*/
import "./dynamic" for Struct
import "./trait" for ByRef
import "./math" for Math, Nums
import "./assert" for Assert
import "./seq" for Stack
import "./fmt" for Fmt
// Note that these are all mutable as we need to pass by reference.
var Point = Struct.create("Point", ["x", "y"])
var QuadSpline = Struct.create("QuadSpline", ["c0", "c1", "c2"]) // non-parametric
var QuadCurve = Struct.create("QuadCurve", ["x", "y"]) // planar parametric
var Workset = Struct.create("Workset", ["p", "q"])
// Subdivision by de Casteljau's algorithm
var subdivideQuadSpline = Fn.new { |q, t, u, v|
var s = 1 - t
var c0 = q.c0
var c1 = q.c1
var c2 = q.c2
u.c0 = c0
v.c2 = c2
u.c1 = s * c0 + t * c1
v.c1 = s * c1 + t * c2
u.c2 = s * u.c1 + t * v.c1
v.c0 = u.c2
}
var subdivideQuadCurve = Fn.new { |q, t, u, v|
subdivideQuadSpline.call(q.x, t, u.x, v.x)
subdivideQuadSpline.call(q.y, t, u.y, v.y)
}
// It is assumed that xa0 <= xa1, ya0 <= ya1, xb0 <= xb1, and yb0 <= yb1.
var rectsOverlap = Fn.new { |xa0, ya0, xa1, ya1, xb0, yb0, xb1, yb1|
return (xb0 <= xa1 && xa0 <= xb1 && yb0 <= ya1 && ya0 <= yb1)
}
// This accepts the point as an intersection if the boxes are small enough.
var testIntersect = Fn.new { |p, q, tol, exclude, accept, intersect|
var pxmin = Nums.min([p.x.c0, p.x.c1, p.x.c2])
var pymin = Nums.min([p.y.c0, p.y.c1, p.y.c2])
var pxmax = Nums.max([p.x.c0, p.x.c1, p.x.c2])
var pymax = Nums.max([p.y.c0, p.y.c1, p.y.c2])
var qxmin = Nums.min([q.x.c0, q.x.c1, q.x.c2])
var qymin = Nums.min([q.y.c0, q.y.c1, q.y.c2])
var qxmax = Nums.max([q.x.c0, q.x.c1, q.x.c2])
var qymax = Nums.max([q.y.c0, q.y.c1, q.y.c2])
exclude.value = true
accept.value = false
if (rectsOverlap.call(pxmin, pymin, pxmax, pymax, qxmin, qymin, qxmax, qymax)) {
exclude.value = false
var xmin = Math.max(pxmin, qxmin)
var xmax = Math.min(pxmax, qxmax)
Assert.ok(xmax >= xmin)
if (xmax - xmin <= tol) {
var ymin = Math.max(pymin, qymin)
var ymax = Math.min(pymax, qymax)
Assert.ok(ymax >= ymin)
if (ymax - ymin <= tol) {
accept.value = true
intersect.x = 0.5 * xmin + 0.5 * xmax
intersect.y = 0.5 * ymin + 0.5 * ymax
}
}
}
}
var findIntersects = Fn.new { |p, q, tol|
var intersects = []
var workload = Stack.new()
workload.push(Workset.new(p, q))
// Quit looking after having found four intersections or emptied the workload.
while (intersects.count != 4 && !workload.isEmpty) {
var work = workload.peek()
workload.pop()
var exclude = ByRef.new(false)
var accept = ByRef.new(false)
var intersect = Point.new(0, 0)
testIntersect.call(work.p, work.q, tol, exclude, accept, intersect)
if (accept.value) {
intersects.add(intersect)
} else if (!exclude.value) {
var p0 = QuadCurve.new(QuadSpline.new(0, 0, 0), QuadSpline.new(0, 0, 0))
var p1 = QuadCurve.new(QuadSpline.new(0, 0, 0), QuadSpline.new(0, 0, 0))
var q0 = QuadCurve.new(QuadSpline.new(0, 0, 0), QuadSpline.new(0, 0, 0))
var q1 = QuadCurve.new(QuadSpline.new(0, 0, 0), QuadSpline.new(0, 0, 0))
subdivideQuadCurve.call(work.p, 0.5, p0, p1)
subdivideQuadCurve.call(work.q, 0.5, q0, q1)
workload.push(Workset.new(p0, q0))
workload.push(Workset.new(p0, q1))
workload.push(Workset.new(p1, q0))
workload.push(Workset.new(p1, q1))
}
}
return intersects
}
var p = QuadCurve.new(QuadSpline.new(-1, 0, 1), QuadSpline.new(0, 10, 0))
var q = QuadCurve.new(QuadSpline.new( 2, -8, 2), QuadSpline.new(1, 2, 3))
var intersects = findIntersects.call(p, q, 0.000001)
for (intersect in intersects) Fmt.print("($ f, $f)", intersect.x, intersect.y)- Output:
( 0.654983, 2.854984) ( 0.881025, 1.118975) (-0.681025, 2.681025) (-0.854984, 1.345017)