Marching Cubes - Triangle Verticies & Faces coordinates too large


I am trying to use skimage.measure.marching_cubes() to find the isosurfaces of constant density in an NxNxN grid.

I have the following data:

a 64x64x64 grid and an estimated density value at each point from kernel density estimation. I have the probability desnity function stored in a 64x64x64 grid as variable pdf. I want to work out the isosurface for a floating point value p(x,y,z) i.e. 0.05 and then plot the surface on the grid. I am using the function as below:

verts, faces, _, _ = measure.marching_cubes(volume=pdf, level=p(x,y,z)) 

and it is returning verts and faces. However, the values of the triangles when indexing verts[faces] have coordinates that are not on the same scale as my grid points. For example, my maximum density value is 0.053177612148756395, my minimum = 3.037898412266252e-11 and the isosurface value = 0.005720253514050692 . For the triangles I am getting coordinates in the order of tens like [ 8. , 25.895016, 35. ].

From my knowledge, I have x, y and z values in F(x,y,z) and the c value (pdf) such that F(x,y,z) = c for each point in my grid and so this approach should work. For reference the limits of my x,y and z grid values are:

x: min = -2.994612860227619, max = 2.6805708407690005
y: min = -3.0461430547999266, max = 3.1709747732901796
z: min = -3.116856591599126, max = 2.9290962417638613

My questions are:

  • Am using the function correctly?
  • How do I scale the triangles into plottable coordinates on my 3D plot?
  • Do I need to use a triangulation lookup table for the 14 unique triangulations?


The coordinates are the coordinates in array coordinates, ie in the range [0, 64). If you want to scale them to the range e.g. (-3, 3), you would do e.g. verts_scaled = verts * 6 / 64 - 3, or if you want to scale each dimension separately:

mins = np.array([-2.99, -3.05, -3.12])
maxes = np.array([2.68, 3.17, 2.93])
ranges = maxes - mins
verts_scaled = verts * ranges / np.array(pdf.shape) - mins

Does that make sense?

To plot marching cubes data see this scikit-image gallery example:

Or you can use napari.view_surface((verts, faces)) — see the surface layer docs.


@jni thanks for getting back so quickly!

returned coordinates being in array coordinates makes a lot of sense. I am using Plotly for the graphing and taking the approach to scale each dimension separately I use go.mesh3d () with the x,y and z vertices, verts_scaled[:, 0], verts_scaled[:, 1] and verts_scaled[:, 2] respectively. It seems to work thereabouts but I think I need to adjust the scaling. For instance, when I use the approach suggested, my axes are [3, 6.5], [3, 7] and [3, 7] for x y and z. Could you elaborate on the scaling you suggested please?

I am trying to compare the isosurface calculated from marchin cubes with the isosurface plotted using the go.Isosurface plot in plotly. Please see the plotly image for both methods below:

Marching Cubes - scaling each dimension as suggested


go.Isosurface plot using the pdf (64,64,64) grid and specifying the p(x,y,z) constant value


I have used the same constant for both isosurfaces. Personally I would prefer to use the marching cubes algo for extensability

Thanks for any help!