From 8bb280ab8540961d72b137ce4aab12122659c181 Mon Sep 17 00:00:00 2001
From: Steve Kelly <kd2cca@gmail.com>
Date: Sat, 29 Jun 2024 23:35:36 -0400
Subject: [PATCH 1/8] new MarchingCubes interface

---
 Project.toml          |   1 +
 docs/src/internals.md |   2 -
 src/lut/mc.jl         |  48 ++++++++---------
 src/marching_cubes.jl | 121 +++++++++++++++++-------------------------
 test/runtests.jl      |  12 ++---
 5 files changed, 81 insertions(+), 103 deletions(-)

diff --git a/Project.toml b/Project.toml
index b44c442..b03a596 100644
--- a/Project.toml
+++ b/Project.toml
@@ -3,6 +3,7 @@ uuid = "e6723b4c-ebff-59f1-b4b7-d97aa5274f73"
 version = "0.7.0"
 
 [deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
diff --git a/docs/src/internals.md b/docs/src/internals.md
index 572bbf2..4bfedca 100644
--- a/docs/src/internals.md
+++ b/docs/src/internals.md
@@ -8,8 +8,6 @@ This is some brief documentation on the basic internal functions.
 ```@docs
 Meshing._DEFAULT_SAMPLES
 Meshing._get_cubeindex
-Meshing._get_cubeindex_pos
-Meshing._determine_types
 ```
 
 ## Marching Cubes
diff --git a/src/lut/mc.jl b/src/lut/mc.jl
index 570bfdd..42f0a59 100644
--- a/src/lut/mc.jl
+++ b/src/lut/mc.jl
@@ -578,30 +578,30 @@ const _mc_eq_mapping = (0x01, 0x01, 0x02, 0x01, 0x03, 0x02, 0x04, 0x01, 0x02, 0x
                         0x15, 0x10, 0x10, 0x03, 0x0c, 0x14, 0x11, 0x01, 0x18, 0x0c, 0x0c,
                         0x11, 0x0c, 0x14, 0x11, 0x01, 0x0c, 0x11, 0x14, 0x01, 0x11, 0x01, 0x01)
 
-const _mc_connectivity = ((0x01, 0x02, 0x03, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x04, 0x02, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x01, 0x04, 0x02, 0x04, 0x05, 0x02, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x03, 0x02, 0x04, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x01, 0x04, 0x02, 0x05, 0x06, 0x07, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x01, 0x03, 0x04, 0x01, 0x04, 0x05, 0x04, 0x03, 0x06, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08, 0x09, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x04, 0x01, 0x03, 0x04, 0x03, 0x05, 0x04, 0x05, 0x06, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x01, 0x04, 0x02, 0x01, 0x05, 0x04, 0x06, 0x02, 0x04, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x04, 0x06, 0x07, 0x06, 0x05, 0x08, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x01, 0x03, 0x04, 0x04, 0x03, 0x05, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x01, 0x03, 0x04, 0x01, 0x04, 0x05, 0x06, 0x04, 0x03, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x05, 0x07, 0x06, 0x05, 0x08, 0x07, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x03, 0x02, 0x04, 0x05, 0x06, 0x07, 0x05, 0x07, 0x08, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x01, 0x03, 0x04, 0x02, 0x05, 0x03, 0x06, 0x03, 0x07, 0x05, 0x07, 0x03),
-                            (0x01, 0x02, 0x03, 0x04, 0x01, 0x03, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08, 0x09, 0x0a, 0x0b, 0x0c, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x01, 0x04, 0x02, 0x01, 0x05, 0x04, 0x06, 0x04, 0x05, 0x07, 0x08, 0x09),
-                            (0x01, 0x02, 0x03, 0x01, 0x03, 0x04, 0x05, 0x06, 0x03, 0x06, 0x04, 0x03, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x01, 0x04, 0x02, 0x01, 0x05, 0x04, 0x04, 0x06, 0x02, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x04, 0x02, 0x01, 0x04, 0x05, 0x02, 0x04, 0x06, 0x05, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x01, 0x04, 0x02, 0x05, 0x06, 0x07, 0x06, 0x08, 0x07, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x02, 0x04, 0x03, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00))
+const _mc_connectivity = ((0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
+                            (0x04, 0x02, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
+                            (0x04, 0x05, 0x06, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
+                            (0x01, 0x04, 0x02, 0x04, 0x05, 0x02, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
+                            (0x03, 0x02, 0x04, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
+                            (0x01, 0x04, 0x02, 0x05, 0x06, 0x07, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
+                            (0x01, 0x03, 0x04, 0x01, 0x04, 0x05, 0x04, 0x03, 0x06, 0x00, 0x00, 0x00),
+                            (0x04, 0x05, 0x06, 0x07, 0x08, 0x09, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
+                            (0x04, 0x01, 0x03, 0x04, 0x03, 0x05, 0x04, 0x05, 0x06, 0x00, 0x00, 0x00),
+                            (0x01, 0x04, 0x02, 0x01, 0x05, 0x04, 0x06, 0x02, 0x04, 0x00, 0x00, 0x00),
+                            (0x04, 0x05, 0x06, 0x04, 0x06, 0x07, 0x06, 0x05, 0x08, 0x00, 0x00, 0x00),
+                            (0x01, 0x03, 0x04, 0x04, 0x03, 0x05, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
+                            (0x01, 0x03, 0x04, 0x01, 0x04, 0x05, 0x06, 0x04, 0x03, 0x00, 0x00, 0x00),
+                            (0x04, 0x05, 0x06, 0x05, 0x07, 0x06, 0x05, 0x08, 0x07, 0x00, 0x00, 0x00),
+                            (0x03, 0x02, 0x04, 0x05, 0x06, 0x07, 0x05, 0x07, 0x08, 0x00, 0x00, 0x00),
+                            (0x01, 0x03, 0x04, 0x02, 0x05, 0x03, 0x06, 0x03, 0x07, 0x05, 0x07, 0x03),
+                            (0x04, 0x01, 0x03, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
+                            (0x04, 0x05, 0x06, 0x07, 0x08, 0x09, 0x0a, 0x0b, 0x0c, 0x00, 0x00, 0x00),
+                            (0x01, 0x04, 0x02, 0x01, 0x05, 0x04, 0x06, 0x04, 0x05, 0x07, 0x08, 0x09),
+                            (0x01, 0x03, 0x04, 0x05, 0x06, 0x03, 0x06, 0x04, 0x03, 0x00, 0x00, 0x00),
+                            (0x01, 0x04, 0x02, 0x01, 0x05, 0x04, 0x04, 0x06, 0x02, 0x00, 0x00, 0x00),
+                            (0x04, 0x02, 0x01, 0x04, 0x05, 0x02, 0x04, 0x06, 0x05, 0x00, 0x00, 0x00),
+                            (0x01, 0x04, 0x02, 0x05, 0x06, 0x07, 0x06, 0x08, 0x07, 0x00, 0x00, 0x00),
+                            (0x02, 0x04, 0x03, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00))
 
 const _mc_edge_list = ((0x01, 0x02), (0x02, 0x03), (0x03, 0x04), (0x04, 0x01),
                         (0x05, 0x06), (0x06, 0x07), (0x07, 0x08), (0x08, 0x05),
diff --git a/src/marching_cubes.jl b/src/marching_cubes.jl
index be4645e..7a934ee 100644
--- a/src/marching_cubes.jl
+++ b/src/marching_cubes.jl
@@ -2,17 +2,15 @@
 #Look up Table
 include("lut/mc.jl")
 
-function isosurface(sdf::AbstractArray{T, 3}, method::MarchingCubes, ::Type{VertType}=SVector{3,Float64}, ::Type{FaceType}=SVector{3, Int};
-                    origin=VertType(-1,-1,-1), widths=VertType(2,2,2)) where {T, VertType, FaceType}
+function isosurface(sdf::AbstractArray{T, 3}, method::MarchingCubes, X=-1:1, Y=-1:1, Z=-1:1) where {T}
     nx, ny, nz = size(sdf)
 
-    # we subtract one from the length along each axis because
-    # an NxNxN SDF has N-1 cells on each axis
-    s = VertType(widths[1]/(nx-1), widths[2]/(ny-1), widths[3]/(nz-1))
+    vts = NTuple{3, float(T)}[]
+    fcs = NTuple{3, Int}[]
 
-    # arrays for vertices and faces
-    vts = VertType[]
-    fcs = FaceType[]
+    xp = LinRange(first(X), last(X), nx)
+    yp = LinRange(first(Y), last(Y), ny)
+    zp = LinRange(first(Z), last(Z), nz)
 
     @inbounds for xi = 1:nx-1, yi = 1:ny-1, zi = 1:nz-1
 
@@ -32,33 +30,32 @@ function isosurface(sdf::AbstractArray{T, 3}, method::MarchingCubes, ::Type{Vert
         # Cube is entirely in/out of the surface
         _no_triangles(cubeindex) && continue
 
-        points = mc_vert_points(xi,yi,zi,s,origin,VertType)
+        points = mc_vert_points(xi+1,yi+1,zi+1,xp,yp,zp)
 
-        # Create the triangle
-        _mc_unique_triangles!(method, points, iso_vals, vts, fcs, cubeindex, FaceType)
+        # process the voxel
+        process_mc_voxel!(vts, fcs, cubeindex, points, method.iso, iso_vals)
     end
     vts,fcs
 end
 
 
-function isosurface(f::Function, method::MarchingCubes, ::Type{VertType}=SVector{3,Float64}, ::Type{FaceType}=SVector{3, Int};
-                    origin=VertType(-1,-1,-1), widths=VertType(2,2,2), samples::NTuple{3,T}=_DEFAULT_SAMPLES) where {T <: Integer, VertType, FaceType}
+function isosurface(f::Function, method::MarchingCubes, X=-1:1, Y=-1:1, Z=-1:1; samples::NTuple{3,T}=_DEFAULT_SAMPLES) where {T <: Integer}
 
     nx, ny, nz = samples[1], samples[2], samples[3]
 
-    # we subtract one from the length along each axis because
-    # an NxNxN SDF has N-1 cells on each axis
-    s = VertType(widths[1]/(nx-1), widths[2]/(ny-1), widths[3]/(nz-1))
+    # find widest type
+    FT = promote_type(eltype(first(X)), eltype(first(Y)), eltype(first(Z)), eltype(T))
 
-    # arrays for vertices and faces
-    vts = VertType[]
-    fcs = FaceType[]
+    vts = NTuple{3, float(FT)}[]
+    fcs = NTuple{3, Int}[]
+
+    xp = LinRange(first(X), last(X), nx)
+    yp = LinRange(first(Y), last(Y), ny)
+    zp = LinRange(first(Z), last(Z), nz)
 
-    zv = zero(eltype(VertType))
-    iso_vals = (zv,zv,zv,zv,zv,zv,zv,zv)
     @inbounds for xi = 1:nx-1, yi = 1:ny-1, zi = 1:nz-1
 
-        points = mc_vert_points(xi,yi,zi,s,origin,VertType)
+        points = mc_vert_points(xi+1,yi+1,zi+1,xp,yp,zp)
 
         iso_vals = (f(points[1]),
                     f(points[2]),
@@ -68,7 +65,6 @@ function isosurface(f::Function, method::MarchingCubes, ::Type{VertType}=SVector
                     f(points[6]),
                     f(points[7]),
                     f(points[8]))
-        # TODO: Cache points
 
         #Determine the index into the edge table which
         #tells us which vertices are inside of the surface
@@ -77,83 +73,66 @@ function isosurface(f::Function, method::MarchingCubes, ::Type{VertType}=SVector
         # Cube is entirely in/out of the surface
         _no_triangles(cubeindex) && continue
 
-        # Create the triangle
-        _mc_unique_triangles!(method, points, iso_vals, vts, fcs, cubeindex, FaceType)
+        process_mc_voxel!(vts, fcs, cubeindex, points, method.iso, iso_vals)
     end
     vts,fcs
 end
 
-"""
-    _mc_unique_triangles!(vts, fcs, vertlist, cubeindex, FaceType)
 
-Create triangles by only adding unique vertices within the voxel.
-Each face may share a reference to a vertex with another face.
-"""
-function _mc_unique_triangles!(method, points, iso_vals, vts, fcs, cubeindex, ::Type{FaceType}) where {FaceType}
+function process_mc_voxel!(vts, fcs, cubeindex, points, iso, iso_vals)
+
+    fct = length(vts)
+
     @inbounds begin
-        fct = length(vts)
+        vert_to_add = _mc_verts[cubeindex]
+        for i = 1:12
+            vt = vert_to_add[i]
+            iszero(vt) && break
+            ed = _mc_edge_list[vt]
+            push!(vts, vertex_interp(iso,points[ed[1]],points[ed[2]],iso_vals[ed[1]],iso_vals[ed[2]]))
+        end
+    end
 
-        find_vertices_interp!(vts, points, iso_vals, cubeindex, method.iso, method.eps)
+    @inbounds begin
 
         offsets = _mc_connectivity[_mc_eq_mapping[cubeindex]]
 
         # There is atleast one face so we can push it immediately
-        push!(fcs, FaceType(fct+offsets[3], fct+offsets[2], fct+offsets[1]))
+        push!(fcs, (fct+3, fct+3, fct+1))
 
-        for i in (4,7,10,13)
+        for i in (1,4,7,10)
             iszero(offsets[i]) && return
-            push!(fcs, FaceType(fct+offsets[i+2], fct+offsets[i+1], fct+offsets[i]))
+            push!(fcs, (fct+offsets[i+2], fct+offsets[i+1], fct+offsets[i]))
         end
     end
 end
 
 
-"""
-    find_vertices_interp(points, iso_vals, cubeindex, iso, eps)
-
-Find the vertices where the surface intersects the cube
-"""
-function find_vertices_interp!(vts, points, iso_vals, cubeindex, iso, eps)
-    @inbounds begin
-        vert_to_add = _mc_verts[cubeindex]
-        for i = 1:12
-            vt = vert_to_add[i]
-            iszero(vt) && break
-            ed = _mc_edge_list[vt]
-            push!(vts, vertex_interp(iso,points[ed[1]],points[ed[2]],iso_vals[ed[1]],iso_vals[ed[2]], eps))
-        end
-    end
-end
-
-"""
-    vertex_interp(iso, p1, p2, valp1, valp2, eps = 0.00001)
+"""    vertex_interp(iso, p1, p2, valp1, valp2)
 
 Linearly interpolate the position where an isosurface cuts
 an edge between two vertices, each with their own scalar value
 """
-function vertex_interp(iso, p1, p2, valp1, valp2, eps = 0.00001)
+function vertex_interp(iso, p1, p2, valp1, valp2)
 
-    abs(iso - valp1) < eps && return p1
-    abs(iso - valp2) < eps && return p2
-    abs(valp1-valp2) < eps && return p1
-    mu = (iso - valp1) / (valp2 - valp1)
-    p = p1 + mu * (p2 - p1)
+    mu = (iso .- valp1) / (valp2 .- valp1)
+    p = p1 .+ mu .* (p2 .- p1)
 
     return p
 end
 
 """
-    mc_vert_points(xi,yi,zi, scale, origin, ::Type{VertType})
+    mc_vert_points(xi,yi,zi,xp,yp,zp)
 
 Returns a tuple of 8 points corresponding to each corner of a cube
 """
-@inline function mc_vert_points(xi,yi,zi, scale, origin, ::Type{VertType}) where VertType
-    (VertType(xi-1,yi-1,zi-1) .* scale .+ origin,
-     VertType(xi  ,yi-1,zi-1) .* scale .+ origin,
-     VertType(xi  ,yi  ,zi-1) .* scale .+ origin,
-     VertType(xi-1,yi  ,zi-1) .* scale .+ origin,
-     VertType(xi-1,yi-1,zi  ) .* scale .+ origin,
-     VertType(xi  ,yi-1,zi  ) .* scale .+ origin,
-     VertType(xi  ,yi  ,zi  ) .* scale .+ origin,
-     VertType(xi-1,yi  ,zi  ) .* scale .+ origin)
+@inline function mc_vert_points(xi,yi,zi,xp,yp,zp)
+    ((xp[xi-1],yp[yi-1],zp[zi-1]),
+    (xp[xi  ],yp[yi-1],zp[zi-1]),   
+    (xp[xi  ],yp[yi  ],zp[zi-1]),
+    (xp[xi-1],yp[yi  ],zp[zi-1]),
+    (xp[xi-1],yp[yi-1],zp[zi  ]),
+    (xp[xi  ],yp[yi-1],zp[zi  ]),
+    (xp[xi  ],yp[yi  ],zp[zi  ]),
+    (xp[xi-1],yp[yi  ],zp[zi  ]))
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index ebe0d91..96b9453 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -46,15 +46,15 @@ end
     resolution = 0.1
 
     for algorithm in (MarchingCubes(0.5),
-                      MarchingTetrahedra(0.5),
-                      MarchingCubes(iso=0.5, reduceverts=false),
-                      MarchingTetrahedra(iso=0.5, reduceverts=false))
+                      #MarchingTetrahedra(0.5),
+                      MarchingCubes(iso=0.5, reduceverts=false))
+                      #MarchingTetrahedra(iso=0.5, reduceverts=false))
         @testset "$(typeof(algorithm))" begin
             # Extract isosurface using a function
-            points, faces = isosurface(f, algorithm, origin=origin, widths=widths, samples=(50, 50, 50))
+            points, faces = isosurface(f, algorithm, samples=(50, 50, 50))
 
             # should be centered on the origin
-            @test mean(points) ≈ [0, 0, 0] atol=0.15*resolution
+            @test mean(collect.(points)) ≈ [0, 0, 0] atol=0.15*resolution
             # and should be symmetric about the origin
             #@test maximum(points) ≈ [0.5, 0.5, 0.5]
             #@test minimum(points) ≈ [-0.5, -0.5, -0.5]
@@ -125,7 +125,7 @@ end
     algo = MarchingCubes()
 
     # Extract isosurfaces using MarchingCubes
-    points_mf, faces_mf = isosurface(sphere_function, algo, origin=SVector(-1, -1, -1), widths=SVector(2, 2, 2), samples=(21, 21, 21))
+    points_mf, faces_mf = isosurface(sphere_function, algo, samples=(21, 21, 21))
 
     # Test the number of vertices and faces
     @test length(points_mf) == 7320

From f6ce852bd503376be3c512c2a8324b49179b473c Mon Sep 17 00:00:00 2001
From: Steve Kelly <kd2cca@gmail.com>
Date: Sun, 30 Jun 2024 09:14:18 -0400
Subject: [PATCH 2/8] add back some optimizations

---
 Project.toml          |  3 --
 src/Meshing.jl        |  2 --
 src/adaptive.jl       |  2 +-
 src/lut/mt.jl         | 20 ++++++-------
 src/marching_cubes.jl | 66 ++++++++++++++++++++++++++++---------------
 test/runtests.jl      | 30 ++++++++++----------
 6 files changed, 68 insertions(+), 55 deletions(-)

diff --git a/Project.toml b/Project.toml
index b03a596..6faa49a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -3,11 +3,8 @@ uuid = "e6723b4c-ebff-59f1-b4b7-d97aa5274f73"
 version = "0.7.0"
 
 [deps]
-BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
-StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
-StaticArrays = "1"
 julia = "1.9"
 
 [extras]
diff --git a/src/Meshing.jl b/src/Meshing.jl
index eb3ff95..cce2316 100644
--- a/src/Meshing.jl
+++ b/src/Meshing.jl
@@ -1,7 +1,5 @@
 module Meshing
 
-using StaticArrays
-
 """
     _DEFAULT_SAMPLES = (24,24,24)
 
diff --git a/src/adaptive.jl b/src/adaptive.jl
index 29b3afe..f2e6619 100644
--- a/src/adaptive.jl
+++ b/src/adaptive.jl
@@ -12,7 +12,7 @@ end
 
 
 function vertices(h::HyperRectangle, ::Type{SV}) where SV
-    z = zero(eltype(SV))
+    z = zero(eltype(h.origin))
     @inbounds begin
         o = SV(h.origin[1],h.origin[2],h.origin[3])
         w = SV(h.widths[1],h.widths[2],h.widths[3])
diff --git a/src/lut/mt.jl b/src/lut/mt.jl
index ae0e280..d02d982 100644
--- a/src/lut/mt.jl
+++ b/src/lut/mt.jl
@@ -22,19 +22,17 @@ Voxel corner and edge indexing conventions
 """
 # this gets vectorized so we want to ensure it is the
 # same type as out vertex
-@inline function voxCrnrPos(::Type{PT}) where {PT}
-    (PT(0, 0, 0),
-    PT(0, 1, 0),
-    PT(1, 1, 0),
-    PT(1, 0, 0),
-    PT(0, 0, 1),
-    PT(0, 1, 1),
-    PT(1, 1, 1),
-    PT(1, 0, 1))
+function voxCrnrPos()
+    ((0, 0, 0),
+    (0, 1, 0),
+    (1, 1, 0),
+    (1, 0, 0),
+    (0, 0, 1),
+    (0, 1, 1),
+    (1, 1, 1),
+    (1, 0, 1))
 end
 
-const voxCrnrPosInt = voxCrnrPos(SVector{3,UInt8})
-
 # the voxel IDs at either end of the tetrahedra edges, by edge ID
 const voxEdgeCrnrs = ((0x01, 0x02),
                       (0x02, 0x03),
diff --git a/src/marching_cubes.jl b/src/marching_cubes.jl
index 7a934ee..1035328 100644
--- a/src/marching_cubes.jl
+++ b/src/marching_cubes.jl
@@ -30,7 +30,7 @@ function isosurface(sdf::AbstractArray{T, 3}, method::MarchingCubes, X=-1:1, Y=-
         # Cube is entirely in/out of the surface
         _no_triangles(cubeindex) && continue
 
-        points = mc_vert_points(xi+1,yi+1,zi+1,xp,yp,zp)
+        points = mc_vert_points(xi,yi,zi,xp,yp,zp)
 
         # process the voxel
         process_mc_voxel!(vts, fcs, cubeindex, points, method.iso, iso_vals)
@@ -39,7 +39,7 @@ function isosurface(sdf::AbstractArray{T, 3}, method::MarchingCubes, X=-1:1, Y=-
 end
 
 
-function isosurface(f::Function, method::MarchingCubes, X=-1:1, Y=-1:1, Z=-1:1; samples::NTuple{3,T}=_DEFAULT_SAMPLES) where {T <: Integer}
+function isosurface(f::F, method::MarchingCubes, X=-1:1, Y=-1:1, Z=-1:1; samples::NTuple{3,T}=_DEFAULT_SAMPLES) where {T <: Integer, F <: Function}
 
     nx, ny, nz = samples[1], samples[2], samples[3]
 
@@ -53,18 +53,40 @@ function isosurface(f::Function, method::MarchingCubes, X=-1:1, Y=-1:1, Z=-1:1;
     yp = LinRange(first(Y), last(Y), ny)
     zp = LinRange(first(Z), last(Z), nz)
 
-    @inbounds for xi = 1:nx-1, yi = 1:ny-1, zi = 1:nz-1
+    # initialize iso_vals so we can cache function evaluations
+    points = mc_vert_points(1,1,1,xp,yp,zp)
+    iso_vals = (f(points[1]),
+                f(points[2]),
+                f(points[3]),
+                f(points[4]),
+                f(points[5]),
+                f(points[6]),
+                f(points[7]),
+                f(points[8]))
 
-        points = mc_vert_points(xi+1,yi+1,zi+1,xp,yp,zp)
+    @inbounds for xi = 1:nx-1, yi = 1:ny-1, zi = 1:nz-1
 
-        iso_vals = (f(points[1]),
-                    f(points[2]),
-                    f(points[3]),
-                    f(points[4]),
-                    f(points[5]),
-                    f(points[6]),
-                    f(points[7]),
-                    f(points[8]))
+        points = mc_vert_points(xi,yi,zi,xp,yp,zp)
+
+        iso_vals = if xi == 1
+            (f(points[1]),
+            f(points[2]),
+            f(points[3]),
+            f(points[4]),
+            f(points[5]),
+            f(points[6]),
+            f(points[7]),
+            f(points[8]))
+        else
+            (iso_vals[2],
+            f(points[2]),
+            f(points[3]),
+            iso_vals[3],
+            iso_vals[6],
+            f(points[6]),
+            f(points[7]),
+            iso_vals[7])
+        end
 
         #Determine the index into the edge table which
         #tells us which vertices are inside of the surface
@@ -114,10 +136,8 @@ Linearly interpolate the position where an isosurface cuts
 an edge between two vertices, each with their own scalar value
 """
 function vertex_interp(iso, p1, p2, valp1, valp2)
-
-    mu = (iso .- valp1) / (valp2 .- valp1)
+    mu = (iso - valp1) / (valp2 - valp1)
     p = p1 .+ mu .* (p2 .- p1)
-
     return p
 end
 
@@ -127,12 +147,12 @@ end
 Returns a tuple of 8 points corresponding to each corner of a cube
 """
 @inline function mc_vert_points(xi,yi,zi,xp,yp,zp)
-    ((xp[xi-1],yp[yi-1],zp[zi-1]),
-    (xp[xi  ],yp[yi-1],zp[zi-1]),   
-    (xp[xi  ],yp[yi  ],zp[zi-1]),
-    (xp[xi-1],yp[yi  ],zp[zi-1]),
-    (xp[xi-1],yp[yi-1],zp[zi  ]),
-    (xp[xi  ],yp[yi-1],zp[zi  ]),
-    (xp[xi  ],yp[yi  ],zp[zi  ]),
-    (xp[xi-1],yp[yi  ],zp[zi  ]))
+    ((xp[xi  ],yp[yi  ],zp[zi  ]),
+     (xp[xi+1],yp[yi  ],zp[zi  ]),   
+     (xp[xi+1],yp[yi+1],zp[zi  ]),
+     (xp[xi  ],yp[yi+1],zp[zi  ]),
+     (xp[xi  ],yp[yi  ],zp[zi+1]),
+     (xp[xi+1],yp[yi  ],zp[zi+1]),
+     (xp[xi+1],yp[yi+1],zp[zi+1]),
+     (xp[xi  ],yp[yi+1],zp[zi+1]))
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index 96b9453..7f4a212 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,7 +1,6 @@
 using Meshing
 using Test
 using ForwardDiff
-using StaticArrays
 using Statistics: mean
 using LinearAlgebra: dot, norm
 using Random
@@ -36,19 +35,23 @@ end
     @test Meshing._get_cubeindex(iso_vals3, iso3) == 0xe0
 end
 
+@testset "vertex interpolation" begin
+    @test Meshing.vertex_interp(0, (0,0,0), (0,1,0), -1, 1) == (0,0.5,0)
+    @test Meshing.vertex_interp(-1, (0,0,0), (0,1,0), -1, 1) == (0,0,0)
+    @test Meshing.vertex_interp(1, (0,0,0), (0,1,0), -1, 1) == (0,1,0)
+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)
-    origin = SVector{3, Float64}(-1, -1, -1)
-    widths = SVector{3, Float64}(2, 2, 2)
     resolution = 0.1
 
     for algorithm in (MarchingCubes(0.5),
-                      #MarchingTetrahedra(0.5),
-                      MarchingCubes(iso=0.5, reduceverts=false))
-                      #MarchingTetrahedra(iso=0.5, reduceverts=false))
+                      MarchingTetrahedra(0.5),
+                      MarchingCubes(iso=0.5, reduceverts=false),
+                      MarchingTetrahedra(iso=0.5, reduceverts=false))
         @testset "$(typeof(algorithm))" begin
             # Extract isosurface using a function
             points, faces = isosurface(f, algorithm, samples=(50, 50, 50))
@@ -56,8 +59,11 @@ end
             # should be centered on the origin
             @test mean(collect.(points)) ≈ [0, 0, 0] atol=0.15*resolution
             # and should be symmetric about the origin
-            #@test maximum(points) ≈ [0.5, 0.5, 0.5]
-            #@test minimum(points) ≈ [-0.5, -0.5, -0.5]
+
+            for i in 1:3
+                @test maximum(map(p -> p[i], collect.(points))) ≈ 0.5 atol=0.001
+                @test minimum(map(p -> p[i], collect.(points))) ≈ -0.5 atol=0.001
+            end
         end
     end
 
@@ -66,7 +72,7 @@ end
     # so a larger tolerance is needed for this one. In addition,
     # quad -> triangle conversion is not functioning correctly
     # see: https://github.com/JuliaGeometry/GeometryTypes.jl/issues/169
-    points, faces = isosurface(f, NaiveSurfaceNets(0.5), origin=origin, widths=widths, samples=(21, 21, 21))
+    points, faces = isosurface(f, NaiveSurfaceNets(0.5), samples=(21, 21, 21))
 
     # should be centered on the origin
     #@test mean(points) ≈ [0, 0, 0] atol=0.15*resolution
@@ -115,12 +121,6 @@ end
     @test length(faces) == 6928
 end
 
-@testset "vertex interpolation" begin
-    @test Meshing.vertex_interp(0, SVector(0,0,0), SVector(0,1,0), -1, 1) == SVector(0,0.5,0)
-    @test Meshing.vertex_interp(-1, SVector(0,0,0), SVector(0,1,0), -1, 1) == SVector(0,0,0)
-    @test Meshing.vertex_interp(1, SVector(0,0,0), SVector(0,1,0), -1, 1) == SVector(0,1,0)
-end
-
 @testset "marching cubes" begin
     algo = MarchingCubes()
 

From da2a92d51f21d52408ceee9cd36e1b5b9bf00f06 Mon Sep 17 00:00:00 2001
From: Steve Kelly <kd2cca@gmail.com>
Date: Sun, 30 Jun 2024 17:45:19 -0400
Subject: [PATCH 3/8] fixup opts

---
 src/lut/mc.jl         | 46 +++++++++++++++++++-------------------
 src/marching_cubes.jl | 51 ++++++++++++-------------------------------
 2 files changed, 37 insertions(+), 60 deletions(-)

diff --git a/src/lut/mc.jl b/src/lut/mc.jl
index 42f0a59..86119c7 100644
--- a/src/lut/mc.jl
+++ b/src/lut/mc.jl
@@ -579,29 +579,29 @@ const _mc_eq_mapping = (0x01, 0x01, 0x02, 0x01, 0x03, 0x02, 0x04, 0x01, 0x02, 0x
                         0x11, 0x0c, 0x14, 0x11, 0x01, 0x0c, 0x11, 0x14, 0x01, 0x11, 0x01, 0x01)
 
 const _mc_connectivity = ((0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
-                            (0x04, 0x02, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
-                            (0x04, 0x05, 0x06, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
-                            (0x01, 0x04, 0x02, 0x04, 0x05, 0x02, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
-                            (0x03, 0x02, 0x04, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
-                            (0x01, 0x04, 0x02, 0x05, 0x06, 0x07, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
-                            (0x01, 0x03, 0x04, 0x01, 0x04, 0x05, 0x04, 0x03, 0x06, 0x00, 0x00, 0x00),
-                            (0x04, 0x05, 0x06, 0x07, 0x08, 0x09, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
-                            (0x04, 0x01, 0x03, 0x04, 0x03, 0x05, 0x04, 0x05, 0x06, 0x00, 0x00, 0x00),
-                            (0x01, 0x04, 0x02, 0x01, 0x05, 0x04, 0x06, 0x02, 0x04, 0x00, 0x00, 0x00),
-                            (0x04, 0x05, 0x06, 0x04, 0x06, 0x07, 0x06, 0x05, 0x08, 0x00, 0x00, 0x00),
-                            (0x01, 0x03, 0x04, 0x04, 0x03, 0x05, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
-                            (0x01, 0x03, 0x04, 0x01, 0x04, 0x05, 0x06, 0x04, 0x03, 0x00, 0x00, 0x00),
-                            (0x04, 0x05, 0x06, 0x05, 0x07, 0x06, 0x05, 0x08, 0x07, 0x00, 0x00, 0x00),
-                            (0x03, 0x02, 0x04, 0x05, 0x06, 0x07, 0x05, 0x07, 0x08, 0x00, 0x00, 0x00),
-                            (0x01, 0x03, 0x04, 0x02, 0x05, 0x03, 0x06, 0x03, 0x07, 0x05, 0x07, 0x03),
-                            (0x04, 0x01, 0x03, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
-                            (0x04, 0x05, 0x06, 0x07, 0x08, 0x09, 0x0a, 0x0b, 0x0c, 0x00, 0x00, 0x00),
-                            (0x01, 0x04, 0x02, 0x01, 0x05, 0x04, 0x06, 0x04, 0x05, 0x07, 0x08, 0x09),
-                            (0x01, 0x03, 0x04, 0x05, 0x06, 0x03, 0x06, 0x04, 0x03, 0x00, 0x00, 0x00),
-                            (0x01, 0x04, 0x02, 0x01, 0x05, 0x04, 0x04, 0x06, 0x02, 0x00, 0x00, 0x00),
-                            (0x04, 0x02, 0x01, 0x04, 0x05, 0x02, 0x04, 0x06, 0x05, 0x00, 0x00, 0x00),
-                            (0x01, 0x04, 0x02, 0x05, 0x06, 0x07, 0x06, 0x08, 0x07, 0x00, 0x00, 0x00),
-                            (0x02, 0x04, 0x03, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00))
+                        (0x04, 0x02, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
+                        (0x04, 0x05, 0x06, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
+                        (0x01, 0x04, 0x02, 0x04, 0x05, 0x02, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
+                        (0x03, 0x02, 0x04, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
+                        (0x01, 0x04, 0x02, 0x05, 0x06, 0x07, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
+                        (0x01, 0x03, 0x04, 0x01, 0x04, 0x05, 0x04, 0x03, 0x06, 0x00, 0x00, 0x00),
+                        (0x04, 0x05, 0x06, 0x07, 0x08, 0x09, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
+                        (0x04, 0x01, 0x03, 0x04, 0x03, 0x05, 0x04, 0x05, 0x06, 0x00, 0x00, 0x00),
+                        (0x01, 0x04, 0x02, 0x01, 0x05, 0x04, 0x06, 0x02, 0x04, 0x00, 0x00, 0x00),
+                        (0x04, 0x05, 0x06, 0x04, 0x06, 0x07, 0x06, 0x05, 0x08, 0x00, 0x00, 0x00),
+                        (0x01, 0x03, 0x04, 0x04, 0x03, 0x05, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
+                        (0x01, 0x03, 0x04, 0x01, 0x04, 0x05, 0x06, 0x04, 0x03, 0x00, 0x00, 0x00),
+                        (0x04, 0x05, 0x06, 0x05, 0x07, 0x06, 0x05, 0x08, 0x07, 0x00, 0x00, 0x00),
+                        (0x03, 0x02, 0x04, 0x05, 0x06, 0x07, 0x05, 0x07, 0x08, 0x00, 0x00, 0x00),
+                        (0x01, 0x03, 0x04, 0x02, 0x05, 0x03, 0x06, 0x03, 0x07, 0x05, 0x07, 0x03),
+                        (0x04, 0x01, 0x03, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
+                        (0x04, 0x05, 0x06, 0x07, 0x08, 0x09, 0x0a, 0x0b, 0x0c, 0x00, 0x00, 0x00),
+                        (0x01, 0x04, 0x02, 0x01, 0x05, 0x04, 0x06, 0x04, 0x05, 0x07, 0x08, 0x09),
+                        (0x01, 0x03, 0x04, 0x05, 0x06, 0x03, 0x06, 0x04, 0x03, 0x00, 0x00, 0x00),
+                        (0x01, 0x04, 0x02, 0x01, 0x05, 0x04, 0x04, 0x06, 0x02, 0x00, 0x00, 0x00),
+                        (0x04, 0x02, 0x01, 0x04, 0x05, 0x02, 0x04, 0x06, 0x05, 0x00, 0x00, 0x00),
+                        (0x01, 0x04, 0x02, 0x05, 0x06, 0x07, 0x06, 0x08, 0x07, 0x00, 0x00, 0x00),
+                        (0x02, 0x04, 0x03, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00))
 
 const _mc_edge_list = ((0x01, 0x02), (0x02, 0x03), (0x03, 0x04), (0x04, 0x01),
                         (0x05, 0x06), (0x06, 0x07), (0x07, 0x08), (0x08, 0x05),
diff --git a/src/marching_cubes.jl b/src/marching_cubes.jl
index 1035328..0196070 100644
--- a/src/marching_cubes.jl
+++ b/src/marching_cubes.jl
@@ -39,7 +39,7 @@ function isosurface(sdf::AbstractArray{T, 3}, method::MarchingCubes, X=-1:1, Y=-
 end
 
 
-function isosurface(f::F, method::MarchingCubes, X=-1:1, Y=-1:1, Z=-1:1; samples::NTuple{3,T}=_DEFAULT_SAMPLES) where {T <: Integer, F <: Function}
+function isosurface(f::Function, method::MarchingCubes, X=-1:1, Y=-1:1, Z=-1:1; samples::NTuple{3,T}=_DEFAULT_SAMPLES) where {T <: Integer}
 
     nx, ny, nz = samples[1], samples[2], samples[3]
 
@@ -53,40 +53,18 @@ function isosurface(f::F, method::MarchingCubes, X=-1:1, Y=-1:1, Z=-1:1; samples
     yp = LinRange(first(Y), last(Y), ny)
     zp = LinRange(first(Z), last(Z), nz)
 
-    # initialize iso_vals so we can cache function evaluations
-    points = mc_vert_points(1,1,1,xp,yp,zp)
-    iso_vals = (f(points[1]),
-                f(points[2]),
-                f(points[3]),
-                f(points[4]),
-                f(points[5]),
-                f(points[6]),
-                f(points[7]),
-                f(points[8]))
-
     @inbounds for xi = 1:nx-1, yi = 1:ny-1, zi = 1:nz-1
 
         points = mc_vert_points(xi,yi,zi,xp,yp,zp)
 
-        iso_vals = if xi == 1
-            (f(points[1]),
-            f(points[2]),
-            f(points[3]),
-            f(points[4]),
-            f(points[5]),
-            f(points[6]),
-            f(points[7]),
-            f(points[8]))
-        else
-            (iso_vals[2],
-            f(points[2]),
-            f(points[3]),
-            iso_vals[3],
-            iso_vals[6],
-            f(points[6]),
-            f(points[7]),
-            iso_vals[7])
-        end
+        iso_vals = (f(points[1]),
+                    f(points[2]),
+                    f(points[3]),
+                    f(points[4]),
+                    f(points[5]),
+                    f(points[6]),
+                    f(points[7]),
+                    f(points[8]))
 
         #Determine the index into the edge table which
         #tells us which vertices are inside of the surface
@@ -106,6 +84,7 @@ function process_mc_voxel!(vts, fcs, cubeindex, points, iso, iso_vals)
     fct = length(vts)
 
     @inbounds begin
+        # Add the vertices
         vert_to_add = _mc_verts[cubeindex]
         for i = 1:12
             vt = vert_to_add[i]
@@ -113,16 +92,14 @@ function process_mc_voxel!(vts, fcs, cubeindex, points, iso, iso_vals)
             ed = _mc_edge_list[vt]
             push!(vts, vertex_interp(iso,points[ed[1]],points[ed[2]],iso_vals[ed[1]],iso_vals[ed[2]]))
         end
-    end
-
-    @inbounds begin
 
+        # Add the faces
         offsets = _mc_connectivity[_mc_eq_mapping[cubeindex]]
 
         # There is atleast one face so we can push it immediately
-        push!(fcs, (fct+3, fct+3, fct+1))
+        push!(fcs, (fct+3, fct+2, fct+1))
 
-        for i in (1,4,7,10)
+        for i in (1, 4,7,10)
             iszero(offsets[i]) && return
             push!(fcs, (fct+offsets[i+2], fct+offsets[i+1], fct+offsets[i]))
         end
@@ -146,7 +123,7 @@ end
 
 Returns a tuple of 8 points corresponding to each corner of a cube
 """
-@inline function mc_vert_points(xi,yi,zi,xp,yp,zp)
+function mc_vert_points(xi,yi,zi,xp,yp,zp)
     ((xp[xi  ],yp[yi  ],zp[zi  ]),
      (xp[xi+1],yp[yi  ],zp[zi  ]),   
      (xp[xi+1],yp[yi+1],zp[zi  ]),

From 419980e379fe38005735c0931e8959317e5277d0 Mon Sep 17 00:00:00 2001
From: Steve Kelly <kd2cca@gmail.com>
Date: Sun, 30 Jun 2024 17:52:58 -0400
Subject: [PATCH 4/8] format

---
 src/marching_cubes.jl | 76 +++++++++++++++++++++----------------------
 1 file changed, 38 insertions(+), 38 deletions(-)

diff --git a/src/marching_cubes.jl b/src/marching_cubes.jl
index 0196070..9208883 100644
--- a/src/marching_cubes.jl
+++ b/src/marching_cubes.jl
@@ -2,11 +2,11 @@
 #Look up Table
 include("lut/mc.jl")
 
-function isosurface(sdf::AbstractArray{T, 3}, method::MarchingCubes, X=-1:1, Y=-1:1, Z=-1:1) where {T}
+function isosurface(sdf::AbstractArray{T,3}, method::MarchingCubes, X=-1:1, Y=-1:1, Z=-1:1) where {T}
     nx, ny, nz = size(sdf)
 
-    vts = NTuple{3, float(T)}[]
-    fcs = NTuple{3, Int}[]
+    vts = NTuple{3,float(T)}[]
+    fcs = NTuple{3,Int}[]
 
     xp = LinRange(first(X), last(X), nx)
     yp = LinRange(first(Y), last(Y), ny)
@@ -14,14 +14,14 @@ function isosurface(sdf::AbstractArray{T, 3}, method::MarchingCubes, X=-1:1, Y=-
 
     @inbounds for xi = 1:nx-1, yi = 1:ny-1, zi = 1:nz-1
 
-        iso_vals = (sdf[xi,yi,zi],
-                    sdf[xi+1,yi,zi],
-                    sdf[xi+1,yi+1,zi],
-                    sdf[xi,yi+1,zi],
-                    sdf[xi,yi,zi+1],
-                    sdf[xi+1,yi,zi+1],
-                    sdf[xi+1,yi+1,zi+1],
-                    sdf[xi,yi+1,zi+1])
+        iso_vals = (sdf[xi, yi, zi],
+            sdf[xi+1, yi, zi],
+            sdf[xi+1, yi+1, zi],
+            sdf[xi, yi+1, zi],
+            sdf[xi, yi, zi+1],
+            sdf[xi+1, yi, zi+1],
+            sdf[xi+1, yi+1, zi+1],
+            sdf[xi, yi+1, zi+1])
 
         #Determine the index into the edge table which
         #tells us which vertices are inside of the surface
@@ -30,24 +30,24 @@ function isosurface(sdf::AbstractArray{T, 3}, method::MarchingCubes, X=-1:1, Y=-
         # Cube is entirely in/out of the surface
         _no_triangles(cubeindex) && continue
 
-        points = mc_vert_points(xi,yi,zi,xp,yp,zp)
+        points = mc_vert_points(xi, yi, zi, xp, yp, zp)
 
         # process the voxel
         process_mc_voxel!(vts, fcs, cubeindex, points, method.iso, iso_vals)
     end
-    vts,fcs
+    vts, fcs
 end
 
 
-function isosurface(f::Function, method::MarchingCubes, X=-1:1, Y=-1:1, Z=-1:1; samples::NTuple{3,T}=_DEFAULT_SAMPLES) where {T <: Integer}
+function isosurface(f::Function, method::MarchingCubes, X=-1:1, Y=-1:1, Z=-1:1; samples::NTuple{3,T}=_DEFAULT_SAMPLES) where {T<:Integer}
 
     nx, ny, nz = samples[1], samples[2], samples[3]
 
     # find widest type
     FT = promote_type(eltype(first(X)), eltype(first(Y)), eltype(first(Z)), eltype(T))
 
-    vts = NTuple{3, float(FT)}[]
-    fcs = NTuple{3, Int}[]
+    vts = NTuple{3,float(FT)}[]
+    fcs = NTuple{3,Int}[]
 
     xp = LinRange(first(X), last(X), nx)
     yp = LinRange(first(Y), last(Y), ny)
@@ -55,16 +55,16 @@ function isosurface(f::Function, method::MarchingCubes, X=-1:1, Y=-1:1, Z=-1:1;
 
     @inbounds for xi = 1:nx-1, yi = 1:ny-1, zi = 1:nz-1
 
-        points = mc_vert_points(xi,yi,zi,xp,yp,zp)
+        points = mc_vert_points(xi, yi, zi, xp, yp, zp)
 
         iso_vals = (f(points[1]),
-                    f(points[2]),
-                    f(points[3]),
-                    f(points[4]),
-                    f(points[5]),
-                    f(points[6]),
-                    f(points[7]),
-                    f(points[8]))
+            f(points[2]),
+            f(points[3]),
+            f(points[4]),
+            f(points[5]),
+            f(points[6]),
+            f(points[7]),
+            f(points[8]))
 
         #Determine the index into the edge table which
         #tells us which vertices are inside of the surface
@@ -75,7 +75,7 @@ function isosurface(f::Function, method::MarchingCubes, X=-1:1, Y=-1:1, Z=-1:1;
 
         process_mc_voxel!(vts, fcs, cubeindex, points, method.iso, iso_vals)
     end
-    vts,fcs
+    vts, fcs
 end
 
 
@@ -90,18 +90,18 @@ function process_mc_voxel!(vts, fcs, cubeindex, points, iso, iso_vals)
             vt = vert_to_add[i]
             iszero(vt) && break
             ed = _mc_edge_list[vt]
-            push!(vts, vertex_interp(iso,points[ed[1]],points[ed[2]],iso_vals[ed[1]],iso_vals[ed[2]]))
+            push!(vts, vertex_interp(iso, points[ed[1]], points[ed[2]], iso_vals[ed[1]], iso_vals[ed[2]]))
         end
 
         # Add the faces
         offsets = _mc_connectivity[_mc_eq_mapping[cubeindex]]
 
         # There is atleast one face so we can push it immediately
-        push!(fcs, (fct+3, fct+2, fct+1))
+        push!(fcs, (fct + 3, fct + 2, fct + 1))
 
-        for i in (1, 4,7,10)
+        for i in (1, 4, 7, 10)
             iszero(offsets[i]) && return
-            push!(fcs, (fct+offsets[i+2], fct+offsets[i+1], fct+offsets[i]))
+            push!(fcs, (fct + offsets[i+2], fct + offsets[i+1], fct + offsets[i]))
         end
     end
 end
@@ -123,13 +123,13 @@ end
 
 Returns a tuple of 8 points corresponding to each corner of a cube
 """
-function mc_vert_points(xi,yi,zi,xp,yp,zp)
-    ((xp[xi  ],yp[yi  ],zp[zi  ]),
-     (xp[xi+1],yp[yi  ],zp[zi  ]),   
-     (xp[xi+1],yp[yi+1],zp[zi  ]),
-     (xp[xi  ],yp[yi+1],zp[zi  ]),
-     (xp[xi  ],yp[yi  ],zp[zi+1]),
-     (xp[xi+1],yp[yi  ],zp[zi+1]),
-     (xp[xi+1],yp[yi+1],zp[zi+1]),
-     (xp[xi  ],yp[yi+1],zp[zi+1]))
+function mc_vert_points(xi, yi, zi, xp, yp, zp)
+    ((xp[xi], yp[yi], zp[zi]),
+        (xp[xi+1], yp[yi], zp[zi]),
+        (xp[xi+1], yp[yi+1], zp[zi]),
+        (xp[xi], yp[yi+1], zp[zi]),
+        (xp[xi], yp[yi], zp[zi+1]),
+        (xp[xi+1], yp[yi], zp[zi+1]),
+        (xp[xi+1], yp[yi+1], zp[zi+1]),
+        (xp[xi], yp[yi+1], zp[zi+1]))
 end

From 7d730698e81be8dcbe2096ffb1bbcbb9c0b28e99 Mon Sep 17 00:00:00 2001
From: Steve Kelly <kd2cca@gmail.com>
Date: Sun, 30 Jun 2024 19:32:53 -0400
Subject: [PATCH 5/8] convert MT to new API

---
 src/lut/mt.jl              |  20 +++----
 src/marching_tetrahedra.jl | 116 ++++++++++++++++++-------------------
 test/runtests.jl           |  46 +++++++--------
 3 files changed, 88 insertions(+), 94 deletions(-)

diff --git a/src/lut/mt.jl b/src/lut/mt.jl
index d02d982..ad54018 100644
--- a/src/lut/mt.jl
+++ b/src/lut/mt.jl
@@ -20,18 +20,14 @@ Voxel corner and edge indexing conventions
   /
  X
 """
-# this gets vectorized so we want to ensure it is the
-# same type as out vertex
-function voxCrnrPos()
-    ((0, 0, 0),
-    (0, 1, 0),
-    (1, 1, 0),
-    (1, 0, 0),
-    (0, 0, 1),
-    (0, 1, 1),
-    (1, 1, 1),
-    (1, 0, 1))
-end
+const voxCrnrPos = ((0, 0, 0),
+                    (0, 1, 0),
+                    (1, 1, 0),
+                    (1, 0, 0),
+                    (0, 0, 1),
+                    (0, 1, 1),
+                    (1, 1, 1),
+                    (1, 0, 1))
 
 # the voxel IDs at either end of the tetrahedra edges, by edge ID
 const voxEdgeCrnrs = ((0x01, 0x02),
diff --git a/src/marching_tetrahedra.jl b/src/marching_tetrahedra.jl
index aacccc5..278b259 100644
--- a/src/marching_tetrahedra.jl
+++ b/src/marching_tetrahedra.jl
@@ -8,7 +8,7 @@ include("lut/mt.jl")
 
 Determines which case in the triangle table we are dealing with
 """
-@inline function tetIx(tIx, cubeindex)
+function tetIx(tIx, cubeindex)
     @inbounds v1 = subTetsMask[tIx][1]
     @inbounds v2 = subTetsMask[tIx][2]
     ix = 0x01 & cubeindex | (0x40 & cubeindex) >> 0x03
@@ -26,8 +26,8 @@ two edges get the same index) and unique (every edge gets the same ID
 regardless of which of its neighboring voxels is asking for it) in order
 for vertex sharing to be implemented properly.
 """
-@inline function vertId(e, x, y, z, nx, ny)
-    dx, dy, dz = voxCrnrPosInt[voxEdgeCrnrs[e][1]]
+function vertId(e, x, y, z, nx, ny)
+    dx, dy, dz = voxCrnrPos[voxEdgeCrnrs[e][1]]
     voxEdgeDir[e]+7*(x-0x01+dx+nx*(y-0x01+dy+ny*(z-0x01+dz)))
 end
 
@@ -39,7 +39,7 @@ occurs.
 eps represents the "bump" factor to keep vertices away from voxel
 corners (thereby preventing degeneracies).
 """
-@inline function vertPos(e, x, y, z, scale, origin, vals::V, iso, eps, ::Type{VertType}) where {V, VertType}
+function vertPos(e, x, y, z, xp, yp, zp, vals::V, iso, eps) where {V}
     T = eltype(vals)
 
     ixs     = voxEdgeCrnrs[e]
@@ -47,10 +47,11 @@ corners (thereby preventing degeneracies).
     tgtVal  = vals[ixs[2]]
     a       = min(max((iso-srcVal)/(tgtVal-srcVal), eps), one(T)-eps)
     b       = one(T)-a
-    c1 = voxCrnrPos(VertType)[ixs[1]]
-    c2 = voxCrnrPos(VertType)[ixs[2]]
+    c1 = voxCrnrPos[ixs[1]]
+    c2 = voxCrnrPos[ixs[2]]
 
-    ((VertType(x,y,z) + c1 .* b + c2.* a) .- 1) .* scale .+ origin
+    d = (xp[x+1], yp[y+1], zp[z+1]) .- (xp[x], yp[y], zp[z])
+    (xp[x],yp[y],zp[z]) .+ (c1 .* b .+ c2 .* a) .* d
 end
 
 """
@@ -63,28 +64,21 @@ end
 Gets the vertex ID, adding it to the vertex dictionary if not already
 present.
 """
-@inline function getVertId(e, x, y, z, nx, ny,
+function getVertId(e, x, y, z, nx, ny,
                            vals, iso::Real,
-                           scale, origin,
+                           xp, yp, zp,
                            vts::Dict,
                            vtsAry::Vector,
-                           eps::Real,
-                           reduceverts::Bool)
+                           eps::Real)
 
-    VertType = eltype(vtsAry)
-    if reduceverts
-        vId = vertId(e, x, y, z, nx, ny)
-        haskey(vts, vId) && return vts[vId]
-    end
+    vId = vertId(e, x, y, z, nx, ny)
+    haskey(vts, vId) && return vts[vId]
 
     # calculate vert position
-    v = vertPos(e, x, y, z, scale, origin, vals, iso, eps, VertType)
+    v = vertPos(e, x, y, z, xp, yp, zp, vals, iso, eps)
     push!(vtsAry, v)
 
-    # if deduplicting, push to dict
-    if reduceverts
-        vts[vId] = length(vtsAry)
-    end
+    vts[vId] = length(vtsAry)
 
     return length(vtsAry)
 end
@@ -94,7 +88,7 @@ end
 Given a sub-tetrahedron case and a tetrahedron edge ID, determines the
 corresponding voxel edge ID.
 """
-@inline function voxEdgeId(subTetIx, tetEdgeIx)
+function voxEdgeId(subTetIx, tetEdgeIx)
     srcVoxCrnr = subTets[subTetIx][tetEdgeCrnrs[tetEdgeIx][1]]
     tgtVoxCrnr = subTets[subTetIx][tetEdgeCrnrs[tetEdgeIx][2]]
     return voxEdgeIx[srcVoxCrnr][tgtVoxCrnr]
@@ -108,11 +102,10 @@ end
 Processes a voxel, adding any new vertices and faces to the given
 containers as necessary.
 """
-function procVox(vals, iso::Real, x, y, z, nx, ny, scale, origin,
+function procVox(vals, iso::Real, x, y, z, nx, ny, xp, yp, zp,
                  vts::Dict, vtsAry::Vector, fcs::Vector,
-                 eps::Real, cubeindex, reduceverts)
-    VertType = eltype(vtsAry)
-    FaceType = eltype(fcs)
+                 eps::Real, cubeindex)
+
     # check each sub-tetrahedron in the voxel
     @inbounds for i = 1:6
         tIx = tetIx(i, cubeindex)
@@ -121,31 +114,31 @@ function procVox(vals, iso::Real, x, y, z, nx, ny, scale, origin,
         e = tetTri[tIx]
 
         # add the face to the list
-        push!(fcs, FaceType(
-                    getVertId(voxEdgeId(i, e[1]), x, y, z, nx, ny, vals, iso, scale, origin, vts, vtsAry, eps, reduceverts),
-                    getVertId(voxEdgeId(i, e[2]), x, y, z, nx, ny, vals, iso, scale, origin, vts, vtsAry, eps, reduceverts),
-                    getVertId(voxEdgeId(i, e[3]), x, y, z, nx, ny, vals, iso, scale, origin, vts, vtsAry, eps, reduceverts)))
+        push!(fcs, (getVertId(voxEdgeId(i, e[1]), x, y, z, nx, ny, vals, iso, xp, yp, zp, vts, vtsAry, eps),
+                    getVertId(voxEdgeId(i, e[2]), x, y, z, nx, ny, vals, iso, xp, yp, zp, vts, vtsAry, eps),
+                    getVertId(voxEdgeId(i, e[3]), x, y, z, nx, ny, vals, iso, xp, yp, zp, vts, vtsAry, eps)))
 
         # bail if there are no more faces
         iszero(e[4]) && continue
-        push!(fcs, FaceType(
-                    getVertId(voxEdgeId(i, e[4]), x, y, z, nx, ny, vals, iso, scale, origin, vts, vtsAry, eps, reduceverts),
-                    getVertId(voxEdgeId(i, e[5]), x, y, z, nx, ny, vals, iso, scale, origin, vts, vtsAry, eps, reduceverts),
-                    getVertId(voxEdgeId(i, e[6]), x, y, z, nx, ny, vals, iso, scale, origin, vts, vtsAry, eps, reduceverts)))
+        push!(fcs, (getVertId(voxEdgeId(i, e[4]), x, y, z, nx, ny, vals, iso, xp, yp, zp, vts, vtsAry, eps),
+                    getVertId(voxEdgeId(i, e[5]), x, y, z, nx, ny, vals, iso, xp, yp, zp, vts, vtsAry, eps),
+                    getVertId(voxEdgeId(i, e[6]), x, y, z, nx, ny, vals, iso, xp, yp, zp, vts, vtsAry, eps)))
     end
 end
 
-function isosurface(sdf::AbstractArray{T, 3}, method::MarchingTetrahedra, ::Type{VertType}=SVector{3,Float64}, ::Type{FaceType}=SVector{3, Int};
-                    origin=VertType(-1,-1,-1), widths=VertType(2,2,2)) where {T, VertType, FaceType}
+function isosurface(sdf::AbstractArray{T, 3}, method::MarchingTetrahedra, X=-1:1, Y=-1:1, Z=-1:1) where {T}
 
     vts    = Dict{Int, Int}()
-    fcs    = FaceType[]
-    vtsAry = VertType[]
+    fcs = NTuple{3,Int}[]
+    vtsAry = NTuple{3,float(T)}[]
 
     # process each voxel
-    scale = widths ./ VertType(size(sdf) .- 1)
     nx::Int, ny::Int, nz::Int = size(sdf)
 
+    xp = LinRange(first(X), last(X), nx)
+    yp = LinRange(first(Y), last(Y), ny)
+    zp = LinRange(first(Z), last(Z), nz)
+
     @inbounds for i = 1:nx-1, j = 1:ny-1, k = 1:nz-1
 
         vals = (sdf[i,  j,  k  ],
@@ -161,36 +154,30 @@ function isosurface(sdf::AbstractArray{T, 3}, method::MarchingTetrahedra, ::Type
 
         _no_triangles(cubeindex) && continue
 
-        procVox(vals, method.iso, i, j, k, nx, ny, scale, origin, vts, vtsAry, fcs, method.eps, cubeindex, method.reduceverts)
+        procVox(vals, method.iso, i, j, k, nx, ny, xp, yp, zp, vts, vtsAry, fcs, method.eps, cubeindex)
     end
 
     vtsAry,fcs
 end
 
-function isosurface(f::Function, method::MarchingTetrahedra,
-                    ::Type{VertType}=SVector{3,Float64}, ::Type{FaceType}=SVector{3, Int};
-                    origin=VertType(-1,-1,-1), widths=VertType(2,2,2),
-                    samples::NTuple{3,T}=_DEFAULT_SAMPLES) where {T, VertType, FaceType}
+function isosurface(f::F, method::MarchingTetrahedra, X=-1:1, Y=-1:1, Z=-1:1;
+                    samples::NTuple{3,T}=_DEFAULT_SAMPLES) where {F, T}
+
+    FT = promote_type(eltype(first(X)), eltype(first(Y)), eltype(first(Z)), eltype(T))
 
     vts    = Dict{Int, Int}()
-    fcs    = FaceType[]
-    vtsAry = VertType[]
+    fcs    = NTuple{3,Int}[]
+    vtsAry = NTuple{3,float(FT)}[]
 
     # process each voxel
-    scale = widths ./ (VertType(samples...) .-1)
-    nx::Int, ny::Int, nz::Int = samples
-    zv = zero(eltype(VertType))
-    vals = (zv,zv,zv,zv,zv,zv,zv,zv)
+    nx, ny, nz = samples[1], samples[2], samples[3]
+    xp = LinRange(first(X), last(X), nx)
+    yp = LinRange(first(Y), last(Y), ny)
+    zp = LinRange(first(Z), last(Z), nz)
 
     @inbounds for i = 1:nx-1, j = 1:ny-1, k = 1:nz-1
-        points = (VertType(i-1,j-1,k-1).* scale .+ origin,
-                  VertType(i-1,j  ,k-1).* scale .+ origin,
-                  VertType(i  ,j  ,k-1).* scale .+ origin,
-                  VertType(i  ,j-1,k-1).* scale .+ origin,
-                  VertType(i-1,j-1,k  ).* scale .+ origin,
-                  VertType(i-1,j  ,k  ).* scale .+ origin,
-                  VertType(i  ,j  ,k  ).* scale .+ origin,
-                  VertType(i  ,j-1,k  ).* scale .+ origin)
+        points = mt_vert_points(i, j, k, xp, yp, zp)
+
         vals = (f(points[1]),
                 f(points[2]),
                 f(points[3]),
@@ -204,8 +191,19 @@ function isosurface(f::Function, method::MarchingTetrahedra,
 
         _no_triangles(cubeindex) && continue
 
-        procVox(vals, method.iso, i, j, k, nx, ny, scale, origin, vts, vtsAry, fcs, method.eps, cubeindex, method.reduceverts)
+        procVox(vals, method.iso, i, j, k, nx, ny, xp, yp, zp, vts, vtsAry, fcs, method.eps, cubeindex)
     end
 
     vtsAry,fcs
 end
+
+function mt_vert_points(i, j, k, xp, yp, zp)
+    ((xp[i], yp[j], zp[k]),
+    (xp[i], yp[j+1], zp[k]),
+    (xp[i+1], yp[j+1], zp[k]),
+    (xp[i+1], yp[j], zp[k]),
+    (xp[i], yp[j], zp[k+1]),
+    (xp[i], yp[j+1], zp[k+1]),
+    (xp[i+1], yp[j+1], zp[k+1]),
+    (xp[i+1], yp[j], zp[k+1]))
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 7f4a212..71fe9e2 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -72,7 +72,7 @@ end
     # so a larger tolerance is needed for this one. In addition,
     # quad -> triangle conversion is not functioning correctly
     # see: https://github.com/JuliaGeometry/GeometryTypes.jl/issues/169
-    points, faces = isosurface(f, NaiveSurfaceNets(0.5), samples=(21, 21, 21))
+    #points, faces = isosurface(f, NaiveSurfaceNets(0.5), samples=(21, 21, 21))
 
     # should be centered on the origin
     #@test mean(points) ≈ [0, 0, 0] atol=0.15*resolution
@@ -80,24 +80,24 @@ end
     #@test maximum(points) ≈ [0.5, 0.5, 0.5] atol=0.2
     #@test minimum(points) ≈ [-0.5, -0.5, -0.5] atol=0.2
 end
-@testset "surface nets" begin
-
-
-    """
-    Test the isosurface function with NaiveSurfaceNets and a function.
-
-    This code tests the `isosurface` function with the `NaiveSurfaceNets` algorithm and two different level set functions: `sphere_function` and `torus_function`. It checks that the number of vertices and faces in the resulting meshes match expected values.
-    """
-    # Test the isosurface function with NaiveSurfaceNets and a function
-    points_sphere, faces_sphere = isosurface(sphere_function, NaiveSurfaceNets(), origin=SVector(-1, -1, -1), widths=SVector(2, 2, 2), samples=(21, 21, 21))
-    points_torus, faces_torus = isosurface(torus_function, NaiveSurfaceNets(), origin=SVector(-2, -2, -2), widths=SVector(4, 4, 4), samples=(81, 81, 81))
-
-    # Add assertions to check the output
-    @test length(points_sphere) == 1832
-    @test length(faces_sphere) == 1830
-    @test length(points_torus) == 5532
-    @test length(faces_torus) == 5532
-end
+#@testset "surface nets" begin
+#
+#
+#    """
+#    Test the isosurface function with NaiveSurfaceNets and a function.
+#
+#    This code tests the `isosurface` function with the `NaiveSurfaceNets` algorithm and two different level set functions: `sphere_function` and `torus_function`. It checks that the number of vertices and faces in the resulting meshes match expected values.
+#    """
+#    # Test the isosurface function with NaiveSurfaceNets and a function
+#    points_sphere, faces_sphere = isosurface(sphere_function, NaiveSurfaceNets(), origin=SVector(-1, -1, -1), widths=SVector(2, 2, 2), samples=(21, 21, 21))
+#    points_torus, faces_torus = isosurface(torus_function, NaiveSurfaceNets(), origin=SVector(-2, -2, -2), widths=SVector(4, 4, 4), samples=(81, 81, 81))
+#
+#    # Add assertions to check the output
+#    @test length(points_sphere) == 1832
+#    @test length(faces_sphere) == 1830
+#    @test length(points_torus) == 5532
+#    @test length(faces_torus) == 5532
+#end
 
 @testset "noisy spheres" begin
     # Produce a level set function that is a noisy version of the distance from
@@ -136,8 +136,8 @@ end
     a2 = MarchingTetrahedra(reduceverts=false)
 
     # Extract isosurfaces using MarchingTetrahedra
-    points_mt1, faces_mt1 = isosurface(v -> sqrt(sum(dot(v, v))) - 1, a1, origin=SVector(-1, -1, -1), widths=SVector(2, 2, 2), samples=(21, 21, 21))
-    points_mt2, faces_mt2 = isosurface(v -> sqrt(sum(dot(v, v))) - 1, a2, origin=SVector(-1, -1, -1), widths=SVector(2, 2, 2), samples=(21, 21, 21))
+    points_mt1, faces_mt1 = isosurface(v -> sqrt(sum(dot(v, v))) - 1, a1, samples=(21, 21, 21))
+    points_mt2, faces_mt2 = isosurface(v -> sqrt(sum(dot(v, v))) - 1, a2, samples=(21, 21, 21))
 
     # Test the number of vertices and faces
     @test length(points_mt1) == 5582
@@ -146,8 +146,8 @@ end
     @test length(faces_mt2) == length(faces_mt1)
 
     # When zero set is not completely within the box
-    points_mt1_partial, faces_mt1_partial = isosurface(v -> v[3] - 50, a1, origin=SVector(0, 0, 40), widths=SVector(50, 50, 60), samples=(2, 2, 2))
-    points_mt2_partial, faces_mt2_partial = isosurface(v -> v[3] - 50, a2, origin=SVector(0, 0, 40), widths=SVector(50, 50, 60), samples=(2, 2, 2))
+    points_mt1_partial, faces_mt1_partial = isosurface(v -> v[3] - 50, a1, 0:50, 0:50, 40:60, samples=(2, 2, 2))
+    points_mt2_partial, faces_mt2_partial = isosurface(v -> v[3] - 50, a2, 0:50, 0:50, 40:60, samples=(2, 2, 2))
 
     # Test the number of vertices and faces
     @test length(points_mt1_partial) == 9

From 2f10b50f8ce831f45509fc9162ae069fa57dbaeb Mon Sep 17 00:00:00 2001
From: Steve Kelly <kd2cca@gmail.com>
Date: Sun, 30 Jun 2024 19:46:56 -0400
Subject: [PATCH 6/8] update tests

---
 src/algorithmtypes.jl | 53 ++++++++++---------------------------------
 test/runtests.jl      | 30 ++++--------------------
 2 files changed, 17 insertions(+), 66 deletions(-)

diff --git a/src/algorithmtypes.jl b/src/algorithmtypes.jl
index d1b81dc..1a00235 100644
--- a/src/algorithmtypes.jl
+++ b/src/algorithmtypes.jl
@@ -13,27 +13,8 @@ abstract type AbstractMeshingAlgorithm end
 abstract type AbstractAdaptiveMeshingAlgorithm end
 
 
-function (::Type{MeshAlgo})(;iso::T1=0.0, eps::T2=1e-3, reduceverts::Bool=true, insidepositive::Bool=false) where {T1, T2, MeshAlgo <: AbstractMeshingAlgorithm}
-    if isconcretetype(MeshAlgo)
-        return MeshAlgo(iso, eps, reduceverts, insidepositive)
-    else
-        return MeshAlgo{promote_type(T1, T2)}(iso, eps, reduceverts, insidepositive)
-    end
-end
-
-function (::Type{MeshAlgo})(iso) where MeshAlgo <: AbstractMeshingAlgorithm
-    MeshAlgo(iso=iso)
-end
-
-function (::Type{MeshAlgo})(iso,eps) where MeshAlgo <: AbstractMeshingAlgorithm
-    MeshAlgo(iso=iso,eps=eps)
-end
-
 """
-    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)
+    MarchingCubes(iso=0.0)
 
 Specifies the use of the Marching Cubes algorithm for isosurface extraction.
 This algorithm provides a good balance between performance and vertex count.
@@ -41,19 +22,14 @@ In contrast to the other algorithms, vertices may be repeated, so mesh size
 may be large and it will be difficult to extract topological/connectivity information.
 
 - `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 slight performance improvement.
 """
-struct MarchingCubes{T} <: AbstractMeshingAlgorithm
-    iso::T
-    eps::T
-    reduceverts::Bool
-    insidepositive::Bool
+Base.@kwdef struct MarchingCubes{T} <: AbstractMeshingAlgorithm
+    iso::T = 0.0
 end
 
 """
-    MarchingTetrahedra(iso=0.0, eps=1e-3, reduceverts=true, insidepositive=false)
-    MarchingTetrahedra(;iso=0.0, eps=1e-3, reduceverts=true, insidepositive=false)
+    MarchingTetrahedra(iso=0.0, eps=1e-3)
+    MarchingTetrahedra(;iso=0.0, eps=1e-3)
     MarchingTetrahedra(iso)
     MarchingTetrahedra(iso,eps)
 
@@ -64,18 +40,15 @@ 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` (default: true) if false, vertices will not be unique and have repeated face references.
 """
-struct MarchingTetrahedra{T} <: AbstractMeshingAlgorithm
-    iso::T
-    eps::T
-    reduceverts::Bool
-    insidepositive::Bool
+Base.@kwdef struct MarchingTetrahedra{T} <: AbstractMeshingAlgorithm
+    iso::T = 0.0
+    eps::T = 1e-3
 end
 
 """
-    NaiveSurfaceNets(iso=0.0, eps=1e-3, reduceverts=true, insidepositive=false)
-    NaiveSurfaceNets(;iso=0.0, eps=1e-3, reduceverts=true, insidepositive=false)
+    NaiveSurfaceNets(iso=0.0, eps=1e-3)
+    NaiveSurfaceNets(;iso=0.0, eps=1e-3)
     NaiveSurfaceNets(iso)
     NaiveSurfaceNets(iso,eps)
 
@@ -96,8 +69,8 @@ struct NaiveSurfaceNets{T} <: AbstractMeshingAlgorithm
 end
 
 """
-    MarchingCubes(iso=0.0, eps=1e-3, reduceverts=true, insidepositive=false)
-    MarchingCubes(;iso=0.0, eps=1e-3, reduceverts=true, insidepositive=false)
+    MarchingCubes(iso=0.0, eps=1e-3)
+    MarchingCubes(;iso=0.0, eps=1e-3)
     MarchingCubes(iso)
     MarchingCubes(iso,eps)
 
@@ -115,7 +88,6 @@ struct AdaptiveMarchingCubes{T} <: AbstractAdaptiveMeshingAlgorithm
     eps::T
     rtol::T
     atol::T
-    reduceverts::Bool
 end
 
 
@@ -124,7 +96,6 @@ struct AdaptiveMarchingTetrahedra{T} <: AbstractAdaptiveMeshingAlgorithm
     eps::T
     rtol::T
     atol::T
-    reduceverts::Bool
 end
 
 
diff --git a/test/runtests.jl b/test/runtests.jl
index 71fe9e2..c02bb76 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -9,17 +9,6 @@ const algos = (MarchingCubes, MarchingTetrahedra, NaiveSurfaceNets)
 
 sphere_function(v) = sqrt(sum(dot(v, v))) - 1
 torus_function(v) = (sqrt(v[1]^2 + v[2]^2) - 0.5)^2 + v[3]^2 - 0.25
-@testset "meshing algorithms" begin
-    for algo in (MarchingCubes, MarchingTetrahedra, NaiveSurfaceNets)
-        @test algo(5) == algo{Float64}(5.0, 0.001, true, false)
-        @test algo(0x00, 0x00) == algo{UInt8}(0x00, 0x00, true, false)
-        @test algo{UInt16}(5, 0x00) == algo{UInt16}(0x0005, 0x0000, true, false)
-        @test algo{Float32}(1) ==  algo{Float32}(1.0f0, 0.001f0, true, false)
-        @test algo{Float32}(1.0, 0x00) == algo{Float32}(1.0f0, 0.0f0, true, false)
-        @test algo{Float32}(iso=1, eps=0x00, insidepositive=true, reduceverts=false) ==
-            algo{Float32}(1.0f0, 0.0f0, false, true)
-    end
-end
 
 @testset "_get_cubeindex" begin
     iso_vals1 = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8]
@@ -48,10 +37,8 @@ end
     f = x -> norm(x)
     resolution = 0.1
 
-    for algorithm in (MarchingCubes(0.5),
-                      MarchingTetrahedra(0.5),
-                      MarchingCubes(iso=0.5, reduceverts=false),
-                      MarchingTetrahedra(iso=0.5, reduceverts=false))
+    for algorithm in (MarchingCubes(iso=0.5),
+                      MarchingTetrahedra(iso=0.5))
         @testset "$(typeof(algorithm))" begin
             # Extract isosurface using a function
             points, faces = isosurface(f, algorithm, samples=(50, 50, 50))
@@ -115,7 +102,7 @@ end
     #
     lambda = N-2*sigma # isovalue
 
-    points, faces = isosurface(distance, MarchingTetrahedra(lambda))
+    points, faces = isosurface(distance, MarchingTetrahedra(iso=lambda))
 
     @test length(points) == 3466
     @test length(faces) == 6928
@@ -133,27 +120,20 @@ end
 end
 @testset "marching tetrahedra" begin
     a1 = MarchingTetrahedra()
-    a2 = MarchingTetrahedra(reduceverts=false)
 
     # Extract isosurfaces using MarchingTetrahedra
     points_mt1, faces_mt1 = isosurface(v -> sqrt(sum(dot(v, v))) - 1, a1, samples=(21, 21, 21))
-    points_mt2, faces_mt2 = isosurface(v -> sqrt(sum(dot(v, v))) - 1, a2, samples=(21, 21, 21))
 
     # Test the number of vertices and faces
     @test length(points_mt1) == 5582
-    @test length(points_mt2) == 33480
     @test length(faces_mt1) == 11160
-    @test length(faces_mt2) == length(faces_mt1)
 
     # When zero set is not completely within the box
     points_mt1_partial, faces_mt1_partial = isosurface(v -> v[3] - 50, a1, 0:50, 0:50, 40:60, samples=(2, 2, 2))
-    points_mt2_partial, faces_mt2_partial = isosurface(v -> v[3] - 50, a2, 0:50, 0:50, 40:60, samples=(2, 2, 2))
 
     # Test the number of vertices and faces
     @test length(points_mt1_partial) == 9
-    @test length(points_mt2_partial) == 24
     @test length(faces_mt1_partial) == 8
-    @test length(faces_mt2_partial) == length(faces_mt1_partial)
 end
 
 @testset "function/array" begin
@@ -161,11 +141,11 @@ end
     for algo in algos
         @testset "$algo" begin
             # Extract isosurface using a function
-            points_fn, faces_fn = isosurface(torus_function, algo(), origin=SVector(-2, -2, -2), widths=SVector(4, 4, 4), samples=(81, 81, 81))
+            points_fn, faces_fn = isosurface(torus_function, algo(), -2:2, -2:2, -2:2, samples=(81, 81, 81))
 
             # Extract isosurface using an array
             torus_array = [torus_function(SVector(x, y, z)) for x in -2:0.05:2, y in -2:0.05:2, z in -2:0.05:2]
-            points_arr, faces_arr = isosurface(torus_array, algo(), origin=SVector(-2, -2, -2), widths=SVector(4, 4, 4))
+            points_arr, faces_arr = isosurface(torus_array, algo(), -2:2, -2:2, -2:2)
 
             # Test that the vertices and faces are the same for both cases
             @test_broken all(points_fn .≈ points_arr)

From 34ada529dee9c4e6eb630d9682291f88f4390f39 Mon Sep 17 00:00:00 2001
From: Steve Kelly <kd2cca@gmail.com>
Date: Mon, 1 Jul 2024 21:24:10 -0400
Subject: [PATCH 7/8] remove NaiveSurfaceNets and cleanup tests (support
 autodiff again)

---
 src/Meshing.jl             |   4 +-
 src/algorithmtypes.jl      |  29 +---
 src/marching_cubes.jl      |   7 +-
 src/marching_tetrahedra.jl |   6 +-
 src/surface_nets.jl        | 270 -------------------------------------
 test/runtests.jl           | 100 ++++++--------
 6 files changed, 51 insertions(+), 365 deletions(-)
 delete mode 100644 src/surface_nets.jl

diff --git a/src/Meshing.jl b/src/Meshing.jl
index cce2316..dbde984 100644
--- a/src/Meshing.jl
+++ b/src/Meshing.jl
@@ -11,13 +11,11 @@ include("algorithmtypes.jl")
 include("common.jl")
 include("marching_tetrahedra.jl")
 include("marching_cubes.jl")
-include("surface_nets.jl")
 include("roots.jl")
 include("adaptive.jl")
 
 export isosurface,
        MarchingCubes,
-       MarchingTetrahedra,
-       NaiveSurfaceNets
+       MarchingTetrahedra
 
 end # module
diff --git a/src/algorithmtypes.jl b/src/algorithmtypes.jl
index 1a00235..57281d4 100644
--- a/src/algorithmtypes.jl
+++ b/src/algorithmtypes.jl
@@ -41,31 +41,9 @@ 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.
 """
-Base.@kwdef struct MarchingTetrahedra{T} <: AbstractMeshingAlgorithm
+Base.@kwdef struct MarchingTetrahedra{T,E} <: AbstractMeshingAlgorithm
     iso::T = 0.0
-    eps::T = 1e-3
-end
-
-"""
-    NaiveSurfaceNets(iso=0.0, eps=1e-3)
-    NaiveSurfaceNets(;iso=0.0, eps=1e-3)
-    NaiveSurfaceNets(iso)
-    NaiveSurfaceNets(iso,eps)
-
-Specifies the use of the Naive Surface Nets algorithm for isosurface extraction.
-This algorithm has a slight performance advantage in some cases over Marching Cubes.
-Each vertex is guaranteed to not be repeated (useful for topological analysis),
-however the algorithm does not 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.
-"""
-struct NaiveSurfaceNets{T} <: AbstractMeshingAlgorithm
-    iso::T
-    eps::T
-    reduceverts::Bool
-    insidepositive::Bool
+    eps::E = 1e-3
 end
 
 """
@@ -83,7 +61,7 @@ may be large and it will be difficult to extract topological/connectivity inform
 - `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 slight performance improvement.
 """
-struct AdaptiveMarchingCubes{T} <: AbstractAdaptiveMeshingAlgorithm
+struct AdaptiveMarchingCubes{T,E} <: AbstractAdaptiveMeshingAlgorithm
     iso::T
     eps::T
     rtol::T
@@ -104,4 +82,3 @@ end
 #
 
 default_face_length(::Union{MarchingCubes,MarchingTetrahedra}) = 3
-default_face_length(::NaiveSurfaceNets) = 4
diff --git a/src/marching_cubes.jl b/src/marching_cubes.jl
index 9208883..4794048 100644
--- a/src/marching_cubes.jl
+++ b/src/marching_cubes.jl
@@ -5,7 +5,10 @@ include("lut/mc.jl")
 function isosurface(sdf::AbstractArray{T,3}, method::MarchingCubes, X=-1:1, Y=-1:1, Z=-1:1) where {T}
     nx, ny, nz = size(sdf)
 
-    vts = NTuple{3,float(T)}[]
+    # find widest type
+    FT = promote_type(eltype(first(X)), eltype(first(Y)), eltype(first(Z)), eltype(T), typeof(method.iso))
+
+    vts = NTuple{3,float(FT)}[]
     fcs = NTuple{3,Int}[]
 
     xp = LinRange(first(X), last(X), nx)
@@ -44,7 +47,7 @@ function isosurface(f::Function, method::MarchingCubes, X=-1:1, Y=-1:1, Z=-1:1;
     nx, ny, nz = samples[1], samples[2], samples[3]
 
     # find widest type
-    FT = promote_type(eltype(first(X)), eltype(first(Y)), eltype(first(Z)), eltype(T))
+    FT = promote_type(eltype(first(X)), eltype(first(Y)), eltype(first(Z)), eltype(T), typeof(method.iso))
 
     vts = NTuple{3,float(FT)}[]
     fcs = NTuple{3,Int}[]
diff --git a/src/marching_tetrahedra.jl b/src/marching_tetrahedra.jl
index 278b259..faffc4a 100644
--- a/src/marching_tetrahedra.jl
+++ b/src/marching_tetrahedra.jl
@@ -128,9 +128,11 @@ end
 
 function isosurface(sdf::AbstractArray{T, 3}, method::MarchingTetrahedra, X=-1:1, Y=-1:1, Z=-1:1) where {T}
 
+    FT = promote_type(eltype(first(X)), eltype(first(Y)), eltype(first(Z)), eltype(T), typeof(method.iso), typeof(method.eps))
+
     vts    = Dict{Int, Int}()
     fcs = NTuple{3,Int}[]
-    vtsAry = NTuple{3,float(T)}[]
+    vtsAry = NTuple{3,float(FT)}[]
 
     # process each voxel
     nx::Int, ny::Int, nz::Int = size(sdf)
@@ -163,7 +165,7 @@ end
 function isosurface(f::F, method::MarchingTetrahedra, X=-1:1, Y=-1:1, Z=-1:1;
                     samples::NTuple{3,T}=_DEFAULT_SAMPLES) where {F, T}
 
-    FT = promote_type(eltype(first(X)), eltype(first(Y)), eltype(first(Z)), eltype(T))
+    FT = promote_type(eltype(first(X)), eltype(first(Y)), eltype(first(Z)), eltype(T), typeof(method.iso), typeof(method.eps))
 
     vts    = Dict{Int, Int}()
     fcs    = NTuple{3,Int}[]
diff --git a/src/surface_nets.jl b/src/surface_nets.jl
deleted file mode 100644
index 5f618d0..0000000
--- a/src/surface_nets.jl
+++ /dev/null
@@ -1,270 +0,0 @@
-#
-# SurfaceNets in Julia
-#
-# Ported from the Javascript implementation by Mikola Lysenko (C) 2012
-# https://github.com/mikolalysenko/mikolalysenko.github.com/blob/master/Isosurface/js/surfacenets.js
-# MIT License
-#
-# Based on: S.F. Gibson, "Constrained Elastic Surface Nets". (1998) MERL Tech Report.
-#
-
-#Look up Table
-include("lut/sn.jl")
-
-function isosurface(sdf::AbstractArray{T, 3}, method::NaiveSurfaceNets, ::Type{VertType}=SVector{3,Float64}, ::Type{FaceType}=SVector{4, Int};
-                    origin=VertType(-1,-1,-1), widths=VertType(2,2,2)) where {T, VertType, FaceType}
-
-    scale = widths ./ VertType(size(sdf) .- 1)  # subtract 1 because an SDF with N points per side has N-1 cells
-
-    dims = size(sdf)
-
-    vertices = VertType[]
-    faces = FaceType[]
-
-    n = 0
-    R = Array{Int}([1, (dims[1]+1), (dims[1]+1)*(dims[2]+1)])
-    buf_no = 1
-
-    buffer = fill(zero(Int),R[3]*2)
-
-    #March over the voxel grid
-    zi = 0
-    @inbounds while zi<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))
-
-        for yi = 0:dims[2]-2
-            for xi = 0:dims[1]-2
-
-                inds = (xi,yi,zi)
-
-                # Read in 8 field values around this vertex and store them in an array
-                # Also calculate 8-bit mask, like in marching cubes, so we can speed up sign checks later
-                @inbounds grid = (sdf[xi+1,yi+1,zi+1],
-                                  sdf[xi+2,yi+1,zi+1],
-                                  sdf[xi+1,yi+2,zi+1],
-                                  sdf[xi+2,yi+2,zi+1],
-                                  sdf[xi+1,yi+1,zi+2],
-                                  sdf[xi+2,yi+1,zi+2],
-                                  sdf[xi+1,yi+2,zi+2],
-                                  sdf[xi+2,yi+2,zi+2])
-
-                mask = _get_cubeindex(grid, method.iso)
-
-                # Check for early termination if cell does not intersect boundary
-                if _no_triangles(mask)
-                    m += 1
-                    continue
-                end
-
-                # iso level correction
-                grid = grid .- method.iso
-
-                #Sum up edge intersections
-                edge_mask = sn_edge_table[mask]
-
-                _sn_add_verts!(inds, vertices, grid, edge_mask, buffer, m, scale, origin, method.eps, VertType)
-
-                #Now we need to add faces together, to do this we just loop over 3 basis components
-                _sn_add_faces!(inds, faces, edge_mask, mask, buffer, m, R, FaceType)
-
-                m += 1
-            end
-            m += 2
-        end
-        zi += 1
-        buf_no = xor(buf_no,1)
-        R[3]=-R[3]
-    end
-
-    vertices, faces
-end
-
-function isosurface(f::Function, method::NaiveSurfaceNets,
-                    ::Type{VertType}=SVector{3,Float64}, ::Type{FaceType}=SVector{4, Int};
-                    origin=VertType(-1,-1,-1), widths=VertType(2,2,2),
-                    samples::NTuple{3,T}=_DEFAULT_SAMPLES) where {T <: Integer, VertType, FaceType}
-
-    scale = widths ./ VertType(samples .- 1)  # subtract 1 because an SDF with N points per side has N-1 cells
-
-    vertices = VertType[]
-    faces = FaceType[]
-
-
-    R = Array{Int}([1, (samples[1]+1), (samples[1]+1)*(samples[2]+1)])
-    buf_no = 1
-
-    buffer = fill(zero(Int),R[3]*2)
-
-    zv = zero(eltype(VertType))
-    zvt = zero(VertType)
-    grid = (zv,zv,zv,zv,zv,zv,zv,zv)
-    points = (zvt,zvt,zvt,zvt,zvt,zvt,zvt,zvt)
-
-    #March over the voxel grid
-    zi = 0
-    @inbounds while zi<samples[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 + (samples[1]+1) * (1 + buf_no * (samples[2]+1))
-
-        for yi = 0:samples[2]-2
-            for xi = 0:samples[1]-2
-
-                inds = (xi,yi,zi)
-
-                if xi == 0
-                    points = (VertType(xi,yi,zi).* scale .+ origin,
-                              VertType(xi+1,yi,zi).* scale .+ origin,
-                              VertType(xi,yi+1,zi).* scale .+ origin,
-                              VertType(xi+1,yi+1,zi).* scale .+ origin,
-                              VertType(xi,yi,zi+1).* scale .+ origin,
-                              VertType(xi+1,yi,zi+1).* scale .+ origin,
-                              VertType(xi,yi+1,zi+1).* scale .+ origin,
-                              VertType(xi+1,yi+1,zi+1).* scale .+ origin)
-                    grid = (f(points[1]),
-                            f(points[2]),
-                            f(points[3]),
-                            f(points[4]),
-                            f(points[5]),
-                            f(points[6]),
-                            f(points[7]),
-                            f(points[8]))
-                else
-                    points = (points[2],
-                              VertType(xi+1,yi,zi).* scale .+ origin,
-                              points[4],
-                              VertType(xi+1,yi+1,zi).* scale .+ origin,
-                              points[6],
-                              VertType(xi+1,yi,zi+1).* scale .+ origin,
-                              points[8],
-                              VertType(xi+1,yi+1,zi+1).* scale .+ origin)
-                    grid = (grid[2],
-                            f(points[2]),
-                            grid[4],
-                            f(points[4]),
-                            grid[6],
-                            f(points[6]),
-                            grid[8],
-                            f(points[8]))
-                end
-
-                # Also calculate 8-bit mask, like in marching cubes, so we can speed up sign checks later
-                mask = _get_cubeindex(grid, method.iso)
-
-                # Check for early termination if cell does not intersect boundary
-                if _no_triangles(mask)
-                    m += 1
-                    continue
-                end
-
-                # iso level correction
-                grid = grid .- method.iso
-
-                #Sum up edge intersections
-                edge_mask = sn_edge_table[mask]
-
-                # add vertices
-                _sn_add_verts!(inds, vertices, grid, edge_mask, buffer, m, scale, origin, method.eps, VertType)
-
-                #Now we need to add faces together, to do this we just loop over 3 basis components
-                _sn_add_faces!(inds, faces, edge_mask, mask, buffer, m, R, FaceType)
-
-                m += 1
-            end
-            m += 2
-        end
-        zi += 1
-        buf_no = xor(buf_no,1)
-        R[3]=-R[3]
-    end
-
-    vertices, faces
-end
-
-@inline function _sn_add_verts!(inds, vertices, grid, edge_mask, buffer, m, scale, origin, eps, ::Type{VertType}) where {VertType}
-    v = zero(VertType)
-    T = eltype(VertType)
-    e_count = 0
-
-    #For every edge of the cube...
-    @inbounds for i=0x00:0x0b
-
-        #Use edge mask to check if it is crossed
-        iszero(edge_mask & (0x0001<<i)) && continue
-
-        #If it did, increment number of edge crossings
-        e_count += 1
-
-        #Now find the point of intersection
-        e0 = cube_edges[(i<<0x01)+1]       #Unpack vertices
-        e1 = cube_edges[(i<<0x01)+2]
-        g0 = grid[e0+1]                 #Unpack grid values
-        g1 = grid[e1+1]
-        t  = g0 - g1                 #Compute point of intersection
-
-        abs(t) <= eps && continue
-
-        t = g0 / t
-
-        #Interpolate vertices and add up intersections (this can be done without multiplying)
-        a1, a2, a3 = 0x01 & e0, 0x02 & e0, 0x04 & e0
-        b1, b2, b3 = 0x01 & e1, 0x02 & e1, 0x04 & e1
-        xj = T(!iszero(a1))
-        yj = T(!iszero(a2))
-        zj = T(!iszero(a3))
-        a1 != b1 && (xj += !iszero(a1) ? -t : t)
-        a2 != b2 && (yj += !iszero(a2) ? -t : t)
-        a3 != b3 && (zj += !iszero(a3) ? -t : t)
-        v = v .+ VertType(xj,yj,zj)
-
-    end # edge check
-
-    #Now we just average the edge intersections and add them to coordinate
-    s = one(T) / e_count
-    v = (VertType(inds...)  .+ s .* v) .* scale .+ origin
-
-    #Add vertex to buffer, store pointer to vertex index in buffer
-    buffer[m+1] = length(vertices)
-    push!(vertices, v)
-end
-
-function _sn_add_faces!(inds, faces, edge_mask, mask, buffer, m, R, ::Type{FaceType}) where {FaceType}
-    for i = 0x00:0x02
-        #The first three entries of the edge_mask count the crossings along the edge
-        iszero(edge_mask & (0x0001<<i)) && continue
-
-        # i = axes we are point along.  iu, iv = orthogonal axes
-        iu = (i+0x01)%0x03
-        iv = (i+0x02)%0x03
-
-        #If we are on a boundary, skip it
-        iszero(inds[iu+1]) || iszero(inds[iv+1]) && continue
-
-        #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 !iszero(mask & 0x01)
-            if length(FaceType) == 4
-                push!(faces,FaceType(buffer[m+1]+1, buffer[m-du+1]+1, buffer[m-du-dv+1]+1, buffer[m-dv+1]+1))
-            elseif length(FaceType) == 3
-                push!(faces,FaceType(buffer[m+1]+1, buffer[m-du+1]+1, buffer[m-du-dv+1]+1))
-                push!(faces,FaceType(buffer[m-du-dv+1]+1, buffer[m-dv+1]+1, buffer[m+1]+1))
-            end
-        else
-            if length(FaceType) == 4
-                push!(faces,FaceType(buffer[m+1]+1, buffer[m-dv+1]+1, buffer[m-du-dv+1]+1, buffer[m-du+1]+1))
-            elseif length(FaceType) == 3
-                push!(faces,FaceType(buffer[m+1]+1, buffer[m-dv+1]+1, buffer[m-du-dv+1]+1))
-                push!(faces,FaceType(buffer[m-du-dv+1]+1, buffer[m-du+1]+1, buffer[m+1]+1))
-            end
-        end
-    end
-end
diff --git a/test/runtests.jl b/test/runtests.jl
index c02bb76..6073c62 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -5,7 +5,7 @@ using Statistics: mean
 using LinearAlgebra: dot, norm
 using Random
 
-const algos = (MarchingCubes, MarchingTetrahedra, NaiveSurfaceNets)
+const algos = (MarchingCubes, MarchingTetrahedra)
 
 sphere_function(v) = sqrt(sum(dot(v, v))) - 1
 torus_function(v) = (sqrt(v[1]^2 + v[2]^2) - 0.5)^2 + v[3]^2 - 0.25
@@ -54,37 +54,8 @@ end
         end
     end
 
-    #TODO: SurfaceNets functional variant is bugged
-    # Naive Surface Nets has no accuracy guarantee, and is a weighted sum
-    # so a larger tolerance is needed for this one. In addition,
-    # quad -> triangle conversion is not functioning correctly
-    # see: https://github.com/JuliaGeometry/GeometryTypes.jl/issues/169
-    #points, faces = isosurface(f, NaiveSurfaceNets(0.5), samples=(21, 21, 21))
-
-    # should be centered on the origin
-    #@test mean(points) ≈ [0, 0, 0] atol=0.15*resolution
-    # and should be symmetric about the origin
-    #@test maximum(points) ≈ [0.5, 0.5, 0.5] atol=0.2
-    #@test minimum(points) ≈ [-0.5, -0.5, -0.5] atol=0.2
 end
-#@testset "surface nets" begin
-#
-#
-#    """
-#    Test the isosurface function with NaiveSurfaceNets and a function.
-#
-#    This code tests the `isosurface` function with the `NaiveSurfaceNets` algorithm and two different level set functions: `sphere_function` and `torus_function`. It checks that the number of vertices and faces in the resulting meshes match expected values.
-#    """
-#    # Test the isosurface function with NaiveSurfaceNets and a function
-#    points_sphere, faces_sphere = isosurface(sphere_function, NaiveSurfaceNets(), origin=SVector(-1, -1, -1), widths=SVector(2, 2, 2), samples=(21, 21, 21))
-#    points_torus, faces_torus = isosurface(torus_function, NaiveSurfaceNets(), origin=SVector(-2, -2, -2), widths=SVector(4, 4, 4), samples=(81, 81, 81))
-#
-#    # Add assertions to check the output
-#    @test length(points_sphere) == 1832
-#    @test length(faces_sphere) == 1830
-#    @test length(points_torus) == 5532
-#    @test length(faces_torus) == 5532
-#end
+
 
 @testset "noisy spheres" begin
     # Produce a level set function that is a noisy version of the distance from
@@ -144,7 +115,7 @@ end
             points_fn, faces_fn = isosurface(torus_function, algo(), -2:2, -2:2, -2:2, samples=(81, 81, 81))
 
             # Extract isosurface using an array
-            torus_array = [torus_function(SVector(x, y, z)) for x in -2:0.05:2, y in -2:0.05:2, z in -2:0.05:2]
+            torus_array = [torus_function((x, y, z)) for x in -2:0.05:2, y in -2:0.05:2, z in -2:0.05:2]
             points_arr, faces_arr = isosurface(torus_array, algo(), -2:2, -2:2, -2:2)
 
             # Test that the vertices and faces are the same for both cases
@@ -159,17 +130,20 @@ end
     samples = randn(10, 10, 10)
 
     # Extract isosurfaces using MarchingTetrahedra
-    points_mt1, faces_mt1 = isosurface(samples, MarchingTetrahedra(), origin=SVector(0, 0, 0), widths=SVector(1, 1, 1))
-    points_mt2, faces_mt2 = isosurface(Float32.(samples), MarchingTetrahedra(), origin=SVector(0, 0, 0), widths=SVector(1, 1, 1))
-
-    @test all(isapprox.(points_mt1, points_mt2, rtol=1e-6))
-    @test all(faces_mt1 == faces_mt2)
+    points_mt1, faces_mt1 = isosurface(samples, MarchingTetrahedra(), 0:1, 0:1, 0:1)
+    points_mt2, faces_mt2 = isosurface(Float32.(samples), MarchingTetrahedra(), 0:1, 0:1, 0:1)
+
+    @test length(points_mt1) == length(points_mt2)
+    #for i in eachindex(points_mt1)
+    #    @test all(points_mt1[i] .≈ points_mt2[i])
+    #end
+    @test length(faces_mt1) == length(faces_mt2)
+    for i in eachindex(faces_mt1)
+        @test all(faces_mt1[i] .≈ faces_mt2[i])
+    end
 
-    #=
     @testset "forward diff" begin
         # Demonstrate that we can even take derivatives through the meshing algorithm
-        origin = SVector(-1, -1, -1)
-        widths = SVector(2, 2, 2)
         resolution = 0.1
 
         for algo in (MarchingCubes, MarchingTetrahedra)
@@ -183,7 +157,7 @@ end
                         # of some quantity you might want to differentiate in a mesh,
                         # and has the benefit for testing of having a trivial expected
                         # solution.
-                        points, _ = isosurface(norm, algo(isovalue,  reduce_vert=false), samples=(21, 21, 21))
+                        points, _ = isosurface(norm, algo(iso=isovalue), samples=(21, 21, 21))
                         mean(norm.(points))
                     end
 
@@ -195,7 +169,7 @@ end
         end
     end
 
-    =#
+
 
     @testset "type stability" begin
         # verify that we don't lose type stability just by mixing arguments
@@ -203,36 +177,38 @@ end
         data = randn(5, 5, 5)
         iso = 0.2
         eps = 1e-4
-        for algo_ty in (MarchingTetrahedra, MarchingCubes,NaiveSurfaceNets)
-            algo1 = algo_ty(iso, eps)
-            algo2 = algo_ty(Float64(iso), Float16(eps))
-            algo3 = algo_ty(Float32(iso), Float64(eps))
-            @inferred(isosurface(data, algo1, SVector{3, Float16}, SVector{3,UInt32}))
-            @inferred(isosurface(Float32.(data), algo2, SVector{3, Float32}, SVector{3,Int16}))
-            @inferred(isosurface(Float64.(data), algo3, SVector{3, Float64}, SVector{3,UInt}))
-        end
+        algo1 = MarchingTetrahedra(iso=iso, eps=eps)
+        algo2 = MarchingTetrahedra(iso=Float16(iso), eps=Float16(eps))
+        algo3 = MarchingTetrahedra(iso=Float32(iso), eps=Float32(eps))
+        @inferred(isosurface(data, algo1))
+        @inferred(isosurface(Float32.(data), algo2))
+        @inferred(isosurface(Float64.(data), algo3))
+
+        algo1 = MarchingCubes(iso=iso)
+        algo2 = MarchingCubes(iso=Float16(iso))
+        algo3 = MarchingCubes(iso=Float32(iso))
+        @inferred(isosurface(data, algo1))
+        @inferred(isosurface(Float32.(data), algo2))
+        @inferred(isosurface(Float64.(data), algo3))
     end
 
     @testset "Float16" begin
         sdf_torus = [((sqrt(x^2 + y^2) - 0.5)^2 + z^2 - 0.25) for x in -2:0.05:2, y in -2:0.05:2, z in -2:0.05:2]
-        points_nsn, faces_nsn = isosurface(Float16.(sdf_torus), NaiveSurfaceNets(Float16(0.0)), origin=SVector(Float16(-2), Float16(-2), Float16(-2)), widths=SVector(Float16(4), Float16(4), Float16(4)))
-        points_mt, faces_mt = isosurface(Float16.(sdf_torus), MarchingTetrahedra(Float16(0.0)), origin=SVector(Float16(-2), Float16(-2), Float16(-2)), widths=SVector(Float16(4), Float16(4), Float16(4)))
-        points_mc, faces_mc = isosurface(Float16.(sdf_torus), MarchingCubes(Float16(0.0)), origin=SVector(Float16(-2), Float16(-2), Float16(-2)), widths=SVector(Float16(4), Float16(4), Float16(4)))
-
-        @test_broken typeof(points_nsn) == Vector{SVector{3, Float16}}
-        @test typeof(faces_nsn) == Vector{SVector{4, Int}}
-        @test_broken typeof(points_mt) == Vector{SVector{3, Float16}}
-        @test typeof(faces_mt) == Vector{SVector{3, Int}}
-        @test_broken typeof(points_mc) == Vector{SVector{3, Float16}}
-        @test typeof(faces_mc) == Vector{SVector{3, Int}}
+        points_mt, faces_mt = isosurface(Float16.(sdf_torus), MarchingTetrahedra(iso=Float16(0.0), eps=Float16(1e-3)), Float16(-2):Float16(2), Float16(-2):Float16(2), Float16(-2):Float16(2))
+        points_mc, faces_mc = isosurface(Float16.(sdf_torus), MarchingCubes(iso=Float16(0.0)), Float16(-2):Float16(2), Float16(-2):Float16(2), Float16(-2):Float16(2))
+
+        @test typeof(points_mt) == Vector{NTuple{3, Float16}}
+        @test typeof(faces_mt) == Vector{NTuple{3, Int}}
+        @test typeof(points_mc) == Vector{NTuple{3, Float16}}
+        @test typeof(faces_mc) == Vector{NTuple{3, Int}}
     end
     @testset "Integers" begin
         A = rand(Int, 10,10,10)
         for algo in algos
             @testset "$algo" begin
                 p, t = isosurface(A, algo())
-                @test typeof(p) <: Vector{SVector{3,Float64}}
-                @test typeof(t) <: Vector{ SVector{ Meshing.default_face_length(algo()), Int } }
+                @test typeof(p) <: Vector{NTuple{3, Float64}}
+                @test typeof(t) <: Vector{NTuple{3, Int}}
             end
         end
     end

From 6929b513d3fb2cceffaa855b8bcea4db37d1b07b Mon Sep 17 00:00:00 2001
From: Steve Kelly <kd2cca@gmail.com>
Date: Tue, 2 Jul 2024 07:20:04 -0400
Subject: [PATCH 8/8] rm cruft

---
 src/algorithmtypes.jl | 12 +-----------
 src/lut/sn.jl         | 35 -----------------------------------
 2 files changed, 1 insertion(+), 46 deletions(-)
 delete mode 100644 src/lut/sn.jl

diff --git a/src/algorithmtypes.jl b/src/algorithmtypes.jl
index 57281d4..02d0bf3 100644
--- a/src/algorithmtypes.jl
+++ b/src/algorithmtypes.jl
@@ -14,7 +14,7 @@ abstract type AbstractAdaptiveMeshingAlgorithm end
 
 
 """
-    MarchingCubes(iso=0.0)
+    MarchingCubes(;iso=0.0)
 
 Specifies the use of the Marching Cubes algorithm for isosurface extraction.
 This algorithm provides a good balance between performance and vertex count.
@@ -28,10 +28,7 @@ Base.@kwdef struct MarchingCubes{T} <: AbstractMeshingAlgorithm
 end
 
 """
-    MarchingTetrahedra(iso=0.0, eps=1e-3)
     MarchingTetrahedra(;iso=0.0, eps=1e-3)
-    MarchingTetrahedra(iso)
-    MarchingTetrahedra(iso,eps)
 
 Specifies the use of the Marching Tetrahedra algorithm for isosurface extraction.
 This algorithm has a roughly 2x performance penalty compared to Marching Cubes,
@@ -75,10 +72,3 @@ struct AdaptiveMarchingTetrahedra{T} <: AbstractAdaptiveMeshingAlgorithm
     rtol::T
     atol::T
 end
-
-
-#
-# Helper functions
-#
-
-default_face_length(::Union{MarchingCubes,MarchingTetrahedra}) = 3
diff --git a/src/lut/sn.jl b/src/lut/sn.jl
deleted file mode 100644
index 610c62f..0000000
--- a/src/lut/sn.jl
+++ /dev/null
@@ -1,35 +0,0 @@
-const cube_edges = (0x00, 0x01, 0x00, 0x02, 0x00, 0x04, 0x01, 0x03, 0x01, 0x05, 0x02, 0x03,
-                    0x02, 0x06, 0x03, 0x07, 0x04, 0x05, 0x04, 0x06, 0x05, 0x07, 0x06, 0x07)
-
-const sn_edge_table = (0x0007, 0x0019, 0x001e, 0x0062, 0x0065, 0x007b, 0x007c, 0x00a8,
-                       0x00af, 0x00b1, 0x00b6, 0x00ca, 0x00cd, 0x00d3, 0x00d4, 0x0304,
-                       0x0303, 0x031d, 0x031a, 0x0366, 0x0361, 0x037f, 0x0378, 0x03ac,
-                       0x03ab, 0x03b5, 0x03b2, 0x03ce, 0x03c9, 0x03d7, 0x03d0, 0x0510,
-                       0x0517, 0x0509, 0x050e, 0x0572, 0x0575, 0x056b, 0x056c, 0x05b8,
-                       0x05bf, 0x05a1, 0x05a6, 0x05da, 0x05dd, 0x05c3, 0x05c4, 0x0614,
-                       0x0613, 0x060d, 0x060a, 0x0676, 0x0671, 0x066f, 0x0668, 0x06bc,
-                       0x06bb, 0x06a5, 0x06a2, 0x06de, 0x06d9, 0x06c7, 0x06c0, 0x0a40,
-                       0x0a47, 0x0a59, 0x0a5e, 0x0a22, 0x0a25, 0x0a3b, 0x0a3c, 0x0ae8,
-                       0x0aef, 0x0af1, 0x0af6, 0x0a8a, 0x0a8d, 0x0a93, 0x0a94, 0x0944,
-                       0x0943, 0x095d, 0x095a, 0x0926, 0x0921, 0x093f, 0x0938, 0x09ec,
-                       0x09eb, 0x09f5, 0x09f2, 0x098e, 0x0989, 0x0997, 0x0990, 0x0f50,
-                       0x0f57, 0x0f49, 0x0f4e, 0x0f32, 0x0f35, 0x0f2b, 0x0f2c, 0x0ff8,
-                       0x0fff, 0x0fe1, 0x0fe6, 0x0f9a, 0x0f9d, 0x0f83, 0x0f84, 0x0c54,
-                       0x0c53, 0x0c4d, 0x0c4a, 0x0c36, 0x0c31, 0x0c2f, 0x0c28, 0x0cfc,
-                       0x0cfb, 0x0ce5, 0x0ce2, 0x0c9e, 0x0c99, 0x0c87, 0x0c80, 0x0c80,
-                       0x0c87, 0x0c99, 0x0c9e, 0x0ce2, 0x0ce5, 0x0cfb, 0x0cfc, 0x0c28,
-                       0x0c2f, 0x0c31, 0x0c36, 0x0c4a, 0x0c4d, 0x0c53, 0x0c54, 0x0f84,
-                       0x0f83, 0x0f9d, 0x0f9a, 0x0fe6, 0x0fe1, 0x0fff, 0x0ff8, 0x0f2c,
-                       0x0f2b, 0x0f35, 0x0f32, 0x0f4e, 0x0f49, 0x0f57, 0x0f50, 0x0990,
-                       0x0997, 0x0989, 0x098e, 0x09f2, 0x09f5, 0x09eb, 0x09ec, 0x0938,
-                       0x093f, 0x0921, 0x0926, 0x095a, 0x095d, 0x0943, 0x0944, 0x0a94,
-                       0x0a93, 0x0a8d, 0x0a8a, 0x0af6, 0x0af1, 0x0aef, 0x0ae8, 0x0a3c,
-                       0x0a3b, 0x0a25, 0x0a22, 0x0a5e, 0x0a59, 0x0a47, 0x0a40, 0x06c0,
-                       0x06c7, 0x06d9, 0x06de, 0x06a2, 0x06a5, 0x06bb, 0x06bc, 0x0668,
-                       0x066f, 0x0671, 0x0676, 0x060a, 0x060d, 0x0613, 0x0614, 0x05c4,
-                       0x05c3, 0x05dd, 0x05da, 0x05a6, 0x05a1, 0x05bf, 0x05b8, 0x056c,
-                       0x056b, 0x0575, 0x0572, 0x050e, 0x0509, 0x0517, 0x0510, 0x03d0,
-                       0x03d7, 0x03c9, 0x03ce, 0x03b2, 0x03b5, 0x03ab, 0x03ac, 0x0378,
-                       0x037f, 0x0361, 0x0366, 0x031a, 0x031d, 0x0303, 0x0304, 0x00d4,
-                       0x00d3, 0x00cd, 0x00ca, 0x00b6, 0x00b1, 0x00af, 0x00a8, 0x007c,
-                       0x007b, 0x0065, 0x0062, 0x001e, 0x0019, 0x0007)