[BUG] Error in computation of Isosurface mesh even in case of no missing data

[aniruddhkb]

Plotly version -- 5.10.0
Environment -- Python 3.8.6 on Windows 10 x64

Below is a code example:

import plotly.graph_objects as go 
import numpy as np 


X = np.linspace(-10,10,100)
Y = np.linspace(-10,10,100)
Z = np.linspace(0,10,20)
coords = [X,Y,Z]
coord_mesh = np.concatenate([i.ravel().reshape(-1,1) for i in np.meshgrid(*coords)], axis=-1)


data = coord_mesh[:,-1].copy()


data_min, data_max = np.nanmin(data), np.nanmax(data)

DATA_REPLACE = 5
for i in range(data.size): 
    currcoords = coord_mesh[i,:]  
    if(currcoords[0]**2 + currcoords[1]**2 < 4):
        data[i] = DATA_REPLACE


ISOVAL = 6


fig = go.Figure(
    go.Isosurface(
        x=coord_mesh[:,0], 
        y= coord_mesh[:,1], 
        z = coord_mesh[:,2], 
        value=data, 
        caps=dict(x_show=False, y_show=False, z_show=False), 
        spaceframe=dict(show=False), 
        surface=dict(count=1), 
        isomin =ISOVAL, 
        isomax=ISOVAL, 
        cmin=data_min,
        cmax=data_max,
        lighting=
        dict(
            ambient=0.3,
            diffuse=1,
            fresnel=0,
            roughness=1,
            specular=1
            ), 
        lightposition=dict(
            x=0,
            y=0,
            z=1.75)
            )
        )

fig.show()

As you can see, the fourth dimension (data) within the cuboidal region is just the z-axis value EXCEPT in a cylindrical region at the center, where the data is replaced with DATA_REPLACE .

The output expected is a plane parallel to x-y with a circular hole in the middle.

The actual output:

The cylindrical boundary rising up is NOT expected behaviour. The value is not "6" or even close to 6 over that cylinder.

This is not simply a case of missing data, or incorrect imputation of missing data. This can affect any plot where the isosurface has a hole.

My intuition is that the code for computing the mesh expects that:

1. The mesh should be a closed surface
   OR
2. Any holes in the mesh should exist at the boundaries of the cuboidal domain.

To me this feels like a generalization of plotly/plotly.js#4192

Thanks,
AKB

PS: If you are already on this / if this is a duplicate issue / it doesn't fit the issue format, sorry for wasting your time.

[aniruddhkb]

I would like to amend my statement ----- Plotly 3D does not handle NULL values properly yet, therefore, when you replace the null values with a default value, artifacts occur which cannot be easily pruned manually. A suggested approach for handling NULL values would be to artificially set a very low value and prune the resulting artifacts.

I'll investigate this further later on.

[AaronStiff]

@alexcjohnson is this indeed a duplicate of the issue @aniruddhkb referenced in plotly.js?

[empet]

The example posted by @aniruddhkb is not conclusive for what he tries to point out,
because it requires to "carve" the isosurface corresponding to the value
6, i.e. a plane with a circular hole. But the algorithm marching cube, behind the go.Isosurface,
is implemented such that to get a curved surface that bounds a subvolume from the initial volume.

Calling the function measure_marching_cubes from scikit image we get the same surface:

import plotly.graph_objects as go 
from skimage import measure
import numpy as np
## Notice how simple it is to define the initial volume
X, Y, Z = np.mgrid[-10:10:100j, -10:10:100j, 0:10:50j] 
vals=Z.copy()
vals[np.where(X**2+Y**2<4)]=5
ISOVAL=6
spacing=(X[1,0, 0]-X[0,0,0], Y[0,1, 0]-Y[0,0,0], Z[0,0, 1]-Z[0,0,0])
verts, faces, _, _ = measure.marching_cubes(vals, ISOVAL, spacing=spacing)
verts =verts - np.array([10, 10, 0])
x, y, z = verts.T
I, J, K = faces.T
fig=go.Figure(go.Mesh3d(x=x,
                   y=y,
                   z=z,
                   intensity= z,
                   i=I,
                   j=J,
                   k=K,
                   name=''
                    ))
fig.update_scenes(zaxis_range=[0, 10.1])
fig.show()

I tried also to get the triangulation of the corresponding isosurface, with the Julia language package, MarchingCubes, https://github.com/JuliaGeometry/MarchingCubes.jl
and to plot the mesh with PlotlyJS.jl, the Julia version of plotly.py,
and got the same result, i.e. the returned isosurface is not a plane with a circular
hole, as the OP expected. And the reason is because the plane with a hole is a boundary
of two infinite volumes, i.e. he wanted the marching cube algorithm to work in conditions for which it was not designed.

L.E.: The plot returned by the above code (using marching_cube from scikit-image):

[alexcjohnson]

If I'm understanding the data correctly, I think the output we give in this case is correct: isosurface is trying to show where the data crosses the value 6; that happens at z=6 outside the hole (where value==z), and it happens at all z>6 on the boundary of the hole because at those locations we transition from z>6 to z=5 so an interpolation must pass through 6. @aniruddhkb is that what you meant by your "amended statement" #3938 (comment)? If so then I think we can close this issue and leave plotly/plotly.js#4192 open to follow up on better handling of null values.

[aniruddhkb]

To keep things in one place, I'll reply at plotly/plotly.js#4192

I agree; this can be closed.