Using numpy / scipy to calculate iso-surface from a 3D array

Question

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

In this particular case, each 1D column = myarray[i, j, :] ( column = myarray[i, j, :] ) of a three-dimensional array can be processed independently. So, what I would like to know is the last position of the index (2D array) where the function is equal to a specific value, for example myvalue .

Some (slow) code to illustrate:

 # myarray = 3D ndarray import numpy as np from scipy import interpolate result = np.zeros(nx, ny) z_values = np.arange(nz) for i in range(nx): for j in range(ny): f = interpolate.interp1d(my_array[i, j], z_values) result[i, j] = f(myvalue) 

I know this can be np.ndenumerate up a bit with np.ndenumerate and other tricks, but I was wondering if there is an even easier way to make this kind of iso-surface. I could not find anything in ndimage or other libraries. I know that mayavi2 and vtk have many tools for working with iso-surfaces, but my goal here is not visualization - I want to perform calculations on these iso-surface values, and not display them. In addition, many vtk isosurface methods seem to include polygons and the like, and I only need a 2D array of positions for each surface value.

Answer 1

Using only numpy , you can get a good solution using argsort , sort , take and proper array manipulation. The following function uses the weighted average to calculate 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 = 1 from 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 in xrange(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 in which there is more than one z corresponding to my_value . An example application for constructing iso-surfaces 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 various interpolation functions. Below is one example that takes the average of the two closest zs :

 def calc_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)/2 return z