data(volcano)
z <- 2 * volcano # Exaggerate the relief
x <- 10 * (1:nrow(z)) # 10 meter spacing (S to N)
y <- 10 * (1:ncol(z)) # 10 meter spacing (E to W)
p <- persp3d(x, y, z, scale = FALSE, axes = FALSE,smooth = T, xlab = '', ylab = '', zlab = '')
r <- range(z)
zvect <- (z[p$indz] - r[1])/(r[2]-r[1])
palette <- color.to.rgb(terrain.colors(50))
rgb <- palette[round(zvect*49,0)+1]
setcolors3d(p$handle, rgb = rgb)


# Plot a random sample and an ellipsoid of concentration corresponding to a 95%
# probability region for a
# trivariate normal distribution with mean 0, unit variances and
# correlation 0.8.
if (requireNamespace("MASS")) {
  Sigma <- matrix(c(10,3,0,3,2,0,0,0,1), 3,3)
  Mean <- 1:3
  x <- MASS::mvrnorm(1000, Mean, Sigma)
  open3d()
  plot3d(x, box=FALSE)
  plot3d( ellipse3d(Sigma, centre=Mean), col="green", alpha=0.5, add = TRUE)
}








data(mtcars)
fit <- lm(mpg ~ disp + cyl , mtcars)
open3d()
plot3d(ellipse3d(fit, level = 0.90), col="blue", alpha=0.5, aspect=TRUE)





x <- c(1:10, 10:1)
y <- rev(c(rep(c(0,2), 5), rep(c(1.5,-0.5),5)))
plot(x, y, type="n")
polygon(x, y)
open3d()
shade3d( extrude3d(x, y), col = "red" )




grid3d(side, at = NULL, col = "gray", lwd = 1, lty = 1, n = 5)


x <- 1:10

y <- 1:10

z <- matrix(outer(x-5,y-5) + rnorm(100), 10, 10)

open3d()

persp3d(x, y, z, col="red", alpha=0.7, aspect=c(1,1,0.5))

grid3d(c("x", "y+", "z"))






#

# a lightsource moving through the scene

#

data(volcano)

z <- 2 * volcano # Exaggerate the relief

x <- 10 * (1:nrow(z)) # 10 meter spacing (S to N)

y <- 10 * (1:ncol(z)) # 10 meter spacing (E to W)

zlim <- range(z)

zlen <- zlim[2] - zlim[1] + 1

colorlut <- terrain.colors(zlen) # height color lookup table

col <- colorlut[ z-zlim[1]+1 ] # assign colors to heights for each point

open3d()

bg3d("gray50")

surface3d(x, y, z, color=col, back="lines")

r <- max(y)-mean(y)

lightid <- spheres3d(1,1,1,alpha=0)

for (a in seq(-pi, pi, length.out=100)) {
  
  save <- par3d(skipRedraw=TRUE)
  
  clear3d(type = "lights")
  
  rgl.pop(id = lightid)
  
  xyz <- matrix(c(r*sin(a)+mean(x), r*cos(a) + mean(y), max(z)), ncol=3)
  
  light3d(x = xyz, diffuse = "gray75",
          
          light3d(diffuse="gray10", specular="gray25")
          
          lightid <- spheres3d(xyz, emission="white", radius=4)
          
          par3d(save)
          
          Sys.sleep(0.02)
          
}







# A 90 degree rotation about the x axis:

rotationMatrix(pi/2, 1, 0, 0)

# Find what happens when you rotate (2,0,0) by 45 degrees about the y axis:

x <- asHomogeneous(c(2,0,0))

y <- x

asEuclidean(y)

# or more simply...

rotate3d(c(2,0,0), pi/4, 0, 1, 0)






vertices <- c(
  
  -1.0, -1.0, 0, 1.0,
  
  1.0, -1.0, 0, 1.0,
  
  1.0, 1.0, 0, 1.0,
  
  -1.0, 1.0, 0, 1.0
  
)

indices <- c( 1, 2, 3, 4 )

open3d()

wire3d( qmesh3d(vertices,indices) )

# render 4 meshes vertically in the current view

open3d()

bg3d("gray")

l0 <- oh3d(tran = par3d("userMatrix"), color = "green" )

shade3d( translate3d( l0, -6, 0, 0 ))

l1 <- subdivision3d( l0 )

shade3d( translate3d( l1 , -2, 0, 0 ), color="red", override = FALSE )

l2 <- subdivision3d( l1 )

shade3d( translate3d( l2 , 2, 0, 0 ), color="red", override = TRUE )

l3 <- subdivision3d( l2 )

shade3d( translate3d( l3 , 6, 0, 0 ), color="red" )

# render all of the Platonic solids

open3d()

shade3d( translate3d( tetrahedron3d(col="red"), 0, 0, 0) )

shade3d( translate3d( cube3d(col="green"), 3, 0, 0) )

shade3d( translate3d( octahedron3d(col="blue"), 6, 0, 0) )

shade3d( translate3d( dodecahedron3d(col="cyan"), 9, 0, 0) )

shade3d( translate3d( icosahedron3d(col="magenta"), 12, 0, 0) )












# (1) The Obligatory Mathematical surface.

# Rotated sinc function.

x <- seq(-10, 10, length= 30)

y <- x

f <- function(x,y) { r <- sqrt(x^2+y^2); 10 * sin(r)/r }

z <- outer(x, y, f)

z[is.na(z)] <- 1

open3d()

bg3d("white")

material3d(col="black")

persp3d(x, y, z, aspect=c(1, 1, 0.5), col = "lightblue",
        
        xlab = "X", ylab = "Y", zlab = "Sinc( r )")

# (2) Add to existing persp plot:

xE <- c(-10,10); xy <- expand.grid(xE, xE)

points3d(xy[,1], xy[,2], 6, col = "red")

lines3d(x, y=10, z= 6 + sin(x), col = "green")

phi <- seq(0, 2*pi, len = 201)

r1 <- 7.725 # radius of 2nd maximum

xr <- r1 * cos(phi)

yr <- r1 * sin(phi)

lines3d(xr,yr, f(xr,yr), col = "pink", lwd = 2)

# (3) Visualizing a simple DEM model

z <- 2 * volcano # Exaggerate the relief

x <- 10 * (1:nrow(z)) # 10 meter spacing (S to N)

y <- 10 * (1:ncol(z)) # 10 meter spacing (E to W)

open3d()

bg3d("slategray")

material3d(col="black")

persp3d(x, y, z, col = "green3", aspect="iso",
        
        axes = FALSE, box = FALSE)

# (4) A cylindrical plot

z <- matrix(seq(0, 1, len=50), 50, 50)

theta <- t(z)

r <- 1 + exp( -pmin( (z - theta)^2, (z - theta - 1)^2, (z - theta + 1)^2 )/0.01 )

x <- r*cos(theta*2*pi)

y <- r*sin(theta*2*pi)

open3d()

persp3d(x, y, z, col="red")

# (5) A globe

lat <- matrix(seq(90,-90, len=50)*pi/180, 50, 50, byrow=TRUE)

long <- matrix(seq(-180, 180, len=50)*pi/180, 50, 50)

r <- 6378.1 # radius of Earth in km

x <- r*cos(lat)*cos(long)

y <- r*cos(lat)*sin(long)

z <- r*sin(lat)

open3d()

persp3d(x, y, z, col="white",
        
        texture=system.file("textures/worldsmall.png",package="rgl"),
        
        specular="black", axes=FALSE, box=FALSE, xlab="", ylab="", zlab="",
        
        normal_x=x, normal_y=y, normal_z=z)

if (!rgl.useNULL())
  
  play3d(spin3d(axis=c(0,0,1), rpm=8), duration=5)

## Not run:

# This looks much better, but is slow because the texture is very big

persp3d(x, y, z, col="white",
        
        texture=system.file("textures/world.png",package="rgl"),
        
        specular="black", axes=FALSE, box=FALSE, xlab="", ylab="", zlab="",
        
        normal_x=x, normal_y=y, normal_z=z)

## End(Not run)
















# Show regression plane with z as dependent variable

x <- rnorm(100)

y <- rnorm(100)

z <- 0.2*x - 0.3*y + rnorm(100, sd=0.3)

fit <- lm(z ~ x + y)

plot3d(x,y,z, type="s", col="red", size=1)

coefs <- coef(fit)

a <- coefs["x"]

b <- coefs["y"]

c <- -1

d <- coefs["(Intercept)"]

planes3d(a, b, c, d, alpha=0.5)

open3d()

plot3d(x,y,z, type="s", col="red", size=1)

clipplanes3d(a, b, c, d)








open3d()

plot3d( cube3d(col="green") )

M <- par3d("userMatrix")

if (!rgl.useNULL())
  
  play3d( par3dinterp( userMatrix=list(M,
                                       
                                       duration=4 )
                       
                       ## Not run:
                       
                       movie3d( spin3d(), duration=5 )
                       
                       ## End(Not run)
                       
                       rotate3d(M, pi/2, 1, 0, 0),
                       
                       rotate3d(M, pi/2, 0, 1, 0) ) ),
duration=4 )

## Not run:

movie3d( spin3d(), duration=5 )

## End(Not run)









open3d()

x <- sort(rnorm(1000))

y <- rnorm(1000)

z <- rnorm(1000) + atan2(x,y)

plot3d(x, y, z, col=rainbow(1000))






# Show 12 random vertices in various ways.

M <- matrix(rnorm(36), 3, 12, dimnames=list(c('x','y','z'),
                                            
                                            rep(LETTERS[1:4], 3)))

# Force 4-tuples to be convex in planes so that quads3d works.

for (i in c(1,5,9)) {
  
  quad <- as.data.frame(M[,i+0:3])
  
  coeffs <- runif(2,0,3)
  
  if (mean(coeffs) < 1) coeffs <- coeffs + 1 - mean(coeffs)
  
  quad$C <- with(quad, coeffs[1]*(B-A) + coeffs[2]*(D-A) + A)
  
  M[,i+0:3] <- as.matrix(quad)
  
}

open3d()

# Rows of M are x, y, z coords; transpose to plot

M <- t(M)

shift <- matrix(c(-3,3,0), 12, 3, byrow=TRUE)

points3d(M)

lines3d(M + shift)

segments3d(M + 2*shift)

triangles3d(M + 3*shift, col='red')

quads3d(M + 4*shift, col='green')

text3d(M + 5*shift, texts=1:12)

# Add labels

shift <- outer(0:5, shift[1,])

shift[,1] <- shift[,1] + 3

text3d(shift,
       
       texts = c('points3d','lines3d','segments3d',
                 
                 'triangles3d', 'quads3d','text3d'),
       
       adj = 0)

rgl.bringtotop()








theta <- seq(0, 4*pi, len=50)

r <- theta + 1

r <- c(r[-50], rev(theta*0.8) + 1)

theta <- c(theta[-50], rev(theta))

x <- r*cos(theta)

y <- r*sin(theta)

plot(x,y,type="n")

polygon(x, y)

polygon3d(x, y, x+y, col = "blue")














x <- c(0,1,0,0)

y <- c(0,0,1,0)

z <- c(0,0,0,1)

labels <- c("Origin", "X", "Y", "Z")

i <- c(1,2,1,3,1,4)

# rgl.* interface

rgl.open()

rgl.texts(x,y,z,labels)

rgl.texts(1,1,1,"rgl.* coordinates")

rgl.lines(x[i],y[i],z[i])

# *3d interface

open3d()

text3d(x,y,z,labels)

text3d(1,1,1,"*3d coordinates")

segments3d(x[i],y[i],z[i])








filename <- tempfile(fileext=".stl")

open3d()

shade3d( icosahedron3d(col="magenta") )

writeSTL(filename)

open3d()

readSTL(filename, col="red")






rgl.open()

rgl.points(rnorm(100), rnorm(100), rnorm(100))

rgl.bbox(color=c("#333377","white"), emission="#333377",
         
         specular="#3333FF", shininess=5, alpha=0.8 )

open3d()

points3d(rnorm(100), rnorm(100), rnorm(100))

bbox3d(color=c("#333377","black"), emission="#333377",
       
       specular="#3333FF", shininess=5, alpha=0.8)













save <- material3d("color")

material3d(color="red")

material3d("color")

material3d(color=save)


# this illustrates the effect of depth_test

x <- c(1:3); xmid <- mean(x)

y <- c(2,1,3); ymid <- mean(y)

z <- 1

open3d()

tests <- c("never", "less", "equal", "lequal", "greater", "notequal", "gequal", "always")


for (i in 1:8) {
  
  triangles3d(x,y,z+i, col=heat.colors(8)[i])
  
  texts3d(xmid,ymid,z+i, paste(i, tests[i], sep=". "), depth_test=tests[i])
  
}











x <- y <- seq(-10,10,length=20)

z <- outer(x,y,function(x,y) x^2 + y^2)

persp3d(x,y,z, col='lightblue')

title3d("Using LaTeX text", col='red', line=3)

rgl.postscript("persp3da.ps","ps",drawText=FALSE)

rgl.postscript("persp3da.pdf","pdf",drawText=FALSE)

rgl.postscript("persp3da.tex","tex")

rgl.pop()

title3d("Using ps/pdf text", col='red', line=3)

rgl.postscript("persp3db.ps","ps")

rgl.postscript("persp3db.pdf","pdf")

rgl.postscript("persp3db.tex","tex",drawText=FALSE)

## Not run:

#

# create a series of frames for an animation

#

rgl.open()

shade3d(oh3d(), color="red")

rgl.viewpoint(0,20)

for (i in 1:45) {
  
  rgl.viewpoint(i,20)
  
  filename <- paste("pic",formatC(i,digits=1,flag="0"),".eps",sep="")
  
  rgl.postscript(filename, fmt="eps")
  
}

## End(Not run)










open3d()

## Not quite right --- this doesn't play well with rescaling

pan3d <- function(button) {
  
  start <- list()
  
  begin <- function(x, y) {
    
    start$userMatrix <<- par3d("userMatrix")
    
    start$viewport <<- par3d("viewport")
    
    start$scale <<- par3d("scale")
    
    start$projection <<- rgl.projection()
    
    start$pos <<- rgl.window2user( x/start$viewport[3], 1 - y/start$viewport[4], 0.5,
                                   projection=start$projection)
                                   
  }
  
  update <- function(x, y) {
    
    xlat <- (rgl.window2user( x/start$viewport[3], 1 - y/start$viewport[4], 0.5,
                              
             projection = start$projection) - start$pos)*start$scale
                              
                              mouseMatrix <- translationMatrix(xlat[1], xlat[2], xlat[3])
                              
                              par3d(userMatrix = start$userMatrix %*% t(mouseMatrix) )
                              
  }
  
  rgl.setMouseCallbacks(button, begin, update)
  
  cat("Callbacks set on button", button, "of rgl device",rgl.cur(),"\n")
  
}

pan3d(3)









shade3d(oh3d(), color="red")

rgl.bringtotop()

rgl.viewpoint(0,20)

setwd(tempdir())

for (i in 1:45) {
  
  rgl.viewpoint(i,20)
  
  filename <- paste("pic",formatC(i,digits=1,flag="0"),".png",sep="")
  
  rgl.snapshot(filename)
  
}

## Now run ImageMagick command:

## convert -delay 10 *.png -loop 0 pic.gif

## End(Not run)
















#

# volcano example taken from "persp"

#

data(volcano)

y <- 2 * volcano # Exaggerate the relief

x <- 10 * (1:nrow(y)) # 10 meter spacing (S to N)

z <- 10 * (1:ncol(y)) # 10 meter spacing (E to W)

ylim <- range(y)

ylen <- ylim[2] - ylim[1] + 1

colorlut <- terrain.colors(ylen) # height color lookup table

col <- colorlut[ y-ylim[1]+1 ] # assign colors to heights for each point

rgl.open()

rgl.surface(x, z, y, color=col, back="lines")










open3d()

points3d(rnorm(100), rnorm(100), rnorm(100))

if (interactive() || !.Platform$OS=="unix") {
  
  # Calculate a square in the middle of the display and plot it
  
  square <- rgl.window2user(c(0.25, 0.25, 0.75, 0.75, 0.25),
                            c(0.25, 0.75, 0.75, 0.25, 0.25), 0.5)
                            par3d(ignoreExtent = TRUE)
                            
                            lines3d(square)
                            
                            par3d(ignoreExtent = FALSE)
                            
}















x <- rnorm(100)

y <- rnorm(100)

z <- rnorm(100)

p <- plot3d(x, y, z, type='s')

rgl.ids()

lines3d(x, y, z)

rgl.ids()

if (interactive()) {
  
  readline("Hit enter to change spheres")
  
  rgl.pop(id = p["data"])
  
  spheres3d(x, y, z, col="red", radius=1/5)
  
  box3d()
  
}







open3d()

z <- 2 * volcano # Exaggerate the relief

x <- 10 * (1:nrow(z)) # 10 meter spacing (S to N)

y <- 10 * (1:ncol(z)) # 10 meter spacing (E to W)

persp3d(x, y, z, col = "green3", aspect="iso")









s <- scene3d()

# Make it bigger

s$par3d$windowRect <- 1.5*s$par3d$windowRect

# and draw it again

plot3d(s)






# Allow the user to select some points, and then redraw them

# in a different color

if (interactive()) {
  
  x <- rnorm(1000)
  
  y <- rnorm(1000)
  
  z <- rnorm(1000)
  
  open3d()
  
  points3d(x,y,z)
  
  f <- select3d()
  
  keep <- f(x,y,z)
  
  rgl.pop()
  
  points3d(x[keep],y[keep],z[keep],color='red')
  
  points3d(x[!keep],y[!keep],z[!keep])
  
}










ids <- plot3d( rnorm(20), rnorm(20), rnorm(20) )

# Click near a point to select it and put a sphere there

selectpoints3d(ids["data"],
               
               multiple = function(x) {
                 
                 spheres3d(x, color = "red", alpha = 0.3, radius = 0.2)
                 
                 TRUE
                 
               })











shapelist3d(icosahedron3d(), x=rnorm(10), y=rnorm(10), z=rnorm(10), col=1:5, size=0.3)







open3d()

spheres3d(rnorm(10), rnorm(10), rnorm(10), radius=runif(10), color=rainbow(10))










# Spin one object

open3d()

plot3d(oh3d(col="lightblue", alpha=0.5))

if (!rgl.useNULL())
  
  play3d(spin3d(axis=c(0,0,1), rpm=10), duration=10)








open3d()

spriteid <- NULL

spin1 <- spin3d(rpm=3 ) # the scene spinner

spin2 <- spin3d(rpm=6 ) # the sprite spinner

f <- function(time) {
  
  par3d(skipRedraw = TRUE) # stops intermediate redraws
  
  on.exit(par3d(skipRedraw=FALSE)) # redraw at the end
  
  rgl.pop(id=spriteid) # delete the old sprite
  
  cubeid <- shade3d(cube3d(), col="red")
  
  spriteid <<- sprites3d(0:1, 0:1, 0:1, shape=cubeid,
                         userMatrix=spin2(time,
                         base=spin1(time)$userMatrix)$userMatrix)
                         spin1(time)
                         
}

if (!rgl.useNULL())
  play3d(f, duration=3)


                 











open3d()

particles3d( rnorm(100), rnorm(100), rnorm(100), color=rainbow(100) )

# is the same as
sprites3d( rnorm(100), rnorm(100), rnorm(100), color=rainbow(100),
           lit=FALSE, alpha=.2,
           textype="alpha", texture=system.file("textures/particle.png", package="rgl") )

sprites3d( rnorm(10)+6, rnorm(10), rnorm(10), shape=shade3d(tetrahedron3d(), col="red") )











open3d()

shade3d( subdivision3d( cube3d(), depth=3 ), color="red", alpha=0.5 )














#

# volcano example taken from "persp"

#

data(volcano)

z <- 2 * volcano # Exaggerate the relief

x <- 10 * (1:nrow(z)) # 10 meter spacing (S to N)

y <- 10 * (1:ncol(z)) # 10 meter spacing (E to W)

zlim <- range(y)

zlen <- zlim[2] - zlim[1] + 1

colorlut <- terrain.colors(zlen) # height color lookup table

col <- colorlut[ z-zlim[1]+1 ] # assign colors to heights for each point

open3d()

surface3d(x, y, z, color=col, back="lines")














open3d()

famnum <- rep(1:4, 8)

family <- c("serif", "sans", "mono", "symbol")[famnum]

font <- rep(rep(1:4, each=4), 2)

cex <- rep(1:2, each=16)

text3d(font, cex, famnum, text=paste(family, font),adj = 0.5,
       
       color="blue", family=family, font=font, cex=cex)

# Not run:

# These FreeType fonts are available from the Amaya project, and are not shipped

# with rgl. You would normally install them to the rgl/fonts directory

# and use fully qualified pathnames, e.g.

# system.file("fonts/FreeSerif.ttf", package= "rgl")


rglFonts(serif=c("FreeSerif.ttf","FreeSerifBold.ttf","FreeSerifItalic.ttf",
                 
                 "FreeSerifBoldItalic.ttf"),
         
         sans =c("FreeSans.ttf", "FreeSansBold.ttf", "FreeSansOblique.ttf",
                 
                 "FreeSansBoldOblique.ttf"),
         
         mono =c("FreeMono.ttf", "FreeMonoBold.ttf", "FreeMonoOblique.ttf",
                 
                 "FreeMonoBoldOblique.ttf"),
         
         symbol=c("ESSTIX10.TTF", "ESSTIX12.TTF", "ESSTIX9_.TTF",
                  
                  "ESSTIX11.TTF"))

## End(Not run)
         
         
         














x <- 1:10

y <- rnorm(10)^2

shade3d(turn3d(x, y), col = "green")










rgl.open()

shade3d(oh3d(), color="red")

start <- proc.time()[3]

while ((i <- 36*(proc.time()[3]-start)) < 360) {
  
  rgl.viewpoint(i,i/4);
  
}








filename <- tempfile(fileext=".obj")

open3d()

shade3d( icosahedron3d() )

writeOBJ(filename)







filename <- tempfile(fileext=".ply")

open3d()

shade3d( icosahedron3d(col="magenta") )

writePLY(filename)











plot3d(rnorm(100), rnorm(100), rnorm(100), type="s", col="red")
# This writes a copy into temporary directory 'webGL', and then displays it
browseURL(paste("file://", writeWebGL(dir=file.path(tempdir(), "webGL"), width=500), sep=""))








library(rgl)
readSTL("CONEX.STL",plot=TRUE,col="red")
rgl.bringtotop(stay = TRUE)


browseURL(paste("file://", writeWebGL(dir=file.path(tempdir(), "webGL"), width=500), sep=""))