Jump to content

Bézier curves/Intersections

From Rosetta Code
Task
Bézier curves/Intersections
You are encouraged to solve this task according to the task description, using any language you may know.

You are given two planar quadratic Bézier curves, having control points and , respectively. They are parabolas intersecting at four points, as shown in the following diagram:


Two intersecting parabolas and their control points.


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

Ada

This is a variant of the Bézier clipping method of Thomas W. Sederberg's lecture notes on "Computer Aided Geometric Design", using symmetric power (s-power) polynomials instead of the usual "control points". I do not know that it is any better than the method of the s-power polynomial flattening, which is done on this page in several languages. The flattening method is considerably simpler. Also, I have doubts about the numerical stability of the clipping method.

On the other hand, the code below introduces a new way of doing Bézier clipping, by transformation of coordinates. And it also hints at a geometric interpretation of s-power curves, in which all the constituent "points" undergo rotations and scalings, but only the endpoints are subject to translations.

Ada is an Algol-like language that is case-insensitive. I have used an all-lowercase style with manual indentation. Subprogram names are overloaded. Subprogram parameters may be passed either positionally or by name. Ada does have a few oddities, among which are that the term for a pointer is not "pointer" but "access".

pragma wide_character_encoding (utf8);
pragma ada_2022;
pragma assertion_policy (check);

--
-- The algorithm here is iterated "Bézier clipping" [1] of quadratic
-- Bézier curves in a plane. However, the curves will be represented
-- in the s-power basis [2,3], and there will be these other
-- differences from the algorithms suggested by [1]:
--
-- * Before clipping, the problem will be moved to a coordinate system
--   in which the "skeleton" of the clipper curve—the line segment
--   between its endpoints—extends from (0,0) and (1,0). In this
--   coordinate system, the "distance function" is the y-coordinate.
--
-- * The fat line will be the best possible parallel to the "skeleton"
--   of the clipper curve. This is easily done in the s-power basis,
--   once one has transformed coordinate systems as described above:
--   the center coefficient of the (y,t)-polynomial is exactly four
--   times the optimum signed fat line width.
--
-- * Intersections with fat line boundaries will be found not by
--   extending the sides of control polygons, but instead by the
--   quadratic formula.
--
-- * Mechanisms to ensure termination of the algorithm, and to improve
--   the chances of finding intersections, may differ.
--
-- Inputs and outputs to the algorithms will be in single
-- precision. However, the internals will be in double precision. This
-- decision helps resolve some minor numerical quandaries.
--
-- References:
--
-- [1] Sederberg, Thomas W., "Computer Aided Geometric Design"
--     (2012). Faculty
--     Publications. 1. https://scholarsarchive.byu.edu/facpub/1
--     http://hdl.lib.byu.edu/1877/2822
--
-- [2] J. Sánchez-Reyes, ‘The symmetric analogue of the polynomial
--     power basis’, ACM Transactions on Graphics, vol 16 no 3, July
--     1997, page 319.
--
-- [3] 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.
--
-- Suggestions for the future:
--
-- * Combination of Bézier clipping and flatness testing. Remember
--   from other Rosetta Code examples, that in the s-power basis it is
--   simple to test the "flatness" of a parametric plane
--   curve. Suppose two curves are running nearly parallel to each
--   other, so they are difficult to reduce by clipping. Nevertheless,
--   they can be made flatter and flatter, so eventually they become
--   "segment-like" and can be tested by a line-intersection routine.
--

with ada.text_io;
use ada.text_io;

with ada.float_text_io;
use ada.float_text_io;

with ada.numerics.elementary_functions;
use ada.numerics.elementary_functions;

with ada.numerics.long_elementary_functions;
use ada.numerics.long_elementary_functions;

with ada.unchecked_deallocation;

procedure bezier_intersections is

  subtype real is float;
  subtype lreal is long_float;

  function between_0_and_1 (x : real)
  return boolean is
  begin
    return 0.0 <= x and x <= 1.0;
  end between_0_and_1;

  function between_0_and_1 (x : lreal)
  return boolean is
  begin
    return 0.0 <= x and x <= 1.0;
  end between_0_and_1;

  procedure quadratic_formula (a, b, c     : lreal;
                               has_smaller : out boolean;
                               smaller     : out lreal;
                               has_larger  : out boolean;
                               larger      : out lreal)
  with post => has_smaller or not has_larger
  is
    --
    -- This routine is based on QUADPL from the Naval Surface Warfare
    -- Center mathematical library nswc.f90 https://archive.ph/2IpWb
    --
    -- We are interested only in real roots. Be aware that, if there
    -- are two roots, they might be a double root.
    --

    bh, e, d        : lreal;

    -- Avogadro’s number, by definition (since 2019):
    arbitrary_value : constant lreal := 6.012214076e23;
  begin
    has_smaller := false;
    has_larger := false;
    smaller := arbitrary_value;
    larger := arbitrary_value;

    if a = 0.0 then
      if b = 0.0 then
        -- The equation is a constant equation. I do not know what to
        -- do with that. There may be no solution, or there may be
        -- infinitely many solutions. Let us treat it as having no
        -- USEFUL solutions.
        null;
      else
        -- The equation is linear. There truly is only one root.
        has_smaller := true;
        smaller := -c / b;
      end if;
    elsif c = 0.0 then
      -- The equation is quadratic, but the constant term
      -- vanishes. One of the roots is trivially zero.
      has_smaller := true;
      has_larger := true;
      smaller := 0.0;
      larger := -b / a;
    else
      -- Compute the discriminant, avoiding overflow.
      bh := b / 2.0;
      if abs (bh) >= abs (c) then
        e := 1.0 - (a / bh) * (c / bh);
        d := sqrt (abs (e)) * abs(bh);
      else
        e := (if c < 0.0 then -a else a);
        e := bh * (bh / abs (c)) - e;
        d := sqrt (abs (e)) * sqrt (abs (c));
      end if;
      if e < 0.0 then
        -- The solutions are complex conjugates, but we are interested
        -- only in real solutions.
        null;
      else
        has_smaller := true;
        has_larger := true;
        if bh >= 0.0 then
          d := -d;
        end if;
        larger := (-bh + d) / a;
        if larger = 0.0 then
          smaller := 0.0;       -- Obviously a double root.
        else
          smaller := (c / larger) / a;
        end if;
      end if;
    end if;
  end quadratic_formula;

  type point is
    record
      x, y : real;
    end record;

  type control_point_curve is
    record
      pt0, pt1, pt2 : point;
    end record;

  type spower_poly is
    record
      c0, c1, c2 : real;
    end record;

  type spower_curve is
    record
      x, y : spower_poly;
    end record;

  function eval (poly : spower_poly;
                 t    : real)
  return real is
    c0, c1, c2 : real;
  begin
    c0 := poly.c0;
    c1 := poly.c1;
    c2 := poly.c2;
    return (if t <= 0.5 then
              c0 + t * ((c2 - c0) + ((1.0 - t) * c1))
            else
              c2 + (1.0 - t) * ((c0 - c2) + (t * c1)));
  end eval;

  function eval (curve : spower_curve;
                 t     : real)
  return point is
  begin
    return (x => eval (curve.x, t),
            y => eval (curve.y, t));
  end eval;

  function to_spower_curve (ctrl : control_point_curve)
  return spower_curve is
    curve    : spower_curve;
    p0x, p0y : long_float;
    p1x, p1y : long_float;
    p2x, p2y : long_float;
  begin
    curve.x.c0 := ctrl.pt0.x;
    curve.y.c0 := ctrl.pt0.y;

    curve.x.c2 := ctrl.pt2.x;
    curve.y.c2 := ctrl.pt2.y;

    p0x := long_float (ctrl.pt0.x);
    p0y := long_float (ctrl.pt0.y);
    p1x := long_float (ctrl.pt1.x);
    p1y := long_float (ctrl.pt1.y);
    p2x := long_float (ctrl.pt2.x);
    p2y := long_float (ctrl.pt2.y);
    curve.x.c1 := float (p1x + p1x - p0x - p2x);
    curve.y.c1 := float (p1y + p1y - p0y - p2y);

    return curve;
  end to_spower_curve;

  type opoint is
    record
      -- A point relative to the origin.
      x, y : lreal;
    end record;

  type dpoint is
    record
      -- A point relative to some other point.
      dx, dy : lreal;
    end record;

  function "-" (pt1, pt2 : opoint)
  return dpoint is
  begin
    return (dx => pt1.x - pt2.x,
            dy => pt1.y - pt2.y);
  end "-";

  function "+" (opt : opoint;
                dpt : dpoint)
  return opoint is
  begin
    return (x => opt.x + dpt.dx,
            y => opt.y + dpt.dy);
  end "+";

  function "+" (pt1, pt2 : dpoint)
  return dpoint is
  begin
    return (dx => pt1.dx + pt2.dx,
            dy => pt1.dy + pt2.dy);
  end "+";

  function "*" (scalar : lreal;
                pt     : dpoint)
  return dpoint is
  begin
    return (dx => scalar * pt.dx,
            dy => scalar * pt.dy);
  end "*";

  type geocurve is
    record
      --
      -- A plane curve in s-power basis, represented geometrically.
      --
      -- pt0 and pt2 form the curve’s ‘skeleton’. Their convex combination,
      -- ((1.0-t)*pt0)+(t*pt1), form the linear backbone of the curve.
      --
      -- (1.0-t)*t*pt1 forms (x,t) and (y,t) parabolas that are
      -- convex-up/down, and which (because this is an s-power curve)
      -- vanish as the curve becomes linear.
      --
      -- The actual geocurve is the locus of (1.0-t)*t*pt1 relative to
      -- the skeleton, respectively for each value of t in [0,1].
      --
      -- If [t0,t1] is not [0.0,1.0], then this is a portion of some
      -- parent curve that is determined by context.
      --
      -- (I feel as if all this can be crafted nicely in terms of a
      -- geometric algebra, with special points not only at infinity
      -- but also at the origin, but let us not worry about that for
      -- now. Likely the calculations would not change much, but the
      -- names of things might.)
      --
      pt0    : opoint;
      pt1    : dpoint;
      pt2    : opoint;
      t0, t1 : lreal;
    end record;

  function to_geocurve (scurve : spower_curve)
  return geocurve is
    gcurve : geocurve;
  begin
    gcurve.pt0.x := lreal (scurve.x.c0);
    gcurve.pt0.y := lreal (scurve.y.c0);
    gcurve.pt1.dx := lreal (scurve.x.c1);
    gcurve.pt1.dy := lreal (scurve.y.c1);
    gcurve.pt2.x := lreal (scurve.x.c2);
    gcurve.pt2.y := lreal (scurve.y.c2);
    gcurve.t0 := 0.0;
    gcurve.t1 := 1.0;
    return gcurve;
  end to_geocurve;

  function eval (gcurve : geocurve;
                 t      : lreal)
  return opoint is
    pt0 : opoint;
    pt1 : dpoint;
    pt2 : opoint;
  begin
    pt0 := gcurve.pt0;
    pt1 := gcurve.pt1;
    pt2 := gcurve.pt2;
    return (if t <= 0.5 then
              pt0 + t * ((pt2 - pt0) + ((1.0 - t) * pt1))
            else
              pt2 + (1.0 - t) * ((pt0 - pt2) + (t * pt1)));
  end eval;

  type affine_projection is
    record
      --
      -- An affine transformation used to project objects into the
      -- coordinate system of what we are calling a "skeleton".
      --
      -- If it is an opoint, first add transpose[dx,dy]. If it is a
      -- dpoint, ignore dx,dy. Then, in either case, premultiply by
      -- [[a,b],[c,d]]. (One can view dx,dy as a dpoint.)
      --
      dx, dy : lreal;
      a, b   : lreal;
      c, d   : lreal;
    end record;

  function project (pt    : opoint;
                    trans : affine_projection)
  return opoint is
    x, y : lreal;
  begin
    -- Apply an affine transformation to an opoint, projecting it from
    -- one coordinate system to another.
    x := pt.x + trans.dx;       -- The origin may have moved.
    y := pt.y + trans.dy;
    return (x => trans.a * x + trans.b * y, -- Rotation, scaling, etc.
            y => trans.c * x + trans.d * y);
  end project;

  function project (pt    : dpoint;
                    trans : affine_projection)
  return dpoint is
    x, y : lreal;
  begin
    -- Apply an affine transformation to a dpoint, projecting it from
    -- one coordinate system to another.
    x := pt.dx;          -- With dpoint, translation is not performed.
    y := pt.dy;
    return (dx => trans.a * x + trans.b * y, -- Rotation, scaling, &c.
            dy => trans.c * x + trans.d * y);
  end project;

  function project (gcurve : geocurve;
                    trans  : affine_projection)
  return geocurve is
  begin
    -- Apply an affine transformation to a geocurve, projecting it
    -- from one coordinate system to another.
    return (pt0 => project (gcurve.pt0, trans),
            pt1 => project (gcurve.pt1, trans),
            pt2 => project (gcurve.pt2, trans),
            t0 => 0.0, t1 => 1.0);
  end project;

  function segment_projection_to_0_1 (pt0, pt1 : opoint)
  return affine_projection is
    a, b  : lreal;
    denom : lreal;
  begin
    -- Return the transformation that projects pt0 to (0,0) and pt1 to
    -- (1,0) without any distortions except translation, rotation, and
    -- scaling.
    a := pt1.x - pt0.x;
    b := pt1.y - pt0.y;
    denom := (a * a) + (b * b);
    pragma assert (denom /= 0.0);
    a := a / denom;
    b := b / denom;
    return (dx => -pt0.x, dy => -pt0.y,
            a => a,  b => b,
            c => -b, d => a);
  end segment_projection_to_0_1;

  procedure test_segment_projection_to_0_1 is
    trans : affine_projection;
    pt    : opoint;
  begin
    -- Some unit testing.

    trans := segment_projection_to_0_1 ((1.0, 2.0), (-3.0, 5.0));
    pt := project ((1.0, 2.0), trans);
    pragma assert (abs (pt.x) < 1.0e-10);
    pragma assert (abs (pt.y) < 1.0e-10);
    pt := project ((-3.0, 5.0), trans);
    pragma assert (abs (pt.x - 1.0) < 1.0e-10);
    pragma assert (abs (pt.y) < 1.0e-10);

    trans := segment_projection_to_0_1 ((0.0, 2.0), (0.0, -50.0));
    pt := project ((0.0, 2.0), trans);
    pragma assert (abs (pt.x) < 1.0e-10);
    pragma assert (abs (pt.y) < 1.0e-10);
    pt := project ((0.0, -50.0), trans);
    pragma assert (abs (pt.x - 1.0) < 1.0e-10);
    pragma assert (abs (pt.y) < 1.0e-10);
  end test_segment_projection_to_0_1;

  function skeleton_projection (clipper : geocurve)
  return affine_projection is
  begin
    -- Return the transformation that projects the "skeleton" of the
    -- clipper into a coordinate system where the skeleton goes from
    -- the origin (0,0) to the point (1,0).
    return segment_projection_to_0_1 (clipper.pt0, clipper.pt2);
  end skeleton_projection;

  procedure fat_line_bounds (projected_clipper : geocurve;
                             dmin, dmax        : out lreal) is
    d       : lreal;
    d1, d2  : lreal;
    padding : constant lreal := lreal (real'model_epsilon) / 1000.0;
  begin
    --
    -- (1-t)*t is bounded on [0,1] by 0 and 1/4. Thus we have our fat
    -- line boundaries. We leave some padding, however, for numerical
    -- safety. The padding is much smaller than our tolerance. (We are
    -- calculating in double precision, but to obtain results in
    -- single precision.)
    --
    -- The fat line bounds will be y=dmin and y=dmax. If you ignore
    -- the padding, the fat line is the best possible among those
    -- parallel to the skeleton. It is likely to be much tighter than
    -- what you would infer from the Bernstein-basis control polygon.
    --

    d := projected_clipper.pt1.dy / 4.0;
    d1 := -padding * d;
    d2 := (1.0 + padding) * d;
    dmin := lreal'min (d1, d2);
    dmax := lreal'max (d1, d2);
  end fat_line_bounds;

  type clipping_places is
    (clip_everywhere,
     clip_nowhere,
     clip_left_and_right);

  type clipping_plan is
    record
      where  : clipping_places;
      t_lft  : lreal;
      t_rgt  : lreal;
    end record;

  function make_clipping_plan (projected_clippehend : geocurve;
                               dmin, dmax           : lreal)
  return clipping_plan is
    ncrossings  : integer range 0 .. 4;
    tcrossings  : array (1 .. 4) of lreal;
    dcrossings  : array (1 .. 4) of lreal;
    has_smaller : boolean;
    has_larger  : boolean;
    smaller     : lreal;
    larger      : lreal;
    c0, c1, c2  : lreal;

    function outside_of_bounds
    return boolean is
      dy, ymin, ymax : lreal;
    begin
      -- (1-t)*t is bounded on [0,1] by 0 and 1/4. Use that fact to
      -- make padding for a simple (although not especially tight)
      -- bounds on y.
      dy := abs (projected_clippehend.pt1.dy) / 4.0;
      ymin := -dy + lreal'min (projected_clippehend.pt0.y,
                               projected_clippehend.pt2.y);
      ymax := dy + lreal'max (projected_clippehend.pt0.y,
                              projected_clippehend.pt2.y);
      return (ymax <= dmin or dmax <= ymin);
    end outside_of_bounds;

    procedure find_crossings (yvalue : lreal) is
    begin
      quadratic_formula (-c1, c1 + c2 - c0, c0 - yvalue,
                         has_smaller, smaller, has_larger, larger);
      if has_smaller then
        if between_0_and_1 (smaller) then
          ncrossings := ncrossings + 1;
          tcrossings(ncrossings) := smaller;
          dcrossings(ncrossings) := yvalue;
        end if;
      end if;
      if has_larger then
        if between_0_and_1 (larger) then
          ncrossings := ncrossings + 1;
          tcrossings(ncrossings) := larger;
          dcrossings(ncrossings) := yvalue;
        end if;
      end if;
    end find_crossings;

    function create_a_clipping_plan
    return clipping_plan is
      plan : clipping_plan;
    begin
      --
      -- For simplicity, there is just one situation in which we will
      -- clip: when there is one crossing of dmin and one crossing of
      -- dmax.
      --
      if ncrossings = 2 and then dcrossings(1) /= dcrossings(2) then
        plan.where := clip_left_and_right;
        plan.t_lft := lreal'min (tcrossings(1), tcrossings(2));
        plan.t_rgt := lreal'max (tcrossings(1), tcrossings(2));
      else
        plan.where := clip_nowhere;
        -- Arbitrary values:
        plan.t_lft := 0.0;
        plan.t_rgt := 0.0;
      end if;
      return plan;
    end create_a_clipping_plan;

    function analyze_and_plan
    return clipping_plan is
    begin
      c0 := projected_clippehend.pt0.y;
      c1 := projected_clippehend.pt1.dy;
      c2 := projected_clippehend.pt2.y;
      ncrossings := 0;
      find_crossings (dmin);
      find_crossings (dmax);
      return create_a_clipping_plan;
    end analyze_and_plan;

  begin
    return (if outside_of_bounds then
              (where => clip_everywhere,
               -- Arbitrary values:
               t_lft => 0.0, t_rgt => 0.0)
            else
              analyze_and_plan);
  end make_clipping_plan;

  function make_portion (gcurve : geocurve;
                         t0, t1 : lreal)
  return geocurve is
    pt0    : opoint;
    pt1    : dpoint;
    pt2    : opoint;
  begin
    pt0 := eval (gcurve, t0);
    pt1 := ((t1 - t0 - t0) * t1 + (t0 * t0)) * gcurve.pt1;
    pt2 := eval (gcurve, t1);
    return (pt0 => pt0, pt1 => pt1, pt2 => pt2,
            t0 => t0, t1 => t1);
  end make_portion;

  procedure segment_parameters (a0, a1 : opoint;
                                b0, b1 : opoint;
                                ta, tb : out lreal) is
    axdiff : lreal;
    aydiff : lreal;
    bxdiff : lreal;
    bydiff : lreal;
    denom  : lreal;
    ta1    : lreal;
    tb1    : lreal;
  begin
    -- Do the line segments (a0,a1) and (b0,b1) intersect?  If so,
    -- return t-parameter values in ta,tb in [0,1] for the point of
    -- intersection, treating them as parametric splines of degree
    -- 1. Otherwise return ta=-1.0 and tb=-1.0.

    ta := -1.0;
    tb := -1.0;

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

    ta1 := ((bxdiff * a0.y) - (bydiff * a0.x)
               + (b0.x * b1.y) - (b1.x * b0.y));
    ta1 := ta1 / denom;
    if between_0_and_1 (ta1) then
      tb1 := -((axdiff * b0.y) - (aydiff * b0.x)
                  + (a0.x * a1.y) - (a1.x * a0.y));
      tb1 := tb1 / denom;
      if between_0_and_1 (tb1) then
        ta := ta1;
        tb := tb1;
      end if;
    end if;
  end segment_parameters;

  type t_pair;
  type t_pair_access is access t_pair;
  type t_pair is
    record
      tp, tq : real;            -- Single precision results.
      next   : t_pair_access;
    end record;
  procedure t_pair_delete is
    new ada.unchecked_deallocation (t_pair, t_pair_access);

  function length (params : t_pair_access)
  return natural is
    n : natural;
    p : t_pair_access;
  begin
    n := 0;
    p := params;
    while p /= null loop
      n := n + 1;
      p := p.next;
    end loop;
    return n;
  end length;

  procedure insert_t_pair (params : in out t_pair_access;
                           tp, tq : real)
  with pre => between_0_and_1 (tp) and between_0_and_1 (tq)
  is
    pair : t_pair_access;
  begin
    -- Test for duplicates and also sort as we insert. This is a
    -- recursive implementation.
    if params = null then
      pair := new t_pair;
      params := pair;
      pair.next := null;
      pair.tp := tp;
      pair.tq := tq;
    elsif abs (params.tp - tp) <= real'model_epsilon then
      null;                     -- A duplicate, to within epsilon.
    elsif params.tp < tp then
      insert_t_pair (params.next, tp, tq);
    else
      pair := new t_pair;
      pair.next := params.next;
      params.next := pair;
      pair.tp := params.tp;
      pair.tq := params.tq;
      params.tp := tp;
      params.tq := tq;
    end if;
  end insert_t_pair;

  procedure insert_t_pair (params : in out t_pair_access;
                           tp, tq : lreal)
  with pre => between_0_and_1 (tp) and between_0_and_1 (tq)
  is
  begin
    insert_t_pair (params, real (tp), real (tq));
  end insert_t_pair;

  type clipping_roles is
    (pportion_clips_qportion,
     qportion_clips_pportion);

  type intersection_job;
  type intersection_job_access is access intersection_job;
  type intersection_job is
    record
      pportion : geocurve;
      qportion : geocurve;
      roles    : clipping_roles;
      next     : intersection_job_access;
    end record;
  procedure intersection_job_delete is
    new ada.unchecked_deallocation (intersection_job,
                                    intersection_job_access);

  function find_intersections (p, q : geocurve)
  return t_pair_access is
    params     : t_pair_access;
    workload   : intersection_job_access;
    pportion   : geocurve;
    qportion   : geocurve;
    roles      : clipping_roles;
    tp, tq     : lreal;

    tolerance : constant lreal := 0.4 * lreal (real'model_epsilon);
    bisection_threshold : constant lreal := 0.2;

    procedure insert_job (pportion : geocurve;
                          qportion : geocurve;
                          roles    : clipping_roles)
    with post => workload /= null
    is
      job : intersection_job_access;
    begin
      job := new intersection_job;
      job.pportion := pportion;
      job.qportion := qportion;
      job.roles := roles;
      job.next := workload;
      workload := job;
    end insert_job;

    procedure insert_job (pportion : geocurve;
                          qportion : geocurve)
    with post => workload /= null
    is
      roles : clipping_roles;
    begin
      if pportion.t1 - pportion.t0 <= qportion.t1 - qportion.t0 then
        roles := pportion_clips_qportion;
      else
        roles := qportion_clips_pportion;
      end if;
      insert_job (pportion, qportion, roles);
    end insert_job;

    procedure choose_job
    with pre => workload /= null
    is
      job : intersection_job_access;
    begin
      job := workload;
      workload := workload.next;
      pportion := job.pportion;
      qportion := job.qportion;
      roles := job.roles;
      intersection_job_delete (job);
    end choose_job;

    function t_parameters_are_within_tolerance
    return boolean is
    begin
      return (pportion.t1 - pportion.t0 <= tolerance
                and qportion.t1 - qportion.t0 <= tolerance);
    end t_parameters_are_within_tolerance;

    procedure insert_intersection_parameters is
    begin
      segment_parameters (pportion.pt0, pportion.pt2,
                          qportion.pt0, qportion.pt2, tp, tq);
      if between_0_and_1 (tp) and between_0_and_1 (tq) then
        tp := (1.0 - tp) * pportion.t0 + tp * pportion.t1;
        tq := (1.0 - tq) * qportion.t0 + tq * qportion.t1;
        insert_t_pair (params, tp, tq);
      end if;
    end insert_intersection_parameters;

  procedure no_intersections is
  begin
    -- Do nothing. The reason to have a procedure for this is so the
    -- code will document itself.
    null;
  end no_intersections;

  procedure bisect_them is
    pportion1  : geocurve;
    pportion2  : geocurve;
    qportion1  : geocurve;
    qportion2  : geocurve;
    t0, t1, th : lreal;
  begin
    --
    -- Sederberg employs heuristics to decide what to bisect and who
    -- should be the clipper. We will simply bisect both curves, and
    -- use some heuristic coded elsewhere to choose the
    -- clipper.
    --
    -- (Bisection as it is normally done in the s-power basis is
    -- equivalent to change of homogeneous coordinates. It can be
    -- applied to p and q directly and, for quadratics, is a simple
    -- operation.)
    --

    t0 := pportion.t0;
    t1 := pportion.t1;
    th := 0.5 * t0 + 0.5 * t1;
    pportion1 := make_portion (p, t0, th);
    pportion2 := make_portion (p, th, t1);

    t0 := qportion.t0;
    t1 := qportion.t1;
    th := 0.5 * t0 + 0.5 * t1;
    qportion1 := make_portion (q, t0, th);
    qportion2 := make_portion (q, th, t1);

    insert_job (pportion1, qportion1);
    insert_job (pportion1, qportion2);
    insert_job (pportion2, qportion1);
    insert_job (pportion2, qportion2);
  end bisect_them;

  procedure clip_and_maybe_bisect (plan : clipping_plan) is
    clipping : geocurve;
    t_lft    : lreal;
    t_rgt    : lreal;
    t0, t1   : lreal;
  begin
    t_lft := plan.t_lft;
    t_rgt := plan.t_rgt;
    if roles = pportion_clips_qportion then
      t0 := (1.0 - t_lft) * qportion.t0 + t_lft * qportion.t1;
      t1 := (1.0 - t_rgt) * qportion.t0 + t_rgt * qportion.t1;
      clipping := make_portion (q, t0, t1);
      if 1.0 - (t_rgt - t_lft) < bisection_threshold then
        qportion := clipping;
        bisect_them;
      else
        insert_job (pportion, clipping, qportion_clips_pportion);
      end if;
    else
      t0 := (1.0 - t_lft) * pportion.t0 + t_lft * pportion.t1;
      t1 := (1.0 - t_rgt) * pportion.t0 + t_rgt * pportion.t1;
      clipping := make_portion (p, t0, t1);
      if 1.0 - (t_rgt - t_lft) < bisection_threshold then
        pportion := clipping;
        bisect_them;
      else
        insert_job (clipping, qportion, pportion_clips_qportion);
      end if;
    end if;
  end clip_and_maybe_bisect;

  procedure perform_clipping is
    clippehend, clipper  : geocurve;
    projected_clippehend : geocurve;
    projected_clipper    : geocurve;
    trans                : affine_projection;
    dmin, dmax           : lreal;
    plan                 : clipping_plan;
  begin
    if roles = pportion_clips_qportion then
      clipper := pportion;
      clippehend := qportion;
    else
      clipper := qportion;
      clippehend := pportion;
    end if;
    trans := skeleton_projection (clipper);
    projected_clipper := project (clipper, trans);
    fat_line_bounds (projected_clipper, dmin, dmax);
    projected_clippehend := project (clippehend, trans);
    plan := make_clipping_plan (projected_clippehend, dmin, dmax);
    case plan.where is
      when clip_everywhere => no_intersections;
      when clip_nowhere => bisect_them;
      when clip_left_and_right => clip_and_maybe_bisect (plan);
    end case;
  end perform_clipping;

  begin
    params := null;
    workload := null;
    insert_job (p, q);
    while workload /= null loop
      choose_job;
      if t_parameters_are_within_tolerance then
        insert_intersection_parameters;
      else
        perform_clipping;
      end if;
    end loop;
    return params;
  end find_intersections;

  function find_intersections (p, q : spower_curve)
  return t_pair_access is
  begin
    return find_intersections (to_geocurve (p), to_geocurve (q));
  end find_intersections;

  pctrl, qctrl : control_point_curve;
  p, q         : spower_curve;
  params, ptr  : t_pair_access;

begin
  test_segment_projection_to_0_1;
  pctrl := ((x => -1.0, y =>  0.0),
            (x =>  0.0, y => 10.0),
            (x =>  1.0, y =>  0.0));
  qctrl := ((x =>  2.0, y =>  1.0),
            (x => -8.0, y =>  2.0),
            (x =>  2.0, y =>  3.0));
  p := to_spower_curve (pctrl);
  q := to_spower_curve (qctrl);
  params := find_intersections (p, q);

  set_col (to => 9);
  put ("convex up");
  set_col (to => 44);
  put_line ("convex left");
  ptr := params;
  while ptr /= null loop
    set_col (to => 3);
    put (ptr.tp, fore => 1, aft => 6, exp => 0);
    set_col (to => 13);
    put ("(");
    put (eval (p, ptr.tp).x, fore => 1, aft => 6, exp => 0);
    put (", ");
    put (eval (p, ptr.tp).y, fore => 1, aft => 6, exp => 0);
    put (")");
    set_col (to => 38);
    put (ptr.tq, fore => 1, aft => 6, exp => 0);
    set_col (to => 48);
    put ("(");
    put (eval (q, ptr.tq).x, fore => 1, aft => 6, exp => 0);
    put (", ");
    put (eval (q, ptr.tq).y, fore => 1, aft => 6, exp => 0);
    put_line (")");
    ptr := ptr.next;
  end loop;
end bezier_intersections;
Output:
        convex up                          convex left
  0.072508  (-0.854983, 1.345017)    0.172508  (-0.854983, 1.345017)
  0.159488  (-0.681025, 2.681026)    0.840512  (-0.681026, 2.681025)
  0.827492  (0.654984, 2.854983)     0.927492  (0.654983, 2.854983)
  0.940512  (0.881025, 1.118975)     0.059488  (0.881025, 1.118975)

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)

C#

Translation of: Java
using System;
using System.Collections.Generic;

public class BezierCurveIntersection
{
    public static void Main(string[] args)
    {
        QuadCurve vertical = new QuadCurve(new QuadSpline(-1.0, 0.0, 1.0), new QuadSpline(0.0, 10.0, 0.0));
        // QuadCurve vertical represents the Bezier curve having control points (-1, 0), (0, 10) and (1, 0) 
        QuadCurve horizontal = new QuadCurve(new QuadSpline(2.0, -8.0, 2.0), new QuadSpline(1.0, 2.0, 3.0));
        // QuadCurve horizontal represents the Bezier curve having control points (2, 1), (-8, 2) and (2, 3) 		

        Console.WriteLine("The points of intersection are:");
        List<Point> intersects = FindIntersects(vertical, horizontal);
        foreach (Point intersect in intersects)
        {
            Console.WriteLine($"( {intersect.X,9:0.000000}, {intersect.Y,9:0.000000} )");
        }
    }

    private static List<Point> FindIntersects(QuadCurve p, QuadCurve q)
    {
        List<Point> result = new List<Point>();
        Stack<QuadCurve> stack = new Stack<QuadCurve>();
        stack.Push(p);
        stack.Push(q);

        while (stack.Count > 0)
        {
            QuadCurve pp = stack.Pop();
            QuadCurve qq = stack.Pop();
            List<object> objects = TestIntersection(pp, qq);
            bool accepted = (bool)objects[0];
            bool excluded = (bool)objects[1];
            Point intersect = (Point)objects[2];

            if (accepted)
            {
                if (!SeemsToBeDuplicate(result, intersect))
                {
                    result.Add(intersect);
                }
            }
            else if (!excluded)
            {
                QuadCurve p0 = new QuadCurve();
                QuadCurve q0 = new QuadCurve();
                QuadCurve p1 = new QuadCurve();
                QuadCurve q1 = new QuadCurve();
                SubdivideQuadCurve(pp, 0.5, p0, p1);
                SubdivideQuadCurve(qq, 0.5, q0, q1);
                stack.Push(p0);
                stack.Push(q0);
                stack.Push(p0);
                stack.Push(q1);
                stack.Push(p1);
                stack.Push(q0);
                stack.Push(p1);
                stack.Push(q1);
            }
        }
        return result;
    }

    private static bool SeemsToBeDuplicate(List<Point> intersects, Point point)
    {
        foreach (Point intersect in intersects)
        {
            if (Math.Abs(intersect.X - point.X) < Spacing && Math.Abs(intersect.Y - point.Y) < Spacing)
            {
                return true;
            }
        }
        return false;
    }

    private static List<object> TestIntersection(QuadCurve p, QuadCurve q)
    {
        double pxMin = Math.Min(Math.Min(p.X.C0, p.X.C1), p.X.C2);
        double pyMin = Math.Min(Math.Min(p.Y.C0, p.Y.C1), p.Y.C2);
        double pxMax = Math.Max(Math.Max(p.X.C0, p.X.C1), p.X.C2);
        double pyMax = Math.Max(Math.Max(p.Y.C0, p.Y.C1), p.Y.C2);

        double qxMin = Math.Min(Math.Min(q.X.C0, q.X.C1), q.X.C2);
        double qyMin = Math.Min(Math.Min(q.Y.C0, q.Y.C1), q.Y.C2);
        double qxMax = Math.Max(Math.Max(q.X.C0, q.X.C1), q.X.C2);
        double qyMax = Math.Max(Math.Max(q.Y.C0, q.Y.C1), q.Y.C2);

        bool accepted = false;
        bool excluded = true;
        Point intersect = new Point(0.0, 0.0);

        if (RectanglesOverlap(pxMin, pyMin, pxMax, pyMax, qxMin, qyMin, qxMax, qyMax))
        {
            excluded = false;
            double xMin = Math.Max(pxMin, qxMin);
            double xMax = Math.Min(pxMax, pxMax);
            if (xMax - xMin <= Tolerance)
            {
                double yMin = Math.Max(pyMin, qyMin);
                double yMax = Math.Min(pyMax, qyMax);
                if (yMax - yMin <= Tolerance)
                {
                    accepted = true;
                    intersect = new Point(0.5 * (xMin + xMax), 0.5 * (yMin + yMax));
                }
            }
        }
        return new List<object> { accepted, excluded, intersect };
    }

    private static bool RectanglesOverlap(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;
    }

    private static 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);
    }

    // de Casteljau's algorithm
    private static void SubdivideQuadSpline(QuadSpline q, double t, QuadSpline u, QuadSpline v)
    {
        double s = 1.0 - t;
        u.C0 = q.C0;
        v.C2 = q.C2;
        u.C1 = s * q.C0 + t * q.C1;
        v.C1 = s * q.C1 + t * q.C2;
        u.C2 = s * u.C1 + t * v.C1;
        v.C0 = u.C2;
    }

    public struct Point
    {
        public double X { get; }
        public double Y { get; }

        public Point(double x, double y)
        {
            X = x;
            Y = y;
        }
    }

    public class QuadSpline
    {
        public double C0 { get; set; }
        public double C1 { get; set; }
        public double C2 { get; set; }

        public QuadSpline(double c0, double c1, double c2)
        {
            C0 = c0;
            C1 = c1;
            C2 = c2;
        }

        public QuadSpline() : this(0.0, 0.0, 0.0) { }
    }

    public class QuadCurve
    {
        public QuadSpline X { get; set; }
        public QuadSpline Y { get; set; }

        public QuadCurve(QuadSpline x, QuadSpline y)
        {
            X = x;
            Y = y;
        }

        public QuadCurve() : this(new QuadSpline(), new QuadSpline()) { }
    }

    private const double Tolerance = 0.000_000_1;
    private const double Spacing = 10 * Tolerance;
}
Output:
The points of intersection are:
(  0.654983,  2.854983 )
( -0.681025,  2.681025 )
(  0.881025,  1.118975 )
( -0.854983,  1.345017 )


C++

#include <cmath>
#include <iomanip>
#include <iostream>
#include <stack>
#include <tuple>
#include <vector>

const double TOLERANCE = 0.000'000'1;
const double SPACING = 10 * TOLERANCE;

typedef std::pair<double, double> point;

class quad_spline {
public:
	quad_spline(double c0, double c1, double c2) : c0(c0), c1(c1), c2(c2) {};
	quad_spline() : c0(0.0), c1(0.0), c2(0.0) {};
	double c0, c1, c2;
};

class quad_curve {
public:
	quad_curve(quad_spline x, quad_spline y) : x(x), y(y) {};
	quad_curve() : x(quad_spline()), y(quad_spline()) {};
	quad_spline x, y;
};

// de Casteljau's algorithm
void subdivide_quad_spline(const quad_spline& q, const double& t, quad_spline& u, quad_spline& v) {
	const double s = 1.0 - t;
	u.c0 = q.c0;
	v.c2 = q.c2;
	u.c1 = s * q.c0 + t * q.c1;
	v.c1 = s * q.c1 + t * q.c2;
	u.c2 = s * u.c1 + t * v.c1;
	v.c0 = u.c2;
}

void subdivide_quad_curve(const quad_curve& q, const double t, quad_curve& u, quad_curve& v) {
	subdivide_quad_spline(q.x, t, u.x, v.x);
	subdivide_quad_spline(q.y, t, u.y, v.y);
}

bool rectangles_overlap(const double& xa0, const double& ya0, const double& xa1, const double& ya1,
					    const double& xb0, const double& yb0, const double& xb1, const double& yb1) {
	return xb0 <= xa1 && xa0 <= xb1 && yb0 <= ya1 && ya0 <= yb1;
}

std::tuple<bool, bool, point> test_intersection(const quad_curve& p, const quad_curve& q) {
	const double px_min = std::min(std::min(p.x.c0, p.x.c1), p.x.c2);
	const double py_min = std::min(std::min(p.y.c0, p.y.c1), p.y.c2);
	const double px_max = std::max(std::max(p.x.c0, p.x.c1), p.x.c2);
	const double py_max = std::max(std::max(p.y.c0, p.y.c1), p.y.c2);

	const double qx_min = std::min(std::min(q.x.c0, q.x.c1), q.x.c2);
	const double qy_min = std::min(std::min(q.y.c0, q.y.c1), q.y.c2);
	const double qx_max = std::max(std::max(q.x.c0, q.x.c1), q.x.c2);
	const double qy_max = std::max(std::max(q.y.c0, q.y.c1), q.y.c2);

	bool accepted = false;
	bool excluded = true;
	point intersect = std::make_pair(0.0, 0.0);

	if ( rectangles_overlap(px_min, py_min, px_max, py_max, qx_min, qy_min, qx_max, qy_max) ) {
		excluded = false;
		const double x_min = std::max(px_min, qx_min);
		const double x_max = std::min(px_max, px_max);
		if ( x_max - x_min <= TOLERANCE ) {
			const double y_min = std::max(py_min, qy_min);
			const double y_max = std::min(py_max, qy_max);
			if ( y_max - y_min <= TOLERANCE ) {
				accepted = true;
				intersect = std::make_pair(0.5 * ( x_min + x_max ), 0.5 * ( y_min + y_max));
			}
		}
	}
	return std::make_tuple(accepted, excluded, intersect);
}

bool seems_to_be_duplicate(const std::vector<point>& intersects, const point& pt) {
	for ( point intersect : intersects ) {
		if ( fabs(intersect.first - pt.first) < SPACING && fabs(intersect.second - pt.second) < SPACING ) {
			return true;
		}
	}
	return false;
}

std::vector<point> find_intersects(const quad_curve& p, const quad_curve& q) {
	std::vector<point> result;
	std::stack<quad_curve> stack;
	stack.push(p);
	stack.push(q);

	while ( ! stack.empty() ) {
		quad_curve pp = stack.top();
		stack.pop();
		quad_curve qq = stack.top();
		stack.pop();

		std::tuple<bool, bool, point> objects = test_intersection(pp, qq);
		bool accepted, excluded;
		point intersect;
		std::tie(accepted, excluded, intersect) = objects;

		if ( accepted ) {
			if ( ! seems_to_be_duplicate(result, intersect) ) {
				result.push_back(intersect);
			}
		} else if ( ! excluded ) {
			quad_curve p0{}, q0{}, p1{}, q1{};
			subdivide_quad_curve(pp, 0.5, p0, p1);
			subdivide_quad_curve(qq, 0.5, q0, q1);
			for ( quad_curve item : { p0, q0, p0, q1, p1, q0, p1, q1 } ) {
				stack.push(item);
			}
		}
	}
	return result;
}

int main() {
	quad_curve vertical(quad_spline(-1.0, 0.0, 1.0), quad_spline(0.0, 10.0, 0.0));
	// QuadCurve vertical represents the Bezier curve having control points (-1, 0), (0, 10) and (1, 0)
	quad_curve horizontal(quad_spline(2.0, -8.0, 2.0), quad_spline(1.0, 2.0, 3.0));
	// QuadCurve horizontal represents the Bezier curve having control points (2, 1), (-8, 2) and (2, 3)

	std::cout << "The points of intersection are:" << std::endl;
	std::vector<point> intersects = find_intersects(vertical, horizontal);
	for ( const point& intersect : intersects ) {
		std::cout << "( " << std::setw(9) << std::setprecision(6) << intersect.first
				  << ", " << intersect.second << " )" << std::endl;
	}
}
Output:
The points of intersection are:
(  0.654983, 2.85498 )
( -0.681025, 2.68102 )
(  0.881025, 1.11898 )
( -0.854983, 1.34502 )

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)

FreeBASIC

FreeBASIC Implementation 1

Translation of: C
' 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.

#define Min(a, b) iif((a) < (b), (a), (b))
#define Max(a, b) iif((a) > (b), (a), (b))

' Note that these are all mutable as we need to pass by reference.
Type Punto
    x As Double
    y As Double
End Type

Type QuadSpline
    c0 As Double
    c1 As Double
    c2 As Double    ' Non-parametric spline
End Type

Type QuadCurve
    x As QuadSpline
    y As QuadSpline     ' Planar parametric spline
End Type

' Subdivision by de Casteljau's algorithm
Sub subdivideQuadSpline(q As QuadSpline, t As Double, u As QuadSpline, v As QuadSpline)
    Dim As Double s = 1.0 - t
    Dim As Double c0 = q.c0
    Dim As Double c1 = q.c1
    Dim As 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
End Sub

Sub subdivideQuadCurve(q As QuadCurve, t As Double, u As QuadCurve, v As QuadCurve)
    subdivideQuadSpline(q.x, t, u.x, v.x)
    subdivideQuadSpline(q.y, t, u.y, v.y)
End Sub

' It is assumed that xa0 <= xa1, ya0 <= ya1, xb0 <= xb1, and yb0 <= yb1.
Function rectsOverlap(xa0 As Double, ya0 As Double, xa1 As Double, ya1 As Double, _
    xb0 As Double, yb0 As Double, xb1 As Double, yb1 As Double) As Boolean
    Return (xb0 <= xa1 And xa0 <= xb1 And yb0 <= ya1 And ya0 <= yb1)
End Function

Function max3(x As Double, y As Double, z As Double) As Double
    Return Max(Max(x, y), z)
End Function

Function min3(x As Double, y As Double, z As Double) As Double
    Return Min(Min(x, y), z)
End Function

' This accepts the point as an intersection if the boxes are small enough.
Sub testIntersect(p As QuadCurve, q As QuadCurve, tol As Double, _
    Byref exclude As Boolean, Byref accept As Boolean, Byref intersect As Punto)
    Dim As Double pxmin = min3(p.x.c0, p.x.c1, p.x.c2)
    Dim As Double pymin = min3(p.y.c0, p.y.c1, p.y.c2)
    Dim As Double pxmax = max3(p.x.c0, p.x.c1, p.x.c2)
    Dim As Double pymax = max3(p.y.c0, p.y.c1, p.y.c2)
    
    Dim As Double qxmin = min3(q.x.c0, q.x.c1, q.x.c2)
    Dim As Double qymin = min3(q.y.c0, q.y.c1, q.y.c2)
    Dim As Double qxmax = max3(q.x.c0, q.x.c1, q.x.c2)
    Dim As 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) Then
        exclude = False
        Dim As Double xmin = Max(pxmin, qxmin)
        Dim As Double xmax = Min(pxmax, qxmax)
        Assert(xmax >= xmin)
        If xmax - xmin <= tol Then
            Dim As Double ymin = Max(pymin, qymin)
            Dim As Double ymax = Min(pymax, qymax)
            Assert(ymax >= ymin)
            If ymax - ymin <= tol Then
                accept = True
                intersect.x = 0.5 * xmin + 0.5 * xmax
                intersect.y = 0.5 * ymin + 0.5 * ymax
            End If
        End If
    End If
End Sub

Function seemsToBeDuplicate(intersects() As Punto, icount As Integer, _
    xy As Punto, spacing As Double) As Boolean
    Dim As Boolean seemsToBeDup = False
    Dim As Integer i = 0
    While Not seemsToBeDup And i <> icount
        Dim As Punto pt = intersects(i)
        seemsToBeDup = Abs(pt.x - xy.x) < spacing And Abs(pt.y - xy.y) < spacing
        i += 1
    Wend
    Return seemsToBeDup
End Function

Sub findIntersects(p As QuadCurve, q As QuadCurve, tol As Double, _
    spacing As Double, intersects() As Punto)
    Dim As Integer numIntersects = 0
    Type workset
        p As QuadCurve
        q As QuadCurve
    End Type
    Dim As workset workload(64)
    Dim As Integer numWorksets = 1
    workload(0) = Type<workset>(p, q)
    
    ' Quit looking after having emptied the workload.
    While numWorksets <> 0
        Dim As workset work = workload(numWorksets-1)
        numWorksets -= 1
        Dim As Boolean exclude, accept
        Dim As Punto intersect
        testIntersect(work.p, work.q, tol, exclude, accept, intersect)
        If accept Then
            ' To avoid detecting the same intersection twice, require some
            ' space between intersections.
            If Not seemsToBeDuplicate(intersects(), numIntersects, intersect, spacing) Then
                intersects(numIntersects) = intersect
                numIntersects += 1
            End If
        Elseif Not exclude Then
            Dim As QuadCurve p0, p1, q0, q1
            subdivideQuadCurve(work.p, 0.5, p0, p1)
            subdivideQuadCurve(work.q, 0.5, q0, q1)
            workload(numWorksets) = Type<workset>(p0, q0)
            numWorksets += 1
            workload(numWorksets) = Type<workset>(p0, q1)
            numWorksets += 1
            workload(numWorksets) = Type<workset>(p1, q0)
            numWorksets += 1
            workload(numWorksets) = Type<workset>(p1, q1)
            numWorksets += 1
        End If
    Wend
End Sub

Dim As QuadCurve p, q
p.x = Type<QuadSpline>(-1.0,  0.0, 1.0)
p.y = Type<QuadSpline>( 0.0, 10.0, 0.0)
q.x = Type<QuadSpline>( 2.0, -8.0, 2.0)
q.y = Type<QuadSpline>( 1.0,  2.0, 3.0)
Dim As Double tol = 0.0000001
Dim As Double spacing = tol * 10.0
Dim As Punto intersects(4)
findIntersects(p, q, tol, spacing, intersects())
For i As Integer = 0 To 3
    Print "("; intersects(i).x; ", "; intersects(i).y; ")"
Next i

Sleep
Output:
( 0.6549834311008453,  2.854983419179916)
( 0.8810249865055084,  1.118975013494492)
(-0.6810249347860431,  2.681024998426437)
(-0.8549834191799164,  1.345016568899155)

FreeBASIC Implementation 2

Translation of: C
' 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. */

Type bernstein_spline
    b0 As Double
    b1 As Double
    b2 As Double
End Type

Type spower_spline
    c0 As Double
    c1 As Double
    c2 As Double
End Type

Type spower_curve
    x As spower_spline
    y As spower_spline
End Type

' Convert a non-parametric spline from Bernstein basis to s-power.
Function bernstein_spline_to_spower(S As bernstein_spline) As spower_spline
    Dim As spower_spline T
    T.c0 = S.b0
    T.c1 = (2 * S.b1) - S.b0 - S.b2
    T.c2 = S.b2
    Return T
End Function

' 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.
Function spower_spline_portion(S As spower_spline, t0 As Double, t1 As Double) As spower_spline
    Dim As Double t0_t0 = t0 * t0
    Dim As Double t0_t1 = t0 * t1
    Dim As Double t1_t1 = t1 * t1
    Dim As Double c2p1m0 = S.c2 + S.c1 - S.c0
    
    Dim As spower_spline T
    T.c0 = S.c0 + (c2p1m0 * t0) - (S.c1 * t0_t0)
    T.c1 = (S.c1 * t1_t1) - (2 * S.c1 * t0_t1) + (S.c1 * t0_t0)
    T.c2 = S.c0 + (c2p1m0 * t1) - (S.c1 * t1_t1)
    Return T
End Function

Function spower_curve_portion(C As spower_curve, t0 As Double, t1 As Double) As spower_curve
    Dim As spower_curve D
    D.x = spower_spline_portion(C.x, t0, t1)
    D.y = spower_spline_portion(C.y, t0, t1)
    Return D
End Function

' Given a parametric curve, is it "flat enough" to have its quadratic
' terms removed?
Function flat_enough(C As spower_curve, tol As Double) As Boolean
  ' 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:
    Dim As Double cx0 = C.x.c0
    Dim As Double cx1 = C.x.c1
    Dim As Double cx2 = C.x.c2
    Dim As Double cy0 = C.y.c0
    Dim As Double cy1 = C.y.c1
    Dim As Double cy2 = C.y.c2
    Dim As Double dx = cx2 - cx0
    Dim As Double dy = cy2 - cy0
    Dim As Double error_squared = 0.125 * ((cx1 * cx1) + (cy1 * cy1))
    Dim As Double length_squared = (dx * dx) + (dy * dy)
    Return (error_squared <= length_squared * tol * tol)
End Function

' 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.
Sub test_line_segment_intersection(ax0 As Double, ax1 As Double, _
    ay0 As Double, ay1 As Double, bx0 As Double, bx1 As Double, _
    by0 As Double, by1 As Double, Byref they_intersect As Boolean, _
    Byref x As Double, Byref y As Double)
    Dim As Double anumer = ((bx1 - bx0) * ay0 - (by1 - by0) * ax0 _
    + bx0 * by1 - bx1 * by0)
    Dim As Double bnumer = -((ax1 - ax0) * by0 - (ay1 - ay0) * bx0 _
    + ax0 * ay1 - ax1 * ay0)
    Dim As Double denom = ((ax1 - ax0) * (by1 - by0) _
    - (ay1 - ay0) * (bx1 - bx0))
    Dim As Double ta = anumer / denom   ' Parameter of segment a.
    Dim As Double tb = bnumer / denom   ' Parameter of segment b.
    they_intersect = (0 <= ta And ta <= 1 And 0 <= tb And tb <= 1)
    If they_intersect Then
        x = ((1 - ta) * ax0) + (ta * ax1)
        y = ((1 - ta) * ay0) + (ta * ay1)
    End If
End Sub

Function too_close(x As Double, y As Double, xs() As Double, ys() As Double, _
    num_points As Integer, spacing As Double) As Boolean
    Dim As Boolean tooclose = False
    Dim As Integer i = 0
    While Not tooclose And i <> num_points
        tooclose = (Abs(x - xs(i)) < spacing And Abs(y - ys(i)) < spacing)
        i += 1
    Wend
    Return tooclose
End Function

Sub recursion(tp0 As Double, tp1 As Double, tq0 As Double, tq1 As Double, _
    P As spower_curve, Q As spower_curve, tol As Double, spacing As Double, _
    max_points As Integer, xs() As Double, ys() As Double, Byref num_points As Integer)
    If num_points = max_points Then
        ' do nothing
    Elseif Not flat_enough(spower_curve_portion(P, tp0, tp1), tol) Then
        Dim As Double tp_half = (0.5 * tp0) + (0.5 * tp1)
        If Not flat_enough(spower_curve_portion(Q, tq0, tq1), tol) Then
            Dim As 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)
        End If
    Elseif Not flat_enough(spower_curve_portion(Q, tq0, tq1), tol) Then
        Dim As 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
        Dim As spower_curve P1 = spower_curve_portion(P, tp0, tp1)
        Dim As spower_curve Q1 = spower_curve_portion(Q, tq0, tq1)
        Dim As Boolean they_intersect
        Dim As 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 And Not too_close(x, y, xs(), ys(), num_points, spacing) Then
        xs(num_points) = x
        ys(num_points) = y
        num_points += 1
    End If
End If
End Sub

Sub find_intersections(P As spower_curve, Q As spower_curve, _
    flatness_tolerance As Double, point_spacing As Double, _
    max_points As Integer, xs() As Double, ys() As Double, _
    Byref num_points As Integer)
    num_points = 0
    recursion(0, 1, 0, 1, P, Q, flatness_tolerance, point_spacing, _
    max_points, xs(), ys(), num_points)
End Sub

Dim As bernstein_spline bPx = Type<bernstein_spline>(-1, 0, 1)
Dim As bernstein_spline bPy = Type<bernstein_spline>(0, 10, 0)
Dim As bernstein_spline bQx = Type<bernstein_spline>(2, -8, 2)
Dim As bernstein_spline bQy = Type<bernstein_spline>(1, 2, 3)

Dim As spower_spline Px = bernstein_spline_to_spower(bPx)
Dim As spower_spline Py = bernstein_spline_to_spower(bPy)
Dim As spower_spline Qx = bernstein_spline_to_spower(bQx)
Dim As spower_spline Qy = bernstein_spline_to_spower(bQy)

Dim As spower_curve P = Type<spower_curve>(Px, Py)
Dim As spower_curve Q = Type<spower_curve>(Qx, Qy)

Dim As Double flatness_tolerance = 0.001
Dim As Double point_spacing = 0.000001 ' Max norm minimum spacing.

Const max_points As Integer = 10
Dim As Double xs(max_points)
Dim As Double ys(max_points)
Dim As Integer num_points

find_intersections(P, Q, flatness_tolerance, point_spacing, _
max_points, xs(), ys(), num_points)

For i As Integer = 0 To num_points - 1
    Print "("; xs(i); ", "; ys(i); ")"
Next i

Sleep

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)

Java

import java.util.ArrayList;
import java.util.List;
import java.util.Stack;

public final class BezierCurveIntersection {

	public static void main(String[] aArgs) {
		QuadCurve vertical = new QuadCurve( new QuadSpline(-1.0, 0.0, 1.0), new QuadSpline(0.0, 10.0, 0.0) );
		// QuadCurve vertical represents the Bezier curve having control points (-1, 0), (0, 10) and (1, 0) 
		QuadCurve horizontal = new QuadCurve( new QuadSpline(2.0, -8.0, 2.0), new QuadSpline(1.0, 2.0, 3.0) );
		// QuadCurve horizontal represents the Bezier curve having control points (2, 1), (-8, 2) and (2, 3) 		

		System.out.println("The points of intersection are:");
		List<Point> intersects = findIntersects(vertical, horizontal);
		for ( Point intersect : intersects ) {
			System.out.println(String.format("%s%9.6f%s%9.6f%s", "( ", intersect.aX, ", ", intersect.aY, " )"));
		}
	}
	
	private static List<Point> findIntersects(QuadCurve aP, QuadCurve aQ) {
		List<Point> result = new ArrayList<Point>();
		Stack<QuadCurve> stack = new Stack<QuadCurve>();
		stack.push(aP);
		stack.push(aQ);
		
		while ( ! stack.isEmpty() ) {
			QuadCurve pp = stack.pop();
		    QuadCurve qq = stack.pop();
		    List<Object> objects = testIntersection(pp, qq);
		    final boolean accepted = (boolean) objects.get(0);
		    final boolean excluded = (boolean) objects.get(1);
		    Point intersect = (Point) objects.get(2);
		    
		    if ( accepted ) {
		    	if ( ! seemsToBeDuplicate(result, intersect) ) {
		    		result.add(intersect);
		    	}
		    } else if ( ! excluded ) {
		    	QuadCurve p0 = new QuadCurve();
		    	QuadCurve q0 = new QuadCurve();
		    	QuadCurve p1 = new QuadCurve();
		    	QuadCurve q1 = new QuadCurve();
		    	subdivideQuadCurve(pp, 0.5, p0, p1);
		    	subdivideQuadCurve(qq, 0.5, q0, q1);
		    	stack.addAll(List.of( p0, q0, p0, q1, p1, q0, p1, q1 ));
		    }
		}
		return result;
	}
	
	private static boolean seemsToBeDuplicate(List<Point> aIntersects, Point aPoint) {
		for ( Point intersect : aIntersects ) {
			if ( Math.abs(intersect.aX - aPoint.aX) < SPACING && Math.abs(intersect.aY - aPoint.aY) < SPACING ) {
				return true;
			}
		}
		return false;
	}
	
	private static List<Object> testIntersection(QuadCurve aP, QuadCurve aQ) {
	    final double pxMin = Math.min(Math.min(aP.x.c0, aP.x.c1), aP.x.c2);
	    final double pyMin = Math.min(Math.min(aP.y.c0, aP.y.c1), aP.y.c2);
	    final double pxMax = Math.max(Math.max(aP.x.c0, aP.x.c1), aP.x.c2);
	    final double pyMax = Math.max(Math.max(aP.y.c0, aP.y.c1), aP.y.c2);

	    final double qxMin = Math.min(Math.min(aQ.x.c0, aQ.x.c1), aQ.x.c2);
	    final double qyMin = Math.min(Math.min(aQ.y.c0, aQ.y.c1), aQ.y.c2);
	    final double qxMax = Math.max(Math.max(aQ.x.c0, aQ.x.c1), aQ.x.c2);
	    final double qyMax = Math.max(Math.max(aQ.y.c0, aQ.y.c1), aQ.y.c2);
	    
	    boolean accepted = false;
	    boolean excluded = true;
	    Point intersect = new Point(0.0, 0.0);
	    
	    if ( rectanglesOverlap(pxMin, pyMin, pxMax, pyMax, qxMin, qyMin, qxMax, qyMax) ) {
		    excluded = false;
		    final double xMin = Math.max(pxMin, qxMin);
		    final double xMax = Math.min(pxMax, pxMax);
		    if ( xMax - xMin <= TOLERANCE ) {
		    	final double yMin = Math.max(pyMin, qyMin);
		    	final double yMax = Math.min(pyMax, qyMax);
		    	if ( yMax - yMin <= TOLERANCE ) {
		    		accepted = true;
		    		intersect = new Point(0.5 * ( xMin + xMax), 0.5 * ( yMin +  yMax));
		    	}
		    }
	    }
	    return List.of( accepted, excluded, intersect );
	}
	
	private static boolean rectanglesOverlap(double aXa0, double aYa0, double aXa1, double aYa1,
										     double aXb0, double aYb0, double aXb1, double aYb1) {
		return aXb0 <= aXa1 && aXa0 <= aXb1 && aYb0 <= aYa1 && aYa0 <= aYb1;
	}
	
	private static void subdivideQuadCurve(QuadCurve aQ, double aT, QuadCurve aU, QuadCurve aV) {
		subdivideQuadSpline(aQ.x, aT, aU.x, aV.x);
		subdivideQuadSpline(aQ.y, aT, aU.y, aV.y);
	}
	
	// de Casteljau's algorithm
	private static void subdivideQuadSpline(QuadSpline aQ, double aT, QuadSpline aU, QuadSpline aV) {
		final double s = 1.0 - aT;
		aU.c0 = aQ.c0;
		aV.c2 = aQ.c2;
		aU.c1 = s * aQ.c0 + aT * aQ.c1;
		aV.c1 = s * aQ.c1 + aT * aQ.c2;
		aU.c2 = s * aU.c1 + aT * aV.c1;
		aV.c0 = aU.c2;
	}
	
	private static record Point(double aX, double aY) {}
	
	private static class QuadSpline {
		
		public QuadSpline(double aC0, double aC1, double aC2) {
			c0 = aC0; c1 = aC1; c2 = aC2;
		}
		
		public QuadSpline() {
			this(0.0, 0.0, 0.0);
		}
		
		private double c0, c1, c2;	
		
	}	
	
	private static class QuadCurve {
		
		public QuadCurve(QuadSpline aX, QuadSpline aY) {
			x = aX; y = aY;
		}
		
		public QuadCurve() {
			this( new QuadSpline(), new QuadSpline() );
		}
		
		private QuadSpline x, y;
		
	}
	
	private static final double TOLERANCE = 0.000_000_1;
	private static final double SPACING = 10 * TOLERANCE;

}
Output:
The points of intersection are:
(  0.654983,  2.854983 )
( -0.681025,  2.681025 )
(  0.881025,  1.118975 )
( -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:

# This program combines the methods of the 2nd C implementation (which
# by itself is inefficient) with those of the Modula-2 implementation,
# and then rearranges the computations to try to achieve greater
# efficiency.
#
# The algorithm actually returns t-parameters for two curves, as a
# pair for each intersection point. This is exactly what one might
# want: for instance, to break a font glyph, made from two or more
# other glyphs, into pieces at the points of intersection of all the
# outlines.
#
# The code below is written to illustrate the algorithm rather than to
# squeeze performance out of Icon. For instance, I use a "set" to
# store the workload, and, when choosing the next workitem-pair to
# work on, do so by random selection. It would be faster, certainly,
# to use an Icon "list", as either a stack or a queue or both.
#
# It is also possible, of course, to cast the algorithm as a recursive
# procedure.
#
#
# 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.
#

record point (x, y)
record spower (c0, c1, c2)
record curve (x, y)
record workitem (P, t0, t1, pt0, pt1)

$define P_controls [point (-1, 0), point ( 0, 10), point ( 1,  0)]
$define Q_controls [point ( 2, 1), point (-8,  2), point ( 2,  3)]

$define DEFAULT_NUMPIECES 2     # Bisection.

# Tolerance of the ratio of a bound on the non-linear component to the
# length of the segment. I use a max norm but you can use your
# favorite norm.
$define DEFAULT_FLATNESS_TOLERANCE 0.0001

# For demonstration, I choose a minimum spacing between intersection
# points equal to several times single precision machine epsilon. I
# measure distance using a max norm, but you can use your favorite
# norm.
$define DEFAULT_MINIMUM_SPACING 1e-6

procedure main ()
  local P, Q
  local intersections, xsect

  P := controls_to_curve ! P_controls
  Q := controls_to_curve ! Q_controls

  intersections := find_intersections (P, Q)

  write ()
  write_tabbed_line ("          convex up\t" ||
                     "          convex left\t")
  every xsect := !intersections do
    write_tabbed_line (xsect[1] || "\t(" ||
                       xsect[2].x || ", " || xsect[2].y || ")\t" ||
                       xsect[3] || "\t(" ||
                       xsect[4].x || ", " || xsect[4].y || ")")
  write ()
end

procedure write_tabbed_line (line)
  write (detab (line, 18, 56, 74))
end

procedure find_intersections (P, Q, tol, spacing)
  # Return a list of [tP, pointP, tQ, pointQ] for the intersections,
  # sorted by tP values.

  local workload, ts
  local tP, ptP
  local tQ, ptQ
  local tbl, intersections

  /tol := DEFAULT_FLATNESS_TOLERANCE
  /spacing := DEFAULT_MINIMUM_SPACING

  workload := initial_workload (P, Q)

  tbl := table ()
  every ts := process_workload (tol, workload) do
  {
    tP := ts[1];  ptP := curve_eval (P, tP)
    tQ := ts[2];  ptQ := curve_eval (Q, tQ)
    not (max (abs ((!tbl)[2].x - ptP.x),
              abs ((!tbl)[2].y - ptP.y)) < spacing) &
      not (max (abs ((!tbl)[4].x - ptQ.x),
                abs ((!tbl)[4].y - ptQ.y)) < spacing) &
      tbl[tP] := [tP, ptP, tQ, ptQ]
  }
  tbl := sort (tbl, 1)
  every put (intersections := [], (!tbl)[2])
  return intersections
end

procedure process_workload (tol, workload)
  # Generate pairs of t-parameters.

  local pair, ts

  while *workload ~= 0 do
  {
    pair := ?workload
    delete (workload, pair)
    if rectangles_overlap (pair[1].pt0, pair[1].pt1,
                           pair[2].pt0, pair[2].pt1) then
    {
      if flat_enough (tol, pair[1]) then
      {
        if flat_enough (tol, pair[2]) then
        {
          if ts := segment_parameters (pair[1].pt0, pair[1].pt1,
                                       pair[2].pt0, pair[2].pt1) then
            suspend [(1 - ts[1]) * pair[1].t0 + ts[1] * pair[1].t1,
                     (1 - ts[2]) * pair[2].t0 + ts[2] * pair[2].t1]
        }
        else
          every insert (workload, [pair[1],
                                   split_workitem (pair[2])])
      }
      else
      {
        if flat_enough (tol, pair[2]) then
          every insert (workload, [split_workitem (pair[1]),
                                   pair[2]])
        else
          every insert (workload, [split_workitem (pair[1]),
                                   split_workitem (pair[2])])
      }
    }
  }
end

procedure split_workitem (W, num_pieces)
  # Split a workitem in pieces and generate the pieces.

  local fraction, t1_t0, ts, pts, i

  /num_pieces := DEFAULT_NUMPIECES

  fraction := 1.0 / num_pieces
  t1_t0 := W.t1 - W.t0

  every put (ts := [],
             W.t0 + (1 to num_pieces - 1) * fraction * t1_t0)
  every put (pts := [], curve_eval (W.P, !ts))
  ts := [W.t0] ||| ts ||| [W.t1]
  pts := [W.pt0] ||| pts ||| [W.pt1]

  every i := 1 to *pts - 1 do
    suspend (workitem (W.P, ts[i], ts[i + 1], pts[i], pts[i + 1]))
end

procedure initial_workload (P, Q)
  # Create an initial workload set, by breaking P and Q at any
  # critical points.

  local workload

  every insert (workload := set (), [break_at_critical_points (P),
                                     break_at_critical_points (Q)])
  return workload
end

procedure break_at_critical_points (P)
  # Generate workitems for the curve P, after breaking it at any
  # critical points.

  local ts, pts, i

  ts := [0] ||| sort (curve_critical_points (P)) ||| [1]
  every put (pts := [], curve_eval (P, !ts))
  every i := 1 to *pts - 1 do
    suspend (workitem (P, ts[i], ts[i + 1], pts[i], pts[i + 1]))
end

procedure flat_enough (tol, P, t0, t1, pt0, pt1)
  # Is the [t0,t1] portion of the curve P flat enough to be treated as
  # a straight line between pt0 and pt1, where pt0 and pt1 are the
  # endpoints of the portion?

  local error, length

  # Let flat_enough be called this way, where W is a workitem:
  # flat_enough(tol,W)
  if /t0 then
  {
    pt1 := P.pt1
    pt0 := P.pt0
    t1 := P.t1
    t0 := P.t0
    P := P.P
  }

  # pt0 and pt1 probably have been computed before and saved, but if
  # necessary they could be computed now:
  /pt0 := curve_eval (P, t0)
  /pt1 := curve_eval (P, t1)

  # 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, which
  # uses "max norm" length measurements. (Substitute your favorite
  # norm.)
  error := 0.25 * max (abs (spower_center_coef (P.x, t0, t1)),
                       abs (spower_center_coef (P.y, t0, t1)))
  length := max (abs (pt1.x - pt0.x), abs (pt1.y - pt0.y))
  ((error <= length * tol) & return) | fail
end

procedure curve_eval (P, t)
  # Return the point that lies on the curve P at parameter value t.
  return point (spower_eval (P.x, t), spower_eval (P.y, t))
end

procedure curve_critical_points (P)
  # Return a set containing parameter values (values of t) for the
  # critical points of curve P.

  local ts

  ts := set ()
  insert (ts, spower_critical_point (P.x))
  insert (ts, spower_critical_point (P.y))
  return ts
end

procedure spower_eval (p, t)
  # Evaluate the s-power spline p at t.
  return (p.c0 + (p.c1 * t)) * (1 - t) + (p.c2 * t)
end

procedure spower_center_coef (p, t0, t1)
  # Return the center coefficient for the [t0,t1] portion of the
  # s-power spline p.
  if /t1 then { t1 := t0[2]; t0 := t0[1] } # Allow a list as t0.
  return p.c1 * ((t1 - t0 - t0) * t1 + (t0 * t0))
end

procedure spower_critical_point (p)
  # Return t in (0,1) where p is at a critical point, else fail.

  local t

  p.c1 = 0 & fail               # The spline is linear
  p.c0 = p.c2 & return 0.5      # The spline is "pulse-like".

  t := (0.5 * (p.c2 + p.c1 - p.c0)) / p.c1 # Root of the derivative.
  0 < t < 1 & return t                    
  fail
end

procedure rectangles_overlap (ptA0, ptA1, ptB0, ptB1)
  # Do the rectangles with corners at (ptA0,ptA1) and (ptB0,ptB1)
  # overlap?
  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

procedure segment_parameters (ptA0, ptA1, ptB0, ptB1)
  # 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.

  local ax0, ax1, ay0, ay1
  local bx0, bx1, by0, by1
  local ax1_ax0, ay1_ay0
  local bx1_bx0, by1_by0
  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

  ax1_ax0 := ax1 - ax0
  ay1_ay0 := ay1 - ay0
  bx1_bx0 := bx1 - bx0
  by1_by0 := by1 - by0

  denom := (ax1_ax0 * by1_by0) - (ay1_ay0 * bx1_bx0)

  anumer :=
    (bx1_bx0 * ay0) - (by1_by0 * ax0) + (bx0 * by1) - (bx1 * by0)
  tA := anumer / denom
  0 <= tA <= 1 | fail

  bnumer :=
    -((ax1_ax0 * by0) - (ay1_ay0 * bx0) + (ax0 * ay1) - (ax1 * ay0))
  tB := bnumer / denom
  0 <= tB <= 1 | fail

  return [tA, tB]
end

procedure controls_to_curve (p0, p1, p2)
  # Convert control points to a curve in s-power basis.
  return curve (spower (p0.x, (2 * p1.x) - p0.x - p2.x, p2.x),
                spower (p0.y, (2 * p1.y) - p0.y - p2.y, p2.y))
end

procedure abs (x)
  return (if x < 0 then -x else x)
end

procedure max (x, y)
  return (if x < y then y else x)
end
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)

Julia

Translation of: 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.
=#

mutable struct Point
    x::Float64
    y::Float64
    Point() = new(0, 0)
end

mutable struct QuadSpline # Non-parametric spline.
    c0::Float64
    c1::Float64
    c2::Float64
    QuadSpline(a = 0.0, b = 0.0, c = 0.0) = new(a, b, c)
end

mutable struct QuadCurve # Planar parametric spline.
    x::QuadSpline
    y::QuadSpline
    QuadCurve(x = QuadSpline(), y = QuadSpline()) = new(x, y)
end

const Workset = Tuple{QuadCurve, QuadCurve}

""" Subdivision by de Casteljau's algorithm. """
function subdivide!(q::QuadSpline, t::Float64, u::QuadSpline, 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
end

function subdivide!(q::QuadCurve, t::Float64, u::QuadCurve, v::QuadCurve)
    subdivide!(q.x, t, u.x, v.x) 
    subdivide!(q.y, t, u.y, v.y)
end

""" It is assumed that xa0 <= xa1, ya0 <= ya1, xb0 <= xb1, and yb0 <= yb1. """
function rectsoverlap(xa0, ya0, xa1, ya1, xb0, yb0, xb1, yb1)
    return xb0 <= xa1 && xa0 <= xb1 && yb0 <= ya1 && ya0 <= yb1
end

""" This accepts the point as an intersection if the boxes are small enough. """
function testintersect(p::QuadCurve, q::QuadCurve, tol::Float64, intersect::Point)
    pxmin = min(p.x.c0, p.x.c1, p.x.c2)
    pymin = min(p.y.c0, p.y.c1, p.y.c2)
    pxmax = max(p.x.c0, p.x.c1, p.x.c2)
    pymax = max(p.y.c0, p.y.c1, p.y.c2)

    qxmin = min(q.x.c0, q.x.c1, q.x.c2)
    qymin = min(q.y.c0, q.y.c1, q.y.c2)
    qxmax = max(q.x.c0, q.x.c1, q.x.c2)
    qymax = max(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 = max(pxmin, qxmin)
        xmax = min(pxmax, pxmax)
        @assert xmin < xmax "Assertion failure: $xmax < $xmin"        
        if xmax-xmin <= tol
            ymin = max(pymin, qymin)
            ymax = min(pymax, qymax)
            @assert ymin < ymax "Assertion failure: $ymax < $ymin"           
            if ymax-ymin <= tol
                accept = true
                intersect.x = 0.5*xmin + 0.5*xmax
                intersect.y = 0.5*ymin + 0.5*ymax
            end
        end
    end
    return exclude, accept
end

""" return true if there seems to be a duplicate of xy in the intersects `isects` """
seemsduplicate(isects, xy, sp) = any(p -> abs(p.x-xy.x) < sp && abs(p.y-xy.y) < sp, isects)

""" find intersections of p and q """
function findintersects(p, q, tol, spacing)
    intersects = Point[]
    workload = Workset[(p, q)]
    while !isempty(workload)
        work = popfirst!(workload)
        intersect = Point()
        exclude, accept = testintersect(first(work), last(work), tol, intersect)
        if accept
            # To avoid detecting the same intersection twice, require some
            # space between intersections.
            if !seemsduplicate(intersects, intersect, spacing)
                pushfirst!(intersects, intersect)
            end
        elseif !exclude
            p0, p1, q0, q1 = QuadCurve(), QuadCurve(), QuadCurve(), QuadCurve()
            subdivide!(first(work), 0.5, p0, p1)
            subdivide!(last(work), 0.5, q0, q1)
            pushfirst!(workload, (p0, q0), (p0, q1), (p1, q0), (p1, q1))
        end
    end
    return intersects
end

""" test the methods """
function testintersections()
    p, q = QuadCurve(), 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
    for intersect in findintersects(p, q, tol, spacing)
        println("($(intersect.x) $(intersect.y))")
    end
end

testintersections()
Output:
(0.6549834311008453 2.8549834191799164)
(0.8810249865055084 1.1189750134944916)
(-0.6810249984264374 2.681024934786043)
(-0.8549834191799164 1.3450165688991547)

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)

Nim

Translation of: Go
import std/strformat

type
  Point = tuple[x, y: float]
  QuadSpline = tuple[c0, c1, c2: float]   # Non-parametric spline.
  QuadCurve = tuple[x, y: QuadSpline]     # Planar parametric spline.


func subdivideQuadSpline(q: QuadSpline; t: float; u, v: var QuadSpline) =
  ## Subdivision by de Casteljau's algorithm.
  let s = 1.0 - t
  u.c0 = q.c0
  v.c2 = q.c2
  u.c1 = s * q.c0 + t * q.c1
  v.c1 = s * q.c1 + t * q.c2
  u.c2 = s * u.c1 + t * v.c1
  v.c0 = u.c2

func subdivideQuadCurve(q: QuadCurve; t: float; u, v: var 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: float): bool =
  xb0 <= xa1 and xa0 <= xb1 and yb0 <= ya1 and ya0 <= yb1

func max(x, y, z: float): float = max(max(x, y), z)
func min(x, y, z: float): float = min(min(x, y), z)

# This accepts the point as an intersection if the boxes are small enough.
func testIntersect(p, q: QuadCurve; tol: float; exclude, accept: var bool; intersect: var Point) =
  let
    pxmin = min(p.x.c0, p.x.c1, p.x.c2)
    pymin = min(p.y.c0, p.y.c1, p.y.c2)
    pxmax = max(p.x.c0, p.x.c1, p.x.c2)
    pymax = max(p.y.c0, p.y.c1, p.y.c2)

    qxmin = min(q.x.c0, q.x.c1, q.x.c2)
    qymin = min(q.y.c0, q.y.c1, q.y.c2)
    qxmax = max(q.x.c0, q.x.c1, q.x.c2)
    qymax = max(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
    let xmin = max(pxmin, qxmin)
    let xmax = min(pxmax, pxmax)
    assert xmin <= xmax, &"Assertion failure: {xmin} <= {xmax}"
    if xmax - xmin <= tol:
      let ymin = max(pymin, qymin)
      let ymax = min(pymax, qymax)
      assert ymin <= ymax, &"Assertion failure: {ymin} <= {ymax}"
      if ymax - ymin <= tol:
          accept = true
          intersect = (0.5 * xmin + 0.5 * xmax, 0.5 * ymin + 0.5 * ymax)

func seemsToBeDuplicate(intersects: openArray[Point]; xy: Point; spacing: float): bool =
  var i = 0
  while not result and i != intersects.len:
    let pt = intersects[i]
    result = abs(pt.x - xy.x) < spacing and abs(pt.y - xy.y) < spacing
    inc i

func findIntersects(p, q: QuadCurve; tol, spacing: float): seq[Point] =
  var workload = @[(p: p, q: q)]

  # Quit looking after having emptied the workload.
  while workload.len > 0:
    var exclude, accept: bool
    var intersect: Point
    let work = workload.pop()
    testIntersect(work.p, work.q, tol, exclude, accept, intersect)
    if accept:
      # To avoid detecting the same intersection twice, require some
      # space between intersections.
      if not seemsToBeDuplicate(result, intersect, spacing):
          result.add intersect
    elif not exclude:
      var p0, p1, q0, q1: QuadCurve
      subdivideQuadCurve(work.p, 0.5, p0, p1)
      subdivideQuadCurve(work.q, 0.5, q0, q1)
      workload.add @[(p0, q0), (p0, q1), (p1, q0), (p1, q1)]

var p, q: QuadCurve
p.x = (-1.0, 0.0, 1.0)
p.y = (0.0, 10.0, 0.0)
q.x = (2.0, -8.0, 2.0)
q.y = (1.0, 2.0, 3.0)

const Tol = 0.0000001
const Spacing = Tol * 10

let intersects = findIntersects(p, q, Tol, Spacing)
for intersect in intersects:
  echo &"({intersect.x: .6f}, {intersect.y: .6f})"
Output:
( 0.654983,  2.854983)
( 0.881025,  1.118975)
(-0.681025,  2.681025)
(-0.854983,  1.345017)

ObjectIcon

Translation of: Python
Translation of: Icon

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

Translation of: Icon

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

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

A Phix implementation of the "rectangle-pruned curve-flattening algorithm"

Translation of: Pascal

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

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

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 *