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

2024-02-10 19:23:34

文章目录

1 多项式曲线拟合函数 polyfit

polyfit — 多项式曲线拟合

主要有 3 种重载方式

NO.1 给定坐标点 ( x , y ) (x,y) (x,y) 和拟合阶数 n n n，返回 n + 1 n+1 n+1 个按将幂排列的最小二乘拟合多项式 p ( x ) p(x) 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(x)=p_1x^n+p_2x^{n-1}+...+p_nx+p_{n+1} p(x)=p1xn+p2xn−1+...+pnx+pn+1，对应的多项式系数为 p = [ p 1 , p 2 , . . . , p n , p n + 1 ] p=[p_1,p_2,...,p_n,p_{n+1}] p=[p1,p2,...,pn,pn+1]

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

[p,S![在这里插入图片描述](https://www.icode9.com/i/ll/?i=3c201e6c05a343bc80b120dd00eb09f4.jpg?,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)

m u mu mu，一个二元列向量，包含中心化值和缩放值。 m u ( 1 ) mu(1) mu(1) 是均值 m e a n ( x ) mean(x) mean(x)， m u ( 2 ) mu(2) mu(2) 是标准差 s t d ( x ) std(x) std(x)。使用这些值时，polyfit 将 x x 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 个位于区间 [ x 1 , x 2 ] [x1,x2] [x1,x2] 等间距的行向量，且间距为 x 2 − x 1 n − 1 \cfrac{x2-x1}{n-1} n−1x2−x1

y = linspace(x1,x2,n)

2.2 多项式与二维点集拟合

代码：

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)));

结果展示：

小注：

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 y=p_1x^n+p_2x^{n-1}+...+p_nx+p_{n+1} y=p1xn+p2xn−1+...+pnx+pn+1
x x x 可以是矩阵或向量。在任一情况下，polyval 计算每个 x 元素处的 p。

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

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