[SN] non-sdf SN implementation

Author: sjkelly

src/surface_nets.jl

@@ -205,6 +205,175 @@ function surface_nets(data::Vector{T}, dims,eps,scale,origin) where {T}
     vertices, faces # faces are quads, indexed to vertices
 end
 
+"""
+Generate a mesh using naive surface nets.
+This takes the center of mass of the voxel as the vertex for each cube.
+"""
+function surface_nets(f::Function, dims::NTuple{3,Int},eps,scale,origin)
+
+    # TODO
+    T = Float64
+
+    vertices = Point{3,T}[]
+    faces = Face{4,Int}[]
+
+    sizehint!(vertices,ceil(Int,maximum(dims)^2/2))
+    sizehint!(faces,ceil(Int,maximum(dims)^2/2))
+
+    n = 0
+    x = [0,0,0]
+    R = Array{Int}([1, (dims[1]+1), (dims[1]+1)*(dims[2]+1)])
+    buf_no = 1
+
+    buffer = fill(zero(Int),R[3]*2)
+
+    v = Vector{T}([0.0,0.0,0.0])
+
+    #March over the voxel grid
+    x[3] = 0
+    @inbounds while x[3]<dims[3]-1
+
+        # m is the pointer into the buffer we are going to use.
+        # This is slightly obtuse because javascript does not have good support for packed data structures, so we must use typed arrays :(
+        # The contents of the buffer will be the indices of the vertices on the previous x/y slice of the volume
+        m = 1 + (dims[1]+1) * (1 + buf_no * (dims[2]+1))
+
+        x[2]=0
+        @inbounds while x[2]<dims[2]-1
+
+            x[1]=0
+            @inbounds while x[1] < dims[1]-1
+
+                # Read in 8 field values around this vertex and store them in an array
+                points = (Point{3,Float64}(x[1],x[2],x[3]).* scale + origin,
+                          Point{3,Float64}(x[1]+1,x[2],x[3]).* scale + origin,
+                          Point{3,Float64}(x[1],x[2]+1,x[3]).* scale + origin,
+                          Point{3,Float64}(x[1]+1,x[2]+1,x[3]).* scale + origin,
+                          Point{3,Float64}(x[1],x[2],x[3]+1).* scale + origin,
+                          Point{3,Float64}(x[1]+1,x[2],x[3]+1).* scale + origin,
+                          Point{3,Float64}(x[1],x[2]+1,x[3]+1).* scale + origin,
+                          Point{3,Float64}(x[1]+1,x[2]+1,x[3]+1).* scale + origin)
+
+                @inbounds grid = map(f, points)
+
+                # Also calculate 8-bit mask, like in marching cubes, so we can speed up sign checks later
+                mask = 0x00
+                signbit(grid[1]) && (mask |= 0x01)
+                signbit(grid[2]) && (mask |= 0x02)
+                signbit(grid[3]) && (mask |= 0x04)
+                signbit(grid[4]) && (mask |= 0x08)
+                signbit(grid[5]) && (mask |= 0x10)
+                signbit(grid[6]) && (mask |= 0x20)
+                signbit(grid[7]) && (mask |= 0x40)
+                signbit(grid[8]) && (mask |= 0x80)
+
+                # Check for early termination if cell does not intersect boundary
+                if mask == 0x00 || mask == 0xff
+                    x[1] += 1
+                    n += 1
+                    m += 1
+                    continue
+                end
+
+                #Sum up edge intersections
+                edge_mask = sn_edge_table[mask+1]
+                v[1] = 0.0
+                v[2] = 0.0
+                v[3] = 0.0
+                e_count = 0
+
+                #For every edge of the cube...
+                @inbounds for i=0:11
+
+                    #Use edge mask to check if it is crossed
+                    if (edge_mask & (1<<i)) == 0
+                        continue
+                    end
+
+                    #If it did, increment number of edge crossings
+                    e_count += 1
+
+                    #Now find the point of intersection
+                    e0 = cube_edges[(i<<1)+1]       #Unpack vertices
+                    e1 = cube_edges[(i<<1)+2]
+                    g0 = grid[e0+1]                 #Unpack grid values
+                    g1 = grid[e1+1]
+                    t  = g0 - g1                 #Compute point of intersection
+                    if abs(t) > eps
+                      t = g0 / t
+                    else
+                      continue
+                    end
+
+                    #Interpolate vertices and add up intersections (this can be done without multiplying)
+                    k = 1
+                    for j = 1:3
+                        a = e0 & k
+                        b = e1 & k
+                        (a != 0) && (v[j] += 1.0)
+                        if a != b
+                            v[j] += (a != 0 ? - t : t)
+                        end
+                        k<<=1
+                    end
+                end # edge check
+
+                #Now we just average the edge intersections and add them to coordinate
+                s = 1.0 / e_count
+                for i=1:3
+                    @inbounds v[i] = (x[i] + s * v[i])# * scale[i] + origin[i]
+                end
+
+                #Add vertex to buffer, store pointer to vertex index in buffer
+                buffer[m+1] = length(vertices)
+                push!(vertices, Point{3,T}(v[1],v[2],v[3]))
+
+                #Now we need to add faces together, to do this we just loop over 3 basis components
+                for i=0:2
+                    #The first three entries of the edge_mask count the crossings along the edge
+                    if (edge_mask & (1<<i)) == 0
+                        continue
+                    end
+
+                    # i = axes we are point along.  iu, iv = orthogonal axes
+                    iu = (i+1)%3
+                    iv = (i+2)%3
+
+                    #If we are on a boundary, skip it
+                    if (x[iu+1] == 0 || x[iv+1] == 0)
+                        continue
+                    end
+
+                    #Otherwise, look up adjacent edges in buffer
+                    du = R[iu+1]
+                    dv = R[iv+1]
+
+                    #Remember to flip orientation depending on the sign of the corner.
+                    if (mask & 0x01) != 0x00
+                        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));
+                    else
+                        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));
+                    end
+                end
+                x[1] += 1
+                n += 1
+                m += 1
+            end
+            x[2] += 1
+            n += 1
+            m += 2
+        end
+        x[3] += 1
+        n+=dims[1]
+        buf_no = xor(buf_no,1)
+        R[3]=-R[3]
+    end
+    #All done!  Return the result
+
+    vertices, faces # faces are quads, indexed to vertices
+end
+
+
 struct NaiveSurfaceNets{T} <: AbstractMeshingAlgorithm
     iso::T
     eps::T
@@ -235,3 +404,18 @@ function (::Type{MT})(sdf::SignedDistanceField, method::NaiveSurfaceNets) where
                             orig)
     MT(vts, fcs)::MT
 end
+
+function (::Type{MT})(f::Function, bounds::HyperRectangle, size::NTuple{3,Int}, method::NaiveSurfaceNets) where {MT <: AbstractMesh}
+    orig = origin(bounds)
+    w = widths(bounds)
+    scale = w ./ Point(size .- 1)  # subtract 1 because an SDF with N points per side has N-1 cells
+
+    # TODO ISO val
+
+    vts, fcs = surface_nets(f,
+                            size,
+                            method.eps,
+                            scale,
+                            orig)
+    MT(vts, fcs)::MT
+end

test/runtests.jl

@@ -9,18 +9,32 @@ using LinearAlgebra: dot, norm
 
 @testset "meshing" begin
     @testset "surface nets" begin
-          sdf_sphere = SignedDistanceField(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.))) do v
-              sqrt(sum(dot(v,v))) - 1 # sphere
-          end
-          sdf_torus = SignedDistanceField(HyperRectangle(Vec(-2,-2,-2.),Vec(4,4,4.)), 0.05) do v
-              (sqrt(v[1]^2+v[2]^2)-0.5)^2 + v[3]^2 - 0.25
-          end
-          sphere = HomogenousMesh(sdf_sphere, NaiveSurfaceNets())
-          torus = HomogenousMesh(sdf_torus, NaiveSurfaceNets())
-          @test length(vertices(sphere)) == 1832
-          @test length(vertices(torus)) == 5532
-          @test length(faces(sphere)) == 1830
-          @test length(faces(torus)) == 5532
+        sdf_sphere = SignedDistanceField(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.))) do v
+            sqrt(sum(dot(v,v))) - 1 # sphere
+        end
+        sdf_torus = SignedDistanceField(HyperRectangle(Vec(-2,-2,-2.),Vec(4,4,4.)), 0.05) do v
+            (sqrt(v[1]^2+v[2]^2)-0.5)^2 + v[3]^2 - 0.25
+        end
+
+        snf_torus = HomogenousMesh(HyperRectangle(Vec(-2,-2,-2.),Vec(4,4,4.)), (81,81,81), NaiveSurfaceNets()) do v
+            (sqrt(v[1]^2+v[2]^2)-0.5)^2 + v[3]^2 - 0.25
+        end
+
+        snf_sphere = HomogenousMesh(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)), (21,21,21), NaiveSurfaceNets()) do v
+            sqrt(sum(dot(v,v))) - 1 # sphere
+        end
+
+        sphere = HomogenousMesh(sdf_sphere, NaiveSurfaceNets())
+        torus = HomogenousMesh(sdf_torus, NaiveSurfaceNets())
+        @test length(vertices(sphere)) == 1832
+        @test length(vertices(torus)) == 5532
+        @test length(faces(sphere)) == 1830
+        @test length(faces(torus)) == 5532
+
+        @test length(vertices(sphere)) == length(vertices(snf_sphere))
+        @test length(vertices(torus)) == length(vertices(snf_torus))
+        @test length(faces(sphere)) == length(faces(snf_sphere))
+        @test length(faces(torus)) == length(faces(snf_torus))
     end
 
 
@@ -69,20 +83,17 @@ using LinearAlgebra: dot, norm
             sqrt(sum(dot(v,v))) - 1 # sphere
         end
 
+        mf = marching_cubes(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)),(21,21,21)) do v
+            sqrt(sum(dot(v,v))) - 1 # sphere
+        end
+
         m = marching_cubes(sdf,0)
         m2 = marching_cubes(sdf)
         @test length(vertices(m)) == 10968
         @test length(faces(m)) == 3656
         @test m == m2
-
-    end
-
-    @testset "marching cubes function" begin
-        m = marching_cubes(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)),(21,21,21)) do v
-            sqrt(sum(dot(v,v))) - 1 # sphere
-        end
-        @test length(vertices(m)) == 10968
-        @test length(faces(m)) == 3656
+        @test length(vertices(m)) == length(vertices(mf))
+        @test length(faces(m)) == length(faces(mf))
     end
 
     @testset "respect origin" begin