CN101779203B - 对相互作用对象集合进行仿真的计算机设备及相应方法 - Google Patents

Abstract

一种对相互作用对象集合进行仿真的计算机设备及相应方法。一种对相互作用对象集合进行仿真的计算机设备，包括：保存有对象树状表示的存储器(8)，其中，每个结点与动力学数据、几何数据、相互作用数据相关联，所述动力学数据、几何数据、相互作用数据取决于子结点数据，并且对于某些结点而言，还取决于局部相互作用数据；仿真控制器(4)，用于重复地启动：相互作用数据分配器(10)，相互作用更新机构，所述相互作用更新机构遍历树状表示，并根据子结点的局部相互作用数据更新结点的相互作用数据，动力学数据更新机构(12)，所述动力学数据更新机构(12)遍历树状表示，并根据有关几何、相互作用数据运行，几何数据更新机构(14)，所述几何数据更新机构(14)为倾向于相互作用的结点遍历树状表示，并根据已更新的动力学数据运行。本发明还涉及一种对相互作用对象集合的进行仿真方法，和一种计算机程序产品。

Description

对相互作用对象集合进行仿真的计算机设备及相应方法

技术领域

本发明涉及对相互作用对象集合的仿真。

背景技术

相互作用对象仿真涉及从多机械体(复杂的或简单的)仿真到分子(例如蛋白质分子)行为仿真的广泛领域。

现有的、最早的方案采用对象归类的原则，即把对象分组汇集以便定义总体行为。这类方案因其较低的计算成本而有意义。

但是，计算成本的降低是以损害仿真的现实性为代价的，因为与分组对象间关系相关的许多方面被有意忽略或简化。

其它方案考虑了相互作用对象，例如关节联结体。关节联结体通常由以关节相互联结的两个体构成。所述两个体自身可以是关节联结体，以此类推。这种方案容许对这类对象集合的机械行为进行逼真建模，同时还保证高性能。

因为对象主要是通过联结它们的关节相互作用的，所以称“机械”行为。为了对对象集合建模，只需要对关节建模。

当对象间的相互作用变得更加复杂时，这些方案很快失效，因为相关计算需求/成本如此高昂而无法实现。

发明内容

本发明的目的是改善这一状况。

为此目的，本发明提出一种用于对相互作用对象集合进行仿真的计算机设备，所述设备包括：存储器，该存储器包含该对象集合的树状表示。在该树状表示中，每个中间结点与动力学数据、几何数据和相互作用数据相关联。

所述设备还包括一个仿真控制器，该仿真控制器根据反复进行的循环启动：

-对象、对象子集相互作用数据的分配器；

-根据有关相互作用数据和几何数据，更新动力学数据的机构，为倾向于相互作用的结点遍历树状表示；

-根据已更新的动力学数据，更新几何数据的机构，为倾向于相互作用的结点遍历树状表示。

为了考虑不依赖于关节的结点间的相互作用，所述存储器还包含与至少一些结点相关联的局部相互作用数据，并且，所述相互作用数据分配器包括为倾向于相互作用的结点遍历树状表示，更新相互作用的机构，以便根据有关子结点的局部相互作用数据来更新结点的相互作用数据。

这种设备优点突出，因为它提供了延伸关节联结体建模方案的可能性，同时考虑了非关节性力，并且没有显著增加计算时间，从而有可能提高仿真的可应用性。

本发明还涉及一种对相互作用对象的行为进行仿真的方法，其中使对象集合保持树状表示。在该树状表示中，每个中间结点与动力学数据、几何数据和相互作用数据相关联，所述动力学数据、几何数据和相互作用数据取决于其子结点数据。

该方法重复实施一个循环，所述循环包括如下步骤：

a.分配一些对象和对象子集的相互作用数据；

b.为倾向于相互作用的结点遍历树状表示，根据有关相互作用数据和几何数据，更新这些结点的动力学数据；和

c.为倾向于相互作用的结点重新遍历树状表示，根据在步骤b中更新的动力学数据，更新几何数据，

更具体地，步骤a包括：

a1.为倾向于相互作用的结点预遍历树状表示，根据有关子结点的局部相互作用数据，更新这些结点的相互作用数据。

附图说明

阅读以下描述之后，可以更好地了解本发明的其它特征和优点。以下描述针对结合如下附图的说明性、非限制性的实例，其中：

-图1是根据本发明的计算机设备的示意图；

-图2是图1所示设备实施的仿真循环的函数执行的方框图；

-图3是图2所示仿真循环的第一个函数的方框图；

-图4是图3所示的一个函数的方框图；

-图5和图6是图4所示的两个函数的示意图；

-图7是图3所示的另一个函数的方框图；

-图8是图2所示仿真循环的第二个函数的方框图；

-图9是图2所示仿真循环的第三个函数的方框图；

-图10是图2所示仿真循环的第四个函数的方框图；

-图11是图2所示仿真循环的第五个函数的方框图；

附图和下文的描述基本上包含了确切特征的元素，因此不仅可以帮助理解本发明，而且在必要时有助于其定义。

具体实施方式

本描述旨在阐明能由著作权和/版权保护的元素。权利人不反对任何人复制包含在官方文件中的本专利文档或其描述。权利人完整保留其余权利。

本发明还包括所描述的、作为产品的软件。所述软件保存在任何可由计算机读取的“介质”(载体)上。术语“可由计算机读取的介质”包括磁的、光的和/或电子的数据存储载体，还包括传输载体或工具，例如模拟或数字信号。

另外，作为对详细描述的补充，附录A给出了本发明中使用的数学公式。该附录被单独列出，目的是为了表述清晰，便于提交。该附录是本说明书的组成部分，因此，不仅帮助理解本发明，而且在必要时有助于其定义。

第一个对关节联结体进行机械仿真的方案在如下文章中提出：Featherstone著的《用于刚体动力学并行o(log(n))计算的关节联结体分治算法，第一部分：基本算法》，国际机器人技术研究杂志18(9)：867-875(”A divide-and-conquer articulated bodyalgorithm forparallel o(log(n))calculation of rigid body dynamics.part 1：BasicAlgorithm”，International Journal of Robotics Research 18(9)：867-875)，下文中记作Featherstone‘99。

Featherstone‘99描述了一个分割算法，将关节联结体分割为二叉树，以便递归地计算带力和位置分量的系统。对于关节联结体，系统包括连接两个体的关节，而这两个体本身可以是关节联结体。

因此，关节联结体的概念不应该局限于其简单的机械含义。还应该从广义上理解关节的概念。例如关节可以包括数量可变的自由度，介于零(两个子系统的刚性耦合)与六(两个子系统相对运动不受限制)之间。

二叉树包括叶结点，这些叶结点表示关节联结对象的最基本对象，并且，每个叶结点都与一个关节(或结点)联结。因此，二叉树的每个结点本身是一个关节联结体，表示一个对象子集。

Redon和Lin所著的文章《用于复杂连接准静态仿真的高效的误差约束近似算法》(美国计算机学会固体与物理建模研讨会汇编，2005年。下文中记作SPM‘05)(”Anefficient，Error-boundedApproximation Algorithm for Simulating Quasi-Staticsof ComplexLinkages”，Proceedings of ACM Symposium on Solid andPhysicalModeling，2005)定义了关节的加速度度量，从而完善了Featherstone‘99的方法。

利用所述度量，可以只根据给定关节的数据评估与该关节关联的关节联结体的加速度，而不需知道叶结点数据。因此可以将树“刚化”：对于加速度度量小于选定阈值的关节联结体，取消其加速度。因此，树的“刚化”容许将树的遍历限制在少量结点中，同时还能控制仿真的逼真度。

在Redon、Gallopo、Lin所著的文章《关节联结体自适应动力学》(美国计算机学会图形学汇刊(SIGGRAPH 2005)，下文中记作SIGGRAPH‘05)(“Adaptative Dynamics ofArticulated Bodies”，ACMTransactions on Graphics(SIGGRAPH 2005))中，仿真的质量得到改进，所采取的手段是还处理动态情况，而不仅仅是是准静态情况，并且增加了考虑刚化的方法，以避免产生没有具体现实性的力。

所描述的实施方式利用了所有这些方法的优点，有时可以作为这些方法的延伸而实施，当然，还可以考虑本发明的大量其它实施方式。

为了仿真，每个树结点都与各种数据相关联。可以区分出三类主要数据，这三类数据将在下文中首先介绍：相互作用数据、动力学数据、几何数据。

相互作用数据可以如此定义：表示施加于某给定结点或该结点表示的对象子集中的某给定对象的力。这些力可以分为两类：

-应力，涉及对象间通过关节的相互作用；

-局部力，包括可能施加于某给定结点的所有其它力。

在本文所描述的实例中，对象是分子和原子，局部力包括：分别施加于关节的外力；所有原子间力和分子间力，例如静电力、范德华(van der waals)力，等等。

动力学数据表示关节联结体动力学的特征复合数据。这些数据容许计算关节的加速度和速度度量，而不必知道全部对象的全部位置数据和力数据。下文中将进一步介绍动力学数据的应用。

几何数据表示关节联结体的只依赖于对象间相对位置的一组特征数据。这些数据尤其容许在局部参照范围内实施大部分计算，从而简化了仿真。这些数据将在下文中阐述。

此处所述例子应用于对象的准静态仿真，也就是说，认为对象的速度在每一时刻都为零。

图1是根据本发明的设备的示意图。

设备2包括一个仿真控制器4、一个仿真模块6、一个存储器8。仿真模块6包括一个相互作用数据分配器10、一个动力学数据更新机构12、一个几何数据更新机构14。

在仿真时，基于存储器8中保存的数据，控制器4根据一个仿真循环(在下文中描述)重复地启动仿真模块6。

设备2可以以多种形式实施。当前的优先实施方式是一个计算机程序，该计算机程序以在计算机上可执行的形式实现控制器4和模块6，所述计算机可以在市面上买到，例如带如下配置：一个主频1.7GHz的微处理器、1GB的RAM存储器、一个容量足够装操作系统和所述程序代码的硬盘。

还可以采用其它实施形式，例如采用专用印刷电路(子卡、ASIC，等等)或专用处理器(FPGA、微芯片，等等)。

图2是模块6实施的仿真循环的示意图。在该循环中，函数Upd_Int()20，Sprd_Int()22，Get_Dyn_Coeff()24，Upd_Dyn()26，Upd_Geo()28相继被调用。

图2所示循环表示给定运行循环中的仿真的主要步骤。这些步骤分别为：更新对象间非关节性相互作用、传播这些相互作用、获得并更新仿真的动力学参数、以及更新仿真的几何参数。

图2所示各步骤对应于实施本发明的程序中的函数，但如下文所示，这些函数本身又包括多个子函数，并且应该被作为功能成分考虑，专业人士完全可以按照自己的意愿实施这些功能成分。

从功能角度来看，可以认为：

◆函数Upd_Int()20和Sprd_Int()22构成数据分配器10的核心；

◆函数Get_Dyn_Coeff()24和Upd_Dyn()26构成动力学数据更新机构12的核心；

函数Upd_Geo()28构成几何数据更新机构14的核心。

图3表示图2所示函数Upd_Int()20的可能实施方式。如图3所示，函数Upd_Int()20包括3个子函数：Upd_BB()30，Upd_Int_L()32，Upd_Int_Val()34。

根据传统方案，一旦一个对象子集被刚化，就暂时停止在仿真中处理该对象子集。因为这类仿真中考虑的力是机械力，所以所述对象子集可以在空间上非常接近另一个对象子集而不产生任何后果。

当希望对蛋白质或分子链的行为进行仿真时，就向应力或外力中加入静电类型的或范德华的原子间类型的局部力。

因此，“机械”仿真不再足够，因为可能有一个刚化的子集与另一个子集相互作用。于是需要对其“解除刚化”，以便考虑这些新的局部力，这些新的局部力只取决于对象间的相对位置。而根据现有方案，无法考虑这些新的局部力，除非放弃现有方案的所有优点。

函数Upd_Int()20容许考虑这些局部力。为此，该函数从仿真给定分子状况开始，容许首先确定哪些是可能相互作用的对象。然后，该函数计算相应局部力，并有序地归纳(intègrer)这些数据。

函数Upd_BB()30是对相互作用检测加以调整以适用于分子仿真。因为原子间力主要是短距离的，所以该函数定义了每个对象周围的潜在相互作用空间。

更精确地说，这些空间被定义为围绕一个对象或一组对象的“有向盒”。所述盒本身是在仿真开始时为二叉树的每个结点定义的。运行函数Upd_BB()30就为每个盒更新参考点位置和取向。

有向盒(在文献中称作“Oriented Bounding Box(有向包围盒)”或OBB)的原理在如下文章中有描述：Gottschalk等人所著的《OBB树：用于快速干涉检测的层次结构》(美国计算机学会图形学汇刊(SIGGRAPH 1996))(“Obbtree：a hierarchical structure forrapidinterference detection”，ACM Translations on Graphics(SIGGRAPH1996))。

OBB只依赖于所围绕的对象的位置，并且在与结点(该结点是与所述OBB相关联的结点)相关的参照系中表达，所以OBB在结点刚化后是不变的。因此，函数Upd_OBB()30只对在前一仿真步骤中位置改变了的结点起作用。该函数从下向上选择性地经历(traverser)二叉树，只经历活动结点，如下文所述。

更新后的OBB形成了函数Upd_Int()32使用的工作基础。该函数(图4显示其实施)自根RN开始，从上向下递归地遍历二叉树。

函数Upd_Int_L()为每个结点维护两个主要清单Int_L和Int_L_Old。在每个结点，清单Int_L(对应地，Int_L_Old)包含一些结点对，这些结点对属于某给定结点确定的对象子集，该对象子集中的对象在当前循环(或前一循环)中通过局部力相互作用。

这些清单容许检测已经变化了的局部相互作用，即在前一循环后，因对象位置而改变了的局部相互作用。

函数Upd_Int()首先检测当前结点的活动性3200。该检测与该当前结点的加速度度量相关(或在需要时与另一适当度量相关)，其实施的形式是计算该度量，或读取一个指示结点活动性与否的特征的数据(通常称为“标记(flag)”)。

在此处所描述的实例中，如果结点的加速度度量值大于用户选定阈值，即如果认为该结点未刚化，那么就认为该结点是活动的。如果结点是刚性的，则意味着与该结点相关联的关节联结体内的整体加速度在仿真中被认为是零。

正如上文所述，这种情况意味着体内各对象的位置没有改变，因为该体的加速度被认为是零了。因此，前一循环的相互作用清单Int_L和Int_Old_L对该结点及其所有子结点仍然有效，函数结束(3202)。

如果结点C是活动的，那么前一循环的清单Iht_L(C)通过覆盖变成当前循环的清单Int_Old_L(C)(3204)，而清单Int_L(C)置零(3206)。

然后，启动一个循环将结点C的子结点间所有可能的相互作用排成队列。因此，在3208，函数CH()取回(recuprer)两个直接子结点A和B；在3210，取回A和B各自的OBB即OBBA和OBBB。

于是一对(OBBA，OBBB)入栈清单Q(3212)，该清单Q用作保存OBB间潜在相互作用的栈。所提到的循环的退出条件为栈Q空(3214)：当该条件满足时，函数Upd_Int_L(C)结束(3216)。

所述循环首先退出清单栈Q，从而可以取回Q中最后一对潜在相互作用(OBBE，OBBF)(3218)。在该循环第一次迭代时，显然取回的是步骤3212的(OBBA，OBBB)对。

在3220，函数Upd_Pos()更新各OBB(即OBBE和OBBF)在现实世界中的位置，以便可以确定二者之间可能的碰撞。这种“现实”位置的更新与函数Upd_OBB()30中实施的更新是分开进行的，函数Upd_OBB()30在OBB的结点的参照系中更新每个OBB的相对位置。

现实位置与相对位置之间的分离可以更有效地、更完整地实施算法。当然，仍然可以直接在函数Upd_OBB()30中更新OBB的现实位置。

在3222，函数Coll()确定OBBE和OBBF是否碰撞，也就是说其有向盒是否相交。如果不相交，则没有局部相互作用，循环从3214重新开始，以便检测下一个潜在相互作用。

如果OBBE和OBBF(分别与结点E和F相关联)碰撞，那么存在四种可能性：

a)E和F都是树的叶，意味着这两个对象之间相互作用，函数结束；

b)F是关节联结对象而E是叶，意味着需要搜索F以便确定其子结点和与E关联的对象之间的相互作用，然后确定F的子结点之间的相互作用；

c)E是关节联结对象而F是叶，意味着需要搜索E以便确定其子结点和与F关联的对象之间的相互作用，然后确定E的子结点之间的相互作用；

d)E和F都是关节联结对象，意味着需要搜索E和F以便确定其各自子结点间的相互作用，然后确定E和F的子结点之间的相互作用。

为此，在3224，关于OBBE调用函数CH_OBB()，以便检测E是否是叶结点。函数CH_OBB()与函数CH()相似，唯一的区别在于CH_OBB()用OBB，而不是用结点，作为输入。

如果E是叶结点，则在3230调用函数DP1()，处理情况a)和b)。图5描述了函数DP1()的一个实施实例。在函数DP1()中，在3232，关于OBBF调用CH_OBB()，以便检测F是否是叶结点。

如果F是叶结点，那么适用情况a)，在3234，清单Int_L(C)接收(E，F)对。在3236，函数DP1()结束，在3214，循环Upd_Int_L(C)重新开始，检测栈Q。

如果F是关节联结对象，那么FA和FB是其子，适应情况b)。为了确定F的子与E的相互作用，在3238取回FA和FB的OBB，在3240和3242分别将两个潜在相互作用(OBBE，OBBFA)和(OBBE，OBBFB)入栈Q。

然后，在3244，关于F调用Upd_Int_L()，以便确定F的各子之间的相互作用。该函数调用容许更新清单Upd_Old_L(F)和Upd_L(F)。最后，在3236，函数DP1()结束。在3214，循环Upd_Int_L(C)重新开始，检测栈Q。

如果E是关节联结对象，那么EA和EB是其子。然后，在3250，调用函数DP2()，该函数DP2()处理情况c)和d)。图6描述了函数DP2()的一个实施例。

在函数DP2()中，在检测步骤3252，关于OBBF调用函数CH_OBB()，以便确定F是否是叶结点。

如果F是叶结点，那么适用情况c)。为了确定E的子与F的相互作用，在3254取回EA和EB的OBB，在3256和3258分别将两个潜在相互作用(OBBEA，OBBF)和(OBBEB，OBBF)入栈Q。

然后，在3260，关于E调用Upd_Int_L()，以便确定E的各子之间的相互作用。该函数调用容许更新清单Upd_Old_L(E)和Upd_L(E)。最后，在3262，函数DP2()结束。在3214，Upd_Int_L(C)重新开始，检测栈Q。

如果F是关节联结对象，那么FA和FB是其子，适用情况d)。为了确定E的子与F的子的相互作用，在3264，关于OBBE和OBBF调用函数Max_Vol()。

函数Max_Vol()确定OBBE和OBBF中哪个的体积最大，以便搜索其相关联的结点的子与另一个结点的相互作用。这一办法可以优化Upd_Int_L()的性能。

结果，得到三个OBB：OBBMA、OBBMB、OBBMin。因此，如果OBBE比OBBF体积大，那么OBBMA包含OBBEA，OBBMB包含OBBEB，OBBMin包含OBBF。在相反的情况下，OBBMA包含OBBFA，OBBMB包含OBBFB，OBBMin包含OBBE。

因此，在3266和3268，分别将两个潜在相互作用(OBBMA，OBBMin)和(OBBMB，OBBMin)入栈Q。然后，在3270，通过关于E调用Upd_Int_L()，确定E的各子之间的相互作用。在3272，通过关于F调用Upd_Int_L()，确定F的各子之间的相互作用。

所述程序调用容许更新清单Upd_Old_L(E)、Upd_L(E)、Upd_Old_L(F)、Upd_L(F)。最后，在3262，函数DP2()结束。在3214，循环Upd_Int_L(C)重新开始，检测栈Q。

从上述描述中可以看出，应用于结点C的函数Upd_Int()带有两个嵌套递归。

第一个递归涉及栈Q的入栈和出栈，以便分析结点C的各子结点(无论其相对于C的深度如何)之间的潜在相互作用力矩(couple)。

第二个递归涉及关于C的子结点选择性地调用函数Upd_Int()，从而可以更新所有活动结点的清单Upd_Old_L()和Upd_L()。

如上文所示，在二叉树中确定局部相互作用主要是基于OBB的属性。不过，尽管OBB是有益的工具，仍然还可以采用其它工具，例如相互作用预测表，或专业人士可能就有关应用设计的其它适当工具。

当清单Upd_Old_L()和Upd_L()针对所有活动结点都更新以后，就可以用函数Upd_Int_Val()34更新相互作用数据了。

图7显示了函数Upd_Int_Val()的实施。如该图所示，该函数以一个递归块700开始。通过关于结点RN调用函数Upd_Int_Val()，块700在二叉树中下行。

块700以702开始，测试当前结点C的活动性，如前所述。如果当前结点是不活动的，则函数于704结束。如果当前结点是活动的，则在706通过函数CH()获得该结点的子女A和B，从而开始在树中的下行。

然后，在708，对C的左子结点A进行活动性检测，如前所述。如果左子结点A是活动的，则在710通过关于A调用Upd_Int_Val()而进行左边下行。

然后，在712，对C的右子结点B进行活动性检测，如前所述。如果A是不活动的，则直接进行712的检测，并且对结点C没有左边下行。

如果B是活动的，则在714关于B调用Upd_Int_Val()，从而进行右边下行。于是递归下行块700结束，函数进行至其有效处理部分。如果B是不活动的，则对结点C没有右边下行，块700结束，函数进行至其有效处理部分。

一旦块700已循环完，并且递归也以进行完，数据处理就于结点C通过在720重新初始化清单Upd()和栈Q而开始。

清单Upd(C)用于从在块700实施递归调用时C的子结点中被更新的数据出发，逐个上溯已更新的相互作用数据。

清单Q用作对象存储栈，其中的相互作用数据需要在当前结点处更新。因此，在722，清单Q接收清单Upd(A)、Upd(B)、Int_Old_L(C)、Int_L(C)。可以注意到，尽管Int_Old_L(C)和Int_L(C)包含的是结点对而不是结点，专业人士可以认为操作722只是将结点相互分离地入栈Q。

在724，相互作用数据上溯循环开始，将清单Q出栈到工作结点R。在726，将工作结点R加入清单Upd(C)，以便保证在块700的递归出栈时关联给该结点的数据的上溯，然后，在728进行检测以便确定R是属于A所代表的对象子集，还是属于B所代表的对象子集。

如果R属于A代表的对象子集，则在730，重新用F(A，R)的值对相互使用数据F(C，R)重新初始化。否则，在732，用F(B，R)的值对F(C，R)重新初始化。最后，在734，对Q进行检测，以便指明是否Q的对象的所有相互作用数据都被重新初始化了。

当Q空后，在736，通过使清单Int_L_t(C)出栈，一个相互作用数据更新循环开始。清单Int_L_t(C)是清单Int_L(C)的一个工作副本，仅仅用于更新循环。清单Int_L_t(C)的出栈结果保存在工作结点对(P，T)中。

在738，通过向在上述循环中重新初始化了的F(C，P)值中加上T对P的相互作用F(P，T)的值，更新相互作用数据F(C，P)。F(P，T)的值是通过函数F_Int()计算出的，给出P与T的原子间相互作用，以P为参照。在740，对T进行同样的运算，通过向在上述循环中重新初始化了的F(C，T)值中加上P对T的相互作用F(T，P)的值，以T为参照。

函数F_Int()基于众所周知的两原子/分子间静电相互作用和范德华(Van derWaals)公式。当然可以用其它公式。要注意的是，F_Int()计算的相互作用是以引用的第一个结点为参照的。

然后，在742进行检测，以便验证Int_L_t(C)是否已经完全出栈。如果确已完全出栈，那么在744，Upd_Int_Val()的有关实例(instance)结束。否则，在736，更新循环重新开始。一旦清单Int_L_t(C)完全出栈，则相互作用数据F(C，X)(其中X是C的一个后代叶结点)表示结点X受到的来自于C的所有其它后代叶结点的相互作用的总和。

鉴于以上所述，显然函数Upd_Int_Val()采用的是自下而上的路径，尽管该函数是关于根结点调用的。该函数还使结点C的每个叶结点受到的C的其它叶结点的所有相互作用逐步上溯。

最后，该函数还保证跟踪通过清单Upd()更新了相互作用的数据的结点。因此，当函数Upd_Int_Val()结束时，根结点RN既包含其叶结点上相互作用的所有数据，又包含了当前循环中更新的相互作用数据的清单，从而确保只访问确实演变了的结点。

一旦函数Upd_Int_Val()结束，函数Upd_Int()也就结束了。如图2所示，在仿真循环中，Upd_Int()的下一个函数是Sprd_Int()。函数Sprd_Int()是关于根结点RN被调用的，如其清单Upd(RN)包含函数Upd_Int()修改过的所有相互作用数据。

图8显示了该函数的一个实施例。函数Sprd_Int()是一个循环，该循环于802开始，检测输入结点的清单Upd(C)。如果该清单是空的，那么在804，函数结束。

如果清单Upd(C)不空，则在806，元素R出栈，在808，相应结点的相互作用数据F(R，R)被更新为F(C，R)值。

然后，在810，函数Set_Flg()记录与R相关联的数据Flg_bph。数据Flg_bph的值指明更新动力学系数b、p、η的必要性。函数Set_Flg()还记录结点R的直至根结点的所有祖先的Flg_bph数据。然后循环重新从802开始，直至清单Upd(C)空。

可以注意到，函数Sprd_Int()可以以多种方式实施。尤其，在这里，更新后的相互作用数据被保存在有关节点中。然而，可以将这些相互作用保存在另一个地方，例如保存在一个分离的表格中。对于函数Set_Flg()保存的数据，也同样如此。

还可以不执行该函数，而只是在必要时适当地调用清单Upd(RN)，尽管这样做会导致不清晰和管理相对复杂。

当相互作用数据更新以后，仿真循环进行到函数Get_Dyn_Coeff()，该函数可以确定对象和树的结点的动力学参数。

图9显示了函数Get_Dyn_Coeff()的一个实施例。函数Get_Dyn_Coeff()是自下而上传播的函数，与函数Upd_Int_Val()一样。为了实现所述传播，在902进行检测，以验证结点C是否有子结点A和B。

如果结点C没有子结点，则在904检测数据Flg_bph(C)是否活动。如果活动，则结点C是需要在906对其调用函数Coeff_bphl()的叶结点，函数将p(C)和η(C)置零，并基于更新的相互作用数据F(C，C)计算b(C)。在步骤906以后，或者如果数据Flg_bph(C)不活动，则函数结束于908。

如果结点C有子结点，则在910和912，函数Get_Coeff_Dyn()分别关于A和B递归调用自己。这种递归调用保证了自下而上地遍历二叉树的必要结点。一旦进行了所有的递归调用，并且所有的必要子结点都被处理了，就在块920中逐个递归循环地进行树的上溯。

在块920中，首先检测当前结点C是否活动。如果活动，则在924，关于C的子结点(也就是A和B)的这些系数值调用函数Coeff_phipsi()，以便计算系数Ф(C)和Ψ(C)。然后，在926，关于C的子结点(也就是A和B)的这些系数值调用函数Coeff_bph()，以便计算系数b(C)、p(C)、η(C)。最后，在928，函数结束。

在当前结点C不活动的情况下，在930，检测数据Flg_bph(C)。如果检测结果为否，则函数直接结束于928。如果检测结果为是，也就是说该数据Flg_bph(C)指明需要更新系数b(C)、p(C)、η(C)，则在932，关于C的子结点(也就是A和B)的这些系数值调用函数Coeff_bph()，以便计算系数b(C)、p(C)、η(C)。最后，在928，函数结束。

附录A的公式(1)至(13)列出了函数Coeff_phipsi()和Coeff_bph()确定这些系数时所用的公式。就本说明书而言，只需要知道：递归地确定的这些系数可以局部地描述关节联结体的行为，并用作仿真中计算度量、加速度、位置的基础。关于更多的信息，文献SPM‘05和SIGGRAPH ‘05阐述了这些方程的精确定义和形成。申请人推荐为此目的查阅这些文献。

尤其可以注意到，系数b(C)、p(C)、η(C)的计算可能直接或间接地产生一个量Q，该量Q表示结点C的关节上的作用力数据。该量Q例如可以用来表示一个依赖于二面角的力矩。因为与结点C相关联的量Q的计算只取决于C的子结点A和B的相对位置，所以只对活动结点更新该量。因此，例如可以就在用函数Coeff_bph()计算b(C)、p(C)、η(C)之前进行所述更新。

当动力学系数被更新以后，函数Upd_Dyn()根据这些系数来更新动力学数据。图10显示了函数Upd_Dyn()的一个实施例。

图10所示的函数Upd_Dyn()是一个递归函数，该函数自上而下地经过二叉树。

该函数采用一个优先级队列，目的在于确定活动结点，如下文所述。为此，在1002对当前结点进行检测，如果当前结点是根结点RN，则在1003调用加速度度量的计算。计算结果在1004用作界限(borne)ε的基础，在1006，界限ε用作结束条件。

加速度度量的计算是由函数Tot_Acc()实施的，函数Tot_Acc()实施附录A的公式21所表达的计算，对应于当前结点的子结点的加速度的平方和。该计算可以定义活动结点和不活动的结点。

事实上，εmax定义了希望得到的、用加速度度量表达的仿真精度。因此，如果ε＜εmax，则意味着达到了希望的精度。于是不计算结点C的加速度qdd(C)，在1010，函数结束。

反之，如果ε＞εmax，则认为仿真还不够精确。在1012进行检测，以便确定结点C是否有子结点。如果没有，则在1010，函数结束。

在1014，调用函数Active_Flg(C)。函数Active_Flg()记录一个包含时间标记的数据，该时间标记与正在进行的仿真循环相关联。就是该时间标记被用在结点活动性检测中。为此，对该给定结点，将其可能关联的时间标记与当前循环的时间标记相比较。如果时间相同，则认为结点活动。可以用类似方法激活数据Flg_bph(C)。

如果结点C有子结点A和B，则在1016，用函数Q_Dot计算结点C的加速度qdd(C)。函数Q_Dot()使用动力学系数Ф(C)、b(C)以及力F(C)。函数Q_Dot()进行的计算由附录A的方程22表述。

力F(C)是一个合成矢量，表示施加在结点C上的关节力的总和。对于根结点RN，有一些特殊的条件可以确定矢量F(C)，如SPM’05的3.1.2节所述。

如文章SPM’05的3.1.1节所述，矢量F(C)有两个分量：F(1，C)和F(2，C)。就本发明而言，结点的矢量F(C)用于计算其子结点F(A)和F(B)的矢量。

事实上，如果考虑C的子结点A和B，则力F(1，C)等于F(1，A)，而力F(2，C)等于F(2，B)。此外，分量F(2，A)和F(1，B)还可以从F(C)得到，如下文所述。因此，F(A)和F(B)分量完全可以从F(C)得到。

一旦结点C的加速度计算出来，则函数Upd_Dyn()进行至1018和1020，用根据附录A的公式23的函数Ch_F()计算矢量F(2，A)和F(1，B)。

然后，在1022和1024，计算结点A和B的加速度度量。在1026，根据Acc(A)和Acc(B)的值，结点A和B进入优先级队列Q_prio。在1028，将加速度qdd(C)的平方值赋给ε，以便更新仿真的精度。最后，在1030，结点R从队列Q_prio出栈，在1032，关于R调用Upd_dyn()，以便实现该函数的递归和树中下行。

队列Q_prio根据结点的加速度度量值将所述结点分类，以便快速确定最重要的结点。因此，鉴于结点的加速度度量表示其子结点的加速度的平方和，函数就在使ε下降最快的结点中逐步下行，从而加快仿真速度。

可以注意到，该函数同时执行两个任务。事实上，结点的活动性是基于加速度度量值的。因此，不计算某些结点的加速度，以及使用数据Active_Flg()还可以将所述结点排除在后续的、考虑活动性的循环之外，而不需要额外的计算成本。

上述的将所述结点排除在后续的、考虑活动性的循环之外之所以能实现，尤其基于如下事实：自下而上遍历树的函数开始于只在活动结点中的递归下行。

因此，以上所有函数通过选择活动结点保证了对树的最小遍历，并选择性地更新某些结点的数据。

还可以注意到，此处仿真精度是通过参数εmax保证的，参数εmax表示加速度的最大误差。也可以用其它手段，例如要仿真的结点的数量规定，或其它度量。

当动力学数据更新以后，就进入仿真循环的最后一个函数Upd_Geo()。事实上，正如上文所提到的，结点中的相互作用数据是相对于每个结点的特定局部方位标(repère)计算出来的。

其实，结点的加速度和推导出来的位置也是相对的，或固有的。因此，为了能显示仿真的结果并能使数据在树的各不同经历中上溯(或再下行)，需要一些过渡矩阵，直到现在还没有讨论这些过渡矩阵。

正是这些过渡矩阵被用于计算图4的函数Upd_Pos()中OBB的实际位置，或者被用于更新图7的步骤732和734中的力F(C，R)。

为了更新这些矩阵，函数Upd_Geo()在1102以一个子函数Upd_Hyb_qdd()开始。函数Upd_Hyb_qdd()计算修正的动力学系数，以便考虑函数Upd_Dyn()实施的刚化。然后，用修正的动力学系数重新计算结点的加速度qdd(C)。关于该函数实施的运算，以及这些运算往qqd(C)计算中的整合，在文章SIGGRAPH’05的第4节和第5节中有描述。

尽管函数Upd_Hyb_qdd()在本发明的优选实施方式中被描述，但该函数并非在所有情况下都必需。作为一种选择，函数Upd_Geo()不包含该函数，直接进入1104。

函数Upd_Geo()在1104执行函数Upd_qpos()，该函数Upd_qpos()根据结点的加速度更新其固有(intrinsèque)位置。这一计算很简单，代表一种基本的时间积分，对应于附录A的方程14。还可以采用其它方法来更新结点的固有位置，例如进行额外的加速度计算。

最后，函数Upd_Geo()在1106执行最后一个函数Upd_CT_Mat()，该函数Upd_CT_Mat()根据结点的已更新的固有位置来更新上文提到的内部矩阵。这些矩阵及其更新在SPM’05的第5节的靠前部分有描述。

在前文中，相互作用数据和几何数据主要是以固有方式计算的，也就是说是相对于每次与有关结点关联的局部方位标计算的。

以固有方式计算可以不必重新计算只依赖于对象相对位置的数据，例如系数Ф(C)、Ψ(C)以及局部相互作用力。因此，只是对于结点的现实位置计算才需要转换到现实世界的坐标。

根据另一实施方式，可以每次都对现实世界的坐标的数据计算，只对不活动结点保留固有数据。

尽管本发明是用二叉树描述的，但也可以用其它类型的树，例如三叉树甚至n叉树。还可以用其它表示类型，只要这些表示类型可以区分活动结点与不活动结点并且遍历成本差不多就行。

上文中的函数体现了实施本发明的主要步骤。所述函数可以由专业人士用作基础，以便编写实施本发明的计算机程序。这种程序可以用C++语言或任何其它的、专业人士认为适当的计算机语言编写。

附录A

节段1

${\mathrm{\Φ}}^C=\left[\begin{array}{cc}{\mathrm{\Φ}}_1^C& {\mathrm{\Φ}}_{12}^C\\ {\mathrm{\Φ}}_{21}^C& {\mathrm{\Φ}}_2^C\end{array}\right]---\left(1\right)$

${\mathrm{\Φ}}_1^C={\mathrm{\Φ}}_1^A-{\mathrm{\Φ}}_{12}^{A}W{\mathrm{\Φ}}_{21}^A$

${\mathrm{\Φ}}_2^C={\mathrm{\Φ}}_2^B-{\mathrm{\Φ}}_{21}^{B}W{\mathrm{\Φ}}_{12}^B---\left(2\right)$

${\mathrm{\Φ}}_{21}^C={\mathrm{\Φ}}_{21}^{B}W{\mathrm{\Φ}}_{21}^A={\left({\mathrm{\Φ}}_{12}^C\right)}^T$

其中

$V={({\mathrm{\Φ}}_2^A+{\mathrm{\Φ}}_1^B)}^{-1}$ (3)

W＝V^-VS(S^TVS)^-1S^TV

其中S是结点的运动子空间，其维数：6×结点自由度数。典型地，S＝(0，0，1，0，0，1)^T

$\mathrm{\Ψ}\left(C\right)=\left[\begin{array}{cc}{\mathrm{\Ψ}}_1^C& {\mathrm{\Ψ}}_{12}^C\\ {\mathrm{\Ψ}}_{21}^C& {\mathrm{\Ψ}}_2^C\end{array}\right]$ (4)

${\mathrm{\Ψ}}_1^C={\mathrm{\Ψ}}_1^A+{\mathrm{\Φ}}_{12}^{A}Y{\mathrm{\Φ}}_{21}^A-({\mathrm{\Φ}}_{12}^{A}W{\mathrm{\Ψ}}_{21}^A+{\mathrm{\Ψ}}_{12}^{A}W{\mathrm{\Φ}}_{21}^A)$

${\mathrm{\Ψ}}_2^C={\mathrm{\Ψ}}_2^B+{\mathrm{\Φ}}_{21}^{B}Y{\mathrm{\Φ}}_{12}^B-({\mathrm{\Φ}}_{21}^{B}W{\mathrm{\Ψ}}_{12}^B+{\mathrm{\Ψ}}_{21}^{B}W{\mathrm{\Φ}}_{12}^B)$

${\mathrm{\Ψ}}_{21}^C=-{\mathrm{\Φ}}_{21}^{B}Y{\mathrm{\Φ}}_{21}^A+{\mathrm{\Φ}}_{21}^{B}W{\mathrm{\Ψ}}_{21}^A+{\mathrm{\Ψ}}_{21}^{B}W{\mathrm{\Φ}}_{21}^A={\left({\mathrm{\Ψ}}_{12}^C\right)}^T$

其中

U＝(S^TVS)^-1S^TV

N＝U^TU (5)

$Y=N+W({\mathrm{\Ψ}}_2^A+{\mathrm{\Ψ}}_1^B)W$

注：矩阵Ф(C)和Ψ(C)是矩阵的矩阵

附录A(续)

节段2

$b\left(C\right)=\left[\begin{array}{c}{b}_1^C\\ b_2^C\end{array}\right]$ (11)

$b_1^C=b_1^A-{\mathrm{\Φ}}_{12}^A\mathrm{\γ}$

$b_2^C=b_2^B+{\mathrm{\Φ}}_{21}^B\mathrm{\γ}$

其中

$\mathrm{\β}=b_2^A-b_1^B$

γ＝Wβ+VS(S^TVS)^-1Q，

其中Q是结点的作用力的向量，其维数：自由度数×1

$p\left(C\right)=\left[\begin{array}{c}{p}_1^C\\ p_2^C\end{array}\right]$

$p_1^C=p_1^A-{\mathrm{\Φ}}_{12}^{A}Z-2{\mathrm{\Ψ}}_{12}^A\mathrm{\γ}---\left(12\right)$

$p_2^C=p_2^B+{\mathrm{\Φ}}_{21}^{B}Z+2{\mathrm{\Ψ}}_{21}^B\mathrm{\γ}$

${\mathrm{\η}}^C={\mathrm{\η}}^A+{\mathrm{\η}}^B+R^{T}R+{\mathrm{\γ}}^T({\mathrm{\Ψ}}_2^A+{\mathrm{\Ψ}}_1^B\mathrm{\γ}-{\mathrm{\γ}}^T(p_2^A-p_1^B)$

其中

R＝(S^TVS)^-1(Q-S^TVβ)(13)

$Z={2U}^{T}R+W(p_2^A-p_1^B)-2W({\mathrm{\Ψ}}_2^A+{\mathrm{\Ψ}}_1^B)\mathrm{\γ}$

附录A(续)

节段3

$\mathrm{Acc}\left(C\right)={\left[\begin{array}{c}{f}_1^A\\ f_2^B\end{array}\right]}^T\left[\begin{array}{cc}{\mathrm{\Ψ}}_1^C& {\mathrm{\Ψ}}_{12}^C\\ {\mathrm{\Ψ}}_{21}^C& {\mathrm{\Ψ}}_2^C\end{array}\right]\left[\begin{array}{c}{f}_1^A\\ f_2^B\end{array}\right]+{\left[\begin{array}{c}{f}_1^A\\ f_2^B\end{array}\right]}^T\left[\begin{array}{c}{p}_1^C\\ p_2^C\end{array}\right]{\mathrm{\η}}^C---\left(21\right)$

$\mathrm{qdd}\left(C\right)={\left(S^T\mathrm{VS}\right)}^{-1}(Q-S^{T}V({\mathrm{\Φ}}_{21}^A{f}_1^A-{\mathrm{\Φ}}_{12}^B{f}_2^B+\mathrm{\β}))---\left(22\right)$

$f_1^B=W{\mathrm{\Φ}}_{21}^A{f}_1^A-W{\mathrm{\Φ}}_{12}^B{f}_2^B\mathrm{\γ}$ (23)

$f_2^A=-f_1^B$

Claims (7)

1.一种用于仿真相互作用的对象集合并确定对象的相互作用力、加速度和位置的相继更新的计算机设备，包括：

存储器(8)，适用于保存对象集合的树状表示，其中所述树状表示包括含有叶结点和中间结点的结点，所述叶结点表示关节联结对象的最基本对象以及所述中间结点每个都是表示一个对象子集的一个关节联结体，其中每个中间结点与动力学数据、几何数据和相互作用数据相关联，所述动力学数据、几何数据和相互作用数据取决于其各子结点的数据，其中所述相互作用数据表示施加于某给定结点或该结点表示的对象子集中的某给定对象的力，以及

仿真控制器(4)，用于按反复进行的循环启动：

对象、对象子集相互作用数据分配器(10)；

根据有关相互作用数据和几何数据更新动力学数据的机构(12)，为倾向于相互作用的结点遍历树状表示；

根据已更新的动力学数据更新几何数据的机构(14)，为倾向于相互作用的结点遍历树状表示，

其特征在于：

存储器(8)还包含与至少一些所述中间结点相关联的局部相互作用数据，其中与一中间结点相关联的局部相互作用数据表示由所述各子结点表示的对象被由所述各子结点表示的其它对象施加的力，而不是由于所述对象和所述其它对象之间通过关节联结的相互作用导致的应力；

相互作用数据分配器(10)包括更新相互作用数据的机构(20)，为倾向于相互作用的结点遍历树状表示，以便根据有关子结点的局部相互作用数据更新中间结点的相互作用数据；局部相互作用数据包括相互作用清单，所述相互作用清单包括指明局部相互作用的结点的对；

更新相互作用数据的机构包括：

更新相互作用清单的数据的函数(32)；

根据至少一些已更新相互作用清单的数据来更新相互作用数据的函数(34)。

2.根据权利要求1所述的计算机设备，其特征在于：

对于至少一些结点，存储器包括有向盒数据，所述有向盒数据表示相关联的结点周围的相互作用空间；

更新相互作用数据的机构还包括更新有向盒的数据的函数(30)，更新相互作用清单的函数根据至少一些已更新的有向盒数据进行操作。

3.根据权利要求1或2所述的计算机设备，其特征在于：

对于至少一些结点，存储器包括动力学系数数据；

更新动力学数据的机构包括：

更新动力学系数数据的函数(26)；

根据已更新动力学系数数据，更新动力学数据的函数(28)。

4.根据权利要求1或2所述的计算机设备，其特征在于：

对于至少一些结点，存储器包括表示活动性标记(Active_Flg( ))的数据；

更新动力学数据的机构(12)还能更新其处理的结点的活动性标记(Active_Flg( ))；

更新相互作用数据的机构、更新动力学数据的机构和更新几何数据的机构中的至少一个被配置为忽略未标记为活动的结点地遍历树状表示。

5.一种仿真相互作用对象集合的行为并确定对象的相互作用力、加速度和位置的相继更新的方法，其中使对象集合保持树状表示，其中所述树状表示包括含有叶结点和中间结点的结点，所述叶结点表示关节联结对象的最基本对象以及所述中间结点每个都是表示一个对象子集的一个关节联结体，其中每个中间结点与动力学数据、几何数据和相互作用数据相关联，所述动力学数据、几何数据和相互作用数据取决于其子结点数据，其中所述相互作用数据表示施加于某给定结点或该结点表示的对象子集中的某给定对象的力，所述方法包括：

a.通过计算机分配一些对象和对象子集的相互作用数据；

b.通过计算机为倾向于相互作用的结点遍历树状表示，根据有关相互作用数据和几何数据，更新这些结点的动力学数据；和

c.通过计算机为倾向于相互作用的结点重新遍历树状表示，根据在步骤b中更新的动力学数据，更新几何数据，

其特征在于：步骤a包括：

a1.为倾向于相互作用的结点预遍历树状表示，根据有关子结点的局部相互作用数据，更新中间结点的相互作用数据，其中与中间结点相关联的局部相互作用数据表示由所述各子结点表示的对象被由所述各子结点表示的其它对象施加的力，而不是由于所述对象和所述其它对象之间通过关节联结的相互作用导致的应力；步骤a1包括：

a1a.更新形成局部相互作用数据的与至少一些结点相关联的相互作用清单，每个清单包括一些有局部相互作用的结点对；

a1b.根据已更新的相互作用清单来更新相互作用数据。

6.根据权利要求5所述的方法，其特征在于：步骤a11包括：

a1a1.更新与至少一些结点相关联的有向盒的数据，有向盒的数据表示相关联的结点周围的相互作用空间；

a1a2.根据已更新的有向盒的数据来更新相互作用清单。

7.根据权利要求5至6中任一项所述的方法，其特征在于：步骤a、b、c中的至少一个步骤根据与至少一些结点相关联的活动性标记并忽略至少一些结点地来遍历树状表示。