## 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)
}
###################################################################