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

clear;

// données (Gander et coll.)

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

function [A] = regression_circulaire(X, Y)
    B = [X.^2 + Y.^2, X, Y, ones(X)];
    [U, S, V]= svd(B);
    u = V(:,4); a = u(1); b =[u(2); u(3)]; c = u(4);
    z = -b/2/a; r = sqrt(norm(z)^2 - c/a);
    A = [z', r];
endfunction

// **********
// programme principal
// **********

// lecture des données

// Xdef, Ydef : vecteurs ligne

Xdef = X0';
Ydef = Y0';

fenetre = scf(0); // création de la fenêtre graphique

// regression
[Aopt] = regression_circulaire(Xdef, Ydef)
centre = Aopt(1:2);
rayon = Aopt(3);
diam = 2*rayon;

print(%io(2), centre)
print(%io(2), rayon)

plot(Xdef, Ydef, 'b+');
plot(centre(1), centre(2), 'r+');
xarc(centre(1)-rayon, centre(2)+rayon, diam, diam, 0, 360*64)
xmin = min(centre(1)-rayon, min(Xdef));
xmax = max(centre(1)+rayon, max(Xdef));
ymin = min(centre(2)-rayon, min(Ydef));
ymax = max(centre(2)+rayon, max(Ydef));
isoview(xmin, xmax, ymin, ymax);