CN104392021B - 一种层次自由曲面的保角参数特性优化方法 - Google Patents

Abstract

本发明公开了一种层次自由曲面的保角参数特性优化方法，初始阶段输入一个自由曲面，使用组合辛普森规则构造其保角能量离散逼近形式；使用里奇流方法构建自由变换的参数形式，使最优自由变换问题转换为最小二乘问题；计算其保角误差，反复迭代可得到满足误差范围内的具有最优保角性的层次自由曲面。该方法在不改变用户给定的自由曲面形状和次数的前提下，通过对自由曲面进行自由变换引入层次自由曲面，实现了对自由曲面的保角特性的优化。

Description

一种层次自由曲面的保角参数特性优化方法

技术领域

本发明涉及一种层次自由曲面的保角参数特性优化方法。

背景技术

国际标准化组织(ISO)于1991年正式确定把NURBS方法作为定义产品形状的唯一数学方法，越来越多的CAD/CAM系统采用NURBS作为其模型，NURBS在工业产品的外形设计中获得广泛的应用。纹理映射、曲面镶嵌、曲面匹配和注册的结果都和曲面的保角性息息相关。

在现代图形学和建筑学中，自由曲面的作用日益重要。曲面渲染(例如纹理映射)、曲面离散、曲面采样等自由曲面应用十分依赖曲面的保角特性。已有的保角映射的成果大多针对离散三角网格曲面，到目前为止还未有三维NURBS自由曲面保角性的相关成熟结果。

设计人员通常先对自由曲面保角特性进行优化，以满足曲面采样、曲面相交、曲率计算等算法的要求。因此基于NURBS曲面保角参数特性优化方法具有很强的实用价值。

发明内容

本发明的目的是为克服上述现有技术的不足，提供一种层次自由曲面的保角参数特性优化方法。该方法通过引入层次自由曲面和自由变换，在不改变原有自由曲面形状的前提下，实现自由曲面保角参数特性的优化

为实现上述目的，本发明采用下述技术方案：

一种层次自由曲面的保角参数特性优化方法，包括如下步骤：

1)输入NURBS空间自由曲面X(u，v)；

2)将自由曲面X(u，v)的保角能量公式化，得到保角能量函数；并对保角能量函数进行离散化处理；

3)采用最小二乘方法极小化离散后的保角能量函数，并计算出相应的优化后的自由曲面；

4)对新的层次自由曲面采样，计算保角性误差，将所得保角性误差与设定的阈值进行比较：如果所述误差小于设定阈值，则所得层次自由曲面具有保角性；如果所述误差大于设定阈值，则所得层次自由曲面不具有保角性，分别增加u，v上采样点数目，提高(u，v)参数域采样精度，返回步骤1)重新计算。

所述步骤1)中NURBS空间自由曲面X(u，v)具体为：

$X(u,v)=\frac{{\mathrm{\Σ}}_{i=0}^{n_u}{\mathrm{\Σ}}_{j=0}^{n_v}{N}_i^{p_1}\left(u\right)N_j^{q_1}\left(v\right){\mathrm{\ω}}_{i,j}{P}_{i,j}}{{\mathrm{\Σ}}_{i=0}^{n_u}{\mathrm{\Σ}}_{j=0}^{n_v}{N}_i^{p_1}\left(u\right)N_j^{q_1}\left(v\right){\mathrm{\ω}}_{i,j}},u,v\∈\left[\mathrm{0,1}\right]$

其中，u和v是曲面X(u，v)的参数，n_u和n_v是曲面X(u，v)所含的控制顶点在u方向和v方向上的个数，P_i，j是曲面(u，v)中序号为(i，j)的控制顶点的坐标值，由用户输入，i＝0，1，...，n_u-1；j＝0，1，...，n_v-1，采用三维直角坐标表示，ω_i，j是曲面X(u，v)中序号为(i，j)的控制顶点的权重，由用户输入，i＝0，1，...，n_u-1；j＝0，1，...，n_v-1，是定义在节点向量序列U之上的序号为i的第p₁阶B-spline基函数，i＝0，1，...，n_u-1；是定义在节点向量序列V之上的序号为j的第q₁阶B-spline基函数，j＝0，1，...，n_v-1，

其中，p，q是曲面X(u，v)在u方向和v方向上的次数，由用户输入，u₀，u₁，...，u_n+p等为节点向量序列U中的节点，其中u₀＝u₁＝…＝u_p＝0，u_n＝u_n+1＝…＝u_n+p＝1，u_p+1，u_n+2，...，u_n-1由用户输入；v₀，v₁，...，v_n+q等为节点向量序列V中的节点，其中v₀＝v₁＝…＝v_q＝0，v_n＝v_n+1＝…＝v_n+q＝1，v_q+1，v_q+2，...，v_n-1由用户输入，上述的B-spline基函数采用Matlab系统样条工具库中B-spline基函数的定义方式。

所述步骤2)中将自由曲面X(u，v)的保角能量公式化，得到保角能量函数的具体为：

$J_c\left(Q_{i,j}\right)={{\∫}_0^1{\∫}_0^1|{\mathrm{\φ}}_t-\left(\begin{array}{cc}-\frac{F}{\sqrt{\mathrm{EG}-F^2}}& -\frac{G}{\sqrt{\mathrm{EG}-F^2}}\\ \frac{E}{\sqrt{\mathrm{EG}-F^2}}& \frac{F}{\sqrt{\mathrm{EG}-F^2}}\end{array}\right){\mathrm{\φ}}_s|}^2\mathrm{dsdt}$

其中，确定曲面X(u，v)的基本形式为：

ds²＝X_u·X_v(du)²+2X_u·X_vdudv+X_v·X_v(dv)²

为曲面X(u，v)的两个偏导数；

E＝X_u·X_v，F＝X_u·X_v，G＝X_v·X_v

其中，E和G代表了两个偏导的长度，F用来度量两个偏导的正交性；

由层次自由曲面变换可得

$X(u,v)=\mathrm{\φ}(s,t)={\mathrm{\Σ}}_{i=0}^{n_s}{\mathrm{\Σ}}_{j=0}^{n_t}{N}_i^{p_2}\left(s\right)N_j^{q_2}\left(t\right)Q_{i,j},$

上述保角能量函数中φ_s和φ_t分别是变换后自由曲面对s和t的偏导数。

所述步骤2)中对保角能量函数利用复合辛普森规则逼近，为了应用该规则，对保角能量函数进行离散化处理的具体方法如下：

在曲面X(u，v)的参数域上进行离散采样，分别将参数域u，v细分为2l和2k个子区间，其均分步长分别为h＝1/(2l)，p＝1/(2k)；细分后参数域如下：

$R_{i,j}=\left(\begin{array}{c}{u}_p+\frac{i\·(u_{n_u+1}-u_p)}{2l}\\ v_q+\frac{i\·(v_{n_v+1}-v_q)}{2k}\end{array}\right),i=0,...2l,j=0,...2k.$

其中，u_p、v_q、均为参数域节点向量中的节点；

使用里奇流方法计算出对应离散化后曲面X(u，v)对应的保角映射M，令

$T={|{\mathrm{\φ}}_t-\left(\begin{array}{cc}-\frac{F}{\sqrt{\mathrm{EG}-F^2}}& -\frac{G}{\sqrt{\mathrm{EG}-F^2}}\\ \frac{E}{\sqrt{\mathrm{EG}-F^2}}& \frac{F}{\sqrt{\mathrm{EG}-F^2}}\end{array}\right){\mathrm{\φ}}_s|}^2$

其中，E＝X_u·X_v，F＝X_u·X_v，G＝X_v·X_v，E和G代表了两个偏导的长度，F用来度量两个偏导的正交性；φ_s和φ_t分别是层次自由曲面变换(u，v)＝φ(s，t)后，对s和t的偏导数；

利用辛普森规则的逼近后的保角能量结果应有如下形式：

$J=\tilde{J}+E_{\mathrm{error}}$

其中，参数域中存在某一个和使得

$\left\{\begin{array}{c}\tilde{J}=\frac{4\mathrm{hp}}{9}\{J_1+{2J}_2+{4J}_3+J_4\}\\ E_{\mathrm{error}}=-\frac{h^4\frac{{\∂}^{4}T(\stackrel{\‾}{s},\stackrel{\‾}{t})}{\∂a^4}+p^4\frac{{\∂}^{4}T(\hat{s},\hat{t})}{\∂t^4}}{180}\\ J_1=T(s_{\mathrm{0,0}},t_{\mathrm{0,0}})+2\underset{i=1}{\overset{l-1}{\mathrm{\Σ}}}T(s_{2i,0},t_{2i,0})+4\underset{i=1}{\overset{l}{\mathrm{\Σ}}}T(s_{2i-\mathrm{1,0}},t_{2i-\mathrm{1,0}})+T(s_{2l,0},t_{2l,0})\\ J_2=\underset{j=1}{\overset{k-1}{\mathrm{\Σ}}}T(s_{0,2j},t_{0,2j})+2\underset{j=1}{\overset{k-1}{\mathrm{\Σ}}}\underset{i=1}{\overset{l-1}{\mathrm{\Σ}}}T(s_{2i,2j},t_{2i,2j})+\\ 4\underset{j=1}{\overset{k-1}{\mathrm{\Σ}}}\underset{i=1}{\overset{l}{\mathrm{\Σ}}}T(s_{2i-1,2j},t_{2i-1,2j})+\underset{j=1}{\overset{k-1}{\mathrm{\Σ}}}T(s_{2l,2j},t_{2l,2j})\\ J_3=\underset{j=1}{\overset{k}{\mathrm{\Σ}}}T(s_{0,2j-1},t_{0,2j-1})+2\underset{j=1}{\overset{k}{\mathrm{\Σ}}}\underset{i=1}{\overset{l-1}{\mathrm{\Σ}}}T(s_{2i,2j-1},t_{2i,2j-1})+\\ 4\underset{j=1}{\overset{k}{\mathrm{\Σ}}}\underset{i=1}{\overset{l}{\mathrm{\Σ}}}T(s_{2i-1,2j-1},t_{2i-1,2j-1})+\underset{j=1}{\overset{k}{\mathrm{\Σ}}}T(s_{2l,2j-1},t_{2l,2j-1})\\ J_4=T(s_{0,2k},t_{0,2k})+2\underset{i=1}{\overset{l-1}{\mathrm{\Σ}}}T(s_{2i,2k},t_{2i,2k})+4\underset{i=1}{\overset{l}{\mathrm{\Σ}}}T(s_{2i-1,2k},t_{2i-1,2k})+T(s_{2l,2k},t_{2l,2k})\end{array}\right.$

其中，J₁，J₂，J₃，J₄是其四个离散项，s_i，j，t_i，j是参数中的采样点u，v坐标，h和p分别是u，v方向的采样间隔，T是保角项；通过增加采样密度l和k来减少E_error项，当l和k均为10时，误差项就变成了如下形式

$E_{\mathrm{error}}=-\frac{\frac{{\∂}^{4}T(\stackrel{\‾}{u},\stackrel{\‾}{v})}{\∂u^4}+\frac{{\∂}^{4}T(\hat{u},\hat{v})}{\∂v^4}}{2.88\×{10}^{11}}$

其中u，v是自由曲面的两个参数，T是保角项；优化其如下离散保角能量：

${\tilde{J}}_p\left(Q_{i,j}\right)=\underset{i=0}{\overset{2l}{\mathrm{\Σ}}}\underset{j=0}{\overset{2k}{\mathrm{\Σ}}}{|R_{i,j}-\mathrm{\φ}(s_{i,j},t_{i,j})|}^2$

R_i，j是u，v参数域中的采样点，s_i，j，t_i，j是st参数域中的两个参数；φ是二维自由曲面变换，Q_i，j是二维变换φ的控制顶点；优化离散保角能量使得生成的层次曲面的参数化逼近保角性。

所述步骤3)中采用最小二乘方法极小化离散后的保角能量函数的具体方法为：

令{R(i，j)}和{C(s_i，j，t_i，j)}，i＝0，…，2l，j＝0，…，2k分别表示采样点集合和其对应满足的保角性约束条件集合；

将{R(i，j)}和{C(s_i，j，t_i，j)}分为无约束部分和有约束部分：使用和来表示无约束的部分；使用和来表示有约束的部分；给有约束部分设定非负权重，有和增加权重可以增加每个数据点逼近其邻近点的逼近程度，减小权重可以减弱每个数据点逼近其邻近点的逼近程度；

令m_u＝r_u+s_u+1，m_c＝r_c+s_c+1，其中m_u，r_u，s_u表示无约束部分元素的数目，m_c，r_c，s_c表示有约束部分元素的数目，且满足m_c＜n，m_c+n＜m_u+1；

令

S_k，k＝0，…，m_u是第k个无约束数据点；

T_k，k＝0，…，m_c是第k个有约束数据点；

ω_k，k＝0，…，mu是第k个无约束数据点的约束项；

定义向量和矩阵如下：

S＝[S_k]，为一个有m_u+1个元素的向量；

T＝[T_k]，为一个有m_c+1个元素的向量；

W＝[ω_k]，为一个(m_u+1)×(m_u+1)的对角矩阵，ω_k在对角线上；

Q＝[Q_i，j]，为一个(n_s+1)×(n_t+1)的未知控制定点向量；

N＝[ND_k]，其中ND_k是一个k次基函数或者是保角项的基函数，N为一个(m_u+1)×(n_s+1)(n_t+1)的标量矩阵；

M＝[MD_k]，其中MD_k是一个k次基函数或者是保角项的基函数，M为一个(m_u+1)×(n_s+1)(n_t+1)的标量矩阵；

令A＝[λ_k]，A是拉格朗日乘子的向量，每个λ_k都是和Q_i，j维度相同的向量；这里，无约束部分有

NQ＝S

有约束部分有

MQ＝T

无约束部分误差是S–NP，我们希望通过极小化这个误差从而使得MP＝T；因此，利用拉格朗日乘子法，等价极小化

(S^T-PN^T)W(S–NP)+A^T(MP-T)

对于未知数A，P；对上面的式子进行微分，令微分为0

-2(S^TWN-P^TN^TWN)+A^TM＝0

MP-T＝0

表示为矩阵形式为

$\left[\begin{array}{cc}{N}^T\mathrm{WN}& M^T\\ M& 0\end{array}\right]\left[\begin{array}{c}P\\ A\end{array}\right]=\left[\begin{array}{c}{N}^T\mathrm{WS}\\ T\end{array}\right]$

通过求解上式我们可得到A和P；解唯一的条件是N^TWN和M(N^TWN)^-1M^T都是可逆的；可得

$\left\{\begin{array}{c}P={\left(N^T\mathrm{WN}\right)}^{-1}{N}^T\mathrm{WS}-{\left(N^T\mathrm{WN}\right)}^{-1}{M}^{T}A\\ A={\left(M{\left(N^T\mathrm{WN}\right)}^{-1}{M}^T\right)}^{-1}({M\left(N^T\mathrm{WN}\right)}^{-1}{N}^T\mathrm{WS}-T)\end{array}\right.$

最终P为

P＝(N^TWN)^-1N^TWS–(N^TWN)^-1M^T(M(N^TWN)^-1M^T)^-1(M(N^TWN)^-1N^TWS-T)

该P就是通过最小二乘法得到的满足要求的解。

所述步骤4)中新的层次自由曲面采样点的最大的保角性误差为：

$\underset{0\≤k\≤m}{\max }{|N\×\frac{\∂X}{\∂s}-\frac{\∂X}{\∂t}|}^2$

其中，N是曲面法向量， 分别是曲面X在s，t两个方向上的切矢。

本发明有益效果：

在实际应用中，本发明方法简单可行。在不改变用户给定的自由曲面形状和次数的前提下，通过对自由曲面进行线性重新参数化，只改变自由曲面的控制顶点、权因子以及节点向量，实现了对自由曲面的保角特性的优化，并且保证了最优保角曲面的唯一性。

附图说明

图1为输入的一个人脸NURBS曲面X(u，v)；

图2为对NURBS曲面X(u，v)进行保角参数优化的最优结果；

图3是整个方法的步骤流程图。

具体实施方式

下面通过具体实例和附图对本发明进行进一步的阐述，应该说明的是，下述说明仅是为了解释本发明，并不对其内容进行限定。

本CAD方法以Windows7中的VisualStudio2008为开发平台，其具体实施方式如下：

图1是输入一个NURBS空间曲面X(u，v)，对于曲面X(u，v)其形式为：

$X(u,v)=\frac{{\mathrm{\Σ}}_{i=0}^{n_u}{\mathrm{\Σ}}_{j=0}^{n_v}{N}_i^{p_1}\left(u\right)N_j^{q_1}\left(v\right){\mathrm{\ω}}_{i,j}{P}_{i,j}}{{\mathrm{\Σ}}_{i=0}^{n_u}{\mathrm{\Σ}}_{j=0}^{n_v}{N}_i^{p_1}\left(u\right)N_j^{q_1}\left(v\right){\mathrm{\ω}}_{i,j}},u,v\∈\left[\mathrm{0,1}\right]$

其中u和v是曲面X(u，v)的参数，n_u和n_v是曲面X(u，v)所含的控制顶点在u方向和v方向上的个数，P_i，j是曲面(u，v)中序号为(i，j)的控制顶点的坐标值，由用户输入，i＝0，1，...，n_u-1；j＝0，1，...，n_v-1，采用三维直角坐标表示，ω_i，j是曲面X(u，v)中序号为(i，j)的控制顶点的权重，由用户输入，i＝0，1，...，n_u-1；j＝0，1，...，n_v-1，是定义在节点向量序列U之上的序号为i的第p₁阶B-spline基函数，i＝0，1，...，n_u-1；是定义在节点向量序列V之上的序号为j的第q₁阶B-spline基函数，j＝0，1，...，n_v-1，

其中p，q是曲面X(u，v)在u方向和v方向上的次数，由用户输入，u₀，u₁，...，u_n+p等为节点向量序列U中的节点，其中u₀＝u₁＝…＝u_p＝0，u_n＝u_n+1＝…＝u_n+p＝1，u_p+1，u_n+2，...，u_n-1由用户输入；v₀，v₁，...，v_n+q等为节点向量序列V中的节点，其中v₀＝v₁＝…＝v_q＝0，v_n＝v_n+1＝…＝v_n+q＝1，v_q+1，v_q+2，...，v_n-1由用户输入，上述的B-spline基函数采用Matlab系统样条工具库中B-spline基函数的定义方式；

S2:保角能量离散化

S21：保角能量公式化

X(u，v)可表现为其基本形式，

ds²＝X_u·X_v(du)²+2X_u·X_vdudv+X_v·X_v(dv)²

曲面X的两个偏导数。第一个基本形式描述了曲面X的度量。设定：

E＝X_u·X_v，F＝X_u·X_v，G＝X_v·X_v

其中的E和G代表了两个偏导的长度，F用来度量两个偏导的正交性。

保角能量公式最终为

$J_c\left(Q_{i,j}\right)={{\∫}_0^1{\∫}_0^1|{\mathrm{\φ}}_t-\left(\begin{array}{cc}-\frac{F}{\sqrt{\mathrm{EG}-F^2}}& -\frac{G}{\sqrt{\mathrm{EG}-F^2}}\\ \frac{E}{\sqrt{\mathrm{EG}-F^2}}& \frac{F}{\sqrt{\mathrm{EG}-F^2}}\end{array}\right){\mathrm{\φ}}_s|}^2\mathrm{dsdt}$

S22：在计算机中，采用如下方法离散保角能量：

使用混合辛普森规则离散化得到的能量进行逼近。为了使用混合辛普森规则，我们需要在参数域上均匀采样，将u，v均匀分为2l和2k份。

$R_{i,j}=\left(\begin{array}{c}{u}_p+\frac{i\·(u_{n_u+1}-u_p)}{2l}\\ v_q+\frac{i\·(v_{n_v+1}-v_q)}{2k}\end{array}\right),i=0,...2l,j=0,...2k.$

这样的细分方法决定了步长分别为h＝1/(2l)，p＝1/(2k)。为了避免非线性问题，我们使用里奇流方法计算出对应离散化后的3DNURBS曲面的对应的保角映射M。假定

$T={|{\mathrm{\φ}}_t-\left(\begin{array}{cc}-\frac{F}{\sqrt{\mathrm{EG}-F^2}}& -\frac{G}{\sqrt{\mathrm{EG}-F^2}}\\ \frac{E}{\sqrt{\mathrm{EG}-F^2}}& \frac{F}{\sqrt{\mathrm{EG}-F^2}}\end{array}\right){\mathrm{\φ}}_s|}^2$

逼近结果形式如下

$J=\tilde{J}+E_{\mathrm{error}}$

其中对于若干参数域中得和有

$\left\{\begin{array}{c}\tilde{J}=\frac{4\mathrm{hp}}{9}\{J_1+{2J}_2+{4J}_3+J_4\}\\ E_{\mathrm{error}}=-\frac{h^4\frac{{\∂}^{4}T(\stackrel{\‾}{s},\stackrel{\‾}{t})}{\∂a^4}+p^4\frac{{\∂}^{4}T(\hat{s},\hat{t})}{\∂t^4}}{180}\\ J_1=T(s_{\mathrm{0,0}},t_{\mathrm{0,0}})+2\underset{i=1}{\overset{l-1}{\mathrm{\Σ}}}T(s_{2i,0},t_{2i,0})+4\underset{i=1}{\overset{l}{\mathrm{\Σ}}}T(s_{2i-\mathrm{1,0}},t_{2i-\mathrm{1,0}})+T(s_{2l,0},t_{2l,0})\\ J_2=\underset{j=1}{\overset{k-1}{\mathrm{\Σ}}}T(s_{0,2j},t_{0,2j})+2\underset{j=1}{\overset{k-1}{\mathrm{\Σ}}}\underset{i=1}{\overset{l-1}{\mathrm{\Σ}}}T(s_{2i,2j},t_{2i,2j})+\\ 4\underset{j=1}{\overset{k-1}{\mathrm{\Σ}}}\underset{i=1}{\overset{l}{\mathrm{\Σ}}}T(s_{2i-1,2j},t_{2i-1,2j})+\underset{j=1}{\overset{k-1}{\mathrm{\Σ}}}T(s_{2l,2j},t_{2l,2j})\\ J_3=\underset{j=1}{\overset{k}{\mathrm{\Σ}}}T(s_{0,2j-1},t_{0,2j-1})+2\underset{j=1}{\overset{k}{\mathrm{\Σ}}}\underset{i=1}{\overset{l-1}{\mathrm{\Σ}}}T(s_{2i,2j-1},t_{2i,2j-1})+\\ 4\underset{j=1}{\overset{k}{\mathrm{\Σ}}}\underset{i=1}{\overset{l}{\mathrm{\Σ}}}T(s_{2i-1,2j-1},t_{2i-1,2j-1})+\underset{j=1}{\overset{k}{\mathrm{\Σ}}}T(s_{2l,2j-1},t_{2l,2j-1})\\ J_4=T(s_{0,2k},t_{0,2k})+2\underset{i=1}{\overset{l-1}{\mathrm{\Σ}}}T(s_{2i,2k},t_{2i,2k})+4\underset{i=1}{\overset{l}{\mathrm{\Σ}}}T(s_{2i-1,2k},t_{2i-1,2k})+T(s_{2l,2k},t_{2l,2k})\end{array}\right.$

我们可以通过增加采样密度l和k来减少E_error项，当l，k均为10时，误差项就变成了如下形式

$E_{\mathrm{error}}=-\frac{\frac{{\∂}^{4}T(\stackrel{\‾}{u},\stackrel{\‾}{v})}{\∂u^4}+\frac{{\∂}^{4}T(\hat{u},\hat{v})}{\∂v^4}}{2.88\×{10}^{11}}$

对于一个给定的NURBS曲面，误差计量是有上界的，当l，k为10时，误差项能够减少7个数量级。等式(5)中得能量函数被离散为如下形式

${\tilde{J}}_c\left(Q_{i,j}\right),i=0...n_s,j=0...n_t.$

除了保角性，我们希望变换可以逼近X^-1οM^-1，这样离散形式

${\tilde{J}}_p\left(Q_{i,j}\right)=\underset{i=0}{\overset{2l}{\mathrm{\Σ}}}\underset{j=0}{\overset{2k}{\mathrm{\Σ}}}{|R_{i,j}-\mathrm{\φ}(s_{i,j},t_{i,j})|}^2$

就可以尽可能地小了。在下一节中给出了一个就加权带约束最小二乘方法用来求最佳变换φ_o，使得层次自由曲面X(φ_o(s，t))收敛于具有保角性。

S3：采用最小二乘方法极小化离散后的保角能量：

为了具有极小化保角能采用了一个加权带约束的最小二乘方法极小化保角能量。首先我们重新描述了这个极小化问题，然后通过拉格朗日乘子法得出了它的解。

假设{R(i，j)}和{(s_i，j，t_i，j)}，i＝0，…，2l，j＝0，…，2k分别是一系列点和要满足的保角性约束。我们将{R(i，j)}和{C(s_i，j，t_i，j)}分为无约束部分和有约束部分。使用和来表示无约束的部分，使用和来表示有约束的部分(即精确插值部分)。我们同行也允许给无约束的部分加上约束，这样我们得到了和通过这些权重，我们增加了对每个数据点逼近其邻近点的逼近程度的额外控制。增加权重可以增加对应项的逼近程度，减小权重可以减弱对应项的逼近程度。

假定m_u＝r_u+s_u+1，m_c＝r_c+s_c+1，我们希望使用一个具有(n_s+1)(n_t+1)个控制顶点的自由变换φ(s，t)来对无约束点逼近，对约束点插值，这里我们要去m_c＜n且m_c+n＜m_u+1。

S_k，k＝0，…，m_u是第k个无约束数据点(一个或者C^u(k)(s_i，j，t_i，j))；

T_k，k＝0，…，m_c是第k个有约束数据点(一个或者C^c(k)(s_i，j，t_i，j))；

ω_k，k＝0，…，mu是第k个无约束数据点的约束项。

定义向量和矩阵如下：

S＝[S_k](一个有m_u+1个元素的向量)；

T＝[T_k](一个有m_c+1个元素的向量)；

W＝[ω_k](一个(m_u+1)×(m_u+1)的对角矩阵，ω_k在对角线上)；

Q＝[Q_i，j](一个(n_s+1)×(n_t+1)的未知控制定点向量)；

N＝[ND_k](其中ND_k是一个k次基函数或者是保角项的基函数，它是一个(m_u+1)×(n_s+1)(n_t+1)的标量矩阵)；

M＝[MD_k](其中MD_k是一个k次基函数或者是保角项的基函数，它是一个(m_u+1)×(n_s+1)(n_t+1)的标量矩阵)；

参数s_i，j和t_i，j可以通过里奇流方法计算得到。自由变换变成了一个基本的带约束、有(n_s+1)×(n_t+1)个未知数、极小化问题，m_c+1个约束的求解极小化问题。基本算法是使用拉格朗日乘子法进行计算。拉格朗日乘子法通过引入m_c+1个附加未知量λ_k(拉格朗日乘子)。因此我们得到了[(n_s+1)(n_t+1)+m_c+1]×[(n_s+1)(n_t+1)+m_c+1]规模的等式矩阵。首先计算出来(m_c+1)个λ_k，然后在计算出(n_s+1)(n_t+1)个Q_i，j。

令

A＝[λ_k]

是兰格朗日乘子的向量，注意每个λ_k都是和Q_i，j维度相同的向量。

我们希望极小化下面的式子

（S^T-PN^T）W（S–NP）+A^T(MP-T)

对于未知数A，P。对上面的式子进行微分，令微分为0

-2(S^TWN-P^TN^TWN)+A^TM＝0

MP-T＝0

表示为矩阵形式为

$\left[\begin{array}{cc}{N}^T\mathrm{WN}& M^T\\ M& 0\end{array}\right]\left[\begin{array}{c}P\\ A\end{array}\right]=\left[\begin{array}{c}{N}^T\mathrm{WS}\\ T\end{array}\right]$

通过求解上式(我们可得到A和P。解唯一的条件是N^TWN和M(N^TWN)^-1M^T都是可逆的。可得

$\left\{\begin{array}{c}P={\left(N^T\mathrm{WN}\right)}^{-1}{N}^T\mathrm{WS}-{\left(N^T\mathrm{WN}\right)}^{-1}{M}^{T}A\\ A={\left(M{\left(N^T\mathrm{WN}\right)}^{-1}{M}^T\right)}^{-1}({M\left(N^T\mathrm{WN}\right)}^{-1}{N}^T\mathrm{WS}-T)\end{array}\right.$

最终P为

P＝(N^TWN)^-1N^TWS–(N^TWN)^-1M^T(M(N^TWN)^-1M^T)^-1(M(N^TWN)^-1N^TWS-T)

S4:与指定精度比较，确定是否继续迭代：

层次自由曲面的保角性可以通过在(s，t)参数域均匀采样，计算采样点的最大的保角性误差度量。

$\underset{0\≤k\≤m}{\max }{|N\×\frac{\∂X}{\∂s}-\frac{\∂X}{\∂t}|}^2$

如果最大的保角误差大于给定的阙值ε_c，则提高(u，v)参数域采样精度，迭代执行S1。否则迭代终止。通过反复迭代，我们最终可以得到一个具有保角性的层次自由曲面。图2为迭代终止曲面。

图3是整个方法的步骤流程图，初始化阶段输入一张NURBS空间曲面S；重定义自由变换形式，提高(u，v)参数域采样精度；使用最小二乘法计算最优变换；计算新的层次自由曲面保角性误差和改善程度；如果最大的保角性误差小于阙值ε_c，或者改善程度小于阙值ε_i，迭代终止；否则层次自由曲面仍不具有保角性，则重新定义自由变换形式，迭代继续

在初始步骤中，设定l₀，k₀，这时保角性误差大于ε_c的点最多。然后提取出iso-等参线作为参数s和t，为它们分别设定ε_c。ε_c可以用来产生新的节点向量。新的节点向量将用于下面的曲面优化步骤。反复迭代，我们最终可以得到一个具有保角性的层次自由曲面

上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。

Claims (5)

1.一种层次自由曲面的保角参数特性优化方法，包括如下步骤：

1)输入NURBS空间自由曲面X(u，v)；

2)将自由曲面X(u，v)的保角能量公式化，得到保角能量函数；并对保角能量函数进行离散化处理；

3)采用最小二乘方法极小化离散后的保角能量函数，并计算出相应的优化后的自由曲面；

4)对新的层次自由曲面采样，计算保角性误差，将所得保角性误差与设定的阈值进行比较：如果所述误差小于设定阈值，则所得层次自由曲面具有保角性；如果所述误差大于设定阈值，则所得层次自由曲面不具有保角性，分别增加u，v上采样点数目，提高(u，v)参数域采样精度，返回步骤1)重新计算；

所述步骤2)中将自由曲面X(u，v)的保角能量公式化，得到保角能量函数的具体为：

$J_c\left(Q_{i,j}\right)={\∫}_0^1{\∫}_0^1|{\mathrm{\φ}}_t-\left(\begin{array}{cc}-\frac{F}{\sqrt{EG-F^2}}& -\frac{G}{\sqrt{EG-F^2}}\\ \frac{E}{\sqrt{EG-F^2}}& \frac{F}{\sqrt{EG-F^2}}\end{array}\right){\mathrm{\φ}}_s{|}^{2}dsdt$

其中，确定曲面X(u，v)的基本形式为：

ds²＝X_u·X_v(du)²+2X_u·X_vdudv+X_v·X_v(dv)²

为曲面X(u，v)的两个偏导数；

E＝X_u·X_v，F＝X_u·X_v，G＝X_v·X_v

其中，E和G代表了两个偏导的长度，F用来度量两个偏导的正交性；Q_i，j是二维变换φ的控制顶点；

由层次自由曲面变换可得

$X(u,v)=\mathrm{\φ}(s,t)={\Σ}_{i=0}^{n_s}{\Σ}_{j=0}^{n_t}{N}_i^{p_2}\left(s\right)N_j^{q_2}\left(t\right)Q_{i,j},$

上述保角能量函数中φ_s和φ_t分别是变换后自由曲面对s和t的偏导数。

2.如权利要求1所述的一种层次自由曲面的保角参数特性优化方法，所述步骤1)中NURBS空间自由曲面X(u，v)具体为：

$X(u,v)=\frac{{\Σ}_{i=0}^{n_u}{\Σ}_{j=0}^{n_v}{N}_i^{p_1}\left(u\right)N_j^{q_1}\left(v\right){\mathrm{\ω}}_{i,j}{P}_{i,j}}{{\Σ}_{i=0}^{n_u}{\Σ}_{j=0}^{n_v}{N}_i^{p_1}\left(u\right)N_j^{q_1}\left(v\right){\mathrm{\ω}}_{i,j}},u,v\∈\[0,1\]$

其中，u和v是曲面X(u，v)的参数，n_u和n_v是曲面X(u，v)所含的控制顶点在u方向和v方向上的个数，P_i，j是曲面X(u，v)中序号为(i，j)的控制顶点的坐标值，由用户输入，i＝0，1，...，n_u-1；j＝0，1，...，n_v-1，采用三维直角坐标表示，ω_i，j是曲面X(u，v)中序号为(i，j)的控制顶点的权重，由用户输入，i＝0，1，...，n_u-1；j＝0，1，...，n_v-1，是定义在节点向量序列U之上的序号为i、次数为p₁的B样条基函数，i＝0，1，...，n_u-1；是定义在节点向量序列V之上的序号为j、次数为q₁的B样条基函数，j＝0，1，...，n_v-1，p1、q1分别为B样条基函数的阶数；

其中，p，q是曲面X(u，v)在u方向和v方向上的次数，由用户输入，u₀，u₁，...，u_n+p等为节点向量序列U中的节点，其中u₀＝u₁＝...＝u_p＝0，u_n＝u_n+1＝...＝u_n+p＝1，u_p+1，u_n+2，...，u_n-1由用户输入；v₀，v₁，...，v_n+q为节点向量序列V中的节点，其中v₀＝v₁＝...＝v_q＝0，v_n＝v_n+1＝...＝v_n+q＝1，v_q+1，v_q+2，...，v_n-1由用户输入，上述的B样条基函数采用Matlab系统样条工具库中B样条基函数的定义方式。

3.如权利要求1所述的一种层次自由曲面的保角参数特性优化方法，所述步骤2)中对保角能量函数利用复合辛普森规则逼近，为了应用该规则，对保角能量函数进行离散化处理的具体方法如下：

在曲面X(u，v)的参数域上进行离散采样，分别将参数域u，v细分为2l和2k个子区间，l和k分别为在参数域u，v上的采样密度，其均分步长分别为h＝1/(21)，p＝1/(2k)；细分后参数域如下：

$R_{i,j}=\left(\begin{array}{c}{u}_p+\frac{i\·(u_{n_u+1}-u_p)}{2l}\\ v_q+\frac{i\·(v_{n_v+1}-v_q)}{2k}\end{array}\right),i=0,...2l,j=0,...2k.$

其中，u_p、v_q、均为参数域节点向量中的节点；

使用里奇流方法计算出对应离散化后曲面X(u，v)对应的保角映射M，令

$T=|{\mathrm{\φ}}_t-\left(\begin{array}{cc}-\frac{F}{\sqrt{EG-F^2}}& -\frac{G}{\sqrt{EG-F^2}}\\ \frac{E}{\sqrt{EG-F^2}}& \frac{F}{\sqrt{EG-F^2}}\end{array}\right){\mathrm{\φ}}_s{|}^2$

其中，E＝X_u·X_v，F＝X_u·X_v，G＝X_v·X_v，E和G代表了两个偏导的长度，F用来度量两个偏导的正交性；φ_s和φ_t分别是层次自由曲面变换(u，v)＝φ(s，t)后，对s和t的偏导数；利用辛普森规则的逼近后的保角能量结果应有如下形式：

$J=\tilde{J}+E_{error}$

其中，参数域中存在某一个和使得

$\left\{\begin{array}{c}\tilde{J}=\frac{4hp}{9}\left\{J_1+2J_2+4J_3+J_4\right\}\\ E_{error}=-\frac{h^4\frac{{\∂}^{4}T\left(\stackrel{\‾}{s},\stackrel{\‾}{t}\right)}{\∂s^4}+p^4\frac{{\∂}^{4}T\left(\hat{s},\hat{t}\right)}{\∂t^4}}{180}\\ J_1=T\left(s_{0,0},t_{0,0}\right)+2\underset{i=1}{\overset{l-1}{\Σ}}T\left(s_{2i,0},t_{2i,0}\right)+4\underset{i=1}{\overset{l}{\Σ}}T\left(s_{2i-1,0},t_{2i-1,0}\right)+T\left(s_{2l,0},t_{2l,0}\right)\\ J_2=\underset{j=1}{\overset{k-1}{\Σ}}T\left(s_{0,2j},t_{0,2j}\right)+2\underset{j=1}{\overset{k-1}{\Σ}}\underset{i=1}{\overset{l-1}{\Σ}}T\left(s_{2i,2j},t_{2i,2j}\right)+\\ 4\underset{j=1}{\overset{k-1}{\Σ}}\underset{i=1}{\overset{l}{\Σ}}T\left(s_{2i-1,2j},t_{2i-1,2j}\right)+\underset{j=1}{\overset{k-1}{\Σ}}T\left(s_{2l,2j},t_{2l,2j}\right)\\ J_3=\underset{j=1}{\overset{k}{\Σ}}T\left(s_{0,2j-1},t_{0,2j-1}\right)+2\underset{j=1}{\overset{k}{\Σ}}\underset{i=1}{\overset{l-1}{\Σ}}T\left(s_{2i,2j-1},t_{2i,2j-1}\right)+\\ 4\underset{j=1}{\overset{k}{\Σ}}\underset{i=1}{\overset{l}{\Σ}}T\left(s_{2i-1,2j-1},t_{2i-1,2j-1}\right)+\underset{j=1}{\overset{k}{\Σ}}T\left(s_{2l,2j-1},t_{2l,2j-1}\right)\\ J_4=T\left(s_{0,2k},t_{0,2k}\right)+2\underset{i=1}{\overset{l-1}{\Σ}}T\left(s_{2i,2k},t_{2i,2k}\right)+4\underset{i=1}{\overset{l}{\Σ}}T\left(s_{2i-1,2k},t_{2i-1,2k}\right)+T\left(s_{2l,2k},t_{2l,2k}\right)\end{array}\right.$

其中，J₁，J₂，J₃，J₄是其四个离散项，s_i，j，t_i，j是参数中的采样点u，v坐标，h和p分别是u，v方向的采样间隔，T是保角项；通过增加采样密度l和k来减少E_error项，当l和k均为10时，误差项就变成了如下形式

$E_{error}=-\frac{\frac{{\∂}^{4}T(\stackrel{\‾}{u},\stackrel{\‾}{v})}{\∂u^4}+\frac{{\∂}^{4}T\left(\hat{u},\hat{v}\right)}{\∂v^4}}{2.88\×{10}^{11}}$

其中u，v是自由曲面的两个参数，T是保角项；优化其如下离散保角能量：

${\tilde{J}}_p\left(Q_{i,j}\right)=\underset{i=0}{\overset{2l}{\Σ}}\underset{j=0}{\overset{2k}{\Σ}}|R_{i,j}-\mathrm{\φ}(s_{i,j},t_{i,j}){|}^2$

R_i，j是u，v参数域中的采样点，s_i，j，t_i，j是st参数域中的两个参数；φ是二维自由曲面变换，Q_i，j是二维变换φ的控制顶点；优化离散保角能量使得生成的层次曲面的参数化逼近保角性。

4.如权利要求1所述的一种层次自由曲面的保角参数特性优化方法，所述步骤3)中采用最小二乘方法极小化离散后的保角能量函数的具体方法为：

令{R(i，j)}和{C(s_i，j，t_i，j)}，i＝0，...，2l，j＝0，...，2k分别表示采样点集合和其对应满足的保角性约束条件集合；

将{R(i，j)}和{C(s_i，j，t_i，j)}分为无约束部分和有约束部分：使用和来表示无约束的部分；使用和来表示有约束的部分；给有约束部分设定非负权重，有和增加权重可以增加每个数据点逼近其邻近点的逼近程度，减小权重可以减弱每个数据点逼近其邻近点的逼近程度；

令m_u＝r_u+s_u+1，m_c＝r_c+s_c+1，其中m_u，r_u，s_u表示无约束部分元素的数目，m_c，r_c，s_c表示有约束部分元素的数目，且满足m_c＜n，m_c+n＜m_u+1；

令

S_k，k＝0，...，m_u是第k个无约束数据点；

T_k，k＝0，...，m_c是第k个有约束数据点；

ω_k，k＝0，...，mu是第k个无约束数据点的约束项；

定义向量和矩阵如下：

S＝[S_k]，为一个有m_u+1个元素的向量；

T＝[T_k]，为一个有m_c+1个元素的向量；

W＝[ω_k]，为一个(m_u+1)×(m_u+1)的对角矩阵，ω_k在对角线上；

Q＝[Q_i，j]，为一个(n_s+1)×(n_t+1)的未知控制定点向量；

N＝[ND_k]，其中ND_k是一个k次基函数或者是保角项的基函数，N为一个(m_u+1)×(n_s+1)(n_t+1)的标量矩阵；

M＝[MD_k]，其中MD_k是一个k次基函数或者是保角项的基函数，M为一个(m_u+1)×(n_s+1)(n_t+1)的标量矩阵；

令A＝[λ_k]，A是拉格朗日乘子的向量，每个λ_k都是和Q_i，j维度相同的向量；这里，无约束部分有

NQ＝S

有约束部分有

MQ＝T

无约束部分误差是S-NP，通过极小化这个误差从而使得MP＝T；因此，利用拉格朗日乘子法，等价极小化

(S^T-PN^T)W(S-NP)+A^T(MP-T)

对于未知数A，P；对上面的式子进行微分，令微分为0

-2(S^TWN-p^TN^TWN)+A^TM＝0

MP-T＝0

表示为矩阵形式为

$\left[\begin{array}{cc}{N}^{T}WN& M^T\\ M& 0\end{array}\right]\left[\begin{array}{c}P\\ A\end{array}\right]=\left[\begin{array}{c}{N}^{T}WS\\ T\end{array}\right]$

通过求解上式得到A和P；解唯一的条件是N^TWN和M(N^TWN)^-1M^T都是可逆的；可得

$\left\{\begin{array}{c}P={\left(N^{T}WN\right)}^{-1}{N}^{T}WS-{\left(N^{T}WN\right)}^{-1}{M}^{T}A\\ A={\left(M{\left(N^{T}WN\right)}^{-1}{M}^T\right)}^{-1}(M{\left(N^{T}WN\right)}^{-1}{N}^{T}WS-T)\end{array}\right.$

最终P为

P＝(N^TWN)-¹N^TWS-(N^TWN)-¹M^T(M(N^TWN)-¹M^T)-¹(M(N^TWN)-¹N^TWS-T)该P就是通过最小二乘法得到的满足要求的解。

5.如权利要求1所述的一种层次自由曲面的保角参数特性优化方法，所述步骤4)中新的层次自由曲面采样点的最大的保角性误差为：

$\underset{0\≤k\≤m}{max}|N\×\frac{\∂X}{\∂s}-\frac{\∂X}{\∂t}{|}^2$

其中，N是曲面法向量，分别是曲面X在s，t两个方向上的切向量。