Search code examples
rshadowraytracingrgl

Raytracing many objects and exporting shaded ground as a raster

I have models consisting of thousands of cylinders which are stored in a data frame with start-coordinates, end-coordinates, lengths and radii. I want to simulate the shading they would create in the real world, given a certain position of a light source. As a result, I would like to have a raster which contains the information whether the ground is shaded or not (on the xy-plane). Is there a way to do this in R? Even when handling several thousand cylinder objects at the same time?

Here a mockup for a single cylinder:

enter image description here

I usually draw my cylinders with the rgl package, but it would be also okay if I have to use another package. I figured I might be able to use the packages rayrender or raytracing, but I don't know to export the shaded ground from the view to an array or raster.

Edit: Code for creating & plotting some cylinders:

library(rgl)

# some cylinders
cylinder <- structure(list(radius = c(0.01, 0.01, 0.01, 0.02, 0.01, 0.01, 
0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 
0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 
0.01, 0.01), length = c(0.07, 0.13, 0.08, 0.08, 0.1, 0.08, 0.09, 
0.08, 0.07, 0.15, 0.02, 0.09, 0.12, 0.12, 0.08, 0.26, 0.1, 0.09, 
0.08, 0.02, 0.12, 0.11, 0.08, 0.06, 0.06, 0.19, 0.05, 0.1, 0.09, 
0.09), start_X = c(0.62, 0.61, 0.62, 0.63, 0.64, 0.65, 0.65, 
0.63, 0.63, 0.63, 0.63, 0.64, 0.63, 0.69, 0.79, 0.81, 0.92, 0.97, 
1.03, 1.07, 1.08, 1.15, 1.24, 1.3, 1.34, 0.61, 0.5, 0.47, 0.4, 
0.37), start_Y = c(0.13, 0.11, 0.09, 0.09, 0.09, 0.08, 0.09, 
0.08, 0.07, 0.08, 0.07, 0.07, 0.07, 0.05, 0.02, 0.04, 0.02, 0.01, 
0.04, 0.05, 0.05, 0.1, 0.15, 0.19, 0.22, 0.07, 0.13, 0.16, 0.17, 
0.26), start_Z = c(361, 361.07, 361.2, 361.29, 361.36, 361.46, 
361.54, 361.62, 361.7, 361.77, 361.9, 361.92, 361.78, 361.88, 
361.98, 362.04, 362.26, 362.35, 362.39, 362.46, 362.48, 362.56, 
362.62, 362.65, 362.66, 361.76, 361.91, 361.95, 362.01, 362.07
), axis_X = c(-0.09, 0.05, 0.12, 0.14, 0.1, -0.03, -0.15, -0.07, 
-0.07, -0.2, 0.52, -0.62, 0.43, 0.54, 0.16, 0.35, 0.53, 0.43, 
0.76, 0.58, 0.63, 0.74, 0.66, 0.56, 0.79, -0.61, -0.65, -0.64, 
-0.33, -0.7), axis_Y = c(-0.12, -0.09, -0.01, -0.08, -0.08, 0.01, 
-0.11, -0.14, -0.04, -0.06, 0.06, 0.59, -0.14, -0.14, 0.38, -0.22, 
0.1, 0, 0.14, 0.15, 0.47, 0.45, 0.46, 0.67, 0.48, 0.28, 0.2, 
0, 0.55, 0.16), axis_Z = c(0.99, 0.99, 0.99, 0.99, 0.99, 1, 0.98, 
0.99, 1, 0.98, 0.85, 0.53, 0.89, 0.83, 0.91, 0.91, 0.84, 0.9, 
0.64, 0.8, 0.61, 0.5, 0.59, 0.48, -0.39, 0.74, 0.74, 0.77, 0.77, 
0.69)), row.names = c(NA, 30L), class = "data.frame")

# calculate end points of cylinders
# (cylinders have starting coordinates, a length and a direction as unit vector)
cylinder$end_X = cylinder$start_X + cylinder$axis_X * cylinder$length
cylinder$end_Y = cylinder$start_Y + cylinder$axis_Y * cylinder$length
cylinder$end_Z = cylinder$start_Z + cylinder$axis_Z * cylinder$length
  
# prepare cylinders
cylinder_list <- lapply(1:nrow(cylinder), function(i) {
  cyl <- cylinder3d(
    center = cbind(
      c(cylinder$start_X[i], cylinder$end_X[i]),
      c(cylinder$start_Y[i], cylinder$end_Y[i]),
      c(cylinder$start_Z[i], cylinder$end_Z[i])),
    radius = cylinder$radius[i],
    closed = -2)
  cyl
})

# plot cylinders
open3d()
par3d(windowRect = c(50,50,650, 650))
shade3d(shapelist3d(cylinder_list, plot = FALSE), color = "blue")

I would like the light source to be a point source from infinite distance, as I would like to simulate sunlight. As the sun is so far away, I would just assume the light beams to be parallel to each other.

Solution

  • I believe this does what you want:

    library(rgl)

# some cylinders
# ... deleted ...

# cylinder$start_Z <- cylinder$start_Z -  361.5

# calculate end points of cylinders
# (cylinders have starting coordinates, a length and a direction as unit vector)
cylinder$end_X = cylinder$start_X + cylinder$axis_X * cylinder$length
cylinder$end_Y = cylinder$start_Y + cylinder$axis_Y * cylinder$length
cylinder$end_Z = cylinder$start_Z + cylinder$axis_Z * cylinder$length

# prepare cylinders
cylinder_list <- lapply(1:nrow(cylinder), function(i) {
  cyl <- cylinder3d(
    center = cbind(
      c(cylinder$start_X[i], cylinder$end_X[i]),
      c(cylinder$start_Y[i], cylinder$end_Y[i]),
      c(cylinder$start_Z[i], cylinder$end_Z[i])),
    radius = cylinder$radius[i],
    closed = -2)
  cyl
})

cylinder_list <- shapelist3d(cylinder_list, plot = FALSE)

# plot cylinders
open3d()
#> glX 
#>   1
par3d(windowRect = c(50,50,650, 650))
# shade3d(cylinder_list, color = "blue")

# Suppose Sun is in direction xyz

sun <- c(0,0,1) #  straight up

# Get a matrix that projects sun to c(0,0,0),
# and leaves anything with z = 0 alone

M <- rbind( cbind( diag(2), c(-sun[1]/sun[3], -sun[2]/sun[3])), 0)

# Replot the shadows of the cylinders

shadows <- transform3d(cylinder_list, t(M))
shade3d(shadows, col = "gray", lit = FALSE)

# decorate3d()   # Add axes

par3d(userMatrix = diag(4)) # Display flat.

lowlevel() # And show it in reprex


    # Use snapshot3d() to get PNG file containing the final image

    Created on 2022-07-12 by the reprex package (v2.0.1)

    It should include your definition of the cylinders dataframe, but I left that out to make it easier to read.

    If you want to show the cylinders as well as the shadows, uncomment the shade3d() call on cylinder_list. You'll probably also want to make the Z values smaller; I subtracted 361.5 to put them near zero. If you don't do this, the graph will be extremely long and thin, and you won't be able to see anything.

    The idea of the code is to project the sun vector to c(0,0,0), while leaving anything with Z==0 alone. This flattens the cylinders into objects that fit in the Z == 0 plane. Setting the userMatrix to the identity matrix then displays the shadow. Use the snapshot3d() function to make this into a raster image (in a PNG file).