3

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?

Armut
  • 969
  • 8
  • 22
  • The function from the tutorial is called `measure.marching_cubes_lewiner` and it takes two arguments. – ImportanceOfBeingErnest Sep 17 '18 at 22:54
  • @ImportanceOfBeingErnest: The second argument is optional to set the values for the isosurfaces. marching_cubes used lewiner algorithm as the default, so both works. But, I edited the question to avoid confusion. – Armut Sep 17 '18 at 22:59
  • My version of skimage does not have a `marching_cubes`. But the code as it is now should work; what problem do you experience? – ImportanceOfBeingErnest Sep 17 '18 at 23:06
  • I want the sphere to be centered at 0,0,0. I don't know how to do that. – Armut Sep 17 '18 at 23:08

1 Answers1

3

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

ImportanceOfBeingErnest
  • 321,279
  • 53
  • 665
  • 712
  • Thank you very much for the answer. Could you elaborate on the operations on the vertices? By the way, `import numpy as np` is missing. – Armut Sep 17 '18 at 23:51