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