【编程实践】点云曲率计算与可视化

基于Open3D和kNN的点云曲率计算与可视化，按下面指标定义的曲率

输入数据

• 输入数据
  点云文件ply或pcd，加载为点云对象

指标定义

• 此处曲率定义：采用基于点云PCA协方差分析的局部曲率估计方法。对于每个点，其曲率定义为：
  $$\kappa =\frac{{\lambda }_1}{{\lambda }_1+{\lambda }_2+{\lambda }_3}$$
  其中${\lambda }_1$是协方差矩阵最小特征值，反映法向方向上的变化程度。该指标不具备严格几何意义的主曲率属性，但在点云中广泛应用于反映局部几何尖锐程度和表面变化趋势。

数学处理过程

1. 邻域搜索（K-Nearest Neighbour）

• 对每个点$p_i$：
• • 使用KD树查找其最近的k个临近点，点本身也被包含在邻域中；
• • 得到临近点集N(pi)={pi1,pi2,...,pik}\mathscr{N}(p_i)=\{p_{i1},p_{i2},...,p_{ik}\}N(pi​)={pi1​,pi2​,...,pik​} ，对于某个点$p_i$有k个邻接点。

2. 协方差矩阵计算
   对于每个邻域点集$\mathcal{N}(p_i)$，计算其三维坐标协方差矩阵：
   $$C_i=Cov(\mathcal{N}(p_i))\in {\mathbb{R}}^{3\times 3}$$
   即：
   $$C_i=\frac{1}{k}\sum _{j=1}^k(p_{ij}-{\overline{p}}_i)(p_{ij}-{\overline{p}}_i{)}^T$$
   其中，${\overline{p}}_i$是邻域点的均值向量。

3. 特征值分级（Eigen Decomposition）
   计算协方差矩阵的特征值${\lambda }_1$、${\lambda }_2$、${\lambda }_3$并排序：
   $${\lambda }_1\le {\lambda }_2\le {\lambda }_3$$
   特征值代表该邻域点云在主轴方向上的分布方差。
   注：此处写的是数学上的标准形式1/k，属于有偏估计；实际适用了numpy的无偏估计1/(k-1).

4. 曲率估计
   使用下面公式估算每个点的局部曲率：
   $${\kappa }_i=\frac{{\lambda }_1}{\sum ({\lambda }_i)+ϵ}$$
   其中，${\lambda }_1$为最小特征值，表示法向方向上的变化；分母为总变化量，归一化；$ϵ$为防止分母为0加的小常数。

曲率可视化（颜色映射）

• 对所有点的曲率值${\kappa }_i$进行归一化：
  $${\kappa }_i^{\text{norm}}=\frac{{\kappa }_i-\min (\kappa )}{\max (\kappa )-\min (\kappa )+ϵ}$$
• 将归一化值使用Matplotlib的 jet 颜色映射为RGB颜色：
  $${\text{color}}_i=jet({\kappa }_i^{\text{norm}})$$
• 将颜色赋值回Open3D点云对象用于可视化。

import open3d as o3d
import numpy as np
import matplotlib.pyplot as plt# 加载点云
pcd = o3d.io.read_point_cloud("Result.ply")
# 构建KD树
pcd_tree = o3d.geometry.KDTreeFlann(pcd)
# ---- 曲率计算函数（基于kNN） ----
def compute_curvature(pcd, knn=30):points = np.asarray(pcd.points)curvatures = []for i in range(len(points)):_, idx, _ = pcd_tree.search_knn_vector_3d(points[i], knn)neighbors = points[idx]covariance = np.cov(neighbors.T)eigvals, _ = np.linalg.eigh(covariance)eigvals = np.sort(eigvals)curvature = eigvals[0] / (np.sum(eigvals) + 1e-10)  # 避免除0curvatures.append(curvature)return np.array(curvatures)
# ---- 计算曲率 ----
curvatures = compute_curvature(pcd, knn=20)
# ---- 映射为颜色（使用matplotlib colormap） ----
norm_curvatures = (curvatures - np.min(curvatures)) / (np.max(curvatures) - np.min(curvatures) + 1e-10)
colors = plt.cm.jet(norm_curvatures)[:, :3]  # 只取 RGB，不要 alpha
pcd.colors = o3d.utility.Vector3dVector(colors)
# ---- 可视化 ----
o3d.visualization.draw_geometries([pcd])# ---- 可选：保存带颜色的点云 ----
# o3d.io.write_point_cloud("curvature_colored.ply", pcd)

计算结果可视化：