In the previous post I mentioned level set method can be used for offset surface. In this post I will use IPython to prototype the general idea.


The main issue of current Blender Solidify modifier is self intersection. It works well for most cases, however, self intersection happens when we offset too much into the shape.


We may expect the geometry could “merge” at the self intersected places. The level set method is good at handling this kind of topology change.

The problem modeling is based on J. A. Sethian’s book1. We first build a signed distance field() for the geometry. And then solve Hamilton-Jacobi equation to propagate the zero level of surface along the normal direction(the gradient) with a normal speed().

Since our geometry is represented in distance field implicitly, the topology change(i.e. merging or melding of surfaces) is handled for free.


The complete source can be found here

I use numpy, matplotlib, and numba(for performance):

import numpy as np
import matplotlib.pyplot as plt
import numba
import math

We first build a meshgrid to place the distance field of geometry.

# Discrete linear array on range [-2, 4) with step 0.001
x = np.arange(-2, 4, 0.001)
y = np.arange(-2, 4, 0.001)

X, Y = np.meshgrid(x, y)

Then we fill the meshgrid with distance field. In my test case I use two touching spheres. The distance field can be computed analytically for spheres. In practice we may traverse all polygons to calculate it for geometry represented in mesh.

def buildSphere(cx, cy, r):
    return lambda x, y: math.sqrt((x - cx)**2 + (y - cy)**2) - r

s0 = buildSphere(0.003, 0.003, 1)
s1 = buildSphere(2 - 0.01 + 0.003, 0.003, 1)

f = lambda x, y: min(s0(x, y), s1(x, y))
F = np.vectorize(f)
Z = F(X, Y)

Next we convect the distance field according to the Hamilton-Jacobi equation:

def Convect(X, Y, Z, dt, speed):
    hx = X[0, 1] - X[0, 0]
    hy = Y[1, 0] - Y[0, 0]
    row, column = Z.shape
    u = np.empty(Z.shape)
    for r in range(row):
        for c in range(column):
            z = Z[r, c]
            dx = (Z[r, (c + 1) % column] - Z[r, (c - 1) % column]) / (2 * hx)
            dy = (Z[(r + 1) % row, c] - Z[(r - 1) % row, c]) / (2 * hy)
            g = math.sqrt(dx * dx + dy * dy)
            u[r, c] = z - dt * speed * g
    return u

Note that for simplicity I use central differencing for spatial derivatives, and do not pay much attention on CFL condition. In a real world implementation some advanced scheme(e.g. upwind or other ones derived from conservation law) has to be used. And some numerical issues should be carefully handled. For now I just want to illustrate the general ideas.

We show the experiment results below. The first one contains no topology change. This happens when we do not offset the surface too much.


The next one has topology change. We can see that the offset surface breaks into two pieces.



We have to perform a mesh-to-distance-field operation in order to use level set method, and also a distance-field-to-mesh operation (e.g. dual contouring, marching cube, etc) in order to represent the offset surface in Blender mesh data structure.

The contouring method should be handled very carefully or non-manifold geometry would be introduced into our mesh.


  1. Sethian, James Albert. Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science. Vol. 3. Cambridge university press, 1999.