CN101706830B - 对刚性材质物体表面网格模型进行钻孔后模型重建的方法 - Google Patents

Abstract

本发明公开了一种对刚性材质物体表面网格模型进行钻孔后模型重建的方法，目的是将三维布尔运算问题简化为二维布尔运算问题，降低模型布尔操作中相交检测过程的计算复杂度，使得模型重建具有实时性。技术方案是首先检测与切割曲面相交的三角形面片，将检测到的相交三角形面片曲面外网格顶点移动至切割曲面上，接着对模型表面网格预处理，实现网格模型表面三角形面片的分类划分，最后动态重建布尔运算几何处理中的网格模型。采用本发明可实现直接采用二维布尔运算处理方法实现三维空间中对网格模型的布尔运算几何处理，省略了对模型表面进行空间划分的过程，简化了相交检测算法，提高了运算效率，达到实时处理的效果。

Description

对刚性材质物体表面网格模型进行钻孔后模型重建的方法

技术领域

本发明涉及一种计算机图形学应用领域对刚性材质物体进行几何处理的方法，尤其是沿固定矢量方向对刚性材质物体表面网格模型进行钻孔后模型重建的方法。

背景技术

计算机图形学中的布尔方法是对刚性材质物体进行几何处理的有效方法之一，广泛应用于计算机辅助设计、虚拟现实、电影特效、动画制作等领域。根据物体模型的物理属性不同，可以分为弹性物体的布尔处理方法和刚性物体的布尔处理方法。前者通常用于对软组织结构构成的生物体的切割特效处理，后者则应用于对刚性材质物体的几何处理及不同模型之间融合特效处理，但经典的布尔处理方法存在表面网格模型重建困难和难以实时处理等问题。

2001年，Biermann等在ACM SIGGRAPH会议上发表的“Approximate Booleanoperations on free-form solids”(《自由形式实体的近似布尔操作方法》)论文中，将物体表面提取为三角形面片网格模型形式，并采用控制网格对表面网格模型进行多级聚类划分，得到不同程度的表面网格的简化模型，在相交检测过程中，逐步迭代求解不同物体模型之间的空间关系，最终得到近似的布尔运算结果。2002年，Museth等在“ACMTransactions on Graphics”上发表的“Level set surface editing operators”(《水平集表面的编辑操作方法》)论文中，采用水平集方法对相交的物体模型表面进行编辑操作，保证了相交表面细节特征和平滑效果。2003年，Bart等在“ACM Transactions on Graphics”上发表的“Interactive Boolean operations on Surfel-bounded solids”(《面元边界实体的交互式布尔运算方法》)论文中，将物体表面提取为矢量点的集合形式，并采用基于八叉树的空间分解方法将矢量点集合进行空间排序，实现了对复杂模型的快速、交互式布尔处理，得到较为精确的结果后，采用点的再采样方法，消除尖锐边缘实现模型表面平滑。2007年，Smith等在“Computer-Aided Design”上发表的“A topologically robust algorithm for Booleanoperations on polyhedral shapes using approximate arithmetic”(《针对多面体形状上布尔运算的一种使用近似运算的拓扑健壮性算法》)论文中，实现了对边界为多面体的物体模型的布尔运算，并使用后续平滑方法对算法得到的结果数据进行平滑处理，以消除模型中存在的失真部分。2009年，刘冬等人在《计算机辅助设计与图形学学报》上发表的《改进的布尔方法及其在虚拟ACL中的应用》中对传统布尔方法进行改进，对模型表面进行空间划分的过程进行了特殊处理，减少了布尔运算中相交检测过程的计算开销。

沿固定矢量方向对刚性材质物体表面网格模型进行钻孔可以看作对一个固定的刚性物体模型和一个不断推进的钻头模型之间进行动态的相减(difference)布尔操作。对于一个待钻孔的刚性物体表面网格模型，它在t时刻可以表示为Ω(t)＝(V(t)，F(t))，其中V(t)为此时刚性物体表面网格模型的顶点集合，F(t)为此时刚性物体表面网格模型的表面三角形面片集合。t时刻的钻头表面网格模型表示为

其中

为此时钻头表面网格模型的顶点集合，

为此时钻头表面网格模型的表面三角形面片集合，则t+1时刻得到的钻孔操作结果模型表示为Ω(t+1)＝Ω(t)-M(t+1)。在现有的布尔操作方法中，动态计算Ω(t)-M(t+1)的基本方法为：在每一个视觉循环中，根据当前模型Ω(t)和M(t+1)中顶点和表面三角形面片在三维空间中的位置，对两个网格模型进行表面网格空间划分的预处理操作，并进行相交检测和布尔运算。具体包括以下几个步骤：

第一步，根据当前时刻刚性物体表面网格模型所处三维空间中的位置，对模型表面网格进行空间划分的预处理操作。

布尔操作的核心在于对两个参与运算的模型进行相交检测。相交检测的效率直接影响整个方法的性能。如果采用蛮力检测方法对Ω(t)和M(t+1)中的每一对面片分别进行计算，则每一次相交检测的计算复杂度为O(m(t)n(t+1))，其中m(t)和n(t+1)分别为模型Ω(t)和M(t+1)中的面片数量。对于数据量大的网格模型，这种蛮力检测的方法无疑是低效的。因此，为节省布尔运算中相交检测的计算开销，各种布尔方法均采用对网格模型表面面片进行空间划分的预处理来提高相交检测过程的效率。目前有两种常用的相交检测模型表面空间划分的预处理方法：一种是基于空间分解的划分方法，一种是基于层次包围盒的划分方法。

第二步，根据预处理生成的空间层次结构，对网格模型进行相交性检测。

为了得到布尔运算的结果模型Ω(t)-M(t+1)，需要在网格模型Ω(t)和M(t+1)之间进行相交性检测。具体的，分为以下两个步骤：

步骤1：检测模型Ω(t)位于模型M(t+1)内的表面，并在Ω(t)中去除这部分表面，得到Ω(t)位于模型M(t+1)之外的表面模型Ω′(t)。

步骤2：检测模型M(t+1)位于模型Ω(t)内的表面，并保留这部分表面，记为M′(t+1)。

第三步，由相交检测结果分别生成当前时刻布尔运算后的结果网格模型。

首先对M′(t+1)中的所有面片进行反向操作，然后将Ω′(t)与M′(t+1)重建为一致的三角形面片网格模型，得到钻孔操作的t+1时刻的结果模型Ω(t+1)＝Ω(t)-M(t+1)。

在方法实施中，这三个步骤循环执行，实现时间的推移。在每一个视觉显示循环中得到布尔操作的当前结果，对刚性物体表面网格模型的钻孔操作进行动态仿真。

现有的对图形几何处理的基本布尔方法主要针对计算机辅助设计中模型编辑方面的应用提出，因此在实时性方面没有很高要求。在模型编辑的一次交互活动中，一般只需要进行一次空间划分的预处理，因此计算复杂度较高的预处理操作不会对方法整体造成大的影响。而在针对钻孔操作的动态三维物体布尔运算中，每一个视觉显示循环中的物体表面网格模型都会发生变化，因此需要重新进行模型空间划分的预处理操作，以便生成更新后的模型空间层次结构。另一方面，对图形几何处理的基本布尔方法第三步中需要将相交检测得到的两部分结果进行合并，重建出一致的三角形面片表面网格模型。在动态三维物体的布尔操作中同样需要在每一个视觉显示循环中进行重建工作，这也将降低整个方法的效率，影响钻孔操作的实时性。

发明内容

本发明要解决的技术问题在于针对现有基于网格模型表面空间划分布尔方法在表面网格模型进行布尔运算的几何处理，尤其是沿特定矢量方向的切割处理时预处理过程复杂和不能实时动态重建表面网格过程的不足，提出一种对刚性材质物体表面网格模型进行钻孔后模型重建的方法，将三维布尔运算问题简化为二维布尔运算问题，降低模型布尔操作中相交检测过程的计算复杂度，使得对刚性材质物体表面网格模型进行钻孔后模型重建具有实时性。

本发明的技术方案是：

对于一个刚性物体表面网格模型Ω＝(V，F)，V为网格顶点v_i的集合，v_i＝(x_i，y_i，z_i)(1≤i≤N_V)，N_V代表顶点集合V中顶点的个数，x_i，y_i，z_i为网格顶点v_i在三维直角坐标系中的三个分量。F为体表面三角形面片f_i网格集合，f_i＝(v_i1，v_i2，v_i3)(1≤i≤N_F)，N_F为F中三角形面片的个数，v_i1，v_i2，v_i3分别为第i个三角形面片按逆时针排序后的三个网格顶点。

第一步，检测与切割曲面相交的三角形面片，并将检测到的相交三角形面片曲面外网格顶点移动至切割曲面上，包括以下三个步骤：

步骤1，检测二维空间中网格顶点相对于切割曲线的位置：

将三维空间中切割曲面C以二维空间中切割曲线C′和切割矢量n形式进行表示，保持刚性物体表面网格模型Ω中顶点相对于y轴的高度参数不变，检测二维空间中三角形面片与C′是否相交。根据C′形状的不同，分别进行切割类布尔运算，若布尔运算结果为真，表示三角形面片与C′相交，若布尔运算结果为假，表示三角形面片与C′不相交：

(1)当C′为封闭的由二维平面中不共线的m(m为正整数)个顶点u_i(0≤i≤m-1)组成的凸多边形时，凸多边形的顶点逆时针排序，顶点u_i到u_i+1的边e_i＝u_i+1-u_i，(其中：u_m＝u₀)，e_i的法向量

均指向外侧。则C′所在平面中的一个点v与切割曲线C′的位置检测可由下面的条件分别判断(·表示向量的点积运算，-为向量减运算，0≤i≤m-1)：

当

表示点v在凸多边形内部，

当

表示点v在凸多边形外部，

当

且

表示点v位于凸多边形边上。

(2)当C′为椭圆形时，则点v与C′的位置检测方法是：

当

表示点在椭圆形内部，

当

表示点在椭圆形内部，

当

表示点位于椭圆形边上。

这里a、b分别为椭圆形长、短轴的半径长度，c为半焦距的长度。

步骤2，将检测到的相交三角形面片f＝(v₁，v₂，v₃)位于切割曲面C外的网格顶点移到C上，使C外的顶点投影到C上，得到新的相交三角形面片f′＝(u₁，u₂，u₃)，其中u₁，u₂，u₃分别为三角形面片的三个顶点，当v_i位于C内部时：u_i＝v_i，(i＝1，2，3)；当v_i位于C外部时：u_i为顶点v_i的投影点，即u_i为位于C上距离v_i最近的顶点。

步骤3，对移动操作后位于模型表面网格的顶点v按以下平滑处理规则进行平滑处理，平滑处理后得到点v′：

1)计算v在切割曲面C上的投影点u，在投影点u位于切割曲面C的邻域内找点v′，方法是判定(v′-u)×n＝0是否成立，若成立，将点u按切割矢量n反向移动到点v′，×表示向量的叉积运算；

2)判定(v′-v)·n_v＝0是否成立，若成立，即判定点v′和点v具有相同矢量方向n_v，n_v为v的矢量。

第二步，模型表面网格预处理，实现网格模型表面三角形面片的分类划分，方法是：

根据相对于切割曲面C及切割矢量n的位置不同，将网格模型表面的三角形面片f进行分类划分。

对刚体几何网格模型Ω＝(V，F)的矢量属性定义如下：

1、矢量(1≤i≤N_F)为三角形面片f_i＝(v_i1，v_i2，v_i3)(1≤i1，i2，i3≤N_V)的法方向：

${\stackrel{\→}{f}}_i=\stackrel{\→}{(v_{i1},v_{i2},v_{i3})}=(v_{i2}-v_{i1})\×(v_{i3}-v_{i1})---\left(1\right)$

2、矢量

(1≤i≤N_V)为网格顶点v_i的矢量方向：

${\stackrel{\→}{v}}_i=\underset{j=1}{\overset{n_s}{\mathrm{\Σ}}}{\stackrel{\→}{f}}_j/n_s,(\∀f_j,v_i\∈f_j)---\left(2\right)$

f_j是三角形面片f_j＝(v_j1，v_j2，v_j3)(1≤j1，j2，j3≤N_V)

其中，n_s为网格顶点相邻的三角形面片个数。

步骤1)，根据模型表面网格矢量相对于切割矢量n的关系，将模型的表面网格划分为两类：沿切割矢量逆方向的上表面网格集合F_i和沿切割矢量方向的下表面网格集合F_k，其中：

$F_l=\left\{f_i\right|\stackrel{\&RightArrow;}{f_i}\&CenterDot;n<0,f_i\&Element;F,1\&le;i\&le;N_F\}---\left(3\right)$

$F_k=\left\{f_i\right|\stackrel{\&RightArrow;}{f_i}\&CenterDot;n\&GreaterEqual;0,f_i\&Element;F,1\&le;i\&le;N_F\}---\left(4\right)$

步骤2)，根据布尔运算几何处理中相对于切割曲面的位置不同，应用第一步步骤1的相交检测方法，提取出布尔运算几何处理过程中切割内部区域表面网格集合F_j：

F_j＝{f_i|N(f_i，C)＜3，f_i∈F}          (5)

其中N(f_i，C)∈{0，1，2，3}，是三角形面片f_i位于C内顶点个数的统计函数；

步骤3)，根据切割矢量n所划分的模型上、下表面，将F_j进一步分为上表面切割和下表面切割内部网格集合F_ji、F_ik，即：F_j＝F_ji∪F_jk  其中：

$F_{\mathrm{ji}}=\left\{f_i\right|\stackrel{\&RightArrow;}{f_i}\&CenterDot;n<0,N(f_i,C)<3,f_i\&Element;F\}---\left(6\right)$

$F_{\mathrm{jk}}=\left\{f_i\right|\stackrel{\&RightArrow;}{f_i}\&CenterDot;n\&GreaterEqual;0,N(f_i,C)<3,f_i\&Element;F\}---\left(7\right)$

第三步，动态重建布尔运算几何处理中的网格模型。

包括以下四个步骤：

步骤(1)，存储几何网格模型Ω的高度信息，将切割矢量n旋转到与y轴正方向平行，并同向旋转初始几何网格模型Ω，得到与y轴正方向平行的切割矢量n′和旋转后的网格模型Ω′，切割方向矢量n旋转过程为：先将n围绕y轴旋转-α，使其位于xy平面上；然后将n围绕z轴旋转-β，使其与y轴正方向平行，其中α为切割矢量n＝(x_n，y_n，z_n)与x轴的夹角，

x_n，y_n，z_n为切割矢量n的三个分量，β为切割矢量n与z轴的夹角，

旋转矩阵M使矢量n变换到y轴正方向n′：nM＝n′，切割矢量n旋转至y轴正方向的矩阵M为：

M＝M_y(-α)M_z(-β)，

$M_y(-\mathrm{\&alpha;})=\left[\begin{array}{ccc}\cos (-\mathrm{\&alpha;})& 0& -\sin (-\mathrm{\&alpha;})\\ 0& 1& 0\\ \sin (-\mathrm{\&alpha;})& 0& \cos (-\mathrm{\&alpha;})\end{array}\right],$ $M_z(-\mathrm{\&beta;})=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos (-\mathrm{\&beta;})& \sin (-\mathrm{\&beta;})\\ 0& -\sin (-\mathrm{\&beta;})& \cos (-\mathrm{\&beta;})\end{array}\right]$

M_y(-α)为切割矢量n围绕y轴转动-α的旋转矩阵，M_z(-β)为切割矢量n围绕z轴转动-β的旋转矩阵；

步骤(2)，动态实时生成模型表面网格，方法是：

引入基于时间的深度控制函数h(t)表示时刻t对模型切割的深度，h(t)＝v×t，v为钻孔速率，在初始时刻切割操作的深度为h(t＝0)＝0，随时间变化的网格顶点的高度参数为：y_i(t)＝y_i(0)-h(t)，则基于时间t的模型顶点v_i(t)＝(x_i(t)，y_i(t)，z_i(t))，1≤i≤N_V具体表述为：

$v_i\left(t\right)=\left\{\begin{array}{c}(x_i\left(0\right),y_i\left(0\right),z_i\left(0\right)),T(v_i\left(0\right),C)\&GreaterEqual;0\\ (x_i\left(0\right),y_i\left(t\right),z_i\left(0\right)),T(v_i\left(0\right),C)<0\end{array}\right.---\left(8\right)$

其中，v_i(0)为钻孔过程t＝0时刻的模型顶点，T(v_i(0)，C)∈{-1，0，1}，为顶点v_i(0)与C的位置检测函数，-1表示v_i(0)位于C内部，0表示v_i(0)位于C上，1表示v_i(0)位于C外部；

基于时间t的模型表面网格f_i(t)，1≤i≤N_F为：

$f_i\left(t\right)=\left\{\begin{array}{c}(v_{i1}\left(0\right),v_{i2}\left(0\right),v_{i3}\left(0\right)),N(f_i\left(0\right),C)=0\\ (v_{i1}\left(t\right),v_{i2}\left(t\right),v_{i3}\left(t\right)),N(f_i\left(0\right),C)>0\end{array}\right.---\left(9\right)$

其中，f_i(0)＝(v_i1(0)，v_i2(0)，v_i3(0))为钻孔过程t＝0时刻的模型表面网格。

基于时间t的模型上表面集合F_i(t)和模型下表面集合F_k(t)为：

$F_i\left(t\right)=\left\{\begin{array}{cc}\left\{f_i\left(0\right)\right|\stackrel{\&RightArrow;}{f_i\left(0\right)}\&CenterDot;n<0\}& (0<t<t_0)\\ \left\{f_i\left(t\right)\right|N(f_i\left(0\right),C)=0,\stackrel{\&RightArrow;}{f_i\left(0\right)}\&CenterDot;n<0\}& (t\&GreaterEqual;t_0)\end{array}\right.---\left(10\right)$

$F_k\left(t\right)=\left\{\begin{array}{cc}\left\{f_i\left(0\right)\right|\stackrel{\&RightArrow;}{f_i\left(0\right)}\&CenterDot;n>0\}& (0<t<t_1)\\ \left\{f_i\left(t\right)\right|N(f_i\left(0\right),C)=0,\stackrel{\&RightArrow;}{f_i\left(0\right)}\&CenterDot;n>0\}& (t\&GreaterEqual;t_1)\end{array}\right.---\left(11\right)$

t₀是开始对模型沿固定矢量方向进行切割操作的时刻，t₁为切割曲面贯穿整个模型的时刻，

为三角形面片f_i(0)的法方向；

已知三维空间中切割曲面C可以由二维空间中切割曲线C′和切割矢量n形式进行表示，C′表示为XZ平面上的二维函数形式。网格模型Ω和三维空间中C的相交检测过程，可以通过同向旋转变换为网格模型Ω′和二维空间中C′的相交检测过程。由于旋转变换后切割矢量n′与y轴正方向平行，且模型Ω′由Ω经同向旋转得到，因此，网格模型Ω′和C′的相交检测与y轴的高度数值无关。

降低维数的布尔运算几何处理即是通过隐含旋转变换后Ω′和C相对于y轴的高度参数来实现的，将三维布尔运算问题转化为二维布尔运算问题。该步骤完成了隐藏模型网格顶点的第三维高度信息，将三维布尔运算简化为二维布尔运算。

步骤(3)，动态生成网格模型布尔运算几何处理的切割内壁网格，方法是：

切割内壁网格顶点集合V_b＝{v_i|T(v_i，C)＝0}，这里T(v_i，C)∈{-1，0，1}，为顶点v_i与C的位置检测函数，-1表示v_i位于C内部，0表示v_i位于C上，1表示v_i位于C外部；切割内壁网格顶点按照所位于模型表面位置不同，可进一步按(12)(13)式将V_b划分为上表面切割内壁网格顶点集V_bi和下表面切割内壁网格顶点集V_bk，即：

V_b＝V_bi∪V_bk   其中：

$V_{\mathrm{bi}}=\left\{v_i\right|T(v_i,C)=0,\stackrel{\&RightArrow;}{v_i}\&CenterDot;n<0\}---\left(12\right)$

$V_{\mathrm{bk}}=\left\{v_i\right|T(v_i,C)=0,\stackrel{\&RightArrow;}{v_i}\&CenterDot;n>0\}---\left(13\right)$

切割内壁的网格动态生成按下述步骤进行：

1)通过基于时间的高度函数h(t)，实现对切割内壁网格顶点集合的计算，来获取特定时刻t的切割内壁网格顶点集合

其中：

${V_b}^t={V_{\mathrm{bi}}}^0\&cup;{V_{\mathrm{bk}}}^t---\left(14\right)$

$V_{\mathrm{bi}}^0=\left\{v_i\left(0\right)\right|T(v_i\left(0\right),C)=0,\stackrel{\&RightArrow;}{v_i\left(0\right)}\&CenterDot;n<0\}---\left(15\right)$

${V_{\mathrm{bk}}}^t=\left\{v_i\left(t\right)\right|T(v_i\left(t\right),C)=0,\stackrel{\&RightArrow;}{v_i\left(t\right)}\&CenterDot;n>0\}---\left(16\right)$

这里矢量为网格顶点v_i(t)按公式(2)计算的矢量方向，V_bi ⁰为0时刻上表面切割内壁网格顶点集，V_bk ^t为t时刻下表面切割内壁网格顶点集；

2)采用常用的前沿推进方法，由上、下表面切割边界网格顶点生成模型布尔操作的内壁网格

$F_{\mathrm{inner}}^t=\left\{f_j\right|f_j=(v_{j1},v_{j2},v_{j3}),v_{j1},v_{j2},v_{j3}\&Element;V_b^t,1\&le;j1,j2,j3\&le;N_V\}---\left(17\right)$

步骤(4)，逆向旋转布尔运算处理后的表面网格结果模型，使该结果模型回到初始位置。

已知初始时切割矢量方向n与x轴的夹角为α、与y轴夹角为β，n′为先将n围绕y轴旋转-α，再围绕z轴旋转-β得到，相应的旋转变换矩阵为：

M＝M_y(-α)M_z(-β)                (18)

将n′旋转至初始矢量n使模型Ω′旋转回到初始位置Ω，按下述两次旋转运算完成：先将n′围绕z轴旋转β、再围绕y轴旋转α即可以得到初始矢量n。

相应的旋转变换矩阵为：

M^-1＝M_z(β)M_y(α)                     (18)

其中：

$M_z\left(\mathrm{\&beta;}\right)=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \left(\mathrm{\&beta;}\right)& \sin \left(\mathrm{\&beta;}\right)\\ 0& -\sin \left(\mathrm{\&beta;}\right)& \cos \left(\mathrm{\&beta;}\right)\end{array}\right],$ $M_y\left(\mathrm{\&alpha;}\right)=\left[\begin{array}{ccc}\cos \left(\mathrm{\&alpha;}\right)& 0& -\sin \left(\mathrm{\&alpha;}\right)\\ 0& 1& 0\\ \sin \left(\mathrm{\&alpha;}\right)& 0& \cos \left(\mathrm{\&alpha;}\right)\end{array}\right]$

将模型的所有表面三角面片f_i和网格顶点v_i都乘以逆向旋转矩阵M^-1，得到逆向旋转布尔运算几何处理的结果网格模型Ω′。

步骤(5)将旋转运算之后几何模型顶点的高度值修改为第三步步骤(1)的旋转运算之前保存的几何模型顶点的高度值，得到特定t时刻对模型进行布尔运算几何处理的结果。

在布尔操作的动态处理过程中，仅需要考虑切割钻孔的深度变化，而不必在每一次网格重建中相对最新状态切割曲面进行额外的初始化操作，简化了布尔运算的预处理过程；同时，通过对初始模型的旋转变换和对结果模型的逆向旋转，隐含了在布尔运算中对几何网格模型第三维高度参数的操作，将三维布尔运算问题简化为二维布尔运算的降维问题，降低了模型布尔操作中相交检测过程的运算复杂度。

采用本发明可获得如下有益效果：

第一，本发明直接采用二维布尔运算处理方法实现三维空间中对几何网格模型进行切割、钻孔类的布尔运算几何处理，与三维空间中布尔方法相比，省略了对模型表面进行空间划分的过程，且简化了相交检测算法，极大提高了系统运算效率，达到实时处理的效果；

第二，本发明采用基于保持网格顶点矢量的表面平滑方法，代替了三维布尔方法中相关的表面平滑等后处理过程，使得每次布尔运算处理的几何网格模型都能够保持平滑的表面及锋锐的切割边界特效；

第三，本发明采用基于时间的深度控制函数，能在视觉的每一帧快速渲染布尔运算几何处理的效果，避免了传统三维布尔方法中每次处理所必然引入的预处理过程，真正实现针对几何网格模型的快速、动态的布尔运算处理效果。

附图说明

图1为本发明总体流程图；

图2为布尔操作中网格模型表面面片划分示意图；

图3为布尔操作中网格顶点平滑处理示意图；

具体实施方式

采用3DS Max软件生成的VRML格式表面网格模型作为输入，弹性物体的表面网格模型表示为Ω＝(V，F)，其中V为网格顶点集合，F为三角形面片集合。如图1所示，本发明具体步骤为：

第一步，检测与切割曲面相交的三角形面片，并将检测到的相交三角形面片曲面外网格顶点移动至切割曲面上。

第二步，模型表面网格预处理，实现网格模型表面三角形面片的分类划分。

第三步，动态重建布尔运算几何处理中的网格模型，包括以下五个步骤：

步骤(1)，沿切割矢量方向逆向旋转网格模型，使切割矢量和固定轴方向一致；

步骤(2)，动态生成模型表面网格；

步骤(3)，动态生成网格模型布尔运算几何处理的切割内壁网格；

步骤(4)，逆向旋转布尔运算处理后的表面网格结果模型，使该结果模型回到初始位置；

步骤(5)将旋转运算之后几何模型顶点的高度值修改为第三步步骤(1)的旋转运算之前保存的几何模型顶点的高度值，得到特定t时刻对模型进行布尔运算几何处理的结果。

图2是模型表面网格划分示意图，Ω为表面网格模型，F_i为沿切割矢量逆方向的上表面网格集合，F_k为沿切割矢量方向的下表面网格集合，n为切割矢量，F_ji、F_ik分别为上表面切割和下表面切割内部网格集合。

图3是模型表面网格的顶点模型表面网格的顶点v平滑处理后得到点v′：的示意图，n为切割矢量，n_v为v的矢量方向。平滑处理的过程是：

1)计算v在切割曲面C上的投影点u，在投影点u位于切割曲面C的邻域内找点v′，方法是判定(v′-u)×n＝0是否成立，若成立，将点u按切割矢量n反向移动到点v′，×表示向量的叉积运算；

2)判定(v′-v)·n_v＝0是否成立，若成立，即判定点v′和点v具有相同矢量方向n_v，n_v为v的矢量。

Claims (1)

1.一种对刚性材质物体表面网格模型进行钻孔后模型重建的方法，其特征在于包括以下步骤：

第一步，检测与切割曲面相交的三角形面片，并将检测到的相交三角形面片曲面外网格顶点移动至切割曲面上，包括以下三个步骤：

步骤1，检测二维空间中网格顶点相对于切割曲线的位置：

将三维空间中切割曲面C以二维空间中切割曲线C′和切割矢量n形式进行表示，保持刚性物体表面网格模型Ω中顶点相对于y轴的高度参数不变，检测二维空间中三角形面片与C′是否相交，根据C′形状的不同，分别进行切割类布尔运算，若布尔运算结果为真，表示三角形面片与C′相交，若布尔运算结果为假，表示三角形面片与C′不相交，刚性物体表面网格模型Ω＝(V，F)，V为网格顶点v_i的集合，v_i＝(x_i，y_i，z_i)，1≤i≤N_V，N_V为网格顶点集合V中顶点的个数，x_i，y_i，z_i为顶点v_i在三维直角坐标系中的三个分量，F为体表面三角形面片f_i网格集合，f_i＝(v_i1，v_i2，v_i3)，1≤i≤N_F，N_F为F中三角形面片的个数，v_i1，v_i2，v_i3分别为第i个三角形面片按逆时针排序后的三个网格顶点：

(1)当C′为封闭的由二维平面中不共线的m个顶点u_i组成的凸多边形时，凸多边形的顶点逆时针排序，顶点u_i到u_i+1的边e_i＝u_i+1-u_i，e_i的法向量

均指向外侧，其中，m为正整数，0≤i≤m-1，u_m＝u₀，则C′所在平面中的一个点v与切割曲线C′的位置检测由下面的条件分别判断：

当

表示点v在凸多边形内部，

当

表示点v在凸多边形外部，

当

且

表示点v位于凸多边形边上；其中，·表示向量的点积运算，-为向量减运算；

(2)当C′为椭圆形时，则点v与C′的位置检测方法是：

当

表示点在椭圆形内部，

当

表示点在椭圆形内部，

当

表示点位于椭圆形边上，a、b分别为椭圆长、短轴的半径长度，c为半焦距的长度；

步骤2，将检测到的相交三角形面片f＝(v₁，v₂，v₃)位于切割曲面C外的网格顶点移到C上，使C外的顶点投影到C上，得到新的相交三角形面片f′＝(u₁，u₂，u₃)，其中u₁，u₂，u₃分别为三角形面片的三个顶点，当v_i位于C内部时：u_i＝v_i，i＝1，2，3；当v_i位于C外部时：u_i为顶点v_i的投影点，即u_i为位于C上距离v_i最近的顶点；

步骤3，对移动操作后位于模型表面网格的顶点v按以下平滑处理规则进行平滑处理，平滑处理后得到点v′：

1)计算v在切割曲面C上的投影点u，在投影点u位于切割曲面C的邻域内找点v′，方法是判定(v′-u)×n＝0是否成立，若成立，将点u按切割矢量n反向移动到点v′，×表示向量的叉积运算；

2)判定(v′-v)·n_v＝0是否成立，若成立，即判定点v′和点v具有相同矢量方向n_v，n_v为v的矢量；

第二步，模型表面网格预处理，实现网格模型表面三角形面片的分类划分，方法是：

步骤1)，根据模型表面网格矢量

相对于切割矢量n的关系，将模型的表面网格划分为两类：沿切割矢量逆方向的上表面网格集合F_l和沿切割矢量方向的下表面网格集合F_k，其中：

$F_l=\left\{f_i\right|\stackrel{\&RightArrow;}{f_i}\&CenterDot;n<0,f_i\&Element;F,1\&le;i\&le;N_F\}---\left(3\right)$

$F_k=\left\{f_i\right|\stackrel{\&RightArrow;}{f_i}\&CenterDot;n\&GreaterEqual;0,f_i\&Element;F,1\&le;i\&le;N_F\}---\left(4\right)$

其中，矢量

为三角形面片f_i＝(v_i1，v_i2，v_i3)的法方向，1≤i≤N_F，1≤i1，i2，i3≤N_V；

${\stackrel{\&RightArrow;}{f}}_i=\stackrel{\&RightArrow;}{(v_{i1},v_{i2},v_{i3})}=(v_{i2}-v_{i1})\&times;(v_{i3}-v_{i1})---\left(1\right)$

步骤2)，根据布尔运算几何处理中相对于切割曲面的位置不同，应用第一步步骤1的相交检测方法，提取出布尔运算几何处理过程中切割内部区域表面网格集合F_j：

F_j＝{f_i|N(f_i，C)＜3，f_i∈F}                    (5)

其中N(f_i，C)∈{0，1，2，3}，是三角形面片f_i位于C内顶点个数的统计函数；

步骤3)，根据切割矢量n所划分的模型上、下表面，将F_j进一步分为与上表面切割和下表面切割内部网格集合F_ji、F_jk，即：

F_j＝F_ji∪F_jk其中：

$F_{\mathrm{ji}}=\left\{f_i\right|\stackrel{\&RightArrow;}{f_i}\&CenterDot;n<0,N(f_i,C)<3,f_i\&Element;F\}---\left(6\right)$

$F_{\mathrm{jk}}=\left\{f_i\right|\stackrel{\&RightArrow;}{f_i}\&CenterDot;n\&GreaterEqual;0,N(f_i,C)<3,f_i\&Element;F\}---\left(7\right)$

第三步，动态重建布尔运算几何处理中的网格模型，方法是：

步骤(1)，存储几何网格模型Ω的高度信息，将切割矢量n旋转到与y轴正方向平行，并同向旋转初始几何网格模型Ω，得到与y轴正方向平行的切割矢量n′和旋转后的网格模型Ω′，切割方向矢量n旋转过程为：先将n围绕y轴旋转-α，使其位于xy平面上；然后将n围绕z轴旋转-β，使其与y轴正方向平行，其中α为切割矢量n＝(x_n，y_n，z_n)与x轴的夹角，

x_n，y_n，z_n为切割矢量n的三个分量，β为切割矢量n与z轴的夹角，

旋转矩阵M使矢量n变换到y轴正方向n′：nM＝n′，切割矢量n旋转至y轴正方向的矩阵M为：

M＝M_y(-α)M_z(-β)，

$M_y(-\mathrm{\&alpha;})=\left[\begin{array}{ccc}\cos (-\mathrm{\&alpha;})& 0& -\sin (-\mathrm{\&alpha;})\\ 0& 1& 0\\ \sin (-\mathrm{\&alpha;})& 0& \cos (-\mathrm{\&alpha;})\end{array}\right],$ $M_z(-\mathrm{\&beta;})=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos (-\mathrm{\&beta;})& \sin (-\mathrm{\&beta;})\\ 0& -\sin (-\mathrm{\&beta;})& \cos (-\mathrm{\&beta;})\end{array}\right]$

M_y(-α)为切割矢量n围绕y轴转动-α的旋转矩阵，M_z(-β)为切割矢量n围绕z轴转动-β的旋转矩阵；

步骤(2)，动态实时生成模型表面网格，方法是：引入基于时间的深度控制函数h(t)表示时刻t对模型切割的深度，h(t)＝v×t，v为钻孔速率，在初始时刻切割操作的深度为h(t＝0)＝0，随时间变化的网格顶点的高度参数为：y_i(t)＝y_i(0)-h(t)，则基于时间t的模型顶点v_i(t)＝(x_i(t)，y_i(t)，z_i(t))，1≤i≤N_V表述为：

$v_i\left(t\right)=\left\{\begin{array}{c}(x_i\left(0\right),y_i\left(0\right),z_i\left(0\right)),T(v_i\left(0\right),C)\&GreaterEqual;0\\ (x_i\left(0\right),y_i\left(t\right),z_i\left(0\right)),T(v_i\left(0\right),C)<0\end{array}\right.---\left(8\right)$

其中，v_i(0)为钻孔过程t＝0时刻的模型顶点，T(v_i(0)，C)∈{-1，0，1}，为顶点v_i(0)与C的位置检测函数，-1表示v_i(0)位于C内部，0表示v_i(0)位于C上，1表示v_i(0)位于C外部；

基于时间t的模型表面网格f_i(t)，1≤i≤N_F为：

$f_i\left(t\right)=\left\{\begin{array}{c}(v_{i1}\left(0\right),v_{i2}\left(0\right),v_{i3}\left(0\right)),N(f_i\left(0\right),C)=0\\ (v_{i1}\left(t\right),v_{i2}\left(t\right),v_{i3}\left(t\right)),N(f_i\left(0\right),C)>0\end{array}\right.---\left(9\right)$

其中，f_i(0)＝(v_i1(0)，v_i2(0)，v_i3(0))为钻孔过程t＝0时刻的模型表面网格；

基于时间t的模型上表面集合F_i(t)和模型下表面集合F_k(t)为：

$F_i\left(t\right)=\left\{\begin{array}{cc}\left\{f_i\left(0\right)\right|\stackrel{\&RightArrow;}{f_i\left(0\right)}\&CenterDot;n<0\}& (0<t<t_0)\\ \left\{f_i\left(t\right)\right|N(f_i\left(0\right),C)=0,\stackrel{\&RightArrow;}{f_i\left(0\right)}\&CenterDot;n<0\}& (t\&GreaterEqual;t_0)\end{array}\right.---\left(10\right)$

$F_k\left(t\right)=\left\{\begin{array}{cc}\left\{f_i\left(0\right)\right|\stackrel{\&RightArrow;}{f_i\left(0\right)}\&CenterDot;n>0\}& (0<t<t_1)\\ \left\{f_i\left(t\right)\right|N(f_i\left(0\right),C)=0,\stackrel{\&RightArrow;}{f_i\left(0\right)}\&CenterDot;n>0\}& (t\&GreaterEqual;t_1)\end{array}\right.---\left(11\right)$

t₀是开始对模型沿固定矢量方向进行切割操作的时刻，t₁为切割曲面贯穿整个模型的时刻，

为三角形面片f_i(0)的法方向；

步骤(3)，动态生成网格模型布尔运算几何处理的切割内壁网格，方法是：

1)通过基于时间的高度函数h(t)，实现对切割内壁网格顶点集合的计算，来获取时刻t的切割内壁表面网格顶点集合

其中：

${V_b}^t={V_{\mathrm{bi}}}^0\&cup;{V_{\mathrm{bk}}}^t---\left(14\right)$

$V_{\mathrm{bi}}^0=\left\{v_i\left(0\right)\right|T(v_i\left(0\right),C)=0,\stackrel{\&RightArrow;}{v_i\left(0\right)}\&CenterDot;n<0\}---\left(15\right)$

${V_{\mathrm{bk}}}^t=\left\{v_i\left(t\right)\right|T(v_i\left(t\right),C)=0,\stackrel{\&RightArrow;}{v_i\left(t\right)}\&CenterDot;n>0\}---\left(16\right)$

这里矢量

为网格顶点v_i(t)按公式(2)计算的矢量方向；

V_b为切割内壁网格顶点集合，V_b＝{v_i|T(v_i，C)＝0}，这里T(v_i，C)∈{-1，0，1}，为顶点v_i与C的位置检测函数，-1表示v_i位于C内部，0表示v_i位于C上，1表示v_i位于C外部；

为时刻t的切割内壁表面网格顶点集合；切割内壁网格顶点按照所位于模型表面位置不同，可进一步按(12)(13)式将V_b划分为上表面切割内壁网格顶点集V_bi和下表面切割内壁网格顶点集V_bk，即：

V_b＝V_bi∪V_bk  其中：

$V_{\mathrm{bi}}=\left\{v_i\right|T(v_i,C)=0,\stackrel{\&RightArrow;}{v_i}\&CenterDot;n<0\}---\left(12\right)$

$V_{\mathrm{bk}}=\left\{v_i\right|T(v_i,C)=0,\stackrel{\&RightArrow;}{v_i}\&CenterDot;n>0\}---\left(13\right)$

其中矢量

为网格顶点v_i的矢量方向，1≤i≤N_v：

${\stackrel{\&RightArrow;}{v}}_i=\underset{j=1}{\overset{n_s}{\mathrm{\&Sigma;}}}{\stackrel{\&RightArrow;}{f}}_j/n_s,(\&ForAll;f_j,v_i\&Element;f_j)---\left(2\right)$

其中，n_s为网格顶点相邻的三角形面片个数；f_j是三角形面片，f_j＝(v_j1，v_j2，v_j3)，1≤j1，j2，j3≤N_V；V_bi ⁰为0时刻上表面切割内壁网格顶点集，V_bk ^t为t时刻下表面切割内壁网格顶点集；

2)采用常用的前沿推进方法，由上、下表面切割内壁网格顶点生成模型布尔操作的内壁网格

$F_{\mathrm{inner}}^t=\left\{f_j\right|f_j=(v_{j1},v_{j2},v_{j3}),v_{j1},v_{j2},v_{j3}\&Element;V_b^t,1\&le;j1,j2,j3\&le;N_V\}---\left(17\right)$

步骤(4)，逆向旋转布尔运算处理后的表面网格结果模型，使该结果模型回到初始位置，即将n′旋转至初始矢量n，使模型Ω′旋转回到初始位置Ω，按下述两次旋转运算完成：先将n′围绕z轴旋转β、再围绕y轴旋转α得到初始矢量n，

相应的旋转变换矩阵为：

M^-1＝M_z(β)M_y(α)                (18)

其中：

$M_z\left(\mathrm{\&beta;}\right)=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \left(\mathrm{\&beta;}\right)& \sin \left(\mathrm{\&beta;}\right)\\ 0& -\sin \left(\mathrm{\&beta;}\right)& \cos \left(\mathrm{\&beta;}\right)\end{array}\right],$ $M_y\left(\mathrm{\&alpha;}\right)=\left[\begin{array}{ccc}\cos \left(\mathrm{\&alpha;}\right)& 0& -\sin \left(\mathrm{\&alpha;}\right)\\ 0& 1& 0\\ \sin \left(\mathrm{\&alpha;}\right)& 0& \cos \left(\mathrm{\&alpha;}\right)\end{array}\right]$

将模型的所有表面三角面片f_i和网格顶点v_i都乘以逆向旋转矩阵M^-1，得到逆向旋转布尔运算几何处理的结果网格模型Ω′；

步骤(5)将旋转运算之后几何模型顶点的高度值修改为第三步步骤(1)的旋转运算之前保存的几何模型顶点的高度值，得到特定t时刻对模型进行布尔运算几何处理的结果。

Images