How to determine the general equation of an ellipse using lsqnonlin

Hello.

Given a set of points (xk,yk) how to determine de general equation of an ellipse: P(x,y) = Ax + Bxy + Cy^2 + Dx + Ey + F = 0

I've tried a few methods (direct least square, Kanatani's hypernormalization, Ellipse fitting with sampson distance). However, this time, I am trying to use de lsqnonlin command. The main difficulty is to guarantee that the coefficients really represent an ellipse, i.e., how do I add the condition: B^2-4*A*C<0 to the lsqnonlin ?

Thank you in advance.

Answer 1

[Matt J]

You would probably have to use fmincon, rather than lsqnonlin. It can accommodate more general nonlinear constraints like B^2-4*A*C<=0.

[António]

One of the problems with using fmincon is that B^2 - 4*A*C <=0 actually doesn't guarantee that the coefficients represent an ellipse. In the case B^2 - 4*A*C =0, I get a parabola.

[Matt J]

There's no way to avoid a parabola if that's what your data wants to be. The asymptotic minimum over B^2-4*A*C<0 will always be the same as the minimum over B^2-4*A*C<=0.

[António]

Correct. However, in my case, I'm not using data that completes a full ellipse, that's why I needed to guarantee the condition B^2-4*A*C<0. In fact, the same amount of points (limited points, of course) can belong to a circle, ellipse, parabola or hyperbola. Thank you for the answer.

[Matt J]

To bound things strictly away from zero, you can use fmincon to impose the constraint B^2-4*A*C+delta<=0 where delta>0 is a strictly positive parameter chosen by you. That's the only way to ensure, by any algorithm, that you don't get a parabola when a more weakly constrained fit would give one.

[António]

Hi. Here is another issue. fmincon asks for a scalar value in the return of the objective function. How do I do that, considering that I have 2 vectors of x and y values, and the objective function is given by P(x,y) ?

The call of fmincon looks like:

tHeta = fmincon(@fit_simp,tHeta0,[],[],[],[],[],[],@condfun,[xData,yData]);

And the objective function is like this:

function diff = fit_simp(tHeta)

      a = tHeta(1);
      B = tHeta(2);
      c = tHeta(3);
      d = tHeta(4);
      e = tHeta(5);
      f = tHeta(6);

      diff = a.*xData.^2 + B.*xData.*yData + c.*yData.^2 + d.*xData + e.*yData + f;

end

Hope you can help me. thanks

[Matt J]

It should look like

 tHeta = fmincon(@(tHeta) fit_simp(tHeta,xData,yData),...
                  tHeta0,[],[],[],[],[],[],@condfun,options);

with objective function

 function val = fit_simp(tHeta,xData,yData)

      a = tHeta(1);
      B = tHeta(2);
      c = tHeta(3);
      d = tHeta(4);
      e = tHeta(5);
      f = tHeta(6);

      diff = a.*xData.^2 + B.*xData.*yData + c.*yData.^2 + d.*xData + e.*yData + f;

      val=norm(diff(:))^2;

   end

[Matt J]

Even better, use more linear algebra ops,

    xData=xData(:);
    yData=yData(:);
    M=[xData.^2,  xData.*yData, yData.^2, xData ,yData ];
     M(:,end+1)=1;

    tHeta = fmincon(@(tHeta) norm(M*tHeta(:))^2,...
                      tHeta0,[],[],[],[],[],[],@condfun,options);

Note that this way, you can get rid of fit_simp.m completely.