MATLAB: Fitting data to a superellipse

curve fittingellipseimage processingmathematicsnonlinear

greeting to the MATLAB community,
xdata=
I have two data sets; xdata and ydata and I need to fit a supperellipse to those data with this line of code.
    %%%%%%%%%%%%%%%%%%%
a=2;
b=5;
p=2;
xc=3;
yc=1;
t=0:pi/20:2*pi;
xdata=xc+a*cos(t); %example

ydata=yc+b*sin(t); %example
a0 = [10 10 2]; %inintial guess
options = optimset('Display','iter');
c = [xc yc]; %given
%
f = @(aa) (((xdata-c(1)))./aa(1)).^aa(3) + (((ydata-c(2)))./aa(2)).^aa(3) -1; %superellipse Equation
af = lsqnonlin(f, a0, [], [], options);
I used the x&y data of an ellipse with a=2,b=5,xc=3,yc=1
but it returns the following values: af= 9.999999837882097 + 0.000000000000000i 9.999999149008248 + 0.000000000000000i 1.999994717969291 – 0.000007892981223i
which should be colse to 2 5 2
any idea? thanks!

Best Answer

There are several problems here, as I expected.
1. An ellipse is not a function. So a sum of squares of residuals does not work well here. Best is to convert to polar coordinates.
2. Are you trying to estimate the center of the ellipse? You have fixed xc and yc, but in any fit to data, you would generally not know the exact center.
3. Most importantly, look carefully at what you are doing. Depending on which side of the center a point lies, the expression (x-xc) may be positive of negative. Then you try to raise the to a non-integer power. For example, what does this evaluate to in MATLAB? TRY IT!
(-1)^2.3
ans =
    0.58779 +    0.80902i
It returns a complex result. Any negative number will do that, when raised to a non-integer power.
So again, it is best to convert to polar coordinates. Express your ellipse as a function of angle and distance from the polar origin. Then everything will be positive, and no complex numbers will be generated. I'll look back in again if you need help with that, but I'll get you started.
Assume that the center is not known.
xc0 = mean(x);
yc0 = mean(y);
% Polar radius
r = sqrt((x - xc0)^2 + (y - yc0)^2);
% polar angle
th = atan2(y - yc0,x - xc0);
Now, you need to re-write your model in polar coordinates, thus r(th), given parameters for the ellipse.
Related Question