Comments (1)
sjkelly
commented on June 12, 2024
steve:~/.julia/v0.4/Meshing(master)$ diff ./src/MarchingTetrahedra.jl ../Meshes/src/isosurface.jl
1,134c1,133
< """
< Marching Tetrahedra is an algorithm for extracting a triangular
< mesh representation of an isosurface of a scalar volumetric
< function sampled on a rectangular grid.
<
< We divide the cube into six tetrahedra. [It is possible to divide
< a cube into five tetrahedra, but not in a way that a translated
< version of the division would share face diagonals. (It reqires a
< reflection.)]
< """
< module MarchingTetrahedra
<
< using GeometryTypes
<
< """
< Voxel corner and edge indexing conventions
<
< Z
< |
<
< 5------5------6 Extra edges not drawn
< /| /| -----------
< 8 | 6 | - face diagonals
< / 9 / 10 - 13: 1 to 3
< 8------7------7 | - 14: 1 to 8
< | | | | - 15: 1 to 6
< | 1------1--|---2 -- Y - 16: 5 to 7
< 12 / 11 / - 17: 2 to 7
< | 4 | 2 - 18: 4 to 7
< |/ |/ - body diagonal
< 4------3------3 - 19: 1 to 7
<
< /
< X
< """
<
< immutable VoxelIndices{T <: Integer}
< voxCrnrPos::NTuple{8,Vec{3,T}}
< voxEdgeCrnrs::NTuple{19, Vec{2,T}}
< voxEdgeDir::NTuple{19,T}
< voxEdgeIx::Matrix{T}
< subTets::Matrix{T}
< tetEdgeCrnrs::Matrix{T}
< tetTri::Matrix{T}
<
< function VoxelIndices()
< VT3 = Vec{3,T}
< VT2 = Vec{2,T}
< voxCrnrPos = (VT3(0, 0, 0),
< VT3(0, 1, 0),
< VT3(1, 1, 0),
< VT3(1, 0, 0),
< VT3(0, 0, 1),
< VT3(0, 1, 1),
< VT3(1, 1, 1),
< VT3(1, 0, 1))
< # the voxel IDs at either end of the tetrahedra edges, by edge ID
< voxEdgeCrnrs = (VT2(1, 2),
< VT2(2, 3),
< VT2(4, 3),
< VT2(1, 4),
< VT2(5, 6),
< VT2(6, 7),
< VT2(8, 7),
< VT2(5, 8),
< VT2(1, 5),
< VT2(2, 6),
< VT2(3, 7),
< VT2(4, 8),
< VT2(1, 3),
< VT2(1, 8),
< VT2(1, 6),
< VT2(5, 7),
< VT2(2, 7),
< VT2(4, 7),
< VT2(1, 7))
<
< # direction codes:
< # 0 => +x, 1 => +y, 2 => +z,
< # 3 => +xy, 4 => +xz, 5 => +yz, 6 => +xyz
< voxEdgeDir = convert(NTuple{19,T}, (1,0,1,0,1,0,1,0,2,2,2,2,3,4,5,3,4,5,6))
<
< # For a pair of corner IDs, the edge ID joining them
< # 0 denotes a pair with no edge
< voxEdgeIx = [0 1 13 4 9 15 19 14;
< 1 0 2 0 0 10 17 0;
< 13 2 0 3 0 0 11 0;
< 4 0 3 0 0 0 18 12;
< 9 0 0 0 0 5 16 8;
< 15 10 0 0 5 0 6 0;
< 19 17 11 18 16 6 0 7;
< 14 0 0 12 8 0 7 0]
<
< # voxel corners that comprise each of the six tetrahedra
< subTets = [ 1 3 2 7;
< 1 8 4 7;
< 1 4 3 7;
< 1 2 6 7;
< 1 5 8 7;
< 1 6 5 7]'
< # tetrahedron corners for each edge (indices 1-4)
< tetEdgeCrnrs = [1 2;
< 2 3;
< 1 3;
< 1 4;
< 2 4;
< 3 4]'
<
< # triangle cases for a given tetrahedron edge code
< tetTri = T[ 0 0 0 0 0 0;
< 1 3 4 0 0 0;
< 1 5 2 0 0 0;
< 3 5 2 3 4 5;
< 2 6 3 0 0 0;
< 1 6 4 1 2 6;
< 1 5 6 1 6 3;
< 4 5 6 0 0 0;
< 4 6 5 0 0 0;
< 1 6 5 1 3 6;
< 1 4 6 1 6 2;
< 2 3 6 0 0 0;
< 3 2 5 3 5 4;
< 1 2 5 0 0 0;
< 1 4 3 0 0 0;
< 0 0 0 0 0 0]'
<
< new(voxCrnrPos,
< voxEdgeCrnrs,
< voxEdgeDir,
< voxEdgeIx,
< subTets,
< tetEdgeCrnrs,
< tetTri)
< end
---
> # *** Marching Tetrahedra ***
> #
> # Marching Tetrahedra is an algorithm for extracting a triangular
> # mesh representation of an isosurface of a scalar volumetric
> # function sampled on a rectangular grid.
> #
> # We divide the cube into six tetrahedra. [It is possible to divide
> # a cube into five tetrahedra, but not in a way that a translated
> # version of the division would share face diagonals. (It reqires a
> # reflection.)]
> #
> # Voxel corner and edge indexing conventions
> #
> # Z
> # |
> #
> # 5------5------6 Extra edges not drawn
> # /| /| -----------
> # 8 | 6 | - face diagonals
> # / 9 / 10 - 13: 1 to 3
> # 8------7------7 | - 14: 1 to 8
> # | | | | - 15: 1 to 6
> # | 1------1--|---2 -- Y - 16: 5 to 7
> # 12 / 11 / - 17: 2 to 7
> # | 4 | 2 - 18: 4 to 7
> # |/ |/ - body diagonal
> # 4------3------3 - 19: 1 to 7
> #
> # /
> # X
> immutable VoxelIndexes{T <: Integer}
> voxCrnrPos ::Vector{NTuple{3, T}}
> voxEdgeCrnrs ::Vector{NTuple{2, T}}
> voxEdgeDir ::Vector{T}
> voxEdgeIx ::Matrix{T}
> subTets ::Matrix{T}
> subTets ::Matrix{T}
> tetEdgeCrnrs ::Matrix{T}
> tetTri ::Matrix{T}
>
> function VoxelIndexes()
> voxCrnrPos = NTuple{3, T}[
> (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
> voxEdgeCrnrs = NTuple{2, T}[
> (1, 2),
> (2, 3),
> (4, 3),
> (1, 4),
> (5, 6),
> (6, 7),
> (8, 7),
> (5, 8),
> (1, 5),
> (2, 6),
> (3, 7),
> (4, 8),
> (1, 3),
> (1, 8),
> (1, 6),
> (5, 7),
> (2, 7),
> (4, 7),
> (1, 7)
> ]
>
> # direction codes:
> # 0 => +x, 1 => +y, 2 => +z,
> # 3 => +xy, 4 => +xz, 5 => +yz, 6 => +xyz
> voxEdgeDir = T[1,0,1,0,1,0,1,0,2,2,2,2,3,4,5,3,4,5,6]
>
> # For a pair of corner IDs, the edge ID joining them
> # 0 denotes a pair with no edge
> voxEdgeIx = [0 1 13 4 9 15 19 14;
> 1 0 2 0 0 10 17 0;
> 13 2 0 3 0 0 11 0;
> 4 0 3 0 0 0 18 12;
> 9 0 0 0 0 5 16 8;
> 15 10 0 0 5 0 6 0;
> 19 17 11 18 16 6 0 7;
> 14 0 0 12 8 0 7 0]
>
> # voxel corners that comprise each of the six tetrahedra
> subTets = [ 1 3 2 7;
> 1 8 4 7;
> 1 4 3 7;
> 1 2 6 7;
> 1 5 8 7;
> 1 6 5 7]'
> # tetrahedron corners for each edge (indices 1-4)
> tetEdgeCrnrs = [1 2;
> 2 3;
> 1 3;
> 1 4;
> 2 4;
> 3 4]'
>
> # triangle cases for a given tetrahedron edge code
> tetTri = T[ 0 0 0 0 0 0;
> 1 3 4 0 0 0;
> 1 5 2 0 0 0;
> 3 5 2 3 4 5;
> 2 6 3 0 0 0;
> 1 6 4 1 2 6;
> 1 5 6 1 6 3;
> 4 5 6 0 0 0;
> 4 6 5 0 0 0;
> 1 6 5 1 3 6;
> 1 4 6 1 6 2;
> 2 3 6 0 0 0;
> 3 2 5 3 5 4;
> 1 2 5 0 0 0;
> 1 4 3 0 0 0;
> 0 0 0 0 0 0]'
>
> new(
> voxCrnrPos,
> voxEdgeCrnrs,
> voxEdgeDir,
> voxEdgeIx,
> subTets,
> subTets,
> tetEdgeCrnrs,
> tetTri)
> end
141,142c140,141
< function hasFaces{T<:Real}(vals::Vector{T}, iso::T)
< if vals[1] < iso
---
> @inline function hasFaces{T<:Real}(vals::Vector{T}, iso::T)
> if vals[1] < iso
144c143
< vals[i] >= iso && return true
---
> vals[i] >= iso && return true
146c145
< else
---
> else
148c147
< vals[i] < iso && return true
---
> vals[i] < iso && return true
150c149
< end
---
> end
155,159c154,158
< function tetIx{T<:Real, IType <: Integer}(tIx::IType, vals::Vector{T}, iso::T, vxidx::VoxelIndices{IType})
< ifelse(vals[vxidx.subTets[1, tIx]] < iso, 1, 0) +
< ifelse(vals[vxidx.subTets[2, tIx]] < iso, 2, 0) +
< ifelse(vals[vxidx.subTets[3, tIx]] < iso, 4, 0) +
< ifelse(vals[vxidx.subTets[4, tIx]] < iso, 8, 0) + 1
---
> @inline function tetIx{T<:Real, IType <: Integer}(tIx::IType, vals::Vector{T}, iso::T, vxidx::VoxelIndexes{IType})
> ifelse(vals[vxidx.subTets[1, tIx]] < iso, 1, 0) +
> ifelse(vals[vxidx.subTets[2, tIx]] < iso, 2, 0) +
> ifelse(vals[vxidx.subTets[3, tIx]] < iso, 4, 0) +
> ifelse(vals[vxidx.subTets[4, tIx]] < iso, 8, 0) + 1
167,168c166,167
< function vertId{IType <: Integer}(e::IType, x::IType, y::IType, z::IType,
< nx::IType, ny::IType, vxidx::VoxelIndices{IType})
---
> @inline function vertId{IType <: Integer}(e::IType, x::IType, y::IType, z::IType,
> nx::IType, ny::IType, vxidx::VoxelIndexes{IType})
178c177
< vals::Vector{T}, iso::T, eps::T, vxidx::VoxelIndices{IType})
---
> vals::Vector{T}, iso::T, eps::T, vxidx::VoxelIndexes{IType})
188,189c187,188
<
< Point{3,T}(
---
>
> Point{3, T}(
198,202c197,203
< function getVertId{T<:Real, IType <: Integer}(e::IType, x::IType, y::IType, z::IType,
< nx::IType, ny::IType,
< vals::Vector{T}, iso::T,
< vts::Dict{IType, Point{3,T}},
< eps::T, vxidx::VoxelIndices{IType})
---
> function getVertId{T<:Real, IType <: Integer}(
> e::IType, x::IType, y::IType, z::IType,
> nx::IType, ny::IType,
> vals::Vector{T}, iso::T,
> vts::Dict{IType, Point{3, T}},
> eps::T, vxidx::VoxelIndexes{IType}
> )
213c214
< function voxEdgeId{IType <: Integer}(subTetIx::IType, tetEdgeIx::IType, vxidx::VoxelIndices{IType})
---
> @inline function voxEdgeId{IType <: Integer}(subTetIx::IType, tetEdgeIx::IType, vxidx::VoxelIndexes{IType})
221,225c222,228
< function procVox{T<:Real, IType <: Integer}(vals::Vector{T}, iso::T,
< x::IType, y::IType, z::IType,
< nx::IType, ny::IType,
< vts::Dict{IType, Point{3,T}}, fcs::Vector{Face{3,IType,0}},
< eps::T, vxidx::VoxelIndices{IType})
---
> function procVox{T<:Real, IType <: Integer}(
> vals::Vector{T}, iso::T,
> x::IType, y::IType, z::IType,
> nx::IType, ny::IType,
> vts::Dict{IType, Point{3, T}}, fcs::Vector{Triangle{IType}},
> eps::T, vxidx::VoxelIndexes{IType}
> )
233c236
< e1 == 0 && break
---
> e1 == 0 && break
238c241
< push!(fcs, Face{3,IType,0}(
---
> push!(fcs, Triangle{IType}(
246,250c249,250
<
< """
< Given a 3D array and an isovalue, extracts a mesh represention of the
< an approximate isosurface by the method of marching tetrahedra.
< """
---
> # Given a 3D array and an isovalue, extracts a mesh represention of the
> # an approximate isosurface by the method of marching tetrahedra.
252,253c252,253
< vts = Dict{indextype, Point{3,T}}()
< fcs = Array(Face{3,indextype,0}, 0)
---
> vts = Dict{indextype, Point{3, T}}()
> fcs = Array(Triangle{indextype}, 0)
256c256
< const vxidx = VoxelIndices{indextype}()
---
> const vxidx = VoxelIndexes{indextype}()
272c272
< function isosurface(lsf, isoval, eps, indextype=Int, index_start=one(Int))
---
> function isosurface(lsf, isoval, eps, indextype=Cuint, index_start=zero(Int))
283c283
< fcAry = Face{3,indextype, index_start-1}[Face{3,indextype, index_start-1}(vtD[f[1]], vtD[f[2]], vtD[f[3]]) for f in fcs]
---
> fcAry = Face{3, indextype, -1}[Face{3, indextype, -1}(vtD[f[1]], vtD[f[2]], vtD[f[3]]) for f in fcs]
286c286
< (vtAry, fcAry)
---
> vtAry, fcAry
295c295
< vts, fcs = isosurface(volume, iso_val, eps_val)
---
> vts, fcs = isosurface(volume, iso_val, eps_val, eltype(facetype(MT)))
298,299d297
<
< end # modulefrom meshing.jl.
Related Issues (20)
- Gracefully handle requests for integer points HOT 1
- Isocaps HOT 2
- Helpful Errors
- Regression in Naive Surface Nets HOT 1
- What is in the name? HOT 1
- Default to Marching Cubes
- varg passed incorrectly in GeometryTypes API
- insidepositive results bugged
- isosurface(x,y,z,V) and rotations
- Makie error when plotting meshes HOT 11
- documentation? HOT 2
- GeometryBasics Normal Calculations/general API for specifying normals
- Example:gyroid does not run HOT 3
- Order of points in triangles HOT 3
- TagBot trigger issue HOT 5
- Naive methods and parallelizations? HOT 1
- Fully de-duplicate vertices in MarchingCube algorithm HOT 2
- GeometryBasics example errors HOT 7
- I could contribute Marching Tetrahedra on Tetmeshes. HOT 5
- Isolines of a 2d function, e.g, (x,y) -> x^2-y^2 HOT 1
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from meshing.jl.