Tuesday, August 26, 2014

Pretty pictures with R

One of the nice things about using the R programming language, apart from the fact that it is free, is that it can produce beautiful graphics. As well as do all sorts of standard and non-standard plots and graphs, R can also generate animated GIF files.
For example, suppose we have some height: \[h= exp( - \alpha r ) cos ( \beta \ r + p)\] where r is the distance from the origin: \[r=\sqrt{ x^2 + y^2 }\] Also we have constants
\(\alpha\) which determines the decay
\(\beta\) determines the frequency
and \(p\) is the phase

We generate a matrix which we've called h. Each element corresponds to an x and a y. We can convert that h matrix to an image and then cycle through values of phase to animate it. We will then have an image such as the one above.
Here is the R code that generates the h matrix:
 # generate an image of a ripple  
 rippleMat <- function( numPx, decay, phase, freq){  
   x <- seq( -1, 1, 2/ (numPx-1))  
   y <- x   
     
   # We have a matrix r[i,j] = sqrt( x[i] * x[i] + y[j] * y[j])  
   r <- sqrt( outer( x*x, y*y, "+"))    
     
   h <- exp( - decay * r ) * cos( phase + freq * r )  
   return (h)  
 }  
   
(The full code used to generate that animated gif is show below)

Alternatively, we could also plot a spiral. The basic equation is: \[r = exp( \phi \alpha + p ) \] Where r is the distance from the origin, \(\phi\) is the angle to the x-axis, p is the phase and \(\alpha\) determines how tightly coilded it is.
We could of course plot such a path using a line chart in R. But it might be a bit more interesting to color it in. Again we could iterate through values of the phase p to generate an animated GIF file:

We could also build an RGB color triangle. Each corner is set to one of { red, green, blue }. That color then fades across the triangle until it is completely absent on the opposite edge.

Just in case you're interested, here's the full source code for the triangle:
# Author: Philip Kinlen, Aug 2014
library(grid) # required for grid.raster(.)

# numPix should be an integer,
# and gapX < (numPix / 2) 
# where gapX is the minimum number of horizontal pixels from the triangle to the edge of the image
RGBTriangle <- function(numPix = 100, gapX = 10, doSmoothing=FALSE) {
  # the verticle gap between the triangle and edge of image
  gapY            <- round(numPix * ( 0.5 - sqrt(3)/4) + gapX * sqrt(3)/2)

  xArr            <- 1:numPix
  yArr            <- numPix:1   # the function call grid.raster(..) below will put the element with 
                                # highest row number at the bottom. 
                                # For our array, the last shall be at the bottom.
  # The x's a constant in the matrix columns, y's are constant along the matrix rows.
  xMat            <- matrix(rep(xArr, numPix), numPix, numPix, F)  
  yMat            <- matrix(rep(yArr, numPix), numPix, numPix, T)  

  m1              <- sqrt(3)                        # slope
  c1              <- gapY - m1 * gapX               # intercept
  
  m2              <- -sqrt(3)                       # slope
  c2              <- gapY - m2 * (numPix - gapX)    # intercept
  
  height          <- numPix - 2 * gapY              # height of triangle in pixels

  red             <- matrix(mapply(getRed,   xMat, yMat, m1, c1, m2, c2, gapY, height), numPix, numPix, T)
  green           <- matrix(mapply(getGreen, xMat, yMat, m1, c1, m2, c2, gapY, height), numPix, numPix, T)
  blue            <- matrix(mapply(getBlue,  xMat, yMat, m1, c1, m2, c2, gapY, height), numPix, numPix, T)
  
  col             <- rgb(red, green, blue)
  dim(col)        <- dim(red)
  grid.raster(col, interpolate=doSmoothing) 
  # One possible source of minor confusion is that when we write a cartesian co-ordinate (x,y)
  # the x comes first.
  # but in a matrix, we specify the row first, which corresponds to y
} 

#################################################################
getRed <- function( x, y, m1, c1, m2, c2, gapY, height){
   if ( isInsideTriangle( x, y, m1, c1, m2, c2, gapY)){
     res <- (y - gapY) / height    # A rounding error may cause res to be outside the range [0,1]
     return (max(0, min(1, res)))  # If left untreated that could cause an error in grid.raster(.)
   } else
     return (0)     
}
#################################################################
getBlue <- function( x, y, m1, c1, m2, c2, gapY, height){
  if ( isInsideTriangle( x, y, m1, c1, m2, c2, gapY)){ 
    res <- sqrt(3) * ((y - c2) / m2 - x) / ( 2 * height) 
    return (max(0, min(1, res))) 
  } else
    return (0)     
}  
#################################################################
getGreen <- function( x, y, m1, c1, m2, c2, gapY, height){
  if ( isInsideTriangle( x, y, m1, c1, m2, c2, gapY)){
    res <- 0.5 * (m1 * x + c1 - y) / height
    return (max(0, min(1, res))) 
  } else
    return (0)     
}  
#################################################################
isInsideTriangle <- function(x,y, m1, c1, m2 ,c2, gapY){
  
  if (  ( y >  gapY) 
      & ( y < ( m1 * x + c1 ) ) 
      & ( y < ( m2 * x + c2 ) ) )
    return(T)
  else 
    return (F)
}

And here is the code used to generate the animated GIFs contains the ripples and the spiral
# Author: Philip Kinlen, Aug 2014
###################################################################
library( caTools)  # library required for write.gif(..)

###################################################################
animateRipples   <- function(){
  
  filePath       <- paste(getwd(), "/animatedRipples.gif", sep="")
  numImages      <- 25
  numSecForCycle <- 1
  
  numPx          <- getNumPixels()
  phases         <- -2 * pi * ( 0:(numImages-1)) / numImages
  
  # The following line does the main work generating the data. 
  # The call to sapply(..) does the work that would have been done by a for-loop in
  # a language such as a Java.
  zeds           <- sapply(phases, rippleMat, simplify="array")    

  delayCS        <- 100 * numSecForCycle / numImages  # in hundreths of a second
  
  write.gif( zeds, filePath, col= getCustomColors(), delay = delayCS )
  
}
###################################################################
animateSpiral     <- function(){
  # for MS Windows might need to set the separator argument to a back-slash
  filePath        <- paste(getwd(), "animatedSpiral.gif", sep="/")
  
  numImagesPerSec <- 20
  numSecForCycle  <- 5
  
  numImages       <- numImagesPerSec * numSecForCycle
  numPx           <- getNumPixels()
  
  phases          <- 2 * pi * ( 0:(numImages-1)) / (numImages * getSpiralBeta())

  # The following line does the main work generating the data. 
  # The call to sapply(..) does the work that would have been done by a for-loop in
  # a language such as a Java.
  zeds            <- sapply( phases, spiralMat, simplify="array")
  
  delayCS         <- 100 * numSecForCycle / numImages  # in hundreths of a second
  
  write.gif( zeds, filePath, col= getCustomColors(), delay = delayCS )
}
###################################################################
getSpiralBeta <- function(){
  return ( 0.2)
}
###################################################################
getNumPixels <- function(){
  return ( 200)
}
###################################################################
getCustomColors   <- function(){
  maxCol     <- 255
  
  rising     <- 0:maxCol
  falling    <- maxCol - rising
  
  myRed      <- c(falling, rising) * 0.4 + 0.2
  myGreen    <- c(falling, rising) * 0.4 + 0.2
  myBlue     <- c(falling, rising) * 0.8 + 0.2
  
  return (rgb(  myRed, myGreen, myBlue,  maxCol, maxColorValue=maxCol ))
}  
###################################################################
# generate a matrix representation of a ripple
rippleMat <- function( phase, numPx = getNumPixels(), decay = 2.5, freq = 12 * pi){
   x <- seq( -1,  1, 2/ (numPx-1))
   y <- x 
   
   # We have a matrix r[i,j] = sqrt( x[i] * x[i] + y[j] * y[j])
   r <- sqrt( outer( x*x, y*y, "+"))   
   
   h <- exp( - decay * r ) * cos( phase + freq * r )
   return (h)
}
###################################################################
spiralMat    <- function(phase, beta = getSpiralBeta() , numPx= getNumPixels(), generateImage= FALSE){

  alpha    <- 0.1   # the main spiral follows: r = exp( alpha * phi)
  
  x        <- seq(-1, 1, 2 / (numPx - 1))
  y        <- x
  matR     <- sqrt(outer(x*x, y*y, "+"))
  
  matX     <- matrix(rep(x, numPx), numPx, numPx, F)  
  matY     <- matrix(rep(y, numPx), numPx, numPx, T)  
  theta    <- ( atan2(matY, matX) + pi) + phase

  # Note that:
  #      r >= exp( (2 * pi * (rev-1) + theta) * alpha)  
  # and  r <= exp( (2 * pi *  rev    + theta) * alpha)
  # where rev is the number of revolutions 
  rev      <- ceiling((log(matR)/ alpha  - theta) /( 2 * pi))

  phi      <- rev * 2 * pi + theta
  z        <- sin(phi * beta)

  if ( generateImage)
    image( z,  col = getCustomColors())
  else 
    return (z)
}
###################################################################