# Merry Christmas 2016 (with R)

[This article was first published on

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

**U.N.A. Matemáticas El Tigre**, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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}\).

**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.

**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

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 | 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.

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

To

**leave a comment**for the author, please follow the link and comment on their blog:**U.N.A. Matemáticas El Tigre**.R-bloggers.com offers

**daily e-mail updates**about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.