A graphing calculator via RGL

Here is a function that I wrote to help me visualize equations of 2 variables (i.e. a 3D surface) that kept cropping up in my Stats 630 class. The core component is fairly simple: I create a 2D grid (a matrix), and for each point on the grid, evaluate the incoming equation or function. The X, Y, and Z coordinates are then shipped off to the RGL package, along with a couple of decorations and RGL does the rest. I’ve also included the empirical marginal distributions for the function on 2 “faces” of the enclosing box. I’ve included a short video of this the code in action at the bottom of this page (below the actual code).


# equationViewer  - a native R tool for visualizing equations in 3D
#  f    - a function such as:  function(x,y) sin(sqrt(x^2 + y^2))
#         or                   function(x,y) 3*x^1.5  + 2*y
#       - in the absence of an actual equation, you can just ship in a data frame or matrix.
#         If it has 2 columns, it will think of them as 'x' and 'y' coordinates, each 
#         representing a single point of data.   The code will then convert that to a 
#         2D kernel density estimate and plot out the results. 
#       - Altenatively, if your dataframe or matrix has more than two columns, it will
#         think of it as a matrix of 'z' values and plot out their heights.
#  xRange   - range over which to evaluate the function. Default is 0 to 1.
#  yRange   - range over which to evaluate the function. Default is 0 to 1.
#  gridSize - I find that just 20 or 30 points per axis is plenty.  Too many and it will eat
#             up your memory and available compute power
#  dropFunction - an extra bit.  You can supply something like function(x,y) { x < y } and the
#                 equationViewer will drop out all points on your grid where x is less than y.
#                 Thus it allows you remove slices from your picture.
#  append   - add another surface to an already existing surface
# Examples:
# equationViewer(f=function(x,y) { sin(x) + sin(y) }, xRange=c(-2,2), yRange=c(-2,2))
# equationViewer(f=data.frame(x=rnorm(1000), y=rt(1000, df=5)))
# equationViewer(f=volcano)  # see ?volcano for details

#  TO DO:  
#   - change xRange to xlim, yRange to ylim, 
#   - args are z, x, y  (account for function in different parts (how to do a 3D spiral)
equationViewer <- function(f, xRange=c(0,1), yRange=c(0,1), gridSize=21, 
                           dropFunction, append=FALSE, xlab=NULL, ylab=NULL) {

  if (!(require(rgl))) {
		cat("\n\n  Looks like you don't have the package 'rgl' installed\n",
		    "  on this machine.   I'll install it now... \n\n")
	xRangeTrue <- NULL
	yRangeTrue <- NULL
  # Cute default function to plot
  if (missing(f)) {   
		f <- function(x,y) sin(sqrt(x^2 + y^2))
		xRange = c(-pi,pi);  yRange=c(-pi,pi);
  if (is.function(f)) {
    x   <- seq(xRange[1], xRange[2], length.out=gridSize)     # set up a 20x20 grid of x's and y's
    y   <- seq(yRange[1], yRange[2], length.out=gridSize)
    g   <- expand.grid(x=x, y=y)
    g$z <- f(g$x, g$y)
    mat <- matrix(g$z, nrow=gridSize, ncol=gridSize, byrow=FALSE)   # put in z values in a matrix
    # example of cutting off the p.d.f. via an "equation"
    if (!missing(dropFunction)) mat[dropFunction(g$x, g$y)] <- NA   

  # Otherwise the user didn't supply a function, rather they supplied a set
  # of data consisting of X and Y coordinates that will be slotted into a kde2d
  # to give the appearance of a density.
  } else if (ncol(f) == 2) {
		if (missing(xlab)) xlab = colnames(data.frame(f))[1]
		if (missing(ylab)) ylab = colnames(data.frame(f))[2]

    dens <- kde2d(f[,1], f[,2], n=gridSize)
    x    <- dens$x
    y    <- dens$y
    mat  <- dens$z
		g    <- expand.grid(x=x, y=y)
    if (!missing(dropFunction)) mat[dropFunction(g$x, g$y)] <- NA   
    # The following few steps are for graphical purposes
    # Artificially inflate the density estimate to make it look nice.
    mat  <- 0.3 * diff(range(c(x,y))) * mat / max(mat, na.rm=TRUE)
    f    <- function(x,y) kde2d(x,y)  
    xRange=range(x);  yRange=range(y);
  # Last possibility is that they supplied an actual data matrix with the
  # cell values as the 'height' values
  } else {
		mat  <- as.matrix(f)
    x    <- 1:nrow(mat)
    y    <- 1:ncol(mat)
		g    <- expand.grid(x=x, y=y)
		if (!missing(xRange)) xRangeTrue <- xRange  # use for the axes in a minute
		if (!missing(yRange)) yRangeTrue <- yRange  # use for the axes in a minute		
		xRange=range(x);  yRange=range(y);
		f <- function(x,y) Input_Matrix    # formalism to keep up with requirements below
    if (!missing(dropFunction)) mat[dropFunction(g$x, g$y)] <- NA     
#	as.Table <- function(x) t(x[nrow(x):1,])
	mat      <- mat[,ncol(mat):1]
  heightScale = diff(range(mat, na.rm=TRUE))  
  colValues = colSums(mat, na.rm=TRUE);                              # marginal sums of mat
  colValues = (colValues - min(colValues)) / diff(range(colValues))  # values from 0..1
  if (!all(is.finite(colValues))) colValues <- rep(1,length(colValues))
  colValues = colValues * heightScale  + min(mat, na.rm=TRUE)   

  rowValues = rowSums(mat, na.rm=TRUE);                              # marginal sums of mat
  rowValues = (rowValues - min(rowValues)) / diff(range(rowValues))  # values from 0..1
  if (!all(is.finite(rowValues))) rowValues <- rep(1,length(rowValues))
  rowValues = rowValues * heightScale  + min(mat, na.rm=TRUE)   
  # A couple of steps to assign colors to the heights in the "z matrix"
  numColors   <- 20
  colorLUT    <- terrain.colors(numColors+1, alpha=0)   
  colorIndex  <- round(numColors*as.vector(mat)/diff(range(mat, na.rm=TRUE)))
  colorIndex  <- colorIndex - min(colorIndex, na.rm=TRUE) + 1
  colorValues <- colorLUT[colorIndex]

  if (!append)  open3d(userMatrix=rotationMatrix(pi/3, pi/3, pi/3, 0))
	# N.B. rgl thinks of the X and Z axes differently than me.
	# It takes a *lot* of horsepower to render partially transparent objects...
	if (length(mat) > 25000) {
		rgl.surface(x=x, z=y, y=mat, color=colorValues)
	# ... but partial alpha transparency is good if you can do it
	} else {
		rgl.surface(x=x, z=y, y=mat, color=colorValues, alpha=0.7, back="lines", specular="grey")          	
#		surface3d(x=x, y=y, z=mat, color=colorValues, alpha=0.7, back="lines", specular="grey")          	
  if (!append) {
    # N.B. rgl thinks of the X and Z axes differently than me
  	xMin   <- xRange[1] - 0.05 * diff(xRange) 
    xSeq   <- signif(seq(from=xRange[1], to=xRange[2], length.out=5), 2)
    yMin   <- yRange[1] - 0.05 * diff(yRange) 
    ySeq   <- signif(seq(from=yRange[1], to=yRange[2], length.out=5), 2)
		zRange <- range(mat, na.rm=TRUE)
    zMin   <- zRange[1] - 0.05 * diff(zRange) 
    zSeq   <- signif(seq(from=zRange[1], to=zRange[2], length.out=5), 2)
		# Some gymnastics are required if we want the x- and y- axes to be on a different
		# scale of numbers than each other or the z- axis.   That is to say, after an 
		# hour or two of messing around, I can't see a way in RGL to define the x-, y- and
		# z- axes to be on different scales without forcing the resulting object into a 
		# 3D cube (e.g. persp3D).
		if (is.null(xRangeTrue)) {
			axis3d('x', pos=c(xMin, zMin, yMin), nticks=5, labels=xSeq, at=xSeq, col='red', alpha=0.7)
		} else {
			xSeqTrue <- signif(seq(from=xRangeTrue[1], to=xRangeTrue[2], length.out=5), 2)
			axis3d('x', pos=c(xMin, zMin, yMin), nticks=5, labels=xSeqTrue, at=xSeq, col='red', alpha=0.7)	
		axis3d('y', pos=c(xMin, zMin, yMin), nticks=5, labels=zSeq, at=zSeq, col='red', alpha=0.7)
		if (is.null(yRangeTrue)) {
			axis3d('z', pos=c(xMin, zMin, yMin), nticks=5, labels=rev(ySeq), at=ySeq, col='red', alpha=0.7)
		} else {
			ySeqTrue <- signif(seq(from=yRangeTrue[1], to=yRangeTrue[2], length.out=5), 2)
			axis3d('z', pos=c(xMin, zMin, yMin), nticks=5, labels=rev(ySeqTrue), at=ySeq, col='red', alpha=0.7)		
    equationName = paste("f(x,y) =", paste(deparse(f)[-1], collapse=""), sep="")
    equationName = gsub("[{}]", "", equationName, perl=TRUE)
    title3d(main=equationName, family="serif", font=1, color="black", cex=0.8, 
            xlab=ifelse(is.null(xlab), 'X axis', xlab), 
						ylab='Z axis', 
						zlab=ifelse(is.null(ylab), 'Y axis', ylab)) 

    # Add on the marginal density for x (ie  f(x))
		segments3d(x=rep(x, each=2),    
		           y=unlist(lapply(rowValues, function(a) c(a,zMin))),  # interdigitate zeros trick
		           z=rep(yMin, 2*length(x)), col='red', alpha=0.3, lwd=0.5)
		lines3d(x=x, y=rowValues, z=rep(yMin, length(x)), col='red', alpha=0.3, lwd=2)
    text3d(x=median(x, na.rm=TRUE), y=median(rowValues, na.rm=TRUE), z=yMin,  
           alpha=0.8, adj=c(0.5,0.5), family="serif", font=4, cex=0.9, col='red', text="f(x)")
    # Add on the marginal density for y (ie  f(y))
		segments3d(z=rep(y, each=2), 
		           y=unlist(lapply(colValues, function(a) c(a,zMin))),  # interdigitate zeros trick
		           x=rep(xMin, 2*length(x)), col='red', alpha=0.3, lwd=0.5)
		lines3d(x=rep(xMin, length(y)), y=colValues, z=y, col='red', alpha=0.3, lwd=2)
    text3d(z=median(y, na.rm=TRUE), x=xMin, y=median(colValues, na.rm=TRUE),
           alpha=0.8, adj=c(0.5,0.5), family="serif", font=4, cex=0.9, col='red', text="f(y)")

    play3d(spin3d(axis=c(0,1,0), rpm=-9), duration=1)  # rpm=-5, duration=3
  invisible(list(x=x, y=y, z=mat))     # return the matrix, mat, without printing it

Leave a Reply

Your email address will not be published. Required fields are marked *