In a break from my usual obsessions and interests here is a guest blog post by Ian Walker. I'm posting it because I think it is rather cool and hope it will be of interest to some of my regular readers. Ian is perhaps best known (in the blogosphere) for his work on transport psychology - particularly cycling - but is also an expert on psychological statistics.

Some time ago, I had some data that lent themselves to a three-dimensional surface plot. The problem was, the plot was quite asymmetrical, and finding the right viewing angle to see it effectively on a computer screen was extremely difficult. I spent ages tweaking angles and every possible view seemed to involve an unacceptable compromise.

Of course, displaying fundamentally three-dimensional items in two dimensions is an ancient problem, as any cartographer will tell you. That night, as I lay thinking in bed, a solution presented itself. I had recently been reading about the work of a fellow University of Bath researcher, Adrian Bowyer, and his RepRap project, to produce an open-source three-dimensional printer. The solution was obvious: I had to find a way to print R data on one of these printers!

I managed to meet up with Adrian back in May 2012, and he explained to me the structure of the STL (stereolithography) files commonly used for three-dimensional printing. These describe an object as a large series of triangles. I decided I'd have a go at writing R code to produce valid STL files.

I'm normally a terrible hacker when it comes to programming; I usually storm in and try to make things work as quickly as possible then fix all the mistakes later. This time, I was much more methodical. As a little lesson to us all, the methodical approach worked: I had the core code producing valid STL files in under 3 hours.

Unfortunately, it then took until September 2012 before I could get hold of somebody with a 3D printer who'd let me test my code. A few days ago the first prototype was produced, as you can see in this photograph:



So now I'd like to share the code under a Creative Commons BY-NC-SA licence, in case anybody else finds it useful. You can download the code here, in a file called r2stl.r. One day, when I learn how, I might try to make this a library, but for now you can just call this code with R's source() command. All that is in the file is the function r2stl(), and having once called the file with source(), you can then use the r2stl function to generate your STL files. The command is:

r2stl(x, y, z, filename='3d-R-object.stl', object.name='r2stl-object', z.expand=FALSE, min.height=0.008, show.persp=FALSE, strict.stl=FALSE)

• x, y and z should be vectors of numbers, exactly as with R's normal persp() plot. x and y represent a flat grid and z represents heights above this grid


• filename is pretty obvious, I hope


• object.name The STL file format requires the object that is being described to have a name specified inside the file. It's unlikely anybody will ever see this, so there's probably no point changing it from the default


• z.expand By default, r2stl() normalizes each axis so it runs from 0 to 1 (this is an attempt to give you an object that is agnostic with regard to how large it will eventually be printed). Normally, the code then rescales the z axis back down so its proportions relative to x and y are what they were originally. If, for some reason, you want your 3D plot to touch all six faces of the imaginary cube that surrounds it, set this parameter to TRUE


• min.height Your printed model would fall apart if some parts of it had z values of zero, as this would mean zero material is laid down in those parts of the plot. This parameter therefore provides a minimum height for the printed material. The default of 0.008 ensures that, when printed, no part of your object is thinner than around 0.5 mm, assuming that it is printed inside a 60 mm x 60 mm x 60 mm cube. Recall that the z axis gets scaled from 0 to 1. If you are printing a 60mm-tall object then a z-value of 1 represents 60mm. The formula is min.height=min.mm/overall.mm, so if we want a minimum printed thickness of 0.5mm and the overall height of your object will be 60mm, 0.5/60 = 0.008, which is the default. If you want the same minimum printed thickness of 0.5mm but want to print your object to 100mm, this parameter would be set to 0.5/100 = 0.005


• show.persp Do you want to see a persp() plot of this object on your screen as the STL is being generated? Default is FALSE


• strict.stl To make files smaller, this code cheats and simply describes the entire rectangular base of your object as two huge triangles. This seems to work fine for printing, but isn't strictly proper STL format. Set this to TRUE if you want the base of your object described as a large number of triangles and don't mind larger files

To view and test your STL files before you print them, you can use various programs. I have had good experiences with the free, open-source Meshlab, which even has iPhone and Android versions so you can let people interact with your data even when you're in the pub. Even if all you ever do is show people your 3D plots using Meshlab, I believe r2stl() still offers a useful service, as it makes viewing data far more interactive than static persp() plots. To actually get your hands on a printer, you might try your local school - apparently lots of schools have got rapid prototypers these days.


Demo


source('r2stl.r')

# Let's do the classic persp() demo plot, as shown in the photograph above

x <- seq(-10, 10, length= 100)

y <- x

f <- function(x,y) { r <- sqrt(x^2+y^2); 10 * sin(r)/r }

z <- outer(x, y, f)

z[is.na(z)] <- 1

r2stl(x, y, z, filename="lovelyfunction.stl", show.persp=TRUE)



# Now let's look at R's Volcano data

z <- volcano

x <- 1:dim(volcano)[1]

y <- 1:dim(volcano)[2]

r2stl(x, y, z, filename="volcano.stl", show.persp=TRUE)


I hope you might find this code useful. Any questions or suggestions, then please get in touch.


September 2012 - Ian Walker, Department of Psychology, University of Bath.