rgl Demos

2023-07-04

Introduction

This vignette holds code that was previously included as “demos” in rgl. Some of the demos require R to be running; those remain available via demo(package = "rgl").

hist3d: 3D histogram using basic building blocks

##########
### 3D HIST EXAMPLE:
##########

################################################################################

##### Required functions 'binplot' and 'hist3d':

binplot.3d<-function(x,y,z,alpha=1,topcol="#ff0000",sidecol="#aaaaaa") {
  save <- par3d(skipRedraw=TRUE)
  on.exit(par3d(save))
    
  x1<-c(rep(c(x[1],x[2],x[2],x[1]),3),rep(x[1],4),rep(x[2],4))
  z1<-c(rep(0,4),rep(c(0,0,z,z),4))
  y1<-c(y[1],y[1],y[2],y[2],rep(y[1],4),rep(y[2],4),rep(c(y[1],y[2],y[2],y[1]),2))
  x2<-c(rep(c(x[1],x[1],x[2],x[2]),2),rep(c(x[1],x[2],rep(x[1],3),rep(x[2],3)),2))
  z2<-c(rep(c(0,z),4),rep(0,8),rep(z,8) )
  y2<-c(rep(y[1],4),rep(y[2],4),rep(c(rep(y[1],3),rep(y[2],3),y[1],y[2]),2) )
  quads3d(x1,z1,y1,col=rep(sidecol,each=4),alpha=alpha)
  quads3d(c(x[1],x[2],x[2],x[1]),rep(z,4),c(y[1],y[1],y[2],y[2]),
              col=rep(topcol,each=4),alpha=1) 
  segments3d(x2,z2,y2,col="#000000")
}

hist3d<-function(x,y=NULL,nclass="auto",alpha=1,col="#ff0000",scale=10) {
  save <- par3d(skipRedraw=TRUE)
  on.exit(par3d(save))
  xy <- xy.coords(x,y)
  x <- xy$x
  y <- xy$y
  n<-length(x)
  if (nclass == "auto") nclass<-ceiling(sqrt(nclass.Sturges(x)))
  breaks.x <- seq(min(x),max(x),length=(nclass+1))
  breaks.y <- seq(min(y),max(y),length=(nclass+1))
  z<-matrix(0,(nclass),(nclass))
  for (i in seq_len(nclass)) {
    for (j in seq_len(nclass)) {
      z[i, j] <- (1/n)*sum(x < breaks.x[i+1] & y < breaks.y[j+1] & 
                            x >= breaks.x[i] & y >= breaks.y[j])
      binplot.3d(c(breaks.x[i],breaks.x[i+1]),c(breaks.y[j],breaks.y[j+1]),
                 scale*z[i,j],alpha=alpha,topcol=col)
    }
  }
}
################################################################################

open3d()
#> null 
#>   15
bg3d(color="gray")
light3d(0, 0)

# Drawing a 'bin' for given coordinates:
binplot.3d(c(-0.5,0.5),c(4.5,5.5),2,alpha=0.6)

# Setting the viewpoint ('theta' and 'phi' have the same meaning as in persp):
view3d(theta=40,phi=40)

##### QUADS FORMING BIN

open3d()
#> null 
#>   16
# Defining transparency and colors:
alpha<-0.7; topcol<-"#ff0000"; sidecol<-"#aaaaaa"

# Setting up coordinates for the quads and adding them to the scene:
y<-x<-c(-1,1) ; z<-4; of<-0.3
x12<-c(x[1],x[2],x[2],x[1]); x11<-rep(x[1],4); x22<-rep(x[2],4)
z00<-rep(0,4); z0z<-c(0,0,z,z); zzz<-rep(z,4); y11<-rep(y[1],4)
y1122<-c(y[1],y[1],y[2],y[2]); y12<-c(y[1],y[2],y[2],y[1]); y22<-rep(y[2],4)

quads3d(c(x12,x12,x11-of,x12,x22+of,x12),
          c(z00-of,rep(z0z,4),zzz+of),
          c(y1122,y11-of,y12,y22+of,y12,y1122),
          col=rep(c(rep(sidecol,5),topcol),each=4),alpha=c(rep(alpha,5),1))

# Setting up coordinates for the border-lines of the quads and drawing them:
yl1<-c(y[1],y[2],y[1],y[2]); yl2<-c(y[1]-of,y[1]-of)
xl<-c(rep(x[1],8),rep(x[1]-of,8),rep(c(x[1],x[2]),8),rep(x[2],8),rep(x[2]+of,8))
zl<-c(0,z,0,z,z+of,z+of,-of,-of,0,0,z,z,0,z,0,z,rep(0,4),rep(z,4),rep(-of,4),
      rep(z+of,4),z+of,z+of,-of,-of,rep(c(0,z),4),0,0,z,z)
yl<-c(yl2,y[2]+of,y[2]+of,rep(c(y[1],y[2]),4),y[1],y[1],y[2],y[2],yl2,
      rep(y[2]+of,4),yl2,y[2],y[2],rep(y[1],4),y[2],y[2],yl1,yl2,y[2]+of,
      y[2]+of,y[1],y[1],y[2],y[2],yl1)

lines3d(xl,zl,yl,col="#000000")
view3d(theta=40,phi=40)

##### COMPLETE HISTOGRAM:

open3d()
#> null 
#>   17
# Setting the rng to a fixed value:
set.seed(1000)

# Drawing a 3d histogramm of 2500 normaly distributed observations:
hist3d(rnorm(2500),rnorm(2500),alpha=0.4,nclass=7,scale=30)
# Choosing a lightgrey background:
bg3d(col="#cccccc")
view3d(theta=40,phi=40)

bivar: Bivariate densities: kernel smoothing using surface3d and alpha-channel (requires MASS package)

# rgl demo: rgl-bivar.r
# author: Daniel Adler

rgl.demo.bivar <- function() {
  if (!requireNamespace("MASS", quietly = TRUE))
    stop("This demo requires MASS")
  
  # parameters:
  n<-50; ngrid<-40
  
  # generate samples:
  set.seed(31415)
  x<-rnorm(n); y<-rnorm(n)
  
  # estimate non-parameteric density surface via kernel smoothing
  denobj <- MASS::kde2d(x, y, n=ngrid)
  den.z <-denobj$z
  
  # generate parametric density surface of a bivariate normal distribution
  xgrid <- denobj$x
  ygrid <- denobj$y
  bi.z <- dnorm(xgrid)%*%t(dnorm(ygrid))
  
  # visualize:
  zscale<-20
  
  # New window
  open3d()
  
  # clear scene:
  clear3d("all")
  
  # setup env:
  bg3d(color="#887777")
  light3d()
  
  # Draws the simulated data as spheres on the baseline
  spheres3d(x,y,rep(0,n),radius=0.1,color="#CCCCFF")
  
  # Draws non-parametric density
  surface3d(xgrid,ygrid,den.z*zscale,color="#FF2222",alpha=0.5)
  
  # Draws parametric density
  surface3d(xgrid,ygrid,bi.z*zscale,color="#CCCCFF",front="lines") 
}

rgl.demo.bivar()

Abundance: Animal abundance, visualization of multi-dimension data using multiple techniques

# RGL-Demo: animal abundance
# Authors: Oleg Nenadic, Daniel Adler

rgl.demo.abundance <- function() {
  open3d()
  clear3d("all")               # remove all shapes, lights, bounding-box, and restore viewpoint
  
  # Setup environment:
  bg3d(col="#cccccc")     # setup background
  light3d()               # setup head-light
  
  # Importing animal data (created with wisp)
  terrain<-dget(system.file("demodata/region.dat",package="rgl"))
  pop<-dget(system.file("demodata/population.dat",package="rgl"))
  
  # Define colors for terrain
  zlim <- range(terrain)
  colorlut <- terrain.colors(82) 
  col1 <- colorlut[9*sqrt(3.6*(terrain-zlim[1])+2)]
  
  # Set color to (water-)blue for regions with zero 'altitude' 
  col1[terrain==0]<-"#0000FF"
  
  # Add terrain surface shape (i.e. population density):
  surface3d( 
      1:100,seq(1,60,length=100),terrain,
      col=col1,spec="#000000", ambient="#333333", back="lines"
  )
  
  # Define colors for simulated populations (males:blue, females:red):
  col2<-pop[,4]
  col2[col2==0]<-"#3333ff"
  col2[col2==1]<-"#ff3333"
  
  # Add simulated populations as sphere-set shape
  spheres3d(
    pop[,1],
    pop[,2],
    terrain[cbind( ceiling(pop[,1]),ceiling(pop[,2]*10/6) )]+0.5,
    radius=0.2*pop[,3], col=col2, alpha=(1-(pop[,5])/10 )
  )

}
rgl.demo.abundance()

lsystem: Plant modelling using a turtle and L-system

# demo: lsystem.r
# author: Daniel Adler

#
# geometry 
#

deg2rad <- function( degree ) {
  return( degree*pi/180 )
}

rotZ.m3x3 <- function( degree ) {
  kc <- cos(deg2rad(degree))
  ks <- sin(deg2rad(degree))
  return( 
    matrix(
      c(
        kc, -ks,   0, 
        ks,  kc,   0,
         0,   0,   1
      ),ncol=3,byrow=TRUE
    ) 
  )
}

rotX.m3x3 <- function( degree ) {
  kc <- cos(deg2rad(degree))
  ks <- sin(deg2rad(degree))
  return(
    matrix(
      c(
         1,   0,   0,
         0,  kc, -ks,
         0,  ks,  kc
      ),ncol=3,byrow=TRUE
    )
  )
}

rotY.m3x3 <- function( degree ) {
  kc <- cos(deg2rad(degree))
  ks <- sin(deg2rad(degree))
  return(
    matrix(
      c(
        kc,   0,   ks,
         0,   1,    0,
       -ks,   0,   kc
      ),ncol=3,byrow=TRUE
    )
  )
}

rotZ <- function( v, degree ) {
  return( rotZ.m3x3(degree) %*% v)
}

rotX <- function( v, degree ) {
  return( rotX.m3x3(degree) %*% v)
}

rotY <- function( v, degree ) {
  return( rotY.m3x3(degree) %*% v)
}


#
# turtle graphics, rgl implementation:
#

turtle.init <- function(pos=c(0,0,0),head=0,pitch=90,roll=0,level=0) {
  clear3d("all")
  bg3d(color="gray")
  light3d()
  return( list(pos=pos,head=head,pitch=pitch,roll=roll,level=level) )
}


turtle.move <- function(turtle, steps, color) {
  
  rm <- rotX.m3x3(turtle$pitch) %*% rotY.m3x3(turtle$head) %*% rotZ.m3x3(turtle$roll)
  
  from <- as.vector( turtle$pos )  
  dir  <- rm %*% c(0,0,-1)
  to   <- from + dir * steps
    
  x <- c( from[1], to[1] )
  y <- c( from[2], to[2] )
  z <- c( from[3], to[3] )
  lines3d(x,y,z,col=color,size=1.5,alpha=0.5)
  turtle$pos <- to
  return(turtle)
}

turtle.pitch <- function(turtle, degree) {
  turtle$pitch <- turtle$pitch + degree
  return(turtle)
}

turtle.head <- function(turtle, degree) {
  turtle$head <- turtle$head + degree
  return(turtle)
}

turtle.roll <- function(turtle, degree) {
  turtle$roll <- turtle$roll + degree
  return(turtle)
}

#
# l-system general
#


lsystem.code <- function( x )
  substitute( x )

lsystem.gen <- function( x, grammar, levels=0 ) {
  code <- eval( substitute( substitute( REPLACE , grammar ), list(REPLACE=x) ) )
  if (levels)
    return( lsystem.gen( code , grammar , levels-1 ) )
  else
    return( code )
}

#
# l-system plot
#

lsystem.plot <- function( expr, level ) {
  turtle <- turtle.init(level=level)
  lsystem.eval( expr, turtle )
}

lsystem.eval <- function( expr, turtle ) {
  if ( length(expr) == 3 ) {
    turtle <- lsystem.eval( expr[[2]], turtle )
    turtle <- lsystem.eval( expr[[3]], turtle )
    turtle <- lsystem.eval( expr[[1]], turtle )
  } else if ( length(expr) == 2 ) {
    saved <- turtle
    turtle <- lsystem.eval( expr[[1]], turtle )
    turtle <- lsystem.eval( expr[[2]], turtle )
    turtle <- saved
  } else if ( length(expr) == 1 ) {
    if ( as.name(expr) == "stem" )      turtle <- turtle.move(turtle, 5, "brown")
    else if ( as.name(expr) == "short") turtle <- turtle.move(turtle, 5, "brown")
    else if ( as.name(expr) == "leaf" ) {
      spheres3d(turtle$pos[1],turtle$pos[2],turtle$pos[3],radius=0.1+turtle$level*0.3,color="green")
      sprites3d(turtle$pos[1],turtle$pos[2],turtle$pos[3],radius=0.5+turtle$level*0.3 ,color="green",texture=system.file("textures/particle.png",package="rgl"),textype="alpha",alpha=0.5)     
    }
    else if ( as.name(expr) == "roll" ) turtle <- turtle.head(turtle, 60)
    else if ( as.name(expr) == "down" ) turtle <- turtle.pitch(turtle,10)
    else if ( as.name(expr) == "up" )   turtle <- turtle.pitch(turtle,-10)
    else if ( as.name(expr) == "left" ) turtle <- turtle.head(turtle, 1)
    else if ( as.name(expr) == "right") turtle <- turtle.head(turtle,-1.5)
    else if ( as.name(expr) == "turnleft") turtle <- turtle.head(turtle,20)
    else if ( as.name(expr) == "turnright") turtle <- turtle.head(turtle,-20)
    else if ( as.name(expr) == "turn") turtle <- turtle.roll(turtle,180)
  }
  return(turtle)
}


#
# example
#

simple <- function(level=0) {
  grammar <- list(
    stem=lsystem.code(
      stem-(up-stem-leaf)-stem-(down-stem-leaf)-stem-leaf
    )
  )
  plant <- lsystem.gen(lsystem.code(stem), grammar, level )
  lsystem.plot(plant,level)
}

rgl.demo.lsystem <- function(level=0) {
  gen   <- list(
    stem=lsystem.code( 
      stem-left-stem-branch( turnleft-down-short-turnleft-down-stem-leaf)-right-right-stem--branch( turnright-up-short-turnright-up-short-turnright-short-stem-leaf)-left-left-left-stem-branch( turnleft-down-short-turnright-down-stem-leaf )-branch( up-turnright-short-up-turnleft-up-stem-leaf ) 
    )
  )    
  plant <- lsystem.gen(lsystem.code(stem), gen, level )
  lsystem.plot(plant,level)  
}

open3d()
#> null 
#>   20
rgl.demo.lsystem(level=1)
#> $pos
#>               [,1]
#> [1,] -2.092747e+00
#> [2,]  7.992083e+01
#> [3,] -4.893740e-15
#> 
#> $head
#> [1] 5
#> 
#> $pitch
#> [1] 90
#> 
#> $roll
#> [1] 0
#> 
#> $level
#> [1] 1
rglwidget()