This is Words and Buttons Online — a collection of interactive #tutorials, #demos , and #quizzes about #mathematics, #algorithms and #programming.

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.

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.

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.

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_xt₀³ + b_xt₀² + c_xt₀ + d_x = x₀
3a_xt₀² + 2b_xt₀ + c_x = dx₀ / dt₀
a_xt₁³ + b_xt₁² + c_xt₁ + d_x = x₁
3a_xt₁² + 2b_xt₁ + c_x = dx₁ / dt₁

a_xt₀³ + b_xt₀² + c_xt₀ + d_x = y₀
3a_xt₀² + 2b_xt₀ + c_x = dy₀ / dt₀
a_xt₁³ + b_xt₁² + c_xt₁ + d_x = y₁
3a_xt₁² + 2b_xt₁ + c_x = dy₁ / dt₁

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

d_x = x₀
c_x = dx₀ / dt₀
a_x + b_x + c_x + d_x = x₁
3a_x + 2b_x + c_x = dx₁ / dt₁

d_x = y₀
c_x = dy₀ / dt₀
a_x + b_x + c_x + d_x = y₁
3a_x + 2b_x + c_x = dy₁ / dt₁

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

// 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.