R中实现数千圆柱体平行光地面阴影栅格生成方法

核心需求 #

我持有存储在data frame中的数千个圆柱体模型数据，每条数据记录了对应圆柱体的起始坐标、终止坐标、长度与半径参数。我需要模拟给定光源条件下这些圆柱体产生的真实阴影效果，最终输出xy平面上标注地面是否被阴影覆盖的栅格（raster）数据，现咨询R语言是否可实现该需求，且能够支撑数千个圆柱体对象的同步处理。

本次模拟的光源为模拟太阳光的无限远距离点光源，假设所有光线互相平行。

现有使用基础 #

• 日常使用rgl包完成圆柱体绘制，也可接受其他R包实现方案
• 曾考虑使用rayrender、raytracing包开展光线追踪计算，但暂未掌握如何将渲染视图中的地面阴影结果导出为数组或栅格格式
• 单圆柱体效果可参考示意图：

现有圆柱体创建与绘制示例代码 #

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")

内容的提问来源于stack exchange，提问作者Zoe