CODEWITHSUNDEEP
pythonmatplotlibplotscipysurface
bob.sacamento
bob.sacamento

Reputation: 6651

matplotlib 2D slice of 3D data

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?

Upvotes: 8

Views: 11652

Answers (3)

petem
petem

Reputation: 820

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

Upvotes: 5

armatita
armatita

Reputation: 13465

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.

Upvotes: 2

michael_j_ward
michael_j_ward

Reputation: 4559

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

Upvotes: 3

Related Questions