packages feed

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)