Bézier curves/Intersections: Difference between revisions

Content added Content deleted

Inline

Revision as of 20:13, 10 June 2023

You are given two planar quadratic Bézier curves, having control points $(-1,0),(0,10),(1,0)$ and $(2,1),(-8,2),(2,3)$ , respectively. They are parabolas intersecting at four points, as shown in the following diagram:

Two intersecting parabolas and their control points.

The task is to write a program that finds all four intersection points and prints their $(x,y)$ coordinates. You may use any algorithm you know of or can think of, including any of those that others have used.

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

Translation of: D

/* 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

Translation of: D

/* 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 $t$ -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.)

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

Works with: GCC version 13.1.1

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)

Pascal

Translation of: Icon

This is the algorithm of the Icon example, recast as a recursive procedure.

Phix

Translation of: D

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}}

Python

Translation of: Icon

Translation of: Pascal

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.)

#!/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 = None):
        '''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)

Wren

Translation of: D

Library: Wren-dynamic

Library: Wren-trait

Library: Wren-math

Library: Wren-assert

Library: Wren-seq

Library: Wren-fmt

/* 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)

@@ Line 1,224: / Line 1,224: @@
 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]]. 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:
+* [[#Pascal|Pascal]]
+* [[#Python|Python]]
 <syntaxhighlight lang="icon">