LXX
读完需要
速读仅需 1 分钟
前言：
高斯滤波利用了高斯函数经傅里叶变换后仍具有高斯函数的特性。令指定区域的权重为高斯分布，从而将高频的噪声点滤除。具体是将某一数据点与其前后各 n 个数据点加权平均，那些远大于操作距离的点被处理成固定的端点，这有助于识别间隙和端点。由于高斯滤波平均效果较小，在滤波的同时，能较好地保持数据原貌，因而常被使用。本文分享对 icesat-2 数据进行高斯滤波的尝试用法。
希望各位同学点个关注，点个小赞，这将是更新的动力，不胜感激❥(^_-)
另外上期中奖的三位同学为：
后续更多活动请大家多多关注！！！
1
对于点云中的每个点P，高斯滤波器会根据预定义的窗口大小，在其邻域内选择一定数量的邻居点，然后根据邻居点与目标点之间的距离的近远关系，为每个邻居点分配一个权重值。这些权重值可以根据高斯分布函数计算得到，通常与邻居点的距离成反比。最后，对邻居点的数值进行加权平均，得到目标点P的新数值。
一般三维点云较多的应用高斯滤波，本文以此为例：点云高斯滤波是一种在三维点云数据上应用高斯滤波器的技术。在点云高斯滤波中，每个点的邻域被加权平均，以产生平滑的点云，权重由一个高斯分布函数计算。在三维空间中函数形式为
式中，d为点云中每个点到邻域点的距离；σ为高斯函数的标准差。
对于每个点，高斯滤波器将其邻域内的所有点乘以一个对应的高斯权重，然后将其加权平均。过程公式为：
式中， p0为当前点的坐标值；pi为邻域内第i个点的坐标值；N为邻域内的点数；G(pi-p0)为第i个点的高斯权重，它与第i个点与当前点的距离pi-p0有关；F（x,y,z）为滤波后点的坐标值。
通过高斯滤波，可以有效地消除点云中的噪声，并平滑点云表面，使得点云数据更具可读性和准确性。接下来，将介绍如何在代码中实现高斯滤波。
PCL 中点云高斯滤波的计算原理（知乎：点云侠）
https://zhuanlan.zhihu.com/p/694782029
2
为了更快运算构建了KD树
import numpy as npfrom scipy.spatial import KDTreefrom scipy.signal import savgol_filterfrom scipy.optimize import curve_fitimport matplotlib.pyplot as plt# 模拟光子点云数据np.random.seed(42)N_total = 1000 # 总光子点数y_values = np.random.normal(loc=0.5, scale=0.1, size=N_total)y_values = np.clip(y_values, 0, 1) # 将 y 值裁剪在合理范围内 (0 到 1)# 生成对应的 x 值x_values = np.random.rand(N_total)# 合并成光子点云数据photon_data = np.column_stack((x_values, y_values))# 参数设置P = 20 # 有效邻域点数radius = np.sqrt(P / (N_total * np.pi)) # 半径计算公式# 构建 KD 树tree = KDTree(photon_data)# 计算每个光子的邻域点数并生成直方图neighbor_counts = []for point in photon_data:neighbors = tree.query_ball_point(point, radius)neighbor_counts.append(len(neighbors))# 直方图统计counts, bin_edges = np.histogram(neighbor_counts, bins=30)# 使用 Savitzky-Golay 滤波平滑直方图smooth_counts = savgol_filter(counts, window_length=5, polyorder=2)# 定义高斯函数，用于拟合双峰分布def gaussian(x, amp, mean, stddev):return amp * np.exp(-((x - mean) ** 2) / (2 * stddev ** 2))def double_gaussian(x, amp1, mean1, stddev1, amp2, mean2, stddev2):return gaussian(x, amp1, mean1, stddev1) + gaussian(x, amp2, mean2, stddev2)# 对平滑后的直方图进行双高斯拟合bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2popt, _ = curve_fit(double_gaussian, bin_centers, smooth_counts, p0=[1, 5, 1, 1, 15, 1])# 计算阈值mean1, stddev1 = popt[1], popt[2]mean2, stddev2 = popt[4], popt[5]threshold = (mean1 + mean2) / 2# 将光子点分类为信号或噪声signal_points = []noise_points = []for i, count in enumerate(neighbor_counts):if count >= threshold:signal_points.append(photon_data[i])else:noise_points.append(photon_data[i])signal_points = np.array(signal_points)noise_points = np.array(noise_points)# 可视化结果plt.figure(figsize=(12, 6))# 原始邻域计数直方图plt.subplot(1, 2, 1)plt.hist(neighbor_counts, bins=30, alpha=0.6, color='red', label='原始直方图')plt.plot(bin_centers, smooth_counts, label='平滑直方图', color='blue')plt.plot(bin_centers, double_gaussian(bin_centers, *popt), label='双高斯拟合', color='red')plt.axvline(threshold, color='green', linestyle='--', label='阈值')plt.legend()plt.rcParams['font.sans-serif'] = ['STSong']plt.rcParams['axes.unicode_minus'] = Falseplt.title("邻近点数量直方图")plt.xlabel("邻近点数量")plt.ylabel("频率")# 光子点分类可视化plt.subplot(1, 2, 2)plt.scatter(noise_points[:, 0], noise_points[:, 1], color='red', label='Noise', alpha=0.5)plt.scatter(signal_points[:, 0], signal_points[:, 1], color='blue', label='Signal', alpha=0.5)plt.legend()plt.rcParams['font.sans-serif'] = ['STSong']plt.rcParams['axes.unicode_minus'] = Falseplt.title("光子点分类")plt.xlabel("X ")plt.ylabel("Y")plt.tight_layout()plt.show()
结果如下:
【参考】
PCL 中点云高斯滤波的计算原理
https://zhuanlan.zhihu.com/p/694782029
