8

I haven't been able to find anything on this, maybe because I don't have the right nomenclature (i.e. I don't know exactly how to ask for it), but anyway, I have a 3D numpy array "a". I would like to identify and plot the 2D surface where a=0. To make clear, the data is double precision floats smoothly varying over 3D space. It is highly likely that the surface a=0 will "thread between" the points of the array, and not exactly lie right on any of them. So I need something that can interpolate to find the a=0 surface and plot it. Does matplotlib have a ready-made routine for doing this?

bob.sacamento
  • 6,283
  • 10
  • 56
  • 115
  • This seems like more of a numpy or scipy thing – kilojoules Jun 17 '16 at 14:36
  • 2
    I believe you want a volume slicer. Neither numpy, scipy or matplotlib (as far as I know) are, by default, prepared to do this. You can calculate slices between layers but you'll need to build the code to do it (I might try presenting a solution if this is acceptable to you). The library closer to matplotlib that does that is Mayavi (I think its only available in Python 2.x for now). You might also consider pyQtGraph, VTK or Vispy. – armatita Jun 17 '16 at 17:06
  • @armatita Thanks for the answer. At least now I know for sure I wasn't asking something trivial! If you want to build a code that can do the volume slicing, I will certainly be interested in seeing it! Thanks. – bob.sacamento Jun 17 '16 at 19:24

3 Answers3

5

With Plotly you can plot both planar and nonlinear slices in volumetric data: http://nbviewer.jupyter.org/github/empet/Plotly-plots/blob/master/Plotly-Slice-in-volumetric-data.ipynb

xecafe
  • 740
  • 6
  • 12
3

To "identify and plot the 2D surface where a=0", you just need to subset the data where a=0 and plot (shown below) If you are looking for a projection of the data onto that plane, then that is a little more complicated.

threeD = np.array([(x,y,z) for x in [0,1] for y in [1,2] for z in [5,6]])
twoD = np.array([(y,z) for (x,y,z) in threeD if x==0])
y,z = zip(*twoD)
plt.plot(y,z, 'o')
plt.xlim(0,3)
plt.ylim(0,7)

enter image description here

michael_j_ward
  • 4,369
  • 1
  • 24
  • 25
2

The volume slicing operation usually relies on a notion of interpolation and the the most typical are: Nearest neighbor, Linear, and Cubic. Notice these methods are adaptable to further number of dimensions (see for example Bilinear or Trilinear interpolation).

In this case you are saying that you have a volume from which you can retrieve slices in X, Y, Z (or an hybrid but lets not consider this case since its a new whole problem and would only bring confusion).

So considering, as example, a slice X=5 and X=6 you want to know how to obtain the X=5.5. Take a look into the example:

def linear_interpolation(p1, p2, x0):
    """
    Function that receives as arguments the coordinates of two points (x,y)
    and returns the linear interpolation of a y0 in a given x0 position. This is the
    equivalent to obtaining y0 = y1 + (y2 - y1)*((x0-x1)/(x2-x1)).
    Look into https://en.wikipedia.org/wiki/Linear_interpolation for more
    information.

    Parameters
    ----------
    p1     : tuple (floats)
        Tuple (x,y) of a first point in a line.
    p2     : tuple (floats)
        Tuple (x,y) of a second point in a line.
    x0     : float
        X coordinate on which you want to interpolate a y0.

    Return float (interpolated y0 value)
    """
    return p1[1] + (p2[1] - p1[1]) * ((x0 - p1[0]) / (p2[0] - p1[0]))

X, Y, Z = np.meshgrid(range(10), range(10), range(10))
vol = np.sqrt((X-5)**2 + (Y-5)**2 + (Z-5)**2)

Slice5dot5 = linear_interpolation((Y[5, :, :], vol[5, :, :]), (Y[6, :, :], vol[6, :, :]), 5.5)

f, (ax1, ax2, ax3) = plt.subplots(1, 3, sharey=True)
ax1.imshow(vol[5, :, :], interpolation='nearest', origin='lower', vmin=vol.min(), vmax=vol.max())
ax1.set_title("vol[5, :, :]")
ax2.imshow(Slice5dot5, interpolation='nearest', origin='lower', vmin=vol.min(), vmax=vol.max())
ax2.set_title("vol[5.5, :, :]")
ax3.imshow(vol[6, :, :], interpolation='nearest', origin='lower', vmin=vol.min(), vmax=vol.max())
ax3.set_title("vol[6, :, :]")
plt.show()

The function appears documented (its part of an old project I did) to be used with numbers but it will work just as well with numpy 2D slices (and it will be much faster than looping all those cells).

The result was this:

Volume slicing matplotlib

You'll notice the color getting paler from left to right. The slice in the middle is fully interpolated and did not exist previously.

armatita
  • 12,825
  • 8
  • 48
  • 49