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

Question

I have a set of 3d data (300 points) that create a surface which looks like two cones or ellipsoids connected to each other. I want a way to find the equation of a best fit ellipsoid or cone to this dataset. The regression method is not important, the easier it is the better. I basically need a way, a code or a matlab function to calculate the constants of the elliptic equation for these data.

Answer 1

The matlab function fit can take arbitrary fit expressions. It takes a bit of figuring out the parameters but it can be done.

You would first create a fittype object that has a string representing your expected form. You'll need to work out the expression yourself that best fits what you're expecting, I'm going to take a cone expression from the Mathworld site for an example and rearrange it for z

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

If it's a simple form matlab can work out which are the variables and which the coefficients but with something more complex like this you'd want to give it a hand.

The 'fitoptions' object holds the configuration for the methods: depending on your dataset you might have to spend some time specifying upper and lower bounds, starting values etc.

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

then get the output

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

Caveat: I've done this plenty with 2D datasets and the docs state it works for three but I haven't done that myself so the above code might not work, check the docs to make sure you've got your syntax right.

It might be worth starting with a simple fit expression, something linear, so that you can get your code working. Then swap the expression out for the cone and play around until you get something that looks like what you're expecting.

After you've got your fit a good trick is that you can use the eval function on the string expression you used in your fit to evaluate the contents of the string as if it was a matlab expression. This means you need to have workspace variables with the same names as the variables and coefficients in your string expression.

Answer 2

You can also try with fminsearch, but to avoid falling on local minima you will need a good starting point given the amount of coefficients (try to eliminate some of them).

Here is an example with a 2D ellipse:

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