[bwarea/fill] Calcul de l'aire de disques

Version imprimable

Bonjour
Je suis une nouvelle utilisatrice de matlab (quelques semaines..). Et grace aux infos que j'ai trouvé sur les forums j'ai reussi à programmer quelque chose (merci) mais c'est assez lent…
J’ai des disques dont je connais les positions des centres et les rayons. Et je voudrais avoir une idée de l’air projetée.

J’utilise la function bwarea de matlab. Le code marche mais j’ai un pb de vitesse et de memoire (segmentation fault) si je fais trop de relances.

Voila mon script :
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 discount %nb de disques %X matrice avec les position des centres pour chaque disque discX = X(1,:); discY = X(2,:); discRadius=r; figure('Position',Pos1,'Visible', 'off') hold on for i=1:discCount c=discX(i); b=discY(i); circle(c,b,discRadius); end axis equal axis off hold off print -f -djpeg 'projection' b = imread('projection.jpg'); delete('projection.jpg') b = im2bw(b); %b=b;%(200:700,352:952); %crop the image intensity=bwarea(b);
et j’appelle circle qui est :
Code:

1 2 3 4 5 6 7 8 function [x y]=circle(a,b,r) theta = linspace(0,2*pi,100000); x = a+r*cos(theta); y = b+r*sin(theta); %plot(x,y); % plot circle fill(x,y,'k') axis('equal');% set equal scale on axes per pixel end
J’utilise la function bwarea de matlab. Mais du coup, il me faut une image binaire.
Est ce que c’est possible de sauvegarder (dans une matrix) les point du cercles qui sont remplie par le fill et de juste les compter (sans consruire une vraie image)?
Ou est ce qu’il y a une fonction comme fill qui donne une matrice au lieu d’une image?
Je ne sais pas si j’ai été assez claire, dites moi s’il vous faut plus d’infos.
Merci beaucoup!

14/08/2009, 17h44
Jerome Briot

Citation:

Envoyé par Nadge

J’ai des disques dont je connais les positions des centres et les rayons. Et je voudrais avoir une idée de l’air projetée.

Je ne comprend pas ce que tu appelles "aire projetée" ?

Pourquoi ne pas simplement calculer la surface des disques dès le départ ?
Pourquoi passer par une image ?

Pour la 1ère partie, c'est vectorisable (a priori ça gagne du temps mais je n'ai pas fait d'énormes tests) :

Code:

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
%discount %nb de disques
%X matrice avec les position des centres pour chaque disque
 
%Tests
% X = [1 2 3 4;
%     1 2 3 4];
% discCount = size(X,2);
% r = 2;
 
fig = figure('Position',Pos1,'Visible', 'off')
 
N = 100000;
theta = linspace(0,2*pi,100000);
x = repmat(X(1,:)',1,N)+repmat(r*cos(theta),discCount,1);
y = repmat(X(2,:)',1,N)+repmat(r*sin(theta),discCount,1);
fill(x',y','k')
axis('equal');% set equal scale on axes per pixel
axis equal
axis off
 
print -f -djpeg projection
 
%** Fermer la figure, car comme elle est invisible tu ne la vois plus mais elle prend des ressources
close(fig);
 
b = imread('projection.jpg');
delete('projection.jpg')
b = im2bw(b);
%b=b;%(200:700,352:952); %crop the image
intensity=bwarea(b);

Pour la suite je ne suis pas assez calée en images pour t'aider.

Bonjour.
Une idée avec polyarea :

Code:

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
clear
X = [1 2 3 4;
    1 2 3 4];
r=2;
N=1000;
Nc=4;%Nombre de cercles
x=zeros(1,N+1)';y=zeros(1,N+1)';
int=0;
for c=1:Nc-1
    for k=c+1:Nc
        int=int+1;
        d(int)=sqrt((X(1,k)-X(1,c))^2+(X(2,k)-X(2,c))^2);%Distance entre centres
        theta(int)=2*acos(d(int)/r/2);
        surf(int)=(r^2*theta(int)/2-r^2*sin(2*theta(int)));%Surfaces communes
        if isreal(surf(int))==0
            surf(int)=0;
        end
    end
end
for c=1:Nc
    L = linspace(0,2*pi,N);
    xv=r*cos(L)'+X(1,c);
    yv=r*sin(L)'+X(2,c);
    xv = [xv ; xv(1)]; yv = [yv ; yv(1)];
    x(:,c)=xv;y(:,c)=yv;
    A(c) = polyarea(xv,yv);
    plot(x(:,c),y(:,c)); title(['Aire = ' num2str(A)]); axis image
    hold on
end
%Surface projetée
S=sum(A)-sum(surf)

polyarea est utilisé pour les polygones, ici on peut évidemment mettre l'aire du cercle :

Code:

A(c)=pi*r^2;

Si les cercles sont projetés, il faut multiplier le résultat par le cosinus de l'angle des plans.

15/08/2009, 15h51
FR119492

Salut!

Citation:

J’ai des disques dont je connais les positions des centres et les rayons.

Dans R^2 ou dans R^3 ?

Citation:

Et je voudrais avoir une idée de l’air projetée.

Projetée sur quoi?
Jean-Marc Blanc
17/08/2009, 09h50
Nadge

Bonjour,

Merci beaucoup pour ces réponses!
Je vais essayer l'idée de Phryte, ça a l'air d'être exactement ce que je cherche!!

Ce que j'appellais aire projetée c'est juste l'aire totale des disques en fait.
Mais comme les disques se recouvrent en partie (distances entre centres < au rayon), je ne veux pas compter plusieurs fois les aires des intersections.

Donc je me suis mal exprimée, ce n'est pas une projection, je suis dans R^2, les disques sont tous dans le plan (sans inclinaisons).

Je suis passée par une image parce que j'avais trouvé la fonction bwarea de matlab qui prend une image et du coup je n'avais pas cherché mieux... Mais effectivement, je dois pouvoir éviter.

Caro-Line.
La fonction repmat c'est pour éviter de faire une boucle for et tracer tous les disques en meme temps?

Citation:

Envoyé par phryte

Code:

1
2
 
Une idée avec polyarea :

Le même code... optimisé à la sauce MATLAB :

Code:

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
X = [1 2 3 4;
    1 2 3 4];
 
 
r = 2;
N = 1000;
Nc = size(X,2); %Nombre de cercles
 
int = 0;
for c = 1:Nc-1
    for k = c+1:Nc
        int = int+1;
        d(int) = sqrt((X(1,k)-X(1,c))^2+(X(2,k)-X(2,c))^2);%Distance entre centres
    end
end
 
theta = 2*acos(d/r/2);
srf = (r^2*theta/2-r^2*sin(2*theta));%Surfaces communes
srf = real(srf);
 
L = linspace(0,2*pi,N).';
 
xv = bsxfun(@plus,r*cos(L),X(1,:));
yv = bsxfun(@plus,r*sin(L),X(2,:));
 
xv = [xv ; xv(1,:)];
yv = [yv ; yv(1,:)];
 
A = polyarea(xv,yv);
 
%Surface projetée
S = sum(A)-sum(srf);

;)

17/08/2009, 11h51
phryte

Bonjour et merci Dut. Plus élégant.
17/08/2009, 12h05
Nadge
Pour matlab, je pense que c'est bon.
Je crois que ce n'est plus qu'un problème de maths maintenant (pour le calcul de srf).
Par exemple quand plusieurs disques se recouvrent complètement (ils ont le même centre), j'obtiens 1.5 fois l'aire du disque.
Et quand il y a des intersections entre plusieurs disques, je peux obtenir des aires totales plus petites que l'aire d'un disque! Parce qu'on supprime plusieurs fois la même zone.

Par exemple avec :
Code:

1 2 3 X = [1 1 2; 1 2 2]; r = 2;

Citation:

Et quand il y a des intersections entre plusieurs disques, je peux obtenir des aires totales plus petites que l'aire d'un disque! Parce qu'on supprime plusieurs fois la même zone.

Effectivement cela n'est pas dans les mêmes conditions.

Une idée avec polyarea (à travailler et vérifier):
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 clear X = [1 1 2; 1 2 2]; r=2; N=100; Nc=size(X,2); %Nombre de cercles for c=1:Nc L = linspace(0,2*pi,N); xv=r*cos(L)'+X(1,c); yv=r*sin(L)'+X(2,c); xv = [xv ; xv(1)]; yv = [yv ; yv(1)]; x(:,c)=xv;y(:,c)=yv; end [lat,lon] = POLYBOOL('union',x(:,1),y(:,1),x(:,2),y(:,2));%Cercles 1 et 2 [lat,lon] = POLYBOOL('union',lat,lon,x(:,3),y(:,3));%Intersection avec cercle 3 plot(lon,lat) grid surface_union=polyarea(lon,lat)

17/08/2009, 14h21
Jerome Briot

Citation:

Envoyé par Nadge

Par exemple quand plusieurs disques se recouvrent complètement (ils ont le même centre), j'obtiens 1.5 fois l'aire du disque.

Ce cas particuliers peut être traité en amont du code en comparant les valeurs dans X et en ne conservant qu'un seul des disques qui se chevauchent dans X
17/08/2009, 16h51
Nadge

Citation:

Ce cas particuliers peut être traité en amont du code en comparant les valeurs dans X et en ne conservant qu'un seul des disques qui se chevauchent dans X

Merci Dut!
J'ai fait ça avec

Code:

X=(unique(X','rows'))'

et ça marche! Un problème de moins!!

Il me reste plus que le pb des "intersections multiples".
Phryte, j'ai essayé ton script mais j'ai un message d'erreur qui apparait.

Code:

??? Undefined function or method 'polybool' for input arguments of type 'char'.

Je suis assez vite perdue avec les initialisations de variables. Et je ne sais pas laquelle devrait être d'un autre type que 'char' et pourquoi.
Et polybool n'apparait pas dans l'aide.

Qu'est ce que c'est exactement? Une fonction qui construit un polygone qui englobe les cercles? Est ce que ça s'écrit bien comme ça? ou est ce que c'est une fonction qu'il faut que j'écrive??
Merci!
17/08/2009, 16h59
Jerome Briot

Citation:

Envoyé par Nadge

Et polybool n'apparait pas dans l'aide.

Il semble que ce soit une fonction de la Mapping Toolbox : http://www.mathworks.com/access/help.../polybool.html
17/08/2009, 17h15
Nadge

Merci! Je n'avais même pas penser à chercher sur internet...
Je vais demander si je peux avoir cette toolbox.
17/08/2009, 17h19
Jerome Briot

Regarde ceci : Polygon Clipper ou ceci : Polygon_Intersection
18/08/2009, 18h19
Nadge

Bonjour,
J'ai utilisé "Polygon Clipper" pour former les nouveaux contours avant de calculer l'aire. Et ça a l'air de bien marcher! Merci beaucoup!!