Code to slice open a Menger sponge

Last month the New York Times ran a story about a sculpture based on cutting open a “Menger sponge,” a shape formed by recursively cutting holes through a cube. All the holes are rectangular, but when you cut the sponge open at an angle, you see six-pointed stars.

Here are some better photos, including both a physical model and a computer animation. Thanks to Mike Croucher for the link.

I’ve written some Python code to take slices of a Menger sponge. Here’s a sample output.

The Menger sponge starts with a unit cube, i.e. all coordinates are between 0 and 1. At the bottom of the code, you specify a plane by giving a point inside the cube and vector normal to the plane. The picture above is a slice that goes through the center of the cube (0.5, 0.5, 0.5) with a normal vector running from the origin to the opposite corner (1, 1, 1).

from math import floor, sqrt
from numpy import empty, array
from matplotlib.pylab import imshow, cm, show

def outside_unit_cube(triple):
    x, y, z = triple
    if x < 0 or y < 0 or z < 0:
        return 1
    if x > 1 or y > 1 or z > 1:
        return 1
    return 0

def in_sponge( triple, level ):
    """Determine whether a point lies inside the Menger sponge
    after the number of iterations given by 'level.' """
    x, y, z = triple
    if outside_unit_cube(triple):
        return 0
    if x == 1 or y == 1 or z == 1:
        return 1
    for i in range(level):
        x *= 3
        y *= 3
        z *= 3

        # A point is removed if two of its coordinates
        # lie in middle thirds.
        count = 0
        if int(floor(x)) % 3 == 1:
            count += 1
        if int(floor(y)) % 3 == 1:
            count += 1
        if int(floor(z)) % 3 == 1:
            count += 1
        if count >= 2:
            return 0

    return 1

def cross_product(v, w):
    v1, v2, v3 = v
    w1, w2, w3 = w
    return (v2*w3 - v3*w2, v3*w1 - v1*w3, v1*w2 - v2*w1)

def length(v):
    "Euclidean length"
    x, y, z = v
    return sqrt(x*x + y*y + z*z)

def plot_slice(normal, point, level, n):
    """Plot a slice through the Menger sponge by
    a plane containing the specified point and having
    the specified normal vector. The view is from
    the direction normal to the given plane."""

    # t is an arbitrary point
    # not parallel to the normal direction.
    nx, ny, nz = normal
    if nx != 0:
        t = (0, 1, 1)
    elif ny != 0:
        t = (1, 0, 1)
    else:
        t = (1, 1, 0)

    # Use cross product to find vector orthogonal to normal
    cross = cross_product(normal, t)
    v = array(cross) / length(cross)

    # Use cross product to find vector orthogonal
    # to both v and the normal vector.
    cross = cross_product(normal, v)
    w = array(cross) / length(cross)

    m = empty( (n, n), dtype=int )
    h = 1.0 / (n - 1)
    k = 2.0*sqrt(3.0)

    for x in range(n):
        for y in range(n):
            pt = point + (h*x - 0.5)*k*v + (h*y - 0.5)*k*w
            m[x, y] = 1 - in_sponge(pt, level)
    imshow(m, cmap=cm.gray)
    show()

# Specify the normal vector of the plane
# cutting through the cube.
normal = (1, 1, 0.5)

# Specify a point on the plane.
point = (0.5, 0.5, 0.5)

level = 3
n = 500
plot_slice(normal, point, level, n)

Related post: A chip off the old fractal block