---
title: Including an interactive 3D `rgl` graphic in a html report with `knitr`
author: Stéphane Laurent
date : 2013-03-08
--- &lead

```{r setup0, echo=FALSE}
opts_chunk$set(fig.path="assets/fig/rglknitr-")
```

I posted [a question on StackOverflow](http://stackoverflow.com/questions/14879210/including-a-interactive-3d-figure-with-knitr): *how to embed an interactive 3D figure created with the `rgl` package into an html report created with `knitr` ?* The next day, Yihui Xie has posted a solution. He has updated `knitr` to include this possibility. 


The first example
--------------------------------

Below is the rendering of the example given by Yihui.

```{r setup}
library(rgl)
knit_hooks$set(webgl = hook_webgl)
```

```{r rgl-firstexample, webgl=TRUE, results='hide'}
x <- sort(rnorm(1000))
y <- rnorm(1000)
z <- rnorm(1000) + atan2(x,y)
open3d()
plot3d(x, y, z, col=rainbow(1000))
```




A real-life motivated example
-------------------------------

As a statistician consultant I sometimes have to write statistical reports about a  response surface  analysis.  Wouldn't it be great to include an interactive  response surface in a statistical report ? 
(At least, for fun --- I do not want to encourage an abusive practice of 3D visualization, even when it is interactive)

Russell V. Lenth's `rsm` package provides convenient functions  to draw the response surface fitted with the `lm()` funtion or the  `rsm()` function,  as shown in the vignette [Surface Plots in the rsm Package](http://cran.r-project.org/web/packages/rsm/vignettes/rsm-plots.pdf). 
Below is the interactive `rgl` version of the first example of the vignette.

```{r persp3dlm, echo=FALSE}
persp3d.lm <- function(x, form, at, bounds, zlim, zlab, xlabs, col = "white", 
    contours = NULL, hook, atpos = 3, decode = TRUE, theta = -25, 
    phi = 20, r = 4, border = NULL, box = TRUE, ticktype = "detailed", 
    ...) 
{
    draw.cont.line = function(line) {
        if (cont.varycol) {
            cont.col = col
            if (length(col) > 1) 
                cont.col = col[cut(c(line$level, dat$zlim), length(col))][1]
        }
        lines(trans3d(line$x, line$y, cont.z, transf), col = cont.col, 
            lwd = cont.lwd)
    }
    plot.data = contour.lm(x, form, at, bounds, zlim, xlabs, 
        atpos = atpos, plot.it = FALSE)
    transf = list()
    if (missing(zlab)) 
        zlab = ""
    facet.col = col
    cont = !is.null(contours)
    if (mode(contours) == "logical") 
        cont = contours
    cont.first = cont
    cont.z = cz = plot.data[[1]]$zlim[1]
    cont.col = 1
    cont.varycol = FALSE
    cont.lwd = 1
    if (is.character(contours)) {
        idx = charmatch(contours, c("top", "bottom", "colors"), 
            0)
        if (idx == 1) {
            cont.first = FALSE
            cont.z = plot.data[[1]]$zlim[2]
        }
        else if (idx == 2) {
        }
        else if (idx == 3) {
            cont.varycol = TRUE
            if (length(col) < 2) 
                col = rainbow(40)
        }
        else cont.col = contours
    }
    else if (is.list(contours)) {
        if (!is.null(contours$z)) 
            cz = contours$z
        if (is.numeric(cz)) 
            cont.z = cz
        else if (cz == "top") {
            cont.first = FALSE
            cont.z = plot.data[[1]]$zlim[2]
        }
        if (!is.null(contours$col)) 
            cont.col = contours$col
        if (!is.null(contours$lwd)) 
            cont.lwd = contours$lwd
        if (charmatch(cont.col, "colors", 0) == 1) {
            cont.varycol = TRUE
            if (length(col) < 2) 
                col = rainbow(40)
        }
    }
    for (i in 1:length(plot.data)) {
        dat = plot.data[[i]]
        cont.lines = NULL
        if (!missing(hook)) 
            if (!is.null(hook$pre.plot)) 
                hook$pre.plot(dat$labs)
        if (cont) 
            cont.lines = contourLines(dat$x, dat$y, dat$z)
        if (cont && cont.first) {
            transf = persp3d(dat$x, dat$y, dat$z, zlim = dat$zlim, 
                theta = theta, phi = phi, r = r, col = NA, border = NA, 
                box = FALSE, ...)
            lapply(cont.lines, draw.cont.line)
            par(new = TRUE)
        }
        if (length(col) > 1) {
            nrz = nrow(dat$z)
            ncz = ncol(dat$z)
            zfacet = dat$z[-1, -1] + dat$z[-1, -ncz] + dat$z[-nrz, 
                -1] + dat$z[-nrz, -ncz]
            zfacet = c(zfacet/4, dat$zlim)
            facet.col = cut(zfacet, length(col))
            facet.col = col[facet.col]
        }
        transf = persp3d(dat$x, dat$y, dat$z, xlab = dat$labs[1], 
            ylab = dat$labs[2], zlab = zlab, zlim = dat$zlim, 
            col = facet.col, border = border, box = box, theta = theta, 
            phi = phi, r = r, ticktype = ticktype, ...)
        if (atpos == 3) 
            #title(sub = dat$labs[5], ...)
        if (cont && !cont.first) 
            lapply(cont.lines, draw.cont.line)
        if (!missing(hook)) 
            if (!is.null(hook$post.plot)) 
                hook$post.plot(dat$labs)
        plot.data[[i]]$transf = transf
    }
#    plot.data
    invisible(plot.data)
}
```


```{r rgl-persp3dlm, message=FALSE, webgl=TRUE, error=FALSE, results='hide'}
library(rsm)
swiss2.lm <- lm(Fertility ~ poly(Agriculture, Education, degree=2), data=swiss)
open3d()
persp3d.lm(swiss2.lm, Education ~ Agriculture, zlab = "Fertility")
```


I have obtained the function `persp3d.lm()` I used above from the `persp.lm()` function of the 
`rsm` package by replacing every occurence of a call to the `persp()` function (`graphics` package)  with a call 
to the `persp3d()` function (`rgl` package). 
Thus this example is nothing but a default output whose aesthetics could be greatly improved 
with the possibilites of Adler & Murdoch's `rgl` package. 


