Curve fitting: general Parabola --> translated and rotated coordinate system

Hi guys, I have a set of points through which I need to fit a general parabola (to find out its vertex, angle of inclination of directrix from x-axis and parameter a). But problem is these points are not well conditioned and they are in a rotated and translated coordinate system, though not inclined set of axes(angle between axes ==90 degree). Does anyone how to do that in matlab?

Points are something like this:

x= [-5.2:.4:0,-5.2:.4:0]';
y=[4,3.6,5.2,6,6.4,6.8,6.8,7.2,7.6,10.4,10.4,11.2,11.6,11.2,3.6,3.2,2.8,2.4,2,1.6,1.2,0.4,0,-.4,-1.2,-1.2,-1.6,-1.6]';

Just do plot(x,y,'.') and you will see what I am talking about.

Any help would be much appreciated.

Thanks, Shubham

Answer 1

[Matt J]

Something like this, maybe:

 function [vertex,theta, a] = myfit(x,y)

    xy=[x(:),y(:)].';

    theta= fminsearch(@(theta) cost(theta,xy), 45);    

    [~,coeffs]=cost(theta,xy);

    [a,b,c]=deal(coeffs(1),coeffs(2), coeffs(3));

        xv=-b/2/a;

     vertex=R(-theta)*[xv;polyval(coeffs,xv)];

 function [Cost,coeffs,xx,yy] = cost(theta,xy)

    Rxy=R(theta)*xy;

    [xx,idx]=sort(Rxy(1,:));
    yy=Rxy(2,idx);

    [coeffs,S]=polyfit(xx,yy,2);

    Cost=S.normr;

 function Rmat=R(theta)

     Rmat=[cosd(theta), -sind(theta); sind(theta), cosd(theta)];

[Shubham]

Hey Matt

Thanks for your reply. Could you please explain your code, what algorithm you are using here? Before this, I was trying to fit general equation of parabola (A*x^2 + 2*B*xy +C*y^2 + 2*g*x + 2*f*y + c) with A= cos(theta), C= sin(theta) and B.^2= 2AC and then try to figure out its parameters. It was fitting a vertical parabola with that. With your code, what is the theta we are getting (78 degrees)?? Is there any way I could visualize fitted parabola.

Regards and Thanks,

Shubham

[Matt J]

Hi Shubham,

In my code, theta is the angle of the directrix to the x-axis or equivalently, the angle of the parabola's axis of symmetry to the y-axis. Positive theta rotates the directrix/symmetry axis clockwise.

My code basically reduces the fit to a minimization over theta. For a given theta, the x,y coordinates are un-rotated and POLYFIT is applied, leading to a polyfit error. The code minimizes over this error with fminsearch (Nelder-Mead). In the code, you will also see

    [a,b,c]=deal(coeffs(1),coeffs(2), coeffs(3));

Here a,b,c are the coefficients of the unrotated parabola. To plot, you could derive coordinates for the unrotated parabola and then rotate by theta clockwise.

Your approach with A*x^2 + 2*B*xy +C*y^2 + 2*g*x + 2*f*y + c looks like it could work, too, but it requires a constrained solver. It would be interesting to see which produces the better fit. If nothing else, my approach might be a good way of deriving an initial parameter guess for yours. A few remarks, though,

1. I think you have a mistake in B.^2= 2AC, which might account for some of your problems. I think it should really be B.^2= 4AC.
2. I don't know how you got to A=cos(theta) and C=sin(theta). If that were true, then in the unrotated case theta=90 you would always reduce to a conventional parabola of the form y=x^2+p*x+q whose second derivative is always fixed at 2. Surely, you want it to reduce to something more general, like y=s*x^2+p*x+q.

[Shubham]

Thanks Matt for such a beautiful and clever solution. It works perfectly.

[Carl Witthoft]

I think there's a bug in the proposed solution: if the original parabola is rotated in the opposite direction, the default starting value (+45) will find a local minimum at a value roughly 90 degrees off, and not very accurately there. Probably the initial value should be selected after determining in which quadrant(s) the majority of the input dataset reside(s). The reason for this problem is that the parabola x=y^2 will yield a minimum more or less for a degenerate parabola, i.e. straight line, along the axis of the parabola.

[Christoph]

First let me thank you very much for your provided code, it helped me a lot. Second this would be an idea how to draw the ellipse afterwards,

a = coeffs(1);
b = coeffs(2);
c = coeffs(3);

t=-1000:1:+1000;
xEllipse = -(a.*t.^2 + b*t + c) * sin(-theta) + cos(-theta)*t;
yEllipse = +(a.*t.^2 + b*t + c) * cos(-theta) + sin(-theta)*t;
plot(xEllipse,yEllipse,'b.');

Example Image or curve fitting, red the original points fitted with a polynomial for each side, plotted in blue.

And my question would be: If I am skipping the fminsearch for minimizing the cost function and finding the best theta angle and just use the 45° theta angle for the coefficient calculation I can also achieve good, sometimes also better cost results as while using the fminseach function. I am not able to understand this at the moment, could you provide a clarification for me. Thank you very much, Christoph