isosurface is interupted

[Kevin-Mattheus-Moerman]

As I mentioned a while back on Slack, I've been comparing this library to MarchingCubes.jl, I found that when one "under samples" the structure e.g. like the below code where nSteps=25. Then I notice the surface can get "crummy" interrupted.
With MarchingCubes.jl the surface remains intact longer it seems (if I go too sparse it will also speak up). Not sure if the "intactness" is wrong or if the "crummy/splits" is wrong but perhaps you can check.

using GeometryBasics
using Meshing
using GLMakie
using MarchingCubes

function padarray(A::Array{T,3}; padAmount = 1, padValue = T[0.0]) where T<:Real
    siz = size(A) # Get size of A 

    # Create padded image 
    B = Array{T,3}(undef,siz .+ (2*padAmount)) 
    B[1:padAmount,:,:]   .= padValue
    B[:,1:padAmount,:]   .= padValue
    B[:,:,1:padAmount]   .= padValue
    B[end-padAmount:end,:,:] .= padValue
    B[:,end-padAmount:end,:] .= padValue
    B[:,:,end-padAmount:end] .= padValue

    # Set centre of B to A
    @inbounds for i in 1:1:siz[1]
        @inbounds for j in 1:1:siz[2]
            @inbounds for k in 1:1:siz[3]
                B[i+padAmount,j+padAmount,k+padAmount] = A[i,j,k]                
            end
        end
    end
    return B
end

function getIsosurface(A; level=0.0, cap=false)    
    if cap == true          
        if isapprox(level,0.0,atol=1e-8)
            padValue=1e10
        else
            padValue = level + 1*10^(round(Int,log10(abs(level)))+10)        
            if level<0                  
                padValue *= 1.0
            end
        end
        println(padValue)
        mc = MarchingCubes.MC(padarray(A; padAmount = 1, padValue = padValue))
    elseif cap==false
        mc = MarchingCubes.MC(A,Int)
    end    
    MarchingCubes.march(mc,level)
    F = [TriangleFace{Int64}(f) for f in mc.triangles]
    V = [Point{3,Float64}(p) for p in mc.vertices]
    return F,V
end

function getIsosurface2(A; level=0.0, cap=false)
    if cap == true        
        if isapprox(level,0.0,atol=1e-8)
            padValue=1e10
        else
            padValue = level + 1*10^(round(Int,log10(abs(level)))+10)        
            if level<0                  
                padValue *= 1.0
            end
        end
        println(padValue)
        vv, ff = isosurface(padarray(A; padAmount = 1, padValue = padValue), Meshing.MarchingCubes(iso=level))
    else  
        vv, ff = isosurface(A, Meshing.MarchingCubes(iso=level))
    end
    F = [TriangleFace{Int64}(f) for f in ff]
    V = [Point{3,Float64}(v) for v in vv]
    return F,V
end

gyroid(v) = cos(v[1])*sin(v[2])+cos(v[2])*sin(v[3])+cos(v[3])*sin(v[1])
gyroid_shell(v) = max(gyroid(v)-0.4,-gyroid(v)-0.4)

nSteps = 25
xr,yr,zr = ntuple(_->range(0,pi*4,nSteps),3)
A = [gyroid_shell((x,y,z)) for x in xr, y in yr, z in zr]

level = 0.0
F1,V1 = getIsosurface(A; level=level, cap=true)
F2,V2 = getIsosurface2(A; level=level, cap=true)

# Visualization
fig = Figure(size=(800,800))

ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "MarchingCubes.jl")
if ~isempty(F1)
    hp1 = poly!(ax1,GeometryBasics.Mesh(V1,F1), strokewidth=1,color=:white,shading=FastShading,transparency=false)
end

ax2 = Axis3(fig[1, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Meshing.MarchingCubes")
if ~isempty(F2)
    hp2 = poly!(ax2,GeometryBasics.Mesh(V2,F2), strokewidth=1,color=:white,shading=FastShading,transparency=false)
end

fig

[Kevin-Mattheus-Moerman]

Here is a more minimal example that does not appear to suffer the interruptions:

using GLMakie
using GeometryBasics
using MarchingCubes
using Meshing

function getIsosurface(A; level=0.0)        
    mc = MarchingCubes.MC(A,Int)
    MarchingCubes.march(mc,level)
    F = [TriangleFace{Int64}(f) for f in mc.triangles]
    V = [Point{3,Float64}(p) for p in mc.vertices]
    return F,V
end

function getIsosurface2(A; level=0.0)
    vv, ff = isosurface(A, Meshing.MarchingCubes(iso=level))    
    F = [TriangleFace{Int64}(f) for f in ff]
    V = [Point{3,Float64}(v) for v in vv]    
    return F,V
end

nSteps = 10 # Use this to set the density 
xr,yr,zr = ntuple(_->range(0,pi*4,nSteps),3)
A = [ cos(x)*sin(y)+cos(y)*sin(z)+cos(z)*sin(x) for x in xr, y in yr, z in zr]

level = 0.0
F1,V1 = getIsosurface(A; level=level)
F2,V2 = getIsosurface2(A; level=level)

# Visualization
fig = Figure(size=(800,800))

ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "MarchingCubes.jl")
if ~isempty(F1)
    hp1 = poly!(ax1,GeometryBasics.Mesh(V1,F1), strokewidth=1,color=:white,shading=FastShading,transparency=false)
end

ax2 = Axis3(fig[1, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Meshing.MarchingCubes")
if ~isempty(F2)
    hp2 = poly!(ax2,GeometryBasics.Mesh(V2,F2), strokewidth=1,color=:white,shading=FastShading,transparency=false)
end

fig

Since with nSteps=10 this is very coarse indeed, I think the error perhaps relates to the shell version. In the above example I used the "sheet" Gyroid where previously I used the "thickened" version (gyroid_shell(v) = max(gyroid(v)-0.4,-gyroid(v)-0.4)).

[Kevin-Mattheus-Moerman]

@sjkelly I've done some more digging it looks like the issue occurs when the spacing between the sheets becomes too small. The value 0.4 in this gyroid_shell(v) = max(gyroid(v)-0.4,-gyroid(v)-0.4) sets the sheet spacing. If I sweep those for a given density (see below video) and go from large to too small we can see that both Meshing and MarchingCubes present with the Swiss cheese issue. My guess is this happens when the sheet spacing is equal to or smaller than the voxel grid spacing. In that case the isosurface can't seem to decide "what side a node should be one" thus the surface connectivity pops in and out of the front/back side. Maybe this is a normal issue that cannot be helped. It does however seem to happen sooner for your library so perhaps there is something to be improved upon. Next I'll check if MATLAB's isosurface function suffers the same Swiss cheese problem.

temp_triply_issue.mp4

[Kevin-Mattheus-Moerman]

Indeed MATLAB's isosurface function suffers the same fate, albeit later like for MarchingCubes, see the below gif animation.

So it seems a normal/expected behaviour, and not really a bug at all. But, like I said, it does happen sooner with Meshes, so perhaps you can do some digging, to explore the issue. Below is my lengthier code to change the spacing with a slider in Julia.

using GLMakie
using GeometryBasics
using MarchingCubes
using Meshing

function getIsosurface(A; level=0.0)        
    mc = MarchingCubes.MC(A,Int)
    MarchingCubes.march(mc,level)
    F = [TriangleFace{Int64}(f) for f in mc.triangles]
    V = [Point{3,Float64}(p) for p in mc.vertices]
    return F,V
end

function getIsosurface2(A; level=0.0)
    vv, ff = isosurface(A, Meshing.MarchingCubes(iso=level))    
    F = [TriangleFace{Int64}(f) for f in ff]
    V = [Point{3,Float64}(v) for v in vv]    
    return F,V
end

gyroid(v) = cos(v[1])*sin(v[2])+cos(v[2])*sin(v[3])+cos(v[3])*sin(v[1])
gyroid_shell(v,s) = max(gyroid(v)-s,-gyroid(v)-s)

function getmeshes(s,nSteps,level)
    xr,yr,zr = ntuple(_->range(0,pi*4,nSteps),3)
    A = [gyroid_shell((x,y,z),s) for x in xr, y in yr, z in zr]
    
    F1,V1 = getIsosurface(A; level=level)
    F2,V2 = getIsosurface2(A; level=level)
    return F1,V1,F2,V2
end

nSteps = 50
level = 0.0
s = 0.15

F1,V1,F2,V2 = getmeshes(s,nSteps,level)

# Visualization
strokewidth = 0
fig = Figure(size=(800,800))

ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "MarchingCubes.jl")
hp1 = poly!(ax1,GeometryBasics.Mesh(V1,F1), strokewidth=strokewidth,color=:white,shading=FastShading,transparency=false)

ax2 = Axis3(fig[1, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Meshing.MarchingCubes")
hp2 = poly!(ax2,GeometryBasics.Mesh(V2,F2), strokewidth=strokewidth,color=:white,shading=FastShading,transparency=false)

stepRange1 = range(0.01,0.8,50)
hSlider1 = Slider(fig[2, :], range = stepRange1, startvalue = 0,linewidth=30)

on(hSlider1.value) do stepVal
    F1,V1,F2,V2 = getmeshes(stepVal,nSteps,level)    
    hp1[1] = GeometryBasics.Mesh(V1,F1)
    hp2[1] = GeometryBasics.Mesh(V2,F2)
end

fig