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

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

// **********
// Constants
// **********

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

// circle
circcenter = [2;2]; // center of the circle
radius = 2; // radius of the circle
sigmafuz = 1; // variance of the gaussian blur (width)
intensity = 128; // 256 levels of grey (1-256)

// background, 256 levels of grey (1-256)
varbkg = 20; // variance of the background noise
meanbkg = 64; // mean of the background

// **********
// Functions
// **********

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

// **********
// Main program
// **********

// coordinates
X = (xmin:step:xmax)';
Y = (ymin:step:ymax)';
sizeX = size(X, '*');
sizeY = size(Y, '*');

// picture background
base = ones(sizeX, sizeY);
background = meanbkg*base + varbkg*rand(base, 'normal')

// circle drawing
circle = zeros(base);
for i = 1:sizeX
    for j = 1:sizeY
        dist_center = sqrt((X(i)-circcenter(1))^2 + (Y(j)-circcenter(2))^2);
        circle(i,j) = intensity*gauss(dist_center, radius, sigmafuz);
    end
end

// total image
image = floor(circle + background);

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

// saving
write('X_noisy_circle.txt', X, '(F4.1)');
write('Y_noisy_circle.txt', Y, '(F4.1)');
write('image_noisy_circle.txt', image, '(I3)');

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

clear;
clf;
chdir('mypath/')

// **********
// Constants
// **********

threshold = 140; // threshold of detection on 256 levels of grey

// **********
// functions
// **********

// Recomposition of the picture matrix
function [matrixpic] = picture_matrix(sizeX, sizeY, vector)
    // sizeX, sizeY : dimensions of the picture, in nb of pts (integers)
    // vector : list of intensities (sizeX*sizeY)
    // matrice : intensities in a sizeX x sizeY matrix
    for i = 1:sizeX
        for j = 1:sizeY
            matrixpic(i,j) = vector((i - 1)*sizeX + j);
        end
    end
endfunction

// Threshold detection
function [X, Y] = detection(Xpicture, Ypicture, intensite, threshold)
    // Xpicture, Ypicture : scales (vectors)
    // intensité : grey levels (matrix)
    // threshold : detection threshold (scalar)
    // X, Y : coordinnates of the points which intensity is above the threshold
    // (vectors)
    sizeX = size(Xpicture,'*');
    sizeY = size(Ypicture,'*');
    k = 1;
    for i = 1:sizeX
        for j = 1:sizeY
            if (intensite(i,j) > = threshold) then
                X(1,k) = Xpicture(i);
                Y(1,k) = Ypicture(j);
                k = k+1;
            end
        end
    end
endfunction

// initial variable substitution

function [d] = varsubst(A)
    // A : m*n matrix
    // n : nb of points ; m : dimension of the space (2 in a plane)
    // each column has the coordinates of the experimental points
    // d : n*1 matrix
    [m, n] = size(A);
    for i = 1:n
        d(i, 1) = norm(A(:,i)).^2;
    end
endfunction

// rLinear circle fitting

function [C, r, R] = circle_fitting(X, Y)
    // X, Y: experimental points (line vectors)
    // C : coordinnates of the center of the circle (vector)
    // r : circle radius (scalar)
    A = [X ; Y];
    d = varsubst(A);
    // xstar = 2*xcenter
    // ystar = 2*ycenter
    // zstar = r^2 - (xcenter^2 + ycenter^2)
    [xystar, zstar, sigma] = reglin(A, d');
    C = 0.5*xystar;
    r = sqrt(zstar + C(1)^2 + C(2)^2);
    R = sigma;
endfunction

// **********
// Main program
// **********

// Data reading

raw_picture = read('image_noisy_circle.txt', -1, 1); // intensity
Xdef = read('X_noisy_circle.txt', -1, 1);
Ydef = read('Y_noisy_circle.txt', -1, 1);
xmin = min(Xdef); xmax = max(Xdef);
ymin = min(Ydef); ymax = max(Ydef);
sizeX = size(Xdef, '*');
sizeY = size(Ydef, '*');

// intensities -> matrixs
picture = picture_matrix(sizeX, sizeY, raw_picture);

// Display of the noisy picture
isoview(xmin, xmax, ymin, ymax);
cmap = [graycolormap(256);
    0, 255, 0; // green
    225, 0, 0]; // red
xset('colormap', cmap);
grayplot(Xdef, Ydef, picture);

// detection of the pixels above the threshold
[Xcercle, Ycercle] = detection(Xdef, Ydef, picture, threshold);
plot(Xcercle, Ycercle, 'g+');

// regression
[center, radius, mismatch] = circle_fitting(Xcercle, Ycercle);

// writing

text = 'C('+string(center(1))+' ; '+string(center(2))+') ;...
    r = '+string(radius);

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

xstring(-0.5, 4.5, text);
txt = gce(); txt2.font_foreground = 1; // black colour (cmap)

// Drawing the model

plot(center(1), center(2), 'r+')
diametre = 2*radius;
xarc(center(1) - radius, center(2) + radius,...
    diametre, diametre,...
    0, 360*64)
b = gce(); // circle
b.foreground = 258; // red colour (cmap)

// Display of parameters

print(%io(2), center)
print(%io(2), radius)
print(%io(2), mismatch)