Source code in the c programming language

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


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