Finding faces with skimage.measure.marching_cubes returns in weird offset based on meshgrid limits

A collegue of mine recently tried to visualise strain moduli using matplotlib. If you do not know what that is: Here's a link if you are interested. https://www.degruyter.com/document/doi/10.1524/zkri.2011.1413/html

It's a 3D visualisation of the anisotropy of the strain tensor for crystal structures when they are compressed under pressure.

For our purposes it is simply an isosurface in 3D space. You compute values for a 3D meshgrid using a function and then the marching cubes algorithm finds vertex points that map out the isosurface at a given value (let's say: every vertex point that has an iso value of 0.4).

Additionally, the marching_cubes algorithm finds all faces that belong to the given iso-surface. In the simplest case, there is no anisotropy and the result will be a perfect sphere.

This is the code:

import matplotlib.pyplot as plt
import numpy as np
from skimage import measure


def func(x, y, z):

    Z1111 = 1 
    Z1133 = 1/3 
    Z3333 = 1 

    return (Z1111 * (x**4 + 2 * x**2 * y**2 + y**4) + 6 * Z1133 * (x**2 * z**2 + y**2 * z**2) + Z3333 * z**4 )


a,b,c = -1,1,.1

X, Y, Z = np.mgrid[a:b:c, a:b:c, a:b:c]
vals = func(X,Y,Z)

iso_val = 0.4
verts, faces, murks, _ = measure.marching_cubes(vals, iso_val)

# shift_zero = np.array((a*10,a*10,a*10))
# verts = verts + shift_zero

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_trisurf(verts[:, 0], verts[:,1], faces, verts[:, 2],
                cmap='jet', lw=1)

ax.plot((0,0),(0,0),(-10,10),'r-',lw=3,zorder=100)
plt.show()

Here's where the unusual behavior comes in:

When you define the 3D-meshgrid and you use points between -1,1 with an interval of 0.1 for all directions (x,y,z), then the sphere is not centered around (0,0,0), but around (10,10,10). If you uncomment the shift_zero and verts = verts + shift_zero lines, then the offset can be compensated.

Is there something wrong with my script? Does the function for the isosurface include an intrinsic shift? Is this an issue with the implementation of marching_cubes?. I could just compute the center point of the vertices using verts.mean(axis=0) and subtract that from the vertices, but that seems like a cop-out.

You expect marching_cubes() to know of your meshgrids. It does not. It thinks in terms of indices into the voxel array you give it.

You have to do the calculations to turn its returned verts into something relative to your meshgrids.

Your meshgrids result from a:b:c, or np.arange(a, b, c) == [-1.0, -0.9, ..., 0.8, 0.9]

Note that this isn't symmetric. to be symmetric, it would have to include index 1.0.

The center of the 20x20x20 voxel volume is index 9.5 == (20-1)/2. The origin of your meshgrid, being (0,0,0), maps to indices (10,10,10) in the voxel grid: (0.0 - a) / c == 10.

The inverse maps from indices into your meshgrid: a + 10.0 * c == 0.0, so a vertex at voxel index 10 is at position 0.0 in your meshgrid.

You can just calculate

verts = a + verts * c

If you then change your plotting to draw a shorter line that doesn't go from +10 to -10, but from +1 to -1:

ax.plot((0,0),(0,0),(-1,+1),'r-',lw=3,zorder=100)

Then you get this: