Merge pull request #14 from JuliaGeometry/algorithms

Fix origin and scaling of meshes, and fix type piracy

Author: SimonDanisch

.travis.yml

@@ -18,4 +18,4 @@ script:
   - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi
   - julia -e 'Pkg.clone(pwd()); Pkg.build("Meshing"); Pkg.test("Meshing"; coverage=true)'
 after_success:
-  - julia -e 'cd(Pkg.dir("Meshing")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(process_folder()); Codecov.submit(process_folder())'
+  - julia -e 'cd(Pkg.dir("Meshing")); Pkg.add("Coverage"); using Coverage; Codecov.submit(process_folder())'

src/Meshing.jl

@@ -4,9 +4,13 @@ module Meshing
 
 using GeometryTypes
 
+abstract type AbstractMeshingAlgorithm end
+
 include("marching_tetrahedra.jl")
 include("marching_cubes.jl")
 
-export marching_cubes
+export marching_cubes,
+       MarchingCubes,
+       MarchingTetrahedra
 
 end # module

src/marching_cubes.jl

@@ -287,12 +287,17 @@ However it may generate non-manifold meshes, while Marching
 Tetrahedra guarentees a manifold mesh.
 """
 function marching_cubes(sdf::SignedDistanceField{3,ST,FT},
-        iso=0.0,
-        MT::Type{M}=SimpleMesh{Point{3,Float64},Face{3,Int}}) where {ST,FT,M<:AbstractMesh}
+                               iso=0.0,
+                               MT::Type{M}=SimpleMesh{Point{3,Float64},Face{3,Int}},
+                               eps=0.00001) where {ST,FT,M<:AbstractMesh}
     nx, ny, nz = size(sdf)
     h = HyperRectangle(sdf)
     w = widths(h)
-    s = Point{3,Float64}(w[1]/nx, w[2]/ny, w[3]/nz)
+    orig = origin(HyperRectangle(sdf))
+
+    # we subtract one from the length along each axis because
+    # an NxNxN SDF has N-1 cells on each axis
+    s = Point{3,Float64}(w[1]/(nx-1), w[2]/(ny-1), w[3]/(nz-1))
 
     # arrays for vertices and faces
     vts = Point{3,Float64}[]
@@ -315,63 +320,63 @@ function marching_cubes(sdf::SignedDistanceField{3,ST,FT},
         # Cube is entirely in/out of the surface
         edge_table[cubeindex] == 0 && continue
 
-        points = (Point{3,Float64}(xi-1,yi-1,zi-1).*s,
-                  Point{3,Float64}(xi,yi-1,zi-1).*s,
-                  Point{3,Float64}(xi,yi,zi-1).*s,
-                  Point{3,Float64}(xi-1,yi,zi-1).*s,
-                  Point{3,Float64}(xi-1,yi-1,zi).*s,
-                  Point{3,Float64}(xi,yi-1,zi).*s,
-                  Point{3,Float64}(xi,yi,zi).*s,
-                  Point{3,Float64}(xi-1,yi,zi).*s)
+        points = (Point{3,Float64}(xi-1,yi-1,zi-1) .* s .+ orig,
+                  Point{3,Float64}(xi,yi-1,zi-1) .* s .+ orig,
+                  Point{3,Float64}(xi,yi,zi-1) .* s .+ orig,
+                  Point{3,Float64}(xi-1,yi,zi-1) .* s .+ orig,
+                  Point{3,Float64}(xi-1,yi-1,zi) .* s .+ orig,
+                  Point{3,Float64}(xi,yi-1,zi) .* s .+ orig,
+                  Point{3,Float64}(xi,yi,zi) .* s .+ orig,
+                  Point{3,Float64}(xi-1,yi,zi) .* s .+ orig)
 
         # Find the vertices where the surface intersects the cube
         if (edge_table[cubeindex] & 1 != 0)
           vertlist[1] =
-             vertex_interp(iso,points[1],points[2],sdf[xi,yi,zi],sdf[xi+1,yi,zi])
+             vertex_interp(iso,points[1],points[2],sdf[xi,yi,zi],sdf[xi+1,yi,zi], eps)
         end
         if (edge_table[cubeindex] & 2 != 0)
           vertlist[2] =
-             vertex_interp(iso,points[2],points[3],sdf[xi+1,yi,zi],sdf[xi+1,yi+1,zi])
+             vertex_interp(iso,points[2],points[3],sdf[xi+1,yi,zi],sdf[xi+1,yi+1,zi], eps)
         end
         if (edge_table[cubeindex] & 4 != 0)
           vertlist[3] =
-             vertex_interp(iso,points[3],points[4],sdf[xi+1,yi+1,zi],sdf[xi,yi+1,zi])
+             vertex_interp(iso,points[3],points[4],sdf[xi+1,yi+1,zi],sdf[xi,yi+1,zi], eps)
         end
         if (edge_table[cubeindex] & 8 != 0)
           vertlist[4] =
-             vertex_interp(iso,points[4],points[1],sdf[xi,yi+1,zi],sdf[xi,yi,zi])
+             vertex_interp(iso,points[4],points[1],sdf[xi,yi+1,zi],sdf[xi,yi,zi], eps)
         end
         if (edge_table[cubeindex] & 16 != 0)
           vertlist[5] =
-             vertex_interp(iso,points[5],points[6],sdf[xi,yi,zi+1],sdf[xi+1,yi,zi+1])
+             vertex_interp(iso,points[5],points[6],sdf[xi,yi,zi+1],sdf[xi+1,yi,zi+1], eps)
         end
         if (edge_table[cubeindex] & 32 != 0)
           vertlist[6] =
-             vertex_interp(iso,points[6],points[7],sdf[xi+1,yi,zi+1],sdf[xi+1,yi+1,zi+1])
+             vertex_interp(iso,points[6],points[7],sdf[xi+1,yi,zi+1],sdf[xi+1,yi+1,zi+1], eps)
         end
         if (edge_table[cubeindex] & 64 != 0)
           vertlist[7] =
-             vertex_interp(iso,points[7],points[8],sdf[xi+1,yi+1,zi+1],sdf[xi,yi+1,zi+1])
+             vertex_interp(iso,points[7],points[8],sdf[xi+1,yi+1,zi+1],sdf[xi,yi+1,zi+1], eps)
         end
         if (edge_table[cubeindex] & 128 != 0)
           vertlist[8] =
-             vertex_interp(iso,points[8],points[5],sdf[xi,yi+1,zi+1],sdf[xi,yi,zi+1])
+             vertex_interp(iso,points[8],points[5],sdf[xi,yi+1,zi+1],sdf[xi,yi,zi+1], eps)
         end
         if (edge_table[cubeindex] & 256 != 0)
           vertlist[9] =
-             vertex_interp(iso,points[1],points[5],sdf[xi,yi,zi],sdf[xi,yi,zi+1])
+             vertex_interp(iso,points[1],points[5],sdf[xi,yi,zi],sdf[xi,yi,zi+1], eps)
         end
         if (edge_table[cubeindex] & 512 != 0)
           vertlist[10] =
-             vertex_interp(iso,points[2],points[6],sdf[xi+1,yi,zi],sdf[xi+1,yi,zi+1])
+             vertex_interp(iso,points[2],points[6],sdf[xi+1,yi,zi],sdf[xi+1,yi,zi+1], eps)
         end
         if (edge_table[cubeindex] & 1024 != 0)
           vertlist[11] =
-             vertex_interp(iso,points[3],points[7],sdf[xi+1,yi+1,zi],sdf[xi+1,yi+1,zi+1])
+             vertex_interp(iso,points[3],points[7],sdf[xi+1,yi+1,zi],sdf[xi+1,yi+1,zi+1], eps)
         end
         if (edge_table[cubeindex] & 2048 != 0)
           vertlist[12] =
-             vertex_interp(iso,points[4],points[8],sdf[xi,yi+1,zi],sdf[xi,yi+1,zi+1])
+             vertex_interp(iso,points[4],points[8],sdf[xi,yi+1,zi],sdf[xi,yi+1,zi+1], eps)
         end
 
         # Create the triangle
@@ -399,3 +404,14 @@ function vertex_interp(iso, p1, p2, valp1, valp2, eps = 0.00001)
 
     return p
 end
+
+struct MarchingCubes{T} <: AbstractMeshingAlgorithm
+     iso::T
+     eps::T
+end
+
+MarchingCubes(iso::T1=0.0, eps::T2=1e-3) where {T1, T2} = MarchingCubes{promote_type(T1, T2)}(iso, eps)
+
+function (::Type{MT})(df::SignedDistanceField, method::MarchingCubes)::MT where {MT <: AbstractMesh}
+     marching_cubes(df, method.iso, MT, method.eps)
+end

src/marching_tetrahedra.jl

@@ -292,16 +292,51 @@ end
 
 isosurface(lsf,isoval) = isosurface(lsf,isoval, convert(eltype(lsf), 0.001))
 
+"""
+The marchingTetrahedra function returns vertices on the (1-based) indices of the
+SDF's data, ignoring its actual bounds. This function adjusts the vertices in
+place so that they correspond to points within the SDF bounds.
+"""
+function _correct_vertices!(vts, sdf::SignedDistanceField)
+    bounds = HyperRectangle(sdf)
+    orig = origin(bounds)
+    w = widths(bounds)
+    s = w ./ Point(size(sdf) .- 1)  # subtract 1 because an SDF with N points per side has N-1 cells
+    for i in eachindex(vts)
+        vts[i] = (vts[i] .- 1) .* s .+ orig  # subtract 1 to fix 1-indexing
+    end
+end
+
 
 function (::Type{MT})(volume::Array{T, 3}, iso_val::Real, eps_val=0.001) where {MT <: AbstractMesh, T}
+    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)
     iso_val = convert(T, iso_val)
     eps_val = convert(T, eps_val)
     vts, fcs = isosurface(volume, iso_val, eps_val)
     MT(vts, fcs)
 end
 
 function (::Type{MT})(df::SignedDistanceField, eps_val=0.001) where MT <: AbstractMesh
+    Base.depwarn("(mesh type)(sdf::SignedDistanceField, eps_val) is deprecated. Please use: (mesh type)(sdf, MarchingTetrahedra(0, eps))", :Meshing_MT_sdf_eps)
     vts, fcs = isosurface(df.data, 0.0, eps_val)
+    _correct_vertices!(vts, df)
     MT(vts, fcs)
 end
 
+struct MarchingTetrahedra{T} <: AbstractMeshingAlgorithm
+    iso::T
+    eps::T
+end
+
+MarchingTetrahedra(iso::T1=0.0, eps::T2=1e-3) where {T1, T2} = MarchingTetrahedra{promote_type(T1, T2)}(iso, eps)
+
+function (::Type{MT})(sdf::SignedDistanceField, method::MarchingTetrahedra)::MT where {MT <: AbstractMesh}
+    vts, fcs = isosurface(sdf.data, method.iso, method.eps)
+    _correct_vertices!(vts, sdf)
+    MT(vts, fcs)
+end
+
+function (::Type{MT}){MT <: AbstractMesh, T}(volume::Array{T, 3}, method::MarchingTetrahedra)::MT
+    vts, fcs = isosurface(volume, convert(T, method.iso), convert(T, method.eps))
+    MT(vts, fcs)
+end

test/runtests.jl

@@ -2,46 +2,100 @@ using Meshing
 using Base.Test
 using GeometryTypes
 
-# Produce a level set function that is a noisy version of the distance from
-# the origin (such that level sets are noisy spheres).
-#
-# The noise should exercise marching tetrahedra's ability to produce a water-
-# tight surface in all cases (unlike standard marching cubes).
-#
-N = 10
-sigma = 1.0
-distance = Float32[ sqrt(Float32(i*i+j*j+k*k)) for i = -N:N, j = -N:N, k = -N:N ]
-distance = distance + sigma*rand(2*N+1,2*N+1,2*N+1)
-
-# Extract an isosurface.
-#
-lambda = N-2*sigma # isovalue
-
-msh = HomogenousMesh(distance,lambda)
-
-s2 = SignedDistanceField(HyperRectangle(Vec(0,0,0.),Vec(1,1,1.))) do v
-    sqrt(sum(dot(v,v))) - 1 # sphere
-end
 
-msh = HomogenousMesh(s2)
-@test length(vertices(msh)) == 973
-@test length(faces(msh)) == 1830
+@testset "meshing" begin
+    @testset "noisy spheres" begin
+        # Produce a level set function that is a noisy version of the distance from
+        # the origin (such that level sets are noisy spheres).
+        #
+        # The noise should exercise marching tetrahedra's ability to produce a water-
+        # tight surface in all cases (unlike standard marching cubes).
+        #
+        N = 10
+        sigma = 1.0
+        distance = Float32[ sqrt(Float32(i*i+j*j+k*k)) for i = -N:N, j = -N:N, k = -N:N ]
+        distance = distance + sigma*rand(2*N+1,2*N+1,2*N+1)
+
+        # Extract an isosurface.
+        #
+        lambda = N-2*sigma # isovalue
+
+        msh = HomogenousMesh(distance, MarchingTetrahedra(lambda))
+
+        s2 = SignedDistanceField(HyperRectangle(Vec(0,0,0.),Vec(1,1,1.))) do v
+            sqrt(sum(dot(v,v))) - 1 # sphere
+        end
+
+        msh = HomogenousMesh(s2, MarchingTetrahedra())
+        @test length(vertices(msh)) == 973
+        @test length(faces(msh)) == 1830
+    end
+
+    @testset "vertex interpolation" begin
+        @test Meshing.vertex_interp(0, Vec(0,0,0), Vec(0,1,0), -1, 1) == Vec(0,0.5,0)
+        @test Meshing.vertex_interp(-1, Vec(0,0,0), Vec(0,1,0), -1, 1) == Vec(0,0,0)
+        @test Meshing.vertex_interp(1, Vec(0,0,0), Vec(0,1,0), -1, 1) == Vec(0,1,0)
+    end
+
+    @testset "marching cubes" begin
+        sdf = SignedDistanceField(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.))) 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 "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)
+        bounds = HyperRectangle(Vec(-1, -1, -1), Vec(2, 2, 2))
+        resolution = 0.1
+        sdf = SignedDistanceField(f, bounds, resolution)
+
+        for algorithm in (MarchingCubes(0.5), MarchingTetrahedra(0.5))
+            mesh = @inferred GLNormalMesh(sdf, algorithm)
+            # should be centered on the origin
+            @test mean(vertices(mesh)) ≈ [0, 0, 0] atol=0.15*resolution
+            # and should be symmetric about the origin
+            @test maximum(vertices(mesh)) ≈ [0.5, 0.5, 0.5]
+            @test minimum(vertices(mesh)) ≈ [-0.5, -0.5, -0.5]
+        end
+    end
+
+    @testset "AbstractMeshingAlgorithm interface" begin
+        f = x -> norm(x) - 0.5
+        bounds = HyperRectangle(Vec(-1, -1, -1), Vec(2, 2, 2))
+        resolution = 0.1
+        sdf = SignedDistanceField(f, bounds, resolution)
+
+        @testset "marching cubes" begin
+            m1 = @test_nowarn marching_cubes(sdf, 0.0, GLNormalMesh)
+            m2 = @test_nowarn GLNormalMesh(sdf, MarchingCubes())
+            @test vertices(m1) == vertices(m2)
+            @test faces(m1) == faces(m2)
+        end
 
-# Vertex interpolation
-@test Meshing.vertex_interp(0, Vec(0,0,0), Vec(0,1,0), -1, 1) == Vec(0,0.5,0)
-@test Meshing.vertex_interp(-1, Vec(0,0,0), Vec(0,1,0), -1, 1) == Vec(0,0,0)
-@test Meshing.vertex_interp(1, Vec(0,0,0), Vec(0,1,0), -1, 1) == Vec(0,1,0)
+        @testset "marching tetrahedra" begin
+            m1 = @test_warn "is deprecated" GLNormalMesh(sdf)
+            m2 = @test_warn "is deprecated" GLNormalMesh(sdf, 1e-3)
+            m3 = @test_nowarn GLNormalMesh(sdf, MarchingTetrahedra())
+            @test vertices(m1) == vertices(m2) == vertices(m3)
+            @test faces(m1) == faces(m2) == faces(m3)
 
-# marching cubes
-sdf = SignedDistanceField(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.))) do v
-    sqrt(sum(dot(v,v))) - 1 # sphere
+            m4 = @test_warn "is deprecated" GLNormalMesh(sdf.data, 0.5)
+            m5 = @test_nowarn GLNormalMesh(sdf.data, MarchingTetrahedra(0.5))
+            @test vertices(m4) == vertices(m5)
+            @test faces(m4) == faces(m5)
+        end
+    end
 end
 
-m = marching_cubes(sdf,0)
-m2 = marching_cubes(sdf)
-@test length(vertices(m)) == 10968
-@test length(faces(m)) == 3656
-@test m == m2
 
 if "--profile" in ARGS
     HomogenousMesh(s2)