Using Numpy/scipy To Calculate Iso-surface From 3d Array

I have a 3D numpy array that contains the values of a given function. I want to calculate a 2D iso-surface, or a set of iso-surfaces that represent certain values of this function.

Solution 1:

Using only numpy you can get a good solution using argsort, sort, take and the proper array manipulation. The function below uses a weighted average to compute the iso-surface:

def calc_iso_surface(my_array, my_value, zs, interp_order=6, power_parameter=0.5):
    if interp_order < 1: interp_order = 1from numpy import argsort, take, clip, zeros
    dist = (my_array - my_value)**2
    arg = argsort(dist,axis=2)
    dist.sort(axis=2)
    w_total = 0.
    z = zeros(my_array.shape[:2], dtype=float)
    for i inxrange(int(interp_order)):
        zi = take(zs, arg[:,:,i])
        valuei = dist[:,:,i]
        wi = 1/valuei
        clip(wi, 0, 1.e6, out=wi) # avoiding overflows
        w_total += wi**power_parameter
        z += zi*wi**power_parameter
    z /= w_total
    return z

This solution does not handle situations where there is more than one z corresponding to my_value. An application example to build the iso-surfaces below is given in the following code:

from numpy import meshgrid, sin, cos, pi, linspace
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
dx = 100; dy =  50; dz = 25
nx = 200; ny = 100; nz = 100
xs = linspace(0,dx,nx)
ys = linspace(0,dy,ny)
zs = linspace(0,dz,nz)
X,Y,Z = meshgrid( xs, ys, zs, dtype=float)
my_array = sin(0.3*pi+0.4*pi*X/dx)*sin(0.3*pi+0.4*pi*Y/dy)*(Z/dz)

fig = plt.figure()
ax = fig.gca(projection='3d')

z = calc_iso_surface( my_array, my_value=0.1, zs=zs, interp_order=6 )
ax.plot_surface(X[:,:,0], Y[:,:,0], z, cstride=4, rstride=4, color='g')

z = calc_iso_surface( my_array, my_value=0.2, zs=zs, interp_order=6 )
ax.plot_surface(X[:,:,0], Y[:,:,0], z, cstride=4, rstride=4, color='y')

z = calc_iso_surface( my_array, my_value=0.3, zs=zs, interp_order=6 )
ax.plot_surface(X[:,:,0], Y[:,:,0], z, cstride=4, rstride=4, color='b')

plt.ion()
plt.show()

You can also play with different interpolation functions. See below one example that takes the average of the two closest zs:

defcalc_iso_surface_2(my_array, my_value, zs):
    '''Takes the average of the two closest zs
    '''from numpy import argsort, take
    dist = (my_array - my_value)**2
    arg = argsort(dist,axis=2)
    z0 = take(zs, arg[:,:,0])
    z1 = take(zs, arg[:,:,1])
    z = (z0+z1)/2return z