如何将椭圆锥拟合到一组数据?

How to fit an elliptic cone to a set of data?

我有一组 3d 数据(300 个点)创建了一个看起来像两个相互连接的圆锥体或椭圆体的表面。我想要一种方法来找到最适合该数据集的椭球体或圆锥体的方程。回归方法不重要,越简单越好。我基本上需要一种方法、代码或 matlab 函数来计算这些数据的椭圆方程的常数。

matlab 函数fit 可以采用任意拟合表达式。计算参数需要一些时间,但可以做到。

您首先要创建一个 fittype 对象,其中包含一个代表您预期形式的字符串。您需要自己计算出最符合您期望的表达式,我将以 a cone expression from the Mathworld site 为例并重新排列 z

ft = fittype('sqrt((x^2 + y^2)/c^2) + z_0', ...
     'independent', {'x', 'y'}, 'coeff', {'c', 'z_0'});

如果它是一个简单的形式,matlab 可以算出哪些是变量,哪些是系数,但是对于像这样更复杂的东西,你会想要帮助它。

'fitoptions' 对象保存方法的配置:根据您的数据集,您可能需要花一些时间指定上限和下限、起始值等。

fo = fitoptions('Upper', [one, for, each, of, your, coeffs, in, the, order, they, appear, in, the, string], ...
     'Lower', [...], `StartPoint', [...]);

然后得到输出

[fitted, gof] = fit([xvals, yvals], zvals, ft, fo);

警告:我已经用 2D 数据集做了很多,docs 声明它适用于三个但我自己没有这样做所以上面的代码可能不起作用,检查文档以确保你的语法正确。

可能值得从一个简单的拟合表达式开始,一些线性的东西,这样您就可以让您的代码正常工作。然后将表达式换成圆锥体并反复尝试,直到得到符合您预期的结果。

在你完成拟合后,一个好技巧是你可以对你在拟合中使用的字符串表达式使用 eval 函数来计算字符串的内容,就好像它是一个 matlab 表达式一样.这意味着您需要具有与字符串表达式中的变量和系数同名的工作区变量。

您也可以尝试使用 fminsearch,但为了避免落入局部最小值,您将需要一个好的起点,给定系数的数量(尝试消除其中的一些系数)。

这是一个二维椭圆的例子:

% implicit equation
fxyc = @(x, y, c_) c_(1)*x.^2 + c_(2).*y.^2 + c_(3)*x.*y + c_(4)*x + c_(5).*y - 1; % free term locked to -1

% solution (ellipse)
c_ = [1, 2, 1, 0, 0]; % x^2, y^2, x*y, x, y (free term is locked to -1)

[X,Y] = meshgrid(-2:0.01:2);
figure(1);
fxy = @(x, y) fxyc(x, y, c_);
c = contour(X, Y, fxy(X, Y), [0, 0], 'b');
axis equal;
grid on;
xlabel('x');
ylabel('y');          
title('solution');          

% we sample the solution to have some data to fit
N = 100; % samples
sample = unique(2 + floor((length(c) - 2)*rand(1, N)));
x = c(1, sample).';
y = c(2, sample).';

x = x + 5e-2*rand(size(x)); % add some noise
y = y + 5e-2*rand(size(y));

fc = @(c_) fxyc(x, y, c_); % function in terms of the coefficients
e = @(c) fc(c).' * fc(c); % squared error function

% we start with a circle
c0 = [1, 1, 0, 0, 0];

copt = fminsearch(e, c0)

figure(2);
plot(x, y, 'rx');
hold on
fxy = @(x, y) fxyc(x, y, copt);
contour(X, Y, fxy(X, Y), [0, 0], 'b');
hold off;
axis equal;
grid on;
legend('data', 'fit');
xlabel('x');                    %# Add an x label
ylabel('y');          
title('fitted solution');