R_christmasTree.R

rm(list=ls())
setwd("image/")

library(plot3D)

trunk <- list(
  x0=0, y0=0, z0=0,
  x1=0, y1=0, z1=6,
  extension=NULL,
  branch1=NULL, branch2=NULL,
  branch3=NULL, branch4=NULL,
  depth=1,
  lwd=1,
  stub_color="brown"
)

draw_tree <- function(tree) {
  #function to draw the tree recursively
  if (is.null(tree)) return()
  with(tree,{
    segments3D(x0,y0,z0,x1,y1,z1,col=stub_color,
               lwd=lwd,add=TRUE)
    draw_tree(branch1)
    draw_tree(branch2)
    draw_tree(branch3)
    draw_tree(extension)
  })
}

extend <- function(tree) {
  #creates an extension of the main branch
  #in the same direction vector u=p1-p0
  with(tree,{
    growth_factor <- runif(1,min=0.9,max=1)
    return(list(x0=x1,y0=y1,z0=z1,
                x1=x1+growth_factor*(x1-x0),
                y1=y1+growth_factor*(y1-y0),
                z1=z1+growth_factor*(z1-z0),
                extension=NULL,
                branch1=NULL, branch2= NULL,
                branch3=NULL, branch4= NULL,
                depth=depth*0.9,
                lwd=1,
                stub_color="darkgreen"))
  })
}

create_branch <- function(tree,c1,c2) {
  #this function returns a branch for tree
  #c1 and c2 are the components for the vector
  #base {v1,v2}. v1 and v2 are two orthonrmal
  #vectors, both perpedicular to the direction
  #vector u=p1-p0 (p1 and p0 are the endpoints
  #of the tree trunk).
  #where_at is a scalar value in (0,1), indicating
  #where along p0=(x0,y0,z0) and p1=(x1,y1,z1)
  #the branch will grow.
  with(tree,{
    #vector u=(x1-x0,y1-y0,z1-z0) gives the 
    #direction of the tree trunk.
    #stub_length is the length of that trunk.
    u <- c(x1-x0,y1-y0,z1-z0)
    stub_length <- sqrt(sum(u*u))
    #growth_factor is how long the branch
    #will be with respect to stub_length
    growth_factor <- runif(1,min=0.7,max=0.8)
    where_at <- runif(1,min=0.4,max=0.8)
    #create two perpendicular vectors to u
    #and set their norms to 1
    i <- which(u!=0)[1]
    oi <- setdiff(1:3,i)
    v1 <- rep(0,3); v1[i] <- -u[oi[1]]; v1[oi[1]] <- u[i]
    v2 <- rep(0,3); v2[i] <- -u[oi[2]]; v2[oi[2]] <- u[i]
    v1 <- v1/sqrt(sum(v1*v1)); v2 <- v2/sqrt(sum(v2*v2))
    #vector v is a linear combination of v1 and v2,
    #hence it lies on the plane perpendicular to u
    v <- v1*c1+v2*c2
    #beta is the angle of separation between the
    #branch and the trunk. The angle is larger
    #for the "deeper" parts of the tree.
    beta <- runif(1,min=0.7*depth,max=0.9*depth)    
    #calculate the scalar k for multiplying with v
    #so that u+kv, the branch, will form a beta
    #angle with the trunk.
    k <- tan(beta)*sqrt(sum(u*u))/sqrt(sum(v*v))
    #new_u is the direction vector of the branch,
    #whose length will be the length of the branch 
    new_u <- u+k*v
    new_u <- new_u / sqrt(sum(new_u*new_u))*
      stub_length*growth_factor
    #the new (x0,y0,z0) point of the branch
    #is where along the main trunk the branch begins.
    new_x0 <- x0+(x1-x0)*where_at
    new_y0 <- y0+(y1-y0)*where_at
    new_z0 <- z0+(z1-z0)*where_at
    #the new (x1,y1,z1) point of the branch is
    #its ending point.
    new_x1 <- new_x0+new_u[1]
    new_y1 <- new_y0+new_u[2]
    new_z1 <- new_z0+new_u[3]
    list(
      x0=new_x0,y0=new_y0,z0=new_z0,
      x1=new_x1,y1=new_y1,z1=new_z1,
      extension=NULL,
      branch1=NULL, branch2= NULL,
      branch3=NULL, branch4= NULL,
      depth=depth*0.8,
      lwd=1,
      stub_color="darkgreen")    
  })
}

grow_tree <- function(tree) {
  if (is.null(tree) ) return(NULL)
  tree$lwd <- tree$lwd*1.4
  if (tree$lwd>2.0) tree$stub_color <- 'brown4'
  if (is.null(tree$extension)) {
    tree$extension <- extend(tree)
    if (tree$depth<0.5) {
      #make two branches at an angle of more than 160 degrees
      rot1 <- runif(1,min=0,max=2*pi)
      rot2 <- (runif(1,min=8*pi/9,max=10*pi/9)+rot1)%%(2*pi)
      tree$branch1 <- create_branch(tree,cos(rot1),sin(rot1))
      tree$branch2 <- create_branch(tree,cos(rot2),sin(rot2))
    } else { 
      #make 3 branches with 120 degree angles between them and
      #rotate the entire branch set.
      rot <- runif(1,min=0,max=2*pi)
      tree$branch1 <- create_branch(tree,-sin(rot),cos(rot))
      tree$branch2 <- create_branch(tree,
                                    sqrt(3)/2*cos(rot)+sin(rot)/2,
                                    sqrt(3)/2*sin(rot)-cos(rot)/2)
      tree$branch3 <- create_branch(tree,
                                    -sqrt(3)/2*cos(rot)+sin(rot)/2,
                                    -sqrt(3)/2*sin(rot)-cos(rot)/2)
    }
  } else {
    tree$extension <- grow_tree(tree$extension)
    tree$branch1 <- grow_tree(tree$branch1)
    tree$branch2 <- grow_tree(tree$branch2)
    if (tree$depth>=0.5) 
      tree$branch3 <- grow_tree(tree$branch3)
  }
  return(tree)
}

create_ornaments <- function(tree) {
  if (is.null(tree)) return()
  po <- (1-tree$depth)^6
  ornament <- sample(c(T,F),size=1,prob=c(po,1-po))
  co <- sample(c("red","darkgoldenrod4"),size=1,
               prob=c(0.6,0.4))
  if (ornament)
    ornaments <<- rbind(ornaments,
                        data.frame(x=tree$x1,y=tree$y1,z=tree$z1,color=co))
  create_ornaments(tree$branch1)
  create_ornaments(tree$branch2)
  create_ornaments(tree$branch3)
  create_ornaments(tree$extension)
}

create_lights1 <- function(tree) {
  if (is.null(tree)) return()
  po <- (1-tree$depth)^4
  light <- sample(c(T,F),size=1,prob=c(po,1-po))
  if (light)
    lights1 <<- rbind(lights1,
                      data.frame(x=tree$x1,y=tree$y1,z=tree$z1))
  create_lights1(tree$branch1)
  create_lights1(tree$branch2)
  create_lights1(tree$branch3)
  create_lights1(tree$extension)
}

create_lights2 <- function(tree) {
  if (is.null(tree)) return()
  po <- (1-tree$depth)^4
  light <- sample(c(T,F),size=1,prob=c(po,1-po))
  if (light)
    lights2 <<- rbind(lights2,
                      data.frame(x=(tree$x1+tree$x0)/2,
                                 y=(tree$y1+tree$y0)/2,
                                 z=(tree$z1+tree$z0)/2))
  create_lights2(tree$branch1)
  create_lights2(tree$branch2)
  create_lights2(tree$branch3)
  create_lights2(tree$extension)
}

draw_ornaments <- function() {
  with(ornaments,
       {points3D(x=x,y=y,z=z,
                 pch=19,cex=1.2,col=as.character(color),
                 colkey=FALSE,add=TRUE)})
}

draw_lights1 <- function() {
  with(lights1,
       {points3D(x=x,y=y,z=z,pch="+",cex=0.8,col="white",
                 colkey=FALSE,add=TRUE)})
}

draw_lights2 <- function() {
  with(lights2,
       {points3D(x=x,y=y,z=z,pch="+",cex=0.8,col="yellow",
                 colkey=FALSE,add=TRUE)})
}

set.seed(20161224)
pine_tree <- trunk
for (i in 1:5) pine_tree <- grow_tree(pine_tree)
ornaments <- NULL;
lights1 <- NULL;
lights2 <- NULL;
create_ornaments(pine_tree)
create_lights1(pine_tree)
create_lights2(pine_tree)

png("tree%02d.png")
for (i in 0:35) {
  perspbox(x=c(-25,25),y=c(-25,25),z=c(0,40),bty="n",
           phi=8,theta=i*10,col="white",alpha=0)
  draw_tree(pine_tree)
  draw_ornaments()
  switch((i%%4)+1,{},
         {draw_lights1()},
         {draw_lights2()},
         {draw_lights1(); draw_lights2()})
}
graphics.off()

system("convert -delay 40 *.png christmasTree.gif")