Problem: Given a sequence of points, connect the points with a smooth path that can be drawn in an SVG file or graphics program:
It turns out there is a simple algorithm that produces a gorgeous path: the smoothest path that can be built with cubic Bézier curves.
In a moment, I’ll provide the optimal algorithm. It accepts a list of points, , as its input, and computes the control points, and , needed to draw a smooth sequence of cubic Bézier curves, , connecting all of the points in :
Now, without further ado, here is the algorithm! This is a working implementation in JavaScript:
const solver = new Solver();
const P = [
new Point(0, 0),
new Point(80, 160),
new Point(160, 80)
];
const [Q, R] = solver.getControlPoints(P);
'use strict';
function Point (x, y) {
this.x = x;
this.y = y;
}
Point.plus = function (p, q) {
return new Point(p.x + q.x, p.y + q.y);
}
Point.minus = function (p, q) {
return new Point(p.x - q.x, p.y - q.y);
}
Point.times = function (c, a) {
return new Point(c * a.x, c * a.y);
}
Point.divide = function (a, c) {
return new Point(a.x / c, a.y / c);
}
'use strict';
function Solver () {}
Solver.prototype.getControlPoints = function (P) {
const Q = this.__getQ(P);
const R = this.__getR(P, Q);
return [Q, R];
}
Solver.prototype.__getQ = function (P) {
const [a, b, c, d] = this.__getProblem(P);
return this.__getSolution(a, b, c, d);
}
Solver.prototype.__getProblem = function (P) {
const a = [];
const b = [];
const c = [];
const d = [];
const n = P.length - 1;
a[0] = 0;
b[0] = 2;
c[0] = 1;
d[0] = Point.plus(P[0], Point.times(2, P[1]));
for (let i = 1; i < n-1; ++i) {
a[i] = 1;
b[i] = 4;
c[i] = 1;
d[i] = Point.plus(Point.times(4, P[i]), Point.times(2, P[i+1]));
}
a[n-1] = 2;
b[n-1] = 7;
c[n-1] = 0;
d[n-1] = Point.plus(Point.times(8, P[n-1]), P[n]);
return [a, b, c, d];
}
Solver.prototype.__getSolution = function (a, b, c, d) {
const x = [];
const n = a.length;
for (let i = 1; i < n; ++i) {
const w = a[i] / b[i-1];
b[i] = b[i] - (w * c[i-1]);
d[i] = Point.minus(d[i], Point.times(w, d[i-1]));
}
x[n-1] = Point.divide(d[n-1], b[n-1]);
for (let i = n - 2; -1 < i; --i) {
x[i] = Point.divide(Point.minus(d[i], Point.times(c[i], x[i+1])), b[i]);
}
return x;
}
Solver.prototype.__getR = function (P, Q) {
const R = [];
const n = Q.length;
for (let i = 0; i < n - 1; ++i) {
R[i] = Point.minus(Point.times(2, P[i+1]), Q[i+1]);
}
R[n-1] = Point.divide(Point.plus(P[n], Q[n-1]), 2);
return R;
}
While the algorithm above is technically all you need, here’s an extra bit of code, so you can see the points and the smooth curve that connects them:
<!DOCTYPE html>
<html lang="en">
<head>
<meta charset="utf-8">
<title></title>
<script src="index.js" defer>
</head>
<body></body>
</html>
'use strict';
window.addEventListener(
const P = [
new Point(0, 0),
new Point(80, 160),
new Point(160, 80)
];
const grapher = new Grapher();
const viewBox = grapher.getViewBox(P);
const svg = grapher.getSvg(viewBox, P);
document.body.appendChild(svg);
console.log(svg.outerHTML);
});
function Point (x, y) {
this.x = x;
this.y = y;
}
Point.plus = function (p, q) {
return new Point(p.x + q.x, p.y + q.y);
}
Point.minus = function (p, q) {
return new Point(p.x - q.x, p.y - q.y);
}
Point.times = function (c, a) {
return new Point(c * a.x, c * a.y);
}
Point.divide = function (a, c) {
return new Point(a.x / c, a.y / c);
}
function Solver () {}
Solver.prototype.getControlPoints = function (P) {
const Q = this.__getQ(P);
const R = this.__getR(P, Q);
return [Q, R];
}
Solver.prototype.__getQ = function (P) {
const [a, b, c, d] = this.__getProblem(P);
return this.__getSolution(a, b, c, d);
}
Solver.prototype.__getProblem = function (P) {
const a = [];
const b = [];
const c = [];
const d = [];
const n = P.length - 1;
a[0] = 0;
b[0] = 2;
c[0] = 1;
d[0] = Point.plus(P[0], Point.times(2, P[1]));
for (let i = 1; i < n-1; ++i) {
a[i] = 1;
b[i] = 4;
c[i] = 1;
d[i] = Point.plus(Point.times(4, P[i]), Point.times(2, P[i+1]));
}
a[n-1] = 2;
b[n-1] = 7;
c[n-1] = 0;
d[n-1] = Point.plus(Point.times(8, P[n-1]), P[n]);
return [a, b, c, d];
}
Solver.prototype.__getSolution = function (a, b, c, d) {
const x = [];
const n = a.length;
for (let i = 1; i < n; ++i) {
const w = a[i] / b[i-1];
b[i] = b[i] - (w * c[i-1]);
d[i] = Point.minus(d[i], Point.times(w, d[i-1]));
}
x[n-1] = Point.divide(d[n-1], b[n-1]);
for (let i = n - 2; -1 < i; --i) {
x[i] = Point.divide(Point.minus(d[i], Point.times(c[i], x[i+1])), b[i]);
}
return x;
}
Solver.prototype.__getR = function (P, Q) {
const R = [];
const n = Q.length;
for (let i = 0; i < n - 1; ++i) {
R[i] = Point.minus(Point.times(2, P[i+1]), Q[i+1]);
}
R[n-1] = Point.divide(Point.plus(P[n], Q[n-1]), 2);
return R;
}
function Grapher () {
this.margin = 5;
this.circleRadius = 3;
this.lineWidth = 2;
}
Grapher.prototype.getViewBox = function (points) {
const [min, max] = this.__getMinMax(points);
min.x -= this.margin;
min.y -= this.margin;
max.x += this.margin;
max.y += this.margin;
for (const p of points) {
p.x -= min.x;
p.y -= min.y;
}
return new Point(max.x - min.x, max.y - min.y);
}
Grapher.prototype.__getMinMax = function (points) {
const min = new Point(points[0].x, points[0].y);
const max = new Point(points[0].x, points[0].y);
for (const point of points) {
if (point.x < min.x) {
min.x = point.x;
} else if (max.x < point.x) {
max.x = point.x;
}
if (point.y < min.y) {
min.y = point.y;
} else if (max.y < point.y) {
max.y = point.y;
}
}
return [min, max];
}
Grapher.prototype.getSvg = function (viewBox, P) {
const parser = new DOMParser();
const svgXml = this.__getSvgXml(viewBox, P);
const fragment = parser.parseFromString(svgXml, 'image/svg+xml');
return fragment.firstChild;
}
Grapher.prototype.__getSvgXml = function (viewBox, P) {
const curveXml = this.__getCurveXml(P);
const pointsXml = this.__getPointsXml(P);
return `<svg viewBox="0 0 ${viewBox.x} ${viewBox.y}" xmlns=
<style>
path {
fill: none;
stroke: #000;
stroke-width: ${this.lineWidth};
}
</style>
${curveXml}
${pointsXml}
</svg>`;
}
Grapher.prototype.__getCurveXml = function (P) {
if (2 < P.length) {
return this.__getSplineXml(P);
} else if (P.length === 2) {
return this.__getLineXml(P[0], P[1]);
} else {
return '';
}
}
Grapher.prototype.__getSplineXml = function (P) {
const solver = new Solver();
const [Q, R] = solver.getControlPoints(P);
const d = this.__getSplineD(P, Q, R);
return this.__getPathXml(d);
}
Grapher.prototype.__getSplineD = function (P, Q, R) {
const segments = [
`M ${P[0].x} ${P[0].y} c`
];
for (let i = 0; i < Q.length; ++i) {
const qx = this.__getNumberText(
const qy = this.__getNumberText(
const rx = this.__getNumberText(
const ry = this.__getNumberText(
const sx = this.__getNumberText(
const sy = this.__getNumberText(
segments.push(`${qx} ${qy} ${rx} ${ry} ${sx} ${sy}`);
}
return segments.join(' ');
}
Grapher.prototype.__getPathXml = function (d) {
return `<path d="${d}"/>`;
}
Grapher.prototype.__getLineXml = function (a, b) {
const d = this.__getLineD(a, b);
return this.__getPathXml(d);
}
Grapher.prototype.__getLineD = function (a, b) {
const ax = this.__getNumberText(
const ay = this.__getNumberText(
const dx = this.__getNumberText(
const dy = this.__getNumberText(
return `M ${ax} ${ay}
}
Grapher.prototype.__getPointsXml = function (points) {
const circles = [];
for (const p of points) {
const circleXml = this.__getCircleXml(p);
circles.push(circleXml);
}
return circles.join("\n");
}
Grapher.prototype.__getCircleXml = function (center) {
const cx = this.__getNumberText(
const cy = this.__getNumberText(
const r = this.__getNumberText(
return `<circle cx="${cx}" cy="${cy}" r="${r}"/>`;
}
Grapher.prototype.__getNumberText = function (number) {
return number.toFixed(3).
}
Why does this work?
In this section, I’ll derive the formulas that make up the majority of the code above.
gives the coordinates of a point as it moves along the th Bézier curve from time to time , starting at and ending at .
So is the velocity of that point, and is its acceleration.
Here’s what we know so far:
Because the spline is :
Because the spline is :
Because the spline is :
This constrains every interior point along the path, but it doesn’t constraint the endpoints.
So we’ll add a constraint at the beginning:
…and at the end:
This linearizes the path near the endpoints. This is a reasonable choice if we don’t know how the curve should bend at the tips, but it’s not the only constraint we could have chosen.
Notice that we have equations, and the same number of unknowns (, , and ), so we’re ready to solve for .
The equation for the th cubic Bézier curve is:
for
We can use this equation to simplify those constraints from a moment ago:
The constraint is:
The constraint becomes (using and the constraint):
The constraint becomes (using and the , constraints):
The beginning constraint becomes (using and the constraint):
The end constraint becomes (using , , with , and with ):
Together, these equations form a system of linear equations, and it’s a type that is particularly easy to solve: The Thomas algorithm was designed to handle systems of equations like this (see the JavaScript implementation above).
After solving for , we’ll use the constraint to find . This gives us , , , and , which is everything we need to draw the smooth path connecting the points in .