CN106991722A - 一种基于斜椭球域影响凸包的几何网格模型变形方法 - Google Patents

Abstract

本发明公开了一种基于斜椭球域影响凸包的几何网格模型变形方法。首先根据感兴趣区域，选定待变形控制顶点，确定从全局网格顶点逻辑正向或反向选定影响区域，建立全局或局部的斜椭球域影响凸包，遍历待变形控制顶点的环形邻域，计算网格顶点处的Laplacian坐标，对待变形区域进行最小二乘约束变形，并将其与非影响凸包区域叠加得到变形模型。本发明用最小体积外接斜椭球域表征控制顶点的影响区域，精准化包络变形顶点的影响范围，可克服Laplacian微分坐标在缩放和旋转变换中局部顶点发生旋转扭曲的弊端，变形后的网格面片法向量或顶点法向量过渡平缓均匀，具有更好的平顺性，有助于实现几何邻近而拓扑不连接区域的几何网格模型变形。

Description

一种基于斜椭球域影响凸包的几何网格模型变形方法

技术领域

本发明涉及先进制造技术领域的复杂几何模型处理方法，尤其是涉及计算机造型设计中几何网格模型的变形方法。

背景技术

网格变形是指在满足几何约束下得到一个新曲面。网格变形技术作为一种重要的造型手段，已经成为计算机图形学领域一个十分活跃的研究热点，在造型设计、变形设计影视动画等领域有着相当广泛的应用。拉格朗日坐标是嵌在质点上，随物体一起运动和变形的坐标，又称物质坐标或随体坐标。欧氏坐标是固定在空间中的坐标，又称空间坐标或固定坐标。网格变形种类有拉伸变形、压缩变形、弯曲变形和扭曲变形，扭曲变形是在弯曲变形的基础上旋扭变换。网格变形技术主要分为：自由变形、基于薄壳能量的变形、基于梯度的网格变形、基于拉普拉斯(Laplacian)坐标的网格变形、多分辨层次网格变形，以及基于解偏微分方程的网格变形等。自由变形是指不直接操作变形物体，而是将物体嵌入一空间，随所嵌空间变形而变形，从某种程度上讲仍然是建立在传统曲线和曲面造型理论基础上的。自由变形技术具有便于用户交互和高效实用的优势，但很难保持模型的几何细节，因此常用于光滑模型的变形。基于薄壳能量的变形是指满足位置约束的网格变形的薄壳能量应最小化，优点是由于位置约束是网格变形最典型、最直观的约束，可直接作用在曲面上，保持变形的细节部分，缺点是要求解大量约束方程，耗时较多。基于梯度的网格变形是指通过求解满足不同边界约束条件的Poisson方程对几何网格进行变换操作，进而约束网格变形，优点是局部微分坐标能表示几何细节，具备简单、健壮和性能优势，缺点是梯度的方向依赖于全局坐标系(Global Coordinate System，GCS)，梯度不是刚性不变量，直接基于位置约束条件求解基于梯度的变形网格难以产生满意的结果。基于Laplacian坐标的网格变形是指基于Laplacian坐标的刻画顶点均值曲率和法向量的能力，用Laplacian坐标代替梯度表示几何细节，优点是把网格变形归结位置约束的优化问题，包含了网格的局部细节特征，因此Laplacian网格变形能够较好地保持网格模型的局部细节，缺点是基于迭代的非线性优化过程计算量大，对于变化幅度较大的网格变形不理想。Laplacian坐标的性质主要包括线性变换、平移不变性和对旋转变化敏感。Laplacian坐标的表示对网格细节特征的保持起到至关重要的作用，因此Laplacian网格变形方法的关键之处在于计算顶点的Laplacian坐标，并根据Laplacian坐标的平移不变性通过求解线性系统来获得变形后网格顶点的欧氏空间坐标。由于Laplacian坐标对旋转敏感，使得网格的局部信息会发生旋转扭曲，特别是对于大尺度变形时，其扭曲尤为严重。要实现网格模型的保特征变形，不能直接使用原网格的Laplacian坐标来重建变形后的网格模型，而应该重新设置微分坐标的方向再重建模型。基于解偏微分方程的网格编辑方法直接作用在网格上，能够有效保持网格曲面的微分特性，能够反映曲面的局部几何细节，具有细节保持的性能。但是基于微分域变形技术需要求解大型稀疏线性方程组或进行非线性优化，时间复杂度较高。椭球作为常见的规则几何体，相对于平面包围体(例如轴对齐包围盒Axis Aligned Bounding Box，AABB，有向包围盒Oriented Bounding Box，OBB)等具有更广泛的拟合性，特别是对于具有弧度的弯曲区域，比单一的包围球或包围盒更加有效灵活。

许多学者提出了一些网格变形算法，例如，1986年Sederberg和Parry在《ACMSiggraph Computer Graphics》(1986，20(4):151-160.)发表论文“Free-formdeformation of solid geometric models”，首次提出自由曲面造型(FFD)方法。Catmull等人于1972年在《ACM Conference》(1972:422-431)发表“A system for computergenerated movies”，提出一种骨架驱动的变形方法。Terzopoulos等人于2015年在《ACMSiggraph Computer Graphics》(2015，21(4):205-214.)发表论文“ElasticallyDeformable Models”，提出了薄壳能量应最小化的网格变形。Alexa于2003年在《VisualComputer》(2003，19(2):105-114.)发表论文“Differential coordinates for localmesh morphing and deformation”，首次提出用均匀Laplacian坐标表示几何细节。Sorkine于2004年发表论文“Laplacian Surface Editing”，提出几何细节是一个表面的一个内在属性，因此，表面编辑是最好的通过操作一个内在的表面表示。提供了基于网格的Laplacian算子，通过编码相对于它的邻域的每个顶点，不断遍历各层环形邻域的顶点，进而表示出网格面片的关系。美国明尼苏达大学Michael Ludwig等人于2015年在《Computers&Graphics》((2015，51:146-156)发表论文“3D shape and texture morphingusing 2D projection and reconstruction”，采用降维的方法把三维变形问题转化为多个二维变形，包括改变拓扑结构、创造出变形时的外观和给用户提供控制界面。美国阿拉巴马大学Chao Peng等人于2016年在《Computers&Graphics》((2016，59:107-118)发表论文“Fast mapping and morphing for genus-zero meshes with cross sphericalparameterization”，提出一种基于交叉式球形参数化对网格的快速映射和变形的方法。荷兰代尔夫特理工大学T.Gillebaart等人于2016年在《Journal of ComputationalPhysics》(2016，321:997-1025)发表论文“Adaptive radial basis function meshdeformation using data reduction”，提出了自适应径向基函数(RBF)网格变形方法，增加网格变形的鲁棒性。

综上所述，现有的网格变形技术都是尽量保证全局变形的光滑性和均匀性，因此难以有效地保持网格模型特征。如何分析和抽取不同类型的网格特征，并在变形中保持这些具有设计和工程意义的特征，满足创意设计、点云重建、快速成型(Rapid Prototyping，RP)中几何网格模型变形需求，是未来需要进一步研究的难题之一。如何精确控制控制顶点的影响区域，是现有技术中所欠缺的。

发明内容

为了解决背景技术中存在的问题，为了表征控制顶点的影响区域，将原网格分解为影响凸包区域和非影响凸包区域，本发明的目的在于提供一种基于斜椭球域影响凸包的几何网格模型变形方法。针对用不同的相对影响凸包处理影响凸包区域和非影响凸包区域，能使变形后的网格过渡平缓均匀，具有更好的平顺性。

为了实现上述目的，如图1所示，本发明采用的技术方案的步骤如下：

第一步：导入待变形的网格模型，根据感兴趣区域(Regions of interest,ROI)作为待变形区域，将原网格模型分解为影响凸包区域和非影响凸包区域，以感兴趣区域作为影响凸包区域，感兴趣区域以外的网格模型为非影响凸包区域，待变形区域中的顶点作为待变形点，获得待变形点集P，从待变形点中选定一个待变形控制顶点V^*；

第二步：根据待变形区域占总体网格的比例，进行影响凸包的正向包络和反向包络，若选择正向包络，则包络区域内的顶点为待变形顶点，若选择反向包络，则包络区域外的顶点为待变形顶点，据此确定从全局网格顶点逻辑正向或反向选定影响区域；

第三步：以待变形控制顶点V^*为球心建立最小体积外接椭球(Minimum VolumeEnclosing Ellipsoids，MVEE)，通过仿射变换将构建的最小体积外接椭球进行斜椭球化，建立斜椭球域影响凸包，最小体积外接斜椭球包络所有待变形点，用最小体积外接斜椭球表征控制顶点的影响区域；

第四步：根据选择的变形种类，对拉伸变形和扭曲变形的不同变形种类采用区分式处理；

第五步：遍历待变形控制顶点的环形邻域，得到具有s个顶点的三角网格模型M，并用Laplacian坐标δ_i描述几何网格模型；

第六步：在最小二乘意义下，对待变形区域进行最小二乘约束变形，无论是拉伸或扭曲变形，均构建能量误差函数E(V′)并求最小值，解得变形后的顶点，从而完成变形。

求解能量误差函数E(V′)的极小值即为求解二次最小化的问题，可以转换为求解稀疏线性系统。

所述第四步中，变形种类的区分和处理具体为：

若变形种类为拉伸变形，则直接进行下一步骤；

若变形种类为扭曲变形，则通过一种矢量组合运算方法完成空间网格模型绕任意动轴的扭曲变形，之后再进行下一步骤。

所述的矢量组合运算方法具体是：

步1：建立世界坐标系OXYZ，以待变形控制顶点V^*为向量的起点并作为旋转基点，各个待变形的顶点为向量的终点，为每一个顶点添加一个变换向量；

步2：以待变形控制顶点V^*的单位法向量n为旋转轴；

步3：构造矢量V_i是第i个待变形点，i为待变形点的序号，过待变形控制顶点V^*取一个法向量为单位法向量n的平面γ；

步4：通过单位法向量n与矢量乘构造矢量v₂；

步5：通过矢量v₂与单位法向量n叉乘构造矢量v₁；由于n为单位向量，v₁与v₂正交且模相等；

步6：构造向量v＝v₁ cosθ+v₂ sinθ，其中θ为旋转角度；由于v₁与v₂正交且模相等，向量v就是向量v₁绕n旋转θ后的新向量；

步7：构造向量其中O为坐标系的原点；

步8：为方便计算，将向量进一步转化为下式：

其中，为旋转基点转化而来，V^* ₁是控制顶点V^*在点V_i旋转平面的投影点，为向量n乘以一个常数ρ，其中ρ为向量n与向量点乘的结果。

向量的末端就是点V_i旋转θ得到的新顶点V′_i，由此获得对任意顶点V_i绕法向量为n的控制顶点V^*旋转θ角度后的顶点V_i′。

所述的矢量组合运算方法是用不同的相对影响凸包处理影响凸包区域和非影响凸包区域，对原始网格模型的Laplacian坐标进行了修正，为影响凸包内的网格顶点隐式地添加变换向量，解决了Laplacian微分坐标对缩放和旋转变换变形中局部信息发生旋转扭曲的弊端。

所述第三步构建最小体积外接椭球的具体步骤为：

步1：给定点集P和包络误差，将点集P划分为多个子点集，对每个子点集生成一个包围椭球，椭球的中心式方程表示如下：

(x-P_s)^T×A×(x-P_s)＝1

其中，P_s是最小体积包围椭球的中心，在三维空间中，A是一个3×3的正定矩阵，且A的特征值为椭球半轴平方的倒数，T表示矩阵转置；

步2：在三维空间中，由椭球方程的定义可知，A是3×3的正定矩阵。对完全包围点集P的椭球方程中的正定矩阵A进行奇异值分解，得到三个矩阵U,Q,D，其中，U,D分别是酉矩阵(Unitary Matrix)，分别为左乘酉矩阵和右乘酉矩阵，U和D的共轭转置等于其逆矩阵U^H＝U^-1，D^H＝U^-1，Q是表示A的特征值的对角矩阵，即表示为椭球方程的标准形式：

其中，x_r表示椭球沿x方向的半轴，y_r表示椭球沿y方向的半轴，z_r表示椭球沿z方向的半轴；

步3：由仿射变换将直笛卡尔坐标系的椭球域进行斜椭球化，令[x,y,z]←[x,y,z]*R，其中R为仿射变换矩阵：

R＝Rz×Ry×Rx

其中，θ_x,θ_y,θ_z分别表示分别绕xyz三个方向的旋转角度；

从而获得包络全局或局部几何网格模型的斜椭球域影响凸包，通过逻辑正向或反向选定来确定影响区域。

所述第五步具体为：

步1：遍历待变形控制顶点V^*的环形邻域，得到由环形邻域中的s个顶点构成的三角网格模型M＝(V，E，F)，V为环形邻域中所有顶点构成的顶点集，E为环形邻域的边集，F为环形邻域的三角面片集合，s表示三角网格模型中的顶点总数；

顶点集V中，对于每个顶点，用传统的笛卡尔坐标V_i表示，记V_i＝(X_i，Y_i，Z_i)；

步2：通过树形数据结构表示顶点的环形邻接关系。循环遍历顶点集V中每个顶点的一阶环形邻域，获得顶点集N_i＝{V_j|(V_i，V_j)∈E}，V_j表示顶点V_i一阶环形邻域中的顶点，(V_i，V_j)表示顶点V_i和顶点V_j之间的连线；

步3：用以下公式表示的Laplacian坐标δ_i计算获得待变形网格模型中的顶点V_i：

其中，δ_i为顶点V_i的Laplacian坐标，L为网格模型的Laplacian算子，w_i,j为V_j点相对于V_i点的权值，V_j表示顶点V_i一阶环形邻域中的顶点。

所述的权值w_i,j可采用以下公式计算：

其中，card(N_i)是顶点集N_i的元素个数，作为顶点V_i的度，即一阶环形邻域中的顶点个数，j表示顶点V_i一阶环形邻域中的顶点的序数。

所述第六步中，能量误差函数为：

其中，V_i′表示变形后的顶点V_i，s表示网格模型中的顶点总数，m表示待变形点总数，V′表示变换后的待变形点，点V_i′和点V′均用欧氏坐标表示；

在能量误差函数中，对待变形点和非待变形点进行区分处理：

若V_i不为待变形点，则V_i′＝u_i，u_i表示顶点V_i的原坐标；

若V_i为待变形点，则采用上述能量误差函数求解获得变形后的顶点V_i′，此时i∈{1,…,m-1}。

所述第六步中，求解能量误差函数E(V′)的极小值具体是转换为求解以下公式的稀疏线性系统，使变形后网格上的待变形点位置逼近于指定的位置，而不变形点位置保持不变，得到影响凸包区域和非影响凸包区域叠加后的网格模型：

其中，L为Laplacian变换的系数矩阵，0表示零矩阵，I为单位矩阵，V′表示待变形点变换后的欧氏坐标，Δ表示变形后网格顶点的Laplacian坐标矩阵，U_i表示非待变形点的集合，U_i＝{u_i|i∈(m…s)}。

求解上述公式的稀疏线性系统是对系数矩阵L进行LU(下三角和上三角)分解，然后采用迭代法求解得到待变形点变换后的欧氏坐标V′。

所述的能量误差函数是指在控制顶点的影响区域内，把待变形点和不变形点区分处理并叠加，将原网格分解为影响凸包区域(待变形网格)和非影响凸包区域(不变形网格)，即得到模型的多模态变形。

本发明具有的有益效果是：

1.本发明提出的方法，可以给定斜椭球必须包围的点集，确定受控制影响变形和不受控制影响变形的顶点区域，确定从全局网格顶点逻辑正向或反向选定影响区域，进一步获得斜椭球域影响凸包的空间表达，有助于实现网格模型上几何邻近而拓扑不连接区域顶点的拉伸变形与扭曲变形的多模态变形。

2.本发明提出的多分辨层次网格变形技术，可对斜椭球域影响凸包确定的变形顶点通过矢量的组合运算完成空间网格模型绕任意动轴的旋转，变形后的网格面片法向量或顶点法向量过渡平缓均匀，具有更好的平顺性，有助于实现几何邻近而拓扑不连接区域的几何网格模型变形。

附图说明

图1是本发明方法流程总图。

图2是本发明的微分域坐标变换的多模态变形示意图。

图3是本发明的扭曲变形矢量构建示意图。

图4是本发明的实例几何网格模型图。

图5是本发明的网格模型的全局与局部斜椭球域影响凸包图。

图6是本发明的全局斜椭球域影响凸包的实例变形图。

图7是本发明的局部斜椭球域影响凸包的实例变形图。

图8是本发明的全局斜椭球域影响凸包变形网格的顶点法向量图。

图9是本发明的局部斜椭球域影响凸包变形网格的顶点法向量图。

具体实施方式

下面结合附图和实施例对本发明作进一步说明。

本发明的实施例及其实施过程如下(进一步完善实施例，核实实施例的完整性，并核查字母名称解释是否统一和缺失)：

第一步：导入待变形的网格模型，图4是本发明的实例几何网格模型图。

实例几何网格模型面片数3789，顶点数1926。在打印坐标系中，变换后网格模型，凸包围盒顶点的坐标极值：[x_min,x_max,y_min,y_max,z_min,z_max]＝[0.0000,11.6360,0.0000,9.3149,0.0000,12.3809]，凸包围盒中心坐标：(5.8180,4.6574,6.1905)，模型形心坐标：(8.0134,4.1450,3.6348)，凸包围盒空间：(11.6360,9.3149,12.3809)，凸包围盒对角线长度为19.3766。顶点法向量与Z轴正方向的夹角，最大值max为176.9595，最小值min为21.4073，平均值mean为98.6952，大于平均值mean的有937个，占48.65％。平均曲率最大值max(Cmean)＝5.1386，最小值min(Cmean)＝0.0308，平均值mean(Cmean)＝0.6843，大于平均平均曲率的顶点有651个，占33.80％。高斯曲率最大值max(Cgaussian)＝16.4826，最小值min(Cgaussian)＝0.0001，平均值mean(Cgaussian)＝0.4136，大于平均高斯曲率的顶点有381个，占19.78％。第一主分量最大值max(λ1)＝3.1806，最小值min(λ1)＝0.0002，平均值mean(λ1)＝0.2420，大于平均λ1的顶点有589个，占30.58％。第二主分量最大值max(λ2)＝8.2885，最小值min(λ2)＝0.0358，平均值mean(λ2)＝1.1265，大于平均λ2的顶点有672个，占34.89％。面片法向量与Z轴正方向的夹角，最大值max为161.5788，最小值min为1.3441，平均值mean为81.5279，大于平均值mean的有1955个，占51.60％。顶点法向量与Z轴正方向的夹角，最大值max为176.9595，最小值min为21.4073，平均值mean为98.6952，大于平均值mean的有937个，占48.65％。网格模型X方向最大值点(11.6360,5.3530,8.4176)，最小值点(0.0000,0.9198,1.3243)，Y方向最大值点(8.5910,9.3149,6.6931)，最小值点(8.0511,0.0000,7.0067)，Z方向最高点(10.9784,6.4637,12.3809)，最低点(2.4196,3.5340,0.0000)。Laplacian微分坐标表示为1926*1926个单元的稀疏矩阵。

根据感兴趣区域(Regions of interest,ROI)作为待变形区域，待变形区域中的顶点作为待变形点，获得待变形点集P，从待变形点中选定一个待变形控制顶点V^*。

第二步：以待变形控制顶点V^*为球心建立最小体积外接椭球，建立全局或局部的斜椭球域影响凸包，最小体积外接椭球包络所有待变形点，用最小体积外接椭球表征控制顶点的影响区域。

第三步：通过仿射变换将构建的最小体积外接椭球进行斜椭球化：根据空间包络关系对最小体积外接椭球进行旋转变换，直至构建的最小体积外接椭球所处的中心满足全部包络指定的顶点。

图5是本实施例的网格模型的全局与局部斜椭球域影响凸包图。全局斜椭球域影响凸包的中心(表示)为(7.1448,3.9422,4.8724)，半轴为(4.3221,5.8075,9.4114)，X方向最大值点(14.2414,7.1413,7.5622)，最小值点(0.0483,0.7430,2.1827)，Y方向最大值点(11.0857,9.7503,7.5111)，最小值点(3.2040,-1.8660,2.2337)，Z方向最高点(10.3725,7.1413,12.3754)，最低点(3.9172,0.7430,-2.6306)。根据ROI变形需求，在网格模型上选定影响凸包需包络的网格顶点为：X方向最大值点(11.5672,5.3529,8.7523)，最小值点(8.9154,4.4288,8.7060)，Y方向最大值点(10.6101,6.8204,11.0361)，最小值点(10.1289,4.1117,8.6023)，Z方向最高点(10.9784,6.4637,12.3809)，最低点(9.5647,5.6849,8.5600)，进一步从全局模型中选择网格顶点21个，构成局部影响域的边界(●表示)，坐标依次为(8.9154,8.9628,9.0919)，(9.2858,9.5647,9.8617)，(10.2120,10.4660,10.7907)，(11.1114,11.3959,11.5589)，(11.5672,11.4277,11.1941)，(10.9256,10.5893,10.1289)，(10.1334,9.5875,9.1684)，(4.4288,4.6009,4.9185)，(5.3072,5.6849,5.9827)，(6.2309,6.2866,6.2593)，(6.1228,5.9109,5.6479)，(5.3529,4.9966,4.5900)，(4.2089,4.1446,4.1117)，(4.3285,4.3516,4.3910)，(8.7060,8.7424,8.8054)，(8.8877,8.5600,8.6599)，(8.8101,8.5869,8.7123)，(8.7943,8.8260,8.7934)，(8.7523,8.7321,8.6865)，(8.5630,8.7802,8.6023)，(9.0327,8.8541,8.7406)，得到局部斜椭球域影响凸包的中心(表示)为(10.3760,5.4256,9.5983)，半轴为(1.1365,1.4408,3.0844)，X方向最大值点(11.8487,5.9893,10.5507)，最小值点(8.9033,4.8619,8.6459)，Y方向最大值点(11.0347,7.1444,11.4363)，最小值点(9.7173,3.7068,7.7603)，Z方向最高点(10.7363,6.6796,12.3653)，最低点(10.0157,4.1716,6.8313)。

最小体积包围椭球(Minimum Volume Enclosing Ellipsoids，MVEE)具有如下性质：

其中，P表示待变形点集，C_h(P)是P的凸包，MVEE(P)表示围绕MVEE(P)的中心对其缩放倍。可用d×N的矩阵代表P，d表示顶点的维数，N表示点集P中含有N个顶点。即P的本质为在d维(d≥3)中包含N个定点。

构建待变形点集P的最小体积外接椭球的具体步骤为：

步1：给定点集P和包络误差，将点集P划分为多个子点集，对每个子点集生成一个包围椭球。椭球的中心式方程表示如下：

(x-P_s)^T×A×(x-P_s)＝1

其中，P_s是最小体积包围椭球的中心，在三维空间中，A是一个3×3正定矩阵，且A的特征值为椭球半轴平方的倒数，T表示矩阵转置；

步2：确定椭球方程的参数A，P_s，得到完全包围点集P的椭球方程。在三维空间中，由椭球方程的定义可知，A是3×3的正定矩阵。对正定矩阵A进行奇异值分解，得到三个矩阵U,Q,D，其中，U,D分别是酉矩阵(Unitary Matrix)，分别为左乘酉矩阵和右乘酉矩阵，其共轭转置等于其逆矩阵U^H＝U^-1，D^H＝D^-1，Q是表示A的特征值的对角矩阵，即椭球方程的标准形式为：

其中，x_r表示椭球沿x方向的半轴，y_r表示椭球沿y方向的半轴，z_r表示椭球沿z方向的半轴；

步3：引入三个方向的旋转角度θ_x,θ_y,θ_z，由仿射变换将直笛卡尔坐标系的椭球域进行斜椭球化，令[x,y,z]←[x,y,z]*R，其中R为仿射变换矩阵：

R＝Rz×Ry×Rx

步4：获得包络全局或局部几何网格模型的斜椭球域影响凸包，通过逻辑正向或反向选定来确定影响区域。

第四步：根据选择的变形种类，对拉伸变形和扭曲变形的不同变形种类采用区分式处理；

对于扭曲变形，常见的计算方法有旋转矩阵法、四元数法等。本发明通过矢量组合运算完成空间网格模型绕任意动轴的旋转，如图3所示。具体步骤如下：

步1：选定网格上的变形控制顶点V^*，设置ROI确定出待变形的顶点数；

步2：建立世界坐标系OXYZ，以选定的控制顶点V^*为旋转基点即向量起点，各个待变形的顶点为向量的终点，为每一个顶点添加一个变换向量；

步3：对于扭曲变形，确定变形控制顶点的单位法向量n为旋转轴；

步4：如图3所示，在世界坐标系OXYZ中，给定单位向量n，一个旋转基点V^*，一个空间的任意点V_i，以构建扭曲变形矢量的方式可重新设置微分坐标的方向；

步5：构造矢量过V^*点取一个法向量为n的平面γ；

步6：通过n与叉乘构造矢量v₂；

步7：通过v₂与n叉乘构造矢量v₁。由于n为单位向量，v₁与v₂正交且模相等；

步8：构造向量v＝v₁ cosθ+v₂sinθ。由于v₁与v₂正交且模相等，向量v就是向量v₁绕n旋转θ后的新向量；

步9：构造向量为方便计算，向量可进一步转化为下式：

其中，为旋转基点转化而来，V^* ₁是控制顶点V^*在点V_i旋转平面的投影点，而为向量n乘以一个常数ρ，其中ρ为向量n与向量点乘的结果。

向量的末端就是点V_i旋转θ得到的新点坐标V′_i，由此获得对任意顶点V_i绕法向量为n的控制顶点V^*旋转θ角度后的顶点坐标V′_i。

第五步：遍历待变形控制顶点的环形邻域，得到具有s个顶点的三角网格模型M，并用Laplacian坐标δ_i描述几何网格模型；

如图2所示，是本实施例构建的3D坐标分层切片平面创建示意图。用Laplacian坐标δ_i描述几何网格模型，其中δ_i的定义如下：

其中，δ_i为顶点V_i的Laplacian坐标，L为网格模型的Laplacian算子，w_i,j为V_j点相对于V_i点的权值，V_j表示顶点V_i一阶环形邻域中的顶点。

第六步：对待变形区域进行最小二乘约束变形，对于拉伸或扭曲变形构建能量误差函数E(V′)并求最小值，解得变形后的顶点，从而完成变形。

由实施例以本发明方法的全局斜椭球域影响凸包的实例变形图如图6所示。选取的控制顶点坐标为(7.8958,4.0050,4.4230)，距离全局网格最近欧氏距离为零，位于网格上，该点处法向量为(0.7735,-0.2506,-0.5821)，该点环形邻域顶点有5个，笛卡尔坐标为(7.8740,8.1180,8.2256)，(7.7321,8.0285,4.0410)，(3.6366,4.0896,3.8121)，(4.2351,4.3792,4.8710)，(4.8514,4.3042,4.5041)。人机交互变形向量为(0.4065,0.8964,0.1089)，人机交互控制顶点偏移量为4.9642，变形后控制顶点坐标(9.9136,8.4550,4.9636)。能量误差函数精度数值设置为1e-6。

实施例以本发明方法的局部斜椭球域影响凸包的实例变形图如图7所示。选取的控制顶点坐标为(7.8958,4.0050,4.4230)，距离局部网格最近欧氏距离为-4.5264，添加新网格顶点(全局第1个)为(8.9154,4.4288,8.7060)可使欧氏距离为零，该新网格顶点处法向量为(0.5221,0.7094,-0.4735)。该新网格顶点处的1环形邻域顶点有4个，笛卡尔坐标为(8.9628,8.9892,9.1684)，(9.1926,4.6009,4.6509)，(4.3910,4.5846,8.7424)，(9.0764,8.7406,9.1042)。人机交互变形向量为(0.4065,0.8964,0.1089)，人机交互控制顶点偏移量为4.9642，变形后控制顶点坐标(9.9136,8.4550,4.9636)。能量误差函数精度数值设置为1e-6。

图8是本发明的全局斜椭球域影响凸包变形网格的顶点法向量图。面片法向量与X轴正方向的夹角，最大值为179.5147°，最小值为0.5961°，平均值为90.6251°，大于平均值的有2052个，占54.16％。面片法向量与Y轴正方向的夹角，最大值为178.8383°，最小值为2.5865°，平均值为91.2862°，大于平均值的有1782个，占47.03％。面片法向量与Z轴正方向的夹角，最大值为169.0496°，最小值为4.3358°，平均值为82.2177°，大于平均值的有1953个，占51.54％。顶点法向量与Z轴正方向的夹角，最大值为178.6990°，最小值为5.0856°，平均值为89.5188°，大于平均值的有882个，占45.79％。顶点法向量与Y轴正方向的夹角，最大值为178.1338°，最小值为0.9345°，平均值为88.6945°，大于平均值的有1031个，占53.53％。顶点法向量与Z轴正方向的夹角，最大值为176.3872°，最小值为16.7448°，平均值为97.9758°，大于平均值的有931个，占48.34％。

图9是本发明的局部斜椭球域影响凸包变形网格的顶点法向量图。面片法向量与X轴正方向的夹角，最大值为179.3977°，最小值为2.1573°，平均值为89.4380°，大于平均值的有2044个，占53.95％。面片法向量与Y轴正方向的夹角，最大值为179.3487°，最小值为3.9884°，平均值为90.9403°，大于平均值的有1762个，占46.50％。面片法向量与Z轴正方向的夹角，最大值为164.3313°，最小值为4.3297°，平均值为81.5884°，大于平均值的有1951个，占51.49％。顶点法向量与Z轴正方向的夹角，最大值为179.4559°，最小值为1.9101°，平均值为90.6369°，大于平均值的有882个，占45.79％。顶点法向量与Y轴正方向的夹角，最大值为175.9330°，最小值为1.8163°，平均值为88.9624°，大于平均值的有1042个，占54.10％。顶点法向量与Z轴正方向的夹角，最大值为175.6531°，最小值为19.7440°，平均值为98.6227°，大于平均值的有928个，占48.18％。通过对比图8-图9可见，本发明使用斜椭球影响凸包，可精准化包络变形顶点的影响范围，变形后的网格面片法向量或顶点法向量过渡平缓均匀，具有更好的平顺性，有助于实现几何邻近而拓扑不连接区域的几何网格模型变形。

由此可见，本发明构建出包围指定边界点集的斜椭球；确定受控制影响变形和不受控制影响变形的顶点区域；能够进一步获得斜椭球域影响凸包的空间表达；可用Laplacian坐标在微分域内变换网格顶点实现网格变形；在ROI内实现网格变形的同时还能够保持模型的局部几何细节；通过矢量的组合运算完成网格绕任意动轴的旋转；能够避免网格的局部信息发生旋转扭曲，保证全局变形的光滑性和均匀性。

上述具体实施方式用来解释说明本发明，而不是对本发明进行限制，在本发明的精神和权利要求的保护范围内，对本发明作出的任何修改和改变，都落入本发明的保护范围。

Claims (6)

1.一种基于斜椭球域影响凸包的几何网格模型变形方法，其特征在于包括以下步骤：

第一步：导入待变形的网格模型，根据感兴趣区域作为待变形区域，待变形区域中的顶点作为待变形点，获得待变形点集P，从网格模型中选定一个待变形控制顶点V*；

第二步：根据待变形区域占总体网格的比例，进行影响凸包的正向包络和反向包络，若选择正向包络，则包络区域内的顶点为待变形顶点，若选择反向包络，则包络区域外的顶点为待变形顶点，据此确定从全局网格顶点逻辑正向或反向选定影响区域；

第三步：以待变形控制顶点V^*为球心建立最小体积外接椭球，通过仿射变换将构建的最小体积外接椭球进行斜椭球化，建立斜椭球域影响凸包，最小体积外接斜椭球包络所有待变形点，用最小体积外接斜椭球表征控制顶点的影响区域；

第四步：根据选择的变形种类，对拉伸变形和扭曲变形的不同变形种类采用区分式处理；

第五步：遍历待变形控制顶点的环形邻域，得到具有s个顶点的三角网格模型M，并用Laplacian坐标δ_i描述几何网格模型；

第六步：对待变形区域进行最小二乘约束变形，构建能量误差函数E(V′)并求最小值，解得变形后的顶点，从而完成变形。

2.根据权利要求1所述的一种基于斜椭球域影响凸包的几何网格模型变形方法，其特征在于：所述第四步中，变形种类的区分和处理具体为：

若变形种类为拉伸变形，则直接进行下一步骤；

若变形种类为扭曲变形，则通过一种矢量组合运算方法完成空间网格模型绕任意动轴的扭曲变形，之后再进行下一步骤。

3.根据权利要求2所述的一种基于斜椭球域影响凸包的几何网格模型变形方法，其特征在于：所述的矢量组合运算方法具体是：

步1：以待变形控制顶点V^*为向量的起点并作为旋转基点，各个待变形的顶点为向量的终点，为每一个顶点添加一个变换向量；

步2：以待变形控制顶点V^*的单位法向量n为旋转轴；

步3：构造矢量V_i是第i个待变形点，i为待变形点的序号，过待变形控制顶点V^*取一个法向量为单位法向量n的平面γ；

步4：通过单位法向量n与矢量叉乘构造矢量v₂；

步5：通过矢量v₂与单位法向量n叉乘构造矢量v₁；

步6：构造向量v＝v₁cosθ+v₂sinθ，其中θ为旋转角度；

步7：构造向量其中O为坐标系的原点；

步8：将向量进一步转化为下式：

其中，为旋转基点转化而来，V^* ₁是控制顶点V^*在点V_i旋转平面的投影点；

向量的末端就是点Vi旋转θ得到的新顶点V_i′，由此获得对任意顶点V_i绕法向量为n的控制顶点V^*旋转θ角度后的顶点V′_i。

4.根据权利要求1所述的一种基于斜椭球域影响凸包的几何网格模型变形方法，其特征在于：所述第三步构建最小体积外接椭球的具体步骤为：

步1：给定点集P和包络误差，将点集P划分为多个子点集，对每个子点集生成一个包围椭球，椭球的中心式方程表示如下：

(x-P_s)^T×A×(x-P_s)＝1

其中，P_s是最小体积包围椭球的中心，在三维空间中，A是一个3×3的正定矩阵，且A的特征值为椭球半轴平方的倒数，T表示矩阵转置；

步2：对完全包围点集P的椭球方程中的正定矩阵A进行奇异值分解，得到三个矩阵U,Q,D，其中，U,D分别是酉矩阵(Unitary Matrix)，分别为左乘酉矩阵和右乘酉矩阵，U和D的共轭转置等于其逆矩阵U^H＝U^-1，D^H＝U^-1，Q是表示A的特征值的对角矩阵，即表示为椭球方程的标准形式：

$\frac{{(x-P_s.x)}^2}{x_r^2}+\frac{{(y-P_s.y)}^2}{y_r^2}+\frac{{(z-P_s.z)}^2}{z_r^2}$

其中，x_r表示椭球沿x方向的半轴，y_r表示椭球沿y方向的半轴，z_r表示椭球沿z方向的半轴；

步3：由仿射变换将直笛卡尔坐标系的椭球域进行斜椭球化，令[x,y,z]←[x,y,z]*R，其中R为仿射变换矩阵：

R＝Rz×Ry×Rx

$Rx=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& cos\left({\mathrm{\θ}}_x\right)& -\sin \left(\mathrm{\θ}x\right)& 0\\ 0& \sin \left(\mathrm{\θ}x\right)& \cos \left(\mathrm{\θ}x\right)& 0\\ 0& 0& 0& 1\end{array}\right]$

$Ry=\left[\begin{array}{cccc}cos\left(\mathrm{\θ}y\right)& 0& \sin \left(\mathrm{\θ}y\right)& 0\\ 0& 1& 0& 0\\ -\sin \left(\mathrm{\θ}y\right)& 0& cos\left(\mathrm{\θ}y\right)& 0\\ 0& 0& 0& 1\end{array}\right]$

$Rz=\left[\begin{array}{cccc}cos\left(\mathrm{\θ}z\right)& -\sin \left(\mathrm{\θ}z\right)& 0& 0\\ \sin \left(\mathrm{\θ}z\right)& \cos \left(\mathrm{\θ}z\right)& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$

其中，θ_x,θ_y,θ_z分别表示分别绕xyz三个方向的旋转角度；

从而获得包络全局或局部几何网格模型的斜椭球域影响凸包，通过逻辑正向或反向选定来确定影响区域。

5.根据权利要求1所述的一种基于斜椭球域影响凸包的几何网格模型变形方法，其特征在于：所述第五步具体为：

步1：遍历变形控制点V^*的环形邻域，得到由环形邻域中的n个顶点构成的三角网格模型M＝(V，E，F)，V为环形邻域中所有顶点构成的顶点集，E为环形邻域的边集，F为环形邻域的三角面片集合，n表示三角网格模型中的顶点总数；

步2：通过树形数据结构表示顶点的环形邻接关系。循环遍历顶点集V中每个顶点的一阶环形邻域，获得环形邻域顶点集N_i＝{j|(i,j)∈E}，顶点j表示顶点i一阶环形邻域中的顶点，(i,j)表示顶点i和顶点j之间的连线；

步3：用以下公式表示的Laplacian坐标δ_i计算获得待变形网格模型中的顶点V_i：

${\mathrm{\δ}}_i=L\left(V_i\right)=V_i-\underset{j\∈N_i}{\Σ}{w}_{i,j}{V}_j$

其中，δ_i为顶点i的Laplacian坐标，L为网格模型的Laplacian算子，w_i,j为j点相对于i点的权值，V_i，V_j分别是顶点i和j的欧式坐标。

6.根据权利要求1所述的一种基于斜椭球域影响凸包的几何网格模型变形方法，其特征在于：所述第六步中求解能量误差函数E(V′)的极小值转换为求解以下公式的稀疏线性系统：

$\left(\frac{L}{0|I}\right)V^{\′}=\left(\begin{array}{c}\mathrm{\Δ}\\ U_m:U_n\end{array}\right)$

其中，L为Laplacian变换的系数矩阵，0表示零矩阵，I为单位矩阵，V′表示待变形点变换后的欧氏坐标，Δ表示变形后网格顶点的Laplacian坐标矩阵，U_m:U_n表示U_i＝{u_i|i∈(m…n)}。