[R] 3-D response surface using wireframe()

[Deepayan Sarkar]

On Fri, Apr 9, 2010 at 1:46 PM, array chip <arrayprofile at yahoo.com> wrote:
> David,
>
> Thanks for the 2 previous posts from Sarkar. Actually, I am now one step closer. I am now able
> to remove the 3 outer lines of the bounding box using par.box argument, even Sarkar said in
> his 2008 post that par.box() does not control different boundaries, so maybe it was fixed.
>
> Replacing "par.box=list(lwd=2)" in my original code with "par.box=list(lwd=2,
> col=c(1,1,1,NA,1,1,NA,NA,1))" will now remove the 3 outer lines of the bounding box. The
> only thing missing here is the 3 inner lines of the box (behind the plot) are dashed lines,
> not solid. And par.box argument only control those 9 visible lines of the bounding box.
>
> As for how to draw grid lines onto the 3 surfaces, I still have no clue. But as you pointed
> out Sarkar indicated in his 2007 post that it might be possible.

Sorry I haven't had time to take a stab at this till now.

As David has pointed out, the trick is to get 3d versions of the grid
lines, and project them using ltransform3dto3d(), just as in the
projected contour lines example. The code given below is tedious but
conceptually straightforward. It may need some manual tweaking to get
the right faces, depending on the value of 'screen'.

 rescale <- function(x, from, to)
{
    ans <- min(to) + abs(diff(to)) * (x - min(from)) / abs(diff(from))
    ans[ans > min(to) & ans < max(to)]
}

panel.3dgrid <-
    function(rot.mat, distance,
             xlim, ylim, zlim,
             xlim.scaled, ylim.scaled, zlim.scaled,
             zero.scaled,
             ...)
{
    xgrid <- rescale(pretty(xlim), xlim, xlim.scaled)
    ygrid <- rescale(pretty(ylim), ylim, ylim.scaled)
    zgrid <- rescale(pretty(zlim), zlim, zlim.scaled)
    ## helper functions
    drawGridSegments <- function(xyz1, xyz2)
    {
        m1 <- ltransform3dto3d(t(xyz1), rot.mat, distance)
        m2 <- ltransform3dto3d(t(xyz2), rot.mat, distance)
        panel.segments(m1[1,], m1[2,], m2[1,], m2[2,], col = "grey")
    }
    drawBox <- function(xyz1, xyz2)
    {
        m1 <- ltransform3dto3d(t(xyz1), rot.mat, distance)
        m2 <- ltransform3dto3d(t(xyz2), rot.mat, distance)
        panel.segments(m1[1,], m1[2,], m2[1,], m2[2,], col = "black")
    }
    ## xgrid * ygrid * min(z) (bottom)
    drawGridSegments(expand.grid(x = xgrid, y = ylim.scaled[1], z =
zlim.scaled[1]),
                     expand.grid(x = xgrid, y = ylim.scaled[2], z =
zlim.scaled[1]))
    drawGridSegments(expand.grid(x = xlim.scaled[1], y = ygrid, z =
zlim.scaled[1]),
                     expand.grid(x = xlim.scaled[2], y = ygrid, z =
zlim.scaled[1]))
    ## max(x) * ygrid * zgrid
    drawGridSegments(expand.grid(x = xlim.scaled[2], y = ygrid, z =
zlim.scaled[1]),
                     expand.grid(x = xlim.scaled[2], y = ygrid, z =
zlim.scaled[2]))
    drawGridSegments(expand.grid(x = xlim.scaled[2], y =
ylim.scaled[1], z = zgrid),
                     expand.grid(x = xlim.scaled[2], y =
ylim.scaled[2], z = zgrid))
    ## xgrid * max(y) * zgrid
    drawGridSegments(expand.grid(x = xgrid, y = ylim.scaled[2], z =
zlim.scaled[1]),
                     expand.grid(x = xgrid, y = ylim.scaled[2], z =
zlim.scaled[2]))
    drawGridSegments(expand.grid(x = xlim.scaled[1], y =
ylim.scaled[2], z = zgrid),
                     expand.grid(x = xlim.scaled[2], y =
ylim.scaled[2], z = zgrid))
    ## boxes around the faces
    drawBox(expand.grid(x = xlim.scaled[1], y = ylim.scaled, z =
zlim.scaled[1]),
            expand.grid(x = xlim.scaled[2], y = ylim.scaled, z =
zlim.scaled[1]))
    drawBox(expand.grid(x = xlim.scaled, y = ylim.scaled[1], z =
zlim.scaled[1]),
            expand.grid(x = xlim.scaled, y = ylim.scaled[2], z =
zlim.scaled[1]))
    drawBox(expand.grid(x = xlim.scaled[2], y = ylim.scaled, z =
zlim.scaled[1]),
            expand.grid(x = xlim.scaled[2], y = ylim.scaled, z =
zlim.scaled[2]))
    drawBox(expand.grid(x = xlim.scaled[2], y = ylim.scaled[1], z =
zlim.scaled),
            expand.grid(x = xlim.scaled[2], y = ylim.scaled[2], z =
zlim.scaled))
    drawBox(expand.grid(x = xlim.scaled[1], y = ylim.scaled[2], z =
zlim.scaled),
            expand.grid(x = xlim.scaled[2], y = ylim.scaled[2], z =
zlim.scaled))
    drawBox(expand.grid(x = xlim.scaled, y = ylim.scaled[2], z =
zlim.scaled[1]),
            expand.grid(x = xlim.scaled, y = ylim.scaled[2], z =
zlim.scaled[2]))
}

cloud(Sepal.Length ~ Petal.Length * Petal.Width,
      data = iris, cex = .8,
      groups = Species,
      par.box = list(col = "transparent"),
      screen = list(z = 20, x = -70, y = 3),
      panel.3d.cloud = function(...) {
          panel.3dgrid(...)
          panel.3dscatter(...)
      },
      par.settings = list(axis.line = list(col = "transparent")),
      scales = list(arrows = FALSE, col = "black"))


wireframe(volcano, shade = TRUE,
          par.box = list(col = "transparent"),
          screen = list(z = 20, x = -70, y = 3),
          panel.3d.wireframe = function(...) {
              panel.3dgrid(...)
              panel.3dwire(...)
          },
          par.settings = list(axis.line = list(col = "transparent")),
          scales = list(arrows = FALSE, col = "black"))