Skip to content

Commit 2c0995c

Browse files
rdeitsrdeits
authored
Merge pull request #20 from JuliaGeometry/sjk/surfacenets2
Add naive surfacenet implementation
2 parents 70e7368 + 5212972 commit 2c0995c

File tree

4 files changed

+327
-3
lines changed

4 files changed

+327
-3
lines changed

‎README.md

Lines changed: 46 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,52 @@ This package provides meshing algorithms for use on distance fields.
88
Including:
99
* [Marching Tetrahedra](https://en.wikipedia.org/wiki/Marching_tetrahedra)
1010
* [Marching Cubes](https://en.wikipedia.org/wiki/Marching_cubes)
11+
* [Naive Surface Nets](https://0fps.net/2012/07/12/smooth-voxel-terrain-part-2/)
12+
13+
## Interface
14+
15+
This package is tightly integrated with [GeometryTypes.jl](https://github.com/JuliaGeometry/GeometryTypes.jl).
16+
17+
All algorithms operate on `SignedDistanceField` and output a concrete `AbstractMesh`. For example:
18+
19+
```
20+
using Meshing
21+
using GeometryTypes
22+
using LinearAlgebra: dot, norm
23+
using FileIO
24+
25+
# generate an SDF of a sphere
26+
sdf_sphere = SignedDistanceField(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.))) do v
27+
sqrt(sum(dot(v,v))) - 1 # sphere
28+
end
29+
30+
m = GLNormalMesh(sdf_sphere, MarchingCubes())
31+
32+
save("sphere.ply",m)
33+
```
34+
35+
The general API is ``(::Type{MT})(sdf::SignedDistanceField, method::AbstractMeshingAlgorithm) where {MT <: AbstractMesh}``
36+
37+
For a full listing of concrete `AbstractMesh` types see [GeometryTypes.jl mesh documentation](http://juliageometry.github.io/GeometryTypes.jl/latest/types.html#Meshes-1).
38+
39+
### Meshing Algorithms
40+
41+
Three meshing algorithms exist:
42+
* `MarchingCubes()`
43+
* `MarchingTetrahedra()``
44+
* `NaiveSurfaceNets()`
45+
46+
Each takes an optional `iso` and `eps` parameter, e.g. `MarchingCubes(0.0,1e-6)`.
47+
48+
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.
49+
50+
Below is a comparison of the algorithms:
51+
52+
| Algorithm | Accurate | Manifold | Performance Penalty | Face Type |
53+
|--------------------|----------|----------|---------------------|-----------|
54+
| MarchingCubes | Yes | No | ~4x | Triangle |
55+
| MarchingTetrahedra | Yes | Yes | ~12x | Triangle |
56+
| NaiveSurfaceNets | No | No | 1x | Quad |
1157

1258
## Credits
1359

‎src/Meshing.jl

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -6,9 +6,11 @@ abstract type AbstractMeshingAlgorithm end
66

77
include("marching_tetrahedra.jl")
88
include("marching_cubes.jl")
9+
include("surface_nets.jl")
910

1011
export marching_cubes,
1112
MarchingCubes,
12-
MarchingTetrahedra
13+
MarchingTetrahedra,
14+
NaiveSurfaceNets
1315

1416
end # module

‎src/surface_nets.jl

Lines changed: 237 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,237 @@
1+
#
2+
# SurfaceNets in Julia
3+
#
4+
# Ported from the Javascript implementation by Mikola Lysenko (C) 2012
5+
# https://github.com/mikolalysenko/mikolalysenko.github.com/blob/master/Isosurface/js/surfacenets.js
6+
# MIT License
7+
#
8+
# Based on: S.F. Gibson, "Constrained Elastic Surface Nets". (1998) MERL Tech Report.
9+
#
10+
11+
12+
13+
# Precompute edge table, like Paul Bourke does.
14+
# This saves a bit of time when computing the centroid of each boundary cell
15+
const cube_edges = ( 0, 1, 0, 2, 0, 4, 1, 3, 1, 5, 2, 3,
16+
2, 6, 3, 7, 4, 5, 4, 6, 5, 7, 6, 7 )
17+
const sn_edge_table = [0, 7, 25, 30, 98, 101, 123, 124, 168, 175, 177, 182, 202,
18+
205, 211, 212, 772, 771, 797, 794, 870, 865, 895, 888,
19+
940, 939, 949, 946, 974, 969, 983, 976, 1296, 1303, 1289,
20+
1294, 1394, 1397, 1387, 1388, 1464, 1471, 1441, 1446,
21+
1498, 1501, 1475, 1476, 1556, 1555, 1549, 1546, 1654,
22+
1649, 1647, 1640, 1724, 1723, 1701, 1698, 1758, 1753,
23+
1735, 1728, 2624, 2631, 2649, 2654, 2594, 2597, 2619,
24+
2620, 2792, 2799, 2801, 2806, 2698, 2701, 2707, 2708,
25+
2372, 2371, 2397, 2394, 2342, 2337, 2367, 2360, 2540,
26+
2539, 2549, 2546, 2446, 2441, 2455, 2448, 3920, 3927,
27+
3913, 3918, 3890, 3893, 3883, 3884, 4088, 4095, 4065,
28+
4070, 3994, 3997, 3971, 3972, 3156, 3155, 3149, 3146,
29+
3126, 3121, 3119, 3112, 3324, 3323, 3301, 3298, 3230,
30+
3225, 3207, 3200, 3200, 3207, 3225, 3230, 3298, 3301,
31+
3323, 3324, 3112, 3119, 3121, 3126, 3146, 3149, 3155,
32+
3156, 3972, 3971, 3997, 3994, 4070, 4065, 4095, 4088,
33+
3884, 3883, 3893, 3890, 3918, 3913, 3927, 3920, 2448,
34+
2455, 2441, 2446, 2546, 2549, 2539, 2540, 2360, 2367,
35+
2337, 2342, 2394, 2397, 2371, 2372, 2708, 2707, 2701,
36+
2698, 2806, 2801, 2799, 2792, 2620, 2619, 2597, 2594,
37+
2654, 2649, 2631, 2624, 1728, 1735, 1753, 1758, 1698,
38+
1701, 1723, 1724, 1640, 1647, 1649, 1654, 1546, 1549,
39+
1555, 1556, 1476, 1475, 1501, 1498, 1446, 1441, 1471,
40+
1464, 1388, 1387, 1397, 1394, 1294, 1289, 1303, 1296,
41+
976, 983, 969, 974, 946, 949, 939, 940, 888, 895, 865,
42+
870, 794, 797, 771, 772, 212, 211, 205, 202, 182, 177,
43+
175, 168, 124, 123, 101, 98, 30, 25, 7, 0]
44+
45+
"""
46+
Generate a mesh using naive surface nets.
47+
This takes the center of mass of the voxel as the vertex for each cube.
48+
"""
49+
function surface_nets(data::Vector{T}, dims,eps,scale,origin) where {T}
50+
51+
vertices = Point{3,T}[]
52+
faces = Face{4,Int}[]
53+
54+
sizehint!(vertices,ceil(Int,maximum(dims)^2/2))
55+
sizehint!(faces,ceil(Int,maximum(dims)^2/2))
56+
57+
n = 0
58+
x = [0,0,0]
59+
R = Array{Int}([1, (dims[1]+1), (dims[1]+1)*(dims[2]+1)])
60+
buf_no = 1
61+
62+
buffer = fill(zero(Int),R[3]*2)
63+
64+
v = Vector{T}([0.0,0.0,0.0])
65+
66+
#March over the voxel grid
67+
x[3] = 0
68+
@inbounds while x[3]<dims[3]-1
69+
70+
# m is the pointer into the buffer we are going to use.
71+
# This is slightly obtuse because javascript does not have good support for packed data structures, so we must use typed arrays :(
72+
# The contents of the buffer will be the indices of the vertices on the previous x/y slice of the volume
73+
m = 1 + (dims[1]+1) * (1 + buf_no * (dims[2]+1))
74+
75+
x[2]=0
76+
@inbounds while x[2]<dims[2]-1
77+
78+
x[1]=0
79+
@inbounds while x[1] < dims[1]-1
80+
81+
# Read in 8 field values around this vertex and store them in an array
82+
# Also calculate 8-bit mask, like in marching cubes, so we can speed up sign checks later
83+
mask = 0x00
84+
@inbounds grid = (data[n+1],
85+
data[n+2],
86+
data[n+dims[1]+1],
87+
data[n+dims[1]+2],
88+
data[n+dims[1]*2+1 + dims[1]*(dims[2]-2)],
89+
data[n+dims[1]*2+2 + dims[1]*(dims[2]-2)],
90+
data[n+dims[1]*3+1 + dims[1]*(dims[2]-2)],
91+
data[n+dims[1]*3+2 + dims[1]*(dims[2]-2)])
92+
93+
signbit(grid[1]) && (mask |= 0x01)
94+
signbit(grid[2]) && (mask |= 0x02)
95+
signbit(grid[3]) && (mask |= 0x04)
96+
signbit(grid[4]) && (mask |= 0x08)
97+
signbit(grid[5]) && (mask |= 0x10)
98+
signbit(grid[6]) && (mask |= 0x20)
99+
signbit(grid[7]) && (mask |= 0x40)
100+
signbit(grid[8]) && (mask |= 0x80)
101+
102+
# Check for early termination if cell does not intersect boundary
103+
if mask == 0x00 || mask == 0xff
104+
x[1] += 1
105+
n += 1
106+
m += 1
107+
continue
108+
end
109+
110+
#Sum up edge intersections
111+
edge_mask = sn_edge_table[mask+1]
112+
v[1] = 0.0
113+
v[2] = 0.0
114+
v[3] = 0.0
115+
e_count = 0
116+
117+
#For every edge of the cube...
118+
@inbounds for i=0:11
119+
120+
#Use edge mask to check if it is crossed
121+
if (edge_mask & (1<<i)) == 0
122+
continue
123+
end
124+
125+
#If it did, increment number of edge crossings
126+
e_count += 1
127+
128+
#Now find the point of intersection
129+
e0 = cube_edges[(i<<1)+1] #Unpack vertices
130+
e1 = cube_edges[(i<<1)+2]
131+
g0 = grid[e0+1] #Unpack grid values
132+
g1 = grid[e1+1]
133+
t = g0 - g1 #Compute point of intersection
134+
if abs(t) > eps
135+
t = g0 / t
136+
else
137+
continue
138+
end
139+
140+
#Interpolate vertices and add up intersections (this can be done without multiplying)
141+
k = 1
142+
for j = 1:3
143+
a = e0 & k
144+
b = e1 & k
145+
(a != 0) && (v[j] += 1.0)
146+
if a != b
147+
v[j] += (a != 0 ? - t : t)
148+
end
149+
k<<=1
150+
end
151+
end # edge check
152+
153+
#Now we just average the edge intersections and add them to coordinate
154+
s = 1.0 / e_count
155+
for i=1:3
156+
@inbounds v[i] = (x[i] + s * v[i]) * scale[i] + origin[i]
157+
end
158+
159+
#Add vertex to buffer, store pointer to vertex index in buffer
160+
buffer[m+1] = length(vertices)
161+
push!(vertices, Point{3,T}(v[1],v[2],v[3]))
162+
163+
#Now we need to add faces together, to do this we just loop over 3 basis components
164+
for i=0:2
165+
#The first three entries of the edge_mask count the crossings along the edge
166+
if (edge_mask & (1<<i)) == 0
167+
continue
168+
end
169+
170+
# i = axes we are point along. iu, iv = orthogonal axes
171+
iu = (i+1)%3
172+
iv = (i+2)%3
173+
174+
#If we are on a boundary, skip it
175+
if (x[iu+1] == 0 || x[iv+1] == 0)
176+
continue
177+
end
178+
179+
#Otherwise, look up adjacent edges in buffer
180+
du = R[iu+1]
181+
dv = R[iv+1]
182+
183+
#Remember to flip orientation depending on the sign of the corner.
184+
if (mask & 0x01) != 0x00
185+
push!(faces,Face{4,Int}(buffer[m+1]+1, buffer[m-du+1]+1, buffer[m-du-dv+1]+1, buffer[m-dv+1]+1));
186+
else
187+
push!(faces,Face{4,Int}(buffer[m+1]+1, buffer[m-dv+1]+1, buffer[m-du-dv+1]+1, buffer[m-du+1]+1));
188+
end
189+
end
190+
x[1] += 1
191+
n += 1
192+
m += 1
193+
end
194+
x[2] += 1
195+
n += 1
196+
m += 2
197+
end
198+
x[3] += 1
199+
n+=dims[1]
200+
buf_no = xor(buf_no,1)
201+
R[3]=-R[3]
202+
end
203+
#All done! Return the result
204+
205+
vertices, faces # faces are quads, indexed to vertices
206+
end
207+
208+
struct NaiveSurfaceNets{T} <: AbstractMeshingAlgorithm
209+
iso::T
210+
eps::T
211+
end
212+
213+
NaiveSurfaceNets(iso::T1=0.0, eps::T2=1e-3) where {T1, T2} = NaiveSurfaceNets{promote_type(T1, T2)}(iso, eps)
214+
215+
function (::Type{MT})(sdf::SignedDistanceField, method::NaiveSurfaceNets) where {MT <: AbstractMesh}
216+
bounds = sdf.bounds
217+
orig = origin(bounds)
218+
w = widths(bounds)
219+
scale = w ./ Point(size(sdf) .- 1) # subtract 1 because an SDF with N points per side has N-1 cells
220+
221+
d = vec(sdf.data)
222+
223+
# Run iso surface additions here as to not
224+
# penalize surface net inner loops, and we are using a copy anyway
225+
if method.iso != 0.0
226+
for i = eachindex(d)
227+
d[i] -= method.iso
228+
end
229+
end
230+
231+
vts, fcs = surface_nets(d,
232+
size(sdf.data),
233+
method.eps,
234+
scale,
235+
orig)
236+
MT(vts, fcs)::MT
237+
end

‎test/runtests.jl

Lines changed: 41 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,22 @@ using LinearAlgebra: dot, norm
88

99

1010
@testset "meshing" begin
11+
@testset "surface nets" begin
12+
sdf_sphere = SignedDistanceField(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.))) do v
13+
sqrt(sum(dot(v,v))) - 1 # sphere
14+
end
15+
sdf_torus = SignedDistanceField(HyperRectangle(Vec(-2,-2,-2.),Vec(4,4,4.)), 0.05) do v
16+
(sqrt(v[1]^2+v[2]^2)-0.5)^2 + v[3]^2 - 0.25
17+
end
18+
sphere = HomogenousMesh(sdf_sphere, NaiveSurfaceNets())
19+
torus = HomogenousMesh(sdf_torus, NaiveSurfaceNets())
20+
@test length(vertices(sphere)) == 1832
21+
@test length(vertices(torus)) == 5532
22+
@test length(faces(sphere)) == 1830
23+
@test length(faces(torus)) == 5532
24+
end
25+
26+
1127
@testset "noisy spheres" begin
1228
# Produce a level set function that is a noisy version of the distance from
1329
# the origin (such that level sets are noisy spheres).
@@ -77,6 +93,16 @@ using LinearAlgebra: dot, norm
7793
@test maximum(vertices(mesh)) ≈ [0.5, 0.5, 0.5]
7894
@test minimum(vertices(mesh)) ≈ [-0.5, -0.5, -0.5]
7995
end
96+
# Naive Surface Nets has no accuracy guarantee, and is a weighted sum
97+
# so a larger tolerance is needed for this one. In addition,
98+
# quad -> triangle conversion is not functioning correctly
99+
# see: https://github.com/JuliaGeometry/GeometryTypes.jl/issues/169
100+
mesh = @inferred GLNormalMesh(sdf, NaiveSurfaceNets(0.5))
101+
# should be centered on the origin
102+
@test mean(vertices(mesh)) ≈ [0, 0, 0] atol=0.15*resolution
103+
# and should be symmetric about the origin
104+
@test maximum(vertices(mesh)) ≈ [0.5, 0.5, 0.5] atol=0.2
105+
@test minimum(vertices(mesh)) ≈ [-0.5, -0.5, -0.5] atol=0.2
80106
end
81107

82108
@testset "AbstractMeshingAlgorithm interface" begin
@@ -96,6 +122,10 @@ using LinearAlgebra: dot, norm
96122
@test_nowarn GLNormalMesh(sdf.data, MarchingTetrahedra(0.5))
97123
@inferred GLNormalMesh(sdf.data, MarchingTetrahedra(0.5))
98124
end
125+
@testset "naive surface nets" begin
126+
@test_nowarn GLNormalMesh(sdf, NaiveSurfaceNets())
127+
@inferred GLNormalMesh(sdf, NaiveSurfaceNets())
128+
end
99129
end
100130

101131
@testset "mixed types" begin
@@ -144,7 +174,16 @@ using LinearAlgebra: dot, norm
144174
@inferred(Meshing.marchingTetrahedra(Float32.(data), Float64(iso), Float16(eps), Int32))
145175
@inferred(Meshing.marchingTetrahedra(Float64.(data), Float32(iso), Float64(eps), Int64))
146176
end
177+
@testset "Float16" begin
178+
sdf_torus = SignedDistanceField(HyperRectangle(Vec{3,Float16}(-2,-2,-2.),
179+
Vec{3,Float16}(4,4,4.)),
180+
0.1, Float16) do v
181+
(sqrt(v[1]^2+v[2]^2)-0.5)^2 + v[3]^2 - 0.25
182+
end
183+
@test typeof(HomogenousMesh(sdf_torus,NaiveSurfaceNets())) ==
184+
PlainMesh{Float16,Face{4,Int}}
185+
m2 = HomogenousMesh(sdf_torus,MarchingTetrahedra())
186+
m3 = HomogenousMesh(sdf_torus,MarchingCubes())
187+
end
147188
end
148189
end
149-
150-

0 commit comments

Comments
 (0)