R/write_stl.R

In rayshader: Create Maps and Visualize Data in 2D and 3D

#' Write STL
#'
#' @param filename The filename
#'
#' @keywords internal
write_stl = function(filename, rotate = FALSE, maxwidth = 100, unit="mm") {
  if(rotate) {
    zrot = matrix(0,3,3)
    zrot[1,1] = 1
    zrot[2,2] = 0
    zrot[3,3] = 0
    zrot[2,3] = 1
    zrot[3,2] = -1
  } else {
    zrot = matrix(0,3,3)
    zrot[1,1] = 1
    zrot[2,2] = 1
    zrot[3,3] = 1
  }
  rot_z_90 = function(vec) {
    vec %*% zrot
  }
  
  con = file(filename, "wb")
  on.exit(close(con))
  #Write STL header
  padding = paste(rep(" ", 80), collapse = "")
  header_val = substr(paste("binary", sprintf("'%s' generated by rayshader\n", filename), padding), 
                  1, 80)
  writeChar(header_val, con, nchars = 80, useBytes = TRUE, 
            eos = NULL)
  writeBin(0L, con, size = 4, endian = "little")
  
  #Write data
  triangles = 0L
  id_val = rgl::tagged3d(c("surface_tris", "base"))
  if(length(id_val) == 0) {
    stop("rayshader: No 3D scene data found")
  }
  min_x = Inf
  min_y = Inf
  min_z = Inf
  max_x = -Inf
  max_y = -Inf
  max_z = -Inf
  
  #Get bounding box
  for(id in id_val) {
    tmp_vert = rgl::rgl.attrib(id, "vertices")
    min_x = min(c(min_x,tmp_vert[,1]),na.rm=TRUE)
    min_y = min(c(min_y,tmp_vert[,2]),na.rm=TRUE)
    min_z = min(c(min_z,tmp_vert[,3]),na.rm=TRUE)
    
    max_x = max(c(max_x,tmp_vert[,1]),na.rm=TRUE)
    max_y = max(c(max_y,tmp_vert[,2]),na.rm=TRUE)
    max_z = max(c(max_z,tmp_vert[,3]),na.rm=TRUE)
  }
  dim1width = abs(min_x-max_x)
  dim2width = abs(min_z-max_z)
  dim3width = abs(min_y-max_y)
  maxdim = max(dim1width,dim2width)
  multiplier = maxwidth/maxdim
  
  for(id in id_val) {
    vertices = rgl::rgl.attrib(id, "vertices")
    indices = rgl::rgl.attrib(id, "indices")
    if(nrow(indices) == 0) {
      indices = matrix(1:nrow(vertices), ncol = 3, nrow = nrow(vertices)/3, byrow = TRUE)
    } else {
      indices = matrix(indices, ncol = 3, nrow = nrow(indices)/3, byrow = TRUE)
    }
    triangles = triangles + nrow(indices)
    for (i in seq_len(nrow(indices))) {
      if(!rotate) {
        vec0 = vertices[indices[i,1],] * multiplier
        vec1 = vertices[indices[i,2],] * multiplier
        vec2 = vertices[indices[i,3],] * multiplier
        normal = unit_vector(cross(as.vector(vec2-vec0), as.vector(vec1-vec0)))
      } else {
        vec0 = rot_z_90(vertices[indices[i,1],] * multiplier) 
        vec1 = rot_z_90(vertices[indices[i,2],] * multiplier) 
        vec2 = rot_z_90(vertices[indices[i,3],] * multiplier) 
        #The inverse transpose of zrot here is the same as zrot, so we can reuse the same function
        normal = rot_z_90(unit_vector(cross(as.vector(vec2-vec0), as.vector(vec1-vec0))))
      }
      writeBin(c(normal, vec0, vec1, vec2), con, size = 4, endian = "little")
      writeBin(0L, con, size = 2, endian = "little")
    }
  }
  seek(con, 80)
  writeBin(as.integer(triangles), con, size = 4, endian = "little")
  if(!rotate) {
    temp1 = dim2width
    dim2width = dim3width
    dim3width = temp1
  }
  if(unit == "mm") {
    message(sprintf("Dimensions of model are: %1.1f mm x %1.1f mm x %1.1f mm",dim1width*multiplier,dim2width*multiplier,dim3width*multiplier))
  } else {
    message(sprintf("Dimensions of model are: %1.2f in x %1.2f in x %1.2f in",dim1width*0.0393*multiplier,dim2width*0.0393*multiplier,dim3width*0.0393*multiplier))
  }
}