Ellipse autour de points

**wiwaxia** · 01/11/2022, 21h00

Envoyé par tbc92

Imaginons un problème un peu différent. On a 100 points.
On cherche une ellipse E. Et pour évaluer cette ellipse E, on regarde les 70 points P pour lesquels la distance de P à l'ellipse soit minimale, et on fait la somme de ces 70 distances.
Et on cherche à minimiser cette somme ...

L'ellipse est un objet géométrique caractérisé par
a) la position (centre C),
b) la forme (rapport b/a),
c) l'étendue (aire),
d) l'orientation (angle θ) ... soit 5 paramètres, quelle que soit la façon dont on s'y prenne.
Comment définir la distance d'un tel objet à un point ?

On va bientôt arriver à la somme (MF1 + MF2) ...

PS: il y a une difficulté rédhibitoire à la construction d'une ellipse sur 4 points quelconques du nuage: c'est instabilité du résultat, un petit déplacement de l'un de ces points pouvant entraîner une déformation énorme de la courbe en rejetant très loin l'un des foyers. Voire conduire à une hyperbole, si l'un des points est à l'intérieur de l'enveloppe convexe

Tandis qu'un petit déplacement de l'un des foyers ne produit qu'une modification modérée de la même courbe ...

**mach1974** · 02/11/2022, 09h23

Bonjour,
L'algorithme de la General Hough Transform permet de calculer une ellipse en fonction des points

**wiwaxia** · 06/11/2022, 00h50

Il est possible de pousser à son maximum la résolution de la recherche et de déterminer au pixel près la position optimale des foyers, tout en ne testant qu'un nombre raisonnable d'ellipses, à condition de procéder en deux étapes:
a) d'abord en utilisant une grille à 26 carreaux couvrant la totalité de l'image,
b) puis en ayant ensuite recours à 4 grilles comportant 28 carreaux, et respectivement centrées sur les positions des foyers déjà repérées.
L'inventaire se limite alors (si j'ose dire

) à 973366 courbes.
Les valeurs retenues sont en rapport avec les dimensions de l'image - elles doivent vérifier dans le cas présent

N₁*N₂ > La , soit ici N₁*N₂ > 700 ≈ (26.46)² .

S'affichent ci-dessous les résultats obtenus au cours des deux étapes successives, qui reprennent le même algorithme déjà décrit:

Nom : 2 tableaux.png
Affichages : 325
Taille : 35,9 Ko

L'ellipse obtenue passe par un point unique (M₀, croix violette); les points les plus proches de la courbe sont signalés en rouge ou en bleu, selon qu'ils sont situés à l'extérieur ou à l'intérieur; ils sont facilement repérables par la comparaison de la somme des distances qui les séparent des foyers, soit ici l'écart

ε = (MF₁ + MF₂) - (M₀F₁ + M₀F₂)

qui permet d'exprimer la condition de sélection: |ε| < 10 .

Nom : Im_3_Km=26=70%_Ne=707281_Se=72670_Etape 2_P rouges et bleus.png
Affichages : 298
Taille : 7,4 Ko

On remarque que la courbe passe très près d'un certain nombre de points situés à l'intérieur. Cela rejoint les observations de tbc92 et Flodelarab, selon lesquels l'ellipse d'aire minimale doit passer par plusieurs points du nuage. Ce qui ne contredit nullement le fait que, dans la recherche de cette courbe limite, le repérage préalable des foyers reste incontournable.

L'algorithme apparaît lourd à mettre en œuvre et difficilement adaptable, bien qu'un plus grand nombre d'étapes soit envisageable et susceptible de réduire substantiellement la quantité de calculs. Il serait théoriquement possible d'affiner la résolution, en attribuant des coordonnées non entières aux foyers; cela impliquerait toutefois d'apporter d'importantes modifications au programme.
Autre défaut: le programme peut ne pas donner la combinaison optimale si l'un des indices se situe sur le bord de l'une des grilles: la comparaison des surfaces ne va pas suffisamment loin pour que soit livré l'extremum:

Nom : Résultats_27_26.png
Affichages : 294
Taille : 23,9 Ko

Le meilleur algorithme me paraît être la recherche du minimum de la fonction qui, aux coordonnées des foyers fait correspondre l'aire minimale de l'ellipse précédemment définies: Amin= F(X1, Y1, X2, Y2).
Bien qu'il faille explorer un domaine de R⁴, la quantité de calcul ne devrait pas dépasser, à performance égale, celle des méthodes précédentes; et l'on pourrait espérer parvenir à la solution, avec tout la précision disponible: 2^-63 ~ 10^-19 (

... en étant très optimiste).

**Guesset** · 06/11/2022, 14h45

Bonjour wiwaxia,

Envoyé par wiwaxia

...la précision disponible (2^-31 ~ 10^-19

... en étant très optimiste).

Tu es sûr que ce n'est pas plutôt (extended) 2^-63 ~ 10^-19 ?

Salut.

**wiwaxia** · 06/11/2022, 16h08

Bonjour Guesset,

Envoyé par Guesset

... Tu es sûr que ce n'est pas plutôt (extended) 2^-63 ~ 10^-19 ? ...

C'est exact, merci de ta remarque. Il commençait à se faire tard.

PS: il m'est revenu le lointain souvenir du théorème de Pascal concernant les ellipses; malheureusement, cela ne donne rien.

Nom : Théorème de Pascal.png
Affichages : 293
Taille : 6,2 Ko

https://www.bibmath.net/dico/index.p...exapascal.html

**Guesset** · 07/11/2022, 12h05

Bonjour wiwaxia,

Envoyé par wiwaxia

...il m'est revenu le lointain souvenir du théorème de Pascal concernant les ellipses; malheureusement, cela ne donne rien...

Bienvenu au club. J'ai souvent des idées. Elles sont souvent fausses. Et les plus séduisantes sont souvent les plus inutiles

Salut

**wiwaxia** · 11/11/2022, 16h07

Ci-dessous le programme source, que chacun pourra traduire dans son langage usuel.
Il est assez lourd, et sa rédaction rapide a sûrement laissé des instructions inutilement longues.
L'affichage permanent des indices et de quelques autres résultats en cours de calcul en allonge la durée d'exécution; il avait pour finalité de visualiser le déroulement correct du programme.

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
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
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
 PROGRAM Modification_Image;
 
(*HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH
 
 Recherche de la plus petite ellipse (aire minimale)
 
 entourant un certain nombre de points du nuage.
 
 Recherche en deux ‚tapes,
 
 TYPE Z_32 = LongInt; Reel = Extended;
 
 Wt(...), We(...), Wr(...): affichage … l'‚cran d'un texte, d'un entier, d'un
                            r‚el
 
 A_: Pause dans l'ex‚cution du programme, et sortie ‚ventuelle
 
 Copie_F1(Chemin + 'Doss_1\F_1'): saisie des caract‚ristiques d'un fichier
                                  bitmap
 
 Creation_F2(Chemin + 'Doss_2\F_2'): cr‚ation du nouveau fichier bitmap
                                     aux mˆmes dimensions
 
HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH*)
 
 USES Crt, E_Texte, U_Math, U_Copie_1F_ter, Math;
 
 CONST Chemin = 'D:\Virtual_Pascal\Fichiers_VP\Z_Dimension_2\Ellipses\';
       N_Point = 100;   Limite = 70;
       Kmax = 26;
       Marge_1 = 0.050; Marge_2 = 1 - Marge_1;   Plage = Marge_2 - Marge_1;
       m = 255;         P000: Pixel = (0, 0, 0); P111: Pixel = (m, m, m);
                        P100: Pixel = (m, 0, 0); P001: Pixel = (0, 0, m);
                        Px01: Pixel = (160, 0, m);
 
 TYPE Ve_2E = RECORD  x, y: Z_32  END;
      Tab_V = ARRAY[1..N_Point] OF Ve_2E;
      T_Ellipse = RECORD  M0, K1, K2, F1, F2: Ve_2E;
                          Rn: Word; Sd, Se: Z_32  END;
 
 VAR Lst_P, Lst_A: Tab_V;
     Ellipse_1, Ellipse_2: T_Ellipse;
 
(*HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH
 
 Trac‚ des points
 
HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH*)
 
 PROCEDURE Croix(Coul: Pixel; Xc, Yc, La, Ha: Z_32; VAR Ma2: Tab_Pix);
   CONST Lc1 = 5; Lc2 = 2;
   VAR x, X1, X2, y, Y1, Y2: Z_32;
   BEGIN
     X1:= Xc - Lc1; X2:= Xc + Lc1;
     Y1:= Yc - Lc1; Y2:= Yc + Lc1;
     FOR x:= X1 TO X2 DO IF (Abs(x - Xc)>Lc2) THEN Ma2[x, Yc]:= Coul;
     FOR y:= Y1 TO Y2 DO IF (Abs(y - Yc)>Lc2) THEN Ma2[Xc, y]:= Coul;
     X1:= Xc - Lc2; X2:= Xc + Lc2;
     Y1:= Yc - Lc2; Y2:= Yc + Lc2;
     FOR x:= X1 TO X2 DO
       FOR y:= Y1 TO Y2 DO
         IF ((Abs(x - Xc)=Lc2) OR (Abs(y - Yc)=Lc2)) THEN Ma2[x, y]:= Coul
                                                     ELSE Ma2[x, y]:= P000
   END;
 
(*HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH
 
 Trac‚ de l'ellipse
 
 T_Ellipse = RECORD  M0, K1, K2, F1, F2: Ve_2E;
                          Rn: Word; Sd, Se: Z_32  END;
 
HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH*)
 
 FUNCTION Dist2(Wa, Wb: Ve_2E): Z_32;
   VAR X2, Y2: Z_32;
   BEGIN
     X2:= Sqr(Wb.x - Wa.x); Y2:= Sqr(Wb.y - Wa.y);
     Result:= X2 + Y2
   END;
 
 FUNCTION Somm_D(W1, W2, W: Ve_2E): Reel;
   VAR p, q, r, Sd2: Reel;
   BEGIN
     p:= Sqrt(Dist2(W1, W));
     q:= Sqrt(Dist2(W2, W));
     Result:= p + q
   END;
 
 PROCEDURE Tr_Points_Proches(L_P: Tab_V; Ell: T_Ellipse);
   CONST Epsilon = 10;
   VAR j: Byte; k, K0: Word; Xm, Ym: Z_32; Lim_1, Lim_2, Lim_3, u: Reel;
       Vf1, Vf2: Ve_2E;
   BEGIN
     WITH Ell DO BEGIN
                   Vf1:= F1; Vf2:= F2; K0:= Rn
                 END;
     Lim_1:= Somm_D(Vf1, Vf2, Lst_P[K0]);
     Lim_2:= Lim_1 + Epsilon; Lim_3:= Lim_1 - Epsilon;
 
     FOR k:= 1 TO N_Point DO
       IF (k<>K0) THEN BEGIN
                         u:= Somm_D(Vf1, Vf2, Lst_P[k]);
                         IF ((Lim_1<u) AND (u<Lim_2)) THEN
                           BEGIN
                             Xm:= Lst_P[k].x; Ym:= Lst_P[k].y;
                             Croix(P100, Xm, Ym, Larg_Image,
                                   Haut_Image, Matrice_2)
                           END;
                         IF ((Lim_1>u) AND (u>Lim_3)) THEN
                           BEGIN
                             Xm:= Lst_P[k].x; Ym:= Lst_P[k].y;
                             Croix(P001, Xm, Ym, Larg_Image,
                                   Haut_Image, Matrice_2)
                           END;
                       END
   END;
 
 PROCEDURE Calc_ElemE(F_1, F_2, V_M: Ve_2E;
                      VAR C_x, C_y, C_t, S_t, D_A, D_B, A_E: Reel);
   VAR Df2, X1, X2, Y1, Y2: Z_32; Df1, p, q, r, Sd, Sd2: Reel;
   BEGIN
     C_x:= 0.5 * (F_1.x + F_2.x);
     C_y:= 0.5 * (F_1.y + F_2.y);
     X1:= F_2.x - F_1.x; X2:= Sqr(X1);
     Y1:= F_2.y - F_1.y; Y2:= Sqr(Y1);
     Df2:= X2 + Y2;      Df1:= Sqrt(Df2);
     IF (Df2=0) THEN BEGIN
                       C_t:= 1; S_t:= 0
                     END
                ELSE BEGIN
                       C_t:= X1 / Df1; S_t:= Y1 / Df1
                     END;
     X2:= Sqr(F_1.x - V_M.x); Y2:= Sqr(F_1.y - V_M.y);
 
     p:= Sqrt(X2 + Y2);
     X2:= Sqr(F_2.x - V_M.x); Y2:= Sqr(F_2.y - V_M.y);
     q:= Sqrt(X2 + Y2);
     Sd:= p + q;              Sd2:= Sqr(Sd);
     IF (Sd2>Df2) THEN r:= Sqrt(Sd2 - Df2)
                  ELSE r:= 0;
     D_A:= 0.5 * Sd;          D_B:= 0.5 * r;           A_E:= Pi * (D_A * D_B)
   END;
 
 PROCEDURE Tr_Ellipse(Ell: T_Ellipse; VAR Ma2: Tab_Pix);
   CONST D_Pi = 2 * Pi; Np_Ell = 1000; I_Np = 1 / Np_Ell; o = 8; Zmax = 5000;
         Pell: Pixel = (0, 255, 0);
         C1 = 4; L1 = 51; L2 = L1 + 2; v = 16; w = 8;
   VAR k: Word; Xm, Ym: Z_32;
       Aell, CosT, DaxA, DaxB, SinT, Xcen, Ycen,
       Cu, p, q, r, Su, Xell, Yell: Reel;
   BEGIN
     WITH Ell DO
       Calc_ElemE(F1, F2, M0, Xcen, Ycen, CosT, SinT, DaxA, DaxB, Aell);
 
       E(0015); Wt(C1, L1, 'Demi-axes:           (DaxA, DaxB) =');
       E(0013); Write(DaxA:v:w, ' ', DaxB:v:w);
       E(0015); Wt(C1, L2, 'Aire de l''ellipse:             Ae =');
       E(0013); Write(Aell:v:w); q:= 4 / Pi;
       E(0012); Write(' ', q * Aell:v:w);
       FOR k:= 0 TO (Np_Ell - 1) DO
         BEGIN
           q:= k * I_Np;     SinCos(D_Pi * q, Su, Cu);
           Xell:= DaxA * Cu; Yell:= DaxB * Su;
           p:= Xell * CosT;  q:= Yell * SinT; r:= p - q; Xm:= Round(r + Xcen);
           p:= Xell * SinT;  q:= Yell * CosT; r:= p + q; Ym:= Round(r + Ycen);
           IF ((Xm>=0) AND (Xm<Zmax)) AND ((Ym>=0) AND (Ym<Zmax))
             THEN Ma2[Xm, Ym]:= Pell
         END;
 
     WITH Lst_P[Ell.Rn] DO BEGIN
                             Xm:= x; Ym:= y
                           END;
     E(0015); Wt(C1, L1 - 2, 'Point Limite:      ');
     E(0010); Write(Ell.Rn:o);
     E(0009); Write(Xm:o, Ym:o);
     Croix(Px01, Xm, Ym, Larg_Image, Haut_Image, Matrice_2); A_
   END;
 
(*HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH
 
 Recherche de l'ellipse minimale
 
HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH*)
 
 PROCEDURE Aff_E(VAR Ell: T_Ellipse);
   CONST C1 = 4; L1 = 38; L2 = L1 + 2; L3 = L2 + 2; L4 = L3 + 2; o = 8;
   BEGIN
     WITH Ell DO
       BEGIN
         E(0015); Wt(C1, L1, 'Indices des foyers:      Ki(x, y) = ');
         E(0011); Write(K1.x:o, K1.y:o, K2.x:o, K2.y:o);
         E(0015); Wt(C1, L2, 'Coordonn‚es des foyers:  Fi(x, y) = ');
         E(0010); Write(F1.x:o, F1.y:o, F2.x:o, F2.y:o);
         E(0015); Wt(C1, L3, 'Aire de l''ellipse (ž4/Pi)      Se = ');
         E(0012); Write(Se:8)
       END
   END;
 
 PROCEDURE Aff_L(L_D: Tab_V);
   CONST Lcol = 20; Hcol = 25; o = 5;
   VAR c, x, y: Byte; k: Word;
   BEGIN
     FOR k:= 1 TO N_Point DO
       BEGIN
         c:= (k - 1) DIV Hcol; x:= c * Lcol; Inc(x);
         y:= (k - 1) MOD Hcol; Inc(y, 51);
         E(0012);              We(x, y, k, o);
         IF (k=Limite) THEN E(0014) ELSE E(0011);
         Write(L_D[k].x:o, L_D[k].y:9)
       END;
     A_
   END;
 
(*HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH
 
 Balayage du domaine et tri des listes selon l'ordre croissant
 
          Si F1F2 = 0 alors D_a = D_b = Sqrt(X2 + Y2) ; Aire = Pi * D_a * D_b
 
          D_a:= 0.5 * Sd; D_b:= 0.5 * r; Se = 4 * Da * Db
          Aire = Pi * D_a * D_b = (Pi/4) * Sd * r
 
HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH*)
 
 PROCEDURE Tri_C(VAR L_A: Tab_V);
   VAR i, j, k: Word; Min: Z_32; Vmin, W: Ve_2E; Liste: Tab_V;
   BEGIN
     Liste:= L_A;
     FOR i:= 1 TO (N_Point - 1) DO
       BEGIN
         k:= i; Min:= Liste[i].y;
         FOR j:= (i + 1) TO N_Point DO
           IF (Min>Liste[j].y) THEN BEGIN
                                      Min:= Liste[j].y; k:= j;
                                      Vmin:= Liste[k]
                                    END;
         IF (k>i) THEN BEGIN
                            W:= Liste[i]; Liste[i]:= Vmin;
                            Liste[k]:= W;
                          END
       END;
     L_A:= Liste
   END;
 
 PROCEDURE Aff_E1(E_1: T_Ellipse);
   CONST o = 8;
   BEGIN
     WITH E_1 DO BEGIN
                   E(0114); We(23, 47, K1.x, o);
                   Write(K1.y:o, K2.x:o, K2.y:o);
                   E(0012); Write(Se:9); E(0010); Write(Rn:o); E(0011)
                 END
   END;
 
 FUNCTION Calc_P(F_1, F_2, V_M: Ve_2E): Z_32;
   VAR Df2, X2, Y2: Z_32; p, q, r, Sd, Sd2: Reel;
   BEGIN
     X2:= Sqr(F_1.x - V_M.x); Y2:= Sqr(F_1.y - V_M.y); p:= Sqrt(X2 + Y2);
     X2:= Sqr(F_2.x - V_M.x); Y2:= Sqr(F_2.y - V_M.y); q:= Sqrt(X2 + Y2);
     Sd:= p + q;              Sd2:= Sqr(Sd);
     X2:= Sqr(F_2.x - F_1.x); Y2:= Sqr(F_2.y - F_1.y); Df2:= X2 + Y2;
     IF (Sd2>Df2) THEN r:= Sqrt(Sd2 - Df2)
                  ELSE r:= 0;
     Result:= Round(Sd * r)
   END;
 
 PROCEDURE Calc_MinMax(Ell_1: T_Ellipse;
                       VAR D_, X1_a, X1_b, X2_a, X2_b,
                       Y1_a, Y1_b, Y2_a, Y2_b: Z_32);
   BEGIN
     WITH Ell_1 DO
       BEGIN
         X1_a:= F1.x - D_; X1_b:= F1.x + D_;
         Y1_a:= F1.y - D_; Y1_b:= F1.y + D_;
         X2_a:= F2.x - D_; X2_b:= F2.x + D_;
         Y2_a:= F2.y - D_; Y2_b:= F2.y + D_
       END
   END;
 
 PROCEDURE Balayage_2(La, Ha: Z_32; Ell_1: T_Ellipse; VAR Ell_2: T_Ellipse);
   CONST I_Km = 1 / Kmax;
         d = 8;   C1 = 23; C2 = C1 + d; C3 = C2 + d; C4 = C3 + d;
         L1 = 45; L2 = L1 + 2; o = 8;
   VAR j: Word;
          Delta, Ha1, Kx1, Kx2, Ky1, Ky2, La1, N_Ell, Smin,
          X1min, X1max, X2min, X2max, Y1min, Y1max, Y2min, Y2max: Z_32;
       w: Reel; Vf1, Vf2: Ve_2E; Ell_: T_Ellipse;
   BEGIN
     La1:= La - 1; Ha1:= Ha - 1; Smin:= Sqr(La + Ha); N_Ell:=0;
     IF (La>Ha) THEN Delta:= Ceil(La / (2 * Kmax))
                ELSE Delta:= Ceil(Ha / (2 * Kmax));
 
     E(0015);      Wt(4, L1, 'Indices courants:');
     Wt(4, L2, 'Valeurs optimales:'); E(0011);
 
     FOR Kx1:= -Delta TO Delta DO
       BEGIN
         Vf1.x:= Ell_1.F1.x + Kx1; We(C1, L1, Kx1, o);
         FOR Ky1:= -Delta TO Delta DO
           BEGIN
             Vf1.y:= Ell_1.F1.y + Ky1; We(C2, L1, Ky1, o);
             FOR Kx2:= -Delta TO Delta DO
               BEGIN
                 Vf2.x:= Ell_1.F2.x + Kx2; We(C3, L1, Kx2, o);
                 FOR Ky2:= -Delta TO Delta DO
                   BEGIN
                     Vf2.y:= Ell_1.F2.y + Ky2; We(C4, L1, Ky2, o);
                     IF ((Vf2.x>Vf1.x) OR ((Vf2.x=Vf1.x) AND (Vf2.y>=Vf1.y)))
                       THEN BEGIN
                              Inc(N_Ell);     Write('    Nell = ', N_Ell:6);
                               FOR j:= 1 TO N_Point DO
                                 BEGIN
                                   Lst_A[j].x:= j;
                                   Lst_A[j].y:= Calc_P(Vf1, Vf2, Lst_P[j])
                                 END;
                              Tri_C(Lst_A); //Aff_L(Lst_A);
                              IF (Smin>Lst_A[Limite].y)
                                THEN BEGIN
                                       Smin:= Lst_A[Limite].y;
                                       WITH Ell_ DO
                                         BEGIN
                                           K1.x:= Kx1; K1.y:= Ky1;
                                           K2.x:= Kx2; K2.y:= Ky2;
                                           F1:= Vf1;   F2:= Vf2;
                                           Se:= Smin;  Rn:= Lst_A[Limite].x;
                                           M0:= Lst_P[Rn];
                                         END;
                                       Aff_E1(Ell_)
                                     END
                            END
                   END
               END
           END
       END;
     Ell_2:= Ell_; Aff_E(Ellipse_2)
   END;
 
 PROCEDURE Balayage_1(La, Ha: Z_32; VAR Ell_1: T_Ellipse);
   CONST I_Km = 1 / Kmax;
         d = 8;   C1 = 23; C2 = C1 + d; C3 = C2 + d; C4 = C3 + d;
         L1 = 45; L2 = L1 + 2; o = 8;
   VAR j: Word; Ha1, Iini, Kx1, Kx2, Ky1, Ky2, La1, N_Ell, Smin: Z_32;
       w: Reel; Vf1, Vf2: Ve_2E; Ell_: T_Ellipse;
   BEGIN
     La1:= La - 1; Ha1:= Ha - 1; Smin:= Sqr(La + Ha);
     N_Ell:=0;     E(0015);      Wt(4, L1, 'Indices courants:');
     Wt(4, L2, 'Valeurs optimales:'); E(0011);
     FOR Kx1:= 0 TO Kmax DO
       BEGIN
         w:= I_Km * Kx1; Vf1.x:= Round(La1 * w); We(C1, L1, Kx1, o);
         FOR Ky1:= 0 TO Kmax DO
           BEGIN
             w:= I_Km * Ky1; Vf1.y:= Round(Ha1 * w); We(C2, L1, Ky1, o);
             FOR Kx2:= Kx1 TO Kmax DO
               BEGIN
                 w:= I_Km * Kx2; Vf2.x:= Round(La1 * w); We(C3, L1, Kx2, o);
                 IF (Kx2=Kx1) THEN Iini:= Ky1
                              ELSE Iini:= 0;
                 FOR Ky2:= Iini TO Kmax DO
                   BEGIN
                     w:= I_Km * Ky2; Vf2.y:= Round(Ha1 * w);
                     We(C4, L1, Ky2, o);
                     Inc(N_Ell);     Write('    Nell = ', N_Ell:6);
                     FOR j:= 1 TO N_Point DO
                       BEGIN
                         Lst_A[j].x:= j;
                         Lst_A[j].y:= Calc_P(Vf1, Vf2, Lst_P[j])
                       END;
                     Tri_C(Lst_A); //Aff_L(Lst_A);
                     IF (Smin>Lst_A[Limite].y)
                       THEN BEGIN
                              Smin:= Lst_A[Limite].y;
                              WITH Ell_ DO
                                BEGIN
                                  K1.x:= Kx1; K1.y:= Ky1; K2.x:= Kx2;
                                  K2.y:= Ky2; F1:= Vf1;   F2:= Vf2;
                                  Se:= Smin;  Rn:= Lst_A[Limite].x;
                                  M0:= Lst_P[Rn];
                                END;
                              Aff_E1(Ell_)
                            END
                   END
               END
           END
       END;
     Ell_1:= Ell_; Aff_E(Ellipse_1); A_; FenI(3700)
   END;
 
(*HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH
 
 Ensemencement du domaine et repr‚sentation du nuage
 
HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH*)
 
 PROCEDURE Ensemencement(La, Ha: Z_32; VAR L_P: Tab_V; VAR Ma2: Tab_Pix);
   CONST Mu = 0.750;
   VAR j: Word; Ha1, Ha2, La1, La2, Xm, Ym: Z_32;
       Lim1, Lim2, p, Ux, Uy: Reel; W: Ve_2E;
 
   BEGIN
     La1:= Round(Plage * La); La2:= Round(Marge_1 * La);
     Ha1:= Round(Plage * Ha); Ha2:= Round(Marge_1 * Ha);
     FOR Xm:= La2 TO (La2 + La1 - 1) DO
       FOR Ym:= Ha2 TO (Ha2 + Ha1 - 1) DO Ma2[Xm, Ym]:= P000;
 
     RandSeed:= 1987654321;
     FOR j:= 1 TO N_Point DO
       BEGIN
         REPEAT
           Ux:= Random;          Uy:= Random; p:= Ux * (1 - Ux);
           Lim1:= Ux - (Mu * p); Lim2:= Ux + (Mu * p);
         UNTIL ((Lim1<Uy) AND (Uy<Lim2));
         W.x:= Round(La1 * Ux); Inc(W.x, La2);
         W.y:= Round(Ha1 * Uy); Inc(W.y, Ha2);
         L_P[j]:= W;
         Croix(P111, W.x, W.y, Larg_Image, Haut_Image, Matrice_2);
       END
   END;
 
 PROCEDURE Ensemencement_01(La, Ha: Z_32; VAR L_P: Tab_V; VAR Ma2: Tab_Pix);
   CONST Mu = 0.25; Nu = 1 - Mu;
   VAR j: Word; Ha1, Ha2, La1, La2, Xm, Ym: Z_32;
       Lim, Ux, Uy: Reel; W: Ve_2E; Test1, Test2: Bool;
   BEGIN
     La1:= Round(Plage * La); La2:= Round(Marge_1 * La);
     Ha1:= Round(Plage * Ha); Ha2:= Round(Marge_1 * Ha);
     FOR Xm:= La2 TO (La2 + La1 - 1) DO
       FOR Ym:= Ha2 TO (Ha2 + Ha1 - 1) DO Ma2[Xm, Ym]:= P000;
 
     RandSeed:= 1333444999;
     FOR j:= 1 TO N_Point DO
       BEGIN
         REPEAT
           Ux:= Random; Uy:= Random; Lim:= Nu * Sqr(Ux);
           Test1:= (Nu<Uy);
           Test2:= (Lim<Uy) AND (Uy<Lim + Mu)
         UNTIL (Test1 OR Test2);
         W.x:= Round(La1 * Ux); Inc(W.x, La2);
         W.y:= Round(Ha1 * Uy); Inc(W.y, Ha2);
         L_P[j]:= W;
         Croix(P111, W.x, W.y, Larg_Image, Haut_Image, Matrice_2);
       END
   END;
 
 PROCEDURE Ensemencement_00(La, Ha: Z_32; VAR L_P: Tab_V; VAR Ma2: Tab_Pix);
   VAR j: Word; Ha1, Ha2, La1, La2, Xm, Ym: Z_32; W: Ve_2E;
   BEGIN
     La1:= Round(Plage * La); La2:= Round(Marge_1 * La);
     Ha1:= Round(Plage * Ha); Ha2:= Round(Marge_1 * Ha);
     FOR Xm:= La2 TO (La2 + La1 - 1) DO
       FOR Ym:= Ha2 TO (Ha2 + Ha1 - 1) DO Ma2[Xm, Ym]:= P000;
 
     RandSeed:= 1333444111;
     FOR j:= 1 TO N_Point DO
       BEGIN
         W.x:= Random(La1); Inc(W.x, La2);
         W.y:= Random(Ha1); Inc(W.y, Ha2);
         L_P[j]:= W;
         Croix(P111, W.x, W.y, Larg_Image, Haut_Image, Matrice_2);
       END
   END;
 
 PROCEDURE Init_M(La, Ha: Z_32; VAR Ma2: Tab_Pix);
   CONST Pcoul: Pixel = (0, 0, 150);
   VAR Xm, Ym: Z_32;
   BEGIN
     FOR Xm:= 0 TO (La - 1) DO
       FOR Ym:= 0 TO (Ha - 1) DO Ma2[Xm, Ym]:= Pcoul
   END;
 
(*HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH
 
 Programme principal
 
HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH*)
 
 BEGIN
   Copie_F1(Chemin + 'Doss_1\F_1');
 
   Init_M(Larg_Image, Haut_Image, Matrice_2);
   Ensemencement(Larg_Image, Haut_Image, Lst_P, Matrice_2);
   Balayage_1(Larg_Image, Haut_Image, Ellipse_1);
   Balayage_2(Larg_Image, Haut_Image, Ellipse_1, Ellipse_2);
 
   Tr_Ellipse(Ellipse_2, Matrice_2);
   Tr_Points_Proches(Lst_P, Ellipse_2);
 
   Creation_F2(Chemin + 'Doss_2\F_2'); A_
 END.

Nom : Image_5.png
Affichages : 279
Taille : 58,4 Ko

Deux détails de calcul un peu délicats:
1°) L'inclinaison (θ) du grand axe de l'ellipse est implicitement définie par les rapports:

Ct = Cos(θ) = (F2.x - F1.x)/Df1 , St = Sin(θ) = (F2.y - F1.y)/Df1 ,
avec Df1 = (Df2)^1/2 = (F₁F₂) ,

ce qui suppose des foyers non confondus; dans le cas contraire (F₁F₂ = 0), la courbe devient circulaire, perd son excentricité et l'angle (θ), alors indéterminé, peut faire l'objet d'un choix arbitraire:

θ = 0 , d'où: Ct = 1 et St = 0 .

2°) Autrement plus sournoise est l'expression du réel

r = ((MF₁ + MF₂)² - (F₁F₂)²)^1/2 = ((p + q)² - Df2)^1/2

incontournable dans le calcul des caractéristiques de l'ellipse: b = r / 2 .
Bien que mathématiquement irréprochable, cette formule où interviennent deux variables (p, q) au format Float constitue, au voisinage de zéro, une source redoutable de plantage aléatoire du programme - d'autant quelle peut être appelée plusieurs centaines de milliers de fois !

D'où la nécessité d'un calcul conditionnel, aprés contrôle du signe de la différence

Les deux écueils sont ainsi contournés lors du calcul des éléments 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
23
24
25
 PROCEDURE Calc_ElemE(F_1, F_2, V_M: Ve_2E;
                      VAR C_x, C_y, C_t, S_t, D_A, D_B, A_E: Reel);
   VAR Df2, X1, X2, Y1, Y2: Z_32; Df1, p, q, r, Sd, Sd2: Reel;
   BEGIN
     C_x:= 0.5 * (F_1.x + F_2.x);
     C_y:= 0.5 * (F_1.y + F_2.y);
     X1:= F_2.x - F_1.x; X2:= Sqr(X1);
     Y1:= F_2.y - F_1.y; Y2:= Sqr(Y1);
     Df2:= X2 + Y2;      Df1:= Sqrt(Df2);
     IF (Df2=0) THEN BEGIN    // Calcul de Ct et St
                       C_t:= 1; S_t:= 0
                     END
                ELSE BEGIN
                       C_t:= X1 / Df1; S_t:= Y1 / Df1
                     END;
     X2:= Sqr(F_1.x - V_M.x); Y2:= Sqr(F_1.y - V_M.y);

     p:= Sqrt(X2 + Y2);
     X2:= Sqr(F_2.x - V_M.x); Y2:= Sqr(F_2.y - V_M.y);
     q:= Sqrt(X2 + Y2);
     Sd:= p + q;              Sd2:= Sqr(Sd);
     IF (Sd2>Df2) THEN r:= Sqrt(Sd2 - Df2)    // Calcul de (r)
                  ELSE r:= 0;
     D_A:= 0.5 * Sd;          D_B:= 0.5 * r;           A_E:= Pi * (D_A * D_B)
   END;

Résultat obtenu à partir d'un nuage coudé comportant 100 points, dont 70 à l'intérieur de la courbe:

Nom : Km=26_Np=100_70%_2 stades_Ens01_N coudé.png
Affichages : 261
Taille : 7,6 Ko

Ce problème en évoque un autre, celui de la recherche de l'ellipse passant par quatre points donnés - en supposant évidemment qu'elle existe.
Ils ont pour point commun la localisation progressive des foyers.

**Sylvain255** · 18/11/2022, 22h14

j'ai trouvé cet article : http://ada-posturologie.fr/EllipseConfiance.htm

mon article fait-il avancer la discussion ?

**tbc92** · 19/11/2022, 18h57

Non, même si tes données vérifient des conditions dont tu n'as pas parlé.

Si je ne m'abuse, cet article traite des données où il y a corrélation linéaire entre les 2 axes (par exemple des individus, pour lesquels on aurait mesuré la taille et le poids, ou des appartements, avec le prix et la surface).
Les points sont donc plus ou moins sur un nuage de forme d'ellipse.
Si tes points sont totalement quelconques (par exemple 2 paquets assez espacés, ou en forme en V, ou en T ...), ça ne marche pas du tout.

De plus, c'est une approche probabiliste. C'est à dire qu'on part de 2 ou 3 indicateurs synthétiques (la moyenne, la variance, la covariance) et on en déduit l'ellipse.
L'ellipse obtenue peut contenir moins de 90% des points, et elle peut aussi être un peu plus grande que l'ellipse minimale.
Il faudrait vraiment un coup de chance énorme (=impossible) pour que l'ellipse obtenue contienne 90% des données, et soit l'ellipse minimale.

**Sylvain255** · 19/11/2022, 18h59

donc la bonne réponse est celle de :
11/11/2022, 16h07
wiwaxia

?

mais je ne comprends pas du tout son algo, ça me dépasse !
moi ce que je veux c'est un algo qui prend en entrée une suite de points (x,y) et le pourcentage de points attendu à l'intérieur de l'ellipse et qui nous donne en sortie les paramètres d'une ellipse
je pense qu'on peut dire que le minimum de points à fournir est trois points

je viens de trouver cet article aussi :
https://www.ngs.noaa.gov/PUBS_LIB/Al...OS107_CGS3.pdf

**Sylvain255** · 19/11/2022, 19h20

Envoyé par Sylvain255

je viens de trouver cet article aussi :
https://www.ngs.noaa.gov/PUBS_LIB/Al...OS107_CGS3.pdf

Non plus, cet article parle de 'confidence' (=confiance). Autrement dit, approche probabiliste : on a toutes raisons de penser que le bateau est à tel endroit, basé sur les mesures GPS, et les imprécisions liées à ces mesures. Et on trace un cercle ou une ellipse.

Envoyé par Sylvain255

donc la bonne réponse est celle de :
11/11/2022, 16h07
wiwaxia

?

mais je ne comprends pas du tout son algo, ça me dépasse !
moi ce que je veux c'est un algo qui prend en entrée une suite de points (x,y) et le pourcentage de points attendu à l'intérieur de l'ellipse et qui nous donne en sortie les paramètres d'une ellipse
je pense qu'on peut dire que le minimum de points à fournir est trois points

Je pense que Wiwaxia est le mieux placé pour donner les limites/les impasses de son programme.
Quoi qu'il en soit, c'est la seule proposition aboutie. Et je pense que tu n'auras pas mieux.

**Sylvain255** · 19/11/2022, 21h37

ok merci

Ellipse autour de points

Algorithmes et structures de données

Discussions similaires

Partager

Partager