%PRACTICA 5. ANALISIS ESTADISTICO DE CLASES

%Para empezar, declaramos una matriz 9x6 (9 clases, 6 bandas) con los valores de la 
%media aritmetica para cada una de las bandas

tiempo=clock;

c=[121.77 75.76 115.78 118.83 118.83 102.95  %clase 1
    102.95 62.32 95.10 109.08 191.84 85.80  %clase 2 ...
    65.69 28.30 24.64 110.19 67.48 18.35
    67.95 30.88 24.64 161.48 85.99 22.59
    108.48 65.64 98.87 103.29 195.23 99.99
    83.60 45.10 63.95 74.36 118.01 52.62
    100.08 62.45 98.40 99.72 199.53 113.31
    89.04 48.44 68.89 69.94 153.83 89.55
    95.34 50.87 67.59 81.77 138.90 74.72];
 

%cargamos las imágenes correspondientes a las diferentes bandas en las que se ha 
%adquirido la imagen por satélite.

imagen1='Tm1.dat';
    fid=fopen(imagen1,'r+b');
    [imagen1,count]=fread(fid,[1024,1024],'uint8');
    fclose(fid);
    
imagen2='Tm2.dat';
    fid=fopen(imagen2,'r+b');
    [imagen2,count]=fread(fid,[1024,1024],'uint8');
    fclose(fid);

imagen3='Tm3.dat';
    fid=fopen(imagen3,'r+b');
    [imagen3,count]=fread(fid,[1024,1024],'uint8');
    fclose(fid); 
    
imagen4='Tm4.dat';
    fid=fopen(imagen4,'r+b');
    [imagen4,count]=fread(fid,[1024,1024],'uint8');
    fclose(fid); 

imagen5='Tm5.dat';
    fid=fopen(imagen5,'r+b');
    [imagen5,count]=fread(fid,[1024,1024],'uint8');
    fclose(fid);

imagen6='Tm7.dat';
    fid=fopen(imagen6,'r+b');
    [imagen6,count]=fread(fid,[1024,1024],'uint8');
    fclose(fid);
    
    
%Dado que la gestión de memoria de Matlab no es excesivamente buena, para evitar problemas
%de desbordamiento, procesaremos la imagen en dos partes de 1024x512, para luego unirlas
%y obtener la imagen definitiva. Para esto, emplearemos bucles anidados que calcularán,
%primero,la distancia de cada píxel a cada una de las medias definidas, es decir, a cada
%una de las clases en las que queremos dividir la imagen.
%Posteriormente, se calculará cual de estas distancias es la menor, para tomar como valor
%final del pixel de la nueva imagen el número de clase más cercano a dicho pixel.

%Este bucle procesa la primera mitad (la izquierda) de la imagen

for(i=1:1024)
    for(j=1:512)
        for(k=1:9)
            for(n=1:6)
               
               im1(k,n)=((imagen1(i,j)-c(k,n))^2);
					im2(k,n)=((imagen2(i,j)-c(k,n))^2);
					im3(k,n)=((imagen3(i,j)-c(k,n))^2);
					im4(k,n)=((imagen4(i,j)-c(k,n))^2);
					im5(k,n)=((imagen5(i,j)-c(k,n))^2);
					im6(k,n)=((imagen6(i,j)-c(k,n))^2);
               
            end
        end
        
        for(n=1:6)
         %Aquí calculamos las distancias de cada píxel a cada clase. 
         d1(i,j)=sqrt(im1(1,n) + im2(1,n) + im3(1,n) + im4(1,n) + im5(1,n) + im6(1,n));
        	d2(i,j)=sqrt(im1(2,n) + im2(2,n) + im3(2,n) + im4(2,n) + im5(2,n) + im6(2,n)); 
        	d3(i,j)=sqrt(im1(3,n) + im2(3,n) + im3(3,n) + im4(3,n) + im5(3,n) + im6(3,n));
        	d4(i,j)=sqrt(im1(4,n) + im2(4,n) + im3(4,n) + im4(4,n) + im5(4,n) + im6(4,n));
        	d5(i,j)=sqrt(im1(5,n) + im2(5,n) + im3(5,n) + im4(5,n) + im5(5,n) + im6(5,n)); 
        	d6(i,j)=sqrt(im1(6,n) + im2(6,n) + im3(6,n) + im4(6,n) + im5(6,n) + im6(6,n));
        	d7(i,j)=sqrt(im1(7,n) + im2(7,n) + im3(7,n) + im4(7,n) + im5(7,n) + im6(7,n));
        	d8(i,j)=sqrt(im1(8,n) + im2(8,n) + im3(8,n) + im4(8,n) + im5(8,n) + im6(8,n)); 
        	d9(i,j)=sqrt(im1(9,n) + im2(9,n) + im3(9,n) + im4(9,n) + im5(9,n) + im6(9,n));
		  end
        
        %Guardamos el valor de mínima distancia de cada píxel a cada clase
        
         dkij=[d1(i,j) d2(i,j) d3(i,j) d4(i,j) d5(i,j) d6(i,j) d7(i,j) d8(i,j) d9(i,j)];
			[Y,dk1(i,j)]=min(dkij);
                    
     end
end

%borramos el workspace para liberar memoria
clear d1 d2 d3 d4 d5 d6 d7 d8 d9 dkij;


%Este bucle procesa la segunda mitad (la derecha) de la imagen de forma análoga al anterior
for(i=1:1024)
    for(j=513:1024)
        for(k=1:9)
            for(n=1:6)
               
               im1(k,n)=((imagen1(i,j)-c(k,n))^2);
					im2(k,n)=((imagen2(i,j)-c(k,n))^2);
					im3(k,n)=((imagen3(i,j)-c(k,n))^2);
					im4(k,n)=((imagen4(i,j)-c(k,n))^2);
					im5(k,n)=((imagen5(i,j)-c(k,n))^2);
					im6(k,n)=((imagen6(i,j)-c(k,n))^2);
               
            end
        end
        
        for(n=1:6)
        	d1(i,j-512)=sqrt(im1(1,n) + im2(1,n) + im3(1,n) + im4(1,n) + im5(1,n) + im6(1,n));
        	d2(i,j-512)=sqrt(im1(2,n) + im2(2,n) + im3(2,n) + im4(2,n) + im5(2,n) + im6(2,n)); 
        	d3(i,j-512)=sqrt(im1(3,n) + im2(3,n) + im3(3,n) + im4(3,n) + im5(3,n) + im6(3,n));
        	d4(i,j-512)=sqrt(im1(4,n) + im2(4,n) + im3(4,n) + im4(4,n) + im5(4,n) + im6(4,n));
        	d5(i,j-512)=sqrt(im1(5,n) + im2(5,n) + im3(5,n) + im4(5,n) + im5(5,n) + im6(5,n)); 
        	d6(i,j-512)=sqrt(im1(6,n) + im2(6,n) + im3(6,n) + im4(6,n) + im5(6,n) + im6(6,n));
        	d7(i,j-512)=sqrt(im1(7,n) + im2(7,n) + im3(7,n) + im4(7,n) + im5(7,n) + im6(7,n));
        	d8(i,j-512)=sqrt(im1(8,n) + im2(8,n) + im3(8,n) + im4(8,n) + im5(8,n) + im6(8,n)); 
        	d9(i,j-512)=sqrt(im1(9,n) + im2(9,n) + im3(9,n) + im4(9,n) + im5(9,n) + im6(9,n));
		  end
        
         dkij=[d1(i,j-512) d2(i,j-512) d3(i,j-512) d4(i,j-512) d5(i,j-512) d6(i,j-512) d7(i,j-512) d8(i,j-512) d9(i,j-512)];
			[Y,dk2(i,j-512)]=min(dkij);
                    
     end
end


%Para acabar, componemos la imagen final a partir de las dos mitades obtenidas y la
%guardamos en un fichero para el posterior procesado con un programa de edición de 
%imágenes, de tal forma que asignaremos a cada clase un color para observar de manera
%rápida y vistosa las diferentes clases en las que hemos dividido la imagen.

dk=[dk1 dk2];

fid=fopen('imagen.dat','w+b');
    count=fwrite(fid,dk,'uint8');
    fclose(fid);
    
tiempo=etime(clock,tiempo)

%FIN