Internals update

Project.toml

@@ -2,6 +2,7 @@ name = "Meshing"
 uuid = "e6723b4c-ebff-59f1-b4b7-d97aa5274f73"
 
 [deps]
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 GeometryTypes = "4d00f742-c7ba-57c2-abde-4428a4b178cb"
 
 [extras]

docs/src/api.md

@@ -2,10 +2,11 @@
 
 ## Quick Start
 
-This package inherits the [mesh types](http://juliageometry.github.io/GeometryTypes.jl/latest/types.html#Meshes-1)
-from [GeometryTypes.jl](https://github.com/JuliaGeometry/GeometryTypes.jl).
+The easiest way to work with Meshing is with GeometryTypes.
+This package extends the [mesh constructors](http://juliageometry.github.io/GeometryTypes.jl/latest/types.html#Meshes-1)
+from [GeometryTypes.jl](https://github.com/JuliaGeometry/GeometryTypes.jl) for convience.
 
-The algorithms operate on a `Function` or a `SignedDistanceField` and output a concrete `AbstractMesh`. For example:
+The algorithms operate on a `Function`, `AbstractArray`, or `SignedDistanceField` and output a concrete `AbstractMesh`. For example:
 
 ```
 using Meshing
@@ -14,7 +15,7 @@ using LinearAlgebra: dot, norm
 using FileIO
 
 # Mesh an equation of sphere in the Axis-Aligned Bounding box starting
-# at -1,-1,-1 and widths of 2,2,2
+# at -1,-1,-1 and widths of 2,2,2 using Marching Cubes
 m = GLNormalMesh(HyperRectangle(Vec(-1,-1,-1.), Vec(2,2,2.)), MarchingCubes()) do v
     sqrt(sum(dot(v,v))) - 1
 end
@@ -23,9 +24,6 @@ end
 save("sphere.ply",m)
 ```
 
-The general API is: ```(::Type{MT})(sdf::Function, method::AbstractMeshingAlgorithm) where {MT <: AbstractMesh}``` or ```(::Type{MT})(sdf::SignedDistanceField, method::AbstractMeshingAlgorithm) where {MT <: AbstractMesh}```
-
-
 For a full listing of concrete `AbstractMesh` types see [GeometryTypes.jl mesh documentation](http://juliageometry.github.io/GeometryTypes.jl/latest/types.html#Meshes-1).
 
 ## Meshing Algorithms
@@ -39,16 +37,31 @@ Each takes an optional `iso` and `eps` parameter, e.g. `MarchingCubes(0.0,1e-6)`
 
 Here `iso` controls the offset for the boundary detection. By default this is set to 0. `eps` is the detection tolerance for a voxel edge intersection.
 
+Users must construct an algorithm type and use it as an argument to a GeometryTypes mesh call or `isosurface` call.
+
 Below is a comparison of the algorithms:
 
-| Algorithm          | Accurate | Manifold | Performance Penalty | Face Type |
-|--------------------|----------|----------|---------------------|-----------|
-| MarchingCubes      | Yes      | No       | ~2x                 | Triangle  |
-| MarchingTetrahedra | Yes      | Yes      | ~3x                 | Triangle  |
-| NaiveSurfaceNets   | No       | Yes      | 1x                  | Quad      |
+|                     | Face Type | Unique Vertices | Performance | Interpolation     |
+|---------------------|-----------|-----------------|-------------|-------------------|
+| Naive Surface Nets  | Quad      | Yes             | ~1x         | Voxel Edge Weight |
+| Marching Cubes      | Triangle  | No/Partial      | 1x          | Linear on Edge    |
+| Marching Tetrahedra | Triangle  | Yes             | 3x          | Linear on Edge    |
 
 ```@docs
 MarchingCubes
 MarchingTetrahedra
 NaiveSurfaceNets
+Meshing.AbstractMeshingAlgorithm
+```
+
+## Isosurface
+
+`isosurface` is the common and generic API for isosurface extraction with any type of abstract vector/vertex/face type.
+
+```@docs
+isosurface
 ```
+
+## GeometryTypes
+
+Meshing extends the mesh types in GeometryTypes for convience and use with visualization tools such as Makie and MeshCat.

docs/src/index.md

@@ -1,8 +1,19 @@
 # Meshing.jl
 
-This package provides a comprehensive suite of meshing (isosurface extraction) algorithms.
+This package provides a suite of meshing (isosurface extraction) algorithms.
 
 Algorithms included:
 * [Marching Tetrahedra](https://en.wikipedia.org/wiki/Marching_tetrahedra)
 * [Marching Cubes](https://en.wikipedia.org/wiki/Marching_cubes)
 * [Naive Surface Nets](https://0fps.net/2012/07/12/smooth-voxel-terrain-part-2/)
+
+## What is isosurface extraction?
+
+Isosurface extraction is a common technique to visualize and analyze scalar fields and functions.
+It takes scalar data that is structured in grids and converts this into a series of points and faces suitable for 3D visualization, GPU computing, 3D Printing, or other analyses.
+
+There are several applications for these techniques such as:
+- Medical Imaging of CT and MRI data
+- Visualization of Functions
+- Solid Modeling
+- Terrain Generation

src/Meshing.jl

@@ -1,16 +1,18 @@
 module Meshing
 
-using GeometryTypes
+using GeometryTypes,
+      StaticArrays
 
-abstract type AbstractMeshingAlgorithm end
+const _DEFAULT_SAMPLES = (24,24,24)
 
 include("algorithmtypes.jl")
+include("geometrytypes_api.jl")
 include("common.jl")
 include("marching_tetrahedra.jl")
 include("marching_cubes.jl")
 include("surface_nets.jl")
 
-export marching_cubes,
+export isosurface,
        MarchingCubes,
        MarchingTetrahedra,
        NaiveSurfaceNets

src/algorithmtypes.jl

@@ -1,4 +1,12 @@
 
+"""
+    AbstractMeshingAlgorithm
+
+Abstract type to specify an algorithm for isosurface extraction.
+    (unexported)
+"""
+abstract type AbstractMeshingAlgorithm end
+
 """
     MarchingCubes(iso=0.0, eps=1e-3, reduceverts=true, insidepositive=false)
     MarchingCubes(;iso=0.0, eps=1e-3, reduceverts=true, insidepositive=false)
@@ -80,3 +88,11 @@ end
 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)
+
+
+#
+# Helper functions
+#
+
+default_face_length(::Union{MarchingCubes,MarchingTetrahedra}) = 3
+default_face_length(::NaiveSurfaceNets) = 4

src/common.jl

@@ -64,4 +64,21 @@ function _determine_types(meshtype, fieldtype=Float64, facelen=3)
         FaceType = Face{facelen, Int}
     end
     VertType, FaceType
-end
+end
+
+#
+# General isosurface docstring
+# TODO FORMATTING?
+
+@doc """
+
+    `function isosurface(sdf::AbstractArray{T, 3}, method::AbstractMeshingAlgorithm,
+                         [ VertType = SVector{3,Float64} ], [ FaceType} = SVector{3, Int} ] ;
+                         origin = SVector(-1.0,-1.0,-1.0), widths = SVector(2.0,2.0,2.0))`
+
+    `function isosurface(f::Function, method::AbstractMeshingAlgorithm,
+                         [ VertType = SVector{3,Float64} ], [FaceType = SVector{3, Int} ] ;
+                         origin = SVector(-1.0,-1.0,-1.0), widths = SVector(2.0,2.0,2.0),
+                         samples=(50,50,50))`
+
+""" isosurface

src/geometrytypes_api.jl

@@ -0,0 +1,28 @@
+function (::Type{MT})(df::SignedDistanceField{3,ST,FT}, method::AbstractMeshingAlgorithm)::MT where {MT <: AbstractMesh, ST, FT}
+    vertex_eltype = promote_type(FT, typeof(method.iso), typeof(method.eps))
+    VertType, FaceType = _determine_types(MT, vertex_eltype, default_face_length(method))
+    h = df.bounds
+    vts, fcs = isosurface(df.data, method, VertType, FaceType, origin=VertType(origin(h)), widths=VertType(widths(h)))
+    MT(vts, fcs)::MT
+end
+
+function (::Type{MT})(f::Function, h::HyperRectangle, samples::NTuple{3,T}, method::AbstractMeshingAlgorithm)::MT where {MT <: AbstractMesh, T <: Integer}
+    vertex_eltype = promote_type(T, typeof(method.iso), typeof(method.eps))
+    VertType, FaceType = _determine_types(MT,vertex_eltype, default_face_length(method))
+    vts, fcs = isosurface(f, method, VertType, FaceType, samples=samples, origin=VertType(origin(h)), widths=VertType(widths(h)))
+    MT(vts, fcs)::MT
+end
+
+function (::Type{MT})(f::Function, h::HyperRectangle, method::AbstractMeshingAlgorithm; samples::NTuple{3,T}=_DEFAULT_SAMPLES)::MT where {MT <: AbstractMesh, T <: Integer}
+    vertex_eltype = promote_type(T, typeof(method.iso), typeof(method.eps))
+    VertType, FaceType = _determine_types(MT,vertex_eltype, default_face_length(method))
+    vts, fcs = isosurface(f, method, VertType, FaceType, samples=samples, origin=VertType(origin(h)), widths=VertType(widths(h)))
+    MT(vts, fcs)::MT
+end
+
+function (::Type{MT})(volume::AbstractArray{T, 3}, method::AbstractMeshingAlgorithm; vargs...) where {MT <: AbstractMesh, T}
+    vertex_eltype = promote_type(T, typeof(method.iso), typeof(method.eps))
+    VertType, FaceType = _determine_types(MT,vertex_eltype, default_face_length(method))
+    vts, fcs = isosurface(volume, method, VertType, FaceType, vargs...)
+    MT(vts, fcs)::MT
+end

src/lut/mt.jl

@@ -20,14 +20,21 @@ Voxel corner and edge indexing conventions
   /
  X
 """
-const voxCrnrPos = (Point(0, 0, 0),
-                Point(0, 1, 0),
-                Point(1, 1, 0),
-                Point(1, 0, 0),
-                Point(0, 0, 1),
-                Point(0, 1, 1),
-                Point(1, 1, 1),
-                Point(1, 0, 1))
+# 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))
+end
+
+const voxCrnrPosInt = voxCrnrPos(SVector{3,Int})
+
 # the voxel IDs at either end of the tetrahedra edges, by edge ID
 const voxEdgeCrnrs = ((1, 2),
                 (2, 3),