function [val] = CritAtt(Imseg,Imori)
% critère de dispersion des attributs de Levine et  Nazif
% function [val] = CritAtt(Imseg,Imori)
% Imseg : matrice correspondant à l'image segmentée  avec des valeurs allant de 1 à n pour les n régions segmentées 
% Imori : matrice correspondant à l'image originale
%
% Utilitaires pour l'évaluation de la segmentation d'images
% Toolbox matlab (version 5.3)
%
% (c) Laboratoire de Vision et Robotique (UPRES EA 2078)
% ENSI de Bourges - Université d'Orléans
%
% Sébastien Chabrier : sebastien.chabrier@ensi-bourges.fr
%
% Si vous utilisez cette toolbox, veuillez citer ce papier svp.
%
%S. Chabrier, B. Emile, C. Rosenberger, H. Laurent, 
%"Unsupervised performance evaluation of image segmentation", 
%Special Issue on Performance Evaluation in Image Processing, 
%EURASIP Journal on Applied Signal Processing, pages 1-12, 2006.


disp('Critère dispersion des attributs')

if (min(min(Imseg))==0)
   Imseg=double(Imseg)+1;
end;


NBclass=double(max(max(Imseg)));

Imseg=double(Imseg);
Imori=double(Imori);

nbre2lig=size(Imori,1);
nbre2col=size(Imori,2);

taille=3;
m1=1;
m2=1;
m3=1;
m4=1;

% statistiques sur l'image originale

deb=(taille-1)/2+1;
milieu=deb-1;
cpt=ones(NBclass)*0;

for i=1:(nbre2lig-taille+1)
   for j=1:(nbre2col-taille+1)
      
      cpt(Imseg(i+milieu,j+milieu))=cpt(Imseg(i+milieu,j+milieu))+1;
      cpt2=0;
      
% calcul de la moyenne     
	aux=0;
   cpt3=0; 
   for l=i:i+taille-1
      for k=j:j+taille-1
         if (Imseg(l,k)==Imseg(i+milieu,j+milieu))
            aux=aux+Imori(l,k);
            cpt3=cpt3+1;
         end;
      end;
   end;
   moyregion=aux/cpt3;
   
   if (m1==1)
      cpt2=cpt2+1;
      param(Imseg(i+milieu,j+milieu),cpt(Imseg(i+milieu,j+milieu)),cpt2)=moyregion;
   end;   
         
% calcul de la variance
      
      if (m2==1)
         cpt2=cpt2+1;
         aux=0;
         cpt3=0;
         for l=i:i+taille-1
            for k=j:j+taille-1
               if (Imseg(l,k)==Imseg(i+milieu,j+milieu))
                  aux2=(Imori(l,k)-moyregion)*2/255;
                  aux=aux+aux2.^2;
                  cpt3=cpt3+1;
               end;
            end;
         end;
         param(Imseg(i+milieu,j+milieu),cpt(Imseg(i+milieu,j+milieu)),cpt2)=aux/cpt3;
      end;
      
% calcul du moment d'ordre 3
      
      if (m3==1)
         cpt2=cpt2+1;
         aux=0;
         cpt3=0;
         for l=i:i+taille-1
            for k=j:j+taille-1
               if (Imseg(l,k)==Imseg(i+milieu,j+milieu))
                  aux2=(Imori(l,k)-moyregion)*2/255;
                  aux=aux+aux2^3;
                  cpt3=cpt3+1;
               end;
            end;
         end;
         param(Imseg(i+milieu,j+milieu),cpt(Imseg(i+milieu,j+milieu)),cpt2)=aux/cpt3;
      end;
      
% calcul du moment d'ordre 4
      
      if (m4==1)
         cpt2=cpt2+1;
         aux=0;
         cpt3=0;
         for l=i:i+taille-1
            for k=j:j+taille-1
               if (Imseg(l,k)==Imseg(i+milieu,j+milieu))
                  aux2=(Imori(l,k)-moyregion)*2/255;
                  aux=aux+aux2^4;
                  cpt3=cpt3+1;
               end;
            end;
         end;
         param(Imseg(i+milieu,j+milieu),cpt(Imseg(i+milieu,j+milieu)),cpt2)=aux/cpt3;
      end;
          
   end;
 end;
 
 % Normalisation et réduction 
 
 l=1;
 for i=1:NBclass
    for j=1:cpt(i)
%       affich=strcat('BOUCLE 2 => i: ',int2str(i),' | j: ',int2str(j),' | l: ',int2str(l),' | k: ',int2str(k));
%       disp(affich);
       param2(l,:)=param(i,j,:);
       l=l+1;
    end;
 end;
  
 barycentre=mean(param2);
 variance=ones(size(param2,2))*0;
    
 for l=1:size(param2,1)
    for k=1:size(param2,2)
%       affich=strcat('BOUCLE 3 => i: ',int2str(i),' | j: ',int2str(j),' | l: ',int2str(l),' | k: ',int2str(k));
%       disp(affich);
       variance(k)=variance(k)+(param2(l,k)-barycentre(1,k))^2;   
    end;
 end;
 
 variance=variance/size(param2,2);
    
 for l=1:size(param2,1)
    for k=1:size(param2,2)   
%       affich=strcat('BOUCLE 1 => i: ',int2str(i),' | j: ',int2str(j),' | l: ',int2str(l),' | k: ',int2str(k));
%       disp(affich);
       param2(l,k)=(param2(l,k)-barycentre(1,k))/variance(k);
    end;
 end;
 
 val=0;
 %val2=0;
 
 % Calcul de la mesure de dispersion
 
 dim=1;
 
 for i=1:NBclass
    
    for k=dim:dim+cpt(i)-1
       for l=1:size(param2,2)
%          affich=strcat('BOUCLE 4 1 => i: ',int2str(i),' | j: ',int2str(j),' | l: ',int2str(l),' | k: ',int2str(k));
%          disp(affich);
          tab(k,l)=param2(k,l);
       end;
       dim=dim+1;
    end;
    
    barycentre=mean(tab);
%    ecart=std(tab)*4/(255*255);
%    val2=val2+sum(ecart)/size(param2,2);
    
    aux=0;
    for j=1:cpt(i)
%       affich=strcat('BOUCLE 4 2 => i: ',int2str(i),' | j: ',int2str(j),' | l: ',int2str(l),' | k: ',int2str(k));
%       disp(affich);
       aux=sqrt(sum((barycentre(1,:)-tab(j,:)).^2))+aux;
    end;
    aux=aux/cpt(i);
    val=val+aux;
 end;
 
 val=val/NBclass;
% val2=val2/NBclass;

 disp(val)
% disp(val2)