Discussion :

[Nadge]

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 :

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
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 :

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

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!

[Jerome Briot]

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 ?

[Caro-Line]

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

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
%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.

[phryte]

Bonjour.
Une idée avec polyarea :

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
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 :

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

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.

[FR119492]

Salut!

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

Dans R^2 ou dans R^3 ?

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

Projetée sur quoi?
Jean-Marc Blanc

[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?

[Jerome Briot]

Code :

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

1
2
 
Une idée avec polyarea :

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

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
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);

[phryte]

Bonjour et merci Dut. Plus élégant.

[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 :

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

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

[phryte]

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 :

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
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)

[Jerome Briot]

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

[Nadge]

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 :

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

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 :

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

 ??? 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!

[Jerome Briot]

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

[Nadge]

Merci! Je n'avais même pas penser à chercher sur internet...
Je vais demander si je peux avoir cette toolbox.

[Jerome Briot]

Regarde ceci : Polygon Clipper ou ceci : Polygon_Intersection

[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!!