Maillage de surfaces de genre zéro sans bord

Sujet proposé par Steve Oudot

Steve.Oudot[at]inria.fr Difficulté : difficile (***).
 

Figure 1: Maillage d'une surface connexe de genre zéro sans bord.


1  Introduction

Le maillage de surfaces est très utilisé en informatique graphique et en CAD pour représenter et manipuler des objets tridimensionnels. L'intérêt des maillages, comparés à d'autres représentations discrètes, est que la plupart des algorithmes de rendu sont câblés dans les puces des cartes graphiques récentes, ce qui permet de manipuler les objets en temps réel. La figure 1 montre un exemple de maillage d'une surface lisse.

L'objectif de ce projet est de mailler des surfaces compactes sans bord. Considéré en toute généralité, ce problème est un sujet de recherche à part entière et n'a trouvé de solutions certifiées que très rcemment [1, 2, 3]. Pour simplifier les choses, nous imposons quelques restrictions sur la surface à mailler, qui sera appelée S par la suite : elle doit être connexe et de genre zéro (i.e. sans anse).

En quoi cela simplifie-t-il les choses ? Du moment que S est de genre zéro sans bord, elle est homéomorphe à la sphère unité S2 et l'on peut (du moins en théorie) trouver un paramétrage sphérique bicontinu de S, i.e. une application Φ:  S2→S bijective, continue et d'inverse continue. Il suffit donc de mailler S2 puis d'envoyer tous les sommets du maillage sur S par Φ, tout en gardant la même connectivité, pour obtenir un maillage de S. Il faut toutefois veiller à mailler S2 assez finement pour s'assurer que le maillage de S ne s'auto-intersecte pas – voir figure 2. Ainsi, la densité locale des sommets du maillage de S2 doit dépendre de S, et donc de Φ. Cette question sera abordée aux parties 3 et 4.


Figure 2: Un cas où le maillage de S (en gras) s'auto-intersecte.


Pour mailler S2, on peut se contenter de mailler le disque unité D dans le plan, puis de remonter les sommets du maillage sur la demi-sphère supérieure par une projection Π. Le choix de Π sera discuté en partie 3. Le maillage de S2 est obtenu en recollant le maillage de la demi-sphère supérieure avec son symétrique par rapport au plan horizontal. Il faut pour cela que les bords des deux maillages concordent dans 3, ce qui peut être assuré en choisissant une projection Π qui laisse le bord du maillage de D dans le plan horizontal. Notre approche consistera à faire en sorte que les sommets du bord du maillage de D se trouvent sur le bord C de D, puis à choisir une fonction Π qui laisse C invariant.

Au final, le travail à efectuer se divise en trois étapes :
  1. construire un maillage triangulaire du disque unité D dans le plan, dont la taille des triangles dépend d'une fonction ϱ:  D→ + donnée en argument (section 2);
  2. relever les sommets du maillage sur la sphère unité S2 (section 3);
  3. projeter les sommets du maillage sur la surface S avec la fonction Φ (section 4).
Les étapes ne sont pas indépendantes : la première est nécessaire pour effectuer la deuxième, qui elle-même est nécessaire pour effectuer la troisième. Chacune contient son lot de difficultés. À noter toutefois que la première consiste essentiellement à implanter un algorithme décrit dans l'énoncé (c'est le minimum à faire), tandis que la deuxième et la troisième sont plus théoriques. Si vous rencontrez des difficultés avec le sujet, n'hésitez pas à me contacter.

2  Maillage de D

2.1  La théorie

Notre objectif ici est de construire un maillage triangulaire de D tel que : Soit E l'ensemble des sommets du maillage. Puisque D est convexe, le bord du maillage est formé des arêtes de Conv(E), l'enveloppe convexe de E. Pour satisfaire le premier point, il suffit donc de faire en sorte que tous les points de E situés sur Conv(E) appartiennent à C.

Pour construire le maillage, nous allons utiliser la triangulation de Delaunay de E.

Définition 1   Une triangulation d'un ensemble fini E de points du plan est un ensemble maximal de triangles reliant les points de E, tel que deux triangles ne peuvent s'intersecter que le long d'un sommet commun ou d'une arête commune.


Définition 2   La triangulation de Delaunay d'un ensemble fini E de points du plan est une triangulation de E telle que tout triangle a un cercle circonscrit vide de points de E. Elle est notée Del(E).


Figure 3: Pour les définitions 1 et 2.


À noter que la triangulation de Delaunay de E existe toujours, et qu'elle est unique si les points de E sont en position générale, i.e. si trois points quelconques ne sont jamais alignés sur le bord de l'enveloppe convexe de E, et si quatre points quelconques ne sont jamais sur un même cercle vide d'autres points de E. Cette hypothèse de généricité sera supposée vérifiée tout au long du projet.

La figure 3 illustre les définitions 1 et 2. L'image de gauche ne représente pas une triangulation car les deux triangles en gras ont des intérieurs qui s'intersectent. L'image centrale montre une triangulation qui n'est pas de Delaunay car le cercle circonscrit au triangle en gras contient un sommet. Enfin, l'image de droite représente une triangulation de Delaunay.

La triangulation de Delaunay est intéressante à deux titres pour le maillage : On donne ci-dessous une variante de l'algorithme de Ruppert [4], basée sur la propriété du cercle vide, qui construit un maillage du disque D :
  1. on choisit 3 points de C qu'on insère dans un ensemble E, dont on calcule la triangulation de Delaunay, Del(E).
  2. tant qu'il reste un triangle t de Del(E) dont le cercle circonscrit (c,r) est tel que r>ϱ(c) :
  3. en sortie de l'algorithme, on retourne Del(E).


Question 1   Montrez que si les points choisis en 1. forment un triangle dont les trois angles sont plus petits que π/2, alors les points de E situés sur Conv(E) appartiennent tous à C, et ce tout au long de l'algorithme.


Question 2   Donnez un contre-exemple dans le cas où les trois points de départ forment un triangle avec au moins un angle plus grand que π/2. Dans ce cas, l'algorithme peut-il se bloquer ? Peut-il boucler indéfiniment ?


En pratique, on prendra les points (1,0), (−1,0) et (0,1) pour initialiser l'algorithme. On déduit immédiatement de la question 1 que le bord du maillage de sortie est conforme à C. Par ailleurs, l'algorithme fait en sorte que les rayons des cercles circonscrits des triangles de Del(E) soient petits, donc le maillage de sortie vérifie les deux points énoncés au début de la section 2.1. Par ailleurs, on admettra que l'algorithme termine toujours. Une preuve de terminaison se trouve dans [4].

2.2  Implantation

L'objectif principal de la section 2 est d'implanter l'algorithme décrit plus haut. Le design de l'implantation est libre, tout comme le langage utilisé.

Pour les utilisateurs de JAVA, nous fournissons une implantation de la triangulation de Delaunay, effectuée par l'un de vos camarades il y a quelques années. Une archive contenant le code source peut être téléchargée à l'adresse suivante : ftp://ftp-sop.inria.fr/geometrica/soudot/cours/IF_Delaunay.tgz. Cette archive contient entre autres un fichier lisez_moi.txt qui donne des instructions pour compiler et exécuter le programme. L'objectif est bien entendu de partir du code fourni pour implanter votre propre code.

Pour les utilisateurs de C++, nous recommandons la bibliothèque CGAL, disponible sur le site http://www.cgal.org et installée sur les machines de l'École. Cette bibliothèque fournit une implantation robuste et efficace de la triangulation de Delaunay.

 

Le stockage des arêtes de Conv(E) ne nécessite pas de structure de données particulière, puisque ces arêtes sont dans la triangulation de Delaunay. En revanche, casser une arête de Conv(E) nécessite d'intersecter sa médiatrice avec C.

La fonction de densité ϱ:  D→+ sera considérée égale à une constante ϱ0 dans un premier temps. Dans les sections 3 et 4, on essaiera de définir une fonction de densité adaptée à S2 puis à S.

Pour la visualisation des maillages 3D, nous recommandons le programme Geomview, installé sur les machines de l'École. Pour lancer le programme, tapez geomview dans un terminal. Geomview prend en entrée un fichier au format OFF, structuré comme suit :

 

 
OFF
nSommets nFaces 0
x0  y0  z0 ## Coordonnées du premier point
x1  y1  z1 ## Coordonnées du deuxième point
⋯
xn−1  yn−1  zn−1 ## Coordonnées du dernier point
 
3  p1  p2  p3 ## 3 pour face triangulaire, suivi des indices (entre 0 et n−1)
3  q1  q2  q3 ## des trois sommets du triangle
⋯


Question 3   Faites tourner votre programme pour diverses valeurs de la constante ϱ0. Comment évolue la densité des sommets du maillage ? Dressez un tableau de l'évolution du nombre de sommets en fonction de la valeur de ϱ0. Quelle loi empirique pouvez-vous extraire de vos données ?


3  Projection de D sur S2

Maintenant que nous disposons d'un maillage M de D conforme à C, il nous faut remonter les sommets de M sur S2 par une certaine projection Π, qui laisse C invariant, afin d'obtenir un maillage de la demi-sphère supérieure qui soit lui aussi conforme à C . Le maillage de S2 est obtenu en recollant le maillage de la demi-sphère supérieure avec son symétrique par rapport au plan horizontal1.

À ce stade, plusieurs fonctions de projection Π sont possibles. La plus simple est Π0: (x , y)↦ ( x , y , √1−x2−y2).

Question 1   Implantez Π0 et utilisez-la pour construire le maillage de S2 à partir de M. Testez votre programme : quel est le problème ? Essayez de modifier la fonction de densité ϱ, par exemple en la faisant dépendre de la distance au centre O de D. Le problème se résout-il ? Expliquez.


On se propose d'utiliser à la place de Π0 une projection stéréographique, définie par :
Πs: (x , y)↦ ⎛
⎜
⎜
⎝
2x
1+x2+y2
 , 
2y
1+x2+y2
 , 
1−x2−y2
1+x2+y2
⎞
⎟
⎟
⎠


Question 2   Remplacez la projection Π0 par Πs dans votre programme et essayez de trouver expérimentalement une fonction ϱ qui résolve le problème évoqué à la question 1.


Remarque Avant de vous lancer dans le codage, vérifiez quand-même par le calcul que Πs est bien à valeurs dans S2...

 

Une fois les sommets du maillage remontés sur S2, certains triangles ont plus grossi que d'autres. C'est un effet de ce qu'on appelle l'étirement. Formellement, étant donné une fonction f: 2→3 et un point p∈2, l'étirement associé à f en p, Efp, est la forme quadratique définie par:
∀ v∈2,  Efp(v )  =  ∂v f(p)2 = 
 
lim
h→ 0+
dist ⎛
⎝		 f(p),f(p+hv) ⎞
⎠		 2
h
    (1)


où ∂v f(p) est la dérivée partielle de f en p le long du vecteur v (∂v f(p) est un vecteur de 3). Efp(v) est l'étirement associé à f en p le long de v.

Puisque Efp est une forme quadratique, sa matrice est diagonalisable dans une base orthonormée. En d'autres termes, il existe deux vecteurs unitaires v1(p) et v2(p), appelés directions principales2, tels que pour n'importe quel vecteur non nul v∈2, Efp(v)/v2 est une combinaison linéaire convexe de Efp(v1) et de Efp(v2).

Question 3   Montrez que l'étirement associé à Πs en tout p∈ D est isotrope, c'est-à-dire qu'il existe un réel λΠs(p) tel que ∀ v∈2, EΠsp(v) = λΠs(p) v2. En déduire que λΠs = (∂Πs/∂ x)2 = (∂Πs/∂ y)2.


Connaître l'étirement permet de fixer la densité ϱ du maillage M de sorte que les cercles circonscrits aux triangles du maillage de S2 aient tous à peu près le même rayon. Dans le cas de Πs, l'étirement est isotrope, i.e. il ne dépend pas de la direction choisie. Nous pouvons donc fixer ϱ(p) = 1/√λΠs(p), où λΠs(p) est défini comme dans la question 3. Ceci nous permet de générer un maillage de la sphère S2, de densité uniforme.

Question 4   Ajoutez à votre code une fonction qui calcule λΠs, et utilisez-la pour construire le maillage de S2. Vérifiez expérimentalement que le maillage que vous obtenez a une densité à peu près uniforme. Pour cela, vous pourrez par exemple tracer l'histogramme de la répartition des rayons des triangles de votre maillage.


4  Projection de S2 sur S

Puisque S est compacte, connexe et sans bord, elle borne un domaine compact Ω de 3. Pour simplifier les choses, on suppose que Ω est étoilé par rapport à O, c'est-à-dire :
∀ p∈ Ω, [Op]⊂Ω


Dès lors, on peut supposer que Φ est une fonction radiale, à savoir :
∀ p∈ S2, Φ(p)=r(p)  p


où r est une fonction de S2 dans +. Voici quelques exemples de fonctions r(p), avec p=(x,y,z) (x2+y2+z2=1 car p∈ S2):

⎧
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎨
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎩
r1(x,y,z) = 2 (sphère de rayon 2)
 
r2(x,y,z) = 0.1 + z2 (cacahuète)
 
r3(x,y,z) = 0.1 + x2 (cacahuète orientée le long de l'axe des abcisses)
 
r4(x,y,z) = 1 + 0.3 ⎛
⎝		 1+z ⎞
⎠		 ⎛
⎝		 1+z ⎞
⎠
(oeuf)
 
r5(x,y,z) =  ⎧
⎪
⎨
⎪
⎩
1
√
x2 + y2
− z
si z≥0
1 si z<0
(toupie)
 
r6(x,y,z) = 1 + (1 − |z|)2 (soucoupe volante)


Vous pouvez implanter ces fonctions dans votre programme pour voir les résultats produits. Un exemple est donné dans la figure 4.


Figure 4: Résultat de l'algorithme pour la fonction radiale r2 : à gauche, le maillage de D ; au centre, le maillage de S2 ; à droite, le maillage de S, de densité uniforme.


Pour générer un maillage uniforme de S, il faut compenser l'étirement lié à Φ ∘ Πs dans le calcul de la fonction de densité ϱ. D'après la formule (1), l'étirement en un point dépend généralement de la direction dans laquelle on regarde. Le cas de la projection stéréographique était simple puisque cette dernière a un étirement isotrope en tout point de D. En revanche, le cas général d'un étirement anisotrope est problématique puisqu'on ne sait pas quelle direction privilégier. Pour s'en sortir, on ne considérera que la composante isotrope de l'étirement, qui est un scalaire. Cette approximation donne d'assez bons résultats en pratique, comme le montre la figure 4.

Plus précisément, puisque Φ est radiale, on peut se placer en coordonnées sphériques (θ, ϕ) pour faire nos calculs. Comme r est à valeurs dans , on peut calculer son gradient ∇ r en tout point p=(θ, ϕ) de S2 :
∇ r = 
1
r
∂ r
∂ θ
 θ + 
1
rsinθ
∂ r
∂ ϕ
 ϕ


où θ et ϕ sont les vecteurs unitaires associés aux coordonnées curvilignes θ et ϕ. Le gradient correspond, dans le plan tangent à S2 en p, à la direction de plus forte variation de r. Or, comme Φ est radiale, son étirement est directement lié à la variation de r. Il résulte que ∇ r/∥∇ r∥ est l'une des deux directions principales de l'étirement de Φ, l'autre direction principale étant donnée par n'importe quel vecteur du plan tangent orthogonal à ∇ r, comme par exemple
∇ r⊥ = −
1
rsinθ
∂ r
∂ ϕ
 θ + 
1
r
∂ r
∂ θ
 ϕ


Les étirements EΦp(∇ r / ∥∇ r∥) et EΦp(∇ r⊥/∥∇ r⊥∥) le long des vecteurs ∇ r/∥∇ r∥ et ∇ r⊥/∥∇ r⊥∥ ne sont autres que les valeurs propres de la forme quadratique EΦp. La composante isotrope λi(p) vaut alors λi(p) = √ EΦp(∇ r/∥∇ r∥)  EΦp(∇ r⊥/∥∇ r⊥∥). Nous pouvons donc fixer ϱ(p) = 1/√λΠs(p) λi(p), où λΠs(p) est défini comme dans la question 3, pour générer un maillage à peu près uniforme de S.

Question 1   Ajoutez à votre code des fonctions qui calculent les composantes isotropes λi des étirements associés aux fonctions r définies plus haut. Vérifiez expérimentalement que les maillages produits par votre programme sont de densité à peu près uniforme.


Remarque Les fonctions r4 (oeuf) et r5 (toupie) ne sont pas symétriques en z, ce qui signifie que l'une des deux moitiés (inférieure ou supérieure) de S ne sera pas maillée uniformément. Pour améliorer votre code, vous pouvez mailler deux fois D, une fois pour chaque moitié de S, puis recoller les deux maillages obtenus. Attention toutefois à vérifier que les deux maillages ont le même bord.

5  Améliorations

Les maillages que nous avons produits jusqu'ici ont un défaut : du fait que l'on recolle les maillages des deux demi-sphères le long de l'équateur, ce dernier apparaît clairement dans le maillage de S2, ainsi que dans celui de S puisque Φ est radiale – voir l'exemple de la figure 4. Une méthode pour pallier ce défaut consiste à perturber aléatoirement les sommets du maillage de D se trouvant sur C. La perturbation doit s'effectuer vers l'intérieur de D, pour que les sommets puissent encore être projetés sur S2 par Πs. Par ailleurs, l'amplitude de la perturbation doit bien évidemment dépendre de la fonction de densité ϱ choisie.

Question 1   Ajoutez à votre code une méthode qui perturbe les sommets du maillage de D se trouvant sur C. Observez le résultat sur les maillages de S2 et de S, pour les fonctions r définies plus haut. Perturber les sommets a-t-il toujours un intérêt ?


On veut maintenant générer des maillages dont la densité soit adaptée à la courbure locale de S. Ceci permettra de réduire l'erreur d'approximation. Puisque notre fonction de densité ϱ nous garantit une densité uniforme, il suffit de la pondérer par la courbure de S, elle-même liée à la différentielle seconde de r.

Question 2   Modifiez votre fonction de densité ϱ pour qu'elle prenne en compte3 la courbure de S. Vérifiez visuellement sur plusieurs exemples que les maillages obtenus ont des densités adaptées à la courbure de S.


6  Extensions

Plusieurs extensions du travail effectué dans le cadre de ce projet sont envisageables :

Références

[1]
Jean-Daniel Boissonnat and Steve Oudot. An effective condition for sampling surfaces with guarantees. In Proc. 9th Symp. on Solid Modeling, pages 101–112, 2004.

[2]
Jean-Daniel Boissonnat, David Cohen Steiner, and Gert Vegter. Isotopic implicit surface meshing. In Proceedings of the thirty-sixth annual ACM symposium on Theory of computing, pages 301 – 309, 2004.

[3]
Siu-Wing Cheng, Tamal K. Dey, Edgar A. Ramos, and Tathagata Ray. Sampling and meshing a surface with guaranteed topology and geometry. In Proc. 20th Annu. ACM Sympos. Comput. Geom., pages 280 – 289, 2004.

[4]
J. Ruppert. A Delaunay refinement algorithm for quality 2-dimensional mesh generation. J. Algorithms, 18:548–585, 1995.
1
Attention à ne pas dédoubler les sommets de M se trouvant sur C !
2
ce sont les vecteurs propres de la matrice de Efp
3
Il y a mille manières de le faire, sans forcément calculer la différentielle seconde de r. Choisissez la méthode qui vous convient le mieux, en justifiant votre choix.
Ce document a été traduit de LATEX par HEVEA