Attempt to vectorize rayshader

shader.R

faster_bilinear <- function (Z, x0, y0){
  i = floor(x0)
  j = floor(y0)
  XT = (x0 - i)
  YT = (y0 - j)
  result = (1 - YT) * (1 - XT) * Z[i, j]
  nx = nrow(Z)
  ny = ncol(Z)
  if(i + 1 <= nx){
    result = result + (1-YT) * XT * Z[i + 1, j]
  }
  if(j + 1 <= ny){
    result = result + YT * (1-XT) * Z[i, j + 1]
  }
  if(i + 1 <= nx && j + 1 <= ny){
    result = result + YT * XT * Z[i + 1, j + 1]
  }
  result
}

bench_rays <- function() {
  volcanoshadow = matrix(1, ncol = ncol(volcano), nrow = nrow(volcano))

  sunangle = 45 / 180*pi
  angle = -90 / 180 * pi
  diffangle = 90 / 180 * pi
  numberangles = 25
  anglebreaks = sapply(seq(angle, diffangle, length.out = numberangles), tan)
  maxdistance = floor(sqrt(ncol(volcano)^2 + nrow(volcano)^2))
  sinsun = sin(sunangle)
  cossun = cos(sunangle)

  for (i in 1:nrow(volcano)) {
    for (j in 1:ncol(volcano)) {
      vij = volcano[i, j]
      for (anglei in anglebreaks) {
        for (k in 1:maxdistance) {

          xcoord = (i + sinsun*k)
          ycoord = (j + cossun*k)

          if(xcoord > nrow(volcano) || 
             ycoord > ncol(volcano) || 
             xcoord < 0 || ycoord < 0) {
            break
          } else {
            tanangheight = vij + 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 < faster_bilinear(volcano, xcoord, ycoord)) {
              volcanoshadow[i, j] =  volcanoshadow[i, j] - 1 / length(anglebreaks)
              break
            }
          }
        }
      }
    }
  }
}

ff_bilin <- function(Z, x, y) {
  i <- as.integer(x)
  j <- as.integer(y)
  XT <- x - i
  YT <- y - j

  # Expand matrix to handle OOB cases; algorithm as written uses zero
  # (implicitly), and so do we, but might be better to just use last row and
  # column values instead.  We have not tested whether the additional allocation
  # for the matrix is slower than trying to explicitly compute each of the
  # special cases.

  Z <- rbind(cbind(Z, 0), 0)

  # Create symbols for re-used vectors

  XT_1 <- 1 - XT
  YT_1 <- 1 - YT
  i1 <- i + 1L
  j1 <- j + 1L

  # compute interpolation in a vectorized manner

  (YT_1) * (XT_1) * Z[cbind(i,  j)] +
  (YT_1) * (XT)   * Z[cbind(i1, j)] +
  (YT)   * (XT_1) * Z[cbind(i,  j1)] +
  (YT)   * (XT)   * Z[cbind(i1,  j1)]
}


faster_bilinear(volcano, 1.5, 2.4)
ff_bilin(volcano, 1.5, 2.4)


bench_rays2 <- function() {
  volcanoshadow = matrix(1, ncol = ncol(volcano), nrow = nrow(volcano))

  sunangle = 45 / 180*pi
  angle = -90 / 180 * pi
  diffangle = 90 / 180 * pi
  numberangles = 25
  anglebreaks = sapply(seq(angle, diffangle, length.out = numberangles), tan)
  maxdistance = floor(sqrt(ncol(volcano)^2 + nrow(volcano)^2))
  sinsun = sin(sunangle)
  cossun = cos(sunangle)

  # for each point in our grid, compute the points towards the light source that
  # are within the grid

  seq <- seq_len(maxdistance)
  xs <- seq_len(maxdistance) * sinsun
  js <-  * cossun

  # big allocation here

  coords <- 







  for (i in 1:nrow(volcano)) {
    for (j in 1:ncol(volcano)) {
      vij = volcano[i, j]
      for (anglei in anglebreaks) {
        for (k in 1:maxdistance) {

          xcoord = (i + sinsun*k)
          ycoord = (j + cossun*k)

          if(xcoord > nrow(volcano) || 
             ycoord > ncol(volcano) || 
             xcoord < 0 || ycoord < 0) {
            break
          } else {
            tanangheight = vij + 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 < faster_bilinear(volcano, xcoord, ycoord)) {
              volcanoshadow[i, j] =  
                volcanoshadow[i, j] - 1 / length(anglebreaks)
              break
            }
          }
        }
      }
    }
  }
}