Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
1517 lines (1269 sloc) 54.1 KB
title: A New Stereoscopic Mount Spectacular
author: ~
date: '2018-12-12'
slug: three-d-pipeline
categories: [r]
tags: [optimization,visualization]
```{r echo=FALSE}
comment = "", fig.align='center'
# Cool, but Not Eye Popping
style='float: left; margin: 5px 15px 5px 0;' width='250'/>
As I finished the [ray shading post](/2018/10/23/do-not-shade-r/), it occurred
to me that it would be pretty cool if we could make these awesome 3D-looking
reliefs actually 3D, or at least stereoscopic. After all, how hard could it
be? We have the elevation data, so some simple trigonometry should allow us to
tilt things left, tilt things right, and voilà. So I did that and it looked
like garbage.
That would have been the right time to walk away and move on to more productive
things. Instead, I started messing around with [`rgl`][12] (fantastic package),
but I struggled to get it to do *exactly* what I wanted programmatically. While
trying to understand how `rgl` reacts to inputs I read up on 3D projections,
and realized that maybe this could all be done outside of `rgl`.
And so began this quixotic challenge[^quixotic] to implement an
"efficient-enough"[^how-fast-is-useful] 3D rendering pipeline using base R only.
I took notes along the way to share with those who like me are interested, but
unlike me have better time management skills. So, come along for a ride through
3D projections, meshes, perspective adjustments, rasterization, and image
processing, all in base R[^base-r-only].
**DISCLAIMER**: I knew nothing about 3D rendering before starting this post, and
now I know just enough to be dangerous. Do not look for 3D rendering best
practices here. My criteria for using things in this post is: "does it _appear_ to work".
# Volcano!
Our subject today will be the much-loved `volcano`, a.k.a. Maunga Whau. We are
going to pretend that the Z-values of the `volcano` matrix are in the same units
as the grid size (i.e. 10m per unit, as per `?volcano`). This is obviously not
the case, but the documentation is just ambiguous enough that I feel comfortable
indulging in this minor fraud for dramatic effect. Here is a height-map of
```{r helper-funs, echo=FALSE}
## Rescale data to a range from 0 to `range` where `range` in (0,1]
rescale <- function(x, range=1, center=0.5)
((x - min(x, na.rm=TRUE)) / diff(range(x, na.rm=TRUE))) * range +
(1 - range) * center
## Prepare a plot with a particular aspect ratio
plot_new <- function(
x, y, xlim=c(0,1), ylim=c(0,1),
par.args=list(mai=numeric(4L), xaxt='n', yaxt='n', xaxs='i', yaxs='i')
) {
if(length(par.args)), par.args)
xlim, ylim, asp=diff(range(y, na.rm=TRUE))/diff(range(x, na.rm=TRUE))
) }
```{r, fig.width=6, fig.height=4}
xv <- row(volcano)
yv <- col(volcano)
plot_new(xv, yv)
points(x=rescale(xv), y=rescale(yv), col=gray(rescale(volcano)), pch=15)
`plot_new` is a thin wrapper around `` that also sets some `par`
values. `rescale` rescales values to between zero and one. Both of these are
[helper functions defined in the appendix](#helper-functions).
As we saw from the [ray shading post](/2018/10/23/do-not-shade-r/), it is
possible to ray shade elevation maps using base R only. For the sake of
exposition I put the code in the [shadow demo package][15], but as far as I'm
concerned I'm still abiding by the base-R-only pledge.
```{r, fig.width=6, fig.height=4}
shadow <- ray_shade2(volcano, seq(-90, 90, length=25), sunangle=180)
plot_new(xv, yv)
points(x=rescale(xv), y=rescale(yv), col=gray(shadow), pch=15)
# 3D Projection
While volcano looks pretty cool from above, we want to provide a more
interesting view. To do this we can use [3D rotation matrices][1]. `shadow`
implements the `rot_*` functions, which are thin wrappers that generate the
rotation matrices described in the [wikipedia article][1]:
round(rot_z(90), 3)
By our convention the Z-axis is always pointed directly at the user, even after
rotations are applied. In other words, it is the model that rotates, not the
viewer[^1]. This also means that the X and Y axes remain in their traditional
In order to use the rotation matrix we need our data in long format. This can
be done easily by taking advantage of the `row` and `col` functions:
mx.example <- matrix((0:3)*10, nrow=2, ncol=2)
rbind(x=c(row(mx.example)), y=c(col(mx.example)), z=c(mx.example))
While this looks rather wide for long format data, I assure you it is long in
spirit. We need it in 3 x n rather n x 3 format for the upcoming matrix
multiplication. As applied to `volcano`:
volc.l <- rbind(x=c(row(volcano)), y=c(col(volcano)), z=c(volcano))
ncol(volc.l) # long data has same number of colums ...
length(volcano) # ... as volcano has elements
The 3D projection is then just a matrix multiply:
rot <- rot_x(-20) %*% rot_z(65) <- rot %*% volc.l
We need to re-order by Z values to ensure that the nearest points are plotted
last so they are not occluded by further points. Recall the Z axis is pointed
directly at you, so points with higher Z values are closer.
```{r, fig.width=6}
zord <- order([3,])
vord <-[, zord] # reorder by Z-values
x <- vord[1,]
y <- vord[2,]
plot_new(x, y)
points(rescale(x), rescale(y), col=gray(shadow[zord]), pch=15, cex=0.5)
The good news is now we have an awesome looking 3D rendition of our volcano...
in the moments after a direct hit from a disruptor cannon.
# Interlude I
As you might have guessed one of the reasons why people don't usually implement
3D rendering pipelines in base R is performance. If we are to have any hope of
creating a pipeline that renders things in _useful_[^how-fast-is-useful] time we
need to be careful about how we structure our data and calculations.
In an earlier blog post we found that [list and list-matrix data structures are
often faster to operate on][2] than the pure-vector equivalent matrices and
arrays. Here we turn our long volcano matrix into list format:
<a name='volcano-list-long'></a>
vl <- lapply(seq_len(nrow(, function(x)[x,])
vl <- c(vl, list(shadow)) # add texture info
names(vl) <- c('x', 'y', 'z', 't') # 't' for texture
Throughout the rest of this blog post we will use techniques explored at length
in the aforementioned post. If you are not familiar with `lapply`, ``,
`Map`, `Reduce`, and similar functional list manipulation facilities, you should
consider reading [that post][2]. We'll wait for you right here =).
# Meshes
## Tiles
As our disintegrating volcano demonstrates, we can't just use elevation data
as point coordinates. One possibility is to turn our point "cloud" into a mesh.
It turns out that making a mesh out of an arbitrary point cloud is one of those
"seems easy but isn't" problems. Worse, the typical implementations
involve some kind of sequential process that is repeated many times until the
mesh is complete. These types of algorithms cannot be implemented efficiently
in R, though there are R packages that implement them in compiled code.
Fortunately for us and our base-R-only challenge, our point cloud is not
arbitrary. The source data contains the information required to build the mesh:
the original x-y grid. We can generate a tile mesh by subsetting a contiguous
portion of our elevation matrix four times, once for each of the four vertices
of each tile. This is what it looks like with an illustrative 4 x 4 grid:
```{r echo=FALSE, fig.width=7, fig.height=1.75}
p.seq <- seq(.1, .9, length.out=4)
points <- expand.grid(x=p.seq, y=p.seq)
len <- length(p.seq)
pad <- .05 <- data.frame(
x=c(pad, p.seq[3] + pad, p.seq[3] + pad, pad),
y=c(-pad + p.seq[2], -pad + p.seq[2], pad + p.seq[4], pad + p.seq[4])
seq.diff <- diff(p.seq[1:2])
#, height=1.5)
old.par <- par(mfrow=c(1,5), mar=c(.5, .25, 2, .25)) <- c(0, 1, 1, 0) <- c(1, 1, 0, 0)
cols <- c(
"#440154FF", "#472D7BFF", "#3B528BFF", "#2C728EFF", "#21908CFF",
"#27AD81FF", "#5DC863FF", "#AADC32FF", "#FDE725FF")
col.basic <- matrix('black', 4, 4)
bg.basic <- matrix(NA_character_, 4, 4)
plot_new(0:1,0:1, par.args=NULL)
tiles <- shadow:::mesh_tile(c(as.list(points), list(1:16)), c(4,4))
x4 <-, c(tiles[,'x'], list(NA)))
y4 <-, c(tiles[,'y'], list(NA)))
polygon(c(0, 1, 1, 0), c(0, 0, 1, 1), col='#EEEEEE', border=NA)
polygon(x4, y4, col=cols, border=NA)
points(points, pch=22, col='black', bg=NA)
for(i in seq_along(p.seq)) {
plot_new(0:1,0:1, par.args=NULL)
title(sprintf('Vertex %d', i))
polygon(c(0, 1, 1, 0), c(0, 0, 1, 1), col='#EEEEEE', border=NA) <- col.basic <- bg.basic[if([i]) -1 else -4, if([i]) -4 else -1] <- cols[if([i]) -1 else -4, if([i]) -4 else -1] <- cols
points(points, pch=22,,
polygon([['x']] + seq.diff *[i],[['y']] - seq.diff *[i]
To form each tile as shown in the first panel we take one vertex from each
of the four vertex panels. Each panel has the 9 of 16 grid points available
for use as vertices enclosed in a black bounding box. For the top left
green tile, we take the top left point from each of the bounding boxes. The
same process is repeated for each tile by picking the corresponding vertex from
each panel. The vertices are colored by the tile they belong to.
The operation we just described is equivalent to alternating dropping the
first/last row/column. Because our projected volcano data is already in long
format, we need to compute which elements in the long data correspond to
the first and last rows. This is actually pretty straightforward if we recall
that our matrix in long format is basically the component vectors "unspooled".
We can "re-spool" an index vector the length of our data, `seq_along(volcano)`,
by making a matrix out of it with the correct dimensions:
nr <- dim(volcano)[1]
nc <- dim(volcano)[2]
idx.raw <- matrix(seq_along(volcano), nr, nc)
Exploiting the underlying vector nature of matrices, and vice versa, is a common
technique when trying to write fast R code.
We can now drop rows and columns in sequence from that index matrix to get the
long format indices for each vertex for each tile:
idx.tile <- list(
v1=c(idx.raw[-nr, -nc]), v2=c(idx.raw[-nr, -1]),
v3=c(idx.raw[ -1, -1]), v4=c(idx.raw[ -1, -nc])
For example, to get the location of vertex data for tile forty two:
sapply(idx.tile, '[', 42)
This means that data for vertex 1 for tile 42 is at position 42 in the long data,
for vertex 2 it is at position 129, and so on.
Recall that the volcano data is in `vl` in [long list
So to get the `x`, `y`, `z` coordinates, and "texture" values for the vertices
of tile forty two we use:
sapply(vl, '[', sapply(idx.tile, '[', 42))
In order to extend this to the full tile set we are going to use some `for`
loops. After all `sapply`, `lapply`, and friends are just for
loops themselves. Additionally in this case we want a specific result structure
that is just easier to fill with for loops. Finally, there will only be sixteen
iterations of the loop, so the overhead from R-level function calls is small.
## Initialize result structure
mesh.tile <- matrix(
list(), nrow=length(idx.tile), ncol=length(vl),
dimnames=list(names(idx.tile), names(vl))
## Fill it with the correctly subset volcano data
for(i in names(idx.tile))
for(j in names(vl))
mesh.tile[[i,j]] <- vl[[j]][idx.tile[[i]]]
This is an instance of the fabled list-matrix. It looks a little weird, but it
is equivalent to a 3D array and is just a list of equal length vectors with a
`dim` attribute:
One problem with this data structure is that it is duplicative. Most vertices
are shared by four tiles, yet we make a copy of each vertex for each tile. This
is one of the trade-offs we accept for improved overall
To extract tile forty two's data from this structure we can use:
apply(mesh.tile, 1:2, function(x) unlist(x)[42])
The required `unlist` makes this a little awkward, but we do not intend to
operate on a single tile at a time. For most of our manipulations, the
list-matrix is convenient. Look how easy it is to display the mesh from raw
<a name='tile-mesh-plot'></a>
zord <- order(Reduce('+', mesh.tile[,'z']))
x <-, c(mesh.tile[,'x'], list(NA)))[,zord]
y <-, c(mesh.tile[,'y'], list(NA)))[,zord]
plot_new(x, y)
polygon(rescale(x), rescale(y), col='#DDDDDD', border='#444444')
```{r echo=FALSE}
x4 <- x # for later re-use
y4 <- y
It may help to look at the data structure of `x` above to understand what's
going on:
`y` shares the same structure. In order to have `polygon` render multiple
polygons, we need to NA-delimit our input vectors. Rather than painstakingly
creating those vectors from the undelimited ones, we `rbind` the four vertex
vectors with one NA vector. That this produces a matrix is besides the point;
the idea is that the vector wrapping nature of the matrix allows us to sandwich
NA values between each of our tiles' data:
Great, but did we _really_ need to subject ourselves to this rigmarole?
Wouldn't `persp`, which would still abide by our base-R-only pledge, do this
just as well and without ridiculous list-matrices? Well, [kind
of](#meshes-with-persp), but more importantly, "this" is not what we're after.
We have far grander ambitions, and "this" is basically about the limit of what
we can do with `persp`.
## Triangles
A tile mesh is not ideal. Among other things, there is no guarantee the
tiles are flat. Triangles on the other hand are guaranteed to be flat, and
additionally benefit from simple shading algorithms. So we will turn our tiles
into triangles:
```{r echo=FALSE}
rescale(c(0, 1, 1, 0), .9), rescale(c(0, 0, 1, 1), .9),
col='#DDDDDD', border='#444444'
rescale(c(0, 1, 1, 0), .95), rescale(c(0, 0, 1, 1), .95),
labels=sprintf('v%d', 1:4)
rescale(c(0, 1, 1), .8) + .0125, rescale(c(0, 0, 1), .8) - .0125,
border='#0000FF', col='#CCCCFF'
rescale(c(0, 1, 1), .7) + .03, rescale(c(0, 0, 1), .7) - .03,
labels=sprintf('v%d', 1:3), col='#0000AA'
rescale(c(1, 0, 0), .8) - .0125, rescale(c(1, 1, 0), .8) + .0125,
border='#00FF00', col='#CCFFCC'
rescale(c(1, 0, 0), .7) - .03, rescale(c(1, 1, 0), .7) + .03,
labels=sprintf('v%d', c(3,4,1)), col='#00AA00'
The list-matrix structure of our data makes the translation from tile to
triangle mesh trivial:
```{r} <- 1:3 <- c(3,4,1)
mesh.tri <- Map('c', mesh.tile[,], mesh.tile[,])
dim(mesh.tri) <- c(3, 4)
dimnames(mesh.tri) <- list(head(rownames(mesh.tile), -1), colnames(mesh.tile))
zord3 <- order(Reduce('+', mesh.tri[,'z']))
x <-, c(mesh.tri[,'x'], list(NA)))[,zord3]
y <-, c(mesh.tri[,'y'], list(NA)))[,zord3]
col <- rep(c('#CCCCFF','#CCFFCC'), each=ncol(x)/2)[zord3]
plot_new(x, y)
polygon(rescale(x), rescale(y), col=col)
And with more realistic colors by shading each facet with the mean of the vertex
texture <- gray(Reduce('+', mesh.tri[,'t']) / 3)[zord3]
plot_new(x, y)
polygon(rescale(x), rescale(y), col=texture, border=texture)
```{r echo=FALSE}
texture3 <- texture # for later use
Now we're cooking!
# Perspective
If I had any perspective I wouldn't be writing this blog post. But that can't
be helped. Our dear friend `volcano` is a different story.
Everything we've done to this point has assumed no perspective effect; objects
of the same size appear the same size irrespective of how far from the observer
they are. This observer, you, has been looking at the object
along a line parallel to the Z axis and passing through the midpoint of the X
and Y extents of the projected object. For expository purposes we are about to,
for the first and only time, rotate Observer 0 (you) 90 degrees about the Y
In the following diagram you are looking along a line parallel to the X axis.
The Z axis now runs from left to right. Observers 1, 2, and 3 are all still
looking at `volcano` along the Z axis:
```{r echo=FALSE, fig.width=7, fig.height=2.33}
rot2 <- rot_x(-20) %*% rot_z(65) %*% rot_x(-90)
volc.l2 <- rbind(x=c(row(volcano)), y=c(col(volcano)), z=c(volcano))
volc.lr2 <- rot2 %*% volc.l
vl2 <- lapply(seq_len(nrow(volc.lr2)), function(x) volc.lr2[x,])
vl2 <- c(vl2, list(shadow)) # add texture info
names(vl2) <- c('x', 'y', 'z', 't') # 't' for texture
mesh.tri2 <- shadow::mesh_tri(vl2, dim(volcano))
zord32 <- order(Reduce('+', mesh.tri2[,'z']))
x32 <-, c(mesh.tri2[,'x'], list(NA)))[,zord32]
y32 <-, c(mesh.tri2[,'y'], list(NA)))[,zord32]
texture32 <- gray((Reduce('+', mesh.tri2[,'t'])/4))[zord32]
yscale <- .8
ylow <- (1 - yscale) * .9
ymid <- .25 + yscale/2
yarrow <- rescale(c(0, .5, 1), yscale * .5, 0) + ylow
yarmax <- max(yarrow)
polygon(rescale(x32, 1/3, 1), rescale(y32, yscale, 1),
col=texture32, border=texture32
x0=c(0, 1/3, 1/2, 2/3),
x1=c(2/3,2/3, 2/3, 1),
y0=c(yarrow, ylow),
y1=c(yarrow, ylow),
col=c(rep('blue', 3), '#00DD00'),
c(0, 0, NA, 1/3, 1/3, NA, 1/2, 1/2, NA, 2/3, 2/3, NA, 1, 1, NA),
c(ylow, ymid, NA,
ylow + (yarmax - ylow) * .5, ymid, NA,
yarmax, ymid, NA,
ylow, .75, NA,
ylow, .83, NA),
points(c(0, 1/3, 1/2), y=rep(ymid, 3), col='blue', pch=16)
x=c(2/3 * .5, 1/3 + 1/6, 1/2 + 1/12, 2/3 + 1/6),
y=c(yarrow - ylow/2, ylow/2),
col=c(rep('blue', 3), '#00DD00'),
label=c('D = ZD * 2', 'D = ZD', 'D = ZD / 2', 'ZD = Z Depth'),
x=c(0, 1/3, 1/2), y=rep(ymid, 3), col='blue', pos=3,
label=sprintf('Obs %d', 1:3)
For convenience we define the quantity `ZD` as the total depth of our object
along the Z axis. Then, we can define the distance of our observers from the
nearest point along the Z axis with the measure `D` as multiples of `ZD`.
Here is our volcano re-rendered from the perspective of each of observers 1-3
pictured above:
```{r echo=FALSE, fig.width=7, fig.height=3}
par(mfrow=c(1,3), mar=c(.5, .25, 2, .25))
D.vals <- c(2,1,.5)
for(i in seq_along(D.vals)) {
vlp <- shadow::persp_rel(vl, D.vals[[i]])
vlpm <- shadow::mesh_tri(vlp, dim=dim(volcano))
xi <-, c(vlpm[,'x'], list(NA)))[,zord3]
yi <-, c(vlpm[,'y'], list(NA)))[,zord3]
# zord <- order(Reduce('+', vlpm[,'z']))
# texture <- Reduce('+', vlpm[,'t'])/3
plot.window(c(0,1), c(0,1), asp=diff(range(vlp['x']))/diff(range(vlp['y'])))
title(sprintf('D = %.01f', D.vals[i]))
polygon(rescale(xi), rescale(yi), col=texture3, border=texture3)
The actual perspective computation is best applied to the original point cloud,
so we go back to the `vl` variable, which is `volcano` in
[long list format](#volcano-list-long), and apply the perspective transformation
for observer 3:
str(vl) # `volcano` in long list format
```{r eval=FALSE}
D <- 0.5 # Observer 3 is half Z-Depth from volcano
z.rng <- range(vl[['z']])
ZD <- diff(z.rng)
## Center observer on x-y extent and at Z=0
vlp <- vl
vlp[c('x','y')] <- lapply(vl[c('x','y')], function(x) x - sum(range(x)) / 2)
vlp[['z']] <- vlp[['z']] - (z.rng[2] + D * ZD)
## Apply perspective transformation
z.factor <- -1 / vlp[['z']]
vlp[c('x','y')] <- lapply(vlp[c('x','y')], '*', z.factor)
```{r echo=FALSE, eval=FALSE}
## just to confirm above works
vlp2 <- sapply(vlp, '[', order(vlp[['z']]), simplify=FALSE)
x=rescale(vlp2[['x']]), y=rescale(vlp2[['y']]), col=gray(vlp2[['t']]),
pch=15, cex=0.5
Now that we have the perspective-transformed point cloud, we can turn it back
into a triangle mesh. This time we use `shadow::mesh_tri`, which is a wrapper
around the calculations described in the [Mesh section](#meshes):
<a name='mesh-tri-persp'></a>
mesh.tri <- mesh_tri(vlp, dim(volcano), order=TRUE)
x <-, c(mesh.tri[,'x'], list(NA)))
y <-, c(mesh.tri[,'y'], list(NA)))
texture <- gray((Reduce('+', mesh.tri[,'t'])/nrow(mesh.tri)))
plot_new(x, y)
polygon(rescale(x), rescale(y), col=texture, border=texture)
`mesh_tri` has the additional benefit that it orders the tiles by Z value so we
do not need to worry about re-ordering before plotting anymore.
Are we there yet? Nope, we're going to drag this baby kicking and screaming
all the way to its painful, bloody conclusion.
# Rasterization
## Overview
So far we have relied on base R plotting facilities, `polygon` mostly, to
actually transform our 3D object representation into an image. There are two
problems with this:
1. `volcano` does not look very natural with all the triangular faces.
2. Speed is a problem with higher polygon counts.
We will resolve these issues by resorting to rasterization, whereby we compute
what each pixel of our final image should look like. This requires determining
which points of our final output pixel grid belong to which triangles. While
this is the type of thing a child could do easily, it is surprisingly difficult
to instruct a computer to do the same, particularly when we need to account for
the need to use internally vectorized[^internally-vectorized] base R code.
Conceptually, the process takes three steps:
```{r raster-overview, echo=FALSE, fig.width=7, fig.height=2.5}
par(mfrow=c(1, 4), mar=c(.5, .25, 2, .25))
xr <- c(0, .8, .2, NA, .8, 1, .2, NA)
yr <- c(0, .2, .8, NA, .2, 1, .8, NA)
triangles <- function() {
c(0, .8, .8, 0, NA, .2, 1, 1, .2, NA),
c(0, 0, .8, .8, NA, .2, .2, 1, 1, NA),
border=c('#0000FF', '#00DD00'), col=NA
plot_new(xr, yr, par.args=list(mar=c(.5, .25, 2, .25)))
title('Mesh Polygons\nto Rasterize')
polygon(xr, yr, col=c('#CCCCFF','#CCFFCC'), border=NA)
plot_new(xr, yr, par.args=list(mar=c(.5, .25, 2, .25)))
title('Step 1:\nBounding Box')
rast_1 <- function(x, y) {
polygon(x, y, col=c('#CCCCFF','#CCFFCC'), border=NA)
rast_1(xr, yr)
plot_new(xr, yr, par.args=list(mar=c(.5, .25, 2, .25)))
title('Step 2:\nCandidate Pixels')
len <- 13
pdat <- as.matrix(
expand.grid(x=seq(0, .8, length.out=len), y=seq(0, .8, length.out=len))
rast_2 <- function(x, y, pdat) {
polygon(x, y, col=c('#CCCCFF','#CCFFCC'), border=NA)
points(pdat, pch=24, col='#0000FF')
points(pdat + .2, pch=25, col='#00DD00')
rast_2(xr, yr, pdat)
plot_new(xr, yr, par.args=list(mar=c(.5, .25, 2, .25)))
title('Step 3:\nActual Pixels')
nr <- nrow(pdat)
mx1 <- matrix(c(lapply(xr[1:3], rep, nr), lapply(yr[1:3], rep, nr)), 3)
colnames(mx1) <- c('x', 'y')
mx1.oob <- Reduce('|', lapply(barycentric(, mx1), '<', 0))
mx2 <- matrix(c(lapply(xr[5:7], rep, nr), lapply(yr[5:7], rep, nr)), 3)
colnames(mx2) <- c('x', 'y')
mx2.oob <- Reduce('|', lapply(barycentric(, mx2), '<', 0))
rast_3 <- function(x, y, pdat, oob1, oob2) {
polygon(x, y, col=c('#CCCCFF','#CCFFCC'), border=NA)
points(pdat[!oob1,], pch=24, bg='#0000FF', col=NA)
points(pdat[!oob2,] + .2, pch=25, bg='#00DD00', col=NA)
rast_3(xr, yr, pdat, mx1.oob, mx2.oob)
<div style='display: none;'>
We use these for the little inset thumbnails.
```{r rast-one, fig.width=3, fig.height=3}
plot_new(xr, yr, par.args=list(mai=c(.1,.1,.1,.1)))
rast_1(xr, yr)
```{r rast-two, fig.width=3, fig.height=3}
plot_new(xr, yr, par.args=list(mai=c(.1,.1,.1,.1)))
rast_2(xr, yr, pdat)
```{r rast-three, fig.width=3, fig.height=3}
plot_new(xr, yr, par.args=list(mai=c(.1,.1,.1,.1)))
rast_3(xr, yr, pdat, mx1.oob, mx2.oob)
## Step 1: Bounding Boxes
style='float: left; margin: 5px 15px 5px 0;' width='250'/>
One thing that is easy to do is to compute the rectangular bounding box for each
triangle. All we need to do is compute the x and y extents of the vertices, and
use those to generate the corners of our bounding box.
Before we do this we need to settle on the resolution of our output grid as
everything is going to be measured in relation to that. This is easiest to do
with the perspective transformed point cloud data rather than the mesh version.
For exposition purposes we keep jumping back and forth between these formats,
although in the final rendering pipeline we will only carry out the mesh
transformation once we no longer have a need for the point data.
<a name='bounding-box-scale'></a>
x.width <- 600L
x.rng <- range(vlp[['x']])
y.rng <- range(vlp[['y']])
asp <- diff(y.rng) / diff(x.rng)
y.width <- as.integer(x.width * asp)
vlpt <- vlp
vlpt[['x']] <- rescale(vlpt[['x']], center = 0) * (x.width - 1) + 1
vlpt[['y']] <- rescale(vlpt[['y']], center = 0) * (y.width - 1) + 1
We can regenerate the mesh, and for each mesh polygon compute the maximum x and y
integer pixel values they will contain:
mesh <- mesh_tri(vlpt, dim(volcano), order=TRUE)
ceil_min <- function(x) list(as.integer(ceiling(, x))))
floor_max <- function(x) list(as.integer(, x)))
mins <- unlist(apply(mesh[,c('x','y')], 2, ceil_min), recursive=FALSE)
maxs <- unlist(apply(mesh[,c('x','y')], 2, floor_max), recursive=FALSE)
str(mins) # `maxs` will have same structure
We use `list(...)` in the ceil/floor functions to prevent `apply` from
simplifying to a numeric matrix, which in turn requires the unlist to get the
desired structure.
This is what bounding boxes look like:
minmax <- unname(c(mins, maxs))
xs <- range(unlist(minmax[c(1,3)]))
ys <- range(unlist(minmax[c(2,4)]))
xss <- lapply(minmax[c(1,3)], function(x) (x - xs[1]) / diff(xs))
yss <- lapply(minmax[c(2,4)], function(y) (y - ys[1]) / diff(ys))
plot_new(xs, ys)
rect(xss[[1]], yss[[1]], xss[[2]], yss[[2]], border='#00000033', asp=1)
## Step 2: Candidate Pixels
style='float: left; margin: 5px 15px 5px 0;' width='250'/>
We will brute force generate pixels to fill each of our bounding boxes. Oddly
this is a little tricky, mostly because we wish to do this fully in internally
vectorized code. And it is particularly important that we do so because we are
about to dramatically increase the amount of data we are dealing with. The plan
of attack is to compute all possible pixel "offset" permutations within each
bounding box, and then add back the min x and y coordinates to get the actual
pixel coordinates.
First we compute the dimensions of the bounding x/y min/max boxes in pixels:
dims.x <- pmax(maxs[['x']] - mins[['x']] + 1, 0)
dims.y <- pmax(maxs[['y']] - mins[['y']] + 1, 0)
dims.xy <- dims.x * dims.y
Then we generate all the offset permutations. This is the type of thing that
one might normally do with `expand.grid`:
```{r eval=FALSE}
Map(function(x, y) expand.grid(x=seq_len(x), y=seq_len(y)), dims.x, dims.y)
Unfortunately, with 10,320 polygons[^2] it becomes computationally burdensome.
We can do a little better with lower overhead variations[^3], but even then we
are left with lists with thousands of elements that are awkward to manage.
Instead we want to generate vectors that contain all the permutations of x and y
offsets we need. We do so here with a toy example for three bounding boxes with
x and y dimensions defined by `dim.x0` and `dim.y0`:
dims.x0 <- c(2, 2, 3)
dims.y0 <- c(2, 3, 2)
dims.xy0 <- dims.x0 * dims.y0 <- rep(seq_along(dims.xy0), dims.xy0)
```{r echo=FALSE}
col <- c('#0000FF', '#00DD00', '#CCCC00')
bb.col <- col[]
bb.html <- sprintf(
'<span style="background-color:%s; color:#EEEEEE;">%d</span>', bb.col,
Here we took advantage of a less known feature of `rep` whereby if the second
parameter is the same length as the first, then each value of the first
parameter is repeated the number of times designated by the corresponding value
of the second. Visualized:
```{r echo=FALSE, results='asis'}
"<pre><code>", paste0(c(" ", bb.html), collapse=" "), "</code></pre>"
) )
This allows us to generate the bounding box ids with internally vectorized code,
despite the bounding boxes having different sizes.
Next we compute the x offsets of each pixel in the bounding box:
## We wish to cycle through the x values in a bounding box
## as many times as there are y values in that bounding box
x.lens <- rep(dims.x0, dims.y0) # special use of `rep` again here
## sequence generates sequences of sequences... <- sequence(x.lens)
Now we have all the possible x offset values for each bounding box:
```{r echo=FALSE, results='asis'}
paste0(c(" ", bb.html), collapse=" "), "\n",
paste0(c(" ",, collapse=" "),
) )
The y offsets are a small variation on the theme:
```{r} <- rep(sequence(dims.y0), x.lens)
Everything together:
```{r echo=FALSE, results='asis'}
paste0(c(" ", bb.html), collapse=" "), "\n",
paste0(c(" ",, collapse=" "), "\n",
paste0(c(" ",, collapse=" "),
) )
And visualized in two dimensions:
```{r echo=FALSE, fig.height=2.5}
col <- c('#0000FF', '#00DD00', '#CCCC00')
par(mfrow=c(1,3), mar=c(4, 4, 4, 0), pty='s')
for(i in unique( {
plot([ == i],[ == i], col=col[i],
xlim=c(.8,3.2), ylim=c(.8,3.2), xaxp=c(1,3,2), yaxp=c(1,3,2),
xlab='', ylab='', asp=1, bty='n', pch=15, cex=3
title(sprintf(" %d", i))
Unfortunately `sequence` is not an internally vectorized function. It is just a
looped call to `seq_len`:
This becomes a liability when we have many sequences as with the full set of
`volcano` polygons. Fortunately it turns out that it is possible to write an
internally vectorized version of `sequence` by taking advantage of the
internally vectorized cumulative functions:
## Start with one continuous sequence
x.lens <- rep(dims.x0, dims.y0)
seq <- seq_len(sum(x.lens))
## Compute the locations at which the sequence should reset
reset <- cumsum(x.lens[-length(x.lens)]) + 1L
## Generate a reset vector
sub.raw <- integer(length(seq))
sub.raw[reset] <- reset - 1L
sub <- cummax(sub.raw)
## Modify the sequence by subtracting the reset vector
seq - cummax(sub)
all.equal(, seq - cummax(sub))
This produces the same result as `sequence`. For convenience we put the above
logic into `shadow::sequence2`. Let's benchmark against the full polygon data:
```{r, eval=FALSE}
x.lens <- rep(dims.x, dims.y) # full volcano bounding boxes
microbenchmark(sequence2(x.lens), sequence(x.lens))
Unit: milliseconds
expr min lq mean median uq max neval
sequence2(x.lens) 13.8 15.9 31.4 22.7 30.8 177 100
sequence(x.lens) 68.8 94.5 139.0 123.2 155.7 363 100
We see a 4-5x improvement despite having to rely on far more function calls.
Normally we wouldn't highlight such internal details, but this is an interesting
use case of the internally vectorized cumulative functions to solve a problem
that doesn't have a built-in vectorized solution.
Now that we know how we can to generate the candidate x/y offsets quickly we can
fill out the bounding boxes.
dims.xy <- dims.x * dims.y <- rep(seq_along(dims.x), dims.xy)
x.lens <- rep(dims.x, dims.y) <- sequence2(x.lens) - 1L <- rep(sequence(dims.y), x.lens) - 1L
## Add back the min values of the bounding boxes
p.x <- + mins[['x']][]
p.y <- + mins[['y']][]
Finally, we're ready to generate our first pixel-level volcano. We will re-use
the textures we computed for our [mesh](#mesh-tri-persp):
p.raster <- matrix(numeric(0), x.width, y.width)
p.raster[cbind(p.x, p.y)] <- rep(texture, dims.xy)
We apply the textures by indexing `bb.raster` with a matrix (`cbind(p.x,
p.y)`). This is a special type of indexing that allows you to directly address
specific cells by providing a matrix with the row and column values you are
There is also a subtle trick deployed here: instead of checking for duplicate
pixels we just wrote all the pixels in increasing Z order. This ensures that
the pixels nearest to the viewer overwrite those further away implicitly
eliminating the duplicates.
Finally, we can plot our matrix as a raster, but due to different conventions
about which directions the dimensions run in we need to flip our matrix:
flip <- function(x) t(x)[rev(seq_len(ncol(x))),]
A volcano to make NovaLogic[^4] proud!
## Step 3: Pixel Trimming
style='float: left; margin: 5px 15px 5px 0;' width='250'/>
We are now worse off than we were when plotting the polygons in our
[mesh](#mesh-tri-persp), but fear not, we can improve the chunky raster volcano.
We need to trim the excess pixels in our bounding boxes. This is as simple as
deciding which pixels are inside the original polygons and which ones are not.
Really something that a three year old could do. Yet, as is often
the case, expressing such tasks in a way that a computer can easily execute is
There are many possible solutions to the problem, but I settled on using
[barycentric coordinates][3]. These have the useful property that points with
negative barycentric coordinates are outside of the triangle referenced by those
coordinates. Additionally, there is a simple arithmetic expression for
barycentric coordinates in triangles that can be internally vectorized.
Finally, there is a substantial side-benefit that we will discuss when we get to
We are going to re-use the [`bary_L`][4] function we implemented as part of our
[previous blog post][5] on list data structures and vertex shading. See that
blog post for a detailed explanation of what the function does. For now, all we
care about is which points end up with negative barycentric coordinates.
`shadow::barycentric` is a copy of that function.
First we need to match up our vertex data to the candidate points so we can use
vectorized operations. Unfortunately, this means that we have to make a copy of
the vertex x-y coordinates for every single point. This is wasteful, but one of
the trade offs often made in R: use more memory in order to make fewer expensive
R function calls.
cols <- c('x','y','t') # don't need Z vals as `` already Z ordered
mesh.all <- mesh[, cols]
## each id in `` reference a mesh triangle
mesh.all[] <- lapply(mesh.all, '[',
Now we have as many polygons as points, which allows us to compute the
barycentric coordinates in one R function call:
bc <- barycentric(list(x=p.x, y=p.y), mesh.all[,c('x','y')])
Of the five points `str` shows, only point 3 will be kept as it is the only one
with positive barycentric coordinates for each vertex.
## determine which points are inbounds and generate raster matrix
inbounds <- Reduce('&', lapply(bc, '>=', 0)) <- p.x[inbounds] <- p.y[inbounds]
p.raster <- matrix(numeric(0), x.width, y.width)
p.raster[cbind(,] <- rep(texture, dims.xy)[inbounds]
plot(as.raster(flip(p.raster)), asp=1)
Great! All that work and we're right back were we started[^5] with our
[triangle mesh](#mesh-tri-persp). But isn't it so satisfying that we built this
thing pixel by pixel, with our bare hands. It's a little bit like that pasta
sauce you might make from scratch, that objectively really isn't that much
better than the stuff in the jar, but it sure does _seem_ to taste better.
Oh, and we're not done yet.
# Shading
The main reason I picked barycentric coordinates to determine which points are
inbounds is that those coordinates also represent the relative weights at each
vertex of the triangle required to "balance" the triangle on the corresponding
points. This will allow us to shade the pixels within polygons based on their
proximity to each vertex and the vertex colors. I cover this at length in the
previously mentioned [blogpost][5].
Since we have the barycentric coordinates, all we have to do is compute weighted
average colors for each pixel:
```{r} <- lapply(bc, '[', inbounds) <- lapply(mesh.all[, 't'], '[', inbounds)
t.shaded <- Reduce('+', Map('*',,
p.raster <- matrix(numeric(0), x.width, y.width)
p.raster[cbind(,] <- t.shaded
And just like that we have a beautifully shaded volcano!
# Interlude II
Hah, you probably thought we were done. Well, we're close. But first,
we need to run some benchmarks. I captured the entire process we went through
in `shadow::render` to make it easier to re-use and benchmark:
```{r, eval=FALSE}
render_elevation(volcano, shadow, rot, resolution=600, d=.5)
Ticks: 3743; Iterations: 10; Time Per: 520.3 milliseconds; Time Total: 5.203 seconds; Time Ticks: 3.743
render_elevation --------------------- : 520.3 - 0.0
rasterize ------------------------ : 518.2 - 8.9
barycentric ------------------ : 132.8 - 132.8
mesh_expand ------------------ : 120.7 - 0.0
| lapply ------------------- : 120.7 - 120.7
candidate_pixels ------------- : 102.7 - 69.5
| sequence2 ---------------- : 33.2 - 12.1
| cummax --------------- : 16.0 - 16.0
| integer -------------- : 4.9 - 4.9
lapply ----------------------- : 55.9 - 55.9
shade ------------------------ : 48.7 - 0.0
| Reduce ------------------- : 48.7 - 0.0
| Map ------------------ : 42.7 - 0.0
| | mapply ----------- : 42.7 - 0.0
| | lapply ------- : 34.1 - 34.1
| | <Anonymous> -- : 8.6 - 8.6
| f -------------------- : 6.0 - 6.0
Reduce ----------------------- : 45.0 - 0.7
f ------------------------ : 23.1 - 23.1
lapply ------------------- : 21.3 - 0.0
FUN ------------------ : 21.1 - 21.1
Almost 2FPS! We're not going to be playing any 3D games with this, but
considering this is done on a 2016 Macbook (yes, the 12" Macbook a.k.a. Apple
Netbook), on CPU, with a single core, in double precision floating
point[^double-precision], and via an interpreted programming language, it's not
too bad.
As for the benchmark details: virtually all the time was spent rasterizing. 3D
projection, meshes, and perspective adjustments take no time, relatively
speaking. So why bother with rasterization? One obvious answer is that it
allows us to implement much better shading. We also want the bitmap data
as an R object for further manipulation, and while it is possible to save
the polygon output to image files, I don't know of an easy way to get that
information back into R with base-R-only. And finally polygon rendering becomes
expensive when the polygon counts increase. Here we illustrate with an
elevation map that has fifty times as many points as `volcano`:
```{r eval=FALSE}
# File originally from
elev2 <- readRDS('static/data/three-d-pipeline-elev-complex.RDS')
# shadow2 <- ray_shade2(elev2, seq(-90, 90, length=25), 45)
shadow2 <- readRDS('static/data/three-d-pipeline-shadow-complex.RDS')
rot2 <- rot_x(-90) %*% rot_z(-90)
rel <- render_elevation(elev2, shadow2, rot2, 1000, 2, zord='pixel', empty=1)
title('Raster Render', cex=1.5)
user system elapsed
1.956 0.649 2.622
```{r eval=FALSE}
mesh2 <- elevation_as_mesh(elev2, shadow2, rot_x(-90) %*% rot_z(-90), d=2)
x <-, c(mesh2[,'x'], list(NA)))
y <-, c(mesh2[,'y'], list(NA)))
texture <- gray((Reduce('+', mesh2[,'t'])/nrow(mesh2)))
polygon(rescale(x), rescale(y), col=texture, border=texture)
title('Polygon Render', cex=5)
user system elapsed
10.389 0.238 10.664
```{r echo=FALSE, eval=FALSE}
'static/images/three-d-pipeline-pixel.png', width=1000, height=601,
title('Raster Render', cex=1.5)
'static/images/three-d-pipeline-polygon.png', width=1000, height=601,
polygon(rescale(x), rescale(y), col=texture, border=texture)
title('Polygon Render', cex=5)
It is unclear to me why the polygon method is slower[^6] given that the
graphical device will rasterize the polygons as well, presumably in compiled
code. With the higher polygon count, the polygon image quality is higher, but
if you compare the two renders closely you look closely you can still see the
polygon artifacts:
<a title='click for full size' href='/images/three-d-pipeline-pixel.png' >
alt='Raster based render of elevation map'
style='max-width: 600px; display: block; margin: auto' />
<a title='click for full size' href='/images/three-d-pipeline-polygon.png' >
alt='Polygon based render of elevation map'
style='max-width: 600px; display: block; margin: auto' />
# Stereoscopy
Finally: the fun stuff. As hinted at in the introduction, we are looking to
make [analygraphs][8] and [stereograms][9]. We can easily render two scenes
with slightly different angles. For convenience, we use
`shadow::mrender_elevation` which is a wrapper to `render_elevation` for
multiple rotations. In addition to supporting multiple rotations,
`mrender_elevation` uses absolute distances and fixed fields of view to ensure
consistent perspective transformations.
rot.l <- rot %*% rot_z(2.5)
rot.r <- rot %*% rot_z(-2.5)
elren <- mrender_elevation(
volcano, shadow, list(rot.l, rot.r), res=1000, d=125, fov=85
left <- flip(elren[[1]])
right <- flip(elren[[2]])
In order to generate an analygraph we need to assign the left image to the red
channel, and the right image to the cyan channel. No cyan channel? No problem:
turns out that cyan is just blue and green added together. This is very easy to
do as the R raster format accepts 3D arrays as inputs, where each value of the
third dimension corresponds to one of the RGB color channels:
```{r analygraph}
analygraph <- function(left, right) {
analygraph <- array(0, dim=c(dim(left), 3))
analygraph[,,1] <- left # red channel
analygraph[,,2:3] <- right # green and blue
elcolor <- analygraph(left, right)
par(bg='black', mai=numeric(4))
shadow::analygraph_glasses(.15) # for the icon h/t @jurbane2
<div style='display: none;'>
```{r analygraph-square, echo=FALSE, fig.width=4, fig.height=4}
# living on the edge here, we should really re-render and set the empty color
# to -1 or something, but we're lazy. This is for the opening thumbnail.
par(bg='transparent', mai=numeric(4))
elcolor[elcolor == 0] <- 1
You can get [analygraph glasses very cheaply][16]. I got ten for five bucks.
One advantage of using red and cyan is that the areas where the images overlap
are white or gray, which looks a lot more natural. If we had used red and green
the overlapping areas would end up yellow.
If you can't be bothered getting the glasses, we can also make stereograms:
```{r stereogram}
eljuxt <- cbind(left, right)
par(bg='black', mai=c(.25, .5, .25, .5))
In order to see the 3D effect on these you need to get your [eyes to
diverge][10] such that each eye is looking at one of the images. It does take
practice to master. One trick that can help is making the image smaller so that
the volcanoes are closer together. If you wear glasses you might also find it
easier to "focus" on the stereogram with them off. The effort is worth it as
the perspective effect is much stronger with stereograms than with analygraphs.
Since nowadays it seems no blog post is complete without some kind of animation,
let's spin up `volcano`. Out of courtesy to those of you who get nauseated
easily, I'm not auto-playing the videos. I also set up the lighting to remain
constant throughout so that it is clear that our frame of reference is fixed.
Tough luck to those visiting `volcano` today.
<video style='display: block; margin: 0 auto;' controls loop>
<source src="/videos/three-d-pipeline-volcano-rot.mp4" type="video/mp4" />
And one more, why not?
<video style='display: block; margin: 0 auto;' controls loop>
src="/videos/three-d-pipeline-volcano-zoom.mp4" type="video/mp4" />
Assembling the video from the component frames with [`ffmpeg` magic](#videos)
and using blogdown to render the post are the only things we resorted to
non base R for.
# So Long, and Thanks For All The Clicks
Hopefully this was somewhat instructive and enjoyable for you. I imagine if you
got this far there is a chance it was. Or maybe you're just looking for my
contact info to flame me on the interwebs. If the latter I'll take some solace
in the fact that you probably overshot and had to scroll back through the
footnotes and appendix.
All the illustrations in this post were generated with base R graphics. You can
see the [source for this post on Github][14].
Questions? Comments? @ me on [Twitter](
# Appendix
## Session Info
## Helper Functions
```{r helper-funs, eval=FALSE}
## Meshes With Persp
With some work I was able to roughly recreate our own [tile
x=seq_len(nrow(volcano)), y=seq_len(ncol(volcano)), z=volcano,
ylim=c(2,52), xlim=c(6,46), theta=-65, phi=70,
col='#DDDDDD', border='#444444', r=1e7, scale=FALSE
I left the box in on purpose to highlight that I had to go to ridiculous lengths
to get `persp` to fit the viewport the way I wanted. Finding the correct `xlim`
and `ylim` values to achieve the desired effect was an exercise in almost futile
stochastic debugging. There probably is a way to do this better.
Notice also that the `theta` and `phi` values are kind of weird given that we
use `rot_x(-20) %*% rot_z(65)`. I wasn't able to figure out what values I
needed until I went through the C source code to figure out what what order
`phi` and `theta` were being applied (it matters), and saw this:
## r-devel@75683:src/library/graphics/src/plot3d.c@1205
XRotate(-90.0); /* rotate x-y plane to horizontal */
YRotate(-theta); /* azimuthal rotation */
XRotate(phi); /* elevation rotation */
I had similar frustrations in understanding `rgl` as well, and similarly had to
dig through sources to figure exactly what was going on and eventually decided
it was more fun to just build the whole thing from scratch to best fit the
specific purpose I had in mind.
As far as I was able to get with `persp` is this:
tile.col <- Reduce('+', mesh.tile[,'t']) / 4
x=seq_len(nrow(volcano)), y=seq_len(ncol(volcano)), z=volcano,
ylim=c(2,52), xlim=c(6,46), theta=-65, phi=70,
col=gray(tile.col), border=NA, r=1e7, scale=FALSE, box=FALSE
This doesn't look great: it's blocky, and there is no [good way to
get rid of the anti-aliasing lines][11]. I also used our `mesh.tile`
list-matrix to compute the shading of the tiles, which requires a fair bit of
work beyond what `persp` does on its own. Finally, while I can get this
out to a bitmap, AFAICT there is no easy base-R-only way to read back a bitmap
into R for manipulation to make the analygraphs.
## Videos
Both videos were made by generating PNG files that I turned into an MP4 with
ffmpeg` with the command:
ffmpeg -pattern_type glob -i '*.png' -r 30 -pix_fmt yuv420p
```{r eval=FALSE}
Analygraph videos:
## Basic video
deg <- 0:119 * 3
deg.stereo <- c(deg + 2.5, deg - 2.5)
rot.360 <- lapply(deg.stereo, function(x) rot %*% rot_z(x))
sun.els <- seq(-30, 120, length=50)
shadows <-
lapply(deg + 180, ray_shade2, anglebreaks=sun.els, heightmap=volcano)
elren <- mrender_elevation(volcano, shadows, rot.360, res=600, d=200)
elren <- lapply(elren, flip)
dim(elren) <- c(length(elren) / 2, 2) # a list matrix again
elrows <- seq_len(nrow(elren))
elren2 <- lapply(elrows, function(i), elren[i,]))
size <- max(dim(elren2[[1]]))
for(i in elrows) {
png(sprintf('tmp/volc%03d.png', i), width=size, height=size, bg='black')
## Coming at ya video
eye.sep <- tan(2.5/180*pi) * 200 # really half eye sep, exaggerate angle
dists <- seq(1400, 50, length.out=90)
deg.d <- atan(eye.sep / dists) * 180 / pi
deg.ds <- c(deg.d, -deg.d)
rot.d <- lapply(deg.ds, function(x) rot %*% rot_z(x))
elrena <- mrender_elevation(volcano, shadow, rot.d, res=600, d=as.list(dists))
elrena <- lapply(elrena, flip)
dim(elrena) <- c(length(elrena) / 2, 2) # a list matrix again
elrows <- seq_len(nrow(elrena))
elrena2 <- lapply(elrows, function(i), elrena[i,]))
size <- max(dim(elrena2[[1]]))
for(i in elrows) {
png(sprintf('tmp1/volc%03d.png', i), width=size, height=size, bg='black')
Stereogram video:
```{r eval=FALSE}
elren3 <- lapply(elrows, function(i), elren[i,]))
width <- ncol(elren3[[1]]) / 2
height <- nrow(elren3[[1]]) / 2
for(i in elrows) {
path <- sprintf('tmp2/volc%03d.png', i)
png(path, width=600, height=600/width*height, bg='black')
[^quixotic]: Why would I subject myself to this? Well, I'm a sucker for
learning new interesting things such as how 3D rendering is actually done, and
there is no better way to learn than to do. I also believe, as an R package
developer, that it is important to know just how far you can push base
R (there are good reasons for this that I will elaborate on in a future blog
post), and there is no better way to know than to push it to do
things it isn't really meant to do.
[^base-r-only]: I use the term loosely to mean the packages that are by default
loaded into the search path of a "vanilla" R session. As of this writing the
packages are `stats`, `graphics`, `grDevices`, `utils`, `datasets`, `methods`,
and the actual `base` package.
[^copies-for-speed]: There is always a tension in R between minimizing R-level
function calls and minimizing memory footprint. Due to the cost of calling
R-level functions it is often though not always useful to duplicate or
partially duplicate data if it avoids many R-level function calls.
[^how-fast-is-useful]: We're not shooting for live rendering of 3D games here;
something that will allow us to render short video clips in single digit
minutes is what we consider useful.
[^internally-vectorized]: R code that carries out looped operations over vectors
in compiled code rather than in R-level code.
[^double-precision]: 3D rendering is typically done in single precision floating
point, which would halve the memory footprint of the process. Granted we only
have on color chanel, but on the whole even if we had three color chanels
we're taking a performance penalty due to R only have one type of floating
point data type.
[^1]: Some 3D packages out there will blatantly lie to you by rotating the
displayed axes along with the object, making it seem like the axes are
rotating, but in reality the underlying axes are not rotating. This lie
made it particularly difficult for me to figure out why things were not
rotating the way I thought they should. I'm still pissed off about it.
[^2]: In reality we only *need* to do this for the polygons that are visible in
the end result, but determining what is visible and not is non-trivial,
particularly for open surfaces at the pixel level, so we'll reserve that for a
possible future low-probability-of-getting-made blog post.
[^3]: A more efficient variation on `expand.grid(x, y)` would be
`list(rep(seq_len(x),y), rep(seq_len(y), each=x))`
[^4]: NovaLogic developed a 3D terrain rendering engine for use in the Comanche
helicopter flight simulations back in the early nineties. It was awesome
then, but looks about as good as our volcano here does right now:
![Comanche Screenshot][7]<br />See [s-macke][6] for an amazing explanation of
their rendering engine.
[^5]: Well, almost, the eagle-eyed reader will notice we've lost the
anti-aliasing, but shhhh, don't tell anyone.
[^6]: Possibly it is because the device still renders the polygon borders,
or perhaps because of anti-aliasing (though turning that off did not seem to
make a difference). I am not curious enough to look at the source of the
device code.
[2]: /2018/11/23/is-your-matrix-running-slow-try-lists#lists-to-the-rescue
[4]: /2018/11/23/is-your-matrix-running-slow-try-lists/#bary_l
[5]: /2018/11/23/is-your-matrix-running-slow-try-lists/#a-case-study