matlab - 如何将椭圆锥拟合到一组数据？

Question

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

Answer 1

matlab 函数fit可以采用任意拟合表达式。需要花点时间弄清楚参数，但可以做到。

您将首先创建一个fittype对象，该对象具有一个代表您预期形式的字符串。您需要自己计算出最符合您期望的表达式，我将从 Mathworld 站点获取一个圆锥表达式作为示例，并将其重新排列为 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 数据集做了很多，文档声明它适用于三个，但我自己没有这样做，所以上面的代码可能不起作用，检查文档以确保你的语法正确。

从一个简单的拟合表达式开始可能是值得的，一些线性的，这样你就可以让你的代码工作。然后将表达式换成锥体并玩弄直到你得到看起来像你期望的东西。

在你适应之后，一个很好的技巧是你可以使用eval你在适应中使用的字符串表达式上的函数来评估字符串的内容，就好像它是一个 matlab 表达式一样。这意味着您需要有与字符串表达式中的变量和系数同名的工作区变量。

Answer 2

您也可以尝试使用fminsearch，但为了避免陷入局部最小值，您需要一个良好的起点，因为您需要考虑系数的数量（尝试消除其中的一些）。

这是一个带有 2D 椭圆的示例：

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