Dispositif informatique pour la simulation d'un ensemble d'objets en interaction et procédé correspondant

L'invention concerne la simulation d'un ensemble d'objets en interaction.

La simulation d'objets en interaction est un domaine vaste qui s'étend de la simulation d'une multitude de corps mécaniques plus ou moins complexes à la simulation du comportement de molécules, par exemple des protéines.

Les premières approches développées ont eu pour principe le rassemblement des objets, c'est-à-dire leur rassemblement par groupe afin de définir un comportement global. Ce type d'approche est intéressant grâce à ses coûts de calcul relativement faibles.

Cependant, la réduction du coût de calcul se fait au détriment du réalisme de la simulation dans la mesure où de nombreux aspects relatifs aux relations des objets groupés sont volontairement ignorés ou réduits à leur plus simple expression.

D'autres approches ont consisté à considérer des objets en interaction comme des corps articulés. Un corps articulé est généralement composé de deux corps reliés entre eux par un joint. Les deux corps peuvent être eux-mêmes des corps articulés et ainsi de suite. Cette approche a permis de modéliser avec fidélité le comportement mécanique de tels ensembles d'objets tout en conservant un niveau de performance élevé.

On parle de comportement "mécanique" car les objets interagissent essentiellement par le biais de l'articulation ou joint qui les rejoint. Il suffit alors de modéliser les articulations pour modéliser l'ensemble d'objets.

Lorsque les interactions entre les objets deviennent plus complexes, ces approches deviennent assez rapidement inefficaces, car les besoins/coûts de calcul qui leur sont associés deviennent prohibitifs.

L'invention vient améliorer la situation.

A cet effet, l'invention propose un dispositif informatique pour la simulation d'un ensemble d'objets en interaction qui comporte une mémoire contenant une représentation arborescente de l'ensemble d'objets. Dans cette représentation arborescente, chacun des intermédiaires est associé à des données dynamiques, des données géométriques et des données d'interaction.

Le dispositif comporte en outre un contrôleur de simulation qui actionne selon un cycle répétitif :

- un distributeur de données d'interaction sur les objets et sous-ensembles d'objets,

- un mécanisme de mise à jour des données dynamiques avec parcours de la représentation arborescente pour noeuds des sujets à interactions en fonction des données d'interaction et des données géométriques concernées, et

- un mécanisme de mise à jour des données géométriques, avec parcours de la représentation arborescente pour des nœuds sujets à interactions, en fonction des données dynamiques mises à jour.

Afin de tenir compte des interactions entre noeuds indépendantes des articulations, la mémoire comporte en outre des données d'interaction locale associées à certains au moins des noeuds, et le distributeur de données d'interactions comprend un mécanisme de mise à jour d'interaction, avec parcours de la représentation arborescente pour des sujets à interaction, pour mettre à jour une donnée d'interaction d'un noeud en fonction des données d'interaction locale de nœuds fils concernés.

Un tel dispositif est particulièrement avantageux car il offre la possibilité de prolonger l'approche de modélisation par corps articulés, tout en tenant compte de forces non exclusivement articulaires, sans augmenter de manière significative les temps de calcul, et d'accroître par conséquent l'exploitabilité de la simulation.

L'invention concerne également un procédé de simulation du comportement d'un ensemble d'objets en interaction dans lequel on maintient une représentation arborescente de l'ensemble d'objets. Dans la représentation, chaque nœud intermédiaire est associé à des données dynamiques, des données géométriques et des données d'interaction dépendantes de données de ses nœuds fils.

Le procédé met en oeuvre répétitivement un cycle comprenant les étapes suivantes :

a. distribuer les données d'interaction sur certains des objets et sous-ensembles d'objets,

b. parcourir la représentation arborescente pour des nœuds sujets à interaction, tout en mettant à jour des données dynamiques de ces noeuds, en fonction des données d'interaction et des données géométriques concernées, et

c. parcourir à nouveau la représentation arborescente pour des nœuds sujets à interaction, tout en mettant à jour des données géométriques, en fonction des données dynamiques telles que mises à jour à l'étape b.

Plus spécifiquement, l'étape a. comprend :

al. un parcours préalable de la représentation arborescente pour des nœuds sujets à interaction, avec la mise à jour des données d'interaction de ces noeuds, en fonction de données d'interaction locale de nœuds fils concernés.

D'autres caractéristiques et avantages de l'invention apparaîtront mieux à la lecture de la description qui suit tirée d'exemples donnés à titre illustratif et non limitatif à partir des dessins sur lesquels :

- la figure 1 représente un diagramme schématique d'un dispositif informatique selon l'invention ;

- la figure 2 représente un diagramme bloc d'une implémentation par des fonctions d'une boucle de simulation effectuée par le dispositif de la figure 1 ;

- la figure 3 représente un diagramme bloc d'une première fonction de la boucle de simulation de la figure 2 ;

- la figure 4 représente un diagramme bloc d'une fonction de la figure 3 ;

- les figures 5 et 6 représentent schématiquement deux fonctions de la figure 4 ;

- la figure 7 représente un diagramme bloc d'une autre fonction de la figure 3 ;

- la figure 8 représente un diagramme bloc d'une deuxième fonction de la boucle de simulation de la figure 2 ;

- la figure 9 représente un diagramme bloc d'une troisième fonction de la boucle de simulation de la figure 2 ;

- la figure 10 représente un diagramme bloc d'une quatrième fonction de la boucle de simulation de la figure 2 ; et

- la figure 11 représente un diagramme bloc d'une cinquième fonction de la boucle de simulation de la figure 2.

Les dessins et la description ci-après contiennent, pour l'essentiel, des éléments de caractère certain. Ils pourront donc non seulement servir à mieux faire comprendre la présente invention, mais aussi contribuer à sa définition, le cas échéant.

La présente description est de nature à faire intervenir des éléments susceptibles de protection par le droit d'auteur et/ou le copyright. Le titulaire des droits n'a pas d'objection à la reproduction à l'identique par quiconque du présent document de brevet

ou de sa description, telle qu'elle apparaît dans les dossiers officiels. Pour le reste, il réserve intégralement ses droits.

L'invention couvre également, en tant que produits, les éléments logiciels décrits, mis à disposition sous tout "médium" (support) lisible par ordinateur. L'expression "médium lisible par ordinateur" comprend les supports de stockage de données, magnétiques, optiques et/ou électroniques, aussi bien qu'un support ou véhicule de transmission, comme un signal analogique ou numérique.

En outre, la description détaillée est augmentée de l'annexe A, qui donne la formulation de certaines formules mathématiques mises en œuvres dans le cadre de l'invention. Cette Annexe est mise à part dans un but de clarification, et pour faciliter les renvois. Elle est partie intégrante de la description, et pourra donc non seulement servir à mieux faire comprendre la présente invention, mais aussi contribuer à sa définition, le cas échéant.

Une première approche de la simulation mécanique de corps articulés a été proposée dans l'article "A divide-and-conquer articulated body algorithm for parallel o(log(n)) calculation ofrigid body dynamics. part 1: Basic Algorithm", International Journal of Robotics Research 18(9):867-875 par Featherstone, noté ci-après Featherstone '99.

Featherstone '99 décrit un algorithme de découpage d'un corps articulé en un arbre binaire, pour calculer l'ensemble des composantes de force et de position de manière récursive. Par corps articulé, on entend l'ensemble qui comprend une articulation, ou joint, qui relie entre eux deux corps pouvant être eux-mêmes des corps articulés.

La notion de corps articulé ne doit donc pas être restreinte à sa simple connotation mécanique. La notion d'articulation doit également être comprise au sens large. Par exemple une articulation peut comporter un nombre variable de degrés de liberté, compris entre zéro (couplage rigide des deux sous-ensembles) et six (mouvement relatif des deux sous-ensembles non contraint).

L'arbre binaire comporte des nœuds feuilles qui représentent les objets les plus élémentaires de l'objet articulé, et qui sont chacun reliés à un joint, ou noeud. Ainsi, chaque nœud de l'arbre binaire est lui-même un corps articulé qui désigne un sous- ensemble d'objets.

L'article "An Efficient, Error-Bounded Approximation Algorithm for Simulating Quasi- Statics o/Complex Linkages", Proceedings of ACM Symposium on Solid and Physical Modeling, 2005, par Redon et Lin, noté ci-après SPM '05, vient compléter cette approche en définissant une métrique d'accélération des joints.

Cette métrique permet d'évaluer l'accélération du corps articulé associé à un joint donné, à partir des seules données de ce joint, sans connaître les données des feuilles. Il est alors possible de "rigidifier" l'arbre, en annulant l'accélération d'un corps articulé pour lequel la métrique d'accélération est inférieure à un seuil choisi. La rigidification de l'arbre permet ainsi de restreindre le parcours de l'arbre à un nombre de noeuds réduit, tout en maîtrisant la fidélité de la simulation.

Dans l'article "Adaptative Dynamics cf Articulated Bodies", ACM Transactions on Graphics (SIGGRAPH 2005), 24(3). par Redon, Gallopo et Lin, noté ci-après SIGGRAPH '05, la qualité de la simulation est améliorée en traitant également le cas dynamique, et plus seulement quasi-statique, et en en ajoutant une méthode de prise en compte de la rigidification, afin d'éviter de créer des forces n'ayant aucune réalité concrète.

Les modes de réalisation décrits ici tirent avantage de l'ensemble de ces méthodes, et pourront parfois être implémentés comme des extensions de celles-ci, bien que de nombreuses autres implémentations de l'invention peuvent être envisagées.

Pour la simulation, chaque nœud de l'arbre est associé à diverses données. On distingue trois types principaux de données qui seront décrites plus avant par la suite : les données d'interaction, les données dynamiques, eι les données géométriques.

Les données d'interaction peuvent être définies comme suit : elles représentent des forces qui s'appliquent à un nœud donné, ou à un objet donné dans le sous-ensemble d'objets que représente ce nœud. Ces forces peuvent être classées en deux classes de forces :

- les forces de contrainte, liées aux interactions des objets par le biais des articulations, et

- les forces locales, qui englobent toutes les autres forces pouvant s'appliquer à un nœud donné.

Dans l'exemple décrit ici, les objets sont des molécules et des atomes, et les forces locales comprennent d'une part les forces externes qui s'appliquent individuellement à des joints, et d'autre part l'ensemble des forces inter-atomiques et inter-moléculaires, telles que les forces électrostatiques, les forces de van der Waals etc.

Les données dynamiques représentent des données composites caractéristiques de la dynamique des corps articulés. Ces données permettent de calculer des métriques d'accélération et de vitesse de joint sans connaître nécessairement l'ensemble des données de position et de force de l'ensemble des objets. Leur application apparaîtra mieux par la suite.

Les données géométriques représentent un ensemble de données caractéristiques des corps articulés qui dépendent exclusivement des positions relatives des objets entre eux. Ces données permettent notamment de réaliser la plupart des calculs dans des repères locaux, ce qui simplifie la simulation. Ici encore, ces données seront explicitées par la suite.

L'exemple décrit ici s'applique à une simulation quasi-statique des objets. Cela signifie que l'on considère que leur vitesse est nulle à chaque instant.

On a représenté sur la figure 1 une représentation schématique d'un dispositif selon l'invention.

Le dispositif 2 comporte un contrôleur de simulation 4, un bloc de simulation 6 et une mémoire 8. Le bloc de simulation 6 comporte un distributeur de données d'interaction 10, un mécanisme de mise à jour de données dynamiques 12 et un mécanisme de mise à jour de données géométriques 14.

Lors d'une simulation, le contrôleur 4 actionne répétitivement le bloc de simulation 6 selon une boucle de simulation qui sera décrite plus bas sur la base de données contenues dans la mémoire 8.

Le dispositif 2 peut être réalisé sous plusieurs formes. Le mode de réalisation préféré actuel est un programme informatique qui implémente le contrôleur 4 et le bloc 6 sous une forme exécutable sur un ordinateur disponible dans le commerce, par exemple un ordinateur muni d'un microprocesseur cadencé à 1.7 GHz, doté de 1 GB de mémoire RAM, et d'un disque dur de taille suffisante pour contenir un système d'exploitation et le code du programme.

Néanmoins, d'autres formes de réalisation sont possibles, comme la réalisation sous la forme d'un circuit imprimé spécialisé (carte fille, ASIC, etc.) ou d'un processeur dédié (FPGA, microchip, etc.).

On a représenté sur la figure 2 un diagramme schématique d'une implémentation de la boucle de simulation réalisée par le bloc 6. Dans cette boucle, des fonctions Upd lnt () 20, Sprd_Int() 22, Get_Dyn_Coeff() 24, Upd_Dyn() 26, et Upd_Geo() 28, sont appelés successivement.

La boucle représentée sur la figure 2 représente les principales étapes de la simulation dans une boucle de fonctionnement donnée. Ces étapes sont respectivement la mise à jour des interactions inter-objets de nature non articulaire, la diffusion de ces

interactions, l'obtention et la mise à jour des paramètres dynamiques de simulation, puis la mise à jour des paramètres géométriques de simulation.

Bien que les différentes étapes représentées sur la figure 2 correspondent à des fonctions dans le programme implémentant l'invention, il apparaîtra par la suite que ces fonctions comportent elles-mêmes plusieurs sous fonctions, et doivent être considérées comme des composants fonctionnels que l'homme du métier sera libre d'implémenter comme il l'entend.

D'un point de vue fonctionnel, on peut considérer que :

* les fonctions Upd_Int() 20 et Sprd_Int() 22 constituent un noyau du distributeur de données 10,

* les fonctions Get_Dyn_Coeff() 24 et Upd_Dyn() 26 constituent un noyau du mécanisme de mise à jour de données dynamiques 12, et

* la fonction Upd_Geo() 28 constitue un noyau du mécanisme de mise à jour de données géométriques 14.

La figure 3 représente une implémentation possible de la fonction Upd_Int() 20 de la figure 2. Comme on peut le voir sur cette figure, la fonction Upd_Int() 20 comporte trois sous fonctions, Upd_BB() 30, Upd_Int_L() 32 et Upd_Int_Val() 34.

Dans l'approche classique, dès lors qu'un sous-ensemble d'objets est rigidifié, celui-ci cesse temporairement d'être traité dans la simulation. Comme les forces considérées dans ce type de simulation sont de type mécanique, ce sous-ensemble peut être très proche dans l'espace d'un autre sous-ensemble d'objets sans que cela n'ait de conséquence.

Lorsque l'on souhaite simuler le comportement d'une protéine ou d'une chaîne moléculaire, des forces locales de type interatomiques du type électrostatique ou de van der Waals viennent s'ajouter aux forces de contrainte et extérieures.

Dès lors, la simulation "mécanique" n'est plus suffisante, car il devient possible qu'un sous-ensemble rigidifié vienne interagir avec un autre sous-ensemble. Il est alors nécessaire de le "dé-rigidifïer" pour tenir compte de ces nouvelles forces locales, qui ne dépendent que de la position relative des objets entre eux. Cela n'est pas possible à réaliser dans l'approche existante, à moins d'en abandonner tous les avantages.

La fonction Upd_Int() 20 permet de tenir compte de ces forces locales. Pour cela, cette fonction part d'une situation moléculaire donnée, et permet de déterminer dans un premier temps quels sont les objets susceptibles d'interagir entre eux. Dans un deuxième temps, cette fonction calcule les forces locales correspondantes, et intègre ces données de manière ordonnée.

La fonction Upd_BB() 30 est une adaptation de la détection d'interaction à la simulation moléculaire. Comme les forces interatomiques sont principalement à courte distance, cette fonction définit des espaces d'interaction potentielle autour de chaque objet.

Plus précisément, ces espaces sont définis comme des "boîtes orientées" entourant un objet ou groupe d'objets. Les boîtes elles-mêmes sont définies pour chaque nœud de l'arbre binaire au début de la simulation, et la fonction Upd_BB() 30 opère en mettant à jour la position d'un point de référence et une direction d'orientation pour chaque boîte.

Les principes des boîtes orientées, appelées "Oriented Bounding Box" ou OBB dans la littérature, sont décrits dans l'article "Obbtree: a hierarchical structure for rapid interférence détection" par Gottschalk et al. dans ACM Translations on Graphics (SIGGRAPH 1996).

Comme les OBB dépendent exclusivement des positions des objets qu'elles entourent, et comme elles sont exprimées dans un référentiel lié au nœud auquel elles sont rattachées,

elles sont inchangées lorsqu'un nœud est rigidifié. La fonction Upd_OBB() 30 opère donc uniquement sur les nœuds dont la position a changé lors de l'étape de simulation précédente. Cette fonction traverse l'arbre binaire de bas en haut de manière sélective, en ne traversant que les noeuds dits actifs, comme cela sera expliqué plus bas.

Les OBB mises à jour forment la base de travail utilisée par la fonction Upd_Int_L() 32. Cette fonction, dont une implémentation est montrée sur la Figure 4 parcourt récursivement l'arbre binaire de haut en bas, depuis la racine RN.

La fonction Upd_Int_L() maintient pour chaque nœud deux listes principales Int_L et Int L Old. A chaque nœud, la liste Int L (respectivement Int L Old) contient des paires de nœuds appartenant au sous-ensemble d'objets désignés par le nœud donné qui interagissent par des forces locales dans la boucle en cours (respectivement la boucle précédente).

Ces listes permettent de détecter les interactions locales qui ont changé, c'est-à-dire les interactions locales qui ont été modifiées par les positions des objets à la suite de la boucle précédente.

La fonction Upd_Int() part d'un test d'activité 3200 du nœud courant. Ce test est lié à la métrique d'accélération du nœud courant (ou à une autre métrique appropriée le cas échéant), et peut être réalisé sous la forme d'un calcul de cette métrique, ou sous la forme de la lecture d'une donnée (communément appelée "flag") qui désigne le caractère actif ou non du noeud, décidé en fonction de la valeur de cette métrique.

Dans l'exemple décrit ici, un nœud est considéré comme actif si la valeur de sa métrique d'accélération est supérieure à un seuil choisi par l'utilisateur, c'est-à-dire s'il est considéré comme non-rigidifïé. Si le nœud est rigide, cela signifie que les accélérations globales au sein du corps articulé associé au nœud sont considérées comme nulle dans le cadre de la simulation.

Comme on l'a vu plus haut, cela signifie que les positions des objets au sein de ce corps n'ont pas changé, puisque les accélérations de ce corps ont été considérées comme nulles. Ainsi, les listes interactions IntJL et Int Old L de la boucle précédente restent valables pour ce nœud et tous ces nœuds enfants, et la fonction s'arrête (3202).

Si le nœud C est actif, la liste Int L(C) de la boucle précédente devient la liste Int Old L(C) pour la boucle courante par écrasement (3204), et la liste M L(C) est remise à zéro (3206).

Ensuite, une boucle est entamée pour mettre en queue toutes les interactions possibles entre les enfants du nœud C. Ainsi, les deux nœuds-fils immédiats A et B sont récupérés par une fonction CH() en 3208, et leurs OBB respectives OBBA et OBBB sont récupérées en 3210.

Le couple (OBBA, OBBB) est alors empilé dans une liste Q (3212), qui sert de pile de stockage des interactions potentielles entre OBB. La boucle mentionnée a pour condition d'arrêt le vidage de la pile Q (3214), et la fonction Upd_Int_L(C) s'arrête lorsque cette condition est remplie (3216).

La boucle commence par un dépilage de la liste Q, qui permet de récupérer le dernier couple d'interaction potentielle (OBBE, OBBF) introduit dans Q (3218). Lors de la première itération de cette boucle, il s'agit évidemment du couple (OBBA, OBBB) de l'étape 3212.

En 3220, une fonction Upd_Pos() met à jour les positions dans le monde réel des OBB OBBE et OBBF, afin de permettre d'établir leur collision éventuelle. Cette mise à jour de position "réelle" est distincte de celle réalisée dans la fonction Upd_OBB() 30, qui met à jour la position relative de chaque OBB dans le référentiel de son nœud.

La distinction entre position réelle et position relative permet d'implémenter l'algorithme de manière plus efficace et plus intégrée. Il serait néanmoins possible de mettre à jour la position réelle des OBB directement dans la fonction Upd OBBQ 30.

En 3222, la fonction Coll() détermine si OBBE et OBBF sont en collision, c'est-à-dire si leurs boîtes orientées présentent une intersection. Si ce n'est pas le cas, alors il n'y a pas interaction locale, et la boucle est réinitialisée en 3214 pour tester l'interaction potentielle suivante.

Dans le cas où OBBE et OBBF, qui sont associées à des nœuds respectifs E et F, sont en collision, quatre possibilités existent :

a) E et F sont tous les deux des feuilles de l'arbre, ce qui signifie que ces objets interagissent entre eux, et que la fonction s'arrête ;

b) F est un objet articulé et E est une feuille, ce qui signifie qu'il faut explorer F pour déterminer les interactions de ses nœuds fils avec l'objet associé à E, puis déterminer les interactions des nœuds fils de F entre eux ;

c) E est un objet articulé et F est une feuille, ce qui signifie qu'il faut explorer E pour déterminer les interactions de ses nœuds fils avec l'objet associé à F, puis déterminer les interactions des nœuds fils de E entre eux ; et

d) E et F sont tous deux des objets articulés, ce qui signifie qu'il faut explorer E et F pour déterminer les interactions entre leurs nœuds fils respectifs, puis déterminer les interactions des nœuds fils de E et de F entre eux.

Pour cela, un test 3224 détermine si E est un nœud feuille, en appelant une fonction CH_OBB() avec OBBE. La fonction CH_0BB() est similaire à la fonction CH(), sauf qu'elle prend en entrée une OBB et pas un nœud.

Si E est un nœud feuille, la fonction DPl () est appelée en 3230, qui gère les cas a) et b). Un exemple d'implémentation de la fonction DP1() est représenté sur la figure 5. Dans cette fonction, un test 3232 détermine si F est un nœud feuille, en appelant la fonction CH_OBB() avec OBBF.

Si F est un nœud feuille, alors le cas a) s'applique, et la liste Int L(C) reçoit la paire (E, F) en 3234, la fonction DP1() s'arrête en 3236, et la boucle de Upd ιnt L(C) reprend en 3214 avec le test de la pile Q.

Si F est un objet articulé, alors FA et FB sont ses fils, et le cas b) s'applique. Pour déterminer les interactions des fils de F avec E, on récupère les OBB de FA et FB en 3238, et on empile alors deux interactions potentielles (OBBE, OBBFA) en 3240 et (OBBE, OBBFB) en 3242 dans la pile Q.

Ensuite, les interactions des fils de F entre eux sont déterminées en appelant Upd_Int_L() avec F en 3244. Cet appel de fonction permet de mettre à jour les listes Upd Old L(F) et Upd L(F). Enfin la fonction DP1() s'arrête en 3236, et la boucle de Upd ιnt L(C) reprend en 3214 avec le test de la pile Q.

Si E est un objet articulé, alors EA et EB sont ses fils. Ensuite une fonction DP2() est appelée en 3250, qui gère les cas c) et d). Un exemple d'implémentation de la fonction DP2() est représenté sur la figure 6.

Dans cette fonction, un test 3252 détermine si F est un nœud feuille, en appelant la fonction CH_OBB() avec OBBF.

Si F est un nœud feuille, alors on est dans le cas c). Pour déterminer les interactions des fils de E avec F, on récupère les OBB de EA et EB en 3254, et on empile alors deux interactions potentielles (OBBEA, OBBF) en 3256 et (OBBEB, OBBF) en 3258 dans la pile Q.

Ensuite, les interactions des fils de E entre eux sont déterminées en appelant Upd_Int_L() avec E en 3260. Cet appel de fonction permet de mettre à jour les listes Upd_01d_L(E) et Upd_L(E). Enfin, la fonction DP2() s'arrête en 3262, et la boucle de Upd_Int_L(C) reprend en 3214 avec le test de la pile Q.

Si F est un objet articulé, alors FA et FB sont ses fils, et le cas d) s'applique. Pour déterminer les interactions des fils de E avec ceux de F, on appelle une fonction Max_Vol() avec OBBE et OBBF en 3264.

La fonction Max_Vol() détermine laquelle de OBBE et OBBF a le volume le plus important, afin d'explorer les interactions des fils du nœud associé avec l'autre nœud. Cela permet d'optimiser les performances de Upd_Int_L().

En résultat, on obtient trois OBB, OBBMA, OBBMB et OBBMin. Ainsi, si OBBE est plus volumineuse que OBBF, alors OBBMA contient OBBEA, OBBMB contient OBBEB et OBBMin contient OBBF. Dans le cas inverse, OBBMA contient OBBFA, OBBMB contient OBBFB et OBBMin contient OBBE.

On empile alors deux interactions potentielles (OBBMA, OBBMin) en 3266 et (OBBMB, OBBMin) en 3268 dans la pile Q. Ensuite, les interactions des fils de E entre eux sont déterminées en appelant Upd_Int_L() avec E en 3270, et celles des fils de F entre eux sont déterminées en appelant Upd_Int_L() avec F en 3272.

Ces appels de fonction permettent de mettre à jour les listes Upd Old L(E), Upd L(E), Upd_01d_L(F) et Upd_L(F). Enfin, la fonction DP2() s'arrête en 3262, et la boucle de Upd ιnt L(C) reprend en 3214 avec le test de la pile Q.

Comme cela apparaît implicitement, la fonction Upd_Int() appliquée à un nœud C présente donc deux récursions imbriquées.

La première récursion concerne l'empilage et le dépilage de la pile Q, de manière à analyser tous les couples d'interaction potentielle entre les fils du nœud C, quelque soit leur niveau de profondeur par rapport à ce nœud.

La deuxième récursion concerne l'appel sélectif de la fonction Upd_Int() avec des nœuds fils de C, qui permet de mettre à jour les listes Upd_01d_L() et Upd_L() de tous les nœuds actifs.

La détermination des interactions locales dans l'arbre binaire est principalement basée sur les propriétés des OBB, comme on l'a vu plus haut. Cependant, bien que les OBB constituent un outil intéressant, d'autres outils peuvent être employés, comme des tables de prédiction d'interaction, ou d'autres outils adaptés que l'homme du métier saura envisager au vu de l'application concernée.

Lorsque les listes Upd_Old_L() et Upd_L() ont été mises à jour pour tous les nœuds actifs, les données d'interaction sont prêtes à être mises à jour avec la fonction Upd_Int_Val() 34.

Une implémentation de la fonction Upd_Int_Val() est montrée sur la Figure 7. Comme on le voit sur cette figure, cette fonction commence par un bloc récursif 700. En appelant la fonction Upd_Int_Val() avec le nœud RN, le bloc 700 réalise une descente dans l'arbre binaire.

Le bloc 700 commence en 702 avec un test d'activité du nœud courant C, comme décrit précédemment. Si le nœud courant est inactif, alors la fonction se termine en 704. Si le nœud courant est actif, la descente dans l'arbre est commencée en obtenant les enfants A et B de ce nœud par la fonction CH() en 706.

Ensuite un test d'activité est réalisé sur le fils gauche A de C, comme décrit précédemment en 708. Si le nœud A est actif, alors la descente à gauche est réalisée en 710, en appelant Upd_Int_Val() avec A.

Ensuite, un test d'activité est réalisé sur le fils droit B de C, comme décrit précédemment, en 712. Si A est inactif, alors le test 712 est réalisé directement, et il n'y a pas de descente à gauche pour le nœud C.

Si le nœud B est actif, alors la descente à droite est réalisée en 714, en appelant Upd_Int_Val() avec B. C'est alors la fin du bloc 700 de descente récursive, et la fonction se poursuit avec sa partie effective de traitement. Si B est inactif, alors il n'y a

pas de descente à droite pour le nœud C, le bloc 700 se termine, et la fonction se poursuit avec sa partie effective de traitement.

Une fois que le bloc 700 a été répété et que les récursions ont été réalisées, le traitement des données commence au niveau du nœud C en réinitialisant une liste Upd(C) et une pile Q en 720.

La liste Upd(C) sert à faire remonter de proche en proche des données d'interaction mises à jour, à partir des données mises à jour dans les fils de C lors des appels récursifs effectués grâce au bloc 700.

La liste Q sert de pile de stockage d'objets dont les données d'interactions doivent à mettre à jour au niveau du nœud courant. Ainsi, la liste Q reçoit des listes Upd(A), Upd(B), Int Old L(C), et Int L(C) en 722. On notera que, bien que Int Old L(C), et IntJL(C) contiennent des paires de nœuds et pas des nœuds, l'homme du métier saura reconnaître que l'opération 722 consiste à empiler dans Q uniquement les nœuds distincts les uns des autres.

Une boucle de remontée de données d'interaction commence alors en 724, avec le dépilage de la liste Q dans un nœud de travail R. Le nœud de travail R est ajouté à la liste Upd(C) en 726, pour assurer la remontée des données qui lui sont associées lors du dépilage des récursions du bloc 700, puis un test 728 est réalisé pour déterminer si R appartient au sous-ensemble d'objets désigné par A, ou celui désigné par B.

Si R appartient au sous-ensemble d'objets désigné par A, alors la donnée d'interaction F(C, R) est réinitialisée avec la valeur de F(A, R) en 730. Sinon, F(C, R) est réinitialisée avec la valeur de F(B, R) en 732. Enfin, un test sur Q en 734 indique si toutes les données d'interaction des objets de Q ont été réinitialisées ou pas.

Lorsque Q est vide, une boucle de mise à jour de données d'interaction commence en 736 avec le dépilage d'une liste Int_L_t(C). La liste Int_L_t(C) est une copie de travail

de la liste Int L(C), et sert uniquement pour la boucle de mise à jour. Le résultat du dépilage de la liste Int L t(C) est stocké dans une paire de nœuds de travail (P, T).

En 738, les données d'interaction F(C, P) sont mises à jour en ajoutant à la valeur de F(C, P) réinitialisée dans la boucle ci-dessus la valeur F(P, T) de l'interaction de T sur P. La valeur de F(P, T) est calculée, au moyen d'une fonction F_Int(), et donne l'intensité de l'interaction inter-atomique de P avec T, dans le repère de P. Cette même opération est réalisée en 740 avec T, en ajoutant à la valeur de F(C, T) réinitialisée dans la boucle ci-dessus la valeur F(T, P) de l'interaction de P sur T, dans le repère de T.

La fonction F_Int() est basée sur les formules bien connues d'interaction électrostatique et de van der Waals entre deux atomes / molécules. D'autres formules peuvent bien sûr être utilisées. Il est important de noter que F_Int() calcule les interactions dans le repère associé au premier nœud invoqué.

Ensuite, en 742 un test vérifie si Int L t(C) a été entièrement dépilée. Si c'est le cas, l'instance concernée de Upd_Int_Val() se termine en 744. Sinon, la boucle de mise à jour reprend en 736. Une fois que la liste Int L t(C) a été entièrement dépilée, les données d'interactions F(C, X), où X est un nœud feuille descendant de C, représentent l'ensemble des interactions subies par le nœud X par tous les autres nœuds feuilles qui sont des descendants de C.

Il apparaît clairement au vu de ce qui précède que la fonction Upd_Int_Val() a une approche du bas vers le haut, bien qu'elle soit invoquée avec le nœud racine. Il apparaît également que cette fonction fait remonter de proche en proche toutes les interactions subies par chaque nœud feuille d'un nœud C par les autres nœuds feuilles de C.

Enfin, cette fonction garde par ailleurs la trace des nœuds dont les données d'interaction ont été mises à jour au moyen de la liste Upd(). Ainsi, lorsque la fonction Upd_Int_Val() se termine, le nœud racine RN contient d'une part toutes les données d'interaction sur ses nœuds feuilles, mais également la liste des données d'interaction

qui ont été mises à jour dans la boucle courante, ce qui assure de n'accéder qu'aux nœuds qui ont effectivement évolué.

Une fois que la fonction Upd_Int_Val() est terminée, la fonction Upd_Int() prend fin. Comme représenté sur la Figure 2, la fonction Upd_Int() est suivie de la fonction Sprd_Int() dans la boucle de simulation. La fonction Sprd_Int() eεt appelée avec le nœud racine RN, comme sa liste Upd(RN) comprend toutes les données d'interaction qui ont été modifiées par la fonction Upd lntQ.

Un exemple d'implémentation de cette fonction est représenté sur la figure 8. La fonction Sprd_Int() est une boucle qui commence par un test 802 sur la liste Upd(C) du nœud en entrée. Si cette liste est vide, la fonction se termine en 804.

Lorsque la liste Upd(C) est non vide, un élément R est dépilé en 806, les données d'interaction F(R, R) du nœud correspondant sont mises à jour à la valeur F(C, R) en 808.

Ensuite, une fonction Set_Flg() enregistre des données Flg bph associées à R en 810. La valeur des données Flg_bph signale la nécessité de mettre à jour les coefficients dynamiques b, p et η. La fonction Set FlgO enregistre également des données Flg bph pour tous les ascendants du nœud R, jusqu'au noeud racine RN. La boucle reprend alors en 802, jusqu'à ce que la liste Upd(C) soit vide.

On notera que la fonction Sprd_Int() peut être implémentée de nombreuses manières. Notamment, les données d'interaction mises à jour sont ici stockées dans le nœud concerné. Cependant, il serait possible de stocker ces interactions à un autre endroit, dans une table séparée par exemple. Il en va de même pour les données enregistrées par la fonction Set_Flg().

II serait également possible de ne pas exécuter cette fonction, et de se contenter d'un appel à la liste Upd(RN) de manière adéquate lorsque cela est nécessaire, bien que cela

aurait pour conséquence un certain manque de clarté et introduirait une relative complexité de gestion.

Lorsque les données d'interaction ont été mises à jour, la boucle de simulation se prolonge avec la fonction Get_Dyn_Coeff(), qui permet de déterminer les paramètres dynamiques des objets et des nœuds de l'arbre.

Un exemple d'implémentation de la fonction Get_Dyn_Coeff() est représenté sur la figure 9. La fonction Get_Dyn_Coeff() est une fonction à propagation de bas en haut, comme la fonction Upd_Int_Val(). Pour obtenir cette propagation, un test 902 vérifie si le nœud C a des nœuds fils A et B.

Si le nœud C n'a pas de nœuds fils, un test 904 vérifie si la donnée Flg bph(C) est activée. Si c'est le cas, alors le nœud C est un nœud feuille pour lequel une fonction Coeff_bphl() est appelée en 906, fonction qui met p(C) et η(C) à 0, et qui calcule b(C) à partir des données d'interaction mises à jour F(C, C). Après l'étape 906, ou si la donnée Flg bph(C) n'est pas activée, la fonction se termine en 908.

Si le nœud C a des fils, la fonction Get_Coeff_Dyn() s'appelle récursivement avec A et B en 910 et 912. Cela assure un parcours des nœuds nécessaires de l'arbre binaire de bas en haut. Une fois que tous les appels récursifs ont été lancés et que tous les nœuds feuilles nécessaires ont été traités, une remontée de l'arbre est opérée, boucle de récursion par boucle de récursion, dans un bloc 920.

Dans le bloc 920, un premier test 922 vérifie si le nœud courant C est actif. Si c'est le cas, alors les coefficients φ(C) et ψ(C) sont calculés en 924 en appelant une fonction Coeff_phipsi() avec les valeurs de ces coefficients pour les fils de C, nommément A et B. Ensuite, les coefficients b(C), p(C) et η(C) sont calculés en 926 en appelant une fonction Coeff_bph(), avec les valeurs de ces coefficients pour les fils de C, nommément A et B. Enfin, la fonction se termine en 928.

Dans le cas où le nœud courant C est inactif, un test est réalisé sur la donnée Flg bph(C) en 930. Dans le cas où le test est négatif, la fonction se termine directement en 928. Si ce test est positif, c'est-à-dire si la donnée indique qu'il est nécessaire de mettre à jour les coefficients b(C), p(C) et η(C), alors ceux-ci sont calculés en 932 en appelant la fonction Coeff_bph() avec les valeurs de ces coefficients pour les fils de C, nommément A et B. Enfin, la fonction se termine en 928.

Les formules (l) à (13) de l'Annexe A donnent la formule du calcul sur la base duquel ces coefficients sont établis dans les fonctions Coeff_phipsi() et Coeff_bph(). Pour les besoins de la présente description, il suffit de savoir que ces coefficients déterminés récursivement permettent de décrire localement le comportement des corps articulés, et servent de base aux calculs de métrique, d'accélération et de position de la simulation. Pour plus d'informations, les documents SPM '05 et SIGGRAPH '05 explicitent la définition précise et l'obtention de ces équations. La demanderesse recommande la consultation de ces documents à cet effet.

On notera en particulier que le calcul des coefficients b(C), p(C) et η(C) peut faire intervenir directement et indirectement une grandeur Q qui représente une donnée de force active sur l'articulation d'un nœud C. Cette grandeur Q peut être utilisée, par exemple, pour représenter un couple dépendant d'un angle dièdre. Dans la mesure où le calcul de la grandeur Q associée au nœud C ne dépend que de la position relative des fils A et B de C, cette grandeur est mise à jour uniquement pour les nœuds actifs. Ainsi, cette mise à jour peut se faire, par exemple, immédiatement avant le calcul des b(C), p(C) et η(C) par la fonction Coeff_bph().

Lorsque les coefficients dynamiques ont été mis à jour, la fonction Upd_Dyn() met à jour les données dynamiques sur la base de ces coefficients. Un exemple d'implémentation de la fonction Upd_Dyn() est représenté sur la Figure 10.

La fonction Upd_Dyn() représentée sur la Figure 10 est une fonction récursive qui traverse l'arbre binaire de haut en bas.

Cette fonction utilise une queue de priorité visant à optimiser la détermination des nœuds actifs, comme cela apparaîtra plus bas. Pour cela, un test 1002 sur le nœud courant appelle un calcul de métrique d'accélération 1003 si le nœud courant est le nœud racine RN. Le résultat sert de base en 1004 à une borne ε qui sert de condition d'arrêt en 1006.

Le calcul de métrique d'accélération est réalisé par une fonction Tot_Acc() qui réalise le calcul décrit avec la formule 21 de l'Annexe A, et qui correspond à la somme des carrés des accélérations des nœuds fils du nœud courant. Cela permet de définir les nœuds actifs et inactifs.

En effet, εmax définit la précision de simulation recherchée en termes de métrique d'accélération. Donc, si ε < εmax, cela signifie que la précision voulue est atteinte. L'accélération du nœud C qdd(C) n'est alors pas calculée, et la fonction se termine en 1010.

En revanche, si ε > εmax, la simulation n'est pas considérée comme encore assez précise. En 1012, un test détermine si le nœud C a des fils ou pas. Si ce n'est pas le cas, la fonction se termine en 1010.

Une fonction Active Flg(C) est alors appelée en 1014. La fonction Active_Flg() enregistre une donnée qui contient un marqueur de temps qui est associé à la boucle de simulation en cours. C'est ce marqueur de temps qui est utilisé dans les tests d'activité d'un nœud. Pour cela, pour le nœud donné, on compare l'éventuel marqueur de temps qui lui est associé au marqueur de temps de la boucle courante. Lorsqu'il y a concomitance, alors le nœud est considéré actif. Une méthode semblable peut être utilisée pour activer les données Flg bph(C).

Si le nœud C a des enfants A et B, l'accélération du nœud C qdd(C) est calculée en 1016 au moyen de la fonction Q_Dot. La fonction Q_Dot() utilise les coefficients dynamiques φ(C) et b(C), ainsi qu'une force F(C). Le calcul réalisé par la fonction Q_Dot() est explicité par l'équation 22 de l'Annexe A.

La force F(C) est un vecteur composite qui représente l'ensemble des forces d'articulation qui s'appliquent sur le nœud C. Pour le nœud racine RN des conditions particulières permettent de déterminer le vecteur F(C), comme décrit dans SPM' 05 section 3.1.2.

Comme décrit dans l'article SPM '05 section 3.1.1, un vecteur F(C) a deux composantes : F(I, C) et F(2, C). Dans le cadre de l'invention, le vecteur F(C) d'un nœud sert à calculer le vecteur de ses nœuds fils F(A) et F(B).

En effet, si on considère les deux fils A et B de C, alors la force F(I, C) est égale à F(I, A), et la force F(2, C) est égale à F(2, B). En outre, les composantes F(2, A) et F(I, B) peuvent être obtenues également à partir de F(C), comme décrit plus bas. Par conséquent les composantes de F(A) et de F(B) peuvent être obtenues totalement à partir de F(C).

Une fois l'accélération du nœud C calculée, la fonction Upd_Dyn() se poursuit en 1018 et en 1020 avec le calcul des vecteurs F(2, A) et F(I, B) au moyen de la fonction Ch_F() selon la formule 23 de l'Annexe A.

Ensuite, les métriques d'accélération des nœuds A et B sont calculées en 1022 et 1024, et les nœuds A et B sont rentrés en 1026 dans une queue de priorité Q_prio en fonction des valeurs de Acc(A) et Acc(B). En 1028, la valeur du carré de l'accélération qdd(C) est retirée à ε pour mettre à jour la précision de la simulation. Enfin, un nœud R est dépilé de la queue Q_prio en 1030, et la fonction Upd_Dyn() est appelée avec R en 1032 pour assurer la récursion de la fonction et la descente dans l'arbre.

La queue Q_prio classe les nœuds en fonction de la valeur de leur métrique d'accélération, afin de déterminer rapidement les noeuds les plus importants. Ainsi, comme la métrique d'accélération d'un nœud représente la somme des carrés des accélérations de ses fils, la fonction descend progressivement dans les nœuds qui feront baisser ε le plus rapidement, ce qui accélère la simulation.

On notera que cette fonction réalise deux tâches simultanément. En effet, l'activité d'un nœud est basée sur la valeur de la métrique d'accélération. Le non calcul des accélérations de certains nœuds, et l'utilisation des données Active_Flg() permet donc également de les écarter des boucles suivantes qui prennent en compte l'activité, sans coût de calcul supplémentaire.

Cela est notamment rendu possible par le fait que les fonctions qui parcourent l'arbre de bas en haut commencent par descendre récursivement dans les nœuds actifs uniquement.

Par conséquent, l'ensemble des fonctions ci-dessus, par sa sélection des nœuds actifs, assure un parcours minimal de l'arbre, et constitue une mise à jour sélective des données de certains des nœuds.

On notera par ailleurs que la précision de la simulation est ici assurée au moyen du paramètre εmax qui représente l'erreur maximale d'accélération. D'autres moyens peuvent être employés, comme la spécification du nombre de nœuds à simuler, ou d'autres métriques.

Lorsque la mise à jour des données dynamiques a été réalisée, la fonction Upd_Geo() termine la boucle de simulation. En effet, comme on l'a mentionné plus haut, les données d'interaction dans les nœuds sont calculées par rapport à un repère local particulier à chaque nœud.

De fait, les accélérations des nœuds et les positions que l'on peut dériver sont également relatives, ou intrinsèques. Pour pouvoir afficher le résultat de la simulation et pour pouvoir faire remonter (ou redescendre) les données dans les diverses traversées de l'arbre, des matrices de passage que l'on n'a pas discuté jusqu'à maintenant sont donc nécessaires. Ce sont ces matrices de passage qui sont utilisées pour calculer la position réelle des OBB dans la fonction Upd_Pos() de la figure 4, ou encore pour mettre à jour la force F(C, R) dans les étapes 732 et 734 de la figure 7.

Pour mettre à jour ces matrices, la fonction Upd_Geo() commence en 1102 avec une sous-fonction Upd_Hyb_qdd(). La fonction Upd_Hyb_qdd() vient calculer des coefficients dynamiques corrigés pour tenir compte de la rigidification réalisée avec la fonction Upd_Dyn(). Ensuite, les accélérations qdd(C) des nœuds sont recalculées avec les coefficients dynamiques corrigés. Les opérations réalisées par cette fonction, ainsi que leur intégration dans le calcul des qdd(C) sont décrites dans l'article SIGGRAPH '05, sections 4 et 5.

Bien qu'elle soit décrite dans le mode de réalisation préféré de l'invention, la fonction Upd_Hyb_qdd()_n'est pas forcément nécessaire dans tous les cas. En option, la fonction Upd_Geo() peut ne pas contenir cette fonction et procéder directement en 1104.

La fonction Upd_Geo() se poursuit en 1104 avec une fonction Upd_qpos(), qui vient mettre à jour les positions intrinsèques des nœuds sur la base de leurs accélérations. Ce calcul est simple et représente une intégration temporelle basique correspondant à l'équation 14 de l'Annexe A. D'autres méthodes peuvent être utilisées pour la mise à jour des positions intrinsèques des nœuds, faisant intervenir par exemple des calculs d'accélérations supplémentaires.

Enfin, la fonction Upd_Geo() se termine en 1106 avec la fonction Upd_CT_Mat(), qui met à jour les matrices internes mentionnées plus haut sur la base des positions intrinsèques mises à jour des nœuds. Ces matrices ainsi que leur mise à jour sont décrites plus avant dans la section 5 de SPM '05. Dans ce qui précède, les données d'interaction et de géométrie sont principalement calculées de manière intrinsèque, c'est-à-dire par rapport à un repère local associé à chaque fois au nœud concerné.

Cela permet de ne pas avoir à recalculer des données qui dépendent seulement des positions relatives des objets, comme les coefficients φ(C) et ψ(C) et les forces d'interaction locale. Ainsi, le passage aux coordonnées du monde réel n'est nécessaire que pour le calcul des positions réelles des nœuds.

Dans d'autres modes de réalisation, il serait possible d'opérer à chaque fois sur des données de coordonnées du monde réel, en ne conservant des données intrinsèques que pour les nœuds qui ne sont pas actifs.

Bien que la présente invention ait été décrite en rapport à la représentation au moyen d'un arbre binaire, il serait possible d'utiliser d'autres types d'arbre, par exemple ternaire ou même n-aire. Il serait également possible d'utiliser d'autres types de représentations, dès lors qu'elles permettent de différencier des nœuds actifs et des nœuds inactifs avec un coût de parcours comparable.

Les fonctions ci-dessus présentent les étapes principales de mise en œuvre de l'invention. Elles peuvent être utilisées par l'homme du métier comme base pour rédiger un programme informatique implémentant l'invention. Un tel programme peut être rédigé en langage C++ ou en tout autre langage informatique approprié tel que saura le reconnaître l'homme du métier.

ANNEXE A

SECTION 1

avec où S est le sous-espac de mouvement du noeud, de dimension 6 x nombre de degrés de liberté du noeud. Typiquement, S = (0, 0, 1, 0, 0, l) ^τ .

avec

Nota : Les matrices φ(C) et φ(C) sont des matrices de matrices.

ANNEXE A (suite)

SECTION 2

avec où Q est le vecteur de dimension nombre de degrés de liberté x 1 des forces actives du noeud.

avec

ANNEXE A (suite)

SECTION 3