[R] persp(): add second plane

[Henrik Bengtsson]

Hmm, very interesting. I guess I am not the only one that thought about
extending Ben Bolkers idea to plot R^3 data by starting with an empty
persp(). Except from the "risk" of plotting outside the bounds of 3d-box
created by persp() I can't see any problems with this approach. I wrote up
some code for a plot3d() and a points3d() method. I have pasted the code at
the end and as you see I made the plot3d() code to be almost identical to
the one for persp(). First an example how the functions can be used and then
the code:

# Simulate three groups of (x,y,z) data
n <- 3*1000
gr <- list()
gr[[1]] <- matrix(rnorm(n, mean=c(0,0,1), sd=c(1,1,1)), ncol=3, byrow=TRUE)
gr[[2]] <- matrix(rnorm(n, mean=c(5,4,2), sd=c(1,2,1)), ncol=3, byrow=TRUE)
gr[[3]] <- matrix(rnorm(n, mean=c(2,2,4), sd=c(0.2,0.5,0.8)), ncol=3,
byrow=TRUE)

# Calculate the overall limits of the axis
xlim <- ylim <- zlim <- NA
for (k in 1:3) {
  gr[[k]] <- as.data.frame(gr[[k]])
  colnames(gr[[k]]) <- c("x", "y", "z")
  xlim <- range(xlim, gr[[k]][,1], na.rm=TRUE)
  ylim <- range(ylim, gr[[k]][,2], na.rm=TRUE)
  zlim <- range(zlim, gr[[k]][,3], na.rm=TRUE)
}

# First plot one group using plot3d()...
plot3d(gr[[1]], phi=50, pch=176, col="red", xlim=xlim, ylim=ylim, zlim=zlim,
xlab="x", ylab="y", zlab="z")

# ...and then use points3d() to plot the other groups
for (k in 2:3)
  points3d(gr[[k]], pch=176, col=c("blue", "green")[k-1])


##########################################################################
# plot3d
##########################################################################
plot3d <- function(x=seq(0, 1, len=nrow(z)), y=seq(0, 1, len=ncol(z)),
    z, xlim=range(x), ylim=range(y), zlim=range(z, na.rm=TRUE),
    xlab=NULL, ylab=NULL, zlab=NULL, main=NULL, sub=NULL, col=NULL,
    theta=0, phi=15, r=sqrt(3), d=1, scale=TRUE, expand=1,
    border=NULL, box=TRUE, axes=TRUE, nticks=5, ticktype="simple", ...) {

  # Sets the labels on the axis
  if (is.null(xlab))
    xlab <- if (!missing(x)) deparse(substitute(x)) else "X"
  if (is.null(ylab))
    ylab <- if (!missing(y)) deparse(substitute(y)) else "Y"
  if (is.null(zlab))
    zlab <- if (!missing(z)) deparse(substitute(z)) else "Z";

  # Extract the (x,y,z) coordinates
  if (missing(z)) {
    if (!missing(x)) {
      if (is.list(x)) {
        z <- x$z;
        y <- x$y;
        x <- x$x;
      } else {
        z <- x;
        x <- seq(0, 1, len=nrow(z));
      }
    } else
      stop("no `z' matrix specified");
  } else if (is.list(x)) {
    y <- x$y;
    x <- x$x;
  }

  # Decide on the axis ranges
  if (is.null(xlim))
    xlim <- range(x, na.rm=TRUE);
  if (is.null(ylim))
    ylim <- range(y, na.rm=TRUE);
  if (is.null(zlim))
    zlim <- range(z, na.rm=TRUE);

  # Some default values for persp(), which are not needed here.
  ltheta <- -135;
  lphi <- 0;
  shade <- NA;

  # Get the ticktype ID
  ticktype <- pmatch(ticktype, c("simple", "detailed"));

  # Create a minimum empty persp plot.
  xdummy <- 0:1;
  ydummy <- 0:1;
  zdummy <- matrix(rep(NA,4), nrow=2, ncol=2);

  r <- .Internal(persp(xdummy, ydummy, zdummy,
         xlim, ylim, zlim, theta, phi, r, d, scale, expand,
         col, border, ltheta, lphi, shade, box, axes, nticks,
         ticktype, xlab, ylab, zlab, ...));

  # Set any main or subtitle of the plot
  if (!is.null(main) || !is.null(sub))
    title(main = main, sub = sub, ...);

  # Plot the (x,y,z) points
  points3d(x,y,z, r, col=col, ...);

  # Remember the 4-by-4 perspective transform matrix
  options(persp.matrix=r);

  # ...and return it.
  invisible(r)
} # plot3d



##########################################################################
# points3d
##########################################################################
points3d <- function(x=seq(0, 1, len=nrow(z)), y=seq(0, 1, len=ncol(z)),
    z, persp.matrix=getOption("persp.matrix"), ...) {
  # Assert that the 4-by-4 perspective transform matrix is available
  if (is.null(persp.matrix))
    stop("Argument persp.matrix must be specified");

  if (!is.matrix(persp.matrix) ||
      nrow(persp.matrix) != 4 || ncol(persp.matrix) != 4) {
    stop("Argument persp.matrix must be a 4x4 matrix as returned by for
instance persp() or plot3d().");
  }

  # Extract the (x,y,z) coordinates
  if (missing(z)) {
    if (!missing(x)) {
      if (is.list(x)) {
        z <- x$z;
        y <- x$y;
        x <- x$x;
      } else {
        z <- x;
        x <- seq(0, 1, len=nrow(z));
      }
    } else
      stop("no `z' matrix specified");
  } else if (is.list(x)) {
    y <- x$y;
    x <- x$x;
  }

  # Ben Bolker's function for projecting (x,y,z) to (x',y')
  trans3d <- function(x,y,z,pmat) {
    tmat <- t((cbind(x,y,z,1)%*% pmat))
    list(x=tmat[1,]/tmat[4,],y=tmat[2,]/tmat[4,])
  }

  # Transform the (x,y,z) in R^3 to (x',y') in R^2
  xy3d <- trans3d(x,y,z, persp.matrix);
  points(xy3d, ...);
} # points3d


Cheers

Henrik Bengtsson

Dept. of Mathematical Statistics @ Centre for Mathematical Sciences
Lund Institute of Technology/Lund University, Sweden (+2h UTC)
Office: P316, +46 46 222 9611 (phone), +46 46 222 4623 (fax)
h b @ m a t h s . l t h . s e, http://www.maths.lth.se/bioinformatics/


> -----Original Message-----
> From: owner-r-help at stat.math.ethz.ch
> [mailto:owner-r-help at stat.math.ethz.ch]On Behalf Of Ben Bolker
> Sent: Thursday, March 14, 2002 10:16 PM
> To: Nina Lieske
> Cc: R help list
> Subject: Re: [R] persp(): add second plane
>
>
>
>   Not directly, but with a little bit of ingenuity you can hack one.
>
>   It's completely undocumented, but persp() returns a 4x4 matrix that
> allows you to compute the perspective transform for any other set of
> (x,y,z) coordinates, and you can then use points/lines/whatever to
> superimpose these on the plot.  (Hidden line removal is a harder problem.)
>
>   Here is some example code (the key is the trans3d() function) that adds
> a point, a line, and a plane (of sorts) to one of the persp() example
> plots.
>
>   This is essentially how the points3d(), plane3d(), xyz.convert()
> functions returned by a call to the scatterplot3d() function (found in its
> own library on CRAN) work ...
>
>
> ## matrix multiply c(3dv,1) by transformation matrix:
> ## plot v[0]/v[3], v[1]/v[3]
>
> x <- seq(-10, 10, length = 50)
> 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
> par(bg = "white")
>
> trans3d <- function(x,y,z,pmat) {
>   tmat <- t((cbind(x,y,z,1)%*% pmat))
>   list(x=tmat[1,]/tmat[4,],y=tmat[2,]/tmat[4,])
> }
>
> pmat <- persp(x, y, z, theta = 30, phi = 30, expand = 0.5,
>               col = "lightblue", xlab = "X", ylab = "Y", zlab = "Z",
>               ticktype="detailed")
> m <- 1e-5
> points(trans3d(m,m,f(m,m),pmat),pch=16)
> z2 <- sapply(1:length(x),function(n)f(x[n],y[n]))
> lines(trans3d(x,y,z2,pmat),col="red",lwd=2)
> lines(trans3d(c(-10,10,10,-10,-10),
>               c(-10,-10,10,10,-10),
>               c(2,2,8,8,2),pmat),col="blue")
>
> On Thu, 14 Mar 2002, Nina Lieske wrote:
>
> > Dear R-users,
> >
> > is it possible to add a second plane to the persp()-plot? I
> couldn't find
> > any hint on that in the news archive...
> > I'm using R.1.4.1 for Windows.
> >
> > Thanks in advance,
> > Nina
> >
> >
> >
> -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.
> -.-.-.-.-.-.-
> > r-help mailing list -- Read
> http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
> > Send "info", "help", or "[un]subscribe"
> > (in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
> >
> _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
> _._._._._._._
> >
>
> --
> 318 Carr Hall                                bolker at zoo.ufl.edu
> Zoology Department, University of Florida    http://www.zoo.ufl.edu/bolker
> Box 118525                                   (ph)  352-392-5697
> Gainesville, FL 32611-8525                   (fax) 352-392-3704
>
> -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.
> -.-.-.-.-.-.-
> r-help mailing list -- Read
http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
_._


-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._