如何将椭圆锥拟合到一组数据？答案

【问题标题】：How to fit an elliptic cone to a set of data?如何将椭圆锥拟合到一组数据？
【发布时间】：2015-07-13 08:49:28
【问题描述】：

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

【问题讨论】：

标签： matlab regression curve-fitting ellipse best-fit-curve

【解决方案1】：

您也可以尝试使用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');

【讨论】：

【解决方案2】：

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 表达式一样。这意味着您需要有与字符串表达式中的变量和系数同名的工作区变量。

【讨论】：