Listing Matlab

[ Précédente ] [ Accueil ]

Voici le listing du programme écrit sous Matlab. Cette version permet de traiter des images couleurs à trois canaux. Il reçoit une image couleur en chargeant séparement ses trois composantes à partir de trois images à échelle de niveaux de gris, chacune représentant un canal spectral. Il s'agit de rentrer les noms des trois images, enregistrées au format .ima de l'ENST, dans les trois champs prévus à cet effet au début du script. Le fichier .m est également disponible dans le répertoire courant.

% PROJET MBAI

% Analyse en composantes principales / Transformation de Karhunen-Loeve

% Joel Budin

% Xavier Collignon

% ENST, 2001

close all, clear all, clc;

% Chargement des images / normalisation entre 0 et 255 / Affichage

I1 = ima2mat('billcat_red');

I2 = ima2mat('billcat_green');

I3 = ima2mat('billcat_blue');

%[I1, I2, I3] = three_little_squares;

n = size(I1,1);

m = size(I1,2);

X1 = I1(:);

X2 = I2(:);

X3 = I3(:);

maxi = [max(max(I1)) max(max(I2)) max(max(I3))];

maxi = max(maxi);

mini = [min(min(I1)) min(min(I2)) min(min(I3))];

mini = min(mini);

I1 = 255*(I1-mini)/(maxi-mini);

I2 = 255*(I2-mini)/(maxi-mini);

I3 = 255*(I3-mini)/(maxi-mini);

figure(1)

subplot(1,3,1)

colormap gray(256);

image(I1);

title('Canal 1');

subplot(1,3,2)

colormap gray(256);

image(I2);

title('Canal 2');

subplot(1,3,3)

colormap gray(256);

image(I3);

title('Canal 3');

% Moyennes, ecarts-type, centrees, pour TKL et reduites pour ACP

l = length(X1);

m1 = mean(X1);

m2 = mean(X2);

m3 = mean(X3);

s1 = std(X1);

s2 = std(X2);

s3 = std(X3);

X1_kl = (X1-m1);

X2_kl = (X2-m2);

X3_kl = (X3-m3);

X1_acp = (X1-m1)/s1;

X2_acp = (X2-m2)/s2;

X3_acp = (X3-m3)/s3;

% Covariance, vecteurs propres et valeurs propres

X_kl = zeros(l,3);

X_kl(:,1) = X1_kl;

X_kl(:,2) = X2_kl;

X_kl(:,3) = X3_kl;

Co_kl = cov(X_kl);

[Vp_kl,Xp_kl] = eig(Co_kl);

E_kl = eig(Co_kl);

P_kl = sum(E_kl);

[mini i3_kl] = min([E_kl(1) E_kl(2) E_kl(3)]); % Classement par ordre decroissant des val.propres

[maxi i1_kl] = max([E_kl(1) E_kl(2) E_kl(3)]); % | -> i1_kl, i2_kl, i3_kl

if (i1_kl+i3_kl==3) i2_kl = 3; % |

elseif (i1_kl+i3_kl==4) i2_kl = 2; % |

elseif (i1_kl+i3_kl==5) i2_kl = 1; % |

end % |

P1_kl = E_kl(i1_kl)/P_kl;

P2_kl = E_kl(i2_kl)/P_kl;

P3_kl = E_kl(i3_kl)/P_kl;

X_acp = zeros(l,3);

X_acp(:,1) = X1_acp;

X_acp(:,2) = X2_acp;

X_acp(:,3) = X3_acp;

Co_acp = cov(X_acp);

[Vp_acp,Xp_acp] = eig(Co_acp);

E_acp = eig(Co_acp);

P_acp = sum(E_acp);

[mini i3_acp] = min([E_acp(1) E_acp(2) E_acp(3)]);

[maxi i1_acp] = max([E_acp(1) E_acp(2) E_acp(3)]);

if (i1_acp+i3_acp==3) i2_acp = 3;

elseif (i1_acp+i3_acp==4) i2_acp = 2;

elseif (i1_acp+i3_acp==5) i2_acp = 1;

end

P1_acp = E_acp(i1_acp)/P_acp;

P2_acp = E_acp(i2_acp)/P_acp;

P3_acp = E_acp(i3_acp)/P_acp;

% Transformations

tX_kl = X_kl*Vp_kl;

Cov_tX_kl = cov(tX_kl);

tX1_kl = tX_kl(:,i1_kl);

tX2_kl = tX_kl(:,i2_kl);

tX3_kl = tX_kl(:,i3_kl);

tX_acp = X_acp*Vp_acp;

Cov_tX_acp = cov(tX_acp);

tX1_acp = tX_acp(:,i1_acp);

tX2_acp = tX_acp(:,i2_acp);

tX3_acp = tX_acp(:,i3_acp);

% Calculs min, max, re-centrage entre 0 et 255, recreation de matrices images

maxi_kl = [max(tX1_kl) max(tX2_kl) max(tX3_kl)];

maxi_kl = max(maxi_kl);

mini_kl = [min(tX1_kl) min(tX2_kl) min(tX3_kl)];

mini_kl = min(mini_kl);

tX1_kl = 255*(tX1_kl-mini_kl)/(maxi_kl-mini_kl);

tX2_kl = 255*(tX2_kl-mini_kl)/(maxi_kl-mini_kl);

tX3_kl = 255*(tX3_kl-mini_kl)/(maxi_kl-mini_kl);

s1_kl = std(tX1_kl);

s2_kl = std(tX2_kl);

s3_kl = std(tX3_kl);

tX1_kl = reshape(tX1_kl,n,m);

tX2_kl = reshape(tX2_kl,n,m);

tX3_kl = reshape(tX3_kl,n,m);

maxi_acp = [max(tX1_acp) max(tX2_acp) max(tX3_acp)];

maxi_acp = max(maxi_acp);

mini_acp = [min(tX1_acp) min(tX2_acp) min(tX3_acp)];

mini_acp = min(mini_acp);

tX1_acp = 255*(tX1_acp-mini_acp)/(maxi_acp-mini_acp);

tX2_acp = 255*(tX2_acp-mini_acp)/(maxi_acp-mini_acp);

tX3_acp = 255*(tX3_acp-mini_acp)/(maxi_acp-mini_acp);

s1_acp = std(tX1_acp);

s2_acp = std(tX2_acp);

s3_acp = std(tX3_acp);

tX1_acp = reshape(tX1_acp,n,m);

tX2_acp = reshape(tX2_acp,n,m);

tX3_acp = reshape(tX3_acp,n,m);

% Affichage des composantes principales et des pourcentages d'inertie correspondants

figure(2)

subplot(2,3,1)

colormap gray(256);

image(tX1_kl);

title(['KL1 : ', num2str(round(1000*P1_kl)/10), '%']);

subplot(2,3,2)

colormap gray(256);

image(tX2_kl);

title(['KL2 : ', num2str(round(1000*P2_kl)/10), '%']);

subplot(2,3,3)

colormap gray(256);

image(tX3_kl);

title(['KL3 : ', num2str(round(1000*P3_kl)/10), '%']);

subplot(2,3,4)

colormap gray(256);

image(tX1_acp);

title(['ACP1 : ', num2str(round(1000*P1_acp)/10), '%']);

subplot(2,3,5)

colormap gray(256);

image(tX2_acp);

title(['ACP2 : ', num2str(round(1000*P2_acp)/10), '%']);

subplot(2,3,6)

colormap gray(256);

image(tX3_acp);

title(['ACP3 : ', num2str(round(1000*P3_acp)/10), '%']);

% Affichage des qualités de représentation sur le premier axe

Rpred_kl_axe = ((Vp_kl(1,i1_kl))^2)/((Vp_kl(1,1))^2+(Vp_kl(2,1))^2+(Vp_kl(3,1))^2)

Rpred_acp_axe = ((Vp_acp(1,i1_acp))^2)/((Vp_acp(1,1))^2+(Vp_acp(2,1))^2+(Vp_acp(3,1))^2)

Rpgreen_kl_axe = ((Vp_kl(2,i1_kl))^2)/((Vp_kl(1,1))^2+(Vp_kl(2,1))^2+(Vp_kl(3,1))^2)

Rpgreen_acp_axe = ((Vp_acp(2,i1_acp))^2)/((Vp_acp(1,1))^2+(Vp_acp(2,1))^2+(Vp_acp(3,1))^2)

Rpblue_kl_axe = ((Vp_kl(3,i1_kl))^2)/((Vp_kl(1,1))^2+(Vp_kl(2,1))^2+(Vp_kl(3,1))^2)

Rpblue_acp_axe = ((Vp_acp(3,i1_acp))^2)/((Vp_acp(1,1))^2+(Vp_acp(2,1))^2+(Vp_acp(3,1))^2)

% Affichage des qualités de représentation sur le premier plan

Rpred_kl_plan = Rpred_kl_axe+((Vp_kl(1,i2_kl))^2)/((Vp_kl(1,2))^2+(Vp_kl(2,2))^2+(Vp_kl(3,2))^2)

Rpred_acp_plan = Rpred_acp_axe+((Vp_acp(1,i2_kl))^2)/((Vp_acp(1,2))^2+(Vp_acp(2,2))^2+(Vp_acp(3,2))^2)

Rpgreen_kl_plan = Rpgreen_kl_axe+((Vp_kl(2,i2_kl))^2)/((Vp_kl(1,2))^2+(Vp_kl(2,2))^2+(Vp_kl(3,2))^2)

Rpgreen_acp_plan = Rpgreen_acp_axe+((Vp_acp(2,i2_kl))^2)/((Vp_acp(1,2))^2+(Vp_acp(2,2))^2+(Vp_acp(3,2))^2)

Rpblue_kl_plan = Rpblue_kl_axe+((Vp_kl(3,i2_kl))^2)/((Vp_kl(1,2))^2+(Vp_kl(2,2))^2+(Vp_kl(3,2))^2)

Rpblue_acp_plan = Rpblue_acp_axe+((Vp_kl(3,i2_acp))^2)/((Vp_acp(1,2))^2+(Vp_acp(2,2))^2+(Vp_acp(3,2))^2)

[ Précédente ] [ Accueil ]