Plotting 3D image form a data in NumPy-array

Question

I have a data file in NumPy array, I would like to view the 3D-image. I am sharing an example, where I can view 2D image of size (100, 100), this is a slice in xy-plane at z = 0.

import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
X, Y, Z = np.mgrid[-10:10:100j, -10:10:100j, -10:10:100j]
T = np.sin(X*Y*Z)/(X*Y*Z)
T=T[:,:,0]
im = plt.imshow(T, cmap='hot')
plt.colorbar(im, orientation='vertical')
plt.show()

How can I view a 3D image of the data T of shape (100, 100, 100)?

Solution 1

I have got a solution to my question. If we have the NumPy data, then we can convert them into TVTK ImageData and then visualization is possible with the help of mlab form Mayavi. The code and its 3D visualization are the following

from tvtk.api import tvtk
import numpy as np
from mayavi import mlab
X, Y, Z = np.mgrid[-10:10:100j, -10:10:100j, -10:10:100j]
data = np.sin(X*Y*Z)/(X*Y*Z)
i = tvtk.ImageData(spacing=(1, 1, 1), origin=(0, 0, 0))
i.point_data.scalars = data.ravel()
i.point_data.scalars.name = 'scalars'
i.dimensions = data.shape
mlab.pipeline.surface(i)
mlab.colorbar(orientation='vertical')
mlab.show()

For another randomly generated data

from numpy import random
data = random.random((20, 20, 20))

The visualization will be

Solution 2

I think the main problem is, that you do have 4 informations for each point, so you are actually interessted in a 4-dimensional object. Plotting this is always difficult (maybe even impossible). I suggest one of the following solutions:

1. You change the question to: I'm not interessted in all combinations of x,y,z, but only the ones, where z = f(x,y)

2. You change the accuracy of you plot a bit, saying that you don't need 100 levels of z, but only maybe 5, then you simply make 5 of the plots you already have.

In case you want to use the first method, then there are several submethods:

A. Plot the 2-dim surface f(x,y)=z and color it with T B. Use any technic that is used to plot complex functions, for more info see here.

The plot given by method 1.A (which I think is the best solution) with z=x^2+y^2 yields:

I used this programm:

import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib as mpl
X, Y = np.mgrid[-10:10:100j, -10:10:100j]
Z = (X**2+Y**2)/10 #definition of f
T = np.sin(X*Y*Z)
norm = mpl.colors.Normalize(vmin=np.amin(T), vmax=np.amax(T))
T = mpl.cm.hot(T) #change T to colors
fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(X, Y, Z, facecolors=T, linewidth=0,
       cstride = 1, rstride = 1)
plt.show()

The second method gives something like:

With the code:

norm = mpl.colors.Normalize(vmin=-1, vmax=1)
X, Y= np.mgrid[-10:10:101j, -10:10:101j]
fig = plt.figure()
ax = fig.gca(projection='3d')
for i in np.linspace(-1,1,5):
    Z = np.zeros(X.shape)+i
    T = np.sin(X*Y*Z)
    T = mpl.cm.hot(T)
    ax.plot_surface(X, Y, Z, facecolors=T, linewidth=0, alpha = 0.5, cstride 
        = 10, rstride = 10)

plt.show()

Note: I changed the function to T = sin(X*Y*Z) because dividing by X*Y*Zmakes the functions behavior bad, as you divide two number very close to 0.