The simplest possible smooth contouring algorithm

Contouring algorithm is something that turns your implicit function into a contour. Marching cubes and dual method are the most noted examples.

There is also a simpler algorithm. It may be not as efficient computationally, but due to it's simplicity, it can be easily improved to produce smooth contours and not just polylines. Here's an interactive plot to show how it works.

Show: distance; expected shape; grid.

Contour: grid cells; linear; smooth.

Undo last Clear

The algorithm accepts a distance function for input. This function denotes the distance to the expected curve for every point on the plane.

The first phase of the algorithm is finding the initial contour. Let's iterate through all the cells of a grid. If the distance function changes its sing between the neighboring cells' centers, we'll add a piece of contour between them. We have to put one horizontal and one vertical piece per cube to form a fully closed contour.

Undo last Clear

function draw_contour_ver_1(context, d_f, w, h, grid_size) { for(var y = 0; y <= h; y += grid_size) for(var x = 0; x <= w; x += grid_size) { var d_in_square = d_f(x + grid_size/2, y + grid_size/2); var d_in_square_on_the_left = d_f(x - grid_size/2, y + grid_size/2); var d_in_square_above = d_f(x + grid_size/2, y - grid_size/2); // draw a horizontal piece if there is a sign change if( Math.sign(d_in_square) != Math.sign(d_in_square_above)) { // draw a contour piece context.moveTo(x, y); context.lineTo(x + grid_size, y); } // draw a vertical piece if there is a sign change if( Math.sign(d_in_square) != Math.sign(d_in_square_on_the_left)) { // draw a contour piece context.moveTo(x, y); context.lineTo(x, y + grid_size); } } }

The second phase is attracting the contour nodes to the expected curve. The distance function shows how far every node is from the expected curve. And the direction to shift our node is also detectable from the distance function's gradient.

All we have to do is to add a gradient vector multiplied by an inverse distance to every node point.

Undo last Clear

function draw_contour_ver_2(context, d_f, w, h, grid_size) { const dx = function(x, y) { return ( d_f(x + gradient_step, y) - d_f(x, y) ) / gradient_step; } const dy = function(x, y) { return ( d_f(x, y + gradient_step) - d_f(x, y) ) / gradient_step; } for(var y = 0; y <= h; y += grid_size) for(var x = 0; x <= w; x += grid_size) for(var piece_no = 0; piece_no <= 1; ++piece_no) { // 1-st piece is horizontal, 2-nd - vertical var x_step = (piece_no == 0) ? grid_size : 0; var y_step = (piece_no == 1) ? grid_size : 0; // check if there is a border crossing var x_in_square = x + grid_size / 2; var y_in_square = y + grid_size / 2; var d_in_square = d_f(x_in_square, y_in_square); var d_in_neighbor = d_f(x_in_square - y_step, y_in_square - x_step); // sign change = border crossing if( Math.sign(d_in_square) != Math.sign(d_in_neighbor)) { // p1 var d_in_p1 = d_f(x, y); var p1x = x - dx(x, y) * d_in_p1; var p1y = y - dy(x, y) * d_in_p1; // p2 var d_in_p2 = d_f(x + x_step, y + y_step); var p2x = x + x_step - dx(x + x_step, y + y_step) * d_in_p2; var p2y = y + y_step - dy(x + x_step, y + y_step) * d_in_p2; // draw horizontal on piece_no = 0 or vertical on piece_no = 1 context.moveTo(p1x, p1y); context.lineTo(p2x, p2y); } } }

The third phase is adding curvature. Once again, we'll measure the function's gradient for every node. We'll use the partial derivatives from it to form a parametric 3-rd degree polynomial.

It probably sounds more complicated than it is. The 3-rd degree polynomials for x and y are the solutions for these systems:

a x t 0 3 + b x t 0 2 + c x t 0 + d x = x 0

3a x t 0 2 + 2b x t 0 + c x = dx 0 / dt 0

a x t 1 3 + b x t 1 2 + c x t 1 + d x = x 1

3a x t 1 2 + 2b x t 1 + c x = dx 1 / dt 1 a x t 0 3 + b x t 0 2 + c x t 0 + d x = y 0

3a x t 0 2 + 2b x t 0 + c x = dy 0 / dt 0

a x t 1 3 + b x t 1 2 + c x t 1 + d x = y 1

3a x t 1 2 + 2b x t 1 + c x = dy 1 / dt 1

If you parametrize the polynomials for t = [0..1], the systems will become:

d x = x 0

c x = dx 0 / dt 0

a x + b x + c x + d x = x 1

3a x + 2b x + c x = dx 1 / dt 1 d x = y 0

c x = dy 0 / dt 0

a x + b x + c x + d x = y 1

3a x + 2b x + c x = dy 1 / dt 1

Undo last Reset all

function draw_contour_ver_3(context, d_f, w, h, grid_size) { // spline related const polynomial_in_x = function(A, x){ var y = 0.0; for(var i = 0; i < A.length; ++i){ y += A[i] * Math.pow(x, i); } return y; } const spline_for = function(p1, p1d, p2, p2d) { // A = [ // [1, 0, 0, 0], // [0, 1, 0, 0], // [1, 1, 1, 1], // [0, 1, 2, 3]]; // B = [p1, p1d, p2, p2d] return [ p1, p1d, 3 * p2 - p2d - 3*p1 - 2*p1d, p2d + p1d - 2*p2 + 2*p1 ]; } const dx = function(x, y) { return ( d_f(x + gradient_step, y) - d_f(x, y) ) / gradient_step; } const dy = function(x, y) { return ( d_f(x, y + gradient_step) - d_f(x, y) ) / gradient_step; } for(var y = 0; y <= h; y += grid_size) for(var x = 0; x <= w; x += grid_size) for(var piece_no = 0; piece_no <= 1; ++piece_no) { // 1-st piece is horizontal, 2-nd - vertical var x_step = (piece_no == 0) ? grid_size : 0; var y_step = (piece_no == 1) ? grid_size : 0; // first square is this square, the second is a neighbor var d_in_square = [d_f(x + grid_size/2, y + grid_size/2), d_f(x + grid_size/2 - y_step, y + grid_size/2 - x_step)]; // sign change = border crossing if( Math.sign(d_in_square[0]) != Math.sign(d_in_square[1])) { // contour piece direction matters for normal calculations var cpd = Math.sign(d_in_square[piece_no]); // p1 var d_in_p1 = d_f(x, y); var p1x = x - dx(x, y) * d_in_p1; var p1y = y - dy(x, y) * d_in_p1; var d1x = dx(p1x, p1y) * grid_size; var d1y = dy(p1x, p1y) * grid_size; // p2 var d_in_p2 = d_f(x + x_step, y + y_step); var p2x = x + x_step - dx(x + x_step, y + y_step) * d_in_p2; var p2y = y + y_step - dy(x + x_step, y + y_step) * d_in_p2; var d2x = dx(p2x, p2y) * grid_size; var d2y = dy(p2x, p2y) * grid_size; // parametric curve var Px = spline_for(p1x, cpd*d1y, p2x, cpd*d2y); var Py = spline_for(p1y, (-cpd)*d1x, p2y, (-cpd)*d2x); // draw contour for(var k = 0; k < grid_size; k += 1) { var t1 = (k + 0) / grid_size; var t2 = (k + 1) / grid_size; var x1 = polynomial_in_x(Px, t1); var y1 = polynomial_in_x(Py, t1); var x2 = polynomial_in_x(Px, t2); var y2 = polynomial_in_x(Py, t2); context.moveTo(x1, y1); context.lineTo(x2, y2); } } } }

The concepts on which the algorithm is built may not be simple: partial derivatives, gradients, parametric polynomials. But the algorithm itself is: cells → lines → curves.