How to generate bounded harmonic functions on regular grids





Generating an height-map  is a good example to demonstrate how useful harmonic functions can be. In (a) user can stroke lines corresponding to different heights. In (b) we associate a height where the strokes lies. In (c) we compute the rest of the values automatically (black is zero, white is one). (d) is a 3D view of the heightmap in (c). It is a 2D function \( f: \mathbb R^2 \rightarrow \mathbb R\) which returns the height. The tutorial explains how to compute the unknown heights and ensure \(f\) is an harmonic function.


In this tutorial I give the basics to be able to compute harmonic functions \( f: \mathbb R^n \rightarrow \mathbb R\) over a regular grid/texture (in 1D, 2D or 3D). All you need is to set a few values on the grid that define a closed region (see figure above), then the remaining values are automatically computed using the so-called laplace operator and solving a linear system of equations. After a quick and "mathless" presentation of an easy way to implement this in C++, I will explain the mathematical jargons. This tutorial is the first of series which will lead to more complicated techniques like Finite Element Methods (FEM).



A practical example in 2D

Our goal is to generate a 2D function \( f: \mathbb R^2 \rightarrow \mathbb R\) discretized and stored in a 2D array array[x][y] = f(x,y):


First we need to define by hand some values of the grid. These values must form a closed boundary, below I give an example of such closed region:


Note that the boundary doesn't necessarily have to be on the border of the grid. For instance we could have defined the closed region with a circle or a rectangle. In the end we will be able to compute values inside the domain defined by the closed boundary. Using the algorithm explained aftwards, this is what \( f \) looks like with these boundary conditions:

harmonic_isos 3D view
values of \( f(x, y) \) are shown with level curves. The curves colors mean: 5 for Pink, below 5 blue, above 5 red and black for 10. The same plot of \(f\) but in 3D (view from top).

Here \(f\) is harmonic: the function respects certain properties and will always have the same shape for a given closed boundary.

The function \( f \) is called differently depending on the context. In mathematics we refer to it as a scalar-function or scalar-field. Indeed the function associates a point in space \( \mathbf p \) (2D, 3D or  otherwise) to a scalar (i.e a single real number). When \( f \) returns a height like the introductory example we may call it a height-field. Sometimes \( f \) is called a weight-map as it maps every points of the plane/volume/etc. to a weight/a real number.

After presenting the algorithm which generates harmonic functions, I will explain their properties. But before I really dive into the technicalities and explanations, lets see how all this is useful in the every day life of a computer graphics nerd!


Being able to compute weight maps by only specifying a few values on a canva is extremely usefull in computer graphics. As described above it can be used to generate heightmaps, but not only. There is tons of applications where you need to generate such smooth transition of values / smooth shade over a surface. One of my favourites is cage deformations, where you enclose an object into a cage and deform it by moving the cage's borders:

A circle in a cage in rest pose. When moving the cage's vertex we deform the circle. The influence of one vertex of the cage. Blue means full influence. This influence is defined as an harmonic function.

Here are more concrete and convincing examples of such deformations:

(These examples actually use bi-harmonic functions but the idea is the same)
Cage deformations from "Bounded Biharmonic Weights for Real-Time Deformation" by Jacobson et Al

If you want to understand how these weight maps are used to deform objects, you can find here a nice essay introducing harmonic coordinates specifically for cage-based deformations (the article doesn't go as far as I plan to do though).

Although I have only presented 2D examples on regular grids, harmonic functions can be computed on various domains. For instance cage deformations can be extended to 3D. On a regular 3D grid we can define an harmonic function \( f: \mathbb R^3 \rightarrow \mathbb R\). This function will tell the influence inside a triangular 3D cage:

We can compute an harmonic function \( f: \mathbb R^3 \rightarrow \mathbb R\) inside a 3D grid A triangular cage defines a closed boundary inside this 3D grid. From there we can compute a weight map (or region of influence) for each cage's vertex. Result of the deformation  with 3D harmonic functions.

But we are not limited to regular grids may it be 2D or 3D. We can also compute these harmonic functions over irregular domains like triangular meshes:

In this case the function is 2D \( f: \mathbb R^2 \rightarrow \mathbb R\) and defines a weight for each vertex of the triangular 3D mesh. Here I represented the boundaries with red and blue lines. The function \( f \) returns 1 for red vertices and decays to zero for blue vertices.This is particularly useful to define what is called skinning weights. These are needed to animate characters ( as a side note computing those skinning weights with harmonic functions is briefly described in context aware deformation section 2.2. I intend to explain this in the series as well)

The diffusion algorithm

Ok, back to our simple 2D problem. We have our 2D rectangular grid and defined by hand some values forming a closed boundary. To compute automatically the values inside we will use a process called diffusion. The algorithm is dead simple:

The idea is that each time you look up the grid and compute the average, the values from the boundarie "leaks out"/diffuse to their direct neighbors. The more you repeat the process the more the values propagate inside the domain untill it stabilizes. I always found amazing the fact that, if you don't change the values at the boundaries then what is inside will eventually converge to a steady state. For the same boundary values, you will always get the same values inside if you iterate enough this averaging process. In other word given a boundary, you always find a unique harmonic function \( f \) (if it exists). Here we found this function using the diffusion algorithm, but I will explain other ways in a couple of sections.


The C++ procedure below compute nb_iter iteration of the diffusion algorithm for the 2D case:

/// (untested code)
bool inside(const Grid_2d& grid, Vec2i idx)
    return idx.x() < grid.size().x() && idx.x() > -1 &&
           idx.y() < grid.size().y() && idx.y() > -1

/// @param grid : 2d grid of float values you want to diffuse
/// @param mask : 2d grid where values set to '0' means don't diffuse this position
void diffuse(Grid_2d& grid, const Grid_2d& mask, int nb_iter)
    assert(grid.size() == mask.size());

    // First we store every grid elements that need to be diffused
    // to 'mask'
    int mask_nb_elts = mask.size().product();
    std::vector idx_list;
    idx_list.reserve( mask_nb_elts );
    for (int x = 0; x < mask.size().x(); ++x) {
        for (int y = 0; y < mask.size().y(); ++y) {
            if (mask(x, y) == 1)
                idx_list.push_back(Vec2i(x, y));

    // Diffuse values with a Gauss-seidel like scheme
    for (int i = 0; i < nb_iter; ++i)
        // Look up grid elements to be diffused
        for (unsigned j = 0; j < idx_list.size(); ++j)
            Vec2i idx         = idx_list[i];
            float sum         = 0.f;
            int nb_neighbours = 0;

            Vec2i above = idx + Vec2i( 0,  1);
            Vec2i below = idx + Vec2i( 0, -1);
            Vec2i left  = idx + Vec2i(-1,  0);
            Vec2i right = idx + Vec2i( 1,  0);
            if( inside(grid, above) ) { nb_neighbours++;  sum += grid(above.x(), above.y()) }
            if( inside(grid, below) ) { nb_neighbours++;  sum += grid(below.x(), below.y()) }
            if( inside(grid, left ) ) { nb_neighbours++;  sum += grid( left.x(),  left.y()) }
            if( inside(grid, right) ) { nb_neighbours++;  sum += grid(right.x(), right.y()) }

            float average = sum / (float)nb_neighbours;
            grid(idx.x(), idx.y()) = average;

The 1D case

Diffusion works in \( nD \) and we can try it out on a \( 1D \) function to better see what is happenning. Begin to store a function \( g : \mathbb R \Rightarrow \mathbb R\) in a one dimensionnal array. The array is intialized to zero except for the first and last  elements (for instance: array[0] = 1.0 and array[size-1] = 3.0). Then we apply the diffusion by looking the array from beginning to end and compute the average ( array[i-1] + array[i+1]) / 2.0:

harmonic 1D begin
relaxation 1d
Plot of \( g \) before diffusion The diffusion process for every iterations until convergence.

Again \( g \) converge to a harmonic solution.

Initial values

The grid must be initialized to some intial value/intial guess before diffusion. The choice of this starting point will impact the rapidity of convergence. If choose wisely this can drastically speed up the diffusion by decreasing the number of iterations needed:

Before diffusion The diffusion converge faster than above

The 3D case

Extension to 3D grids is trivial, the boundary must define a closed volume (as opposed to a closed area in 2D) and the diffusion now average the 6th neighbors instead of the 4th neighbors in 2D. The number of iterations before it stabilizes will be extremly important as there is usually more grid's nodes in 3D than in 2D. But the diffusion algorithm is the slowest method to compute Harmonic functions. I will introduce more efficient methods afterwards.

More examples

Here I show more examples of harmonic functions for different boundary conditions:

3d view

3D view
Boundary conditions and associated harmonic function \( f: \mathbb R^2 \rightarrow \mathbb R\) plotted with level-curves. The harmonic function in 3D view

In the second example notice that we can add fixed values inside the boundary as well.

Evaluate \(f\)

So far I have described how to compute a harmonic function at the location of the nodes of a grid, but what about in between those nodes? For instance in the 1D case you know the values \(f(1), f(2), \dots, f(n)\), but how do you compute \( f(0.6) \)? Well you need to interpolate the values of the nodes. You may choose the nearest value to \(f(0.6)\) or even do a linear interpolation. This is already an advanced tutorial leading to more complicated ones, so I don't feel explaining interpolation methods here is appropriate. I might write a separate post or find good links.


We have seen how to compute harmonic functions over regular grids in 1D, 2D or 3D with the diffusion process. I also gave some examples of use like cage-based deformations or heightmaps. In the next tutorial we will talk about their mathematical properties. Knowing and understanding these properties will be very useful to use them correctly or compute them more efficiently.

No comments

All html tags except <b> and <i> will be removed from your comment. You can make links by just typing the url or mail-address.
Spam bot question: