CN112862957A - 一种基于约束投影的gpu并行试衣仿真方法 - Google Patents

Abstract

本发明提供了一种基于约束投影的GPU并行试衣仿真方法，通过Morton码对人体模型并行构建二叉基数树形成BVH，每次碰撞检测时直接用衣物顶点遍历BVH进行测试；衣物仿真利用Galerkin多重网格法加速基于速度求解的Projective Dynamics方法，对衣物顶点采用线性蒙皮坐标对衣物顶点进行自由度扩充，在每一帧迭代的局部求解中，对衣物顶点施加各种几何约束投影，并对与人体碰撞的衣物顶点进行Signorini‑Coulomb约束的投影，进而在不改变全局求解中的系统矩阵的情况下，仅通过局部求解迭代中的约束投影实现对试衣中衣物的碰撞响应与摩擦物理效果的仿真控制。

Description

一种基于约束投影的GPU并行试衣仿真方法

技术领域

本发明属于计算机图形学领域与物理仿真领域，尤其涉及一种基于约束投影的GPU(Graphic Processing Unit，图形处理器)并行试衣仿真方法。

背景技术

用计算机仿真试衣中的人衣交互具有广泛的应用场景，如虚拟试衣、游戏场景、VR应用等，人衣交互主要涉及到衣物的仿真以及衣物人体碰撞检测与响应两个方面的技术问题，基于物理的方法对该场景进行仿真是目前主流的方法，而该方法的核心就是求解牛顿力学指导下的运动方程。

目前主流的关于衣物的物理仿真的方法主要包含以下两类：

一、基于力学的仿真方法：该方法直接利用牛顿第二定律的相关公式作为系统数值积分求解的目标方程，常见的模型有质点弹簧模型。

二、基于能量优化的仿真方法：该方法通过将隐式欧拉积分重新进行构造，得到一个优化问题，通过求解一个能量最小化问题来得到系统在对应时间步的关于惯性势能和弹性势能之间的稳态位置。基于位置的动力学方法(PBD\XPBD)、基于投影的动力学方法(Projective Dynamics)、基于交替方向乘子的方法(ADMM)等方法都属于该方法的范畴。

衣物和人体的碰撞检测和响应是仿真出合理的试衣效果必不可少的步骤，目前碰撞检测通用的方法是对场景中的物体分别构建BVH(Bounding volume hierarchy，层次包围盒)，然后通过遍历BVH进行碰撞检测。而碰撞响应方法主要包含两类：

一、离散碰撞响应，离散碰撞响应仅仅基于当前的状态，不需要额外的历史信息。距离场方法、相交轮廓最小化方法等方法都属于离散碰撞响应的范畴。

二、连续碰撞响应，连续碰撞响应需要上一时刻无碰撞的历史状态作为输入进行参考。排斥力冲量法、碰撞区法等方法都属于连续碰撞响应的范畴。

目前关于试衣的物理仿真方法的实现，基于CPU(Central Processing Unit，中央处理器)的实现往往远远达不到实时仿真的要求。随着显卡的广泛应用，并行计算的应用普及，为基于GPU的试衣仿真得到实时的仿真结果提供了可能。

对于衣物的物理仿真的方法中，1)基于力学的仿真方法；2)基于能量优化的仿真方法。方法1较为简单，但是在仿真的稳定性以及仿真物理效果的可扩展性较弱；方法2是目前业界采用的方法，可扩展性强，但目前对实时人衣交互中如人衣摩擦、衣物材质等物理效果的仿真探索缺少GPU实现上的扩展。

人衣之间碰撞检测和碰撞响应是人衣交互中限制仿真速度提升的最大瓶颈，在每一个时间步中对变化的人体与衣物重新构建BVH并进行碰撞检测、响应时间开销很大。

发明内容

发明目的：本发明主要针对当前计算机对3D(三维)试衣仿真中对衣物人体之间交互缺乏物理效果仿真的问题，提供一种基于约束投影的GPU并行试衣仿真方法，包括如下步骤：

步骤1，在每一帧中利用输入的新的人体模型以及上一帧输入的人体模型的信息进行层次包围盒的并行构建或者更新；本发明中，人体和衣物模型可以从法国国家信息与自动化研究所INRIA申请许可证然后下载得到，为Obj(3D模型文件格式)文件，许多Obj文件构成一个连续的人体运动的输入数据。而衣物模型就是一个Obj文件，通过运算进行模型更新。帧是指一个动画帧，也就是计算结果绘制一次算一帧。

步骤2，并行建立需要构造Signorini-Coulomb(西尼奥里尼-库伦)约束的点集；

步骤3，使用Galerkin(伽辽金)多重网格法加速基于速度求解的ProjectiveDynamics(基于投影的动力学)方法，在局部求解的每次迭代中对步骤2得到的点集中的点进行Signorini-Coulomb约束投影，对所有衣物顶点进行通用的几何约束投影，最后经过迭代后得到当前帧衣物与人体的无碰撞冲突的位置。

步骤1包括：用输入的相邻两帧的人体模型同一个编号的面片在相邻两帧的时间间隔内形成的线性轨迹，计算所述线性轨迹对应的AABB(轴对齐)包围盒中心作为Morton(莫顿)码编码的依据，通过Morton码进行编码，对Morton码进行去重并按照升序排序，从而将三维的空间信息映射到一维的空间；

根据获得的Morton编码并行构造BRT(Binary Radix Tree，二叉前缀树)，所述二叉前缀树BRT作为人体模型的层次包围盒BVH的组织层次。

步骤1中，所述二叉前缀树BRT包括内部节点和叶节点，对二叉前缀树BRT的内部节点按照深度优先的遍历顺序进行重新排序；

每一个叶节点对应于一个Morton码，叶节点是升序排列的；设定有N个叶节点，对于二叉前缀树BRT则有N-1个内部节点；每一个内部节点涵盖的范围就是对应子节点所涵盖的Morton码的取值范围。

步骤1中，只在初始化人体模型或者碰撞检测速度下降为原有记录的最快检测速度的2倍时进行二叉前缀树BRT的构建或重新构建。

步骤2包括：并行地将每一个衣物顶点与人体模型进行碰撞检测，并用几何映射到面的方式对有冲突的顶点进行碰撞响应，通过几何的方法消除由于新的人体模型更新所产生的少量碰撞冲突。然后根据衣物顶点所受外力，包括重力、风力等，对顶点下一时刻的位置进行预测，根据预测的衣物顶点位置再次与新的人体位置进行碰撞检测，将所有新的产生碰撞冲突的衣物顶点纳入需要构造Signorini-Coulomb约束的点集当中；

步骤2包括：通过对预测的衣物顶点位置与人体位置之间进行碰撞检测，对于检测到的每一个产生碰撞冲突的衣物顶点都会对应于一个库伦摩擦锥，并且所述衣物顶点需要被添加入Signorini-Coulomb约束点集，e代表摩擦锥接触的法向，r代表摩擦锥的局部接触力，u代表摩擦锥中碰撞的衣物顶点与对应发生碰撞的人体三角面片之间的相对速度；局部接触力和全局接触力、相对速度和速度之间的关系用公式1a、公式1b来描述，其中P_j表示第j个碰撞冲突的旋转矩阵，用于改变第j个碰撞对的局部接触力r_j与相对速度u_j的基底，将它们转换成第i个顶点的全局接触力ξ_i和全局速度v_i：

ξ_i＝P_jr_j (公式1a)

v_i＝P_ju_j (公式1b)

步骤3包括：参考文献(Bouaziz，S.，Martin，S.，Liu，T.，Kavan，L.，&Pauly，M.(2014).Projective dynamics：Fusing constraint projections for fastsimulation.ACM transactions on graphics(TOG)，33(4)，1-11.)中由隐式欧拉积分推导得到的Projective Dynamics方法，其全局求解的目标方程为公式2，其中M代表顶点质量矩阵，h代表时间步长，q_n+1、q_n分别代表在时刻n+1顶点的位置向量和时刻n顶点的位置向量，A_i、B_i是关联于第i个顶点的系数矩阵，S_i是关联于第i个顶点的筛选矩阵，ω_i是关联于第i个顶点权重系数，p_i是Projective Dynamics方法引入的关联于第i个顶点的辅助变量：

公式2此时的Proiective Dynamics全局求解的目标方程的未知量为时刻n+1的顶点的位置组成的列向量q_n+1，将目标方程的未知量——位置向量q_n+1，转换为时刻n+1的速度向量v_n+1，则此时基于位置求解的Projective Dynamics方法被转换成基于速度求解的ProjectiveDynamics方法，从而得到公式3，根据公式4、5，对公式3后续公式进行简写，其中b(p)是公式3中等式右边部分的简写，A_sys是公式3中等式左边括号内表达式的简写：

步骤3中，使用层次数为3的Galerkin多重网格法加速公式3的求解，具体包括：

步骤3-1，根据公式4、5将公式3简写为A_sysx＝b(p)的形式，其中用x代表未知量v_n+1，对A_sysx＝b(p)的x向量进行初始化，设定为x₁；

步骤3-2，对A_sysx＝b(p)进行平滑，更新x₁；

步骤3-3，计算残差r₁＝b(p)-A_sysx₁；

步骤3-4，利用Galerkin多重网格法第1层对应的约束矩阵R₁对残差向量r₁和系统矩阵A_sys进行下采样，得到第2层次的系统矩阵

目标向量b₂＝R₁r₁；

步骤3-5，对A₂e₂＝b₂进行平滑得到修正量e₂；

步骤3-6，计算残差r₂＝b₂-A₂e₂；

步骤3-7，利用Galerkin多重网格法第2层对应的约束矩阵R₂对残差向量r₂和系统矩阵A₂进行下采样，得到第3层次的系统矩阵

目标向量b₃＝R₂r₂；

步骤3-8，利用共轭梯度法求解A₃e₃＝b₃，得到第3层次的修正量e₃；

步骤3-9，利用Galerkin多重网格法第2层对应的插值矩阵

对修正量e₃进行插值，得到

步骤3-10，对步骤3-5中的e₂进行误差更新：e₂＝e₂+e′₂；

步骤3-11，利用Galerkin多重网格法第1层对应的插值矩阵

对修正量e₂进行插值，

步骤3-12，对步骤3-2中的x₁进行误差更新：x₁＝x₁+e₁；

步骤3-13，对A_sysx＝b(p)进行平滑，得到最后的x₁值作为结果返回。

步骤3中，利用Galerkin多重网格法在相邻层次进行约束和插值时，能够使用线性蒙皮坐标进行顶点自由度的扩充，而线性蒙皮坐标可以用仿射矩阵表示，如公式6所示第i个顶点的坐标x_i，等价于由多重网格相邻低层次所有控制点的仿射矩阵乘其自身齐次坐标X_i的叠加，其中T_j表示第j个控制点的仿射矩阵，ω_ij是第j个下采样点与第i个顶点共同决定的权重系数；对比于基于位置求解的Projective Dynamics方法，在基于速度的Projective Dynamics方法中，对应方程组中的蒙皮坐标信息将会转换成蒙皮坐标的变化率，将每个仿射矩阵行用vec()函数向量化得到公式7中的q向量，并用其表示蒙皮坐标，此时原始坐标和蒙皮坐标以及原始速度与蒙皮坐标下速度的关系如公式7所示，其中U是顶点坐标以及蒙皮坐标的转换矩阵，它将作为使用蒙皮坐标的多重网格法第1层到第2层的插值矩阵，即

根据Galerkin多重网格法的性质，在对应层次之间的约束矩阵为U的转置，即R₁＝U^T，k为第一层下采样的顶点数；公式8给出了U矩阵的计算方式，U矩阵是一个分块矩阵，U_ij是U矩阵中第i个顶点与第j个下采样点所索引的块矩阵，其中n为衣物的顶点数，

表示克罗内克乘积，I₃是维度为3的单位矩阵，ω_ij是第j个下采样点与第i个顶点共同决定的权重系数：

x_i＝∑_jω_iiT_jX_i，T_j∈R^3×4，X_i∈R^4×1 (公式6)

当多重网格法使用大于等于3个层次时，多重网格法中第k层次的修正量e_k下采样到第k+1(k≥2)层次的目标向量b_k+1变换关系及其插值矩阵U_k的计算方式分别如公式9、10所示，其中k≥2，原始网格为第1层，I₁₂是维度为12的单位矩阵，U_k+1矩阵是一个分块矩阵，U_k|ij是U_k矩阵中第k层的第i个顶点与第k+1层的第j个下采样点所索引的块矩阵，ω_k|ij是第k层的第i个顶点与第k+1层的第j个下采样点共同决定的权重系数：

e_k＝U_kb_k+1 (公式9)

U_k|ij＝ω_k|ij I₁₂ (公式10)

步骤3包括：为了在基于速度的Projective Dynamics方法中对产生Signorini-Coulomb约束的点进行对应的约束投影，由公式3、5得到公式11，其中∑_iξ_i(r)是Signorini-Coulomb约束的投影项，ξ_i(r)是一个关于第i个顶点的局部接触力r的函数，具体的物理意义为第i个顶点所受的全局接触力：

A_SYSv_n+1＝b(p)+∑_iξ_i(r) (公式11)

对公式3进行变换得到公式12、13、14，其中公式12通过雅可比风格的变换，得到公式13中关于碰撞脉冲系统的近似力f；通过力f的方向与大小来近似碰撞冲突产生的力，第j个碰撞冲突的局部接触力r_j、相对速度u_j与与该碰撞相关的第i个衣物顶点力f_i的局部关系由公式14表示，其中

是一个旋转矩阵，将力f_i的全局方向改为局部接触中的相对方向：

关于局部接触力r_j以及相对速度u_j的取值，根据

与对应的库伦摩擦锥的接触法向e之间的关系，分以下三种情况进行讨论，库仑摩擦锥示意图如图3所示，图中用A、B代表库伦摩擦锥相关的两个物体：

脱离：当

向量在接触法向e方向的分量

时，此时根据公式9a的脉冲物理含义，法向相对速度u_i|N必然大于零，A、B物体在下一帧将会脱离接触，此时局部接触力r_j将会置零；

滑动：当

向量在接触法向e方向的分量

且切向分量

时，此时A、B物体将会产生紧贴滑动，此时局部接触力r_j的法向分量的值为

切向分量的值为

法向相对速度u_j|N＝0，切向相对速度方向为局部接触力r_j切向分量的反方向；

粘滞：当

向量在接触法向e方向的分量

且切向分量

时，此时A、B物体将会处于粘滞状态，此时局部接触力r_j为

相对速度u_j＝0。

利用GPU在每一帧中并行构建BVH，然后并行进行碰撞检测与响应是提高仿真速度的有力方法，针对现有GPU实时解决方案中缺少对试衣中人衣交互中物理效果的具体仿真，本发明提供了一种基于GPU的并行试衣实时仿真方案，通过Morton码对人体模型并行构建二叉基数树形成BVH，并对BVH内部节点按照深度优先顺序进行重新排列以提高遍历BVH时GPU的缓存命中率，每次碰撞检测时直接用衣物顶点遍历BVH进行测试；同时衣物仿真利用Galerkin多重网格法加速基于速度求解的ProiectiveDynamics方法，对衣物顶点采用线性蒙皮坐标方法进行自由度扩充，在每一帧迭代的局部求解中，对衣物顶点施加各种几何约束投影，并对与人体碰撞的衣物顶点进行Signorini-Coulomb约束的投影，进而在不改变全局求解中的系统矩阵的情况下，仅通过局部求解迭代中的约束投影实现对试衣中衣物的碰撞响应与摩擦物理效果的仿真控制。

有益效果：本提案提供了一套完整的GPU并行试衣技术解决流程，通过Galerkin多重网格法加速非线性系统求解，可以满足实时仿真的需求；通过并行构建的重排序BVH树以及并行的局部约束投影以对试衣中人衣交互的部分物理效果进行仿真，从而对虚拟试衣的真实感仿真做出贡献。

附图说明

下面结合附图和具体实施方式对本发明做更进一步的具体说明，本发明的上述和/或其他方面的优点将会变得更加清楚。

图1是本发明方法的渲染管线及处理流程示意图。

图2是基于Morton码构造的BRT示意图。

图3是库伦摩擦锥示意图。

具体实施方式

本发明提供了一种基于约束投影的GPU并行试衣仿真方法，图1是本发明方法的渲染管线及处理流程示意图。渲染管线中，包括：

(1)首先在每一帧中利用新输入的人体模型与上一帧输入的人体模型的信息进行层次包围盒的并行构建或者更新。

(2)在上一帧的衣物和人体经过渲染管线的一次完整仿真处于无碰撞冲突状态，而衣物当前帧还未进行物理步进时，单纯由于人体的运动必然会产生新的冲突状态，因此在渲染管线中，先处理由于人体模型变化引发的碰撞冲突。在该步骤中并行地将每一个衣物顶点与人体模型进行碰撞检测，并用几何映射到面的方式对有冲突的顶点进行碰撞响应，通过几何的方法消除由于新的人体模型更新所产生的少量碰撞冲突。然后根据衣物顶点所受外力，包括重力、风力等对衣物顶点在当前帧的位置利用简单的牛顿第二定律对下一时刻衣物顶点的位置进行预测，然后再次进行碰撞检测，将所有新的产生碰撞冲突的衣物顶点纳入需要构造Signorini-Coulomb约束的点集当中。

(3)对于新的碰撞冲突，其碰撞响应的处理融于时间积分的求解当中。使用Galerkin多重网格法加速基于速度求解的Projective Dynamics方法，并在多重网格法中使用线性蒙皮坐标提高层次之间传递的顶点自由度，在局部求解的每次迭代中对步骤(2)得到的点集中的点进行Signorini-Coulomb约束投影，对所有衣物顶点进行通用的几何约束投影，最后经过迭代后得到当前帧衣物与人体的无碰撞冲突的所有顶点的位置。

以下渲染管线中涉及到的核心算法的详细阐述：

一、并行构建层次包围盒：

由于在仿真中涉及到的人体模型是多个连续Obj文件组成的的动作序列，其中每一帧的人体模型状态都是由某两个相邻的Obj文件插值得到的，传统的并行构建层次包围盒需要将人体模型中所有三角面片求取重心进行Morton码编码，而对于动态人体模型这种仿真情况，则需要取相邻两帧的人体模型对应面片的线性轨迹形成的AABB包围盒中心作为Morton码编码的依据。

通过Morton码进行编码，并对Morton码进行去重并按照升序排序，从而将三维的空间信息映射到一维的空间。根据获得的Morton编码并行构造二叉前缀树BRT，该BRT将作为人体模型的层次包围盒BVH的组织层次；然后对BRT的内部节点按照深度优先的遍历顺序进行重新排序，以提高后期在进行碰撞检测中遍历BRT时的GPU缓存命中率。

BRT由内部节点和叶节点构成，每一个叶节点对应于一个Morton码，因此叶节点也是升序排列的；设定有N个叶节点，N取值为自然数，对于BRT来说则有N-1个内部节点。每一个内部节点涵盖的范围就是对应子节点的Morton码的取值范围。如图2所示，这是一个BRT的示意图，内部节点和叶节点的编号代表了它们各自存储位置的顺序。黑色方框内的编号是对BRT内部节点重新排序后的编号情况。

由于人体模型不像其他可变形体会产生很大的形变，因此可以只在初始化人体模型或者碰撞检测速度下降为原有记录的最快检测速度的2倍时进行二叉前缀树BRT的构建或重新构建，而正常在每一帧的物理仿真中只需依据BRT树的层次结构进行更新对应的包围盒信息即可。

二、多重网格法加速基于速度求解的Projective Dynamics目标方程：

参考文献(Bouaziz，S.，Martin，S.，Liu，T.，Kavan，L.，&Pauly，M.(2014).Projective dynamics：Fusing constraint projections for fast simulation.ACMtransactions on graphics(TOG)，33(4)，1-11.)中由隐式欧拉积分推导得到的ProjectiveDynamics方法，其全局求解的目标方程为公式1，其中M代表顶点质量矩阵，h代表时间步长，q_n+1、q_n分别代表在时刻n+1顶点的位置向量和时刻n顶点的位置向量，A_i、B_i是关联于第i个顶点的系数矩阵，S_i是关联于第i个顶点的筛选矩阵，ω_i是关联于第i个顶点权重系数，p_i是ProjectiveDynamics方法引入的关联于第i个顶点的辅助变量：

公式1此时的Projective Dynamics全局求解的目标方程的未知量为时刻n+1的顶点的位置组成的列向量q_n+1，将目标方程的未知量——位置向量q_n+1，转换为时刻n+1的速度向量v_n+1，则此时基于位置求解的Projective Dynamics方法被转换成基于速度求解的ProjectiveDynamics方法，从而得到公式2，根据公式3、4，可以对公式2后续公式进行简写：

为了在每一帧求解公式2这样一个线性方程组，本方案在GPU中使用了Galerkin多重网格法加速线性方程组的求解。多重网格法通过对原始问题进行采样构建层次，逐级削减每一个层级的问题的规模，从而加速整个线性方程组的求解与收敛。像公式2这样的线性方程组可以简写为形如A_sysx＝b(p)的形式，使用层次数为3的Galerkin多重网格法加速公式3求解的加速算法流程如下所示：

1)根据公式3、4将公式2简写为A_sysx＝b(p)的形式，其中用x代表未知量v_n+1，对A_sysx＝b(p)的x向量进行初始化，设定为x₁；

2)利用雅可比迭代法或者高斯赛德尔迭代法对A_sysx＝b(p)进行平滑，更新x₁；

3)计算残差r₁＝b(p)-A_sysx₁；

4)利用Galerkin多重网格法第1层对应的约束矩阵R₁对残差向量r₁和系统矩阵A_sys进行下采样，得到第2层次的系统矩阵

目标向量b₂＝R₁r₁；；

5)对A₂e₂＝b₂进行平滑得到修正量e₂；

6)计算残差r₂＝b₂-A₂e₂；

7)利用Galerkin多重网格法第2层对应的约束矩阵R₂对残差向量r₂和系统矩阵A₂进行下采样，得到第3层次的系统矩阵

目标向量b₃＝R₂r₂；

8)利用共轭梯度法求解A₃e₃＝b₃，得到第3层次的修正量e₃；

9)利用Galerkin多重网格法第2层对应的插值矩阵

对修正量e₃进行插值，得到

10)对步骤3-5中的e₂进行误差更新：e₂＝e₂+e′₂；

11)利用Galerkin多重网格法第1层对应的插值矩阵

对修正量e₂进行插值，

12)对步骤3-2中的x₁进行误差更新：x₁＝x₁+e₁；

13)对A_sysx＝b(p)进行平滑，得到最后的x₁值作为结果返回。

总之，这是一个简单的3层次结构的多重网格，子问题求解上一个问题经下采样得到的残差得到修正量，通过插值修正高层次问题的解，从而提高线性方程组的求解速度与收敛速度。

为了进一步提高多重网格法的收敛性，利用Galerkin多重网格法在相邻层次进行约束和插值时，可以使用线性蒙皮坐标进行顶点自由度的扩充，而线性蒙皮坐标可以用仿射矩阵表示，如公式5所示第i个顶点的坐标x_i，等价于由多重网格相邻低层次所有控制点的仿射矩阵乘其自身齐次坐标X_i的叠加，其中T_j表示第j个控制点的仿射矩阵，ω_ij是第j个下采样点与第i个顶点共同决定的权重系数；对比于基于位置求解的Projective Dynamics方法，在基于速度的Projective Dynamics方法中，对应方程组中的蒙皮坐标信息将会转换成蒙皮坐标的变化率，将每个仿射矩阵行用vec()函数向量化得到公式6中的q向量，并用其表示蒙皮坐标，此时原始坐标和蒙皮坐标以及原始速度与蒙皮坐标下速度的关系如公式6所示，其中U是顶点坐标以及蒙皮坐标的转换矩阵，它将作为使用蒙皮坐标的多重网格法第1层到第2层的插值矩阵，也即

根据Galerkin多重网格法的性质，在对应层次之间的约束矩阵为U的转置，也即R₁＝U^T，k为第一层下采样的顶点数；公式7给出了U矩阵的计算方式，U矩阵是一个分块矩阵，U_ij是U矩阵中第i个顶点与第j个下采样点所索引的块矩阵，其中n为衣物的顶点数，

表示克罗内克乘积，I₃是维度为3的单位矩阵，ω_ij是第j个下采样点与第i个顶点共同决定的权重系数：

x_i＝∑_jω_ijT_jX_i，T_j∈R^3×4，X_i∈R^4×1 (公式5)

当多重网格法使用大于等于3个层次时，多重网格法中第k(k≥2)层次的修正量e_k下采样到第k+1(k≥2)层次的目标向量b_k+1变换关系及其插值矩阵U_k的计算方式分别如公式8、9所示，其中原始网格为第1层，I₁₂是维度为12的单位矩阵，U_k+1矩阵是一个分块矩阵，U_k|ij是U_k矩阵中第k层的第i个顶点与第k+1层的第j个下采样点所索引的块矩阵，ω_k|ij是第k层的第i个顶点与第k+1层的第j个下采样点共同决定的权重系数：

e_k＝U_kb_k+1 (公式8)

U_k|ij＝ω_k|ijI₁₂ (公式9)

三、Signorini-Coulomb约束的施加与投影：

Signorini-Coulomb定律定义了库仑摩擦锥之间三种独立的情况：脱离、滑动、粘滞三种情况，我们方案中对此进行了简单的抽象，只有和人体模型产生了碰撞冲突的衣物顶点才会利用该定律进行Signorini-Coulomb约束施加。

如图3所示，这是一个局部的库伦摩擦锥示意图，通过对预测的衣物顶点位置与人体位置之间进行碰撞检测，检测到的每一个产生碰撞冲突的衣物顶点都会被添加入Signorini-Coulomb约束点集，并对应于这样一个库伦摩擦锥。图中e代表摩擦锥接触的法向，r代表摩擦锥的局部接触力，u代表摩擦锥中碰撞的衣物顶点与对应发生碰撞的人体三角面片之间的相对速度。局部接触力和全局接触力、相对速度和速度之间的关系用公式10、公式11来描述，其中P_j表示第j个碰撞冲突的旋转矩阵，用于改变第j个碰撞对的局部接触力r_j与相对速度u_j的基底，将它们转换成第i个顶点的全局接触力ξ_i和全局速度v_i：

ξ_i＝P_jr_j (公式10)

v_i＝P_ju_j (公式11)

由于碰撞会涉及到速度跳变，用加速度对力的描述在处理该情况存在局限性，为了更好地描述这样一个脉冲的情况，对于碰撞力的描述，方案中直接利用关于速度的公式进行描述，从而由公式2、4得到公式12，其中∑_iξ_i(r)是Signorini-Coulomb约束的投影项，ξ_i(r)是一个关于第i个顶点的局部接触力r的函数，具体的物理意义为第i个顶点所受的全局接触力：

A_SYSv_n+1＝b(p)+∑_iξ_i(r) (公式12)

若把r看作未知量，那么在求解公式10时将不得不采用复杂的求解方法，因此本方案将公式10中的r看作常量，在局部迭代求解中不断调整近似值，从而简化了系统的求解。在局部求解中对r不断进行近似就是Signorini-Coulomb约束在局部进行投影的具体行为。为了对r进行近似，对公式2进行变换得到公式13、14、15，其中公式13通过雅可比风格的变换，得到公式14中关于碰撞脉冲系统的近似力f；通过力f的方向与大小来近似碰撞冲突产生的力，第j个碰撞冲突的局部接触力r_j、相对速度u_j与与该碰撞相关的第i个衣物顶点力f_i的局部关系由公式15表示，其中

是一个旋转矩阵，将力f_i的全局方向改为局部接触中的相对方向：

关于局部接触力r_j以及相对速度u_j的取值可以根据

与对应的库伦摩擦锥的接触法向e之间的关系，分以下三种情况进行讨论，库仑摩擦锥示意图如图3所示，图中用A、B代表库伦摩擦锥相关的两个物体：

脱离：当

向量在接触法向e方向的分量

时，此时根据公式9a的脉冲物理含义，法向相对速度u_j|N必然大于零，A、B物体在下一帧将会脱离接触，此时局部接触力r_j将会置零；

滑动：当

向量在接触法向e方向的分量

且切向分量

时，此时A、B物体将会产生紧贴滑动，此时局部接触力r_j的法向分量的值为

切向分量的值为

法向相对速度u_j|N＝0，切向相对速度方向为局部接触力r_j切向分量的反方向；

粘滞：当

向量在接触法向e方向的分量

且切向分量

时，此时A、B物体将会处于粘滞状态，此时局部接触力r_j为

相对速度u_j＝0。

通过以上规则，可以在局部迭代中确定每一个涉及碰撞冲突的衣物顶点的局部接触力r，从而将得到的近似值其当作常量，直接用于线性方程组的求解，最后通过迭代求解线性方程组，解决仿真中所有的碰撞冲突。

本发明提供了一种基于约束投影的GPU并行试衣仿真方法，具体实现该技术方案的方法和途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。本实施例中未明确的各组成部分均可用现有技术加以实现。

Claims (10)

1.一种基于约束投影的GPU并行试衣仿真方法，其特征在于，包括如下步骤：

步骤1，在每一帧中利用输入的新的人体模型以及上一帧输入的人体模型的信息进行层次包围盒的并行构建或者更新；

步骤2，并行建立需要构造Signorini-Coulomb约束的点集；

步骤3，使用Galerkin多重网格法加速基于速度求解的Projective Dynamics方法，在局部求解的每次迭代中对步骤2得到的点集中的点进行Signorini-Coulomb约束投影，对所有衣物顶点进行通用的几何约束投影，最后经过迭代后得到当前帧衣物与人体的无碰撞冲突的位置。

2.根据权利要求1所述的方法，其特征在于，步骤1包括：利用输入的相邻两帧的人体模型同一个编号的面片在相邻两帧的时间间隔内形成的线性轨迹，计算所述线性轨迹对应的AABB轴对齐包围盒中心作为Morton码编码的依据，通过Morton码进行编码，对Morton码进行去重并按照升序排序，从而将三维的空间信息映射到一维的空间；

根据获得的Morton编码并行构造二叉前缀树BRT，所述二叉前缀树BRT作为人体模型的层次包围盒BVH的组织层次。

3.根据权利要求2所述的方法，其特征在于，步骤1中，所述二叉前缀树BRT包括内部节点和叶节点，对二叉前缀树BRT的内部节点按照深度优先的遍历顺序进行重新排序；

每一个叶节点对应于一个Morton码，叶节点是升序排列的；设定有N个叶节点，对于二叉前缀树BRT则有N-1个内部节点；每一个内部节点涵盖的范围就是对应子节点所涵盖的Morton码的取值范围。

4.根据权利要求3所述的方法，其特征在于，步骤1中，只在初始化人体模型或者碰撞检测速度下降为原有记录的最快检测速度的2倍时进行二叉前缀树BRT的构建或重新构建。

5.根据权利要求4所述的方法，其特征在于，步骤2包括：并行地将每一个衣物顶点与人体模型进行碰撞检测，并用几何映射到面的方式对有冲突的顶点进行碰撞响应，通过几何的方法消除由于新的人体模型更新所产生的碰撞冲突，然后根据衣物顶点所受外力，对顶点下一时刻的位置进行预测，根据预测的衣物顶点位置再次与新的人体位置进行碰撞检测，将所有新的产生碰撞冲突的衣物顶点纳入需要构造Signorini-Coulomb约束的点集当中。

6.根据权利要求5所述的方法，其特征在于，步骤2中，通过对预测的衣物顶点位置与人体位置之间进行碰撞检测，对于检测到的每一个产生碰撞冲突的衣物顶点都会对应于一个库伦摩擦锥，并且所述衣物顶点需要被添加入Signorini-Coulomb约束点集，e代表摩擦锥接触的法向，r代表摩擦锥的局部接触力，u代表摩擦锥中碰撞的衣物顶点与对应发生碰撞的人体三角面片之间的相对速度；局部接触力和全局接触力、相对速度和速度之间的关系用公式1a、公式1b来描述，其中P_j表示第j个碰撞冲突的旋转矩阵，用于改变第j个碰撞对的局部接触力r_j与相对速度u_j的基底，将它们转换成第i个顶点的全局接触力ξ_i和全局速度v_i：

ξ_i＝P_jr_j (公式1a)

v_i＝P_ju_j (公式1b)。

7.根据权利要求6所述的方法，其特征在于，步骤3包括：基于ProiectiveDynamics方法，其全局求解的目标方程为公式2，其中M代表顶点质量矩阵，h代表时间步长，q_n+1、q_n分别代表在时刻n+1顶点的位置向量和时刻n顶点的位置向量，A_i、B_i是关联于第i个顶点的系数矩阵，S_i是关联于第i个顶点的筛选矩阵，ω_i是关联于第i个顶点权重系数，p_i是ProjectiveDynamics方法引入的关联于第i个顶点的辅助变量：

公式2此时的Projective Dynamics全局求解的目标方程的未知量为时刻n+1的顶点的位置组成的列向量q_n+1，将q_n+1转换为时刻n+1的速度向量v_n+1，则此时基于位置求解的Projective Dynamics方法被转换成基于速度求解的Projective Dynamics方法，从而得到公式3，根据公式4、5，对公式3后续公式进行简写，其中b(p)是公式3中等式右边部分的简写，A_sys是公式3中等式左边括号内表达式的简写：

8.根据权利要求7所述的方法，其特征在于，步骤3中，使用层次数为3的Galerkin多重网格法加速公式3的求解，具体包括：

步骤3-1，根据公式4、5将公式3简写为A_sysx＝b(p)的形式，其中用x代表未知量v_n+1，对A_sysx＝b(p)的x向量进行初始化，设定为x₁；

步骤3-2，对A_sysx＝b(p)进行平滑，更新x₁；

步骤3-3，计算残差r₁＝b(p)-A_sysx₁；

步骤3-4，利用Galerkin多重网格法第1层对应的约束矩阵R₁对残差向量r₁和系统矩阵A_sys进行下采样，得到第2层次的系统矩阵

目标向量b₂＝R₁r₁；

步骤3-5，对A₂e₂＝b₂进行平滑得到修正量e₂；

步骤3-6，计算残差r₂＝b₂-A₂e₂；

步骤3-7，利用Galerkin多重网格法第2层对应的约束矩阵R₂对残差向量r₂和系统矩阵A₂进行下采样，得到第3层次的系统矩阵

目标向量b₃＝R₂r₂；

步骤3-8，利用共轭梯度法求解A₃e₃＝b₃，得到第3层次的修正量e₃；

步骤3-9，利用Galerkin多重网格法第2层对应的插值矩阵

对修正量e₃进行插值，得到

步骤3-10，对步骤3-5中的e₂进行误差更新：e₂＝e₂+e′₂；步骤3-11，利用Galerkin多重网格法第1层对应的插值矩阵

对修正量e₂进行插值，

步骤3-12，对步骤3-2中的x₁进行误差更新：x₁＝x₁+e₁；

步骤3-13，对A_sysx＝b(p)进行平滑，得到最后的x₁值作为结果返回。

9.根据权利要求8所述的方法，其特征在于，步骤3中，利用Galerkin多重网格法在相邻层次进行约束和插值时，使用线性蒙皮坐标进行顶点自由度的扩充，而线性蒙皮坐标能够用仿射矩阵表示，如公式6所示第i个顶点的坐标x_i，等价于由多重网格相邻低层次所有控制点的仿射矩阵乘其自身齐次坐标X_i的叠加，其中T_j表示第j个控制点的仿射矩阵，ω_ij是第j个下采样点与第i个顶点共同决定的权重系数；对比于基于位置求解的ProjectiveDynamics方法，在基于速度的Projective Dynamics方法中，对应方程组中的蒙皮坐标信息将会转换成蒙皮坐标的变化率，将每个仿射矩阵行用vec()函数向量化得到公式7中的q向量，并用其表示蒙皮坐标，此时原始坐标和蒙皮坐标以及原始速度与蒙皮坐标下速度的关系如公式7所示，其中U是顶点坐标以及蒙皮坐标的转换矩阵，它将作为使用蒙皮坐标的多重网格法第1层到第2层的插值矩阵，即

根据Galerkin多重网格法的性质，在对应层次之间的约束矩阵为U的转置，即R₁＝U^T，k为第一层下采样的顶点数；公式8给出了U矩阵的计算方式，U矩阵是一个分块矩阵，U_ij是U矩阵中第i个顶点与第j个下采样点所索引的块矩阵，其中n为衣物的顶点数，

表示克罗内克乘积，I₃是维度为3的单位矩阵，ω_ij是第j个下采样点与第i个顶点共同决定的权重系数：

x_i＝∑_jω_ijT_jX_i，T_j∈R^3×4，X_i∈R^4×1 (公式6)

当多重网格法使用大于等于3个层次时，多重网格法中第k层次的修正量e_k下采样到第k+1层次的目标向量b_k+1变换关系及其插值矩阵U_k的计算方式分别如公式9、10所示，其中k≥2，原始网格为第1层，I₁₂是维度为12的单位矩阵，U_k+1矩阵是一个分块矩阵，U_k|ij是U_k矩阵中第k层的第i个顶点与第k+1层的第j个下采样点所索引的块矩阵，ω_k|ij是第k层的第i个顶点与第k+1层的第j个下采样点共同决定的权重系数：

e_k＝U_kb_k+1 (公式9)

U_k|ij＝ω_k|ijI₁₂ (公式10)。

10.根据权利要求9所述的方法，其特征在于，步骤3包括：为了在基于速度的Projective Dynamics方法中对产生Signorini-Coulomb约束的点进行对应的约束投影，由公式3、5得到公式11，其中∑_iξ_i(r)是Signorini-Coulomb约束的投影项，ξ_i(r)是一个关于第i个顶点的局部接触力r的函数：

A_SYSv_n+1＝b(p)+∑_iξ_i(r) (公式11)

对公式3进行变换得到公式12、13、14，其中公式12通过雅可比风格的变换，得到公式13中关于碰撞脉冲系统的近似力f；通过力f的方向与大小来近似碰撞冲突产生的力，第j个碰撞冲突的局部接触力r_j、相对速度u_j与与该碰撞相关的第i个衣物顶点力f_i的局部关系由公式14表示，其中

是一个旋转矩阵，将力f_i的全局方向改为局部接触中的相对方向：

关于局部接触力r_j以及相对速度u_j的取值，根据

与对应的库伦摩擦锥的接触法向e之间的关系，分以下三种情况进行讨论，用A、B代表库伦摩擦锥相关的两个物体：

脱离：当

向量在接触法向e方向的分量

时，此时根据公式9a的脉冲物理含义，法向相对速度u_j|N必然大于零，A、B物体在下一帧将会脱离接触，此时局部接触力r_j将会置零；

滑动：当

向量在接触法向e方向的分量

且切向分量

时，此时A、B物体将会产生紧贴滑动，此时局部接触力r_j的法向分量的值为

切向分量的值为

法向相对速度u_j|N＝0，切向相对速度方向为局部接触力r_j切向分量的反方向；

粘滞：当

向量在接触法向e方向的分量

且切向分量

时，此时A、B物体将会处于粘滞状态，此时局部接触力r_j为

相对速度u_j＝0。

Images