在Lidr R中计算LAS文件的EigenVector12指标及提速问询

开发者阿泽

2026-7-2

问题：计算ALS点云的EigenVector12指标并加速处理

我需要计算TreeLS中FastPointMetrics函数的EigenVector12指标（定义为“3D特征向量系数，第二特征向量的第一载荷”），但该函数仅支持TLS数据，无法处理我的ALS点云。于是我用point_metrics()写了自定义函数，为每个点的20个最近邻计算该指标，再生成20m像素范围的均值栅格。

我的代码如下：

EV12_metrics = function(x,y,z, th1 = 20) {
  xyz <- cbind(x,y,z)         # create XYZ matrix
  cov_m <- cov(xyz)           # compute covariance of the point locations
  e = eigen(cov_m)            # get eigenvalues and eigenvectors from the matrix
  vect12 = e$vectors[4]       # Select the 1st load of the second eigenvector
  return(list(EV12 = vect12)) # return with a list of values 
}

M = point_metrics(laz, ~EV12_metrics(X,Y,Z), k = 20)   # apply the function to the point cloud 
laz = add_attribute(laz, M$EV12, "EV12")               # Add the attribute to the LAS points
gridmets = grid_metrics(laz, mean(EV12), res = 20)     # Create raster of the mean of the values at 20m resolution

代码能正常运行，但处理速度太慢。后来发现point_eigenvalues()函数，但它只返回特征值和其他我不需要的指标，想知道能不能用它实现需求，或者有没有办法加速自定义函数？

解答

一、先纠正代码中的索引歧义

你代码里e$vectors[4]的写法虽然结果正确（R中矩阵按列优先存储，3x3矩阵的第4个元素就是第1行第2列），但用行列索引e$vectors[1, 2]更清晰，明确指向第二特征向量的第一载荷，避免后续维护时混淆。

二、关于`point_eigenvalues()`的可用性

point_eigenvalues()是基于C++底层实现的高效函数，但它仅返回特征值（最大、中等、最小）和衍生的形状指标（如线性度、平面度），没有暴露特征向量的结果，因此无法直接用它获取EigenVector12指标。

三、加速自定义函数的方案

1. 利用对称矩阵的快速特征分解

协方差矩阵是对称矩阵，给eigen()函数加上symmetric = TRUE参数，会调用专门针对对称矩阵的高效分解算法，能大幅提升计算速度：

EV12_metrics = function(x,y,z) {
  xyz <- cbind(x,y,z)
  cov_m <- cov(xyz)
  # 针对对称矩阵启用快速特征分解
  e = eigen(cov_m, symmetric = TRUE)
  # 明确索引：第二特征向量的第一载荷
  vect12 = e$vectors[1, 2]
  return(list(EV12 = vect12))
}

另外，函数中未使用的th1参数可直接删除，减少冗余。

2. 启用并行计算

lidR的point_metrics()支持多核心并行处理，通过cores参数指定核心数（根据你的CPU核心数调整），能显著提升大数据量下的处理速度：

# 示例：使用4核心并行计算
M = point_metrics(laz, ~EV12_metrics(X,Y,Z), k = 20, cores = 4)

3. 跳过点属性添加，直接栅格化

如果不需要将EV12作为点云属性保存，可以直接在grid_metrics()中嵌套计算逻辑，减少一次数据拷贝的内存开销，进一步提升效率：

gridmets = grid_metrics(laz, {
  xyz <- cbind(X,Y,Z)
  cov_m <- cov(xyz)
  e = eigen(cov_m, symmetric = TRUE)
  # 直接计算栅格内的EV12均值
  mean(e$vectors[1, 2])
}, res = 20, k = 20, cores = 4)

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