## use in slides:

with(iris,
     ecloud(x = jitter(Petal.Length),
            y = jitter(Petal.Width),
            z = jitter(Sepal.Length),
            xlab = "Petal Length",
            ylab = "Petal Width",
            zlab = "Sepal Length",
            alpha.density = 0.6,
            col = rainbow(3)[Species]))

vp <-
    structure(c(0.97834062576294, -0.170206964015961, 0.117810525000095, 
                0, 0.204019293189049, 0.696589708328247, -0.687851071357727, 
                0, 0.0350112244486809, 0.696988046169281, 0.716227293014526, 
                0, 0, 0, 0, 1), .Dim = c(4, 4))

view3d(userMatrix = vp)
rgl.snapshot("iris3d.png", fmt = "png")
rgl.postscript("iris3d.pdf", fmt = "pdf")





with(iris,
     ecloud(x = jitter(Petal.Length),
            y = jitter(Petal.Width),
            z = jitter(Sepal.Length),
            xlab = "Petal Length",
            ylab = "Petal Width",
            zlab = "Sepal Length",
            aspect = c(0.5, 2),
            alpha.density = 0.8,
            col = rainbow(3)[Species]))

with(quakes,
     ecloud(x = long,
            y = lat,
            z = depth,
            radius = 0.001))

for (i in seq(1, 360, length = 36)) rgl.viewpoint(i,i/4);



ecloud <-
    function(x, y, z,
             aspect = c(1, 1),
             xlab = deparse(substitute(x)),
             ylab = deparse(substitute(y)),
             zlab = deparse(substitute(z)),
             xlim = extend(range(x)),
             ylim = extend(range(y)),
             zlim = extend(range(z)),
             bg = "white",
             col = "pink",
             col.box = "turquoise",
             col.density = "grey",
             col.axis = "black",
             alpha.box = 0.1,
             alpha.density = 0.8,
             radius = 0.01,
             dh = 0.2,
             new = FALSE)
{
    stopifnot(require(rgl))
    aspect <- rep(aspect, length = 2)
    if (new) open3d()
    clear3d()
    if (!is.null(bg)) bg3d(color = bg)
    extend <- function(r) r + 0.05 * c(-1, 1) * diff(r)
    tx <- (x - xlim[1]) / diff(xlim)
    ty <- aspect[1] * (y - ylim[1]) / diff(ylim)
    tz <- aspect[2] * (z - zlim[1]) / diff(zlim)
    spheres3d(x = tx, y = ty, z = tz,
              radius = radius,
              color = col)
    bbox3d(color = c(col.box, col.axis),
           xat = (pretty(xlim, 3) - xlim[1]) / diff(xlim),
           yat = aspect[1] * (pretty(ylim, 3) - ylim[1]) / diff(ylim),
           zat = aspect[2] * (pretty(zlim, 3) - zlim[1]) / diff(zlim),
           xlab = format(pretty(xlim, 3)),
           ylab = format(pretty(ylim, 3)),
           zlab = format(pretty(zlim, 3)),
           alpha = c(alpha.box, 1))
    text3d(x = c(0.5,   0,   0),
           y = c(  0, 0.5,   0),
           z = c(  0,   0, 0.5),
           text = c(xlab, ylab, zlab),
           color = "black")
    require(MASS)
    dxy <- kde2d(tx, ty)
    dxy$z[] <- aspect[2] + dh * dxy$z / max(dxy$z)
    rgl.surface(x = dxy$x,
                y = dxy$z,
                z = dxy$y,
                coords = c(1, 3,  2),
                alpha = alpha.density)
    dyz <- kde2d(ty, tz)
    dyz$z[] <- 1 + dh * dyz$z / max(dyz$z)
    rgl.surface(x = dyz$x,
                y = dyz$z,
                z = dyz$y,
                coords = c(2, 1,  3),
                alpha = alpha.density)
    dzx <- kde2d(tz, tx)
    dzx$z[] <- aspect[1] + dh * dzx$z / max(dzx$z)
    rgl.surface(x = dzx$x,
                y = dzx$z,
                z = dzx$y,
                coords = c(3, 2,  1),
                alpha = alpha.density)
    invisible()
}