最小二乘法拟合球及其相关代码实现

最新推荐文章于 2024-07-08 13:42:21 发布

程序员小白学开发

最新推荐文章于 2024-07-08 13:42:21 发布

阅读量6.5k

点赞数 4

分类专栏： c++ 程序设计

c++ 同时被 2 个专栏收录

28 篇文章 1 订阅

订阅专栏

程序设计

8 篇文章 0 订阅

订阅专栏

当我们手中握有大量的数据时，对于二维的数据，我们会对他们进行直线拟合、对数拟合，圆曲线的拟合等等。这些拟合的方法都是运用的了非常古老而又非常有效的方法，即最小二乘法。
今天给大家介绍一种三维球面数据的拟合方法，该方法也是运用的最小二乘的方法。旨在使拟合的半径在均方意义下误差达到最小。

公式推导

设拟合后的球面的球心为(x_0,y_0,z_0)及半径r。
对于每一点拟合后估计的值与实际值的差值为:

则误差的平方和为:

注意E是x_0,y_0,z_0,r的函数。因此令E分别对x_0,y_0,z_0,r的偏导数等于0，即可求出x_0,y_0,z_0,r，有:

令：

则有：

由(1)-(4)得

由(2)-(4)得

由(3)-(4)得

写成矩阵的形式：

求解该矩阵得到x_0,y_0,z_0值，然后带入(4)式中得到r的值。

Matlab仿真

首先生成若干个球面数据，然后加入一定能量的噪声。最后利用上述的公式计算拟合后的球心坐标和球面半径，下面给出的是matlab仿真代码:

%最小二乘的方法进行拟合
clear all;
close all
clc;
R = 2;         %球面半径
x0 = 100;      %球心x坐标
y0 = 1;        %球心y坐标
z0 = 76;       %球心z坐标
%********************************生成随机球面数据************************************
alfa = 0:pi/50:pi;
sita = 0:pi/50:2*pi;
num_alfa = length(alfa);
num_sita = length(sita);
x = zeros(num_alfa,num_sita);
y = zeros(num_alfa,num_sita);
z = zeros(num_alfa,num_sita);
for i = 1:num_alfa
    for j = 1:num_sita
        x(i,j) = x0+R*sin(alfa(i))*cos(sita(j));
        y(i,j) = y0+R*sin(alfa(i))*sin(sita(j));
        z(i,j) = z0+R*cos(alfa(i));
    end
end

x = reshape(x,num_alfa*num_sita,1);
y = reshape(y,num_alfa*num_sita,1);
z = reshape(z,num_alfa*num_sita,1);
figure;
plot3(x,y,z,'*');
title('生成的没有噪声的球面数据');
%加入均值为0的高斯分布噪声 
amp = 0.1;
x = x + amp*rand(num_alfa*num_sita,1);
y = y + amp*rand(num_alfa*num_sita,1);
z = z + amp*rand(num_alfa*num_sita,1);
figure;
plot3(x,y,z,'*');
title('加入噪声的球面数据');
%*******************************************************************************************
%球面拟合算法
num_points = length(x);
x_avr = sum(x)/num_points;
y_avr = sum(y)/num_points;
z_avr = sum(z)/num_points;

xx_avr = sum(x.*x)/num_points;
yy_avr = sum(y.*y)/num_points;
zz_avr = sum(z.*z)/num_points;
xy_avr = sum(x.*y)/num_points;
xz_avr = sum(x.*z)/num_points;
yz_avr = sum(y.*z)/num_points;

xxx_avr = sum(x.*x.*x)/num_points;
xxy_avr = sum(x.*x.*y)/num_points;
xxz_avr = sum(x.*x.*z)/num_points;
xyy_avr = sum(x.*y.*y)/num_points;
xzz_avr = sum(x.*z.*z)/num_points;
yyy_avr = sum(y.*y.*y)/num_points;
yyz_avr = sum(y.*y.*z)/num_points;
yzz_avr = sum(y.*z.*z)/num_points;
zzz_avr = sum(z.*z.*z)/num_points;
%计算求解线性方程的系数矩阵
A = [xx_avr - x_avr*x_avr,xy_avr - x_avr*y_avr,xz_avr - x_avr*z_avr;
     xy_avr - x_avr*y_avr,yy_avr - y_avr*y_avr,yz_avr - y_avr*z_avr;
     xz_avr - x_avr*z_avr,yz_avr - y_avr*z_avr,zz_avr - z_avr*z_avr];
b = [xxx_avr - x_avr*xx_avr + xyy_avr - x_avr*yy_avr + xzz_avr - x_avr*zz_avr;
     xxy_avr - y_avr*xx_avr + yyy_avr - y_avr*yy_avr + yzz_avr - y_avr*zz_avr;
     xxz_avr - z_avr*xx_avr + yyz_avr - z_avr*yy_avr + zzz_avr - z_avr*zz_avr];
b = b/2;

resoult = inv(A)*b;

x00 = resoult(1);     %拟合出的x坐标
y00 = resoult(2);     %拟合出的y坐标
z00 = resoult(3);     %拟合出的z坐标
r = sqrt(xx_avr-2*x00*x_avr+x00*x00 + yy_avr-2*y00*y_avr+y00*y00 + zz_avr-2*z00*z_avr+z00*z00);   %拟合出的球半径r

 
 1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77