Python package for building isosurfaces from 3D scalar fields
conda install -c ifilot pyqint
pip install pytessel
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()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()
python3 setup.py build_ext --inplace bdist