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]]
-> Bool
-> Bool
-> Maybe FilePath
-> 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
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')
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)
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)
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)
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)
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))
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
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)
where
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)
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)
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
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
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