如何高效对3D数组每行执行异尺寸数组运算以加速高光谱校正？

问题背景 #

我有一个spatRaster格式的大型3D高光谱立方体，还有两个1行像素的数组（由spatRaster转换而来）：暗参考图像dark_current和白参考图像white_reference，以及几个校正常数。需要对样本图像的每一行执行如下运算：

校正值 = R_ref × int_time × ((样本行像素值 - 暗参考行像素值) / 白参考行像素值)

目前用for循环实现但速度很慢，尝试apply未成功，求最快的加速方案。

原代码如下：

#hyperspectral cube
nlines_s <- 1412  #number of lines of pixels in the image
npixels_s <- 1024 #number of pixels within each line
nbands_s <- 448   #number of spectral bands within each pixel

sample_image <- terra::rast(array(runif(1),dim = c(nlines_s ,npixels_s ,nbands_s )))

# spatRasters already transformed to array type
dark_current <- array(runif(1),dim = c(1024,448)) 
white_reference <- array(runif(1),dim = c(1024,448))

#integration time constants
t_s <- 13
t_w <- 13

int_time <- t_w/t_s

#reflectance of reference material
R_ref <- 0.99 

#Process that needs to be speed up
reflectance_image <- array(0, c(nlines_s,npixels_s,nbands_s ))

for (fr in 1:nlines_s) {
      
      reflectance_image[fr, , ] <- as.matrix(R_ref * int_time  * ((sample_image[fr, ,]-dark_current)/white_reference) )

print(paste0("Calculating reflectance: ",round((fr/nlines_s)*100,2), "%"))      
}

ri_r <- terra::rast(reflectance_image)
names(ri_r) <- names(sample_image)

最优加速方案 #

方案1：直接使用terra栅格的向量化运算 #

terra包底层基于C++实现，支持栅格的向量化元素级运算，无需手动循环，效率远高于R原生循环。核心是将暗/白参考转换为与样本图像维度匹配的栅格，利用自动广播完成行维度的匹配。

# 将暗参考数组转换为spatRaster，并复制到与样本图像相同的行数
dark_rast <- terra::rast(dark_current, nrow = 1, ncol = npixels_s, nlyr = nbands_s)
dark_rast <- terra::rep(dark_rast, nlines_s)

# 同理处理白参考
white_rast <- terra::rast(white_reference, nrow = 1, ncol = npixels_s, nlyr = nbands_s)
white_rast <- terra::rep(white_rast, nlines_s)

# 直接执行向量化校正运算，terra会自动优化计算（支持多线程）
ri_r <- R_ref * int_time * ((sample_image - dark_rast) / white_rast)

# 保留原光谱波段名称
names(ri_r) <- names(sample_image)

优势：

• 避免数组与栅格的频繁转换，减少内存开销
• 底层C++运算比R循环快数倍至数十倍
• 默认开启多线程，可通过terra::terraOptions(nthreads = 你的线程数)手动调整

方案2：数组级广播运算 #

若需基于数组操作，利用R的自动广播特性（R 4.0+支持）扩展暗/白参考的行维度，直接执行向量化运算，跳过循环。

# 将暗/白参考扩展为与样本图像相同的三维数组（自动广播行维度）
dark_broadcast <- array(dark_current, dim = c(nlines_s, npixels_s, nbands_s))
white_broadcast <- array(white_reference, dim = c(nlines_s, npixels_s, nbands_s))

# 将spatRaster转为三维数组
sample_array <- terra::values(sample_image, mat = FALSE)

# 执行向量化校正运算
reflectance_array <- R_ref * int_time * ((sample_array - dark_broadcast) / white_broadcast)

# 转换回spatRaster
ri_r <- terra::rast(reflectance_array, nrow = nlines_s, ncol = npixels_s)
names(ri_r) <- names(sample_image)

注意事项：

• 提前检查white_reference中是否存在0值，避免除以0错误：white_reference[white_reference == 0] <- 1e-8
• 若内存不足，可按波段或行块分批次处理，但上述方案已是内存允许下的最优选择

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