raytracing volcano dataset in R

gistfile1.txt

# Code taken from article: http://www.tylermw.com/throwing-shade/
# Author: Tyler Morgan-Wall
# License: GPLv3

volcanoshadow = matrix(1, ncol = ncol(volcano), nrow = nrow(volcano))

volc = list(x=1:nrow(volcano), y=1:ncol(volcano), z=volcano)
sunangle = 45 / 180*pi
angle = -90 / 180 * pi
diffangle = 90 / 180 * pi
numberangles = 25
anglebreaks = seq(angle, diffangle, length.out = numberangles)
maxdistance = floor(sqrt(ncol(volcano)^2 + nrow(volcano)^2))

for (i in 1:nrow(volcano)) {
  for (j in 1:ncol(volcano)) {
    for (anglei in anglebreaks) {
      for (k in 1:maxdistance) {
        xcoord = (i + sin(sunangle)*k)
        ycoord = (j + cos(sunangle)*k)
        if(xcoord > nrow(volcano) || 
           ycoord > ncol(volcano) || 
           xcoord < 0 || ycoord < 0) {
          break
        } else {
          tanangheight = volcano[i, j] + tan(anglei) * k 
          if (all(c(volcano[ceiling(xcoord), ceiling(ycoord)],
                    volcano[floor(xcoord), ceiling(ycoord)],
                    volcano[ceiling(xcoord), floor(ycoord)],
                    volcano[floor(xcoord), floor(ycoord)]) < tanangheight)) next
          if (tanangheight < akima::bilinear(volc$x,volc$y,volc$z,x0=xcoord,y0=ycoord)$z) {
            volcanoshadow[i, j] =  volcanoshadow[i, j] - 1 / length(anglebreaks)
            break
          }
        }
      }
    }
  }
}