当前位置：首页 > news >正文

MLS平滑滤波

news 2025/8/5 11:22:15

1.前言

最近在学习，因此查阅相关资料，该怎么表述感觉有些困难

2.代码

2.1代码1

使用全局坐标系

参考：python点云移动最小二乘法(Moving Least Squares)平滑_移动最小二乘法python-CSDN博客

def Moving_Least_Squares_Smoothing_v1_explain( points, radius, tau, weights_name='gaussian'):"""如果输入open3d点云，则points = np.asarray(inlier_cloud.points)功能：使用Moving Least Squares Smooth进行点云的平滑必要参数：points: 点云数据[[x1,y1,z1],[x2,y2,z2],...,[xn,yn,zn]]radius: 邻居搜索半径tau: 权重系数可选参数：weights_name: 权重函数，可选：gaussian:高斯权重函数。根据距离计算权重，权重随epanechnikov:Epanechnikov核权重函数。在一定范围内提供权重，超出这个范围权重急剧下降。cubic_spline:三次样条权重函数。使用三次样条函数计算权重，提供平滑的权重变化。uniform:均匀权重函数。为所有点提供相同的权重。inverse_distance:距离反比权重函数。根据距离的反比例来计算权重。wendland:Wendland权重函数。使用Wendland函数计算局部光滑的权重。shepard_weight:Shepard权重函数。使用距离的负幂计算权重，权重随距离增大而减小。"""x = points[:, 0]y = points[:, 1]z = points[:, 2]# 构建KD树以快速找到邻居tree = KDTree(np.vstack((x, y)).T)smoothed_points = np.zeros((len(x), 3))for i in range(len(x)):# 使用KD树找到指定半径内的所有邻居idx = tree.query_ball_point([x[i], y[i]], r=radius)"""这是 SciPy KDTree 中的一个关键方法，用于查找指定半径内的所有邻近点。其核心功能是：以点 [x[i], y[i]] 为中心，在 radius 半径范围内，查找KD树中所有满足条件的邻近点.返回值返回类型：list of int内容：包含所有邻域点在原始点云数组中的索引值示例：[3, 5, 8, 12] 表示原始点云中索引为3,5,8,12的点在当前点的邻域内"""# 获取邻居点x_neigh = x[idx]y_neigh = y[idx]z_neigh = z[idx]# 如果邻域内没有足够的点，则跳过if len(x_neigh) < 3:continue# 计算邻居点到当前点的距离distances = np.sqrt((x_neigh - x[i]) ** 2 + (y_neigh - y[i]) ** 2)# 计算权重weights = utils.choose_weight(weights_name, distances, tau)"""weights = np.exp(-0.5 * (distance / bandwidth) ** 2)高斯核函数输出一个介于 0 到 1 的值，表示两点之间的相似度（距离越近，值越接近1；距离越远，值越接近0）。参数解释distance：两点之间的欧氏距离（或其他距离度量）。bandwidth：控制高斯曲线的宽度（平滑参数）。bandwidth越大：曲线越宽，远距离的点仍有较高权重（更平滑）。bandwidth越小：曲线越窄，只有近距离的点才有显著权重（更敏感）。计算步骤将距离归一化：distance / bandwidth。平方后乘以 -0.5，使结果符合高斯分布的形式。通过 np.exp 计算指数，得到最终权重。"""# 二次曲面拟合: z = ax^2 + by^2 + cxy + dx + ey + f# 构建加权最小二乘问题A = np.vstack([x_neigh ** 2 * weights, y_neigh ** 2 * weights, x_neigh * y_neigh * weights, x_neigh * weights,y_neigh * weights, weights]).Tb = z_neigh * weights   #形状：(n_neighbors,)"""设计矩阵AA = np.vstack([x_neigh ** 2 * weights,  # 二次项 x²y_neigh ** 2 * weights,  # 二次项 y²x_neigh * y_neigh * weights,  # 交叉项 xyx_neigh * weights,       # 线性项 xy_neigh * weights,       # 线性项 yweights                  # 常数项]).T形状：(n_neighbors, 6)numpy.vstack(tup)功能：沿第一个轴（垂直方向）拼接多个数组输入：元组 tup 包含要堆叠的数组序列输出：堆叠后的新数组stacked = np.array([[x1²w1, x2²w2, ..., xm²wm],  # 行1[y1²w1, y2²w2, ..., ym²wm],  # 行2[x1y1w1, x2y2w2, ..., xmymwm],  # 行3[x1w1, x2w2, ..., xmwm],     # 行4[y1w1, y2w2, ..., ymwm],     # 行5[w1, w2, ..., wm]            # 行6])结果形状：(6, m)m 是邻近点数量转置操作：A = stacked.T转置后形状：(m, 6)pythonA = [[x1²w1, y1²w1, x1y1w1, x1w1, y1w1, w1],  # 点1[x2²w2, y2²w2, x2y2w2, x2w2, y2w2, w2],  # 点2...,[xm²wm, ym²wm, xmymwm, xmwm, ymwm, wm]  # 点m]"""# 解线性方程组得到拟合参数fit_params = np.linalg.lstsq(A, b, rcond=None)[0]     #返回六个参数"""A：加权设计矩阵 (m_neighbors, 6)b：加权观测 z 值 (m_neighbors,)返回：二次曲面参数 [a, b, c, d, e, f]用最小二乘法求解线性方程组 A⋅x≈b 的最优解，即找到参数 x 使误差平方和最小。np.linalg.lstsq	NumPy 提供的最小二乘解函数（least squares）。A	输入的设计矩阵（形状为 (m, n)），每行是一个样本，每列是一个特征。b	目标值向量（形状为 (m,) 或 (m, k)，对应回归问题中的观测值）。rcond=None	禁用秩的截断阈值（默认行为），避免旧版本 NumPy 的警告。[0]	lstsq 返回多个值（解、残差、秩、奇异值），这里只取最优参数 x数学意义最小化以下目标函数：见b本如果 A 是满秩的（即列线性无关），解是唯一的。如果 A 秩亏（列相关），返回最小范数解（Moore-Penrose 伪逆）"""# 计算当前点的平滑值smoothed_z = fit_params[0] * x[i] ** 2 + fit_params[1] * y[i] ** 2 + fit_params[2] * x[i] * y[i] + \fit_params[3] * x[i] + fit_params[4] * y[i] + fit_params[5]smoothed_points[i] = [x[i], y[i], smoothed_z]"""这段代码使用拟合的二次曲面参数，计算当前点的平滑高度值：$(x[i], y[i])$：当前点的原始坐标{fit_params} = [a, b, c, d, e, f]$：拟合参数参数含义参数索引	符号	几何意义	对z值的影响0	a	x²系数	控制x方向的曲率1	b	y²系数	控制y方向的曲率2	c	xy系数	控制扭曲/扭转3	d	x系数	控制x方向的坡度4	e	y系数	控制y方向的坡度5	f	常数项	基准高度偏移几何解释曲率修正：$ax^2 + by^2$ 补偿局部曲率扭曲校正：$cxy$ 消除非正交畸变坡度调整：$dx + ey$ 修正倾斜基准校准：$f$ 确保高度连续性平滑效果分析项	平滑作用	频率响应x²/y²	抑制高频噪声	低通滤波xy	消除斜向噪声	方向滤波x/y	补偿线性趋势	带通滤波常数项	保持基准高度	直流保持""""""这段代码是在 用二次多项式曲面（二次曲面）来平滑（拟合）一组三维点，并把平滑后的 z 值存回 smoothed_points[i]。fit_params[0] 到 fit_params[5] 是 二次曲面参数,这些参数是之前通过最小二乘法（如 np.linalg.lstsq）拟合得到的把原始 x[i]、y[i] 和平滑后的 z 组成一个新点，存回 smoothed_points。"""result =  smoothed_pointspcd_mls = o3d.geometry.PointCloud()pcd_mls.points = o3d.utility.Vector3dVector(result)return pcd_mls

2.2代码2

参考：python——移动最小二乘拟合曲线/曲面_python 移动最小二乘法曲线拟合-CSDN博客

#2025.8.3
#点云侠  MLSimport numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.spatial import KDTreedef cubic_spline_3d(r):"""三维权函数（紧支域[0,1]）"""r = np.clip(r, 0, 1)return np.where(r <= 0.5,2 / 3 - 4 * r ** 2 + 4 * r ** 3,4 / 3 - 4 * r + 4 * r ** 2 - 4 / 3 * r ** 3)def mls_surface_fit(points, grid_x, grid_y, s_max=0.5, degree=2):"""MLS曲面拟合核心函数:param points: 点云数据 (n,3) [x,y,z]:param grid_x, grid_y: 拟合网格坐标:param s_max: 支持域半径:param degree: 多项式阶数:return: 拟合曲面Z值"""# 构建KD树加速邻域搜索tree = KDTree(points[:, :2])# 初始化拟合曲面Z_fit = np.zeros_like(grid_x)grid_points = np.column_stack([grid_x.ravel(), grid_y.ravel()])# 确定基函数 (二次曲面: 6个基函数)basis_funcs = [lambda x, y: 1,lambda x, y: x,lambda x, y: y,lambda x, y: x * y,lambda x, y: x ** 2,lambda x, y: y ** 2]if degree == 1:  # 线性曲面basis_funcs = basis_funcs[:3]for idx, (x, y) in enumerate(grid_points):# 搜索邻域点dists, idxs = tree.query([x, y], k=50, distance_upper_bound=s_max)"""用于查找最近的k个点，同时设置最大搜索距离限制。核心功能：以点 [x, y] 为中心，查找最近的50个点，但不超过 s_max 距离.dists	array	邻居点到中心点的距离（长度≤k）idxs	array	邻居点在原始数组中的索引（长度≤k）"""## 创建有效点掩码（过滤无穷远点）valid = dists < float('inf')# 提取有效邻域点local_points = points[idxs[valid]]if len(local_points) < len(basis_funcs) + 1:  #至少大于四个， 曲面至少大于3个点， 曲线至少大于2个点continue  # 邻域点不足# 计算权重r = dists[valid] / s_maxweights = cubic_spline_3d(r)# 构造设计矩阵A和观测向量bA = np.zeros((len(local_points), len(basis_funcs)))  #矩阵大小 mx6b = local_points[:, 2]for i, func in enumerate(basis_funcs):A[:, i] = [func(pt[0] - x, pt[1] - y) for pt in local_points]#A[:i]"""A[:, i] 是 NumPy 数组的切片语法，用于取出矩阵（二维数组）A 的第 i 列符号	含义A	一个二维 NumPy 数组（矩阵），形状为 (行数, 列数)。:	表示“所有行”。i	表示“第 i 列”。""""""这段代码用于构建标准移动最小二乘（MLS）的设计矩阵和观测向量，核心是：初始化设计矩阵 A准备观测向量 b用局部坐标系中的基函数填充 AA = np.zeros((len(local_points), len(basis_funcs)))  # 矩阵大小 m x n    功能：初始化设计矩阵  维度：m = len(local_points)：邻域点数n = len(basis_funcs)：基函数个数示例：二次曲面：6个基函数 → 6列线性平面：3个基函数 → 3列b = local_points[:, 2]功能：提取观测值向量内容：所有邻域点的 z 坐标形状：(m,) 向量for i, func in enumerate(basis_funcs):  #（x,y)当前点， pt近邻点A[:, i] = [func(pt[0] - x, pt[1] - y) for pt in local_points]功能：填充设计矩阵的每一列关键操作：遍历每个基函数 func对每个邻域点：计算局部坐标：Δx = pt[0]-x, Δy = pt[1]-y计算基函数值：func(Δx, Δy)将结果存入当前列"""# 加权最小二乘求解W = np.diag(weights)  ## 将权重向量转换为对角矩阵AWA = A.T @ W @ A"""矩阵	维度	说明A	(n, p)	    n 个观测，p 个特征W	(n, n)	    对角权重矩阵，对角线元素为 w₁, w₂, …, wₙAᵀ	(p, n)	    A 的转置Aᵀ @ W @ A	(p, p)	最终结果""""""维度变化：$\mathbf{A}^T$: n×m$\mathbf{W}$: m×m$\mathbf{A}$: m×n结果: n×n 矩阵物理意义：加权信息矩阵（Fisher信息矩阵）"""AWb = A.T @ W @ b"""维度变化：$\mathbf{A}^T$: n×m$\mathbf{W}$: m×m$\mathbf{b}$: m×1结果: n×1 向量物理意义：加权观测向量"""try:theta = np.linalg.solve(AWA, AWb)except np.linalg.LinAlgError:theta = np.linalg.lstsq(AWA, AWb, rcond=None)[0]"""方法	              原理	         适用场景	        优点     	缺点np.linalg.solve	直接求解（LU分解）	矩阵满秩且条件数好	速度快、精度高	对病态矩阵不稳定np.linalg.lstsq	最小二乘（SVD分解）	秩亏或病态矩阵	数值稳定、提供最小范数解	计算量大当正规方程矩阵病态时：解空间为超平面而非单点lstsq 返回最短的解向量对应最平滑的曲面"""# 计算当前点拟合高度z_fit = sum(theta[i] * func(0, 0) for i, func in enumerate(basis_funcs))"""theta[i]：第 i 个基函数的拟合系数func(0, 0)：在局部原点 (0,0) 计算基函数值sum：将所有基函数的贡献相加，  ：由于（x,y)=(0,0),因此 6个系数---仅常数项不为0"""Z_fit.flat[idx] = z_fit"""Z_fit：网格高度矩阵（二维数组）.flat[idx]：通过展平索引访问元素idx：当前网格点在展平数组中的索引Z_fit：二维网格矩阵（形状同 grid_x/grid_y）.flat：展平视图（一维索引访问）idx：当前网格点索引等效操作：Z_fit[row, col] = z_fit"""return Z_fit
#数学证明
# 定理：局部坐标系下 MLS 拟合曲面在原点处高度等于常数项系数"""
方法1：使用全局坐标系， 6个系数全部使用
方法2：点云侠： 使用局部坐标系，  6个系数仅使用常数项局部坐标系的优势
简化计算：只需常数项系数
几何直观：$\theta_0$ 直接表示中心点高度
数值稳定：避免高阶项贡献
物理一致：与曲面定义原点一致"""#
# # 示例：拟合鞍形曲面
# np.random.seed(42)
#
# # 生成带噪声的鞍形曲面数据
# x = np.random.uniform(-2, 2, 500)
# y = np.random.uniform(-2, 2, 500)
# z = x ** 2 - y ** 2 + np.random.normal(0, 0.1, 500)
# points = np.column_stack([x, y, z])
#
# # 创建拟合网格
# grid_x, grid_y = np.meshgrid(np.linspace(-2, 2, 50),
#                              np.linspace(-2, 2, 50))
#
# # MLS曲面拟合
# Z_fit = mls_surface_fit(points, grid_x, grid_y, s_max=1.0, degree=2)
#
# # 可视化
# # 设置中文显示
# plt.rcParams['font.sans-serif'] = ['LiSu']  # 显示中文
# plt.rcParams['axes.unicode_minus'] = False  # 正常显示负号
# fig = plt.figure(figsize=(12, 6))
# ax1 = fig.add_subplot(121, projection='3d')
# ax1.scatter(x, y, z, c=z, cmap='viridis', s=5)
# ax1.set_title('原始点云数据')
#
# ax2 = fig.add_subplot(122, projection='3d')
# ax2.plot_surface(grid_x, grid_y, Z_fit, cmap='plasma', alpha=0.8)
# ax2.set_title('MLS曲面拟合结果')
# plt.tight_layout()
# plt.show()