{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE FlexibleContexts #-}
module MarchingCubes.MarchingCubes
( Voxel
, XYZ
, marchingCubes
, makeVoxel
) where
import Data.Array.Unboxed ( IArray
, UArray
, listArray
)
import qualified Data.Foldable as F
import Data.List ( elemIndices
, findIndices
)
import Data.Matrix ( (<->)
, Matrix(ncols)
, getCol
, getRow
, mapCol
, submatrix
)
import qualified Data.Matrix as M
import Data.Maybe ( catMaybes
, mapMaybe
)
import qualified Data.Vector as V
import Data.Vector.Unboxed ( (!)
, Unbox
)
import qualified Data.Vector.Unboxed as UV
import MarchingCubes.Internal ( calPoints
, faces7
, facesNo7
, getBasic1
, getBasic2
, getPoints
, getR
, getTcase
, levCells
)
import MarchingCubes.Tables ( edgePoints
, edgesLengths
, edgesTable
, edgesTable2
, facesTable
, specialInd
, specialName
, specialNedge
, specialNface
, specialPos
)
import MarchingCubes.Utils ( cbind
, jthColumn
, matrix2listMinusFirstColumn
, replicateEach
, replicateEach'
, subMatrix
, vector2matrix
)
type Bounds a = ((a, a), (a, a), (a, a))
type Dims = (Int, Int, Int)
type Voxel a = ((UArray Dims a, a), (Bounds a, Dims))
type XYZ a = (a, a, a)
makeVoxel
:: (RealFloat a, IArray UArray a)
=> (XYZ a -> a)
-> Bounds a
-> Dims
-> Voxel a
makeVoxel :: forall a.
(RealFloat a, IArray UArray a) =>
(XYZ a -> a) -> Bounds a -> Dims -> Voxel a
makeVoxel XYZ a -> a
fun bds :: Bounds a
bds@((a
xm, a
xM), (a
ym, a
yM), (a
zm, a
zM)) dims :: Dims
dims@(Int
nx, Int
ny, Int
nz) =
((forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
(i, i) -> [e] -> a i e
listArray ((Int
0, Int
0, Int
0), (Int
nx forall a. Num a => a -> a -> a
- Int
1, Int
ny forall a. Num a => a -> a -> a
- Int
1, Int
nz forall a. Num a => a -> a -> a
- Int
1)) [a]
values, a
mxmm), (Bounds a
bds, Dims
dims))
where
x_ :: [a]
x_ = [ a
xm forall a. Num a => a -> a -> a
+ (a
xM forall a. Num a => a -> a -> a
- a
xm) forall a. Num a => a -> a -> a
* forall {a} {a}. (Fractional a, Real a) => a -> a
fracx Int
i | Int
i <- [Int
0 .. Int
nx forall a. Num a => a -> a -> a
- Int
1] ]
fracx :: a -> a
fracx a
p = forall a b. (Real a, Fractional b) => a -> b
realToFrac a
p forall a. Fractional a => a -> a -> a
/ (forall a b. (Real a, Fractional b) => a -> b
realToFrac Int
nx forall a. Num a => a -> a -> a
- a
1)
y_ :: [a]
y_ = [ a
ym forall a. Num a => a -> a -> a
+ (a
yM forall a. Num a => a -> a -> a
- a
ym) forall a. Num a => a -> a -> a
* forall {a} {a}. (Fractional a, Real a) => a -> a
fracy Int
i | Int
i <- [Int
0 .. Int
ny forall a. Num a => a -> a -> a
- Int
1] ]
fracy :: a -> a
fracy a
p = forall a b. (Real a, Fractional b) => a -> b
realToFrac a
p forall a. Fractional a => a -> a -> a
/ (forall a b. (Real a, Fractional b) => a -> b
realToFrac Int
ny forall a. Num a => a -> a -> a
- a
1)
z_ :: [a]
z_ = [ a
zm forall a. Num a => a -> a -> a
+ (a
zM forall a. Num a => a -> a -> a
- a
zm) forall a. Num a => a -> a -> a
* forall {a} {a}. (Fractional a, Real a) => a -> a
fracz Int
i | Int
i <- [Int
0 .. Int
nz forall a. Num a => a -> a -> a
- Int
1] ]
fracz :: a -> a
fracz a
p = forall a b. (Real a, Fractional b) => a -> b
realToFrac a
p forall a. Fractional a => a -> a -> a
/ (forall a b. (Real a, Fractional b) => a -> b
realToFrac Int
nz forall a. Num a => a -> a -> a
- a
1)
values :: [a]
values = [ XYZ a -> a
fun (a
x, a
y, a
z) | a
x <- [a]
x_, a
y <- [a]
y_, a
z <- [a]
z_ ]
mxmm :: a
mxmm = forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (forall a. (a -> Bool) -> [a] -> [a]
filter (Bool -> Bool
not forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. RealFloat a => a -> Bool
isNaN) [a]
values)
rescale :: Fractional a => (a, a) -> Int -> a -> a
rescale :: forall a. Fractional a => (a, a) -> Int -> a -> a
rescale (a
minmm, a
maxmm) Int
n a
w = a
minmm forall a. Num a => a -> a -> a
+ (a
maxmm forall a. Num a => a -> a -> a
- a
minmm) forall a. Num a => a -> a -> a
* a
w forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
n forall a. Num a => a -> a -> a
- Int
1)
rescaleMatrix :: Fractional a => Matrix a -> Bounds a -> Dims -> Matrix a
rescaleMatrix :: forall a. Fractional a => Matrix a -> Bounds a -> Dims -> Matrix a
rescaleMatrix Matrix a
mtrx ((a, a)
xbds, (a, a)
ybds, (a, a)
zbds) (Int
nx, Int
ny, Int
nz) = Matrix a
mtrx'''
where
mtrx' :: Matrix a
mtrx' = forall a. (Int -> a -> a) -> Int -> Matrix a -> Matrix a
mapCol (\Int
_ a
w -> forall a. Fractional a => (a, a) -> Int -> a -> a
rescale (a, a)
xbds Int
nx (a
w forall a. Num a => a -> a -> a
- a
1)) Int
1 Matrix a
mtrx
mtrx'' :: Matrix a
mtrx'' = forall a. (Int -> a -> a) -> Int -> Matrix a -> Matrix a
mapCol (\Int
_ a
w -> forall a. Fractional a => (a, a) -> Int -> a -> a
rescale (a, a)
ybds Int
ny (a
w forall a. Num a => a -> a -> a
- a
1)) Int
2 Matrix a
mtrx'
mtrx''' :: Matrix a
mtrx''' = forall a. (Int -> a -> a) -> Int -> Matrix a -> Matrix a
mapCol (\Int
_ a
w -> forall a. Fractional a => (a, a) -> Int -> a -> a
rescale (a, a)
zbds Int
nz (a
w forall a. Num a => a -> a -> a
- a
1)) Int
3 Matrix a
mtrx''
marchingCubes
:: (RealFloat a, Unbox a, IArray UArray a) => Voxel a -> a -> Matrix a
marchingCubes :: forall a.
(RealFloat a, Unbox a, IArray UArray a) =>
Voxel a -> a -> Matrix a
marchingCubes ((UArray Dims a
voxel, a
mx), (Bounds a
bds, Dims
dims)) a
level = forall a. Fractional a => Matrix a -> Bounds a -> Dims -> Matrix a
rescaleMatrix
(forall b a. b -> (a -> b) -> Maybe a -> b
maybe Matrix a
triangles1 (Matrix a
triangles1 forall a. Matrix a -> Matrix a -> Matrix a
<->) Maybe (Matrix a)
triangles2)
Bounds a
bds
Dims
dims
where
ijkt :: Matrix Int
ijkt = forall a.
(Real a, IArray UArray a) =>
UArray Dims a -> a -> a -> Matrix Int
levCells UArray Dims a
voxel a
level a
mx
vt :: Vector Int
vt = forall a. Int -> Matrix a -> Vector a
getRow Int
4 Matrix Int
ijkt
tcase :: Vector Int
tcase = Vector Int -> Vector Int
getTcase Vector Int
vt
r :: Vector Int
r = Vector Int -> Vector Int
getR Vector Int
tcase
nR :: Int
nR = forall a. Unbox a => Vector a -> Int
UV.length Vector Int
r
ijk :: Matrix Int
ijk = forall a. Int -> Int -> Int -> Int -> Matrix a -> Matrix a
submatrix Int
1 Int
3 Int
1 (forall a. Matrix a -> Int
ncols Matrix Int
ijkt) Matrix Int
ijkt
vivjvk :: Matrix Int
vivjvk = forall a. Matrix a -> Matrix a
M.transpose Matrix Int
ijk
cubeco :: Matrix Int
cubeco = Vector Int -> Matrix Int -> Matrix Int
getBasic1 Vector Int
r Matrix Int
vivjvk
values :: Vector a
values = forall a.
(Num a, Unbox a, IArray UArray a) =>
UArray Dims a -> a -> Matrix Int -> Vector a
getBasic2 UArray Dims a
voxel a
level Matrix Int
cubeco
p1 :: [Int]
p1 = [ Int
8 forall a. Num a => a -> a -> a
* Int
i forall a. Num a => a -> a -> a
+ Int
1 | Int
i <- [Int
0 .. Int
nR forall a. Num a => a -> a -> a
- Int
1] ]
cases :: Vector Int
cases = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map (\Int
j -> Vector Int
vt forall a. Vector a -> Int -> a
V.! Int
j forall a. Num a => a -> a -> a
- Int
1) Vector Int
r
edgeslengths :: Vector Int
edgeslengths = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map (Vector Int
edgesLengths forall a. Unbox a => Vector a -> Int -> a
UV.!) Vector Int
cases
p1rep :: [Int]
p1rep = forall (t :: * -> *) a. Foldable t => t a -> [a]
F.toList forall a b. (a -> b) -> a -> b
$ forall a. [a] -> [Int] -> Seq a
replicateEach [Int]
p1 (forall a. Unbox a => Vector a -> [a]
UV.toList Vector Int
edgeslengths)
edges :: [Int]
edges =
[ (Vector (Vector Int)
edgesTable forall a. Vector a -> Int -> a
V.! (Vector Int
cases forall a. Unbox a => Vector a -> Int -> a
! Int
i)) forall a. Unbox a => Vector a -> Int -> a
! Int
j
| Int
i <- [Int
0 .. Int
nR forall a. Num a => a -> a -> a
- Int
1]
, Int
j <- [Int
0 .. Vector Int
edgeslengths forall a. Unbox a => Vector a -> Int -> a
! Int
i forall a. Num a => a -> a -> a
- Int
1]
]
epoints :: [Vector Int]
epoints = forall a b. (a -> b) -> [a] -> [b]
map (\Int
j -> Vector (Vector Int)
edgePoints forall a. Vector a -> Int -> a
V.! (Int
j forall a. Num a => a -> a -> a
- Int
1)) [Int]
edges
x1 :: [Int]
x1 = forall a b. (a -> b) -> [a] -> [b]
map (forall a. Unbox a => Vector a -> Int -> a
! Int
1) [Vector Int]
epoints
x2 :: [Int]
x2 = forall a b. (a -> b) -> [a] -> [b]
map (forall a. Unbox a => Vector a -> Int -> a
! Int
2) [Vector Int]
epoints
points :: Matrix a
points = forall a.
(Unbox a, RealFrac a) =>
Matrix Int -> Vector a -> [Int] -> [Int] -> [Int] -> Matrix a
getPoints Matrix Int
cubeco Vector a
values [Int]
p1rep [Int]
x1 [Int]
x2
triangles1 :: Matrix a
triangles1 = forall a. Fractional a => Matrix a -> Matrix a
calPoints Matrix a
points
specialCases :: [Vector Int]
specialCases = forall a b. (a -> b) -> [a] -> [b]
map (forall a. (Unbox a, Eq a) => a -> Vector a -> Vector Int
`UV.elemIndices` Vector Int
tcase) [Int]
specialName
r3s :: [Vector Int]
r3s = forall a. (a -> Bool) -> [a] -> [a]
filter (Bool -> Bool
not forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Unbox a => Vector a -> Bool
UV.null) [Vector Int]
specialCases
cs :: [Int]
cs = forall a. (a -> Bool) -> [a] -> [Int]
findIndices (Bool -> Bool
not forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Unbox a => Vector a -> Bool
UV.null) [Vector Int]
specialCases
setOfTriangles :: [Matrix a]
setOfTriangles = forall a. [Maybe a] -> [a]
catMaybes (forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Int -> Vector Int -> Maybe (Matrix a)
special [Int]
cs [Vector Int]
r3s)
special :: Int -> Vector Int -> Maybe (Matrix a)
special Int
c Vector Int
r3 = if forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Matrix a]
newtriangles
then forall a. Maybe a
Nothing
else forall a. a -> Maybe a
Just forall a b. (a -> b) -> a -> b
$ forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldl1 forall a. Matrix a -> Matrix a -> Matrix a
(<->) [Matrix a]
newtriangles
where
nR3 :: Int
nR3 = forall a. Unbox a => Vector a -> Int
UV.length Vector Int
r3
cubeco3 :: Matrix Int
cubeco3 = Vector Int -> Matrix Int -> Matrix Int
getBasic1 Vector Int
r3 Matrix Int
vivjvk
values3 :: Vector a
values3 = forall a.
(Num a, Unbox a, IArray UArray a) =>
UArray Dims a -> a -> Matrix Int -> Vector a
getBasic2 UArray Dims a
voxel a
level Matrix Int
cubeco3
p13 :: Vector Int
p13 = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
UV.map (\Int
i -> Int
8 forall a. Num a => a -> a -> a
* Int
i forall a. Num a => a -> a -> a
+ Int
1) (forall a. (Unbox a, Num a) => a -> Int -> Vector a
UV.enumFromN Int
0 Int
nR3)
cases3 :: [Int]
cases3 = [ Vector Int
vt forall a. Vector a -> Int -> a
V.! (Vector Int
r3 forall a. Unbox a => Vector a -> Int -> a
! Int
i) forall a. Num a => a -> a -> a
- Int
1 | Int
i <- [Int
0 .. Int
nR3 forall a. Num a => a -> a -> a
- Int
1] ]
nedge :: Int
nedge = Vector Int
specialNedge forall a. Unbox a => Vector a -> Int -> a
! Int
c
faces3 :: Vector Int
faces3 = forall a. Unbox a => [Vector a] -> Vector a
UV.concat forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (Vector (Vector Int)
facesTable forall a. Vector a -> Int -> a
V.!) [Int]
cases3
index3 :: [Int]
index3 = case Int
c of
Int
0 -> forall a.
(Num a, Ord a, Unbox a) =>
Vector Int -> Vector Int -> Vector a -> Int -> Int -> [Int]
facesNo7 Vector Int
faces3 Vector Int
p13 Vector a
values3 Int
nR3 Int
1
Int
1 -> forall a.
(RealFloat a, Unbox a) =>
Vector Int -> Vector Int -> Vector a -> Int -> Int -> [Int]
faces7 Vector Int
faces3 Vector Int
p13 Vector a
values3 Int
nR3 Int
1
Int
_ -> [Int]
index3'
where
nface :: Int
nface = Vector Int
specialNface forall a. Unbox a => Vector a -> Int -> a
! Int
c
facelast :: Vector Int
facelast = Vector Int -> Int -> Int -> Vector Int
jthColumn Vector Int
faces3 Int
nface (Int
nface forall a. Num a => a -> a -> a
- Int
1)
index3' :: [Int]
index3' = Int -> [Int] -> [Int]
loop Int
0 (forall a.
(RealFloat a, Unbox a) =>
Vector Int -> Vector Int -> Vector a -> Int -> Int -> [Int]
faces7 Vector Int
facelast Vector Int
p13 Vector a
values3 Int
nR3 Int
nface)
loop :: Int -> [Int] -> [Int]
loop Int
j ![Int]
idx
| Int
j forall a. Eq a => a -> a -> Bool
== Int
nface forall a. Num a => a -> a -> a
- Int
1
= [Int]
idx
| Bool
otherwise
= let facej :: Vector Int
facej = Vector Int -> Int -> Int -> Vector Int
jthColumn Vector Int
faces3 Int
nface Int
j
in let temp :: [Int]
temp = forall a.
(Num a, Ord a, Unbox a) =>
Vector Int -> Vector Int -> Vector a -> Int -> Int -> [Int]
facesNo7 Vector Int
facej Vector Int
p13 Vector a
values3 Int
nR3 (Int
j forall a. Num a => a -> a -> a
+ Int
1)
in Int -> [Int] -> [Int]
loop (Int
j forall a. Num a => a -> a -> a
+ Int
1) (forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
(+) [Int]
idx [Int]
temp)
edges3' :: [Int]
edges3' = forall a. Unbox a => Vector a -> [a]
UV.toList forall a b. (a -> b) -> a -> b
$ forall a. Unbox a => [Vector a] -> Vector a
UV.concat forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (Vector (Vector Int)
edgesTable2 forall a. Vector a -> Int -> a
V.!) [Int]
cases3
edges3 :: Matrix Int
edges3 = forall a. [a] -> Int -> Matrix a
vector2matrix [Int]
edges3' Int
nedge
edgesp1index :: Matrix Int
edgesp1index = forall a. Matrix a -> [a] -> [a] -> Matrix a
cbind Matrix Int
edges3 (forall a. Unbox a => Vector a -> [a]
UV.toList Vector Int
p13) [Int]
index3
ind3 :: Vector Int
ind3 = Vector (Vector Int)
specialInd forall a. Vector a -> Int -> a
V.! Int
c
newtriangles :: [Matrix a]
newtriangles = forall a b. (a -> Maybe b) -> [a] -> [b]
mapMaybe Int -> Maybe (Matrix a)
f [Int
0 .. forall a. Unbox a => Vector a -> Int
UV.length Vector Int
ind3 forall a. Num a => a -> a -> a
- Int
1]
f :: Int -> Maybe (Matrix a)
f Int
j = Maybe (Matrix a)
triangles3
where
wrows :: [Int]
wrows = forall a. Eq a => a -> [a] -> [Int]
elemIndices (Vector Int
ind3 forall a. Unbox a => Vector a -> Int -> a
! Int
j) [Int]
index3
triangles3 :: Maybe (Matrix a)
triangles3 = if forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Int]
wrows then forall a. Maybe a
Nothing else forall a. a -> Maybe a
Just forall a b. (a -> b) -> a -> b
$ forall a. Fractional a => Matrix a -> Matrix a
calPoints Matrix a
points3
where
wcols :: Vector Int
wcols = forall a. Unbox a => a -> Vector a -> Vector a
UV.cons Int
nedge ((Vector (Vector (Vector Int))
specialPos forall a. Vector a -> Int -> a
V.! Int
c) forall a. Vector a -> Int -> a
V.! Int
j)
ed :: Matrix Int
ed = forall a. Matrix a -> [Int] -> Vector Int -> Matrix a
subMatrix Matrix Int
edgesp1index [Int]
wrows Vector Int
wcols
col0ed :: [Int]
col0ed = forall a. Vector a -> [a]
V.toList forall a b. (a -> b) -> a -> b
$ forall a. Int -> Matrix a -> Vector a
getCol Int
1 Matrix Int
ed
col0edrep :: [Int]
col0edrep = forall (t :: * -> *) a. Foldable t => t a -> [a]
F.toList forall a b. (a -> b) -> a -> b
$ forall a. [a] -> Int -> Seq a
replicateEach' [Int]
col0ed (forall a. Unbox a => Vector a -> Int
UV.length Vector Int
wcols forall a. Num a => a -> a -> a
- Int
1)
edge1 :: [Int]
edge1 = forall a. Matrix a -> [a]
matrix2listMinusFirstColumn Matrix Int
ed
epoints' :: [Vector Int]
epoints' = forall a b. (a -> b) -> [a] -> [b]
map (\Int
k -> Vector (Vector Int)
edgePoints forall a. Vector a -> Int -> a
V.! (Int
k forall a. Num a => a -> a -> a
- Int
1)) [Int]
edge1
x1' :: [Int]
x1' = forall a b. (a -> b) -> [a] -> [b]
map (forall a. Unbox a => Vector a -> Int -> a
! Int
1) [Vector Int]
epoints'
x2' :: [Int]
x2' = forall a b. (a -> b) -> [a] -> [b]
map (forall a. Unbox a => Vector a -> Int -> a
! Int
2) [Vector Int]
epoints'
points3 :: Matrix a
points3 = forall a.
(Unbox a, RealFrac a) =>
Matrix Int -> Vector a -> [Int] -> [Int] -> [Int] -> Matrix a
getPoints Matrix Int
cubeco3 Vector a
values3 [Int]
col0edrep [Int]
x1' [Int]
x2'
triangles2 :: Maybe (Matrix a)
triangles2 = if forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Matrix a]
setOfTriangles
then forall a. Maybe a
Nothing
else forall a. a -> Maybe a
Just forall a b. (a -> b) -> a -> b
$ forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldl1 forall a. Matrix a -> Matrix a -> Matrix a
(<->) [Matrix a]
setOfTriangles