######### ZOBRAZENI 3D DAT ######################################################################################################
# nahodne vytvorime 3 vektory souradnic x, y, z
x <- c(3,5,7,9,8,2,1,1,1,4)
y <- c(5,4,1,1,0,6,9,5,7,2)
z <- c(2,1,2,0,2,1,0,4,2,4)
# install.packages("rgl")
library(rgl)
rgl.open()                                        # otevreni prazdneho RGL pole (sam. okna); note, nejsou videt osy
rgl.spheres(1,1,1, r = 0.2, color="violetred4")   # vykresleni samostatneho bodu na souradnicich x = 1, y = 1, z = 1
rgl.bbox(color="darkblue")                        # dokreslime box
# note: je vykreslen bod - sphere ("kulicka") - se stredem v bode S = 1, 1, 1 a polomeru r = 0.2
rgl.spheres(x,y,z, r = 0.2, color="lightgray")    # prikreslime body o souradnicich x, y, z (10 novych bodu)
# grafy lze natacet a koleckem misi priblizovat event. oddalovat

# moznost zobrazeni bez vyvolavani okna
plot3d(x,y,z, type = "s", r = 0.2, col = "lightgray")
### parametr type:  'p' for points, 's' for spheres, 'l' for lines, 'h' for line segments from z = 0, and 'n' for nothing

# dva soubory dat:
x1 <- c(3,5,7,9,8,2,1,1,1,4)
y1 <- c(5,4,1,1,0,6,9,5,7,2)
z1 <- c(2,1,2,0,2,1,0,4,2,4)
z2 <- c(3,5,7,9,8,2,1,1,1,4)
x2 <- c(5,4,1,1,0,6,9,5,7,2)
y2 <- c(2,1,2,0,2,1,0,4,2,4)
library(rgl)
rgl.open()
rgl.spheres(x1,y1,z1, r = 0.2, color="lightgray")
rgl.spheres(x2,y2,z2, r = 0.2, color="tomato")
rgl.bbox(color="darkblue")

# ZAKLADNI PARAMETRY
rgl.open()                                                # new RGL device
rgl.bg(color = "aliceblue")                               # definuje pozadi otevreneho okna
rgl.viewpoint(theta = 20, phi = 45, zoom = 3)             # zorny uhel (argumenty theta a phi jsou polarni souradnice); zoom - priblizeni (preddefinovano 1)
rgl.spheres(x,y,z, r = 0.2, color="blueviolet")           # vykresleni bodu o souradnicich x, y, z; polomer bodu 0.2 + barva bodu 
# vykreslime box
rgl.bbox(color = c("darkolivegreen3", "black"), xlen=10, ylen=10, zlen=10, marklen=15.9, shininess=8, alpha=0.8)
# vektor barev - prvni jake budou steny, druha popisky
# xlen, ylen, zlen - pocet znacek (tickmarks) na prislusnych osach
# marklen - udava jak dlouhe znacky (tickmarks) budou
# shininess - lesk boxu (0 = zcela bily)
# alpha - udava prusvitnost boxu, od 0.0 (komplene prusvitne) do 1.0 (nepruhledne)

# pokud nechceme zobrazovat box, muzeme zobrazit osy - pridanim rgl.lines o patricnych souradnicich 
rgl.open()                                                
rgl.bg(color = "aliceblue")                               
rgl.viewpoint(theta = 20, phi = 45, zoom = 3)             
rgl.spheres(x,y,z, r = 0.2, color="blueviolet")  
rgl.lines(c(0, 10), c(0, 0), c(0, 0), color = "midnightblue")           # osa x
rgl.lines(c(0, 0), c(0,10), c(0, 0), color = "coral4")                  # osa y
rgl.lines(c(0, 0), c(0, 0), c(0,10), color = "seagreen")                # osa z
# NOTE: obecne se zobrazuje osa x doprava, osa y smerujici dopredu a osa z smerujici nahoru
# zde budeme spise pouzivat zobrazeni y nahoru a z dopredu
# pokud chceme osy i do zapornych cisel:
rgl.open()                                                
rgl.bg(color = "aliceblue")                               
rgl.viewpoint(theta = 20, phi = 45, zoom = 1)             
rgl.spheres(x,y,z, r = 0.2, color="blueviolet")
rgl.lines(c(-10, 10), c(0, 0), c(0, 0), color = "midnightblue", lwd = 2)           # osa x
rgl.lines(c(0, 0), c(-10,10), c(0, 0), color = "coral4", lwd = 2)                  # osa y
rgl.lines(c(0, 0), c(0, 0), c(-10,10), color = "seagreen", lwd = 2)                # osa z

# parametr par3d, podobne jako ve 2D, nastavime velikost okna
# vyzkousime jinou moznost vykresleni os: pomoci funkce axis3d
rgl.open()                                                
rgl.bg(color = "aliceblue")
par3d(windowRect=c(100,100,600,600))
rgl.viewpoint(theta = 20, phi = 45, zoom = 1)             
rgl.spheres(x,y,z, r = 0.2, color="blueviolet")
axis3d("x", pos=c( NA, 0, 0 ), col = "darkgrey")
axis3d("y", pos=c( 0, NA, 0 ), col = "darkgrey")
axis3d("z", pos=c( 0, 0, NA ), col = "darkgrey")
rgl.close()                                                                        # zavreni aktivnio rgl okna

# pokud se nam nelibi body v podobe kulicek rgl ma moznost sesti zakladnich tvaru
open3d()
shapelist3d(tetrahedron3d(col = "slateblue4"), 0, 0, 0)
shapelist3d(octahedron3d(col = "slateblue3"), 3, 0, 0)
shapelist3d(cube3d(col = "slateblue1"), 6, 0, 0)
shapelist3d(dodecahedron3d(col = "plum2"), 9, 0, 0)
shapelist3d(icosahedron3d(col = "palevioletred3"), 12, 0, 0)
shapelist3d(cuboctahedron3d(col = "salmon"), 15, 0, 0)
rgl.snapshot("typy_bodu.png", fmt = "png")
rgl.postscript("typy_bodu.pdf", fmt = "pdf")
rgl.close() 

# bonusovy point - tvar
shapelist3d(oh3d(col = "aquamarine"),18,0,0)

# pouziti napr. dodecahedronu
rgl.open()                                                
rgl.bg(color = "aliceblue")
par3d(windowRect=c(100,100,600,600))
rgl.viewpoint(theta = 20, phi = 45, zoom = 1)             
shapelist3d(dodecahedron3d(col = "royalblue4"), x, y, z, size = 1.2, alpha = 0.5)   # size = velikost bodu, alpha = pruhlednost
axis3d("x", pos=c( NA, 0, 0 ), col = "darkgrey")
axis3d("y", pos=c( 0, NA, 0 ), col = "darkgrey")
axis3d("z", pos=c( 0, 0, NA ), col = "darkgrey")
rgl.snapshot("body.png", fmt = "png")
rgl.close() 

####### "REGRESNI" ROVINA - prolozeni rovinou

library(rgl)
open3d()
x <- rnorm(100)                                        # nahodne souradnice x
y <- rnorm(100)                                        # nahodne souradnice y
z <- 0.2*x - 0.3*y + rnorm(100, sd = 0.3)              # souradnice z (nahodne ale s jistou zavislosti na x a y)
mod <- lm(z ~ x + y)                                   # regresni model (z v zavislosti na x a y)
plot3d(x, y, z, type = "s", col = "tomato", size = 1)
coefs <- coef(mod)
# funkce coef je genericka funkce, ktera extrahuje koeficienty modelu z objektu, ktere vraci modelove funkce
# zde vybere koeficienty z objektu mod, ktery vygenerovala modelova funkce lm (linear models)
# podivame se na koeficienty (pozor, nezapominame ze, obecna rovnice roviny je: ax + by + cz + d = 0)
# tzn. -cz = ax + by + d, tedy c = -1 v tomto pripade
coefs
# (Intercept)           x           y 
# -0.03352896  0.26614672 -0.30761018 
# tzn.:
a <- coefs["x"]
b <- coefs["y"]
c <- -1
d <- coefs["(Intercept)"]
# a pomoci funkce plane3d rovinu zobrazime
planes3d(a, b, c, d, alpha = 0.5)
# pro zobrazeni roviny lze take pouzit funkci
rgl.planes(a, b, c, d, alpha = 0.5)

# pokud vezmeme nase puvodni "nahodna data" - predpokladali jsme moznost prolozeni a zopakujeme postup
x <- c(3,5,7,9,8,2,1,1,1,4)
y <- c(5,4,1,1,0,6,9,5,7,2)
z <- c(2,1,2,0,2,1,0,4,2,4)
mod <- lm(z ~ x + y)                                   # regresni model
plot3d(x, y, z, type = "s", col = "tomato", size = 1)
coefs <- coef(mod)
coefs
a <- coefs["x"]
b <- coefs["y"]
c <- -1
d <- coefs["(Intercept)"]
# planes3d(a, b, c, d, alpha = 0.5)
rgl.planes(a, b, c, d, alpha = 0.5, col = "tomato")

# pokud nechceme pocitat vse rucne, moznost c.2 napr. pomoci balicku car (pozor balicek car je trochu vetsi)
# install.packages("car")
library(rgl)
library(car)
# nacteme data
DATA <- read.csv("07_data.csv", header = TRUE)
str(DATA)
attach(DATA)
# mame kvantifikovane vzdelani (Education), "sluzebni vek" (Seniority) a prijem (Income)
# chceme se podivat jak zavisi prijem na vzdelani a sluzebnim veku
scatter3d(Income ~ Seniority + Education, data=DATA, fit="smooth")
# pokud bychom chteli rovinu
scatter3d(Income ~ Seniority + Education, data=DATA)
# muzeme odstranit mrizku
scatter3d(Income ~ Seniority + Education, data=DATA, fit="smooth", grid = F, point.col = "purple")

# v pripade pouziti na data s kategorickou promenou
# note: vyuzijeme implementovana data iris (obsahujici rozmery kvetu tri druhu rostlin)
data("iris")
head(iris)
str(iris)
# zde nezadavame formuli, pouze souradnice x, y a z
scatter3d(x = iris$Sepal.Length, y = iris$Petal.Length, z = iris$Sepal.Width)
# mame ale kategoricka data
scatter3d(x = iris$Sepal.Length, y = iris$Petal.Length, z = iris$Sepal.Width, 
          groups = iris$Species, fit = "smooth")
# muzeme upravit barvy
scatter3d(x = iris$Sepal.Length, y = iris$Petal.Length, z = iris$Sepal.Width, 
          groups = iris$Species, fit = "smooth", surface.col = c("#4a148c", "#607d8b", "#00897b"))
# muzeme odstranit roviny a nechat jenom body
scatter3d(x = iris$Sepal.Length, y = iris$Petal.Length, z = iris$Sepal.Width, 
          groups = iris$Species, fit = "smooth", surface = FALSE)

# dokumentace k: https://cran.r-project.org/web/packages/scatterplot3d/vignettes/s3d.pdf

################## ZOBRAZENI ANTROPOLOGICKYCH DAT VE 3D
DATA3D <- read.delim("07_data3d.txt", header = TRUE)
attach(DATA3D)
str(DATA3D)
dim(DATA3D)

DATA3D_f <- DATA3D[sex == "f",]
DATA3D_m <- DATA3D[sex == "m",]

# vidime, ze mame souradnice dvanacti bodu na lebce pro 58 jedincu
# nejprve zobrazime glabellu (G) - a idealne zvlast pro zeny a zvlast pro muze
library(rgl)
library(car)
# rgl.open()
open3d()
rgl.spheres(DATA3D_f$x.G,DATA3D_f$y.G,DATA3D_f$z.G, data = DATA3D_f, r = 2, color="palevioletred1")
rgl.spheres(DATA3D_m$x.G,DATA3D_m$y.G,DATA3D_m$z.G, data = DATA3D_m, r = 2, color="lightskyblue")
rgl.bbox(color="darkblue")

# pro zobrazeni muzeme vyuzit i scatter3d (na grupovana data)
scatter3d(x = DATA3D$x.G, y = DATA3D$y.G, z = DATA3D$z.G, data=DATA3D, groups = DATA3D$sex)
scatter3d(x = DATA3D$x.G, y = DATA3D$y.G, z = DATA3D$z.G, data=DATA3D, groups = DATA3D$sex, fit="smooth")

# ale mame mracno bodu, bylo by vhodne najit prumer tohoto mracna
# zvlast pro muze a zeny!
prumF_x <- mean(DATA3D_f$x.G)
prumF_y <- mean(DATA3D_f$y.G)
prumF_z <- mean(DATA3D_f$z.G)

prumM_x <- mean(DATA3D_m$x.G)
prumM_y <- mean(DATA3D_m$y.G)
prumM_z <- mean(DATA3D_m$z.G)

# rgl.open()
open3d()
rgl.spheres(DATA3D_f$x.G,DATA3D_f$y.G,DATA3D_f$z.G, data = DATA3D_f, r = 2, color="palevioletred4")
rgl.spheres(DATA3D_m$x.G,DATA3D_m$y.G,DATA3D_m$z.G, data = DATA3D_m, r = 2, color="royalblue")
rgl.spheres(prumF_x,prumF_y,prumF_z, r = 4, color="violetred4")
rgl.spheres(prumM_x,prumM_y,prumM_z, r = 4, color="mediumblue")

# chceme zobrazit vsechny body
# nachystame si souradnice do objektu
DATA3D_x <- DATA3D[,c(1,2,3,6,9,12,15,18,21,24,27,30,33,36)]
DATA3D_y <- DATA3D[,c(1,2,4,7,10,13,16,19,22,25,28,31,34,37)]
DATA3D_z <- DATA3D[,c(1,2,5,8,11,14,17,20,23,26,29,32,35,38)]

# funkce rgl.spheres ale vyzaduje zadavani souradnic jako vektory 
# nacteme knihovnu balicku reshape a pomoci funkce melt "spojime" pod sebe vsechny sloupce danych objektu
# install.packages("reshape")
library(reshape)
x_s <- melt(DATA3D_x, id.vars= c("id", "sex"),var = "point")
y_s <- melt(DATA3D_y, id.vars= c("id", "sex"),var = "point")
z_s <- melt(DATA3D_z, id.vars= c("id", "sex"),var = "point")

# zobrazime vsechny body
# rgl.open()
open3d()
rgl.spheres(x_s[sex == "f","value"],y_s[sex == "f","value"],z_s[sex == "f","value"], r = 2, color="palevioletred4")
rgl.spheres(x_s[sex == "m","value"],y_s[sex == "m","value"],z_s[sex == "m","value"], r = 2, color="royalblue")

# jedna se o pomerne neprehledne pole bodu, vypocitame souradnice prumernych bodu
# zvlast pro muze a zeny
prumerF_x <- colMeans(DATA3D[sex == "f",c(3,6,9,12,15,18,21,24,27,30,33,36)])
prumerF_y <- colMeans(DATA3D[sex == "f",c(4,7,10,13,16,19,22,25,28,31,34,37)])
prumerF_z <- colMeans(DATA3D[sex == "f",c(5,8,11,14,17,20,23,26,29,32,35,38)])

prumerM_x <- colMeans(DATA3D[sex == "m",c(3,6,9,12,15,18,21,24,27,30,33,36)])
prumerM_y <- colMeans(DATA3D[sex == "m",c(4,7,10,13,16,19,22,25,28,31,34,37)])
prumerM_z <- colMeans(DATA3D[sex == "m",c(5,8,11,14,17,20,23,26,29,32,35,38)])

# rgl.open()
open3d()
rgl.spheres(x_s[sex == "f","value"],y_s[sex == "f","value"],z_s[sex == "f","value"], r = 2, color="palevioletred4")
rgl.spheres(x_s[sex == "m","value"],y_s[sex == "m","value"],z_s[sex == "m","value"], r = 2, color="royalblue")
rgl.spheres(prumerF_x,prumerF_y,prumerF_z, r = 4, color="violetred4")
rgl.spheres(prumerM_x,prumerM_y,prumerM_z, r = 4, color="mediumblue")

# muzeme zobrazit pouze prumery:
# rgl.open()
open3d()
rgl.viewpoint(theta = 0, phi = -90)
rgl.spheres(prumerF_x,prumerF_y,prumerF_z, r = 4, color="violetred4")
rgl.spheres(prumerM_x,prumerM_y,prumerM_z, r = 4, color="mediumblue")
# chceme k bodum priradit nazvy bodu
colnames(DATA3D)
# vidime, ze se jedna o body: "G","P","FMOdx","FMOsin","JUdx","JUsin","B","L","OP","BA","POdx","POsin"
rgl.texts(prumerF_x,prumerF_y,prumerF_z, text = c("G","P","FMOdx","FMOsin","JUdx","JUsin","B","L","OP","BA","POdx","POsin"), 
          adj = 2.5, color = "gray15")
# chceme spojit bod G s bodem B
rgl.lines(prumerM_x[c(1,7)],prumerM_y[c(1,7)],prumerM_z[c(1,7)], color="mediumblue", lwd = 2)
rgl.lines(prumerM_x[c(7,8)],prumerM_y[c(7,8)],prumerM_z[c(7,8)], color="mediumblue", lwd = 2)

# Priklad: souradnice bodu glabella - pridame elipsoidy pro 95% interval spolehlivosti
rgl.open()
par3d(windowRect=c(20, 30, 800, 800))
rgl.viewpoint(theta = 20, phi = 45, zoom = 1) 
rgl.spheres(DATA3D_f$x.G,DATA3D_f$y.G,DATA3D_f$z.G, data = DATA3D_f, r = 2, color="plum2")
rgl.spheres(DATA3D_m$x.G,DATA3D_m$y.G,DATA3D_m$z.G, data = DATA3D_m, r = 2, color="paleturquoise2")
rgl.spheres(prumF_x,prumF_y,prumF_z, r = 4, color="palevioletred1")
rgl.spheres(prumM_x,prumM_y,prumM_z, r = 4, color="paleturquoise4")
# vypocet elipsoidu, funkce elipse3d(): 
# parametr x = korelacni nebo kovariancni matice mezi x, y a z
# parametr centre = pozice stredu elipsoidu - v nasem pripade v bode o souradnicich: prumM_x,prumM_y,prumM_z (prumery)
# parametr level = interval spolehlivost
# korelacni elipsu3d vykreslujeme pomoci finkce wire3d() nebo napr shade3d()

ellips_F <- ellipse3d(cov(cbind(DATA3D_f$x.G,DATA3D_f$y.G,DATA3D_f$z.G)), centre=c(prumF_x,prumF_y,prumF_z), level = 0.95)
wire3d(ellips_F, col = "plum2",  lit = FALSE)
ellips_M <- ellipse3d(cov(cbind(DATA3D_m$x.G,DATA3D_m$y.G,DATA3D_m$z.G)), centre=c(prumM_x,prumM_y,prumM_z), level = 0.95)
wire3d(ellips_M, col = "paleturquoise2",  lit = FALSE)
# chceme ulozit nahled jako png
rgl.snapshot("G_ellipse3d.png",fmt = "png")
# chceme ulozit nahled jako pdf (note - obcas komplikace v pripade polo-pruhlednych objektu)
# v pripade komplikovanejsich objektu trva prikaz dlouho (zejm. na pomalejsich pocitacich!)
rgl.postscript("G_ellipse3d.pdf", fmt = "pdf")

# vidime jeden bod u glabelly nespada do 95% intervalu => chceme zjistit ktery, pridame text. ID
rgl.open()
par3d(windowRect=c(20, 30, 800, 800))
rgl.viewpoint(theta = 20, phi = 45, zoom = 1) 
rgl.spheres(DATA3D_f$x.G,DATA3D_f$y.G,DATA3D_f$z.G, data = DATA3D_f, r = 2, color="plum2")
rgl.spheres(DATA3D_m$x.G,DATA3D_m$y.G,DATA3D_m$z.G, data = DATA3D_m, r = 2, color="paleturquoise2")
rgl.spheres(prumF_x,prumF_y,prumF_z, r = 4, color="palevioletred1")
rgl.spheres(prumM_x,prumM_y,prumM_z, r = 4, color="paleturquoise4")
rgl.texts(DATA3D_m$x.G,DATA3D_m$y.G,DATA3D_m$z.G, text = DATA3D_m$id, adj = 2.5, color = "ivory")
# vidime, ze se jedna o pripad 24, otazkou je, zda jsou spatne superponovana data, zda se jedna o chybu mereni, extremni hodnotu,
# patologickeho jedince nebo jedince z jine populace

##### DOMACI UKOL 9 ###################################
