From 64ccb81dff2203647a09a7b88e58e0c718363394 Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Thu, 4 Jul 2024 09:19:43 -0400 Subject: [PATCH 1/3] remove function sampling API --- bench/Project.toml | 4 - bench/iso_masjk.jl | 27 ---- src/Meshing.jl | 10 +- src/adaptive.jl | 316 ------------------------------------- src/isosurface.jl | 35 ++++ src/marching_cubes.jl | 41 +---- src/marching_tetrahedra.jl | 37 ----- src/roots.jl | 73 --------- test/runtests.jl | 238 +++++++++++----------------- 9 files changed, 133 insertions(+), 648 deletions(-) delete mode 100644 bench/Project.toml delete mode 100644 bench/iso_masjk.jl delete mode 100644 src/adaptive.jl create mode 100644 src/isosurface.jl delete mode 100644 src/roots.jl diff --git a/bench/Project.toml b/bench/Project.toml deleted file mode 100644 index 2dececf..0000000 --- a/bench/Project.toml +++ /dev/null @@ -1,4 +0,0 @@ -[deps] -BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" -JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b" -Meshing = "e6723b4c-ebff-59f1-b4b7-d97aa5274f73" diff --git a/bench/iso_masjk.jl b/bench/iso_masjk.jl deleted file mode 100644 index e05939b..0000000 --- a/bench/iso_masjk.jl +++ /dev/null @@ -1,27 +0,0 @@ -using BenchmarkTools - -@inline function _get_cubeindex(iso_vals, iso) - cubeindex = iso_vals[1] < iso ? 0x01 : 0x00 - iso_vals[2] < iso && (cubeindex |= 0x02) - iso_vals[3] < iso && (cubeindex |= 0x04) - iso_vals[4] < iso && (cubeindex |= 0x08) - iso_vals[5] < iso && (cubeindex |= 0x10) - iso_vals[6] < iso && (cubeindex |= 0x20) - iso_vals[7] < iso && (cubeindex |= 0x40) - iso_vals[8] < iso && (cubeindex |= 0x80) - cubeindex -end - -a = rand(8) - -@benchmark _get_cubeindex(a, 0.5) - -@inline function _get_cubeindex2(iso_vals, iso) - cubeindex = 0 - for i in 1:8 - cubeindex |= (iso_vals[i] < iso ? 1 : 0) << (i - 1) - end - cubeindex -end - -@benchmark _get_cubeindex2(a, 0.5) diff --git a/src/Meshing.jl b/src/Meshing.jl index dbde984..530ffe6 100644 --- a/src/Meshing.jl +++ b/src/Meshing.jl @@ -1,18 +1,10 @@ module Meshing -""" - _DEFAULT_SAMPLES = (24,24,24) - -Global default sampling count for functions. -""" -const _DEFAULT_SAMPLES = (24,24,24) - include("algorithmtypes.jl") include("common.jl") include("marching_tetrahedra.jl") include("marching_cubes.jl") -include("roots.jl") -include("adaptive.jl") +include("isosurface.jl") export isosurface, MarchingCubes, diff --git a/src/adaptive.jl b/src/adaptive.jl deleted file mode 100644 index f2e6619..0000000 --- a/src/adaptive.jl +++ /dev/null @@ -1,316 +0,0 @@ -# Code derived from Robin Deits' RegionTrees and AdaptiveDistanceFields -# The RegionTrees.jl package is licensed under the MIT "Expat" License -# The AdaptiveDistanceFields.jl package is licensed under the MIT "Expat" License - -# this is designed to avoid allocation of a full octtree and generate a mesh -# while sampling - -struct HyperRectangle{N,T} - origin::NTuple{N,T} - widths::NTuple{N,T} -end - - -function vertices(h::HyperRectangle, ::Type{SV}) where SV - z = zero(eltype(h.origin)) - @inbounds begin - o = SV(h.origin[1],h.origin[2],h.origin[3]) - w = SV(h.widths[1],h.widths[2],h.widths[3]) - (o, - o.+SV(w[1],z,z), - o.+SV(w[1],w[2],z), - o.+SV(z,w[2],z), - o.+SV(z,z,w[3]), - o.+SV(w[1],z,w[3]), - o.+w, - o.+SV(z,w[2],w[3])) - end -end - -function vertices_mt(h::HyperRectangle, ::Type{SV}) where SV - z = zero(eltype(SV)) - @inbounds begin - o = SV(h.origin[1],h.origin[2],h.origin[3]) - w = SV(h.widths[1],h.widths[2],h.widths[3]) - (o, - o.+SV( z,w[2], z), - o.+SV(w[1],w[2], z), - o.+SV(w[1], z, z), - o.+SV( z, z,w[3]), - o.+SV( z,w[2],w[3]), - o.+w, - o.+SV(w[1], z,w[3])) - end -end - -function interpolate_mt(c) - @inbounds begin - (sum(c)*0.125, # center - (c[1]+c[2]+c[3]+c[4])*0.25, - (c[1]+c[4]+c[5]+c[8])*0.25, - (c[1]+c[2]+c[5]+c[6])*0.25, - (c[5]+c[6]+c[7]+c[8])*0.25, - (c[2]+c[3]+c[6]+c[7])*0.25, - (c[3]+c[4]+c[7]+c[8])*0.25) - end -end - -function face_center_vertices(h::HyperRectangle{N,T}) where{N,T} - SV = SVector{3,T} - z = zero(T) - @inbounds begin - o = SV(h.origin[1],h.origin[2],h.origin[3]) - w = SV(h.widths[1],h.widths[2],h.widths[3]) - hw = w.*0.5 - c = center(h) - (SV(c[1],c[2],c[3]), - o.+SV(hw[1],hw[2], z), - o.+SV(hw[1], z,hw[3]), - o.+SV(z ,hw[2],hw[3]), - o.+SV(hw[1],hw[2], w[3]), - o.+SV(hw[1], w[2],hw[3]), - o.+SV( w[1],hw[2],hw[3])) - end -end - -center(rect::HyperRectangle) = rect.origin + 0.5 * rect.widths - -function octsplit(h::HyperRectangle) - ET = eltype(h.origin) - VT = Vec{3,ET} - z = zero(ET) - hw = h.widths - nw = h.widths.*ET(0.5) - no = h.origin - @inbounds begin - (HyperRectangle(no, nw), - HyperRectangle(no.+nw, nw), - HyperRectangle(no.+VT(nw[1],z,z), nw), - HyperRectangle(no.+VT(z,nw[2],z), nw), - HyperRectangle(no.+VT(z,z,nw[3]), nw), - HyperRectangle(no.+VT(nw[1],nw[2],z), nw), - HyperRectangle(no.+VT(nw[1],z,nw[3]), nw), - HyperRectangle(no.+VT(z,nw[2],nw[3]), nw)) - end -end - -function get_iso_vals(f, val_store, points::NTuple{N,T}) where {N, T} - @inbounds begin - return ntuple(N) do i - get_iso_vals(f,val_store,points[i]) - end - end -end - -function get_iso_vals(f, val_store, pt) - # tuples are faster to hash than SVec - key_tup = pt.data # tuple in SVectors/GeometryTypes, probably shoudl make more generic - if haskey(val_store,key_tup) - val_store[key_tup] - else - val_store[key_tup] = f(pt)::valtype(val_store) - end -end - -@inline function _get_interpindex(iso_vals, iso) - cubeindex = iso_vals[1] < iso ? 0x01 : 0x00 - iso_vals[2] < iso && (cubeindex |= 0x02) - iso_vals[3] < iso && (cubeindex |= 0x04) - iso_vals[4] < iso && (cubeindex |= 0x08) - iso_vals[5] < iso && (cubeindex |= 0x10) - iso_vals[6] < iso && (cubeindex |= 0x20) - iso_vals[7] < iso && (cubeindex |= 0x40) - cubeindex -end - -""" - _get_cubeindex_pos(iso_vals, iso) - -given `iso_vals` and iso, return an 8 bit value corresponding -to each corner of a cube. In each bit position, -0 indicates in the isosurface and 1 indicates outside the surface, -where the sign convention indicates positive inside the surface -""" -@inline function _get_interpindex_pos(iso_vals, iso) - cubeindex = iso_vals[1] > iso ? 0x01 : 0x00 - iso_vals[2] > iso && (cubeindex |= 0x02) - iso_vals[3] > iso && (cubeindex |= 0x04) - iso_vals[4] > iso && (cubeindex |= 0x08) - iso_vals[5] > iso && (cubeindex |= 0x10) - iso_vals[6] > iso && (cubeindex |= 0x20) - iso_vals[7] > iso && (cubeindex |= 0x40) - cubeindex -end - - -function isosurface(f::Function, method::AdaptiveMarchingCubes, ::Type{VertType}=SVector{3,Float64}, ::Type{FaceType}=SVector{3, Int}; - origin=VertType(-1,-1,-1), widths=VertType(2,2,2)) where {VertType, FaceType} - - ET = eltype(VertType) - - # arrays for vertices and faces - vts = VertType[] - fcs = FaceType[] - - # refinement queue - refinement_queue = HyperRectangle{3,ET}[HyperRectangle{3,ET}(origin,widths)] - - val_store = Dict{NTuple{3,ET},ET}(); - - @inbounds while true - - if isempty(refinement_queue) - break - end - - cell = pop!(refinement_queue) - points = vertices(cell, VertType) - - iso_vals = get_iso_vals(f,val_store,points) - - #Determine the index into the edge table which - #tells us which vertices are inside of the surface - cubeindex = method.insidepositive ? _get_cubeindex_pos(iso_vals, method.iso) : _get_cubeindex(iso_vals, method.iso) - #interpindex = method.insidepositive ? _get_interpindex_pos(iso_vals, method.iso) : _get_interpindex(iso_vals, method.iso) - - value_interp = sum(iso_vals)*0.125 - value_true = get_iso_vals(f, val_store, center(cell)) - - if (cubeindex == 0xff && value_true < 0) || (iszero(cubeindex) && value_true > 0) - continue - elseif minimum(cell.widths) > method.atol && !isapprox(value_interp, value_true, rtol=method.rtol, atol=method.atol) - append!(refinement_queue, octsplit(cell)) - else - # Find the vertices where the surface intersects the cube - # The underlying space is non-linear so there will be error otherwise - vertlist = find_vertices_interp(points, iso_vals, cubeindex, method.iso, method.eps) - - # Create the triangle - method.reduceverts && _mc_unique_triangles!(vts, fcs, vertlist, cubeindex, FaceType) - !method.reduceverts && _mc_create_triangles!(vts, fcs, vertlist, cubeindex, FaceType) - end - end - vts,fcs -end - - -@inline function vertPos(f, e, width, origin, ::Type{VertType}) where {VertType} - T = eltype(eltype(VertType)) - - ixs = voxEdgeCrnrs[e] - c1 = voxCrnrPos(VertType)[ixs[1]].*width .+ origin - c2 = voxCrnrPos(VertType)[ixs[2]].*width .+ origin - a = brent(f, c1, c2.-c1) # find root using brents method - - c1 + a.*(c2.-c1) -end - - -@inline function getVertId(f, e, width, vals, iso::Real, origin, vtsAry::Vector, vertex_store, eps::Real) - - VertType = eltype(vtsAry) - - # calculate vert position - v = vertPos(f, e, width, origin, VertType) - vt_key = v.data - if haskey(vertex_store, vt_key) - return vertex_store[vt_key] - else - push!(vtsAry, v) - return vertex_store[vt_key] = length(vtsAry) - end -end - -slivers(x,y,z) = x == y || y == z || x == z - -""" - procVox(vals, iso::Real, x, y, z, nx, ny, - vts::Dict, vtsAry::Vector, fcs::Vector, - eps::Real) - -Processes a voxel, adding any new vertices and faces to the given -containers as necessary. -""" -function procVox(f, vals, iso::Real, width, origin, vtsAry::Vector, vertex_store, fcs::Vector, - eps::Real, cubeindex) - VertType = eltype(vtsAry) - FaceType = eltype(fcs) - # check each sub-tetrahedron in the voxel - @inbounds for i = 1:6 - tIx = tetIx(i, cubeindex) - (tIx == 0x00 || tIx == 0x0f) && continue - - e = tetTri[tIx] - v1 = getVertId(f, voxEdgeId(i, e[1]), width, vals, iso, origin, vtsAry, vertex_store, eps) - v2 = getVertId(f, voxEdgeId(i, e[2]), width, vals, iso, origin, vtsAry, vertex_store, eps) - v3 = getVertId(f, voxEdgeId(i, e[3]), width, vals, iso, origin, vtsAry, vertex_store, eps) - - # add the face to the list - !slivers(v1,v2,v3) && push!(fcs, FaceType(v1,v2,v3)) - - # bail if there are no more faces - iszero(e[4]) && continue - v4 = getVertId(f, voxEdgeId(i, e[4]), width, vals, iso, origin, vtsAry, vertex_store, eps) - v5 = getVertId(f, voxEdgeId(i, e[5]), width, vals, iso, origin, vtsAry, vertex_store, eps) - v6 = getVertId(f, voxEdgeId(i, e[6]), width, vals, iso, origin, vtsAry, vertex_store, eps) - !slivers(v4,v5,v6) && push!(fcs, FaceType(v4,v5,v6)) - - end -end - -function isosurface(f::Function, method::AdaptiveMarchingTetrahedra, ::Type{VertType}=SVector{3,Float64}, ::Type{FaceType}=SVector{3, Int}; - origin=VertType(-1,-1,-1), widths=VertType(2,2,2)) where {VertType, FaceType} - - ET = eltype(VertType) - - # arrays for vertices and faces - vts = VertType[] - fcs = FaceType[] - - # refinement queue - # initialize with depth two - refinement_queue = HyperRectangle{3,ET}[] - for h in octsplit(HyperRectangle{3,ET}(origin,widths)) - append!(refinement_queue,octsplit(h)) - end - - val_store = Dict{NTuple{3,ET},ET}(); - vertex_store = Dict{NTuple{3,ET},ET}(); - - @inbounds while true - - if isempty(refinement_queue) - break - end - - cell = pop!(refinement_queue) - points = vertices_mt(cell, VertType) - interp_points = face_center_vertices(cell) - - iso_vals = get_iso_vals(f,val_store,points) - true_vals = get_iso_vals(f,val_store,interp_points) - - #Determine the index into the edge table which - #tells us which vertices are inside of the surface - cubeindex = method.insidepositive ? _get_cubeindex_pos(iso_vals, method.iso) : _get_cubeindex(iso_vals, method.iso) - interpindex = method.insidepositive ? _get_interpindex_pos(true_vals, method.iso) : _get_interpindex(true_vals, method.iso) - - if (cubeindex == 0xff && interpindex == 0x7f) || (iszero(cubeindex) && iszero(interpindex)) - continue - else - value_interp = interpolate_mt(iso_vals) - accurate = true - for i = 1:7 - if !isapprox(true_vals[i], value_interp[i], rtol=method.rtol, atol=method.atol) - append!(refinement_queue, octsplit(cell)) - accurate = false - break - end - end - if accurate - procVox(f, iso_vals, method.iso, cell.widths, cell.origin, vts, vertex_store, fcs, method.eps, cubeindex) - end - end - end - vts,fcs -end diff --git a/src/isosurface.jl b/src/isosurface.jl new file mode 100644 index 0000000..c643633 --- /dev/null +++ b/src/isosurface.jl @@ -0,0 +1,35 @@ + +# +# General isosurface docstring +# + +@doc """ + isosurface(V) + isosurface(V, method::AbstractMeshingAlgorithm, X, Y, Z) + +`isosurface` is the general interface to all isosurface extraction algorithms. + +Returns: (Vector, Vector{FaceType}) + +Defaults: +`method` must be an instance of an `AbstractMeshingAlgorithm`, e.g.: +- MarchingCubes() +- MarchingTetrahedra() + +If `isosurface` is called without a specified algorithm, it will default to MarchingCubes. + +If a subtype of `AbstractArray` is specified, the mesh will by default be centered at the origin between +(-1,1) in each axis. + + +See also: +- [MarchingCubes](@ref) +- [MarchingTetrahedra](@ref) + +""" +function isosurface(A, args...) + isosurface(A, MarchingCubes(), args...) +end + + + diff --git a/src/marching_cubes.jl b/src/marching_cubes.jl index 4794048..3ec4308 100644 --- a/src/marching_cubes.jl +++ b/src/marching_cubes.jl @@ -2,6 +2,7 @@ #Look up Table include("lut/mc.jl") + function isosurface(sdf::AbstractArray{T,3}, method::MarchingCubes, X=-1:1, Y=-1:1, Z=-1:1) where {T} nx, ny, nz = size(sdf) @@ -42,46 +43,6 @@ function isosurface(sdf::AbstractArray{T,3}, method::MarchingCubes, X=-1:1, Y=-1 end -function isosurface(f::Function, method::MarchingCubes, X=-1:1, Y=-1:1, Z=-1:1; samples::NTuple{3,T}=_DEFAULT_SAMPLES) where {T<:Integer} - - nx, ny, nz = samples[1], samples[2], samples[3] - - # 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}[] - - xp = LinRange(first(X), last(X), nx) - yp = LinRange(first(Y), last(Y), ny) - zp = LinRange(first(Z), last(Z), nz) - - @inbounds for xi = 1:nx-1, yi = 1:ny-1, zi = 1:nz-1 - - points = mc_vert_points(xi, yi, zi, xp, yp, zp) - - iso_vals = (f(points[1]), - f(points[2]), - f(points[3]), - f(points[4]), - f(points[5]), - f(points[6]), - f(points[7]), - f(points[8])) - - #Determine the index into the edge table which - #tells us which vertices are inside of the surface - cubeindex = _get_cubeindex(iso_vals, method.iso) - - # Cube is entirely in/out of the surface - _no_triangles(cubeindex) && continue - - process_mc_voxel!(vts, fcs, cubeindex, points, method.iso, iso_vals) - end - vts, fcs -end - - function process_mc_voxel!(vts, fcs, cubeindex, points, iso, iso_vals) fct = length(vts) diff --git a/src/marching_tetrahedra.jl b/src/marching_tetrahedra.jl index faffc4a..0b685de 100644 --- a/src/marching_tetrahedra.jl +++ b/src/marching_tetrahedra.jl @@ -162,43 +162,6 @@ function isosurface(sdf::AbstractArray{T, 3}, method::MarchingTetrahedra, X=-1:1 vtsAry,fcs end -function isosurface(f::F, method::MarchingTetrahedra, X=-1:1, Y=-1:1, Z=-1:1; - samples::NTuple{3,T}=_DEFAULT_SAMPLES) where {F, T} - - FT = promote_type(eltype(first(X)), eltype(first(Y)), eltype(first(Z)), eltype(T), typeof(method.iso), typeof(method.eps)) - - vts = Dict{Int, Int}() - fcs = NTuple{3,Int}[] - vtsAry = NTuple{3,float(FT)}[] - - # process each voxel - nx, ny, nz = samples[1], samples[2], samples[3] - xp = LinRange(first(X), last(X), nx) - yp = LinRange(first(Y), last(Y), ny) - zp = LinRange(first(Z), last(Z), nz) - - @inbounds for i = 1:nx-1, j = 1:ny-1, k = 1:nz-1 - points = mt_vert_points(i, j, k, xp, yp, zp) - - vals = (f(points[1]), - f(points[2]), - f(points[3]), - f(points[4]), - f(points[5]), - f(points[6]), - f(points[7]), - f(points[8])) - - cubeindex = _get_cubeindex(vals, method.iso) - - _no_triangles(cubeindex) && continue - - procVox(vals, method.iso, i, j, k, nx, ny, xp, yp, zp, vts, vtsAry, fcs, method.eps, cubeindex) - end - - vtsAry,fcs -end - function mt_vert_points(i, j, k, xp, yp, zp) ((xp[i], yp[j], zp[k]), (xp[i], yp[j+1], zp[k]), diff --git a/src/roots.jl b/src/roots.jl deleted file mode 100644 index 3c3d629..0000000 --- a/src/roots.jl +++ /dev/null @@ -1,73 +0,0 @@ -# derived from Modesto Mas implementation: -# https://mmas.github.io/brent-julia -function brent(f, origin, width, xtol=1e-3, ytol=1e-3, maxiter=50) - ET = eltype(origin) - EPS = eps(ET) - - x0 = zero(ET) - x1 = one(ET) - y0 = f(x0.*width .+ origin)::ET - y1 = f(x1.*width .+ origin)::ET - if abs(y0) < abs(y1) - # Swap lower and upper bounds. - x0, x1 = x1, x0 - y0, y1 = y1, y0 - end - x2 = x0 - y2 = y0 - x3 = x2 - bisection = true - for _ in 1:maxiter - # x-tolerance. - if abs(x1-x0) < xtol - return x1 - end - - # Use inverse quadratic interpolation if f(x0)!=f(x1)!=f(x2) - # and linear interpolation (secant method) otherwise. - if abs(y0-y2) > ytol && abs(y1-y2) > ytol - x = x0*y1*y2/((y0-y1)*(y0-y2)) + - x1*y0*y2/((y1-y0)*(y1-y2)) + - x2*y0*y1/((y2-y0)*(y2-y1)) - else - x = x1 - y1 * (x1-x0)/(y1-y0) - end - - # Use bisection method if satisfies the conditions. - delta = abs(2EPS*abs(x1)) - min1 = abs(x-x1) - min2 = abs(x1-x2) - min3 = abs(x2-x3) - if (x < (3x0+x1)/4 && x > x1) || - (bisection && min1 >= min2/2) || - (!bisection && min1 >= min3/2) || - (bisection && min2 < delta) || - (!bisection && min3 < delta) - x = (x0+x1)/2 - bisection = true - else - bisection = false - end - - y = f(x.*width .+ origin)::ET - # y-tolerance. - if abs(y) < ytol - return x - end - x3 = x2 - x2 = x1 - if signbit(y0) != signbit(y) - x1 = x - y1 = y - else - x0 = x - y0 = y - end - if abs(y0) < abs(y1) - # Swap lower and upper bounds. - x0, x1 = x1, x0 - y0, y1 = y1, y0 - end - end - error("Max iteration exceeded") -end diff --git a/test/runtests.jl b/test/runtests.jl index 6073c62..6e50216 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -9,6 +9,8 @@ const algos = (MarchingCubes, MarchingTetrahedra) sphere_function(v) = sqrt(sum(dot(v, v))) - 1 torus_function(v) = (sqrt(v[1]^2 + v[2]^2) - 0.5)^2 + v[3]^2 - 0.25 +const norm_sdf =[norm((x,y,z)) for x in -1:0.01:1, y in -1:0.01:1, z in -1:0.01:1] +const sphere_sdf =[sphere_function((x,y,z)) for x in -1:0.01:1, y in -1:0.01:1, z in -1:0.01:1] @testset "_get_cubeindex" begin iso_vals1 = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8] @@ -30,21 +32,111 @@ end @test Meshing.vertex_interp(1, (0,0,0), (0,1,0), -1, 1) == (0,1,0) end +@testset "type stability" begin + # verify that we don't lose type stability just by mixing arguments + # of different types + data = randn(5, 5, 5) + iso = 0.2 + eps = 1e-4 + algo1 = MarchingTetrahedra(iso=iso, eps=eps) + algo2 = MarchingTetrahedra(iso=Float16(iso), eps=Float16(eps)) + algo3 = MarchingTetrahedra(iso=Float32(iso), eps=Float32(eps)) + @inferred(isosurface(data, algo1)) + @inferred(isosurface(Float32.(data), algo2)) + @inferred(isosurface(Float64.(data), algo3)) + + algo1 = MarchingCubes(iso=iso) + algo2 = MarchingCubes(iso=Float16(iso)) + algo3 = MarchingCubes(iso=Float32(iso)) + @inferred(isosurface(data, algo1)) + @inferred(isosurface(Float32.(data), algo2)) + @inferred(isosurface(Float64.(data), algo3)) +end + +@testset "Float16" begin + sdf_torus = [((sqrt(x^2 + y^2) - 0.5)^2 + z^2 - 0.25) for x in -2:0.05:2, y in -2:0.05:2, z in -2:0.05:2] + points_mt, faces_mt = isosurface(Float16.(sdf_torus), MarchingTetrahedra(iso=Float16(0.0), eps=Float16(1e-3)), Float16(-2):Float16(2), Float16(-2):Float16(2), Float16(-2):Float16(2)) + points_mc, faces_mc = isosurface(Float16.(sdf_torus), MarchingCubes(iso=Float16(0.0)), Float16(-2):Float16(2), Float16(-2):Float16(2), Float16(-2):Float16(2)) + + @test typeof(points_mt) == Vector{NTuple{3, Float16}} + @test typeof(faces_mt) == Vector{NTuple{3, Int}} + @test typeof(points_mc) == Vector{NTuple{3, Float16}} + @test typeof(faces_mc) == Vector{NTuple{3, Int}} +end +@testset "Integers" begin + A = rand(Int, 10,10,10) + for algo in algos + @testset "$algo" begin + p, t = isosurface(A, algo()) + @test typeof(p) <: Vector{NTuple{3, Float64}} + @test typeof(t) <: Vector{NTuple{3, Int}} + end + end +end +@testset "Defaults" begin + A = rand(10,10,10) + @test isosurface(A) == isosurface(A, MarchingCubes()) +end + +@testset "mixed types" begin + # https://github.com/JuliaGeometry/Meshing.jl/issues/9 + samples = randn(10, 10, 10) + + # Extract isosurfaces using MarchingTetrahedra + points_mt1, faces_mt1 = isosurface(samples, MarchingTetrahedra(), 0:1, 0:1, 0:1) + points_mt2, faces_mt2 = isosurface(Float32.(samples), MarchingTetrahedra(), 0:1, 0:1, 0:1) + + @test length(points_mt1) == length(points_mt2) + #for i in eachindex(points_mt1) + # @test all(points_mt1[i] .≈ points_mt2[i]) + #end + @test length(faces_mt1) == length(faces_mt2) + for i in eachindex(faces_mt1) + @test all(faces_mt1[i] .≈ faces_mt2[i]) + end +end +@testset "forward diff" begin + # Demonstrate that we can even take derivatives through the meshing algorithm + resolution = 0.1 + + for algo in (MarchingCubes, MarchingTetrahedra) + @testset "$algo" begin + let + function surface_distance_from_origin(isovalue) + # Compute the mean distance of each vertex in the isosurface + # mesh from the origin. This should return a value equal to the + # isovalue and should have a derivative of 1.0 w.r.t. the + # isovalue. This function is just meant to serve as an example + # of some quantity you might want to differentiate in a mesh, + # and has the benefit for testing of having a trivial expected + # solution. + rng = -2:0.01:2 + sdf = [norm((x,y,z)) for x in rng, y in rng, z in rng] + points, _ = isosurface(sdf, algo(iso=isovalue), rng, rng, rng) + mean(norm.(points)) + end + + isoval = 0.85 + @test surface_distance_from_origin(isoval) ≈ isoval atol=1e-2 + @test ForwardDiff.derivative(surface_distance_from_origin, isoval) ≈ 1 atol=1e-2 + end + end + end +end + @testset "respect origin" begin # verify that when we construct a mesh, that mesh: # a) respects the origin of the SDF # b) is correctly spaced so that symmetric functions create symmetric meshes - f = x -> norm(x) - resolution = 0.1 for algorithm in (MarchingCubes(iso=0.5), MarchingTetrahedra(iso=0.5)) @testset "$(typeof(algorithm))" begin # Extract isosurface using a function - points, faces = isosurface(f, algorithm, samples=(50, 50, 50)) + points, faces = isosurface(norm_sdf, algorithm) # should be centered on the origin - @test mean(collect.(points)) ≈ [0, 0, 0] atol=0.15*resolution + @test mean(collect.(points)) ≈ [0, 0, 0] atol=0.015 # and should be symmetric about the origin for i in 1:3 @@ -79,141 +171,3 @@ end @test length(faces) == 6928 end -@testset "marching cubes" begin - algo = MarchingCubes() - - # Extract isosurfaces using MarchingCubes - points_mf, faces_mf = isosurface(sphere_function, algo, samples=(21, 21, 21)) - - # Test the number of vertices and faces - @test length(points_mf) == 7320 - @test length(faces_mf) == 3656 -end -@testset "marching tetrahedra" begin - a1 = MarchingTetrahedra() - - # Extract isosurfaces using MarchingTetrahedra - points_mt1, faces_mt1 = isosurface(v -> sqrt(sum(dot(v, v))) - 1, a1, samples=(21, 21, 21)) - - # Test the number of vertices and faces - @test length(points_mt1) == 5582 - @test length(faces_mt1) == 11160 - - # When zero set is not completely within the box - points_mt1_partial, faces_mt1_partial = isosurface(v -> v[3] - 50, a1, 0:50, 0:50, 40:60, samples=(2, 2, 2)) - - # Test the number of vertices and faces - @test length(points_mt1_partial) == 9 - @test length(faces_mt1_partial) == 8 -end - -@testset "function/array" begin - - for algo in algos - @testset "$algo" begin - # Extract isosurface using a function - points_fn, faces_fn = isosurface(torus_function, algo(), -2:2, -2:2, -2:2, samples=(81, 81, 81)) - - # Extract isosurface using an array - torus_array = [torus_function((x, y, z)) for x in -2:0.05:2, y in -2:0.05:2, z in -2:0.05:2] - points_arr, faces_arr = isosurface(torus_array, algo(), -2:2, -2:2, -2:2) - - # Test that the vertices and faces are the same for both cases - @test_broken all(points_fn .≈ points_arr) - @test_broken faces_fn == faces_arr - end - end -end - -@testset "mixed types" begin - # https://github.com/JuliaGeometry/Meshing.jl/issues/9 - samples = randn(10, 10, 10) - - # Extract isosurfaces using MarchingTetrahedra - points_mt1, faces_mt1 = isosurface(samples, MarchingTetrahedra(), 0:1, 0:1, 0:1) - points_mt2, faces_mt2 = isosurface(Float32.(samples), MarchingTetrahedra(), 0:1, 0:1, 0:1) - - @test length(points_mt1) == length(points_mt2) - #for i in eachindex(points_mt1) - # @test all(points_mt1[i] .≈ points_mt2[i]) - #end - @test length(faces_mt1) == length(faces_mt2) - for i in eachindex(faces_mt1) - @test all(faces_mt1[i] .≈ faces_mt2[i]) - end - - @testset "forward diff" begin - # Demonstrate that we can even take derivatives through the meshing algorithm - resolution = 0.1 - - for algo in (MarchingCubes, MarchingTetrahedra) - @testset "$algo" begin - let - function surface_distance_from_origin(isovalue) - # Compute the mean distance of each vertex in the isosurface - # mesh from the origin. This should return a value equal to the - # isovalue and should have a derivative of 1.0 w.r.t. the - # isovalue. This function is just meant to serve as an example - # of some quantity you might want to differentiate in a mesh, - # and has the benefit for testing of having a trivial expected - # solution. - points, _ = isosurface(norm, algo(iso=isovalue), samples=(21, 21, 21)) - mean(norm.(points)) - end - - isoval = 0.85 - @test surface_distance_from_origin(isoval) ≈ isoval atol=1e-2 - @test ForwardDiff.derivative(surface_distance_from_origin, isoval) ≈ 1 atol=1e-2 - end - end - end - end - - - - @testset "type stability" begin - # verify that we don't lose type stability just by mixing arguments - # of different types - data = randn(5, 5, 5) - iso = 0.2 - eps = 1e-4 - algo1 = MarchingTetrahedra(iso=iso, eps=eps) - algo2 = MarchingTetrahedra(iso=Float16(iso), eps=Float16(eps)) - algo3 = MarchingTetrahedra(iso=Float32(iso), eps=Float32(eps)) - @inferred(isosurface(data, algo1)) - @inferred(isosurface(Float32.(data), algo2)) - @inferred(isosurface(Float64.(data), algo3)) - - algo1 = MarchingCubes(iso=iso) - algo2 = MarchingCubes(iso=Float16(iso)) - algo3 = MarchingCubes(iso=Float32(iso)) - @inferred(isosurface(data, algo1)) - @inferred(isosurface(Float32.(data), algo2)) - @inferred(isosurface(Float64.(data), algo3)) - end - - @testset "Float16" begin - sdf_torus = [((sqrt(x^2 + y^2) - 0.5)^2 + z^2 - 0.25) for x in -2:0.05:2, y in -2:0.05:2, z in -2:0.05:2] - points_mt, faces_mt = isosurface(Float16.(sdf_torus), MarchingTetrahedra(iso=Float16(0.0), eps=Float16(1e-3)), Float16(-2):Float16(2), Float16(-2):Float16(2), Float16(-2):Float16(2)) - points_mc, faces_mc = isosurface(Float16.(sdf_torus), MarchingCubes(iso=Float16(0.0)), Float16(-2):Float16(2), Float16(-2):Float16(2), Float16(-2):Float16(2)) - - @test typeof(points_mt) == Vector{NTuple{3, Float16}} - @test typeof(faces_mt) == Vector{NTuple{3, Int}} - @test typeof(points_mc) == Vector{NTuple{3, Float16}} - @test typeof(faces_mc) == Vector{NTuple{3, Int}} - end - @testset "Integers" begin - A = rand(Int, 10,10,10) - for algo in algos - @testset "$algo" begin - p, t = isosurface(A, algo()) - @test typeof(p) <: Vector{NTuple{3, Float64}} - @test typeof(t) <: Vector{NTuple{3, Int}} - end - end - end - @testset "Defaults" begin - A = rand(10,10,10) - @test isosurface(A) == isosurface(A, MarchingCubes()) - end -end \ No newline at end of file From 3c55752ba71609491a33a525ec760188bf78d155 Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Thu, 4 Jul 2024 09:31:05 -0400 Subject: [PATCH 2/3] docs updates --- docs/make.jl | 3 +- docs/src/api.md | 32 +-------------------- docs/src/examples.md | 65 ++++++++++--------------------------------- docs/src/index.md | 1 - docs/src/internals.md | 58 -------------------------------------- 5 files changed, 17 insertions(+), 142 deletions(-) delete mode 100644 docs/src/internals.md diff --git a/docs/make.jl b/docs/make.jl index e482fc6..03ca8e5 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -7,8 +7,7 @@ makedocs( modules = [Meshing], pages = ["Index" => "index.md", "API" => "api.md", - "Examples" => "examples.md", - "Internals" => "internals.md"] + "Examples" => "examples.md"] ) deploydocs( diff --git a/docs/src/api.md b/docs/src/api.md index 663d155..54ecf7e 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -22,24 +22,6 @@ A = rand(50,50,50) # 3D Matrix points,faces = isosurface(A, MarchingCubes(iso=1)) ``` -Alternatively, we can use the `isosurface` API to sample a function, avoiding allocations: - -```julia -using Meshing -using LinearAlgebra -using StaticArrays - -points, faces = isosurface(origin=SVector(-1,-1,-1.), widths = SVector(2,2,2.), samples = (40,40,40)) do v - sqrt(sum(dot(v,v))) - 1 -end - -# by default MarchingCubes() is used, but we may specify a different algorithm as follows - -points, faces = isosurface(MarchingTetrahedra(), origin=SVector(-1,-1,-1.), widths = SVector(2,2,2.), samples = (40,40,40)) do v - sqrt(sum(dot(v,v))) - 1 -end -``` - ## Isosurface `isosurface` is the common and generic API for isosurface extraction with any type of abstract vector/vertex/face type. @@ -48,30 +30,19 @@ end isosurface ``` - ## Meshing Algorithms Three meshing algorithms exist: * `MarchingCubes()` * `MarchingTetrahedra()` -* `NaiveSurfaceNets()` -Each takes optional `iso`, `eps`, and `insidepositive` parameters, e.g. `MarchingCubes(iso=0.0,eps=1e-6,insidepositive=false)`. +Each takes optional `iso` and `eps` parameters, e.g. `MarchingCubes(iso=0.0,eps=1e-6)`. Here `iso` controls the offset for the boundary detection. By default this is set to 0. `eps` is the detection tolerance for a voxel edge intersection. -`insidepositive` sets the sign convention for inside/outside the surface (default: false). Users must construct an algorithm type and use it as an argument to a GeometryTypes mesh call or `isosurface` call. -Below is a comparison of the algorithms: - -| Algorithm Type | Face Type | Unique Vertices | Performance | Interpolation | -|---------------------|-----------|-----------------|-------------|-------------------| -| Naive Surface Nets | Quad | Yes | ~1x | Voxel Edge Weight | -| Marching Cubes | Triangle | No/Partial | 1x | Linear on Edge | -| Marching Tetrahedra | Triangle | Yes | 3x | Linear on Edge | - Visual Comparison: From left: Marching Cubes, Naive Surface Nets, Marching Tetrahedra @@ -80,6 +51,5 @@ From left: Marching Cubes, Naive Surface Nets, Marching Tetrahedra ```@docs MarchingCubes MarchingTetrahedra -NaiveSurfaceNets Meshing.AbstractMeshingAlgorithm ``` diff --git a/docs/src/examples.md b/docs/src/examples.md index e1f2394..543fe10 100644 --- a/docs/src/examples.md +++ b/docs/src/examples.md @@ -5,31 +5,29 @@ The file for this example can be found here: [http://www.slicer.org/slicerWiki/images/0/00/CTA-cardio.nrrd](http://www.slicer.org/slicerWiki/images/0/00/CTA-cardio.nrrd) ```julia - +using Meshing using FileIO using NRRD -using Meshing -using MeshIO +using WGLMakie +using Downloads using GeometryBasics +nrrd = Downloads.download("http://www.slicer.org/slicerWiki/images/0/00/CTA-cardio.nrrd") + # load the file as an AxisArray -ctacardio = load("CTA-cardio.nrrd") +ctacardio = load(nrrd) # use marching cubes with isolevel at 100 -algo = MarchingCubes(iso=100, insidepositive=true) +#algo = MarchingCubes(iso=100) # use marching tetrahedra with iso at 100 -# algo = MarchingTetrahedra(iso=100, insidepositive=true) -# use Naive Surface Nets with iso at 100 -# algo = NaiveSurfaceNets(iso=100, insidepositive=true) + algo = MarchingTetrahedra(iso=100) + # generate the mesh using marching cubes -mc = Mesh(ctacardio, algo) +vts, fcs = isosurface(ctacardio, algo) -# we can call isosurface to get a vector of points and vector of faces indexing to the points -# vertices, faces = isosurface(ctacardio, algo, Point{3,Float32}, TriangleFace{Int}) +WGLMakie.mesh(vts, map(v -> GeometryBasics.TriangleFace(v...), fcs)) -# save the file as a PLY file (change extension to save as STL, OBJ, OFF) -save("ctacardio_mc.ply", mc) ``` ![cta cardio](./img/ctacardio.png) @@ -60,44 +58,11 @@ Makie.mesh(gy_mesh, color=[norm(v) for v in coordinates(gy_mesh)]) ![gyroid](./img/gyroid.png) -```julia -using Meshing -using GeometryBasics -using LinearAlgebra: dot, norm -using FileIO - -# Mesh an equation of sphere in the Axis-Aligned Bounding box starting -# at -1,-1,-1 and widths of 2,2,2 using Marching Cubes -m = GLNormalMesh(Rect(Vec(-1,-1,-1.), Vec(2,2,2.)), MarchingCubes()) do v - sqrt(sum(dot(v,v))) - 1 -end - -# save the Sphere as a PLY file -save("sphere.ply",m) -``` - -For a full listing of concrete `AbstractMesh` types see [GeometryBasics.jl mesh documentation](http://juliageometry.github.io/GeometryBasics.jl/latest/types.html#Meshes-1). - -Alternatively, we can use the `isosurface` API to sample a function: - -```julia -using Meshing -using LinearAlgebra -using StaticArrays - -points, faces = isosurface(origin=SVector(-1,-1,-1.), widths = SVector(2,2,2.), samples = (40,40,40)) do v - sqrt(sum(dot(v,v))) - 1 -end - -# by default MarchingCubes() is used, but we may specify a different algorithm as follows - -points, faces = isosurface(MarchingTetrahedra(), origin=SVector(-1,-1,-1.), widths = SVector(2,2,2.), samples = (40,40,40)) do v - sqrt(sum(dot(v,v))) - 1 -end -``` - ## Isocaps +We do not provide an equivalent to `isocaps` in Matlab, though +a similar result may be achieved by setting the boundary to a large value: + ```julia using GeometryTypes @@ -120,4 +85,4 @@ gy_mesh = GLNormalMesh(A, MarchingCubes()) import Makie using LinearAlgebra Makie.mesh(gy_mesh, color=[norm(v) for v in gy_mesh.vertices]) -``` \ No newline at end of file +``` diff --git a/docs/src/index.md b/docs/src/index.md index 5ed2183..8907943 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -6,7 +6,6 @@ Algorithms included: * [Marching Tetrahedra](https://en.wikipedia.org/wiki/Marching_tetrahedra) * [Marching Cubes](https://en.wikipedia.org/wiki/Marching_cubes) -* [Naive Surface Nets](https://0fps.net/2012/07/12/smooth-voxel-terrain-part-2/) ## What is isosurface extraction? diff --git a/docs/src/internals.md b/docs/src/internals.md deleted file mode 100644 index 4bfedca..0000000 --- a/docs/src/internals.md +++ /dev/null @@ -1,58 +0,0 @@ -# Internals - -The internals of Meshing have been optimized for performance and to be generic. -This is some brief documentation on the basic internal functions. - -## Common - -```@docs -Meshing._DEFAULT_SAMPLES -Meshing._get_cubeindex -``` - -## Marching Cubes - -```@docs -Meshing._mc_create_triangles! -Meshing._mc_unique_triangles! -Meshing.find_vertices_interp -Meshing.vertex_interp -Meshing.mc_vert_points -``` - - -## Marching Tetrahedra - -``` -Voxel corner and edge indexing conventions - - Z - | - - 5------5------6 Extra edges not drawn - /| /| ----------- - 8 | 6 | - face diagonals - / 9 / 10 - 13: 1 to 3 - 8------7------7 | - 14: 1 to 8 - | | | | - 15: 1 to 6 - | 1------1--|---2 -- Y - 16: 5 to 7 - 12 / 11 / - 17: 2 to 7 - | 4 | 2 - 18: 4 to 7 - |/ |/ - body diagonal - 4------3------3 - 19: 1 to 7 - - / - X -``` - -```@docs -Meshing._correct_vertices! -Meshing.procVox -Meshing.voxEdgeId -Meshing.getVertId -Meshing.vertPos -Meshing.vertId -Meshing.tetIx -``` - -## Surface Nets \ No newline at end of file From 8831634fbf7d9337eddce1ab3ca82df0c45db2f1 Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Fri, 5 Jul 2024 09:49:56 -0400 Subject: [PATCH 3/3] cleanup docs --- src/algorithmtypes.jl | 39 +++------------------------------------ src/marching_cubes.jl | 21 +++++++++++++++++++++ 2 files changed, 24 insertions(+), 36 deletions(-) diff --git a/src/algorithmtypes.jl b/src/algorithmtypes.jl index 02d0bf3..ac7b9c6 100644 --- a/src/algorithmtypes.jl +++ b/src/algorithmtypes.jl @@ -4,14 +4,11 @@ Abstract type to specify an algorithm for isosurface extraction. See: -- MarchingCubes -- MarchingTetrahedra -- NaiveSurfaceNets +- [MarchingCubes](@ref) +- [MarchingTetrahedra](@ref) """ abstract type AbstractMeshingAlgorithm end -abstract type AbstractAdaptiveMeshingAlgorithm end - """ MarchingCubes(;iso=0.0) @@ -31,8 +28,7 @@ end MarchingTetrahedra(;iso=0.0, eps=1e-3) Specifies the use of the Marching Tetrahedra algorithm for isosurface extraction. -This algorithm has a roughly 2x performance penalty compared to Marching Cubes, -and generates more faces. However, each vertex is guaranteed to not be repeated, +This algorithm generates more faces. However, each vertex is guaranteed to not be repeated, making this algorithm useful for topological analysis. - `iso` specifies the iso level to use for surface extraction. @@ -43,32 +39,3 @@ Base.@kwdef struct MarchingTetrahedra{T,E} <: AbstractMeshingAlgorithm eps::E = 1e-3 end -""" - MarchingCubes(iso=0.0, eps=1e-3) - MarchingCubes(;iso=0.0, eps=1e-3) - MarchingCubes(iso) - MarchingCubes(iso,eps) - -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. -- `eps` (default: 1e-3) is the tolerence around a voxel corner to ensure manifold mesh generation. -- `reduceverts` (default: true) if true will merge vertices within a voxel to reduce mesh size by around 30% and with slight performance improvement. -""" -struct AdaptiveMarchingCubes{T,E} <: AbstractAdaptiveMeshingAlgorithm - iso::T - eps::T - rtol::T - atol::T -end - - -struct AdaptiveMarchingTetrahedra{T} <: AbstractAdaptiveMeshingAlgorithm - iso::T - eps::T - rtol::T - atol::T -end diff --git a/src/marching_cubes.jl b/src/marching_cubes.jl index 3ec4308..673e578 100644 --- a/src/marching_cubes.jl +++ b/src/marching_cubes.jl @@ -2,6 +2,27 @@ #Look up Table include("lut/mc.jl") +#= +Voxel corner and edge indexing conventions + + Z + | + + 5------8------8 + /| /| + 5 | 7 | + / 9 / 12 + 6------6------7 | + | | | | + | 1------4--|---4 -- Y + 10 / 11 / + | 1 | 3 + |/ |/ + 2------2------3 + + / + X +=# function isosurface(sdf::AbstractArray{T,3}, method::MarchingCubes, X=-1:1, Y=-1:1, Z=-1:1) where {T} nx, ny, nz = size(sdf)