Search code examples
pythonmatplotlibscikit-imagemplot3dmarching-cubes

How to set the origin for the mesh with mplot3d?


Following the example in scikit-image doc I generate a spherical surface mesh with marching cubes algorithm. I want to center the unit sphere shell at the origin defined by the x,y,z grid. However, I cannot do that, since I don't know how to put x,y,z info with mpl_toolkits.mplot3d.art3d.Poly3DCollection. Here is the code:

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
from skimage import measure
import numpy as np

x, y, z = np.ogrid[-4:4:20j, -4:4:20j, -4:4:20j]
r = np.sqrt(x ** 2 + y ** 2 + z ** 2)
verts, faces, normals, values = measure.marching_cubes_lewiner(r,level=1)
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111, projection='3d')
mesh = Poly3DCollection(verts[faces])
mesh.set_edgecolor('k')
ax.add_collection3d(mesh)
plt.show()

The problem is that the marching_cubes_lewiner function does not take x,y,z into account. How can I center the resulting sphere at 0,0,0 as implied by the grid?


Solution

  • The measure.marching_cubes_lewiner takes the indices of the points in the grid for calculating the topology. It does not seem to have a way to specify the actual grid and neither any offset to it.

    Hence you may manipulate the resulting verts in the desired way. I.e. one can first multiply by the difference between the grid points, effectively scaling the output, and then add the offset of the grid. In this case the transformation would be newverts = 0.42105 * oldverts - 4.

    import numpy as np
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d.art3d import Poly3DCollection
    from skimage import measure
    
    x, y, z = np.ogrid[-4:4:20j, -4:4:20j, -4:4:20j]
    r = np.sqrt(x ** 2 + y ** 2 + z ** 2)
    
    verts, faces, normals, values = measure.marching_cubes_lewiner(r, level=1)
    
    verts *= np.array([np.diff(ar.flat)[0] for ar in [x,y,z]])
    verts += np.array([x.min(),y.min(),z.min()])
    
    fig = plt.figure(figsize=(10, 10))
    ax = fig.add_subplot(111, projection='3d')
    mesh = Poly3DCollection(verts[faces])
    mesh.set_edgecolor('k')
    ax.add_collection3d(mesh)
    ax.set_xlim(-2, 2) 
    ax.set_ylim(-2, 2)
    ax.set_zlim(-2, 2)
    plt.show()
    

    enter image description here