% Calcul du veritable mean square deviation
% En fonction du seuil mais aussi calcul de la piece wise constant pour ne
% pas perdre de temps a la recalculer apres
function [msd,MC,C] = msd_relier2(l, c, N, R, k)

% On appelle les fonctions handles

ed = @edginess;
cons = @byaverage;
n = @norme;
moy = @moyenne;
bds = @boundaries;
link = @relier;

% Voila la premiere boucle pour trouver
% les edges et la matrice qui va avec
% En plus on cree la matrice temoin M

B = zeros(size(N));
B(l,:) = N(l,:);
B(:,c) = N(:,c);
M = zeros([l c]);
% M est la matrice temoin, mais on commence par enlever 
% tous les points qui sont des edges
% B est la matrice couleur qui a aussi les edges en blanc
    [M, B] = feval(bds, l, c, N, R, k,'no');
% Maintenant on applique a M RELIER, pour
% ameliorer l'etude
    [M, B]= feval(link, M, B, l, c,'no');
    %[M, B]= feval(link, M, B, l, c,'no');
    %[M, B]= feval(link, M, B, l, c,'no');
    % On sauvegarde la matrice contour apres link
    MC = M;
% on va limiter l'etude a la matrice temoin
M(l,:) = 1;
M(:,c) = 1;
M(1,:) = 1;
M(:,1) = 1;
t = 0;
% Ensuite on calcule la reelle mean square 
% deviation al'aide du "fill" programme
couple =[[-1 0];[1 0];[0 -1];[0 1]];
C = B;
for i = 1:l
    for j = 1:c
        if M(i,j) ~= 1
            M(i,j) = 1;
            % On cree la matrice AP qui represente 
            % les coordonnees des derniers points crees         
            % initialisation
            AP = [i j];
            [m,n] = size(AP);
            % On va essayer de tout mettre a la moyenne sur tous les points
            s = 0;
            So = 0;
            % d va representer le nombre de lignes dans l'ancienne matrice A
            while m > 0
                AV = AP;
                AP = [];                   
                % u va representer le nombre de lignes dans la nouvelle matrice A
                u = 0;
                for k = 1:m
                    for r = 1:4           
                        if M(AV(k,1) + couple(r,1), AV(k,2) + couple(r,2)) ~= 1
                            M(AV(k,1) + couple(r,1), AV(k,2) + couple(r,2)) = 1;
                            % Ici on agrandit l'after matrice
                            u = u + 1;
                            % Ici on agrandit le vecteur contenant les coordonnes des points
                            % dont on veut la moyenne
                            s = s + 1;
                            % On transforme le vecteur avant en vecteur apres
                            AP(u,1) = AV(k,1) + couple(r,1);
                            AP(u,2) = AV(k,2) + couple(r,2);
                            % Et on ajoute un vecteur dans le vecteur moyenne
                            AM(s,1) = AV(k,1) + couple(r,1);
                            AM(s,2) = AV(k,2) + couple(r,2);
                            So = So + B(AV(k,1) + couple(r,1), AV(k,2) + couple(r,2));
                        end
                    end
                end
                [m n] = size(AP);
            end
            % On calcule le mean square deviation
            w = 0;
            for e = 1:s
                w = (B(AM(e,1),AM(e,2)) - (So)/s)^2 + w;
                C(AM(e,1),AM(e,2)) = (So)/s;
            end
            % Et on donne la valeur a tous les points concernes
            for e = 1:s
                msd(AM(e,1),AM(e,2)) = w;                
            end
        else continue    
        end
    end
end