Skip to content

Instantly share code, notes, and snippets.

@brodieG
Created October 19, 2018 16:56
Show Gist options
  • Save brodieG/7afcc22c0a218e316442281f96fdb9d8 to your computer and use it in GitHub Desktop.
Save brodieG/7afcc22c0a218e316442281f96fdb9d8 to your computer and use it in GitHub Desktop.
Attempt to vectorize rayshader
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
}
}
}
}
}
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment