packages feed

qhull-0.1.0.1: src/Voronoi/R.hs

module Voronoi.R
  where
import           ConvexHull         (convexHull, ConvexHull (..))
import qualified Data.IntMap.Strict as IM
import           Data.List          (intercalate, transpose)
import           Data.List.Index    (imap, iconcatMap)
import           Data.Maybe
import           Delaunay.R
import           Delaunay.Types     (Tesselation)
import           Qhull.Shared
import           Voronoi2D
import           Voronoi3D

voronoi2ForR :: Voronoi2 -> Maybe Tesselation -> String
voronoi2ForR v d =
  (if isJust d then dcode else "")
  ++ "colors <- rainbow(" ++ show (length boundedCells +1) ++ ")\n"
  ++ unlines (imap boundedCellForR boundedCells)
  ++ unlines (map cellForR v)
  where
    dcode = delaunay2ForR (fromJust d) True
    cellForR :: ([Double], Cell2) -> String
    cellForR (site, edges) =
      point ++ "\n" ++ unlines (map f edges)
      where
        point =
          "points(" ++ show (site!!0) ++ ", " ++ show (site!!1) ++
                    ", pch=19, col=\"blue\")"
        f :: Edge2 -> String
        f edge = case edge of
          Edge2 ((x0,y0),(x1,y1)) ->
            "segments(" ++ intercalate "," (map show [x0,y0,x1,y1]) ++
                        ", col=\"black\", lty=1, lwd=1)"
          IEdge2 ((x0,y0),(x1,y1)) ->
            "segments(" ++ intercalate "," (map show [x0,y0,x0+x1,y0+y1]) ++
                        ", col=\"red\", lty=1, lwd=1)"
          TIEdge2 ((x0,y0),(x1,y1)) ->
            "segments(" ++ intercalate "," (map show [x0,y0,x1,y1]) ++
                        ", col=\"red\", lty=1, lwd=1)"
    boundedCells = map snd (restrictVoronoi2 v)
    boundedCellForR :: Int -> Cell2 -> String
    boundedCellForR i cell =
      "polygon(c(" ++ intercalate "," (map show x) ++ "), c(" ++
      intercalate "," (map show y) ++ "), col=colors[" ++ show (i+1) ++ "])\n"
      where
        [x, y] = transpose (cell2Vertices' cell)


voronoi3ForRgl :: Voronoi3 -> Maybe Int -> Maybe Tesselation -> String
voronoi3ForRgl v n d =
  let v' = if isJust n then roundVoronoi3 (fromJust n) v else v in
  let code = unlines $ map cellForRgl v' in
  "library(rgl)\n" ++
  if isJust d
    then code ++ "\n" ++ "# Delaunay:\n" ++
         delaunay3rgl (fromJust d) True True True True (Just 0.9)
    else code
  where
    cellForRgl :: ([Double], Cell3) -> String
    cellForRgl (site, cell) = plotpoint ++ unlines (map f cell)
      where
        plotpoint = "spheres3d(" ++ intercalate "," (map show site) ++ ", radius=0.1, color=\"red\")\n"
        f :: Edge3 -> String
        f edge = case edge of
          Edge3 (x,y) ->
            "segments3d(rbind(c" ++ show x ++ ", \n\tc" ++ show y ++ "))"
          TIEdge3 (x,y) ->
            "segments3d(rbind(c" ++ show x ++ ", \n\tc" ++ show y ++ "), col=c(\"red\",\"red\"))"
          IEdge3 (x,y) ->
            "segments3d(rbind(c" ++ show x ++ ", \n\tc" ++ show (sumTriplet x y) ++ "), col=c(\"red\",\"red\"))"
        sumTriplet (a,b,c) (a',b',c') = (a+a',b+b',c+c')

-- | plot with facets
voronoi3ForRgl' :: Voronoi3 -> Maybe Int -> Maybe Tesselation -> IO String
voronoi3ForRgl' v n d = do -- faudrait un argument approx
  let code1 = voronoi3ForRgl v n d
      v' = restrictVoronoi3' v
      v'' = if isJust n then roundVoronoi3 (fromJust n) v' else v'
      boundedCells = map (cell3Vertices . snd) (restrictVoronoi3' v'')
--      boundedCells' = map (nub . map (map (approx 13))) boundedCells
  hulls <- mapM (\cell -> convexHull cell True False Nothing) boundedCells
  let triangles = map (map verticesCoordinates . IM.elems . _hfacets) hulls
      code_colors = "colors <- rainbow(" ++ show (length triangles +1) ++ ")\n"
      code2 = iconcatMap (\i x -> concatMap (rglTriangle i) x ++ "\n") triangles
  return $ code_colors ++ code1 ++ code2
  where
    -- approx :: RealFrac a => Int -> a -> a
    -- approx n x = fromInteger (round $ x * (10^n)) / (10.0^^n)
    asTriplet p = (p!!0, p!!1, p!!2)
    rglTriangle :: Int -> [[Double]] -> String
    rglTriangle i threepoints =
      "triangles3d(rbind(c" ++ show p1 ++ ", \n\tc" ++ show p2 ++
      ", \n\tc" ++ show p3 ++ "), color=colors[" ++ show (i+1) ++ "]" ++
      ", alpha=0.75)\n"
      where
        p1 = asTriplet $ threepoints!!0
        p2 = asTriplet $ threepoints!!1
        p3 = asTriplet $ threepoints!!2