splines.mjs

// First, a little vector math library:
//
export const vec = (x, y) => ({ x, y });
export const mat = (a, b, c, d, e, f) => ({ a, b, c, d, e, f });

export const V = Object.freeze({
  init: vec,
  zero: vec(0, 0),
  add: (a, b) => vec(a.x + b.x, a.y + b.y),
  sub: (a, b) => vec(a.x - b.x, a.y - b.y),
  mul: (a, b) => vec(a.x * b.x, a.y * b.y),
  div: (a, s) => vec(a.x / s, a.y / s),
  scale: (a, s) => vec(a.x * s, a.y * s),
  neg: (a) => vec(-a.x, -a.y),
  dot: (a, b) => a.x * b.x + a.y * b.y,
  norm: (a) => Math.sqrt(V.dot(a, a)),
  normed: (a) => V.scale(a, 1.0 / V.norm(a)),
  normSq: (a) => V.dot(a, a),
  clockwise: (a) => vec(a.y, -a.x),
  // z value of the 3d cross product
  cross: (a, b) => a.x * b.y - a.y * b.x,
  isZero: (a) => a.x === 0 && a.y === 0,
  rot90: (a) => vec(-a.y, a.x),
});

// Affine matrix, representing a 2d rotate/scale/translate transform:
// a c e
// b d f
// 0 0 1
export const M = Object.freeze({
  init: mat,
  zero: mat(0, 0, 0, 0, 0, 0),
  identity: mat(1, 0, 0, 1, 0, 0),
  // Matrix-matrix multiplication
  mmul: (a, b) => ({
    a: a.a * b.a + a.c * b.b,
    b: b.a * a.b + b.b * a.d,
    c: a.a * b.c + a.c * b.d,
    d: a.b * b.c + a.d * b.d,
    e: a.a * b.e + a.c * b.f + a.e,
    f: a.b * b.e + a.d * b.f + a.f,
  }),
  // Matrix-vector multiplication. Treats B as a 3-vector whose last component is 1
  mul: (a, b) => ({
    x: a.a * b.x + a.c * b.y + a.e,
    y: a.b * b.x + a.d * b.y + a.f,
  }),
  // Determinant
  det: (a, b) => a.a * a.d - a.b * a.c,
  // Specialized matrix constructors
  translate: (x, y) => mat(1, 0, 0, 1, x, y),
  scale: (xScale, yScale) => mat(xScale, 0, 0, yScale, 0, 0),
  rotate: (angle) => {
    const c = Math.cos(angle);
    const s = Math.sin(angle);
    return mat(c, s, -s, c, 0, 0);
  },
});

// And now, the main event.

// For efficiency, this implementation contains improvements that
// take advantage of the shared sub-computations between contiguous curve patches
// to reduce the number of computations that need to be performed when rendering
// incrementally along a path of points. 
export function CurveToQuadraticBeziers() {
  //
}

// Init can be called multiple times to reinitialize the curve
// object with a new starting point. Use the same point as p0
// and p1 if you want the curve to go through it.
CurveToQuadraticBeziers.prototype.init = function (p0, p1, p2) {
  this.p0 = p0;
  this.p1 = p1;
  this.p2 = p2;

  const sub10 = V.sub(p1, p0);
  const sub12 = V.sub(p1, p2);

  this.norm01 = V.norm(sub10);
  this.norm12 = V.norm(sub12);
  this.sqrtNorm01 = Math.sqrt(this.norm01);
  this.sqrtNorm12 = Math.sqrt(this.norm12);

  return p1;
};

// Emit two quadratic curves that interpolate smoothly from p1 to p2,
// using p0 and p3 as control points. To go through the start and end
// points, one option is to emit the first and last points twice (ie.
// with p0 === p1 in init, and two calls to point() with the same p3.
CurveToQuadraticBeziers.prototype.point = function (p3) {
  const { p0, p1, p2, norm01, norm12, sqrtNorm01, sqrtNorm12 } = this;
  const sub10 = V.sub(p1, p0);
  const sub12 = V.sub(p1, p2);
  const sub23 = V.sub(p2, p3);
  const norm23 = V.norm(sub23);
  const sqrtNorm23 = Math.sqrt(norm23);
  const sub21 = V.neg(sub12);

  // Calculate q01, avoiding divide-by-zero when the denominator is zero
  // This code would be easier to follow if we had method chaining or |>.
  // See the accompanying Mathematica notebook in the static subdirectory.
  const q01 =
    norm01 === 0
      ? p1
      : V.add(
        p1,
        V.div(V.add(V.scale(sub21, norm01), V.scale(sub10, norm12)), 4 * sqrtNorm01 * (sqrtNorm01 + sqrtNorm12)),
      );

  // Calculate q11, avoiding divide-by-zero when the denominator is zero
  const q11 =
    norm23 === 0
      ? p2
      : V.add(
        p2,
        V.div(V.add(V.scale(sub23, norm12), V.scale(sub12, norm23)), 4 * sqrtNorm23 * (sqrtNorm12 + sqrtNorm23)),
      );

  // The shared common point is the mean of the two control points.
  // The value for the first control point of the second curve, q10,
  // is the same as q02.
  const q02 = V.scale(V.add(q01, q11), 0.5);

  // Update points that depend on p0, p1, and p2 for the next iteration.
  this.norm01 = norm12;
  this.norm12 = norm23;
  this.sqrtNorm01 = sqrtNorm12;
  this.sqrtNorm12 = sqrtNorm23;
  this.p0 = p1;
  this.p1 = p2;
  this.p2 = p3;

  // Emit two quadratic curves with a shared control point:
  //   {q00 = p1 = p0, q01, q02}
  //   {q10 = q02, q11, q12 = p2 = p1}
  // We don't return q00 any more, since it's initially the the
  // same as the point returned by init (p1) and then the same as
  // q12 from the previous point call. Meaning: 
  // moveTo(curve.init(…)); quadraticCurveTo(q01, q02); quadraticCurveTo(q11, q12);
  return { q01, q02, q11, q12: p2 };
};

// Approximate a cubic bezier segment by two quadratic bezier segments.
// See the `Quad Approximation of Cubic Curves` paper and the accompanying
// Mathematica notebook.
export function cubicBezierToQuadraticBeziers(cp0, cp1, cp2, cp3) {
  const cp1mul3 = V.scale(cp1, 3);
  const cp2mul3 = V.scale(cp2, 3);

  const cp0_plus_cp1mul3 = V.add(cp0, cp1mul3);
  const cp3_plus_cp2mul3 = V.add(cp3, cp2mul3);

  // const q00 = cp0;
  const q01 = V.scale(cp0_plus_cp1mul3, 0.25);
  const q02 = V.scale(V.add(cp0_plus_cp1mul3, cp3_plus_cp2mul3), 0.125);

  // const q10 = q02;
  const q11 = V.scale(cp3_plus_cp2mul3, 0.25);
  const q12 = cp3;

  // Emit two quadratic curves with a shared control point:
  //   {q00 = cp0, q01, q02}
  //   {q10 = q02, q11, q12}
  return { q00: cp0, q01, q02, q11, q12 };
}