commit 937f37c7d3f92f1f17c4d8048232319d16ddeb46 Author: Laureηt Date: Sun Jun 25 16:38:01 2023 +0200 init diff --git a/TP1/E_estimee.mat b/TP1/E_estimee.mat new file mode 100644 index 0000000..42b887b Binary files /dev/null and b/TP1/E_estimee.mat differ diff --git a/TP1/Fontaine/1.png b/TP1/Fontaine/1.png new file mode 100644 index 0000000..0932fe5 Binary files /dev/null and b/TP1/Fontaine/1.png differ diff --git a/TP1/Fontaine/2.png b/TP1/Fontaine/2.png new file mode 100644 index 0000000..2aaf51c Binary files /dev/null and b/TP1/Fontaine/2.png differ diff --git a/TP1/Fontaine/R_sol.mat b/TP1/Fontaine/R_sol.mat new file mode 100644 index 0000000..8230981 Binary files /dev/null and b/TP1/Fontaine/R_sol.mat differ diff --git a/TP1/Fontaine/matK.mat b/TP1/Fontaine/matK.mat new file mode 100644 index 0000000..164ba76 Binary files /dev/null and b/TP1/Fontaine/matK.mat differ diff --git a/TP1/Gargouille/.DS_Store b/TP1/Gargouille/.DS_Store new file mode 100644 index 0000000..5008ddf Binary files /dev/null and b/TP1/Gargouille/.DS_Store differ diff --git a/TP1/Gargouille/1.JPG b/TP1/Gargouille/1.JPG new file mode 100644 index 0000000..08381a5 Binary files /dev/null and b/TP1/Gargouille/1.JPG differ diff --git a/TP1/Gargouille/2.JPG b/TP1/Gargouille/2.JPG new file mode 100644 index 0000000..3e79508 Binary files /dev/null and b/TP1/Gargouille/2.JPG differ diff --git a/TP1/Gargouille/R_sol.mat b/TP1/Gargouille/R_sol.mat new file mode 100644 index 0000000..6926abb Binary files /dev/null and b/TP1/Gargouille/R_sol.mat differ diff --git a/TP1/Gargouille/matK.mat b/TP1/Gargouille/matK.mat new file mode 100644 index 0000000..60d01a6 Binary files /dev/null and b/TP1/Gargouille/matK.mat differ diff --git a/TP1/affichage_3D.m b/TP1/affichage_3D.m new file mode 100644 index 0000000..8d3b794 --- /dev/null +++ b/TP1/affichage_3D.m @@ -0,0 +1,28 @@ +function affichage_3D(I_gauche,points_gauche_apparies,points_3D,t,R,numero_figure) + +% Récupération de la couleur des points 3D : +[nb_lignes_gauche,nb_colonnes_gauche,nb_canaux_gauche] = size(I_gauche); +nb_pixels_gauche = nb_lignes_gauche*nb_colonnes_gauche; +couleurs_image_gauche = reshape(I_gauche,[nb_pixels_gauche,nb_canaux_gauche]); +i_points_gauche_apparies = round(points_gauche_apparies(:,1)); +j_points_gauche_apparies = round(points_gauche_apparies(:,2)); +indices_couleurs = sub2ind([nb_lignes_gauche,nb_colonnes_gauche],j_points_gauche_apparies,i_points_gauche_apparies); +couleurs_3D = couleurs_image_gauche(indices_couleurs,:); + +% Affichage des caméras : +taille_camera = 0.2; +if nargin == 6 + subplot(2,2,numero_figure); +end +plotCamera('Size',taille_camera,'Color','r','Label','1','Opacity',0); +hold on; +grid on; +position = single(-R'*t); +orientation = single(R'); +absPose = rigidtform3d(orientation,position); +plotCamera('AbsolutePose',absPose,'Size',taille_camera,'Color','b','Label','2','Opacity',0); +axis off + +% Création et affichage du nuage de points 3D : +nuage_points_3D = pointCloud(points_3D,'Color',couleurs_3D); +pcshow(nuage_points_3D,'VerticalAxis','y','VerticalAxisDir','down','MarkerSize',45); diff --git a/TP1/donnees_appariees.mat b/TP1/donnees_appariees.mat new file mode 100644 index 0000000..1511c18 Binary files /dev/null and b/TP1/donnees_appariees.mat differ diff --git a/TP1/estimation_E.m b/TP1/estimation_E.m new file mode 100644 index 0000000..6147676 --- /dev/null +++ b/TP1/estimation_E.m @@ -0,0 +1,20 @@ +function E = estimation_E(w1, w2) + + A = zeros(size(w1, 1), 9); + for i=1:size(w1, 1) + A(i, :) = kron(w1(i,:), w2(i,:)); + end + + % on chope e à partir des valeurs propres de A'A + [V, D] = eig(A' * A); + e = V(:, 1); + + % on obtient E à partir de e + E = reshape(e, 3, 3); + + % on force les valeurs singulières de E + [U, S, V] = svd(E); + + E = U * [1 0 0; 0 1 0; 0 0 0] * V'; + +end \ No newline at end of file diff --git a/TP1/estimation_E_robuste.m b/TP1/estimation_E_robuste.m new file mode 100644 index 0000000..83b64a7 --- /dev/null +++ b/TP1/estimation_E_robuste.m @@ -0,0 +1,27 @@ +function [E, w1_new, w2_new] = estimation_E_robuste(w1, w2) + + nb_min = 8; + nb_paires = size(w1, 1); + + mediane = +inf; + S = 5e-5; + + while mediane > S + + selection = randperm(nb_paires, nb_min); + w1_new = w1(selection, :); + w2_new = w2(selection, :); + + E = estimation_E(w1_new, w2_new); + + d1 = w2 * E; + d2 = w1 * E'; + + dist1 = ( d1(:,1) .* w1(:, 1) + d1(:,2) .* w1(:, 2) + d1(:,3) ).^2 ./ ( d1(:,1).^2 + d1(:,2).^2 ); + dist2 = ( d2(:,1) .* w2(:, 1) + d2(:,2) .* w2(:, 2) + d2(:,3) ).^2 ./ ( d2(:,1).^2 + d2(:,2).^2 ); + + mediane = mean([dist1 ; dist2]); + + end + +end \ No newline at end of file diff --git a/TP1/estimation_t.m b/TP1/estimation_t.m new file mode 100644 index 0000000..3e39d66 --- /dev/null +++ b/TP1/estimation_t.m @@ -0,0 +1,7 @@ +function t_estimee = estimation_t(E_estimee, R_sol) + + T = E_estimee / R_sol; + + t_estimee = [ T(3,2) ; T(1,3) ; T(2,1) ]; + +end \ No newline at end of file diff --git a/TP1/exercice_0.m b/TP1/exercice_0.m new file mode 100644 index 0000000..9528855 --- /dev/null +++ b/TP1/exercice_0.m @@ -0,0 +1,116 @@ +clear; +close all; +taille_ecran = get(0,'ScreenSize'); +L = taille_ecran(3); +H = taille_ecran(4); + +% Choix des images : +% dossierImages = 'Gargouille/'; +dossierImages = 'Fontaine/'; +im1 = 1; +im2 = 2; +% ext = '.JPG'; +ext = '.png'; + +% Redimensionnement des images : +facteur = 1; + +% Images à la verticale : +% verticale = 1; +verticale = 0; + +% Utilisation des régions d'intérêt : +ROI = 0; + +% Options : +affichage = 1; % 0 pour éviter les affichages +nb_p_affiches = 400; + +% Noms des images : +nom_1 = strcat(dossierImages,int2str(im1),ext); +nom_2 = strcat(dossierImages,int2str(im2),ext); + +% Lecture des images : +I_1 = imresize(imread(nom_1),facteur); +I_2 = imresize(imread(nom_2),facteur); + +if verticale + I_1 = permute(I_1, [2,1,3]); % On inverse lignes et colonnes + I_1 = I_1(end:-1:1,:,:); % On remet l'image à l'endroit + I_2 = permute(I_2, [2,1,3]); % On inverse lignes et colonnes + I_2 = I_2(end:-1:1,:,:); % On remet l'image à l'endroit +end + +% Détection des points d'intérêt : +[nb_lignes,nb_colonnes,nb_canaux] = size(I_1); +I_1_nvg = rgb2gray(I_1); +I_2_nvg = rgb2gray(I_2); + +if ROI + RI_1 = selection_RI(I_1,L,H); + p_1 = detectMinEigenFeatures(I_1_nvg,'ROI',RI_1,'MinQuality',0.0001); + + RI_2 = selection_RI(I_2,L,H); + p_2 = detectMinEigenFeatures(I_2_nvg,'ROI',RI_2,'MinQuality',0.0001); +else + p_1 = detectMinEigenFeatures(I_1_nvg,'MinQuality',0.0001); + p_2 = detectMinEigenFeatures(I_2_nvg,'MinQuality',0.0001); +end + +% Affichage des images et des points d'intérêt : +if affichage + figure('Name','Detection des points d''interet','Position',[0.1*L,0.4*H,0.8*L,0.5*H]); + + subplot(1,2,1); + imagesc(I_1); + title('Image gauche','FontSize',20); + axis equal; + axis off; + hold on; + plot(p_1.selectStrongest(nb_p_affiches)); + + subplot(1,2,2); + imagesc(I_2); + title('Image droite','FontSize',20); + axis equal; + axis off; + hold on; + plot(p_2.selectStrongest(nb_p_affiches)); +end + +% Mise en correspondance : +traqueur = vision.PointTracker('MaxBidirectionalError',1,'NumPyramidLevels',5); +p_1 = p_1.Location; +initialize(traqueur,p_1,I_1); +[p_2,indices_droite] = step(traqueur,I_2); +p_1_apparies = p_1(indices_droite,:); +p_2_apparies = p_2(indices_droite,:); + +nb_paires = size(p_1_apparies,1); + +% Affichage des appariements : +if affichage + figure('Name','Mise en correspondance'); + pas_affichage = max(floor(nb_paires/nb_p_affiches),1); + showMatchedFeatures(I_1,I_2,p_1_apparies(1:pas_affichage:end,:),p_2_apparies(1:pas_affichage:end,:)); +end + +% Paramètres intrinsèques de la caméra : +nomFichier = strcat(dossierImages,'matK.mat'); +load(nomFichier); + +if verticale + u0 = K(1,3); + v0 = K(2,3); + + K(1,3) = v0; + K(2,3) = u0; +end + +if facteur~=1 + K = K*facteur; + K(3,3) = 1; +end +inverse_K = inv(K); + +save('donnees_appariees.mat','K','inverse_K','p_1_apparies','p_2_apparies','nb_paires','I_1','I_2','nb_lignes','nb_colonnes'); diff --git a/TP1/exercice_1.m b/TP1/exercice_1.m new file mode 100644 index 0000000..9e5920e --- /dev/null +++ b/TP1/exercice_1.m @@ -0,0 +1,26 @@ +clear; +close all; + +load donnees_appariees; +nb_min = 8; + +% Coordonnées homogènes des pixels appariés : +p_1_tilde = [p_1_apparies ones(nb_paires,1)]; +p_2_tilde = [p_2_apparies ones(nb_paires,1)]; + +% Points 3D w : +w_1 = transpose(inverse_K*p_1_tilde'); +w_2 = transpose(inverse_K*p_2_tilde'); + +% Tirage aléatoire de nb_min points : +selection = randperm(nb_paires,nb_min); +w_1_select = w_1(selection,:); +w_2_select = w_2(selection,:); + +% Estimation de E : +E_estimee = estimation_E(w_1_select,w_2_select); + +% Tracé des droites épipolaires passant par les points tirés aléatoirement : +trace_epipoles; + +save('E_estimee.mat','E_estimee'); diff --git a/TP1/exercice_2.m b/TP1/exercice_2.m new file mode 100644 index 0000000..104097f --- /dev/null +++ b/TP1/exercice_2.m @@ -0,0 +1,20 @@ +clear; +close all; + +load donnees_appariees; + +% Coordonnées homogènes dans le repère pixels : +p_1_tilde = [p_1_apparies ones(nb_paires,1)]; +p_2_tilde = [p_2_apparies ones(nb_paires,1)]; + +% Points 3D w : +w_1 = transpose(inverse_K*p_1_tilde'); +w_2 = transpose(inverse_K*p_2_tilde'); + +% Estimation robuste de E : +[E_estimee,w_1_select,w_2_select] = estimation_E_robuste(w_1,w_2); + +% Tracé des droites épipolaires passant par les paires de points sélectionnés : +trace_epipoles; + +save('E_estimee.mat','E_estimee'); diff --git a/TP1/exercice_3.m b/TP1/exercice_3.m new file mode 100644 index 0000000..076c3b1 --- /dev/null +++ b/TP1/exercice_3.m @@ -0,0 +1,23 @@ +clear; +close all; +taille_ecran = get(0,'ScreenSize'); +L = taille_ecran(3); +H = taille_ecran(4); + +load('./donnees_appariees.mat'); +load('./E_estimee.mat'); +load('./Fontaine/R_sol.mat'); + +p_1_tilde = [p_1_apparies ones(nb_paires,1)]; +p_2_tilde = [p_2_apparies ones(nb_paires,1)]; + +w_1 = transpose(inverse_K*p_1_tilde'); +w_2 = transpose(inverse_K*p_2_tilde'); + +% Estimation de la translation : +t_estimee = estimation_t(E_estimee,R_sol); + +% Calcul et affichage des points 3D reconstruits : +Q = reconstruction_3D(w_1,w_2,t_estimee,R_sol); +figure('Name','Reconstruction 3D','Position',[0.1*L,0,L,H]); +affichage_3D(I_1,p_1_apparies,Q,t_estimee,R_sol); diff --git a/TP1/reconstruction_3D.m b/TP1/reconstruction_3D.m new file mode 100644 index 0000000..72d14b2 --- /dev/null +++ b/TP1/reconstruction_3D.m @@ -0,0 +1,18 @@ +function Q = reconstruction_3D(w1, w2, t_estimee, R_sol) + + d = size(w1, 1); + Q = zeros(3, d); + + for i=1:d + A = [ -R_sol * w1(i,:)' , w2(i,:)', -t_estimee ]; + + [V, D] = eig(A' * A); + Z = V(:,1); + Z = Z / Z(3); + + Q(:,i) = 1/2 * (Z(2) * w2(i,:)' + R_sol * Z(1) * w1(i,:)' + t_estimee); + end + + Q = Q'; + +end \ No newline at end of file diff --git a/TP1/selection_RI.m b/TP1/selection_RI.m new file mode 100644 index 0000000..343c8ad --- /dev/null +++ b/TP1/selection_RI.m @@ -0,0 +1,23 @@ +function RI = selection_RI(image,L,H) + +[nb_lignes,nb_colonnes,nb_canaux] = size(image); +figure('Name','Selection de la region d''interet','Position',[0,0,0.5*L,0.5*H]); +imagesc(image); +axis equal; +axis off; + +disp('Cliquez pour determiner la region d''interet'); +[x_r,y_r] = ginput(2); +i_r = min(max(round(y_r),1),nb_lignes); +j_r = min(max(round(x_r),1),nb_colonnes); +j_min = min(j_r(:)); +j_max = max(j_r(:)); +i_min = min(i_r(:)); +i_max = max(i_r(:)); +line([j_min j_max],[i_min,i_min],'Color','r','LineWidth',2); +line([j_min j_max],[i_max,i_max],'Color','r','LineWidth',2); +line([j_min j_min],[i_min,i_max],'Color','r','LineWidth',2); +line([j_max j_max],[i_min,i_max],'Color','r','LineWidth',2); +drawnow; + +RI = [j_min,i_min,j_max-j_min,i_max-i_min]; diff --git a/TP1/sujet_TP1.pdf b/TP1/sujet_TP1.pdf new file mode 100644 index 0000000..92c0baa Binary files /dev/null and b/TP1/sujet_TP1.pdf differ diff --git a/TP1/trace_epipoles.m b/TP1/trace_epipoles.m new file mode 100644 index 0000000..09009b3 --- /dev/null +++ b/TP1/trace_epipoles.m @@ -0,0 +1,58 @@ +taille_ecran = get(0,'ScreenSize'); +L = taille_ecran(3); +H = taille_ecran(4); + +decalage_bord = 500; +I_bord = 240*ones(size(I_1,1),decalage_bord,3); +decalage_milieu = 1000; +I_milieu = 240*ones(size(I_1,1),decalage_milieu,3); + +figure('Name','Detection des points d''interet','Position',[0.1*L,H,0.8*L,H]); +I_1_2 = cat(2,I_bord,I_1,I_milieu,I_2,I_bord) ; +imagesc(I_1_2); +axis equal; +axis image off; +hold on; + +F = inverse_K'*E_estimee*inverse_K; + +p_1_select = transpose(K*w_1_select(1:8,:)'); +p_2_select = transpose(K*w_2_select(1:8,:)'); + +for i = 1:size(p_1_select,1) + + % Tracé coloré des pixels sélectionnés : + x_g = p_1_select(i,1)+decalage_bord; + y_g = p_1_select(i,2); + scatter(x_g,y_g,'r','filled'); + + x_d = p_2_select(i,1)+decalage_bord+decalage_milieu+size(I_1,2); + y_d = p_2_select(i,2); + scatter(x_d,y_d,'g','filled'); + + % Équations cartésiennes des droites épipolaires : + p1 = [x_g-decalage_bord ; y_g ; 1]; + D1 = F*p1; + p2 = [x_d-decalage_bord-nb_colonnes-decalage_milieu ; y_d ; 1]; + D2 = F'*p2; + + % Tracé de la droite épipolaire droite : + x_trace = -(decalage_bord/2):nb_colonnes+decalage_bord/2; + y_trace = -D1(1)/D1(2)*x_trace-D1(3)/D1(2); + indices = find(y_trace>0 & y_trace0 & y_trace S + + selection = randperm(nb_paires, nb_min); + w1_new = w1(selection, :); + w2_new = w2(selection, :); + + E = estimation_E(w1_new, w2_new); + + d1 = w2 * E; + d2 = w1 * E'; + + dist1 = ( d1(:,1) .* w1(:, 1) + d1(:,2) .* w1(:, 2) + d1(:,3) ).^2 ./ ( d1(:,1).^2 + d1(:,2).^2 ); + dist2 = ( d2(:,1) .* w2(:, 1) + d2(:,2) .* w2(:, 2) + d2(:,3) ).^2 ./ ( d2(:,1).^2 + d2(:,2).^2 ); + + mediane = mean([dist1 ; dist2]); + + end + +end \ No newline at end of file diff --git a/TP2/estimation_pose.m b/TP2/estimation_pose.m new file mode 100644 index 0000000..58824d8 --- /dev/null +++ b/TP2/estimation_pose.m @@ -0,0 +1,23 @@ +function [t, R] = estimation_pose(t4, R4, w1, w2) + + d = size(w1, 1); + Q = zeros(3, d); + + counts = zeros(4, 1); + for j=1:4 + for i=1:d + A = [ -R4(:,:,j) * w1(i,:)' , w2(i,:)']; + b = t4(:,j); + + Z = A \ b; + if Z(1) > 0 && Z(2) > 0 + counts(j) = counts(j) + 1; + end + end + end + + [~, indice] = max(counts); + t = t4(:, indice); + R = R4(:, :, indice); + +end \ No newline at end of file diff --git a/TP2/exercice_1.m b/TP2/exercice_1.m new file mode 100644 index 0000000..9925e4a --- /dev/null +++ b/TP2/exercice_1.m @@ -0,0 +1,22 @@ +clear; +close all; + +load E_estimee; + +% Estimation de la pose (4 solutions) : +[t_4,R_4] = estimation_4_poses(E_estimee); + +% Vérification : +verification = zeros(1,4); +for i = 1:4 + t = t_4(:,i); + t_chapeau = [0 -t(3) t(2) ; t(3) 0 -t(1) ; -t(2) t(1) 0]; + R = R_4(:,:,i); + matrice_verif = (t_chapeau*R)./E_estimee; + verification(i) = abs(abs(prod(matrice_verif(:)))-1)<1e-2; +end +if prod(verification) + fprintf('Vous avez tout bon : passez a l''exercice 2 !\n'); +else + fprintf('Erreur : revoyez votre code, svp !\n'); +end diff --git a/TP2/exercice_2.m b/TP2/exercice_2.m new file mode 100644 index 0000000..ceb3cf9 --- /dev/null +++ b/TP2/exercice_2.m @@ -0,0 +1,35 @@ +clear; +close all; + +load donnees_appariees; +load E_estimee; + +% Coordonnées homogènes des pixels appariés : +p_1_tilde = [p_1_apparies ones(nb_paires,1)]; +p_2_tilde = [p_2_apparies ones(nb_paires,1)]; +w_1 = transpose(inverse_K*p_1_tilde'); +w_2 = transpose(inverse_K*p_2_tilde'); + +% Estimation de la pose (4 solutions) : +[t_4,R_4] = estimation_4_poses(E_estimee); + +% Reconstruction 3D (4 solutions) : +figure('Name','Reconstruction 3D : 4 solutions'); +for i = 1:4 + t = t_4(:,i); + R = R_4(:,:,i); + Q = reconstruction_3D(w_1,w_2,t,R); + affichage_3D(I_1,p_1_apparies,Q,t,R,i); +end + +input('Tapez Entree pour afficher la bonne solution !'); + +% Détermination de la "bonne" pose : +[t,R] = estimation_pose(t_4,R_4,w_1,w_2); + +% Reconstruction 3D de la "bonne" solution : +Q = reconstruction_3D(w_1,w_2,t,R); + +% Affichage de la "bonne" solution : +figure('Name','Reconstruction 3D'); +affichage_3D(I_1,p_1_apparies,Q,t,R); diff --git a/TP2/exercice_3.m b/TP2/exercice_3.m new file mode 100644 index 0000000..e65ec71 --- /dev/null +++ b/TP2/exercice_3.m @@ -0,0 +1,78 @@ +clear; +close all; + +load donnees_SfM; + +% Liste de couleurs : +couleurs = uint8([ 0 255 255 ; 255 0 0 ; 0 255 0 ; 0 0 255 ; 255 255 0 ; 0 0 0 ]); + +% Initialisations : +t_1 = zeros(3,1); +R_1 = eye(3); +liste_t = t_1; +liste_R = R_1; +Q = []; +couleurs_points = []; + +% SfM à n > 2 images : +for i = 1:n-1 + + % Estimation du couple (t,R) : + E_estimee = matrices_E(:,:,i); + p_1_tilde = p_1_tilde_complet{i}; + p_2_tilde = p_2_tilde_complet{i}; + w_1 = transpose(inverse_K*p_1_tilde'); + w_2 = transpose(inverse_K*p_2_tilde'); + [t_4,R_4] = estimation_4_poses(E_estimee); + [t,R] = estimation_pose(t_4,R_4,w_1,w_2); + + % Estimation de l'échelle : + if i>1 + % Paires de points d'intérêt associées à la paire d'images {I_{i-1},I_i}} : + p_1_tilde_i_moins_1 = p_1_tilde_complet{i-1}; + [~,indices_i_moins_1,indices_i] = intersect(round(p_1_tilde_i_moins_1),round(p_1_tilde),'rows'); + + % Estimation du facteur d'échelle : + alpha_estime = estimation_echelle(Q_i_moins_1(indices_i_moins_1,:),w_2(indices_i,:),t,R); + t = alpha_estime*t; + end + + % Reconstruction 3D (méthode "du point milieu", cf. TP1) : + nouveaux_Q = reconstruction_3D(w_1,w_2,t,R); + + % Points 3D reconstruits, exprimés dans le repère de la caméra i+1 : + if i==1 + Q = nouveaux_Q; + else + Q = [ (R*Q')'+repmat(t',size(Q,1),1) ; nouveaux_Q ]; + end + Q_i_moins_1 = nouveaux_Q; + + % La couleur d'un point 3D est sa couleur dans l'image gauche : + I_1 = images{i}; + for j = 1:size(p_1_tilde,1) + currentColor = I_1(round(p_1_tilde(j,2)),round(p_1_tilde(j,1)),:); + couleurs_points = [ couleurs_points ; ipermute(currentColor,[1,3,2])]; + end + + % Concaténation des changements de pose : + liste_t = cat(2,liste_t,t); + liste_R = cat(3,liste_R,R); +end + +% Réorganisation des caméras pour l'affichage : +couleurs_cameras = [ couleurs(1:n,:) ; 255 255 255 ]; +for i = 1:size(liste_t,2) + if i==1 + liste_t_affichage = t_1; + liste_R_affichage = R_1; + else + liste_R_affichage(:,:,i) = liste_R(:, :,size(liste_t,2)-i+2)'*liste_R_affichage(:,:,i-1); + liste_t_affichage(:,i) = -liste_R_affichage(:,:,i)*liste_t(:,size(liste_t,2)... + -i+2)+liste_t_affichage(:,i-1); + end +end + +% Affichage des nuages de points 3D colorés : +figure('Name','Reconstruction 3D par SfM'); +affichage_3D_melange(Q,liste_t_affichage,liste_R_affichage,couleurs_points,couleurs_cameras); diff --git a/TP2/reconstru.png b/TP2/reconstru.png new file mode 100644 index 0000000..bd3b0e4 Binary files /dev/null and b/TP2/reconstru.png differ diff --git a/TP2/reconstruction_3D.m b/TP2/reconstruction_3D.m new file mode 100644 index 0000000..f2969ae --- /dev/null +++ b/TP2/reconstruction_3D.m @@ -0,0 +1,36 @@ +function Q = reconstruction_3D(w1, w2, t_estimee, R_sol) + + d = size(w1, 1); + Q = zeros(3, d); + + for i=1:d + A = [ -R_sol * w1(i,:)' , w2(i,:)']; + b = t_estimee; + + Z = A \ b; + + Q(:,i) = 1/2 * (Z(2) * w2(i,:)' + R_sol * Z(1) * w1(i,:)' + t_estimee); + end + + Q = Q'; + +end + +% function Q = reconstruction_3D(w1, w2, t_estimee, R_sol) +% +% d = size(w1, 1); +% Q = zeros(3, d); +% +% for i=1:d +% A = [ -R_sol * w1(i,:)' , w2(i,:)', -t_estimee ]; +% +% [V, D] = eig(A' * A); +% Z = V(:,1); +% Z = Z / Z(3); +% +% Q(:,i) = 1/2 * (Z(2) * w2(i,:)' + R_sol * Z(1) * w1(i,:)' + t_estimee); +% end +% +% Q = Q'; +% +% end \ No newline at end of file diff --git a/TP2/selection_RI.m b/TP2/selection_RI.m new file mode 100644 index 0000000..7b9e042 --- /dev/null +++ b/TP2/selection_RI.m @@ -0,0 +1,23 @@ +function RI = selection_RI(image) + +[nb_lignes,nb_colonnes,nb_canaux] = size(image); +figure('Name','Selection de la region d''interet'); +imagesc(image); +axis equal; +axis off; + +disp('Cliquez pour determiner la region d''interet'); +[x_r,y_r] = ginput(2); +i_r = min(max(round(y_r),1),nb_lignes); +j_r = min(max(round(x_r),1),nb_colonnes); +j_min = min(j_r(:)); +j_max = max(j_r(:)); +i_min = min(i_r(:)); +i_max = max(i_r(:)); +line([j_min j_max],[i_min,i_min],'Color','r','LineWidth',2); +line([j_min j_max],[i_max,i_max],'Color','r','LineWidth',2); +line([j_min j_min],[i_min,i_max],'Color','r','LineWidth',2); +line([j_max j_max],[i_min,i_max],'Color','r','LineWidth',2); +drawnow; + +RI = [j_min,i_min,j_max-j_min,i_max-i_min]; diff --git a/TP2/sfm.m b/TP2/sfm.m new file mode 100644 index 0000000..833d610 --- /dev/null +++ b/TP2/sfm.m @@ -0,0 +1,135 @@ +clear; +close all; + +% Choix des images : +dossierImages = 'Fontaine/'; +extension = '.png'; + +% Paramètres intrinsèques de la caméra : +nomFichier = strcat(dossierImages,'matK.mat'); +load(nomFichier); +inverse_K = inv(K); + +% Nombre d'images utilisées : +n = 4; + +% Utilisation des régions d'intérêt : +ROI = 0; + +% Lecture des images : +for i = 1:n + nom = strcat(dossierImages,int2str(i),extension); + images{i} = imread(nom); +end + +% Estimation des n-1 matrices essentielles : +matrices_E = []; +for i = 1:n-1 + + % Lecture des images i et i+1 : + I_1 = images{i}; + I_2 = images{i+1}; + + % Conversion en niveaux de gris : + I_1_nvg = rgb2gray(I_1); + I_2_nvg = rgb2gray(I_2); + + % Détection des points d'intérêt dans l'image gauche : + if ROI + RI_1 = selection_RI(I_1); + p_1 = detectMinEigenFeatures(I_1_nvg,'ROI',RI_1,'MinQuality',0.0001); + else + p_1 = detectMinEigenFeatures(I_1_nvg,'MinQuality',0.0001); + end + + % Mise en correspondance avec l'image droite : + traqueur = vision.PointTracker('MaxBidirectionalError',1,'NumPyramidLevels',5); + p_1 = p_1.Location; + initialize(traqueur,p_1,I_1); + [p_2,indices_droite] = step(traqueur,I_2); + p_1_apparies = p_1(indices_droite,:); + p_2_apparies = p_2(indices_droite,:); + + % Coordonnées homogènes des points image dans les deux repères : + nb_paires = size(p_1_apparies,1); + p_1_tilde = [p_1_apparies ones(nb_paires,1)]; + p_2_tilde = [p_2_apparies ones(nb_paires,1)]; + w_1 = transpose(inverse_K*p_1_tilde'); + w_2 = transpose(inverse_K*p_2_tilde'); + + % S'il y a trop de paires, on n'en conserve que 10000 : + nb_paires_majore = min(nb_paires,10000); + indices = randperm(nb_paires); + indices_select = indices(1:nb_paires_majore); + w_1_select = w_1(indices_select,:); + w_2_select = w_2(indices_select,:); + + % Estimation robuste de E : + E_estimee = estimation_E_robuste(w_1_select,w_2_select); + + matrices_E = cat(3,matrices_E,E_estimee); + p_1_tilde_complet{i} = p_1_tilde; + p_2_tilde_complet{i} = p_2_tilde; + +end + +save donnees_SfM; + +input('Tapez Entree pour lancer le SfM !'); + +% Liste de couleurs : +couleurs = uint8([ 0 255 255 ; 255 0 0 ; 0 255 0 ; 0 0 255 ; 255 255 0 ; 0 0 0 ]); + +% Initialisations : +t_1 = zeros(3,1); +R_1 = eye(3); +liste_t = t_1; +liste_R = R_1; +Q = []; +couleurs_points = []; + +% SfM à n > 2 images : +for i = 1:n-1 + + % Estimation de (t,R) : + E_estimee = matrices_E(:,:,i); + p_1_tilde = p_1_tilde_complet{i}; + p_2_tilde = p_2_tilde_complet{i}; + w_1 = transpose(inverse_K*p_1_tilde'); + w_2 = transpose(inverse_K*p_2_tilde'); + [t_4,R_4] = estimation_4_poses(E_estimee); + [t,R] = estimation_pose(t_4,R_4,w_1,w_2); + + % Reconstruction 3D (méthode "du point milieu", cf. TP1) : + nouveaux_Q = reconstruction_3D(w_1,w_2,t,R); + + % Points 3D reconstruits, exprimés dans le repère de la caméra i+1 : + if i==1 + Q = nouveaux_Q; + else + Q = [ (R*Q')'+repmat(t',size(Q,1),1) ; nouveaux_Q ]; + end + + % Une couleur différente est affectée à chaque nouveau nuage de points 3D : + couleurs_points = [ couleurs_points ; repmat(couleurs(i+1,:),size(nouveaux_Q,1),1) ]; + + liste_t = cat(2,liste_t,t); + liste_R = cat(3,liste_R,R); +end + +% Réorganisation des caméras pour l'affichage : +couleurs_cameras = [ couleurs(1:n,:) ; 255 255 255 ]; +for i = 1:size(liste_t,2) + if i==1 + liste_t_affichage = t_1; + liste_R_affichage = R_1; + else + liste_R_affichage(:,:,i) = liste_R(:,:,size(liste_t,2)-i+2)'*liste_R_affichage(:,:,i-1); + liste_t_affichage(:,i) = -liste_R_affichage(:,:,i)*liste_t(:,size(liste_t,2)... + -i+2)+liste_t_affichage(:,i-1); + end +end + +% Affichage des nuages de points 3D colorés : +figure('Name','Reconstruction 3D par SfM'); +affichage_3D_melange(Q,liste_t_affichage,liste_R_affichage,couleurs_points,couleurs_cameras); diff --git a/TP2/sujet_TP2.pdf b/TP2/sujet_TP2.pdf new file mode 100644 index 0000000..07a600e Binary files /dev/null and b/TP2/sujet_TP2.pdf differ diff --git a/TP3/MVS.m b/TP3/MVS.m new file mode 100644 index 0000000..721f301 --- /dev/null +++ b/TP3/MVS.m @@ -0,0 +1,61 @@ +function erreur_reprojection = MVS(I, k_ref, masque_ref, K, R, t, valeurs_z) + + % Indices des pixels du masque de ref + ind_masque = transpose(find(masque_ref > 0)); + [i_masque, j_masque] = ind2sub(size(masque_ref), ind_masque); + + % Passage du repère pixels au repère image : + u_masque = j_masque; + v_masque = i_masque; + p_tilde = [ u_masque ; v_masque ; ones(1, length(ind_masque)) ]; + + erreur_reprojection = zeros(size(p_tilde, 2), size(valeurs_z, 2)); + + % construction de la liste des indices des images sur lesquelles on va reprojeter + valeurs_k = [1:k_ref-1, k_ref+1:size(t, 2)]; + + % on va tester plusieurs valeurs de déprojection + for j = 1:size(valeurs_z, 2) + z = valeurs_z(j); + + % calcul de w à partir de p_tilde + Q = K \ p_tilde * z; + + % on déprojete le point dans le repère monde + Qm = R(:,:,k_ref)' * (Q - t(:, k_ref)); + + % pour chacune des autres image (non référence) + for i = 1:size(valeurs_k, 2) + k = valeurs_k(i); + + % on reprojète le point (repère monde) dans le repère de l'image k + Qk = R(:,:,k) * Qm + t(:,k); + wk = Qk ./ Qk(3,:); + + % on transforme Q en pixels + pk_tilde = K * wk; + + % on arrondi les coords des pixels + pk_tilde = round(pk_tilde); + + for o = 1:size(pk_tilde, 2) + i_pixel_reproj = pk_tilde(2, o); % v + j_pixel_reproj = pk_tilde(1, o); % u + i_pixel_ref = p_tilde(2, o); % v + j_pixel_ref = p_tilde(1, o); % u + + invisible = i_pixel_reproj <= 0 || i_pixel_reproj > 540 || j_pixel_reproj <= 0 || j_pixel_reproj > 540; + if invisible + erreur_reprojection(o, k) = Inf; + else + erreur = I(i_pixel_ref, j_pixel_ref, k_ref) - I(i_pixel_reproj, j_pixel_reproj, k); + erreur = erreur^2; + erreur_reprojection(o, j) = erreur; + end + end + + end + + end + +end \ No newline at end of file diff --git a/TP3/affichage_nuage.m b/TP3/affichage_nuage.m new file mode 100644 index 0000000..4167395 --- /dev/null +++ b/TP3/affichage_nuage.m @@ -0,0 +1,37 @@ +function affichage_nuage(masque,K,z,I) + +% Dimensions du masque : +[nb_lignes,nb_colonnes] = size(masque); +nb_pixels = nb_lignes*nb_colonnes; + +% Indices des pixels du masque : +ind_masque = transpose(find(masque>0)); +[i_masque,j_masque] = ind2sub(size(masque),ind_masque); + +% Passage du repère pixels au repère image : +u_masque = j_masque; +v_masque = i_masque; +p_tilde = [ u_masque ; v_masque ; ones(1,length(ind_masque)) ]; +w = inv(K)*p_tilde; + +% Lissage médian de la fonction de profondeur : +z_image = zeros(nb_pixels,1); +z_image(ind_masque) = z; +z_image = reshape(z_image,nb_lignes,nb_colonnes); +z_image_lissee = medfilt2(z_image,[3 3]); +z_lissee = z_image_lissee(ind_masque); +z_lissee = z_lissee(:); + +% Calcul des points 3D : +Q = NaN*ones(nb_pixels,3); +Q(ind_masque,:) = w'.*repmat(z_lissee,1,3); +Q = reshape(Q,[nb_lignes nb_colonnes 3]); + +% Affichage des points 3D coloriés par l'image de référence : +couleurs_3D = reshape(uint8(255*I),[nb_pixels,3]); +nuage_points_3D = pointCloud(reshape(Q,[nb_pixels,3]),'Color',couleurs_3D); +pcshow(nuage_points_3D,'VerticalAxis','y','VerticalAxisDir','down','MarkerSize',45); +axis equal; +axis off; +rotate3d; +zoom(1.6); diff --git a/TP3/affichage_surface.m b/TP3/affichage_surface.m new file mode 100644 index 0000000..dbf3d4e --- /dev/null +++ b/TP3/affichage_surface.m @@ -0,0 +1,31 @@ +function affichage_surface(masque,K,z) + +% Dimensions du masque : +[nb_lignes,nb_colonnes] = size(masque); +nb_pixels = nb_lignes*nb_colonnes; + +% Indices des pixels du masque : +ind_masque = transpose(find(masque>0)); +[i_masque,j_masque] = ind2sub(size(masque),ind_masque); + +% Passage du repère pixels au repère image : +u_masque = j_masque; +v_masque = i_masque; +p_tilde = [ u_masque ; v_masque ; ones(1,length(ind_masque)) ]; +w = inv(K)*p_tilde; + +% Calcul des points 3D : +Q = NaN*ones(nb_pixels,3); +Q(ind_masque,:) = w'.*repmat(z,1,3); +Q = reshape(Q,[nb_lignes nb_colonnes 3]); + +% Affichage des points 3D : +surfl(Q(:,:,1),Q(:,:,2),-Q(:,:,3),[0 90]); +shading flat; +colormap gray; +axis ij; +axis tight; +axis equal; +axis off; +rotate3d; +zoom(1.4); diff --git a/TP3/exercice_1.m b/TP3/exercice_1.m new file mode 100644 index 0000000..09a67da --- /dev/null +++ b/TP3/exercice_1.m @@ -0,0 +1,46 @@ +clear; +close all; + +load lapin; +valeurs_z = 1.5:0.002:2.5; +nb_images = size(I,3); + +% Masque de l'image de référence : +k_ref = 1; % Indice de l'image de référence +masque_ref = masque(:,:,k_ref); % Masque de l'image de référence + +% Affichage des images : +figure('Name','Images'); +for k = 1:nb_images + subplot(2,ceil(nb_images/2),k); + imagesc(I(:,:,k)); + colormap gray; + axis equal; + axis off; + if k=9; + +% Calcul de l'erreur de reprojection : +erreur_reprojection = MVS_2(I,k_ref,masque_ref,K,R,t,valeurs_z); + +% Affichage de l'erreur de reprojection en fonction de la profondeur : +figure('Name','Erreur de reprojection'); +imagesc(erreur_reprojection); +xlabel('Valeur de la profondeur','FontSize',20); +ylabel('Indice du pixel','FontSize',20); + +input('Tapez Entree pour afficher le relief reconstruit !'); + +% Affichage de la reconstruction 3D : +figure('Name','Reconstruction 3D'); +[~,indices_min] = min(erreur_reprojection,[],2); +z = transpose(valeurs_z(indices_min)); +affichage_surface(masque_ref,K,z); diff --git a/TP3/lapin.mat b/TP3/lapin.mat new file mode 100644 index 0000000..40e3833 Binary files /dev/null and b/TP3/lapin.mat differ diff --git a/TP3/sujet_TP3.pdf b/TP3/sujet_TP3.pdf new file mode 100644 index 0000000..01caa37 Binary files /dev/null and b/TP3/sujet_TP3.pdf differ diff --git a/TP4/sujet_TP4.pdf b/TP4/sujet_TP4.pdf new file mode 100644 index 0000000..37a3400 Binary files /dev/null and b/TP4/sujet_TP4.pdf differ diff --git a/TP5/conversion.m b/TP5/conversion.m new file mode 100644 index 0000000..7e9ce1b --- /dev/null +++ b/TP5/conversion.m @@ -0,0 +1,6 @@ +function [theta, phi] = conversion(s) + + phi = atan2(s(1,:), s(2,:)); + theta = acos(s(3,:)); + +end \ No newline at end of file diff --git a/TP5/donnees.mat b/TP5/donnees.mat new file mode 100644 index 0000000..e80ee4f Binary files /dev/null and b/TP5/donnees.mat differ diff --git a/TP5/eclairages.mat b/TP5/eclairages.mat new file mode 100644 index 0000000..bf048e4 Binary files /dev/null and b/TP5/eclairages.mat differ diff --git a/TP5/etalonnage.m b/TP5/etalonnage.m new file mode 100644 index 0000000..9fc3a9a --- /dev/null +++ b/TP5/etalonnage.m @@ -0,0 +1,40 @@ +function [directions] = etalonnage(spheres, exterieur_masque, centre, rayon) + + taille_noyau = round(rayon/20); + triangle = @(n)(1+n-abs(-n:n))/n; + + [height, width, N_photo, N_sphere] = size(spheres); + + directions = zeros(3, N_photo); + for p=1:N_photo + for s=1:N_sphere + % recup info sphere + image_sphere = spheres(:,:,p,s); + + % filtrage + image_filtree = conv2(triangle(taille_noyau), triangle(taille_noyau), image_sphere, 'same'); + image_filtree = image_filtree .* ( 1 - exterieur_masque); + + % recup point brillant + [~, I] = max(image_filtree, [], 'all'); + + % recup des coords + xc = centre(2); + yc = centre(1); + [x, y] = ind2sub([height, width], I); + + % calcul de la normale + normal = [ x-xc, y-yc, 0 ]; + normal(3) = sqrt(rayon^2 - normal(1)^2 - normal(2)^2); + normal = normal / norm(normal'); + + % calcul de la direction + axe_optique = [0, 0, 1]; + direction = 2 * (axe_optique * normal') * normal - axe_optique; + + directions(:, p) = directions(:, p) + direction'; + end + directions(:, p) = directions(:, p) / N_sphere; + directions(:, p) = directions(:, p) / norm(directions(:, p)); + end +end \ No newline at end of file diff --git a/TP5/exercice_1.m b/TP5/exercice_1.m new file mode 100644 index 0000000..167521b --- /dev/null +++ b/TP5/exercice_1.m @@ -0,0 +1,24 @@ +clear; +close all; + +load spheres; + +nb_lignes = size(spheres,1); +nb_colonnes = size(spheres,2); + +centre = [nb_lignes/2 nb_colonnes/2]; +rayon = min(nb_colonnes,nb_lignes)/2; + +[X,Y] = meshgrid(1:nb_colonnes,1:nb_lignes); +marge_rayon = 5; +exterieur_masque = (X-rayon).^2+(Y-rayon).^2>(rayon-marge_rayon)^2; + +s = etalonnage(spheres,exterieur_masque,centre,rayon); + +[theta,phi] = conversion(s); +plot(theta,phi,'o','Color','r','LineWidth',4,'MarkerSize',10); +xlabel('$\theta$','Interpreter','Latex','FontSize',20); +ylabel('$\phi$','Interpreter','Latex','FontSize',20); +axis([0,pi/2,-pi,pi]); + +save eclairages theta phi; diff --git a/TP5/exercice_2.m b/TP5/exercice_2.m new file mode 100644 index 0000000..de5bc76 --- /dev/null +++ b/TP5/exercice_2.m @@ -0,0 +1,42 @@ +clear; +close all; + +load donnees; +load eclairages; + +% Calcul des n fonctions d'interpolation : +[A,I] = interpolation(images,theta,phi); + +% Tirage aléatoire d'un pixel: +n = size(images,1)*size(images,2); +ind_test = randi(n); + +% Affichage des échantillons de ce pixel : +plot3(theta,phi,I(:,ind_test),'x','Color','r','LineWidth',4,'MarkerSize',10); +hold on; + +% Affichage de la fonction d'interpolation de ce pixel : +[phi_affichage,theta_affichage] = meshgrid(-pi:0.1:pi,0:0.05:(pi/2)); +z = A(1,ind_test)+A(2,ind_test)*theta_affichage+A(3,ind_test)*phi_affichage+... + A(4,ind_test)*theta_affichage.^2+A(5,ind_test)*theta_affichage.*phi_affichage+... + A(6,ind_test)*phi_affichage.^2; +surf(theta_affichage,phi_affichage,z); +xlabel('$\theta$','Interpreter','Latex','FontSize',30); +ylabel('$\phi$','Interpreter','Latex','FontSize',30); +zlabel('$I$','Interpreter','Latex','FontSize',30); +axis([0,pi/2,-pi,pi,0,1]); + +input('Tapez un caractere pour lancer la simulation !'); + +% Simulation d'un éclairage tournant : +close; +theta = 1.0; +valeurs_phi = -pi:0.2:pi; +for k = 1:length(valeurs_phi) + phi = valeurs_phi(k); + l = [1 theta phi theta^2 theta*phi phi^2]; + image_test = max(l*A,0); + image_test = reshape(image_test,[size(images,1),size(images,2)]); + imshow(image_test); + pause(0.01); +end diff --git a/TP5/exercice_3.m b/TP5/exercice_3.m new file mode 100644 index 0000000..b043642 --- /dev/null +++ b/TP5/exercice_3.m @@ -0,0 +1,50 @@ +clear; +close all; + +load donnees; +load eclairages; + +[~,I] = interpolation(images,theta,phi); +N = normales(I,theta,phi); + +x_affichage = 255*(N(:,1) + 1)/2; +y_affichage = 255*(N(:,2) + 1)/2; +z_affichage = 255*(N(:,3)); + +interieur = find(masque>0); % Interieur du domaine de reconstruction +exterieur = find(masque==0); % Exterieur du domaine de reconstruction + +N_affichage = reshape([x_affichage,y_affichage,z_affichage],[size(masque,1),size(masque,2),3]); + +% Intégration du champ de normales : +N(exterieur,3) = 1; % Pour éviter les divisions par 0 +p_estime = reshape(-N(:,1)./N(:,3),size(masque)); +p_estime(exterieur) = 0; +q_estime = reshape(-N(:,2)./N(:,3),size(masque)); +q_estime(exterieur) = 0; +z_estime = integration_SCS(-q_estime,p_estime); + +% Ambiguïté concave/convexe : +if (z_estime(floor(size(masque,1)/2),floor(size(masque,2)/2))y +% | +% x +% +% The problem comes to solve a linear system Ax=b, where A is a bloc +% Toplitz matrix. Fast solution is provided by Discrete Cosine Transform + +% Compute div(p,q): +px = 0.5*(p([2:end end],:)-p([1 1:end-1],:)); +qy = 0.5*(q(:,[2:end end])-q(:,[1 1:end-1])); + +% Div(p,q) + Boundary Condition: +f = px+qy; +f(1,2:end-1) = 0.5*(p(1,2:end-1)+p(2,2:end-1)); +f(end,2:end-1) = 0.5*(-p(end,2:end-1)-p(end-1,2:end-1)); +f(2:end-1,1) = 0.5*(q(2:end-1,1)+q(2:end-1,2)); +f(2:end-1,end) = 0.5*(-q(2:end-1,end)-q(2:end-1,end-1)); + +f(1,1) = 0.5*(p(1,1)+p(2,1)+q(1,1)+q(1,2)); +f(end,1) = 0.5*(-p(end,1)-p(end-1,1)+q(end,1)+q(end,2)); +f(1,end) = 0.5*(p(1,end)+p(2,end)-q(1,end)-q(1,end-1)); +f(end,end) = 0.5*(-p(end,end)-p(end-1,end)-q(end,end)-q(end,end-1)); + +% Cosine Transform of f: +fsin = dct2(f); + +% Denominator: +[x,y] = meshgrid(0:size(p,2)-1,0:size(p,1)-1); +denom = (2*cos(pi*x/(size(p,2)))-2) + (2*cos(pi*y/(size(p,1))) - 2); +Z = fsin./(denom); +Z(1,1) = 0.5*Z(1,2)+0.5*Z(2,1); % Or whatever... + +% Inverse Cosine Transform: +U = idct2(Z); + +% Ground truth available: +if (nargin>2) + moyenne_ecarts = mean(U(:)-gt(:)); + U = U-moyenne_ecarts; + npix = size(p,1)*size(p,2); + rmse = sqrt((sum((U(:)-gt(:)).^2))/npix); +else + U = U-min(U(:)); +end + +end diff --git a/TP5/spheres.mat b/TP5/spheres.mat new file mode 100644 index 0000000..e91b610 Binary files /dev/null and b/TP5/spheres.mat differ diff --git a/TP5/sujet_TP5.pdf b/TP5/sujet_TP5.pdf new file mode 100644 index 0000000..8d0c284 Binary files /dev/null and b/TP5/sujet_TP5.pdf differ diff --git a/TP6/Beethoven.png b/TP6/Beethoven.png new file mode 100644 index 0000000..afb881d Binary files /dev/null and b/TP6/Beethoven.png differ diff --git a/TP6/Beethoven_masque.png b/TP6/Beethoven_masque.png new file mode 100644 index 0000000..6db0d00 Binary files /dev/null and b/TP6/Beethoven_masque.png differ diff --git a/TP6/Lena.png b/TP6/Lena.png new file mode 100644 index 0000000..e7fd2f2 Binary files /dev/null and b/TP6/Lena.png differ diff --git a/TP6/delta_z.m b/TP6/delta_z.m new file mode 100644 index 0000000..1f51838 --- /dev/null +++ b/TP6/delta_z.m @@ -0,0 +1,11 @@ +function [d_z] = delta_z(I, voisins, distances) + + I_max = max(I, [], "all"); + abs_grad_z = sqrt( (I_max ./ I).^2 - 1 ); + + p1 = voisins(:, 1); + p2 = voisins(:, 2); + + d_z = distances .* ( abs_grad_z(p1) + abs_grad_z(p2) ) / 2; + +end \ No newline at end of file diff --git a/TP6/exercice_1.m b/TP6/exercice_1.m new file mode 100644 index 0000000..8037f99 --- /dev/null +++ b/TP6/exercice_1.m @@ -0,0 +1,71 @@ +clear; +close all; +taille_ecran = get(0,'ScreenSize'); +L = taille_ecran(3); +H = taille_ecran(4); + +I_max = 255; +pas_decimation = 1; +% nom_image = 'Beethoven.png'; +nom_image = "Lena.png" +nom_masque = 'Beethoven_masque.png'; +nb_iterations = 500; +epsilon = 1e-6; + +% Lecture et affichage de l'image : +I = imread(nom_image); +I = I(1:pas_decimation:end,1:pas_decimation:end,:); % Décimation de l'image +I = double(rgb2gray(I))/I_max; +[nb_lignes,nb_colonnes] = size(I); +figure('Name','Resolution iterative du SfS','Position',[0,0,0.5*L,0.7*H]); +subplot(2,2,1); +imagesc(I); +axis equal; +axis off; +hold on; +colormap gray; +title('Image','Interpreter','Latex','Fontsize',20); +drawnow; + +% Lecture du masque : +% masque = imread(nom_masque); +% masque = masque(1:pas_decimation:end,1:pas_decimation:end); % Décimation du masque +% masque = masque>0; +masque = ones(size(I)); +subplot(2,2,2); +imagesc(masque); +axis equal; +axis off; +hold on; +colormap gray; +title('Masque','Interpreter','Latex','Fontsize',20); +drawnow; + +% Membre droit de l'équation eikonale (attention aux troncations) ; +f = sqrt(1./(min(0.95,max(I,0.05)).^2)-1); + +% Initialisation de la solution : +z = zeros(nb_lignes,nb_colonnes); + +% Schéma itératif : +for k = 1:nb_iterations + + % Sauvegarde de l'itération courante : + z_k_moins_1 = z; + + % Pas de l'itération : + z = lax_friedrichs(z_k_moins_1,f,masque); + + % Affichage du résultat courant : + relief_et_image(z,k); + pause(0.001); + + % Test de convergence : + if (norm(z_k_moins_1(masque>0)-z(masque>0))./norm(z(masque>0))0); % Intérieur du domaine de reconstruction +exterieur = find(masque==0); % Extérieur du domaine de reconstruction +[n,p] = size(I); % n = nombre d'images, p = nombre de pixels + +% Affichage des images : +figure('Name','Donnees','Position',[0.8*L,0,0.2*L,H]); +if n==3 + n_l = n; + n_c = 1; +else + n_l = 7; + n_c = floor(n/7); +end +for i = 1:n_l*n_c + img = reshape(I(i,:),nb_lignes,nb_colonnes); + subplot(n_l,n_c,i); + imshow(uint8(img)); + hold on; + axis image; + axis off; + title(['$I_{' num2str(i,'%2d') '}$'],'Interpreter','Latex','FontSize',20); +end +colormap gray; +drawnow; + +% Estimation de M, rho et N : +[rho_estime,N_estime] = estimation(I,S,masque); + +% Intégration du champ de normales : +N_estime(3,exterieur) = 1; % Pour éviter les divisions par 0 +p_estime = reshape(-N_estime(1,:)./N_estime(3,:),size(masque)); +p_estime(exterieur) = 0; +q_estime = reshape(-N_estime(2,:)./N_estime(3,:),size(masque)); +q_estime(exterieur) = 0; +z_estime = integration_SCS(q_estime,p_estime); + +% Ambiguïté concave/convexe : +if (z_estime(floor(nb_lignes/2),floor(nb_colonnes/2))0); % Indices of non-masked pixels + +B_ = B(Iinds,:); + +% Compute dB/dx and dB/dy: +smooth_kernel = fspecial('gaussian',[25 25],2.5); +%~ smooth_kernel = fspecial('disk',3.0); +%~ smooth_kernel = 1; +for i = 1:3 + B_smooth = imfilter(reshape(B(:,i).*mask_integrability(:),Isize),smooth_kernel,'same'); + %~ B_smooth = reshape(B(:,i).*mask_integrability(:),Isize); + [dBidx dBidy] = gradient(B_smooth); + dBidy = -dBidy; + dB_dx(:,i) = dBidx(Iinds); + dB_dy(:,i) = dBidy(Iinds); +end + +% Solve the following system of equations for x_ij and y_ij: +% sum_{j=2:3} sum_{i=1:j-1} ( a_ij*x_ij - b_ij*y_ij) = 0 +% where, +% a_ij = e(i)*dedx(j)-e(j)*dedx(i) +% c_ij = P_3i*P_2j - P_2i*P_3j +% b_ij = e(i)*dedy(j)-e(j)*dedy(i) +% d_ij = P_3i*P_1j - P_1i*P_3j +a12 = B_(:,1).*dB_dx(:,2) - B_(:,2).*dB_dx(:,1); +a13 = B_(:,1).*dB_dx(:,3) - B_(:,3).*dB_dx(:,1); +a23 = B_(:,2).*dB_dx(:,3) - B_(:,3).*dB_dx(:,2); +b12 = B_(:,1).*dB_dy(:,2) - B_(:,2).*dB_dy(:,1); +b13 = B_(:,1).*dB_dy(:,3) - B_(:,3).*dB_dy(:,1); +b23 = B_(:,2).*dB_dy(:,3) - B_(:,3).*dB_dy(:,2); + +% Solve the system A*x = 0: +% A = [a12 a13 a23 -b12 -b13 -b23] +% x = [c12 c13 c23 d12 d13 d23] +A = [a12 a13 a23 -b12 -b13 -b23 ]; +[AU,AS,AV] = svd(A(:,1:6),0); +x = AV(:,size(AV,2)); + +% Construct the "co-factor matrix" (eqn 12) from the paper: +% A = [-c23 d23 ???; +% c13 -d13 ???; +% -c12 d12 ???]; +mean_x = mean(abs(x(1:6))); +Ainv = [-x(3) x(6) 1.2*mean_x; ... + x(2) -x(5) 0.5*mean_x; ... + -x(1) x(4) -mean_x]'; + +% Invert Ainv to get A: +A = inv(Ainv); + +% Compute the B and S matrices from P3 and Q3: +B = B*A; + +if (mean(B(:,3)<0)) + B = -B; +end + +M1 = transpose(B); + +end diff --git a/TP7/integration_SCS.m b/TP7/integration_SCS.m new file mode 100644 index 0000000..11f0f53 --- /dev/null +++ b/TP7/integration_SCS.m @@ -0,0 +1,62 @@ +function [U,rmse] = integration_SCS(p,q,gt) +% Ref : Direct Analytical Methods for Solving Poisson Equation in +% Computer Vision problems - PAMI 1990 +% +% example : +% p = ones(100,100); +% q = ones(100,100); +% u = simchony(p,q); +% +% This performs the least square solution to \nabla u = [p,q], i.e. : +% min \int_\Omega \| \nablua U - [p,q] \|^2 +% where \Omega is square and the natural Neumann boundary condition +% \mu \cdotp (\nabla u -[p,q]) = 0 is used (2nd order derivatives +% are neglected), where \mu is the outer +% normal to \Omega. +% +% % Axis : O->y +% | +% x +% +% The problem comes to solve a linear system Ax=b, where A is a bloc +% Toplitz matrix. Fast solution is provided by Discrete Cosine Transform + +% Compute div(p,q): +px = 0.5*(p([2:end end],:)-p([1 1:end-1],:)); +qy = 0.5*(q(:,[2:end end])-q(:,[1 1:end-1])); + +% Div(p,q) + Boundary Condition: +f = px+qy; +f(1,2:end-1) = 0.5*(p(1,2:end-1)+p(2,2:end-1)); +f(end,2:end-1) = 0.5*(-p(end,2:end-1)-p(end-1,2:end-1)); +f(2:end-1,1) = 0.5*(q(2:end-1,1)+q(2:end-1,2)); +f(2:end-1,end) = 0.5*(-q(2:end-1,end)-q(2:end-1,end-1)); + +f(1,1) = 0.5*(p(1,1)+p(2,1)+q(1,1)+q(1,2)); +f(end,1) = 0.5*(-p(end,1)-p(end-1,1)+q(end,1)+q(end,2)); +f(1,end) = 0.5*(p(1,end)+p(2,end)-q(1,end)-q(1,end-1)); +f(end,end) = 0.5*(-p(end,end)-p(end-1,end)-q(end,end)-q(end,end-1)); + +% Cosine Transform of f: +fsin = dct2(f); + +% Denominator: +[x,y] = meshgrid(0:size(p,2)-1,0:size(p,1)-1); +denom = (2*cos(pi*x/(size(p,2)))-2) + (2*cos(pi*y/(size(p,1))) - 2); +Z = fsin./(denom); +Z(1,1) = 0.5*Z(1,2)+0.5*Z(2,1); % Or whatever... + +% Inverse Cosine Transform: +U = idct2(Z); + +% Ground truth available: +if (nargin>2) + moyenne_ecarts = mean(U(:)-gt(:)); + U = U-moyenne_ecarts; + npix = size(p,1)*size(p,2); + rmse = sqrt((sum((U(:)-gt(:)).^2))/npix); +else + U = U-min(U(:)); +end + +end diff --git a/TP7/modele.png b/TP7/modele.png new file mode 100644 index 0000000..1b7acfe Binary files /dev/null and b/TP7/modele.png differ diff --git a/TP7/sujet_TP7.pdf b/TP7/sujet_TP7.pdf new file mode 100644 index 0000000..56712fc Binary files /dev/null and b/TP7/sujet_TP7.pdf differ diff --git a/TP7/visages.png b/TP7/visages.png new file mode 100644 index 0000000..a34dd65 Binary files /dev/null and b/TP7/visages.png differ