From 09b07011cc5e7c8d6bcec9e24c268f2b831985d7 Mon Sep 17 00:00:00 2001 From: SimonDanisch Date: Wed, 30 Sep 2026 14:21:41 +0200 Subject: [PATCH 1/4] MarchingCubes: shared vertices, gradient normals, smooth_sdf --- src/Meshing.jl | 2 + src/algorithmtypes.jl | 6 ++- src/common.jl | 31 ++++++++++++ src/marching_cubes.jl | 107 ++++++++++++++++++++++++++++++++++++++++-- test/runtests.jl | 37 ++++++++++++++- 5 files changed, 177 insertions(+), 6 deletions(-) diff --git a/src/Meshing.jl b/src/Meshing.jl index 530ffe6..dac27b0 100644 --- a/src/Meshing.jl +++ b/src/Meshing.jl @@ -7,6 +7,8 @@ include("marching_cubes.jl") include("isosurface.jl") export isosurface, + isosurface_normals, + smooth_sdf, MarchingCubes, MarchingTetrahedra diff --git a/src/algorithmtypes.jl b/src/algorithmtypes.jl index ac7b9c6..4a7dad5 100644 --- a/src/algorithmtypes.jl +++ b/src/algorithmtypes.jl @@ -11,7 +11,7 @@ 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. @@ -19,10 +19,14 @@ 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) diff --git a/src/common.jl b/src/common.jl index 5ff2bfb..8d25ec1 100644 --- a/src/common.jl +++ b/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 diff --git a/src/marching_cubes.jl b/src/marching_cubes.jl index 673e578..4fbef78 100644 --- a/src/marching_cubes.jl +++ b/src/marching_cubes.jl @@ -25,6 +25,19 @@ 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 @@ -32,11 +45,23 @@ function isosurface(sdf::AbstractArray{T,3}, method::MarchingCubes, X=-1:1, Y=-1 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 diff --git a/test/runtests.jl b/test/runtests.jl index 6e50216..4eda759 100644 --- a/test/runtests.jl +++ b/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) @@ -171,3 +171,36 @@ end @test length(faces) == 6928 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 From fdb9490b772db9d071e5bbd943ef6400ecb8b18a Mon Sep 17 00:00:00 2001 From: SimonDanisch Date: Wed, 30 Sep 2026 15:23:37 +0200 Subject: [PATCH 2/4] Test noisy spheres for watertightness, not RNG-dependent counts --- test/runtests.jl | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 4eda759..93e7d00 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -167,8 +167,14 @@ 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 From 98820f87719e55d8ee2df3ff9cb8a16d0a427ddb Mon Sep 17 00:00:00 2001 From: SimonDanisch Date: Wed, 30 Sep 2026 15:28:44 +0200 Subject: [PATCH 3/4] CI: replace the retired actions/cache@v1, test the current release --- .github/workflows/CI.yml | 20 ++++++-------------- 1 file changed, 6 insertions(+), 14 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 570f2b8..0f6cef1 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -11,6 +11,7 @@ jobs: matrix: version: - '1.9' + - '1' - 'nightly' os: - ubuntu-latest @@ -18,29 +19,20 @@ jobs: 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: | From 1120cf778aab4dd983014cbcbf41615234993429 Mon Sep 17 00:00:00 2001 From: SimonDanisch Date: Wed, 30 Sep 2026 16:50:30 +0200 Subject: [PATCH 4/4] CI: test Julia 1.10 to 1.13, require 1.10 --- .github/workflows/CI.yml | 7 ++++--- Project.toml | 2 +- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 0f6cef1..9e198f4 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -10,9 +10,10 @@ jobs: fail-fast: false matrix: version: - - '1.9' - - '1' - - 'nightly' + - '1.10' + - '1.11' + - '1.12' + - '1.13' os: - ubuntu-latest - macOS-latest diff --git a/Project.toml b/Project.toml index 6faa49a..43185e8 100644 --- a/Project.toml +++ b/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"