{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE TupleSections    #-}
module MarchingCubes.Mesh
  ( makeMesh
  , Mesh (..)
  ) where
import           Data.Bifunctor                 ( bimap )
import           Data.Foldable                  ( foldl' )
import           Data.IntMap.Strict             ( IntMap
                                                , elems
                                                , unionWith
                                                )
import qualified Data.IntMap.Strict            as IM
import           Data.List.Split                ( chunksOf )
import           Data.Map.Strict                ( insertLookupWithKey )
import qualified Data.Map.Strict               as M
import           Data.Matrix                    ( Matrix
                                                , nrows
                                                , toLists
                                                )
import           Data.Tuple.Extra               ( (***)
                                                , curry3
                                                , second
                                                , thd3
                                                )
import           Data.Vector.Unboxed            ( (!)
                                                , Unbox
                                                , Vector
                                                , fromList
                                                )
import           Linear                         ( V3(..)
                                                , (^+^)
                                                , (^-^)
                                                , cross
                                                , signorm
                                                )
import           MarchingCubes.MarchingCubes    ( Voxel
                                                , marchingCubes
                                                )


data Mesh a = Mesh 
  { 
    forall a. Mesh a -> Vector (a, a, a)
_vertices :: Vector (a, a, a)
  , forall a. Mesh a -> [(Int, Int, Int)]
_faces :: [(Int, Int, Int)] 
  , forall a. Mesh a -> Vector (a, a, a)
_normals :: Vector (a, a, a)
  } 
  deriving Int -> Mesh a -> ShowS
forall a. (Show a, Unbox a) => Int -> Mesh a -> ShowS
forall a. (Show a, Unbox a) => [Mesh a] -> ShowS
forall a. (Show a, Unbox a) => Mesh a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Mesh a] -> ShowS
$cshowList :: forall a. (Show a, Unbox a) => [Mesh a] -> ShowS
show :: Mesh a -> String
$cshow :: forall a. (Show a, Unbox a) => Mesh a -> String
showsPrec :: Int -> Mesh a -> ShowS
$cshowsPrec :: forall a. (Show a, Unbox a) => Int -> Mesh a -> ShowS
Show

matrix2listOfTriples :: Matrix a -> [(a, a, a)]
matrix2listOfTriples :: forall a. Matrix a -> [(a, a, a)]
matrix2listOfTriples Matrix a
mtrx =
  forall a b. (a -> b) -> [a] -> [b]
map (\[a]
row -> ([a]
row forall a. [a] -> Int -> a
!! Int
0, [a]
row forall a. [a] -> Int -> a
!! Int
1, [a]
row forall a. [a] -> Int -> a
!! Int
2)) (forall a. Matrix a -> [[a]]
toLists Matrix a
mtrx)

uniqueWithIndices :: (Unbox a, Ord a) => [a] -> (Vector a, Vector Int)
uniqueWithIndices :: forall a. (Unbox a, Ord a) => [a] -> (Vector a, Vector Int)
uniqueWithIndices [a]
list = (forall a. Unbox a => [a] -> Vector a
fromList forall a a' b b'. (a -> a') -> (b -> b') -> (a, b) -> (a', b')
*** forall a. Unbox a => [a] -> Vector a
fromList) forall a b. (a -> b) -> a -> b
$ forall {a} {t}. (Ord a, Num t) => t -> Map a t -> [a] -> ([a], [t])
go Int
0 forall k a. Map k a
M.empty [a]
list
 where
  go :: t -> Map a t -> [a] -> ([a], [t])
go t
_ Map a t
_ []       = ([], [])
  go t
i Map a t
m (a
x : [a]
xs) = case forall k a.
Ord k =>
(k -> a -> a -> a) -> k -> a -> Map k a -> (Maybe a, Map k a)
insertLookupWithKey (forall a b c d. ((a, b, c) -> d) -> a -> b -> c -> d
curry3 forall a b c. (a, b, c) -> c
thd3) a
x t
i Map a t
m of
    (Maybe t
Nothing, Map a t
m') -> forall (p :: * -> * -> *) a b c d.
Bifunctor p =>
(a -> b) -> (c -> d) -> p a c -> p b d
bimap (a
x forall a. a -> [a] -> [a]
:) (t
i forall a. a -> [a] -> [a]
:) (t -> Map a t -> [a] -> ([a], [t])
go (t
i forall a. Num a => a -> a -> a
+ t
1) Map a t
m' [a]
xs)
    (Just t
j , Map a t
m') -> forall b b' a. (b -> b') -> (a, b) -> (a, b')
second (t
j forall a. a -> [a] -> [a]
:) (t -> Map a t -> [a] -> ([a], [t])
go t
i Map a t
m' [a]
xs)

undupMesh
  :: (Unbox a, Ord a) => ([(a, a, a)], [[Int]]) -> (Vector (a, a, a), [[Int]])
undupMesh :: forall a.
(Unbox a, Ord a) =>
([(a, a, a)], [[Int]]) -> (Vector (a, a, a), [[Int]])
undupMesh ([(a, a, a)]
vertices, [[Int]]
faces) = (Vector (a, a, a)
vertices', [[Int]]
faces')
 where
  (Vector (a, a, a)
vertices', Vector Int
idx) = forall a. (Unbox a, Ord a) => [a] -> (Vector a, Vector Int)
uniqueWithIndices [(a, a, a)]
vertices
  faces' :: [[Int]]
faces'           = forall a b. (a -> b) -> [a] -> [b]
map (forall a b. (a -> b) -> [a] -> [b]
map (Vector Int
idx forall a. Unbox a => Vector a -> Int -> a
!)) [[Int]]
faces

normal :: Floating a => (a, a, a) -> (a, a, a) -> (a, a, a) -> V3 a
normal :: forall a. Floating a => (a, a, a) -> (a, a, a) -> (a, a, a) -> V3 a
normal (a, a, a)
v1 (a, a, a)
v2 (a, a, a)
v3 = forall (f :: * -> *) a. (Metric f, Floating a) => f a -> f a
signorm forall a b. (a -> b) -> a -> b
$ forall a. Num a => V3 a -> V3 a -> V3 a
cross (V3 a
v3' forall (f :: * -> *) a. (Additive f, Num a) => f a -> f a -> f a
^-^ V3 a
v1') (V3 a
v2' forall (f :: * -> *) a. (Additive f, Num a) => f a -> f a -> f a
^-^ V3 a
v1') -- (NaN,NaN,NaN) if degenerate face
 where
  toV3 :: (a, a, a) -> V3 a
toV3 (a
x, a
y, a
z) = forall a. a -> a -> a -> V3 a
V3 a
x a
y a
z
  v1' :: V3 a
v1' = forall {a}. (a, a, a) -> V3 a
toV3 (a, a, a)
v1
  v2' :: V3 a
v2' = forall {a}. (a, a, a) -> V3 a
toV3 (a, a, a)
v2
  v3' :: V3 a
v3' = forall {a}. (a, a, a) -> V3 a
toV3 (a, a, a)
v3

faceNormals
  :: (Floating a, Unbox a) => Vector (a, a, a) -> [Int] -> IntMap (V3 a)
faceNormals :: forall a.
(Floating a, Unbox a) =>
Vector (a, a, a) -> [Int] -> IntMap (V3 a)
faceNormals Vector (a, a, a)
vs [Int]
face = forall a. [(Int, a)] -> IntMap a
IM.fromList [(Int, V3 a)]
nrmls
 where
  nrml :: V3 a
nrml  = forall a. Floating a => (a, a, a) -> (a, a, a) -> (a, a, a) -> V3 a
normal (Vector (a, a, a)
vs forall a. Unbox a => Vector a -> Int -> a
! ([Int]
face forall a. [a] -> Int -> a
!! Int
0)) (Vector (a, a, a)
vs forall a. Unbox a => Vector a -> Int -> a
! ([Int]
face forall a. [a] -> Int -> a
!! Int
1)) (Vector (a, a, a)
vs forall a. Unbox a => Vector a -> Int -> a
! ([Int]
face forall a. [a] -> Int -> a
!! Int
2))
  nrmls :: [(Int, V3 a)]
nrmls = forall a b. (a -> b) -> [a] -> [b]
map (, V3 a
nrml) [Int]
face

normalsV3
  :: (Floating a, Unbox a) => (Vector (a, a, a), [[Int]]) -> IntMap (V3 a)
normalsV3 :: forall a.
(Floating a, Unbox a) =>
(Vector (a, a, a), [[Int]]) -> IntMap (V3 a)
normalsV3 (Vector (a, a, a)
vs, [[Int]]
faces) = forall a b. (a -> b) -> IntMap a -> IntMap b
IM.map
  forall (f :: * -> *) a. (Metric f, Floating a) => f a -> f a
signorm
  (forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' (forall a. (a -> a -> a) -> IntMap a -> IntMap a -> IntMap a
unionWith forall (f :: * -> *) a. (Additive f, Num a) => f a -> f a -> f a
(^+^)) forall a. IntMap a
IM.empty (forall a b. (a -> b) -> [a] -> [b]
map (forall a.
(Floating a, Unbox a) =>
Vector (a, a, a) -> [Int] -> IntMap (V3 a)
faceNormals Vector (a, a, a)
vs) [[Int]]
faces))

normals
  :: (Floating a, Unbox a) => (Vector (a, a, a), [[Int]]) -> Vector (a, a, a)
normals :: forall a.
(Floating a, Unbox a) =>
(Vector (a, a, a), [[Int]]) -> Vector (a, a, a)
normals (Vector (a, a, a), [[Int]])
mesh = forall a. Unbox a => [a] -> Vector a
fromList forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map forall {c}. V3 c -> (c, c, c)
fromV3 (forall a. IntMap a -> [a]
elems forall a b. (a -> b) -> a -> b
$ forall a.
(Floating a, Unbox a) =>
(Vector (a, a, a), [[Int]]) -> IntMap (V3 a)
normalsV3 (Vector (a, a, a), [[Int]])
mesh)
  where fromV3 :: V3 c -> (c, c, c)
fromV3 (V3 c
x c
y c
z) = (c
x, c
y, c
z)

degenerateFace :: (Unbox a, Eq a) => Vector (a, a, a) -> [Int] -> Bool
degenerateFace :: forall a. (Unbox a, Eq a) => Vector (a, a, a) -> [Int] -> Bool
degenerateFace Vector (a, a, a)
vertices [Int]
face = (a, a, a)
v1 forall a. Eq a => a -> a -> Bool
== (a, a, a)
v2 Bool -> Bool -> Bool
|| (a, a, a)
v1 forall a. Eq a => a -> a -> Bool
== (a, a, a)
v3 Bool -> Bool -> Bool
|| (a, a, a)
v2 forall a. Eq a => a -> a -> Bool
== (a, a, a)
v3
 where
  v1 :: (a, a, a)
v1 = Vector (a, a, a)
vertices forall a. Unbox a => Vector a -> Int -> a
! ([Int]
face forall a. [a] -> Int -> a
!! Int
0)
  v2 :: (a, a, a)
v2 = Vector (a, a, a)
vertices forall a. Unbox a => Vector a -> Int -> a
! ([Int]
face forall a. [a] -> Int -> a
!! Int
1)
  v3 :: (a, a, a)
v3 = Vector (a, a, a)
vertices forall a. Unbox a => Vector a -> Int -> a
! ([Int]
face forall a. [a] -> Int -> a
!! Int
2)

-- | Make mesh of isosurface.
makeMesh :: 
     Voxel Double -- ^ voxel obtained with `makeVoxel`
  -> Double       -- ^ isovalue
  -> IO (Mesh Double)
makeMesh :: Voxel Double -> Double -> IO (Mesh Double)
makeMesh Voxel Double
voxel Double
level = do
  Matrix Double
mtrx <- Voxel Double -> Double -> Bool -> IO (Matrix Double)
marchingCubes Voxel Double
voxel Double
level Bool
True
  let vertices :: [(Double, Double, Double)]
vertices            = forall a. Matrix a -> [(a, a, a)]
matrix2listOfTriples Matrix Double
mtrx
      faces :: [[Int]]
faces               = forall e. Int -> [e] -> [[e]]
chunksOf Int
3 [Int
0 .. forall a. Matrix a -> Int
nrows Matrix Double
mtrx forall a. Num a => a -> a -> a
- Int
1]
      (Vector (Double, Double, Double)
vertices', [[Int]]
faces') = forall a.
(Unbox a, Ord a) =>
([(a, a, a)], [[Int]]) -> (Vector (a, a, a), [[Int]])
undupMesh ([(Double, Double, Double)]
vertices, [[Int]]
faces)
      faces'' :: [[Int]]
faces''             = forall a. (a -> Bool) -> [a] -> [a]
filter (Bool -> Bool
not forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. (Unbox a, Eq a) => Vector (a, a, a) -> [Int] -> Bool
degenerateFace Vector (Double, Double, Double)
vertices') [[Int]]
faces'
      toTriplet :: [c] -> (c, c, c)
toTriplet [c]
ikj       = ([c]
ikj forall a. [a] -> Int -> a
!! Int
0, [c]
ikj forall a. [a] -> Int -> a
!! Int
2, [c]
ikj forall a. [a] -> Int -> a
!! Int
1)
      faces''' :: [(Int, Int, Int)]
faces'''            = forall a b. (a -> b) -> [a] -> [b]
map forall {c}. [c] -> (c, c, c)
toTriplet [[Int]]
faces''
      mesh :: (Vector (Double, Double, Double), [[Int]])
mesh                = (Vector (Double, Double, Double)
vertices', [[Int]]
faces'')
  forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ Mesh 
    {
      _vertices :: Vector (Double, Double, Double)
_vertices = Vector (Double, Double, Double)
vertices', _faces :: [(Int, Int, Int)]
_faces = [(Int, Int, Int)]
faces''', _normals :: Vector (Double, Double, Double)
_normals = forall a.
(Floating a, Unbox a) =>
(Vector (a, a, a), [[Int]]) -> Vector (a, a, a)
normals (Vector (Double, Double, Double), [[Int]])
mesh
    }