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]]
-> Bool
-> Bool
-> Maybe Double
-> IO 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
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))
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
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)
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)
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))
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
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
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
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)