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

clear;
clf;
chdir('monchemin/');

// **********
// Constantes
// **********

// image
xmin = -1;
xmax = 5;
ymin = -1;
ymax = 5;
pas=0.1;

// cercle
centre = [2;2]; // centre du cercle
rayon = 2; // rayon du cercle
sigmaflou = 1; // variance du flou gaussien (largeur)
intensite = 128; // 256 niveau de gris (0-255)

// fond, 256 niveau de gris (0-255)
varfond = 20; // variance du bruit de fond
moyfond = 64; // moyenne du fond

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

function [y] = gauss(x, esperance, ecart_type)
    y = 1/(ecart_type*sqrt(2*%pi))*exp(-((x - esperance)/ecart_type)^2/2);
endfunction

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

// coordonnées
X = (xmin:pas:xmax)';
Y = (ymin:pas:ymax)';
tailleX = size(X, '*');
tailleY = size(Y, '*');

// fond d'image
base = ones(tailleX, tailleY);
fond = moyfond*base + varfond*rand(base, 'normal')

// tracé du cercle
cercle = zeros(base);
for i = 1:tailleX
    for j = 1:tailleY
        dist_centre = sqrt((X(i)-centre(1))^2 + (Y(j)-centre(2))^2);
        cercle(i,j) = intensite*gauss(dist_centre, rayon, sigmaflou);
    end
end

// image globale
image = floor(cercle + fond);

// Affichage
isoview(xmin, xmax, ymin, ymax)
cmap = graycolormap(256);
xset('colormap', cmap);
grayplot(X, Y, image);

// enregistrement
write('X_cercle_bruite.txt', X, '(F4.1)');
write('Y_cercle_bruite.txt', Y, '(F4.1)');
write('image_cercle_bruite.txt', image, '(I3)');

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

clear;
clf;
chdir('monchemin/')

// **********
// Constantes
// **********

seuil = 140; // seuil de détection sur 256 niveaux de gris

// **********
// fonctions
// **********

// Recomposition de la matrice image
function [matrice]=matrice_image(tailleX, tailleY, vecteur)
    // tailleX, tailleY : dimensions de l'image en nb de pts (entiers)
    // vecteur : liste des intensités (tailleX*tailleY)
    // matrice : intensités sous forme de matrice tailleX x tailleY
    for i = 1:tailleX
        for j = 1:tailleY
            matrice(i,j) = vecteur((i - 1)*tailleX + j);
        end
    end
endfunction

// Détection de seuil
function [X, Y] = detection(Ximage, Yimage, intensite, seuil)
    // Ximage, Yimage : échelles (vecteurs)
    // intensité : niveaux de gris (matrice)
    // seuil : seuil de détection (scalaire)
    // X, Y : coordonnées des points dont l'intensité dépasse le seuil
    // (vecteurs)
    tailleX = size(Ximage,'*');
    tailleY = size(Yimage,'*');
    k = 1;
    for i = 1:tailleX
        for j = 1:tailleY
            if (intensite(i,j) > = seuil) then
                X(1,k) = Ximage(i);
                Y(1,k) = Yimage(j);
                k = k+1;
            end
        end
    end
endfunction

// changement de variable initial

function [d]=chgtvar(A)
    // A : matrice m*n
    // n : nb de points ; m : dimension de l'espace (2 dans le plan)
    // chaque colonne contient les coordonées des points exp
    // d : matrice n*1
    [m, n] = size(A);
    for i = 1:n
        d(i, 1) = norm(A(:,i)).^2;
    end
endfunction

// régression circulaire linéaire

function [C, r, R]=regression_circulaire(X, Y)
    // X, Y : points expérimentaux (vecteurs ligne)
    // C : coordonnées du centre du cercle (vecteur)
    // r : rayon du cercle (scalaire)
    A = [X ; Y];
    d = chgtvar(A);
    // xetoile = 2*xcentre
    // yetoile = 2*ycentre
    // zetoile = r^2 - (xcentre^2 + ycentre^2)
    [xyetoile, zetoile, sigma] = reglin(A, d');
    C = 0.5*xyetoile;
    r = sqrt(zetoile + C(1)^2 + C(2)^2);
    R = sigma;
endfunction

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

// lecture des données

image_brute = read('image_cercle_bruite.txt', -1, 1); // intensité
Xdef = read('X_cercle_bruite.txt', -1, 1);
Ydef = read('Y_cercle_bruite.txt', -1, 1);
xmin = min(Xdef); xmax = max(Xdef);
ymin = min(Ydef); ymax = max(Ydef);
tailleX = size(Xdef, '*');
tailleY = size(Ydef, '*');

// mise en matrice des intensités
image = matrice_image(tailleX, tailleY, image_brute);

// Affichage de l'image bruitée
isoview(xmin, xmax, ymin, ymax);
cmap = [graycolormap(256);
    0, 255, 0; // vert
    225, 0, 0]; // rouge
xset('colormap', cmap);
grayplot(Xdef, Ydef, image);

// détection des pixels dépassant le seuil
[Xcercle, Ycercle] = detection(Xdef, Ydef, image, seuil);
plot(Xcercle, Ycercle, 'g+');

// regression
[centre, rayon, ecart] = regression_circulaire(Xcercle, Ycercle);

// écriture

texte = 'C('+string(centre(1))+' ; '+string(centre(2))+') ;...
    r = '+string(rayon);

// fond blanc (4x le même texte décalé)
xstring(-0.48, 4.5, texte);
txt = gce(); txt.font_foreground = 256; // couleur blanc (cmap)
xstring(-0.52, 4.5, texte);
txt = gce(); txt.font_foreground = 256;
xstring(-0.5, 4.52, texte);
txt = gce(); txt.font_foreground = 256;
xstring(-0.5, 4.48, texte);
txt = gce(); txt.font_foreground = 256;

xstring(-0.5, 4.5, texte);
txt = gce(); txt2.font_foreground = 1; // couleur noir (cmap)

// tracé du modèle

plot(centre(1), centre(2), 'r+')
diametre = 2*rayon;
xarc(centre(1) - rayon, centre(2) + rayon,...
    diametre, diametre,...
    0, 360*64)
b = gce(); // cercle
b.foreground = 258; // couleur rouge (cmap)

// Affichage des paramètres

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