Skip to content

Commit 3969a98

Browse files
sjkellysjkelly
committed
add naive surfacenet implmentation
1 parent 70e7368 commit 3969a98

File tree

3 files changed

+257
-3
lines changed

3 files changed

+257
-3
lines changed

‎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: 229 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,229 @@
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, dims,eps,scale,origin)
50+
51+
vertices = Point{3,Float64}[]
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(Int32),R[3]*2)
63+
64+
v = [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+
94+
(grid[1] < 0.0) && (mask |= 0x01)
95+
(grid[2] < 0.0) && (mask |= 0x02)
96+
(grid[3] < 0.0) && (mask |= 0x04)
97+
(grid[4] < 0.0) && (mask |= 0x08)
98+
(grid[5] < 0.0) && (mask |= 0x10)
99+
(grid[6] < 0.0) && (mask |= 0x20)
100+
(grid[7] < 0.0) && (mask |= 0x40)
101+
(grid[8] < 0.0) && (mask |= 0x80)
102+
103+
# Check for early termination if cell does not intersect boundary
104+
if mask == 0x00 || mask == 0xff
105+
x[1] += 1
106+
n += 1
107+
m += 1
108+
continue
109+
end
110+
111+
#Sum up edge intersections
112+
edge_mask = sn_edge_table[mask+1]
113+
v[1] = 0.0
114+
v[2] = 0.0
115+
v[3] = 0.0
116+
e_count = 0
117+
118+
#For every edge of the cube...
119+
@inbounds for i=0:11
120+
121+
#Use edge mask to check if it is crossed
122+
if (edge_mask & (1<<i)) == 0
123+
continue
124+
end
125+
126+
#If it did, increment number of edge crossings
127+
e_count += 1
128+
129+
#Now find the point of intersection
130+
e0 = cube_edges[(i<<1)+1] #Unpack vertices
131+
e1 = cube_edges[(i<<1)+2]
132+
g0 = grid[e0+1] #Unpack grid values
133+
g1 = grid[e1+1]
134+
t = g0 - g1 #Compute point of intersection
135+
if abs(t) > eps
136+
t = g0 / t
137+
else
138+
continue
139+
end
140+
141+
#Interpolate vertices and add up intersections (this can be done without multiplying)
142+
j = 0
143+
k = 1
144+
while j<3
145+
a = e0 & k
146+
b = e1 & k;
147+
if a != b
148+
v[j+1] += (a != 0 ? 1.0 - t : t)
149+
else
150+
v[j+1] += (a != 0 ? 1.0 : 0)
151+
end
152+
j+=1
153+
k<<=1
154+
end
155+
end # edge check
156+
157+
#Now we just average the edge intersections and add them to coordinate
158+
s = 1.0 / e_count
159+
for i=1:3
160+
@inbounds v[i] = (x[i] + s * v[i]) * scale[i] + origin[i]
161+
end
162+
163+
#Add vertex to buffer, store pointer to vertex index in buffer
164+
buffer[m+1] = length(vertices)
165+
push!(vertices, Point{3,Float64}(v[1],v[2],v[3]));
166+
167+
#Now we need to add faces together, to do this we just loop over 3 basis components
168+
for i=0:2
169+
#The first three entries of the edge_mask count the crossings along the edge
170+
if (edge_mask & (1<<i)) == 0
171+
continue
172+
end
173+
174+
# i = axes we are point along. iu, iv = orthogonal axes
175+
iu = (i+1)%3
176+
iv = (i+2)%3
177+
178+
#If we are on a boundary, skip it
179+
if (x[iu+1] == 0 || x[iv+1] == 0)
180+
continue
181+
end
182+
183+
#Otherwise, look up adjacent edges in buffer
184+
du = R[iu+1]
185+
dv = R[iv+1]
186+
187+
#Remember to flip orientation depending on the sign of the corner.
188+
if (mask & 0x01) != 0x00
189+
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));
190+
else
191+
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));
192+
end
193+
end
194+
x[1] += 1
195+
n += 1
196+
m += 1
197+
end
198+
x[2] += 1
199+
n += 1
200+
m += 2
201+
end
202+
x[3] += 1
203+
n+=dims[1]
204+
buf_no = xor(buf_no,1)
205+
R[3]=-R[3]
206+
end
207+
#All done! Return the result
208+
209+
vertices, faces # faces are quads, indexed to vertices
210+
end
211+
212+
struct NaiveSurfaceNets{T} <: AbstractMeshingAlgorithm
213+
eps::T
214+
end
215+
216+
NaiveSurfaceNets() = NaiveSurfaceNets(1e-6)
217+
218+
function (::Type{MT})(sdf::SignedDistanceField, method::NaiveSurfaceNets) where {MT <: AbstractMesh}
219+
bounds = sdf.bounds
220+
orig = origin(bounds)
221+
w = widths(bounds)
222+
scale = w ./ Point(size(sdf) .- 1) # subtract 1 because an SDF with N points per side has N-1 cells
223+
vts, fcs = surface_nets(reshape(sdf.data,(1,length(sdf.data))),
224+
size(sdf.data),
225+
method.eps,
226+
scale,
227+
orig)
228+
MT(vts, fcs)::MT
229+
end

‎test/runtests.jl

Lines changed: 25 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,31 @@ 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+
@test for vt in vertices(sphere)
25+
d = sqrt(sum(vt .^ 2))
26+
if 0.001 < (d-1) < 0.001
27+
continue
28+
else
29+
return false
30+
end
31+
true
32+
end
33+
end
34+
35+
1136
@testset "noisy spheres" begin
1237
# Produce a level set function that is a noisy version of the distance from
1338
# the origin (such that level sets are noisy spheres).
@@ -146,5 +171,3 @@ using LinearAlgebra: dot, norm
146171
end
147172
end
148173
end
149-
150-

0 commit comments

Comments
 (0)