// Configuration du mode d'affichage dans la console
lines(0);

// Augmentation de la taille de la mémoire allouée à Scilab
stacksize("max")

// Contenu du fichier readimage.sce
function m = readimage(fil)
    lignecourante=[]
    m=[]
    u=file("open",fil,"old")
    head=read(u,1,1,"(a)")
    niveau=1;
    if head<>"P2"
        error("Not a PGM file")
    end

    c="#";
    while c=="#",
        lignecourante=read(u,1,1,"(a)");
        c=sscanf(lignecourante,"%c");
    end;
    [sz1,sz2]=sscanf(lignecourante,"%d%d");
    m=matrix(read(u,1,sz1*sz2),sz1,sz2);

    file("close",u)
    m = m';
endfunction

function N = dummy_convol2d(K, P)
    N = zeros(P);
    for x=1:size(P, "r")
        for y=1:size(P, "c")
            r = 0;
            for i = -1:1
                for j = -1:1
                    if (x + i > 0 & x + i <= size(P, "r") & y + j > 0 & y + j <= size(P, "c"))
                        r = r + K(i+2, j+2) * P(x+i, y+j)
                    end
                end
            end
            N(x,y) = r;
        end
    end
endfunction

function N = convol2d(K, P)
    N = zeros(P);
    P2 = [zeros(1,size(P,"c")); P; zeros(1,size(P,"c"))]
    P2 = [zeros(size(P2, "r"), 1), P2, zeros(size(P2, "r"), 1)]
    for x=2:(size(P2, "r") - 1)
        for y=2:(size(P2, "c") - 1)
            r = 0;
            N(x-1,y-1) = sum(K .* P2(x-1:x+1, y-1:y+1));
        end
    end
endfunction

function N = dummy_blurr(P)
    N = zeros(P);
    K = 1/159 * [2  4  5  4 2
    4  9 12  9 4
    5 12 15 12 5
    4  9 12  9 4
    2  4  5  4 2];
    for x=1:size(P, "r")
        for y=1:size(P, "c")
            r = 0;
            for i = -2:2
                for j = -2:2
                    if (x + i > 0 & x + i <= size(P, "r") & y + j > 0 & y + j <= size(P, "c"))
                        r = r + K(i+3, j+3) * P(x+i, y+j)
                    end
                end
            end
            N(x,y) = r;
        end
    end
endfunction

function N = blurr(P)
    N = zeros(P);
    P2 = [zeros(1,size(P,"c")); zeros(1,size(P,"c")); P; zeros(1,size(P,"c")); zeros(1,size(P,"c"))]
    P2 = [zeros(size(P2, "r"), 1), zeros(size(P2, "r"), 1), P2, zeros(size(P2, "r"), 1), zeros(size(P2, "r"), 1)]
    K = 1/159 * [2  4  5  4 2
    4  9 12  9 4
    5 12 15 12 5
    4  9 12  9 4
    2  4  5  4 2];

    for x=3:(size(P2, "r") - 2)
        for y=3:(size(P2, "c") - 2)
            r = 0;
            N(x-2,y-2) = sum(K .* P2(x-2:x+2, y-2:y+2));
        end
    end
endfunction

function []=showImage(imageMatrix)
    // Normalisation de l'image
    maxImage = max(imageMatrix);
    imageNormalisee = floor((imageMatrix ./ maxImage) .* 255);
    // Affichage
    f = figure();
    f.color_map = graycolormap(255);
    f.background = -2;
    Matplot(imageNormalisee);
    a = gca();
    a.axes_visible = ["off" "off" "off"]
endfunction


function edgeDemo(filename)
    // Ouvre l'image PGM filename
    disp("Ouverture "+filename+" ...");
    gray_m = readimage(filename);

    // Affiche l'image en niveaux de gris.
    disp("Affiche l''image en niveaux de gris...");
    showImage(gray_m);

    // Lissage
    disp("Lissage de l''image ...");
    blurred_m = blurr(gray_m);
    showImage(blurred_m);


    // Sobel
    disp("Application de la détection de contour Canny (Sobel) sur "+filename+"...");
    GX = convol2d([-1 0 1 ; -2 0 2 ; -1 0 1], blurred_m);
    GY = convol2d([-1 -2 -1 ; 0 0 0 ; 1 2 1], blurred_m);
    G = sqrt(GX.^2+GY.^2);

    showImage(G);

    // Prewitt
    disp("Application de la détection de contour Canny (Prewitt) sur "+filename+"...");
    GX = convol2d([-1 0 1 ; -1 0 1 ; -1 0 1], blurred_m);
    GY = convol2d([-1 -1 -1 ; 0 0 0 ; 1 1 1], blurred_m);
    G = sqrt(GX.^2+GY.^2);

    showImage(G)

    // Scharr
    disp("Application de la détection de contour Canny (Scharr) sur "+filename+"...");
    GX = convol2d([3 0 -3 ; 10 0 -10 ; 3 0 -3], blurred_m);
    GY = convol2d([3 10 3 ; 0 0 0 ; -3 -10 -3], blurred_m);
    G = sqrt(GX.^2+GY.^2);

    showImage(G)
endfunction

edgeDemo("Scilab.pgm")
