qhull-0.1.0.1: src/Delaunay/R.hs
module Delaunay.R
where
-- import qualified Data.HashMap.Strict.InsOrd as H
import qualified Data.IntMap.Strict as IM
import qualified Data.IntSet as IS
import Data.List
import Data.List.Index (iconcatMap)
-- import Data.Maybe
import Delaunay
import Text.Printf
-- | R code to plot a 2D Delaunay tesselation
delaunay2ForR :: Tesselation -> Bool -> String
delaunay2ForR tess colors =
let tiles = IM.elems (_tiles tess) in
"plot(0, 0, type=\"n\", xlim=c(0,5), ylim=c(0,5)) # please set the limits\n" ++
(if colors
then printf "colors <- heat.colors(%d, alpha=0.5)\n" (length tiles + 1)
else "\n") ++
concatMap triangle (zip [1 .. length tiles] tiles)
where
triangle :: (Int, Tile) -> String
triangle (i, tile) =
let pts = map (\p -> [p!!0,p!!1,p!!2]) (verticesCoordinates tile)
in
printf "polygon(c(%f,%f,%f), c(%f,%f,%f), border=\"black\", "
(pts!!0!!0) (pts!!1!!0) (pts!!2!!0)
(pts!!0!!1) (pts!!1!!1) (pts!!2!!1) ++
-- "polygon(c(" ++ show (pts!!0!!0) ++ ", " ++ show (pts!!1!!0) ++
-- ", " ++ show (pts!!2!!0) ++ "), "
-- ++ "c(" ++ show (pts!!0!!1) ++ ", " ++ show (pts!!1!!1) ++
-- ", " ++ show (pts!!2!!1) ++
-- "), border=\"black\", " ++
(if colors
then printf "col=colors[%d])\n" i
else "col=\"lightblue\")\n")
-- | R code to plot a 3D Delaunay tesselation
delaunay3rgl :: Tesselation -> Bool -> Bool -> Bool -> Bool -> Maybe Double -> String
delaunay3rgl tess interior exterior segments colors alpha =
let allridges = IM.elems (_tilefacets tess) in
let ridges | exterior && interior = allridges
| interior = filter sandwichedFacet allridges
| exterior = filter (not . sandwichedFacet) allridges
in
"library(rgl)\n" ++
(if colors
then printf "colors <- topo.colors(%d, alpha=0.5)\n" (length ridges + 1)
else "\n") ++
concatMap rglRidge ridges ++
(if segments
then if exterior && interior
then concatMap rglSegment (edgesCoordinates tess)
else concatMap rglSegment (tilefacetsEdges ridges)
else "")
where
rglRidge :: TileFacet -> String
rglRidge ridge =
let i = 1 + head (IS.elems $ _facetOf ridge) in
printf "\ntriangles3d(rbind(c%s,c%s,c%s"
(show (pts!!0)) (show (pts!!1)) (show (pts!!2)) ++
(if colors
then
printf "), color=colors[%d]" i
else
"), color=\"blue\"") ++
maybe ")\n" (printf "alpha=%f)\n") alpha
where
pts = map (\p -> (p!!0,p!!1,p!!2)) (verticesCoordinates ridge)
rglSegment :: ([Double], [Double]) -> String
rglSegment (p1, p2) =
printf "segments3d(rbind(c%s,c%s), color=\"black\")\n"
(show p1') (show p2')
where
p1' = (p1!!0, p1!!1, p1!!2)
p2' = (p2!!0, p2!!1, p2!!2)
tilefacetsEdges :: [TileFacet] -> [([Double], [Double])]
tilefacetsEdges tilefacets = foldl' union [] (map tilefacetEdges tilefacets)
where
tilefacetEdges :: TileFacet -> [([Double], [Double])]
tilefacetEdges tilefacet = [(v!!0,v!!1),(v!!1,v!!2),(v!!2,v!!0)]
where
v = verticesCoordinates tilefacet
delaunaySpheres :: Tesselation -> String
delaunaySpheres tess =
let tiles = IM.elems (_tiles tess) in
"library(rgl)\n" ++
printf "colors <- rainbow(%d)\n" (length tiles) ++
iconcatMap rglSphere tiles
where
rglSphere :: Int -> Tile -> String
rglSphere i tile =
printf "spheres3d(%f, %f, %f, radius=%f, color=colors[%d], alpha=0.75)\n"
(c!!0) (c!!1) (c!!2) r (i+1)
where
c = _center tile
r = _circumradius (_simplex tile)