Discussion :

[Azzo2]

Bonjour à tous,

J'essaie de déterminer le centre d'une spirale logarithmique en sélectionnant 3 points (pas trop loin les uns des autres) à un endroit aléatoire de ma spirale. Je trace la tangente à tous les points de ma spirale. Je connais le coefficient directeur des 3 tan aux 3 points choisis et je veux trouver: un angle de rotation, xc et yc (coordonnées du centre de la spirale). Pour ça, je résous un système de 3 équations (3 produits scalaires). Noter que l'angle entre un rayon vecteur de la spirale et une tangente en un point est toujours le même.

x_square correspond aux abscisses de la spirale, y_square aux ordonnées de la spirale, true_tan1 2 et 3 les tan aux points P1 P2 et P3, t1, t2 et t3 sont les vecteurs directeurs des tan, r1 r2 et r3 les vecteurs directeurs des rayons aux points P1 P2 et P3.

Parfois le système se résout parfaitement: l'erreur entre le vrai centre et celui déterminé est de l'ordre de 10e-4.
D'autre fois l'erreur est de l'ordre de 10e26 et je ne comprends pas pourquoi.
De plus, l'angle déterminé varie entre 270 degrés (assez rare) et 3000-4000 degrés (voir plus). Pourquoi ?

Ci-dessous le code:

Code MATLAB :

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
 
syms xc yc ang
 
t1 = [x_square(2)-x_square(1) true_tan1(2)-true_tan1(1) 0];
t2 = [x_square(2)-x_square(1) true_tan2(2)-true_tan2(1) 0];
t3 = [x_square(2)-x_square(1) true_tan3(2)-true_tan3(1) 0];
 
r1 = [xc-x_square(P1) yc-y_square(P1) 0];
r2 = [xc-x_square(P2) yc-y_square(P2) 0];
r3 = [xc-x_square(P3) yc-y_square(P3) 0];
 
scal1 = t1.*r1;
scal2 = t2.*r2;
scal3 = t3.*r3;
 
eq = [(scal1(1)+scal1(2))/(sqrt((t1(1))^2+(t1(2))^2)*sqrt((r1(1))^2+(r1(2))^2)) == cos(ang), (scal2(1)+scal2(2))/(sqrt((t2(1))^2+(t2(2))^2)*sqrt((r2(1))^2+(r2(2))^2)) == cos(ang), (scal3(1)+scal3(2))/(sqrt((t3(1))^2+(t3(2))^2)*sqrt((r3(1))^2+(r3(2))^2)) == cos(ang)];
centre = vpasolve(eq, [xc yc ang]);

Noter également que les tan sont calculées comme suit:

Code MATLAB :

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

1
2
3
4
5
6
7
 
dy = diff(y_square);
dx = diff(x_square); 
 
true_tan1 = (dy(P1)*(x_square-x_square(P1)))/dx(P1) + y_square(P1);
true_tan2 = (dy(P2)*(x_square-x_square(P2)))/dx(P2) + y_square(P2);
true_tan3 = (dy(P3)*(x_square-x_square(P3)))/dx(P3) + y_square(P3);

Merci d'avance pour votre aide,

Azzo2.

[tbc92]

A partir de 3 points, il y a (en général) beaucoup de solutions (une infinité ?)

Exemple : à partir des points (10,1),(10,2),(9,5), il y a la spirale 'naturelle' avec le centre près de (0,0) et qui passe par nos 3 points en 1 dizième de tour. Je ne sais pas s'il y a une seule ou plusieurs solutions ici, Notons cette famille la famille 1, avec éventuellement un seul élément.
On a aussi plein de spirales dont le centre va être vers (10, -100) et qui vont passer par le 1er point, puis après un tour complet, par le 2ème point, et après 3 ou 4 tours complets, par le 3ème point. Famille 2.
On peut aussi trouver plein de spirales avec le centre vers (10, 100) Famille 3, et d'autres avec le centre vers (100, 100) ou bien au delà , famille 4.

Je ne comprends rien au code que tu as posté, mais à mon avis, ton code n'a aucune raison de choisir un centre de la famille 1 plutôt que 2 ou 3 ou 4.

[Azzo2]

Sauf que mon code utilise des tangentes donc le passage de la spirale aux points 1 2 et 3 est "imposé" donc la solution aussi.

Pour mieux comprendre je poste l'intégralité du code, un petit copier coller et F5 dans Matlab permettra de mieux comprendre.

Code MATLAB :

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
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
% ***************************************************************************
% spiral and square
% ***************************************************************************
 
a=1;
b=0.05;
n=1;t=0;l=0;
X = randi([-10 10],1,1);
Y = randi([-10 10],1,1);
x(n)=a*exp(b*t)*cos(t);
y(n)=a*exp(b*t)*sin(t);
long = 50; % 50 is good
while t<long 
    n=n+1;
    while l < 0.00001 % ~ distance between two points
        xx=a*exp(b*t)*cos(t);
        yy=a*exp(b*t)*sin(t);
        l=sqrt((xx-x(n-1))^2+(yy-y(n-1))^2);
        t=t+1e-5; % reduire t pour augmenter le nombre de points (distance entre 2 pts diminuera)
    end
    x(n)=a*exp(b*t)*cos(t);
    y(n)=a*exp(b*t)*sin(t);
    d(n)=sqrt((x(n)-x(n-1))^2+(y(n)-y(n-1))^2); % distance entre points
    l=0;
    dy(n) = a*b*exp(b*t)*sin(t)+a*exp(b*t)*cos(t);
    dx(n) = a*b*exp(b*t)*cos(t)-a*exp(b*t)*sin(t);
end
 
extremite = sqrt((x(n)).^2+(y(n)).^2);
square_length = 5;
 
% ***************************************************************************
% conserve only points into the square
% ***************************************************************************
 
x_square = zeros(1,2);
y_square = zeros(1,2);
t = 1;
for k = 1:length(x)
    if X < x(k) && x(k) < X+square_length && Y < y(k) && y(k) < Y+square_length
        x_square(t) = x(k);
        y_square(t) = y(k);
        t = t + 1;
    end
end
 
% ***************************************************************************
% points selection          RED = 1         BLUE = 2        BLACK = 3
% ***************************************************************************
 
% central point
distance = sqrt(2*(square_length^2));
point_number = 0;
for k = 1:t-1
    if x_square(k) < X+2.5 && y_square(k) < Y+2.5
        di = sqrt((X+2.5-x_square(k)).^2+(Y+2.5-y_square(k)).^2);
    else if x_square(k) > X+2.5 && y_square(k) > Y+2.5
        di = sqrt((x_square(k)-X-2.5).^2+(y_square(k)-Y-2.5).^2);
    else if x_square(k) < X+2.5 && y_square(k) > Y+2.5
        di = sqrt((X+2.5-x_square(k)).^2+(y_square(k)-Y-2.5).^2);
    else if x_square(k) > X+2.5 && y_square(k) < Y+2.5
        di = sqrt((x_square(k)-X-2.5).^2+(Y+2.5-y_square(k)).^2);
        end
        end
        end
    end
    if di < distance
        distance = di;
        point_number = k;
    end
end
 
% tangent points
eloignement = 100;
P1 = point_number + eloignement;
P2 = point_number - eloignement;
P3 = point_number;
 
% ***************************************************************************
% tangents          RED = 1         BLUE = 2        BLACK = 3
% ***************************************************************************
 
dy = diff(y_square);
dx = diff(x_square);    
 
true_tan1 = (dy(P1)*(x_square-x_square(P1)))/dx(P1) + y_square(P1);
true_tan2 = (dy(P2)*(x_square-x_square(P2)))/dx(P2) + y_square(P2);
true_tan3 = (dy(P3)*(x_square-x_square(P3)))/dx(P3) + y_square(P3);
 
coeff1 = (true_tan1(2)-true_tan1(1))/(x_square(2)-x_square(1));
xtan1 = -extremite:0.1:extremite;
ytan1 = coeff1*xtan1 + y_square(P1) - coeff1*x_square(P1);
 
coeff2 = (true_tan2(2)-true_tan2(1))/(x_square(2)-x_square(1));
xtan2 = -extremite:0.1:extremite;
ytan2 = coeff2*xtan2 + y_square(P2) - coeff2*x_square(P2);
 
coeff3 = (true_tan3(2)-true_tan3(1))/(x_square(2)-x_square(1));
xtan3 = -extremite:0.1:extremite;
ytan3 = coeff3*xtan3 + y_square(P3) - coeff3*x_square(P3);
 
% ***************************************************************************
% find centre
% ***************************************************************************
 
syms xc yc ang
 
t1 = [x_square(2)-x_square(1) true_tan1(2)-true_tan1(1) 0];
t2 = [x_square(2)-x_square(1) true_tan2(2)-true_tan2(1) 0];
t3 = [x_square(2)-x_square(1) true_tan3(2)-true_tan3(1) 0];
 
r1 = [xc-x_square(P1) yc-y_square(P1) 0];
r2 = [xc-x_square(P2) yc-y_square(P2) 0];
r3 = [xc-x_square(P3) yc-y_square(P3) 0];
 
scal1 = t1.*r1;
scal2 = t2.*r2;
scal3 = t3.*r3;
 
eq = [(scal1(1)+scal1(2))/(sqrt((t1(1))^2+(t1(2))^2)*sqrt((r1(1))^2+(r1(2))^2)) == cos(ang), (scal2(1)+scal2(2))/(sqrt((t2(1))^2+(t2(2))^2)*sqrt((r2(1))^2+(r2(2))^2)) == cos(ang), (scal3(1)+scal3(2))/(sqrt((t3(1))^2+(t3(2))^2)*sqrt((r3(1))^2+(r3(2))^2)) == cos(ang)];
centre = vpasolve(eq, [xc yc ang]);
 
xc = centre.xc;
yc = centre.yc;
ang = centre.ang;
 
r1 = [xc-x_square(P1) yc-y_square(P1) 0];
r2 = [xc-x_square(P2) yc-y_square(P2) 0];
r3 = [xc-x_square(P3) yc-y_square(P3) 0];
 
figure;
axis equal;
hold on
plot(x,y,'k');
plot([X X+square_length X+square_length X X],[Y Y Y+square_length Y+square_length Y],'g');
plot(0,0,'Xk');
plot(x_square,y_square,'.g')
plot(x_square(P1),y_square(P1),'Or')
plot(x_square(P2),y_square(P2),'Ob')
plot(x_square(P3),y_square(P3),'O','color',[0.5,0.7,0.7])
% plot tan
plot(xtan1,xtan1*(t1(2)/t1(1))+y_square(P1)-(t1(2)/t1(1))*x_square(P1),'r');
plot(xtan1,xtan1*(t2(2)/t2(1))+y_square(P2)-(t2(2)/t2(1))*x_square(P2),'b');
plot(xtan1,xtan1*(t3(2)/t3(1))+y_square(P3)-(t3(2)/t3(1))*x_square(P3),'color',[0.5,0.7,0.7]);
% plot radius
plot(xtan1,xtan1*(r1(2)/r1(1))+y_square(P1)-(r1(2)/r1(1))*x_square(P1),'r');
plot(xtan1,xtan1*(r2(2)/r2(1))+y_square(P2)-(r2(2)/r2(1))*x_square(P2),'b');
plot(xtan1,xtan1*(r3(2)/r3(1))+y_square(P3)-(r3(2)/r3(1))*x_square(P3),'color',[0.5,0.7,0.7]);
 
plot(xc,yc,'x-','color',[0.5,0.7,0.7]);
title('Error');
hold off
 
['xc' xc 
 'yc' yc
 'ang' ang*180/pi
 'erreur' sqrt((xc)^2+(yc)^2)]

Si le programme est lancé plusieurs fois, tu t'apercevras que l'erreur varie (parfois énormément).

[tbc92]

1. Je n'ai pas Matlab, donc je ne peux pas tester.
2. même avec le code complet, ce n'est pas clair.
3. Tu dis que tu utilises la fonction tangente etc etc ... et donc tout devrait aller comme il faut.

Avec les fonctions sinusoïdales, il y a un autre problème qui vient se greffer.
Si tes 3 points sont ceux que je donnais au départ (10,1),(10,2) et (9,5) , et que le vrai centre est vers (0,0), alors ok, les 3 angles valent 0°, 10° et 25° , et la fonction tangente est continue sur l'intervalle en question.

Mais si tes 3 points sont (1,10), (-1,10) et (-2,9) , le centre de la spirale est vers (0,0) là encore. Les 3 angles valent 85°, 95° et 115°... et là, patatras, la fonction tangente n'est pas continue sur cet intervalle. Et tu auras le même probème avec des fonctions Arcsinus ou que sais-je encore.

Peut-être que les jeux de données qui font planter ton programme sont dans cette configuration ?

[Azzo2]

Je parle de l'angle PSI dans le lien ci-dessous, pas de l'angle de mon point avec une quelconque référence. Cet angle est toujours le même.

https://www.mathcurve.com/courbes2d/...arithmic.shtml

[wiwaxia]

Bonjour,

... x_square correspond aux abscisses de la spirale, y_square aux ordonnées de la spirale, true_tan1 2 et 3 les tan aux points P1 P2 et P3, t1, t2 et t3 sont les vecteurs directeurs des tan, r1 r2 et r3 les vecteurs directeurs des rayons aux points P1 P2 et P3.

Parfois le système se résout parfaitement: l'erreur entre le vrai centre et celui déterminé est de l'ordre de 10e-4.
D'autre fois l'erreur est de l'ordre de 10e26 et je ne comprends pas pourquoi.
De plus, l'angle déterminé varie entre 270 degrés (assez rare) et 3000-4000 degrés (voir plus). Pourquoi ? ...

1°) Si j'ai bien compris, les données sont:
a) les positions des 3 points envisagés (P1, P2, P3);
b) l'orientation des tangentes au graphe en chacun de ces points
Sait-on si la spirale est croissante ou décroissante ?

2°) Le problème devrait normalement se ramener au point de concours (approché ?) de trois droites, donc à un problème de triangulation
a) Que se passe-t-il dans les cas de tangentes quasi-parallèles ? Le centre risque d'être rejeté à l'infini ... et le nombre de tours deviendra sans doute lui aussi très grand.
b) Cela s'apparente d'une certaine manière à la recherche du centre d'un cercle passant par 3 points; ne pourrait-on envisager l'intervention d'une méthode des moindres carrés ?

[Azzo2]

Bonjour
1) a) oui b) oui.
On ne sait pas si la spirale est croissante ou décroissante car on regarde seulement une partie en supposant qu'on ne connait pas le reste.

2) a) Dans le cas de tangentes quasi parallèles les résultats sont quasi similaires (erreur de 10e-4, 10e-3 parfois; et erreur de 10e27 à 10e30 à d'autres moments). Néanmoins j'évite ce problème car "l'éloignement" des points P1, P2 et P3 peut être fixé.
b) Oui, je l'ai envisagé mais sans vraiment réussir à le mettre en forme.

[wiwaxia]

Il y a curieusement une solution analytique assez simple.

Soient r₁ = CA₁ , r₂ = CA₂ les positions des deux premiers points,
t₁ et t₂ les vecteurs unitaires tangents correspondants.

Les deux premiers déterminent l'angle orienté

(r₁, r₂) = (r₁, t₁) + (t₁, t₂) + (t₂, r₂) ;

or dans le cas d'une spirale logarithmique, les rayons vecteurs font avec la tangente orientée un angle constant en relation directe avec son équation polaire:

(r₁, t₁) = (r₂, t₂) = a ;

on en tire:

(r₁, r₂) = a + (t₁, t₂) - a = (t₁, t₂) .

Le centre (C) appartient donc au cercle (C₁₂) sous lequel on voit (A₁A₂) sous l'angle constant b₁₂ = (t₁, t₂);
ce cercle, centré sur la médiatrice du segment, admet pour rayon R₁₂ = A₁A₂/(2.Sin(b₁₂)) .

Le point recherché se situe donc à l'intersection commune de trois cercles, dont la détermination devient problématique dans le cas de tangentes quasi-parallèles (b_ij ~ 0 ou ±Pi) .

Mais peut-être pourrais-tu suivre cette piste: exprimer pour un point quelconque (M) les 3 produits scalaires
q_i = (MA_i│t_i) = (MA_i.u_i│t_i) = r_i*Cos(u_i, t_i) = r_i*Cos(a_i).

Sachant que dans le cas où (M) vient à se confondre avec (C), on obtient a₁ = a₂ = a₃ (= a), il y a une propriété à exploiter par le calcul d'une moyenne et d'un écart-type ? ... Je n'ai pas le temps de poursuivre.

[Azzo2]

N'est-ce pas ce que j'ai fait ligne 121 ?

Je résous 3 produits scalaires avec xc et yc les inconnues de mon point (M par exemple) et ang l'angle dont il est question. J'égalise les 3 produits scalaires et j'ai mes solutions.

A moins d'avoir manqué une étape de votre raisonnement... ?

[wiwaxia]

Je voulais parler d'un procédé statistique, qui évite si possible la manipulation de la fonction Arccos(x) . La résolution des équations cartésiennes conduit régulièrement au piège des points très (ou infiniment) éloignés.
Il faut si possible dégager une recherche des moindres carrés portant sur des distances.

On a pour tout point (M) autre que le centre (C):

(r_i, r_j) ≠ (t_i, t_j) ,

soit en passant au cosinus, puis en multipliant chaque terme par le produit des distances (r_i*r_j) = (MA_i*MA_j):

(r_i*r_j)*Cos(r_i, r_j) ≠ (r_i*r_j)*Cos(t_i, t_j) .

Le terme de gauche n'est autre que le produit scalaire (MA_i|MA_j) , et si l'on note (u_{i, j}) l'angle (t_i, t_j), il vient:

(MA_i|MA_j) ≠ (r_i*r_j)*Cos(u_{i, j}) .

Il découle tout naturellement de ceci que la somme des carrés (développée pour les 3 couples d'entiers 1, 2, 3)

S = ∑_j≠i((MA_i|MA_j) - (r_i*r_j)*Cos(u_{i, j}))²

s'annule lorsque (M) se confond avec (C) et seulement avec celui-ci.

La localisation du centre (C) se réduit ainsi à la recherche du minimum nul de la fonction S(M).
Si les termes sont aisément dérivables, on peut employer la méthode du gradient.

[Michel]

Bonjour à tous,
je n'ai pas la solution mais si on examine les 3 tangentes, elles forment un angle constant avec le rayon vecteur. Cette spirale est parfois appelée équiangulaire pour cette raison. (C'est peut-être une piste...)
Par contre, à première vue, on ne peut rien dire sur l'intersection des 3 perpendiculaires au 3 tangentes comme on peut le voir sur ce petit exemple fait avec geogebra

[wiwaxia]

J'ai repris par curiosité la fonction du point S(M) décrite au message #10, de laquelle je n'avais qu'une vague idée du comportement.

... Il découle tout naturellement de ceci que la somme des carrés (développée pour les 3 couples d'entiers 1, 2, 3)

S = ∑_j≠i((MA_i|MA_j) - (r_i*r_j)*Cos(u_{i, j}))²

s'annule lorsque (M) se confond avec (C) et seulement avec celui-ci.

La localisation du centre (C) se réduit ainsi à la recherche du minimum nul de la fonction S(M) ...

Il importe en effet de s'assurer de la présence d'un minimum nul (et d'un seul) dans une région contenant les 3 points de la spirale recherchée.
La fonction étant homogène à la quatrième puissance d'une distance, il est probable que l'on obtienne plus d'informations sur le relief de la surface étudiée (z = F(x, y)) en considérant l'équation F(x, y) = S(M)^1/4 .

Voici ce que l'on obtient, par exemple, à partir de 3 points donnés associés à un vecteur unitaire tangent d'inclinaison (u):

Code :

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

1
2
3
4
        P1        P2        P3    
Xi      9.50     13.0      13.0
Yi      6.50     10.0      17.0
Ui(°)  21.20     66.20    111.20

1°) sur le domaine [-10 ; +30]×[-10 ; +30], en envisageant les 3 fonctions:

F₁(x, y) = S(M) , F₂(x, y) = S(M)^1/2 , F₃(x, y) = S(M)^1/4 :

On constate que la première image ne livre aucune information sur le relief de la partie centrale; les détails n'apparaissent que sur la suivante, et surtout la dernière.

2°) sur un domaine plus étroit [0 ; +20]×[0 ; +20], en reprenant les mêmes fonctions:

F₁(x, y) = S(M) , F₂(x, y) = S(M)^1/2 , F₃(x, y) = S(M)^1/4 :

Les détails y présentent désormais une finesse comparable; il apparaît dans tous les cas une cuvette dont le point le plus bas représente le centre de la spirale, et qui se situe à proximité d'un maximum secondaire.

On relève très simplement, lors de l'exécution du programme de synthèse de l'image, les coordonnées du point recherché:

La précision, évidemment médiocre (~ 1 %), est liée à la résolution de l'image, et à ses dimensions; elle est de l'ordre de 1/(La*Ha)^1/2 .
Le fait essentiel est l'existence du minimum d'ordonnée nulle; un procédé itératif donnera donc rapidement la position du centre, avec une précision beaucoup plus grande.

Matlab présente probablement des procédures appropriées pour la recherche envisagée.

[wiwaxia]

Quelque chose ne va pas dans cet énoncé: le tracé de la spirale passe systématiquement à côté des points sur laquelle la courbe est théoriquement construite. Les écarts observés excèdent largement les approximations faites, et modifier la valeur de l'exposant (N) n'y change rien.
Voici les figures obtenues lorsque l'on prend N = 1 , puis N = 1/4:

Par ailleurs, la valeur minimale de la fonction S(M) n'est en rien négligeable devant la plus grande valeur observée en l'un des sommets: S_min^1/4 représente environ 25 % de S_max^1/4, alors que l'on devait trouver pour la fonction étudiée un minimum nul:

Ce résultat est sans appel: il n'y a pas de solution rigoureuse au problème tel qu'il est formulé, et qui se résume à l'égalité S(M) = 0 .

Cela vient de ce que l'énoncé conduit à un système surdéterminé, comportant plus d'équations que de variables ajustables, et ne laissant comme échappatoire que les procédés de recherche des moindres carrés.

Sur trois points, on peut construire un cercle unique, et une infinité de spirales logarithmiques; pour ne retenir qu'une seule de ces dernières, il suffit de fixer l'orientation d'une tangente en l'un des points.

La résolution rigoureuse du système à 3 points et 3 tangentes peut faire illusion lorsque les orientations sont choisies avec soin: le calcul peut conduire à l'une des 3 spirales possibles (il y en a sûrement d'autres !), relativement proches les unes des autres.

[Azzo2]

Merci beaucoup pour vos réponses wiwaxia !