function [ParG, code] = AtoG (ParA);

%  Conversion of Algebraic parameters of a conic to its Geometric parameters

%   Algebraic parameters are coefficients A,B,C,D,E,F in the algebraic
%   equation     Ax^2 + 2Bxy + Cy^2 + 2Dx + 2Ey + F = 0

%   Geometric parameters depend on the type of the conic (ellipse, etc.)

%   Input:  ParA = (A,B,C,D,E,F)' is the vector of Algebraic parameters

%   Output: code is the code of the conic type (see below)
%           ParG is the vector of Geometric parameters (see below)

%   code:   1 - ellipse  2 - hyperbola  3 - parabola
%           4 - intersecting lines  5 - parallel lines
%           6 - coincident lines    7 -single line
%           8 - single point        9 - imaginary ellipse
%          10 - imaginary parallel lines
%          11 - "impossible" equation, 1=0 or -1=0 (no solutions)

%           Geometric parameters are determined only for 
%               the most common types of conics:

%   1. Ellipses:  canonical equation  x^2/a^2 + y^2/b^2 = 1
%                 a>b or a=b are the major and minor semi-axes

%       ParG = [Xcenter, Ycenter, a, b, AngleOfTilt]'

%   2. Hyperbolas:  canonical equation  x^2/a^2 - y^2/b^2 = 1
%                 a  is the distance from the center to each vertex
%                 b  is the vertical distance from each vertex to asymptotes

%       ParG = [Xcenter, Ycenter, a, b, AngleOfTilt]'

%   3. Parabolas:  canonical equation  y^2 = 2px
%                 p  is the distance from the focus to the directix

%       ParG = [Xcenter, Ycenter, p, AngleOfTilt]'

%   4. Intersecting lines:  canonical equation  
%                (cos(u)*x+sin(u)*y+d)*(cos(v)*x+sin(v)*y+e) = 0
%                   u,v are directional angles of the lines
%                   d,e are distances from the origin (0,0) to the lines

%       ParG = [u,d,v,e]'

%        Nikolai Chernov,  February 2012

ParA=ParA/norm(ParA);    %  normalize the given algebraic parameter vector

ParG = -1;   %  this will be returned for imaginary or degenerate conics

if (1+ParA(1)==1) && (1+ParA(2)==1) && (1+ParA(3)==1)  %  if no quadratic part...
    if (1+ParA(4)==1) && (1+ParA(5)==1)
        code = 11;     %  the "pole", extreme singularity
    else
        code = 7;      %  single line (no quadratic part, only linear part)
    end
    return;
end

M33 = [ParA(1) ParA(2) ParA(4);
       ParA(2) ParA(3) ParA(5);
       ParA(4) ParA(5) ParA(6)];   %    big, 3x3 matrix
M22 = [ParA(1) ParA(2);
       ParA(2) ParA(3)];           %  small, 2x2 matrix

det3x3 = det(M33);                 %    3x3 determinant
det2x2 = det(M22);                 %    2x2 determinant

if (1+det3x3 == 1)                   %  if the big matrix is singular...
    if (1+det2x2 == 1)
        dettwo = ParA(1)*ParA(6) - ParA(4)^2 + ParA(3)*ParA(6) - ParA(5)^2;
        if (dettwo > 0),  code = 10;  end   %  imaginary parallel lines
        if (dettwo < 0),  code =  5;  end   %  parallel lines
        if (1+dettwo == 1),  code = 6;  end   %  coincident lines
        return;
    end
    if (det2x2 > 0)
        code = 8;                %  single point
    else
        code = 4;                            %  intersecting lines
        [Q D] = eig(M33);    %  eigendecomposition
        [Dmax,Imax] = max(diag(D));
        [Dmin,Imin] = min(diag(D));
        Qmax = Q(1:3,Imax)*sqrt(abs(Dmax));
        Qmin = Q(1:3,Imin)*sqrt(abs(Dmin));
        Q1 = Qmax + Qmin;
        Q2 = Qmax - Qmin;
        theta1 = atan2(Q1(2),Q1(1));
        d1 = Q1(3)/norm(Q1(1:2));
        theta2 = atan2(Q2(2),Q2(1));
        d2 = Q2(3)/norm(Q2(1:2));
        ParG = [theta1; d1; theta2; d2];
    end
    return;
end

%          Next: non-degenrate types of conics

[Q D] = eig(M22);    %  eigendecomposition

U  = Q'*[ParA(4); ParA(5)];       %  orthogonal transformation

if (1+D(1,1)==1)||(1+D(2,2)==1)
    code = 3;                             %   parabola
    if (abs(D(1,1))>abs(D(2,2)))
        Uc1 = -U(1)/D(1,1);
        Uc2 = -(U(1)*Uc1 + ParA(6))/2/U(2);
        Center = Q*[Uc1; Uc2];
        p = -U(2)/D(1,1);
        Angle = atan2(Q(2,2), Q(1,2));
        ParG = [Center; p; Angle];
    else
        Uc2 = -U(2)/D(2,2);
        Uc1 = -(U(2)*Uc2 + ParA(6))/2/U(1);
        Center = Q*[Uc1; Uc2];
        p = -U(1)/D(2,2);
        Angle = atan2(Q(2,1), Q(1,1));
        ParG = [Center; p; Angle];
    end
    return;
end

Uc = -U./diag(D);
Center = Q*Uc;
H = -U'*Uc - ParA(6);
if (D(1,1)*D(2,2) < 0)
    code = 2;                            %   hyperbola
    if (D(1,1)*H > 0)
        a = sqrt(H/D(1,1));
        b = sqrt(-H/D(2,2));
        Angle = atan2(Q(2,1), Q(1,1));
        if (Angle < 0),  Angle = Angle + pi;    end
    else
        a = sqrt(H/D(2,2));
        b = sqrt(-H/D(1,1));
        Angle = atan2(Q(2,2), Q(1,2));
        if (Angle < 0),  Angle = Angle + pi;    end
    end
    ParG = [Center; a; b; Angle];
else
    if (H*D(1,1) <= 0)
        code = 9;                     %  imaginary ellipse
    else
        code = 1;                        %   ellipse
        a = sqrt(H/D(1,1));
        b = sqrt(H/D(2,2));
        Angle = atan2(Q(2,1), Q(1,1));
        if (Angle < 0),  Angle = Angle + pi;    end
        if a < b      %  making sure that a is major, b is minor
            ab = a; a = b; b = ab;
            Angle = Angle - pi/2;
            if (Angle < 0),  Angle = Angle + pi;  end
        end
        ParG = [Center; a; b; Angle];
    end
end
    
end    %  end of function AtoG