Support positive sign conventions in MarchingCubes

docs/src/examples.md

@@ -3,7 +3,7 @@
 
 # NRRD Data
 
-The file for this example can be found here: http://www.slicer.org/slicerWiki/images/0/00/CTA-cardio.nrrd
+The file for this example can be found here: [http://www.slicer.org/slicerWiki/images/0/00/CTA-cardio.nrrd](http://www.slicer.org/slicerWiki/images/0/00/CTA-cardio.nrrd)
 
 ```
 
@@ -13,25 +13,20 @@ using Meshing
 using MeshIO
 using GeometryTypes
 
-here = dirname(@__FILE__)
-
 # load the file as an AxisArray
-ctacardio = load(here*"/../data/CTA-cardio.nrrd")
-q = -100 # isolevel
-
-
-# flip sign on nrrd data since we need inside to be negative
-# TODO
-for i in eachindex(ctacardio.data)
-    ctacardio.data[i] = -ctacardio.data[i]
-end
+ctacardio = load("CTA-cardio.nrrd")
 
 # convert AxisArray to SignedDistanceField
 ctasdf = SignedDistanceField(HyperRectangle(Vec(0,0,0), Vec(10,10,10)),ctacardio.data)
 
+# use marching cubes with iso at 100
+algo = MarchingCubes(iso=100, insidepositive=true)
+
 # generate the mesh using marching cubes
-mc = HomogenousMesh{Point{3,Float32},Face{3,Int}}(ctasdf, MarchingCubes(q))
+mc = HomogenousMesh{Point{3,Float32},Face{3,Int}}(ctasdf, algo)
 
 # save the file as a PLY file
 save("ctacardio_mc.ply", mc)
 ```
+
+![cta cardio](./img/ctacardio.png)

docs/src/internals.md

@@ -8,6 +8,7 @@ This is some brief documentation on the basic internal functions.
 
 ```@docs
 Meshing._get_cubeindex
+Meshing._get_cubeindex_pos
 Meshing._determine_types
 ```
 

src/algorithmtypes.jl

@@ -1,7 +1,7 @@
 
 """
-    MarchingCubes(iso=0.0, eps=1e-3, reduceverts=true)
-    MarchingCubes(;iso=0.0, eps=1e-3, reduceverts=true)
+    MarchingCubes(iso=0.0, eps=1e-3, reduceverts=true, insidepositive=false)
+    MarchingCubes(;iso=0.0, eps=1e-3, reduceverts=true, insidepositive=false)
     MarchingCubes(iso)
     MarchingCubes(iso,eps)
 
@@ -10,24 +10,26 @@ This algorithm provides a good balance between performance and vertex count.
 In contrast to the other algorithms, vertices may be repeated, so face counts
 may be large.
 
-- `iso` specifies the iso level to use for surface extraction.
-- `eps` is the tolerence around a voxel corner to ensure manifold mesh generation.
-- `reduceverts` if true will merge vertices within a voxel to reduce mesh size by around 30% and with no performance penalty.
+- `iso` (default: 0.0) specifies the iso level to use for surface extraction.
+- `eps` (default: 1e-3) is the tolerence around a voxel corner to ensure manifold mesh generation.
+- `reduceverts` (default: true) if true will merge vertices within a voxel to reduce mesh size by around 30% and with no performance penalty.
+- `insidepositive` (default: false) set true if the sign convention inside the surface is positive, common for NRRD and DICOM data
 """
 struct MarchingCubes{T} <: AbstractMeshingAlgorithm
     iso::T
     eps::T
     reduceverts::Bool
+    insidepositive::Bool
 end
 
-MarchingCubes(;iso::T1=0.0, eps::T2=1e-3, reduceverts::Bool=true) where {T1, T2} = MarchingCubes{promote_type(T1, T2)}(iso, eps, reduceverts)
+MarchingCubes(;iso::T1=0.0, eps::T2=1e-3, reduceverts::Bool=true, insidepositive::Bool=false) where {T1, T2} = MarchingCubes{promote_type(T1, T2)}(iso, eps, reduceverts, insidepositive)
 MarchingCubes(iso) = MarchingCubes(iso=iso)
 MarchingCubes(iso,eps) = MarchingCubes(iso=iso,eps=eps)
 
 
 """
-    MarchingTetrahedra(iso=0.0, eps=1e-3, reduceverts=true)
-    MarchingTetrahedra(;iso=0.0, eps=1e-3, reduceverts=true)
+    MarchingTetrahedra(iso=0.0, eps=1e-3, reduceverts=true, insidepositive=false)
+    MarchingTetrahedra(;iso=0.0, eps=1e-3, reduceverts=true, insidepositive=false)
     MarchingTetrahedra(iso)
     MarchingTetrahedra(iso,eps)
 
@@ -39,20 +41,22 @@ making this algorithm useful for topological analysis.
 - `iso` specifies the iso level to use for surface extraction.
 - `eps` is the tolerence around a voxel corner to ensure manifold mesh generation.
 - `reduceverts` reserved for future use.
+- `insidepositive` reserved for future use.
 """
 struct MarchingTetrahedra{T} <: AbstractMeshingAlgorithm
     iso::T
     eps::T
     reduceverts::Bool
+    insidepositive::Bool
 end
 
-MarchingTetrahedra(;iso::T1=0.0, eps::T2=1e-3, reduceverts::Bool=true) where {T1, T2} = MarchingTetrahedra{promote_type(T1, T2)}(iso, eps, reduceverts)
+MarchingTetrahedra(;iso::T1=0.0, eps::T2=1e-3, reduceverts::Bool=true, insidepositive::Bool=false) where {T1, T2} = MarchingTetrahedra{promote_type(T1, T2)}(iso, eps, reduceverts, insidepositive)
 MarchingTetrahedra(iso) = MarchingTetrahedra(iso=iso)
 MarchingTetrahedra(iso,eps) = MarchingTetrahedra(iso=iso,eps=eps)
 
 """
-    NaiveSurfaceNets(iso=0.0, eps=1e-3, reduceverts=true)
-    NaiveSurfaceNets(;iso=0.0, eps=1e-3, reduceverts=true)
+    NaiveSurfaceNets(iso=0.0, eps=1e-3, reduceverts=true, insidepositive=false)
+    NaiveSurfaceNets(;iso=0.0, eps=1e-3, reduceverts=true, insidepositive=false)
     NaiveSurfaceNets(iso)
     NaiveSurfaceNets(iso,eps)
 
@@ -64,13 +68,15 @@ guarantee accuracy and generates quad faces.
 - `iso` specifies the iso level to use for surface extraction.
 - `eps` is the tolerence around a voxel corner to ensure manifold mesh generation.
 - `reduceverts` reserved for future use.
+- `insidepositive` reserved for future use.
 """
 struct NaiveSurfaceNets{T} <: AbstractMeshingAlgorithm
     iso::T
     eps::T
     reduceverts::Bool
+    insidepositive::Bool
 end
 
-NaiveSurfaceNets(;iso::T1=0.0, eps::T2=1e-3, reduceverts::Bool=true) where {T1, T2} = NaiveSurfaceNets{promote_type(T1, T2)}(iso, eps, reduceverts)
+NaiveSurfaceNets(;iso::T1=0.0, eps::T2=1e-3, reduceverts::Bool=true, insidepositive::Bool=false) where {T1, T2} = NaiveSurfaceNets{promote_type(T1, T2)}(iso, eps, reduceverts, insidepositive)
 NaiveSurfaceNets(iso) = NaiveSurfaceNets(iso=iso)
 NaiveSurfaceNets(iso,eps) = NaiveSurfaceNets(iso=iso,eps=eps)

src/common.jl

@@ -4,7 +4,8 @@
 
 given `iso_vals` and iso, return an 8 bit value corresponding
 to each corner of a cube. In each bit position, 
-0 indicates in the isosurface and 1 indicates outside the surface
+0 indicates in the isosurface and 1 indicates outside the surface,
+where the sign convention indicates negative inside the surface
 """
 @inline function _get_cubeindex(iso_vals, iso)
     cubeindex = iso_vals[1] < iso ? 0x01 : 0x00
@@ -18,6 +19,26 @@ to each corner of a cube. In each bit position,
     cubeindex
 end
 
+"""
+    _get_cubeindex_pos(iso_vals, iso)
+
+given `iso_vals` and iso, return an 8 bit value corresponding
+to each corner of a cube. In each bit position, 
+0 indicates in the isosurface and 1 indicates outside the surface,
+where the sign convention indicates positive inside the surface
+"""
+@inline function _get_cubeindex_pos(iso_vals, iso)
+    cubeindex = iso_vals[1] > iso ? 0x01 : 0x00
+    iso_vals[2] > iso && (cubeindex |= 0x02)
+    iso_vals[3] > iso && (cubeindex |= 0x04)
+    iso_vals[4] > iso && (cubeindex |= 0x08)
+    iso_vals[5] > iso && (cubeindex |= 0x10)
+    iso_vals[6] > iso && (cubeindex |= 0x20)
+    iso_vals[7] > iso && (cubeindex |= 0x40)
+    iso_vals[8] > iso && (cubeindex |= 0x80)
+    cubeindex
+end
+
 """
     _determine_types(meshtype, fieldtype=Float64, facelen=3)
 

src/marching_cubes.jl

@@ -22,7 +22,7 @@ end
 function marching_cubes(sdf::SignedDistanceField{3,ST,FT}, ::Type{VertType}, ::Type{FaceType},
                                iso=0.0,
                                MT::Type{M}=SimpleMesh{Point{3,Float64},Face{3,Int}},
-                               eps=0.00001, reduceverts=true) where {ST,FT,M<:AbstractMesh, VertType, FaceType}
+                               eps=0.00001, reduceverts=true, insidepositive=Val(false)) where {ST,FT,M<:AbstractMesh, VertType, FaceType}
     nx, ny, nz = size(sdf)
     h = HyperRectangle(sdf)
     w = widths(h)
@@ -54,7 +54,7 @@ function marching_cubes(sdf::SignedDistanceField{3,ST,FT}, ::Type{VertType}, ::T
 
         #Determine the index into the edge table which
         #tells us which vertices are inside of the surface
-        cubeindex = _get_cubeindex(iso_vals, iso)
+        cubeindex = insidepositive === Val(true) ? _get_cubeindex_pos(iso_vals, iso) : _get_cubeindex(iso_vals, iso)
 
         # Cube is entirely in/out of the surface
         (cubeindex == 0x00 || cubeindex == 0xff) && continue
@@ -85,7 +85,7 @@ function marching_cubes(f::Function,
                         samples::NTuple{3,Int}=(256,256,256),
                         iso=0.0,
                         MT::Type{M}=SimpleMesh{Point{3,Float64},Face{3,Int}},
-                        eps=0.00001, reduceverts=true) where {M<:AbstractMesh, VertType, FaceType}
+                        eps=0.00001, reduceverts=true, insidepositive=Val(false)) where {M<:AbstractMesh, VertType, FaceType}
 
     FT = eltype(VertType)
 
@@ -134,7 +134,7 @@ function marching_cubes(f::Function,
 
         #Determine the index into the edge table which
         #tells us which vertices are inside of the surface
-        cubeindex = _get_cubeindex(iso_vals, iso)
+        cubeindex = insidepositive === Val(true) ? _get_cubeindex_pos(iso_vals, iso) : _get_cubeindex(iso_vals, iso)
 
         # Cube is entirely in/out of the surface
         (cubeindex == 0x00 || cubeindex == 0xff) && continue
@@ -313,15 +313,15 @@ end
 
 function (::Type{MT})(df::SignedDistanceField{3,ST,FT}, method::MarchingCubes)::MT where {MT <: AbstractMesh, ST, FT}
     VertType, FaceType = _determine_types(MT, FT)
-    marching_cubes(df, VertType, FaceType, method.iso, MT, method.eps, method.reduceverts)
+    marching_cubes(df, VertType, FaceType, method.iso, MT, method.eps, method.reduceverts, Val(method.insidepositive))
 end
 
 function (::Type{MT})(f::Function, h::HyperRectangle, size::NTuple{3,Int}, method::MarchingCubes)::MT where {MT <: AbstractMesh}
     VertType, FaceType = _determine_types(MT)
-    marching_cubes(f, h, VertType, FaceType, size, method.iso, MT, method.eps, method.reduceverts)
+    marching_cubes(f, h, VertType, FaceType, size, method.iso, MT, method.eps, method.reduceverts, Val(method.insidepositive))
 end
 
 function (::Type{MT})(f::Function, h::HyperRectangle, method::MarchingCubes; size::NTuple{3,Int}=(128,128,128))::MT where {MT <: AbstractMesh}
     VertType, FaceType = _determine_types(MT)
-    marching_cubes(f, h, VertType, FaceType, size, method.iso, MT, method.eps, method.reduceverts)
+    marching_cubes(f, h, VertType, FaceType, size, method.iso, MT, method.eps, method.reduceverts, Val(method.insidepositive))
 end

test/runtests.jl

@@ -99,6 +99,7 @@ using LinearAlgebra: dot, norm
     @testset "marching cubes" begin
 
         algo = MarchingCubes()
+        algo_pos = MarchingCubes(insidepositive=true)
 
         sdf = SignedDistanceField(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.))) do v
             sqrt(sum(dot(v,v))) - 1 # sphere
@@ -108,6 +109,10 @@ using LinearAlgebra: dot, norm
             sqrt(sum(dot(v,v))) - 1 # sphere
         end
 
+        mf_pos = SimpleMesh(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)),algo_pos, size=(21,21,21)) do v
+            -sqrt(sum(dot(v,v))) + 1 # sphere positive inside
+        end
+
         mfrv = SimpleMesh(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)),MarchingCubes(reduceverts=false), size=(21,21,21)) do v
             sqrt(sum(dot(v,v))) - 1 # sphere
         end
@@ -126,8 +131,10 @@ using LinearAlgebra: dot, norm
         @test length(vertices(m)) == 7320
         @test length(faces(m)) == 3656
         @test length(faces(mf)) == length(faces(mfrv))
+        @test length(faces(mf)) == length(faces(mf_pos))
         @test m == m2
         @test length(vertices(m)) == length(vertices(mf))
+        @test length(vertices(mf)) == length(vertices(mf_pos))
         @test length(faces(m)) == length(faces(mf))
     end