Version 0.3.0: Big Honkin' Update

-Added `sphere_shade` function to apply spherical texture maps to elevation matrices.
-Added `create_texture` function to allow user to programmatically create texture maps by specifying colors.
-Added `ambient_shade` function to calculate ambient occlusion shadow maps.
-Added `add_shadow` function to add shadow to existing shadow matrices/texture arrays.
-Added `calculate_normal` function to pre-calculate normal vectors of elevation matrix for faster spherical shading.
-Added `detect_water` function to automatically detect bodies of water (of a specified size) on an elevation matrix.
-Added `write_png` function to save hillshaded maps to file.
-Added `plot_map` function to parse and plot resulting hillshaded maps.
-Added the ability to add a mask to raytracing and ambient occlusion so only certain areas are raytraced.

-Changed names of functions to be snake_case.

-Fixed bugs.

.Rbuildignore

@@ -1,3 +1,7 @@
 ^.*\.Rproj$
 ^\.Rproj\.user$
-^README.Rmd$
+^README.Rmd$
+^\.RDataTmp$
+^README_cache$
+^README_files$
+^tests$

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: rayshader
 Type: Package
 Title: Raytraces Elevation Matrices to Produce Global Illumination Relief Maps
-Version: 0.2.1
-Date: 2018-05-19
+Version: 0.3.0
+Date: 2018-06-28
 Author: Tyler Morgan-Wall
 Maintainer: Tyler Morgan-Wall <[email protected]>
 Description: Uses ray tracing to produce global illumination maps of elevation matrices.
@@ -11,6 +11,12 @@ Imports: suncalc,
          doParallel,
          foreach,
          Rcpp,
-         progress
+         progress,
+         raster,
+         scales,
+         jpeg,
+         png,
+         akima,
+         magrittr
 LinkingTo: Rcpp, progress
 RoxygenNote: 6.0.1

NAMESPACE

@@ -1,10 +1,21 @@
 # Generated by roxygen2: do not edit by hand
 
-export(lambshade)
-export(rayshade)
+export(add_shadow)
+export(add_water)
+export(ambient_shade)
+export(calculate_normal)
+export(create_texture)
+export(detect_water)
+export(lamb_shade)
+export(plot_map)
+export(ray_shade)
+export(sphere_shade)
+export(write_png)
 import(doParallel)
 import(foreach)
 import(parallel)
 import(progress)
 importFrom(Rcpp,evalCpp)
+importFrom(grDevices,col2rgb)
+importFrom(grDevices,rainbow)
 useDynLib(rayshader, .registration = TRUE)

R/RcppExports.R

@@ -1,15 +1,35 @@
 # Generated by using Rcpp::compileAttributes() -> do not edit by hand
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
+calculate_normal_cpp <- function(heightmap, progbar) {
+    .Call(`_rayshader_calculate_normal_cpp`, heightmap, progbar)
+}
+
+construct_matrix <- function(image_reference, number_rows, number_cols, index_x, index_y) {
+    .Call(`_rayshader_construct_matrix`, image_reference, number_rows, number_cols, index_x, index_y)
+}
+
+fill_find_groups <- function(max_z_matrix) {
+    .Call(`_rayshader_fill_find_groups`, max_z_matrix)
+}
+
+find_groups_cpp <- function(max_z_matrix) {
+    .Call(`_rayshader_find_groups_cpp`, max_z_matrix)
+}
+
+interpolate_color <- function(color_nw, color_ne, color_se, color_sw) {
+    .Call(`_rayshader_interpolate_color`, color_nw, color_ne, color_se, color_sw)
+}
+
 lambshade_cpp <- function(heightmap, rayvector) {
     .Call(`_rayshader_lambshade_cpp`, heightmap, rayvector)
 }
 
-rayshade_multicore <- function(sunangle, anglebreaks, heightmap, zscale, maxsearch, row) {
-    .Call(`_rayshader_rayshade_multicore`, sunangle, anglebreaks, heightmap, zscale, maxsearch, row)
+rayshade_multicore <- function(sunangle, anglebreaks, heightmap, zscale, maxsearch, row, cache_mask) {
+    .Call(`_rayshader_rayshade_multicore`, sunangle, anglebreaks, heightmap, zscale, maxsearch, row, cache_mask)
 }
 
-rayshade_cpp <- function(sunangle, anglebreaks, heightmap, zscale, maxsearch) {
-    .Call(`_rayshader_rayshade_cpp`, sunangle, anglebreaks, heightmap, zscale, maxsearch)
+rayshade_cpp <- function(sunangle, anglebreaks, heightmap, zscale, maxsearch, cache_mask, progbar) {
+    .Call(`_rayshader_rayshade_cpp`, sunangle, anglebreaks, heightmap, zscale, maxsearch, cache_mask, progbar)
 }
 

R/add_shadow.R

@@ -0,0 +1,31 @@
+#'@title Add Shadow
+#'
+#'@description Multiplies a texture array or shadow map by a shadow map. 
+#'
+#'@param hillshade A three-dimensional RGB array or 2D matrix of shadow intensities. 
+#'@param shadowmap A matrix that incidates the intensity of the shadow at that point. 0 is full darkness, 1 is full light.
+#'@param max_darken Default 0.7. The lower limit for how much the image will be darkened. 0 is completely black,
+#'1 means the shadow map will have no effect.
+#'@return Shaded texture map.
+#'@export
+#'@examples
+#'#TBD
+add_shadow = function(hillshade, shadowmap, max_darken = 0.7) {
+  if(class(shadowmap) == "array" && class(hillshade) == "matrix") {
+    tempstore = hillshade
+    hillshade = shadowmap
+    shadowmap = tempstore
+  }
+  shadowmap = t(shadowmap[,ncol(shadowmap):1])
+  if(class(hillshade) == "array") {
+    hillshade[,,1] = hillshade[,,1] * scales::rescale(shadowmap,c(max_darken,1))
+    hillshade[,,2] = hillshade[,,2] * scales::rescale(shadowmap,c(max_darken,1))
+    hillshade[,,3] = hillshade[,,3] * scales::rescale(shadowmap,c(max_darken,1))
+    hillshade
+  } else if (class(hillshade) == "matrix") {
+    if(any(hillshade > 1 | hillshade < 0)) {
+      stop("Error: Not a shadow matrix. Intensities must be between 0 and 1. Pass your elevation matrix to ray_shade/lamb_shade/ambient_shade/sphere_shade first.")
+    }
+    hillshade * t(scales::rescale(shadowmap[nrow(shadowmap):1,],c(max_darken,1)))
+  }
+}

R/add_water.R

@@ -0,0 +1,68 @@
+#'@title add_water
+#'
+#'@description Adds a layer of water to a map.
+#'
+#'@param hillshade A three-dimensional RGB array.
+#'@param watermap Matrix indicating whether water was detected at that point. 1 indicates water, 0 indicates no water.
+#'@param color Default `imhof1`. The water fill color. A hexcode or recognized color string. 
+#'Also includes built-in colors to match the palettes included in sphere_shade: 
+#'(`imhof1`,`imhof2`,`imhof3`,`imhof4`, `desert`, `bw`, and `unicorn`). 
+#'@importFrom grDevices col2rgb rainbow
+#'@export
+#'@examples
+#'library(magrittr)
+#'#Here we even out a portion of the volcano dataset to simulate water:
+#'island_volcano = volcano
+#'island_volcano[island_volcano < mean(island_volcano)] = mean(island_volcano)
+#'
+#'#Setting a minimum area avoids classifying small flat areas as water:
+#'island_volcano %>%
+#'  sphere_shade(texture="imhof3") %>%
+#'  add_water(detect_water(island_volcano, min_area = 400),color="imhof3") %>%
+#'  plot_map()
+add_water = function(hillshade, watermap, color="imhof1") {
+  flipud = function(x) {
+    x[nrow(x):1,]
+  }
+  watermap = flipud(t(watermap))
+  if (color == "imhof1") {
+    color = col2rgb("#e9f9ee")
+  } else if (color == "imhof2") {
+    color = col2rgb("#337c73")
+  } else if (color == "imhof3") {
+    color = col2rgb("#4e7982")
+  } else if (color == "imhof4") {
+    color = col2rgb("#638d99")
+  } else if (color == "desert") {
+    color = col2rgb("#b2f4ff")
+  } else if (color == "bw") {
+    color = col2rgb("#ffffff")
+  } 
+  if(class(hillshade) != "array") {
+    stop("Argument `hillshade` must be a 3-layer array.")
+  }
+  if(missing(watermap)) {
+    stop("User must provide matrix indicating locations of bodies of water")
+  }
+  if(all(dim(watermap) != dim(hillshade)[1:2])) {
+    stop("`hillshade` and `watermap` dimensions must be the same; hillshade is ",
+         paste0(dim(hillshade)[1:2],collapse="x"),", watermap is ", paste0(dim(watermap)[1:2],collapse="x"))
+  }
+  for(i in 1:3) {
+    tempmat = hillshade[,,i]
+    if(color[1] != "unicorn") {
+      tempmat[watermap >= 1] = color[i]/256
+    } else {
+      unicolors = col2rgb(rainbow(256))
+      for(row in 1:nrow(watermap)) {
+        for(col in 1:ncol(watermap)) {
+          if(watermap[row,col] >= 1) {
+            tempmat[row,col] = unicolors[i,col %% 256+1]/256
+          }
+        }
+      }
+    }
+    hillshade[,,i] = tempmat
+  }
+  hillshade
+}

R/ambient_shade.R

@@ -0,0 +1,45 @@
+#'@title Ambient Occlusion
+#'
+#'@description Calculates Ambient Occlusion Shadow Map
+#'
+#'@param heightmap  two-dimensional matrix, where each entry in the matrix is the elevation at that point. All points are assumed to be evenly spaced.
+#'@param anglebreaks The angle(s), in degrees, as measured from the horizon from which the light originates.
+#'@param sunbreaks Default 12. Number of rays to be sent out in a circle, evenly spaced, around the point being tested.
+#'@param maxsearch Default 20. The maximum distance that the system should propogate rays to check for .
+#'@param multicore Default FALSE. If TRUE, multiple cores will be used to compute the shadow matrix. By default, this uses all cores available, unless the user has
+#'set `options("cores")` in which the multicore option will only use that many cores.
+#'@param zscale Default 1. The ratio between the x and y spacing (which are assumed to be equal) and the z axis. 
+#'@param remove_edges Default `TRUE`. Slices off artifacts on the edge of the shadow matrix.
+#'@param cache_mask Default `NULL`. A matrix of 1 and 0s, indicating which points on which the raytracer will operate.
+#'@param shadow_cache Default `NULL`. The shadow matrix to be updated at the points defined by the argument `cache_mask`.
+#'@param progbar Default `TRUE`. If `FALSE`, turns off progress bar.
+#'@return Shaded texture map.
+#'@export
+#'@examples
+#'#Here we produce a ambient occlusion map of the `volcano` elevation map.
+#'amb = ambient_shade(heightmap = volcano, 
+#'    sunbreaks = 15, 
+#'    maxsearch = 100)
+ambient_shade = function(heightmap, anglebreaks = seq(0,45,15), sunbreaks = 12, 
+                        maxsearch=20, multicore=FALSE, zscale=1, remove_edges=TRUE, cache_mask = NULL, shadow_cache=NULL, progbar=TRUE) {
+  if(sunbreaks < 3) {
+    stop("sunbreaks needs to be at least 3")
+  }
+  if(remove_edges) {
+    shademat = matrix(0,nrow=nrow(heightmap)-2,ncol = ncol(heightmap)-2)
+  } else {
+    shademat = matrix(0,nrow=nrow(heightmap),ncol = ncol(heightmap))
+  }
+  for(angle in seq(0,360,length.out = sunbreaks+1)[-(sunbreaks+1)]) {
+    shademat = shademat + ray_shade(heightmap,anglebreaks=anglebreaks,sunangle = angle, maxsearch = maxsearch, zscale=zscale,
+                         multicore=multicore,  remove_edges = remove_edges,
+                         lambert=FALSE, cache_mask = cache_mask, progbar = progbar)
+  }
+  shademat = shademat/sunbreaks
+  if(!is.null(shadow_cache)) {
+    cache_mask = (cache_mask)
+    shadow_cache[cache_mask == 1] = shademat[cache_mask == 1]
+    shademat = matrix(shadow_cache,nrow=nrow(shademat),ncol=ncol(shademat))
+  }
+  return(shademat)
+}

R/calculate_normal.R

@@ -0,0 +1,29 @@
+#'@title Calculate Normal
+#'
+#'@description Calculates the normal unit vector for every point on the grid.
+#'
+#'@param heightmap A two-dimensional matrix, where each entry in the matrix is the elevation at that point. All points are assumed to be evenly spaced.
+#'@param zscale Default 1.
+#'@param remove_edges Default TRUE. Slices off artifacts on the edge of the shadow matrix.
+#'@param progbar Default `TRUE`. If `FALSE`, turns off progress bar.
+#'@return Matrix of light intensities at each point.
+#'@export
+#'@examples
+#'#Here we produce a light intensity map of the `volcano` elevation map.
+calculate_normal = function(heightmap, zscale=1, remove_edges = FALSE, progbar=TRUE) {
+  heightmap = heightmap / zscale
+  heightmap = heightmap[1:nrow(heightmap),ncol(heightmap):1]
+  matrices = calculate_normal_cpp(heightmap=heightmap, progbar = progbar)
+  returnnormal = list()
+  if(remove_edges) {
+    returnnormal[["x"]] = t(matrices$x[c(-1,-nrow(heightmap)),c(-1,-ncol(heightmap))])
+    returnnormal[["y"]] = t(matrices$y[c(-1,-nrow(heightmap)),c(-1,-ncol(heightmap))])
+    returnnormal[["z"]] = t(matrices$z[c(-1,-nrow(heightmap)),c(-1,-ncol(heightmap))])
+    return(returnnormal)
+  } else {
+    returnnormal[["x"]] = t(matrices$x)
+    returnnormal[["y"]] = t(matrices$y)
+    returnnormal[["z"]] = t(matrices$z)
+    return(returnnormal)
+  }
+}

R/create_texture.R

@@ -0,0 +1,58 @@
+#'@title create_texture
+#'
+#'@description Creates a texture map based on 5 user-supplied colors.
+#'
+#'@param lightcolor The main highlight color. Corresponds to the top center of the texture map.
+#'@param shadowcolor The main shadow color. Corresponds to the bottom center of the texture map. This color represents slopes directed 
+#'directly opposite to the main highlight color.
+#'@param leftcolor The left fill color. Corresponds to the left center of the texture map. This color represents slopes directed 
+#'90 degrees to the left of the main highlight color.
+#'@param rightcolor The right fill color. Corresponds to the right center of the texture map. This color represents slopes directed 
+#'90 degrees to the right of the main highlight color.
+#'@param centercolor The center color. Corresponds to the center of the texture map. This color represents flat areas.
+#'@param cornercolors Default `NULL`. The colors at the corners, in this order: NW, NE, SW, SE. If this vector isn't present (or
+#'all corners are specified), the mid-points will just be interpolated from the main colors.
+#'@export
+#'@examples
+#'#Here we produce a light intensity map of the `volcano` elevation map.
+create_texture = function(lightcolor,shadowcolor,leftcolor,rightcolor,centercolor,cornercolors=NULL) {
+  lightrgb = col2rgb(lightcolor)
+  shadowrgb = col2rgb(shadowcolor)
+  leftrgb = col2rgb(leftcolor)
+  rightrgb = col2rgb(rightcolor)
+  centerrgb = col2rgb(centercolor)
+  if(is.null(cornercolors) || length(cornercolors) != 4) {
+    nw_corner = (lightrgb+leftrgb)/2
+    ne_corner = (lightrgb+rightrgb)/2
+    sw_corner = (shadowrgb+leftrgb)/2
+    se_corner = (shadowrgb+rightrgb)/2
+  } else {
+    cornercolorsrgb = lapply(cornercolors,col2rgb)
+    nw_corner = cornercolorsrgb[[1]]
+    ne_corner = cornercolorsrgb[[2]]
+    se_corner = cornercolorsrgb[[3]]
+    sw_corner = cornercolorsrgb[[4]]
+  }
+  colorarray = array(0,dim=c(3,3,3))
+  for(i in 1:3) {
+    #center
+    colorarray[2,2,i] = centerrgb[i]
+
+    #edges
+    colorarray[1,2,i] = lightrgb[i]
+    colorarray[2,1,i] = leftrgb[i]
+    colorarray[2,3,i] = rightrgb[i]
+    colorarray[3,2,i] = shadowrgb[i]
+
+    #corners
+    colorarray[1,1,i] = nw_corner[i]
+    colorarray[1,3,i] = ne_corner[i]
+    colorarray[3,1,i] = se_corner[i]
+    colorarray[3,3,i] = sw_corner[i]
+  }
+  returnarray = array(0,dim=c(512,512,3))
+  for(i in 1:3) {
+    returnarray[,,i] = floor(akima::bilinear.grid(x=c(0,0.5,1),y=c(0,0.5,1),z=colorarray[,,i],nx=512,ny=512)$z)/256
+  }
+  returnarray
+}