Draw a smooth curve through a sequence of points

Problem: Given a sequence of points, connect the points with a smooth path that can be drawn in an SVG file or graphics program:

A smooth curve through the sequence of points — After.

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 :

In this diagram, the two cubic Bézier curves, and , form a smooth path connecting the three points , , . Each cubic Bézier curve, , is completely determined by its control points, , , , and —which are exactly what you need to draw the curve in an SVG or graphics program.

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

}

/**

* This code is in the public domain.

* @author Spencer Mortensen

* @see https://spencermortensen.com/articles/draw-a-smooth-curve/

* @license https://creativecommons.org/public-domain/cc0/ Public domain

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

// 0*Q[-1] + 2*Q[0] + 1*Q[1] = P[0] + 2*P[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) {

// 1*Q[i-1] + 4*Q[i] + 1*Q[i+1] = 4*P[i] + 2*P[i+1]:

a[i] = 1;

b[i] = 4;

c[i] = 1;

d[i] = Point.plus(Point.times(4, P[i]), Point.times(2, P[i+1]));

}

// 2*Q[n-3] + 7*Q[n-2] + 0*Q[n-1] = 8*P[n-2] + P[n-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];

}

// The Thomas algorithm solves a tridiagonal system of equations:

// a[i]*x[i-1] + b[i]*x[i] + c[i]*x[i+1] = d[i]

// See: https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm

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;

}

/**

* This code is in the public domain.

* @author Spencer Mortensen

* @see https://spencermortensen.com/articles/draw-a-smooth-curve/

* @license https://creativecommons.org/public-domain/cc0/ Public domain

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>

<head>

</head>

</html>

'use strict';

window.addEventListener('DOMContentLoaded', function () {

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

// How can I see this?

console.log(svg.outerHTML);

});

// Point

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

}

// Solver

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;

// 0*Q[-1] + 2*Q[0] + 1*Q[1] = P[0] + 2*P[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) {

// 1*Q[i-1] + 4*Q[i] + 1*Q[i+1] = 4*P[i] + 2*P[i+1]:

a[i] = 1;

b[i] = 4;

c[i] = 1;

d[i] = Point.plus(Point.times(4, P[i]), Point.times(2, P[i+1]));

}

// 2*Q[n-3] + 7*Q[n-2] + 0*Q[n-1] = 8*P[n-2] + P[n-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];

}

// The Thomas algorithm solves a tridiagonal system of equations:

// a[i]*x[i-1] + b[i]*x[i] + c[i]*x[i+1] = d[i]

// See: https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm

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;

}

// Grapher

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="http://www.w3.org/2000/svg">

<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(Q[i].x - P[i].x);

const qy = this.__getNumberText(Q[i].y - P[i].y);

const rx = this.__getNumberText(R[i].x - P[i].x);

const ry = this.__getNumberText(R[i].y - P[i].y);

const sx = this.__getNumberText(P[i+1].x - P[i].x);

const sy = this.__getNumberText(P[i+1].y - P[i].y);

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(a.x);

const ay = this.__getNumberText(a.y);

const dx = this.__getNumberText(b.x - a.x);

const dy = this.__getNumberText(b.y - a.y);

return `M ${ax} ${ay} l ${dx} ${dy}`;

}

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(center.x);

const cy = this.__getNumberText(center.y);

const r = this.__getNumberText(this.circleRadius);

return `<circle cx="${cx}" cy="${cy}" r="${r}"/>`;

}

Grapher.prototype.__getNumberText = function (number) {

// return Math.round(number).toString();

return number.toFixed(3).replace(/\.?0+$/, '');

}

/**

* This code is in the public domain.

* @author Spencer Mortensen

* @see https://spencermortensen.com/articles/draw-a-smooth-curve/

* @license https://creativecommons.org/public-domain/cc0/ Public domain

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 :

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 .