Discussion :

[Giansolo]

Bonjour à toutes et à tous,

J'ai un ensemble d'enveloppes convexes qui se chevauchent sur une figure matlab et je souhaiterais calculer l'aire que représente l'union de ces enveloppes convexes. il peut y avoir plusieurs "patches" éloignées regroupant des ensemble d'enveloppes convexes, ce qui complique la chose.

Un exemple en images (je cherche l'aire couverte par la zone noire totale) dans les fichiers attachés à ce topic.

J'ai pensé aux méthodes d'imagerie: patcher les enveloppes convexes (comme dans l'image2) et ensuite estimer l'aire, problème: lent si beaucoup d'images différentes (ce qui est mon cas) et qui revient à une estimation dépendante de la résolution de l'image.

Existe-t'il des solutions plus simples et moins couteuses en temps de calcul? La triangulation de Delaunay peut-elle m'être utile dans ce cas?

Merci par avance,

Gian

Edit: je possède l'aire de chaque enveloppe convexe, et j'ai aussi pensé a la fonction polyarea() qui pourrait peut-être m'aider?

[Giansolo]

Pour polybool() et polyarea(), un exemple simple montre que cette fonction n'est pas adaptée à mon problème (un polygone distant renvoie 'NaN' :-/):

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
 
%définition de 3 carrés:
x = [0,2; 0,0; 2,0 ;2,2; 0,2];
y = x + 1;
z = x + 5;
 
[x(:,1), x(:,2)] = poly2cw(x(:,1),x(:,2));
[y(:,1), y(:,2)] = poly2cw(y(:,1),y(:,2));
[z(:,1), z(:,2)] = poly2cw(z(:,1),z(:,2));
 
%Unions:
[X,Y] = polybool('union',x(:,1),x(:,2),y(:,1),y(:,2));
[X2,Y2] = polybool('union',X,Y,z(:,1),z(:,2));
s = polyarea(X2,Y2);
 
plot(x(:,1),x(:,2),'b-'); hold on;
plot(y(:,1),y(:,2),'r-');
plot(z(:,1),z(:,2),'r-');
axis equal;

[Jean Dumoncel]

Bonjour,

Après union tu auras 2 polygones disjoints. Pour différencier ces 2 polygones, polybool les sépare avec des NaN dans ta liste (X2,Y2). Il faut que tu appliques polyarea à chaque liste séparée par des NaN.

[Giansolo]

Ok, merci pour l'info, j'essaye ca!

Je suis également en train d'implémenter une autre solution en partant sur l'idée des algorithmes d'imagerie et ca donne ca:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
 
%definir un carré:
x = [0,2; 0,0; 2,0 ;2,2; 0,2];
y = x + 1;
z = x + 5;
 
[x(:,1), x(:,2)] = poly2cw(x(:,1),x(:,2));
[y(:,1), y(:,2)] = poly2cw(y(:,1),y(:,2));
[z(:,1), z(:,2)] = poly2cw(z(:,1),z(:,2));
 
%Unions:
[X,Y] = polybool('union',x(:,1),x(:,2),y(:,1),y(:,2));
[X2,Y2] = polybool('union',X,Y,z(:,1),z(:,2));
s = polyarea(X2,Y2);
 
h = figure;
plot(x(:,1),x(:,2),'b-'); hold on;
plot(y(:,1),y(:,2),'r-');
plot(z(:,1),z(:,2),'r-');
axis equal;
%patch:
patch(x(:,1),x(:,2),'k')
patch(y(:,1),y(:,2),'k')
patch(z(:,1),z(:,2),'k')
 
 
%Export (sans les axes):
saveas(h, 'toto.jpg', 'jpg');
close
%Reimport:
BW = imread('toto.jpg','jpg');
 
%Binarise:
ZZ = BW > 255/2;
 
STATS = regionprops(double(ZZ(:,:,1)),'Area');

Seul problème, à quoi correspond l'aire? c'est l'aire des pixels noirs je suppose il faut donc que je fasse une opération vis à vis de la résolution et que j'estime la taille d'un pixel? (il faut aussi que j'exporte l'image correctement sans les axes).

Gian

[Jean Dumoncel]

Non, pour regionprops, une région c'est une zone de pixel blanc. Dans ton exemple, il faut donc inverser les pixels blancs/noirs.

Sinon, tu peux éviter la sauvegarde de la figure en jpg :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
%definir un carré:
x = [0,2; 0,0; 2,0 ;2,2; 0,2];
y = x + 1;
z = x + 5;

[x(:,1), x(:,2)] = poly2cw(x(:,1),x(:,2));
[y(:,1), y(:,2)] = poly2cw(y(:,1),y(:,2));
[z(:,1), z(:,2)] = poly2cw(z(:,1),z(:,2));

%Unions:
[X,Y] = polybool('union',x(:,1),x(:,2),y(:,1),y(:,2));
[X2,Y2] = polybool('union',X,Y,z(:,1),z(:,2));
s = polyarea(X2,Y2);

h = figure('Color',[1 1 1]);
plot(x(:,1),x(:,2),'b-'); hold on;
plot(y(:,1),y(:,2),'r-');
plot(z(:,1),z(:,2),'r-');
axis equal;
%patch:
patch(x(:,1),x(:,2),'k')
patch(y(:,1),y(:,2),'k')
patch(z(:,1),z(:,2),'k')

axis off

f = getframe;

ZZ = rgb2gray(f.cdata) > 255/2;

Puis le calcul de l'aire :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

STATS = regionprops(~ZZ(:,:,1),'Area');

tu peux aussi éviter l'inversion en changeant la couleur de tes patchs en blancs et changeant le fond de la figure en noir.

[Giansolo]

Merci pour cette correction, c'est en effet plus propre, cependant tel quel ça ne fonctionne pas (regionprops veut une 'label matrix'). J'ai essayé de faire un cast en double():

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
ZZ = double(rgb2gray(f.cdata) > 255/2);
STATS = regionprops(~ZZ,'Area');

puis en utilisant bwlabel:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
L = bwlabel(ZZ,8);
STATS = regionprops(~L,'Area');

mais regionprops() ne veut rien savoir.

Ensuite question idiote, mais comment récupère-ton le vecteur unité de la figure pour connaitre la surface que représente un pixel?

Merci encore pour ton aide.

[Jean Dumoncel]

Ben si, normalement tel quel le code aurait dû fonctionner, regionprops accepte les tableaux de valeurs logiques, et ZZ est un tableau de valeurs logiques :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
ZZ = rgb2gray(f.cdata) > 255/2;
STATS = regionprops(~ZZ,'Area');

Quelle version de matlab as-tu?


Pour la résolution, l'image jpg aurait peut-être permis d'aller plus rapidement. Tu peux avoir les intervalles des axes de ton image :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

inter = axis;

puis tu divises les intervalles par la taille de l'image :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
dx = (inter(2)-inter(1))/size(ZZ,1);
dy = (inter(4)-inter(3))/size(ZZ,2);

Avec le code du dessus, normalement, dx=dy.

[Invité]

Bonjour,

S'il demande une 'label matrix', alors je dirais

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
ZZ = rgb2gray(f.cdata) > 255/2;
STATS = regionprops(bwlabel(~ZZ),'Area');

[Giansolo]

Avec cette dernière version ça marche correctement (en appliquant le bwlabel()), mais:
-j'obtiens 2 aires par regionprops()??
-mon estimation finale est très loin du résultat (j’obtiens une surface de 7 contre 11 en vrai), a quoi est du cette mauvaise estimation? faut-il augmenter la résolution obtenue avec getframe() ? Si oui, comment procède t-on?

Je mets le code complet:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
 
%definir trois carrés:
x = [0,2; 0,0; 2,0 ;2,2; 0,2];
y = x + 1;
z = x + 5;
 
h = figure('Color',[1 1 1]);
plot(x(:,1),x(:,2),'b-'); hold on;
plot(y(:,1),y(:,2),'r-');
plot(z(:,1),z(:,2),'r-');
axis equal;
%axis square;
%patch:
patch(x(:,1),x(:,2),'k')
patch(y(:,1),y(:,2),'k')
patch(z(:,1),z(:,2),'k')
 
axis off;
 
f = getframe;
 
%Binarise:
ZZ = rgb2gray(f.cdata) > 255/2;
 
STATS = regionprops(bwlabel(~ZZ),'Area');
aire = STATS.Area;
 
inter = axis;
dx = (inter(2)-inter(1))/size(ZZ,1);
dy = (inter(4)-inter(3))/size(ZZ,2);
 
airenoire = aire*(dx*dy)

Edit: Ah, j'ai une piste, en faisant un cast en double(), plutot qu'en passant par bwlabel():

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
 
STATS = regionprops(double(~ZZ),'Area');

j'obtiens une estimation plus correcte (=11.1368) et regionprops() ne retourne qu'une valeur d'aire. D'autre part, en jouant sur les axes, et la résolution on influ sur la solution (mais on le savait déjà! ). Toutefois je ne m'explique pas un telle différence avec bwlabel()... si vuos avez des idées pour optimiser la solution, je suis preneur!
Edit2: Pour les deux areas de bwlabel() j'ai compris, il faut les sommer, ca correspond aux aires "disjointes"