définir des points dans une ellipse

**gegon** · 20/12/2012, 11h38

Bonjour,

J'ai défini une ellipse ainsi que des points dont je connais les coordonnées.
Je voudrais connaitre le nombre de points qui se trouve dans l'ellipse

définition de l'ellipse :

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
    N = 50;
    theta = 0:1/N:2*pi+1/N;
 
    % Parametric equation of the ellipse
    %----------------------------------------
    state(1,:) = arand/2*cos(theta); 
    state(2,:) = brand/2*sin(theta);
 
    % Coordinate transform (since your ellipse is axis aligned)
    %----------------------------------------
    X = state;
    X(1,:) = X(1,:) + C(1);
    X(2,:) = X(2,:) + C(2);
 
    % Plot
    %----------------------------------------
    plot(X(1,:),X(2,:));
     axis ( [ 0 25 0 25] )
    hold on;
   % plot(C(1),C(2),'r*');
    grid;
-----------------------------------------------------------

Les points sont définis par le vecteur P .

Je suis partie sur une boucle en définissant d'abord les point dont les x sont bien compris dans l'ellipse et je fais deux cas, si le point est dans la partie supérieure de l'ellipse ou inférieur.
Puis je compare l'ordonnée du point par rapport à l'ordonnée de l'ellipse pour la même abscisse

Pour le y de l'ellipse, il faut d'abord retrouver l'angle en utilisant la fonction arcos:

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
s=0;
for t=1
    %angle pour que le point dont l'abscisse est P(t,1) soit sur l'ellipse
    thetat=arccos(P(t,1)*2/brand);
    y= brand/2*sin(thetat);
    % si le l'absicsse du point est compris dnas le diametre de l'ellipse,
    % si l'ordonnée est dans la moitié supérieure alors on regarde si cette
    % ordonnée est superieure à l'ordonnée de l'ellipse pour le meme x
    if ( C(1)-arand/2 < P(t,1) < C(1)+ arand/2 & (0 < P(t,2) < C(2)+ brand/2)& (P(t,2)<=y) )
    elseif ( C(1)-arand/2 < P(t,1) < C(1)+ arand/2 & (C(2)- brand/2 < P(t,2) < 0)& (P(t,2)>=y) )
                s=s+1;
    end
end

C est le centre de l'ellipse et arand et brand les diamètres.

mais j'obtiens un message d'erreur suivant :

??? Undefined function or method 'arccos' for input arguments of type 'double'.

Est ce que quelqu'un pourrait m'aider à ce sujet en espérant avoir été compréhensible.

Merci d'avance à vous

**FLB** · 20/12/2012, 11h56

La fonction arc cosinus en Matlab se note acos

Code :

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

doc acos

Une autre solution pour savoir si un point est dans une ellipse est de passer par les coordonnées polaires du point et l'équation polaire de l'ellipse :

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
%Coordonnées du point par rapport au centre de l'ellipse
X = Xpoint-XcentreEllipse;
Y = Ypoint-YcentreEllipse;
%L'axe focale de ton ellipse est l'axe de tes coordonnées, donc pas besoin de tourner les points X Y
%Coordonnées polaires de ton point
[TH,R] = cart2pol(X,Y);
%Valeur du rayon de l'ellipse à cette angle - a et b sont les demi-grand et demi-petit axes
e = sqrt(a*a-b*b)/a;
rEllipse = sqrt(b*b/(1-e*e*cos(TH)*cos(TH)));
if (R<=rEllipse)
disp('Dedans!');
else
disp('Dehors!');
end

**gegon** · 20/12/2012, 15h58

Merci beaucoup pour cette réponse rapide.

Je voudrais tout d'abord essayer de finaliser mon idée et j'essaierai avec les coordonnées polaires juste après.

(Je suis désolée d'afficher le code ainsi mais la balise pour afficher le code ne fonctionne pas).

j'ai donc essayé ceci:

Code :

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

1
2
3
4
5
6
7
8
9
10
s=0;
for t=1:5
    %angle pour que le point dont l'abscisse est P(t,1) soit sur l'ellipse
    thetat=acos(P(t,1)*2/brand)*(2*pi); 
    y= brand/2*sin(thetat);
 
    if (( (C(1)-arand/2) < P(t,1) < (C(1)+ arand/2)) && (C(2) < P(t,2) < (C(2)+ brand/2))&& (P(t,2)<=y) )|| (( (C(1)-arand/2) < P(t,1) < (C(1)+ arand/2)) && ((C(2)- brand/2) < P(t,2) < C(2))&& (P(t,2)>=y) )
    s=s+1;
    end
end

Je m'explique : l'ellipse est centré autour de C la grand diamètre est arand et le petit diamètre est brand, je définis 2 cas , si le point est sous l'axe principale de l'ellipse ou au dessus, je met tout cela dans le même if, si l'une ou l'autre des conditions est avérée alors j’incrémente s.
Cependant, cela ne m'affiche jamais la bonne valeur. SI quelqu'un voit l'erreur...

Merci d'avance

**gegon** · 20/12/2012, 17h07

J'ai essayé avec les coordonnées polaires, ce qui semblent bien plus simple en y rajoutant l'incrémentation mais les résultats que j'obtiens sont comme dans l'autre méthode pas en accord avec ce que je vois sur la figure,

le s indique 5 alors qu'il n'y a que 4 points dans l'ellipse par exemple.

la balise de code ne fonctionne toujours pas :

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
s=0;
for t=1:5
%Coordonnées du point par rapport au centre de l'ellipse
X = P(t,1)-C(1);
Y = P(t,2)-C(2);
%L'axe focale de ton ellipse est l'axe de tes coordonnées, donc pas besoin de tourner les points X Y
%Coordonnées polaires de ton point
[TH,R] = cart2pol(X,Y);
%Valeur du rayon de l'ellipse à cette angle - arand et brand sont les demi-grand et demi-petit axes
f = sqrt(arand*arand-brand*brand)/arand;
rEllipse = sqrt(brand*brand/(1-f*f*cos(TH)*cos(TH)));
if (R<=rEllipse)
    s=s+1;
end
end
 
s

Merci de votre aide et bonne soirée!

**FLB** · 21/12/2012, 09h20

Tu peux nous donner tes valeurs?
J'ai fait un test et je trouve le bon résultat :

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
 
P = [0 0;1 1; 2 1; -1 0; 4 1];
C = [0.5 0.5]
arand = 2;
brand = 1;
figure;
rectangle('Curvature', [1 1],'Position', [C(1)-arand C(2)-brand 2*arand 2*brand]);
hold on; 
plot(P(:,1),P(:,2),'+');
axis equal;
 
%On s'attend à avoir s = 4
%je lance ton script :
s =
     4
%Ca marche!

**gegon** · 21/12/2012, 10h04

Oui il est vrai que ça marche mais pas tout le temps. Quand je le répète plusieurs fois, ca ne donne pas obligatoirement la bonne valeurs et je ne comprends pas pourquoi.
J'efface toujours les anciennes valeurs pour être sure que cela ne prenne pas en compte les anciennes.

et voilà tout mon code:

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
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
 
 
 
% Definition d'un noyau en forme d'ellipse a =grd diametre, b = petit diametre,
% exemple (a=17, b=13 et C = [10,10])
%function plotEllipse(a,b,C)
%function plotEllipse
a=17.64; aerror=4.77 ;
b=12.34; berror = 2.62;
C=[0,0];
e = 5; % coté du carré des points
r = 4; % erreur pour chaque particule
arand = a - aerror + 2*aerror*rand(1);
brand = b - berror + 2*berror*rand(1);
 
    % range to plot over
    %------------------------------------
    N = 50;
    theta = 0:1/N:2*pi+1/N;
 
    % Parametric equation of the ellipse
    %----------------------------------------
    state(1,:) = arand/2*cos(theta); 
    state(2,:) = brand/2*sin(theta);
 
    % Coordinate transform (since your ellipse is axis aligned)
    %----------------------------------------
    X = state;
    X(1,:) = X(1,:) + C(1);
    X(2,:) = X(2,:) + C(2);
 
    % Plot
    %----------------------------------------
    plot(X(1,:),X(2,:));
     axis ( [ -12.5 12.5 -12.5 12.5] )
    hold on;
   % plot(C(1),C(2),'r*');
    grid;
 
% placement de 5 points dans le noyau
%celui du milieu doit etre placé au centre de l'ellipse
thetar =rand (1,5)*(2*pi); 
XP = [];
YP = [];
for i =1:5
ang = thetar(i);
XP= [XP ;  C(1)+ r*cos(ang)];
YP= [YP ;  C(2)+ r*sin(ang)];
end
P = [XP YP];
 
P(1,1) = P(1,1);
P(1,2) = P(1,2);
 
P(2,1) =P(2,1)+ e/2;
P(2,2) =P(2,2)+ e/2;
 
P(3,1) =P(3,1)+ e/2;
P(3,2) =P(3,2)- e/2;
 
P(4,1) =P(4,1)- e/2;
P(4,2) =P(4,2)+ e/2;
 
P(5,1) =P(5,1)- e/2;
P(5,2) =P(5,2)- e/2;
 
 
plot(P(1,1),P(1,2),'r*');hold on;
plot(P(2,1),P(2,2),'r*');hold on;
plot(P(3,1),P(3,2),'r*');hold on;
plot(P(4,1),P(4,2),'r*');hold on;
plot(P(5,1),P(5,2),'r*');hold on;
 
%--------------------
%s=0;
%for t=1:5
%angle pour que le point dont l'abscisse est P(t,1) soit sur l'ellipse
% thetat=((acos(P(t,1)*2/brand)));
% y=( brand/2*sin(thetat));
% si le l'absicsse du point est compris dnas le diametre de l'ellipse,
% si l'ordonnée est dans la moitié supérieure alors on regarde si cette
% ordonnée est superieure à l'ordonnée de l'ellipse pour le meme x
 %if ((( (C(1)-arand/2) < P(t,1) < (C(1)+ arand/2)) && (C(2) < P(t,2) < (C(2)+ brand/2))&& (P(t,2)<=y) ) || (( (C(1)-arand/2) < P(t,1) < (C(1)+ arand/2)) && ((C(2)- brand/2) < P(t,2) < C(2))&& (P(t,2)>=y) ))
    %s=s+1;
   % else
    %s=s;
    %end
%end
%s
 
%-----------------------------------
 
s=0;
for t=1:5
%Coordonnées du point par rapport au centre de l'ellipse
X = P(t,1)-C(1);
Y = P(t,2)-C(2);
%L'axe focale de ton ellipse est l'axe de tes coordonnées, donc pas besoin de tourner les points X Y
%Coordonnées polaires de ton point
[TH,R] = cart2pol(X,Y);
%Valeur du rayon de l'ellipse à cette angle - a et b sont les demi-grand et demi-petit axes
f = sqrt(arand*arand-brand*brand)/arand;
rEllipse = sqrt(brand*brand/(1-f*f*cos(TH)*cos(TH)));
if (R<=rEllipse)
    s=s+1;
end
end
 
s

Merci pour tout

**FLB** · 21/12/2012, 10h17

Tu saurais trouver juste un set de points qui ne marche pas?
Au passage, tu peux vectoriser le code et te passer de la boucle for :

Code :

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

1
2
3
4
5
6
7
8
9
10
11
s=0;
%Coordonnées du point par rapport au centre de l'ellipse
X = P(:,1)-C(1);
Y = P(:,2)-C(2);
%L'axe focale de ton ellipse est l'axe de tes coordonnées, donc pas besoin de tourner les points X Y
%Coordonnées polaires de ton point
[TH,R] = cart2pol(X,Y);
%Valeur du rayon de l'ellipse à cette angle - arand et brand sont les demi-grand et demi-petit axes
f = sqrt(arand*arand-brand*brand)/arand;
rEllipse = sqrt(brand*brand./(1-f*f*cos(TH).*cos(TH)));
s = sum(R<=rEllipse);

**gegon** · 21/12/2012, 10h33

Oui j'ai voulu joindre une figure ou seul 3 points sont au sein de l'ellipse et où s vaut 5 mais toutes les balises pour ne fonctionne pas y compris pour joindre une image.

En lançant plusieurs fois le code et en effaçant bien au préalable le résultats précédent, j'obtiens bien à un moment une ellipse dont seuls 3 points sont dedans :

les coordonnées des points sont :

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
>> P
 
P =
 
    2.8258   -2.8311
    2.0164    6.4707
    3.4920   -6.3750
   -6.2267    1.0469
    1.1989   -0.9773
 
>> arand
 
arand =
 
   21.3769
 
>> brand
 
brand =
 
   11.4956

voilà

**FLB** · 21/12/2012, 11h12

Quand je trace l'ellipse (avec le bout de code précédent) je vois bien tous les points dans l'ellipse!
Ton ellipse est bien centrée en 0,0 ?
arand et brand correspondent aux demi grand axes et demi petit axes?
=> quand je prend arand et brand comme étant les grands axes alors en effet 2 points sont dehors !

Les points servent à faire une multiplication element à element entre 2 matrices, et non un produit matriciel:

Code :

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

1
2
3
4
5
6
7
8
>> [1 1;1 1].*[2 2;3 3]
ans =
     2     2
     3     3
>> [1 1;1 1]*[2 2;3 3]
ans =
     5     5
     5     5

**gegon** · 21/12/2012, 11h22

Alors arand est bien le grand dimatre et brand est bien le petit diametre, ce sont les grand axes (c'est pourquoi je les divise par 2 dans :

Code :

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

1
2
    state(1,:) = arand/2*cos(theta); 
    state(2,:) = brand/2*sin(theta);

je viens de refaire un exemple où s=1 :

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
>> s
 
s =
 
     1
 
>> arand
 
arand =
 
   15.5449
 
>> brand
 
brand =
 
   10.3370
 
>> P
 
P =
 
   -1.0303   -3.8650
   -1.0483    0.6534
    5.2809    0.3752
   -0.4615    5.9416
   -3.2016   -6.4380

la piece jointe fonctionne enfin donc voilà .

**FLB** · 21/12/2012, 11h39

Le problème vient donc de là :
dans le code

Code :

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

1
2
3
4
5
6
7
8
9
10
%Coordonnées du point par rapport au centre de l'ellipse
X = P(:,1)-C(1);
Y = P(:,2)-C(2);
%L'axe focale de ton ellipse est l'axe de tes coordonnées, donc pas besoin de tourner les points X Y
%Coordonnées polaires de ton point
[TH,R] = cart2pol(X,Y);
%Valeur du rayon de l'ellipse à cette angle - arand et brand sont les demi-grand et demi-petit axes
f = sqrt(arand*arand-brand*brand)/arand;
rEllipse = sqrt(brand*brand./(1-f*f*cos(TH).*cos(TH)));
sum(R<=rEllipse)

arand et brand correspondent aux demi diamètres, donc il faut faire la modification suivante :

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
 
demiGrandAxe = arand/2;
demiPetitAxe = brand/2;
%Coordonnées du point par rapport au centre de l'ellipse
X = P(:,1)-C(1);
Y = P(:,2)-C(2);
%L'axe focale de ton ellipse est l'axe de tes coordonnées, donc pas besoin de tourner les points X Y
%Coordonnées polaires de ton point
[TH,R] = cart2pol(X,Y);
%Valeur du rayon de l'ellipse à cette angle - arand et brand sont les demi-grand et demi-petit axes
f = sqrt(demiGrandAxe *demiGrandAxe -demiPetitAxe *demiPetitAxe )/demiGrandAxe ;
rEllipse = sqrt(demiPetitAxe*demiPetitAxe ./(1-f*f*cos(TH).*cos(TH)));
sum(R<=rEllipse)

**gegon** · 21/12/2012, 11h48

Super!

MErci beaucoup.

définir des points dans une ellipse [Débutant]

MATLAB

Discussions similaires

Partager

Partager