布料仿真中的自适应各向异性应变限制研究①

尤祖寰,董兰芳

(中国科学技术大学 计算机科学与技术学院,合肥 230027)

布料仿真一直是计算机图形学与虚拟现实领域研究的热门课题.布料仿真技术,旨在通过计算机构建虚拟的布料模型,模拟逼真的褶皱及折痕,从而获得高度的真实感与动态效果.近年来,随着虚拟现实的日益火热,人们更加注重真实感与代入感.同时伴随着硬件性能的进步,各种大型虚拟现实应用成为可能.因此布料仿真技术获得了学者们的广泛关注和进一步研究.

在进行布料仿真的过程中,布料的物理性质通过其本构模型得以体现.本构模型旨在尽可能真实地反映现实世界中布料所受应力与应变之间的关系.通常,将布料视为弹性材料,虽然能够较好地仿真布料的形变特性,但容易造成过度拉伸或压缩,即“超弹性”现象.应变限制算法一般指对布料质点应变情况进行一定限制,在增强真实感的同时防止其发生“超弹性”现象.

另外,在布料仿真的一些实际应用中,例如虚拟试衣场合,用户可能对服装进行局部拉伸操作,如果不加以限制,服装有可能因所受应力过大而发生非真实的过度形变.因此,需要对布料合理地施加应变限制约束,增强布料仿真时的非线性塑性特征,防止“超弹性”现象的发生.

1 仿真流程概述

布料的基本仿真流程如图1所示.

1)读入模型:读入布料模型与障碍物模型,并对读取后发生穿透的部分进行简单的分离处理;

2)添加物理属性:根据布料的本构模型及相应力学性质为布料模型添加物理属性;

3)模型运动:根据运动时各节点受力情况计算各个节点的速度与位置变化,并更新其状态;

4)碰撞检测与响应:根据模型包围盒的距离、位置等信息来计算是否发生碰撞,并进行响应处理.

图1 布料仿真流程

一般情况下,只要实现了布料的基本仿真流程,即使不做任何优化,也可以得到具有一定真实感的仿真布料,其仿真精度与读入布料的面片数量、所使用的物理模型、碰撞检测精度与响应方法均密切相关.

除此之外,在上述所有条件都已经确定了的情况下,还有一些可进一步丰富布料细节,增强布料真实感的仿真步骤,如网格重划、应变限制、塑性/脆性形变处理等等.

2 相关工作

2.1 网格重划

网格的分辨率大小调节在布料仿真中相当重要,因为逼真的布料褶皱细节需要通过高分辨率的网格进行仿真,然而仿真高分辨率的网格又会大大降低计算效率.因此在布料仿真中,网格重划无疑是兼顾仿真的实时性和真实感的有效方法.通过网格重划,仿真系统可以自适应地改变布料模型的局部分辨率.在仿真过程中,表现真实的褶皱细节时就需要用小而密的三角网格区域,而在布料较为平坦之处即可使用低分辨率的大块三角面片.

网格重划大大提高了仿真效率,同时又能最大程度地保证了应有的布料细节.目前已经存在几种不同的方法,通过使用自适应地网格重划来进行布料仿真.Hutchinson等人是最早提出网格重划算法的,其通过曲率作为质点-弹簧模型的网格重划标准[1].Villard等人将类似的方法推广到了常规四边形网格的质点弹簧模型中[2].二者当时并未考虑布料褶皱的形成.Simnett等人在前者基础上增加了压缩和碰撞作为额外细化标准,进一步地提高了细化效率[3].Lee等人则是结合质点弹簧模型与Loop细分的方法对三角网格模型进行网格重划,通过预先计算细分步骤,得到多分辨率的层次结构,从而自适应地降低线性方程组的维数,更好地进行隐式积分求解[4].Narain等人则通过动态张量场重划网格,使用边的分割,翻转以及塌陷这三种基本操作,已经能够较好地展现出布料的褶皱效果[5].

总而言之,网格重划主要根据节点法向、曲率等参数,对碰撞响应结束之后的网格进行重划分,以满足展现布料微观形变的网格精度要求.

2.2 应变限制

应变限制是对布料应变量的约束,旨在防止过度拉伸或压缩的情况发生.应变限制算法目前可分为直接约束方法与额外约束方法.

早期的研究者们普遍通过使用直接约束方法对应变进行限制.Provot最早提出将应变限制作为构建刚性弹簧的方法,对每个环节施加最大最小允许应变限制约束[6].Hauth等人通过增大布料刚度,从而尽可能地减少过度拉伸[7],但由于布料刚度的增大,布料数值求解的精度将无法保证.

综上所述,当布料受力较小时,采用直接约束处理的方法限制应变可以得到较为精确的结果,而当布料受力较大时,直接约束方法的效率将大幅降低.因此,有学者提出采用额外约束的方法对应变进行限制.Thomaszewski等人提出基于连续介质的方法,将应变限制应用到有限元方法之中,独立地约束应变张量的三个部分[8].Narain等人提出上述方法均不适用于各向异性的高度不均匀网格,通过将应变限制视为一个不等式约束非线性规划问题,对三角网格进行各向同性的应变限制约束[5].Ma等人在Narain等人的工作基础上,直接使用布料平面坐标与世界坐标的映射关系提取出经纬方向的应变,以此近似奇异值分解而得到的两个主方向的应变,得到了各向异性的应变限制约束[9].

3 应变限制

3.1 应变限制定义

假设布料厚度及其交织模式相比其弹性行为对布料形变的影响要小得多,因此本文将布料视为一个二维连续体.将布料离散为三角形网格并进行分析,设三角形顶点i的平面坐标为(ui,vi),世界坐标为xi,如图2所示.

图2 三角形网格示意图

不妨设w(u,v)为将平面坐标映射到世界坐标系的一个线性函数,这里我们并不关心w(u,v)的具体形式.考虑其偏导数,其大小决定了在对应方向上顶点的拉压程度.根据上述假设,给定任意一个三角形元素,关于其三个顶点i,j,k,定义Δx1=xj-xi,Δx2=xk-xi,Δu1=uj-ui,Δu2=uk-ui,Δv1=vj-vi,Δv2=vk-vi,则可得到Δx1=wuΔu1-wvΔv1,Δx2=wuΔu2-wvΔv2,即:

式(1)可简写为:

其中,X=[xixjxk],Δ为 3×2 的矩阵差分算子,β-1为该三角形有限元基础矩阵.利用矩阵乘法结合律,我们可以得到:

一般情况下,平面坐标(ui,vi)是不变的,仅在网格重划之后发生改变.因此,可以将视为世界坐标xi的函数.

对式(4)进行奇异值分解:

其中,U的列向量即为两个应变主方向,Σ中的两个对角线元素η1、η2即为其主应变率.通过对η1、η2进行限制,即可控制布料模型的应变情况.

3.2 约束非线性规划[3]

传统的求解应变限制问题的方法是通过雅可比迭代或高斯赛德尔迭代进行求解的,但当方程组较大时收敛慢是不可避免的一大问题.另外针对高度不规则与各向异性网格时,这种求解方法很难保证其稳定性.因此,将应变限制问题转化为约束非线性规划问题,可以更为有效地处理上述问题,同时提高计算效率.定义最优化模型如下:

其中,mn为顶点n的质量(即其所在三角形的平均值),为未受约束时所处的位置(即顶点的当前世界坐标).关于两个主应变率ηj(j=1,2)的约束条件c(x)为:

其中,a为材质平面空间下三角形的面积,tj,min为压缩应变阈值,tj,max为拉伸应变阈值.又ηj关于xn的偏导数为:

则xi方向上的应变梯度为:

3.3 基于经纬方向的应变限制[6]

Ma等人借鉴了如下拉压能量状态函数[7]

通过将3.2节中由SVD得到的两个主应变率η1、η2分别用‖wu(x)‖、‖wv(x)‖进行近似,则可以得到关于经纬方向的最优化约束条件

又定义剪切向量为ws=wu+wv,类似地,则可以增加两个关于剪切向量的约束条件

该方法一方面直接使用经纬方向的应变率近似SVD得到的主方向应变率,从而避免了SVD的开销,提高了计算效率,另一方面通过定义剪切向量保留了剪切变形,限制了拉伸应变,从而使得布料能够不表现出各向同性的行为.

通过该方法施加的应变限制约束,一方面在形变后经纬方向无法保证正交,另一方面在每次网格重划后,三角形网格的经纬方向都将发生变化,如图3所示.这样得到的应变限制虽然大体上是各向异性的,但是局部约束方向并不一定总能沿着所希望的方向.为此,本文提出一种方法,即通过内积进行映射,使得布料每个三角形网格的局部应变限制在网格重划之后仍然能够和全局应变限制保持一致,从而提高了布料仿真的准确性和真实感.

4 内积倒数映射

考虑到布料的针织显然是具有方向性的,因此经向与纬向的应变一般并不相同,然而在每一帧对网格进行重划分之后,重划分的网格中每个三角形的主方向均发生了变化.Ma等人的工作虽然实现了各向异性的应变限制,但无法较好地处理上述情景,因此针对网格中每个三角面片的主方向不固定的特点,需要定义一种新的应变限制来处理这些问题.为此我们提出了内积倒数映射的概念.

为了满足需求,考虑满足如下性质的一个映射关系:

1)若三角形面片的主方向与布料全局应变限制方向相同,则约束值保持不变;

基于上述性质,定义映射函数:

其中,ei表示布料的全局应变方向,eu表示三角面片内的局部应变单位方向,εi表示ei方向上的全局应变,εu,i表示由ei方向上的全局应变εi映射到三角面片eu方向上的局部应变.

下面确定求解约束非线性规划问题中所需的约束条件c(x)及其关于xn的梯度向量∇xnc,分别就关于SVD主应变方向与关于布料经纬方向的应变限制约束进行讨论.

4.1 约束条件的确定

首先由式(13),我们可以得到三角面片应变阈值为

其中,εmin为压缩应变阈值,εmax为拉伸应变阈值.结合式(7)和式(14)可得

4.2 梯度向量的确定

要求∇xnc,则需要先确定∇xnηj与∇xneu.由式(9)可得:

为确定∇xneu,首先需要确定eu与ei的关系,基于式(5),我们可以得到如下关系式:

其中,V的列向量即是全局应变方向到局部应变主方向的旋转矩阵.

血栓形成可以导致严重心脑血管方面疾病，其中肺栓塞、大面积脑梗死可以直接致死和致残，此类患者往往预后很差。静脉血栓方面疾病具有一定的发病率，在住院患者中尤其危重病患者亦有较高的发病率，可增加死亡率[1] 。深静脉血栓形成（deep venous thromboembolism，DVT）是最常见的静脉血栓疾病，在大手术后及危重病患者人群中，其发病率较高，往往可能导致较严重的并发症。有研究显示，可导致肺动脉高压[2] 等并发症。甚至出现意外死亡和严重伤残，这对该病的治疗及预防提出了巨大挑战。下面详细阐述国内外研究进展情况。

因此,从ei到eu的旋转映射为:

借鉴Papadopoulo等人[10,11]的工作,由SVD得到的矩阵V对原矩阵元素wij的偏导数为:

由链式法则,可以得到:

5 改进乘子法

经过上一节的详细推导,我们得到了本文提出的自适应各向异性应变限制算法的目标函数、约束条件及其关于自变量的梯度向量.下面则需要对算法得到的约束非线性规划问题进行求解.

求解约束非线性规划问题,首先想到的是罚函数法和拉格朗日乘子法.考虑到二种求解算法各有各的优点,因此本文将二者相结合,提出通过改进乘子法来处理该最优化规划问题.结合后的改进方法相当于在目标函数中同时引入了二次惩罚项与线性逼近项.迭代过程中仍是从小到大地依次增大惩罚因子序列,这样在算法过程早期当惩罚因子很小时,拉格朗日乘子的更新可以很大,即通过引入二次惩罚项使得算法的收敛速度更快,同时又因为拉格朗日乘子的加速更新,一般可以在惩罚因子趋于无穷之前就收敛到最优解,从而避免了之前提到的目标函数海塞矩阵随惩罚因子增大而失效的情况.

下面具体介绍求解方法,结合式(6)得到如下约束非线性规划问题:

记可行域D={x|x∈Rn,hi(x)≤0,i=1,2,...,m}.若存在不等式约束gj(x)≤0,可加入松弛变量将其转换为等式约束,即:

下面将罚函数法和拉格朗日乘子法的思想分别引入,即通过拉格朗日乘子和惩罚函数的组合,将约束条件加入到目标函数当中,构造函数:

其中,拉格朗日乘子λ∈Rm,惩罚因子M∈R.整个过程分两步:

1)对于固定的λ和M,通过最小化LA更新x.

2)使用前文提及的方法更新拉格朗日乘子:

与惩罚因子M.

对于足够大的M,迭代一定会收敛,且收敛时λ→λ∗,x→x∗.

6 实验结果

本文的实验仿真是在Intel(R)Core(TM)i3-2120 CPU @ 3.30 GHz,4 GB内存的机器下进行的,采用C++语言进行编写,并通过OpenGL对所有模型进行渲染与显示.

为了验证算法,本文设计了两个对比实验,即分别模拟两端悬挂布料以及仿真服装,通过观察每组实验对象的形变情况,进行算法比较与分析.

6.1 两端悬挂布料实验

实验分为三组,模拟两端悬挂布料在稳定状态下的形变情况,分别实现了基于SVD主应变方向的应变限制、基于布料经纬方向的应变限制以及本文提出的基于内积倒数映射的应变限制三种方法,具体实验构成参数如表1所示.其中,图3(a)(b)(c)垂直方向阈值为[0.99,1.01],水平方向阈值为[0.95,1.05],显然此时垂直方向要比水平方向具有更高的刚性;图3(d)、图3(e)、图3(f)则与图3(a)、图3(b)、图3(c)相反,垂直方向阈值为[0.95,1.05],水平方向阈值为[0.99,1.01],此时水平方向较垂直方向更具有刚性.

图3(a)、图3(d)施加的是基于SVD主应变方向的各向同性应变限制,即使规定的两个主应变方向下的限制不同,展现出来的整体布料形态仍十分相似,仅因为阈值不同而使得细微褶皱上有所区别,即阈值越小,布料应变越小,其褶皱也越少.

表1 布料实验参数

图3 布料对比实验结果图

图3(b)、图3(e)展示的是基于经纬方向的应变限制,因为是各向异性的应变限制,所以当经向与纬向的应变限制约束阈值不同时,影响的不仅仅是细节的褶皱部分,其展现出来的整体布料形态也不相同,即布料应变方向上的形变量受到其阈值的限制.

图3(c)、图3(f)施加的是本文所定义的应变限制,与图3(b)、图3(e)展示的整体布料形态比较相似,因为二者均受到各向异性的应变限制.但仔细观察两种方法的实验结果可以发现,在相同限制阈值的情况下,图3c、图3(f)相对于图3(b)、图3(e)而言展现了更加丰富的褶皱,以及更大的形变状态.究其原因在于本文所定义的应变限制考虑到整体阈值到局部的映射,因此图3(c)、图3(f)在垂直方向与水平方向上的效果差异要比前面两种定义更加明显.

6.2 仿真服装实验

将应变限制运用到服装仿真中,通过三组实验进一步地验证了本文方法较基于经纬方向的方法限制效果更加准确,实验参数如表2所示.如图4展示了通过三种不同的应变限制约束方法下仿真人偶在静止状态下其身着T恤的形变情况.其中,前两组水平、垂直方向阈值均为±5%,不同之处在于第一组施加的是基于布料经纬方向的应变限制,第二组施加的是本文提出的基于内积倒数映射的应变限制;为进一步地控制实验参数变量,第三组施加基于布料经纬方向的应变限制,并将其水平、垂直方向阈值设置为±3%.

表2 服装实验参数

显然地,细微褶皱应变的变化程度较为剧烈,而平坦或较大范围褶皱区域的应变变化程度则较缓慢.因此在限制了T恤的应变后,最终应该得到较为平整且仅保留较大范围褶皱的结果.可以看到,第一组使用本文方法,在应变限制阈值为±5%的情况下得到了预期结果,即细微褶皱被如期限制,如图4(a)所示.

而第二组使用基于经纬方向的应变限制方法,在相同的应变限制阈值(±5%)下并未取得预期结果.如图4(b)所示,其T恤的胸口、领口处仍存在理应受到限制的细微褶皱.

第三组表明若在施加基于经纬方向的应变限制的情况下希望得到上述效果,则需要将应变限制阈值进一步缩小.如图4(c)所示,在阈值限定为±3%的情况下,基于经纬方向的应变限制方法才能较好地保证网格中每个三角形的应变都得到了如期的限制.

图4 服装对比实验结果图

综上所述,同为自适应各向异性应变限制方法,本文所定义的内积倒数映射充分考虑到了整体阈值到局部阈值的映射,使得布料在预期方向上的应变限制效果差异要比基于经纬方向的方法更加准确.

7 结语

本文在基于网格重划的布料仿真流程中,对布料模型的应变限制添加了一种新的约束,解决了网格重划后各个三角形形变主轴的局部应变限制与宏观上规定的全局应变限制不一致的情况,并使用改进乘子法对约束非线性规划问题进行求解.同时,利用OpenGL图形库进行显示与交互,成功地完成了该布料仿真实验,取得了较其他工作更为科学准确的仿真结果.

在本文的实验案例中,内积倒数映射的作用还没有得到完全地展现.在之后的工作中,可以针对具体的布料材质可以由力学实验测得其相应的极限应变,例如布料横纹(即布料经向)、直纹(布料纬向)与斜纹(布料经纬45度方向)的拉伸刚度就具有很大的差异,并通过材料力学的相关理论确定其不同方向的应变限制阈值,从而取代人为设定的阈值.再结合本文所提出的基于内积倒数映射的应变限制,可以精确地仿真出各向异性的布料材质,从而进一步地提高布料仿真实验的准确性与真实感.

1Hutchinson D,Preston M,Hewitt T.Adaptive refinement for mass/spring simulations.Boulic R,Hégron G.Computer Animation and Simulation’96.Vienna,Austria:Springer,1996.31-45.

2Villard J,Borouchaki H.Adaptive meshing for cloth animation.Engineering with Computers,2005,20(4):333-341.[doi:10.1007/s00366-005-0302-1]

3Simnett TJR,Laycock SD,Day AM.An edge-based approach to adaptively refining a mesh for cloth deformation.Proceedings of 2009 EG UK Theory and Practice of Computer Graphics.Cardiff University,United Kingdom.2009.77-84.

4Lee Y,Yoon SE,Oh S,et al.Multi-resolution cloth simulation.Computer Graphics Forum,2010,29(7):2225-2232.[doi:10.1111/cgf.2010.29.issue-7]

5Narain R,Samii A,O’Brien JF.Adaptive anisotropic remeshing for cloth simulation.ACM Transactions on Graphics (TOG),2012,31(6):152.

6Provot X.Deformation constraints in a mass-spring model to describe rigid cloth behaviour.Graphics interface.Canadian Information Processing Society.Toronto,Canada.1995.147-147.

7Hauth M,Etzmuß O,Straßer W.Analysis of numerical methods for the simulation of deformable models.The Visual Computer,2003,19(7-8):581-600.[doi:10.1007/s00371-003-0206-2]

8Thomaszewski B,Pabst S,Straßer W.Continuum-based strain limiting.Computer Graphics Forum,2009,28(2):569-576.[doi:10.1111/cgf.2009.28.issue-2]

9Ma GH,Ye JT,Li JT,et al.Anisotropic strain limiting for quadrilateral and triangular cloth meshes.Computer Graphics Forum,2015,35(1):89-99.

10Baraff D,Witkin A.Large steps in cloth simulation.Proceedings of the 25th Annual Conference on Computer Graphics and Interactive Techniques.New York,NY,USA.1998.43-54.

11Papadopoulo T,Lourakis MIA.Estimating the Jacobian of the singular value decomposition:Theory and applications.Proceedings of European Conference on Computer Vision.Berlin Heidelberg,Germany.2000.554-570.