implicit-0.0.1: Graphics/Implicit/Export/MarchingCubes.hs
-- Implicit CAD. Copyright (C) 2011, Christopher Olah (chris@colah.ca)
-- Released under the GNU GPL, see LICENSE
module Graphics.Implicit.Export.MarchingCubes (getMesh, getMesh2) where
import Graphics.Implicit.Definitions
import Control.Parallel (par, pseq)
-- | getMesh gets a triangle mesh describe the boundary of your 3D
-- object.
-- There are many getMesh functions in this file. THis one is the
-- simplest and should be least bug prone. Use it for debugging.
getMesh :: ℝ3 -> ℝ3 -> ℝ -> Obj3 -> TriangleMesh
getMesh (x1, y1, z1) (x2, y2, z2) res obj =
let
-- How many steps will we take on each axis?
nx = fromIntegral $ ceiling $ (x2 - x1) / res
ny = fromIntegral $ ceiling $ (y2 - y1) / res
nz = fromIntegral $ ceiling $ (y2 - y1) / res
-- Divide it up and compute the polylines
triangles :: [TriangleMesh]
triangles = [getCubeTriangles
(x1 + (x2 - x1)*mx/nx, y1 + (y2 - y1)*my/ny, z1 + (z2 - z1)*mz/nz)
(x1 + (x2 - x1)*(mx+1)/nx, y1 + (y2 - y1)*(my+1)/ny, z1 + (z2 - z1)*(mz+1)/nz)
obj
| mx <- [0.. nx-1], my <- [0..ny-1], mz <- [0..nz-1] ]
in
concat $ triangles
getMesh2 (x1,y1,z1) (x2,y2,z2) res obj =
let
dx = abs $ x2 - x1
dy = abs $ y2 - y1
dz = abs $ z2 - z1
d = maximum [dx, dy, dz]
ffloor = fromIntegral . floor
fceil = fromIntegral . ceiling
in
if (abs.obj) ( (x1 + x2)/2, (y1 + y2)/2, (z1 + z2)/2) > d*0.9 then []
else
if d <= res
then getCubeTriangles (x1,y1,z1) (x1+res,y1+res,z1+res) obj
else let
xs = if dx <= res then [(x1, x2)] else [(x1,xm), (xm, x2)]
where xm = x1 + res * fceil ( ffloor (dx/res) / 2.0)
ys = if dy <= res then [(y1, y2)] else [(y1,xm), (xm, y2)]
where xm = y1 + res * fceil ( ffloor (dy/res) / 2.0)
zs = if dz <= res then [(z1, z2)] else [(z1,xm), (xm, z2)]
where xm = z1 + res * fceil ( ffloor (dz/res) / 2.0)
partitions = [getMesh (x1', y1', z1') (x2', y2', z2') res obj
| (x1',x2') <- xs, (y1', y2') <- ys, (z1',z2') <- zs ]
in
concat partitions
getMesh3 (x1,y1,z1) (x2,y2,z2) res obj =
let
dx = abs $ x2 - x1
dy = abs $ y2 - y1
dz = abs $ z2 - z1
d = maximum [dx, dy, dz]
ffloor = fromIntegral . floor
fceil = fromIntegral . ceiling
in
if (abs.obj) ( (x1 + x2)/2, (y1 + y2)/2, (z1 + z2)/2) > d*0.9 then []
else
if d <= res
then getCubeTriangles (x1,y1,z1) (x1+res,y1+res,z1+res) obj
else let
xs = if dx <= res then [(x1, x2)] else [(x1,xm), (xm, x2)]
where xm = x1 + res * fceil ( ffloor (dx/res) / 2.0)
ys = if dy <= res then [(y1, y2)] else [(y1,xm), (xm, y2)]
where xm = y1 + res * fceil ( ffloor (dy/res) / 2.0)
zs = if dz <= res then [(z1, z2)] else [(z1,xm), (xm, z2)]
where xm = z1 + res * fceil ( ffloor (dz/res) / 2.0)
partitions = [getMesh (x1', y1', z1') (x2', y2', z2') res obj
| (x1',x2') <- xs, (y1', y2') <- ys, (z1',z2') <- zs ]
in
foldr1 par partitions `pseq` concat partitions
-- | This monstrosity of a function gives triangles to divde negative interior
-- regions and positive exterior ones inside a cube, based on its vertices.
-- It is based on the linearly-interpolated marching cubes algorithm.
getCubeTriangles :: ℝ3 -> ℝ3 -> Obj3 -> [Triangle]
getCubeTriangles (x1, y1, z1) (x2, y2, z2) obj =
let
(x,y,z) = (x1, y1, z1)
x1y1z1 = obj (x1, y1, z1)
x2y1z1 = obj (x2, y1, z1)
x1y2z1 = obj (x1, y2, z1)
x2y2z1 = obj (x2, y2, z1)
x1y1z2 = obj (x1, y1, z2)
x2y1z2 = obj (x2, y1, z2)
x1y2z2 = obj (x1, y2, z2)
x2y2z2 = obj (x2, y2, z2)
dx = x2 - x1
dy = y2 - y1
dz = z2 - z1
--{- Linearly interpolated
x1y1 = (x, y, z+dz*x1y1z1/(x1y1z1-x1y1z2))
x1y2 = (x, y+dy, z+dz*x1y2z1/(x1y2z1-x1y2z2))
x2y1 = (x+dx, y, z+dz*x2y1z1/(x2y1z1-x2y1z2))
x2y2 = (x+dx, y+dy, z+dz*x2y2z1/(x2y2z1-x2y2z2))
x1z1 = (x, y+dy*x1y1z1/(x1y1z1-x1y2z1), z)
x1z2 = (x, y+dy*x1y1z2/(x1y1z2-x1y2z2), z+dz)
x2z1 = (x+dx, y+dy*x2y1z1/(x2y1z1-x2y2z1), z)
x2z2 = (x+dx, y+dy*x2y1z2/(x2y1z2-x2y2z2), z+dz)
y1z1 = (x+dx*x1y1z1/(x1y1z1-x2y1z1), y, z)
y1z2 = (x+dx*x1y1z2/(x1y1z2-x2y1z2), y, z+dz)
y2z1 = (x+dx*x1y2z1/(x1y2z1-x2y2z1), y+dy, z)
y2z2 = (x+dx*x1y2z2/(x1y2z2-x2y2z2), y+dy, z+dz)
--}
{- Non-linearly interpolated
x1y1 = (x, y, z+dz/2)
x1y2 = (x, y+dy, z+dz/2)
x2y1 = (x+dx, y, z+dz/2)
x2y2 = (x+dx, y+dy, z+dz/2)
x1z1 = (x, y+dy/2, z)
x1z2 = (x, y+dy/2, z+dz)
x2z1 = (x+dx, y+dy/2, z)
x2z2 = (x+dx, y+dy/2, z+dz)
y1z1 = (x+dx/2, y, z)
y1z2 = (x+dx/2, y, z+dz)
y2z1 = (x+dx/2, y+dy,z)
y2z2 = (x+dx/2, y+dy,z+dz)
--}
-- Convenience function
square a b c d = [(a,b,c),(d,a,c)]
in case
-- whether the vertices are "in" or "out" form the topological
-- basis of our triangles constructions. We must consider every
-- possible case.
-- We arrange the vertices in a human readable way
-- BOTTOM LAYER TOP LAYER
(x1y2z1<=0, x2y2z1<=0, x1y2z2<=0, x2y2z2<=0,
x1y1z1<=0, x2y1z1<=0, x1y1z2<=0, x2y1z2<=0)
of
-- There are 256 cases to implement.
-- Only about half are, but they're the most common ones.
-- In practice, this has no issues redering reasonable objects.
-- Yes, there's some symetries that could reduce the amount of code...
-- But I don't think they're worth exploiting...
-- In particular, since we're not implementing any case,
-- it would make catching the ones we don't implement... problematic.
-- Uniform cases = empty
(False,False, False,False,
False,False, False,False) -> []
(True, True, True, True,
True, True, True, True ) -> []
-- 2 uniform layers
(True, True, False,False,
True, True, False,False) -> square x1y1 x2y1 x2y2 x1y2
(False,False, True, True,
False,False, True, True ) -> square x1y1 x2y1 x2y2 x1y2
(True, True, True, True,
False,False, False,False) -> square x1z1 x2z1 x2z2 x1z2
(False,False, False,False,
True, True, True, True ) -> square x1z1 x2z1 x2z2 x1z2
(False,True, False,True,
False,True, False,True ) -> square y1z1 y2z1 y2z2 y1z2
(True, False, True, False,
True, False, True, False) -> square y1z1 y2z1 y2z2 y1z2
-- z single column
(True, False, True, False,
False,False, False,False) -> square x1z1 y2z1 y2z2 x1z2
(False,True, False,True,
False,False, False,False) -> square x2z1 y2z1 y2z2 x2z2
(False,False, False,False,
True, False, True, False) -> square x1z1 y1z1 y1z2 x1z2
(False,False, False,False,
False,True, False,True ) -> square y1z1 x2z1 x2z2 y1z2
(False,True, False,True,
True, True, True, True ) -> square x1z1 y2z1 y2z2 x1z2
(True, False, True, False,
True, True, True, True ) -> square x2z1 y2z1 y2z2 x2z2
(True, True, True, True,
False,True, False,True ) -> square x1z1 y1z1 y1z2 x1z2
(True, True, True, True,
True, False, True, False) -> square y1z1 x2z1 x2z2 y1z2
-- single y column
(True, False, False,False,
True, False, False,False) -> square x1y1 y1z1 y2z1 x1y2
(False,True, False,False,
False,True, False,False) -> square x2y1 y1z1 y2z1 x2y2
(False,False, True, False,
False,False, True, False) -> square x1y1 y1z2 y2z2 x1y2
(False,False, False,True,
False,False, False,True ) -> square x2y1 y1z2 y2z2 x2y2
(False,True, True, True,
False,True, True, True) -> square x1y1 y1z1 y2z1 x1y2
(True, False, True, True,
True, False, True, True) -> square x2y1 y1z1 y2z1 x2y2
(True, True, False, True,
True, True, False, True) -> square x1y1 y1z2 y2z2 x1y2
(True, True, True, False,
True, True, True, False) -> square x2y1 y1z2 y2z2 x2y2
-- since x column
(True, True, False,False,
False,False, False,False) -> square x1y2 x1z1 x2z1 x2y2
(False,False, False,False,
True, True, False,False) -> square x1y1 x1z1 x2z1 x2y1
(False,False, True, True,
False,False, False,False) -> square x1y2 x1z2 x2z2 x2y2
(False,False, False,False,
False,False, True, True ) -> square x1y1 x1z2 x2z2 x2y1
(False,False, True, True,
True, True, True, True ) -> square x1y2 x1z1 x2z1 x2y2
(True, True, True, True,
False,False, True, True ) -> square x1y1 x1z1 x2z1 x2y1
(True, True, False,False,
True, True, True, True ) -> square x1y2 x1z2 x2z2 x2y2
(True, True, True, True,
True, True, False,False) -> square x1y1 x1z2 x2z2 x2y1
-- lone points
(True, False, False,False,
False,False, False,False) -> [(x1z1, y2z1, x1y2)]
(False,True, False,False,
False,False, False,False) -> [(x2z1, y2z1, x2y2)]
(False,False, False,False,
True, False, False,False) -> [(x1z1, y1z1, x1y1)]
(False,False, False,False,
False,True, False,False) -> [(x2z1, y1z1, x2y1)]
(False,False, True, False,
False,False, False,False) -> [(x1z2, y2z2, x1y2)]
(False,False, False,True,
False,False, False,False) -> [(x2z2, y2z2, x2y2)]
(False,False, False,False,
False,False, True, False) -> [(x1z2, y1z2, x1y1)]
(False,False, False,False,
False,False, False,True ) -> [(x2z2, y1z2, x2y1)]
(False,True, True, True,
True, True, True, True ) -> [(x1z1, y2z1, x1y2)]
(True, False, True, True,
True, True, True, True ) -> [(x2z1, y2z1, x2y2)]
(True, True, True, True,
False,True, True, True ) -> [(x1z1, y1z1, x1y1)]
(True, True, True, True,
True, False, True, True ) -> [(x2z1, y1z1, x2y1)]
(True, True, False,True,
True, True, True, True ) -> [(x1z2, y2z2, x1y2)]
(True, True, True, False,
True, True, True, True ) -> [(x2z2, y2z2, x2y2)]
(True, True, True, True,
True, True, False,True ) -> [(x1z2, y1z2, x1y1)]
(True, True, True, True,
True, True, True, False) -> [(x2z2, y1z2, x2y1)]
-- z flat + 1
(False,False, True, False,
False,False, True, True) -> [(x1y1,x2y1,x2z2), (x1y1,x2z2,y2z2), (x1y1,y2z2,x1y2)]
(True, True, False,True,
True, True, False,False) -> [(x1y1,x2y1,x2z2), (x1y1,x2z2,y2z2), (x1y1,y2z2,x1y2)]
(False,False, False,True,
False,False, True, True) -> [(x2y1,x1y1,x1z2), (x2y1,x1z2,y2z2), (x2y1,y2z2,x2y2)]
(True, True, True, False,
True, True, False,False) -> [(x2y1,x1y1,x1z2), (x2y1,x1z2,y2z2), (x2y1,y2z2,x2y2)]
(False,False, True, True,
False,False, True, False) -> [(x1y2,x2y2,x2z2), (x1y2,x2z2,y1z2), (x1y2,y1z2,x1y1)]
(True, True, False,False,
True, True, False,True ) -> [(x1y2,x2y2,x2z2), (x1y2,x2z2,y1z2), (x1y2,y1z2,x1y1)]
(False,False, True, True,
False,False, False,True) -> [(x2y2,x1y2,x1z2), (x2y2,x1z2,y1z2), (x2y2,y1z2,x2y1)]
(True, True, False,False,
True, True, True, False) -> [(x2y2,x1y2,x1z2), (x2y2,x1z2,y1z2), (x2y2,y1z2,x2y1)]
(True, False, False,False,
True, True, False,False) -> [(x1y1,x2y1,x2z1), (x1y1,x2z1,y2z1), (x1y1,y2z1,x1y2)]
(False,True, True, True,
False,False, True, True) -> [(x1y1,x2y1,x2z1), (x1y1,x2z1,y2z1), (x1y1,y2z1,x1y2)]
(False,True, False,False,
True, True, False,False) -> [(x2y1,x1y1,x1z1), (x2y1,x1z1,y2z1), (x2y1,y2z1,x2y2)]
(True, False, True, True,
False,False, True, True) -> [(x2y1,x1y1,x1z1), (x2y1,x1z1,y2z1), (x2y1,y2z1,x2y2)]
(True, True, False,False,
True, False, False,False) -> [(x1y2,x2y2,x2z1), (x1y2,x2z1,y1z1), (x1y2,y1z1,x1y1)]
(False,False, True, True,
False,True, True, True) -> [(x1y2,x2y2,x2z1), (x1y2,x2z1,y1z1), (x1y2,y1z1,x1y1)]
(True, True, False,False,
False,True, False,False) -> [(x2y2,x1y2,x1z1), (x2y2,x1z1,y1z1), (x2y2,y1z1,x2y1)]
(False,False, True, True,
True, False, True, True) -> [(x2y2,x1y2,x1z1), (x2y2,x1z1,y1z1), (x2y2,y1z1,x2y1)]
-- y flat + 1
(True, False, True, True,
True, False, True, False) -> [(y2z1,x2y2,x2z2),(y2z1,x2z2,y1z1),(y1z1,x2z2,y1z2)]
(False,True, False,False,
False,True, False,True ) -> [(y2z1,x2y2,x2z2),(y2z1,x2z2,y1z1),(y1z1,x2z2,y1z2)]
(True, False, True, False,
True, False, True, True ) -> [(y1z1,x2y1,x2z2),(y1z1,x2z2,y2z1),(y2z1,x2z2,y2z2)]
(False,True, False,True,
False,True, False,False) -> [(y1z1,x2y1,x2z2),(y1z1,x2z2,y2z1),(y2z1,x2z2,y2z2)]
(False,True, True, True,
False,True, False,True ) -> [(y2z1,x1y2,x1z2),(y2z1,x1z2,y1z1),(y1z1,x1z2,y1z2)]
(True, False, False,False,
True, False, True, False) -> [(y2z1,x1y2,x1z2),(y2z1,x1z2,y1z1),(y1z1,x1z2,y1z2)]
(False,True, False,True,
False,True, True, True ) -> [(y1z1,x1y1,x1z2),(y1z1,x1z2,y2z1),(y2z1,x1z2,y2z2)]
(True, False, True, False,
True, False, False,False) -> [(y1z1,x1y1,x1z2),(y1z1,x1z2,y2z1),(y2z1,x1z2,y2z2)]
(True, True, True, False,
True, False, True, False) -> [(y2z2,x2y2,x2z1),(y2z2,x2z1,y1z2),(y1z2,x2z1,y1z1)]
(False,False, False,True,
False,True, False,True ) -> [(y2z2,x2y2,x2z1),(y2z2,x2z1,y1z2),(y1z2,x2z1,y1z1)]
(True, False, True, False,
True, True, True, False) -> [(y1z2,x2y1,x2z1),(y1z2,x2z1,y2z2),(y2z2,x2z1,y2z1)]
(False,True, False,True,
False,False, False,True) -> [(y1z2,x2y1,x2z1),(y1z2,x2z1,y2z2),(y2z2,x2z1,y2z1)]
(True, True, False,True,
False,True, False,True) -> [(y2z2,x1y2,x1z1),(y2z2,x1z1,y1z2),(y1z2,x1z1,y1z1)]
(False,False, True, False,
True, False, True, False) -> [(y2z2,x1y2,x1z1),(y2z2,x1z1,y1z2),(y1z2,x1z1,y1z1)]
(False,True, False,True,
True, True, False,True) -> [(y1z2,x1y1,x1z1),(y1z2,x1z1,y2z2),(y2z2,x1z1,y2z1)]
(True, False, True, False,
False,False, True, False) -> [(y1z2,x1y1,x1z1),(y1z2,x1z1,y2z2),(y2z2,x1z1,y2z1)]
-- x flat +1
(True, True, True, True,
False,False, True, False) -> [(x1z1,x2z1,x1y1),(x1y1,x2z1,x2z2),(x1y1,x2z2,y1z2)]
(False,False, False,False,
True, True, False,True ) -> [(x1z1,x2z1,x1y1),(x1y1,x2z1,x2z2),(x1y1,x2z2,y1z2)]
(False,False, True, False,
True, True, True, True) -> [(x1z1,x2z1,x1y2),(x1y2,x2z1,x2z2),(x1y2,x2z2,y2z2)]
(True, True, False,True,
False,False, False,False) -> [(x1z1,x2z1,x1y2),(x1y2,x2z1,x2z2),(x1y2,x2z2,y2z2)]
(True, True, True, True,
False,False, False,True) -> [(x2z1,x1z1,x2y1),(x2y1,x1z1,x1z2),(x2y1,x1z2,y1z2)]
(False,False, False,False,
True, True, True, False) -> [(x2z1,x1z1,x2y1),(x2y1,x1z1,x1z2),(x2y1,x1z2,y1z2)]
(False,False, False,True,
True, True, True, True) -> [(x2z1,x1z1,x2y2),(x2y2,x1z1,x1z2),(x2y2,x1z2,y2z2)]
(True, True, True, False,
False,False, False,False) -> [(x2z1,x1z1,x2y2),(x2y2,x1z1,x1z2),(x2y2,x1z2,y2z2)]
(True, True, True, True,
True, False, False,False) -> [(x1z2,x2z2,x1y1),(x1y1,x2z2,x2z1),(x1y1,x2z1,y1z1)]
(False,False, False,False,
False,True, True, True ) -> [(x1z2,x2z2,x1y1),(x1y1,x2z2,x2z1),(x1y1,x2z1,y1z1)]
(True, False, False,False,
True, True, True, True ) -> [(x1z2,x2z2,x1y2),(x1y2,x2z2,x2z1),(x1y2,x2z1,y2z1)]
(False,True, True, True,
False,False, False,False) -> [(x1z2,x2z2,x1y2),(x1y2,x2z2,x2z1),(x1y2,x2z1,y2z1)]
(True, True, True, True,
False,True, False,False) -> [(x2z2,x1z2,x2y1),(x2y1,x1z2,x1z1),(x2y1,x1z1,y1z1)]
(False,False, False,False,
True, False, True, True) -> [(x2z2,x1z2,x2y1),(x2y1,x1z2,x1z1),(x2y1,x1z1,y1z1)]
(False,True, False,False,
True, True, True, True) -> [(x2z2,x1z2,x2y2),(x2y2,x1z2,x1z1),(x2y2,x1z1,y2z1)]
(True, False, True, True,
False,False, False,False) -> [(x2z2,x1z2,x2y2),(x2y2,x1z2,x1z1),(x2y2,x1z1,y2z1)]
(True, True, True, False,
True, False, False,False) -> [(x1y1,x1z2,y1z1),(y1z1,x1z2,y2z2),(y1z1,y2z2,x2z1),(x2z1,y2z2,x2y2)]
(False,False, False,True,
False,True, True, True ) -> [(x1y1,x1z2,y1z1),(y1z1,x1z2,y2z2),(y1z1,y2z2,x2z1),(x2z1,y2z2,x2y2)]
(True, True, False,True,
False,True, False,False) -> [(x2y1,x2z2,y1z1),(y1z1,x2z2,y2z2),(y1z1,y2z2,x1z1),(x1z1,y2z2,x1y2)]
(False,False, True, False,
True, False, True, True ) -> [(x2y1,x2z2,y1z1),(y1z1,x2z2,y2z2),(y1z1,y2z2,x1z1),(x1z1,y2z2,x1y2)]
(True, False, False,False,
True, True, True, False) -> [(x1y2,x1z2,y2z1),(y2z1,x1z2,y1z2),(y2z1,y1z2,x2z1),(x2z1,y1z2,x2y1)]
(False,True, True, True,
False,False, False,True ) -> [(x1y2,x1z2,y2z1),(y2z1,x1z2,y1z2),(y2z1,y1z2,x2z1),(x2z1,y1z2,x2y1)]
(False,True, False,False,
True, True, False,True ) -> [(x2y2,x2z2,y2z1),(y2z1,x2z2,y1z2),(y2z1,y1z2,x1z1),(x1z1,y1z2,x1y1)]
(True, False, True, True,
False,False, True, False) -> [(x2y2,x2z2,y2z1),(y2z1,x2z2,y1z2),(y2z1,y1z2,x1z1),(x1z1,y1z2,x1y1)]
_ -> []