## R code to reproduce plots in Dai Feng, Luke Tierney (2008):
## "Computing and Displaying Isosurfaces in R," Journal of Statistical
## Software 28(1), URL http://www.jstatsoft.org/v28/i01/.

library("misc3d")
library("rgl")

## Figure 1: Locations of earthquake epicenters rendered using rgl.
## Locations:
open3d()
points3d(quakes$long/22,quakes$lat/28,-quakes$depth/640, size=2)
box3d(col="gray")
title3d(xlab="long",ylab="lat",zlab="depth")
## Add a density contour:
de <- kde3d(quakes$long, quakes$lat, -quakes$depth, n = 40)
contour3d(de$d, level = exp(-12), x = de$x/22, y = de$y/28, z = de$z/640,
          color = "green", color2 = "gray", add = TRUE)


## Figure 2: Density contour surface for variables Sepal.Length,
## Sepal.Width, and Petal.Length from the iris data set, with false
## color showing the levels of the fourth variable, Petal.Width,
## predicted by aloess fit.
de <- kde3d(iris[,1],iris[,2],iris[,3],n=40)
fit <- loess(Petal.Width ~ Sepal.Length + Sepal.Width + Petal.Length,
             data = iris, control = loess.control(surface = "direct"))
fitCols <- function(x, y, z) {
    d <- data.frame(Sepal.Length = x, Sepal.Width = y, Petal.Length = z)
    p <- predict(fit, d)
    k <- 32
    terrain.colors(k)[cut(p, k,levels=FALSE)]
}
contour3d(de$d, 0.1, de$x, de$y, de$z, color = fitCols)
box3d(col="gray")
title3d(xlab="Sepal Length",ylab="Sepal Width",zlab="Petal Length")


## Figure 3: Multiple isosurfaces of the density of a mixture of three
## tri-variate normal distributions rendered with partial
## transparency.
nmix3 <- function(x, y, z) {
    m <- 0.5
    s <- 0.5
    0.4 * dnorm(x, m, s) * dnorm(y, m, s) * dnorm(z, m, s) +
    0.3 * dnorm(x, -m, s) * dnorm(y, -m, s) * dnorm(z, -m, s) +
    0.3 * dnorm(x, m, s) * dnorm(y, -1.5 * m, s) * dnorm(z, m, s)
}
n <- 40; k <- 5; alo <- 0.1; ahi <- 0.5; cmap = heat.colors
lev <- seq(0.05, 0.2, len = k)
col <- rev(cmap(length(lev)))
al <- seq(alo, ahi, len = length(lev))
x <- seq(-2, 2, len = n)
## Rendered with rgl
contour3d(nmix3, lev, x = x, y = x, z = x, color = col, alpha = al)
## Rendered with grid graphics
alo <- 0.05; ahi <- 0.3;
al <- seq(alo, ahi, len = length(lev))
contour3d(nmix3, lev, x = x, y = x, z = x, color = col,
          alpha = al, engine = "grid")


## Figure 4: Isosurfaces of a normal mixture density rendered by
## standard graphics using a cutaway strategy to show the nested
## contours.
cmap <- rainbow
col <- rev(cmap(length(lev)))
m <- function(x,y,z) x > .25 | y < -.3
contour3d(nmix3, lev, x = x, y = x, z = x,
          color = col, color2 = "lightgray",
          mask = m, engine = "standard",
          scale = FALSE, screen=list(z = 130, x = -80))


## Figure 5: Two isosurfaces of a CT scan of an engine block.  Data
## ara available at
##     http://www.sph.sc.edu/comd/rorden/render.html#Sample
## in
##     http://www.sph.sc.edu/comd/rorden/engine.zip
library("AnalyzeFMRI")
Engine <- f.read.analyze.volume("engine.hdr")[,,,1]
contour3d(Engine, c(120, 200), color = c("lightblue", "red"),
          alpha=c(0.1, 1))


## Figure 6: Isosurfaces of a brain and two intensity differences
## between two tasks in a PET experiment.  The data (not publically
## available) in ANALYZE format consist of an MRI image
## "template.{hdr,img}" and an image of average differnces in PET
## intensities between the tasks, "tmap1-8.{hdr,img}".
template <- f.read.analyze.volume("template.img")
template <- aperm(template,c(1,3,2,4))[,,95:1,1]
brain <- contour3d(template, level = 10000, alpha = 0.3, draw = FALSE)
tm <- f.read.analyze.volume("tmap1-8.img")[,,95:1]
brainmask <- template > 10000
activ <- contour3d(ifelse(brainmask, tm, 0),level = 4:5,
                   color=c("red", "yellow"), alpha=c(0.3, 1),
                   draw = FALSE)
drawScene.rgl(c(list(brain), activ))


## Figure 12: Contour surface of kernel density estimate for the first
## three variables in the iris data set rendered using the four
## standard material settings.
xlim <- c(4, 8)
ylim <- c(1, 5)
zlim <- c(0, 7)
de <- kde3d(iris[,1], iris[,2], iris[,3], n = 40,
            lims = c(xlim, ylim, zlim))
opar <- par(mar=c(1,1,4,1), mfrow=c(2,2))
for (m in c("default", "dull", "metal", "shiny")) {
    contour3d(de$d, 0.1, de$x, de$y, de$z,color="lightblue",
              engine="standard", material = m)
    title(paste('material = "', m, '"', sep = ""))
}
par(opar)


## Figure 13: Contour surface of kernel density estimate for the first
## three variables in the iris data set rendered using no additional
## shading and using two iterations of shading.
opar <- par(mar=c(1,1,4,1))
## no shading
contour3d(de$d, 0.1, de$x, de$y, de$z, color="lightblue",
          engine="standard", smooth = 0)
## shading with smooth=2
contour3d(de$d, 0.1, de$x, de$y, de$z, color="lightblue",
          engine="standard", smooth = 2)
par(opar)


## Figure 14: Surface plot of the Maunga Whau volcano with and without
## depth cuing and using smooth = 3 level shading.
z <- 2 * volcano
x <- 10 * (1:nrow(z))
y <- 10 * (1:ncol(z))
vtri <- surfaceTriangles(x, y, z, color =
                         function(x,y,z) {
                             cols <- terrain.colors(diff(range(z)))
                             cols[z - min(z) + 1]
                         })
opar <- par(mar = rep(0, 4))
## no depth cuing
drawScene(updateTriangles(vtri, material = "default", smooth = 3),
          screen = list(x = 40, y = -40, z = -135), scale = FALSE)
## with depth cuing
drawScene(updateTriangles(vtri, material = "default", smooth = 3),
          screen = list(x = 40, y = -40, z = -135), scale = FALSE,
          depth = 0.3)
par(opar)


## Figure 15: A scene combining several triangular mesh objects.
Volcano <- vtri
## Create the Quakes contour and teapot objects:
de <- kde3d(quakes$long, quakes$lat, -quakes$depth, n = 40)
Quakes <- contour3d(de$d, level = exp(-12),
                    x = de$x/22, y = de$y/28, z = de$z/640,
                    color = "green", color2 = "gray", draw = FALSE)
data("teapot")
Teapot <- makeTriangles(teapot$vertices, teapot$edges,
                        color = "red", color2 = "lightblue")
## Convenience function for rotation of an object.
## lattice::ltransform3dMatrix() is used here; misc3d:::trans3dMat()
## could be used instead but is not public at this point.
r <- function(tris, screen) {
    M <- lattice::ltransform3dMatrix(screen)
    transformTriangles(tris, M)
}
## Reposition the Volcano and Teapot objects:
Volcano <- scaleTriangles(Volcano,x = 0.002, y = 0.002, z = 0.002)
Volcano <- translateTriangles(Volcano, z = -1.5, x = 7)
Teapot <- r(Teapot, list(x = 90, z = 180))
Teapot <- scaleTriangles(Teapot, x = 0.2, y = 0.2, z = 0.2)
Teapot <- translateTriangles(Teapot, x = 9.5, y = 0, z = -1)
## Draw the scene
opar <- par(mar = rep(0, 4))
drawScene(list(Volcano, Quakes, Teapot),scale = FALSE,
          screen = list(z = -40, x = -70), perspective = TRUE, depth = 0.3)
par(opar)


## Figure 16: Earthquake epicenters and contour surface of kernel
## density estimate with points represented by small tetrahedra and
## smooth=2 shading of the contour surface.
de <- kde3d(quakes$long, quakes$lat, -quakes$depth, n = 40)
v <- contour3d(de$d, exp(-12),de$x/22, de$y/28, de$z/640,
               color="green", color2="gray", draw=FALSE, smooth = 2)
p <- pointsTetrahedra(quakes$long/22, quakes$lat/28, -quakes$depth/640,
                      size = 0.005)
drawScene(list(v, p))


## Figure 17 Using wireframe() or persp() to provide axes and labels
## for a contour surface of a kernel density estimate for the iris
## data.
## wireframe
xlim <- c(4, 8)
ylim <- c(1, 5)
zlim <- c(0, 7)
de <- kde3d(iris[,1], iris[,2], iris[,3], n = 40,
            lims = c(xlim, ylim, zlim))
v <- contour3d(de$d, 0.1, de$x, de$y, de$z, color="lightblue", draw = FALSE)
library("lattice")
w <- wireframe(matrix(zlim[1], 2, 2) ~ rep(xlim,2)*rep(ylim,each=2),
               xlim = xlim, ylim = ylim, zlim = zlim,
               aspect=c(diff(ylim) / diff(xlim), diff(zlim) / diff(xlim)),
               xlab = "Sepal Length", ylab = "Sepal Width",
               zlab = "Petal Width", scales = list(arrows = FALSE),
          panel.3d.wireframe = function(x, y, z, rot.mat, distance,
                                        xlim.scaled, ylim.scaled,
                                        zlim.scaled, ...) {
              scale <- c(diff(xlim.scaled) / diff(xlim),
                         diff(ylim.scaled) / diff(ylim),
                         diff(zlim.scaled) / diff(zlim))
              shift <- c(mean(xlim.scaled) - mean(xlim) * scale[1],
                         mean(ylim.scaled) - mean(ylim) * scale[2],
                         mean(zlim.scaled) - mean(zlim) * scale[3])
	      P <- rbind(cbind(diag(scale), shift), c(0, 0, 0, 1))
              rot.mat <- rot.mat %*% P
              drawScene(v, screen = NULL, R.mat = rot.mat,
                        distance = distance, add = TRUE, scale = FALSE,
                        light = c(.5,0,1), engine = "grid")
          })
print(w)
## persp
M <- persp(xlim, ylim, matrix(zlim[1], 2, 2), theta = 30, phi = 30,
	   col = "lightgray", zlim = zlim, ticktype="detailed",
           scale=FALSE, d=4, xlab = "Sepal Length",
           ylab = "Sepal Width", zlab = "Petal Width")
drawScene(v, screen = NULL, R.mat = t(M), add = TRUE, scale = FALSE,
          light = c(.5, 0, 1))
trans3dz <- function (x, y, z, pmat) {
    tr <- cbind(x, y, z, 1) %*% pmat
    list(x = tr[, 1]/tr[, 4], y = tr[, 2]/tr[, 4], z = tr[, 3]/tr[, 4])
}
g <- as.matrix(expand.grid(x = 1:2, y = 1:2, z = 1:2))
b <- trans3dz(xlim[g[,1]], ylim[g[,2]], zlim[g[,3]], M)
fci <- which.max(b$z)
fc <- g[fci,]
coord2index <- function(coord)
    as.integer((coord - 1) %*% c(1, 2, 4) + 1)
v1i <- coord2index(c(if (fc[1] == 1) 2 else 1, fc[2], fc[3]))
v2i <- coord2index(c(fc[1], if (fc[2] == 1) 2 else 1, fc[3]))
v3i <- coord2index(c(fc[1], fc[2], if (fc[3] == 1) 2 else 1))
segments(b$x[fci], b$y[fci], b$x[v1i], b$y[v1i], lty = 3)
segments(b$x[fci], b$y[fci], b$x[v2i], b$y[v2i], lty = 3)
segments(b$x[fci], b$y[fci], b$x[v3i], b$y[v3i], lty = 3)