# Empty graph using implicit_plot3d with contour-option

Hello everybody,

I wanted to make an isosurface plot from a 3d matrix, which contains random values from 0 to 1 at an equally spaced 3d grid.

Therefore I generated random numbers, and reshaped them in a 3d matrix. I interpolated the 3d matrix linearly with the RegularGridInterpolator from scipy. To make a 3D plot of it I am using the implicit_plot3d function of sage with a given contour value.

I get no errors, but in the end the graph is empty, which should not be in my opinion.

Here is my code:

from scipy.interpolate import griddata
import numpy as np
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
from matplotlib import cm
import scipy.interpolate

numbers=np.random.random_sample((1000,)) #generate random numbers
#print numbers

data = np.reshape(numbers, (10, 10, 10))  #reshape the numbers on a 3D matrix
#print data

xi = yi = zi = np.linspace(1, 10, 10)
interp = scipy.interpolate.RegularGridInterpolator((xi,yi,zi), data) #interpolate the 3d matrix with a function

#print xi

var('x,y,z')
#test = implicit_plot3d(interp,(x,1,10),(y,1,10),(z,1,10),contour=0.5)
test = implicit_plot3d(interp==0.5,(x,1,10), (y,1,10),(z,1,10),plot_points=60, color='seagreen')
test.show() #plot the the function at a certain value

#a=(5.,5.,2.)
#interp(a)


Any ideas on that?

Thank you very much!

edit retag close merge delete

Sort by » oldest newest most voted

The following works for me:

import numpy
import scipy.interpolate

n = 10
numbers = numpy.random.random_sample( ( n**3, ) )
data    = numpy.reshape( numbers, (n,n,n) )
xi = yi = zi = np.linspace(1, n, n)
interp = scipy.interpolate.RegularGridInterpolator((xi,yi,zi), data)
var('x,y,z')
test = implicit_plot3d( lambda x,y,z: interp( (x,y,z) )
, (x,1,n)
, (y,1,n)
, (z,1,n)
, plot_points=60, color='seagreen', contour=0.5)
test.show()


My favorite values for n were in the tests $3$ or $5$, to make the machine less angry. And for the numbers, it was numbers = [ 1 for _ in range(n**3) ] .

The difference is the use of lambda x,y,z: interp( (x,y,z) ) instead of interp.

more

Thank you very much for your answer, I wasn't familiar with lambda functions. However, I am still a little bit confused, why it takes so long to make the implicit_plot3d. In mathematica the in my opinion equivalent function would be ContourPlot3D (http://reference.wolfram.com/language...). I already used it to plot interpolated functions to 50x50x50 matrices and it seems to be considerably faster, or in other words with implicit_plot3D it doesn't seem to be possible to do that for n=50. Any ideas on that? Ps. I think it should be numpy.linspace instead of np.linspace (line 7)

Yes, thanks, that np.linspace(1, n, n) should be numpy.linspace(1, n, n) . (Too much copy+paste.)

For the performance... maybe the reason is the "complicated" type of the result of interpolation. With the above notations, we have for instance (this time):

sage: v = interp( (1,1,1) )
sage: v
array(0.4639585570135113)


and it is not obvious to make a number out of v. The trick is to call "v of nothing", v[ () ] to get:

sage: v[()]
0.4639585570135113


and the implicit plot also has to find this or his way. So i tried

test = implicit_plot3d( lambda x,y,z: interp( (x,y,z) )[()]
, (x,1,n), (y,1,n), (z,1,n)
, plot_points=60, color='seagreen', contour=0.5 )


Or maybe the interpolation kills the performance.

Thank you for your ideas! I do not think it is because of the interpolation; commenting out the lines test =... test.show() and starting sage takes about 2 seconds for the whole process (n=2). The trick with the brackets [()] doesn't seem to have any impact, for n=2 with and without brackets sage needs about 2 minutes on my computer. Later I will try the matplotlib contour3D function (not sure if they are suited to this problem), maybe they are faster... (http://matplotlib.org/mpl_toolkits/mp...)

this is probably slow not of implicit_plot3d itself, but because it goes through the Python function lambda x,y,z: interp( (x,y,z) ) (which is not statically typed). the plot data could be efficiently generated by a single call to interp, which accepts array arguments, eg. try interp(np.array([[1.0, 2.0, 3.0], [3.0, 5.0, 6.0]])). but then i don't know how to interface with Sage's implicit_plot3d.