% Calcul du veritable mean square deviation
% En fonction du seuil
function matrice_msd = msd(l, c, N, k)

% On appelle les fonctions handles

ed = @edginess;
cons = @byaverage;
n = @norme;
moy = @moyenne;
bds = @boundaries;
% 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]);

for i = 1:l-1
    for j = 1:c-1
        e(i,j) = feval( ed, N, i, j);
        d(i,j) = feval( cons, N ,i, j);
        R(i,j) = e(i,j)*d(i,j);
        if R(i,j) > k
            B(i,j) = 1;
            msd(i,j) = 0;
            M(i,j) = 1;
        else B(i,j) = N(i,j);    
        end
    end     
end

% on va limiter l'etude de la matrice temoin

M(l,:) = 1;
M(:,c) = 1;
M(1,:) = 1;
M(:,1) = 1;
t = 0;
% Ensuite on va fill les holes
couple =[[-1 0];[1 0];[0 -1];[0 1]];
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;                
            end
            % Et on donne la valeur a tous les points concernes
            for e = 1:s
                matrice_msd(AM(e,1),AM(e,2)) = w;                
            end
        else continue    
        end
    end
end  
%Images_contours = feval(bds, l, c, N, k)