module Geometry.Delaunay.Delaunay
  ( delaunay
  , vertexNeighborFacets
  , sandwichedFacet
  , facetOf
  , facetFamilies
  , facetCenters
  , facetOf'
  , facetFamilies'
  , facetCenters'
  , getDelaunayTiles
  ) 
  where
import           Control.Monad               ( unless, when )
import           Data.IntMap.Strict          ( IntMap )
import qualified Data.IntMap.Strict          as IM
import qualified Data.IntSet                 as IS
import           Data.List.Unique            ( allUnique )
import           Data.Maybe                  ( fromMaybe )
import           Geometry.Delaunay.CDelaunay ( c_tessellation
                                             , cTessellationToTessellation 
                                             )
import           Geometry.Delaunay.Types     ( Tessellation(_tilefacets, _sites, _tiles)
                                             , Tile (..)
                                             , Simplex(_vertices')
                                             , TileFacet(_facetOf)
                                             , Site(_neighfacetsIds) 
                                             )
import           Foreign.C.Types             ( CDouble, CUInt )
import           Foreign.Marshal.Alloc       ( free, mallocBytes )
import           Foreign.Marshal.Array       ( pokeArray )
import           Foreign.Storable            ( peek, sizeOf )
import           Geometry.Qhull.Types        ( HasCenter(_center) 
                                             , HasFamily(_family)
                                             , Family
                                             , Index 
                                             )

delaunay :: [[Double]]      -- ^ sites (vertex coordinates)

         -> Bool            -- ^ whether to add a point at infinity

         -> Bool            -- ^ whether to include degenerate tiles

         -> Maybe Double    -- ^ volume threshold

         -> IO Tessellation -- ^ Delaunay tessellation

delaunay :: [[Double]] -> Bool -> Bool -> Maybe Double -> IO Tessellation
delaunay [[Double]]
sites Bool
atinfinity Bool
degenerate Maybe Double
vthreshold = do
  let n :: Int
n     = forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
sites
      dim :: Int
dim   = forall (t :: * -> *) a. Foldable t => t a -> Int
length (forall a. [a] -> a
head [[Double]]
sites)
  forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
dim forall a. Ord a => a -> a -> Bool
< Int
2) forall a b. (a -> b) -> a -> b
$
    forall a. HasCallStack => [Char] -> a
error [Char]
"dimension must be at least 2"
  forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
n forall a. Ord a => a -> a -> Bool
<= Int
dimforall a. Num a => a -> a -> a
+Int
1) forall a b. (a -> b) -> a -> b
$
    forall a. HasCallStack => [Char] -> a
error [Char]
"insufficient number of points"
  forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
unless (forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (forall a. Eq a => a -> a -> Bool
== Int
dim) (forall a b. (a -> b) -> [a] -> [b]
map forall (t :: * -> *) a. Foldable t => t a -> Int
length (forall a. [a] -> [a]
tail [[Double]]
sites))) forall a b. (a -> b) -> a -> b
$
    forall a. HasCallStack => [Char] -> a
error [Char]
"the points must have the same dimension"
  forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
unless (forall a. Ord a => [a] -> Bool
allUnique [[Double]]
sites) forall a b. (a -> b) -> a -> b
$
    forall a. HasCallStack => [Char] -> a
error [Char]
"some points are duplicated"
  let vthreshold' :: Double
vthreshold' = forall a. a -> Maybe a -> a
fromMaybe Double
0 Maybe Double
vthreshold 
  Ptr CDouble
sitesPtr <- forall a. Int -> IO (Ptr a)
mallocBytes (Int
n forall a. Num a => a -> a -> a
* Int
dim forall a. Num a => a -> a -> a
* forall a. Storable a => a -> Int
sizeOf (forall a. HasCallStack => a
undefined :: CDouble))
  forall a. Storable a => Ptr a -> [a] -> IO ()
pokeArray Ptr CDouble
sitesPtr (forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap (forall a b. (a -> b) -> [a] -> [b]
map forall a b. (Real a, Fractional b) => a -> b
realToFrac) [[Double]]
sites)
  Ptr CUInt
exitcodePtr <- forall a. Int -> IO (Ptr a)
mallocBytes (forall a. Storable a => a -> Int
sizeOf (forall a. HasCallStack => a
undefined :: CUInt))
  Ptr CTessellation
resultPtr <- Ptr CDouble
-> CUInt
-> CUInt
-> CUInt
-> CUInt
-> CDouble
-> Ptr CUInt
-> IO (Ptr CTessellation)
c_tessellation Ptr CDouble
sitesPtr
               (forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
dim) (forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n)
               (forall a b. (Integral a, Num b) => a -> b
fromIntegral forall a b. (a -> b) -> a -> b
$ forall a. Enum a => a -> Int
fromEnum Bool
atinfinity)
               (forall a b. (Integral a, Num b) => a -> b
fromIntegral forall a b. (a -> b) -> a -> b
$ forall a. Enum a => a -> Int
fromEnum Bool
degenerate)
               (forall a b. (Real a, Fractional b) => a -> b
realToFrac Double
vthreshold') Ptr CUInt
exitcodePtr
  CUInt
exitcode <- forall a. Storable a => Ptr a -> IO a
peek Ptr CUInt
exitcodePtr
  forall a. Ptr a -> IO ()
free Ptr CUInt
exitcodePtr
  forall a. Ptr a -> IO ()
free Ptr CDouble
sitesPtr
  if CUInt
exitcode forall a. Eq a => a -> a -> Bool
/= CUInt
0
    then do
      forall a. Ptr a -> IO ()
free Ptr CTessellation
resultPtr
      forall a. HasCallStack => [Char] -> a
error forall a b. (a -> b) -> a -> b
$ [Char]
"qhull returned an error (code " forall a. [a] -> [a] -> [a]
++ forall a. Show a => a -> [Char]
show CUInt
exitcode forall a. [a] -> [a] -> [a]
++ [Char]
")"
    else do
      CTessellation
result <- forall a. Storable a => Ptr a -> IO a
peek Ptr CTessellation
resultPtr
      Tessellation
out <- [[Double]] -> CTessellation -> IO Tessellation
cTessellationToTessellation [[Double]]
sites CTessellation
result
      forall a. Ptr a -> IO ()
free Ptr CTessellation
resultPtr
      forall (m :: * -> *) a. Monad m => a -> m a
return Tessellation
out

-- | tile facets a vertex belongs to, vertex given by its index;

-- the output is the empty map if the index is not valid

vertexNeighborFacets :: Tessellation -> Index -> IntMap TileFacet
vertexNeighborFacets :: Tessellation -> Int -> IntMap TileFacet
vertexNeighborFacets Tessellation
tess Int
i = forall a. IntMap a -> IntSet -> IntMap a
IM.restrictKeys (Tessellation -> IntMap TileFacet
_tilefacets Tessellation
tess) IntSet
ids
  where
    ids :: IntSet
ids = forall b a. b -> (a -> b) -> Maybe a -> b
maybe IntSet
IS.empty Site -> IntSet
_neighfacetsIds (forall a. Int -> IntMap a -> Maybe a
IM.lookup Int
i (Tessellation -> IndexMap Site
_sites Tessellation
tess))

-- | whether a tile facet is sandwiched between two tiles

sandwichedFacet :: TileFacet -> Bool
sandwichedFacet :: TileFacet -> Bool
sandwichedFacet TileFacet
tilefacet = IntSet -> Int
IS.size (TileFacet -> IntSet
_facetOf TileFacet
tilefacet) forall a. Eq a => a -> a -> Bool
== Int
2

-- | the tiles a facet belongs to

facetOf :: Tessellation -> TileFacet -> IntMap Tile
facetOf :: Tessellation -> TileFacet -> IntMap Tile
facetOf Tessellation
tess TileFacet
tilefacet = forall a. IntMap a -> IntSet -> IntMap a
IM.restrictKeys (Tessellation -> IntMap Tile
_tiles Tessellation
tess) (TileFacet -> IntSet
_facetOf TileFacet
tilefacet)

-- | the families of the tiles a facet belongs to

facetFamilies :: Tessellation -> TileFacet -> IntMap Family
facetFamilies :: Tessellation -> TileFacet -> IntMap Family
facetFamilies Tessellation
tess TileFacet
tilefacet = forall a b. (a -> b) -> IntMap a -> IntMap b
IM.map forall m. HasFamily m => m -> Family
_family (Tessellation -> TileFacet -> IntMap Tile
facetOf Tessellation
tess TileFacet
tilefacet)

-- | the circumcenters of the tiles a facet belongs to

facetCenters :: Tessellation -> TileFacet -> IntMap [Double]
facetCenters :: Tessellation -> TileFacet -> IntMap [Double]
facetCenters Tessellation
tess TileFacet
tilefacet =
  forall a b. (a -> b) -> IntMap a -> IntMap b
IM.map forall m. HasCenter m => m -> [Double]
_center (Tessellation -> TileFacet -> IntMap Tile
facetOf Tessellation
tess TileFacet
tilefacet)

funofFacetToFunofInt :: (Tessellation -> TileFacet -> IntMap a)
                     -> (Tessellation -> Int -> IntMap a)
funofFacetToFunofInt :: forall a.
(Tessellation -> TileFacet -> IntMap a)
-> Tessellation -> Int -> IntMap a
funofFacetToFunofInt Tessellation -> TileFacet -> IntMap a
f Tessellation
tess Int
i =
  forall b a. b -> (a -> b) -> Maybe a -> b
maybe forall a. IntMap a
IM.empty (Tessellation -> TileFacet -> IntMap a
f Tessellation
tess) (forall a. Int -> IntMap a -> Maybe a
IM.lookup Int
i (Tessellation -> IntMap TileFacet
_tilefacets Tessellation
tess))

-- | the tiles a facet belongs to, facet given by its id

facetOf' :: Tessellation -> Int -> IntMap Tile
facetOf' :: Tessellation -> Int -> IntMap Tile
facetOf' = forall a.
(Tessellation -> TileFacet -> IntMap a)
-> Tessellation -> Int -> IntMap a
funofFacetToFunofInt Tessellation -> TileFacet -> IntMap Tile
facetOf

-- | the families of the tiles a facet belongs to, facet given by its id

facetFamilies' :: Tessellation -> Int -> IntMap Family
facetFamilies' :: Tessellation -> Int -> IntMap Family
facetFamilies' = forall a.
(Tessellation -> TileFacet -> IntMap a)
-> Tessellation -> Int -> IntMap a
funofFacetToFunofInt Tessellation -> TileFacet -> IntMap Family
facetFamilies

-- | the circumcenters of the tiles a facet belongs to, facet given by its id

facetCenters' :: Tessellation -> Int -> IntMap [Double]
facetCenters' :: Tessellation -> Int -> IntMap [Double]
facetCenters' = forall a.
(Tessellation -> TileFacet -> IntMap a)
-> Tessellation -> Int -> IntMap a
funofFacetToFunofInt Tessellation -> TileFacet -> IntMap [Double]
facetCenters

-- | list of the maps of vertices for all tiles

getDelaunayTiles :: Tessellation -> [IntMap [Double]]
getDelaunayTiles :: Tessellation -> [IntMap [Double]]
getDelaunayTiles Tessellation
tess = forall a. IntMap a -> [a]
IM.elems forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> IntMap a -> IntMap b
IM.map (Simplex -> IntMap [Double]
_vertices' forall b c a. (b -> c) -> (a -> b) -> a -> c
. Tile -> Simplex
_simplex) (Tessellation -> IntMap Tile
_tiles Tessellation
tess)