{-# 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)

-- | Make the voxel. 
makeVoxel
  :: (RealFloat a, IArray UArray a)
  => (XYZ a -> a) -- ^ the function defining the isosurface
  -> Bounds a     -- ^ bounds of the grid
  -> Dims         -- ^ numbers of subdivisions of the grid
  -> 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
  -- ~~## special cases ##~~
  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