CN104157013A - 任意四边形网格上光滑曲面的重建方法 - Google Patents

Abstract

任意四边形网格上光滑曲面的重建方法，包括以下步骤：1)生成四边形区域边界曲线；2)对每个四边形区域，由边界曲线和内部数据点拟合生成B样条曲面；3)在正则顶点处进行C²相容性调整，奇异顶点处G¹相容性调整；4)正则边界上进行C²连续调整；5)非正则边界上进行近似G¹连续调整。本发明提供了一种简化求解、提高效率、提升重建效果的任意四边形网格上光滑B样条曲面的重建方法。

Description

任意四边形网格上光滑曲面的重建方法

技术领域

本发明属于曲面造型领域，涉及一种在经过四边形划分后网格上重建光滑拼接曲面片的方法。

背景技术

在对诸如人体器官、骨骼、玩具公仔等以自由曲面为主的复杂模型逆向工程建模中，很难用单一曲面片对模型进行表述，通常首先需要对原始模型进行四边形划分，然后在每个四边形区域上进行曲面重建，并逼近内部数据点，同时曲面片之间应该实现光滑拼接。衡量光滑程度通常用两种连续性：一是参数连续性，记为Cⁿ；二是几何连续性，记为Gⁿ，n指连续阶数。参数连续是指在曲面切平面的任一方向上两曲面片的方向导数相等，Cⁿ连续指0-n阶方向导数相等。几何连续中的G⁰是指位置连续，G¹指切平面连续，G²指主曲率和主方向连续。

四边形网格的顶点根据顶点的度数(与顶点相连的边的条数)分为正则顶点和奇异顶点：网格内部的正则顶点度数为4，奇异顶点度数为3、5及以上；网格边界上的正则顶点度数为3，奇异顶点度数为2、4及以上。如果一条边的两个顶点都是正则顶点，则该边为正则边，否则为非正则边。

多张曲面间的连续拼接问题，已有不少学者对此进行了研究，出现了多种连续拼接曲面的重建方法。Milroy等[1]提出了一种近似G¹连续的分片B样条曲面拟合方法，将G¹连续条件作为约束条件，通过求解非线性最小二乘问题的解来保证曲面间的连续性。该方法计算量比较大，而且只能得到近似的G¹连续。Eck和Hoppe^[2]利用Peter^[3]提出的曲面样条构造方法来生成双三次B样条曲面片，对四边形网格进行两次Doo-Sabin细分^[5]后得到一个生成B样条曲面片的控制网格，B样条曲面片的控制顶点由控制网格顶点的线性组合得到，构造过程自动保证了B样条曲面片间的G¹连续。曲面的拟合转化为对控制网格顶点位置的最小二乘求解。为了得到保证精度的曲面，曲面构造需要迭代进行，整个过程非常耗时。Peter^[4]通过对任意四边形网格进行Catmull-Clark细分^[6]，在每个四边形面片上，得到一张双三次B样条曲面，并且曲面片间除了奇异点处C¹连续，其他各处达到了C²连续，但由于Catmull-Clark细分的收缩性，重建的曲面模型会有较大的收缩。施锡泉等在[7]中给出了在经过四边形区域划分后的模型上重建收敛G¹连续的双三次B样条曲面片的方法，随着曲面控制顶点的增加，曲面片间趋向于G¹连续。随后，在[8]中讨论了内部没有重节点的双五次B样条曲面片间的G¹连续性条件，给出了生成严格的G¹连续的双五次B样条曲面的局部构造方法。局部构造方法首先拟合得到一个G⁰连续的曲面片网格，然后再根据G¹连续性条件对曲面的控制顶点进行局部修改使曲面片间连续。

已有的分片连续曲面重建方法得到的曲面大部分只有G¹连续，而且有的方法求解复杂，效率不高，或者重建的曲面不能较好的保持原始模型的几何特征，并不实用。

发明内容

为了克服已有分片连续曲面重建方法的求解复杂、效率不高、重建效果较差的不足,本发明提供了一种简化求解、提高效率、提升重建效果的任意四边形网格上光滑B样条曲面的重建方法。

本发明解决其技术问题所采用的技术方案是：

任意四边形网格上光滑曲面的重建方法，所述重建方法包括下述步骤：

1)生成四边形区域边界曲线

用三次B样条曲线拟合四边形区域边界上的采样点得到边界曲线，而且对所有的边界曲线采用相同的节点向量，使相对的两条边的节点向量相容；

在生成边界曲线时，首先找出四边形网格上的长边，长边由依次相连的多条边组成，长边的首末端点为奇异顶点或网格边界上的正则顶点；得到长边后，取出组成长边的每条边上的采样点，将它们合在一起，拟合出一条长曲线，该曲线插值长边上的顶点，逼近每条边上的采样点；再将长曲线在顶点处打断，得到每条边上对应的边界曲线，使边界曲线在正则顶点处C²连续，奇异顶点处G⁰连续；

2)对每个四边形区域，由边界曲线和内部数据点拟合生成B样条曲面；

3)在正则顶点处进行C²相容性调整，再进行奇异顶点处G¹相容性调整，过程如下：

3.1)正则顶点处的C²相容性调整按下述方法进行：

计算与正则顶点P相连的4个曲面{S_i(u_i,v_i),i＝1,...,4}在P处的混合偏导 $V_i=\frac{{\∂}^2{S}_i(u_i,v_i)}{{\∂u}_i{\∂v}_i}{|}_{u_i=v_i=0},X_i=\frac{{\∂}^3{S}_i(u_i,v_i)}{{\∂u}_i^2{\∂v}_i}{|}_{u_i=v_i=0},Y_i=\frac{{\∂}^3{S}_i(u_i,v_i)}{{\∂u}_i{\∂v}_i^2}{|}_{u_i=v_i=0},Z_i=\frac{{\∂}^4{S}_i(u_i,v_i)}{{\∂u}_i^2{\∂v}_i^2}{|}_{u_i=v_i=0};$ 对它们进行平均，将平均后的混合偏导作为曲面在P点处的新的混合偏导，并调整曲面相应的控制顶点，则曲面在P点处有公共的2至4阶混合偏导，在P点处C²相容；

3.2)奇异顶点处G¹相容性调整按下述方法进行：

记与奇异顶点P相连的曲面片为{S_i(u_i,v_i),i＝1,...,n}，边界曲线为{C_i,i＝1,...,n}；要使G¹连续的曲面片{S_i(u_i,v_i),i＝1,...,n}在顶点P处扭矢相容，要求边界曲线{C_i}在顶点P处达到G²连续；在顶点P点处构造一张C²连续的基曲面S，其中S在P点的一阶、二阶导数分别为S_u、S_v、S_uu、S_uv、S_vv；边界曲线{C_i,i＝1,...,n}在P点处G²连续，则P点处新的切矢和曲率矢满足下述关系式：

$T_i^{*}={\mathrm{\α}}_i{S}_u+{\mathrm{\β}}_i{S}_v,$

$H_i^{*}={\mathrm{\α}}_i^2{S}_{\mathrm{uu}}+{2\mathrm{\α}}_i{\mathrm{\β}}_i{S}_{\mathrm{uv}}+{\mathrm{\β}}_i^2{S}_{\mathrm{vv}}+{\mathrm{\λ}}_i{S}_u+{\mathrm{\γ}}_i{S}_v,$

其中α_i,β_i,λ_i,γ_i为常数。将曲面片{S_i(u_i,v_i),i＝1,...,n}在P点处的扭矢取为

$V_i^{*}={\mathrm{\α}}_{i+1}{\mathrm{\α}}_i{S}_{\mathrm{uu}}+({\mathrm{\α}}_{i+1}{\mathrm{\β}}_i+{\mathrm{\α}}_i{\mathrm{\β}}_{i+1})S_{\mathrm{uv}}+{\mathrm{\β}}_{i+1}{\mathrm{\β}}_i{S}_{\mathrm{vv}}+{\mathrm{\ξ}}_i{S}_u+{\mathrm{\η}}_i{S}_v;$

则{S_i(u_i,v_i),i＝1,...,n}在P点处G¹相容；

首先在P点处求出一个公共切平面∏，将T_i投影到∏，得到新的切矢然后用最小二乘法求出曲率矢和扭矢在法向N上的投影调整后的曲率矢和扭矢为

$H_i^{*}=H_i+(<H_i^{*},N>-<H_i,N>)N,V_i^{*}=V_i+(<V_i^{*},N>-<V_i,N>)N.$

求得新的曲率矢和扭矢后，调整曲面片的相应控制顶点，使{S_i,i＝1,...,n}在顶点P处G²连续并且G1相容；

4)正则边界上进行C²连续调整

计算与正则边C_i(v)相连的两曲面S_i-1、S_i的一阶、二阶跨界导矢曲线 对它们进行平均，得到新的跨界导矢曲线 $D_i^{+*}\left(v\right),{\mathrm{DD}}_i^{+*}\left(v\right):$

$D_i^{-*}\left(v\right)=D_i^{+*}\left(v\right)=\frac{1}{2}(D_i^{-}\left(v\right)+D_i^{+}\left(v\right)),{\mathrm{DD}}_i^{-*}\left(v\right)={\mathrm{DD}}_i^{+*}\left(v\right)=\frac{1}{2}({\mathrm{DD}}_i^{-}\left(v\right)+{\mathrm{DD}}_i^{+}\left(v\right)).$

再根据新的跨界导矢曲线调整曲面S_i-1、S_i靠近边界C_i(v)的第二、第三排控制顶点，使S_i-1、S_i在边界C_i(v)上C²连续；

5)非正则边界上进行近似G¹连续调整

首先计算两曲面在公共边界上内部节点处的一阶跨界导矢 根据切平面矫正调整这些跨界导矢，然后求出新的跨界导矢曲线使插值 $\{D_i^{-}\left(v_k\right),k=4,...,n\},\{D_i^{+}\left(v_k\right),k=4,...,n\},$ 且满足如下端点条件：

$D_i^{-}\left(0\right)=T_{i-\mathrm{1,0}},D_i^{-}\left(1\right)=T_{i-\mathrm{1,1}},D_i^{{-}^{\&prime;}}\left(0\right)=V_{i-\mathrm{1,0}},D_i^{{-}^{\&prime;}}\left(1\right)=V_{i-\mathrm{1,1}};$

$D_i^{+}\left(0\right)=T_{i+\mathrm{1,0}},D_i^{+}\left(1\right)=T_{i+\mathrm{1,1}},D_i^{{+}^{\&prime;}}\left(0\right)=V_{i,0},D_i^{{+}^{\&prime;}}\left(1\right)=V_{i,1}.$

接着由新的跨界导矢曲线调整两个曲面靠近公共边界的第二排控制顶点，检查两曲面在公共边界上法矢的夹角，在夹角大于给定的误差允许值的地方插入节点，重新计算跨界导矢曲线；整个过程迭代进行，直到各个节点区间的法矢夹角都小于给定的误差允许值或插入的节点数量达到一个上限值。

进一步，所述步骤2)中，给定四边形区域的四条边界曲线及内部数据点集，通过最小二乘拟合重构一张双三次B样条曲面，重构曲面插值四条边界曲线，逼近内部点。

再进一步，所述步骤2)中，首先由4条边界曲线创建一张双线性Coons曲面，然后将散乱点P_k投影到该基曲面计算得到点P_k的参数值(u_k,v_k)；曲面S(u,v)边界上的控制顶点由4条边界曲线的控制顶点确定，内部的控制顶点通过最小二乘拟合如下的目标函数F得到，

F＝F_lsq+λF_fair，其中F_lsq为最小二乘项，F_fair为光顺项，λ为光顺项系数。

更进一步，所述步骤5)中，具体步骤如下：

Step0:设定初始值：两个曲面在公共边界上法矢夹角θ的最大误差为θ_max，在公共边界上允许插入的最大节点个数为N_max，当前插入的节点个数N_inserted＝0，边界曲线初始的节点向量为V＝{0,0,0,0,v₄,v₅,...,v_n,1,1,1,1}；

Step1:在节点{v_k,k＝4,...,n}处，计算曲线C_i(v)的切矢{C′_i(v_k),k＝4,...,n}，S_i-1、S_i沿着边界C_i(v)的跨界导矢S_i-1、S_i在点C_i(v_k)处的法矢为 $N_k^{-}=D_i^{-}\left(v_k\right)\&times;C_i^{\&prime;}\left(v_k\right),N_k^{+}=C_i^{\&prime;}\left(v_k\right)\&times;D_i^{+}\left(v_k\right)$ (k＝4,...,n)；

Step2:矫正采样点C_i(v_k)处的跨界导矢，使共面，满足G¹连续条件。首先确定两个曲面在点C_i(v_k)处的公共法向

$N_k=(N_k^{-}+N_k^{+})/\left|\right|N_k^{-}+N_k^{+}\left|\right|;$

矫正后的跨界导矢为

$D_i^{-}{\left(v_k\right)}^{*}=D_i^{-}\left(v_k\right)-<D_i^{-}\left(v_k\right),N_k>*N_k,$

$D_i^{+}{\left(v_k\right)}^{*}=D_i^{+}\left(v_k\right)-<D_i^{+}\left(v_k\right),N_k>*N_k.$

Step3:以当前边界曲线C_i(v)的节点向量V作为跨界导矢曲线的节点向量，拟合新的跨界导矢曲线使插值 $D^{-}\{T_{i-\mathrm{1,0}},D_i^{-}{\left(v_4\right)}^{*},...,D_i^{-}{\left(v_n\right)}^{*},T_{i-\mathrm{1,1}}\}$ 且 $D_i^{-\&prime;}\left(0\right)=V_{i-\mathrm{1,0}},D_i^{-\&prime;}\left(1\right)=V_{i-\mathrm{1,1}},$ 插值 $D^{+}\{T_{i+\mathrm{1,0}},D_i^{+}{\left(v_4\right)}^{*},...,D_i^{+}{\left(v_n\right)}^{*},T_{i+\mathrm{1,1}}\}$ 且 $D_i^{+\&prime;}\left(0\right)=V_{i,0},D_i^{+\&prime;}\left(1\right)=V_{i,1};$

Step4:根据矫正后的跨界导矢曲线调整S_i-1、S_i靠近边界C_i(v)的第二排控制顶点；

Step5:经过上述调整，S_i-1、S_i在节点{v_k,k＝3,...,n+1}处已经达到G¹连续，在每个节点区间[v_k,v_k+1](k＝3,...,n)内再取m个采样点

$v_{k,j}=v_k+\frac{j}{m+1}(v_{k+1}-t_k),(j=1,...,m);$

计算S_i-1、S_i在节点{v_k,j,k＝3,...,n,j＝1,...,m}处的跨界导矢法矢求出两个法向的夹角{θ_k,j}；

Step6:找出每个区间内夹角最大处的节点值{v_k,j,k＝3,...,n,j＝p_k}，如果该处的夹角大于θ_max，则将加入到待插入的节点集合,记待插入节点集合为 $V^{*}=\{v_l^{*},l=0,...,s\};$

Step7:如果，将节点分别插入到曲线C_i(v)及曲面S_i-1、S_i中去，得到新的节点向量V＝V∪V^*，计算当前已插入的节点数N_inserted＝N_inserted+s+1；矫正处的跨界导矢，得到新的跨界导矢 $D^{-*}=\{D_i^{-}\left(v_l^{*}\right),l=0,...,s\},D^{+*}=\{D_i^{+}\left(v_l^{*}\right),l=0,...,s\},$ 新的待拟合的跨界导矢集合为D^-＝D^-∪D^-*、D⁺＝D⁺∪D^+*；

Step8:重复Step3－Step7，直到或者N_inserted＞N_max，算法结束。

本发明的技术构思为：在任意拓扑四边形网格上重建G¹/C²连续双三次B样条曲面片的方法。四边形网格的顶点根据顶点的度数(与顶点相连的边的条数)分为正则顶点和奇异顶点：网格内部的正则顶点度数为4，奇异顶点度数为3、5及以上；网格边界上的正则顶点度数为3，奇异顶点度数为2、4及以上。如果一条边的两个顶点都是正则顶点，则该边为正则边，否则为非正则边。整个算法包括下述步骤：1)生成四边形区域边界曲线；2)对每个四边形区域，由边界曲线和内部数据点拟合生成B样条曲面；3)在正则顶点处进行C²相容性调整，奇异顶点处G¹相容性调整；4)正则边界上进行C²连续调整；5)非正则边界上进行近似G¹连续调整。

本发明的有益效果主要表现在：1、简化求解、提高效率；2、得到的曲面片具有较好的连续性，并能保持原始模型的几何特征。

附图说明

图1是正则顶点处相连的曲面图。

图2是奇异顶点处相连的曲面图。

图3是公共边界为正则边的两张曲面图。

图4是公共边界为非正则边的两张曲面图。

图5是封闭网格模型曲面重建示例，其中，(a)是经过四边形划分的牙齿模型，(b)是重建的光滑曲面模型。

图6是非封闭网格模型曲面重建示例，其中，(a)是经过四边形划分的米老鼠模型，(b)是重建的光滑曲面模型。

具体实施方式

下面结合附图对本发明作进一步描述。

参照图1～图6，任意四边形网格上光滑曲面的重建方法，包括以下步骤：

1)四边形区域边界曲线的生成

用三次B样条曲线拟合四边形区域边界上的采样点得到边界曲线，而且对所有的边界曲线采用相同的节点向量，使相对的两条边的节点向量相容，以减少整个曲面网格的控制顶点数目。在生成曲线时，首先找出四边形网格上的长边，长边由依次相连的多条边组成，长边的首末端点为奇异顶点或网格边界上的正则顶点。得到长边后，取出组成长边的每条边上的采样点，将它们合在一起，拟合出一条长曲线，该曲线插值长边上的顶点，逼近每条边上的采样点。再将长曲线在顶点处打断，得到每条边上对应的边界曲线，使边界曲线在正则顶点处C²连续，奇异顶点处G⁰连续。

2)B样条曲面片的生成

给定四边形区域的4条边界曲线{C_i,i＝1,...,4}及内部数据点集{P_k,k＝1,...,l}，通过最小二乘拟合重构一张双三次B样条曲面S(u,v)，使重构曲面在逼近{P_k,k＝1,...,l}的同时，插值4条边界曲线。

散乱数据点参数化采用经典的基面投影法，首先由4条边界曲线创建一张双线性Coons曲面，然后将散乱点P_k投影到该基曲面计算得到点P_k的参数值(u_k,v_k)。曲面S(u,v)边界上的控制顶点由4条边界曲线的控制顶点确定，内部的控制顶点通过最小二乘拟合如下的目标函数F得到，

F＝F_lsq+λF_fair，其中F_lsq为最小二乘项，F_fair为光顺项，λ为光顺项系数。

3)正则顶点处的C²相容性矫正，奇异顶点处的相容性调整；

如图1所示，P是正则顶点，与P相连的4条边界曲线{C_i,i＝1,...,4}两两C²连续，与P相连的4个曲面{S_i(u_i,v_i),i＝1,...,4}具有一致的u、v方向。计算4个曲面{S_i,i＝1,...,4}在P处的混合偏导 $V_i=\frac{{\&PartialD;}^2{S}_i(u_i,v_i)}{{\&PartialD;u}_i{\&PartialD;v}_i}{|}_{u_i=v_i=0},X_i=\frac{{\&PartialD;}^3{S}_i(u_i,v_i)}{{\&PartialD;u}_i^2{\&PartialD;v}_i}{|}_{u_i=v_i=0},$ $Y_i=\frac{{\&PartialD;}^3{S}_i(u_i,v_i)}{{\&PartialD;u}_i{\&PartialD;v}_i^2}{|}_{u_i=v_i=0},Z_i=\frac{{\&PartialD;}^4{S}_i(u_i,v_i)}{{\&PartialD;u}_i^2{\&PartialD;v}_i^2}{|}_{u_i=v_i=0};$ 对它们进行平均，将平均后的混合偏导作为曲面在P点处的新的混合偏导，并调整曲面相应的控制顶点，则曲面在P点处有公共的2至4阶混合偏导，在P点处C²相容。

如图2所示，P为奇异顶点，与P相连的曲面片为{S_i,i＝1,...,n}，边界曲线为{C_i,i＝1,...,n}。要使G¹连续的曲面片{S_i}在顶点P处扭矢相容，要求边界曲线{C_i}在顶点P处达到G²连续。我们在顶点P点处构造一张C²连续的基曲面S，其中S在P点的一阶、二阶导数分别为S_u、S_v、S_uu、S_uv、S_vv。边界曲线{C_i,i＝1,...,n}在P点处G²连续，则P点处新的切矢和曲率矢满足下述关系式：

$T_i^{*}={\mathrm{\&alpha;}}_i{S}_u+{\mathrm{\&beta;}}_i{S}_v,---\left(1\right)$

$H_i^{*}={\mathrm{\&alpha;}}_i^2{S}_{\mathrm{uu}}+{2\mathrm{\&alpha;}}_i{\mathrm{\&beta;}}_i{S}_{\mathrm{uv}}+{\mathrm{\&beta;}}_i^2{S}_{\mathrm{vv}}+{\mathrm{\&lambda;}}_i{S}_u+{\mathrm{\&gamma;}}_i{S}_v.---\left(2\right)$

其中α_i,β_i,λ_i,γ_i为常数。如果曲面片{S_i}在P点处的扭矢取为

$V_i^{*}={\mathrm{\&alpha;}}_{i+1}{\mathrm{\&alpha;}}_i{S}_{\mathrm{uu}}+({\mathrm{\&alpha;}}_{i+1}{\mathrm{\&beta;}}_i+{\mathrm{\&alpha;}}_i{\mathrm{\&beta;}}_{i+1})S_{\mathrm{uv}}+{\mathrm{\&beta;}}_{i+1}{\mathrm{\&beta;}}_i{S}_{\mathrm{vv}}+{\mathrm{\&xi;}}_i{S}_u+{\mathrm{\&eta;}}_i{S}_v,---\left(3\right)$

则{S_i}在顶点P处G¹相容。

为使(1)式成立，首先调整切矢T_i，使它们共面。每一块曲面片在P点处的法向N_i＝T_i×T_i+1，由{N_i,i＝1,...,n}估算曲面片在P点处的公共法向

$N=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{N}_i/\left|\right|\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{N}_i\left|\right|.$

由P和N可以确定公共切平面∏，将T_i投影到∏，得到新的切矢

$T_i^{*}=T_i-<T_i,N>*N(i=1,...,n);$

再根据这些新的切矢调整曲面片边界上对应的控制顶点。

不失一般性，取由下式可求得式(1)中每个所对应的系数(α_i,β_i)

${\mathrm{\&alpha;}}_i=<T_i^{*},S_u>/{\left|\right|S_u\left|\right|}^2,{\mathrm{\&beta;}}_i=<T_i^{*},S_v>/{\left|\right|S_v\left|\right|}^2.$

由(2)式，在法向N的投影为

$<H_i^{*},N>={\mathrm{\&alpha;}}_i^2<S_{\mathrm{uu}},N>+{2\mathrm{\&alpha;}}_i{\mathrm{\&beta;}}_i<S_{\mathrm{uv}},N>+{\mathrm{\&beta;}}_i^2<S_{\mathrm{vv}},N>,$

记s₀＝＜S_uu,N＞，s₁＝＜S_uv,N＞，s₂＝＜S_vv,N＞，则

$<H_i^{*},N>={\mathrm{\&alpha;}}_i^2{s}_0+{2\mathrm{\&alpha;}}_i{\mathrm{\&beta;}}_i{s}_1+{\mathrm{\&beta;}}_i^2{s}_2.---\left(4\right)$

用最小二乘法拟合求出(s₀,s₁,s₂)，目标函数为

$\mathrm{Min}\left(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(<H_i^{*},N>-<H_i,N>)}^2\right).$

求得后，由(4)式可以算出调整后的由下式求出

$H_i^{*}=H_i+(<H_i^{*},N>-<H_i,N>)N.$

同样的，有

$<V_i^{*},N>={\mathrm{\&alpha;}}_{i-1}{\mathrm{\&alpha;}}_i{s}_0+({\mathrm{\&alpha;}}_{i-1}{\mathrm{\&beta;}}_i+{\mathrm{\&alpha;}}_i{\mathrm{\&beta;}}_{i-1})s_1+{\mathrm{\&beta;}}_{i-1}{\mathrm{\&beta;}}_i{s}_2.$ 新的扭矢为

$V_i^{*}=V_i+(<V_i^{*},N>-<V_i,N>)N.$

最后根据这些新的曲率矢和扭矢调整曲面片的对应控制顶点，使{S_i,i＝1,...,n}在顶点P处G²连续并且G¹相容。

4)正则边界上的C²连续调整

如图3所示，P₀、P₁是两个正则顶点，曲面S_i-1与S_i的公共边界C_i(v)是一条正则边，S_i-1、S_i沿着公共边界方向具有相同的节点向量，如不一样，也可通过节点插入使它们相同。计算S_i-1、S_i的一阶、二阶跨界导矢曲线 对它们进行平均，得到新的跨界导矢曲线 ${\mathrm{DD}}_i^{+*}\left(v\right):$

$D_i^{-*}\left(v\right)=D_i^{+*}\left(v\right)=\frac{1}{2}(D_i^{-}\left(v\right)+D_i^{+}\left(v\right)),$

$D_i^{-*}\left(v\right)={\mathrm{DD}}_i^{+*}\left(v\right)=\frac{1}{2}({\mathrm{DD}}_i^{-}\left(v\right)+{\mathrm{DD}}_i^{+}\left(v\right)).$

再根据新的跨界导矢曲线调整曲面S_i-1、S_i靠近边界C_i(v)的第二、第三排控制顶点，使S_i-1、S_i在边界C_i(v)上C²连续。

5)非正则边界上的G¹连续调整

如图4所示，S_i-1、S_i是与非正则边界C_i(v)相连的两张双三次B样条曲面，沿着公共边界方向具有相同的节点向量，记节点向量为V＝{0,0,0,0,v₄,v₅,...,v_n,1,1,1,1}，其中v_k＜v_k+1(k＝3,...,n)，n为边界曲线C_i(v)的控制顶点个数。S_i-1、S_i在顶点P₀、P₁处G²连续且G¹相容。在这里，通过一种迭代逼近方法来调整曲面S_i-1、S_i沿着公共边界C_i(v)的一阶跨界导矢，使S_i-1、S_i在C_i(v)上达到近似的G¹连续。在工程上，当两相连曲面在公共边界上的法矢的夹角的最大值小于某个允许的误差值ε(比如ε＝0.5°)时，两曲面可看成是G¹连续。

该算法首先计算两曲面在公共边界上内部节点处的一阶跨界导矢根据切平面矫正调整这些跨界导矢，然后求出新的跨界导矢曲线使插值 $\{D_i^{-}\left(v_k\right),k=4,...,n\},$ 且满足如下端点条件：

$D_i^{-}\left(0\right)=T_{i-\mathrm{1,0}},D_i^{-}\left(1\right)=T_{i-\mathrm{1,1}},D_i^{{-}^{\&prime;}}\left(0\right)=V_{i-\mathrm{1,0}},D_i^{{-}^{\&prime;}}\left(1\right)=V_{i-\mathrm{1,1}};$

$D_i^{+}\left(0\right)=T_{i+\mathrm{1,0}},D_i^{+}\left(1\right)=T_{i+\mathrm{1,1}},D_i^{{+}^{\&prime;}}\left(0\right)=V_{i,0},D_i^{{+}^{\&prime;}}\left(1\right)=V_{i,1}.$

接着由新的跨界导矢曲线调整两个曲面靠近公共边界的第二排控制顶点。检查两曲面在公共边界上法矢的夹角，在夹角大于给定的误差允许值的地方插入节点，重新计算跨界导矢曲线。整个过程迭代进行，具体步骤如下：

Step0:设定初始值：两个曲面在公共边界上法矢夹角θ的最大误差为θ_max(比如0.5°)，在公共边界上允许插入的最大节点个数为N_max，当前插入的节点个数N_inserted＝0，边界曲线初始的节点向量为V＝{0,0,0,0,v₄,v₅,...,v_n,1,1,1,1}。

Step1:在节点{v_k,k＝4,...,n}处，计算曲线C_i(v)的切矢{C′_i(v_k),k＝4,...,n}，S_i-1、S_i沿着边界C_i(v)的跨界导矢S_i-1、S_i在点C_i(v_k)处的法矢为 $N_k^{-}=D_i^{-}\left(v_k\right)\&times;C_i^{\&prime;}\left(v_k\right),N_k^{+}=C_i^{\&prime;}\left(v_k\right)\&times;D_i^{+}\left(v_k\right)$ (k＝4,...,n)。

Step2:矫正采样点C_i(v_k)处的跨界导矢，使共面，满足G¹连续条件。首先确定两个曲面在点C_i(v_k)处的公共法向

$N_k=(N_k^{-}+N_k^{+})/\left|\right|N_k^{-}+N_k^{+}\left|\right|;$

矫正后的跨界导矢为

$D_i^{-}{\left(v_k\right)}^{*}=D_i^{-}\left(v_k\right)-<D_i^{-}\left(v_k\right),N_k>*N_k,$

$D_i^{+}{\left(v_k\right)}^{*}=D_i^{+}\left(v_k\right)-<D_i^{+}\left(v_k\right),N_k>*N_k.$

Step3:以当前边界曲线C_i(v)的节点向量V作为跨界导矢曲线的节点向量，拟合新的跨界导矢曲线使插值 $D^{-}\{T_{i-\mathrm{1,0}},D_i^{-}{\left(v_4\right)}^{*},...,D_i^{-}{\left(v_n\right)}^{*},T_{i-\mathrm{1,1}}\}$ 且 $D_i^{-\&prime;}\left(0\right)=V_{i-\mathrm{1,0}},D_i^{-\&prime;}\left(1\right)=V_{i-\mathrm{1,1}},$ 插值 $D^{+}\{T_{i+\mathrm{1,0}},D_i^{+}{\left(v_4\right)}^{*},...,D_i^{+}{\left(v_n\right)}^{*},T_{i+\mathrm{1,1}}\}$ 且 $D_i^{+\&prime;}\left(0\right)=V_{i,0},D_i^{+\&prime;}\left(1\right)=V_{i,1}.$

Step4:根据矫正后的跨界导矢曲线调整S_i-1、S_i靠近边界C_i(v)的第二排控制顶点。

Step5:经过上面的调整，S_i-1、S_i在节点{v_k,k＝3,...,n+1}处已经达到G¹连续。在每个节点区间[v_k,v_k+1](k＝3,...,n)内再取m个采样点

$v_{k,j}=v_k+\frac{j}{m+1}(v_{k+1}-t_k),(j=1,...,m);$

计算S_i-1、S_i在节点{v_k,j,k＝3,...,n,j＝1,...,m}处的跨界导矢法矢求出两个法向的夹角{θ_k,j}。

Step6:找出每个区间内夹角最大处的节点值{v_k,j,k＝3,...,n,j＝p_k}，如果该处的夹角大于θ_max，则将加入到待插入的节点集合,记待插入节点集合为 $V^{*}=\{v_l^{*},l=0,...,s\}.$

Step7:如果，将节点分别插入到曲线C_i(v)及曲面S_i-1、S_i中去，得到新的节点向量V＝V∪V^*，计算当前已插入的节点数N_inserted＝N_inserted+s+1；矫正处的跨界导矢，得到新的跨界导矢 新的待拟合的跨界导矢集合为D^-＝D^-∪D^-*、D⁺＝D⁺∪0^+*。

Step8:重复Step3－Step7，直到或者N_inserted＞N_max，算法结束。

至此，曲面S_i-1、S_i沿着公共边界C_i(v)在节点处达到精确的G¹连续，在其它地方达到近似的G¹连续，且在拟合跨界曲线时将矫正后的扭矢作为插值条件，解决了顶点处的相容性问题。

本实施例的经过四边形划分的网格模型上重建连续拼接B样条曲面片的方法，首先构造在正则顶点处C²连续的曲线网格，然后由边界曲线和内部数据点拟合得到一个G⁰连续的双三次B样条曲面片网格，再通过跨界导矢矫正法使曲面片间连续拼接。该算法较好地解决了曲面片在顶点处的相容性，并使相连的曲面片在正则边界上保持C²连续，在非正则边界上保持近似的G¹连续。算法实现简单，速度快，对四边形网格的拓扑结构没有限制，重建的曲面片在边界上保持光滑拼接的同时，内部逼近原始模型上的数据点，适用于逆向工程中对经过任意四边形划分的模型进行光滑B样条曲面重建。图5、6给出了利用该算法对四边形划分的三角网格模型进行B样条曲面重建的两个实例。图5a是个封闭的三角网格模型，有62856个面片，31432个顶点，划分为222个四边形区域；图5b是重建的曲面模型，曲面重建过程共耗时4.5s，测试环境为Inteli53.0GHz,4GB内存。图6a是个非封闭的三角网格模型，有30355个面片，15304个顶点，划分为367个四边形区域，图6b是重建的曲面模型，曲面重建过程共耗时5.6s。从图中可以看出，重建的曲面片在边界上具有较好的连续性，同时在曲面内部保持了模型的几何特征。

Claims (4)

1.任意四边形网格上光滑曲面的重建方法，其特征在于：所述重建方法包括下述步骤：

1)生成四边形区域边界曲线

用三次B样条曲线拟合四边形区域边界上的采样点得到边界曲线，而且对所有的边界曲线采用相同的节点向量，使相对的两条边的节点向量相容；

在生成边界曲线时，首先找出四边形网格上的长边，长边由依次相连的多条边组成，长边的首末端点为奇异顶点或网格边界上的正则顶点；得到长边后，取出组成长边的每条边上的采样点，将它们合在一起，拟合出一条长曲线，该曲线插值长边上的顶点，逼近每条边上的采样点；再将长曲线在顶点处打断，得到每条边上对应的边界曲线，使边界曲线在正则顶点处C²连续，奇异顶点处G⁰连续；

2)对每个四边形区域，由边界曲线和内部数据点拟合生成B样条曲面；

3)在正则顶点处进行C²相容性调整，再进行奇异顶点处G¹相容性调整，过程如下：

3.1)正则顶点处的C²相容性调整按下述方法进行：

计算与正则顶点P相连的4个曲面{S_i,i＝1,...,4}在P处的混合偏导 $V_i=\frac{{\&PartialD;}^2{S}_i(u_i,v_i)}{{\&PartialD;u}_i{\&PartialD;v}_i}{|}_{u_i=v_i=0},X_i=\frac{{\&PartialD;}^3{S}_i(u_i,v_i)}{{\&PartialD;u}_i^2{\&PartialD;v}_i}{|}_{u_i=v_i=0},Y_i=\frac{{\&PartialD;}^3{S}_i(u_i,v_i)}{{\&PartialD;u}_i{\&PartialD;v}_i^2}{|}_{u_i=v_i=0},Z_i=\frac{{\&PartialD;}^4{S}_i(u_i,v_i)}{{\&PartialD;u}_i^2{\&PartialD;v}_i^2}{|}_{u_i=v_i=0};$ 对它们进行平均，将平均后的混合偏导作为曲面在P点处的新的混合偏导，并调整曲面相应的控制顶点，则曲面在P点处有公共的2至4阶混合偏导，在P点处C²相容；

3.2)奇异顶点处G¹相容性调整按下述方法进行：

记与奇异顶点P相连的曲面片为{S_i(u_i,v_i),i＝1,...,n}，边界曲线为{C_i,i＝1,...,n}；要使G¹连续的曲面片{S_i(u_i,v_i),i＝1,...,n}在顶点P处扭矢相容，要求边界曲线{C_i}在顶点P处达到G²连续；在顶点P点处构造一张C²连续的基曲面S，其中S在P点的一阶、二阶导数分别为S_u、S_v、S_uu、S_uv、S_vv；边界曲线{C_i,i＝1,...,n}在P点处G²连续，则P点处新的切矢和曲率矢满足下述关系式：

$T_i^{*}={\mathrm{\&alpha;}}_i{S}_u+{\mathrm{\&beta;}}_i{S}_v,$

$H_i^{*}={\mathrm{\&alpha;}}_i^2{S}_{\mathrm{uu}}+{2\mathrm{\&alpha;}}_i{\mathrm{\&beta;}}_i{S}_{\mathrm{uv}}+{\mathrm{\&beta;}}_i^2{S}_{\mathrm{vv}}+{\mathrm{\&lambda;}}_i{S}_u+{\mathrm{\&gamma;}}_i{S}_v,$

其中α_i,β_i,λ_i,γ_i为常数，将曲面片{S_i(u_i,v_i),i＝1,...,n}在P点处的扭矢取为

$V_i^{*}={\mathrm{\&alpha;}}_{i+1}{\mathrm{\&alpha;}}_i{S}_{\mathrm{uu}}+({\mathrm{\&alpha;}}_{i+1}{\mathrm{\&beta;}}_i+{\mathrm{\&alpha;}}_i{\mathrm{\&beta;}}_{i+1})S_{\mathrm{uv}}+{\mathrm{\&beta;}}_{i+1}{\mathrm{\&beta;}}_i{S}_{\mathrm{vv}}+{\mathrm{\&xi;}}_i{S}_u+{\mathrm{\&eta;}}_i{S}_v;$

则{S_i(u_i,v_i),i＝1,...,n}在P点处G¹相容；

首先在P点处求出一个公共切平面Π，将T_i投影到Π，得到新的切矢然后用最小二乘法求出曲率矢和扭矢在法向N上的投影调整后的曲率矢和扭矢为

$H_i^{*}=H_i+(<H_i^{*},N>-<H_i,N>)N,V_i^{*}=V_i+(<V_i^{*},N>-<V_i,N>)N;$

求得新的曲率矢和扭矢后，调整曲面片的相应控制顶点，使{S_i,i＝1,...,n}在顶点P处G²连续并且G¹相容；

4)正则边界上进行C²连续调整

计算与正则边C_i(v)相连的两曲面S_i-1、S_i的一阶、二阶跨界导矢曲线 对它们进行平均，得到新的跨界导矢曲线 $D_i^{+*}\left(v\right),{\mathrm{DD}}_i^{+*}\left(v\right):$

$D_i^{-*}\left(v\right)=D_i^{+*}\left(v\right)=\frac{1}{2}(D_i^{-}\left(v\right)+D_i^{+}\left(v\right)),{\mathrm{DD}}_i^{-*}\left(v\right)={\mathrm{DD}}_i^{+*}\left(v\right)=\frac{1}{2}({\mathrm{DD}}_i^{-}\left(v\right)+{\mathrm{DD}}_i^{+}\left(v\right));$

再根据新的跨界导矢曲线调整曲面S_i-1、S_i靠近边界C_i(v)的第二、第三排控制顶点，使S_i-1、S_i在边界C_i(v)上C²连续；

5)非正则边界上进行近似G¹连续调整

首先计算两曲面在公共边界上内部节点处的一阶跨界导矢 根据切平面矫正调整这些跨界导矢，然后求出新的跨界导矢曲线使插值 $\{D_i^{-}\left(v_k\right),k=4,...,n\},\{D_i^{+}\left(v_k\right),k=4,...,n\},$ 且满足如下端点条件：

$D_i^{-}\left(0\right)=T_{i-\mathrm{1,0}},D_i^{-}\left(1\right)=T_{i-\mathrm{1,1}},D_i^{{-}^{\&prime;}}\left(0\right)=V_{i-\mathrm{1,0}},D_i^{{-}^{\&prime;}}\left(1\right)=V_{i-\mathrm{1,1}};$

$D_i^{+}\left(0\right)=T_{i+\mathrm{1,0}},D_i^{+}\left(1\right)=T_{i+\mathrm{1,1}},D_i^{{+}^{\&prime;}}\left(0\right)=V_{i,0},D_i^{{+}^{\&prime;}}\left(1\right)=V_{i,1}$

接着由新的跨界导矢曲线调整两个曲面靠近公共边界的第二排控制顶点，检查两曲面在公共边界上法矢的夹角，在夹角大于给定的误差允许值的地方插入节点，重新计算跨界导矢曲线；整个过程迭代进行，直到各个节点区间的法矢夹角都小于给定的误差允许值或插入的节点数量达到一个上限值。

2.如权利要求1所述的任意四边形网格上光滑曲面的重建方法，其特征在于：所述步骤2)中，给定四边形区域的四条边界曲线及内部数据点集，通过最小二乘拟合重构一张双三次B样条曲面，重构曲面插值四条边界曲线，逼近内部点。

3.如权利要求2所述的任意四边形网格上光滑曲面的重建方法，其特征在于：所述步骤2)中，首先由4条边界曲线创建一张双线性Coons曲面，然后将散乱点P_k投影到该基曲面计算得到点P_k的参数值(u_k,v_k)；曲面S(u,v)边界上的控制顶点由4条边界曲线的控制顶点确定，内部的控制顶点通过最小二乘拟合如下的目标函数F得到，

F＝F_lsq+λF_fair，其中F_lsq为最小二乘项，F_fair为光顺项，λ为光顺项系数。

4.如权利要求1～3之一所述的任意四边形网格上光滑曲面的重建方法，其特征在于：所述步骤5)的具体步骤如下：

Step0:设定初始值：两个曲面在公共边界上法矢夹角θ的最大误差为θ_max，在公共边界上允许插入的最大节点个数为N_max，当前插入的节点个数N_inserted＝0，边界曲线初始的节点向量为V＝{0,0,0,0,v₄,v₅,...,v_n,1,1,1,1}；

Step1:在节点{v_k,k＝4,...,n}处，计算曲线C_i(v)的切矢{C′_i(v_k),k＝4,...,n}，S_i-1、S_i沿着边界C_i(v)的跨界导矢S_i-1、S_i在点C_i(v_k)处的法矢为 $N_k^{-}=D_i^{-}\left(v_k\right)\&times;C_i^{\&prime;}\left(v_k\right),N_k^{+}=C_i^{\&prime;}\left(v_k\right)\&times;D_i^{+}\left(v_k\right)$ (k＝4,...,n)；

Step2:矫正采样点C_i(v_k)处的跨界导矢，使C′_i(v_k)、共面，满足G¹连续条件，首先确定两个曲面在点C_i(v_k)处的公共法向

$N_k=(N_k^{-}+N_k^{+})/\left|\right|N_k^{-}+N_k^{+}\left|\right|;$

矫正后的跨界导矢为

$D_i^{-}{\left(v_k\right)}^{*}=D_i^{-}\left(v_k\right)-<D_i^{-}\left(v_k\right),N_k>*N_k,$

$D_i^{+}{\left(v_k\right)}^{*}=D_i^{+}\left(v_k\right)-<D_i^{+}\left(v_k\right),N_k>*N_k$

Step3:以当前边界曲线C_i(v)的节点向量V作为跨界导矢曲线的节点向量，拟合新的跨界导矢曲线使插值 $D^{-}\{T_{i-\mathrm{1,0}},D_i^{-}{\left(v_4\right)}^{*},...,D_i^{-}{\left(v_n\right)}^{*},T_{i-\mathrm{1,1}}\}$ 且 $D_i^{-\&prime;}\left(0\right)=V_{i-\mathrm{1,0}},D_i^{-\&prime;}\left(1\right)=V_{i-\mathrm{1,1}},$ 插值 $D^{+}\{T_{i+\mathrm{1,0}},D_i^{+}{\left(v_4\right)}^{*},...,D_i^{+}{\left(v_n\right)}^{*},T_{i+\mathrm{1,1}}\}$ 且 $D_i^{+\&prime;}\left(0\right)=V_{i,0},D_i^{+\&prime;}\left(1\right)=V_{i,1};$

Step4:根据矫正后的跨界导矢曲线调整S_i-1、S_i靠近边界C_i(v)的第二排控制顶点；

Step5:经过上述调整，S_i-1、S_i在节点{v_k,k＝3,...,n+1}处已经达到G¹连续，在每个节点区间[v_k,v_k+1](k＝3,...,n)内再取m个采样点

$v_{k,j}=v_k+\frac{j}{m+1}(v_{k+1}-t_k),(j=1,...,m);$

计算S_i-1、S_i在节点{v_k,j,k＝3,...,n,j＝1,...,m}处的跨界导矢法矢求出两个法向的夹角{θ_k,j}；

Step6:找出每个区间内夹角最大处的节点值{v_k,j,k＝3,...,n,j＝p_k}，如果该处的夹角大于θ_max，则将加入到待插入的节点集合,记待插入节点集合为 $V^{*}=\{v_l^{*},l=0,...,s\};$

Step7:如果将节点分别插入到曲线C_i(v)及曲面S_i-1、S_i中去，得到新的节点向量V＝V∪V^*，计算当前已插入的节点数N_inserted＝N_inserted+s+1；矫正处的跨界导矢，得到新的跨界导矢 $D^{-*}=\{D_i^{-}\left(v_l^{*}\right),l=0,...,s\},D^{+*}=\{D_i^{+}\left(v_l^{*}\right),l=0,...,s\},$ 新的待拟合的跨界导矢集合为D^-＝D^-∪D^-*、D⁺＝D⁺∪D^+*；

Step8:重复Step3－Step7，直到或者N_inserted＞N_max，算法结束。