Recherche des triangles dans lesquels se trouvent les points

16

Supposons que j'ai un maillage 2D constitué de triangles non chevauchants , et un ensemble de points { p i } M i = 1N k = 1 T K . Quelle est la meilleure façon de déterminer dans quel triangle chacun des points se trouve?{Tk}k=1N{pi}i=1Mk=1NTK

Par exemple, dans l'image suivante, nous avons , p 2T 4 , p 3T 2 , donc je voudrais une fonction f qui retourne la liste f ( p 1 , p 2 , p 3 ) = [ 2 , 4 , 2 ] .p1T2p2T4p3T2ff(p1,p2,p3)=[2,4,2]

entrez la description de l'image ici

Matlab a la fonction pointlocation qui fait ce que je veux pour les maillages Delaunay, mais elle échoue pour les maillages généraux.

Ma première pensée (stupide) est que, pour tous les nœuds , parcourez tous les triangles pour savoir dans quel triangle p i se trouve. Cependant, cela est extrêmement inefficace - vous pourriez avoir à parcourir chaque triangle pour chaque point, donc cela pourrait prendre du travail O ( N M ) .pipiO(NM)

Ma prochaine pensée est, pour tous les points , de trouver le nœud de maillage le plus proche via la recherche du plus proche voisin, puis de regarder à travers les triangles attachés au nœud le plus proche. Dans ce cas, le travail serait O ( a M l o g ( N ) ) , où a est le nombre maximal de triangles attachés à n'importe quel nœud du maillage. Il y a quelques problèmes résolubles mais ennuyeux avec cette approche,piO(aMlog(N))a

  • Cela nécessite la mise en œuvre d'une recherche efficace du plus proche voisin (ou la recherche d'une bibliothèque qui en dispose), ce qui pourrait être une tâche non triviale.
  • Cela nécessite de stocker une liste des triangles attachés à chaque nœud, pour lequel mon code n'est actuellement pas configuré - pour le moment, il n'y a qu'une liste de coordonnées de nœud et une liste d'éléments.

Dans l'ensemble, cela semble inélégant, et je pense qu'il devrait y avoir une meilleure façon. Cela doit être un problème qui se pose souvent, donc je me demandais si quelqu'un pouvait recommander la meilleure façon d'aborder la recherche des triangles dans lesquels se trouvent les nœuds, théoriquement ou en termes de bibliothèques disponibles.

Merci!

Nick Alger
la source

Réponses:

16

La méthode habituelle de saut de bord aléatoire devrait fonctionner. Fondamentalement, commencez par n'importe quel triangle du maillage, puis déterminez de quels bords se situe le point cible du côté opposé. C'est-à-dire, déterminer lequel des bords, lorsqu'il s'étend sur une ligne, sépare le point de l'intérieur du triangle. Lorsqu'il y a deux possibilités, choisissez-en une au hasard, considérez le triangle adjacent à cette arête partagée et répétez. La randomisation devrait faire converger cette méthode avec la probabilité 1 pour les triangulations de Delaunay, et je ne vois pas de raison pour laquelle cela ne fonctionnerait pas pour les triangulations arbitraires.

Edit : je dois ajouter que le saut de bord doit être avec une constante raisonnable pour un seul point, donc ce serait O ( M log N ) pour M points. Cependant, si vous triez vos points par localité (comme en utilisant d'abord un ordre de courbe de Hilbert), vous pouvez initialiser chaque nouvelle requête avec le triangle de la requête précédente, pour réduire davantage le temps d'exécution (je ne suis pas un théoricien CS donc je peux '' t vous dire ce que le big-O serait là).O(JournalN)O(MJournalN)M

Edit2 : Trouvé ce PDF décrivant un tel schéma de "marche" qui est garanti de se terminer, et passe en revue les approches les plus naïves.

Une autre alternative à l'utilisation d'arbres quadratiques consiste à utiliser une hiérarchie de triangulation. Voir Olivier Devillers. Amélioration de la triangulation Delaunay randomisée incrémentielle. Dans Proc. 14th Annu. ACM Sympos. Comput. Geom., Pages 106-115, 1998. Il fonctionne mieux pour les triangulations de Delaunay, mais peut également fonctionner pour les non-Delaunay.

Fondamentalement, quoi que vous fassiez pour accélérer la localisation des points, il faudra construire une structure de données auxiliaire. Dans le cas de quadtrees ou d'une autre subdivision spatiale, vous devez créer l'arborescence de subdivision. Dans le cas du saut de bord, vous devez construire la structure topologique adjacente au triangle. La hiérarchie de triangulation nécessite également de construire un arbre de triangulations plus grossières.

Victor Liu
la source
Victor - connaissez-vous un code open source qui implémente l'approche de saut de bord? IT look slike coudl être une très bonne solution pour mon cas. (modèle de suivi des particules entraîné par les champs actuels dans un maillage de traingualr) -Merci
Chris Barker
J'ai un code pour cela et je peux vous l'envoyer; c'est en C / C ++. Je n'ai pas encore eu le temps de le nettoyer et de le publier sur Github. J'ai dû écrire cela au moins deux fois dans ma vie, une fois avec une structure de données halfedge, encore une fois avec un quadedge, mais il peut facilement être utilisé lorsque ceux-ci ne sont pas disponibles et que vous devez créer vous-même une structure topologique. Regardez sur ma page de profil pour mon site Web, où vous pouvez trouver des informations de contact. Nous pouvons en discuter davantage hors ligne.
Victor Liu
Je suis sur le point de terminer l'implémentation de cela dans matlab en utilisant l'ordre des courbes de Hilbert et la marche aléatoire du triangle. C'est du code de recherche: non optimisé, non documenté, etc., mais toujours assez rapide - je peux vous donner le code si vous êtes intéressé.
Nick Alger
2
A propos: "" "le saut de bord doit être O (logN)" "" Je ne vois pas cela. Par exemple, dans le cas pathologique d'une grande bande triangulaire longue (comme un canal étroit uniquement sur un triangle large), dans le pire des cas, vous devrez sauter d'un triangle au suivant jusqu'à la fin. Dans le cas moyen, à mi-chemin. Donc, si vous doublez le nombre de triangles, ce serait O (N) Dans le cas plus normal d'un arrangement carré de triangles, je m'attendrais à O (sqrt (N)). Ou est-ce que je manque quelque chose? -Chris
Chris Barker
@Chris - Bienvenue sur scicomp! Dans le cadre de l'entretien ménager de scicomp, j'ai migré vos réponses et la conversation qui a suivi pour commenter la réponse de Victor. Nous nous réjouissons de votre participation sur le site.
Aron Ahmadia
8

Je ne suis pas convaincu que votre solution soit correcte. Considérez la situation où vous avez ces nœuds:

  • R: (-3, 1)
  • B: (0, 2)
  • C: (3, 1)
  • D: (0, -5)

Il existe des triangles ABC et ACD. Maintenant, B est le point le plus proche de l'origine, mais l'origine est dans le triangle ACD, qui ne contient pas B.

Pour le reste, cette réponse fournit une solution qui peut dans des cas dégénérés être aussi lente que O(NM), mais sera généralement plus rapide - en particulier pour les triangulations de Delaunay et les triangulations proches de cela dans un certain sens.

J'envisagerais la possibilité de construire un quadtree contenant les triangles eux-mêmes. C'est-à-dire que vous avez un arbre quaternaire qui stocke dans chaque nœud (ce qui correspond à une boîte englobante):

  • Les coordonnées auxquelles la boîte est divisée, ou bien les boîtes englobantes des quatre sous-arbres;
  • Pointeurs vers les sous-arbres;
  • L'ensemble des triangles qui tombent complètement dans la zone de délimitation de ce rectangle, mais pas complètement dans l'un des quatre sous-arbres. En d'autres termes, les triangles qui coupent l'un des deux segments de ligne de subdivision du quadtree.

Lorsqu'on lui donne un point P, parcourez tous les nœuds sur le chemin de la racine du quadtree à la plus petite boîte contenant P. Vous devrez examiner tous les triangles que vous rencontrez dans ces nœuds. Pour une triangulation «bien conduite», il ne devrait y avoir que quelque chose commen triangles qui doivent être examinés au niveau d'un nœud qui contient n triangles dans son sous-arbre, et la profondeur doit être limitée par Journaln. Pour une triangulation «mal comportée», vous pouvez obtenir leO(NM) travail.

Erik P.
la source
Hmm tu as raison. En revanche, si la triangulation était Delaunay, je pense que le plus proche voisin fonctionnerait. C'est trop restrictif pour ce que j'essaie de faire, mais dans le cas Delaunay, considérons le diagramme de Voronoi double - les cellules Voronoi sont l'ensemble des points les plus proches d'un nœud, et les bords des triangles delaunay rencontrent tous les bords du Voronoi cellules à angle droit, donc tout point doit être dans un triangle connecté à son nœud le plus proche. Je me demande si c'est ainsi que fonctionne la fonction pointLocation de matlab sous le capot ..?
Nick Alger
1

CGAL gère les triangulations et l'emplacement des points (et bien plus encore). En particulier, le module de triangulation 2D peut faire ce que vous voulez.

Ronaldo Carpio
la source