Skip to content

Commit e0810f2

Browse files
rdeitsrdeits
committed
reorganize into meshing algorithms and fix long-standing vertices bug
1 parent 1d24905 commit e0810f2

File tree

4 files changed

+163
-59
lines changed

4 files changed

+163
-59
lines changed

‎src/Meshing.jl

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -4,9 +4,13 @@ module Meshing
44

55
using GeometryTypes
66

7+
abstract type AbstractMeshingAlgorithm end
8+
79
include("marching_tetrahedra.jl")
810
include("marching_cubes.jl")
911

10-
export marching_cubes
12+
export marching_cubes,
13+
MarchingCubes,
14+
MarchingTetrahedra
1115

1216
end # module

‎src/marching_cubes.jl

Lines changed: 39 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -287,12 +287,17 @@ However it may generate non-manifold meshes, while Marching
287287
Tetrahedra guarentees a manifold mesh.
288288
"""
289289
function marching_cubes(sdf::SignedDistanceField{3,ST,FT},
290-
iso=0.0,
291-
MT::Type{M}=SimpleMesh{Point{3,Float64},Face{3,Int}}) where {ST,FT,M<:AbstractMesh}
290+
iso=0.0,
291+
MT::Type{M}=SimpleMesh{Point{3,Float64},Face{3,Int}},
292+
eps=0.00001) where {ST,FT,M<:AbstractMesh}
292293
nx, ny, nz = size(sdf)
293294
h = HyperRectangle(sdf)
294295
w = widths(h)
295-
s = Point{3,Float64}(w[1]/nx, w[2]/ny, w[3]/nz)
296+
orig = origin(HyperRectangle(sdf))
297+
298+
# we subtract one from the length along each axis because
299+
# an NxNxN SDF has N-1 cells on each axis
300+
s = Point{3,Float64}(w[1]/(nx-1), w[2]/(ny-1), w[3]/(nz-1))
296301

297302
# arrays for vertices and faces
298303
vts = Point{3,Float64}[]
@@ -315,63 +320,63 @@ function marching_cubes(sdf::SignedDistanceField{3,ST,FT},
315320
# Cube is entirely in/out of the surface
316321
edge_table[cubeindex] == 0 && continue
317322

318-
points = (Point{3,Float64}(xi-1,yi-1,zi-1).*s,
319-
Point{3,Float64}(xi,yi-1,zi-1).*s,
320-
Point{3,Float64}(xi,yi,zi-1).*s,
321-
Point{3,Float64}(xi-1,yi,zi-1).*s,
322-
Point{3,Float64}(xi-1,yi-1,zi).*s,
323-
Point{3,Float64}(xi,yi-1,zi).*s,
324-
Point{3,Float64}(xi,yi,zi).*s,
325-
Point{3,Float64}(xi-1,yi,zi).*s)
323+
points = (Point{3,Float64}(xi-1,yi-1,zi-1) .* s .+ orig,
324+
Point{3,Float64}(xi,yi-1,zi-1) .* s .+ orig,
325+
Point{3,Float64}(xi,yi,zi-1) .* s .+ orig,
326+
Point{3,Float64}(xi-1,yi,zi-1) .* s .+ orig,
327+
Point{3,Float64}(xi-1,yi-1,zi) .* s .+ orig,
328+
Point{3,Float64}(xi,yi-1,zi) .* s .+ orig,
329+
Point{3,Float64}(xi,yi,zi) .* s .+ orig,
330+
Point{3,Float64}(xi-1,yi,zi) .* s .+ orig)
326331

327332
# Find the vertices where the surface intersects the cube
328333
if (edge_table[cubeindex] & 1 != 0)
329334
vertlist[1] =
330-
vertex_interp(iso,points[1],points[2],sdf[xi,yi,zi],sdf[xi+1,yi,zi])
335+
vertex_interp(iso,points[1],points[2],sdf[xi,yi,zi],sdf[xi+1,yi,zi], eps)
331336
end
332337
if (edge_table[cubeindex] & 2 != 0)
333338
vertlist[2] =
334-
vertex_interp(iso,points[2],points[3],sdf[xi+1,yi,zi],sdf[xi+1,yi+1,zi])
339+
vertex_interp(iso,points[2],points[3],sdf[xi+1,yi,zi],sdf[xi+1,yi+1,zi], eps)
335340
end
336341
if (edge_table[cubeindex] & 4 != 0)
337342
vertlist[3] =
338-
vertex_interp(iso,points[3],points[4],sdf[xi+1,yi+1,zi],sdf[xi,yi+1,zi])
343+
vertex_interp(iso,points[3],points[4],sdf[xi+1,yi+1,zi],sdf[xi,yi+1,zi], eps)
339344
end
340345
if (edge_table[cubeindex] & 8 != 0)
341346
vertlist[4] =
342-
vertex_interp(iso,points[4],points[1],sdf[xi,yi+1,zi],sdf[xi,yi,zi])
347+
vertex_interp(iso,points[4],points[1],sdf[xi,yi+1,zi],sdf[xi,yi,zi], eps)
343348
end
344349
if (edge_table[cubeindex] & 16 != 0)
345350
vertlist[5] =
346-
vertex_interp(iso,points[5],points[6],sdf[xi,yi,zi+1],sdf[xi+1,yi,zi+1])
351+
vertex_interp(iso,points[5],points[6],sdf[xi,yi,zi+1],sdf[xi+1,yi,zi+1], eps)
347352
end
348353
if (edge_table[cubeindex] & 32 != 0)
349354
vertlist[6] =
350-
vertex_interp(iso,points[6],points[7],sdf[xi+1,yi,zi+1],sdf[xi+1,yi+1,zi+1])
355+
vertex_interp(iso,points[6],points[7],sdf[xi+1,yi,zi+1],sdf[xi+1,yi+1,zi+1], eps)
351356
end
352357
if (edge_table[cubeindex] & 64 != 0)
353358
vertlist[7] =
354-
vertex_interp(iso,points[7],points[8],sdf[xi+1,yi+1,zi+1],sdf[xi,yi+1,zi+1])
359+
vertex_interp(iso,points[7],points[8],sdf[xi+1,yi+1,zi+1],sdf[xi,yi+1,zi+1], eps)
355360
end
356361
if (edge_table[cubeindex] & 128 != 0)
357362
vertlist[8] =
358-
vertex_interp(iso,points[8],points[5],sdf[xi,yi+1,zi+1],sdf[xi,yi,zi+1])
363+
vertex_interp(iso,points[8],points[5],sdf[xi,yi+1,zi+1],sdf[xi,yi,zi+1], eps)
359364
end
360365
if (edge_table[cubeindex] & 256 != 0)
361366
vertlist[9] =
362-
vertex_interp(iso,points[1],points[5],sdf[xi,yi,zi],sdf[xi,yi,zi+1])
367+
vertex_interp(iso,points[1],points[5],sdf[xi,yi,zi],sdf[xi,yi,zi+1], eps)
363368
end
364369
if (edge_table[cubeindex] & 512 != 0)
365370
vertlist[10] =
366-
vertex_interp(iso,points[2],points[6],sdf[xi+1,yi,zi],sdf[xi+1,yi,zi+1])
371+
vertex_interp(iso,points[2],points[6],sdf[xi+1,yi,zi],sdf[xi+1,yi,zi+1], eps)
367372
end
368373
if (edge_table[cubeindex] & 1024 != 0)
369374
vertlist[11] =
370-
vertex_interp(iso,points[3],points[7],sdf[xi+1,yi+1,zi],sdf[xi+1,yi+1,zi+1])
375+
vertex_interp(iso,points[3],points[7],sdf[xi+1,yi+1,zi],sdf[xi+1,yi+1,zi+1], eps)
371376
end
372377
if (edge_table[cubeindex] & 2048 != 0)
373378
vertlist[12] =
374-
vertex_interp(iso,points[4],points[8],sdf[xi,yi+1,zi],sdf[xi,yi+1,zi+1])
379+
vertex_interp(iso,points[4],points[8],sdf[xi,yi+1,zi],sdf[xi,yi+1,zi+1], eps)
375380
end
376381

377382
# Create the triangle
@@ -399,3 +404,14 @@ function vertex_interp(iso, p1, p2, valp1, valp2, eps = 0.00001)
399404

400405
return p
401406
end
407+
408+
struct MarchingCubes{T} <: AbstractMeshingAlgorithm
409+
iso::T
410+
eps::T
411+
end
412+
413+
MarchingCubes(iso::T1=0.0, eps::T2=1e-3) where {T1, T2} = MarchingCubes{promote_type(T1, T2)}(iso, eps)
414+
415+
function (::Type{MT})(df::SignedDistanceField, method::MarchingCubes) where {MT <: AbstractMesh}
416+
marching_cubes(df, method.iso, MT, method.eps)
417+
end

‎src/marching_tetrahedra.jl

Lines changed: 35 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -292,16 +292,51 @@ end
292292

293293
isosurface(lsf,isoval) = isosurface(lsf,isoval, convert(eltype(lsf), 0.001))
294294

295+
"""
296+
The marchingTetrahedra function returns vertices on the (1-based) indices of the
297+
SDF's data, ignoring its actual bounds. This function adjusts the vertices in
298+
place so that they correspond to points within the SDF bounds.
299+
"""
300+
function _correct_vertices!(vts, sdf::SignedDistanceField)
301+
bounds = HyperRectangle(sdf)
302+
orig = origin(bounds)
303+
w = widths(bounds)
304+
s = w ./ Point(size(sdf) .- 1) # subtract 1 because an SDF with N points per side has N-1 cells
305+
for i in eachindex(vts)
306+
vts[i] = (vts[i] .- 1) .* s .+ orig # subtract 1 to fix 1-indexing
307+
end
308+
end
309+
295310

296311
function (::Type{MT})(volume::Array{T, 3}, iso_val::Real, eps_val=0.001) where {MT <: AbstractMesh, T}
312+
Base.depwarn("(mesh type)(volume::Array, iso_val, eps_val) is deprecated. Please use: (mesh type)(volume, MarchingTetrahedra(iso_val, eps_val))", :Meshing_MT_array_eps)
297313
iso_val = convert(T, iso_val)
298314
eps_val = convert(T, eps_val)
299315
vts, fcs = isosurface(volume, iso_val, eps_val)
300316
MT(vts, fcs)
301317
end
302318

303319
function (::Type{MT})(df::SignedDistanceField, eps_val=0.001) where MT <: AbstractMesh
320+
Base.depwarn("(mesh type)(sdf::SignedDistanceField, eps_val) is deprecated. Please use: (mesh type)(sdf, MarchingTetrahedra(0, eps))", :Meshing_MT_sdf_eps)
304321
vts, fcs = isosurface(df.data, 0.0, eps_val)
322+
_correct_vertices!(vts, df)
305323
MT(vts, fcs)
306324
end
307325

326+
struct MarchingTetrahedra{T} <: AbstractMeshingAlgorithm
327+
iso::T
328+
eps::T
329+
end
330+
331+
MarchingTetrahedra(iso::T1=0.0, eps::T2=1e-3) where {T1, T2} = MarchingTetrahedra{promote_type(T1, T2)}(iso, eps)
332+
333+
function (::Type{MT})(sdf::SignedDistanceField, method::MarchingTetrahedra)::MT where {MT <: AbstractMesh}
334+
vts, fcs = isosurface(sdf.data, method.iso, method.eps)
335+
_correct_vertices!(vts, sdf)
336+
MT(vts, fcs)
337+
end
338+
339+
function (::Type{MT}){MT <: AbstractMesh, T}(volume::Array{T, 3}, method::MarchingTetrahedra)::MT
340+
vts, fcs = isosurface(volume, convert(T, method.iso), convert(T, method.eps))
341+
MT(vts, fcs)
342+
end

‎test/runtests.jl

Lines changed: 84 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -2,46 +2,95 @@ using Meshing
22
using Base.Test
33
using GeometryTypes
44

5-
# Produce a level set function that is a noisy version of the distance from
6-
# the origin (such that level sets are noisy spheres).
7-
#
8-
# The noise should exercise marching tetrahedra's ability to produce a water-
9-
# tight surface in all cases (unlike standard marching cubes).
10-
#
11-
N = 10
12-
sigma = 1.0
13-
distance = Float32[ sqrt(Float32(i*i+j*j+k*k)) for i = -N:N, j = -N:N, k = -N:N ]
14-
distance = distance + sigma*rand(2*N+1,2*N+1,2*N+1)
15-
16-
# Extract an isosurface.
17-
#
18-
lambda = N-2*sigma # isovalue
19-
20-
msh = HomogenousMesh(distance,lambda)
21-
22-
s2 = SignedDistanceField(HyperRectangle(Vec(0,0,0.),Vec(1,1,1.))) do v
23-
sqrt(sum(dot(v,v))) - 1 # sphere
24-
end
255

26-
msh = HomogenousMesh(s2)
27-
@test length(vertices(msh)) == 973
28-
@test length(faces(msh)) == 1830
6+
@testset "meshing" begin
7+
@testset "noisy spheres" begin
8+
# Produce a level set function that is a noisy version of the distance from
9+
# the origin (such that level sets are noisy spheres).
10+
#
11+
# The noise should exercise marching tetrahedra's ability to produce a water-
12+
# tight surface in all cases (unlike standard marching cubes).
13+
#
14+
N = 10
15+
sigma = 1.0
16+
distance = Float32[ sqrt(Float32(i*i+j*j+k*k)) for i = -N:N, j = -N:N, k = -N:N ]
17+
distance = distance + sigma*rand(2*N+1,2*N+1,2*N+1)
18+
19+
# Extract an isosurface.
20+
#
21+
lambda = N-2*sigma # isovalue
22+
23+
msh = HomogenousMesh(distance, MarchingTetrahedra(lambda))
24+
25+
s2 = SignedDistanceField(HyperRectangle(Vec(0,0,0.),Vec(1,1,1.))) do v
26+
sqrt(sum(dot(v,v))) - 1 # sphere
27+
end
28+
29+
msh = HomogenousMesh(s2, MarchingTetrahedra())
30+
@test length(vertices(msh)) == 973
31+
@test length(faces(msh)) == 1830
32+
end
33+
34+
@testset "vertex interpolation" begin
35+
@test Meshing.vertex_interp(0, Vec(0,0,0), Vec(0,1,0), -1, 1) == Vec(0,0.5,0)
36+
@test Meshing.vertex_interp(-1, Vec(0,0,0), Vec(0,1,0), -1, 1) == Vec(0,0,0)
37+
@test Meshing.vertex_interp(1, Vec(0,0,0), Vec(0,1,0), -1, 1) == Vec(0,1,0)
38+
end
39+
40+
@testset "marching cubes" begin
41+
sdf = SignedDistanceField(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.))) do v
42+
sqrt(sum(dot(v,v))) - 1 # sphere
43+
end
44+
45+
m = marching_cubes(sdf,0)
46+
m2 = marching_cubes(sdf)
47+
@test length(vertices(m)) == 10968
48+
@test length(faces(m)) == 3656
49+
@test m == m2
50+
51+
end
52+
@testset "respect origin" begin
53+
# verify that when we construct a mesh, that mesh:
54+
# a) respects the origin of the SDF
55+
# b) is correctly spaced so that symmetric functions create symmetric meshes
56+
f = x -> norm(x)
57+
bounds = HyperRectangle(Vec(-1, -1, -1), Vec(2, 2, 2))
58+
resolution = 0.1
59+
sdf = SignedDistanceField(f, bounds, resolution)
60+
61+
for algorithm in (MarchingCubes(0.5), MarchingTetrahedra(0.5))
62+
mesh = GLNormalMesh(sdf, algorithm)
63+
# should be centered on the origin
64+
@test mean(vertices(mesh)) ≈ [0, 0, 0] atol=0.15*resolution
65+
# and should be symmetric about the origin
66+
@test maximum(vertices(mesh)) ≈ [0.5, 0.5, 0.5]
67+
@test minimum(vertices(mesh)) ≈ [-0.5, -0.5, -0.5]
68+
end
69+
end
70+
71+
@testset "AbstractMeshingAlgorithm interface" begin
72+
f = x -> norm(x) - 0.5
73+
bounds = HyperRectangle(Vec(-1, -1, -1), Vec(2, 2, 2))
74+
resolution = 0.1
75+
sdf = SignedDistanceField(f, bounds, resolution)
2976

30-
# Vertex interpolation
31-
@test Meshing.vertex_interp(0, Vec(0,0,0), Vec(0,1,0), -1, 1) == Vec(0,0.5,0)
32-
@test Meshing.vertex_interp(-1, Vec(0,0,0), Vec(0,1,0), -1, 1) == Vec(0,0,0)
33-
@test Meshing.vertex_interp(1, Vec(0,0,0), Vec(0,1,0), -1, 1) == Vec(0,1,0)
77+
@testset "marching cubes" begin
78+
m1 = @test_nowarn marching_cubes(sdf, 0.0, GLNormalMesh)
79+
m2 = @test_nowarn GLNormalMesh(sdf, MarchingCubes())
80+
@test vertices(m1) == vertices(m2)
81+
@test faces(m1) == faces(m2)
82+
end
3483

35-
# marching cubes
36-
sdf = SignedDistanceField(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.))) do v
37-
sqrt(sum(dot(v,v))) - 1 # sphere
84+
@testset "marching tetrahedra" begin
85+
m1 = @test_warn "is deprecated" GLNormalMesh(sdf)
86+
m2 = @test_warn "is deprecated" GLNormalMesh(sdf, 1e-3)
87+
m3 = @test_nowarn GLNormalMesh(sdf, MarchingTetrahedra())
88+
@test vertices(m1) == vertices(m2) == vertices(m3)
89+
@test faces(m1) == faces(m2) == faces(m3)
90+
end
91+
end
3892
end
3993

40-
m = marching_cubes(sdf,0)
41-
m2 = marching_cubes(sdf)
42-
@test length(vertices(m)) == 10968
43-
@test length(faces(m)) == 3656
44-
@test m == m2
4594

4695
if "--profile" in ARGS
4796
HomogenousMesh(s2)

0 commit comments

Comments
 (0)