{- HLINT ignore "Use head" -}
module Geometry.ConvexHull.ConvexHull
  where
import           Control.Monad                   ( unless, when )
import           Geometry.ConvexHull.CConvexHull ( c_convexhull, cConvexHullToConvexHull )
import           Geometry.ConvexHull.Types       ( Ridge
                                                 , Facet(
                                                     _fvertices
                                                   , _fridges
                                                   , _fedges
                                                   , _normal'
                                                   , _centroid
                                                   , _area
                                                   )
                                                 , ConvexHull(
                                                     _hfacets
                                                   , _hridges
                                                   , _dimension
                                                   , _simplicial
                                                   ) 
                                                 )
import           Data.Function                   ( on )
import           Data.Graph                      ( flattenSCCs, stronglyConnComp )
import qualified Data.HashMap.Strict.InsOrd      as H
import           Data.IntMap.Strict              ( IntMap )
import qualified Data.IntMap.Strict              as IM
import qualified Data.IntSet                     as IS
import           Data.List                       ( nubBy, findIndices, groupBy, partition, intercalate )
import           Data.List.Index                 ( imap )
import           Data.List.Unique                ( allUnique, count_ )
import           Data.Tuple.Extra                ( both )
import           Foreign.C.String                ( newCString )
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.Shared           ( verticesIds, sameFamily, nVertices, nEdges, edgesIds )
import           Geometry.Qhull.Types            ( IndexPair(Pair)
                                                 , IndexMap
                                                 , Index
                                                 , HasVertices(_vertices)
                                                 , HasNormal(_normal)
                                                 , HasFamily(_family)
                                                 , HasEdges(_edges)
                                                 , Family(None)
                                                 , EdgeMap )
import           Text.Printf                     ( printf )
import           Text.Regex                      ( mkRegex, subRegex )

convexHull :: [[Double]]     -- ^ vertices

           -> Bool           -- ^ whether to triangulate

           -> Bool           -- ^ whether to print output to stdout

           -> Maybe FilePath -- ^ write summary to a file

           -> IO ConvexHull
convexHull :: [[Double]] -> Bool -> Bool -> Maybe String -> IO ConvexHull
convexHull [[Double]]
points Bool
triangulate Bool
stdout Maybe String
file = do
  let n :: Int
n   = forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
points
      dim :: Int
dim = forall (t :: * -> *) a. Foldable t => t a -> Int
length (forall a. [a] -> a
head [[Double]]
points)
  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 => String -> a
error String
"dimension must be at least 2"
  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 b c a. (b -> c) -> (a -> b) -> a -> c
. forall (t :: * -> *) a. Foldable t => t a -> Int
length) (forall a. [a] -> [a]
tail [[Double]]
points)) forall a b. (a -> b) -> a -> b
$
    forall a. HasCallStack => String -> a
error String
"the points must have the same dimension"
  forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
n forall a. Ord a => a -> a -> Bool
<= Int
dim) forall a b. (a -> b) -> a -> b
$
    forall a. HasCallStack => String -> a
error String
"insufficient number of points"
  forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
unless (forall a. Ord a => [a] -> Bool
allUnique [[Double]]
points) forall a b. (a -> b) -> a -> b
$
    forall a. HasCallStack => String -> a
error String
"some points are duplicated"
  Ptr CDouble
pointsPtr <- 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
pointsPtr (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]]
points)
  Ptr CUInt
exitcodePtr <- forall a. Int -> IO (Ptr a)
mallocBytes (forall a. Storable a => a -> Int
sizeOf (forall a. HasCallStack => a
undefined :: CUInt))
  CString
summaryFile <- forall b a. b -> (a -> b) -> Maybe a -> b
maybe (String -> IO CString
newCString []) String -> IO CString
newCString Maybe String
file
  Ptr CConvexHull
resultPtr <- Ptr CDouble
-> CUInt
-> CUInt
-> CUInt
-> CUInt
-> CString
-> Ptr CUInt
-> IO (Ptr CConvexHull)
c_convexhull Ptr CDouble
pointsPtr (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
triangulate)
               (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
stdout) CString
summaryFile 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
pointsPtr
  if CUInt
exitcode forall a. Eq a => a -> a -> Bool
/= CUInt
0
    then do
      forall a. Ptr a -> IO ()
free Ptr CConvexHull
resultPtr
      forall a. HasCallStack => String -> a
error forall a b. (a -> b) -> a -> b
$ String
"qhull returned an error (code " forall a. [a] -> [a] -> [a]
++ forall a. Show a => a -> String
show CUInt
exitcode forall a. [a] -> [a] -> [a]
++ String
")"
    else do
      ConvexHull
result <- forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
(>>=) (forall a. Storable a => Ptr a -> IO a
peek Ptr CConvexHull
resultPtr) CConvexHull -> IO ConvexHull
cConvexHullToConvexHull
      forall a. Ptr a -> IO ()
free Ptr CConvexHull
resultPtr
      forall (m :: * -> *) a. Monad m => a -> m a
return ConvexHull
result


-- | convex hull summary

hullSummary :: ConvexHull -> String
hullSummary :: ConvexHull -> String
hullSummary ConvexHull
hull =
  String
"Convex hull:\n" forall a. [a] -> [a] -> [a]
++
  forall r. PrintfType r => String -> r
printf String
"%d vertices\n" (forall a. IntMap a -> Int
IM.size IndexMap [Double]
vertices) forall a. [a] -> [a] -> [a]
++
  forall r. PrintfType r => String -> r
printf String
"%d facets (%s)\n" Int
nfacets String
families forall a. [a] -> [a] -> [a]
++
  forall r. PrintfType r => String -> r
printf String
"%d ridges\n" Int
nridges forall a. [a] -> [a] -> [a]
++
  forall r. PrintfType r => String -> r
printf String
"%d edges\n" Int
nedges forall a. [a] -> [a] -> [a]
++
  (if Int
dim forall a. Ord a => a -> a -> Bool
> Int
2
    then forall r. PrintfType r => String -> r
printf String
"number of vertices per facet: %s\n" (forall a. Show a => a -> String
show [(Int, Int)]
counts_vertices) forall a. [a] -> [a] -> [a]
++
         (if Int
dim forall a. Ord a => a -> a -> Bool
> Int
3
           then forall r. PrintfType r => String -> r
printf String
"number of edges per facet: %s\n" (forall a. Show a => a -> String
show [(Int, Int)]
counts_edges) forall a. [a] -> [a] -> [a]
++
                forall r. PrintfType r => String -> r
printf String
"number of ridges per facet: %s\n" (forall a. Show a => a -> String
show [(Int, Int)]
counts_ridges)
           else String
"")
    else String
"")
  where
    vertices :: IndexMap [Double]
vertices = forall m. HasVertices m => m -> IndexMap [Double]
_vertices ConvexHull
hull
    nedges :: Int
nedges = forall a. HasEdges a => a -> Int
nEdges ConvexHull
hull
    nridges :: Int
nridges = forall a. IntMap a -> Int
IM.size (ConvexHull -> IntMap Ridge
_hridges ConvexHull
hull)
    facets :: IntMap Facet
facets = ConvexHull -> IntMap Facet
_hfacets ConvexHull
hull
    facets' :: [Facet]
facets' = forall a. IntMap a -> [a]
IM.elems IntMap Facet
facets
    nfacets :: Int
nfacets = forall a. IntMap a -> Int
IM.size IntMap Facet
facets
    (Int
nf1,Int
nf2) = forall a b. (a -> b) -> (a, a) -> (b, b)
both forall (t :: * -> *) a. Foldable t => t a -> Int
length forall a b. (a -> b) -> a -> b
$
                forall a. (a -> Bool) -> [a] -> ([a], [a])
partition (Family
None forall a. Eq a => a -> a -> Bool
==) (forall a. (a -> a -> Bool) -> [a] -> [a]
nubBy Family -> Family -> Bool
sameFamily (forall a b. (a -> b) -> [a] -> [b]
map forall m. HasFamily m => m -> Family
_family [Facet]
facets'))
    families :: String
families = forall a. Show a => a -> String
show Int
nf1 forall a. [a] -> [a] -> [a]
++ String
" single, " forall a. [a] -> [a] -> [a]
++
               forall a. Show a => a -> String
show Int
nf2 forall a. [a] -> [a] -> [a]
++ if Int
nf2 forall a. Ord a => a -> a -> Bool
> Int
1 then String
" families" else String
" family"
    dim :: Int
dim = forall (t :: * -> *) a. Foldable t => t a -> Int
length forall a b. (a -> b) -> a -> b
$ forall a. [a] -> a
head (forall a. IntMap a -> [a]
IM.elems IndexMap [Double]
vertices)
    counts_vertices :: [(Int, Int)]
counts_vertices = forall a. Ord a => [a] -> [(a, Int)]
count_ (forall a b. (a -> b) -> [a] -> [b]
map forall a. HasVertices a => a -> Int
nVertices [Facet]
facets')
    counts_edges :: [(Int, Int)]
counts_edges = forall a. Ord a => [a] -> [(a, Int)]
count_ (forall a b. (a -> b) -> [a] -> [b]
map forall a. HasEdges a => a -> Int
nEdges [Facet]
facets')
    counts_ridges :: [(Int, Int)]
counts_ridges = forall a. Ord a => [a] -> [(a, Int)]
count_ (forall a b. (a -> b) -> [a] -> [b]
map (IntSet -> Int
IS.size forall b c a. (b -> c) -> (a -> b) -> a -> c
. Facet -> IntSet
_fridges) [Facet]
facets')

-- | volume of the convex hull (area in dimension 2, volume in dimension 3, 

-- hypervolume in higher dimension)

hullVolume :: ConvexHull -> Double
hullVolume :: ConvexHull -> Double
hullVolume ConvexHull
hull = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (forall a b. (a -> b) -> [a] -> [b]
map Facet -> Double
fvolume [Facet]
facets) forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral (ConvexHull -> Int
_dimension ConvexHull
hull)
  where
    fvolume :: Facet -> Double
    fvolume :: Facet -> Double
fvolume Facet
fct = 
      forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
(*) (Facet -> [Double]
_centroid Facet
fct) (Facet -> [Double]
_normal' Facet
fct)) forall a. Num a => a -> a -> a
* Facet -> Double
_area Facet
fct
    facets :: [Facet]
facets = forall a. IntMap a -> [a]
IM.elems (ConvexHull -> IntMap Facet
_hfacets ConvexHull
hull)

-- | facets ids an edge belongs to

edgeOf :: ConvexHull -> (Index, Index) -> [Int]
edgeOf :: ConvexHull -> (Int, Int) -> [Int]
edgeOf ConvexHull
hull (Int
v1,Int
v2) = forall a. IntMap a -> [Int]
IM.keys forall a b. (a -> b) -> a -> b
$ forall a. (a -> Bool) -> IntMap a -> IntMap a
IM.filter (forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
elem (Int -> Int -> IndexPair
Pair Int
v1 Int
v2)) IntMap [IndexPair]
facetsEdges
  where
    facetsEdges :: IntMap [IndexPair]
facetsEdges = forall a b. (a -> b) -> IntMap a -> IntMap b
IM.map forall a. HasEdges a => a -> [IndexPair]
edgesIds (ConvexHull -> IntMap Facet
_hfacets ConvexHull
hull)

-- | ridges of a facet

facetRidges :: ConvexHull -> Facet -> IntMap Ridge
facetRidges :: ConvexHull -> Facet -> IntMap Ridge
facetRidges ConvexHull
hull Facet
facet = forall a. IntMap a -> IntSet -> IntMap a
IM.restrictKeys (ConvexHull -> IntMap Ridge
_hridges ConvexHull
hull) (Facet -> IntSet
_fridges Facet
facet)

-- | vertices ids of all facets

facetsVerticesIds :: ConvexHull -> [[Index]]
facetsVerticesIds :: ConvexHull -> [[Int]]
facetsVerticesIds ConvexHull
hull = forall a b. (a -> b) -> [a] -> [b]
map forall a. HasVertices a => a -> [Int]
verticesIds (forall a. IntMap a -> [a]
IM.elems forall a b. (a -> b) -> a -> b
$ ConvexHull -> IntMap Facet
_hfacets ConvexHull
hull)

-- | vertices ids of all ridges

ridgesVerticesIds :: ConvexHull -> [[Index]]
ridgesVerticesIds :: ConvexHull -> [[Int]]
ridgesVerticesIds ConvexHull
hull = forall a b. (a -> b) -> [a] -> [b]
map forall a. HasVertices a => a -> [Int]
verticesIds (forall a. IntMap a -> [a]
IM.elems (ConvexHull -> IntMap Ridge
_hridges ConvexHull
hull))

-- | group facets of the same family

groupedFacets :: ConvexHull -> [(Family, [IndexMap [Double]], [EdgeMap])]
groupedFacets :: ConvexHull -> [(Family, [IndexMap [Double]], [EdgeMap])]
groupedFacets ConvexHull
hull =
  forall a b c. [a] -> [b] -> [c] -> [(a, b, c)]
zip3 (forall a b. (a -> b) -> [a] -> [b]
map forall a. [a] -> a
head [[Family]]
families) [[IndexMap [Double]]]
verticesGroups [[EdgeMap]]
edgesGroups
  where
    facets :: [Facet]
facets         = forall a. IntMap a -> [a]
IM.elems (ConvexHull -> IntMap Facet
_hfacets ConvexHull
hull)
    facetsGroups :: [[Facet]]
facetsGroups   = forall a. (a -> a -> Bool) -> [a] -> [[a]]
groupBy (Family -> Family -> Bool
sameFamily forall b c a. (b -> b -> c) -> (a -> b) -> a -> a -> c
`on` forall m. HasFamily m => m -> Family
_family) [Facet]
facets
    edgesGroups :: [[EdgeMap]]
edgesGroups    = forall a b. (a -> b) -> [a] -> [b]
map (forall a b. (a -> b) -> [a] -> [b]
map Facet -> EdgeMap
_fedges) [[Facet]]
facetsGroups
    verticesGroups :: [[IndexMap [Double]]]
verticesGroups = forall a b. (a -> b) -> [a] -> [b]
map (forall a b. (a -> b) -> [a] -> [b]
map Facet -> IndexMap [Double]
_fvertices) [[Facet]]
facetsGroups
    families :: [[Family]]
families       = forall a b. (a -> b) -> [a] -> [b]
map (forall a b. (a -> b) -> [a] -> [b]
map forall m. HasFamily m => m -> Family
_family) [[Facet]]
facetsGroups

-- | group facets of the same family and merge vertices and edges

groupedFacets' :: ConvexHull -> [(Family, IndexMap [Double], EdgeMap)]
groupedFacets' :: ConvexHull -> [(Family, IndexMap [Double], EdgeMap)]
groupedFacets' ConvexHull
hull =
  forall a b. (a -> b) -> [a] -> [b]
map (\(Family
f,[IndexMap [Double]]
v,[EdgeMap]
e) -> (Family
f, forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr forall a. IntMap a -> IntMap a -> IntMap a
IM.union forall a. IntMap a
IM.empty [IndexMap [Double]]
v, forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr EdgeMap -> EdgeMap -> EdgeMap
delta forall k v. InsOrdHashMap k v
H.empty [EdgeMap]
e))
      (ConvexHull -> [(Family, [IndexMap [Double]], [EdgeMap])]
groupedFacets ConvexHull
hull)
  -- zip3 (map head families) (map (foldr IM.union IM.empty) verticesGroups)

  --      (map (foldr delta H.empty) edgesGroups)

  where
    -- facets         = IM.elems (_hfacets hull)

    -- facetsGroups   = groupBy (sameFamily `on` _family) facets

    -- edgesGroups    = map (map _fedges) facetsGroups

    -- verticesGroups = map (map _fvertices) facetsGroups

    -- families       = map (map _family) facetsGroups

    delta :: EdgeMap -> EdgeMap -> EdgeMap
    delta :: EdgeMap -> EdgeMap -> EdgeMap
delta EdgeMap
e1 EdgeMap
e2 = forall k v w.
(Eq k, Hashable k) =>
InsOrdHashMap k v -> InsOrdHashMap k w -> InsOrdHashMap k v
H.difference (forall k v.
(Eq k, Hashable k) =>
InsOrdHashMap k v -> InsOrdHashMap k v -> InsOrdHashMap k v
H.union EdgeMap
e1 EdgeMap
e2) (forall k v w.
(Eq k, Hashable k) =>
InsOrdHashMap k v -> InsOrdHashMap k w -> InsOrdHashMap k v
H.intersection EdgeMap
e1 EdgeMap
e2)

-- data Vertex3 = Vertex3 Double Double Double

--   deriving Show

--

-- toVertex3 :: [Double] -> Vertex3

-- toVertex3 xs = Vertex3 (xs!!0) (xs!!1) (xs!!2)


-- | for 3D only, orders the vertices of the facet (i.e. provides a polygon);

-- also returns a Boolean indicating the orientation of the vertices

facetToPolygon :: Facet -> ([(Index, [Double])], Bool)
facetToPolygon :: Facet -> ([(Int, [Double])], Bool)
facetToPolygon Facet
facet = ([(Int, [Double])]
polygon, Double
dotProduct forall a. Ord a => a -> a -> Bool
> Double
0)
  where
  vs :: [(Int, [Double])]
vs = forall a. IntMap a -> [(Int, a)]
IM.toList forall a b. (a -> b) -> a -> b
$ forall m. HasVertices m => m -> IndexMap [Double]
_vertices Facet
facet
  x :: [((Int, [Double]), Int, [Int])]
x = forall a b. (Int -> a -> b) -> [a] -> [b]
imap (\Int
i (Int, [Double])
v -> ((Int, [Double])
v, Int
i, forall a. (a -> Bool) -> [a] -> [Int]
findIndices ((Int, [Double]) -> (Int, [Double]) -> Bool
connectedVertices (Int, [Double])
v) [(Int, [Double])]
vs)) [(Int, [Double])]
vs
    where
    connectedVertices :: (Index, [Double]) -> (Index, [Double]) -> Bool
    connectedVertices :: (Int, [Double]) -> (Int, [Double]) -> Bool
connectedVertices (Int
i,[Double]
_) (Int
j,[Double]
_) = Int -> Int -> IndexPair
Pair Int
i Int
j forall k a. (Eq k, Hashable k) => k -> InsOrdHashMap k a -> Bool
`H.member` forall m. HasEdges m => m -> EdgeMap
_edges Facet
facet
  polygon :: [(Int, [Double])]
polygon = forall a. [SCC a] -> [a]
flattenSCCs (forall key node. Ord key => [(node, key, [key])] -> [SCC node]
stronglyConnComp [((Int, [Double]), Int, [Int])]
x)
  vertices :: [[Double]]
vertices = forall a b. (a -> b) -> [a] -> [b]
map forall a b. (a, b) -> b
snd [(Int, [Double])]
polygon
  v1 :: [Double]
v1 = [[Double]]
vertices forall a. [a] -> Int -> a
!! Int
0
  v2 :: [Double]
v2 = [[Double]]
vertices forall a. [a] -> Int -> a
!! Int
1
  v3 :: [Double]
v3 = [[Double]]
vertices forall a. [a] -> Int -> a
!! Int
2
  normal :: [Double]
normal = forall {a}. Num a => [a] -> [a] -> [a]
crossProd (forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
subtract [Double]
v1 [Double]
v2) (forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
subtract [Double]
v1 [Double]
v3)
    where
    crossProd :: [a] -> [a] -> [a]
crossProd [a]
u [a]
v = [ [a]
uforall a. [a] -> Int -> a
!!Int
1 forall a. Num a => a -> a -> a
* [a]
vforall a. [a] -> Int -> a
!!Int
2 forall a. Num a => a -> a -> a
- [a]
uforall a. [a] -> Int -> a
!!Int
2 forall a. Num a => a -> a -> a
* [a]
vforall a. [a] -> Int -> a
!!Int
1
                    , [a]
uforall a. [a] -> Int -> a
!!Int
2 forall a. Num a => a -> a -> a
* [a]
vforall a. [a] -> Int -> a
!!Int
0 forall a. Num a => a -> a -> a
- [a]
uforall a. [a] -> Int -> a
!!Int
0 forall a. Num a => a -> a -> a
* [a]
vforall a. [a] -> Int -> a
!!Int
2
                    , [a]
uforall a. [a] -> Int -> a
!!Int
0 forall a. Num a => a -> a -> a
* [a]
vforall a. [a] -> Int -> a
!!Int
1 forall a. Num a => a -> a -> a
- [a]
uforall a. [a] -> Int -> a
!!Int
1 forall a. Num a => a -> a -> a
* [a]
vforall a. [a] -> Int -> a
!!Int
0 ]
  dotProduct :: Double
dotProduct = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum forall a b. (a -> b) -> a -> b
$ forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
(*) [Double]
normal (forall m. HasNormal m => m -> [Double]
_normal Facet
facet)

-- | for 3D only, orders the vertices of the facet (i.e. provides a polygon)

-- in anticlockwise orientation

facetToPolygon' :: Facet -> [(Index, [Double])]
facetToPolygon' :: Facet -> [(Int, [Double])]
facetToPolygon' Facet
facet = if Bool
test then [(Int, [Double])]
polygon else forall a. [a] -> [a]
reverse [(Int, [Double])]
polygon
  where
  ([(Int, [Double])]
polygon, Bool
test) = Facet -> ([(Int, [Double])], Bool)
facetToPolygon Facet
facet

-- -- | like `facetToPolygon`, but returns the vertices indices

-- facetToPolygon' :: Facet -> [Index]

-- facetToPolygon' facet = map fst $ flattenSCCs (stronglyConnComp x)

--   where

--     vs = IM.toList $ _vertices facet

--     x = imap (\i v -> (v, i, findIndices (connectedVertices v) vs)) vs

--     connectedVertices :: (Index, [Double]) -> (Index, [Double]) -> Bool

--     connectedVertices (i,_) (j,_) = Pair i j `H.member` _edges facet


-- | for 4D only, orders the vertices of a ridge (i.e. provides a polygon)

ridgeToPolygon :: Ridge -> [(Index, [Double])]
ridgeToPolygon :: Ridge -> [(Int, [Double])]
ridgeToPolygon Ridge
ridge = forall a. [SCC a] -> [a]
flattenSCCs (forall key node. Ord key => [(node, key, [key])] -> [SCC node]
stronglyConnComp [((Int, [Double]), Int, [Int])]
x)
  where
  vs :: [(Int, [Double])]
vs = forall a. IntMap a -> [(Int, a)]
IM.toList forall a b. (a -> b) -> a -> b
$ forall m. HasVertices m => m -> IndexMap [Double]
_vertices Ridge
ridge
  x :: [((Int, [Double]), Int, [Int])]
x = forall a b. (Int -> a -> b) -> [a] -> [b]
imap (\Int
i (Int, [Double])
v -> ((Int, [Double])
v, Int
i, forall a. (a -> Bool) -> [a] -> [Int]
findIndices ((Int, [Double]) -> (Int, [Double]) -> Bool
connectedVertices (Int, [Double])
v) [(Int, [Double])]
vs)) [(Int, [Double])]
vs
    where
    connectedVertices :: (Index, [Double]) -> (Index, [Double]) -> Bool
    connectedVertices :: (Int, [Double]) -> (Int, [Double]) -> Bool
connectedVertices (Int
i,[Double]
_) (Int
j,[Double]
_) = Int -> Int -> IndexPair
Pair Int
i Int
j forall k a. (Eq k, Hashable k) => k -> InsOrdHashMap k a -> Bool
`H.member` forall m. HasEdges m => m -> EdgeMap
_edges Ridge
ridge

-- | for 3D only, convert the convex hull to a STL file;

-- the STL format is valid for triangle meshes only

hullToSTL :: ConvexHull -> FilePath -> IO ()
hullToSTL :: ConvexHull -> String -> IO ()
hullToSTL ConvexHull
chull String
filename = do
  forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (ConvexHull -> Int
_dimension ConvexHull
chull forall a. Eq a => a -> a -> Bool
/= Int
3) forall a b. (a -> b) -> a -> b
$
    forall a. HasCallStack => String -> a
error String
"dimension is not 3"
  forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
unless (ConvexHull -> Bool
_simplicial ConvexHull
chull) forall a b. (a -> b) -> a -> b
$
    forall a. Show a => a -> IO ()
print String
"* are you sure this convex hull is triangulated?\n"
  let facets :: [Facet]
facets = forall a. IntMap a -> [a]
IM.elems (ConvexHull -> IntMap Facet
_hfacets ConvexHull
chull)
      normals :: [[Double]]
normals = forall a b. (a -> b) -> [a] -> [b]
map forall m. HasNormal m => m -> [Double]
_normal [Facet]
facets
      facetNormals :: [String]
facetNormals = 
        forall a b. (a -> b) -> [a] -> [b]
map (\[Double]
n -> String
"facet normal  " forall a. [a] -> [a] -> [a]
++ forall a. Show a => String -> [a] -> String
stringify String
" " [Double]
n forall a. [a] -> [a] -> [a]
++ String
"\nouter loop\n")
            [[Double]]
normals
      polygons :: [String]
polygons = 
        forall a b. (a -> b) -> [a] -> [b]
map (\Facet
f -> String
"vertex  " forall a. [a] -> [a] -> [a]
++
              (forall a. Show a => String -> [a] -> String
stringify String
"\nvertex  " forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (a -> b) -> [a] -> [b]
map forall a b. (a, b) -> b
snd forall b c a. (b -> c) -> (a -> b) -> a -> c
. Facet -> [(Int, [Double])]
facetToPolygon') Facet
f) [Facet]
facets
      vertices :: [String]
vertices = forall a b. (a -> b) -> [a] -> [b]
map (forall a. [a] -> [a] -> [a]
++ String
"\nendloop\nendfacet\n") [String]
polygons
      vertices' :: [String]
vertices' = 
        forall a b. (a -> b) -> [a] -> [b]
map (\String
v -> Regex -> String -> String -> String
subRegex (String -> Regex
mkRegex String
"\\]") (Regex -> String -> String -> String
subRegex (String -> Regex
mkRegex String
"\\[") String
v String
"") String
"") 
            [String]
vertices
      out :: String
out = Regex -> String -> String -> String
subRegex 
        (String -> Regex
mkRegex String
",") (forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [String
x forall a. [a] -> [a] -> [a]
++ String
y | String
x <- [String]
facetNormals, String
y <- [String]
vertices']) String
" "
  String -> String -> IO ()
writeFile String
filename (String
"solid " forall a. [a] -> [a] -> [a]
++ String
filename forall a. [a] -> [a] -> [a]
++ String
" produced by Haskell\n" forall a. [a] -> [a] -> [a]
++ String
out)
  where
    stringify :: Show a => String -> [a] -> String
    stringify :: forall a. Show a => String -> [a] -> String
stringify String
sep = forall a. [a] -> [[a]] -> [a]
intercalate String
sep forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (a -> b) -> [a] -> [b]
map forall a. Show a => a -> String
show