MarchingCubes: shared vertices, gradient normals, smooth_sdf

.github/workflows/CI.yml

@@ -10,37 +10,30 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.9'
-          - 'nightly'
+          - '1.10'
+          - '1.11'
+          - '1.12'
+          - '1.13'
         os:
           - ubuntu-latest
           - macOS-latest
         arch:
           - x64
     steps:
-      - uses: actions/checkout@v2
-      - uses: julia-actions/setup-julia@v1
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: actions/cache@v1
-        env:
-          cache-name: cache-artifacts
-        with:
-          path: ~/.julia/artifacts
-          key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
-          restore-keys: |
-            ${{ runner.os }}-test-${{ env.cache-name }}-
-            ${{ runner.os }}-test-
-            ${{ runner.os }}-
+      - uses: julia-actions/cache@v2
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
   docs:
     name: Documentation
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v2
-      - uses: julia-actions/setup-julia@v1
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
         with:
           version: '1'
       - run: |

Project.toml

@@ -5,7 +5,7 @@ version = "0.7.0"
 [deps]
 
 [compat]
-julia = "1.9"
+julia = "1.10"
 
 [extras]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"

src/Meshing.jl

@@ -7,6 +7,8 @@ include("marching_cubes.jl")
 include("isosurface.jl")
 
 export isosurface,
+       isosurface_normals,
+       smooth_sdf,
        MarchingCubes,
        MarchingTetrahedra
 

src/algorithmtypes.jl

@@ -11,18 +11,22 @@ abstract type AbstractMeshingAlgorithm end
 
 
 """
-    MarchingCubes(;iso=0.0)
+    MarchingCubes(;iso=0.0, reduceverts=false)
 
 Specifies the use of the Marching Cubes algorithm for isosurface extraction.
 This algorithm provides a good balance between performance and vertex count.
 In contrast to the other algorithms, vertices may be repeated, so mesh size
 may be large and it will be difficult to extract topological/connectivity information.
 
 - `iso` (default: 0.0) specifies the iso level to use for surface extraction.
+- `reduceverts` (default: false) shares vertices between neighbouring voxels, giving a
+  connected mesh with about a quarter of the vertices.
 """
 Base.@kwdef struct MarchingCubes{T} <: AbstractMeshingAlgorithm
     iso::T = 0.0
+    reduceverts::Bool = false
 end
+MarchingCubes(iso) = MarchingCubes(iso, false)  # the positional form from before `reduceverts`
 
 """
     MarchingTetrahedra(;iso=0.0, eps=1e-3)

src/common.jl

@@ -27,3 +27,34 @@ Called after `_get_cubeindex`. Determines if a voxel index has triangles.
 @inline function _no_triangles(cubeindex::UInt8)
     cubeindex == 0x00 || cubeindex == 0xff
 end
+
+"""
+    smooth_sdf(sdf; sigma=0.7)
+
+`sdf` blurred by a Gaussian of `sigma` voxels, edges clamped. Suppresses the banding a
+grid-aligned level set gives the extracted surface. `sigma < 0.3` returns a copy.
+"""
+function smooth_sdf(sdf::AbstractArray{T,3}; sigma::Real=0.7) where {T}
+    sigma < 0.3 && return copy(sdf)
+    r = ceil(Int, 3sigma)
+    w = [exp(-(k / sigma)^2 / 2) for k in -r:r]
+    w = float(T).(w ./ sum(w))
+    a, b = similar(sdf, float(T)), similar(sdf, float(T))
+    blur!(a, sdf, w, 1)
+    blur!(b, a, w, 2)
+    blur!(a, b, w, 3)
+end
+
+function blur!(dst, src, w, axis)
+    r = length(w) ÷ 2
+    n = size(src, axis)
+    @inbounds for I in CartesianIndices(src)
+        acc = zero(eltype(dst))
+        for k in -r:r
+            J = CartesianIndex(ntuple(a -> a == axis ? clamp(I[a] + k, 1, n) : I[a], Val(3)))
+            acc += w[k+r+1] * src[J]
+        end
+        dst[I] = acc
+    end
+    dst
+end

src/marching_cubes.jl

@@ -25,18 +25,43 @@ Voxel corner and edge indexing conventions
 =#
 
 function isosurface(sdf::AbstractArray{T,3}, method::MarchingCubes, X=-1:1, Y=-1:1, Z=-1:1) where {T}
+    vts, fcs, _ = mc_isosurface(sdf, method, X, Y, Z, Val(false))
+    vts, fcs
+end
+
+"""
+    isosurface_normals(sdf, method::MarchingCubes, X=-1:1, Y=-1:1, Z=-1:1) -> (vertices, faces, normals)
+
+Like [`isosurface`](@ref), with a unit normal per vertex from the gradient of `sdf`.
+"""
+isosurface_normals(sdf::AbstractArray{T,3}, method::MarchingCubes, X=-1:1, Y=-1:1, Z=-1:1) where {T} =
+    mc_isosurface(sdf, method, X, Y, Z, Val(true))
+
+function mc_isosurface(sdf::AbstractArray{T,3}, method::MarchingCubes, X, Y, Z, ::Val{normals}) where {T,normals}
     nx, ny, nz = size(sdf)
 
     # find widest type
     FT = promote_type(eltype(first(X)), eltype(first(Y)), eltype(first(Z)), eltype(T), typeof(method.iso))
 
     vts = NTuple{3,float(FT)}[]
     fcs = NTuple{3,Int}[]
+    nms = normals ? NTuple{3,float(FT)}[] : nothing
 
     xp = LinRange(first(X), last(X), nx)
     yp = LinRange(first(Y), last(Y), ny)
     zp = LinRange(first(Z), last(Z), nz)
 
+    # the vertex of each crossed grid edge, when vertices are shared
+    ids = method.reduceverts ? Dict{Int,Int}() : nothing
+    mc_voxels!(vts, nms, fcs, ids, sdf, xp, yp, zp, method.iso)
+    vts, fcs, nms
+end
+
+function mc_voxels!(vts, nms, fcs, ids, sdf, xp, yp, zp, iso)
+    nx, ny, nz = size(sdf)
+    h = (step(xp), step(yp), step(zp))
+    idx = Vector{Int}(undef, 12)
+
     @inbounds for xi = 1:nx-1, yi = 1:ny-1, zi = 1:nz-1
 
         iso_vals = (sdf[xi, yi, zi],
@@ -50,17 +75,69 @@ function isosurface(sdf::AbstractArray{T,3}, method::MarchingCubes, X=-1:1, Y=-1
 
         #Determine the index into the edge table which
         #tells us which vertices are inside of the surface
-        cubeindex = _get_cubeindex(iso_vals, method.iso)
+        cubeindex = _get_cubeindex(iso_vals, iso)
 
         # Cube is entirely in/out of the surface
         _no_triangles(cubeindex) && continue
 
         points = mc_vert_points(xi, yi, zi, xp, yp, zp)
+        grads = mc_vert_grads(nms, sdf, CartesianIndex(xi, yi, zi), h)
 
         # process the voxel
-        process_mc_voxel!(vts, fcs, cubeindex, points, method.iso, iso_vals)
+        mc_voxel!(vts, nms, fcs, ids, idx, cubeindex, (xi, yi, zi), (nx, ny), points, grads, iso, iso_vals)
     end
-    vts, fcs
+end
+
+# no shared vertices and no normals: the original path
+mc_voxel!(vts, ::Nothing, fcs, ::Nothing, idx, cubeindex, voxel, dims, points, grads, iso, iso_vals) =
+    process_mc_voxel!(vts, fcs, cubeindex, points, iso, iso_vals)
+
+function mc_voxel!(vts, nms, fcs, ids, idx, cubeindex, voxel, dims, points, grads, iso, iso_vals)
+    @inbounds begin
+        vert_to_add = _mc_verts[cubeindex]
+        for i = 1:12
+            vt = vert_to_add[i]
+            iszero(vt) && break
+            idx[i] = mc_vertex!(vts, nms, ids, vt, voxel, dims, points, grads, iso, iso_vals)
+        end
+
+        offsets = _mc_connectivity[_mc_eq_mapping[cubeindex]]
+        push!(fcs, (idx[3], idx[2], idx[1]))
+        for i in (1, 4, 7, 10)
+            iszero(offsets[i]) && return
+            push!(fcs, (idx[offsets[i+2]], idx[offsets[i+1]], idx[offsets[i]]))
+        end
+    end
+end
+
+# the index of edge `vt`'s vertex: a new one, or the one its grid edge already has
+mc_vertex!(vts, nms, ::Nothing, vt, voxel, dims, points, grads, iso, iso_vals) =
+    push_mc_vertex!(vts, nms, vt, points, grads, iso, iso_vals)
+mc_vertex!(vts, nms, ids::Dict, vt, voxel, dims, points, grads, iso, iso_vals) =
+    get!(() -> push_mc_vertex!(vts, nms, vt, points, grads, iso, iso_vals), ids, mc_edge_id(voxel, vt, dims))
+
+function push_mc_vertex!(vts, nms, vt, points, grads, iso, iso_vals)
+    a, b = _mc_edge_list[vt]
+    push!(vts, vertex_interp(iso, points[a], points[b], iso_vals[a], iso_vals[b]))
+    push_mc_normal!(nms, vertex_interp(iso, grads[a], grads[b], iso_vals[a], iso_vals[b]))
+    length(vts)
+end
+function push_mc_vertex!(vts, nms, vt, points, ::Nothing, iso, iso_vals)
+    a, b = _mc_edge_list[vt]
+    push!(vts, vertex_interp(iso, points[a], points[b], iso_vals[a], iso_vals[b]))
+    length(vts)
+end
+
+push_mc_normal!(nms, g) = push!(nms, all(iszero, g) ? (zero(g[1]), zero(g[1]), one(g[1])) : g ./ sqrt(sum(abs2, g)))
+
+# each cube edge as (offset of its lower grid point, axis), so voxels sharing an edge give it one id
+const mc_edge_offsets = ((0, 0, 0, 1), (1, 0, 0, 2), (0, 1, 0, 1), (0, 0, 0, 2),
+                         (0, 0, 1, 1), (1, 0, 1, 2), (0, 1, 1, 1), (0, 0, 1, 2),
+                         (0, 0, 0, 3), (1, 0, 0, 3), (1, 1, 0, 3), (0, 1, 0, 3))
+
+@inline function mc_edge_id((xi, yi, zi), edge, (nx, ny))
+    dx, dy, dz, axis = mc_edge_offsets[edge]
+    axis + 3 * ((xi + dx - 1) + nx * ((yi + dy - 1) + ny * (zi + dz - 1)))
 end
 
 
@@ -118,3 +195,27 @@ function mc_vert_points(xi, yi, zi, xp, yp, zp)
         (xp[xi+1], yp[yi+1], zp[zi+1]),
         (xp[xi], yp[yi+1], zp[zi+1]))
 end
+
+# corner offsets, in `mc_vert_points` order
+const mc_corners = (CartesianIndex(0, 0, 0), CartesianIndex(1, 0, 0), CartesianIndex(1, 1, 0), CartesianIndex(0, 1, 0),
+                    CartesianIndex(0, 0, 1), CartesianIndex(1, 0, 1), CartesianIndex(1, 1, 1), CartesianIndex(0, 1, 1))
+
+mc_vert_grads(::Nothing, sdf, I, h) = nothing
+mc_vert_grads(nms, sdf, I, h) = map(c -> sdf_grad(sdf, I + c, h), mc_corners)
+
+sdf_grad(sdf, I, h) = ntuple(a -> sdf_deriv(sdf, I, a, h[a]), Val(3))
+
+# 4th-order central difference, 2nd-order next to the border, one-sided on it
+@inline function sdf_deriv(sdf, I, a, h)
+    d = CartesianIndex(ntuple(b -> Int(b == a), Val(3)))
+    i, n = I[a], size(sdf, a)
+    @inbounds if 2 < i < n - 1
+        (sdf[I-2d] - 8sdf[I-d] + 8sdf[I+d] - sdf[I+2d]) / 12h
+    elseif 1 < i < n
+        (sdf[I+d] - sdf[I-d]) / 2h
+    elseif i == 1
+        (sdf[I+d] - sdf[I]) / h
+    else
+        (sdf[I] - sdf[I-d]) / h
+    end
+end

test/runtests.jl

@@ -1,8 +1,8 @@
 using Meshing
 using Test
 using ForwardDiff
-using Statistics: mean
-using LinearAlgebra: dot, norm
+using Statistics: mean, std
+using LinearAlgebra: dot, norm, normalize
 using Random
 
 const algos = (MarchingCubes, MarchingTetrahedra)
@@ -167,7 +167,46 @@ end
 
     points, faces = isosurface(distance, MarchingTetrahedra(iso=lambda))
 
-    @test length(points) == 3466
-    @test length(faces) == 6928
+    # watertight: every edge in exactly two faces. Not counts or topology: `MersenneTwister(0)`
+    # draws a different stream since Julia 1.11, and the noise decides both.
+    edges = Dict{Tuple{Int,Int},Int}()
+    for f in faces, e in ((f[1], f[2]), (f[2], f[3]), (f[3], f[1]))
+        edges[minmax(e...)] = get(edges, minmax(e...), 0) + 1
+    end
+    @test length(faces) > 1000
+    @test all(==(2), values(edges))
+end
+
+@testset "MarchingCubes reduceverts" begin
+    @test MarchingCubes(0.5) == MarchingCubes(iso=0.5)
+    p, f = isosurface(sphere_sdf, MarchingCubes())
+    pr, fr = isosurface(sphere_sdf, MarchingCubes(reduceverts=true))
+    # face for face the same triangles (a shared edge is interpolated from either end)
+    @test length(fr) == length(f)
+    @test maximum(i -> maximum(k -> maximum(abs, pr[fr[i][k]] .- p[f[i][k]]), 1:3), eachindex(f)) < 1e-12
+    # one vertex per crossed grid edge
+    crossed(d) = count(I -> checkbounds(Bool, sphere_sdf, I + d) && (sphere_sdf[I] < 0) != (sphere_sdf[I+d] < 0),
+                       CartesianIndices(sphere_sdf))
+    @test length(pr) == sum(crossed, (CartesianIndex(1, 0, 0), CartesianIndex(0, 1, 0), CartesianIndex(0, 0, 1)))
 end
 
+@testset "isosurface_normals" begin
+    for reduceverts in (false, true)
+        method = MarchingCubes(; reduceverts)
+        p, f, n = isosurface_normals(sphere_sdf, method)
+        @test (p, f) == isosurface(sphere_sdf, method)
+        @test all(v -> norm(v) ≈ 1, n)
+        @test minimum(i -> dot(collect(n[i]), normalize(collect(p[i]))), eachindex(p)) > 0.999
+    end
+    # z spacing doubled: an ellipsoid, whose normal at (x, y, z) is along (x, y, z/4)
+    p, f, n = isosurface_normals(sphere_sdf, MarchingCubes(), -1:1, -1:1, -2:2)
+    @test minimum(i -> dot(collect(n[i]), normalize([p[i][1], p[i][2], p[i][3] / 4])), eachindex(p)) > 0.999
+end
+
+@testset "smooth_sdf" begin
+    @test smooth_sdf(sphere_sdf; sigma=0.1) == sphere_sdf
+    noisy = sphere_sdf .+ 0.02 .* randn(MersenneTwister(0), size(sphere_sdf))
+    s = smooth_sdf(noisy; sigma=1.0)
+    @test std(s .- sphere_sdf) < std(noisy .- sphere_sdf) / 2
+    @test mean(s) ≈ mean(noisy) atol=1e-3
+end