EP2904587A2

EP2904587A2 - Procédé de traitement de données définissant un élément dans un espace e de dimensions d, programme d'ordinateur associé

Info

Publication number: EP2904587A2
Application number: EP13774392.8A
Authority: EP
Inventors: Bruno Levy
Original assignee: Institut National de Recherche en Informatique et en Automatique INRIA
Current assignee: Institut National de Recherche en Informatique et en Automatique INRIA
Priority date: 2012-10-05
Filing date: 2013-10-03
Publication date: 2015-08-12
Also published as: US10460006B2; CN104937640B; WO2014053606A2; CN104937640A; FR2996669A1; US20150254209A1

Description

Procédé de traitement de données définissant un élément dans un espace E de dimensions d, programme d'ordinateur associé

La présente invention concerne les procédés de traitement de données définissant un élément dans un espace E de dimensions d, comprenant les étapes suivantes :

al on définit un jeu de n points dans l'espace E ;

b/ on calcule l'intersection entre ledit élément et des cellules de Voronoï déterminées chacune à partir d'un point respectivement associé parmi lesdits n points définis.

De telles techniques sont utilisées, entre autres, pour générer un maillage d'un élément tel qu'une surface ou un volume.

Les maillages de surface sont utilisés dans de nombreuses applications techniques, par exemple dans la simulation numérique de l'écoulement pour la conception aéronautique ou pour l'exploration pétrolière, dans la mécanique des structures pour la conception de charpentes, ponts et autres ouvrages d'art, dans la mécanique des déformations pour la simulation de tests d'impact automobiles etc.

Le maillage consiste à découper l'élément considéré en un jeu de cellules de base. Chaque cellule de base et le jeu de cellules de base dans son ensemble doivent satisfaire des conditions de validité prédéfinies, portant par exemple sur la forme géométrique des cellules de base (carré, triangle, tétraèdre, hexaèdre, cube, etc.) ou sur des valeurs d'angle aux sommets (respect de valeurs d'angles minimales/maximales).

Ces conditions prédéfinies dépendent en grande partie de l'application à laquelle le maillage est destiné.

Il existe de nombreux procédés pour générer un maillage d'un domaine.

On connaît certains procédés basés sur l'optimisation d'une fonction « objectif » dépendant des coordonnées aux sommets des cellules. De tels procédés sont parfois désignés "procédés variationnels".

Des procédés utilisent le pavage barycentrique de Voronoï (en anglais "Centroidal

Voronoi Tessellation", ou "CVT") et la triangulation optimale de Delaunay (Optimal Delaunay Triangulation" ou "ODT' en anglais).

Ils permettent de générer efficacement et de manière robuste des maillages isotropes. Ces procédés se révèlent notamment particulièrement efficaces lorsqu'il s'agit de générer un maillage isotrope à partir de cellules de base triangulaires (maillage d'une surface) ou tétraédriques (maillage d'un volume).

Par ailleurs, certaines applications peuvent nécessiter la génération d'un maillage anisotrope, c'est-à-dire qui présente des éléments qui diffèrent les uns des autres en taille et/ou en orientation, dans des zones prédéfinies de l'élément. Une matrice d'anisotropie est définie pour chaque cellule de base spécifiant ces contraintes d'orientation et de taille.

Le document FR 2962582 décrit une technique de génération de maillage anisotrope.

Des techniques sont apparues, selon lesquelles on remplace l'anisotropie par des dimensions supplémentaires de l'espace : un élément anisotrope de dimension dO, par exemple 3D, est représenté par un élément isotrope correspondant de dimension supérieure d, par exemple 6D. On parle alors de plongement de l'élément de dimension dO en un élément correspondant de dimension d.

A titre d'illustration, en référence aux figures 1 et 2, pour générer un maillage anisotrope 2D gouverné par la métrique anisotrope 2D représentée en partie A de la figure 1 , une surface 3D isotrope correspondante est générée (cf. partie A de la figure 2).

Sur la partie A de la figure 1 , chaque cercle déformé représente les points équidistants du centre du cercle déformé en terme de distance anisotrope définie par la métrique anisotrope considérée.

Selon ces techniques, un diagramme de Voronoï isotrope est alors généré pour cette surface 3D, représenté en partie B de la figure 2. Puis un maillage isotrope par triangulation de Delaunay, représenté en partie C de la figure 2, est déduit de ce diagramme de Voronoï.

Les surfaces 3D isotropes représentées sur les parties A, B, C de la figure 2 sont projetées sur le plan (xOy) pour générer les surfaces 2D anisotropes correspondantes des parties A, B, C de la figure 1 .

Certaines de ces techniques sont notamment exposées dans les documents ci- après, nommés références_1 :

- J. F. Nash, The imbedding problem for riemannian manifolds. Annals of Math- ematics, 63:20-63, 1956 ;

- F. Labelle and J.-R. Shewchuk, Anisotropic Voronoi diagrams and guaranteed- quality anisotropic mesh génération, In SCG '03: Proceedings of the nineteenth annual symposium on Computational geometry, pages 191 -200, 2003.

Ces techniques donnent de bons résultats pour la prise en compte de l'anisotropie, mais le plongement dans des dimensions d supérieures donne lieu à un accroissement très important du volume de calculs nécessaire.

En effet, le calcul d'une triangulation de Delaunay a une complexité proportionnelle à d!, où d est la dimension de l'espace de calcul. Généralement, l'espace de plongement est de dimension 6D, voire 10D, ce qui rend prohibitif le volume de calcul. La présente invention vise à proposer une solution pour réduire la charge de calcul dans un espace de dimension d.

A cet effet, suivant un premier aspect, l'invention propose un procédé de traitement de données du type précité caractérisé en ce que dans l'étape b/, on détermine le résultat de l'intersection entre ledit élément et une cellule de Voronoï déterminée à partir du point qui lui associé x,, en mettant en œuvre un traitement itératif selon lequel à un pas d'itération courant k, on sélectionne un point supplémentaire x_jk du jeu de points autre que le point x, et on calcule un résultat actualisé d'intersection égal à l'intersection entre d'une part la moitié de l'espace E délimitée par l'hyperplan médiateur du segment (x,, x_Jk ) et qui contient x,, et d'autre part le résultat d'intersection calculé au pas d'itération précédent,

et selon lequel la sélection d'au moins un point supplémentaire lors des pas d'itération est fonction d'une comparaison entre la distance entre le point associé x, et ledit point supplémentaire et le double de la distance maximale existant entre le point associé x, et un point d'un résultat d'intersection calculé lors du traitement itératif.

L'invention permet de réduire le volume de calculs nécessaires au calcul des intersections des cellules de Voronoï avec l'élément considéré, en utilisant un test qui permet d'éviter des calculs inutiles. Elle permet en outre de paralléliser les traitements pour des cellules de Voronoï distinctes.

Dans des modes de réalisation, le procédé suivant l'invention comporte en outre une ou plusieurs des caractéristiques suivantes :

- dans l'étape b/, les n-1 autres points du jeu classés par ordre de distance croissante par rapport au point x,, sont libellés _y ,- - - x_{Jn 1} , et le traitement itératif est mis en œuvre en sélectionnant au pas d'itération courant k, le point x_jk du jeu de points qui succède dans ledit ordre de distance croissante le point x_{Jk 1} sélectionné au pas d'itération k-1 ;

- ledit traitement itératif est stoppé en fonction d'une comparaison entre la distance entre le point associé x, et le point _ik+i et le double de la distance maximale existant entre le point associé x, et un point du résultat d'intersection calculé au pas d'itération k ;

- dans une étape d'initialisation, l'indice k est fixé à 1 et le résultat d'intersection est fixé égal à l'élément ;

- on réitère l'ensemble des étapes al et b/ en définissant dans chaque nouvelle étape al un jeu de n points tirés d'au moins le jeu précédent ;

- on réitère l'ensemble des étapes al et b/ ; - lors de l'étape de calcul d'un résultat actualisé d'intersection égal à l'intersection entre d'une part la moitié de l'espace E délimitée par un hyperplan médiateur d'un segment (χ,, x_jk ) et qui contient χ,, et d'autre part le résultat d'intersection calculé au pas d'itération précédent, on met en œuvre les étapes suivantes :

- on détermine si un sommet q du résultat actualisé d'intersection est l'intersection entre l'hyperplan médiateur de [x,, x_p] et une arête [qi, q₂] entre deux sommets q1 , q2 de l'élément, où x_p est un point du jeu de n points ;

- et si ledit sommet q est déterminé comme étant l'intersection entre l'hyperplan médiateur de [x,, x_p] et une arête [qi, q₂] entre deux sommets q1, q2 de l'élément, les coordonnées de q sont déterminées à l'aide des formules suivantes :

q = A^_! + A₂q₂, avec Κ = Δ^" (-||χ_ρ||² + 2<q₂, x_p>)

λ₂ = Δ¹(||χ_ρ||² - 2<q_l5 x_p>)

où Δ = -2<qi, x_p> + 2<q₂, x_p>, et

<.,.> représente la fonction produit scalaire ;

- lors de l'étape de calcul d'un résultat actualisé d'intersection égal à l'intersection entre d'une part la moitié de l'espace E délimitée par un hyperplan médiateur d'un segment (x,, x_jk ) et qui contient x,, et d'autre part le résultat d'intersection calculé au pas d'itération précédent, on met en œuvre les étapes suivantes :

- on détermine si un sommet q du résultat actualisé d'intersection est l'intersection entre l'hyperplan médiateur de [x,, x_p], l'hyperplan médiateur de [χ,, x,], et un triangle de sommets ç q₂, q₃ de l'élément, où x_p et x, sont deux points du jeu de n points ;

- et si q est déterminé comme étant l'intersection entre l'hyperplan médiateur de [x,, x_p], l'hyperplan médiateur de [x,, X|], et un triangle de sommets qi, q₂, q₃ de l'élément, les coordonnées de q sont déterminées à l'aide des formules suivantes : q = q₁ + λ₂ q₂ + λ₃ q₃, avec

Κ =A^" [(a₂₃-a₂₂)||x_p||² + (a₁₂-a₁₃)||xi||² + c₃₁ ] ;

A₂ = A^" [(a₂₁ - ¾3)||xp||² + (ai3-aii)||xi||² + c₃2] ;

λ₃ = Δ^"1 [(a₂₂ - a₂₁ ) Il x_p II ² + (ai 1 - a₁₂) || x, || ² + c₃₃ ] où an= - 2<q_1; x_p> ; a₁₂= - 2<q₂, x_p> ; a₁₃= - 2<q₃, x_p>

a₂₁= - 2<q_1; X|> ; a₂₂= - 2<q₂, x,> ; a₂₃= - 2<q₃, x,>.

c₃i = a₂₃ ai₂- a₂₂ ai₃ ; c₃₂= a₂i a-i₃- a₂₃ an ; c₃₃= a₂₂ an - a₂i a-i₂

<.,.> représente la fonction produit scalaire ; - lors de l'étape de calcul d'un résultat actualisé d'intersection égal à l'intersection entre d'une part la moitié de l'espace E délimitée par un hyperplan médiateur d'un segment (χ,, x_jk ) et qui contient χ,, et d'autre part le résultat d'intersection calculé au pas d'itération précédent, comprenant les étapes suivantes :

- on détermine si un sommet q du résultat actualisé d'intersection est un sommet de l'élément ;

- si ledit sommet q est déterminé comme un sommet de l'élément, la présence dudit sommet q dans la moitié de l'espace E délimitée par l'hyperplan médiateur du segment (x,, x_Jk ) et qui contient x, est déterminée en fonction du signe de la fonction orient par : orient ( Y[⁺(i, j_k) ,q) = signe ( || x_Jk || ² -2<q, x_Jk >), où

<.,.> représente la fonction produit scalaire et signe (x) est la fonction fournissant le signe de la variable x ;

- lors de l'étape de calcul d'un résultat actualisé d'intersection égal à l'intersection entre d'une part la moitié de l'espace E délimitée par un hyperplan médiateur d'un segment (x,, x_jk ) et qui contient x,, et d'autre part le résultat d'intersection calculé au pas d'itération précédent :

- on détermine si un sommet q du résultat actualisé d'intersection est l'intersection entre l'hyperplan médiateur de [x,, x_p] et une arête [qi , q₂] entre deux sommets q1 , q2 de l'élément, où x_p est un point du jeu de n points ;

- si q est déterminé comme l'intersection entre l'hyperplan médiateur de [x,, x_k] et une arête [qi , q₂] entre deux sommets q1 , q2 de l'élément, la présence dudit sommet q dans la moitié de l'espace E délimitée par l'hyperplan médiateur du segment Xi,x_jk) et qui contient x, est déterminée en fonction du signe de la fonction orient orient ( f[⁺(/, yJ ,q) = signe (Δ || x_jk || ² -2< A q, x_jk >). signe (Δ), avec Δ = -2<qi , x_p> + 2<q₂, x_p>, et où <.,.> représente la fonction produit scalaire et signe (x) est la fonction fournissant le signe de la variable x ;

- lors de l'étape de calcul d'un résultat actualisé d'intersection égal à l'intersection entre d'une part la moitié de l'espace E délimitée par un hyperplan médiateur d'un segment (x,, x_Jk ) et qui contient x,, et d'autre part le résultat d'intersection calculé au pas d'itération précédent : - on détermine si un sommet q du résultat actualisé d'intersection est l'intersection entre l'hyperplan médiateur de [χ,, x_p], l'hyperplan médiateur de [χ,, X|], et un triangle de sommets qi , q₂, q3 de l'élément, où x_p et X| sont deux points du jeu de n points ;

- si ledit sommet q du résultat actualisé d'intersection est déterminé comme l'intersection entre l'hyperplan médiateur de [x,, x_p], l'hyperplan médiateur de [x,, X|], et un triangle de sommets qi , q₂, q3 de l'élément, la présence dudit sommet q dans la moitié de l'espace E délimitée par l'hyperplan médiateur du segment (x,, x_h- ) et qui contient x, est déterminée en fonction du signe de la fonction orient

( Y[⁺(i _k) .q) définie par : orient ( Y[⁺(i, j_k) ,q) = signe (Δ || x_k || ² -2< Δ q, x_k >). signe (Δ), où

Δ = C31 + C32 + C33 ;

C31 = a₂3 3i2 - a₂2 3i 3 ; C32 = 821 3i3 " ¾3 3i 1 ; C33 = 322 an - 321 3i2 ;

3i i = - 2<qi , x_p> ; 3i2= - 2<q₂, x_P> ; 3i₃= - 2<q₃, x_P>

3₂i= - <qi , X|> ; 3₂₂= - 2<q₂, Xi> ; 3₂₃= - 2<q₃, Xi> et où

<.,.> représente la fonction produit scalaire et signe (x) est la fonction fournissant le signe de la variable x.

Suivant un deuxième aspect, la présente invention propose un programme d'ordinateur de traitement de données définissant un élément dans un espace E de dimensions d, ledit programme comportant des instructions pour mettre en œuvre les étapes d'un procédé suivant le premier aspect de l'invention, lors d'une exécution du programme par des moyens de traitement.

Ces caractéristiques et avantages de l'invention apparaîtront à la lecture de la description qui va suivre, donnée uniquement à titre d'exemple, et faite en référence aux dessins annexés, sur lesquels :

- la figure 1 représente en partie A un champ d'anisotropie prescrit sur une surface 2D, en partie B un diagramme de Voronoï résultant et en partie C le résultat de la triangulation correspondante de Delaunay ;

- la figure 2 représente en partie A une surface isotrope D, en partie B un diagramme de Voronoï résultant et en partie C le résultat de la triangulation correspondante de Delaunay ;

- la figure 3 représente en partie A un diagramme de Voronoï d'un ensemble de points, en partie B la configuration obtenue après une itération d'un algorithme de relaxation de Lloyd, en partie C la configuration obtenue après 100 itérations d'un algorithme de relaxation de Lloyd et en partie D le résultat de la triangulation de Delaunay résultant de la configuration représentée en partie C ;

- la figure 4 est une vue d'un dispositif de traitement de données dans un mode de réalisation de l'invention ;

- la figure 5 représente des étapes d'un procédé dans un mode de réalisation de l'invention ;

- la figure 6 représente des étapes d'un procédé dans un mode de réalisation de l'invention ;

- la figure 7 représente des étapes d'un procédé dans un mode de réalisation de l'invention.

Diagramme de Voronoï et triangulation de Delaunay :

Tout d'abord, les définitions d'un diagramme de Voronoï et de la triangulation de Delaunay sont rappelées.

Etant donné un ensemble X de points tel que X = {xi , x₂,..., x_n} avec x, 6 ¾^d , où 9î est l'espace réel et d la dimension considérée, le diagramme de Voronoï Vor(X) est la collection des cellules de Voronoï Vor(x,), i= 1 à n.

La cellule de Voronoï Vor(x,) associée au point x, est définie par :

Vor(Xi)= {x 6 ¾^d /d(x,Xi) < d(x,Xj) , V j=1 à n} où d(.,.) indique la distance euclidienne ¾^d .

Le dual du diagramme de Voronoï Vor(X), i.e. le complexe simplicial abstrait déduit de la combinatoire de Vor(X), est appelé la triangulation de Delaunay.

Chaque paire de cellules de Voronoï Vor(x,), Vor(Xj) , qui a une intersection non nulle, définit un côté (x,, Xj) dans la triangulation de Delaunay, et chaque triplet de cellules de Voronoï Vor(x,), Vor(Xj), Vor(x_k) qui a une intersection non nulle définit un triangle (x,, Xj, Xk) -

La triangulation de Delaunay présente plusieurs propriétés géométriques intéressantes et est utilisée dans de nombreuses applications, par exemple, mais non seulement, les traitements de génération de maillage (cf. Jean-Daniel Boissonnat and Mariette Yvinec. Algorithmic Geometry. Cambridge University Press, 1998).

Le diagramme de Voronoï barycentrique ou CVT (en anglais « Centroidal Voronoï Tessallation ») est un diagramme de Voronoï dans lequel le point x,, pour i= 1 à n, est le centre de la cellule de Voronoï associée à x,. Une CVT peut être réalisée par un algorithme, dit de relaxation de Lloyd, qui déplace itérativement chaque point x, au centre de Vor(Xi), et fournit des triangles isotropes (cf. Stuart P. Lloyd. Least squares quantization in PCM. IEEE Transactions on Information Theory, 28(2):129-137, 1982).

En référence à la figure 3, un diagramme de Voronoï en 2D associé aux points x,

(correspondant aux points non gras) est représenté en partie A. Les centres des cellules de Voronoï sont en points gras. En partie B de la figure 3, la configuration obtenue après une itération de l'algorithme de relaxation de Lloyd est représentée. En partie C de la figure 3, la configuration obtenue après 100 itérations de l'algorithme de relaxation de Lloyd est représentée. En partie D de la figure 3, la triangulation de Delaunay déduite de la configuration de la partie C de la figure 3 est représentée.

A présent, les notions de diagramme de Voronoï restreint et de triangulation de Delaunay restreinte sont définies. Soit un domaine S inclus dans ¾^d .

Par exemple, S est une surface ou un volume, ou tout autre domaine.

Le diagramme de Voronoï de l'ensemble X restreint à S, libellé Vor(X)|_s, est l'ensemble des cellules de Voronoï Vor(x,), i= 1 à n, restreintes à S, libellées Vor(Xi)|_s.

La cellule de Voronoï restreinte à S, Vor(Xi)|_s, est égale à Vor(Xi)|_s = Vor(x,) Π S. Le dual du diagramme de Voronoï Vor(X), i.e. le complexe simplicial abstrait déduit de la combinatoire de Vor(X) |_s, est appelé la triangulation de Delaunay restreinte. Chaque triangle d'une triangulation de Delaunay restreinte correspond à trois cellules de Voronoï restreintes ayant une intersection non vide.

Conséquence directe de leur définition, les cellules de Voronoï restreintes peuvent aussi être définies par Vor(Xi)|_s = S Π Q Π (', /) ,

y=1 à n,j≠i où ]^⁺(/, y) = {x|d(x, x, ) < c/(x, x_y )} est le demi-espace dans 9î^d limité par l'hyperplan médiateur du segment [x,, Xj] et qui contient x,.

On voit ainsi qu'une cellule de Voronoï restreinte Vor(Xi)|_s peut être obtenue en partant de l'espace 3i^d et en effectuant itérativement des découpes à l'aide des hyperplans Π⁺('>/) pour j = 1 à n et j≠i.

Appliquer une telle itération en pratique n'a pas d'intérêt car l'algorithme qui en résulterait aurait une complexité superquadratique. C'est pourquoi les algorithmes existants sont basés sur d'autres considérations. On connaît par exemple les méthodes de calcul dediagrammes de Voronoi fondés sur les propriétés de la triangulation de Delaunay décrits dans l'ouvrage de référence suivan t: "Algorithmic Geometry", Boissonnat et Yvinec, ISBN-13: 978-0521565295. Ces algorithmes classiques calculent directement la triangulation de Delaunay, en insérant les points un par un et en corrigeant itérativement le maillage de manière à vérifier la propriété dite de la "sphère vide".

Principe du rayon de sécurité :

Un des aspects de l'invention est de permettre de déterminer efficacement, parmi les hyperplans médiateurs définis par les segments [x,, x,], lesquels sont contributeurs à la détermination de cellule de Voronoï Vor(x,) et lesquels sont non contributeurs, i.e. Vor(x,)

Considérons x_h- ,... x_{Jn 1} les n-1 points de l'ensemble X autres que x, et classés par ordre croissant de distance à x,.

Soit V_k(Xi) l'intersection des k premiers hyperplans médiateurs situés entre x, et chacun de ces k premiers points, et R_k son rayon centré sur x,, soit :

V_k(Xi) = ΠΐΓ(^/', 7/ ) ^{et Rk = max} t ^d(^Xi'^x) ^{1 6 ν}*( ) >^■

/=1

Le théorème du rayon de sécurité selon l'invention est le suivant :

Pour tout j tel que d(x„Xj) > 2.R_k, l'hyperplan j) est non contributeur à la construction de la cellule de Voronoï Vor(x,), i.e. V_k(x,) c ·

Preuve : en considérant x 6 V_k(x,) et Xj tel que d(Xi,x_jy) > 2 R_k,

par définition de R_k, d(x,x,) < R_k.

On a d(x,,x) + d(x,Xj)> d(x_{i ;}Xj) (inégalité triangulaire)

par conséquent, d(x,Xj)> R_k > d(x,x,) et x C ]^[⁺(/^', y) .

On notera que ce théorème a des applications pratiques avantageuses lorsqu'on l'applique à des cellules de Voronoï bornées, et par conséquent à des cellules de Voronoï restreintes.

Conséquence directe du théorème du rayon de sécurité :

si d(Xj,x_jk+1) > 2.R_k alors V_k(Xi) = Vor(Xi) . On appelle rayon de sécurité la première valeur de Rk (i.e. la plus grande) rencontrée, en parcourant les x_h- ,... χ_{ί n i} dans cet ordre, qui satisfait cette condition

« d(x_hx_jk+1) > 2.R_k » . Exemple d'application du principe du rayon de sécurité

Dans ce qui suit, on expose un exemple d'application du théorème du rayon de sécurité à un dispositif et un procédé de traitement de données, ces données définissant un domaine S0.

Dans le cas considéré, S0 est une surface s'étendant dans un espace de dimension dO, par exemple dO = 3.

Dans d'autres cas, S0 pourra être un volume délimité par une surface ou tout autre type de domaine.

Dans le cas considéré, le traitement réalisé a pour but la génération d'un maillage anisotrope du domaine S0 conformément à un champ d'anisotropie prescrit.

A cette fin, on considère un dispositif 10 de traitement de données représenté en figure 4.

Un tel dispositif de traitement 10 comporte une mémoire 1 1 , un microcalculateur

12 et une interface homme-machine 13, comprenant notamment un écran d'affichage sur lequel afficher un maillage généré pour le domaine S0.

La mémoire 10 comporte notamment des données numériques de définition du domaine S0 et un programme d'ordinateur P.

Le programme P comprend des instructions logicielles, qui lorsqu'elles sont exécutées par le microcalculateur 12, mettent en œuvre les étapes indiquées ci-dessous en référence aux figures 5-7.

Soit le domaine S0 de dimension dO, pour lequel un maillage anisotrope gouverné par une métrique d'anisotropie donnée doit être généré.

Dans une étape 100, on fait correspondre au domaine S0 de dimension dO, un domaine S de dimensions d > dO conformément aux documents nommés références_1 plus haut. Par exemple, d = 6, ou 10.

Pour dO = 3 et d= 6, à tout point [x,y,z] de S0, on fait par exemple correspondre le point [x,y,z,s.N_x, s.N_y, s.N_z] de S, où le facteur s indique le degré d'anisotropie souhaité

(une petite valeur de s génère un maillage isotrope et une grande valeur de s, un maillage anisotrope). N_x, N_y, N_z sont les vecteurs unitaires normaux à la surface S0 en le point

[x,y,

Puis un traitement itératif 101 est mis en œuvre pour déterminer un maillage isotrope du domaine S. Ainsi dans une étape 102, la valeur n du nombre de cellules de Voronoï étant fixée, un ensemble X de n points de l'espace ¾^d , {xi , x₂,..., x_n} avec x, 6 ¾^d , est déterminé pour l'itération courante.

Lors de la première itération, l'ensemble X de points est par exemple choisi aléatoirement.

Puis lors des itérations suivantes, l'ensemble X est déterminé en fonction des résultats de la dernière itération réalisée pour le traitement 101 .

Dans une étape 103, une étape de détermination du diagramme de Voronoï restreint à S, i.e. Vor(Xj) |_s, est réalisée.

Dans une étape 104, une condition de sortie de boucle est testée. Cette condition de sortie de boucle comprend par exemple :

- la comparaison entre la valeur courante k d'itération et un seuil maximal fixé ; et/ou

- la norme du gradient d'une fonction « objectif » est inférieure à un certain seuil. Par exemple, la fonction « objectif » représente la puissance du bruit de l'échantillonnage

(« quantization noise power » en anglais), à savoir la somme des moments d'inertie des cellules de Voronoi.

Si la condition de sortie de boucle 105 n'est pas satisfaite, la valeur de k est augmentée de 1 et une itération supplémentaire du traitement 101 est mise en œuvre.

Si la condition de sortie de boucle 104 est satisfaite, le traitement itératif 101 est stoppé. Une étape 105 de triangulation de Delaunay restreinte est alors mise en œuvre, pour déterminer le dual du diagramme de Voronoï fourni en sortie du traitement itératif 101 , par déduction directe de la combinatoire du diagramme de Voronoi restreint, à savoir pour chaque sommet du diagramme de Voronoi restreint, on génère le triangle de Delaunay correspondant.

Des opérations diverses peuvent ensuite être mises en œuvre, notamment la projection sur l'espace de dimension dO initiale du résultat de la triangulation de Delaunay obtenue dans l'espace de dimension d à l'issue du traitement itératif 101 , de manière à obtenir le maillage anisotrope souhaité du domaine S0.

Dans le mode de réalisation considéré, l'étape 103 de détermination du Vor(Xj)|_s, est réalisée de la manière indiquée ci-après en référence à la figure 6 :

Dans une étape 103_1 , pour chaque x,, avec i= 1 à n, les autres points x,, j= 1 à n et j≠ i, sont triés par ordre croissant de distance à x, ce qui résulte en une liste ordonnée Xj , . . . Xj ^■ Pour ce faire, l'outil ANN est par exemple utilisé (David M. Mount and Sunil

Arya, ANN: A library for approximate nearest neighbor searching. In Proceedings CGC Workshop on Computational Geometry, pages 33-40, 1997).

Dans une étape 103_2 d'initialisation, on initialise à la valeur 0 la variable m, et chacune des coordonnées de points g, de ¾^d , i= 1 à n.

Dans une étape 103_3, un découpage du domaine S en sous-domaines est effectué, ceci afin de pouvoir paralléliser les étapes réalisées sur des sous-domaines distincts (triangles, tétraèdres etc).

Dans le cas considéré, les sous-domaines f de la surface S sont des triangles.

Dans une étape 103_4, pour chaque point x,, avec i= 1 à n, et pour chaque sous- domaine f telle que l'intersection entre le sous-domaine f et la cellule de Voronoï associée à X, est non vide, Vor(Xi)|_f, la cellule de Voronoï associée à x, et restreinte au sous- domaine f, est déterminée, de la manière indiquée ci-dessous en référence à la figure 7.

Puis on affecte à une variable m la valeur de la « masse » de Vor(Xj) |_f, on ajoute à la valeur courante de m, la valeur de m, et on actualise la valeur des coordonnées du point g, suivant la formule : g, = g, + m.centre Vor(Xi)|_f, où centre Vor(Xi) |, est le centre de gravité de Vor(Xj)|_f.

La « masse » de Vor(Xi) |, est égale à . J cx . Ainsi la « masse » de Vor(Xi) |, est xeVor

1 r

Le centre de gravite de Vor(Xi)|, est 1— x.dx .

masse de Vor(_Xi )\ _{f Vor}l_{x )l f}

Dans un mode de réalisation, la détermination de l'ensemble de points X dans une étape 102 pour l'itération k+1 du traitement 101 , est par exemple telle que le point x, = g,. 1/rrii, où i= 1 à n et g, est le point obtenu à l'issue de l'étape 103_4 mise en œuvre pour l'itération k du traitement 101 .

Selon l'invention, l'étape 103_4 de détermination de la cellule de Voronoï Vor(Xi)|, associée à x, et restreinte au sous-domaine f, est mise en œuvre en utilisant le théorème du rayon de sécurité, et en considérant les points x_h x_{jn 1} ordonnés selon une distance à X, croissante.

Dans un exemple de mode de réalisation, en référence à la figure 7, l'étape 103_4 comprend ainsi les opérations suivantes. Dans une étape 103_41 d'initialisation, on considère un domaine V égal au sous- domaine f considéré, une valeur t égale à 1 et une valeur R égale à max {d(x,, x)/ x 6V}.

Dans une étape 103_42, on itère les étapes suivantes tant que d(Xi, _i( ) < 2R et que t <n :

v = v n []⁺(/, i₍) ;

t = t+1 .

Quand d(x,, _y )≥ 2R ou que t = n, l'étape 103_42 est stoppée.

Et la valeur courante de V est alors égale à la cellule de Voronoï Vor(Xi)|, associée à x, et restreinte au sous-domaine f (poursuivre les itérations de l'étape 103_42 n'est plus utile puisque les hyperplans médiateurs non encore pris en compte ne sont pas contributeurs).

Pour le calcul effectif de l'intersection V y^' _f ) , l'algorithme de fenêtrage réentrant de Sutherland& Hodgman pourra par exemple être utilisé (Ivan Sutherland and Gary W. Hodgman, Reentrant polygon clipping. Communications of the ACM, 17:32-42, 1974).

L'utilisation du théorème du rayon de sécurité permet de ne déterminer que les sections de la cellule de Voronoï qui sont utiles pour calculer l'intersection entre cette cellule et le sous-domaine f considéré, réduisant ainsi amplement le volume de calcul nécessaire, ce qui est très appréciable notamment quand la valeur de la dimension d de l'espace considéré augmente.

Par ailleurs, elle autorise le calcul en parallèle et de façon indépendante, des intersections pour des cellules de Voronoï respectives. Dans un mode de réalisation, un procédé selon l'invention met en œuvre les étapes indiquées ci-dessous, pour définir l'intersection entre un sous-domaine f et une cellule de Voronoï, dans un espace de dimensions d, en utilisant les produits scalaires et la combinaison linéaire de vecteurs, le volume de calcul mis en œuvre étant indépendant de la dimension d.

On notera que ces dispositions peuvent être mises en œuvre indépendamment de l'utilisation du test dit du rayon de sécurité.

Dans l'étape de calcul de V y^' _i ) à l'aide l'algorithme de fenêtrage ré-entrant de Sutherland & Hodgman, il est nécessaire de déterminer les coordonnées des nouveaux sommets du domaine résultant de cette intersection, et également de déterminer de quel côté un sommet se situe par rapport à un demi-espace ]^[⁺(/^', y) .

La fonction permettant de déterminer de quel côté un sommet q se situe par rapport au demi-espace Π⁺('>/) ^est nommée orient ( ]^[⁺(/, y) ,q). Si le résultat de orient ( YY(i, y) ,q) est positif, le sommet q est situé dans le demi-espace Π (/, y)■

Sinon, le sommet q est à l'extérieur du demi-espace Π⁺('>/)■

Des sommets q du domaine résultant dans certains cas ont déjà été déterminé lors d'itération(s) précédente(s) ou n'apparaissent que dans l'itération courante.

Trois configurations différentes existent pour un sommet q :

a- q est un sommet du sous-domaine f considéré ;

b- q est l'intersection entre l'hyperplan médiateur de [x,, x_k] et une arête [qi , q₂] entre deux sommets qi , q₂ du sous-domaine f ;

c- q est l'intersection entre deux hyperplans médiateurs (par exemple l'hyperplan médiateur de [x,, x_k] et celui de [χ,, x,], et le sous-domaine f ; dans le cas présent, le sous-domaine est le triangle de sommets qi , q₂, q₃.

Sans perte de généralité, on suppose x, à l'origine de l'espace considéré. On s'y ramène en appliquant une translation.

Dans le cas a/, les coordonnées de q sont connues, et

orient ( ]^[⁺(/, y) ,q) = signe ( || x_j || ² -2<q, Xj>), où <,> représente la fonction produit scalaire.

Dans le cas b/, les coordonnées barycentriques (λ_{1 ;} λ₂) de q par rapport à qi , q₂ sont déterminées : q = qi + λ₂ q₂ .

En nommant Δ le déterminant de la matrice 2x2 ci-dessus :

Δ = -2<q_{1 ;} x_k> + 2<q₂, x_k>

= Δ ¹ (- Il x_k H ² + 2<q₂, x_k>) et A₂ = A- ( || x_k || ² - 2<q₁ , x_k>)

et orient ( Π }) ,q) = signe (Δ || Xj || ² -2< Δ q, Xj>). signe (Δ).

Dans le cas c/, les coordonnées barycentriques (λι , λ₂, λ₃) de q par rapport à qi , q₂, q₃ sont déterminées : q = λι qi + λ₂ q₂ + λ₃ q₃.

an= - 2<q_1; x_k> ; a₁₂= - 2<q₂, x_k> ; a₁₃= - 2<q₃, x_k>

a₂₁= - 2<q_1; X|> ; a₂₂= - 2<q₂, x,> ; a₂₃= - 2<q₃, x,>.

Les cofacteurs suivants sont calculés :

c₃i = a₂₃3i₂- a₂₂ ai₃ ; c₃₂= a₂i a-i₃ - a₂₃ an ; c₃₃= a₂₂ an - a₂i a-i₂

Δ étant le développement du déterminant par rapport à la dernière ligne de la matrice 3x3 : Δ = c₃i + c₃₂ + c₃₃.

Il en découle :

λι =A^" [(a₂₃-a₂₂) ||x_k||²+ (a₁₂-a₁₃) ||xi ||² + c₃₁ ] ;

A₂ = A^" [(a₂₁ -a₂₃) ||x_k||²+ (a₁₃-a₁₁) ||xi ||² + c₃₂] ;

λ₃ = Δ^"1 [(a₂₂ - a₂₁ ) Il x_k II ² + (a! -a^) || xi || ² + c₃₃ ] ; et orient ( Π (', j) ,q) = signe (Δ || x || ² -2< Δ q, Xj>). signe (Δ). La formulation de q et de la fonction orient est donc indépendante de la dimension d, en ceci que seuls le produit scalaire <.,.> et des combinaisons linéaires interviennent, et où la dimension des systèmes linéaires est indépendante de d.

En outre, elle se prête bien à une implémentation en arithmétique à précision arbitraire et à un filtrage (par arithmétique d'intervalle ou filtres quasi-statiques).

Le recentrage (translation de x, à l'origine) améliore les performances de filtrage.

L'algorithme mettant en œuvre ces étapes de calcul des coordonnées de q et de la fonction orient se paramétrise par un noyau géométrique, définissant les types de points et de vecteurs, le produit scalaire et les combinaisons linéaires de vecteurs. Ceci permet d'avoir une implémentation fonctionnant quelle que soit la dimension d.

Dans l'exemple décrit ci-dessus, on a considéré que le sous-domaine f est un triangle. Toutefois, les formules indiquées ci-dessus sont valables également pour tout sous-domaine f, par exemple un polygone arbitraire ou tout objet de dimension supérieure, par exemple un tétraèdre (ces formules sont donc applicables dans le cas d'un maillage volumique).

L'utilisation du théorème du rayon de sécurité pour diminuer les calculs d'un diagramme de Voronoï restreint ou d'une triangulation de Delaunay restreinte a été décrite ci-dessus dans le cadre d'une application au maillage d'un élément. Bien sûr, elle peut intervenir dans d'autres domaines mettant en œuvre de tels calculs et procure les mêmes avantages.

Il en va de même pour les principes de calcul des sommets résultant des intersections indépendamment de la dimension d, et de la fonction orient.

Ces dispositions sont très avantageuses, et ce pour des valeurs de d égales à 2, 3, ou supérieures ou égales à 4.

Les applications d'une technique selon l'invention autre que le maillage comprennent par exemple :

- des applications de reconstruction 3D, notamment les étapes de passage d'un nuage de points à une surface ;

- des applications de robotique : par exemple dans la représentation d'espaces de configuration, ou encore la gestion de flottilles de robots autonomes ou la planification de trajectoires ;

- des applications astrophysiques : par exemple, calculs dans des espaces des phases ;

- des applications d'intelligence artificielle, telles que des recherches de documents (type "google", qui utilise aussi des espaces de grandes dimensions).

Claims

REVENDICATIONS

1 . Procédé de traitement de données définissant un élément dans un espace E de dimensions d, ledit procédé comprenant les étapes suivantes :

al on définit un jeu de n points dans ledit espace ;

b/ on calcule l'intersection entre ledit élément et des cellules de Voronoï déterminées chacune à partir d'un point respectivement associé parmi lesdits n points définis ;

ledit procédé étant caractérisé en ce que dans l'étape b/, on détermine le résultat de l'intersection entre ledit élément et une cellule de Voronoï déterminée à partir du point qui lui associé χ,, en mettant en œuvre un traitement itératif selon lequel

à un pas d'itération courant k, on sélectionne un point supplémentaire x_h- du jeu de points autre que le point x, et on calcule un résultat actualisé d'intersection égal à l'intersection entre d'une part la moitié de l'espace E délimitée par l'hyperplan médiateur du segment (χ,, x_jk ) et qui contient χ,, et d'autre part le résultat d'intersection calculé au pas d'itération précédent,

2. Procédé selon la revendication 1 , selon lequel dans l'étape b/, les n-1 autres points du jeu classés par ordre de distance croissante par rapport au point x,, sont libellés _y- _{in i} , et le traitement itératif est mis en œuvre en sélectionnant au pas d'itération courant k, le point x_h- du jeu de points qui succède dans ledit ordre de distance croissante le point x_h- ₁ sélectionné au pas d'itération k-1 .

3. Procédé selon la revendication 2, selon lequel ledit traitement itératif est stoppé en fonction d'une comparaison entre la distance entre le point associé x, et le point _ik+i et le double de la distance maximale existant entre le point associé x, et un point du résultat d'intersection calculé au pas d'itération k.

4. Procédé selon l'une des revendications précédentes, selon lequel dans une étape d'initialisation, l'indice k est fixé à 1 et le résultat d'intersection est fixé égal à l'élément.

5. Procédé selon l'une des revendications précédentes, selon lequel on réitère l'ensemble des étapes al et b/ en définissant dans chaque nouvelle étape a/ un jeu de n points tirés d'au moins le jeu précédent.

6. Procédé selon l'une des revendications précédentes, selon lequel on réitère l'ensemble des étapes al et b/.

7. Procédé selon l'une des revendications précédentes, selon lequel, lors de l'étape de calcul d'un résultat actualisé d'intersection égal à l'intersection entre d'une part la moitié de l'espace E délimitée par un hyperplan médiateur d'un segment (x,, x_h- ) et qui contient

Xi , et d'autre part le résultat d'intersection calculé au pas d'itération précédent, on met en œuvre les étapes suivantes :

- et si ledit sommet q est déterminé comme étant l'intersection entre l'hyperplan médiateur de [x,, x_p] et une arête [qi, q₂] entre deux sommets q1 , q2 de l'élément, les coordonnées de q sont déterminées à l'aide des formules suivantes :

q = λι qi + A₂ q₂ , avec λι = Δ^" (- || χ_ρ || ² + 2<q₂, x_p>)

λ₂ = Δ ¹ ( Il x_P II ² - 2<q_! , x_p>)

où Δ = -2<q _; x_p> + 2<q₂, x_p>, et

<.,.> représente la fonction produit scalaire.

8. Procédé selon l'une des revendications précédentes, selon lequel, lors de l'étape de calcul d'un résultat actualisé d'intersection égal à l'intersection entre d'une part la moitié de l'espace E délimitée par un hyperplan médiateur d'un segment (x,, x_h- ) et qui contient

- on détermine si un sommet q du résultat actualisé d'intersection est l'intersection entre l'hyperplan médiateur de [x,, x_p], l'hyperplan médiateur de [x,, X|], et un triangle de sommets qi , q₂, q₃ de l'élément, où x_p et X| sont deux points du jeu de n points ;

- et si q est déterminé comme étant l'intersection entre l'hyperplan médiateur de [x,, x_p], l'hyperplan médiateur de [χ,, x,], et un triangle de sommets qi , q₂, q₃ de l'élément, les coordonnées de q sont déterminées à l'aide des formules suivantes : q = λι qi + λ₂ q₂ + λ₃ q₃, avec

λι = A^"1 [(a₂₃ - a₂₂) || x_P || ² + (ai2 - a₁₃) || xi || ² + c₃i ] ;

A₂ = A^~ [(a₂₁ - a₂₃) || xp || ² + (a₁3 - a_{1 1}) || xi || ² + c₃₂ ] ;

λ₃ = Δ^"1 [(a₂₂ - a₂₁ ) || x_P || ² + (ai ₁ - a₁₂) || x, || ² + c₃₃ ] où an= - 2<qi , x_p> ; a₁₂= - 2<q₂, x_p> ; a₁₃= - 2<q₃, x_p>

a₂₁= - 2<qi , xi> ; a₂₂= - 2<q₂, x,> ; a₂₃= - 2<q₃, x,>.

C31 = ¾3 3ΐ₂ " ¾₂ 3.13 ! C₃₂ = 3₂₁ 3i3 " 3₂₃ &_{1 1} ^', C₃₃ = a₂₂ a-n - a₂i Ά₁₂

<.,.> représente la fonction produit scalaire.

9. Procédé selon l'une des revendications précédentes, selon lequel, lors de l'étape de calcul d'un résultat actualisé d'intersection égal à l'intersection entre d'une part la moitié de l'espace E délimitée par un hyperplan médiateur d'un segment (x,, x_h- ) et qui contient

Xi , et d'autre part le résultat d'intersection calculé au pas d'itération précédent, on met en oeuvre les étapes suivantes :

- si ledit sommet q est déterminé comme un sommet de l'élément, la présence dudit sommet q dans la moitié de l'espace E délimitée par l'hyperplan médiateur du segment (x,, x_h- ) et qui contient x, est déterminée en fonction du signe de la fonction orient ( ]^[⁺(/, y_fc ) ,q) définie par : orient ( Y[⁺(i, j_k ) ,q) = signe ( || x_jk || ² -2<q, x_k >), où

<.,.> représente la fonction produit scalaire, et signe (x) est la fonction fournissant le signe de la variable x.

10. Procédé selon l'une des revendications précédentes, selon lequel, lors de l'étape de calcul d'un résultat actualisé d'intersection égal à l'intersection entre d'une part la moitié de l'espace E délimitée par un hyperplan médiateur d'un segment (x,, x_h- ) et qui contient x,, et d'autre part le résultat d'intersection calculé au pas d'itération précédent, on met en oeuvre les étapes suivantes :

- on détermine si un sommet q du résultat actualisé d'intersection est l'intersection entre l'hyperplan médiateur de [x,, x_p] et une arête [qi , q₂] entre deux sommets q1 , q2 de l'élément, où x_p est un point du jeu de n points ; - si q est déterminé comme l'intersection entre l'hyperplan médiateur de [χ,, x_k] et une arête [qi , q₂] entre deux sommets q1 , q2 de l'élément, la présence dudit sommet q dans la moitié de l'espace E délimitée par l'hyperplan médiateur du segment (Xi,x_jk) et qui contient x, est déterminée en fonction du signe de la fonction orient ( ]^[⁺(/, y_fc) ,q) définie par : orient ( f[⁺(/, y_A) ,q) = signe (Δ || x_Jk || ² -2< A q, x_Jk >). signe (Δ), avec Δ = -2<q_{1 ;} x_p> + 2<q₂, x_p>, et où <.,.> représente la fonction produit scalaire, et signe (x) est la fonction fournissant le signe de la variable x.

1 1 . Procédé selon l'une des revendications précédentes, selon lequel lors de l'étape de calcul d'un résultat actualisé d'intersection égal à l'intersection entre d'une part la moitié de l'espace E délimitée par un hyperplan médiateur d'un segment (x,, x_jk ) et qui contient x,, et d'autre part le résultat d'intersection calculé au pas d'itération précédent, on met en oeuvre les étapes suivantes :

- si ledit sommet q du résultat actualisé d'intersection est déterminé comme l'intersection entre l'hyperplan médiateur de [x,, x_p], l'hyperplan médiateur de [χ,, x,], et un triangle de sommets qi , q₂, c ₃ de l'élément, la présence dudit sommet q dans la moitié de l'espace E délimitée par l'hyperplan médiateur du segment (x,, x_jk ) et qui contient x, est déterminée en fonction du signe de la fonction orient

Δ = C31 + C3₂ + C33 ;

C31 = a₂3 ai₂ - a₂₂ ai 3 ; C3₂ = a₂i a-i3 - a₂3 an ; 033 = a₂₂ an - a₂i a-i₂ ;

an = - 2<q_{1 ;} x_p> ; a₁₂= - 2<q₂, x_p> ; a₁₃= - 2<q₃, x_p>

a₂₁= - 2<q_{1 ;} X|> ; a₂₂= - 2<q₂, x,> ; a₂₃= - 2<q₃, Xi> et où

12. Programme d'ordinateur de traitement de données définissant un élément dans un espace E de dimensions d, ledit programme comportant des instructions pour mettre en œuvre les étapes d'un procédé selon l'une des revendications précédentes, lors d'une exécution du programme par des moyens de traitement.