Skip to content

Commit

Permalink
export lists
Browse files Browse the repository at this point in the history
  • Loading branch information
stla committed Nov 15, 2023
1 parent 71c0b7f commit 92014e2
Show file tree
Hide file tree
Showing 6 changed files with 65 additions and 40 deletions.
2 changes: 1 addition & 1 deletion README.md
Expand Up @@ -161,7 +161,7 @@ Finally, the output of `delaunay` has the `_edges'` field, providing the
edges:

```haskell
_edges' d
> _edges' d
fromList
[ ( Pair 0 1 , ( [ -5.0 , -5.0 , 16.0 ] , [ -5.0 , 8.0 , 3.0 ] ) )
, ( Pair 0 3 , ( [ -5.0 , -5.0 , 16.0 ] , [ 4.0 , -5.0 , 7.0 ] ) )
Expand Down
37 changes: 14 additions & 23 deletions src/Geometry/Delaunay/CDelaunay.hs
@@ -1,8 +1,10 @@
-- {-# LINE 1 "delaunay.hsc" #-}
{-# LANGUAGE ForeignFunctionInterface #-}
module Geometry.Delaunay.CDelaunay
( cTesselationToTesselation
, c_tesselation )
(
cTessellationToTessellation
, c_tessellation
)
where
import Control.Monad ((<$!>))
import qualified Data.HashMap.Strict.InsOrd as H
Expand Down Expand Up @@ -95,7 +97,7 @@ cSiteToSite sites csite = do
where
filterAscList :: Int -> [Int] -> [Int]
filterAscList n list =
let i = findIndex (>n) list in
let i = findIndex (> n) list in
if isJust i
then drop (fromJust i) list
else []
Expand Down Expand Up @@ -144,11 +146,9 @@ cSimplexToSimplex sites simplexdim csimplex = do
dim = length (head sites)
sitesids <- (<$!>) (map fromIntegral)
(peekArray simplexdim (__sitesids csimplex))
-- putStrLn "cSimplexToSimplex - peek sitesids"
let points = fromAscList
(zip sitesids (map ((!!) sites) sitesids))
center <- (<$!>) (map cdbl2dbl) (peekArray dim (__center csimplex))
-- putStrLn "cSimplexToSimplex - peek center"
return Simplex { _vertices' = points
, _circumcenter = center
, _circumradius = radius
Expand Down Expand Up @@ -216,7 +216,6 @@ cSubTiletoTileFacet points csubtile = do
offset = realToFrac $ __offset csubtile
simplex <- cSimplexToSimplex points dim subsimplex
normal <- (<$!>) (map realToFrac) (peekArray dim (__normal csubtile))
-- putStrLn "cSubTiletoTileFacet - peek normal"
return (id', TileFacet { _subsimplex = simplex
, _facetOf = IS.fromAscList ridgeOf
, _normal' = normal
Expand Down Expand Up @@ -292,14 +291,11 @@ cTileToTile points ctile = do
family = __family ctile
orient = __orientation ctile
dim = length (head points)
-- putStrLn $ "tile id: " ++ show id'
simplex <- cSimplexToSimplex points (dim+1) csimplex
neighbors <- (<$!>) (map fromIntegral)
(peekArray nneighbors (__neighbors ctile))
-- putStrLn "cTileToTile - peek neighbors"
ridgesids <- (<$!>) (map fromIntegral)
(peekArray nridges (__ridgesids ctile))
-- putStrLn "cTileToTile - peek ridges"
return (id', Tile { _simplex = simplex
, _neighborsIds = IS.fromAscList neighbors
, _facetsIds = IS.fromAscList ridgesids
Expand All @@ -308,15 +304,15 @@ cTileToTile points ctile = do
else Family (fromIntegral family)
, _toporiented = orient == 1 })

data CTesselation = CTesselation {
data CTessellation = CTessellation {
__sites :: Ptr CSite
, __tiles :: Ptr CTile
, __ntiles :: CUInt
, __subtiles :: Ptr CSubTile
, __nsubtiles :: CUInt
}

instance Storable CTesselation where
instance Storable CTessellation where
sizeOf __ = (40)
-- {-# LINE 246 "delaunay.hsc" #-}
alignment __ = 8
Expand All @@ -332,14 +328,14 @@ instance Storable CTesselation where
-- {-# LINE 252 "delaunay.hsc" #-}
nsubtiles' <- (\hsc_ptr -> peekByteOff hsc_ptr 32) ptr
-- {-# LINE 253 "delaunay.hsc" #-}
return CTesselation {
return CTessellation {
__sites = sites'
, __tiles = tiles'
, __ntiles = ntiles'
, __subtiles = subtiles'
, __nsubtiles = nsubtiles'
}
poke ptr (CTesselation r1 r2 r3 r4 r5)
poke ptr (CTessellation r1 r2 r3 r4 r5)
= do
(\hsc_ptr -> pokeByteOff hsc_ptr 0) ptr r1
-- {-# LINE 263 "delaunay.hsc" #-}
Expand All @@ -352,36 +348,31 @@ instance Storable CTesselation where
(\hsc_ptr -> pokeByteOff hsc_ptr 32) ptr r5
-- {-# LINE 267 "delaunay.hsc" #-}

foreign import ccall unsafe "tesselation" c_tesselation
foreign import ccall unsafe "tessellation" c_tessellation
:: Ptr CDouble -- sites
-> CUInt -- dim
-> CUInt -- nsites
-> CUInt -- 0/1, point at infinity
-> CUInt -- 0/1, include degenerate
-> CDouble -- volume threshold
-> Ptr CUInt -- exitcode
-> IO (Ptr CTesselation)
-> IO (Ptr CTessellation)

cTesselationToTesselation :: [[Double]] -> CTesselation -> IO Tesselation
cTesselationToTesselation vertices ctess = do
cTessellationToTessellation :: [[Double]] -> CTessellation -> IO Tessellation
cTessellationToTessellation vertices ctess = do
let ntiles = fromIntegral $ __ntiles ctess
nsubtiles = fromIntegral $ __nsubtiles ctess
nsites = length vertices
sites'' <- peekArray nsites (__sites ctess)
-- putStrLn "peek sites"
tiles'' <- peekArray ntiles (__tiles ctess)
-- putStrLn "peek tiles"
subtiles'' <- peekArray nsubtiles (__subtiles ctess)
-- putStrLn "peek ridges"
sites' <- mapM (cSiteToSite vertices) sites''
let sites = fromAscList (map (fst3 &&& snd3) sites')
edgesIndices = concatMap thd3 sites'
edges = map (toPair &&& both (_point . ((!) sites))) edgesIndices
tiles' <- mapM (cTileToTile vertices) tiles''
-- putStrLn "mapped cTileToTile"
subtiles' <- mapM (cSubTiletoTileFacet vertices) subtiles''
-- putStrLn "mapped cSubTiletoTileFacet"
return Tesselation
return Tessellation
{ _sites = sites
, _tiles = fromAscList tiles'
, _tilefacets = fromAscList subtiles'
Expand Down
24 changes: 12 additions & 12 deletions src/Geometry/Delaunay/Delaunay.hs
Expand Up @@ -28,7 +28,7 @@ delaunay :: [[Double]] -- ^ sites (vertex coordinates)
-> Bool -- ^ whether to add a point at infinity
-> Bool -- ^ whether to include degenerate tiles
-> Maybe Double -- ^ volume threshold
-> IO Tesselation -- ^ Delaunay tessellation
-> IO Tessellation -- ^ Delaunay tessellation
delaunay sites atinfinity degenerate vthreshold = do
let n = length sites
dim = length (head sites)
Expand All @@ -44,7 +44,7 @@ delaunay sites atinfinity degenerate vthreshold = do
sitesPtr <- mallocBytes (n * dim * sizeOf (undefined :: CDouble))
pokeArray sitesPtr (concatMap (map realToFrac) sites)
exitcodePtr <- mallocBytes (sizeOf (undefined :: CUInt))
resultPtr <- c_tesselation sitesPtr
resultPtr <- c_tessellation sitesPtr
(fromIntegral dim) (fromIntegral n)
(fromIntegral $ fromEnum atinfinity)
(fromIntegral $ fromEnum degenerate)
Expand All @@ -58,13 +58,13 @@ delaunay sites atinfinity degenerate vthreshold = do
error $ "qhull returned an error (code " ++ show exitcode ++ ")"
else do
result <- peek resultPtr
out <- cTesselationToTesselation sites result
out <- cTessellationToTessellation sites result
free resultPtr
return 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 :: Tesselation -> Index -> IntMap TileFacet
vertexNeighborFacets :: Tessellation -> Index -> IntMap TileFacet
vertexNeighborFacets tess i = IM.restrictKeys (_tilefacets tess) ids
where
ids = maybe IS.empty _neighfacetsIds (IM.lookup i (_sites tess))
Expand All @@ -74,31 +74,31 @@ sandwichedFacet :: TileFacet -> Bool
sandwichedFacet tilefacet = IS.size (_facetOf tilefacet) == 2

-- | the tiles a facet belongs to
facetOf :: Tesselation -> TileFacet -> IntMap Tile
facetOf :: Tessellation -> TileFacet -> IntMap Tile
facetOf tess tilefacet = IM.restrictKeys (_tiles tess) (_facetOf tilefacet)

-- | the families of the tiles a facet belongs to
facetFamilies :: Tesselation -> TileFacet -> IntMap Family
facetFamilies :: Tessellation -> TileFacet -> IntMap Family
facetFamilies tess tilefacet = IM.map _family (facetOf tess tilefacet)

-- | the circumcenters of the tiles a facet belongs to
facetCenters :: Tesselation -> TileFacet -> IntMap [Double]
facetCenters :: Tessellation -> TileFacet -> IntMap [Double]
facetCenters tess tilefacet =
IM.map _center (facetOf tess tilefacet)

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

-- | the tiles a facet belongs to, facet given by its id
facetOf' :: Tesselation -> Int -> IntMap Tile
facetOf' :: Tessellation -> Int -> IntMap Tile
facetOf' = funofFacetToFunofInt facetOf

-- | the families of the tiles a facet belongs to, facet given by its id
facetFamilies' :: Tesselation -> Int -> IntMap Family
facetFamilies' :: Tessellation -> Int -> IntMap Family
facetFamilies' = funofFacetToFunofInt facetFamilies

-- | the circumcenters of the tiles a facet belongs to, facet given by its id
facetCenters' :: Tesselation -> Int -> IntMap [Double]
facetCenters' :: Tessellation -> Int -> IntMap [Double]
facetCenters' = funofFacetToFunofInt facetCenters
15 changes: 11 additions & 4 deletions src/Geometry/Delaunay/Types.hs
@@ -1,4 +1,11 @@
module Geometry.Delaunay.Types
(
Site (..)
, Simplex (..)
, TileFacet (..)
, Tile (..)
, Tessellation (..)
)
where
import Data.IntMap.Strict (IntMap)
import qualified Data.IntMap.Strict as IM
Expand Down Expand Up @@ -68,18 +75,18 @@ instance HasVolume Tile where
instance HasCenter Tile where
_center = _circumcenter . _simplex

data Tesselation = Tesselation {
data Tessellation = Tessellation {
_sites :: IndexMap Site
, _tiles :: IntMap Tile
, _tilefacets :: IntMap TileFacet
, _edges' :: EdgeMap
} deriving Show

instance HasEdges Tesselation where
instance HasEdges Tessellation where
_edges = _edges'

instance HasVertices Tesselation where
instance HasVertices Tessellation where
_vertices tess = IM.map _point (_sites tess)

instance HasVolume Tesselation where
instance HasVolume Tessellation where
_volume tess = sum (IM.elems $ IM.map (_volume' . _simplex) (_tiles tess))
13 changes: 13 additions & 0 deletions src/Geometry/Qhull/Shared.hs
@@ -1,4 +1,17 @@
module Geometry.Qhull.Shared
(
sameFamily
, verticesIds
, verticesCoordinates
, nVertices
, edgesIds
, edgesIds'
, edgesCoordinates
, nEdges
, isEdge
, toPoints
, toPoints'
)
where
import qualified Data.HashMap.Strict.InsOrd as H
import qualified Data.IntMap.Strict as IM
Expand Down
14 changes: 14 additions & 0 deletions src/Geometry/Qhull/Types.hs
@@ -1,4 +1,18 @@
module Geometry.Qhull.Types
(
Index
, IndexMap
, IndexSet
, IndexPair (..)
, EdgeMap
, Family (..)
, HasCenter (..)
, HasEdges (..)
, HasFamily (..)
, HasNormal (..)
, HasVertices (..)
, HasVolume (..)
)
where
import Data.Hashable
import Data.HashMap.Strict.InsOrd (InsOrdHashMap)
Expand Down

0 comments on commit 92014e2

Please sign in to comment.