// **********
// Initialisation
// **********

clear;

// **********
// Données
// **********

X0 = [1, 2, 5, 7, 9, 6, 3, 8];
Y0 = [7, 6, 8, 7, 5, 7, 2, 4];

// **********
// Fonctions
// **********

function [a] = regression_elliptique(X, Y) // Fitzgibbon
    // méthode de la distance algébrique
    // X, Y : points expérimentaux, matrices colonnes réelles
    // a : coefficients de la formule quadratique (matrice colonne réelle)
    D = [X.*X, X.*Y, Y.*Y, X, Y, ones(X)]; // matrice de conception (design m.)
    S = D'*D; // matrice de dispersion (scatter m.)
    C = zeros(6,6);
    C(1,3) = 2; C(2,2) = -1; C(3,1) = 2; // matrice de contrainte
    [vecpropres, valpropres] = spec(inv(S)*C); // détermination du
    // système propre
    if imag(vecpropres) <> 0 then
        error('Les vecteurs propres contiennent des valeurs complexes')
    end
    if imag(valpropres) <> 0 then
        error('Les valeurs propres contiennent des valeurs complexes')
    end
    vecpropres = real(vecpropres); // complexes -> réels
    valpropres = real(valpropres);
    [PosLigne, PosColonne] = find((valpropres > 0 & ~isinf(valpropres)));
    // recherche les indices des valeurs propres positives
    a = vecpropres(:, PosLigne); // vecteur propre correspondant
endfunction

function [phi]=trouve_rotation(A)
    // A : coefficients de la formule quadratique (matrice colonne réelle)
    // phi : angle que fait un axe de l'ellipse avec x (radians)
    delta = 1 - 1/(1 + (A(3) - A(1))^2/A(2)^2);
    absphi = acos(sqrt((1 + sqrt(delta))/2));
    signephi = sign(A(2)*(cos(absphi)^2 - sin(absphi)^2)/(A(1) - A(3)));
    phi = signephi*absphi;
endfunction

function [x,y]=trouve_centre(A)
    // A : coefficients de la formule quadratique (matrice colonne réelle)
    // x, y : coordonées du centre de l'ellipse (réels)
    delta = A(2)^2 - 4*A(1)*A(3);
    x = (2*A(3)*A(4) - A(2)*A(5))/delta;
    y = (2*A(1)*A(5) - A(2)*A(4))/delta;
endfunction

function [rx, ry]=trouve_rayons(a, phi, xc, yc)
    // a : coefficients de la formule quadratique (matrice colonne réelle)
    // phi : angle que fait un axe de l'ellipse avec x
    // xc, yc : coordonnées du centre de l'ellipse
    // rx, ry : rayons (grand et petit demi-grands axes) de l'ellipse
    A = [a(1), a(2)/2 ; a(2)/2, a(3)];
    Q = rotate([1,0;0,1], phi); // matrice de rotation
    t = [xc;yc]; // matrice de translation
    Abar = Q'*A*Q;
    b = [a(4);a(5)];
    bbar = (2*t'*A + b')*Q;
    c = a(6);
    cbar = t'*A*t + b'*t + c;
    rx = sqrt(-cbar/Abar(1,1));
    ry = sqrt(-cbar/Abar(2,2));
endfunction

function [] = trace_ellipse(xc, yc, a, b, phi)
    // trace l'ellipse de centre (xc, yc)
    // de rayons a et b et tournée de phi
    pas = 0.1;
    t = 0:pas:%pi/2;
    X = a*cos(t);
    Y = b*sin(t);
    n = 4*size(X,'*');
    XY1 = [X, -flipdim(X,2), -X, flipdim(X,2);...
        Y, flipdim(Y,2), -Y, -flipdim(Y,2)];
    XY = rotate(XY1, phi) + [xc*ones(1,n);yc*ones(1,n)];
    xpoly(XY(1,:), XY(2,:));
endfunction

// **********
// Programme principal
// **********

// lecture des données

Xdef = X0';
Ydef = Y0';

// Régression
aopt = regression_elliptique(Xdef, Ydef);

// affichage des paramètres
disp(aopt)

phi = trouve_rotation(aopt);
phideg = phi*180/%pi;
[xc, yc] = trouve_centre(aopt);
[a, b] = trouve_rayons(aopt, phi, xc, yc);
disp('phi = '+string(phi)+' rad = '+string(phideg)+'°.');
disp('C('+string(xc)+' ; '+string(yc)+').');
disp('a = '+string(a)+' ; b = '+string(b)+'.');

// tracé
clf;

plot(Xdef, Ydef, 'b+')
isoview(0, 10, 1, 9);
plot(xc, yc, 'r+')
trace_ellipse(xc, yc, a, b, phi);
ell = gce();
ell.foreground = 5;

function [a] = regression_elliptique(X, Y) // Halir
    // méthode de la distance algébrique
    // X, Y : points expérimentaux, matrices colonnes réelles
    // a : coefficients de la formule quadratique (matrice colonne réelle)
    D1 = [X.*X, X.*Y, Y.*Y];
    D2 = [X, Y, ones(X)];
    // matrices de conception (design m.)
    S1 = D1'*D1; 
    S2 = D1'*D2;
    S3 = D2'*D2;
    // matrices de dispersion (scatter m.)
    T = -inv(S3)*S2';
    N = S1+ S2*T;
    M = [0.5*N(3, :) ; -N(2,:) ; 0.5*N(1, :)]; // mult par inv(C1) à gauche
    // matrice de dispersion réduite
    [vecpropres, valpropres] = spec(M);
    vep = real(vecpropres);
    // détermination du système propre
    condition = 4*vep(1, :).*vep(3, :) - vep(2, :).^2;
    // évaluation de a'Ca
    a1 = vep(:, find(condition > 0));
    a = [a1 ; T*a1]; // vecteur propre correspondant à la solution
endfunction