init

TP1/calcul_parametres.m

@@ -0,0 +1,16 @@
+## Author: Laurent <laurent@pop-os>
+## Created: 2020-11-02
+
+function [r, a, b] = calcul_parametres (X, Y)
+
+  moyX = mean(X)
+  moyY = mean(Y)
+  sigmaX = sqrt( (X - moyX).'*(X - moyX) / size(X)(1) )
+  sigmaY = sqrt( (Y - moyY).'*(Y - moyY) / size(Y)(1) )
+  covariance = (X - moyX).'*(Y - moyY) / size(X)(1)
+  
+  r = covariance/(sigmaX*sigmaY)
+  a = covariance/(sigmaX^2)
+  b = moyY - moyX*a
+
+endfunction

TP1/ecriture_RVB.m

@@ -0,0 +1,15 @@
+## Author: Laurent <laurent@pop-os>
+## Created: 2020-11-02
+
+function image_RVB = ecriture_RVB(image_originale)
+
+  [H, L] = size(image_originale);
+
+  image_RVB = zeros(H/2, L/2);
+  
+  image_RVB(:,:,1) = image_originale( 1:2:end, 2:2:end ); ## rouge
+  image_RVB(:,:,3) = image_originale( 2:2:end, 1:2:end ); ## bleu
+  image_RVB(:,:,2) = (  image_originale( 1:2:end, 1:2:end )
+                      + image_originale( 2:2:end, 2:2:end ) )/2; ## vert
+
+endfunction

TP1/exercice_1.m

@@ -0,0 +1,49 @@
+clear;
+close all;
+taille_ecran = get(0,'ScreenSize');
+L = taille_ecran(3);
+H = taille_ecran(4);
+
+% Lecture d'une image interne a Matlab et conversion en doubles :
+I_max = 255;
+I_min = 0;
+I = rgb2gray(imread('autumn.tif'));
+I = double(I);
+
+% Affichage de l'image :
+figure('Name','Exemple d''image interne a Matlab','Position',[0,0,0.33*L,0.3*L]);
+imagesc(I);
+axis off;
+axis equal;
+colormap gray;
+
+% Affichage de l'histogramme de l'image :
+histogramme = hist(I(:),[I_min:I_max]);
+figure('Name','Histogramme des niveaux de gris','Position',[0.33*L,0,0.33*L,0.3*L]);
+plot([I_min:I_max],histogramme);
+axis([I_min I_max 0 1.1*max(histogramme)]);
+xlabel('Niveau de gris','FontSize',15);
+ylabel('Frequence','FontSize',15);
+set(gca,'FontSize',20);
+
+% Calcul de vecteurs contenant les niveaux de gris de gauche (de droite) d'une paire de pixels :
+[X,Y] = vectorisation(I);
+
+% Affichage des paires de niveaux de gris sous la forme d'un nuage de points :
+figure('Name','Mise en evidence de la correlation entre pixels voisins','Position',[0.66*L,0,0.33*L,0.3*L]);
+plot(X,Y,'b+');
+axis([I_min I_max I_min I_max]);
+xlabel('Niveau de gris du pixel de gauche','FontSize',15);
+ylabel('Niveau de gris du pixel de droite','FontSize',15);
+hold on;
+set(gca,'FontSize',20);
+axis equal;
+
+% Calcul des parametres de la droite de regression :
+[r,a,b] = calcul_parametres(X,Y);
+fprintf('Coefficient de correlation lineaire : %.4f\n',r);
+
+% Affichage de la droite de regression (d'equation y = a*x+b) :
+x = I_min:I_max;
+y = a*x+b;
+plot(x,y,'r','LineWidth',3);

TP1/exercice_2.m

@@ -0,0 +1,22 @@
+clear;
+close all;
+taille_ecran = get(0,'ScreenSize');
+L = taille_ecran(3);
+H = taille_ecran(4);
+
+% Lecture d'une image interne a Matlab et conversion en doubles :
+I_max = 255;
+I_min = 0;
+I = rgb2gray(imread('autumn.tif'));
+I = double(I);
+
+% décorrélation
+I(:,2:end) = I(:,2:end) - I(:,1:end-1);
+
+% Affichage de l'image :
+figure('Name','Exemple d''image interne a Matlab','Position',[0,0,0.33*L,0.3*L]);
+imagesc(I);
+axis off;
+axis equal;
+colormap gray;
+

TP1/exercice_Matlab.m

@@ -0,0 +1,25 @@
+clear;
+close all;
+taille_ecran = get(0,'ScreenSize');
+L = taille_ecran(3);
+H = taille_ecran(4);
+
+% Lecture et affichage de l'image originale :
+image_originale = imread('image_originale.png');
+figure('Name','Image originale','Position',[0,0,0.33*L,0.3*L]);
+imagesc(image_originale);
+axis off;
+axis equal;
+colormap gray;
+
+% Conversion des entiers non signés en doubles :
+image_originale = double(image_originale);
+
+% Ecriture de image_RVB :
+image_RVB = ecriture_RVB(image_originale);
+
+% Affichage de l'image RVB convertie en entiers non signes :
+figure('Name','Image originale convertie en RVB','Position',[[0.33*L,0,0.33*L,0.3*L]]);
+imagesc(uint8(image_RVB));
+axis off;
+axis equal;

TP1/vectorisation.m

@@ -0,0 +1,12 @@
+## Author: Laurent <laurent@pop-os>
+## Created: 2020-11-02
+
+function [X,Y] = vectorisation(I)
+
+  I_g = I(:, 1:end-1);
+  I_d = I(:, 2:end  );
+  
+  X = reshape(I_g,[],1);
+  Y = reshape(I_d,[],1);
+  
+endfunction

TP2/calcul_frequences.m

@@ -0,0 +1,14 @@
+## Author: Laurent
+## Created: 2020-11-10
+
+function freq = calcul_frequences (texte, alphabet)
+  
+  freq = [];
+
+  for lettre = alphabet.'
+    freq = [freq; length(find(texte==lettre))/length(texte)];
+  endfor
+  
+  freq = freq.';
+  
+endfunction

TP2/coeff_compression_image.m

@@ -0,0 +1,7 @@
+## Created: 2020-11-10
+
+function coeff_compression = coeff_compression_image(histogramme,dico)
+
+  
+
+endfunction

TP2/coeff_compression_texte.m

@@ -0,0 +1,7 @@
+## Created: 2020-11-10
+
+function coeff_compression = coeff_compression_texte(texte, texte_encode)
+
+  coeff_compression = length(texte)*8/length(texte_encode);
+
+endfunction

TP2/exercice_1.m

@@ -0,0 +1,24 @@
+clear;
+close all;
+load exercice_Matlab;
+
+% Encodage du texte :
+alphabet_cell = num2cell(selection_alphabet);			% Format 'cell' requis par la fonction huffmandict
+dico = huffmandict(alphabet_cell,selection_frequences.')		% dico{j,2} = code binaire du caractere ASCII dico{j,1}
+texte_encode = huffmanenco(num2cell(texte), dico);
+fprintf('Codage binaire du texte : %s\n',strrep(num2str(texte_encode),' ',''));
+
+% Calcul du coefficient de compression obtenu par le codage de Huffman :
+coeff_compression = coeff_compression_texte(texte,texte_encode);
+fprintf('Coefficient de compression : %.4f\n',coeff_compression);
+
+% Decodage du texte : 
+texte_decode = huffmandeco(texte_encode,dico);
+
+% Verification de l'integrite du texte decode :
+integrite_decodage = ~strcmp(texte,texte_decode);
+if integrite_decodage
+	fprintf('Integrite du decodage OK\n');
+else
+	fprintf('Probleme d''integrite du decodage !\n');
+end

TP2/exercice_2.m

@@ -0,0 +1,28 @@
+clear;
+close all;
+
+% Codage de Huffman de l'image initiale :
+load exercice_1_TP1;
+frequences = histogramme/sum(histogramme);
+dico = huffmandict(I_min:I_max,frequences);
+
+% Calcul du coefficient de compression obtenu par le codage de Huffman :
+coeff_compression_avant_decorrelation = coeff_compression_image(histogramme,dico);
+fprintf('Coefficient de compression avant decorrelation : %.4f\n',coeff_compression_avant_decorrelation);
+
+## 1.13
+
+% Codage de Huffman de l'image decorrelee :
+load exercice_2_TP1;
+set(0,'RecursionLimit',550);	% Plus de 500 appels recursifs (nombre d'entiers a coder : 2*I_max+1 = 511 > 500)
+frequences = histogramme/sum(histogramme);
+dico = huffmandict(I_min:I_max,frequences);
+
+% Calcul du coefficient de compression obtenu par decorrelation prealable au codage de Huffman :
+coeff_compression_apres_decorrelation = coeff_compression_image(histogramme,dico);
+fprintf('Coefficient de compression apres decorrelation : %.4f\n',coeff_compression_apres_decorrelation);
+
+## 1.57
+
+% Calcul du gain en compression :
+fprintf('Gain en compression : %.4f\n',coeff_compression_apres_decorrelation/coeff_compression_avant_decorrelation);

TP2/exercice_3.m

@@ -0,0 +1,54 @@
+clear;
+close all;
+
+% Definition d'un alphabet ASCII avec les lettres A..Z, a..z, les chiffres 1..9 et la ponctuation :
+taille_alphabet = 128;
+codes_ASCII = transpose(0:taille_alphabet-1);		% Vecteur (colonne) de codes ASCII
+alphabet = char(codes_ASCII);				% Vecteur de caracteres ASCII
+
+% Texte a encoder :
+texte = 'Taratatatsointsoin';
+
+% Affichage du texte :
+fprintf('Texte : %s\n',texte);
+
+% Frequences d'apparition des lettres dans le texte :
+frequences = calcul_frequences(texte,alphabet);
+
+% Selection des caracteres de frequences non nulles :
+[selection_frequences,selection_alphabet] = selection(frequences,alphabet);
+
+% Partitionnement de l'intervalle [0,1] :
+bornes = partitionnement(selection_frequences);
+
+% Codage arithmetique :
+[borne_inf,borne_sup] = codage_arithmetique(texte,selection_alphabet,bornes);
+
+% Choix (arbitraire) d'un nombre compris dans l'intervalle [borne_inf,borne_sup] pour encoder le texte :
+k = 1;
+while (floor(k*borne_inf)==floor(k*borne_sup))
+	k = 10*k;	
+end
+texte_encode = floor(k*borne_sup)/k;
+fprintf('Texte encode : %.15f\n',texte_encode);
+
+% Decodage du texte encode :
+texte_decode = '';
+texte_encode_courant = texte_encode;
+for i = 1:length(texte)
+	c = find((texte_encode_courant>bornes(1,:))&(texte_encode_courant<bornes(2,:)));	% Identification du caractere
+	texte_decode = [texte_decode selection_alphabet(c)];				% Concatenation du caractere
+	texte_encode_courant = (texte_encode_courant-bornes(1,c))/selection_frequences(c);	% Soustraction de la borne inferieure et division par la probabilite pour obtenir le caractere suivant
+end
+fprintf('Texte decode : %s\n',texte_decode);
+
+% Calcul du nombre de bits necessaires pour encoder la partie decimale :
+dec_texte_encode = texte_encode*k;					% Partie decimale
+nb_bits_codage_arithmetique = length(dec2bin(dec_texte_encode));	% Conversion de la partie decimale en binaire (dec2bin)
+fprintf('Nombre de bits du codage arithmetique : %d\n',nb_bits_codage_arithmetique);
+
+% Comparaison avec le codage de Huffman :
+dico = huffmandict(num2cell(selection_alphabet),selection_frequences);
+texte_encode = huffmanenco(texte,dico);
+nb_bits_codage_Huffman = length(texte_encode);
+fprintf('Nombre de bits du codage de Huffman : %d\n',nb_bits_codage_Huffman);

TP2/exercice_Matlab.m

@@ -0,0 +1,36 @@
+clear;
+close all;
+taille_ecran = get(0,'ScreenSize');
+L = taille_ecran(3);
+H = taille_ecran(4);
+
+% Definition d'un alphabet ASCII avec les lettres A..Z, a..z, les chiffres 1..9 et la ponctuation :
+taille_alphabet = 128;
+codes_ASCII = transpose(0:taille_alphabet-1);		% Vecteur (colonne) de codes ASCII
+alphabet = char(codes_ASCII);				% Vecteur de caracteres ASCII
+
+% Texte ecrit en francais :
+texte = 'Ceci est un texte en Francais. Il contient des lettres (majuscules et minuscules), des chiffres de zero (0) a neuf (9) et quelques caracteres de ponctuation comme la virgule (,), le point (.), les deux points (:) etc. Par contre il ne contient pas de caractere accentue, de guillemets, dapostrophe ou de cedille, qui ne sont pas des caracteres ASCII standards. Il sert simplement a illustrer la loi de probabilite de lalphabet francais. En fait nimporte quel texte convient, du moment quil soit assez long pour illustrer la langue francaise. Du coup on met bout a bout des phrases sans reel lien logique, mais qui permettront dobtenir un histogramme assez representatif de la langue francaise. Normalement le caractere le plus present est lespace, suivi du e, du s, a, i, t... Ce genre dhistogramme est tres utile en cryptographie : chaque langue possede son propre histogramme, ce qui permet de reconnaitre la langue dun texte code meme si les lettres ont ete changees. Evidemment, un bon outil de cryptographie fera en sorte de ne pas chiffrer deux occurrences dune meme lettre de la meme facon, sinon ce serait trop simple ! Pour finir ce texte, jajoute un k et un w, parce quil ny en avait pas encore.';
+
+% Affichage du texte :
+fprintf('Texte : %s\n',texte);
+
+% Frequences relatives des lettres dans le texte :
+frequences = calcul_frequences(texte,alphabet);
+
+% Selection des caracteres de frequences non nulles :
+[selection_frequences,selection_alphabet] = selection(frequences,alphabet);
+
+% Tri des caracteres par frequences decroissantes :
+[frequences_triees,indices_frequences_triees] = tri(selection_frequences);
+
+% Histogramme des frequences des differents caracteres :
+figure('Name','Histogramme des caracteres du texte','Position',[0.33*L,0,0.67*L,H]);
+bar(frequences_triees);
+set(gca,'XTick',[1:length(alphabet)]);
+set(gca,'XTickLabel',selection_alphabet(indices_frequences_triees));
+xlabel('Alphabet','FontSize',15);
+ylabel('Frequence','FontSize',15);
+set(gca,'FontSize',20);
+
+save exercice_Matlab;

TP2/selection.m

@@ -0,0 +1,10 @@
+## Author: Laurent
+## Created: 2020-11-10
+
+function [selection_frequences, selection_alphabet] = selection(frequences, alphabet)
+  
+  indices = find(frequences>0.0).';
+  selection_frequences = frequences(indices).';
+  selection_alphabet = alphabet(indices);
+  
+endfunction

TP2/test_calcul_frequences.m

@@ -0,0 +1,7 @@
+frequences = calcul_frequences('aabc',transpose('abc'))
+resultat = (frequences==[0.5 ; 0.25 ; 0.25]);
+if resultat~=1
+	fprintf('Il y a un probleme !\n');
+else
+	fprintf('La fonction calcul_frequences est correcte !\n');
+end

TP2/tri.m

@@ -0,0 +1,7 @@
+## Created: 2020-11-10
+
+function [frequences_triees, indices_frequences_triees] = tri(selection_frequences)
+
+  [frequences_triees, indices_frequences_triees] = sort(selection_frequences, 'descend');
+
+endfunction

TP3/calcul_proba.m

@@ -0,0 +1,5 @@
+function [x_min,x_max,probabilite] = calcul_proba(E_nouveau_repere,p);
+
+    %todo
+
+end

TP3/carres_noirs.m

@@ -0,0 +1,72 @@
+clear;
+close all;
+taille_ecran = get(0,'ScreenSize');
+L = taille_ecran(3);
+H = taille_ecran(4);
+
+% Affichage d'images aleatoires contenant au moins un carre noir :
+figure('Name','Images binaires aleatoires','Position',[0.3*L,0,0.7*L,H]);
+numero_image = 1;
+T = 64;
+p = 0.5;
+liste_valeurs_t = 2:5;
+for t = liste_valeurs_t
+
+	% Tirage d'images aleatoires tant qu'aucun carre noir n'a ete detecte :
+	aucun_carre_noir_detecte = 1;
+	cpt_tirages = 0;
+	while aucun_carre_noir_detecte
+
+		% Tirage aleatoire d'une image :
+		I = rand(T);
+		I = I>p;
+
+		% Detection des carres de cote t entierement noirs :
+		J = zeros(T-t+1,T-t+1,t*t);
+		for i = 1:t
+			for j = 1:t
+				k = (i-1)*t+j;
+				J(:,:,k) = I(i:end-t+i,j:end-t+j);
+			end
+		end
+		somme = sum(J,3);
+		indices_carres_noirs = find(somme==0);
+		if length(indices_carres_noirs)>0
+			aucun_carre_noir_detecte = 0;
+		end
+		cpt_tirages = cpt_tirages+1;
+	end
+
+	% Affichage de l'image :
+	subplot(2,2,numero_image);
+	numero_image = numero_image+1;
+	imagesc(I);
+	axis off;
+	axis equal;
+	colormap gray;
+	set(gca,'FontSize',20);
+	nb_carres = length(indices_carres_noirs);
+	if nb_carres==1
+		title(['Il y a 1 carre noir ' num2str(t) '$\times$' num2str(t) ' : le voyez-vous ?'],'Interpreter','Latex');
+	else
+		title(['Il y a ' num2str(nb_carres) ' carres noirs ' num2str(t) '$\times$' num2str(t) ' : les voyez-vous ?'],'Interpreter','Latex');
+	end
+	hold on;
+	fprintf('Tapez Entree pour surligner les carres noirs\n');
+	pause;
+
+	% Surlignage des carres de cote t entierement noirs :
+	for k = 1:length(indices_carres_noirs)
+		l = indices_carres_noirs(k);
+		x = ceil(l/(T-t+1))-0.5;
+		y = mod(l-1,T-t+1)+1-0.5;
+		line([x,x],[y,y+t],'Color','r','Linewidth',2);
+		line([x+t,x+t],[y,y+t],'Color','r','Linewidth',2);
+		line([x,x+t],[y,y],'Color','r','Linewidth',2);
+		line([x,x+t],[y+t,y+t],'Color','r','Linewidth',2);
+	end
+	if t<liste_valeurs_t(end)
+		fprintf('Tapez Entree pour afficher une nouvelle image\n');
+		pause;
+	end
+end

TP3/exercice_1.m

@@ -0,0 +1,93 @@
+clear;
+close all;
+taille_ecran = get(0,'ScreenSize');
+L = taille_ecran(3);
+H = taille_ecran(4);
+
+% Parametres :
+set(0,'RecursionLimit',1000);				% Pour eviter les erreurs a l'execution
+alpha = pi/8;						% Seuil sur l'orientation du gradient (en radians)
+cos_alpha = cos(alpha);
+seuil_norme = 2/sin(alpha);				% Seuil sur la norme du gradient (decoule d'une etude theorique)
+
+% Lecture de l'image :
+I = imread('Images/Piree.png');
+% I = imread('Images/chaises.png');
+% I = imread('Images/Morlaix.png');
+if size(I,3)==3
+	I = rgb2gray(I);
+end
+I = imresize(I,0.8);					% Permet de limiter l'effet de crenelage
+[nb_lignes,nb_colonnes] = size(I);
+
+% Parametres :
+diagonale = floor(sqrt(nb_lignes^2+nb_colonnes^2));
+card_max = floor(diagonale/2);				% Cardinal max d'un ensemble E
+card_min = floor(card_max/20);				% Cardinal min d'un ensemble E
+
+% Affichage de l'image :
+figure('Name','Ensembles candidats','Position',[0,0,L,0.67*H]);
+subplot(1,2,1);
+imagesc(I);
+axis equal;
+axis off;
+colormap gray;
+hold on;
+
+% Gradient du niveau de gris (x vers la droite, y vers le bas) :
+I = double(I);
+[G_x,G_y] = gradient(I);
+G_norme = sqrt(G_x.^2+G_y.^2);
+
+% Selection des pixels de contour :
+contour = G_norme>seuil_norme;
+
+% Pas de pixel de contour sur le bord de l'image (traitement simplifie !) :
+contour([1,nb_lignes],:) = 0;
+contour(:,[1,nb_colonnes]) = 0;
+
+% Affichage d'une petite fleche sur les pixels de contour :
+[i,j] = find(contour); % i, j vecteurs d'indices ou les valeurs de contour sont non nulles
+quiver(j,i,G_x(contour),G_y(contour),'r');		% Attention : x = j et y = i
+axis equal;
+title('Pixels de contour','FontSize',30);
+drawnow;
+
+% Boucle de traitement des pixels de contour :
+voisins = [ -1 -1 ; -1 0 ; -1 1 ; 0 1 ;...
+		1 1 ; 1 0 ; 1 -1 ; 0 -1 ];		% Positions relatives des 8 plus proches voisins
+liste_E = [];
+while any(contour(:))
+
+	% Pixel de contour dont la norme du gradient est la plus elevee = "germe" :
+	indices_contour = find(contour); % indices linéaires des valeurs non nulles de contour
+	[~,indice_max] = max(G_norme(indices_contour)); % on choppe l'indice de la valeur maximale de la norme du gradient des contours
+	indice_germe = indices_contour(indice_max); % on se sert de cette indice comme de notre germe
+	[i_germe,j_germe] = ind2sub(size(I),indice_germe); % on converti notre indice linéaire en coordonnées
+    E = [i_germe,j_germe]; % on rempli E avec cette première germe
+	G_somme = [G_x(i_germe,j_germe),G_y(i_germe,j_germe)];
+
+	% Construction recursive de l'ensemble E :
+	[E,contour] = recursion(E,contour,G_somme,i_germe,j_germe,voisins,G_x,G_y,card_max,cos_alpha);
+
+	% Ajout de l'ensemble E a la liste d'ensembles E, si son cardinal est superieur a card_min :
+	if size(E,1)>=card_min
+		liste_E = [liste_E,{E}];		% Concatenation d'une cellule {E} a liste_E
+	end
+end
+
+% Affichage des ensembles E (la fonction label2rgb donne a chaque ensemble E une couleur differente) :
+etiquette = zeros(size(I));
+for k = 1:length(liste_E)
+	E = liste_E{k};
+	indices_E = sub2ind(size(I),E(:,1),E(:,2));
+	etiquette(ind2sub(size(I),indices_E)) = k;
+end
+I_resultat = label2rgb(etiquette,'parula',[0.5,0.5,0.5]);
+subplot(1,2,2);
+imagesc(I_resultat);
+axis equal;
+axis off;
+title('Ensembles candidats','FontSize',30);
+
+save exercice_1;

TP3/exercice_2.m

@@ -0,0 +1,65 @@
+clear;
+close all;
+
+load exercice_1;
+
+% Parametres :
+p = alpha/pi;
+epsilon = 1e-20;
+
+% Affichage de l'image :
+figure('Name','Detection des alignements','Position',[0,0,L,0.67*H]);
+subplot(1,2,1);
+imagesc(I);
+axis equal;
+axis off;
+colormap gray;
+title('L''image','FontSize',30);
+drawnow;
+
+% Affichage d'un fond gris pour l'esquisse :
+subplot(1,2,2);
+imagesc(120*ones(size(I)),[0,255]);
+axis equal;
+axis off;
+title('Son esquisse','FontSize',30);
+hold on;
+
+% Boucle sur les ensembles E :
+nb_elimines = 0;
+for k = 1:length(liste_E)
+
+	% Calcul de la matrice d'inertie de l'ensemble candidat numero k :
+	E = liste_E{k};
+	E_x = E(:,2);					% Attention : x = j
+	E_y = E(:,1);					% Attention : y = i
+	indices_E = sub2ind(size(I),E_y,E_x);
+	G_norme_E = G_norme(indices_E);
+	[C_x,C_y,M] = matrice_inertie(E_x,E_y,G_norme_E);
+
+	% Calcul des valeurs/vecteurs propres de M :
+	[V,D] = eig(M);
+
+	% Tri des valeurs propres de M par ordre decroissant :
+	[~,indices] = sort(diag(D),'descend');
+
+	% Matrice de passage du repere xy vers les axes propres du nuage de points :
+	P = V(:,indices);
+
+	% Coordonnees des points du nuage apres rotation dans le repere des axes principaux :
+	E_nouveau_repere = [E_x-C_x,E_y-C_y]*P;
+
+	% Calcul de probabilite (loi binomiale) :
+	[x_min,x_max,probabilite] = calcul_proba(E_nouveau_repere,p);
+
+	% Test de saillance visuelle :
+	if probabilite<epsilon				% Si l'evenement est suffisamment rare
+		% Affichage du segment :
+		extremites = P*[x_min x_max ; 0 0]+[C_x C_x ; C_y C_y];
+		line(extremites(1,:),extremites(2,:),'Color','w','LineWidth',2);
+		drawnow;
+	else
+		nb_elimines = nb_elimines+1;
+	end
+end
+fprintf( 'Sur %d ensembles candidats, %d ont ete elimines\n',length(liste_E),nb_elimines);

TP3/matrice_inertie.m

@@ -0,0 +1,12 @@
+function [x_bar, y_bar, M] = matrice_inertie(E_x, E_y, G_norme)
+
+    PI = sum(G_norme);
+    
+    x_bar = (G_norme.'*E_x)/PI;
+    y_bar = (G_norme.'*E_y)/PI;
+    
+    M = [ (G_norme'*(E_x-x_bar).^2)/PI, (G_norme'*((E_x-x_bar).*(E_y-y_bar)))/PI
+          (G_norme'*((E_x-x_bar).*(E_y-y_bar)))/PI, (G_norme'*(E_y-y_bar).^2)/PI ];
+
+end
+

TP3/recursion.m

@@ -0,0 +1,29 @@
+function [E,contour,G_somme] = recursion(E,contour,G_somme,i,j,voisins,G_x,G_y,card_max,cos_alpha)
+% Fonction recursive permettant de construire un ensemble candidat E
+
+    contour(i,j) = 0;
+    G_somme_normalise = G_somme/norm(G_somme);
+    nb_voisins = size(voisins,1);
+    k = 1;
+    while k<=nb_voisins & size(E,1)<=card_max
+        i_voisin = i + voisins(k, 1);
+        j_voisin = j + voisins(k, 2);
+        if contour(i_voisin,j_voisin)
+            
+            G_ij = [G_x(i_voisin, j_voisin), G_y(i_voisin, j_voisin)];
+            G_ij_normalise = G_ij/norm(G_ij);
+            
+            prod_scal = G_somme_normalise*G_ij_normalise';
+            
+            if prod_scal >= cos_alpha
+               
+                E = [E; i_voisin j_voisin];
+                G_somme = [G_somme; G_ij_normalise];
+                recursion(E, contour, G_somme, i_voisin, j_voisin, voisins, G_x, G_y, card_max, cos_alpha);
+                
+            end
+            
+        end
+        k = k+1;
+    end
+end

TP3/test_matrice_inertie.m

@@ -0,0 +1,13 @@
+[J,I] = meshgrid(1:5,1:5);
+G_norme = I+J;
+E = [ 1 2 ; 2 2 ; 2 3 ; 3 3 ; 3 4 ; 3 5 ];
+indices_E = sub2ind(size(G_norme),E(:,1),E(:,2));
+
+[C_x,C_y,M] = matrice_inertie(E(:,2),E(:,1),G_norme(indices_E));
+
+eps = 0.0001;
+if abs(C_x-3.4848)<eps & abs(C_y-2.5455)<eps & abs(M(1,1)-1.1589)<eps & abs(M(1,2)-0.5234)<eps & abs(M(2,2)-0.4298)<eps
+	fprintf('La fonction matrice_inertie semble correcte !\n');
+else
+	fprintf('Attention : probleme !!!\n');
+end