ellipseConic2param

% estimates values of parametric represenation of a 2D ellipse

 Usage:
   p = ellipseConic2param(C)

   C: double 3x3 ellipse-conic, det(Chh)!>0
   p: double 5x1, [x0,y0,a,b,phi]

 10.12.13 Susanne Wenzel
 wenzel@igg.uni-bonn.de

0001 %% estimates values of parametric represenation of a 2D ellipse
0002 %
0003 % Usage:
0004 %   p = ellipseConic2param(C)
0005 %
0006 %   C: double 3x3 ellipse-conic, det(Chh)!>0
0007 %   p: double 5x1, [x0,y0,a,b,phi]
0008 %
0009 % 10.12.13 Susanne Wenzel
0010 % wenzel@igg.uni-bonn.de
0011 
0012 function p = ellipseConic2param(C)
0013 
0014 % help values
0015 Chh = C(1:2,1:2);
0016 
0017 if det(Chh)<=0
0018     keyboard
0019 end
0020 
0021 Chhinv = inv(Chh);
0022 ch0 = C(1:2,  3);
0023 c00 = C(  3,  3);
0024 
0025 % centre point
0026 x0 = -Chhinv*ch0;
0027 
0028 % eigenvalue decomposition of normalized homogeneous part of the conic
0029 Chhe = -Chh/(c00-ch0'*Chhinv*ch0); %#ok<MINV>
0030 [R,Lambda] = eig( Chhe );
0031 
0032 l = diag(Lambda);
0033 [minlambda,idx_min] = min(l);
0034 maxlambda = max(l);
0035 
0036 % main radii
0037 a = sqrt(1/minlambda);
0038 b = sqrt(1/maxlambda);
0039 
0040 if  abs(minlambda-maxlambda)<eps
0041     % its a circle, orientation 0
0042     phi = 0;
0043 else
0044     phi = atan2(R(2,idx_min),R(1,idx_min));
0045     if phi<0
0046         phi = phi+2*pi;
0047     end
0048 end
0049 
0050  p = [x0',a,b,phi];