Cercles et contour

[lyceen95]

Oui et non. Je disais précédemment qu'il faut 'imposer' un sens de rotation.
- Tu as un jeu de n points.
- Etape n°1 : déterminer si les points tournent dans le sens des aiguilles d'une montre, ou dans le sens opposé.
- Etape n°2 : pour tous les triplets de 3 segments, rechercher le cercle qui va bien. Le cercle qui va bien, c'est celui qui est 'à l'intérieur du polygone' : si les points du polygone tournent dans le sens des aiguilles d'une montre, c'est le cercle qui est sur la droite des 3 segments.

Pour optimiser le truc, et trouver assez vite des cercles qui devraient être optimums, on peut éventuellement classer les segments. Statistiquement, les 3 segments qui vont être retenus pour le cercle optimal, ce sont 3 segments relativement longs.
On peut donc classer les segments du plus long au plus court. On peut aussi rechercher un point proche du centre de gravité du polygone, et mesurer l'écart angulaire de chaque segment (je ne sais pas si je suis clair), et les 3 segments qui vont nous intéresser seront très probablement parmi les segments qui ont un """écart angulaire""" les plus élevés.

A suivre ...

[sylvain231]

et comment tu détermine le sens de tour de points, sachant que le contour n'est pas convexe ?
qu'est ce que tu appelle le cercle qui est sur la droite des 3 segments ?
qu'est ce que tu appelle écart angulaire d'un segment au centre de gravité ?
désolé de toutes ces questions mais j'essaie de comprendre

[lyceen95]

Si le polygone n'est pas convexe, tout ce que j'ai pu raconter jusque là tombe à l'eau, ou presque. En effet, imaginons une forme en étoile plus ou moins régulière. Je disais que le plus grand cercle inscrit serait tangent à 3 segments, et là, on se retrouve avec une configuration où le plus grand cercle inscrit s'appuie sur 3 sommets.
Du coup, on parlait d'une boucle triple sur tous les segments, et ça devient une boucle triple sur tous les segments, et tous les sommets '<orientés vers l'intérieur>'.

Comment savoir si les points sont orientés dans le sens des aiguilles d'une montre, ou dans l'autre sens.
On a n sommets. La somme des angles fait (n-2)*180, ou bien (2-n)*180.
Selon qu'on obtient (n-2)*180 ou (2-n)*180, ça veut dire que les points sont dans un sens ou dans l'autre.
Ca reste valable même si le polygone n'est pas convexe.

Le cercle sur la droite des 3 segments... Effectivement, en relisant, je me suis dit que ça ne voulait rien dire.... mais si. Ici, le mot droite est utilisé dans le sens 'droite ou gauche'.
si les points du polygone tournent dans le sens des aiguilles d'une montre, c'est le cercle qui est sur la droite (= quand on tourne les yeux vers la droite) des 3 segments. Et si les points du polygone sont dans le sens inverse des aiguilles d'une montre, le cercle qui nous intéresse c'est celui qui est sur la gauche des 3 segments.

Ecart angulaire : On a notre polygone avec par exemple 100 points, on a son centre de gravité (exact ou approximatif, peu importe).
Pour tout couple de points consécutifs, on peut calculer l'angle
En moyenne, cet angle vaut 3.6° (cf PS), puisqu'on a 100 points et que la somme donne 360°. Parfois, si le polygone n'est pas convexe, cet angle peut être négatif. Cet angle , c'est ce que j'appelais l'écart angulaire. Et plus cet écart angulaire est grand, plus on a de chance de retrouver ce segment dans les segments tangents au cercle optimal.
Par exemple, avec ces 4 points (0,0) (2,70)(-62,70)(-60,0), les 4 angles sont proches de 90°, le segment (-60,0)(0,0) est celui qui donne l'angle le plus petit, et les 3 autres sont effectivement tangents au cercle recherché.

PS : si les points sont dans le sens des aiguilles d'une montre, cette moyenne n'est pas de +3.6°, mais de -3.6° ... Inversons tous les points pour retomber sur nos pieds.

[sylvain231]

On a beaucoup discuté et je suis un peu perdu, peux-tu stp me redire l'algo complet à jour ?
D'ailleurs on ne sait toujours pas comment trouver si un cercle est à l'intérieur d'un polygone car la solution trop simpliste du forum est fausse

[sylvain231]

j'ai trouvé ça :
https://stackoverflow.com/questions/427 ... ex-polygon

[sylvain231]

grâce à ce forum j'ai fait ce code en C++ avec la librairie OpenCV et ça marche et est quasi instantané :

Code: Tout sélectionner

Circle cercleCirconscrit(std::vector<cv::Point> contour) {
   /*bool br;
   Circle res = Circle(0, Point(-1, -1));
   double minRadius = 9999999999;
   for (int i = 0; i < contour.size(); i += precision) {
      cout << "i=" << i + 1 << "/" << contour.size() << endl;
      for (int j = i + 1; j < contour.size(); j += precision) {
         for (int k = j + 1; k < contour.size(); k += precision) {
            br = false;
            Circle c(contour[i], contour[j], contour[k]);
            //vérif si tous les points du contour sont dans le cercle
            Point P = c.getCenter();
            double R = c.getRadius();
            if (R < minRadius) {
               for (int l = 0; l < contour.size(); l += precision) {
                  if (distanceP(P, contour[l]) > R)
                  {
                     br = true;
                     break;
                  }
               }
               if (!br) {
                  res = c;
                  minRadius = R;
               }
            }
         }
      }
   }*/
   Point2f center;
   float radius;
   minEnclosingCircle(contour, center, radius);
   return Circle(radius,center);
}

Circle cercleInscrit(std::vector<cv::Point> contour,Size imgSize) {
   Mat dist_map;
   double* minVal=new double;
   double* maxVal=new double;
   Point* minLoc=new Point;
   Point* maxLoc=new Point;
   //création d'une image binaire avec le contour rempli de blanc
   Mat src_map = Mat::zeros(imgSize, CV_8UC1);
   //remplissage
   vector<vector<Point>> polygones;
   polygones.push_back(contour);
   fillPoly(src_map, polygones, Scalar(255));
   imwrite("src_map.jpg", src_map);
   distanceTransform(src_map, dist_map, DIST_L2, DIST_MASK_PRECISE);
   minMaxLoc(dist_map, minVal, maxVal, minLoc, maxLoc);
   return Circle(*maxVal, *maxLoc);
}

[sylvain231]

l'algorithme de transformation de distances utilisé est décrit dans cet article : http://www.theoryofcomputing.org/articl ... 08a019.pdf

[sylvain231]

mais cet article est très compliqué, la notion de transformation de distance est plus claire dans cet article de Wikipédia : https://en.wikipedia.org/wiki/Distance_transform

[lyceen95]

Je n'ai pas d'algo complet à jour. D'autant plus que je partais plutôt sur un polygone convexe.
L'idée était celle exposée au début : Dans le cas d'un polygone convexe, le cercle solution est tangent à 3 segments.
Du coup, on recense tous les triplets de 3 segments, et pour chaque triplet, on recherche le cercle qui est tangent à ces 3 segments, et qui est dans la portion de plan '<Intérieur du polygone>'.
On vérifie que ce cercle est bien intérieur au polygone (aucun sommet n'est à l'intérieur du cercle), et si ce cercle est valide, on regarde si ce cercle est mieux que le cercle actuellement optimum.
Si ce cercle est mieux, on le retient, jusqu'à ce qu'on trouve mieux.

On pourrait s'en tenir là, mais c'est en o(n^4). On va donc chercher à optimiser, en recherchant les segments les plus prometteurs.
Si on sait identifier les segments les plus prometteurs, on va trier tous les segments du plus prometteur au moins prometteur, et on va faire une boucle comme celle ci :

Code: Tout sélectionner

for ijk = 6 to 3n-3
  for i=0 to least(ijk,n-2)
    for j = i+1 to least(ijk,n-1)
      k=ijk-i-j 
      if k>j and k <= n
        // on regarde les 3 segments i,j,k
      end if
    next j
  next i
next ijk

Avec cette façon 'particulière' de boucler sur tous les triplets, on évite les segments peu prometteurs, on les traitera tout à la fin, sauf si on a des raisons de ne pas les traiter.

Et on peut trouver des raisons de ne pas les traiter.
Revenons au polygone initial (supposé convexe), et prenons 4 points consécutifs ABCD. Pour le segment BC, en fonction de la longueur BC et des 2 angles ABC et BCD, on sait trouver le rayon maximum d'un cercle tangent à BC, et qui ne coupe pas les demi-droites BA ni CD. Et donc on sait déterminer un RayonMax pour chaque segment. On peut donc trier les segments selon ce RayonMax décroissant, et faire ainsi :

Code: Tout sélectionner

RayonBest=0
for ijk = 6 to 3n-3
  for i=0 to least(ijk,n-2)
   if RayonMax(i)> RayonBest then
   for j = i+1 to least(ijk,n-1)
      if RayonMax(j)> RayonBest then
      k=ijk-i-j 
      if k>j and k <= n
        if RayonMax(k)> RayonBest then
        // on regarde les 3 segments i,j,k
        // et si ce cercle est valide et mieux que RayonBest, on mémorise ce cercle.
      end if
      end if
      end if
    next j
    end if
  next i
next ijk

Dès qu'on a trouvé un cercle valide à peu près optimum, tous les segments qui ne peuvent pas donner un rayon supérieur à ce rayon optimum sont exclus de l'analyse. Et on s'est débrouillé pour que tous ces segments ne soient pas analysés au début.

[sylvain231]

tu n'as pas regardé mes messages précédents !
je t'explique l'algo :
on dessine le polygone rempli de blanc sur une image noire.
on calcule la carte des distances, en fait chaque pixel de cette carte contient la distance du pixel blanc correspondant sur l'image de départ au plus près des pixels noirs toujours sur l'image de départ
ensuite on cherche le pixel de cette carte ayant la plus grande valeur, c'est le centre du cercle, on prend alors la valeur de ce pixel : c'est le rayon du cercle et voilou
évidemment cette solution est approximative car un pixel c'est un quantum et n'est pas continu mais c'est largement suffisant pour mon application

[sylvain231]

j'ai un peu élaboré mon algo de cercle inscrit pour éviter les fuites de mémoire :

Code: Tout sélectionner

Circle cercleInscrit(std::vector<cv::Point> contour,Size imgSize) {
   Mat dist_map;
   double* minVal=new double;
   double* maxVal=new double;
   Point* minLoc=new Point;
   Point* maxLoc=new Point;
   //création d'une image binaire avec le contour rempli de blanc
   Mat mask = Mat::zeros(imgSize, CV_8UC1);
   //remplissage
   vector<vector<Point>> polygones;
   polygones.push_back(contour);
   fillPoly(mask, polygones, Scalar(255));
   //imwrite("mask.jpg", mask);
   distanceTransform(mask, dist_map, DIST_L2, DIST_MASK_PRECISE);
   minMaxLoc(dist_map, minVal, maxVal, minLoc, maxLoc);
   Circle res= Circle(*maxVal, *maxLoc);
   free(maxVal);
   free(maxLoc);
   free(minVal);
   free(minLoc);
   return res;
}

[lyceen95]

Je t'ai exposé l'algorithme que j'avais en tête dans l'après-midi, sur la base d'un polygone convexe. Cet algorithme donne un point extrêmement précis (valeur réelle, et non approchée) et avec les optimisations proposées, on devrait avoir une performance raisonnable.

[sylvain231]

oui mais il nous reste des inconnues :
comment savoir si un cercle est dans un polygone ?
de plus je suppose que n est le nombre de points dans ton algo non ?
je ne comprends pas ton algo de tri des segments, tu sembles les trier par indices alors que je croyais qu'on les triait par taille ? La distance entre deux points consécutifs n'est en effet pas constante
enfin "// on regarde les 3 segments i,j,k" est un peu vague, les regarder c'est-à-dire ?

[lyceen95]

J'avais tapé un très long message, mais je l'ai perdu au moment de l'envoyer.... Je remets ce ue je disais sur une des questions.

Comment sait-on si un cercle est à l'intérieur d'un polygone ?
Ca revient à se demander , pour chaque segment P1 P2, ce segment P1 P2 est-il à l'extérieur du cercle de centre P0 et de rayon R ?
Si une des 2 extrémités est dans le cercle, c'est mort, fin du traitement.
On recherche le point H, projection de P0 sur la droite P1 P2.
On a 4 équations :
H est sur la droite P1 P2, donc il existe un réel k tel que Hx = P1x+ k(P2x-P1x) et Hy = P1y+ k(P2y-P1y)
Et par ailleurs, H est sur la droite passant par P0, et perpendiculaire à P1P2, donc il existe un réel m tel que
Hx=P0x+m(P2y-P1y) et Hy=P0y-m(P2x-P1x)
On a donc ces4 équations, et 4 inconnues Hx,Hy,k et m. On peut résoudre ce système. La valeur de m ne nous intéresse pas vraiment.
La valeur de k nous intéresse. Si k n'est pas entre 0 et 1, ce point H est à l'extérieur du segment P1 P2. et donc pas de problème.
Si k est entre 0 et 1 (autrement dit si ce point H est entre P1 et P2), on calcule la distance P0 H, et si cette distance est inférieure à R, ça veut dire qu'une partie du segment P1 P2 est dans le cercle.

Si tu veux, je peux mettre ça en fonction dans la soirée. Ce ne sera pas du C mais du GO, mais tu sauras traduire.

[sylvain231]

oui je veux bien et pour mes autres questions STP ?

[lyceen95]

J'ai un jour de retard, normal

Voici un bout de code, je reviens en fin de soirée :

Code: Tout sélectionner

type ppoint struct {
   xx int
   yy int
}

func fdistance(P0 ppoint, P1 ppoint) float64 {
   return math.Sqrt(float64((P1.xx-P0.xx)*(P1.xx-P0.xx) + (P1.yy-P0.yy)*(P1.yy-P0.yy)))
}

func fDistanceProjection(P0 ppoint, P1 ppoint, P2 ppoint, rayon float64) bool {
   // On recherche la projection de P0 sur la droite P1 P2
   // Mais on ne s'intéresse à ce point H que s'il est entre P1 et P2
   // Si H est en dehors du segment P1 P2, ou si la distance de P0 à H est bien supérieure à rayon, alors return True.
   var H ppoint
   var k float64
   deltaX := float64(P2.xx - P1.xx)
   deltaY := float64(P2.yy - P1.yy)
   if deltaX == 0 && deltaY == 0 {
                // Les points P1 et P2 sont superposés.
      return true
   }
   k = (float64(P0.xx-P1.xx)*deltaX + float64(P0.yy-P1.yy)*deltaY) / (deltaX*deltaX + deltaY*deltaY)
   if k < 0 || k > 1 {
      // H n'est pas entre P1 et P2, donc ok.
      return true
   }
   H.xx = P1.xx + int(k*float64(P2.xx-P1.xx))
   H.yy = P1.yy + int(k*float64(P2.yy-P1.yy))
   fmt.Println(" Point H", H.xx, H.yy)
   if fdistance(P0, H) < rayon {
      return false
   }
   return true
}

func okCercle(rayon int, CCentre ppoint, P2 ppoint, P3 ppoint) bool {
   // Vérifie si le segment P1 P2 est bien tout entier à l'extérieur du cercle de centre CCentre et rde rayon rayon
   if fdistance(CCentre, P2) < float64(rayon) {
      //fmt.Println(" Point 1 dans le cercle ", rayon)
      return false
   }
   if fdistance(CCentre, P3) < float64(rayon) {
      //fmt.Println(" Point 2 dans le cercle ", rayon)
      return false
   }
   if fDistanceProjection(CCentre, P2, P3, float64(rayon)) == false {
      //fmt.Println(" Proj dans le cercle ", rayon)
      return false
   }
   return true
}
func main() {
   var P1, P2, P3 ppoint
   P1.xx = 0
   P1.yy = 0
   P2.xx = 800
   P2.yy = 100
   P3.xx = 600
   P3.yy = 630
   for i := 700; i < 900; i++ {
      ok := okCercle(i, P1, P2, P3)
      if ok {
         fmt.Println(i)
      }
   }
}