Merry Christmas 2016 (with R)

I'd like to wish all my readers a Merry Christmas 2016- R style! Behold my 3d Christmas tree created using the plot3D R package:

While this might seem like yet another Christmas decoration done in R, it is unique in that the tree is rendered in 3d perspective. I myself wrote some code for a 2d Christmas tree a while ago (see this entry in Spanish), in turn based on the code in the Wiekvoet blog. This time around though, I took it a bit further- the challenge now being to create a 3d model of a Christmas tree that I could rotate and animate in perspective. Heck, if others ~~with enough leisure/boredom time to spare~~ were to join in and showcase their own Christmas R creations for this time of the year, we might be witnessing the start of a new Christmas tradition for all R geekdom across the planet.

Mathematical model of a Christmas tree

Long before the birth of Jesus, most pagan cultures in Europe were sun worshipers and in December they celebrated the rebirth of the sun and the rich harvests. This is the significance of the Christmas tree in these festivities of fire and light, abounding in evergreens and holly branches. In Spanish America, instead of Christmas trees, people would place Nativity scenes (Nacimientos) in their homes, but with globalization and the prevalence of the Anglo-Saxon culture, the Christmas tree eventually found its way into our homes. At any rate, I for one have always found the lights and the bright decorations of the Christmas tree captivating and reminiscent of happier moments when friends are remembered and grudges are forgotten. So in a sort of weird tribute to the spirit of these festivities, I decided I would make a mathematical model of the Christmas tree. How geeky is that?

Mathematically, a Christmas tree is nothing but a fractal. A fractal consisting of a stub or segment of a tree trunk with two or three branches that separate outward from the trunk, and an extension, or outward projection of the trunk in the same direction as the trunk. The branches can in turn be considered as trunk segments, each comprising more branches and an extension outward in the same direction. The extension or outward projection of a trunk is in turn another stub with more branching off points and so on. Consequently, the tree construction process is recursive. Being a descendant of LISP, base R supports the list as the data structure of choice for implementing this recursive definition of a tree.

There are, to be sure, some considerations to be made so that this recursively-built tree ends up looking like a pine tree and not a palm tree, a cherry blossom tree, a mango tree or a bonzai tree. For one thing, pine trees are characteristically conical. This means that the angle of separation of the branches from the trunk is initially larger but "closes up" as the pine tree extends out or up. To control the degree of "outwardness" or "upness" associated to a particular tree-trunk stub there will be a depth parameter. Initially, the pine tree starts as one little green trunk stub of depth 1. As this stub branches off and extends upward, each of the children stubs has a lower depth than the parent stub. The most outward branches and leafs of the tree have a depth of almost 0. This will be useful for the dressing the tree with lights and decorations, as these lights are usually found on the more superficial or outward parts of the tree.

Another aspect of our tree growth model worth mentioning is how the tree trunk becomes thicker as the tree grows out an up. This is done via a width parameter that will eventually determine how thick the lines representing the trunk and the branches should be drawn. This width parameter also determines when the tree stem changes color from dark green to brown. In each growth cycle, the growth algorithm increases the width of each stub by a factor of 1.4, it then goes through the branches and the extension looking for the ending stubs (those whose branches and extensions are NULL or nonexistent). It then adds the branches and the extension and backtracks down to find the remaining end stubs. The width of the new branches and extensions is initialized to 1, but as can be seen, the width of the older stubs has been increased more.

The Christmas tree data structure is a list whose components are:

The starting point of the trunk segment, given as a vector with coordinates in $\mathbb{R}^3$ as $(x_0,y_0,z_0)$.
The ending point of the trunk segment, as given by the vector with coordinates $(x_1,y_1,z_1)$. Together, the start and end points determine the direction vector of the tree stem as $\vec{u}=(x_1-x_0,y_1-y_0,z_1-z_0)$. This direction vector will be useful for creating the extension stub, since the extension stub grows in the same direction as the parent stub. It is also used for determining where along the stub the branches start off and in what direction those branches are created.
The width parameter lwd which is also the thickness with which the tree stem is drawn as a segment when plotted.
The depth, which indicates how outward a branch or tree stem is, as already explained above.
Slots for three branches and one extension (branch1, branch2, branch3 and extension), which are nothing but lists recursively defined like this one. When a branch or extension is created, these slots are initialized to NULL.

Creating the extension of a tree stem poses no major problem. One simply takes as starting point of the extension the ending point of the parent tree stem. The ending point of the extension is determined by adding to its starting point a multiple of the direction vector $\vec{u}=(x_1-x_0,y_1-y_0,z_1-z_0)$. This scalar multiple is such that the resulting length of the extension stem is slightly shorter than that of the parent stem. The most difficult issue in creating the Christmas tree lies with the creation of the branches. How can we create the branches so that they branch out from the stem in apparently uniform angles around the stem and not have all the branches branch out from the same side of the stem? How can we obtain the direction vectors for those branches so that they extend outwards from the stem? A little trigonometry and vector geometry is of use to us here.

If the $\vec{u}$ vector is the direction of the tree stem being considered, then what we need is a vector perpendicular to $\vec{u}$ - let's call it $\vec{v}$ - so that we can obtain the direction vector of the branch $\vec{b}$ by adding $\vec{u}+\vec{v}$ and then multiplying $\vec{b}$ by a scalar to set its length to the desired magnitude (see Fig. 1). However, considering we're dealing with the $\mathbb{R}^3$ metric space, there are infinitely many such vectors. In fact, there is an entire two-dimensional space - a plane $\mathcal{V}$ comprised of vectors that are all perpendicular to $\vec{u}$ (see Fig. 2). So what we need to obtain $\vec{v}$ is a orthonormal base of two vectors $\{\vec{v_1},\vec{v_2}\}$ that will permit us to express $\vec{v}$ as a linear combination of the two base vectors of $\mathcal{V}$.

Geometrical relationship between a tree stem and one of its branches

Fig. 1 Geometrical relationship between the main tree stem ($\vec{u}$) and one of its branches ($\vec{b}$). $\theta$ is the angle at which $\vec{b}$ branches outward from the main stem.

This is easy enough. Two vectors are said to be orthogonal if their dot-product is zero. So if we have a vector, say $(x,y,z)$ all we need to do is pick any non-zero coordinate and switch it with one of the remaining coordinates, changing one of the signs to negative and setting the remaining third coordinate to zero. So for example $(-z,0,x)$, $(0,z,-y)$, and $(-y,x,0)$ are all perpendicular to $(x,y,z)$, provided they are all non-null vectors. In short, to find an orthonormal base for the plane perpendicular to a vector $\vec{u}$, pick any non-zero component of this vector, pick a second component to interchange it with while changing the sign of any of the two and set the third component to zero. This will give you a perpendicular vector to $\vec{u}$. To find the other perpendicular vector, just repeat the above process, but interchanging the chosen non-zero component of $\vec{u}$ with the component you had set to zero for the first vector. Finally, once you have your two vectors $\vec{v_1}$ and $\vec{v_2}$, you must multiply them by the inverse of their norms to set their norms to one. The procedure is better illustrated in the R code on this post.

Vector v as a linear combination of v1 and v2

Fig. 2 Finding vector $\vec{v}$ as a linear combination of an orthonormal base of perpendicular plane $\mathcal{V}$ and then using it to find the direction vector $\vec{b}$ of the branch.

So once you have the base of the perpendicular plane to a given tree stem $\vec{u}$, all you need is two scalar components $c_1$ and $c_2$ to define the perpendicular vector $v$ used for constructing the branch, as a linear combination of the base. How are we to go about this? Well, let's think about what we really need. For some parts of the tree - the deeper parts - we need three branches to a stem. For the more outward parts, two branches to a stem will suffice. We want these branches to be evenly distributed about the stem: in the case of 3 branches, we would like the angle of separation between the branches to be $120^\circ$ (like the Mercedes-Benz star symbol), in the case of 2 branches, an angle of separation between them of about $180^\circ\pm 20^\circ$ should be adequate.

What this means is that we want to choose the $\vec{v}$ vectors in $\mathcal{V}$ to have the sort of layout shown in Fig. 3. This is to ensure the pine tree will have a nice conical shape and not have its branches all projecting to the same side.


Fig. 3a 3 branch defining vectors	Fig. 3b 2 branch defining vectors
Fig. 3 Distribution of two and three branch generating vectors on $\mathcal{V}$

Our problem now is to generate 2 or 3 vectors on a plane with such angles of separations. The good news is that we can define the vectors on the regular 2d plane with the canonical base $\{\hat{i},\hat{j}\}$, and, because linear algebra is so cool, we can use the components of those vectors as $c_1$ and $c_2$ - the scalar coeficients used to define the $\vec{v}$ vectors as linear combinations of $\vec{v_1}$ and $\vec{v_2}$. This is so because $\{\vec{v_1},\vec{v_2}\}$ is an orthonormal base of a plane in the 3d space. So for example if the angle between $(0,1)$ and $(-\tfrac{1}{2},\tfrac{\sqrt{3}}{2})$ is $120^\circ$, then the angle between $\vec{v_2}$ and $-\tfrac{1}{2}\vec{v_1} + \tfrac{\sqrt{3}}{2}\vec{v_2}$ will also be $120^\circ$. Not only that, but the two vectors will lie on the $\mathcal{V}$ plane in three dimensional space. Thus, we can translate the 2 dimensional vectors into vectors located on the $\mathcal{V}$ plane perpendicular to the tree stem $\vec{u}$.

Choosing the 2 or 3 vectors on the regular 2d plane later translating them as vectors on the $\mathcal{V}$ plane is a question of doing the following:

2 branches: We first consider $(1,0)$ as the first vector and we determine the second unit-length vector by choosing a random angle between $160^\circ$ and $200^\circ$. We choose another random angle between $0^\circ$ and $360^\circ$ to rotate the entire set of two vectors.
3 branches: Our three vectors will initially be $(1,0)$, $(\tfrac{\sqrt{3}}{2},-\tfrac{1}{2})$, and $(-\tfrac{\sqrt{3}}{2},-\tfrac{1}{2})$. We then choose a random angle between $0^\circ$ and $360^\circ$ to rotate the entire set of three vectors.

The reason we subsequently rotate the two or three vector set is to randomize the orientation of the branches so that our trees branches will project in all sorts of different directions. In either case, to rotate a 2d vector by an angle of $\theta$ about the origin, we use the following linear transformation:

\[ (x,y)\quad\mapsto (x\,cos(\theta)-y\,sin(\theta),x\,sin(\theta)+y\,cos(\theta)) \]

There are a few remaining details to explain before we go into the R code. As we have already mentioned, the more superficial stems of the three have two branches whereas the deeper parts of the tree (those closest to the original thick stub from which the entire Christmas tree branches out) have three branches. This is controlled by the depth parameter. Whenever $depth>0.8$, the tree stem will branch out into three branches, otherwise it will have only two. The depth parameter also controls where along the stem the branches are chosen. Stems of greater depth tend to branch off closer to its ending point $(x_1,y_1,z_1)$ whereas the more superficial stems branch off at points closer to its starting point $(x_0,y_0,z_0)$. Needless to say, the branching off points are chosen randomly. Finally, the placement of the lights and the tree decorations is also controlled by the depth parameter. The more superficial parts of the tree have a higher likelihood of containing lights and decorations. Again, the placement of these decorations is done stochastically.

R script for the 3d Christmas tree

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()

Bibliography

Apostol, T. (1985). Calculus, Vol. II (2nd ed). (F. Veléz trans.). Caracas: Editorial Reverté.
carlop. (April 26, 2012). Answer to "How to find vector perpendicular to another vector?" on math.stackexchange. [Retrieved 12/20/2016 from http://math.stackexchange.com/questions/137362/how-to-find-perpendicular-vector-to-another-vector].
Cooney, B. (1967). Christmas. New York: Thomas Y. Crowell Company.
Duineveld, K. (2013, December). "Merry Christmas". Blog post. [Retrieved 12/20/2016 from http://wiekvoet.blogspot.com/2013/12/merry-christmas.html].
R DEVELOPMENT CORE TEAM (2009). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0. http://www.R-project.org.
Soetaert, K. (2016). plot3D: Plotting multi-dimensional data. R package version 1.1. https://cran.r-project.org/web/packages/plot3D/index.html.

Dear Reader:

If you found this post interesting or useful, please share it on Google+, Facebook or Twitter so that others may find it too.

Below you will find a link to download a printable pdf version of this article for off-line reading and the R source code featured in this post, but essentially it's the same information I'm sharing with you in this blog entry. It will take you to a Bitcoin payment gateway which is very secure, anonymous and hassle free (I will not collect your sensitive personal information). The fee is small - 0.0005 BTC - or roughly $0.45 at today's rate. Think of this as a donation, not a business transaction that will help me continue writing content and survive in this hell-hole that Venezuela has become. Thanks in advance!