title: Using numpy/scipy to calculate iso-surface from 3D array
body: <p>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. </p>
  
  <p>In this particular case, each 1D column (<code>column = myarray[i, j, :]</code>) of the 3D array can be treated independently. So what I would like to know are the last index positions (2D array) where the function is equal to a certain value, say <code>myvalue</code>.</p>
  
  <p>Some (slow) code to exemplify:</p>
  
  <pre><code># 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)
  </code></pre>
  
  <p>I know this can be sped up a bit with <code>np.ndenumerate</code> and other tricks, but was wondering if there is already a simpler way of doing this kind of iso-surface. I couldn't find anything in <code>ndimage</code> or other libraries. I know that mayavi2 and vtk have a lot of tools to deal with iso-surfaces, but my aim here is not visualisation -- I want to perform calculations on those iso-surface values, not display them. Plus, a lot of the iso-surface methods of vtk seem to involve polygons and the like, and what I need is just a 2D array of positions for each iso-surface value.</p>