CN101655832B - 一种基于标量场梯度的物理变形方法 - Google Patents

Abstract

一种基于标量场梯度的物理变形方法，步骤为：(1)预处理：确定计算域，预计算模型的有向距离场，建立节点应变标量函数；(2)受力计算：利用应变标量函数的梯度计算能量，根据能量计算节点受力；(3)时间积分：先根据受力更新速度，其次根据速度首先处理非对流项，然后再处理对流项；(4)处理拓扑变化；(5)计算域更新；(6)表面提取与绘制。其中步骤(2)至(6)为一个仿真循环。本发明首先将不规则计算域转化为规则计算域从而简化计算。其次，本发明显式记录局部形变状态而无需保留原始模型的信息，在大变形与拓扑改变时也无需进行重采样与重新网格化。再次，本发明基于欧拉计算框架，对处理拓扑变化有较好的灵活性。

Description

一种基于标量场梯度的物理变形方法

技术领域

本发明属于计算机虚拟现实和计算机图形学技术领域，特别是涉及计算机图形学基于物理的变形仿真方法。

背景技术

物理变形仿真技术如今已广泛应用于各领域：数字娱乐、工业制造、教育、医学手术等。模拟复杂的变形过程(大变形、用户与物体之间进行交互、物体与物体间交互、变形物体拓扑发生变化)是当前物理变形研究的趋势。传统的物理变形方法：质点弹簧系统、有限元方法、粒子系统方法，其计算模型都是分布在空间的不规则网格或散乱的点。在复杂的物理变形中，变形的实现涉及求解大量偏微分方程，从这个角度说，规则的计算域是更适合的。首先，规则计算域有利于简化计算、提高求解的稳定性。其次，规则计算域有利于并行算法的实现。

任何基于物理的变形方法，如何定量衡量材料局部所发生的变形(应变)是一个关键性问题，因为它是后续计算能量、受力的前提。传统方法是通过“对比”的方式进行计算的，即在整个变形过程中记录原始非变形模型(参考模型)的信息，然后用变形模型与参考模型做对比从而计算出局部形变(应变)。具体地，质点弹簧系统通过对比弹簧长度来计算沿弹簧的一维应变。基于连续介质力学的有限元方法与无网格方法，通过比较对应点变形前后的位置求得位移场，然后估计位移场的梯度计算应变张量。然而在处理大变形或拓扑改变时，当前的计算模型要发生变化：如大变形时增加采样点，拓扑改变时进行重新网格化。另外，参考模型的结构也必须与当前计算模型一致，因此还必须进一步更新参考模型。这些都给计算带来不小的开销与潜在的出错几率。

与上述基于拉格朗日框架的变形方法相比，基于欧拉的与标量场相结合的方法具有灵活处理拓扑变化的优势；其次，标量场在空间上可以定义在非常规则的栅格上。已有的基于标量场的变形技术在隐式雕刻(implicit sculpting)应用中已有广泛研究。这种技术又称“虚拟橡皮泥”(virtual clay techniques)，它是一种人机交互的建模方式，一般提供力反馈，使用户有如同操作真实材料进行模型雕刻与塑造的体验。然而这种技术并没有提供足够的基于物理的机制来计算应变、定义能量等，因此一般只适用于静态的建模而不适用于动态的变形仿真。Hua Jing，Hong Qin曾提出一种“密度弹簧系统”(参见Hua Jing，Hong Qin.Haptics-based Volumetric Modeling Using DynamicSpline-based Implicit Functions.2002.)，其中他们把三维弹簧系统扩充到四维，第四维是标量场的值。这种方法在标量场与材料力学属性上建立了映射关系，从而可以实现简单的变形动画。但是这种四维弹簧机制并不是物理上真实的，因此不宜直接应用到全局变形，其所能实现的变形仿真应用是有限的。

发明内容

本发明要解决的技术问题：首先将不规则计算域转化为规则计算域。其次，将变形问题转化为处理规则域上的标量场，建立一种物理上合理的标量场到形变的映射关系，用变量函数显式记录形变状态，无需记录与更新原始模型的信息，避免在大变形与拓扑改变时更新参考模型、重采样与重新网格化。

本发明采用的技术方案：一种基于标量场梯度的物理变形方法，其特点在于步骤如下：

(1)预处理

初始化材料的力学参数，确定计算域，预计算模型的有向距离场，建立节点应变标量函数；

(2)受力计算

首先通过应变标量函数计算能量，并计算能量关于节点位移的梯度，其次根据能量关于节点位移的梯度计算节点受力；

(3)时间积分

根据受力更新速度，其次根据速度首先处理非对流项，然后再处理对流项；

(4)处理拓扑变化

首先检测有拓扑改变的计算节点，其次通过局部修改节点的距离场值实现表面的分离或融合；

(5)计算域更新

将变形后物体新覆盖的栅格加入计算域，将未覆盖的栅格从计算域中删除；

(6)表面提取与绘制

利用Marching Cubes算法从距离场中提取变形体的表面并绘制。

其中步骤(2)至(6)为一个仿真循环。

所述步骤(2)中根据应变标量函数计算能量的步骤如下：

(1)用应变标量函数关于空间的梯度计算应变能；

(2)用应变标量函数对应矩阵的特征值计算保体积能。

所述步骤(3)中处理非对流项的步骤如下：

(1)对于节点的每个邻居，沿相对速度的反向找到下一时刻将到达该邻居的点并计算该点的应变标量函数值；

(2)用计算所得的应变标量函数值更新邻居当前的值并重新拟合节点的应变标量函数。

本发明与现有技术相比的有益效果在于：首先将传统方法的不规则计算域(空间分布不规则的网格或粒子)转化为规则计算栅格，将物体的几何属性与物理变形属性均记录为栅格上的标量场(分别对应于有向距离场与应变标量函数)，规则结构能够简化计算、提高计算的稳定性，有利于并行算法的实现。其次，本发明提出一种新的计算变形的机制，初始时在每个计算节点通过应变标量函数定义一族间隔均匀(梯度的大小恒为1)的等值球面，变形过程中等值面受到扰动成为椭球面而不再均匀(梯度大小发生变化)，因此可以通过等值面的梯度衡量局部变形。这种方法的物理背景在于，这些变形椭球与该点的应变张量相对应。这种机制显式记录局部形变状态(应变标量函数)而无需保存变形体的初始信息(参考模型)计算应变，并且该形变状态能保证随着变形物体一起运动刻画局部形变，因此在大变形与拓扑改变时无需重建参考模型，也无需进行重采样与重新网格化。再次，本发明基于欧拉计算框架(计算栅格固定，只更新栅格节点上的各属性值)，对处理拓扑变化有较好的灵活性。最后，通过对应变标量函数的直接控制，可以实现不同程度的塑性变形并且避免计算的不稳定性。

附图说明

图1为本发明的整体流程示意图；

图2、图3为本发明定义的节点应变标量函数示意图；

图4为本发明的节点应变标量函数非对流项更新示意图；

图5为本发明的节点应变标量函数对流项更新示意图。

具体实施方式

下面结合附图与实施例对本发明作进一步详细描述：

本发明实施过程包括六个主要步骤：预处理、受力计算、时间积分、处理拓扑变化、更新计算域、表面提取与绘制。如图1所示。

步骤一：预处理，主要分成三个阶段。

第一个阶段：初始化材料属性

这一阶段的主要任务是初始化计算区域，设置计算节点材料的密度、刚度等属性。本发明中的标量场定义在正交、等间距的规则栅格上，其中计算区域是所有与模型相交的栅格，计算节点为计算域中栅格的节点。密度与刚度均匀分布，在各节点的值相同，其值根据具体需要仿真的材料而定。计算节点的体积规定为一个立方栅格的体积，节点的质量可由密度与体积算出。另外这一阶段还应预计算公式(11)中的矩阵B。

第二个阶段：预计算距离场

本发明中变形体的表面定义为有向距离场下的一个零等值面，因此在预处理的这个阶段必须为计算域中的所有节点计算初始模型的有向距离场值u。

第三个阶段：建立节点应变标量函数

节点应变标量函数的建立是本发明思想的核心内容之一，其目的是在每个计算节点上建立一种标量场与局部变形(应变)之间的映射关系。本发明中应变标量函数的建立基于这一直观认识：材料的局部形变可以通过标量场局部梯度的变化来刻画。以一个单向变形为例，如图2所示，左图中，在物体内部定义了一个标量场，其等值面为一组等距(梯度的大小恒为1)的水平面。右图中材料在竖直方向局部受压，将变形映射到标量场上，则该变形可以直观地从标量场的等值面上看到：在等值面密的区域(梯度的大小大于1)材料沿梯度方向受压。同理可得：在等值面疏的区域(梯度的大小小于1)材料沿梯度方向受拉。因此可以将材料的局部形变映射为标量场在该点梯度大小的变化量。本发明根据这一思想，在所有计算节点处都初始化定义一个应变标量函数，这个标量函数在该点的邻域梯度大小恒为1。当该点及周围材料发生变形时，将变形映射到标量函数(详细内容见步骤三中对应变标量函数的更新)，此时通过计算标量函数的梯度大小变化，就能得到该点的应变。

具体地，对任一节点P，考虑到要能够表示任意方向的变形，本发明将上述节点应变标量函数定义为：

f_P(x)＝(x^TTx)^1/2-r₀，     (1)

其中x＝(x，y，z)^T为P的邻域任意一点相对P点的坐标，r₀是栅格的间距，T是一个3×3的对称、正定矩阵，在预处理阶段设为单位矩阵。如图3(左)所示，该函数实际上是点P的有向距离场函数，其等值面是以P为球心的等距同心球面，并且在P的去心领域内有：

$|\▿f_P\left(x\right)|\≡1.---\left(2\right)$

在发生变形时，T不再是单位阵，如图3(右)所示，其等值面成为椭球面，这些椭球可以认为是P点应变张量在空间上的一种可视化。根据公式(1)可知，P的标量函数只与T有关，由于对称性，它只有6个独立的分量，在计算过程中，该函数记录为以这6个分量为元素的一个系数向量：s(a，b，c，d，e，f)。初始化时，该向量为s(1，1，1，0，0，0)。

步骤二：受力计算，分为以下两个阶段。

第一个阶段：计算能量

本发明在这一阶段定义两种能量：应变能与保体积能。应变能的大小反映了计算节点周围材料的变形大小，而保体积能反映了节点周围材料体积变化的大小。

对于应变能，根据步骤一，应变可以表示成节点应变标量函数的梯度大小的变化，从而P的应变能可以定义为：

$E_T=\frac{1}{n}{k}_{T}v\underset{j}{\mathrm{\Σ}}{\left(\right|{\▿f}_P\left(x_j\right)|-1)}^2,---\left(3\right)$

其中，k_T为材料刚度，v为节点的体积(定义为一个立方栅格所占空间的体积)，j是P的任意一个邻居节点，n是邻居的数量。

对于体积变化，可以描述成变形椭球的体积与原始球体体积的差，椭球的体积正比于D＝(det(T))^-1/2，表示椭球三个半轴长度的乘积，由于初始时T为单位矩阵，则D＝1，因此P的保体积能可以定义为：

E_V＝k_Vv(D-1)²，                (4)

其中，k_V为保体积常数，由用户设定，v为节点的体积。

第二阶段：计算节点力

对于任意节点P，首先计算P对每个邻居节点j的受力f_j，根据弹性力学，力可以表示成能量关于位移的负梯度：

$f_j=-({\▿E}_T+{\▿E}_V),---\left(5\right)$

这个力将累加到节点j的受力中。根据牛顿第三定律，P应受到邻居节点的反作用力：

$f_P=-\underset{j}{\mathrm{\Σ}}{f}_j,---\left(6\right)$

这个力将累加到节点P的受力中。公式(5)中梯度的计算都可以转化为计算节点标量函数系数s(a，b，c，d，e，f)关于位移的梯度，这可以由公式(11)对邻居X_j的位移求偏导数得出(其中偏导数计算只涉及公式(11)中g(X_j ^t)部分，其余为常量)。

步骤三：时间积分，分为以下三个阶段。

第一个阶段：用节点力更新节点速度

根据牛顿第二定律，更新节点的速度：

v←v+(F/m)·Δt，           (7)

其中m为质量，F为节点所受合力，Δt为每次时间积分的时间步长。

第二个阶段：非对流项的更新

这个阶段是对节点应变标量函数(其系数向量s)非对流项的更新。事实上，应变标量函数随时间的变化可以用以下微分方程描述：

$\frac{\mathrm{ds}}{\mathrm{dt}}=s_{t,\mathrm{neighbor}}+s_{t,\mathrm{convect}},---\left(8\right)$

等式右边的第一项即为非对流项，第二项为对流项(对流项将在下一个阶段更新)。非对流项的更新可以理解为：节点P处的物质与其邻居节点j处的物质发生相对运动，使P处应变状态发生改变，从而P的应变标量函数发生变化。如图4所示，节点P邻居位置的物质相对P位置的物质发生向外的相对运动，从而使P处的材料在相应方向发生拉伸，而P的应变标量函数随之变化。对于这一变化，本发明根据以下原理进行更新：对于P邻域内，位置随时间改变的任一动点，在任意不同时刻，该动点始终停留在应变标量函数的同一个等值面上，即对任意两个不同时刻t1、t2，有

$f_P^{t1}\left(x^{t1}\right)\≡f_P^{t2}\left(x^{t2}\right),---\left(9\right)$

其中x^t1、x^t2分别为动点在这两个时刻相对P的位置，f_P ^t1、f_P ^t2分别为这两个时刻P的标量函数。具体的更新策略分为以下步骤：

(1)如图4所示，在当前时刻t，对P的所有邻居X_j，j＝1，2，...，n，计算X_j相对P的速度v_XP，并从X_j出发沿这一相对速度反向延长v_XP·Δt的距离找到下一时刻将要到达X_j的点X_j′；

(2)计算f_P ^t(X_j′)，根据公式(9)，则有：

$f_P^{t+\mathrm{\Δt}}\left(X_j\right)=f_P^t\left({X_j}^{\′}\right);---\left(10\right)$

(3)根据n个点(X_j，f_P ^t(X_j′))，j＝1，2，...，n，用最小二乘法重新拟合P在新时刻r+Δt的应变标量函数f_P ^t+Δt(系数向量)：

$s={\left[\underset{j}{\mathrm{\Σ}}b\left(X_j\right)b{\left(X_j\right)}^T\right]}^{-1}\underset{j}{\mathrm{\Σ}}b\left(X_j\right)g\left({X_j}^{\′}\right),---\left(11\right)$

其中，

$g\left({X_j}^{\′}\right)={(f_P^t\left({X_j}^{\′}\right)+r_0)}^2.---\left(12\right)$

矩阵 $B={\left[\underset{j}{\mathrm{\Σ}}b\left(X_j\right)b{\left(X_j\right)}^T\right]}^{-1}$ 仅与P和其邻居的相对位置有关，由于本发明采用规则栅格，这个相对位置在空间与时间上都是不变的，因此B为常量，可以且仅需预计算一次。

第三阶段：对流项的更新

本阶段实际是求解微分方程(8)等式右边的第二项。对流项的更新可以理解为：由于本发明基于欧拉框架，即栅格的节点固定，而仅仅更新节点上的属性，因此一个节点上的属性值会随着速度“移动”到其它节点，从而引起节属性的改变。本阶段需要更新的属性有节点应变标量函数的系数向量s，有向距离场值u，以及速度本身。如图(5)所示(该图示以二维情况为例，实际中为三维)，以对s的更新为例对节点P进行更新，利用semi-Lagrangian方法，具体步骤如下：

(1)从P出发沿P速度的反向找到下一时刻将要到达P的点P′；

(2)用P′所在体元节点P，Q，R，S的s插值得到P′的s；

(3)将P的s更新为P′的s。

对速度也应进行上述更新，方法一样，不再赘述。对于有向距离场值u，利用level set方法进行更新(level set方法参见Osher，S.，and Sethian，J.1988.Frontspropagating with curvature-dependent speed：Algorithms based on Hamilton-Jacobiformulations)。

步骤四：处理拓扑变化。

由于物体的表面定义为有向距离场的零等值面，因此可以通过局部修改距离场值实现拓扑改变。具体地，首先找出拓扑发生变化的节点，如应变超过某个阈值或在用户定义的切割工具作用范围内的节点，若距离场值小于零，则将其修改为一个很小的正数，这样便能实现表面的分离。表面融合的过程是类似的，仅需将相应节点距离场值由正改为负。另外为准确起见，可以在提取等值面后根据提取的表面重新计算一次距离场。该过程无需进行任何重采样与重新网格化，充分体现了标量场的优势。

步骤五：计算域更新。

计算域只包含那些与物体相交的栅格，因此在物体变形与移动的过程中计算域是动态变化的。这个过程只需处理边界上的栅格体元：如果一个边界体元各顶点的距离场值均为正，则将其从计算域中删除；如果一个边界体元的某顶点距离场值在本次更新中由正变负，则对于那些以该节点为顶点且不在计算域内的体元，将其加入计算域，并通过外插值计算这些体元顶点的各属性值。

步骤六：表面提取并绘制。

在该步骤中，首先利用Marching Cubes算法从有向距离场中提取三角面片，然后对这些面片进行绘制。

Claims (4)

1.一种基于标量场梯度的物理变形的仿真方法，其特征在于步骤如下：

(I)预处理：初始化材料的力学参数，确定计算域，为计算域中的所有节点计算初始模型的有向距离场值，建立节点的应变标量函数；所述标量场定义在正交、等间距的规则栅格上，其中所述计算域是所有与初始模型相交的栅格，所述节点为计算域中栅格的节点；变形体的表面定义为有向距离场下的一个零等值面，为计算域中的所有节点计算初始模型的有向距离场值u；对任一节点P，应变标量函数定义为：

，       (1)

其中

为P的邻域内任意一点相对于P点的坐标；r₀是栅格的间距，T是一个3×3的对称、正定矩阵，在预处理阶段设为单位矩阵；P的应变标量函数是点P的有向距离场函数，P的应变标量函数的等值面是以P为球心的等距同心球面，并且在P的去心邻域内有：

在发生变形时，T不再是单位阵，应变标量函数的等值面成为椭球面，这些椭球面认为是P点应变张量在空间上的一种可视化；根据公式(1)，P的应变标量函数与T有关，由于T的对称性，T只有6个独立的分量，在计算过程中，将P的应变标量函数记录为以这6个分量为元素的一个系数向量：s(a，b，c，d，e，f)；初始化时，该向量为s(1，1，1，0，0，0)；

(II)受力计算：首先通过应变标量函数计算应变能和保体积能，应变能反映计算节点周围材料的变形大小，保体积能反映节点周围材料体积变化的大小；并计算应变能和保体积能关于节点位移的梯度，其次根据应变能和保体积能关于节点位移的梯度计算节点受力；

(III)时间积分：分为三个阶段，第一个阶段为用节点力更新节点速度；第二阶段是对节点应变标量函数非对流项的更新，其更新节点应变标量函数的系数向量s；第三阶段是对节点应变标量函数对流项的更新，其更新节点的应变标量函数的系数向量s、节点的有向距离场值u、以及节点的速度v；

(IV)处理拓扑变化：首先检测有拓扑改变的计算节点，其次通过局部修改节点的距离场值实现表面的分离或融合； 

(V)计算域更新：将变形后物体新覆盖的栅格加入计算域，将未覆盖的栅格从计算域中删除；

(VI)表面提取与绘制：利用Marching Cubes算法从距离场中提取变形体的三角面片并绘制；

其中步骤(II)至(VI)为一个仿真循环。

2.根据权利要求1所述的基于标量场梯度的物理变形的仿真方法，其特征在于：所述步骤(I)中建立节点应变标量函数：该标量函数的等值面为一族椭球面，反映了该节点周围材料变形的分布。

3.根据权利要求1所述的基于标量场梯度的物理变形的仿真方法，其特征在于：所述步骤(II)中通过应变标量函数计算应变能和保体积能的步骤如下：

(II-1)用应变标量函数关于空间的梯度计算应变能；

(II-2)用应变标量函数对应矩阵的特征值计算保体积能。

4.根据权利要求1所述的基于标量场梯度的物理变形的仿真方法，其特征在于：所述步骤(III)中处理非对流项的步骤如下：

(III-1)对于节点的每个邻居，沿相对速度的反向找到下一时刻将到达该邻居的节点并计算该节点的应变标量函数值；

(III-2)用计算所得的应变标量函数值更新邻居当前的值并重新拟合节点的应变标量函数。 

Images