pytessel

PyTessel

Purpose

Python package for building isosurfaces from 3D scalar fields

Installation

Anaconda

conda install -c ifilot pyqint

PyPi

pip install pytessel

Getting started

In the script below, the isosurface of a threedimensional Gaussian function is constructed. The isosurface is written to test.ply and can, for example, be opened using ctmviewer (Linux) or 3D viewer (Windows).

from pytessel import PyTessel
import numpy as np

def main():
    pytessel = PyTessel()

    # generate some data
    x = np.linspace(0, 10, 50)
    # the grid is organized with z the slowest moving index and x the fastest moving index
    grid = np.flipud(np.vstack(np.meshgrid(x, x, x, indexing='ij')).reshape(3,-1)).T

    R = [5,5,5]
    scalarfield = np.reshape(np.array([gaussian(r,R) for r in grid]), (len(x),len(x),len(x)))
    unitcell = np.diag(np.ones(3) * 10.0)

    vertices, normals, indices = pytessel.marching_cubes(scalarfield.flatten(), scalarfield.shape, unitcell.flatten(), 0.1)

    pytessel.write_ply('test.ply', vertices, normals, indices)

def gaussian(r, R):
    return np.exp(-(r-R).dot((r-R)))

if __name__ == '__main__':
    main()

Isosurface quality

In the script below, 6 different images are created of an icosahedral "metaball" using a grid of 10x10x10,20x20x20,25x25x25,50x50x50,100x100x100, and 200x200x200 points. The resulting .ply files are rendered using Blender.

from pytessel import PyTessel
import numpy as np

def main():
    """
    Build 6 .ply files of increasing quality
    """
    pytessel = PyTessel()

    for nrpoints in [10,20,25,50,100,200]:
        sz = 3

        x = np.linspace(-sz, sz, nrpoints)
        y = np.linspace(-sz, sz, nrpoints)
        z = np.linspace(-sz, sz, nrpoints)

        xx, yy, zz, field = icosahedron_field(x,y,z)

        unitcell = np.diag(np.ones(3) * sz * 2)
        pytessel = PyTessel()
        isovalue = 3.75
        vertices, normals, indices = pytessel.marching_cubes(field.flatten(), field.shape, unitcell.flatten(), isovalue)

        pytessel.write_ply('icosahedron_%03i.ply' % nrpoints, vertices, normals, indices)

def icosahedron_field(x,y,z):
    """
    Produce a scalar field for the icosahedral metaballs
    """
    phi = (1 + np.sqrt(5)) / 2
    vertices = [
        [0,1,phi],
        [0,-1,-phi],
        [0,1,-phi],
        [0,-1,phi],
        [1,phi,0],
        [-1,-phi,0],
        [1,-phi,0],
        [-1,phi,0],
        [phi,0,1],
        [-phi,0,-1],
        [phi,0,-1],
        [-phi,0,1]
    ]

    xx,yy,zz = np.meshgrid(x,y,z)
    field = np.zeros_like(xx)
    for v in vertices:
        field += f(xx,yy,zz,v[0], v[1],v[2])

    return xx,yy,zz,field

def f(x,y,z,X0,Y0,Z0):
    """
    Single metaball function
    """
    return 1 / ((x - X0)**2 + (y - Y0)**2 + (z - Z0)**2)

if __name__ == '__main__':
    main()

Local compilation (Linux)

python3 setup.py build_ext --inplace bdist