troussil / convergents Goto Github PK

Dans la fonction : Value operator()(const Point& aPoint)
On teste les trois positions du point par rapport à la droite : appartient, en dessous, en dessus.

Je me demande qu'elle est le comportement de la fonction si myX = 0 et que la droite est verticale ?

Souhaites-tu que j'implémente une position appartient, à gauche ou à droite de la droite si myX = 0.
(qui pose un problème symétrique si myY = 0)

Je me pose une question sur la relation entre les deux ensembles.

Dans le cas de l'enveloppe convexe, si plusieurs points appartiennent à un côté de l'enveloppe convexe, s'il appartiennent à une même droite, on choisit de ne prendre que le dernier point dans l'enveloppe.

Par contre, dans le cas de l'alpha-shape, si plusieurs points sont dans le même cas de figure, par définition, on prend l'ensemble des points de la droite dans l'alpha-shape.

Du coup, même avec un très grand rayon, si ce cas de figure se présente, les ensembles ne pourront pas être égaux.

Sometime the convex hull computation does not finish.

https://github.com/homeostasie/convergents/blob/alphaShape/testConvexHull.cpp

I modified the test files in order to recover the three vertices whicj lies on the circle and I make a few more (200) test in order to make the bug happen quicker.

I've got the bug for this value :

Disk[ Center : (25.5891, 130.797 ), Radius : 34.799 ] | Points : (26,96)(32,165)(-7,143) - First vertex : (25,96)

Disk[ Center : (-171.251, -11.357 ), Radius : 149.643 ] | Points : (-162,138)(-57,-108)(-64,93) - First vertex : (-172,-161)

Disk[ Center : (70.8515, 30.5562 ), Radius : 153.556 ] | Points : (-51,124)(197,-57)(71,-123) - First vertex : (70,-123)

...

Works fine with g++, but there are warnings with clang++ and segmentation fault with c++.

L'algo pour construire l'alpha-shade d'une droite discrète semble assez simple sur le papier :

Il est récursif sur un segment [a,b]

• On calcule le point de Bézout minimal du segment [a, b] : B1.

  Je pense qu'on peut l'extraire en prenant le dernier convergent en dessous du segment [a,b]. On doit pouvoir également passer par ailleurs, mais j'avoue ne pas avoir encore trop réfléchi à la question.

• On calcule le disque de rayon R = -1/alpha et passant par a et b. On appelle C son (futur) centre.

  Dans l'idée on cherche le milieu de [a,b] et on translate orthogonalement tant que [a,C] ne mesure pas R. Ça doit pouvoir se faire avec un peu de Pythagore.

• On regarde si B1 appartient à ce disque.

• * Si non : a et b sont des alpha-neighbours

• * Si oui : on fait tourner la procédure sur [a, B1] et [B1, b].

Je pense avoir réussi les différentes opérations.

Je ne sais pas trop si tu préfères que je fasse des pull request souvent quand une fonction est implémentée et fonctionnelle ou bien quand j'ai un peu plus avancé sur le sujet en terminant une "partie" plus complète de code.

Bonne soirée.

Je crois qu'il y a un cas particulier non traité.

Le cas par exemple du segment de droite formé par les points (0,0) et (1,3).

On trouve les convergents p-2 = (1,0); p-1 (0,1).
On trouve p0 = (1,3)

Du coup, on ne cherche jamais si on a un sommet faisant partie de l'alpha-shape.

Je suis entrain de retravailler tooldisplay pour pouvoir sortir plus facilement les différentes enveloppes alpha-shape (>0 et <0) et convexe par le même exécutable.

Je me pose la question de la pertinence de continuer à garder la méthode de Graham pour les cas convexe et alpha < 0.

Faut-il mieux que j'évite de multiplier les return.

Par exemple dans ce genre de cas :

if (k % 2 != 0)
{
    *aAlphaShapeHull++ = pConvM1;
    return pConvM1;
}
else
{
    *alphaShapeHull++ = pConv;    
    return pConv;
}

qui reste assez facilement échangeable par un code de ce genre :

if (k % 2 != 0)
{
    *aAlphaShapeHull++ = pConvM1;
    pNext = pConvM1;
}
else
{
    *alphaShapeHull++ = pConv;
    pNext = pConv;
}
return pNext;

Can we refactor the next method in OutputSensitiveAlphaShape in order to call the same procedure (dichotomic search step) when k is odd and k is even ?

J'ai posé sur le dépôt de ma branch convexHull le code.

https://github.com/homeostasie/convergents/tree/convexHull

J'ai commencé à rédiger correctement la partie test pour le calcul de l'enveloppe convexe.
Je teste dans un premier temps sur un triangle simple, puis sur le même triangle orienté différemment.

Puis,je complique un peu la chose en réalisant plus de tests à partir de 3 points tirés aléatoirement.

Néanmoins, j'ai un problème au niveau de la vérification.

Le premier point n'est pas toujours le même dans l'algo de Melkman. Le point de getbottom() n'est pas toujours présent.
Sinon, tous les suivants se retrouvent bien dans les deux listes.

Je ne sais pas si cela vient du fait que je ne suis pas sûr de réinitialiser correctement les paramètres d'itérations ou si cela provient d'ailleurs...

Exemple de résultat produit :

 - 3.4 - Convex hull on a random circle
Disk[ (-1.81902, 11.2607 ), 21.2765 ] | aStartingPoint : (-1,-10)
Get
(-1,-10), (4,-9), (7,-8), (9,-7), (13,-4), (16,0), (18,4), (19,7), (19,15), 
(18,19), (17,21), (15,24), (11,28), (8,30), (6,31), (2,32), (-6,32), 
(-9,31), (-13,29), (-16,27), (-18,25), (-20,22),(-22,18), (-23,13), 
(-23,10), (-22,5), (-19,-1), (-14,-6), (-8,-9), (-2,-10), 
Melkman's convex hull of the boundary
(-1,-10), (4,-9), (7,-8), (9,-7), (13,-4), (16,0), (18,4), (19,7), (19,15), 
(18,19), (17,21), (15,24), (11,28), (8,30), (6,31), (2,32), (-6,32), 
(-9,31), (-13,29), (-16,27), (-18,25), (-20,22), (-22,18), (-23,13), 
(-23,10), (-22,5), (-19,-1), (-14,-6), (-8,-9), (-2,-10), 
(4 tests passed / 4 tests)

 - 3.3 - Convex hull on a random circle
Disk[ (40.5263, 24.8421 ), 41.7037 ] | aStartingPoint : (48,-16)
Get
(48,-16), (52,-15), (58,-13), (60,-12), (66,-8), (74,0), (76,3), (79,9), 
(81,15), (82,21), (82,29), (81,34), (80,38), (78,43), (77,45), (73,51), 
(67,57), (61,61), (57,63), (51,65), (47,66), (34,66), (30,65), (24,63), 
(20,61), (14,57), (9,52), (6,48), (3,43), (1,38), (0,34), (-1,28), (-1,21),
 (0,16), (1,12), (2,9), (5,3), (8,-1), (15,-8), (21,-12), (23,-13), (29,-15), 
(33,-16), 
Melkman's convex hull of the boundary
(52,-15), (58,-13), (60,-12), (66,-8), (74,0), (76,3), (79,9), (81,15), 
(82,21), (82,29), (81,34), (80,38), (78,43), (77,45), (73,51), (67,57), 
(61,61), (57,63), (51,65), (47,66), (34,66), (30,65), (24,63), (20,61), 
(14,57), (9,52), (6,48), (3,43), (1,38), (0,34), (-1,28), (-1,21), (0,16), 
(1,12), (2,9), (5,3), (8,-1), (15,-8), (21,-12), (23,-13), (29,-15), (33,-16), 
(4 tests passed / 5 tests)

J'ai vu que l'alpha-shape dans le cas du disque (testAlphaShape.cpp) est comparé avec la méthode de graham en utilisant openTracking et openGrahamScan au lieu de closedTracking et closedGrahamScan.

C'est problablement un mauvais copié-collé de testAlphaShapeStraightLine.cpp.

J'ai commencé l'investigation de notre problème quand la taille des rayons augmentent. Pour un cercle de rayon 1025 (première valeur entière à posée problème), les calculs de Delta sont déjà assez grand :

Entrée : (0,-1025)(1,0)(0,1)-2101250 // 4305461250
Delta : 108332956690448384 ; aS, aD : (1,-1025)(0,1)4307562500
x1, x2 : 1103.32 , 946.68
0946(0,-1025)(1,0)(0,1)(1,-79)
Delta : 8650722642933972048 ; aS, aD : (0,-1024)(1,946)4070978560000
x1, x2 : 1.08323 , 1.08167

Soit un delta entre 10¹⁸ et 10¹⁹.

Je vais continuer mon investigation dans ce sens pour voir s'il est possible de séquencer un peu le calcul pour éviter d'avoir une variable aussi grande.

@troussil

Il faut reprendre étape par étape la méthode de calcul des alpha-shapes pour alpha > 0 :

• Utiliser closedGrahamScan
• Tester son utilisation.
• Le déplier en interne jusqu'à l'utilisation de la méthode updateConvexHull.

Enfin, dans un second temps, il faudra essayer de rendre la méthode incrémentale ;

• Utilisation des méthodes next et updateConvexHull pour chaque nouveau sommet.

Il va également falloir retravailler un peu l’exécutable toolAlphaShape pour prendre en compte les alpha-shape avec alpha > 0 pour les sorties.

Un travail pour refactoriser la fonction principale serait également bénéfique.

• compter le nombre de sommets de l'alpha-shape pour un 1/alpha variant en O(R).
• compter le nombre de sommets de l'enveloppe convexe (avec et sans les points sur le bord).
• générer 100 cercles avec centres rationels dans [0;1]x[0;1](paramétrisation a,b,c,d) pour des rayons variant de 2^5 à 2^15 (puis 2^20).

Il y a bien un problème pour alpha = 0.Tous les sommets sur le bord sont pris en compte.

Le problème provient entre autres de la ligne 220 du fichier inc/OutputSensitiveAlphaShape.h. Je ne sais pas encore trop comment régler le problème, mais je m'en occupe dès que je rentre des présentations.

/**
* If pConv lie on the alpha-shape, we reach a vertex of the
* convex-shape.
*/
            if (myShape(pConv) == 0)
            {...

Pour un exemple de problème mise en valeur : (46,-56), (48,-55), (50,-54), (52,-53), (54,-52), (56,-51) à la place de juste (46,-56), (56,-51).

-- Disk[ Center : (33.4208, -18.505 ), Radius : 39.5791 ] | 
Points : (1025,0)(0,1025)(-1025,0) - First vertex : (36,-58)

Radius predicate : Num2 / Den2 : 10/2

#1 - alpha-shape of the boundary using OpenGrahamScan
(36,-58), (42,-57), (46,-56), (56,-51), (65,-42), (67,-39), (70,-33),
 (72,-27), (72,-10), (70,-4), (67,2), (65,5), (56,14), (46,19), 
(42,20), (35,21), (31,21), (25,20), (21,19), (16,17), (11,14), 
(6,10), (5,9), (1,4), (-1,1), (-2,-1), (-4,-6), (-5,-9), (-6,-15), 
(-6,-22), (-5,-28), (-4,-31), (-2,-36), (-1,-38), (1,-41), (5,-46), 
(6,-47), (11,-51), (16,-54), (21,-56), (25,-57), (31,-58), (36,-58), 

#2 - alpha-shape
(36,-58), (42,-57), (46,-56), (48,-55), (50,-54), (52,-53), (54,-52), 
(56,-51), (57,-50), (58,-49), (59,-48), (60,-47), (61,-46), (62,-45), 
(63,-44), (64,-43), (65,-42), (67,-39), (68,-37), (69,-35), (70,-33), 
(71,-30), (72,-27), (72,-26), (72,-25), (72,-24), (72,-23), (72,-22), 
(72,-21), (72,-20), (72,-19), (72,-18), (72,-17), (72,-16), (72,-15), 
(72,-14), (72,-13), (72,-12), (72,-11), (72,-10), (71,-7), (70,-4), 
(69,-2), (68,0), (67,2), (65,5), (64,6), (63,7), (62,8), (61,9), (60,10), 
(59,11), (58,12), (57,13), (56,14), (54,15), (52,16), (50,17), (48,18), 
(46,19), (42,20), (35,21), (34,21), (33,21), (32,21), (31,21), (25,20),
(21,19), (16,17), (11,14), (6,10), (5,9), (1,4), (-1,1), (-2,-1), (-4,-6),
(-5,-9), (-6,-15), (-6,-16), (-6,-17), (-6,-18), (-6,-19), (-6,-20), (-6,-21),
(-6,-22), (-5,-28), (-4,-31), (-2,-36), (-1,-38), (1,-41), (5,-46), (6,-47),
(11,-51), (16,-54), (21,-56), (25,-57), (31,-58), (32,-58), (33,-58), 
(34,-58), (35,-58), (36,-58), 

Dans OutputSensitiveConvexHull, remplacer méthode .det par un prédicat exact (calculs intermédiaires avec des variables de type BigInteger).

Dans l'algorithme qui regarde si un rayon intersecte une droite, je suis amené à utiliser un point et son vecteur correspondant qui se trouvent dans R².

L'algorithme suit cette idée :

J'ai crée une droite passant par aS de direction aD.

• Je vérifie qu'elle ne soit pas parallèle et que aS n'appartienne pas à d pour pouvoir continuer.
• Je calcule le point d'intersection en résolvant le système linéaire à deux équations, deux inconnus. (j'ai pris en compte le cas ou la droite est parallèle).

Je note Pinter ce point. Il se trouve dans R².

• Je me demande si le rayon coupe la droite en regardant le signe du produit scalaire de aD et du vecteur formé par (Pinter - aS).
• S'il est positif et qu'ainsi le rayon coupe la droite, je normalise mon vecteur ainsi formé par aD et en retire les informations utiles : aQ et aC et je retourne true.

Difficultés

L'une des difficultés rencontrée est que le point d'intersection se trouve être dans le plan R² et non Z². Je n'ai pas vraiment réussi à ne pas m'en servir (sans pour autant rendre la chose complexe).
Du coup, je me retrouve avec un point dans R² puis un vecteur dans R²... De plus, je me vois obligé d'utiliser le produit scalaire et la normL2 pour les réels alors que je les ai créées pour les entiers...

Je me pose la question de savoir si j'aurai pu tirer partie plus facilement du fait d'avoir écrit le programme de manière générique et de ce que j'aurai dû faire pour utiliser ces fonctions dans R² de manière plus transparente.

Apporter une première pull request avec un squelette du rapport avec un plan détaillé et un nouveau patron pour éviter les pages blabnches.

Détaillé dans le sens, mettre tout ce qui va être dans le rapport.

• ExactRayIntersectableCircle
• InexactRayIntersectableCircle
• OutputSensitiveConvexHull ok
• IncrementalNegativeAlphaShape (actuellement OutputSensitiveAlphaShape)
• RecursiveNegativeAlphaShape
• BottomUpPositiveAlphaShape (actuellement PositiveAlphaShape)
• TopDownPositiveAlphaShape (non prioritaire)
• OutputSensitiveAlphaShape (chapeau combinant négatif, 0, positive alpha-shape)

Dans notre cas avec alpha négatif, je sais que les sommets de l'enveloppe convexe sont également des points de l'alpha-shape. Par conséquent, chercher l'alpha-shape d'un convexe revient à chercher l'alpha-shape des segments de droite reliant les sommets.

Je me pose quand même la question de savoir si ce résultat est une généralité ou si on peut également l'étendre a cas où alpha est positif.

• Tous les points à l'intérieur d'un cercle sont-ils également forcément dans l'apha-shape ?

Il faut permettre de choisir la forme intersectable souhaité : Cercle exact et inexact pour le moment. Cette dernière forme doit être testée.

Refactoriser ce qui peut l'être.

Convex Hull, Alpha-shape pour alpha <0 et >0.

Réduire la taille de la fonction toolmean.

Je me retrouve avec une erreur à la compilation que j'ai du mal à comprendre :

testConvergents.cpp:20:30: erreur: missing template arguments before ‘DroiteRatio’
testConvergents.cpp:20:30: erreur: expected ‘;’ before ‘DroiteRatio’
testConvergents.cpp:23:3: erreur: ‘DroiteRatio’ was not declared in this scope

Il n'arrive pas à créer l'objet RayIntersectableStraightLine DroiteRation(p1, p2).

Ai-je oublié quelque chose de particulier ?

Le fichier concerné étant : testConvergents.cpp. Il est disponible dans la branch dev de mon dépôt. Voici l'extrait probablement compromettant.

#include <iostream>
#include <vector>
#include <iterator>

#include "PointVector2D.h"
#include "RayIntersectableStraightLine.h"

template <typename Point, typename OutputIterator>
void geometricConvergents(const Point& direction,
                          OutputIterator res)
{
  Point pS(0,0);
  Point pm2(0,1);
  Point pm1(1,0);
  Point pconv;
  int qk;
  RayIntersectableStraightLine DroiteRatio(pS, direction);

  DroiteRatio.dray(pm2, pm1, qk, pconv);

  while( direction != pconv )
 {
pm2 = pm1;
pm1 = pconv;
*res++ = pconv;
DroiteRatio.dray(pm2, pm1, qk, pconv);
 }
}