[java] Triangulation de Delaunay (incrémentale)

[Rumpel]

Code :

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
	bool QuadEdge::inCircle(Point& a, Point& b, Point& c, Point& d) {
		/*
		 if "d" is strictly INSIDE the circle, then
 
		     |d² dx dy 1|
             |a² ax ay 1|
		 det |b² bx by 1| < 0
		     |c² cx cy 1|
 
		*/
		long a2 = a.x*a.x + a.y*a.y;
		long b2 = b.x*b.x + b.y*b.y;
		long c2 = c.x*c.x + c.y*c.y;
		long d2 = d.x*d.x + d.y*d.y;
 
		long det44 = 0;
		det44 += d2  * det33( a.x,a.y, 1,  b.x,b.y, 1,  c.x,c.y, 1 );
		det44 -= d.x * det33( a2 ,a.y, 1,  b2 ,b.y, 1,  c2 ,c.y, 1 );
		det44 += d.y * det33( a2 ,a.x, 1,  b2 ,b.x, 1,  c2 ,c.x, 1 );
		det44 -= 1   * det33( a2,a.x,a.y,  b2,b.x,b.y,  c2,c.x,c.y );
 
		if (det44<0) return true;
		return false;
	}
 
	long QuadEdge::det33(long m0, long m1, long m2, long m3, long m4, long m5, long m6, long m7, long m8) {
		long det33=0;
		det33 += m0 * (m4*m8 - m5*m7);
		det33 -= m1 * (m3*m8 - m5*m6);
		det33 += m2 * (m3*m7 - m4*m6);
		return det33;
	}

Et moi pour ces points il me renvoie det -471 078 406 donc true...

[pseudocode]

J'obtiens cela:

a2 = 9325
b2 = 133832
c2 = 175012
d2 = 7642

det44 = 0 + (412362320) - (5266563450) + (-222207280) - (-8900297300)

--> det44 = 3823888890

[Rumpel]

Ah oui, mais en java les long sont sur 8 octets et en C++ les long sont seulement sur 4 octets, j'ai dépassé les 2 147 483 647 de valeur max...

Comment vais je résoudre ce problème ...
(Merci pseudoCode)

[pseudocode]

Envoyé par Rumpel

Ah oui, mais en java les long sont sur 8 octets et en C++ les long sont seulement sur 4 octets, j'ai dépassé les 2 147 483 647 de valeur max...

Arf ! essaye "long long", ou au pire "float" / "double" ?

[Rumpel]

Merci pseudocode, je connaissais pas les long long.
Bon tout semble rentrer dans l'ordre maintenant

[Rumpel]

Ah oui encore une chose,

Code :

this.quadEdge.remove(e.onext());

Ne faut il pas enlever aussi l'arête symétrique ?
Ca me paraitrait logique...

[pseudocode]

Envoyé par Rumpel

Ah oui encore une chose,

Code :

this.quadEdge.remove(e.onext());

Ne faut il pas enlever aussi l'arête symétrique ?
Ca me paraitrait logique...

Heu oui. Tout à fait.

Gros oubli dans mon code que je vais aller corriger de ce pas. Merci !

[souviron34]

Envoyé par pseudocode

Moi ca m'intéresse, spécialement la partie construction de la liste des triangles à partir des quadedges.

Dans l'implémentation Java que j'ai postée, cette partie là est une création originale. Je me suis assez cassé la tête dessus pour avoir quelque chose d'efficace, mais je n'ai jamais démontré rigoureusement son exactitude. J'ai juste lancé plusieurs milliers de triangulations et comparé les résultats avec d'autre implémentations.

Pour info, j'ai posté le code C ici

[pseudocode]

Envoyé par souviron34

Pour info, j'ai posté le code C ici

Oui, j'ai vu. Tu travailles directement sur les triangles donc pas de problème pour les énumérer . Moi je travaille avec la structure QuadEdge donc l'extraction des triangles/régions nécessite une opération supplémentaire.

[souviron34]

Envoyé par pseudocode

Oui, j'ai vu. Tu travailles directement sur les triangles donc pas de problème pour les énumérer . Moi je travaille avec la structure QuadEdge donc l'extraction des triangles/régions nécessite une opération supplémentaire.

quelle horreur les maths

( c'est vrai, je vais toujours au plus simple... Dès que j'ai besoin d'une notion supplémentaire par rapport au problème, ça m'énerve et je tente d'éviter)



PS : par contre j'ai réfléchi au handicap que j'avais trouvé, et en fait il est général : dès qu'un triangle est très aplati, la méthode IsInCircle ne suffit pas... J'avais fait un "patch" qui était de vérifier par rapport aux sommets. On pourrait aussi vérifier en prenant le cercle dont le centre est le plus proche du point considéré, ou de rayon le plus petit, ou dont la "corde" est la plus petite...

PS2 : toutes les méthodes que j'ai vu se basent sur le fait que l'enveloppe extérieure est convexe. J'avais là aussi développé un "patch" (en utilisant ma méthode d'enveloppe non-convexe) pour des enveloppes concaves, mais je ne sais pas si il existe des méthodes "générales" traitant ce problème..

[Rumpel]

Hello,

Je regarde ça en fin de semaine, histoire de ne pas me mélanger avec mon propre code . Merci à toi souviron, je te dirais ce que j'en pense bientôt...

[Rumpel]

Code :

		if ( p.x<=bbox.minx || p.x>=bbox.maxx || p.y<=bbox.miny || p.y>=bbox.maxy ) {

C'est mieux, parce que si un point arrive juste sur la ligne de la Bounding box, tu vas avoir mal, c'est arrivé juste pendant ma démo

Sinon encore des problème avec inCircle
a (50,90)
b (70,50)
c (50,10)
d (60,30)

A vu de nez d est dans le cercle, pourtant ça renvoit faux...
Et je suis passé en double cette fois ci...

[pseudocode]

Envoyé par Rumpel

Code :

		if ( p.x<=bbox.minx || p.x>=bbox.maxx || p.y<=bbox.miny || p.y>=bbox.maxy ) {

C'est mieux, parce que si un point arrive juste sur la ligne de la Bounding box, tu vas avoir mal, c'est arrivé juste pendant ma démo

Dans mon cas, la surrounding-box est beaucoup plus grande que la Bounding box. Donc pas de problèmes.

Sinon encore des problème avec inCircle
A vu de nez d est dans le cercle, pourtant ça renvoit faux...
Et je suis passé en double cette fois ci...

La routine "inCircle" ne fonctionne que si les points (a,b,c) sont dans le sens trigo:

Code java :

1
2
if (QuadEdge.isAtRightOf(e,t.dest()) && QuadEdge.inCircle(e.orig(), t.dest(), e.dest(), p)) {
	// flip triangles

[Rumpel]

Ah ok merci Pseudocode, j'avais pas noté, ça...

[Rumpel]

Dans locate :

Code :

1
2
			else
				return e;

J'ai modifié comme ça :

Code :

1
2
3
4
5
6
7
8
		else
		{
			if (QuadEdge::isOnEdge(*(e->getOnext()),p))
				e = e->getOnext();
			if (QuadEdge::isOnEdge(*(e->dprev()),p))
				e = e->dprev();
			return e;
		}

A mettre en relation avec insertPoint

Code :

		if (QuadEdge::isOnEdge(*e,*p)) {

Sinon p peut être sur une des 3 arêtes du triangle, pas forcément sur e...

[pseudocode]

heu... si tu le dis. J'avoue que je ne vois pas trop où est le problème dans le code original...

locate() renvoie juste l'un des 3 quadedges du triangle contenant le point p.

Ou alors je n'ai pas compris le problème.

[Rumpel]

Oui, mais c'est pour traiter le cas où le point P est sur l'une des 3 edges du triangle (sur le triangle, pas à l'intérieur strictement). Forcément dans ce cas, tu dois détruire l'edge sur lequel P est tombé. Mais dans ton code tu as 2 chances sur 3 de te tromper
Toi ça ne doit pas t'arriver souvent, mais moi dans mon problème j'ai souvent des points qui doivent tomber pile sur les edges...

[Commander Keen]

Cela a peut-être été évoqué avant, mais on peut remplacer la ligne
"if (QuadEdge.isOnEdge(e,p)) {" de la fonction insertPoint

par "if ((p.x-e.orig().x)*(p.y-e.dest().y)==(p.y-e.orig().y)*(p.x-e.dest().x)){"
et supprimer les 3 fonction nommées ci-dessus. J'ai fait des tests avec des 200, 500 et 1000 points aléatoires, et ça a l'air de fonctionner.
En plus, cela retire le "d>0.01"

La relation ci-dessus signifie que lorsque A, B et C sont sur la même droite, le produit vectoriel de AB et AC vaut zéro. On augmente artificiellement l'espace d'une dimension pour calculer le produit vectoriel, mais on peut fixer la composante z de chaque point à 0. On ne s'intéresse qu'à "B appartient au segment [AC]", par exemple, mais dans notre cas, on peut se contenter de "B est sur la droite (AC)" car locate nous confirme déjà que si "B est sur (AC)" alors il est forcément sur le segment.

[bastiien]

Bonjour pseudocode,
J'essaie d'utiliser ton algo java pour déterminer la zone la plus proche des villes sur la France. J'ai du l'adapter un peu pour avoir plus de précision donc j'ai passé tous les Point en Point2D.Double et tous les int en double.

J'obtiens 2 problèmes :
- Le premier, un peu incompréhensible, c'est que sur certains ajouts de points, l'aglo reste indéfiniment bloqué dans l'insertPoint car les valeurs next des QuadEdge rebouclent...

- Le 2ème problème, je vais l'illustrer sur les villes du Rhône. Sur ce département, aucun ajout de point ne bloque par contre j'obtiens quelques zones complètement incohérentes puisqu'elles se superposent à d'autres zones correctes.
J'ai joint l'image obtenu sous google earth grâce à ma génération.
J'ai également joint à ce message le positionnement de mes villes.

Pensez vous que le passage en Point2D puisse provoquer ces erreurs ? Est-ce qu'il pourrait y avoir une autre anomalie ?

Merci
Bastien

[bastiien]

Rebonjour,

J'avais pas fait attention au dernier poste de Commander Keen et je viens d'ajouter sa modification à l'algo... Et magie, plus de blocage dans la génération et plus d'erreur sur celle-ci non plus!

Ha la journée commence bien!

Bastien