Discussion :

[tbc92]

@Guesset : on veut l'ellipse la plus petite possible (= de surface minimale). Donc une ellipse qui contiendrait tous les points ne sera a priori pas optimale, si on vise 60% ou 80%.

Trouver une droite qui passe par un seul point : oui, c'est toujours possible. Passant par un point, la liste des directions interdites est un nombre fini, et on a une infinité de directions possibles, et donc, c'est une certitude, on peut toujours construire une droite passant par un point et un seul.

[wiwaxia]

Soit un nuage de (N) points, de coordonnées barycentriques

x_G = Σ_j=0^N-1x_j , y_G = Σ_j=0^N-1y_j ,

à partir desquelles on peut calculer les sommes:

S_x2 = Σ_j=0^N-1(x_j - x_G)² , S_y2 = Σ_j=0^N-1(y_j - y_G)² , S_xy = Σ_j=0^N-1(x_j - x_G)(y_j - y_G) .

L'étalement du nuage de points selon une direction donnée de vecteur unitaire u = (Cos(θ), Sin(θ)) est mesuré par la somme des carrés des produits scalaires:

S(θ) = Σ_j=0^N-1(u.GM_j)² = Σ_j=0^N-1(Cos(θ)*(x_j - x_G) + Sin(θ)*(y_j - y_G))²
= Cos(θ)²S_x2 + Sin(θ)²S_y2 + 2Sin(θ)Cos(θ)S_xy .

Cette somme, qui dépend de (θ), admet pour dérivée:

S'(θ) = 2Sin(θ)Cos(θ)(S_y2 - S_x2) + 2(Cos(θ)² - Sin(θ)²)S_xy = Sin(2θ)(S_y2 - S_x2) + 2Cos(2θ)S_xy ,

et ses extremums sont caractérisés par l'annulation de cette dernière fonction, l'égalité S'(θ) = 0 impliquant:

Tan(2θ) = 2S_xy/(S_x2 - S_y2) = T (par convention),

soit encore

2θ = Arctan(T) + kπ (avec k = 0 ou 1) ,

ce qui conduit à deux directions orthogonales:

θ₁ = (1/2)Arctan(T) et θ₂ = (1/2)Arctan(T) + π/2 .

L'égalité S_y2 = S_x2 , peu probable quoique non exclue, correspondrait à un nuage de forme circulaire, dépourvu de tout allongement.

Dans le repère (XGY) utilisant les deux axes ainsi définis, on peut associer au nuage de points l'ellipse d'équation

X²/S_x2 + Y²/S_y2 = 1 :

le retour aux anciennes coordonnées s'effectue par l'intermédiaire des équations:

x = x_G + X.Cos(θ) - Y.Sin(θ) , y = y_G + X.Sin(θ) + Y.Cos(θ) (avec θ = θ₁) ,

d'où l'expression des nouvelles:

X = (x - x_G)Cos(θ) + (y - y_G).Sin(θ) , Y = -(x - x_G).Sin(θ) + (y - y_G).Cos(θ) .

Ainsi toute famille d'ellipses coaxiales à la précédente et de même forme admettra pour équation générale

X²/S_x2 + Y²/S_y2 = λ² ;

leurs demi-diamètres restant toujours dans le même rapport:

B'/A' = λB/λA = B/A = (S_y2/S_x2)^1/2 .

Une chose est certaine, il y a 2 phases indépendantes :
Fixer la frontière, c'est-à-dire le nombre de points (qu'il soit "strictement égal" ou "au moins égal"). Cette phase est purement arithmétique. Et "facile".
Déterminer l'ellipse et les points intérieurs. Moins facile.

Il suffira de consigner toutes les valeurs de (λ) ou (λ²) dans un tableau, et de trier celui-ci selon l'ordre croissant: la réponse sera immédiate.

Calculs à vérifier.

[wiwaxia]

Bonjour Guesset,

... Je pense que tu voulais mettre 2.k.n.Pi/N sinon cela ne parcourt qu'environ un tiers d'ellipse ...

C'était bien la bonne expression 2.k.Pi/N

mais la graphie du caractère grec n'est vraiment pas réussie, et prête à confusion.

La solution graphique que tu proposes, bien que conforme à la lettre de l'énoncé, me paraît aberrante: une parenté de forme entre le nuage de points et la frontière cherchée me semble aller de soi, on peut sinon obtenir n'importe quoi. Mais c'est peut-être là un point de vue trop restreint.

[Guesset]

Bonjour,

Je n'ai pas dit que l'ellipse passe par un seul point, j'ai juste montré qu'il était possible de partager le plan avec une droite qui, au plus a un point de l'ensemble. Le nombre de points entourés est le nombre de points dans l'un des demi-plan (+1 s'il y a un point sur la droite de clivage).

L'objectif n'était pas de trouver la plus petite ellipse mais de répondre à la question qui restait en suspend : y a-t il une ellipse qui entoure un nombre donné de points ? Cela semblait un préalable avant de chercher la plus petite.

J'étais toujours sur un nombre n donné. Si c'est n et plus, la contrainte est moins forte et il y a évidemment une solution voire plusieurs.

L'isobarycentre est tentant mais se fait assez rapidement piégé par une distribution en haltère, de même la densité n'est pas sans risque. Illustration rapide et très perfectible :

Salutations

[Guesset]

Bonjour wiwaxia,

C'était bien la bonne expression 2.k.Pi/N

Même en gros, j'ai du mal à reconnaître PI mais j'aurais du y penser car, à défaut, il y aurait eu deux variables k et n. Mea culpa .

...La solution graphique que tu proposes, bien que conforme à la lettre de l'énoncé, me paraît aberrante: une parenté de forme entre le nuage de points et la frontière cherchée me semble aller de soi, on peut sinon obtenir n'importe quoi...

Si je comprends bien, tu partirais de l'ellipse minimale des 100% que tu réduirais progressivement sans la déplacer ni la tourner (sinon tu pourrais retomber sur le genre de solution du dessin).

Cela montre que le problème mériterait d'être précisé. Si chacun cherche une solution à des problèmes voisins mais subtilement différents, on peut discuter et s'amuser longtemps en vain

Salut

[wiwaxia]

... Si je comprends bien, tu partirais de l'ellipse minimale des 100% que tu réduirais progressivement sans la déplacer ni la tourner (sinon tu pourrais retomber sur le genre de solution du dessin) ...

C'est effectivement cela. Et dans le cas d'un nuage "normal" en forme de cigare, dépourvu de lacunes, le procédé a l'avantage de conduire rapidement à une solution.
Et rien n'empêche de se limiter d'emblée aux (N') points les plus proches du barycentre, pour un meilleur ajustement.

J'ai peut-être tort d'envisager le problème sous le biais de la statistique; cependant envisager les configurations les plus inattendues conduit inéluctablement à un balayage multiple, une ellipse quelconque étant caractérisée par les coordonnées de son centre, l'inclinaison de son grand axe et son excentricité - soit en tout 4 variables.

Intervient de plus la question du critère de comparaison: comment distinguer la plus petite de deux ellipses ?
Par son périmètre (à rejeter d'emblée !), son aire ? Ou par la somme des distances aux foyers (MF1 + MF2) ?
L'énoncé aurait pu effectivement donner plus de précisions ...

[Guesset]

Bonjour wiwaxia,

...Et dans le cas d'un nuage "normal" en forme de cigare, dépourvu de lacunes, le procédé a l'avantage de conduire rapidement à une solution...

Si je conçois bien des lacunes dans une une surface continue, j'ai un peu plus de mal avec une collection de points qui n'est pas une surface.

Tout triangle entre points proches (i.e. sans point à l'intérieur du triangle, triangulations de Delaunay par exemple) pourrait, me semble-t-il, être assimilé à une lacune à moins de fixer un seuil de densité peut être ?

Plus petite ellipse : j'aurais tendance à la définir comme étant celle qui a la plus petite surface et, en cas d'égalité, celle qui a la distance inter-foyer la plus faible (ce qui, sauf erreur, est équivalent au plus petit périmètre).

Salut

[wiwaxia]

Bonjour Guesset,

Si je conçois bien des lacunes dans une une surface continue, j'ai un peu plus de mal avec une collection de points qui n'est pas une surface.
Tout triangle entre points proches (i.e. sans point à l'intérieur du triangle, triangulations de Delaunay par exemple) pourrait, me semble-t-il, être assimilé à une lacune à moins de fixer un seuil de densité peut être ? ...

Toute structure granulaire comporte nécessairement des interstices ! Ce qui est en cause, c'est l'existence de domaines placés à l'intérieur de l'enveloppe convexe du nuage, vides ou sinon de densité très inférieure à celles observée sur d'autres zones, et dont les dimensions dépassent largement la distance moyenne séparant les points les plus proches.
C'est par exemple le cas de la distribution pseudo-aléatoire ci-dessous (comme celui de quelques exemples cités auparavant), et dont voit mal le rapport de forme avec une ellipse:

Le procédé de Delaunay fera apparaître quelques triangles anormalement grands, d'aire et de dimensions très supérieures à celles des triangles voisins, largement majoritaires et présents sur le reste de la structure.
Ta remarque me donne d'ailleurs de nouvelles idées concernant des images mathématiques

... Plus petite ellipse : j'aurais tendance à la définir comme étant celle qui a la plus petite surface et, en cas d'égalité, celle qui a la distance inter-foyer la plus faible (ce qui, sauf erreur, est équivalent au plus petit périmètre) ...

Cela paraît un bon choix.

[wiwaxia]

Bonjour,

Voici le tracé de l'ellipse d'aire minimale entourant 70 des 100 points d'un nuage donné. La courbe passe par exactemnt par le point tracé en rouge.
Le résultat est légèrement affecté par le nombre de courbes testées (7381 dans le premier cas, 1413721 dans le suivant):

Le résultat devient moins prévisible dans le cas d'un nuage de forme plus complexe:

[tbc92]

Tu as mis un seul point rouge par dessin, mais j'imagine qu'il y en a au moins 2.

[wiwaxia]

Non (sauf cas peu probable), parce que les positions des foyers sont prédéterminées.

Je donnerai ultérieurement plus de détails.

[tbc92]

J'ai dû rater un épisode. Disons que dans ton processus, tu imposes les foyers, et ensuite tu réduis l'ellipse autant que possible. Et dans ce cas, évidemment, tu as un seul point de contact, sauf cas exceptionnel.

[wiwaxia]

J'ai repris le problème sous une autre perspective, que j'avais d'abord délaissée parce qu'elle paraissait conduire à une impasse.
Le procédé donne néanmoins des résultats, malgré la lourdeur des calculs et des difficultés inattendues.

[Flodelarab]

Ou rien du tout. Je n'ai pas connaissance de véritables résultats de convergence d'un algo génétique, même sur un problème convexe (ce qui est le cas ici), c'est-à-dire extrêmement facile à résoudre (d'un point de vue théorie de la complexité, pas forcément en pratique ) ; le mieux que j'ai, c'est des choses comme https://ls11-www.cs.tu-dortmund.de/p...papers/CCJ.pdf, avec des résultats d'une utilité douteuse. L'approche que j'ai citée, à base de SDP, aura la solution optimale en temps polynomial (en pratique, je suppose qu'il faudra un très bon solveur comme présenté sur http://plato.asu.edu/ftp/sparse_sdp.html).

• Sur la forme, je n'aime pas beaucoup les pouces baissés. On n'est pas là pour travailler à la place des internautes, mais leur ouvrir des horizons. L'idée de l'algo génétique n'a pas à être moinssée.
• Sur le fond, tu dénies la capacité des algorithmes génétiques à arriver à l'optimal. Pourtant, ils sont créés pour cela.

Cela m'a suffisamment énerver pour que je le fasse.


J'ai donc considéré 4 gènes :

• abscisse du centre de l'ellipse
• ordonnée du centre de l'ellipse
• angle de l'axe de l'ellipse avec l'horizontal
• grand rayon de l'ellipse

Le petit rayon n'est pas un gène. Il est fixé automatiquement par la quantité optimale de points qu'il doit englober. Si au moins n points sont dans le cercle, de centre le centre de l'ellipse, et de rayon le grand rayon de l'ellipse, alors il existe une ellipse solution optimale contenant exactement n points. Je considère les points suffisamment aléatoires pour que 2 points ne trouvent pas exactement le même petit rayon d'ellipse.
J'obtiens ça (avec 10000 individus et 100 brassages (par pourcentage recherché):

J'ai repris le problème sous une autre perspective, que j'avais d'abord délaissée parce qu'elle paraissait conduire à une impasse.
Le procédé donne néanmoins des résultats, malgré la lourdeur des calculs et des difficultés inattendues.

@wiwaxia : ta méthode, comme la mienne, s'avère bancale. Le fait qu'un seul point serve de contact à l'extérieur montre que le résultat est faux. Il suffirait de déplacer légèrement le centre pour réduire l'ellipse. Les points de contact extérieurs sont forcément multiples. Et là, je soupire.

[tbc92]

Ici, on a un domaine de validité : 'sont valides toutes les ellipses qui contiennent au moins n points', et une fonction objectif : 'minimiser la surface'.
Et forcément, la solution optimale est à la frontière du domaine de validité.
Et quelque chose me dit que les outils divers et variés sont efficaces pour rechercher un optimum quand cet optimum est au 'centre' du domaine de validité, mais que ça marche moins bien dès qu'on doit chercher à la frontière.
Ici, non seulement la solution optimale passe par au moins 2 points de notre jeu de données, mais je pense qu'elle passe même par 3 points. On est donc plus qu'à la frontière.

Il y a peut-être un moyen de chercher parmi toutes les ellipses qui passent par 3 pts... possible. Pour chaque triplet de point, il y a une infinité d'ellipse qui passe par ces 3 points, mais ça paraît """facile""" de trouver la plus petite, qui contiendrait n pts du nuage.
Et dans un second temps, trouver la plus petite des ellipses parmi les n(n-1)(n-2)/6 ellipses qu'on vient de trouver, c'est facile aussi.

Reste à valider que la plus petite ellipse passe forcément par 3 pts.

[Flodelarab]

Il faut 5 points pour définir une conique. Cela peut nécessiter 5 points frontière.

[wiwaxia]

La recherche de l'ellipse d'aire minimale enlaçant un nombre donné d'un nuage de points peut s'effectuer comme suit:

On suppose les coordonnées des points du nuage stockées dans un tableau de vecteurs à composantes entières au format LongInt:

Code Pascal :

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

1
2
3
 TYPE Ve_2E = RECORD  x, y: Z_32  END;
      Tab_V = ARRAY[1..N_Point] OF Ve_2E;
 VAR Lst_P: Tab_V;

1°) Procéder à un balayage du domaine de l'image à l'aide de quatre indices définissant une distribution uniforme des diverses positions des deux foyers (F₁, F₂), en partant sommairement de quatre boucles imbriquées:

Code Pascal :

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

1
2
3
4
FOR Kx1:= 0 TO Kmax DO
  FOR Ky1:= 0 TO Kmax DO
    FOR Kx2:= 0 TO Kmax DO
      FOR Ky2:= 0 TO Kmax DO <...>

au sein desquelles les coordonnées des foyers s'expriment simplement en fonction des dimensions (La, Ha) de l'image:

Code :

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

F1.x:= Round((La - 1) * Kx1 / Kmax) ; F1.y:= Round((Ha - 1) * Ky1 / Kmax) ... etc

Afin d'éviter tout doublon résultant d'un échange entre les deux foyers, on conviendra que (F₁) représente le foyer le plus à gauche: (Kx1 < Kx2)
ou sinon le plus bas, hors superposition: (Ky1 ≤ Ky2 si Kx1 = Kx2).

L'énumération de tous les couples résultera par conséquent des instructions plus strictes:

Code Pascal :

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

1
2
3
4
5
6
7
FOR Kx1:= 0 TO Kmax DO
  FOR Ky1:= 0 TO Kmax DO
    FOR Kx2:= Kx1 TO Kmax DO
      BEGIN
        IF (Kx2=Kx1) THEN Kini:= Ky1 ELSE Kini:= 0;
        FOR Ky2:= Kini TO Kmax DO < ... >        
      END

mais fournira cependant un nombre rapidement croissant de combinaisons; le nombre d'ellipses ainsi envisagées est:

Nell = (1/2)*(K_max + 1)²(1 + (K_max + 1)²) ;

soit par exemple:

Code :

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

1
2
Kmax =   10		15		20		30		40		50
Nell =   7381		32896		97461		462241		1413721		3383901

2°) Pour chaque paire de foyers, lister dans un tableau de même type l'aire de l'ellipse passant par le point correspondant (M_j):

Code Pascal :

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
 VAR Lst_A: Tab_V;
 ... / ...
 FUNCTION Aire_E(F_1, F_2, V_M: Ve_2E): Z_32;
   VAR Df2, X2, Y2: Z_32; D_2, D_a, D_b, p, q, r, Sd, Se: Reel;
   BEGIN
     X2:= Sqr(F_2.x - F_1.x); Y2:= Sqr(F_2.y - F_1.y); Df2:= X2 + Y2;
     IF (Df2=0) THEN BEGIN
                       X2:= Sqr(V_M.x - F_1.x); Y2:= Sqr(V_M.y - F_1.y);
                       // D_a = D_b = Sqrt(X2 + Y2) ; Aire = Pi * D_a * D_b
                       Se:= 4 * (X2 + Y2)
                     END
                ELSE BEGIN
                       X2:= Sqr(V_M.x - F_1.x);  Y2:= Sqr(V_M.y - F_1.y);
                       p:= Sqrt(X2 + Y2);
                       X2:= Sqr(V_M.x - F_2.x);  Y2:= Sqr(V_M.y - F_2.y);
                       q:= Sqrt(X2 + Y2);        Sd := p + q; D_2:= Sqr(Sd);
                       // D_a:= 0.5 * Sd; D_b:= 0.5 * r; Se = 4 * Da * Db
                       r:= Sqrt(Abs(D_2 - Df2)); Se:= Sd * r
                     END;
     Result:= Round(Se)       // (Aire = Pi * D_a * D_b = (Pi/4) * Sd * r
   END;
... / ...
FOR j:= 1 TO N_Point DO
                       BEGIN
                         Lst_A[j].x:= j;
                         Lst_A[j].y:= Aire_E(Vf1, Vf2, Lst_P[j])
                       END;

Toutes les données de l'ellipse se déduisent des positions des foyers, ainsi que de celle du point considéré; la courbe est en effet caractérisée par une constante, la somme des distances: Sd = MF₁ + MF₂;
on obtient pour les demi-diamètres et la surface:

a = CA = Sd/2 ; b = CB = (1/2)(Sd² - F₁F₂²)^1/2 ; Aell = π*a*b .

Comme seule importe la comparaison des résultats, on allégera les calculs en s'en tenant au produit

Sell = 4*a*b = Sd * (Sd² - F₁F₂²)^1/2 .

3°) Trier chaque tableau selon l'ordre croissant des surfaces; examiner la surface de la plus grande ellipse provisoirement retenue, de rang

Limite = Round(70%*100) = 70 (pour l'exemple choisi),

retenir la valeur (ainsi que le point correspondant) quand celle-ci apparaît inférieure à celles précédemment obtenues:

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
                     IF (Smin>Lst_A[Limite].y)
                       THEN BEGIN
                              Smin:= Lst_A[Limite].y;
                              WITH Ell_1 DO BEGIN
                                              E(0114);   We(23, L2, Kx1, o);
                                              Write(Ky1:o, Kx2:o, Ky2:o);
                                              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];
                                              Write(Smin:9, Rn:o); E(0011)
                                            END
                            END

Toutes les données utiles seront mémorisées dans la variable Ellipse:

Code Pascal :

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

1
2
3
 TYPE T_Ellipse = RECORD  M0, K1, K2, F1, F2: Ve_2E; Rn: Word; Sd, Se: Z_32  END;
 
 VAR Ellipse: T_Ellipse;

Le tracé final de l'ellipse fait appel au repère local construit sur les axes de symétrie:

Xell = a*Cos(u) ; Yell = b*Sin(u) ; u = 2π(k/Npell) avec 0 ≤ k < Npell et Npell = 1000 ;
Xm = Round(Xell * CosT - Yell * SinT + Xcen) ;
Ym = Round(Xell * SinT + Yell * CosT + Ycen) .

Voir la figure plus haut.

[wiwaxia]

Il faut 5 points pour définir une conique. Cela peut nécessiter 5 points frontière.

Un point de l'ellipse et ses deux foyers suffisent à la construction de la courbe (voir les calculs précédents).

@wiwaxia : ta méthode, comme la mienne, s'avère bancale. Le fait qu'un seul point serve de contact à l'extérieur montre que le résultat est faux. Il suffirait de déplacer légèrement le centre pour réduire l'ellipse. Les points de contact extérieurs sont forcément multiples.

C'est une affaire de résolution, liée au pas de la grille constituée de l'ensemble des positions envisagées pour (F₁, F₂).

Et là, je soupire.

On respire profondément, sans angoisse parce que l'algorithme s'avère lourd à mettre en œuvre, le temps d'exécution apparaissant proportionnel au produit Kmax⁴*N_Point² .
Il n'est pas possible, par exemple, de parvenir directement à la localisation optimale des foyers au pixel près. On peut toutefois procéder par approximations successives, en utilisant à la nouvelle étape une grille de pas plus faible centrée sur les positions précédemment obtenues.
Dans le cas d'un nuage de 100 points, deux étapes devraient suffire; mais s'il y en a 1000, elles seraient forcément plus nombreuses.

Autre idée: considérer la surface de la plus petite ellipse comme une fonction des variables: Sell = F(X1, Y1, X2, Y2). La recherche du minimum dans un espace à 4 dimensions peut être envisagée par la comparaison de la valeur centrale à celle des 16 points voisins.
Cependant, pour éviter le piège des minimums secondaires, il vaut mieux partir de la solution approchée donnée par le programme précédent.

Et là, je n'ai vraiment pas le temps de m'investir dans un nouveau programme. Moi aussi, je soupire

[wiwaxia]

Ici, on a un domaine de validité : 'sont valides toutes les ellipses qui contiennent au moins n points', et une fonction objectif : 'minimiser la surface'.
Et forcément, la solution optimale est à la frontière du domaine de validité.
Et quelque chose me dit que les outils divers et variés sont efficaces pour rechercher un optimum quand cet optimum est au 'centre' du domaine de validité, mais que ça marche moins bien dès qu'on doit chercher à la frontière.
Ici, non seulement la solution optimale passe par au moins 2 points de notre jeu de données, mais je pense qu'elle passe même par 3 points. On est donc plus qu'à la frontière.
... / ...

Le blocage de la discussion vient de l'a priori selon lequel la recherche de l'ellipse devrait se faire exclusivement sur quelques points du nuage, et rien d'autre. Pourquoi ?
Outre le fait qu'il faut prendre 4 points (déjà 3 dans le cas d'un cercle) et que la construction de la courbe est beaucoup plus difficile, qu'est-ce qui interdit la recherche de plus en plus fine de la position des foyers ? On sait déjà qu'ils se trouvent nécessairement à l'intérieur du nuage, et qu'on en déduit aisément les caractéristiques recherchées ... Pourquoi s'en priver ?

[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.
A priori la complexité est comparable à notre exercice. Mais je pense qu'en fait ce nouveau problème est beaucoup plus simple. Des algorithmes à base de descente du gradient / monte-carlo vont être efficaces pour ce nouvel exercice, parce que l'ellipse solution n'est pas à la frontière du domaine de validité.