Bézier curves/Intersections: Difference between revisions
Content added Content deleted
(→{{header|Scheme}}: Assume the C definition of epsilon.) |
|||
| Line 3,715: | Line 3,715: | ||
=={{header|Scheme}}== |
=={{header|Scheme}}== |
||
===The algorithm of the [[#Icon|Icon]], [[#Python|Python]], [[#Pascal|Pascal]], etc.=== |
|||
{{trans|ObjectIcon}} |
{{trans|ObjectIcon}} |
||
Revision as of 15:11, 15 June 2023
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
(* One might be curious why "t@ype" instead of "type". The answer is:
the notation "type" is restricted to types that take up the same
space as a C void-pointer, which includes ATS pointers, "boxed"
types, etc. A "t@ype" can take up any amount of space, and so
includes any type there is (except for linear types, which is a
whole other subject). For instance, "int", "double", unboxed
records, unboxed tuples, and so on. *)
fun {a, b : t@ype} (* A polymorphic template function. *)
list_any (pred : (a, b) -<cloref1> bool,
obj : a,
lst : List0 b)
: bool =
(* Does pred(obj, item) return true for any list item? Here the
<cloref1> notation means that pred is a CLOSURE of the ordinary
garbage-collected kind, such as functions tend implicitly to be
in Lisps, MLs, Haskell, etc. *)
case+ lst of
| NIL => false
| hd :: tl =>
if pred (obj, hd) then
true
else
list_any (pred, obj, tl)
fun
find_intersection_parameters
(px : @(double, double, double),
py : @(double, double, double),
qx : @(double, double, double),
qy : @(double, double, double),
tol : double,
spacing : double)
: List0 double =
let
val px = bernstein2spower_degree2 px
and py = bernstein2spower_degree2 py
fun
within_spacing (t_candidate : double,
t_in_list : double)
:<cloref1> bool =
abs (t_candidate - t_in_list) < spacing
fun
loop {n : nat}
(params : list (double, n),
n : int n,
workload : List0 (@(double, double)))
: List0 double =
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 &&
~list_any (within_spacing, root0, params) then
begin
if 0.0 <= root1 && root1 <= 1.0 &&
~list_any (within_spacing, root1, params) then
loop (root0 :: root1 :: params, n + 2, tl)
else
loop (root0 :: params, n + 1, tl)
end
else if 0.0 <= root1 && root1 <= 1.0 &&
~list_any (within_spacing, root1, params) 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 spacing = 0.0001 (* Min. spacing between parameters. *)
val t_list = find_intersection_parameters (px, py, qx, qy,
tol, spacing)
(* For no particular reason, sort the intersections so they go
from top to bottom. *)
val t_list = list_vt2t (list_vt_reverse (list_mergesort t_list))
val () = println! ("From top to bottom:")
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:
From top to bottom: (0.654983, 2.854983) (-0.681024, 2.681025) (-0.854982, 1.345016) (0.881023, 1.118975)
C
C implementation 1
/* 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.
*/
#include <stdio.h>
#include <stdbool.h>
#include <math.h>
#include <assert.h>
typedef struct {
double x;
double y;
} point;
typedef struct {
double c0;
double c1;
double c2;
} quadSpline; // Non-parametric spline.
typedef struct {
quadSpline x;
quadSpline y;
} quadCurve; // Planar parametric spline.
// Subdivision by de Casteljau's algorithm.
void subdivideQuadSpline(quadSpline q, double t, quadSpline *u, quadSpline *v) {
double s = 1.0 - t;
double c0 = q.c0;
double c1 = q.c1;
double 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 subdivideQuadCurve(quadCurve q, double t, quadCurve *u, quadCurve *v) {
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.
bool rectsOverlap(double xa0, double ya0, double xa1, double ya1,
double xb0, double yb0, double xb1, double yb1) {
return (xb0 <= xa1 && xa0 <= xb1 && yb0 <= ya1 && ya0 <= yb1);
}
double max3(double x, double y, double z) {
return fmax(fmax(x, y), z);
}
double min3(double x, double y, double z) {
return fmin(fmin(x, y), z);
}
// This accepts the point as an intersection if the boxes are small enough.
void testIntersect(quadCurve p, quadCurve q, double tol,
bool *exclude, bool *accept, point *intersect) {
double pxmin = min3(p.x.c0, p.x.c1, p.x.c2);
double pymin = min3(p.y.c0, p.y.c1, p.y.c2);
double pxmax = max3(p.x.c0, p.x.c1, p.x.c2);
double pymax = max3(p.y.c0, p.y.c1, p.y.c2);
double qxmin = min3(q.x.c0, q.x.c1, q.x.c2);
double qymin = min3(q.y.c0, q.y.c1, q.y.c2);
double qxmax = max3(q.x.c0, q.x.c1, q.x.c2);
double 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;
double xmin = fmax(pxmin, qxmin);
double xmax = fmin(pxmax, qxmax);
assert(xmax >= xmin);
if (xmax - xmin <= tol) {
double ymin = fmax(pymin, qymin);
double ymax = fmin(pymax, qymax);
assert(ymax >= ymin);
if (ymax - ymin <= tol) {
*accept = true;
intersect->x = 0.5 * xmin + 0.5 * xmax;
intersect->y = 0.5 * ymin + 0.5 * ymax;
}
}
}
}
bool seemsToBeDuplicate(point intersects[], int icount, point xy, double spacing) {
bool seemsToBeDup = false;
int i = 0;
while (!seemsToBeDup && i != icount) {
point pt = intersects[i];
seemsToBeDup = fabs(pt.x - xy.x) < spacing && fabs(pt.y - xy.y) < spacing;
++i;
}
return seemsToBeDup;
}
void findIntersects(quadCurve p, quadCurve q, double tol, double spacing, point intersects[]) {
int numIntersects = 0;
typedef struct {
quadCurve p;
quadCurve q;
} workset;
workset workload[64];
int numWorksets = 1;
workload[0] = (workset){p, q};
// Quit looking after having emptied the workload.
while (numWorksets != 0) {
workset work = workload[numWorksets-1];
--numWorksets;
bool exclude, accept;
point intersect;
testIntersect(work.p, work.q, tol, &exclude, &accept, &intersect);
if (accept) {
// To avoid detecting the same intersection twice, require some
// space between intersections.
if (!seemsToBeDuplicate(intersects, numIntersects, intersect, spacing)) {
intersects[numIntersects++] = intersect;
assert(numIntersects <= 4);
}
} else if (!exclude) {
quadCurve p0, p1, q0, q1;
subdivideQuadCurve(work.p, 0.5, &p0, &p1);
subdivideQuadCurve(work.q, 0.5, &q0, &q1);
workload[numWorksets++] = (workset){p0, q0};
workload[numWorksets++] = (workset){p0, q1};
workload[numWorksets++] = (workset){p1, q0};
workload[numWorksets++] = (workset){p1, q1};
assert(numWorksets <= 64);
}
}
}
int main() {
quadCurve p, q;
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};
double tol = 0.0000001;
double spacing = tol * 10.0;
point intersects[4];
findIntersects(p, q, tol, spacing, intersects);
int i;
for (i = 0; i < 4; ++i) {
printf("(% f, %f)\n", intersects[i].x, intersects[i].y);
}
return 0;
}
- Output:
( 0.654983, 2.854983) ( 0.881025, 1.118975) (-0.681025, 2.681025) (-0.854983, 1.345017)
C implementation 2
Unfortunately two of us were writing C implementations at the same time. Had I known this, I would have written the following in a different language.
This implementation uses a recursive function rather than a "workload container". (Another approach to adaptive bisection requires no extra storage at all. That is to work with integer values in a clever way.)
// If you are using GCC, compile with -std=gnu2x because there may be
// C23-isms: [[attributes]], empty () instead of (void), etc.
/* In this program, both of the curves are adaptively "flattened":
that is, converted to a piecewise linear approximation. Then the
problem is reduced to finding intersections of line segments.
How efficient or inefficient the method is I will not try to
answer. (And I do sometimes compute things "too often", although a
really good optimizer might fix that.)
I will use the symmetric power basis that was introduced by
J. Sánchez-Reyes:
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.
Flattening a quadratic that is represented in this basis has a few
advantages, which I will not go into here. */
#include <stdio.h>
#include <stdbool.h>
#include <math.h>
static inline void
do_nothing ()
{
}
struct bernstein_spline
{
double b0;
double b1;
double b2;
};
struct spower_spline
{
double c0;
double c1;
double c2;
};
typedef struct bernstein_spline bernstein_spline;
typedef struct spower_spline spower_spline;
struct spower_curve
{
spower_spline x;
spower_spline y;
};
typedef struct spower_curve spower_curve;
// Convert a non-parametric spline from Bernstein basis to s-power.
[[gnu::const]] spower_spline
bernstein_spline_to_spower (bernstein_spline S)
{
spower_spline T =
{
.c0 = S.b0,
.c1 = (2 * S.b1) - S.b0 - S.b2,
.c2 = S.b2
};
return T;
}
// Compose (c0, c1, c2) with (t0, t1). This will map the portion
// [t0,t1] onto [0,1]. (To get these expressions, I did not use the
// general-degree methods described by Sánchez-Reyes, but instead used
// Maxima, some while ago.)
//
// This method is an alternative to de Casteljau subdivision, and can
// be done with the coefficients in any basis. Instead of breaking the
// spline into two pieces at a parameter value t, it gives you the
// portion lying between two parameter values. In general that
// requires two applications of de Casteljau subdivision. On the other
// hand, subdivision requires two applications of the following.
[[gnu::const]] inline spower_spline
spower_spline_portion (spower_spline S, double t0, double t1)
{
double t0_t0 = t0 * t0;
double t0_t1 = t0 * t1;
double t1_t1 = t1 * t1;
double c2p1m0 = S.c2 + S.c1 - S.c0;
spower_spline T =
{
.c0 = S.c0 + (c2p1m0 * t0) - (S.c1 * t0_t0),
.c1 = (S.c1 * t1_t1) - (2 * S.c1 * t0_t1) + (S.c1 * t0_t0),
.c2 = S.c0 + (c2p1m0 * t1) - (S.c1 * t1_t1)
};
return T;
}
[[gnu::const]] inline spower_curve
spower_curve_portion (spower_curve C, double t0, double t1)
{
spower_curve D =
{
.x = spower_spline_portion (C.x, t0, t1),
.y = spower_spline_portion (C.y, t0, t1)
};
return D;
}
// Given a parametric curve, is it "flat enough" to have its quadratic
// terms removed?
[[gnu::const]] inline bool
flat_enough (spower_curve C, double tol)
{
// The degree-2 s-power polynomials are 1-t, t(1-t), t. We want to
// remove the terms in t(1-t). The maximum of t(1-t) is 1/4, reached
// at t=1/2. That accounts for the 1/8=0.125 in the following:
double cx0 = C.x.c0;
double cx1 = C.x.c1;
double cx2 = C.x.c2;
double cy0 = C.y.c0;
double cy1 = C.y.c1;
double cy2 = C.y.c2;
double dx = cx2 - cx0;
double dy = cy2 - cy0;
double error_squared = 0.125 * ((cx1 * cx1) + (cy1 * cy1));
double length_squared = (dx * dx) + (dy * dy);
return (error_squared <= length_squared * tol * tol);
}
// Given two line segments, do they intersect? One solution to this
// problem is to use the implicitization method employed in the Maxima
// example, except to do it with linear instead of quadratic
// curves. That is what I do here, with the the roles of who gets
// implicitized alternated. If both ways you get as answer a parameter
// in [0,1], then the segments intersect.
void
test_line_segment_intersection (double ax0, double ax1,
double ay0, double ay1,
double bx0, double bx1,
double by0, double by1,
bool *they_intersect,
double *x, double *y)
{
double anumer = ((bx1 - bx0) * ay0 - (by1 - by0) * ax0
+ bx0 * by1 - bx1 * by0);
double bnumer = -((ax1 - ax0) * by0 - (ay1 - ay0) * bx0
+ ax0 * ay1 - ax1 * ay0);
double denom = ((ax1 - ax0) * (by1 - by0)
- (ay1 - ay0) * (bx1 - bx0));
double ta = anumer / denom; /* Parameter of segment a. */
double tb = bnumer / denom; /* Parameter of segment b. */
*they_intersect = (0 <= ta && ta <= 1 && 0 <= tb && tb <= 1);
if (*they_intersect)
{
*x = ((1 - ta) * ax0) + (ta * ax1);
*y = ((1 - ta) * ay0) + (ta * ay1);
}
}
bool
too_close (double x, double y, double xs[], double ys[],
size_t num_points, double spacing)
{
bool too_close = false;
size_t i = 0;
while (!too_close && i != num_points)
{
too_close = (fabs (x - xs[i]) < spacing
&& fabs (y - ys[i]) < spacing);
i += 1;
}
return too_close;
}
void
recursion (double tp0, double tp1, double tq0, double tq1,
spower_curve P, spower_curve Q,
double tol, double spacing, size_t max_points,
double xs[max_points], double ys[max_points],
size_t *num_points)
{
if (*num_points == max_points)
do_nothing ();
else if (!flat_enough (spower_curve_portion (P, tp0, tp1), tol))
{
double tp_half = (0.5 * tp0) + (0.5 * tp1);
if (!(flat_enough (spower_curve_portion (Q, tq0, tq1), tol)))
{
double tq_half = (0.5 * tq0) + (0.5 * tq1);
recursion (tp0, tp_half, tq0, tq_half, P, Q, tol,
spacing, max_points, xs, ys, num_points);
recursion (tp0, tp_half, tq_half, tq1, P, Q, tol,
spacing, max_points, xs, ys, num_points);
recursion (tp_half, tp1, tq0, tq_half, P, Q, tol,
spacing, max_points, xs, ys, num_points);
recursion (tp_half, tp1, tq_half, tq1, P, Q, tol,
spacing, max_points, xs, ys, num_points);
}
else
{
recursion (tp0, tp_half, tq0, tq1, P, Q, tol,
spacing, max_points, xs, ys, num_points);
recursion (tp_half, tp1, tq0, tq1, P, Q, tol,
spacing, max_points, xs, ys, num_points);
}
}
else if (!(flat_enough (spower_curve_portion (Q, tq0, tq1), tol)))
{
double tq_half = (0.5 * tq0) + (0.5 * tq1);
recursion (tp0, tp1, tq0, tq_half, P, Q, tol,
spacing, max_points, xs, ys, num_points);
recursion (tp0, tp1, tq_half, tq1, P, Q, tol,
spacing, max_points, xs, ys, num_points);
}
else
{
spower_curve P1 = spower_curve_portion (P, tp0, tp1);
spower_curve Q1 = spower_curve_portion (Q, tq0, tq1);
bool they_intersect;
double x, y;
test_line_segment_intersection (P1.x.c0, P1.x.c2,
P1.y.c0, P1.y.c2,
Q1.x.c0, Q1.x.c2,
Q1.y.c0, Q1.y.c2,
&they_intersect, &x, &y);
if (they_intersect &&
!too_close (x, y, xs, ys, *num_points, spacing))
{
xs[*num_points] = x;
ys[*num_points] = y;
*num_points += 1;
}
}
}
void
find_intersections (spower_curve P, spower_curve Q,
double flatness_tolerance,
double point_spacing,
size_t max_points,
double xs[max_points],
double ys[max_points],
size_t *num_points)
{
*num_points = 0;
recursion (0, 1, 0, 1, P, Q, flatness_tolerance, point_spacing,
max_points, xs, ys, num_points);
}
int
main ()
{
bernstein_spline bPx = { .b0 = -1, .b1 = 0, .b2 = 1 };
bernstein_spline bPy = { .b0 = 0, .b1 = 10, .b2 = 0 };
bernstein_spline bQx = { .b0 = 2, .b1 = -8, .b2 = 2 };
bernstein_spline bQy = { .b0 = 1, .b1 = 2, .b2 = 3 };
spower_spline Px = bernstein_spline_to_spower (bPx);
spower_spline Py = bernstein_spline_to_spower (bPy);
spower_spline Qx = bernstein_spline_to_spower (bQx);
spower_spline Qy = bernstein_spline_to_spower (bQy);
spower_curve P = { .x = Px, .y = Py };
spower_curve Q = { .x = Qx, .y = Qy };
double flatness_tolerance = 0.001;
double point_spacing = 0.000001; /* Max norm minimum spacing. */
const size_t max_points = 10;
double xs[max_points];
double ys[max_points];
size_t num_points;
find_intersections (P, Q, flatness_tolerance, point_spacing,
max_points, xs, ys, &num_points);
for (size_t i = 0; i != num_points; i += 1)
printf ("(%f, %f)\n", xs[i], ys[i]);
return 0;
}
- Output:
(-0.854982, 1.345017) (-0.681024, 2.681024) (0.881023, 1.118977) (0.654983, 2.854982)
D
This program subdivides both curves by de Casteljau's algorithm, until only very small subdivisions with overlapping control polygons remain. (You could use recursion instead of the workload container. With the container it is easier to terminate early, and also the program then uses only constant stack space.)
- 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 seemsToBeDuplicate(intersects []point, xy point, spacing float64) bool {
seemsToBeDup := false
i := 0
for !seemsToBeDup && i != len(intersects) {
pt := intersects[i]
seemsToBeDup = math.Abs(pt.x-xy.x) < spacing && math.Abs(pt.y-xy.y) < spacing
i++
}
return seemsToBeDup
}
func findIntersects(p, q quadCurve, tol, spacing float64) []point {
var intersects []point
type workset struct {
p, q quadCurve
}
workload := []workset{workset{p, q}}
// Quit looking after having emptied the workload.
for 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 {
// To avoid detecting the same intersection twice, require some
// space between intersections.
if !seemsToBeDuplicate(intersects, intersect, spacing) {
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}
tol := 0.0000001
spacing := tol * 10
intersects := findIntersects(p, q, tol, spacing)
for _, intersect := range intersects {
fmt.Printf("(% f, %f)\n", intersect.x, intersect.y)
}
}
- Output:
( 0.654983, 2.854983) ( 0.881025, 1.118975) (-0.681025, 2.681025) (-0.854983, 1.345017)
Icon
This implementation combines rectangle overlap (as in the Modula-2) and curve flattening (as in one of the C implementations), and tries to arrange the calculations efficiently. I think it may be a new algorithm worth serious consideration. Curve flattening alone requires too many line-segment-intersection checks, but rectangle overlap "prunes the tree".
Furthermore, the algorithm returns -parameter pairs, which is ideal.
Icon is a very good language in which to express this algorithm, if the audience can read Icon. But Icon is not designed to be readable by the programming public. See instead the Pascal. Pascal is meant to be readable. There also instead of a "workload" set there is a recursive procedure, which probably more programmers will recognize as "how to do adaptive bisection". (But both methods are common.)
See also:
- Output:
For each estimated intersection point, the program prints out the -parameters and the corresponding values of the curves.
convex up convex left 0.07250828117 (-0.8549834377, 1.345016607) 0.1725082997 (-0.8549837251, 1.345016599) 0.1594875309 (-0.6810249382, 2.681025168) 0.8405124691 (-0.681025168, 2.681024938) 0.8274917003 (0.6549834005, 2.854983725) 0.9274917188 (0.6549833933, 2.854983438) 0.9405124667 (0.8810249334, 1.118975334) 0.05948753331 (0.8810246662, 1.118975067)
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.
In theory, doing just the above could find an intersection that lies outside the parameter range of the implicitized curve. But we know all four points lie in that range. One could doublecheck by reversing the roles of the two curves.
/*
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 crude (but likely 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 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)
ObjectIcon
This program is mostly a translation of the Python.
Note that that the Icon program could itself be run as an Object Icon program, after a few necessary modifications. ("import io" would have to be added, the implementations of "abs" and "max" would have to be removed, and maybe a few other minor changes would be needed.)
See also:
# This is the algorithm of the Icon and Python implementations.
# Primarily it is a translation of the Python. Snippets are taken from
# the Icon.
#
# References on the symmetric power basis:
#
# 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.
#
import ipl.printf(printf)
procedure main ()
local flatness_tolerance
local minimum_spacing
local p, q
local tbl, params, tp, tq, ppoint, qpoint
local keys, k
flatness_tolerance := 0.0001
minimum_spacing := 0.000001
p := Curve.from_controls (Point (-1.0, 0.0),
Point (0.0, 10.0),
Point (1.0, 0.0))
q := Curve.from_controls (Point (2.0, 1.0),
Point (-8.0, 2.0),
Point ( 2.0, 3.0))
tbl := table ()
every params := find_intersections (p, q, flatness_tolerance) do
{
tp := params[1]; ppoint := p.eval(tp)
tq := params[2]; qpoint := q.eval(tq)
not (length ((!tbl)[2].x - ppoint.x,
(!tbl)[2].y - ppoint.y) < minimum_spacing) &
not (length ((!tbl)[4].x - qpoint.x,
(!tbl)[4].y - qpoint.y) < minimum_spacing) &
tbl[tp] := [tp, ppoint, tq, qpoint]
}
every put (keys := [], key (tbl))
keys := sort (keys)
printf ("\n")
printf (" convex up " ||
" convex left\n")
every k := !keys do
{
tp := tbl[k][1]
ppoint := tbl[k][2]
tq := tbl[k][3]
qpoint := tbl[k][4]
printf (" %11.8r (%11.8r, %11.8r) " ||
"%11.8r (%11.8r, %11.8r)\n",
tp, ppoint.x, ppoint.y, tq, qpoint.x, qpoint.y)
}
printf ("\n")
end
# Generate t-parameter pairs detected as corresponding to intersection
# points of p and q. There may be duplicate detections. It is assumed
# those will be weeded out by later processing. The tol parameter
# specifies the "flatness tolerance" and is a relative tolerance.
procedure find_intersections (p, q, tol)
local i, j
local tp, tq, pportion, qportion
local workload, workpair, params
# The initial workload is the cartesian product of the monotonic
# portions of p and q, respectively.
tp := [0.0] ||| sort (p.critical_points()) ||| [1.0]
tq := [0.0] ||| sort (q.critical_points()) ||| [1.0]
workload := set ()
every i := 1 to *tp - 1 do
every j := 1 to *tq - 1 do
{
pportion := Portion (p, tp[i], tp[i + 1],
p.eval(tp[i]), p.eval(tp[i + 1]))
qportion := Portion (q, tq[j], tq[j + 1],
q.eval(tq[j]), q.eval(tq[j + 1]))
insert (workload, [pportion, qportion])
}
while *workload ~= 0 do
{
workpair := ?workload
delete (workload, workpair)
pportion := workpair[1]
qportion := workpair[2]
if rectangles_overlap (pportion.endpt0, pportion.endpt1,
qportion.endpt0, qportion.endpt1) then
{
if pportion.flat_enough(tol) then
{
if qportion.flat_enough(tol) then
{
if params := segment_parameters(pportion.endpt0,
pportion.endpt1,
qportion.endpt0,
qportion.endpt1) then
{
tp := params[1]
tq := params[2]
tp := (1 - tp) * pportion.t0 + tp * pportion.t1
tq := (1 - tq) * qportion.t0 + tq * qportion.t1
suspend [tp, tq]
}
}
else
every insert (workload, [pportion, qportion.split()])
}
else
{
if qportion.flat_enough(tol) then
every insert (workload, [pportion.split(), qportion])
else
every insert (workload, [pportion.split(),
qportion.split()])
}
}
}
end
class Point ()
public const x, y
public new (x, y)
self.x := x
self.y := y
return
end
end
class SPower () # Non-parametric spline in s-power basis.
public const c0, c1, c2
private critpoints
public new (c0, c1, c2)
self.c0 := c0
self.c1 := c1
self.c2 := c2
return
end
# Evaluate at t.
public eval (t)
return (c0 + (c1 * t)) * (1.0 - t) + (c2 * t)
end
# Return the center coefficient for the [t0,t1] portion. (The other
# coefficients can be found with the eval method.)
public center_coef (t0, t1)
return c1 * ((t1 - t0 - t0) * t1 + (t0 * t0))
end
# Return a set of independent variable values for the critical
# points that lie in (0,1). (This is a memoizing implementation.)
public critical_points ()
local t
if /critpoints then
{
critpoints := set ()
if c1 ~= 0.0 then # If c1 is zero then the spline is linear.
{
if c0 = c2 then
insert (critpoints, 0.5) # The spline is "pulse-like".
else
{
# The root of the derivative is the critical point.
t = (0.5 * (c2 + c1 - c0)) / c1
0 < t < 1 | fail
insert (critpoints, t)
}
}
}
return critpoints
end
end
class Curve () # Parametric spline in s-power basis.
public const x, y
public new (x, y)
self.x := x
self.y := y
return
end
public static from_controls (ctl0, ctl1, ctl2)
local c0x, c0y, c1x, c1y, c2x, c2y
c0x := ctl0.x
c0y := ctl0.y
c1x := (2.0 * ctl1.x) - ctl0.x - ctl2.x
c1y := (2.0 * ctl1.y) - ctl0.y - ctl2.y
c2x := ctl2.x
c2y := ctl2.y
return Curve (SPower (c0x, c1x, c2x),
SPower (c0y, c1y, c2y))
end
# Evaluate at t.
public eval (t)
return Point (x.eval(t), y.eval(t))
end
# Return a set of t-parameter values for the critical points that
# lie in (0,1).
public critical_points ()
return (x.critical_points() ++ y.critical_points())
end
end
class Portion () # The [t0,t1]-portion of a parametric spline.
public const curve, t0, t1, endpt0, endpt1
public static default_num_pieces
private static init ()
default_num_pieces := 2
end
public new (curve, t0, t1, endpt0, endpt1)
self.curve := curve
self.t0 := t0
self.t1 := t1
self.endpt0 := endpt0
self.endpt1 := endpt1
return
end
# Is the Portion close enough to linear to be treated as a line
# segment?
public flat_enough (tol)
local xcentercoef, ycentercoef
local xlen, ylen
# The degree-2 s-power polynomials are 1-t, t(1-t), t. We want to
# remove the terms in t(1-t). The maximum of t(1-t) is 1/4,
# reached at t=1/2. That accounts for the 1/4=0.25 in the
# following.
xcentercoef := curve.x.center_coef(t0, t1)
ycentercoef := curve.y.center_coef(t0, t1)
xlen := endpt1.x - endpt0.x
ylen := endpt1.y - endpt0.y
return compare_lengths (0.25 * xcentercoef,
0.25 * ycentercoef,
tol * xlen, tol * ylen) <= 0
end
# Generate num_pieces (or default_num_pieces) sections of the
# Portion.
public split (num_pieces)
local k1, k, ts, vals, i
/num_pieces := default_num_pieces
k1 := 1.0 / num_pieces
every put (ts := [], (k := k1 * (1 to num_pieces - 1) &
(1.0 - k) * t0 + k * t1))
every put (vals := [], curve.eval(!ts))
ts := [t0] ||| ts ||| [t1]
vals := [endpt0] ||| vals ||| [endpt1]
every i := 1 to *ts - 1 do
suspend Portion (curve, ts[i], ts[i + 1], vals[i], vals[i + 1])
end
end
# Do the rectangles with corners at (ptA0,ptA1) and (ptB0,ptB1)
# overlap?
procedure rectangles_overlap (ptA0, ptA1, ptB0, ptB1)
local ax0, ay0, ax1, ay1
local bx0, by0, bx1, by1
ax0 := ptA0.x; ax1 := ptA1.x
bx0 := ptB0.x; bx1 := ptB1.x
if ax1 < ax0 then ax0 :=: ax1
if bx1 < bx0 then bx0 :=: bx1
bx1 < ax0 & fail
ax1 < bx0 & fail
ay0 := ptA0.y; ay1 := ptA1.y
by0 := ptB0.y; by1 := ptB1.y
if ay1 < ay0 then ay0 :=: ay1
if by1 < by0 then by0 :=: by1
by1 < ay0 & fail
ay1 < by0 & fail
return
end
# Return the respective [0,1] parameters of line segments
# (ptA0,ptA1) and (ptB0,ptB1), for their intersection point. Fail if
# there are not such parameters.
procedure segment_parameters (ptA0, ptA1, ptB0, ptB1)
local ax0, ax1, ay0, ay1
local bx0, bx1, by0, by1
local axdiff, aydiff
local bxdiff, bydiff
local anumer, bnumer, denom
local tA, tB
ax0 := ptA0.x; ax1 := ptA1.x
ay0 := ptA0.y; ay1 := ptA1.y
bx0 := ptB0.x; bx1 := ptB1.x
by0 := ptB0.y; by1 := ptB1.y
axdiff := ax1 - ax0
aydiff := ay1 - ay0
bxdiff := bx1 - bx0
bydiff := by1 - by0
denom := (axdiff * bydiff) - (aydiff * bxdiff)
anumer :=
(bxdiff * ay0) - (bydiff * ax0) + (bx0 * by1) - (bx1 * by0)
tA := anumer / denom
0 <= tA <= 1 | fail
bnumer :=
-((axdiff * by0) - (aydiff * bx0) + (ax0 * ay1) - (ax1 * ay0))
tB := bnumer / denom
0 <= tB <= 1 | fail
return [tA, tB]
end
# Length according to some norm, where (ax,ay) is a "measuring stick"
# vector. Here I use the max norm.
procedure length (a_x, a_y)
return max (abs (a_x), abs (a_y))
end
# Having a compare_lengths function lets one compare lengths in the
# euclidean metric by comparing the squares of the lengths, and thus
# avoiding the square root. The following, however, is a general
# implementation.
procedure compare_lengths (a_x, a_y, b_x, b_y)
local len_a, len_b, cmpval
len_a := length (a_x, a_y)
len_b := length (b_x, b_y)
if len_a < len_b then
cmpval := -1
else if len_a > len_b then
cmpval := 1
else
cmpval := 0
return cmpval
end- Output:
convex up convex left 0.07250828 (-0.85498344, 1.34501661) 0.17250830 (-0.85498373, 1.34501660) 0.15948753 (-0.68102494, 2.68102517) 0.84051247 (-0.68102517, 2.68102494) 0.82749170 ( 0.65498340, 2.85498373) 0.92749172 ( 0.65498339, 2.85498344) 0.94051247 ( 0.88102493, 1.11897533) 0.05948753 ( 0.88102467, 1.11897507)
Pascal
This is the algorithm of the Icon example, recast as a recursive procedure.
See also:
{$mode ISO} { Tell Free Pascal Compiler to use "ISO 7185" mode. }
{
This is the algorithm of the Icon example, recast as a recursive
procedure.
The "var" notation in formal parameter lists means pass by
reference. All other parameters are implicitly passed by value.
Pascal is case-insensitive.
In the old days, when Pascal was printed as a means to express
algorithms, it was usually in a fashion similar to Algol 60 reference
language. It was printed mostly in lowercase and did not have
underscores. Reserved words were in boldface and variables, etc., were
in italics. The effect was like that of Algol 60 reference language.
Code entry practices for Pascal were another matter. It may have been
all uppercase, with ML-style comment braces instead of squiggly
braces. It may have had uppercase reserved words and "Pascal case"
variables, etc., as one sees also in Modula-2 and Oberon-2 code.
Here I have deliberately adopted an all-lowercase style.
References on the s-power basis:
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.
}
program bezierintersections;
const
flatnesstolerance = 0.0001;
minimumspacing = 0.000001;
maxintersections = 10;
type
point =
record
x, y : real
end;
spower = { non-parametric spline in s-power basis }
record
c0, c1, c2 : real
end;
curve = { parametric spline in s-power basis }
record
x, y : spower
end;
portion = { portion of a parametric spline in [t0,t1] }
record
curv : curve;
t0, t1 : real;
endpt0, endpt1 : point { pre-computed for efficiency }
end;
intersectionscount = 0 .. maxintersections;
intersectionsrange = 1 .. maxintersections;
tparamsarray = array [intersectionsrange] of real;
coordsarray = array [intersectionsrange] of point;
var
numintersections : intersectionscount;
tparamsp : tparamsarray;
coordsp : coordsarray;
tparamsq : tparamsarray;
coordsq : coordsarray;
pglobal, qglobal : curve;
i : integer;
{ Minimum of two real. }
function rmin (x, y : real) : real;
begin
if x < y then rmin := x else rmin := y
end;
{ Maximum of two real. }
function rmax (x, y : real) : real;
begin
if x < y then rmax := y else rmax := x
end;
{ Insertion sort of an array of real. }
procedure realsort ( n : integer;
var a : array of real);
var
i, j : integer;
x : real;
done : boolean;
begin
i := low (a) + 1;
while i < n do
begin
x := a[i];
j := i - 1;
done := false;
while not done do
begin
if j + 1 = low (a) then
done := true
else if a[j] <= x then
done := true
else
begin
a[j + 1] := a[j];
j := j - 1
end
end;
a[j + 1] := x;
i := i + 1
end
end;
{ "Length" according to some definition. Here I use a max norm. The
"distance" between two points is the "length" of the differences of
the corresponding coordinates. (The sign of the difference should be
immaterial.) }
function length (ax, ay : real) : real;
begin
length := rmax (abs (ax), abs (ay))
end;
{ Having a "comparelengths" function makes it possible to use a
euclidean norm for "length", and yet avoid square roots. One
compares the squares of lengths, instead of the lengths
themselves. However, here I use a more general implementation. }
function comparelengths (ax, ay, bx, by : real) : integer;
var lena, lenb : real;
begin
lena := length (ax, ay);
lenb := length (bx, by);
if lena < lenb then
comparelengths := -1
else if lena > lenb then
comparelengths := 1
else
comparelengths := 0
end;
function makepoint (x, y : real) : point;
begin
makepoint.x := x;
makepoint.y := y
end;
function makeportion (curv : curve;
t0, t1 : real;
endpt0, endpt1 : point) : portion;
begin
makeportion.curv := curv;
makeportion.t0 := t0;
makeportion.t1 := t1;
makeportion.endpt0 := endpt0;
makeportion.endpt1 := endpt1;
end;
{ Convert from control points (that is, Bernstein basis) to the
symmetric power basis. }
function controlstospower (ctl0, ctl1, ctl2 : point) : curve;
begin
controlstospower.x.c0 := ctl0.x;
controlstospower.y.c0 := ctl0.y;
controlstospower.x.c1 := (2.0 * ctl1.x) - ctl0.x - ctl2.x;
controlstospower.y.c1 := (2.0 * ctl1.y) - ctl0.y - ctl2.y;
controlstospower.x.c2 := ctl2.x;
controlstospower.y.c2 := ctl2.y
end;
{ Evaluate an s-power spline at t. }
function spowereval (spow : spower;
t : real) : real;
begin
spowereval := (spow.c0 + (spow.c1 * t)) * (1.0 - t) + (spow.c2 * t)
end;
{ Evaluate a curve at t. }
function curveeval (curv : curve;
t : real) : point;
begin
curveeval.x := spowereval (curv.x, t);
curveeval.y := spowereval (curv.y, t)
end;
{ Return the center coefficient for the [t0,t1] portion of an s-power
spline. (The endpoint coefficients can be found with spowereval.) }
function spowercentercoef (spow : spower;
t0, t1 : real) : real;
begin
spowercentercoef := spow.c1 * ((t1 - t0 - t0) * t1 + (t0 * t0))
end;
{ Return t in (0,1) where spow is at a critical point, else return
-1.0. }
function spowercriticalpt (spow : spower) : real;
var t : real;
begin
spowercriticalpt := -1.0;
if spow.c1 <> 0.0 then { If c1 is zero, then the spline is linear. }
begin
if spow.c1 = spow.c2 then
spowercriticalpt := 0.5 { The spline is "pulse-like". }
else
begin
{ t = root of the derivative }
t := (spow.c2 + spow.c1 - spow.c0) / (spow.c1 + spow.c1);
if (0.0 < t) and (t < 1.0) then
spowercriticalpt := t
end
end
end;
{ Bisect a portion and pre-compute the new shared endpoint. }
procedure bisectportion ( port : portion;
var port1, port2 : portion);
begin
port1.curv := port.curv;
port2.curv := port.curv;
port1.t0 := port.t0;
port1.t1 := 0.5 * (port.t0 + port.t1);
port2.t0 := port1.t1;
port2.t1 := port.t1;
port1.endpt0 := port.endpt0;
port1.endpt1 := curveeval (port.curv, port1.t1);
port2.endpt0 := port1.endpt1;
port2.endpt1 := port.endpt1;
end;
{ Do the rectangles with corners at (a0,a1) and (b0,b1) overlap at
all? }
function rectanglesoverlap (a0, a1, b0, b1 : point) : boolean;
begin
rectanglesoverlap := ((rmin (a0.x, a1.x) <= rmax (b0.x, b1.x))
and (rmin (b0.x, b1.x) <= rmax (a0.x, a1.x))
and (rmin (a0.y, a1.y) <= rmax (b0.y, b1.y))
and (rmin (b0.y, b1.y) <= rmax (a0.y, a1.y)))
end;
{ Set the respective [0,1] parameters of line segments (a0,a1) and
(b0,b1), for their intersection point. If there are not two such
parameters, set both values to -1.0. }
procedure segmentparameters ( a0, a1, b0, b1 : point;
var ta, tb : real);
var
anumer, bnumer, denom : real;
axdiff, aydiff, bxdiff, bydiff : real;
begin
axdiff := a1.x - a0.x;
aydiff := a1.y - a0.y;
bxdiff := b1.x - b0.x;
bydiff := b1.y - b0.y;
denom := (axdiff * bydiff) - (aydiff * bxdiff);
anumer := ((bxdiff * a0.y) - (bydiff * a0.x)
+ (b0.x * b1.y) - (b1.x * b0.y));
ta := anumer / denom;
if (ta < 0.0) or (1.0 < ta) then
begin
ta := -1.0;
tb := -1.0
end
else
begin
bnumer := -((axdiff * b0.y) - (aydiff * b0.x)
+ (a0.x * a1.y) - (a1.x * a0.y));
tb := bnumer / denom;
if (tb < 0.0) or (1.0 < tb) then
begin
ta := -1.0;
tb := -1.0
end
end
end;
{ Is a curve portion flat enough to be treated as a line segment
between its endpoints? }
function flatenough (port : portion;
tol : real) : boolean;
var
xcentercoef, ycentercoef : real;
begin
{ The degree-2 s-power polynomials are 1-t, t(1-t), t. We want to
remove the terms in t(1-t). The maximum of t(1-t) is 1/4, reached
at t=1/2. That accounts for the 1/4=0.25 in the following. }
{ The "with" construct here is a shorthand to implicitly use fields
of the "port" record. Thus "curv.x" means "port.curv.x", etc. }
with port do
begin
xcentercoef := spowercentercoef (curv.x, t0, t1);
ycentercoef := spowercentercoef (curv.y, t0, t1);
flatenough := comparelengths (0.25 * xcentercoef,
0.25 * ycentercoef,
tol * (endpt1.x - endpt0.x),
tol * (endpt1.y - endpt0.y)) <= 0
end
end;
{ If the intersection point corresponding to tp and tq is not already
listed, insert it into the arrays, sorted by the value of tp. }
procedure insertintersection (p : curve;
tp : real;
q : curve;
tq : real);
var
ppoint, qpoint : point;
lenp, lenq : real;
i : intersectionscount;
insertionpoint : intersectionscount;
begin
if numintersections <> maxintersections then
begin
ppoint := curveeval (p, tp);
qpoint := curveeval (q, tq);
insertionpoint := numintersections + 1; { Insert at end. }
i := 0;
while (0 < insertionpoint) and (i <> numintersections) do
begin
i := i + 1;
lenp := length (coordsp[i].x - ppoint.x,
coordsp[i].y - ppoint.y);
lenq := length (coordsq[i].x - qpoint.x,
coordsq[i].y - qpoint.y);
if (lenp < minimumspacing) and (lenq < minimumspacing) then
insertionpoint := 0 { The point is already listed. }
else if tp < tparamsp[i] then
begin
insertionpoint := i; { Insert here instead of at end. }
i := numintersections
end
end;
if insertionpoint <> numintersections + 1 then
for i := numintersections + 1 downto insertionpoint + 1 do
begin
tparamsp[i] := tparamsp[i - 1];
coordsp[i] := coordsp[i - 1];
tparamsq[i] := tparamsq[i - 1];
coordsq[i] := coordsq[i - 1]
end;
tparamsp[insertionpoint] := tp;
coordsp[insertionpoint] := ppoint;
tparamsq[insertionpoint] := tq;
coordsq[insertionpoint] := qpoint;
numintersections := numintersections + 1
end
end;
{ Find intersections between portions of two curves. }
procedure findportionintersections (pportion, qportion : portion);
var
tp, tq : real;
pport1, pport2 : portion;
qport1, qport2 : portion;
begin
if rectanglesoverlap (pportion.endpt0, pportion.endpt1,
qportion.endpt0, qportion.endpt1) then
begin
if flatenough (pportion, flatnesstolerance) then
begin
if flatenough (qportion, flatnesstolerance) then
begin
segmentparameters (pportion.endpt0, pportion.endpt1,
qportion.endpt0, qportion.endpt1,
tp, tq);
if 0.0 <= tp then
begin
tp := (1.0 - tp) * pportion.t0 + tp * pportion.t1;
tq := (1.0 - tq) * qportion.t0 + tq * qportion.t1;
insertintersection (pportion.curv, tp,
qportion.curv, tq)
end
end
else
begin
bisectportion (qportion, qport1, qport2);
findportionintersections (pportion, qport1);
findportionintersections (pportion, qport2)
end
end
else
begin
bisectportion (pportion, pport1, pport2);
if flatenough (qportion, flatnesstolerance) then
begin
findportionintersections (pport1, qportion);
findportionintersections (pport2, qportion)
end
else
begin
bisectportion (qportion, qport1, qport2);
findportionintersections (pport1, qport1);
findportionintersections (pport1, qport2);
findportionintersections (pport2, qport1);
findportionintersections (pport2, qport2)
end
end
end
end;
{ Find intersections in [0,1]. }
procedure findintersections (p, q : curve);
var
tpx, tpy, tqx, tqy : real;
tp, tq : array [1 .. 4] of real;
ppoints, qpoints : array [1 .. 4] of point;
np, nq, i, j : integer;
pportion, qportion : portion;
procedure pfindcriticalpts;
var i : integer;
begin
tp[1] := 0.0;
tp[2] := 1.0;
np := 2;
tpx := spowercriticalpt (p.x);
tpy := spowercriticalpt (p.y);
if (0.0 < tpx) and (tpx < 1.0) then
begin
np := np + 1;
tp[np] := tpx
end;
if (0.0 < tpy) and (tpy < 1.0) and (tpy <> tpx) then
begin
np := np + 1;
tp[np] := tpy
end;
realsort (np, tp);
for i := 1 to np do
ppoints[i] := curveeval (p, tp[i])
end;
procedure qfindcriticalpts;
var i : integer;
begin
tq[1] := 0.0;
tq[2] := 1.0;
nq := 2;
tqx := spowercriticalpt (q.x);
tqy := spowercriticalpt (q.y);
if (0.0 < tqx) and (tqx < 1.0) then
begin
nq := nq + 1;
tq[nq] := tqx
end;
if (0.0 < tqy) and (tqy < 1.0) and (tqy <> tqx) then
begin
nq := nq + 1;
tq[nq] := tqy
end;
realsort (nq, tq);
for i := 1 to nq do
qpoints[i] := curveeval (q, tq[i])
end;
begin
{ Break the curves at critical points, so one can assume the portion
between two endpoints is monotonic along both axes. }
pfindcriticalpts;
qfindcriticalpts;
{ Find intersections in the cartesian product of portions of the two
curves. (If you would like to compare with the Icon code: In the
Icon, goal-directed evaluation is inserting such cartesian
products into the "workload" set. However, to do this requires
only one "every" construct instead of two, and there is no need
for loop/counter variables.) }
for i := 1 to np - 1 do
for j := 1 to nq - 1 do
begin
pportion := makeportion (p, tp[i], tp[i + 1],
ppoints[i], ppoints[i + 1]);
qportion := makeportion (q, tq[j], tq[j + 1],
qpoints[j], qpoints[j + 1]);
findportionintersections (pportion, qportion);
end
end;
begin
pglobal := controlstospower (makepoint (-1.0, 0.0),
makepoint ( 0.0, 10.0),
makepoint ( 1.0, 0.0));
qglobal := controlstospower (makepoint ( 2.0, 1.0),
makepoint (-8.0, 2.0),
makepoint ( 2.0, 3.0));
numintersections := 0;
findintersections (pglobal, qglobal);
writeln;
writeln (' convex up ',
' convex left');
for i := 1 to numintersections do
writeln (' ',
tparamsp[i]:11:8, ' (',
coordsp[i].x:11:8, ', ',
coordsp[i].y:11:8, ') ',
tparamsq[i]:11:8, ' (',
coordsq[i].x:11:8, ', ',
coordsq[i].y:11:8, ')');
writeln
end.
- Output:
convex up convex left 0.07250828 (-0.85498344, 1.34501661) 0.17250830 (-0.85498373, 1.34501660) 0.15948753 (-0.68102494, 2.68102517) 0.84051247 (-0.68102517, 2.68102494) 0.82749170 ( 0.65498340, 2.85498373) 0.92749172 ( 0.65498339, 2.85498344) 0.94051247 ( 0.88102493, 1.11897533) 0.05948753 ( 0.88102467, 1.11897507)
Phix
Aside: at long last found my first ever real-world use of sq_atom()... and o/c it had a silly bug!
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 {c1,c2,c3} = q, s = 1 - t, u1 = (s * c1) + (t * c2), v1 = (s * c2) + (t * c3), m = (s * u1) + (t * v1) return {{c1,u1,m},{m,v1,c3}} 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 xa1, ya1, xa2, ya2, xb1, yb1, xb2, yb2) assert(xa1<=xa2 and ya1<=ya2 and xb1<=xb2 and yb1<=yb2) return xb1 <= xa2 and xa1 <= xb2 and yb1 <= ya2 and ya1 <= yb2 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), ymin = max(pymin,qymin), ymax = min(pymax,qymax) assert(xmax >= xmin and ymax >= ymin) if xmax-xmin <= tolerance and ymax-ymin <= tolerance then -- we found a suitable intersection! return {(xmin+xmax)/2,(ymin+ymax)/2} end if return true -- accept/further subdivide end if return false -- exclude end function function seems_to_be_a_duplicate(sequence intersections, xy, atom spacing) for pt in intersections do if abs(pt[X]-xy[X])<spacing and abs(pt[Y]-xy[Y])<spacing then return true end if end for return false end function function find_intersections(quadratic_curve p, q, atom tolerance) sequence insects = {}, todo = {{p,q}} while length(todo) do {{p,q},todo} = {todo[1],todo[2..$]} object insect = test_intersection(p, q, tolerance) if sequence(insect) then if not seems_to_be_a_duplicate(insects,insect,tolerance*10) then insects &= {insect} end if elsif insect then sequence {p1,p2} = subdivide_quadratic_curve(p,0.5), {q1,q2} = subdivide_quadratic_curve(q,0.5) todo &= {{p1,q1},{p1,q2},{p2,q1},{p2,q2}} end if end while insects = sort_columns(insects,{-Y,X}) return insects 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) printf(1,"Intersections from top to bottom:\n") pp(intersections,{pp_Nest,1,pp_FltFmt,"%9.6f"})
- Output:
Intersections from top to bottom:
{{ 0.654983, 2.854984},
{-0.681025, 2.681025},
{-0.854984, 1.345017},
{ 0.881025, 1.118975}}
A Phix implementation of the "rectangle-pruned curve-flattening algorithm"
This is the recursive version of the algorithm that appeared first (in non-recursive form) in the Icon example.
-- -*- mode: indented-text; tab-width: 2; -*-
enum X, Y
enum C0, C1, C2
enum CURV, T0, T1, ENDPT0, ENDPT1
function normlength (atom a_x, a_y)
-- A "length" according to our chosen norm (which happens to be the
-- max norm). Here (a_x,a_y) is a vector used as measuring rod.
return max (abs (a_x), abs (a_y))
end function
function compare_normlengths (atom a_x, a_y, b_x, b_y)
-- Here is a general implementation for comparison of two
-- "normlength". For the euclidean norm, you might wish to skip
-- taking square roots.
atom len_a = normlength (a_x, a_y),
len_b = normlength (b_x, b_y),
cmpval = 0
if len_a < len_b then
cmpval = -1
elsif len_a > len_b then
cmpval = 1
end if
return cmpval
end function
function controls_to_spower (sequence controls)
-- Convert from control points (that is, Bernstein basis) to the
-- symmetric power basis.
sequence pt0 = controls[1],
pt1 = controls[2],
pt2 = controls[3]
return {{pt0[X], (2 * pt1[X]) - pt0[X] - pt2[X], pt2[X]},
{pt0[Y], (2 * pt1[Y]) - pt0[Y] - pt2[Y], pt2[Y]}}
end function
function spower_eval (sequence spow, atom t)
-- Evaluate an s-power spline at t.
return (spow[C0] + (spow[C1] * t)) * (1 - t) + (spow[C2] * t)
end function
function curve_eval (sequence curv, atom t)
-- Evaluate a curve at t.
return {spower_eval (curv[X], t),
spower_eval (curv[Y], t)}
end function
function spower_center_coef (sequence spow, atom t0, t1)
-- Return the center coefficient for the [t0,t1] portion of an
-- s-power spline. (The endpoint coefficients can be found with
-- spower_eval.) }
return spow[C1] * ((t1 - t0 - t0) * t1 + (t0 * t0))
end function
function spower_critical_pt (sequence spow)
-- Return t in (0,1) where p is at a critical point, else return -1.
atom crit_pt = -1
if spow[C1] != 0 then -- If c1 is zero, then the spline is linear.
if spow[C1] = spow[C2] then
crit_pt = 0.5 -- The spline is "pulse-like".
else
-- The critical point is the root of the derivative.
atom t = (0.5 * (spow[C2] + spow[C1] - spow[C0])) / spow[C1]
if (0 < t) and (t < 1) then
crit_pt = t
end if
end if
end if
return crit_pt
end function
function bisect_portion (sequence port)
-- Bisect a portion and pre-compute the new shared endpoint.
atom t_mid = 0.5 * (port[T0] + port[T1])
sequence pt_mid = curve_eval (port[CURV], t_mid)
return {{port[CURV], port[T0], t_mid, port[ENDPT0], pt_mid},
{port[CURV], t_mid, port[T1], pt_mid, port[ENDPT1]}}
end function
function rectangles_overlap (sequence a0, a1, b0, b1)
-- Do the rectangles with corners at (a0,a1) and (b0,b1) overlap at
-- all?
return ((min (a0[X], a1[X]) <= max (b0[X], b1[X]))
and (min (b0[X], b1[X]) <= max (a0[X], a1[X]))
and (min (a0[Y], a1[Y]) <= max (b0[Y], b1[Y]))
and (min (b0[Y], b1[Y]) <= max (a0[Y], a1[Y])))
end function
function segment_parameters (sequence a0, a1, b0, b1)
-- Return the respective [0,1] parameters of line segments (a0,a1)
-- and (b0,b1), for their intersection point. If there are not two
-- such parameters, return -1 for both values.
atom axdiff = a1[X] - a0[X],
aydiff = a1[Y] - a0[Y],
bxdiff = b1[X] - b0[X],
bydiff = b1[Y] - b0[Y],
denom = (axdiff * bydiff) - (aydiff * bxdiff),
anumer = ((bxdiff * a0[Y]) - (bydiff * a0[X])
+ (b0[X] * b1[Y]) - (b1[X] * b0[Y])),
ta = anumer / denom,
tb = -1
if (ta < 0.0) or (1.0 < ta) then
ta = -1
else
atom bnumer = -((axdiff * b0[Y]) - (aydiff * b0[X])
+ (a0[X] * a1[Y]) - (a1[X] * a0[Y]))
tb = bnumer / denom
if (tb < 0.0) or (1.0 < tb) then
ta = -1
tb = -1
end if
end if
return {ta, tb}
end function
function flat_enough (sequence port, atom tol)
-- Is a curve portion flat enough to be treated as a line segment
-- between its endpoints?
-- The degree-2 s-power polynomials are 1-t, t(1-t), t. We want to
-- remove the terms in t(1-t). The maximum of t(1-t) is 1/4, reached
-- at t=1/2. That accounts for the 1/4=0.25 in the following.
atom xcentercoef = spower_center_coef (port[CURV][X], port[T0],
port[T1]),
ycentercoef = spower_center_coef (port[CURV][Y], port[T0],
port[T1]),
xlen = port[ENDPT1][X] - port[ENDPT0][X],
ylen = port[ENDPT1][Y] - port[ENDPT0][Y]
return (compare_normlengths (0.25 * xcentercoef,
0.25 * ycentercoef,
tol * xlen, tol * ylen) <= 0)
end function
function find_portion_intersections (sequence pportion, qportion,
atom tol)
-- Find intersections between portions of two curves. Return pairs
-- of t-parameters. There may be duplicates and (due to truncation
-- error) near-intersections.
sequence intersections = {}
if rectangles_overlap (pportion[ENDPT0], pportion[ENDPT1],
qportion[ENDPT0], qportion[ENDPT1]) then
if flat_enough (pportion, tol) then
if flat_enough (qportion, tol) then
atom tp, tq
{tp, tq} = segment_parameters (pportion[ENDPT0],
pportion[ENDPT1],
qportion[ENDPT0],
qportion[ENDPT1])
if 0 <= tp then
tp = (1 - tp) * pportion[T0] + tp * pportion[T1]
tq = (1 - tq) * qportion[T0] + tq * qportion[T1]
intersections &= {{tp, tq}}
end if
else
sequence qport1, qport2
{qport1, qport2} = bisect_portion (qportion)
intersections &=
find_portion_intersections (pportion, qport1, tol)
intersections &=
find_portion_intersections (pportion, qport2, tol)
end if
else
if flat_enough (qportion, tol) then
sequence pport1, pport2
{pport1, pport2} = bisect_portion (pportion)
intersections &=
find_portion_intersections (pport1, qportion, tol)
intersections &=
find_portion_intersections (pport2, qportion, tol)
else
sequence pport1, pport2
sequence qport1, qport2
{pport1, pport2} = bisect_portion (pportion)
{qport1, qport2} = bisect_portion (qportion)
intersections &=
find_portion_intersections (pport1, qport1, tol)
intersections &=
find_portion_intersections (pport1, qport2, tol)
intersections &=
find_portion_intersections (pport2, qport1, tol)
intersections &=
find_portion_intersections (pport2, qport2, tol)
end if
end if
end if
return intersections
end function
function find_intersections (sequence p, q, atom tol)
-- Find intersections in [0,1]. Return pairs of t-parameters.
-- There may be duplicates and (due to truncation error)
-- near-intersections.
-- Break the curves at critical points, so one can assume the
-- portion between two endpoints is monotonic along both axes.
sequence tp = {0, 1}
atom tpx = spower_critical_pt (p[X]),
tpy = spower_critical_pt (p[Y])
if 0 < tpx then
tp &= {tpx}
end if
if 0 < tpy and tpy != tpx then
tp &= {tpy}
end if
tp = sort (tp)
sequence tq = {0, 1}
sequence pvals = {}
for t in tp do
pvals &= {curve_eval (p, t)}
end for
atom tqx = spower_critical_pt (q[X]),
tqy = spower_critical_pt (q[Y])
if 0 < tqx then
tq &= {tqx}
end if
if 0 < tqy and tqy != tqx then
tq &= {tqy}
end if
tq = sort (tq)
sequence qvals = {}
for t in tq do
qvals &= {curve_eval (q, t)}
end for
-- Find intersections in the cartesian product of the monotonic
-- portions of the two curves.
sequence intersections = {}
for i = 1 to length (tp) - 1 do
for j = 1 to length (tq) - 1 do
sequence pportion = {p, tp[i], tp[i + 1],
pvals[i], pvals[i + 1]},
qportion = {q, tq[j], tq[j + 1],
qvals[j], qvals[j + 1]}
intersections &=
find_portion_intersections (pportion, qportion, tol)
end for
end for
return intersections
end function
sequence pcontrols = {{-1, 0}, {0, 10}, {1, 0}},
qcontrols = {{2, 1}, {-8, 2}, {2, 3}},
p = controls_to_spower (pcontrols),
q = controls_to_spower (qcontrols)
atom flatness_tolerance = 0.0001
--
-- With this algorithm and its extension to higher degrees:
--
-- I suspect (albeit VERY, VERY, VERY weakly) merely removing exact
-- duplicates from the returned pairs of t-parameters will suffice to
-- eliminate repeated detections, because (aside from intersections
-- with multiplicities) these SHOULD occur only at endpoints, which
-- adjacent portions share, and where the t-parameters are exact zeros
-- and ones.
--
-- In any case, comparing t-parameters is certainly an alternative to
-- comparing point distances, especially if you want to let
-- intersections have multiplicity (as can easily happen with
-- cubics). Scheme's SRFI-1 has "delete-duplicates", which lets one
-- specify an equivalence predicate other than the default "equal?"--
-- the predicate can be defined as a closure to test closeness to
-- within some tolerance. That is how the code below SHOULD be
-- written.
--
-- But I do not know how to do the same thing so simply in Phix, and
-- thus will merely say "unique" here and let others update the code
-- if they wish. :)
--
sequence t_pairs =
unique (find_intersections (p, q, flatness_tolerance))
printf (1, "\n")
printf (1, " convex up ")
printf (1, " convex left\n")
for tt in t_pairs do
atom tp, tq
{tp, tq} = tt
sequence ppoint = curve_eval (p, tp),
qpoint = curve_eval (q, tq)
printf
(1, " %11.8f (%11.8f, %11.8f) %11.8f (%11.8f, %11.8f)\n",
{tp, ppoint[X], ppoint[Y], tq, qpoint[X], qpoint[Y]})
end for
printf (1, "\n")
- Output:
convex up convex left
0.07250828 (-0.85498344, 1.34501661) 0.17250830 (-0.85498373, 1.34501660)
0.15948753 (-0.68102494, 2.68102517) 0.84051247 (-0.68102517, 2.68102494)
0.82749170 ( 0.65498340, 2.85498373) 0.92749172 ( 0.65498339, 2.85498344)
0.94051247 ( 0.88102493, 1.11897533) 0.05948753 ( 0.88102467, 1.11897507)
Python
This is mostly based on the Icon, but with aspects of the Pascal. In particular, using a set instead of recursion is like the Icon, but having subprograms for computing and comparing lengths is borrowed from the Pascal. (Perhaps that enhancement will be retrofitted into the Icon at some time.)
See also:
#!/bin/env python3
#
# * * *
#
# This is the algorithm that was introduced with the Icon example, and
# perhaps is new (at least in its details). It works by putting both
# curves into the symmetric power basis, then first breaking them at
# their critical points, then doing an adaptive flattening process
# until the problem is reduced to the intersection of two
# lines. Possible lines of inquiry are pruned by looking for overlap
# of the rectangles formed by the endpoints of curve portions.
#
# Unlike Icon, Python does not have goal-directed evaluation
# (GDE). What Python DOES have are "iterables" and
# "comprehensions". Where you see "yield" and comprehensions in the
# Python you will likely see "suspend" and "every" in the Icon.
#
# To put it another way: In Python, there are objects that "can be
# iterated over". In Icon, there are objects that "can produce values
# more than once". In either case, the objects are equivalent to a set
# (albeit an ordered set), and really what THIS algorithm deals with
# is (unordered) sets.
#
# Another thing about Icon to be aware of, when comparing this
# algorithm's implementations, is that Icon does not have boolean
# expressions. It has "succeed" and "fail". An Icon expression either
# "succeeds" and has a value or it "fails" and has no value. An "if"
# construct tests whether an expression succeeded, not what the
# expression's value is. (Booleans are easily "faked", of course, if
# you want them. The language variant Object Icon explicitly
# introduces &yes and &no as boolean values.)
#
# * * *
#
# References on the symmetric power basis:
#
# 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.
#
def length(ax, ay):
'''Length according to some norm, where (ax,ay) is a "measuring
stick" vector. Here I use the max norm.'''
assert isinstance(ax, float)
assert isinstance(ay, float)
return max(abs(ax), abs(ay))
def compare_lengths(ax, ay, bx, by):
'''Having a compare_lengths function lets one compare lengths in
the euclidean metric by comparing the squares of the lengths, and
thus avoiding the square root. The following, however, is a
general implementation.'''
assert isinstance(ax, float)
assert isinstance(ay, float)
assert isinstance(bx, float)
assert isinstance(by, float)
len_a = length(ax, ay)
len_b = length(bx, by)
if len_a < len_b:
cmpval = -1
elif len_a > len_b:
cmpval = 1
else:
cmpval = 0
return cmpval
def rectangles_overlap(a0, a1, b0, b1):
'''Do the rectangles with corners at (a0,a1) and (b0,b1) overlap
at all?'''
assert isinstance(a0, Point)
assert isinstance(a1, Point)
assert isinstance(b0, Point)
assert isinstance(b1, Point)
return ((min(a0.x, a1.x) <= max(b0.x, b1.x))
and (min(b0.x, b1.x) <= max(a0.x, a1.x))
and (min(a0.y, a1.y) <= max(b0.y, b1.y))
and (min(b0.y, b1.y) <= max(a0.y, a1.y)))
def segment_parameters(a0, a1, b0, b1):
'''Do the line segments (a0,a1) and (b0,b1) intersect? If so,
return a tuple of their t-parameter values for the point of
intersection, treating them as parametric splines of degree
1. Otherwise return None.'''
assert isinstance(a0, Point)
assert isinstance(a1, Point)
assert isinstance(b0, Point)
assert isinstance(b1, Point)
retval = None
axdiff = a1.x - a0.x
aydiff = a1.y - a0.y
bxdiff = b1.x - b0.x
bydiff = b1.y - b0.y
denom = (axdiff * bydiff) - (aydiff * bxdiff)
anumer = ((bxdiff * a0.y) - (bydiff * a0.x)
+ (b0.x * b1.y) - (b1.x * b0.y))
ta = anumer / denom
if 0.0 <= ta and ta <= 1.0:
bnumer = -((axdiff * b0.y) - (aydiff * b0.x)
+ (a0.x * a1.y) - (a1.x * a0.y))
tb = bnumer / denom
if 0.0 <= tb and tb <= 1.0:
retval = (ta, tb)
return retval
class Point:
def __init__(self, x, y):
assert isinstance(x, float)
assert isinstance(y, float)
self.x = x
self.y = y
class SPower:
'''Non-parametric spline in s-power basis.'''
def __init__(self, c0, c1, c2):
assert isinstance(c0, float)
assert isinstance(c1, float)
assert isinstance(c2, float)
self.c0 = c0
self.c1 = c1
self.c2 = c2
def val(self, t):
'''Evaluate at t.'''
assert isinstance(t, float)
return (self.c0 + (self.c1 * t)) * (1.0 - t) + (self.c2 * t)
def center_coef(self, t0, t1):
'''Return the center coefficient for the [t0,t1] portion. (The
other coefficients can be found with the val method.)'''
assert isinstance(t0, float)
assert isinstance(t1, float)
return self.c1 * ((t1 - t0 - t0) * t1 + (t0 * t0))
def critical_points(self):
'''Return a set of independent variable values for the
critical points that lie in (0,1).'''
critpoints = set()
if self.c1 != 0: # If c1 is zero then the spline is linear.
if self.c0 == self.c2:
critpoints = {0.5} # The spline is "pulse-like".
else:
# The root of the derivative is the critical point.
t = (0.5 * (self.c2 + self.c1 - self.c0)) / self.c1
if 0.0 < t and t < 1.0:
critpoints = {t}
return critpoints
class Curve:
'''Parametric spline in s-power basis.'''
def __init__(self, x, y):
assert isinstance(x, SPower)
assert isinstance(y, SPower)
self.x = x
self.y = y
@staticmethod
def from_controls(ctl0, ctl1, ctl2):
assert isinstance(ctl0, Point)
assert isinstance(ctl1, Point)
assert isinstance(ctl2, Point)
c0x = ctl0.x
c0y = ctl0.y
c1x = (2.0 * ctl1.x) - ctl0.x - ctl2.x
c1y = (2.0 * ctl1.y) - ctl0.y - ctl2.y
c2x = ctl2.x
c2y = ctl2.y
return Curve(SPower(c0x, c1x, c2x),
SPower(c0y, c1y, c2y))
def val(self, t):
'''Evaluate at t.'''
assert isinstance(t, float)
return Point(self.x.val(t), self.y.val(t))
def critical_points(self):
'''Return a set of t-parameter values for the critical points
that lie in (0,1).'''
return (self.x.critical_points() | self.y.critical_points())
class Portion:
'''Portion of a parametric spline in [t0,t1].'''
default_num_pieces = 2
def __init__(self, curve, t0, t1, endpt0, endpt1):
assert isinstance(curve, Curve)
assert isinstance(t0, float)
assert isinstance(t1, float)
assert isinstance(endpt0, Point)
assert isinstance(endpt1, Point)
self.curve = curve
self.t0 = t0
self.t1 = t1
self.endpt0 = endpt0
self.endpt1 = endpt1
def flat_enough(self, tol):
'''Is the Portion close enough to linear to be treated as a
line segment?'''
# The degree-2 s-power polynomials are 1-t, t(1-t), t. We want
# to remove the terms in t(1-t). The maximum of t(1-t) is 1/4,
# reached at t=1/2. That accounts for the 1/4=0.25 in the
# following.
xcentercoef = self.curve.x.center_coef(self.t0, self.t1)
ycentercoef = self.curve.y.center_coef(self.t0, self.t1)
xlen = self.endpt1.x - self.endpt0.x
ylen = self.endpt1.y - self.endpt0.y
return compare_lengths(0.25 * xcentercoef,
0.25 * ycentercoef,
tol * xlen, tol * ylen) <= 0
def split(self, num_pieces = default_num_pieces):
'''Generate num_pieces sections of the Portion.'''
assert isinstance(num_pieces, int)
assert 2 <= num_pieces
k = 1.0 / num_pieces
ts = [(1.0 - (k * i)) * self.t0 + (k * i) * self.t1
for i in range(1, num_pieces)]
vals = [self.curve.val(t) for t in ts]
ts = [self.t0] + ts + [self.t1]
vals = [self.endpt0] + vals + [self.endpt1]
for i in range(len(ts) - 1):
yield Portion(self.curve, ts[i], ts[i + 1],
vals[i], vals[i + 1])
def find_intersections(p, q, tol):
'''Generate t-parameter pairs detected as corresponding to
intersection points of p and q. There may be duplicate
detections. It is assumed those will be weeded out by later
processing. The tol parameter specifies the "flatness tolerance"
and is a relative tolerance.'''
assert isinstance(p, Curve)
assert isinstance(q, Curve)
# The initial workload is the cartesian product of the monotonic
# portions of p and q, respectively.
tp = [0.0] + sorted(p.critical_points()) + [1.0]
tq = [0.0] + sorted(q.critical_points()) + [1.0]
workload = {(Portion(p, tp[i], tp[i + 1],
p.val(tp[i]), p.val(tp[i + 1])),
Portion(q, tq[j], tq[j + 1],
q.val(tq[j]), q.val(tq[j + 1])))
for i in range(len(tp) - 1)
for j in range(len(tq) - 1)}
while len(workload) != 0:
(pportion, qportion) = workload.pop()
if rectangles_overlap(pportion.endpt0, pportion.endpt1,
qportion.endpt0, qportion.endpt1):
if pportion.flat_enough(tol):
if qportion.flat_enough(tol):
params = segment_parameters(pportion.endpt0,
pportion.endpt1,
qportion.endpt0,
qportion.endpt1)
if params is not None:
(tp, tq) = params
tp = (1 - tp) * pportion.t0 + tp * pportion.t1
tq = (1 - tq) * qportion.t0 + tq * qportion.t1
yield (tp, tq)
else:
workload |= {(pportion, qport)
for qport in qportion.split()}
else:
if qportion.flat_enough(tol):
workload |= {(pport, qportion)
for pport in pportion.split()}
else:
workload |= {(pport, qport)
for pport in pportion.split()
for qport in qportion.split()}
if __name__ == '__main__':
flatness_tolerance = 0.0001
minimum_spacing = 0.000001
p = Curve.from_controls(Point(-1.0, 0.0),
Point(0.0, 10.0),
Point(1.0, 0.0))
q = Curve.from_controls(Point(2.0, 1.0),
Point(-8.0, 2.0),
Point( 2.0, 3.0))
intersections = dict()
for (tp, tq) in find_intersections(p, q, flatness_tolerance):
pval = p.val(tp)
qval = q.val(tq)
if all([(minimum_spacing <=
length(pval.x - intersections[t][1].x,
pval.y - intersections[t][1].y))
and (minimum_spacing <=
length(qval.x - intersections[t][3].x,
qval.y - intersections[t][3].y))
for t in intersections]):
intersections[tp] = (tp, pval, tq, qval)
print()
print(' convex up ',
' convex left');
for k in sorted(intersections):
(tp, pointp, tq, pointq) = intersections[k]
print((" %11.8f (%11.8f, %11.8f) " +
"%11.8f (%11.8f, %11.8f)")
%(tp, pointp.x, pointp.y, tq, pointq.x, pointq.y))
print()
- Output:
convex up convex left 0.07250828 (-0.85498344, 1.34501661) 0.17250830 (-0.85498373, 1.34501660) 0.15948753 (-0.68102494, 2.68102517) 0.84051247 (-0.68102517, 2.68102494) 0.82749170 ( 0.65498340, 2.85498373) 0.92749172 ( 0.65498339, 2.85498344) 0.94051247 ( 0.88102493, 1.11897533) 0.05948753 ( 0.88102467, 1.11897507)
Scheme
The algorithm of the Icon, Python, Pascal, etc.
For R7RS-small, plus SRFI-132 and SRFI-144, both adopted as part of R7RS-large. SRFI-144 is used only to get the value of fl-epsilon.
This is the algorithm of the Icon. The implementation here has support for three different norms and for both relative and absolute tolerance.
;;
;; For R7RS-small plus very limited R7RS-large. The implementation
;; must support libraries, Unicode identifiers, IEEE arithmetic, etc.
;;
;; These will work:
;;
;; gosh -r7 file.scm
;; csc -X r7rs -R r7rs file.scm
;;
(define-library (spower quadratic)
(export spower spower?
spower-c0 spower-c1 spower-c2
plane-curve plane-curve?
plane-curve-x plane-curve-y
control-points->plane-curve
spower-eval plane-curve-eval
center-coef critical-points)
(import (scheme base)
(scheme case-lambda)
(srfi 132) ;; R7RS-large name: (scheme sort)
)
(begin
(define-record-type <spower>
(spower c0 c1 c2)
spower?
(c0 spower-c0)
(c1 spower-c1)
(c2 spower-c2))
(define-record-type <plane-curve>
(plane-curve x y)
plane-curve?
(x plane-curve-x)
(y plane-curve-y))
(define (point->values point)
;; A bit of playfulness on my part: a few ways to write an
;; (x,y)-point: (cons x y), (list x y), (vector x y), and their
;; equivalents.
(cond ((and (pair? point) (real? (car point)))
(cond ((real? (cdr point))
(values (car point) (cdr point)))
((and (real? (cadr point)) (null? (cddr point)))
(values (car point) (cadr point)))
(else (error "not a plane point" point))))
((and (vector? point) (= (vector-length point) 2))
(values (vector-ref point 0) (vector-ref point 1)))
(else (error "not a plane point" point))))
(define control-points->plane-curve
;; A bit more playfulness: the control points can be given
;; separately, as a list, or as a vector.
(case-lambda
((pt0 pt1 pt2)
(let-values
(((cx0 cy0) (point->values pt0))
((cx1 cy1) (point->values pt1))
((cx2 cy2) (point->values pt2)))
(plane-curve (spower cx0 (- (+ cx1 cx1) cx0 cx2) cx2)
(spower cy0 (- (+ cy1 cy1) cy0 cy2) cy2))))
((pt-list-or-vector)
(apply control-points->plane-curve
(if (vector? pt-list-or-vector)
(vector->list pt-list-or-vector)
pt-list-or-vector)))))
(define (spower-eval poly t)
;; Evaluate the polynomial at t.
(let ((c0 (spower-c0 poly))
(c1 (spower-c1 poly))
(c2 (spower-c2 poly)))
(+ (* (+ c0 (* c1 t)) (- 1 t)) (* c2 t))))
(define (plane-curve-eval curve t)
;; Evaluate the plane curve at t, returning x and y as multiple
;; values.
(values (spower-eval (plane-curve-x curve) t)
(spower-eval (plane-curve-y curve) t)))
(define (center-coef poly t0 t1)
;; Return the center coefficient for the [t0,t1] portion. (The
;; other coefficients can be found with the spower-eval
;; procedured.)
(let ((c1 (spower-c1 poly)))
(* c1 (+ (* (- t1 t0 t0) t1) (* t0 t0)))))
(define (critical-points poly-or-curve)
;; Return a sorted list of independent variable values for the
;; critical points that lie in (0,1). Duplicates are removed.
(cond ((plane-curve? poly-or-curve)
(let ((X (plane-curve-x poly-or-curve))
(Y (plane-curve-y poly-or-curve)))
(list-delete-neighbor-dups
=
(list-sort < (append (critical-points X)
(critical-points Y))))))
((spower? poly-or-curve)
(let ((c0 (spower-c0 poly-or-curve))
(c1 (spower-c1 poly-or-curve))
(c2 (spower-c2 poly-or-curve)))
(cond ((zero? c1) '()) ; The spline is linear.
((= c0 c2) '(1/2)) ; The spline is "pulse-like".
(else
;; The root of the derivative is the critical
;; point.
(let ((t (/ (- (+ c2 c1) c0) (+ c1 c1))))
(if (and (positive? t) (< t 1))
(list t)
'()))))))
(else (error "not an spower polynomial or plane-curve"
poly-or-curve))))
)) ;; end library (spower quadratic)
(define-library (spower quadratic intersections)
;; Parameters. (That is, variables whose value settings have
;; "dynamic" rather than lexical scope.)
(export *tolerance-norm*
*flatness-tolerance*
*absolute-tolerance*)
(export find-intersection-parameters)
(import (scheme base)
(scheme case-lambda)
(scheme inexact)
(spower quadratic)
(srfi 144) ;; R7RS-large name: (scheme flonum)
)
(begin
(define-record-type <portion>
(make-portion curve t0 t1 endpt0 endpt1)
portion?
(curve curve@)
(t0 t0@)
(t1 t1@)
(endpt0 endpt0@)
(endpt1 endpt1@))
(define (curve-eval curve t)
(call-with-values (lambda () (plane-curve-eval curve t))
cons))
(define (0<=x<=1 x)
(and (<= 0 x) (<= x 1)))
(define (lerp t a b) ; "linear interpolation"
(+ (* (- 1 t) a) (* t b)))
(define (bisect-portion portion)
;; Bisect portion and return the two new portions as multiple
;; values.
(let ((curve (curve@ portion))
(t0 (t0@ portion))
(t1 (t1@ portion))
(pt0 (endpt0@ portion))
(pt1 (endpt1@ portion)))
(let* ((t½ (* 1/2 (+ t0 t1)))
(pt½ (curve-eval curve t½)))
(values (make-portion curve t0 t½ pt0 pt½)
(make-portion curve t½ t1 pt½ pt1)))))
(define (check-norm x)
(cond ((and (positive? x) (infinite? x)) x)
((= x 1) (exact x))
((= x 2) (exact x))
(else (error "not a norm we can handle" x))))
(define (check-flatness-tolerance x)
(cond ((zero? x) x)
((positive? x) x)
(else
(error
"not a flatness (relative) tolerance we can handle"
x))))
(define (check-absolute-tolerance x)
(cond ((zero? x) x)
((positive? x) x)
(else (error "not an absolute tolerance we can handle"
x))))
(define *tolerance-norm*
;; To be fancier, we could take strings such as "taxicab",
;; "euclidean", "max", etc., as arguments to the
;; constructor. But here we are taking only the value of p in a
;; p-norm (= pth root of the sum of |x| raised p and |y| raised
;; p), and then only for p = 1, 2, +inf.0.
(make-parameter +inf.0 check-norm)) ;; Default is the max norm.
;; A relative tolerance. This setting seems to me rather strict
;; for a lot of applications. But you can override it.
(define *flatness-tolerance*
(make-parameter 0.0001 check-flatness-tolerance))
(define *absolute-tolerance*
(make-parameter (* 50 fl-epsilon) check-absolute-tolerance))
(define (compare-lengths norm ax ay bx by)
(define (compare-lengths-1-norm)
(let ((‖a‖ (+ (abs ax) (abs ay)))
(‖b‖ (+ (abs bx) (abs by))))
(cond ((< ‖a‖ ‖b‖) -1)
((> ‖a‖ ‖b‖) 1)
(else 0))))
(define (compare-lengths-2-norm)
(let ((‖a‖² (* ax ay))
(‖b‖² (* bx by)))
(cond ((< ‖a‖² ‖b‖²) -1)
((> ‖a‖² ‖b‖²) 1)
(else 0))))
(define (compare-lengths-inf-norm)
(let ((‖a‖ (max (abs ax) (abs ay)))
(‖b‖ (max (abs bx) (abs by))))
(cond ((< ‖a‖ ‖b‖) -1)
((> ‖a‖ ‖b‖) 1)
(else 0))))
(cond ((= norm 1) (compare-lengths-1-norm))
((= norm 2) (compare-lengths-2-norm))
(else (compare-lengths-inf-norm))))
(define (compare-to-tol norm ax ay tol)
(define (compare-to-tol-1-norm)
(let ((‖a‖ (+ (abs ax) (abs ay))))
(cond ((< ‖a‖ tol) -1)
((> ‖a‖ tol) 1)
(else 0))))
(define (compare-to-tol-2-norm)
(let ((‖a‖² (* ax ay))
(tol² (* tol tol)))
(cond ((< ‖a‖² tol²) -1)
((> ‖a‖² tol²) 1)
(else 0))))
(define (compare-to-tol-inf-norm)
(let ((‖a‖ (max (abs ax) (abs ay))))
(cond ((< ‖a‖ tol) -1)
((> ‖a‖ tol) 1)
(else 0))))
(cond ((= norm 1) (compare-to-tol-1-norm))
((= norm 2) (compare-to-tol-2-norm))
(else (compare-to-tol-inf-norm))))
(define (flat-enough? portion norm rtol atol)
;; Is the portion flat enough or small enough to be treated as
;; if it were a line segment?
;; The degree-2 s-power polynomials are 1-t, t(1-t), t. We
;; want to remove the terms in t(1-t). The maximum of t(1-t)
;; is 1/4, reached at t=1/2. That accounts for the 1/4 in the
;; following.
(let ((curve (curve@ portion))
(t0 (t0@ portion))
(t1 (t1@ portion))
(pt0 (endpt0@ portion))
(pt1 (endpt1@ portion)))
(let ((X (plane-curve-x curve))
(Y (plane-curve-y curve)))
(let ((errx (* 1/4 (center-coef X t0 t1)))
(erry (* 1/4 (center-coef Y t0 t1)))
(testx (* rtol (- (car pt1) (car pt0))))
(testy (* rtol (- (cdr pt1) (cdr pt0)))))
(or (<= (compare-lengths norm errx erry testx testy) 0)
(<= (compare-to-tol norm errx erry atol) 0))))))
(define (rectangles-overlap a0 a1 b0 b1)
;;
;; Do the rectangles with corners at (a0,a1) and (b0,b1) have
;; overlapping interiors or edges? The corner points must be
;; represented as cons-pairs, with x as the car and y as the
;; cdr.
;;
;; (A side note: I had thought for a few hours that checking
;; only for overlapping interiors offered some advantages, but I
;; changed my mind.)
;;
(let ((a0x (min (car a0) (car a1)))
(a1x (max (car a0) (car a1)))
(b0x (min (car b0) (car b1)))
(b1x (max (car b0) (car b1))))
(cond ((< b1x a0x) #f)
((< a1x b0x) #f)
(else
(let ((a0y (min (cdr a0) (cdr a1)))
(a1y (max (cdr a0) (cdr a1)))
(b0y (min (cdr b0) (cdr b1)))
(b1y (max (cdr b0) (cdr b1))))
(cond ((< b1y a0y) #f)
((< a1y b0y) #f)
(else #t)))))))
(define (segment-parameters a0 a1 b0 b1)
;; Return (as multiple values) the respective [0,1] parameters
;; of line segments (a0,a1) and (b0,b1), for their intersection
;; point. Return #f and #f if there is no intersection. The
;; endpoints must be represented as cons-pairs, with x as the
;; car and y as the cdr.
(let ((a0x (car a0)) (a0y (cdr a0))
(a1x (car a1)) (a1y (cdr a1))
(b0x (car b0)) (b0y (cdr b0))
(b1x (car b1)) (b1y (cdr b1)))
(let ((axdiff (- a1x a0x))
(aydiff (- a1y a0y))
(bxdiff (- b1x b0x))
(bydiff (- b1y b0y)))
(let* ((denom (- (* axdiff bydiff) (* aydiff bxdiff)))
(anumer (+ (* bxdiff a0y)
(- (* bydiff a0x))
(* b0x b1y)
(- (* b1x b0y))))
(ta (/ anumer denom)))
(if (not (0<=x<=1 ta))
(values #f #f)
(let* ((bnumer (- (+ (* axdiff b0y)
(- (* aydiff b0x))
(* a0x a1y)
(- (* a1x a0y)))))
(tb (/ bnumer denom)))
(if (not (0<=x<=1 tb))
(values #f #f)
(values ta tb))))))))
(define (initial-workload P Q)
;; Generate an initial workload from plane curves P and Q. One
;; does this by splitting the curves at their critical points
;; and constructing the Cartesian product of the resulting
;; portions. The workload is a list of cons-pairs of
;; portions. (The list represents a set, but the obvious place,
;; from which to get an arbitrary pair to work on next, is the
;; head of the list. Thus the list effectively is a stack.)
(define (t->point curve) (lambda (t) (curve-eval curve t)))
(let* ((pparams0 `(0 ,@(critical-points P) 1))
(pvalues0 (map (t->point P) pparams0))
(qparams0 `(0 ,@(critical-points Q) 1))
(qvalues0 (map (t->point Q) qparams0)))
(let loop ((pparams pparams0) (pvalues pvalues0)
(qparams qparams0) (qvalues qvalues0)
(workload '()))
(cond ((null? (cdr pparams)) workload)
((null? (cdr qparams))
(loop (cdr pparams) (cdr pvalues)
qparams0 qvalues0 workload))
(else
(let ((pportion
(make-portion P (car pparams) (cadr pparams)
(car pvalues) (cadr pvalues)))
(qportion
(make-portion Q (car qparams) (cadr qparams)
(car qvalues) (cadr qvalues))))
(loop pparams pvalues (cdr qparams) (cdr qvalues)
(cons (cons pportion qportion)
workload))))))))
(define (params≅? a b)
(and (≅? (car a) (car b))
(≅? (cdr a) (cdr b))))
(define (≅? x y)
;; I do not know this test is any better, for THIS algorithm,
;; than an exact equality test. But it is no worse.
(<= (abs (- x y)) (* 0.5 fl-epsilon (max (abs x) (abs y)))))
(define (include-params tp tq lst)
(let ((params (cons tp tq)))
(if (not (member params lst params≅?))
(cons params lst)
lst)))
(define find-intersection-parameters
(case-lambda
((P Q) (find-intersection-parameters P Q #f #f #f))
((P Q norm) (find-intersection-parameters P Q norm #f #f))
((P Q norm rtol) (find-intersection-parameters
P Q norm rtol #f))
((P Q norm rtol atol)
(let ((norm (check-norm
(or norm (*tolerance-norm*))))
(rtol (check-flatness-tolerance
(or rtol (*flatness-tolerance*))))
(atol (check-absolute-tolerance
(or atol (*absolute-tolerance*)))))
(%%find-intersection-parameters P Q norm rtol atol)))))
(define (%%find-intersection-parameters P Q norm rtol atol)
(let loop ((workload (initial-workload P Q))
(params '()))
(if (null? workload)
params
(let ((pportion (caar workload))
(qportion (cdar workload))
(workload (cdr workload)))
(if (not (rectangles-overlap (endpt0@ pportion)
(endpt1@ pportion)
(endpt0@ qportion)
(endpt1@ qportion)))
(loop workload params)
(if (flat-enough? pportion norm rtol atol)
(if (flat-enough? qportion norm rtol atol)
(let-values
(((tp tq)
(segment-parameters
(endpt0@ pportion)
(endpt1@ pportion)
(endpt0@ qportion)
(endpt1@ qportion))))
(if tp
(let* ((tp0 (t0@ pportion))
(tp1 (t1@ pportion))
(tq0 (t0@ qportion))
(tq1 (t1@ qportion))
(tp (lerp tp tp0 tp1))
(tq (lerp tq tq0 tq1)))
(loop workload (include-params
tp tq params)))
(loop workload params)))
(let-values (((qport1 qport2)
(bisect-portion qportion)))
(loop `(,(cons pportion qport1)
,(cons pportion qport2)
. ,workload)
params)))
(if (flat-enough? qportion norm rtol atol)
(let-values (((pport1 pport2)
(bisect-portion pportion)))
(loop `(,(cons pport1 qportion)
,(cons pport2 qportion)
. ,workload)
params))
(let-values (((pport1 pport2)
(bisect-portion pportion))
((qport1 qport2)
(bisect-portion qportion)))
(loop `(,(cons pport1 qport1)
,(cons pport1 qport2)
,(cons pport2 qport1)
,(cons pport2 qport2)
. ,workload)
params)))))))))
)) ;; end library (spower quadratic intersections)
(import (scheme base)
(scheme write)
(spower quadratic)
(spower quadratic intersections)
(srfi 132) ;; R7RS-large name: (scheme sort)
)
(define P (control-points->plane-curve '((-1 0) (0 10) (1 0))))
(define Q (control-points->plane-curve '((2 1) (-8 2) (2 3))))
;; Sort the results by the parameters for P. Thus the intersection
;; points will be sorted left to right.
(define params (list-sort (lambda (tpair1 tpair2)
(< (car tpair1) (car tpair2)))
(find-intersection-parameters P Q)))
(for-each
(lambda (pair)
(let ((tp (car pair))
(tq (cdr pair)))
(let-values (((ax ay) (plane-curve-eval P tp))
((bx by) (plane-curve-eval Q tq)))
(display (inexact tp)) (display "\t(")
(display (inexact ax)) (display ", ")
(display (inexact ay)) (display ") \t")
(display (inexact tq)) (display "\t(")
(display (inexact bx)) (display ", ")
(display (inexact by)) (display ")\n"))))
params)
- Output:
There is no standard for formatted output in R7RS Scheme, so in the following I do not bother. (This does not mean there are not excellent unstandardized ways to do formatted output in Schemes.) This is the output from Gauche Scheme.
0.07250828117222911 (-0.8549834376555417, 1.3450166066735616) 0.17250829973466453 (-0.8549837251463935, 1.345016599469329) 0.15948753092173137 (-0.6810249381565372, 2.6810251680444233) 0.8405124690782686 (-0.6810251680444231, 2.681024938156537) 0.8274917002653355 (0.654983400530671, 2.8549837251463934) 0.9274917188277709 (0.6549833933264384, 2.854983437655542) 0.9405124666929737 (0.8810249333859473, 1.118975333761435) 0.059487533307026316 (0.881024666238565, 1.1189750666140525)
The Scheme implementation from above, plus methods from the ATS
It is possible to reduce the number of bisections by solving a quadratic equation as soon as just one of the curve portions is flat enough to be treated as if it were a line segment. That refinement is implemented here.
;;
;; For R7RS-small plus very limited R7RS-large. The implementation
;; must support libraries, Unicode identifiers, IEEE arithmetic, etc.
;;
;; These will work:
;;
;; gosh -r7 file.scm
;; csc -X r7rs -R r7rs file.scm
;;
(define-library (spower quadratic)
(export spower spower?
spower-c0 spower-c1 spower-c2
plane-curve plane-curve?
plane-curve-x plane-curve-y
control-points->plane-curve
spower-eval plane-curve-eval
center-coef critical-points)
(import (scheme base)
(scheme case-lambda)
(srfi 132) ;; R7RS-large name: (scheme sort)
)
(begin
(define-record-type <spower>
(spower c0 c1 c2)
spower?
(c0 spower-c0)
(c1 spower-c1)
(c2 spower-c2))
(define-record-type <plane-curve>
(plane-curve x y)
plane-curve?
(x plane-curve-x)
(y plane-curve-y))
(define (point->values point)
;; A bit of playfulness on my part: a few ways to write an
;; (x,y)-point: (cons x y), (list x y), (vector x y), and their
;; equivalents.
(cond ((and (pair? point) (real? (car point)))
(cond ((real? (cdr point))
(values (car point) (cdr point)))
((and (real? (cadr point)) (null? (cddr point)))
(values (car point) (cadr point)))
(else (error "not a plane point" point))))
((and (vector? point) (= (vector-length point) 2))
(values (vector-ref point 0) (vector-ref point 1)))
(else (error "not a plane point" point))))
(define control-points->plane-curve
;; A bit more playfulness: the control points can be given
;; separately, as a list, or as a vector.
(case-lambda
((pt0 pt1 pt2)
(let-values
(((cx0 cy0) (point->values pt0))
((cx1 cy1) (point->values pt1))
((cx2 cy2) (point->values pt2)))
(plane-curve (spower cx0 (- (+ cx1 cx1) cx0 cx2) cx2)
(spower cy0 (- (+ cy1 cy1) cy0 cy2) cy2))))
((pt-list-or-vector)
(apply control-points->plane-curve
(if (vector? pt-list-or-vector)
(vector->list pt-list-or-vector)
pt-list-or-vector)))))
(define (spower-eval poly t)
;; Evaluate the polynomial at t.
(let ((c0 (spower-c0 poly))
(c1 (spower-c1 poly))
(c2 (spower-c2 poly)))
(+ (* (+ c0 (* c1 t)) (- 1 t)) (* c2 t))))
(define (plane-curve-eval curve t)
;; Evaluate the plane curve at t, returning x and y as multiple
;; values.
(values (spower-eval (plane-curve-x curve) t)
(spower-eval (plane-curve-y curve) t)))
(define (center-coef poly t0 t1)
;; Return the center coefficient for the [t0,t1] portion. (The
;; other coefficients can be found with the spower-eval
;; procedured.)
(let ((c1 (spower-c1 poly)))
(* c1 (+ (* (- t1 t0 t0) t1) (* t0 t0)))))
(define (critical-points poly-or-curve)
;; Return a sorted list of independent variable values for the
;; critical points that lie in (0,1). Duplicates are removed.
(cond ((plane-curve? poly-or-curve)
(let ((X (plane-curve-x poly-or-curve))
(Y (plane-curve-y poly-or-curve)))
(list-delete-neighbor-dups
=
(list-sort < (append (critical-points X)
(critical-points Y))))))
((spower? poly-or-curve)
(let ((c0 (spower-c0 poly-or-curve))
(c1 (spower-c1 poly-or-curve))
(c2 (spower-c2 poly-or-curve)))
(cond ((zero? c1) '()) ; The spline is linear.
((= c0 c2) '(1/2)) ; The spline is "pulse-like".
(else
;; The root of the derivative is the critical
;; point.
(let ((t (/ (- (+ c2 c1) c0) (+ c1 c1))))
(if (and (positive? t) (< t 1))
(list t)
'()))))))
(else (error "not an spower polynomial or plane-curve"
poly-or-curve))))
)) ;; end library (spower quadratic)
(define-library (spower quadratic intersections)
;; Parameters. (That is, variables whose value settings have
;; "dynamic" rather than lexical scope.)
(export *tolerance-norm*
*flatness-tolerance*
*absolute-tolerance*)
(export find-intersection-parameters)
(import (scheme base)
(scheme case-lambda)
(scheme inexact)
(spower quadratic)
(srfi 144) ;; R7RS-large name: (scheme flonum)
)
;; REMOVE THIS FOR A PRACTICAL APPLICATION.
(import (scheme write))
(begin
(define-record-type <portion>
(make-portion curve t0 t1 endpt0 endpt1)
portion?
(curve curve@)
(t0 t0@)
(t1 t1@)
(endpt0 endpt0@)
(endpt1 endpt1@))
(define (curve-eval curve t)
(call-with-values (lambda () (plane-curve-eval curve t))
cons))
(define (0<=x<=1 x)
(and (<= 0 x) (<= x 1)))
(define (lerp t a b) ; "linear interpolation"
(+ (* (- 1 t) a) (* t b)))
(define (bisect-portion portion)
;; Bisect portion and return the two new portions as multiple
;; values.
(let ((curve (curve@ portion))
(t0 (t0@ portion))
(t1 (t1@ portion))
(pt0 (endpt0@ portion))
(pt1 (endpt1@ portion)))
(let* ((t½ (* 1/2 (+ t0 t1)))
(pt½ (curve-eval curve t½)))
(values (make-portion curve t0 t½ pt0 pt½)
(make-portion curve t½ t1 pt½ pt1)))))
(define (check-norm x)
(cond ((and (positive? x) (infinite? x)) x)
((= x 1) (exact x))
((= x 2) (exact x))
(else (error "not a norm we can handle" x))))
(define (check-flatness-tolerance x)
(cond ((zero? x) x)
((positive? x) x)
(else
(error
"not a flatness (relative) tolerance we can handle"
x))))
(define (check-absolute-tolerance x)
(cond ((zero? x) x)
((positive? x) x)
(else (error "not an absolute tolerance we can handle"
x))))
(define *tolerance-norm*
;; To be fancier, we could take strings such as "taxicab",
;; "euclidean", "max", etc., as arguments to the
;; constructor. But here we are taking only the value of p in a
;; p-norm (= pth root of the sum of |x| raised p and |y| raised
;; p), and then only for p = 1, 2, +inf.0.
(make-parameter +inf.0 check-norm)) ;; Default is the max norm.
;; A relative tolerance. This setting seems to me rather strict
;; for a lot of applications. But you can override it.
(define *flatness-tolerance*
(make-parameter 0.0001 check-flatness-tolerance))
(define *absolute-tolerance*
(make-parameter (* 50 fl-epsilon) check-absolute-tolerance))
(define (compare-lengths norm ax ay bx by)
(define (compare-lengths-1-norm)
(let ((‖a‖ (+ (abs ax) (abs ay)))
(‖b‖ (+ (abs bx) (abs by))))
(cond ((< ‖a‖ ‖b‖) -1)
((> ‖a‖ ‖b‖) 1)
(else 0))))
(define (compare-lengths-2-norm)
(let ((‖a‖² (* ax ay))
(‖b‖² (* bx by)))
(cond ((< ‖a‖² ‖b‖²) -1)
((> ‖a‖² ‖b‖²) 1)
(else 0))))
(define (compare-lengths-inf-norm)
(let ((‖a‖ (max (abs ax) (abs ay)))
(‖b‖ (max (abs bx) (abs by))))
(cond ((< ‖a‖ ‖b‖) -1)
((> ‖a‖ ‖b‖) 1)
(else 0))))
(cond ((= norm 1) (compare-lengths-1-norm))
((= norm 2) (compare-lengths-2-norm))
(else (compare-lengths-inf-norm))))
(define (compare-to-tol norm ax ay tol)
(define (compare-to-tol-1-norm)
(let ((‖a‖ (+ (abs ax) (abs ay))))
(cond ((< ‖a‖ tol) -1)
((> ‖a‖ tol) 1)
(else 0))))
(define (compare-to-tol-2-norm)
(let ((‖a‖² (* ax ay))
(tol² (* tol tol)))
(cond ((< ‖a‖² tol²) -1)
((> ‖a‖² tol²) 1)
(else 0))))
(define (compare-to-tol-inf-norm)
(let ((‖a‖ (max (abs ax) (abs ay))))
(cond ((< ‖a‖ tol) -1)
((> ‖a‖ tol) 1)
(else 0))))
(cond ((= norm 1) (compare-to-tol-1-norm))
((= norm 2) (compare-to-tol-2-norm))
(else (compare-to-tol-inf-norm))))
(define (flat-enough? portion norm rtol atol)
;; Is the portion flat enough or small enough to be treated as
;; if it were a line segment?
;; The degree-2 s-power polynomials are 1-t, t(1-t), t. We
;; want to remove the terms in t(1-t). The maximum of t(1-t)
;; is 1/4, reached at t=1/2. That accounts for the 1/4 in the
;; following.
(let ((curve (curve@ portion))
(t0 (t0@ portion))
(t1 (t1@ portion))
(pt0 (endpt0@ portion))
(pt1 (endpt1@ portion)))
(let ((X (plane-curve-x curve))
(Y (plane-curve-y curve)))
(let ((errx (* 1/4 (center-coef X t0 t1)))
(erry (* 1/4 (center-coef Y t0 t1)))
(testx (* rtol (- (car pt1) (car pt0))))
(testy (* rtol (- (cdr pt1) (cdr pt0)))))
(or (<= (compare-lengths norm errx erry testx testy) 0)
(<= (compare-to-tol norm errx erry atol) 0))))))
(define (rectangles-overlap a0 a1 b0 b1)
;;
;; Do the rectangles with corners at (a0,a1) and (b0,b1) have
;; overlapping interiors or edges? The corner points must be
;; represented as cons-pairs, with x as the car and y as the
;; cdr.
;;
;; (A side note: I had thought for a few hours that checking
;; only for overlapping interiors offered some advantages, but I
;; changed my mind.)
;;
(let ((a0x (min (car a0) (car a1)))
(a1x (max (car a0) (car a1)))
(b0x (min (car b0) (car b1)))
(b1x (max (car b0) (car b1))))
(cond ((< b1x a0x) #f)
((< a1x b0x) #f)
(else
(let ((a0y (min (cdr a0) (cdr a1)))
(a1y (max (cdr a0) (cdr a1)))
(b0y (min (cdr b0) (cdr b1)))
(b1y (max (cdr b0) (cdr b1))))
(cond ((< b1y a0y) #f)
((< a1y b0y) #f)
(else #t)))))))
(define (segment-parameters a0 a1 b0 b1)
;; Return (as multiple values) the respective [0,1] parameters
;; of line segments (a0,a1) and (b0,b1), for their intersection
;; point. Return #f and #f if there is no intersection. The
;; endpoints must be represented as cons-pairs, with x as the
;; car and y as the cdr.
(let ((a0x (car a0)) (a0y (cdr a0))
(a1x (car a1)) (a1y (cdr a1))
(b0x (car b0)) (b0y (cdr b0))
(b1x (car b1)) (b1y (cdr b1)))
(let ((axdiff (- a1x a0x))
(aydiff (- a1y a0y))
(bxdiff (- b1x b0x))
(bydiff (- b1y b0y)))
(let* ((denom (- (* axdiff bydiff) (* aydiff bxdiff)))
(anumer (+ (* bxdiff a0y)
(- (* bydiff a0x))
(* b0x b1y)
(- (* b1x b0y))))
(ta (/ anumer denom)))
(if (not (0<=x<=1 ta))
(values #f #f)
(let* ((bnumer (- (+ (* axdiff b0y)
(- (* aydiff b0x))
(* a0x a1y)
(- (* a1x a0y)))))
(tb (/ bnumer denom)))
(if (not (0<=x<=1 tb))
(values #f #f)
(values ta tb))))))))
(define (initial-workload P Q)
;; Generate an initial workload from plane curves P and Q. One
;; does this by splitting the curves at their critical points
;; and constructing the Cartesian product of the resulting
;; portions. The workload is a list of cons-pairs of
;; portions. (The list represents a set, but the obvious place,
;; from which to get an arbitrary pair to work on next, is the
;; head of the list. Thus the list effectively is a stack.)
(define (t->point curve) (lambda (t) (curve-eval curve t)))
;;
;; There are endless ways to write the loop or recursion to
;; compute the Cartesian product. The original Scheme did it one
;; way. Here is a completely different way. It involves having
;; two pointers go through a list at the same time, spaced one
;; node apart. This is done on more than one list at a
;; time. Also, this time the algorithm is done "procedurally"
;; rather than "functionally".
;;
(let* ((pparams0 `(0 ,@(critical-points P) 1))
(pvalues0 (map (t->point P) pparams0))
(qparams0 `(0 ,@(critical-points Q) 1))
(qvalues0 (map (t->point Q) qparams0))
(workload '()))
(do ((ppi pparams0 (cdr ppi))
(ppj (cdr pparams0) (cdr ppj))
(pvi pvalues0 (cdr pvi))
(pvj (cdr pvalues0) (cdr pvj)))
((null? ppj))
(do ((qpi qparams0 (cdr qpi))
(qpj (cdr qparams0) (cdr qpj))
(qvi qvalues0 (cdr qvi))
(qvj (cdr qvalues0) (cdr qvj)))
((null? qpj))
(let ((pportion (make-portion P (car ppi) (car ppj)
(car pvi) (car pvj)))
(qportion (make-portion Q (car qpi) (car qpj)
(car qvi) (car qvj))))
(set! workload `((,pportion . ,qportion)
. ,workload)))))
workload))
(define (params≅? a b)
(and (≅? (car a) (car b))
(≅? (cdr a) (cdr b))))
(define (≅? x y)
;; PERHAPS IT WOULD BE BETTER TO HAVE THIS DEFINITION BE A
;; PARAMETER.
(<= (abs (- x y)) (* 0.5 fl-epsilon (max (abs x) (abs y)))))
(define (include-params tp tq lst)
(let ((params (cons tp tq)))
(if (not (member params lst params≅?))
(cons params lst)
lst)))
(define find-intersection-parameters
(case-lambda
((P Q) (find-intersection-parameters P Q #f #f #f))
((P Q norm) (find-intersection-parameters P Q norm #f #f))
((P Q norm rtol) (find-intersection-parameters
P Q norm rtol #f))
((P Q norm rtol atol)
(let ((norm (check-norm
(or norm (*tolerance-norm*))))
(rtol (check-flatness-tolerance
(or rtol (*flatness-tolerance*))))
(atol (check-absolute-tolerance
(or atol (*absolute-tolerance*)))))
(%%find-intersection-parameters P Q norm rtol atol)))))
;; REMOVE THIS FOR A PRACTICAL APPLICATION.
(define NUM-ITERATIONS 0)
(define (%%find-intersection-parameters P Q norm rtol atol)
;;
;; Among other things that this version of the program
;; demonstrates: you can safely break up a standard Scheme loop
;; into a mutual tail recursion. There will be no "stack
;; blow-up". (At least if you do not count as "stack blow-up"
;; the very strange way CHICKEN Scheme works.)
;;
;; (It is interesting, by the way, to know in what languages one
;; can do this sort of thing, to what degree. In standard
;; Scheme, it is possible without any restrictions. In ATS, one
;; can do it safely as long as an "fnx" construct is possible,
;; because this is precisely what "fnx" is for. But I have tried
;; a very, very simple mutual recursion in Standard ML, and had
;; it work fine in Poly/ML but blow up the stack in MLton, even
;; though MLton is overall the more aggressively optimizing
;; compiler.)
;;
(define (start-here workload params)
;; REMOVE THIS FOR A PRACTICAL APPLICATION.
(set! NUM-ITERATIONS (+ NUM-ITERATIONS 1))
(if (null? workload)
(begin
;; REMOVE THIS FOR A PRACTICAL APPLICATION.
(display NUM-ITERATIONS)
(display "\n")
params)
(let ((pportion (caar workload))
(qportion (cdar workload))
(workload (cdr workload)))
(if (not (rectangles-overlap (endpt0@ pportion)
(endpt1@ pportion)
(endpt0@ qportion)
(endpt1@ qportion)))
(start-here workload params)
(if (flat-enough? pportion norm rtol atol)
(if (flat-enough? qportion norm rtol atol)
(linear-linear pportion qportion
workload params)
(linear-quadratic pportion qportion
workload params #f))
(if (flat-enough? qportion norm rtol atol)
(linear-quadratic qportion pportion
workload params #t)
(bisect-both pportion qportion
workload params)))))))
(define (linear-linear pportion qportion workload params)
;; Both portions are flat enough to treat as linear. Treat
;; them as line segments and see if the segments intersect.
;; Include the intersection if they do.
;; REMOVE THIS FOR A PRACTICAL APPLICATION.
(display "linear-linear\n")
(let-values (((tp tq)
(segment-parameters (endpt0@ pportion)
(endpt1@ pportion)
(endpt0@ qportion)
(endpt1@ qportion))))
(if tp
(let ((tp (lerp tp (t0@ pportion) (t1@ pportion)))
(tq (lerp tq (t0@ qportion) (t1@ qportion))))
(start-here workload (include-params tp tq params)))
(start-here workload params))))
(define (linear-quadratic pportion qportion workload params
reversed-roles?)
;; Only pportion is flat enough to treat as linear. Find its
;; intersections with the quadratic qportion, and include
;; either or both if they are within the correct
;; intervals. (Use the qportion argument instead of Q
;; directly, so the roles can be reversed for
;; quadratic-linear.)
;;
;; The following Maxima commands will supply formulas for
;; finding values of t for a quadratic in s-power basis
;; plugged into a line in s-power (or Bernstein) basis:
;;
;; /* The line. */
;; xp(t) := px0*(1-t) + px1*t$
;; yp(t) := py0*(1-t) + py1*t$
;;
;; /* The quadratic (s-power basis). */
;; xq(t) := qx0*(1-t) + qx1*t*(1-t) + qx2*t$
;; yq(t) := qy0*(1-t) + qy1*t*(1-t) + qy2*t$
;;
;; /* Implicitize and plug in. */
;; impl(t) := resultant(xq(t)-xp(tau), yq(t)-yp(tau), tau)$
;; expand(impl(t));
;;
;; REMOVE THIS FOR A PRACTICAL APPLICATION.
(if reversed-roles?
(display "quadratic-linear\n")
(display "linear-quadratic\n"))
(let* ((p0 (endpt0@ pportion))
(p1 (endpt1@ pportion))
(QX (plane-curve-x (curve@ qportion)))
(QY (plane-curve-y (curve@ qportion)))
(px0 (car p0)) (py0 (cdr p0))
(px1 (car p1)) (py1 (cdr p1))
(qx0 (spower-c0 QX)) (qy0 (spower-c0 QY))
(qx1 (spower-c1 QX)) (qy1 (spower-c1 QY))
(qx2 (spower-c2 QX)) (qy2 (spower-c2 QY))
(px0×py1 (* px0 py1))
(px1×py0 (* px1 py0))
(px0×qy0 (* px0 qy0))
(px0×qy1 (* px0 qy1))
(px0×qy2 (* px0 qy2))
(px1×qy0 (* px1 qy0))
(px1×qy1 (* px1 qy1))
(px1×qy2 (* px1 qy2))
(py0×qx0 (* py0 qx0))
(py0×qx1 (* py0 qx1))
(py0×qx2 (* py0 qx2))
(py1×qx0 (* py1 qx0))
(py1×qx1 (* py1 qx1))
(py1×qx2 (* py1 qx2))
;; Construct the equation A×t² + B×t + C = 0 and solve
;; it by the quadratic formula.
(A (+ px1×qy1 (- px0×qy1) (- py1×qx1) py0×qx1))
(B (+ (- px1×qy2) px0×qy2 (- px1×qy1) px0×qy1
px1×qy0 (- px0×qy0) py1×qx2 (- py0×qx2)
py1×qx1 (- py0×qx1) (- py1×qx0) py0×qx0))
(C (+ (- px1×qy0) px0×qy0 py1×qx0 (- py0×qx0)
(- px0×py1) px1×py0))
(discr (- (* B B) (* 4 A C))))
(define (invert-param tq)
;;
;; If one of the t-parameter solutions to the quadratic is
;; in [0,1], invert the corresponding (x,y)-coordinates,
;; to find the corresponding t-parameter solution for the
;; linear. If it is in [0,1], then the (x,y)-coordinates
;; are an intersection point. Return that corresponding
;; t-parameter. Otherwise return #f.
;;
(and (0<=x<=1 tq)
(let ((dx (- px1 px0))
(dy (- py1 py0)))
(if (>= (abs dx) (abs dy))
(let* ((x (spower-eval QX tq))
(tp (/ (- x px0) dx)))
(and (0<=x<=1 tp)
(lerp tp (t0@ pportion)
(t1@ pportion))))
(let* ((y (spower-eval QY tq))
(tp (/ (- y py0) dy)))
(and (0<=x<=1 tp)
(lerp tp (t0@ pportion)
(t1@ pportion))))))))
(unless (negative? discr)
(let* ((rootd (sqrt discr))
(tq1 (/ (+ (- B) rootd) (+ A A)))
(tq2 (/ (- (- B) rootd) (+ A A)))
(tp1 (invert-param tq1))
(tp2 (invert-param tq2)))
(when tp1
(set! params (if reversed-roles?
(include-params tq1 tp1 params)
(include-params tp1 tq1 params))))
(when tp2
(set! params (if reversed-roles?
(include-params tq2 tp2 params)
(include-params tp2 tq2 params))))))
(start-here workload params)))
(define (bisect-both pportion qportion workload params)
;; Neither portion is flat enough to treat as linear. Bisect
;; them and add the four new portion-pairs to the workload.
(let-values (((pport1 pport2) (bisect-portion pportion))
((qport1 qport2) (bisect-portion qportion)))
(start-here `((,pport1 . ,qport1) (,pport1 . ,qport2)
(,pport2 . ,qport1) (,pport2 . ,qport2)
. ,workload)
params)))
(start-here (initial-workload P Q) '()))
)) ;; end library (spower quadratic intersections)
(import (scheme base)
(scheme write)
(spower quadratic)
(spower quadratic intersections)
(srfi 132) ;; R7RS-large name: (scheme sort)
)
(define P (control-points->plane-curve '((-1 0) (0 10) (1 0))))
(define P1 (control-points->plane-curve '((-1 0) (-1/2 5) (33 50))))
(define P2 (control-points->plane-curve '((0 5) (1/2 5) (1 0))))
(define Q (control-points->plane-curve '((2 1) (-8 2) (2 3))))
;; Sort the results by the parameters for P. Thus the intersection
;; points will be sorted left to right.
(define params (list-sort (lambda (tpair1 tpair2)
(< (car tpair1) (car tpair2)))
(find-intersection-parameters P Q)))
(for-each
(lambda (pair)
(let ((tp (car pair))
(tq (cdr pair)))
(let-values (((ax ay) (plane-curve-eval P tp))
((bx by) (plane-curve-eval Q tq)))
(display (inexact tp)) (display "\t(")
(display (inexact ax)) (display ", ")
(display (inexact ay)) (display ") \t")
(display (inexact tq)) (display "\t(")
(display (inexact bx)) (display ", ")
(display (inexact by)) (display ")\n"))))
params)
(newline)
(define params (list-sort (lambda (tpair1 tpair2)
(< (car tpair1) (car tpair2)))
(find-intersection-parameters P1 Q)))
(for-each
(lambda (pair)
(let ((tp (car pair))
(tq (cdr pair)))
(let-values (((ax ay) (plane-curve-eval P1 tp))
((bx by) (plane-curve-eval Q tq)))
(display (inexact tp)) (display "\t(")
(display (inexact ax)) (display ", ")
(display (inexact ay)) (display ") \t")
(display (inexact tq)) (display "\t(")
(display (inexact bx)) (display ", ")
(display (inexact by)) (display ")\n"))))
params)
(newline)
(define params (list-sort (lambda (tpair1 tpair2)
(< (car tpair1) (car tpair2)))
(find-intersection-parameters Q P1)))
(for-each
(lambda (pair)
(let ((tp (car pair))
(tq (cdr pair)))
(let-values (((ax ay) (plane-curve-eval Q tp))
((bx by) (plane-curve-eval P1 tq)))
(display (inexact tp)) (display "\t(")
(display (inexact ax)) (display ", ")
(display (inexact ay)) (display ") \t")
(display (inexact tq)) (display "\t(")
(display (inexact bx)) (display ", ")
(display (inexact by)) (display ")\n"))))
params)
(newline)
(define params (list-sort (lambda (tpair1 tpair2)
(< (car tpair1) (car tpair2)))
(find-intersection-parameters P2 Q)))
(for-each
(lambda (pair)
(let ((tp (car pair))
(tq (cdr pair)))
(let-values (((ax ay) (plane-curve-eval P2 tp))
((bx by) (plane-curve-eval Q tq)))
(display (inexact tp)) (display "\t(")
(display (inexact ax)) (display ", ")
(display (inexact ay)) (display ") \t")
(display (inexact tq)) (display "\t(")
(display (inexact bx)) (display ", ")
(display (inexact by)) (display ")\n"))))
params)
- Output:
It turns out that the only bisection branch being executed, for the posed problem, was the one where both portions are bisected. Thus, in the new implementation here, the quadratic formula is never employed. However, fiddling with the numbers turns up problems for which this is not the case.
Also I have added a fourth problem, in which the first the convex-up parabola is cut in half. Interestingly, it takes more iterations to solve that problem than the original! Note, though, that I have not categorized iterations according to their kinds.
$ gosh -r7 bezier_intersections_2.scm linear-linear linear-linear linear-linear linear-linear linear-linear linear-linear 217 0.07250828117222911 (-0.8549834376555417, 1.3450166066735616) 0.17250829973466453 (-0.8549837251463935, 1.345016599469329) 0.15948753092173137 (-0.6810249381565372, 2.6810251680444233) 0.8405124690782686 (-0.6810251680444231, 2.681024938156537) 0.8274917002653355 (0.654983400530671, 2.8549837251463934) 0.9274917188277709 (0.6549833933264384, 2.854983437655542) 0.9405124666929737 (0.8810249333859473, 1.118975333761435) 0.059487533307026316 (0.881024666238565, 1.1189750666140525) quadratic-linear quadratic-linear quadratic-linear quadratic-linear quadratic-linear quadratic-linear 404 0.09485068481613748 (-0.6082597856508847, 1.3083729445649852) 0.15418647228249246 (-0.6082600809514531, 1.3083729445649848) 0.1670105772077877 (0.0874641628839754, 2.7858070880490136) 0.892903544024507 (0.08746389814035438, 2.7858070880490144) linear-quadratic linear-quadratic linear-quadratic linear-quadratic linear-quadratic linear-quadratic 591 0.15418647228249246 (-0.6082600809514531, 1.3083729445649848) 0.09485068481613748 (-0.6082597856508847, 1.3083729445649852) 0.892903544024507 (0.08746389814035438, 2.7858070880490144) 0.1670105772077877 (0.0874641628839754, 2.7858070880490136) linear-linear linear-linear linear-linear 702 0.654983400530671 (0.654983400530671, 2.8549837251463934) 0.9274917188277709 (0.6549833933264384, 2.854983437655542) 0.8810249333859473 (0.8810249333859473, 1.118975333761435) 0.059487533307026316 (0.881024666238565, 1.1189750666140525)
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 seemsToBeDuplicate = Fn.new { |intersects, xy, spacing|
var seemsToBeDup = false
var i = 0
while (!seemsToBeDup && i != intersects.count) {
var pt = intersects[i]
seemsToBeDup = (pt.x - xy.x).abs < spacing && (pt.y - xy.y).abs < spacing
i = i + 1
}
return seemsToBeDup
}
var findIntersects = Fn.new { |p, q, tol, spacing|
var intersects = []
var workload = Stack.new()
workload.push(Workset.new(p, q))
// Quit looking after having emptied the workload.
while (!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) {
// To avoid detecting the same intersection twice, require some
// space between intersections.
if (!seemsToBeDuplicate.call(intersects, intersect, spacing)) {
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 tol = 0.0000001
var spacing = 10 * tol
var intersects = findIntersects.call(p, q, tol, spacing)
for (intersect in intersects) Fmt.print("($ f, $f)", intersect.x, intersect.y)- Output:
( 0.654983, 2.854983) ( 0.881025, 1.118975) (-0.681025, 2.681025) (-0.854983, 1.345017)