Plot 3D isosurface in python

Question

I am trying to plot the following function in python using plotly or matplotlib for a given value of omega:

omega = (1/(12*np.pi*r**3)*((3*np.cos(THETA)**2-1)+1.5*np.sin(THETA)**2*np.cos(2*PHI)))

To do this I specify the value of omega, calculate r and then convert from polar to cartesian coordinates.

import numpy as np
import plotly.graph_objects as go

omega = 10
theta, phi = np.linspace(0,2*np.pi, 400), np.linspace(0,np.pi, 400)
THETA,PHI = np.meshgrid(theta,phi)

#Calculate R for a given value of omega
R1 = (1/(12*np.pi*omega)*((3*np.cos(THETA)**2-1)+1.5**np.sin(THETA)**2*np.cos(2*PHI)))
R = np.sign(R1) * (np.abs(R1))**(1/3)

#Convert to cartesians
X = R * np.sin(PHI) * np.cos(THETA)
Y = R * np.sin(PHI) * np.sin(THETA)
Z = R * np.cos(PHI)

#Plot isosurface plot
fig = go.Figure(data=[go.Surface(z=Z, x=X, y=Y)])
fig.show()

However, by doing this I lose the information about the positive and negative lobes, i.e. here both are plotted.

I can get round this by setting negative values of omega to NaN. This switches off the negative lobes but results in rendering artefacts for the graphs.

import numpy as np
import plotly.graph_objects as go

omega = 10
theta, phi = np.linspace(0,2*np.pi, 400), np.linspace(0,np.pi, 400)
THETA,PHI = np.meshgrid(theta,phi)

#Calculate R for a given value of omega
R1 = (1/(12*np.pi*omega)*((3*np.cos(THETA)**2-1)+1.5**np.sin(THETA)**2*np.cos(2*PHI)))
R = np.sign(R1) * (np.abs(R1))**(1/3)

#Remove negative lobes
R[R1 < 0.] = np.NaN

#Convert to cartesians
X = R * np.sin(PHI) * np.cos(THETA)
Y = R * np.sin(PHI) * np.sin(THETA)
Z = R * np.cos(PHI)

#Plot isosurface plot
fig = go.Figure(data=[go.Surface(z=Z, x=X, y=Y)])
fig.show()

I'm not sure how to overcome this issue - if I increase the number of points for THETA and PHI, the graph renders very slowly, and it's still not possible to increase the number of points sufficiently to remove the artefacts completely. Ideally I would pass r, theta and phi values to the plotting function and plot an isosurface at a given value of omega, but this is possible only in cartesians. Converting the whole function to cartesians would lead to f(x,y,z,r) which I also can't plot.

Stéphane Laurent · Answer 1 · 2021-10-15T19:39:10.447

You will find here a built-in way to plot an isosurface with pyvista.

Or you can use my implementation of the marching cubes (original code from the R package misc3d):

from math import cos, sin
def f(ρ, θ, ϕ):
    return (1/(12*np.pi*ρ**3)*((3*cos(θ)**2-1)+1.5*sin(θ)**2*cos(2*ϕ)))

ρmax = 0.3
θmax = 4*np.pi
ϕmax = 2*np.pi 
nρ = 200
nθ = 200
nϕ = 200

mesh = marchingCubes(
    f, 10, 0.001, ρmax, 0, θmax, 0, ϕmax, nρ, nθ, nϕ
)
print(mesh)

def sph2cart(sph):
    ρ = sph[0]
    θ = sph[1]
    ϕ = sph[2]
    return [
        ρ * cos(θ) * sin(ϕ),
        ρ * sin(θ) * sin(ϕ),
        ρ * cos(ϕ)
    ]
mesh.points = np.apply_along_axis(sph2cart, 1, mesh.points)

mesh.plot(smooth_shading=True, specular=5, color="purple")

I found that theta must go to 4pi and phi must go to 2pi, otherwise the surface is not closed. Here is the result:

Observe the hole at the center. It corresponds to the minimal value of rho (one can't set it to 0 otherwise there would be a division by zero).

Now, if you want to get this surface in plotly, you can extract the vertices and the faces (triangles) from the mesh object, and use Mesh3d in plotly:

vertices = mesh.points
triangles = mesh.faces.reshape(-1, 4)

import plotly.graph_objects as go
import plotly.io as pio


fig = go.Figure(data=[
    go.Mesh3d(
        x=vertices[:,0],
        y=vertices[:,1],
        z=vertices[:,2],
        colorbar_title='z',
        colorscale=[[0, 'gold'],
                    [0.5, 'mediumturquoise'],
                    [1, 'magenta']],
        # i, j and k give the vertices of triangles
        i=triangles[:,1],
        j=triangles[:,2],
        k=triangles[:,3],
        name='y',
        showscale=True
    )
])

#fig.show()

pio.write_html(fig, file='SO.html', auto_open=True)

To get colors:

distances = np.linalg.norm(vertices, axis=1)
distances = distances / (distances.max())

fig = go.Figure(data=[
    go.Mesh3d(
        x=vertices[:,0],
        y=vertices[:,1],
        z=vertices[:,2],
        colorbar_title='z',
        colorscale=[[0, 'gold'],
                    [0.5, 'mediumturquoise'],
                    [1, 'magenta']],
        # Intensity of each vertex, which will be interpolated and color-coded
        intensity=distances,
        # i, j and k give the vertices of triangles
        i=triangles[:,1],
        j=triangles[:,2],
        k=triangles[:,3],
        name='y',
        showscale=True
    )
])

Plot 3D isosurface in python

1 Answers1