ENSEEIHT/TP-modelisation-images
1
0
Fork 0
Code Issues Pull requests Packages Projects Releases Wiki Activity

style: autoformat

Browse source
This commit is contained in:
Laurent Fainsin 2022-04-06 22:17:11 +02:00
parent 74365d944c
commit 2cf841e339
2 changed files with 55 additions and 58 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

97
TP_maillage.m
View file

@ -2,12 +2,12 @@ clear;
close all; close all;
load("mask.mat"); load("mask.mat");
im_mask = im_mask(:,:,:) == 0; im_mask = im_mask(:, :, :) == 0;
im_mask(1:8,:,:) = 0; im_mask(1:8, :, :) = 0;
im_mask(end-5:end,:,:) = 0; im_mask(end - 5:end, :, :) = 0;
im_mask(:,1:5,:) = 0; im_mask(:, 1:5, :) = 0;
im_mask(:,end-30:end,:) = 0; im_mask(:, end - 30:end, :) = 0;
nb_images = 36; % Nombre d'images nb_images = 36; % Nombre d'images
@ -100,14 +100,14 @@ germes = [germes W];
%% %%
figure(4); figure(4);
tiledlayout(2,2, 'Padding', 'none', 'TileSpacing', 'compact'); tiledlayout(2, 2, 'Padding', 'none', 'TileSpacing', 'compact');
bw_img = reshape(germes(image_labelise, 6), size(im, 1), []); bw_img = reshape(germes(image_labelise, 6), size(im, 1), []);
nexttile; nexttile;
imshow(bw_img); imshow(bw_img);
CC = bwconncomp(bw_img,4); CC = bwconncomp(bw_img, 4);
bw_img = zeros(size(bw_img)); bw_img = zeros(size(bw_img));
bw_img(CC.PixelIdxList{1}) = 1; bw_img(CC.PixelIdxList{1}) = 1;
@ -116,7 +116,7 @@ imshow(bw_img);
bw_img = ~bw_img; bw_img = ~bw_img;
CC = bwconncomp(bw_img,4); CC = bwconncomp(bw_img, 4);
bw_img = zeros(size(bw_img)); bw_img = zeros(size(bw_img));
bw_img(CC.PixelIdxList{1}) = 1; bw_img(CC.PixelIdxList{1}) = 1;
@ -125,31 +125,27 @@ bw_img = ~bw_img;
nexttile; nexttile;
imshow(bw_img); imshow(bw_img);
bw_img = im_mask(:,:,ind_img); bw_img = im_mask(:, :, ind_img);
nexttile; nexttile;
imshow(bw_img); imshow(bw_img);
figure(5); figure(5);
tiledlayout(2,2, 'Padding', 'none', 'TileSpacing', 'compact'); tiledlayout(2, 2, 'Padding', 'none', 'TileSpacing', 'compact');
nexttile; nexttile;
imshow(bw_img); imshow(bw_img);
hold on hold on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A FAIRE SI VOUS UTILISEZ LES MASQUES BINAIRES FOURNIS % % A FAIRE SI VOUS UTILISEZ LES MASQUES BINAIRES FOURNIS %
% Chargement des masques binaires % % Chargement des masques binaires %
% de taille nb_lignes x nb_colonnes x nb_images % % de taille nb_lignes x nb_colonnes x nb_images %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pixel_b = find(bw_img == 1); pixel_b = find(bw_img == 1);
[r, c] = ind2sub(size(bw_img), pixel_b(1)); [r, c] = ind2sub(size(bw_img), pixel_b(1));
contour = bwtraceboundary(bw_img, [r c], 'W', 8); contour = bwtraceboundary(bw_img, [r c], 'W', 8);
plot(contour(:,2), contour(:,1), 'g', 'LineWidth', 3); plot(contour(:, 2), contour(:, 1), 'g', 'LineWidth', 3);
% r = delaunay(contour); % r = delaunay(contour);
% barycentres = (contour(r(:,1),:) + contour(r(:,2),:) + contour(r(:,3),:)) / 3; % barycentres = (contour(r(:,1),:) + contour(r(:,2),:) + contour(r(:,3),:)) / 3;
@ -157,18 +153,17 @@ plot(contour(:,2), contour(:,1), 'g', 'LineWidth', 3);
% % triplot(r, contour(:,2), contour(:,1)); % % triplot(r, contour(:,2), contour(:,1));
T = 1; T = 1;
[vx, vy] = voronoi(contour(1:T:end,1), contour(1:T:end,2)); [vx, vy] = voronoi(contour(1:T:end, 1), contour(1:T:end, 2));
plot(vy, vx, 'b'); plot(vy, vx, 'b');
% Selection des segments qui ont leurs extrémités dans l'image % Selection des segments qui ont leurs extrémités dans l'image
ok = vx(1,:) > 1 & vx(1,:) < size(bw_img, 1) & ... ok = vx(1, :) > 1 & vx(1, :) < size(bw_img, 1) & ...
vx(2,:) > 1 & vx(2,:) < size(bw_img, 1) & ... vx(2, :) > 1 & vx(2, :) < size(bw_img, 1) & ...
vy(1,:) > 1 & vy(1,:) < size(bw_img, 2) & ... vy(1, :) > 1 & vy(1, :) < size(bw_img, 2) & ...
vy(2,:) > 1 & vy(2,:) < size(bw_img, 2); vy(2, :) > 1 & vy(2, :) < size(bw_img, 2);
vx = floor(vx(:,ok)); vx = floor(vx(:, ok));
vy = floor(vy(:,ok)); vy = floor(vy(:, ok));
% subplot(2,2,2); % subplot(2,2,2);
nexttile; nexttile;
@ -178,15 +173,15 @@ plot(vy, vx, 'b');
% Selection des segments avec les extremités dans la forme % Selection des segments avec les extremités dans la forme
ind1 = sub2ind(size(bw_img), vx(1,:), vy(1,:)); ind1 = sub2ind(size(bw_img), vx(1, :), vy(1, :));
ok1 = bw_img(ind1) > 0; ok1 = bw_img(ind1) > 0;
ind2 = sub2ind(size(bw_img), vx(2,:), vy(2,:)); ind2 = sub2ind(size(bw_img), vx(2, :), vy(2, :));
ok2 = bw_img(ind2) > 0; ok2 = bw_img(ind2) > 0;
ok = ok1 & ok2; ok = ok1 & ok2;
vx = vx(:,ok); vx = vx(:, ok);
vy = vy(:,ok); vy = vy(:, ok);
% subplot(2,2,3); % subplot(2,2,3);
nexttile; nexttile;
@ -194,59 +189,56 @@ imshow(bw_img);
hold on hold on
plot(vy, vx, 'b'); plot(vy, vx, 'b');
% Remise en forme de vx et vy % Remise en forme de vx et vy
vx_ = vx'; vx_ = vx';
vx_ = [vx_(:,1) ; vx_(:,2)]; vx_ = [vx_(:, 1); vx_(:, 2)];
vy_ = vy'; vy_ = vy';
vy_ = [vy_(:,1) ; vy_(:,2)]; vy_ = [vy_(:, 1); vy_(:, 2)];
V_ = [vx_ vy_]; V_ = [vx_ vy_];
% Calcule des rayons % Calcule des rayons
contour_ = contour'; contour_ = contour';
R = complex(V_(:,1), V_(:,2)) - complex(contour_(1,:), contour_(2,:)); R = complex(V_(:, 1), V_(:, 2)) - complex(contour_(1, :), contour_(2, :));
R = abs(R); R = abs(R);
R = min(R, [], 2); R = min(R, [], 2);
R = R'; R = R';
vx_vy = [vy(1,:) vy(2,:); vx(1,:) vx(2,:)]'; vx_vy = [vy(1, :) vy(2, :); vx(1, :) vx(2, :)]';
viscircles(vx_vy(1:1:end,:), R(1:1:end)); viscircles(vx_vy(1:1:end, :), R(1:1:end));
% Filtrage naif % Filtrage naif
% truc = find(R < 20); % truc = find(R < 20);
% truc = mod(truc - 1, length(vx)) + 1; % truc = mod(truc - 1, length(vx)) + 1;
% %
% vx(:,truc) = []; % vx(:,truc) = [];
% vy(:,truc) = []; % vy(:,truc) = [];
% Filtrage scalaire % Filtrage scalaire
R_scaled = 1.05 * R; R_scaled = 1.05 * R;
dist = abs(complex(V_(:,1), V_(:,2)) - transpose(complex(V_(:,1), V_(:,2)))); dist = abs(complex(V_(:, 1), V_(:, 2)) - transpose(complex(V_(:, 1), V_(:, 2))));
R_vertical = ones(length(R_scaled),1) * R_scaled; R_vertical = ones(length(R_scaled), 1) * R_scaled;
R_horizontal = R_scaled' * ones(1,length(R_scaled)); R_horizontal = R_scaled' * ones(1, length(R_scaled));
[~, c] = ind2sub(size(dist), find(dist + R_vertical < R_horizontal)); [~, c] = ind2sub(size(dist), find(dist + R_vertical < R_horizontal));
vx(:, mod(c - 1, length(vx)) + 1) = []; vx(:, mod(c - 1, length(vx)) + 1) = [];
vy(:, mod(c - 1, length(vy)) + 1) = []; vy(:, mod(c - 1, length(vy)) + 1) = [];
R = [R(1:length(R)/2); R(length(R)/2+1:end)]; R = [R(1:length(R) / 2); R(length(R) / 2 + 1:end)];
R(:, mod(c - 1, length(R)) + 1) = []; R(:, mod(c - 1, length(R)) + 1) = [];
R = [R(1,:) R(2,:)]; R = [R(1, :) R(2, :)];
% subplot(2,2,4); % subplot(2,2,4);
nexttile; nexttile;
imshow(bw_img); imshow(bw_img);
hold on hold on
plot(vy, vx, 'b'); plot(vy, vx, 'b');
vx_vy = [vy(1,:) vy(2,:); vx(1,:) vx(2,:)]'; vx_vy = [vy(1, :) vy(2, :); vx(1, :) vx(2, :)]';
viscircles(vx_vy(1:1:end,:), R(1:1:end)); viscircles(vx_vy(1:1:end, :), R(1:1:end));
%% %%
@ -307,7 +299,7 @@ view(80, -10);
% A COMPLETER % A COMPLETER
% Tetraedrisation de Delaunay % Tetraedrisation de Delaunay
T = delaunayTriangulation(X(1,:)', X(2,:)', X(3,:)'); T = delaunayTriangulation(X(1, :)', X(2, :)', X(3, :)');
% A DECOMMENTER POUR AFFICHER LE MAILLAGE % A DECOMMENTER POUR AFFICHER LE MAILLAGE
% fprintf('Tetraedrisation terminee : %d tetraedres trouves. \n',size(T,1)); % fprintf('Tetraedrisation terminee : %d tetraedres trouves. \n',size(T,1));
@ -315,13 +307,13 @@ T = delaunayTriangulation(X(1,:)', X(2,:)', X(3,:)');
% figure; % figure;
% tetramesh(T); % tetramesh(T);
% A DECOMMENTER ET A COMPLETER % A DECOMMENTER ET A COMPLETER
% Calcul des barycentres de chacun des tetraedres % Calcul des barycentres de chacun des tetraedres
poids = [1 1 1 1] / 4; poids = [1 1 1 1] / 4;
nb_barycentres = size(T.ConnectivityList, 1); nb_barycentres = size(T.ConnectivityList, 1);
for i = 1:size(T,1)
for i = 1:size(T, 1)
% Calcul des barycentres differents en fonction des poids differents % Calcul des barycentres differents en fonction des poids differents
% En commencant par le barycentre avec poids uniformes % En commencant par le barycentre avec poids uniformes
C_g(1:3, i) = poids * T.Points(T.ConnectivityList(i, :), :); C_g(1:3, i) = poids * T.Points(T.ConnectivityList(i, :), :);
@ -350,18 +342,21 @@ tri = T.ConnectivityList;
% ne se trouvent pas dans au moins un des masques des images de travail % ne se trouvent pas dans au moins un des masques des images de travail
% Pour chaque barycentre % Pour chaque barycentre
to_save = []; to_save = [];
for k = 1:nb_barycentres for k = 1:nb_barycentres
valide = 0; valide = 0;
for i = 1:nb_images for i = 1:nb_images
o = P{i}*C_g(:,k); o = P{i} * C_g(:, k);
o = o / o(3); o = o / o(3);
x = floor(o(1)); x = floor(o(1));
y = floor(o(2)); y = floor(o(2));
if im_mask(x,y,i) == 0 if im_mask(x, y, i) == 0
valide = 1; valide = 1;
break break
end end
end end
if valide if valide
@ -371,15 +366,15 @@ for k = 1:nb_barycentres
to_save = [to_save k]; to_save = [to_save k];
end end
triBis = tri(to_save,:); triBis = tri(to_save, :);
nb_barycentres = length(to_save); nb_barycentres = length(to_save);
% A DECOMMENTER POUR AFFICHER LE MAILLAGE RESULTAT % A DECOMMENTER POUR AFFICHER LE MAILLAGE RESULTAT
% Affichage des tetraedres restants % Affichage des tetraedres restants
fprintf('Retrait des tetraedres exterieurs a la forme 3D termine : %d tetraedres restants. \n',size(T,1)); fprintf('Retrait des tetraedres exterieurs a la forme 3D termine : %d tetraedres restants. \n', size(T, 1));
figure(7); figure(7);
tetramesh(triBis, T.Points); tetramesh(triBis, T.Points);
view(80, -10); view(80, -10);
% Sauvegarde des donnees % Sauvegarde des donnees
save donnees; save donnees;

16
TP_maillage_suite.m
View file

@ -3,17 +3,19 @@ load donnees;
% Calcul des faces du maillage à garder % Calcul des faces du maillage à garder
FACES = [sort(triBis(:, [2 3 4]), 2); sort(triBis(:, [1 3 4]), 2); sort(triBis(:, [1 2 4]), 2); sort(triBis(:, [1 2 3]), 2)]; FACES = [sort(triBis(:, [2 3 4]), 2); sort(triBis(:, [1 3 4]), 2); sort(triBis(:, [1 2 4]), 2); sort(triBis(:, [1 2 3]), 2)];
FACES = sortrows(FACES); FACES = sortrows(FACES);
rep = sum(FACES(1:end-1, :) == FACES(2:end, :), 2) == 3; rep = sum(FACES(1:end - 1, :) == FACES(2:end, :), 2) == 3;
FACES([0; rep] | [rep; 0], :) = []; FACES([0; rep] | [rep; 0], :) = [];
fprintf('Calcul du maillage final termine : %d faces. \n',size(FACES,1)); fprintf('Calcul du maillage final termine : %d faces. \n', size(FACES, 1));
% Affichage du maillage final % Affichage du maillage final
figure; figure;
hold on hold on
for i = 1:size(FACES,1)
plot3([X(1,FACES(i,1)) X(1,FACES(i,2))],[X(2,FACES(i,1)) X(2,FACES(i,2))],[X(3,FACES(i,1)) X(3,FACES(i,2))],'r'); for i = 1:size(FACES, 1)
plot3([X(1,FACES(i,1)) X(1,FACES(i,3))],[X(2,FACES(i,1)) X(2,FACES(i,3))],[X(3,FACES(i,1)) X(3,FACES(i,3))],'r'); plot3([X(1, FACES(i, 1)) X(1, FACES(i, 2))], [X(2, FACES(i, 1)) X(2, FACES(i, 2))], [X(3, FACES(i, 1)) X(3, FACES(i, 2))], 'r');
plot3([X(1,FACES(i,3)) X(1,FACES(i,2))],[X(2,FACES(i,3)) X(2,FACES(i,2))],[X(3,FACES(i,3)) X(3,FACES(i,2))],'r'); plot3([X(1, FACES(i, 1)) X(1, FACES(i, 3))], [X(2, FACES(i, 1)) X(2, FACES(i, 3))], [X(3, FACES(i, 1)) X(3, FACES(i, 3))], 'r');
plot3([X(1, FACES(i, 3)) X(1, FACES(i, 2))], [X(2, FACES(i, 3)) X(2, FACES(i, 2))], [X(3, FACES(i, 3)) X(3, FACES(i, 2))], 'r');
end; end;
view(80, -10);
view(80, -10);