{-# LANGUAGE ForeignFunctionInterface #-}
{-# LANGUAGE FlexibleContexts         #-}
module MarchingCubes.MarchingCubes
  ( Voxel
  , XYZ
  , marchingCubes
  , makeVoxel
  ) where
import           Control.Monad                 (when)
import           Data.List                     (transpose)
import           Foreign.C.Types
import           Foreign.Marshal.Alloc         (free, mallocBytes)
import           Foreign.Marshal.Array         (peekArray, pokeArray)
import           Foreign.Ptr                   (Ptr)
import           Foreign.Storable              (peek, sizeOf)
import           Data.Matrix                   ( Matrix
                                               , mapCol
                                               , fromLists )

type Bounds a = ((a, a), (a, a), (a, a))
type Dims = (Int, Int, Int)
type Voxel a = ([a], (Bounds a, Dims))
type XYZ a = (a, a, a)

-- | Make the voxel. 
makeVoxel
  :: RealFloat 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 =>
(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) =
  ([a]
values, (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_ ]

voxelMax :: RealFloat a => Voxel a -> a
voxelMax :: forall a. RealFloat a => Voxel a -> a
voxelMax ([a]
values, (Bounds a, Dims)
_) = 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''

foreign import ccall unsafe "computeContour3d" c_computeContour3d
    :: Ptr CDouble
    -> CUInt
    -> CUInt
    -> CUInt
    -> CDouble
    -> CDouble
    -> Ptr CSize
    -> IO (Ptr (Ptr CDouble))

computeContour3d :: 
  Voxel Double -> Double -> IO (Ptr (Ptr CDouble), Int)
computeContour3d :: Voxel Double -> Double -> IO (Ptr (Ptr CDouble), Int)
computeContour3d Voxel Double
voxel Double
level = do 
  let voxmax :: Double
voxmax = forall a. RealFloat a => Voxel a -> a
voxelMax Voxel Double
voxel
      (Int
nx, Int
ny, Int
nz) = forall a b. (a, b) -> b
snd forall a b. (a -> b) -> a -> b
$ forall a b. (a, b) -> b
snd Voxel Double
voxel 
      voxel' :: [CDouble]
voxel' = forall a b. (a -> b) -> [a] -> [b]
map forall a b. (Real a, Fractional b) => a -> b
realToFrac (forall a b. (a, b) -> a
fst Voxel Double
voxel)
  Ptr CSize
nrowsPtr <- forall a. Int -> IO (Ptr a)
mallocBytes (forall a. Storable a => a -> Int
sizeOf (forall a. HasCallStack => a
undefined :: CSize))
  String -> IO ()
putStrLn String
"Allocating voxel"
  Ptr CDouble
voxelPtr <- forall a. Int -> IO (Ptr a)
mallocBytes (Int
nxforall a. Num a => a -> a -> a
*Int
nyforall a. Num a => a -> a -> a
*Int
nz forall a. Num a => a -> a -> a
* forall a. Storable a => a -> Int
sizeOf (forall a. HasCallStack => a
undefined :: CDouble))
  String -> IO ()
putStrLn String
"Poking voxel..."
  forall a. Storable a => Ptr a -> [a] -> IO ()
pokeArray Ptr CDouble
voxelPtr [CDouble]
voxel'
  String -> IO ()
putStrLn String
"Run C code..."
  Ptr (Ptr CDouble)
result <- Ptr CDouble
-> CUInt
-> CUInt
-> CUInt
-> CDouble
-> CDouble
-> Ptr CSize
-> IO (Ptr (Ptr CDouble))
c_computeContour3d Ptr CDouble
voxelPtr
            (forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
nx) (forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
ny) (forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
nz)
            (forall a b. (Real a, Fractional b) => a -> b
realToFrac Double
voxmax) (forall a b. (Real a, Fractional b) => a -> b
realToFrac Double
level) Ptr CSize
nrowsPtr
  CSize
nrows <- forall a. Storable a => Ptr a -> IO a
peek Ptr CSize
nrowsPtr
  forall a. Ptr a -> IO ()
free Ptr CSize
nrowsPtr
  forall a. Ptr a -> IO ()
free Ptr CDouble
voxelPtr
  forall (m :: * -> *) a. Monad m => a -> m a
return (Ptr (Ptr CDouble)
result, forall a b. (Integral a, Num b) => a -> b
fromIntegral CSize
nrows)

marchingCubes :: Voxel Double -> Double -> Bool -> IO (Matrix Double)
marchingCubes :: Voxel Double -> Double -> Bool -> IO (Matrix Double)
marchingCubes Voxel Double
voxel Double
level Bool
summary = do
  let (bds :: Bounds Double
bds@((Double, Double)
xbds, (Double, Double)
ybds, (Double, Double)
zbds), dims :: Dims
dims@(Int
nx, Int
ny, Int
nz)) = forall a b. (a, b) -> b
snd Voxel Double
voxel
  (Ptr (Ptr CDouble)
ppCDouble, Int
nrows) <- Voxel Double -> Double -> IO (Ptr (Ptr CDouble), Int)
computeContour3d Voxel Double
voxel Double
level
  [[CDouble]]
cpoints <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM (forall a. Storable a => Int -> Ptr a -> IO [a]
peekArray Int
3) forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< forall a. Storable a => Int -> Ptr a -> IO [a]
peekArray Int
nrows Ptr (Ptr CDouble)
ppCDouble
  let toDouble :: CDouble -> Double
      toDouble :: CDouble -> Double
toDouble = forall a b. (Real a, Fractional b) => a -> b
realToFrac
      points :: [[Double]]
points = forall a b. (a -> b) -> [a] -> [b]
map (forall a b. (a -> b) -> [a] -> [b]
map CDouble -> Double
toDouble) [[CDouble]]
cpoints
  forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when Bool
summary forall a b. (a -> b) -> a -> b
$ do
    let tpoints :: [[Double]]
tpoints = forall a. [[a]] -> [[a]]
transpose [[Double]]
points
        x :: [Double]
x = [[Double]]
tpoints forall a. [a] -> Int -> a
!! Int
0
        y :: [Double]
y = [[Double]]
tpoints forall a. [a] -> Int -> a
!! Int
1
        z :: [Double]
z = [[Double]]
tpoints forall a. [a] -> Int -> a
!! Int
2
        xm :: Double
xm = forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum [Double]
x forall a. Num a => a -> a -> a
- Double
1
        xM :: Double
xM = forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Double]
x forall a. Num a => a -> a -> a
- Double
1
        ym :: Double
ym = forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum [Double]
y forall a. Num a => a -> a -> a
- Double
1
        yM :: Double
yM = forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Double]
y forall a. Num a => a -> a -> a
- Double
1
        zm :: Double
zm = forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum [Double]
z forall a. Num a => a -> a -> a
- Double
1
        zM :: Double
zM = forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Double]
z forall a. Num a => a -> a -> a
- Double
1
        xmin :: Double
xmin = forall a. Fractional a => (a, a) -> Int -> a -> a
rescale (Double, Double)
xbds Int
nx Double
xm
        xmax :: Double
xmax = forall a. Fractional a => (a, a) -> Int -> a -> a
rescale (Double, Double)
xbds Int
nx Double
xM
        ymin :: Double
ymin = forall a. Fractional a => (a, a) -> Int -> a -> a
rescale (Double, Double)
ybds Int
ny Double
ym
        ymax :: Double
ymax = forall a. Fractional a => (a, a) -> Int -> a -> a
rescale (Double, Double)
ybds Int
ny Double
yM
        zmin :: Double
zmin = forall a. Fractional a => (a, a) -> Int -> a -> a
rescale (Double, Double)
zbds Int
nz Double
zm
        zmax :: Double
zmax = forall a. Fractional a => (a, a) -> Int -> a -> a
rescale (Double, Double)
zbds Int
nz Double
zM
    String -> IO ()
putStrLn String
"Prebounds:"
    forall a. Show a => a -> IO ()
print ((Double
xm, Double
ym, Double
zm), (Double
xM, Double
yM, Double
zM))
    String -> IO ()
putStrLn String
"Bounds:"
    forall a. Show a => a -> IO ()
print ((Double
xmin, Double
ymin, Double
zmin), (Double
xmax, Double
ymax, Double
zmax))
  forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ forall a. Fractional a => Matrix a -> Bounds a -> Dims -> Matrix a
rescaleMatrix (forall a. [[a]] -> Matrix a
fromLists [[Double]]
points) Bounds Double
bds Dims
dims