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

clear;
chdir('monchemin/');

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

phiinit = 30/180*%pi;
xcinit = 4; ycinit = 3;
ainit = 4; binit = 2;
varbruit = 0.2;

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

function [XY] = cree_points(xc, yc, a, b, phi, var)
    // trace l'ellipse de centre (xc, yc)
    // de rayons a et b et tournée de phi
    pas = 0.5;
    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)];
    bruit = var*rand(XY1, 'normal')
    XY = rotate(XY1, phi) + [xc*ones(1,n);yc*ones(1,n)] + bruit;
endfunction

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

// Génération des données

XYdef = cree_points(xcinit, ycinit, ainit, binit, phiinit, varbruit);
Xdef = XYdef(1,:)';
Ydef = XYdef(2,:)';

// enregistrement

write('ellipse_bruitee.txt', XYdef')

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

clear;
chdir('monchemin/');

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

function [a] = regression_elliptique(X, Y) 
    // 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

XYdef = read('ellipse_bruitee.txt', -1, 2);
Xdef = XYdef(:,1);
Ydef = XYdef(:,2);

// 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, 8, -0.5, 6);
plot(xc, yc, 'r+');
trace_ellipse(xc, yc, a, b, phi);
ell = gce();
ell.foreground = 5;