Calcul de squelette topologique et reconstruction de maillage-Resultats

De Ensiwiki
Révision de 17 juin 2010 à 09:09 par Hetroyf (discussion | contributions)

(diff) ← Version précédente | Voir la version courante (diff) | Version suivante → (diff)
Aller à : navigation, rechercher


Calcul de squelette topologique et reconstruction de maillage
Projet Projet de spécialité 2A
Sous-projet Image
Étudiants Emeline BANTEGNIE (MMIS-IRVM), Marion NOT (SLE), Laura PEYTHIEUX (MMIS-IRVM), Edouard SIEGEL (MMIS-IRVM)
Promo 2011
Tuteur Franck HETROY franck.hetroy@grenoble-inp.fr

Etudiants

Bantegnie, Emeline (MMIS IRVM)

Not, Marion (SLE)

Peythieux, Laura (MMIS IRVM)

Siegel, Edouard (MMIS IRVM)

Introduction

Le scanner acquis par l'Ensimag

Auparavant réservés à des spécialistes du domaine de la modélisation en 3D, les scanners 3D se démocratisent aujourd'hui de plus en plus. On assiste à l'apparition sur le marché de différents modèles allant du scanner pour bâtiments au scanner de bureau. De fait, beaucoup de non spécialistes ont récemment découvert cette technologies et se sont rendus compte des avantages qu'elle pouvait leur apporter. Modéliser de façon fiable le réel sur un ordinateur facilite calculs, traitements, simulations etc.

Il existe cependant des limites à leur utilisation. La sortie d'un scanner est obtenue sous forme d'un maillage, c'est à dire un ensemble de faces triangulaires adjacentes formant une surface fermée. Mais, de par le principe même du scanner (lancement et réception de rayons lasers), les maillages obtenus sont incomplets, les rayons ne pouvant attendre toutes les zones de l'objet à scanner (occlusions, formes complexes…). De nombreuses recherches se sont donc axées sur des méthodes de reconstruction de ces maillages incomplets. Ce projet a pour but d'étudier une de ces méthodes, issue de l'article de Tagliasacchi et al.[1], l'implémenter et en souligner les avantages et les failles.

Eléments de prérequis

Bibliothèque OpenMesh - Tutoriel

De nombreuses bibliothèques existent pour créer et manipuler des maillages 3D. Nous avons choisis d'utiliser OpenMesh car, bien que doté d'une documentation peu claire, elle est assez facile de compréhension et d'implémentation dans un code en C++. Elle permet de plus de traiter les fichiers .off et .obj obtenus en sortie de notre scanner.

Nous avons cependant eu des difficultés à installer les bibliothèques nécessaires à son utilisation ainsi que celles liées à l'interface graphique fournie. Pour éviter ce type de problème, qui nous a beaucoup coûté en temps, aux générations suivantes qui travailleront sur OpenMesh, nous fournissons un Tutoriel expliquant comment installer ces bibliothèques et configurer l'IDE (ici eclipse) pour les utiliser. Nous proposons également un lien vers une page de la documentation d'OpenMesh qui fournit un ensemble de fonctions très pratiques, dans l'optique d'une prise en main rapide de la bibliothèque.

Maillage

Un maillage tel que ceux que nous avons manipulés est un ensemble de sommets décrits par leurs coordonnées 3D associé à un ensemble de faces reliant entre elles ces sommets. Nous avons exclusivement utilisé des maillages à faces triangulaires dans ce projet, leur manipulation étant géométriquement moins complexe (en particulier, une face triangulaire est toujours plane).

La manipulation des maillages est facilitée par la structure proposée par OpenMesh. En particulier, les points sont reliés entre eux par des arêtes constituées de deux demi-arêtes orientées (edges et halfedges), ce qui facilite l'itération sur les éléments. La librairie fournit des itérateurs sur tous les types d'objets et des sous-itérateurs adaptés (itération sur tous les sommets voisins d'un sommet par exemple).

Un maillage peut être décrit par un fichier .off (Object File Format). Nous avons utilisé ce format tout au long de notre projet.

Éléments de maillages

Le squelette topologique

Skeleton.png

Le squelette topologique est, par définition, le lieu des centres des boules de rayon maximal incluses dans l'objet.

Il est assimilable à l'axe médian d'un objet.

Ici, comme nous étudions des objets à géométrie cylindrique, le squelette sera constitué de branches (axes de symétrie des cylindres) qu'il s'agira de relier.

Algorithme de Tagliasacchi

Cet algorithme vise à extraire le squelette topologique d'un maillage incomplet, avec pour objectif à terme de reconstituer le maillage complet.

Il se base sur l'hypothèse que les objets traités peuvent être découpés en sections globalement cylindriques, appelées "branches", reliées entre elles par des sections non cylindriques appelées "jointures". Pour chaque branche, la partie du squelette correspondante sera l'axe de rotation de la section cylindrique. On appliquera ensuite un traitement spécial aux jointures afin de reconstituer un squelette unidimensionnel pertinent.

L'algorithme de Tagliasacchi se veut très résistant au manque de données, même dans le cas où celui-ci est très important. En contrepartie, il ne peut s'appliquer à des modèles aux formes sphériques ou planes, ni à des modèles trop détaillés ou complexes.

Travail Effectué

Calcul du squelette

Le squelette est obtenu à partir des points du ROSA (ROtationnal Symmetry Axis), qui correspondent localement à l'axe de symétrie rotationnelle d'une section cylindrique du modèle traité. On peut obtenir un point du ROSA pour chacun des points du modèle : En pratique, on constate qu'en ne traitant qu'un point sur 10 ou 20, on obtient suffisamment de points ROSA pour reconstituer un squelette pertinent.

Pour chaque point traité :

  • On choisit un plan initial passant par ce point
  • On calcule le voisinage du point relativement à ce plan
  • Grâce au voisinage, on met à jour le plan
  • On itère les deux actions précédentes jusqu'à obtenir un plan de coupe satisfaisant, c'est-à dire le plus perpendiculaire possible à la direction de la section cylindrique traitée.
  • On extrait le point du ROSA du voisinage

Choix du plan initial

Plan initial-Tore

L'article de Tagliasacchi et al. mentionne que l'on peut choisir comme plan initial tout plan contenant la normale au point traité. Cependant, lors de l'implémentation, on constate que l'on ne converge vers le plan de coupe souhaité que si le plan initial présente une orientation "favorable", c'est-à-dire peu éloignée de la direction souhaitée. Des orientations trop différentes conduisent soit à un crash de l'application lorsqu'elle essaie de déterminer l'orientation du plan suivant, soit la convergence vers un plan non désiré.

Afin de pallier à ce problème, nous adaptons le choix du plan initial à la géométrie du modèle. On définit une "direction interdite" qui ne peut en aucun cas être la direction de la normale du plan (dans le cas du tore par exemple, la direction interdite est la verticale), puis on effectue le produit vectoriel entre cette direction et la normale du point (flèches rouges sur le schéma). La direction obtenue (flèche noire) définit le plan initial.

Cette méthode entraîne une limitation : Dans le cas de modèles trop complexes (avec des sections cylindriques dans toutes les directions de l'espace), le calcul du plan initial devient impossible.

Calcul du voisinage

Le calcul du voisinage est basé sur la distance de Mahalanobis. Nous nous référons pour cela au travail de Lehtinen et al.[2].

La distance de Mahalanobis entre deux points permet de prendre en compte à la fois la distance euclidienne entre eux, mais y ajoute la comparaison entre normales. En effet, si l'on utilise la seule distance géométrique entre deux points, si deux portions cylindriques sont poches, le voisinage d'un point comprendra des points des deux branches et non seulement de celle à laquelle il appartient. Ainsi deux points proches géométriquement mais de normales opposées seront loin au sens de Mahalanobis, ce qui permet alors de distinguer deux branches cylindriques proches. La distance de Mahalanobis entre deux points P(x,n) et Pi(xi,ni), où x représente la position du point et n sa normale est:

 \displaystyle{dist_{j}(x) = \| x - x_{j} + \mathsf{F}_{squash} <x - x_{j},x - n_{j}> n_{j}\|}


Nous avons choisi de suivre Lehtinen et de choisir un Fsquash=2, les résultats étant satisfaisants.


Un voisinage d'un point est définit par un plan dit plan de coupe («Cutting plane»). De manière pratique, le calcul du voisinage d'un point se fait de la façon suivante:

  • On élimine tous les points situés à une distance supérieure à une constante delta du plan de coupe. Delta, en accord avec l'article, est égal à 2,5% de la diagonale de la bounding box du modèle.
  • Dans le nuage de point restant, on ne conserve que les points situés à une distance de Mahalanobis inférieure à un \epsilon_{Mah} constant.


Concernant le \epsilon_{Mah}, l'article ne précisait rien sur sa détermination. Nous nous sommes aperçu que sa valeur dépendait fortement du modèle utilisé, et même de la branche du modèle traitée. Pour cela, nous avons décidé d'adapter le \epsilon_{Mah} au lieu de le garder constant pour un même modèle.

On fixe un nombre de sommets maximal et minimal nécessaires et on adapte le \epsilon_{Mah} pour obtenir un voisinage pertinent (souvent plus de 50 et moins de 400 sommets).

Calcul du plan de coupe optimal

Ayant déterminé un premier voisinage, il va nous permettre de mettre à jour notre plan de coupe pour se rapprocher d'un plan perpendiculaire à l'axe du cylindre.

Pour cela, sur un voisinage donné, on souhaite trouver la normale du plan qui minimise la variance de l'angle entre les normales des points et cette normale. Ceci revient donc, à l'itération t à trouver la normale v_{i} correspondant au voisinage N_{i} qui vérifie:

 v_{i}^{t+1} = \underset{v\in\mathbb{R}^{3}, \|v\| = 1}{argmin} var\{<v,n(p_{j})> : p_{j} \in N_{i}^{(t)}\}, t \geqslant 0 


Ceci revient à minimiser la forme quadratique Quad.png où v est de norme 1 et

Mat.png

Comme proposé dans l'article, on a résolu ce problème de minimisation en utilisant la SVD (Singular Value Decomposition).

On s'est muni pour cela de la librairie AlgLib[3].

La résolution par SVD revient à factoriser la matrice sous la forme:

 \mathsf{M} = \mathsf{U}\mathsf{S}\mathsf{V}^{\mathsf{T}}

Le vecteur minimisant notre problème étant la dernière colonne de la matrice V ("right singular value").

Les détails de la résolution mathématique sont dans les articles référencés en [4].

Calcul du point de ROSA

Point du ROSA

Une fois que l'on a obtenu un voisinage pertinent, c'est-à-dire une portion globalement cylindrique, on recherche la position du point du ROSA correspondant. Dans le cas idéal, c'est-à-dire une portion réellement cylindrique, ce point serait l'unique point d'intersection de toutes les normales aux points situés dans le plan de coupe.

Dans notre cas, on cherche le point qui minimise la somme de ses projections sur les normales aux points du voisinage par la formule :

 r_{i}^{*} = \underset{r\in\mathbb{R}^{3}}{argmin} \sum_{i=1}^{N}{\|(r - p_{j}) * n(p_{j})\|^2}

Résultats

Maillage complet simple
Maillage complet complexe
Maillage incomplet complexe

Gestion des jointures

Ceci se fait en 3 étapes successives successives.

Lissage laplacien

Dans un premier temps, on souhaite augmenter la cohérence spatiale des points de rosa au niveau des jointures. L'absence de symétrie rotationnelle a crée des nuages de points sur ces jointures. On va appliquer un lissage laplacien à chacun des points du squelette. Le lissage laplacien est un opérateur assez simple. Chaque point va se voir remplacer par le barycentre de tout ses points voisins. On a choisi de ne pas prendre en compte la surface des faces pour pondérer ce barycentre. Cette opération est itérée un certain nombre de fois, déterminé expérimentalement selon le modèle. 5 à 10 itérations ont donné de bons résultats lors de nos tests.

Etape du lissage Laplacien

Analyse des composantes principales (PCA) et Méthode des moindres carrés (Moving Least Squares- MLS)

Après le lissage laplacien, l'article nous disait, pour affiner encore le squelette au niveau des branches, d'appliquer des MLS (moving least squares) en dimension 1 à certains points.

Etape d'application des MLS


Pour savoir sur quels points appliquer ce principe, on a eu besoin d'effectuer une analyse en composantes principales sur le voisinage (calculé précédemment) du point considéré.

La première étape a été de calculer la matrice de covariance des points du voisinage selon la formule suivante, tirée de l'article de Mederos et al.[5]

 C = \left( p_{i} - \overline{p} \right) ^{T} \left( p_{i} - \overline{p} \right) 

Les p_{i} sont les points du voisinage, et \overline{p} = \displaystyle{\frac{1}{N}\sum_{i = 1}^{N}{p_{i}}}.


Une fois que l'on avait cette matrice, on a pu calculer la fonction \psi

 \displaystyle{\psi (r_{i}) = \frac{\lambda_{i}^{1}}{\lambda_{i}^{1} + \lambda_{i}^{2} + \lambda_{i}^{3}}}

dans laquelle les \lambda_{i} sont les valeurs propres de la matrice de covariance au point r_{i} avec  \lambda_{i}^{1} \leqslant \lambda_{i}^{2} \leqslant \lambda_{i}^{3} .

Cette valeur, si inférieure à un \epsilon_{MLS} (égal à 0,4 dans l'article), indique que le point considéré est trop éloigné de la branche, et que l'on doit donc appliquer les MLS sur ce point.

L'étape suivante a été plus laborieuse, en particulier parce que l'article était très flou concernant cette partie. La formule des Moving Least Squares[6] est donnée par :

 \displaystyle{\sum_{i=1}^{N}{(p(x) - f_{i})^{2} \theta(\|x - x_{i}\|)}}


Le principe est en fait de trouver le p (dans notre cas sous forme d'une fonction polynômiale) qui minimise la formule ci-dessus. Le problème qui s'est alors posé est que les f_{i} sont des fonctions f:\mathbb{R}^{3} \rightarrow \mathbb{R} des x_{i}.

Sachant, qu'en entrée, on avait le voisinage du point x, on connaissait les x_{i}, mais il manquait les valeurs correspondantes des f_{i}.

Comme le but est de rapprocher le point de la branche, on a considéré que les valeurs des f_{i} correspondaient à la distance entre le point x_{i} et sa projection sur la branche. La branche correspond en fait à la direction principale du voisinage, que l'on a donc calculée à partir de la matrice de covariance (il s'agit d'un vecteur propre associé à la plus grande valeur propre de la matrice de covariance)

Ensuite, on a résolu le problème des moindres carrés en ajoutant des poids (les \theta(\|x - x_{i}\|)) puis on a calculé la valeur du p trouvé en x qui nous a permis de déplacer le point à la distance voulue de la direction principale (c'est à dire la branche).



C'est la partie de l'article que l'on a implémentée mais qui ne fonctionne pas comme il faudrait. Le problème vient en fait principalement des voisinages, en effet, dans l'article il était question d'échantillons du maillage sans plus de précisions. En théorie, si les échantillons sont bien choisis, la méthode devrait donner des résultats cohérents, seulement si les échantillons ne sont pas adaptés la direction principale trouvée ne correspond pas du tout à celle de la branche et les résultats obtenus sont totalement inattendus.

Recentrage et extraction du squelette

Après avoir appliqué la méthode des moindres carrés, les branches se retrouvent décentrées, et les nuages de points bien que plus compacts, ne permettent pas encore d'obtenir un squelette unidimensionnel. Dans cette dernière étape, on réapplique un calcul de point de ROSA, avec un traitement différent selon que l'on s'occupe d'une branche ou d'une jointure.

Pour une branche, on applique simplement le calcul du point, pour le recentrer. Pour une jointure par contre, on applique la formule au voisinage global des points de la jointure. C'est à dire à l'union des voisinages de chacun des points de la jointure.

Etape de recentrage

Reste alors à relier le point de la jointure à chacune de ses branches adjacentes. L'implémentation choisie consiste juste à rajouter des points entre celles-ci, selon un pas assez petit pour se rapprocher d'une ligne continue. Le discrétisation du squelette est alors, quasi immédiate.

Reconstruction du maillage

Cette partie n'a pas été implémentée lors de notre projet, du fait d'un manque de temps principalement. En substance, l'objectif est de raffiner le nuage de points à partir du squelette précédemment construit. Pour ça, on souhaite compléter les trous des zones localement cylindriques sous l'hypothèse d'objets scannés lisses. On discrétise alors radialement le voisinage de chaque point du squelette pour résoudre une minimisation de fonction implicite.

Le nuage de points ainsi complété se voit appliquer une reconstruction de poisson, d'après les sources de Kazhdan et al.[7]. Code et exécutable sont disponibles[8]. On obtient alors un modèle par surfaces implicites, auquel on applique un traitement classique de Marching Cubes ou Marching Triangles afin de le mailler.

Bilan

Une partie de l'algorithme est implémentée et fonctionne. Les parties de recentrage et de lissage peuvent cependant être améliorées, en créant un lien plus pertinent entre un point du squelette et ses voisins.

Il reste enfin à implémenter l'étape de complétion du nuage de point et de maillage de la surface implicite.


Ce projet nous a beaucoup apporté à tous. Etant très dense, il faisait appel à beaucoup de principes et méthodes bien connues du domaine mais qui nous étaient inconnues. Cela nous a donc permis de les aborder et de les comprendre. Cela a également développé notre esprit critique vis à vis d'un algorithme, d'en trouver les failles et réfléchir à comment les combler.

Nous déplorons simplement que le temps nous ait manqué pour terminer l'implémentation de l'algorithme global.

Références

Documents additionnels