MATLAB点云处理（十七）：最小二乘多项式曲线拟合

Posted 2021-08-30 没事就要敲代码

文章目录

• • 1 多项式曲线拟合函数 polyfit
  • 2 代码实现
  • • 2.1 多项式与三角函数拟合
    • 2.2 多项式与二维点集拟合

1 多项式曲线拟合函数 polyfit

polyfit — 多项式曲线拟合

主要有 3 种重载方式

NO.1 给定坐标点 $(x,y)$ 和拟合阶数 $n$，返回 $n+1$ 个按将幂排列的最小二乘拟合多项式 $p(x)$ 的系数

p = polyfit(x,y,n)

对于多项式 $p(x)=p_1{x}^n+p_2{x}^{n-1}+...+p_{n}x+p_{n+1}$，对应的多项式系数为 $p=[p_1,p_2,...,p_n,p_{n+1}]$

NO.2 给定坐标点 $(x,y)$ 和拟合阶数 $n$，返回 $n+1$ 个按将幂排列的最小二乘拟合多项式 $p(x)$ 的系数，并且返回一个结构体 S，可用作 polyval 的输入来获取拟合误差。

[p,S![在这里插入图片描述](https://img-blog.csdnimg.cn/3c201e6c05a343bc80b120dd00eb09f4.jpg?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80NjA5ODU3Nw==,size_16,color_FFFFFF,t_70#pic_center)
= polyfit(x,y,n)

NO.3 采样中心化和缩放变换，改善多项式和拟合算法的数值属性。

[p,S,mu] = polyfit(x,y,n)

$mu$，一个二元列向量，包含中心化值和缩放值。$mu(1)$ 是均值 $mean(x)$，$mu(2)$ 是标准差 $std(x)$。使用这些值时，polyfit 将 $x$ 的中心置于零值处，并缩放为具有单位标准差

2 代码实现

2.1 多项式与三角函数拟合

代码：

clc;
clear;

% 在区间 [0,4*pi] 中沿正弦曲线生成 10 个等间距的点
x = linspace(0,4*pi,10);
y = sin(x);

% 使用 polyfit 将一个 7 次多项式与这些点拟合，并输出拟合系数
p = polyfit(x,y,7)

% 在更精细的网格上计算多项式并绘制结果图
x1 = linspace(0,4*pi);
y1 = polyval(p,x1);
figure;
plot(x,y,'o');
hold on;
plot(x1,y1);
hold off;
title('最小二乘拟合结果');
xlabel('X(m)');
ylabel('Y(m)');

输出结果：

p =

   -0.0001    0.0028   -0.0464    0.3702   -1.3808    1.9084   -0.1141    0.0002

小注：

linspace — 生成线性间距向量

有两种语法

• 生成 x1 和 x2 之间的 100 个等间距点 的行向量。

y = linspace(x1,x2)

• 生成 n 个位于区间 $[x1,x2]$ 等间距的行向量，且间距为 $\frac{x2-x1}{n-1}$

y = linspace(x1,x2,n)

2.2 多项式与二维点集拟合

示例1：简单的线性回归（直线拟合）

代码：

clc;
clear;

% 加载点云
ptCloud = pcread('test.pcd');

% 获取(x,y)坐标
xyPoints = ptCloud.Location;
x = xyPoints(:,1);
y = xyPoints(:,2);

% 最小二乘多项式拟合
p1 = polyfit(x,y,1);

% 更精细网格上计算多项式拟合的结果
x1 = linspace(min(x),max(x));
y1 = polyval(p1,x1); % 等价于 y1 = p1(1)*x1 + p1(2);

% 在区间 [min(x),max(x)] 中绘制多项式拟合结果
figure;
plot(x,y,'o');
hold on;
plot(x1,y1);
hold off;
title('最小二乘拟合结果');
xlabel('X(m)');
ylabel('Y(m)');
legend(sprintf('数据点'),sprintf('yn = (%f)x + (%f)',p1(1),p1(2)));

结果展示：

示例2：二次曲线拟合

代码：

clc;
clear;

% 加载点云
ptCloud = pcread('test2.pcd');

xyPoints = ptCloud.Location;

x = xyPoints(:,1);
y = xyPoints(:,2);

% 最小二乘多项式拟合
p1 = polyfit(x,y,2);

% 更精细网格上计算原始函数和多项式拟合
x1 = linspace(min(x),max(x));
y1 = polyval(p1,x1); 

% 在区间 [min(x),max(x)] 中绘制函数值和多项式拟合
figure;
plot(x,y,'o');
hold on;
plot(x1,y1);
hold off;
title('最小二乘拟合结果');
xlabel('X(m)');
ylabel('Y(m)');
legend(sprintf('数据点'),sprintf('yn = (%f)x^2 + (%f)x + (%f)',p1(1),p1(2),p1(3)));

输出结果：

小注：

polyval — 多项式计算

y = polyval(p,x)

返回在 x 处计算的 n 次多项式的值。输入参数 p 是长度为 n+1 的向量，其元素是按降幂排序的多项式系数。
$$y=p_1{x}^n+p_2{x}^{n-1}+...+p_{n}x+p_{n+1}$$
$x$ 可以是矩阵或向量。在任一情况下，polyval 计算每个 x 元素处的 p。

示例： 计算多项式 $p(x)=2x^2+3x+5$ 在 $x$ = 1、4 和 8 时的值：

p = [2 3 5];
x = [5 7 9];
y = polyval(p,x)

输出结果：

y =

    70   124   194