[Open3d系列]--点云曲线拟合

Open3d: 点云曲线拟合

因为项目需要分析点云数据，此文总结其中拟合点云的部分。

曲线拟合

首先定一个曲线方程:

def func(x, a, b, c):
    return a * x**2 + b * x + c

然后将点云数据结构转换为numpy数组：

points = np.asarray(pcd.points)

读取点数组中，x轴、y轴的数组：

xy_points = points[:, :2]
x = xy_points[:,x_axis_idx]
y = xy_points[:,y_axis_idx]

调用scipy.optimize中的curve_fit进行点拟合, 得到各项系数：

popt, pcov = curve_fit(func,x, y)

以上结果系数在返回值popt中，为一个元组

为了评估其拟合的好坏，此处计算了残差, 然后分析其分布情况:

residuals = y - func(x, *popt)

# 通过标准差， 判断其离散程度
np.sqrt(np.sum(residuals**2)/x.shape[0]) < std_threshold

拟合评估

R方，决定系数，又称拟合优度，通常用来描述数据对模型拟合程度的好坏，表示自变量X对因变量Y的解释程度。R方的取值在[0,1]之间，越接近1，说明回归拟合效果越好。比如R方=0.5，那么说明自变量可以解释因变量50%的变化原因。

R方计算公式：

R^{2} = 1 - \frac{\sum_{i}^{n} (y_{i} - \hat{y_{i}})^{2}}{\sum_{i}^{n} (y_{i} - \bar{y_{i}})^{2}} = 1 - S S E / S S T = S S R / S S T

其中 $\hat{y_{i}}$ 为残差， $\bar{y_{i}}$ 为平均值。

$S S E = \sum_{i}^{n} (y_{i} - \hat{y_{i}})^{2}$ : 残差平方和，SSE越大，模型不能解释的部分越多，则拟合效果越差
$S S R = \sum_{i}^{n} (\hat{y_{i}} - \bar{y_{i}})^{2}$ : 回归平方和，是反映自变量与因变量之间的相关程度的偏差平方和
$S S T = \sum_{i}^{n} (y_{i} - \bar{y_{i}})^{2}$ : 总体平方和，它是观测值与平均值之间差值的平方和。它的大小描述了数据集中的数的分散程度

在多元线性回归中，R方有个致命的问题，那就是随着自变量X的个数增加，R方会越来越大，R方越来越大就会认为模型拟合越来越好，但是实际上可能是由于自变量个数的增加导致的R方增大。

调整后的R方来解决这个问题,将自变量个数考虑进公式中，就得到了调整后R方:

R^{2} = 1 - (1 - R^{2}) \frac{n - 1}{n - k - 1}

其中n为样本数， k为自变量个数。

可视化

import matplotlib.pyplot as plt
import matplotlib
# 保证中文显示正常
matplotlib.rcParams['axes.unicode_minus'] = False
matplotlib.rcParams['font.family'] = 'SimHei'

x_fit = np.linspace(np.min(x), np.max(x), x.shape[0])
y_fit = func(x_fit, a_fit, b_fit, c_fit)
# 创建两个子图
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 6))
ax1.scatter(x,y, color='blue', label='Point Cloud')
ax1.plot(x_fit,y_fit, color='red', linewidth=2, label='Fitted Curve')
ax1.legend()
ax1.set_xlabel('X')
ax1.set_ylabel('Y')
ax1.set_title(f"曲线拟合(Fitted Curve)\n$y={a_fit:.3f}*x^2+{b_fit:.3f}*x+{c_fit:.3f}$")

# 绘制残差图
ax2.scatter(x, residuals)
ax2.axhline(y=0, color='r', linestyle='--')
ax2.set_xlabel('X')
ax2.set_ylabel('Residuals')
ax2.set_title('残差(Residuals)')

# 计算R-squared
mean_y = np.mean(y)
ss_total = np.sum((y - mean_y) ** 2)
ss_residual = np.sum(residuals ** 2)
r_squared = 1 - (ss_residual / ss_total)

# 值越大拟合效果越好
ax2.text(0.02, 0.95, f"R平方:{r_squared:.3f},残差的标准差:{np.sqrt(np.sum(residuals**2)/x.shape[0]):.3f}", transform=ax2.transAxes, fontsize=14,verticalalignment='top', bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))

# 0.06为观测阈值
fig.suptitle("点云拟合分析", fontsize=20) 
plt.tight_layout()
plt.show()