矩阵权NURBS曲面的拟合与光顺

doi:10.3785/j.issn.1008-9497.2023.06.017

矩阵权NURBS曲面的拟合与光顺

董国鑫^,, 杨勋年^,^,

浙江大学数学科学学院，浙江杭州 310058

Fitting and fairing quad-meshes by matrix weighted NURBS surfaces

DONG Guoxin^,, YANG Xunnian^,^,

School of Mathematical Sciences，Zhejiang University，Hangzhou 310058，China

通讯作者: ORCID：https：//orcid.org/0000-0001-7501-1378，E-mail：yxn@zju.edu.cn.

收稿日期: 2023-06-21 修回日期: 2023-07-19 接受日期: 2023-07-26

基金资助:

国家自然科学基金资助项目. 12171429

Corresponding authors: zdxb_l@zju.edu.cn

Received: 2023-06-21 Revised: 2023-07-19 Accepted: 2023-07-26

作者简介 About authors

董国鑫（1998—），ORCID：https：//orcid.org/0000-0001-9112-4022，男，硕士研究生，主要从事计算机辅助几何设计研究. 。

摘要

介绍了利用矩阵权非均匀有理 B-spline（non-uniform rational B-spline，NURBS）曲面拟合与光顺带法向的网格数据。通过输入带法向的规则四边形网格，以网格的顶点为曲面的控制顶点，用法向信息计算每个控制顶点对应的矩阵权，构造矩阵权NURBS曲面。与传统NURBS曲面相比，矩阵权NURBS曲面具有拟柱面精度，当数据均匀采样自光滑曲面时，构造的矩阵权NURBS曲面具有较好的光顺性且能很好地拟合网格模型；当输入的网格数据掺杂噪声时，通过迭代在已有的矩阵权NURBS曲面上重新采样顶点和计算法向量构造新的矩阵权NURBS曲面，最后得到具有较好光顺性且仍能逼近原始网格数据的拟合曲面。

关键词： 矩阵权NURBS ; 曲面拟合 ; 光顺

Abstract

This paper proposes to employ matrix weighted NURBS surfaces to fit and fair quad meshes. For a quadrilateral mesh with given or estimated unit normals at vertices, a matrix weighted NURBS surface can be constructed by choosing the mesh vertices as control points and employing normals at each vertex for computing matrix weights. Compared with traditional NURBS surfaces, matrix weighted NURBS surfaces have quasi-cylindrical accuracy. When the input data is uniformly sampled from a smooth surface, the constructed matrix weighted NURBS surface has good smoothness and fits the mesh model well; if the input grid data contain noise, a fair fitting surface that approximates the original grid well can still be obtained by resampling control vertices on current fitting surfaces and re-calculating vertex normals based on the new quad meshes iteratively.

Keywords： matrix weighted NURBS ; surface fitting ; fairing

PDF (1541KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

董国鑫, 杨勋年. 矩阵权NURBS曲面的拟合与光顺. 浙江大学学报(理学版)[J], 2023, 50(6): 820-828 doi:10.3785/j.issn.1008-9497.2023.06.017

DONG Guoxin, YANG Xunnian. Fitting and fairing quad-meshes by matrix weighted NURBS surfaces. Journal of Zhejiang University(Science Edition)[J], 2023, 50(6): 820-828 doi:10.3785/j.issn.1008-9497.2023.06.017

0　引言

曲面拟合与光顺在几何建模、计算机辅助设计、三维物体重建、数控加工等领域有广泛应用^［1-4］。通过参数曲面拟合一系列离散点，例如非均匀有理B样条（non-uniform rational B-spline，NURBS）曲面拟合，可在节省存储空间的同时获得较为光顺的曲面。本文主要研究用矩阵权NURBS曲面拟合带法向的规则四边形网格数据。矩阵权NURBS曲面是NURBS曲面的自然的推广^［5］，与传统NURBS曲面相比，矩阵权NURBS曲面具有更多的自由度控制与设计曲面形状。通过合理定义矩阵权并调整矩阵权中的参数，使得构造的矩阵权NURBS曲面能很好地逼近给定数据且曲面具有很好的光顺性，避免了传统B样条曲面拟合需要求解线性方程组的问题。

曲面拟合是指通过构造特定曲面逼近或插值给定的几何数据集，B样条曲面拟合主要分以下几步^［1］：首先根据输入的数据点计算得到合适的参数，然后根据参数生成B样条的节点向量，最后通过构造线性方程组求解B样条曲面的控制顶点。目前已有不少对B样条曲面拟合的研究。ROGERS等^［6］研究了带约束的B样条曲线和曲面的拟合问题；MA等^［7］研究了逆向工程中NURBS曲线和曲面的拟合，采用两步线性方法拟合NURBS曲线和曲面；LIN等^［8］提出在不求解线性方程组的情况下迭代构造有序点集的拟合曲面；SAPPA等^［9］提出基于隐式B样条的方法对点云数据进行曲面重建，用隐式B样条表示可避免参数化问题；LIU等^［10］提出了一种新的渐进迭代逼近方法，并将其用于B样条曲面拟合，与通常的拟合方法相比，该方法能得到低次数的B样条曲面和有效平滑的噪声数据，但求解效率不高。JIANG等^［11］通过迭代的节点插入算法构造插值且全局光滑的B样条曲面，构造的曲面具有较少的控制顶点。

曲面光顺^［12］是指对曲面进行平滑处理，其目的是获得尽可能光滑的曲面形状。如何度量曲面的光顺性取决于具体的应用，一般来说光顺曲面具有简单的形状，且不含多余的细节和振荡。通常光顺曲面通过求解包含光顺项的函数最小化构造^［13］，例如LOTT等^［14］构造了关于曲面曲率和距离的目标函数，基于最小化约束求解控制顶点得到光顺的B样条曲面；GREINER^［15］提出通过简化的薄板能量泛函构造光顺曲面；DIETZ^［16］通过模拟金属薄板的拉伸过程，结合薄板能量对离散的点云数据进行计算得到光顺的逼近曲面。除这一类曲面光顺方法外，也有基于几何迭代的拟合与光顺方法，即通过几何信息迭代更新控制顶点等数据。MAEKAWA等^［17］提出了一种用几何算法构造光滑的曲面插值给定的网格；KINERI等^［18］提出了一种基于几何迭代的算法，得到拟合点云数据的B样条曲面；KINERI等^［19］研究了基于曲率编辑的曲面设计方法，根据几何规则调整控制顶点位置，使得B样条曲面能够插值特定位置和节点处的曲率。

YANG^［5］将传统的NURBS推广至矩阵权NURBS，通过输入的顶点和法向数据构造可插值控制顶点的曲线和曲面。与传统的NURBS曲面相比，矩阵权NURBS曲面具有更优的性质和拟柱面再生等优点。对于光滑曲面，可以在曲面上合理采样控制顶点和法向信息并利用矩阵权NURBS曲面拟合该曲面。此外，对于同样的控制顶点和法向量，高阶矩阵权NURBS曲面的光顺性通常比低阶矩阵权NURBS曲面的更好，甚至高阶矩阵权NURBS曲面可有效地光顺噪声数据。由于矩阵权NURBS曲面丰富的性质，本文拟利用矩阵权NURBS曲面的构造与表达对规则网格数据曲面进行拟合与光顺。

1　矩阵权NURBS曲面

假设给定控制顶点 ${P_{i, j}}$ ， $0 \leq i \leq m$ ， $0 \leq j \leq n$ ，节点向量分别设为 $U = {u_{0}, u_{1}, \dots, u_{m + k}}$ 和 $V = {v_{0}, v_{1}, \dots, v_{n + l}}$ ，传统的NURBS曲面定义为

S (u, v) = \frac{\sum_{i = 0}^{m} \sum_{j = 0}^{n} w_{i, j} P_{i, j} N_{i, k} (u) N_{j, l} (v)}{\sum_{i = 0}^{m} \sum_{j = 0}^{n} w_{i, j} N_{i, k} (u) N_{j, l} (v)}

，

其中， $N_{i, k} (u)$ 表示k阶B样条基函数：

N_{i, 1} (u) = \{\begin{array}{l} 1, u_{i} \leq u < u_{i + 1} ， \\ 0, 其他 ； \end{array}

\begin{array}{l} N_{i, k} (u) = \frac{u - u_{i}}{u_{i + k - 1} - u_{i}} N_{i, k - 1} (u) + \\ \frac{u_{i + k} - u}{u_{i + k} - u_{i + 1}} N_{i + 1, k - 1} (u), k \geq 2 。 \end{array}

作为传统NURBS的推广，YANG^［5］提出了矩阵权NURBS，即将NURBS表达式中的实数权系数推广为矩阵。矩阵权NURBS曲面定义为

S (u, v) = {[\sum_{i = 0}^{m} \sum_{j = 0}^{n} M_{i, j} N_{i, k} (u) N_{j, l} (v)]}^{- 1} \times \sum_{i = 0}^{m} \sum_{j = 0}^{n} M_{i, j} P_{i, j} N_{i, k} (u) N_{j, l} (v),

（1）

其中， $u \in [u_{k - 1}, u_{m + 1}],$ $v \in [v_{l - 1}, v_{n + 1}],$ $M_{i, j} \in$ $R^{d \times d}$ $（ 0 \leq i \leq m,$ $0 \leq j \leq n ）$ 为给定的一系列权矩阵。矩阵权可以有不同的选取方法，但首先要保证曲面 $S (u, v)$ 定义的合理性，即 $\sum_{i = 0}^{m} \sum_{j = 0}^{n} M_{i, j} N_{i, k} (u) N_{j, l} (v)$ 为可逆矩阵。记 $I$ 为3阶单位矩阵， $n_{i, j} \in R^{3}$ 为单位向量，矩阵权定义为

\begin{array}{l} M_{i, j} = ω_{i, j} (I + μ_{i, j} n_{i, j} n_{i, j}^{T}), ω_{i, j} > 0, μ_{i, j} > - 1, \\ i = 0,1, \dots, m, j = 0,1, \dots, n 。 \end{array}

（2）

这样选取的矩阵权能够保证 $M (u, v) = \sum_{i = 0}^{m} \sum_{j = 0}^{n} M_{i, j} N_{i, k} (u) N_{j, l} (v)$ 为可逆矩阵^［20］，使矩阵权NURBS曲面是良定义的。

定理1　设 $n_{i, j}$ 为单位向量，矩阵 $M_{i, j} = ω_{i, j} (I + μ_{i, j} n_{i, j} n_{i, j}^{T}),$ $0 \leq i \leq m, 0 \leq j \leq n,$ 对任意的正实数 $s_{i, j}, 0 \leq i \leq m, 0 \leq j \leq n,$ 当 $ω_{i, j} > 0$ 且 $μ_{i, j} > - 1$ 时， $\sum_{i = 0}^{m} \sum_{j = 0}^{n} s_{i, j} M_{i, j}$ 为可逆矩阵。

证明　由矩阵 $M_{i, j}$ 的形式，可知 $M_{i, j}$ 为对称矩阵。下证 $M_{i, j}$ 为正定矩阵。为此只需考虑矩阵 $I + μ_{i, j} n_{i, j} n_{i, j}^{T}$ 的特征值。由于该矩阵的特征多项式为

| λ I - (I + μ_{i, j} n_{i, j} n_{i, j}^{T}) | = | (λ - 1) I - μ_{i, j} n_{i, j} n_{i, j}^{T} |

，

因此，可知该矩阵的特征值为1或 $1 + μ_{i, j}$ ，均大于0，所以 $M_{i, j}$ 为正定矩阵，又因为 $s_{i, j} > 0$ ，因此 $\sum_{i = 0}^{m} \sum_{j = 0}^{n} s_{i, j} M_{i, j}$ 为正定矩阵，故其为可逆矩阵。证毕。

假设控制顶点 $P_{i, j}$ 处对应法向量 $n_{i, j}$ 的切平面为 $F_{i, j}$ ，对于空间中任意一点 $Q$ ，假设 $Q_{i, j}$ 为 $Q$ 点在 $F_{i, j}$ 上的投影，于是从 $Q_{i, j}$ 到 $Q$ 的向量可表示为

Q - Q_{i, j} = (n_{i, j} n_{i, j}^{T}) (Q - P_{i, j}) = A_{i, j} (Q - P_{i, j})

。

容易发现 $A_{i, j}^{T} = A_{i, j}$ 以及 $A_{i, j}^{2} = A_{i, j}$ ，于是从 $Q_{i, j}$ 到 $Q$ 的距离为

| | Q - Q_{i, j} | |^{2} = [Q - P_{i, j}]^{T} A_{i, j}^{2} [Q - P_{i, j}]

。

下面构造能量函数：

\begin{array}{l} F (Q (u, v)) = \sum_{i = 0}^{m} \sum_{j = 0}^{n} ω_{i, j} N_{i, k} (u) N_{j, l} (v) [Q - P_{i, j}]^{2} + \\ \sum_{i = 0}^{m} \sum_{j = 0}^{n} ω_{i, j} μ_{i, j} N_{i, k} (u) N_{j, l} (v) \times \\ [Q - P_{i, j}]^{T} A_{i, j}^{2} [Q - P_{i, j}] ， \end{array}

（3）

可以证明 $F (Q (u, v))$ 有唯一的最小二乘解，且该解就是所构造的矩阵权NURBS曲面。

性质1　当参数 $ω_{i, j} > 0$ 且 $μ_{i, j} > - 1$ 时，能量函数式（3）取最小值时的解 $Q (u, v)$ 为矩阵权NURBS曲面。

证明　首先证明 $F (Q (u, v)) \geq 0$ 。由于 $ω_{i, j} > 0, N_{i, k} (u) \geq 0, N_{j, l} (v) \geq 0$ ，只需证明

[Q - P_{i, j}]^{2} + μ_{i, j} [Q - P_{i, j}]^{T} A_{i, j}^{2} [Q - P_{i, j}] \geq 0

。

由于

\begin{array}{l} [Q - P_{i, j}]^{2} + μ_{i, j} [Q - P_{i, j}]^{T} A_{i, j}^{2} [Q - P_{i, j}] = \\ [Q - P_{i, j}]^{T} (I + μ_{i, j} A_{i, j}^{2}) [Q - P_{i, j}], \end{array}

而矩阵 $I + μ_{i, j} A_{i, j}^{2} = I + μ_{i, j} n_{i, j} n_{i, j}^{T}$ 为正定矩阵，因此上式非负，从而 $F (Q (u, v)) \geq 0$ ，且当 $F (Q)$ 取最小值时，有

\begin{array}{l} \frac{1}{2} \frac{\partial F (Q)}{\partial Q} = \sum_{i = 0}^{m} \sum_{j = 0}^{n} ω_{i, j} N_{i, k} (u) N_{j, l} (v) \times \\ [Q - P_{i, j}] + \sum_{i = 0}^{m} \sum_{j = 0}^{n} ω_{i, j} μ_{i, j} N_{i, k} (u) \times \\ N_{j, l} (v) A_{i, j}^{T} A_{i, j} [Q - P_{i, j}] = \\ \sum_{i = 0}^{m} \sum_{j = 0}^{n} ω_{i, j} N_{i, k} (u) N_{j, l} (v) {[Q - P_{i, j}] + \\ μ_{i, j} A_{i, j} [Q - P_{i, j}]} = 0 。 \end{array}

其解为

\begin{array}{l} Q (u, v) = {[\sum_{i = 0}^{m} \sum_{j = 0}^{n} M_{i, j} N_{i, k} (u) N_{j, l} (v)]}^{- 1} \times \\ \sum_{i = 0}^{m} \sum_{j = 0}^{n} M_{i, j} P_{i, j} N_{i, k} (u) N_{j, l} (v), \end{array}

（4）

其中，

\begin{array}{l} M_{i, j} = ω_{i, j} (I + μ_{i, j} n_{i, j} n_{i, j}^{T}), ω_{i, j} > 0, μ_{i, j} > - 1, \\ i = 0,1, \dots, m, j = 0,1, \dots, n 。 \end{array}

（5）

性质2　设 $Q (u, v)$ 为矩阵权NURBS曲面，若参数 $μ_{i, j}$ 趋向于无穷大，则曲面片上的点 $Q (u, v), u \in (u_{i}, u_{i + k}), v \in (v_{j}, v_{j + l})$ 几乎处于 $P_{i, j}$ 的切平面上；同样，若参数 $ω_{i, j}$ 趋向于无穷大，则曲面片上的点 $Q (u, v), u \in (u_{i}, u_{i + k}), v \in (v_{j}, v_{j + l})$ 趋向于 $P_{i, j}$ 。

证明　由式（3）可知，当参数 $μ_{i, j}$ 趋于无穷大时， $N_{i, k} (u) N_{j, l} (v) [Q - P_{i, j}]^{T} A_{i, j}^{2} [Q - P_{i, j}]$ 必须趋于0才能使能量函数取最小值，此时 $[Q - P_{i, j}]^{T} A_{i, j}^{2} [Q - P_{i, j}]$ 趋于0，即 $Q (u, v)$ ， $u \in (u_{i}, u_{i + k}),$ $v \in (v_{j}, v_{j + l})$ 趋向于 $P_{i, j}$ 的切平面；同理，当参数 $ω_{i, j}$ 趋向于无穷大时， $[Q - P_{i, j}]^{2}$ 趋于0，即曲面片上的点 $Q (u, v), u \in (u_{i}, u_{i + k}), v \in (v_{j}, v_{j + l})$ 趋向于 $P_{i, j}$ 。

由性质2可知，本文的拟合和光顺算法可根据具体需要调整相应参数，以满足特定点以及特定点处法向逼近的要求。例如若希望顶点 $P_{i, j}$ 处的法向量逼近 $n_{i, j}$ ，则可适当增大对应的 $μ_{i, j}$ 值。此外，类似于矩阵权NURBS曲线具有螺线精度，矩阵权NURBS曲面也有相似的性质。方便起见，下文只考虑B样条函数的阶数为偶数且B样条函数的节点向量均匀选取的情形，其他情形可类似讨论。

性质3　假设点集 ${P_{i, j}}$ 和法向量集合 ${n_{i, j}}, i = 0,1, \dots ， m, j = 0,1, \dots, n$ ，是从柱面上均匀采样得到的，设 $d_{1} = \frac{k - 1}{2}, d_{2} = \frac{l - 1}{2}$ ，正整数 $q$ 和 $p$ 满足 $d_{1} \leq q \leq m - d_{1}, d_{2} \leq p \leq n - d_{2}$ ，令 $λ_{q + i} = N_{q + i, k} (u_{q + \frac{k}{2}})$ ， $ξ_{p + j} = N_{p + j, l} (v_{p + \frac{l}{2}})$ ，其中 $N_{i, k} (u)$ 和 $N_{j, l} (v)$ 是k阶和 $l$ 阶B样条基函数， $0 < | i | \leq d_{1}$ ， $0 < | j | \leq d_{2}$ ， $u_{q + \frac{k}{2}} = \frac{1}{2} (u_{q} + u_{q + k})$ ， $v_{p + \frac{l}{2}} = \frac{1}{2} (v_{p} + v_{p + l})$ 。参数

μ_{q, p} = \frac{\sum_{| i | = 0}^{d_{1}} \sum_{| j | = 0}^{d_{2}} λ_{q + i} ξ_{p + j} n_{q, p}^{T} (P_{q, p} - P_{q + i, p + j})}{\sum_{| i | = 0}^{d_{1}} \sum_{| j | = 0}^{d_{2}} λ_{q + i} ξ_{p + j} n_{q, p}^{T} A_{q + i, p + j} (P_{q + i, p + j} - P_{q, p})},

（6）

其中， $A_{q + i, p + j} = n_{q + i, p + j} n_{q + i, p + j}^{T}$ ，如果矩阵权中的参数 $w_{i, j} = 1, 0 \leq i \leq m, 0 \leq j \leq n$ ，且对于固定的 $q, p, μ_{q + i, p + j} = μ_{q, p}, | i | \leq d_{1}, | j | \leq d_{2}$ ，则定义的矩阵权NURBS曲面在参数 $(u_{q + \frac{k}{2}}, v_{p + \frac{l}{2}})$ 处取到顶点 $P_{q, p}$ 。

证明　由性质1的证明可知，在矩阵权有理B样条曲面 $Q (u, v)$ 上的任意一点均满足

\sum_{i = 0}^{m} \sum_{j = 0}^{n} ω_{i, j} N_{i, k} (u) N_{j, l} (v) [Q - P_{i, j}] + \sum_{i = 0}^{m} \sum_{j = 0}^{n} ω_{i, j} N_{i, k} (u) N_{j, l} (v) μ_{i, j} A_{i, j} [Q - P_{i, j}] = 0 。

假设 $Q (u, v)$ 在参数 $(u_{q + \frac{k}{2}}, v_{p + \frac{l}{2}})$ 处取到顶点 $P_{q, p}$ ，由条件 $w_{i, j} = 1,$ $0 \leq i \leq m, 0 \leq j \leq n,$ $μ_{q + i, p + j} = μ_{q, p},$ $| i | \leq d_{1}, | j | \leq d_{2}$ ，结合B样条基函数的性质，化简可得

\begin{array}{l} \sum_{| i | = 0}^{d_{1}} \sum_{| j | = 0}^{d_{2}} λ_{q + i} ξ_{p + j} [P_{q, p} - P_{q + i, p + j}] + \sum_{| i | = 0}^{d_{1}} \sum_{| j | = 0}^{d_{2}} λ_{q + i} ξ_{p + j} μ_{q, p} A_{q + i, p + j} \times \\ [P_{q, p} - P_{q + i, p + j}] = 0 \end{array}

，（7）

等式两边点乘 $n_{q, p}$ ，即可得到 $μ_{q, p}$ 的表达式。

下证当 $μ_{q, p}$ 的表达式为式（6）时， $Q (u, v)$ 在参数 $(u_{q + \frac{k}{2}}, v_{p + \frac{l}{2}})$ 处取到顶点 $P_{q, p}$ 。

首先，当 $μ_{q, p}$ 的表达式为式（6）时，有

\begin{array}{l} n_{q, p}^{T} \{\sum_{| j | = 0}^{d_{2}} λ_{q + i} ξ_{p + j} [P_{q, p} - P_{q + i, p + j}] + \sum_{| i | = 0}^{d_{1}} \sum_{| j | = 0}^{d_{2}} λ_{q + i} ξ_{p + j} μ_{q, p} A_{q + i, p + j} \times \\ [P_{q, p} - P_{q + i, p + j} \overset{\overset{}{}}{\underset{}{]}}\} = 0 。 \end{array}

（8）

令 $Q_{q + i, p + j}$ 为 $P_{q, p}$ 在 $P_{q + i, p + j}$ 的切平面上的投影点，则 $A_{q + i, p + j} (P_{q, p} - P_{q + i, p + j}) = P_{q, p} - Q_{q + i, p + j},$ 并且根据B样条基函数的性质， $λ_{q + i} = λ_{q - i}, ξ_{p + j} = ξ_{p - j}$ ，于是式（8）可化简为

\begin{array}{l} n_{q, p}^{T} \{\sum_{j = 0}^{d_{2}} λ_{q + i} ξ_{p + j} [4 P_{q, p} - P_{q + i, p + j} - \\ P_{q - i, p + j} - P_{q + i, p - j} - P_{q - i, p - j}] + \\ \sum_{i = 0}^{d_{1}} \sum_{j = 0}^{d_{2}} λ_{q + i} ξ_{p + j} μ_{q, p} [4 P_{q, p} - Q_{q + i, p + j} - \\ Q_{q - i, p + j} - Q_{q + i, p - j} - Q_{q - i, p - j} \overset{\overset{}{}}{\underset{}{]}}\} = 0 。 \end{array}

（9）

分别考虑式（9）的前后两项，以前一项为例，

\begin{array}{l} 4 P_{q, p} - P_{q + i, p + j} - P_{q - i, p + j} - P_{q + i, p - j} - P_{q - i, p - j} = \\ 2 (2 P_{q, p} - P_{q, p + j} - P_{q, p - j}) + (2 P_{q, p + j} - P_{q + i, p + j} - \\ P_{q - i, p + j}) + (2 P_{q, p - j} - P_{q + i, p - j} - P_{q - i, p - j}), \end{array}

由于 ${P_{i, j}}$ 和 ${n_{i, j}} (i = 0,1, \dots ， m, j = 0,1, \dots, n)$ 是从柱面上均匀采样得到的，不妨设 $P_{i, j} = (r c o s θ_{i}, r s i n θ_{i}, z_{j})$ ， $n_{i, j} = (c o s θ_{i}, s i n θ_{i}, 0)$ ，因此 $2 P_{q, p} - P_{q, p + j} - P_{q, p - j}, 2 P_{q, p} - P_{q + i, p} - P_{q - i, p},$ $2 P_{q, p} - Q_{q, p + j} - Q_{q, p - j}, 2 P_{q, p} - Q_{q + i, p} - Q_{q - i, p}$ 均平行于 $P_{q, p}$ 处的法向量 $n_{q, p}$ ，从而 $4 P_{q, p} - P_{q + i, p + j} - P_{q - i, p + j} - P_{q + i, p - j} - P_{q - i, p - j}$ 平行于 $n_{q, p}$ ，同理可得 $4 P_{q, p} - Q_{q + i, p + j} - Q_{q - i, p + j} - Q_{q + i, p - j} - Q_{q - i, p - j}$ 平行于 $n_{q, p}$ ，于是由式（8）可推得式（7），从而可知当 $μ_{q, p}$ 的表达式为式（6）时，有 $Q (u_{q + \frac{k}{2}}, v_{p + \frac{l}{2}}) = P_{q, p}$ 。得证。

图（1）为不同阶数的矩阵权NURBS曲面拟合加噪声的柱面网格。给定一个柱面，其参数表示为

\{\begin{array}{l} x (u, v) = c o s u, \\ y (u, v) = s i n u, u \in [0, π], v \in [0, π] 。 \\ z (u, v) = v, \end{array}

从该柱面上均匀地采样： $u_{0} = 0, u_{i} = u_{i - 1} + 0.05 π,$ $v_{0} = 0, v_{j} = v_{j - 1} + 0.05 π,$ $i, j = 1,2, \dots, 20$ ，随后添加随机函数产生高斯噪声，将这些采样点作为控制顶点分别构造不同次数的矩阵权NURBS曲面，其中矩阵权和参数 $μ_{i, j}$ 分别由式（2）和式（6）计算，曲面的颜色深浅代表该点处的Gauss曲率。从图1中可以看到，对于柱面采样网格，构造的矩阵权NURBS曲面经过所有的控制顶点，且曲面上所有点的Gauss曲率均较均匀；对于添加噪声的网格，双三次矩阵权NURBS曲面的Gauss曲率连续但不光滑，而提高次数后曲面的光顺性明显得到改善。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 矩阵权NURBS曲面拟合采样自柱面的网格

注数对（i， j）分别表示参数u，v方向的B样条基函数的次数。（a）和（b）为双三次矩阵权NURBS曲面分别拟合无噪声和有噪声的柱面网格；（c）和（d）分别为双五次和双七次矩阵权NURBS曲面的拟合结果。

Fig.1 The matrix weighted NURBS surface is used to fit the mesh sampled from a cylinder

Note The number pair （i， j） represents the degree of the B-spline basis function in the parameter u and v directions respectively. （a） and （b） are bicubic matrix weighted NURBS surfaces fitted to cylindrical meshes with and without noise， respectively；（c） and （d） are the fitting results of matrix weighted NURBS surfaces with bi-degree five or bi-degree seven， respectively.

2　算法实现

2.1　矩阵权NURBS曲面的拟合

假设 ${P_{i, j}}$ 和 ${n_{i, j}}$ $(0 \leq i \leq m, 0 \leq j \leq n)$ 分别为给定的采样顶点和法向量，由式（6）选取参数 $μ_{i, j}$ ，并且取 $ω_{i, j} = 1$ ，于是每个控制顶点对应的矩阵权 $M_{i, j}$ 均可直接计算，从而可以构造矩阵权NURBS曲面，但直接构造的矩阵权NURBS曲面不能拟合边界控制顶点，为了构造能拟合所有初始顶点的矩阵权NURBS曲面，需要在初始控制网格的边界添加边界对称点，获得扩展的控制网格。

假设要构造的矩阵权NURBS曲面的阶数为 $(k, l)$ ，设 $d_{1} = \frac{k - 1}{2}, d_{2} = \frac{l - 1}{2}$ ，对于固定的 $i, 0 \leq i \leq m$ ，添加 $P_{i, j}$ 关于 $P_{i, 0}$ 的对称点 $P_{i, - j}, j = 1,2, \dots, d_{2}$ ，其中，

\begin{matrix} P_{i, - j} = 2 P_{i, 0} - P_{i, j} + 2 a_{i, j} n_{i, 0}, \\ a_{i, j} = (P_{i, j} - P_{i, 0}) n_{i, 0}, j = 1,2, \dots, d_{2} 。 \end{matrix}

同样，添加 $P_{i, n - j}$ 关于 $P_{i, n}$ 的对称点 $P_{i, n + j}, j = 1,2, \dots, d_{2}$ ，类似地，这些对称点的单位法向量可通过计算对应点对称的法向量得到。因为控制网格的边界有4条，还需用同样方法对另外2条边界添加对称点，例如对于固定的 $j, - d_{2} \leq i \leq n + d_{2}$ ，添加 $P_{i, j}$ 关于 $P_{0, j}$ 的对称点 $P_{- i, j}, i = 1,2, \dots, d_{1}$ ，其中，

\begin{matrix} P_{- i, j} = 2 P_{0, j} - P_{i, j} + 2 a_{i, j} n_{0, j}, \\ a_{i, j} = (P_{i, j} - P_{0, j}) n_{0, j}, i = 1,2, \dots, d_{1} 。 \end{matrix}

利用 $(P_{i, j}, n_{i, j}),$ $i = - d_{1}, - d_{1} + 1, \dots,$ $m + d_{1},$ $j = - d_{2}, - d_{2} + 1, \dots,$ $n + d_{2}$ 构造矩阵权NURBS曲面

Q (u, v) = {[\sum_{i = - d_{1}}^{m + d_{1}} \sum_{j = - d_{2}}^{n + d_{2}} M_{i, j} N_{i, k} (u) N_{j, l} (v)]}^{- 1} \times \sum_{i = - d_{1}}^{m + d_{1}} \sum_{j = - d_{2}}^{n + d_{2}} M_{i, j} P_{i, j} N_{i, k} (u) N_{j, l} (v) 。

由性质3可知，由扩展后的控制网格构造的矩阵权NURBS曲面 $Q (u, v)$ 可以插值初始网格边界的控制顶点。为描述方便，对扩展后的顶点和法向量下标重新编号，分别记为 ${P_{i, j}^{'}}$ 和 ${{n_{i, j}}^{'}},$ $i = 0,1, \dots, m^{'},$ $j = 0,1, \dots, n^{'}$ ，其中， $m^{'} = m + 2 d_{1}, n^{'} = n + 2 d_{2}$ 。矩阵权NURBS曲面的拟合算法步骤如下：

输入　点集 ${P_{i, j}}$ 和法向量集 ${n_{i, j}}$ 以及 u， v 两个方向B样条基函数的阶数。

输出　矩阵权NURBS曲面。

Step1　计算并添加边界对称点，设添加后的点集和法向量集为 ${P_{i, j}^{'}}$ 和 ${{n_{i, j}}^{'}}$ ；

Step2　选取节点向量 $U = {u_{0}, u_{1}, \dots ， u_{m^{'} + k}}$ 和 $V = {v_{0}, v_{1}, \dots ， v_{n^{'} + l}}$ ；

Step3　计算参数 $μ_{i, j}$ 和矩阵权 $M_{i, j}$ ；

Step4　通过点集 ${{P_{i, j}}^{'}}$ 和矩阵权集合 ${M_{i, j}}$ 构造矩阵权NURBS曲面；

Step5　输出矩阵权NURBS曲面 $Q (u, v)$ 。

2.2　矩阵权NURBS曲面的光顺

由前面的讨论可知，对于无噪声的数据，矩阵权NURBS曲面能很好地拟合且光顺性也很好，然而添加噪声后构造的矩阵权NURBS曲面可能并不光顺。为提高构造曲面的光顺性，参考YANG^［5］关于矩阵权NURBS曲线光顺的思路，将其推广应用至曲面情形。首先构造初始的矩阵权NURBS曲面 $S (u, v)$ ，然后在该曲面上重新采样控制顶点，并重新计算每个顶点处的法向量，随后参照矩阵权NURBS曲面拟合算法构造曲面 $Q (u, v)$ ，根据给定的光顺迭代次数重复这一过程，得到更加光顺的曲面并且该曲面仍然逼近原始的控制网格。

假设初始构造的矩阵权NURBS曲面为

S (u, v) = {[\sum_{i = 0}^{m^{'}} \sum_{j = 0}^{n^{'}} M_{i, j} N_{i, k} (u) N_{j, l} (v)]}^{- 1} \times \sum_{i = 0}^{m^{'}} \sum_{j = 0}^{n^{'}} M_{i, j} P_{i, j} N_{i, k} (u) N_{j, l} (v) ，

其中，节点向量 $U = {u_{0}, u_{1}, \dots, u_{m^{'} + k}}$ ， $V = {v_{0}, v_{1}, \dots, v_{n^{'} + l}}$ ，在以下参数处采样： $(u_{i}, v_{j}),$ $i = k - 1, k, \dots, m^{'} + 1, j = l - 1, l, \dots, n^{'} + 1$ ，获得新的控制顶点网格。类似地，对此网格做边界扩展，在边界处添加同样多的边界对称点。对法向量也进行迭代更新，顶点 $P_{i, j}$ 处的法向量由该顶点1-领域内三角面片的法向量面积加权平均得到，将网格中的四边形视为由2个三角形拼接而成。当然原始的法向量数据可能包含噪声，因此在计算得到新的法向量后可考虑采用法向滤波技术^［21］对法向量场进行降噪处理。矩阵权NURBS曲面的光顺算法步骤如下：

输入　点集 ${P_{i, j}}$ 、法向量集 ${n_{i, j}}$ 、光顺迭代次数 $K$ 以及 u， v 两个方向B样条基函数的阶数。

输出　矩阵权NURBS曲面 $Q (u, v)$ 。

Step1　利用输入数据构造初始的矩阵权NURBS曲面 $S (u, v)$ ；

Step2　在曲面 $S (u, v)$ 的节点位置重新采样并将其作为新的控制顶点；

Step3　在曲面的边界添加边界对称点；

Step4　重新计算每个顶点对应的法向量；

Step5　由新的顶点和对应的法向量构造矩阵权NURBS曲面 $Q (u, v)$ ；

Step6　用 $Q (u, v)$ 代替 $S (u, v)$ ，重复step2~step5，直到满足迭代次数 $K$ ；

Step7　输出矩阵权NURBS曲面 $Q (u, v)$ 。

图2为利用带有噪声的控制网格构造的双五次矩阵权NURBS曲面，顶点和法向量采样自一个锥面，且人为地添加噪声。可以看到，初始构造的矩阵权NURBS曲面光顺性较差，经过10次顶点重采样、法向量计算以及曲面重建的迭代，最终得到了光顺的矩阵权NURBS曲面。

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 矩阵权NURBS曲面

Fig.2 The matrix weighted NURBS surface

3　结果与讨论

给出几个具体的用矩阵权NURBS曲面拟合和光顺给定网格数据的例子，并且将实验结果与B样条曲面拟合结果^［1］进行对比。

首先为平面例子，给定一个包含噪声的平面网格数据，通过本文算法构造双三次矩阵权NURBS曲面并分别用光顺算法迭代10次、100次，所得结果如图（3）所示。可以看到，B样条曲面虽然能较好地拟合平面网格，但曲面的光顺性对数据噪声比较敏感；相比之下虽然图3（b）中初始构造的矩阵权NURBS曲面光顺性非常差，但随着迭代次数的增加，曲面逐渐变得平滑，最后得到的双三次矩阵权NURBS曲面不仅逼近初始平面网格，而且有效克服了噪声的影响。其中初始的双三次矩阵权NURBS曲面以及迭代10次和100次的曲面与原网格的最大距离偏差分别为0.063 9%，0.049 2%和0.060 0%，平均距离偏差分别为0.014 2%，0.008 2%和0.005 9%。

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 平面噪声网格的拟合与光顺

Fig.3 Fitting and fairing a noise plane

图4为用双三次B样条曲面和双三次矩阵权NURBS曲面拟合人耳模型得到的结果。与B样条曲面拟合得到的结果相比，光顺迭代5次的矩阵权NURBS曲面的光顺性相对较好，且随着迭代次数的增加，曲面越来越光滑。双三次B样条曲面与原网格的最大距离偏差为0.611 9%，双三次矩阵权NURBS曲面光顺迭代5次和20次后与原网格的最大距离偏差分别为0.474 3%和0.761 2%，平均距离偏差分别为0.054 5%和0.062 6%。

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 矩阵权NURBS曲面拟合与光顺噪声耳朵网格模型

Fig.4 Fitting and fairing an ear mesh model by matrix weighted NURBS surfaces

图5和图6分别为含有噪声的人头模型和牙齿模型通过双三次矩阵权NURBS曲面拟合得到的结果。可以看到，经过多次光顺迭代后构造的矩阵权NURBS曲面令人满意。其中人头模型迭代5次和50次的最大距离偏差分别为0.556 8%和1.227 9%，平均距离偏差分别为0.052 0%和0.045 8%。牙齿模型迭代3次和30次的最大距离偏差分别为0.443 8%和0.616 1%，平均距离偏差分别为0.057 7%和0.070 1%。光顺迭代次数的增多可能会增大构造的矩阵权NURBS曲面与原网格数据的偏差，但也会增加其光顺性。

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 矩阵权NURBS曲面拟合与光顺人头模型

Fig.5 Fitting and fairing a head model by matrix weighted NURBS surfaces

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 矩阵权NURBS曲面拟合与光顺牙齿模型

Fig.6 Fitting and fairing a tooth model by matrix weighted NURBS surface

4　结语

提出了利用矩阵权NURBS曲面拟合和光顺规则的四边形网格数据的方法，对于给定的顶点和法向量，直接选取这些顶点作为矩阵权NURBS曲面的控制顶点，利用法向数据计算矩阵权。当数据从光滑曲面上均匀采样时，构造的矩阵权NURBS曲面能很好地逼近顶点数据；对于包含噪声的数据，通过光顺算法，能直接且高效地构造光顺的矩阵权NURBS曲面。如果网格模型不包含法向数据，也可通过给定的顶点数据估计每个点处的法向量。本文方法适用于规则的四边形网格数据，不能直接用于任意拓扑三角网格模型或点云数据，并且本文只考虑有边界的单片NURBS曲面的拟合，对于可用NURBS表示的封闭网格可类似处理。其次，本文中B样条函数的节点向量是均匀选取的，并没有研究节点向量的不同选取对结果的影响。此外，如果光顺迭代次数过大，模型的尖锐特征也会逐渐被平滑，有待进一步研究。

http://dx.doi.org/10.3785/j.issn.1008-9497.2023.06.017

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

PIEGL

， TILLER

. The Nurbs Book［M］. 2nd ed. New York：Springer-Verlag， 1997. doi:10.1007/978-3-642-59223-2

[本文引用: 3]

[2]

Z J

， MA

， WANG

Z J

， et al.

Surface modeling method for aircraft engine blades by using speckle patterns based on the virtual stereo vision system

［J］. Optics Communications， 2018， 411： 33-39. DOI：10.1016/j.optcom.2017.10.064

[3]

YANG

Z Y

， SHEN

L Y

， YUAN

C M

， et al.

Curve fitting and optimal interpolation for CNC machining under confined error using quadratic B-splines

［J］. Computer-Aided Design， 2015， 66： 62-72. DOI：10.1016/j.cad.2015.04.010

[4]

CATMULL

， CLARK

Recursively generated B-spline surfaces on arbitrary topological meshes

［J］. Computer-Aided Design， 1978， 10（6）： 350-355. DOI：10.1016/0010-4485（78）90110-0