Get the point(s) at which a line segment intersects a cubic bezier curve.

getSegmentCubicBezierIntersections.ts

// Adapted from https://github.com/w8r/bezier-intersect

/**
 * Get the point(s) at which a line segment intersects a cubic bezier curve.
 * @param p0x The x-coordinate of the curve's starting point.
 * @param p0y The y-coordinate of the curve's starting point.
 * @param p1x The x-coordinate of the curve's first control point.
 * @param p1y The y-coordinate of the curve's first control point.
 * @param p2x The x-coordinate of the curve's second control point.
 * @param p2y The y-coordinate of the curve's second control point.
 * @param p3x The x-coordinate of the curve's ending point.
 * @param p3y The y-coordinate of the curve's ending point.
 * @param x0 The x-coordinate of the segment's starting point.
 * @param y0 The x-coordinate of the segment's starting point.
 * @param x1 The x-coordinate of the segment's ending point.
 * @param y1 The y-coordinate of the segment's ending point.
 */
export function getSegmentCubicBezierIntersections(
  p0x: number,
  p0y: number,
  p1x: number,
  p1y: number,
  p2x: number,
  p2y: number,
  p3x: number,
  p3y: number,
  x0: number,
  y0: number,
  x1: number,
  y1: number
) {
  const result: number[][] = [];

  let ax: number,
    ay: number,
    bx: number,
    by: number,
    cx: number,
    cy: number,
    dx: number,
    dy: number, // temporary variables
    c3x: number,
    c3y: number,
    c2x: number,
    c2y: number,
    c1x: number,
    c1y: number,
    c0x: number,
    c0y: number, // coefficients of cubic
    cl: number, // c coefficient for normal form of line
    nx: number,
    ny: number; // normal for normal form of line

  // used to determine if point is on line segment
  const minx = Math.min(x0, x1),
    miny = Math.min(y0, y1),
    maxx = Math.max(x0, x1),
    maxy = Math.max(y0, y1);

  // Calculate the coefficients
  ax = p0x * -1;
  ay = p0y * -1;
  bx = p1x * 3;
  by = p1y * 3;
  cx = p2x * -3;
  cy = p2y * -3;
  dx = ax + bx + cx + p3x;
  dy = ay + by + cy + p3y;
  c3x = dx;
  c3y = dy; // vec

  ax = p0x * 3;
  ay = p0y * 3;
  bx = p1x * -6;
  by = p1y * -6;
  cx = p2x * 3;
  cy = p2y * 3;
  dx = ax + bx + cx;
  dy = ay + by + cy;
  c2x = dx;
  c2y = dy; // vec

  ax = p0x * -3;
  ay = p0y * -3;
  bx = p1x * 3;
  by = p1y * 3;
  cx = ax + bx;
  cy = ay + by;
  c1x = cx;
  c1y = cy; // vec

  c0x = p0x;
  c0y = p0y;

  // Convert line to normal form: ax + by + c = 0
  // Find normal to line: negative inverse of original line's slope
  nx = y0 - y1;
  ny = x1 - x0;

  // Determine new c coefficient
  cl = x0 * y1 - x1 * y0;

  // Rotate each cubic coefficient using line for new coordinate system?
  // Find roots of rotated cubic
  var roots = getPolynomialRoots(
    // dot products => x * x + y * y
    nx * c3x + ny * c3y,
    nx * c2x + ny * c2y,
    nx * c1x + ny * c1y,
    nx * c0x + ny * c0y + cl
  );

  // Any roots in closed interval [0,1] are intersections on Bezier, but
  // might not be on the line segment.
  // Find intersections and calculate point coordinates
  for (var i = 0; i < roots.length; i++) {
    var t = roots[i];

    if (0 <= t && t <= 1) {
      // We're within the Bezier curve
      // Find point on Bezier
      // lerp: y1 + (x2 - y1) * t
      var p5x = p0x + (p1x - p0x) * t;
      var p5y = p0y + (p1y - p0y) * t; // lerp(p1, p2, t);

      var p6x = p1x + (p2x - p1x) * t;
      var p6y = p1y + (p2y - p1y) * t;

      var p7x = p2x + (p3x - p2x) * t;
      var p7y = p2y + (p3y - p2y) * t;

      var p8x = p5x + (p6x - p5x) * t;
      var p8y = p5y + (p6y - p5y) * t;

      var p9x = p6x + (p7x - p6x) * t;
      var p9y = p6y + (p7y - p6y) * t;

      // candidate
      var p10x = p8x + (p9x - p8x) * t;
      var p10y = p8y + (p9y - p8y) * t;

      // See if point is on line segment
      if (x0 === x1) {
        // vertical
        if (miny <= p10y && p10y <= maxy) {
          result.push([p10x, p10y]);
        }
      } else if (y0 === y1) {
        // horizontal
        if (minx <= p10x && p10x <= maxx) {
          result.push([p10x, p10y]);
        }
      } else if (p10x >= minx && p10y >= miny && p10x <= maxx && p10y <= maxy) {
        result.push([p10x, p10y]);
      }
    }
  }
  return result;
}

// Helpers -----------------------------

const POLYNOMIAL_TOLERANCE = 1e-6;
const TOLERANCE = 1e-12;

export function getPolynomialRoots(...C: number[]) {
  var degree = C.length - 1;
  var n = degree;
  var results: number[] = [];
  for (var i = 0; i <= degree; i++) {
    if (Math.abs(C[i]) <= TOLERANCE) degree--;
    else break;
  }

  switch (degree) {
    case 1:
      getLinearRoots(C[n], C[n - 1], results);
      break;
    case 2:
      getQuadraticRoots(C[n], C[n - 1], C[n - 2], results);
      break;
    case 3:
      getCubicRoots(C[n], C[n - 1], C[n - 2], C[n - 3], results);
      break;
    default:
      break;
  }

  return results;
}

export function getLinearRoots(C0: number, C1: number, results: number[] = []) {
  if (C1 !== 0) results.push(-C0 / C1);
  return results;
}

export function getQuadraticRoots(
  C0: number,
  C1: number,
  C2: number,
  results: number[] = []
) {
  var a = C2;
  var b = C1 / a;
  var c = C0 / a;
  var d = b * b - 4 * c;

  if (d > 0) {
    var e = Math.sqrt(d);

    results.push(0.5 * (-b + e));
    results.push(0.5 * (-b - e));
  } else if (d === 0) {
    results.push(0.5 * -b);
  }

  return results;
}

export function getCubicRoots(
  C0: number,
  C1: number,
  C2: number,
  C3: number,
  results: number[] = []
) {
  var c3 = C3;
  var c2 = C2 / c3;
  var c1 = C1 / c3;
  var c0 = C0 / c3;

  var a = (3 * c1 - c2 * c2) / 3;
  var b = (2 * c2 * c2 * c2 - 9 * c1 * c2 + 27 * c0) / 27;
  var offset = c2 / 3;
  var discrim = (b * b) / 4 + (a * a * a) / 27;
  var halfB = b / 2;
  var tmp, root;

  if (Math.abs(discrim) <= POLYNOMIAL_TOLERANCE) discrim = 0;

  if (discrim > 0) {
    var e = Math.sqrt(discrim);

    tmp = -halfB + e;
    if (tmp >= 0) root = Math.pow(tmp, 1 / 3);
    else root = -Math.pow(-tmp, 1 / 3);

    tmp = -halfB - e;
    if (tmp >= 0) root += Math.pow(tmp, 1 / 3);
    else root -= Math.pow(-tmp, 1 / 3);

    results.push(root - offset);
  } else if (discrim < 0) {
    var distance = Math.sqrt(-a / 3);
    var angle = Math.atan2(Math.sqrt(-discrim), -halfB) / 3;
    var cos = Math.cos(angle);
    var sin = Math.sin(angle);
    var sqrt3 = Math.sqrt(3);

    results.push(2 * distance * cos - offset);
    results.push(-distance * (cos + sqrt3 * sin) - offset);
    results.push(-distance * (cos - sqrt3 * sin) - offset);
  } else {
    if (halfB >= 0) tmp = -Math.pow(halfB, 1 / 3);
    else tmp = Math.pow(-halfB, 1 / 3);

    results.push(2 * tmp - offset);
    // really should return next root twice, but we return only one
    results.push(-tmp - offset);
  }

  return results;
}