在自由曲面造型领域往往需要拼接多张曲面片，在这个过程中容易产生一些 “洞”，因此曲面 “补洞”问题是计算机辅助几何设计(CAGD)研究的一项重要工作，解决曲面 “补洞” 问题的首要任务是构造插值边界曲线的曲面片^[1]. 目前已有很多方法处理插值3条或4条边界曲线的问题，如Coons等^[2-3]提出的方法. 对于插值5条及5条以上边界曲线的问题，目前有以下5种主要的方法. 1)中心分割法，把参数域从中心分割成三角形域或四边形域，再用相应的面片来插值边界曲线，如Gregory等^[4-6]提出的方法. 2)细分法，利用细分产生的极限曲面来补洞，如Lai等^[7]提出的方法. 3)控制网格法，利用边界条件和连续性条件计算控制网格，生成多边面片，如Várady等^[8-9]提出的方法. 4)裁剪曲面法，利用三角形或四边形面片裁剪成多边面片，如Hamann等^[10]提出的方法. 5)超限插值曲面法，利用边界曲线及其各阶导矢和混合函数构造多边面片，如Várady等^[11-13]提出的方法.

超限曲面插值最早可以追溯到20世纪60年代，Coons^[2]构造出插值型曲面，称为Coons曲面. Gordon^[14]利用三次Hermite基函数构造出双三次Coons面片，将Coons曲面解释为布尔和曲面，提出超限插值的概念. Gordon^[15]将该方法推广到三角形域上，提出三角形域上Coons面片的构造方法. Gregory^[16]提出在正多边形域上构造多边面片的角混合方法，本质上是对多边形的每个角 (相邻的两条边) 构造插值这个角的面片，结合混合函数构造多边面片. Kato^[12]提出边混合方法，对每一条边构造插值该边所对应边界曲线的面片，计算效率比角混合方法高. Várady等^[17]将角混合和边混合方法中规则的正多边形域推广到不规则域. Salvi等^[9]提出三边混合方法，核心是对每一组相连的3条边界曲线构造插值这3条边界曲线及其跨界导矢的面片，结合相应的混合函数得到多边面片，该方法使得多边面片在边界处的曲率过渡更自然. Várady等^[8]提出GB面片，这是结合超限插值曲面法和控制网格法构造的多边Bezier面片，GB面片在边界处可以达到$ {G^1} $连续(切线方向连续). Salvi等^[18]将GB面片的凸多边形参数域推广到一般多边形参数域. Várady等^[19]将一般多边形参数域推广到曲边多边形参数域，这种参数域更接近输入的三维边界曲线围成的形状，可以解决输入复杂边界曲线构造多边面片时失真的问题. Qin等^[20]提出BB面片，这是张量积Bezier曲面在多边形域上的推广，在边界处能够达到$ {G^2} $连续 (曲率连续).

通常相邻面片在边界曲线处至少需要达到切线$ {G^1} $连续，但在某些应用中只需要位置连续 ($ {G^0} $)，或者在边界曲线处法向量和跨界导矢信息缺失的情况下只能达到$ {G^0} $连续^[21]. Salvi^[21]提出简单的三边混合方法，构造边界处$ {G^0} $连续的多边Coons面片. 本文继承了双线性Coons曲面的布尔和构造方法，提出构造简单、计算高效的$ {G^0} $连续多边形域Coons面片.

给定一张四边面片，4条边界曲线分别为 ${\boldsymbol{P}}\left( {u,0} \right)、 {\boldsymbol{P}}\left( {u,1} \right)、{\boldsymbol{P}}(0,v)、{\boldsymbol{P}}(1,v),\; u、v \in \left[ {0,1.0} \right],$ 4个角点分别为$ {\boldsymbol{P}}(0,0)、{\boldsymbol{P}}(0,1)、{\boldsymbol{P}}(1,0)、{\boldsymbol{P}}(1,1)，$ 如图1所示.

Farin^[3]给出双线性Coons面片的定义和构造方式. 给定2条边界曲线$ {\boldsymbol{P}}(u，0) $和$ {\boldsymbol{P}}(u，1) $，可以构造线性插值2条边界曲线的直纹面$ {{\boldsymbol{Q}}_1}(u，v) $，定义为

(1)$ {{\boldsymbol{Q}}_1}(u,v) = {\boldsymbol{P}}(u,0)(1 - v)+{\boldsymbol{P}}(u,1)v.{\text{ }} $

如图2(a)所示，横向直线表示直纹面$ {{\boldsymbol{Q}}_1}(u,v) $的等参线，横向曲线表示未插值的边界曲线，箭头表示参数u变大的方向. 类似地，给定另外2条边界曲线$ {\boldsymbol{P}}(0,v) $和$ {\boldsymbol{P}}(1,v) $，可以构造线性插值2条边界曲线的直纹面$ {{\boldsymbol{Q}}_2}(u,v) $，定义为

(2)$ {{\boldsymbol{Q}}_2}(u,v) = {\boldsymbol{P}}(0,v)(1 - u)+{\boldsymbol{P}}(1,v)u.{\text{ }} $

如图2(b)所示，纵向直线表示直纹面$ {{\boldsymbol{Q}}_2}(u，v) $的等参线，纵向曲线表示未插值的边界曲线，箭头表示参数v变大的方向. 简单迭加上述2张直纹面所得的曲面${\boldsymbol{S}}(u，v) = {{\boldsymbol{Q}}_1}(u，v)+{{\boldsymbol{Q}}_2}(u，v)$，不能同时插值2组边界曲线. 从曲面${\boldsymbol{S}}(u，v)$的每条边界来看，迭加导致多出了连接边界两端点的直边. 为了得到同时插值4条边界曲线的曲面，迭加所得的曲面必须减去插值4个角点的双线性曲面${{\boldsymbol{Q}}_3}(u，v)$，如图2(c)所示，定义为

(3)$ \begin{split} {{\boldsymbol{Q}}_3}(u,v) = &[{\boldsymbol{P}}(0,0)(1 - u)+{\boldsymbol{P}}(1,0)u](1 - v)+ \\ &{\text{ [}}{\boldsymbol{P}}(0,1)(1 - u)+{\boldsymbol{P}}(1,1)u{\text{]}}v{\text{.}} \end{split}$

插值给定4条边界曲线的双线性Coons面片${\boldsymbol{Q}}(u,v)$等参线图，如图2(d)所示. 双线性Coons曲面${\boldsymbol{Q}}(u,v)$的参数表达式可以简写为如下的布尔和形式：

(4)$ {\boldsymbol{Q}}(u,v) = {{\boldsymbol{Q}}_1}(u,v)+{{\boldsymbol{Q}}_2}(u,v) - {{\boldsymbol{Q}}_3}(u,v).{\text{ }} $

Sederberg等^[22-23]研究双线性映射及其逆映射. 双线性坐标可以定义为：设四边形的4个顶点分别是$ {{\boldsymbol{P}}_1}、{{\boldsymbol{P}}_2}、{{\boldsymbol{P}}_3}、{{\boldsymbol{P}}_4} $，则四边形内的点${\boldsymbol{P}}$关于4个顶点的双线性坐标表示为

$ {\boldsymbol{P}} = {\lambda _1}{{\boldsymbol{P}}_1}+{\lambda _2}{{\boldsymbol{P}}_2}+{\lambda _3}{{\boldsymbol{P}}_3}+{\lambda _4}{{\boldsymbol{P}}_4}. $

式中：

(5)$ \left.\begin{aligned} {\lambda _1} =& (1 - u)(1 - v),{\text{ }}{\lambda _2} = (1 - u)v, \\ {\lambda _3} =& uv,{\text{ }}{\lambda _4} = u(1 - v),{\text{ }}u、v \in [0,1.0].\end{aligned} \right\}$

由式(5)可以得出参数$u、v、1 - u、1 - v$和双线性坐标的关系:

(6)$ \left. \begin{aligned} {\lambda _1}+{\lambda _2} =& 1 - u,{\text{ }}{\lambda _2}+{\lambda _3} = v, \\ {\lambda _3}+{\lambda _4} =& u,{\text{ }}{\lambda _4}+{\lambda _1} = 1 - v.\end{aligned}\right\} $

设$ P \subset {{\bf{R}}^2} $是凸多边形，顶点为$ {{\boldsymbol{P}}_1}，{{\boldsymbol{P}}_2}， \cdots ，{{\boldsymbol{P}}_n}， n \geqslant 3 $， 以逆时针方向排列， 如图3(a)所示为$n = 5$的例子.

函数${\lambda _i}({\boldsymbol{x}}) \geqslant 0,{\text{ }}i = 1, \cdots ,n,{\text{ }}{\boldsymbol{x}} \in P$满足以下条件：

$ \sum\limits_{i = 1}^n {{\lambda _i}({\boldsymbol{x}}) = 1} ,{\text{ }}\sum\limits_{i = 1}^n {{\lambda _i}({\boldsymbol{x}}){{\boldsymbol{P}}_i} = {\boldsymbol{x}}.} {\text{ }} $

称函数${\lambda _i}:P \to {\bf{R}},{\text{ }}i = 1，2，\cdots ,n$为（广义）重心坐标. 当$n = 3$时，函数${\lambda _1}、{\lambda _2}、{\lambda _3}$是唯一确定的，且函数值是x点关于其3个顶点$ {{\boldsymbol{P}}_1}、{{\boldsymbol{P}}_2}、{{\boldsymbol{P}}_3} $的常规的重心坐标(面积坐标).当$ n \geqslant 4 $时，${\lambda _1}，{\lambda _2}，\cdots ,{\lambda _n}$的选择不再唯一.

Wachspress坐标是广义重心坐标中的一种，由Wachspress^[24]提出，Warren^[25]改进. Wachspress坐标定义为

(7)$ {\lambda _i}({\boldsymbol{x}}) = \frac{{{w_i}({\boldsymbol{x}})}}{{\displaystyle \sum\nolimits_{j = 1}^n {{w_j}({\boldsymbol{x}})} }}.{\text{ }} $

式中：

$ {w_i}({\boldsymbol{x}}) = \frac{{A({{\boldsymbol{P}}_{i - 1}},{{\boldsymbol{P}}_i},{{\boldsymbol{P}}_{i+1}})}}{{A({\boldsymbol{x}},{{\boldsymbol{P}}_{i - 1}},{{\boldsymbol{P}}_i})A({\boldsymbol{x}},{{\boldsymbol{P}}_i},{{\boldsymbol{P}}_{i+1}})}}{\text{. }} $

其中$ A({{\boldsymbol{x}}_1},{{\boldsymbol{x}}_2},{{\boldsymbol{x}}_3}) $表示以$ {{\boldsymbol{x}}_1}、{{\boldsymbol{x}}_2}、{{\boldsymbol{x}}_3} $为顶点的三角形的有向面积，如图3(b)所示.

利用式(7)，可以求出任意一个凸多边形的Wachspress坐标，图4给出正五边形的Wachspress坐标等高线图.

利用式 (6)，可以将式 (1)~(3) 改写成双线性坐标形式，如下所示：

(8)$ \begin{split} {{\boldsymbol{Q}}_1}(u,v) = &{\boldsymbol{P}}(u,0)(1 - v)+{\boldsymbol{P}}(u,1)v = \\ &{\boldsymbol{P}}({\lambda _3}+{\lambda _4},0)({\lambda _4}+{\lambda _1})+{\boldsymbol{P}}({\lambda _3}+{\lambda _4},1)({\lambda _2}+{\lambda _3}),{\text{ }}\end{split}$

(9)$ \begin{split} {{\boldsymbol{Q}}_2}(u,v) =& {\boldsymbol{P}}(0,v)(1 - u)+{\boldsymbol{P}}(1,v)u = \\ &{\boldsymbol{P}}(0,{\lambda _2}+{\lambda _3})({\lambda _1}+{\lambda _2})+{\boldsymbol{P}}(1,{\lambda _2}+{\lambda _3})({\lambda _3}+{\lambda _4}),\end{split} $

(10)$ \begin{split} {{\boldsymbol{Q}}_3}(u,v) = &[{\boldsymbol{P}}(0,0)(1 - u)+{\boldsymbol{P}}(1,0)u](1 - v)+ \\ &[{\boldsymbol{P}}(0,1)(1 - u)+{\boldsymbol{P}}(1,1)u]v = \\ &{\lambda _1}{{\boldsymbol{P}}_1}+{\lambda _2}{{\boldsymbol{P}}_2}+{\lambda _3}{{\boldsymbol{P}}_3}+{\lambda _4}{{\boldsymbol{P}}_4}.\end{split} $

为了使双线性Coons面片的表达式更加简洁，将式 (8)~(10) 中的边界曲线$ {\boldsymbol{P}}(0，v)、{\boldsymbol{P}}(u，1)、 {\boldsymbol{P}}(1，v) $和$ {\boldsymbol{P}}(u，0) $分别记作$ {{\boldsymbol{C}}_1}、{{\boldsymbol{C}}_2}、{{\boldsymbol{C}}_3}、{{\boldsymbol{C}}_4} $，规定逆时针方向为正方向，如图5所示.

如图2(d)所示，双线性Coons面片$ {\boldsymbol{Q}}(u,v) $中参数$ u、v $的正方向分别是纵向、横向箭头所指的方向.边界曲线$ {\boldsymbol{P}}(1,v) $和$ {\boldsymbol{P}}(u,0) $的正方向分别与$ {{\boldsymbol{C}}_3} $、$ {{\boldsymbol{C}}_4} $的方向相反，通过修改边界曲线的参数，使得边界曲线参数正方向与$ {{\boldsymbol{C}}_3} $、$ {{\boldsymbol{C}}_4} $的方向一致，即

$ {{\boldsymbol{C}}_3} = {\boldsymbol{P}}(1,1 - v) = {\boldsymbol{P}}(1,{\lambda _4}+{\lambda _1}),{\text{ }} $

$ {{\boldsymbol{C}}_4} = {\boldsymbol{P}}(1 - u,0) = {\boldsymbol{P}}({\lambda _1}+{\lambda _2},0).{\text{ }} $

$ {{\boldsymbol{Q}}_1}(u,v) $和$ {{\boldsymbol{Q}}_2}(u,v) $的表达式可以改写如下：

(11)$ {{\boldsymbol{Q}}_1}(u,v) = ({\lambda _4}+{\lambda _1}){{\boldsymbol{C}}_4}({\lambda _1}+{\lambda _2})+ {\text{ }}({\lambda _2}+{\lambda _3}){{\boldsymbol{C}}_2}({\lambda _3}+{\lambda _4}), $

(12)$ {{\boldsymbol{Q}}_2}(u,v) = ({\lambda _1}+{\lambda _2}){{\boldsymbol{C}}_1}({\lambda _2}+{\lambda _3})+ {\text{ }}({\lambda _3}+{\lambda _4}){{\boldsymbol{C}}_3}({\lambda _4}+{\lambda _1}). $

将式 (10)~(12) 代入式 (4)，整理合并可得

(13)$ {\boldsymbol{Q}}({\lambda _1},{\lambda _2},{\lambda _3},{\lambda _4}) = \sum\limits_{i = 1}^4 {({\lambda _i}+{\lambda _{i+1}}){{\boldsymbol{C}}_i}({\lambda _{i+1}}+{\lambda _{i+2}})} - {\text{ }}\sum\limits_{i = 1}^4 {{\lambda _i}{{\boldsymbol{P}}_i}} {\text{. }}$

式中：$ {\lambda _i} $和$ {{\boldsymbol{P}}_i} $为循环索引，即$ {\lambda _5} = {\lambda _1}，{{\boldsymbol{P}}_5} = {{\boldsymbol{P}}_1}$，$ {{\boldsymbol{C}}_i} $为边界曲线的函数，对应关系如图5所示.

类似地，将四边形域上的双线性Coons面片推广到多边形域上，可以得到多边形域Coons面片的定义.

定义2.1　给定n条 (n > 3) 定义在[0,1.0]上的边界曲线$ {{\boldsymbol{C}}_1},{{\boldsymbol{C}}_2},\cdots ,{{\boldsymbol{C}}_n} $，它们首尾相连构成封闭曲线，其n个角点为$ {{\boldsymbol{P}}_1},{{\boldsymbol{P}}_2},\cdots ,{{\boldsymbol{P}}_n} $，参数域是凸多边形，则多边形域Coons面片的定义为

(14)$ {\boldsymbol{Q}}({\lambda _1}, \cdots ,{\lambda _n}) = \sum\limits_{i = 1}^n {({\lambda _i}+{\lambda _{i+1}}){{\boldsymbol{C}}_i}({\lambda _{i+1}}+{\lambda _{i+2}})} {\text{ }} - \sum\limits_{i = 1}^n {{\lambda _i}{{\boldsymbol{P}}_i}} {\text{. }} $

式中：$ \left\{ {{\lambda _i}} \right\}(i = 1, 2, \cdots ,n )$为多边形中的点关于多边形顶点${{\boldsymbol{P}}_i}$的Wachspress坐标，$ {{\boldsymbol{C}}_i} $为相应的边界曲线函数，对应关系如图6所示.

多边形参数域可以由n个角点投影到二维平面后顺次连接得到，但该方法在某些情况下可能会生成凹多边形参数域. Várady等^[17]针对这种可能出现的情况，提出基于边长度/角生成凸多边形参数域的方法，利用该方法解决了参数域可能出现非凸多边形的问题.

定理2.1　多边形域Coons面片具有边界插值性.

由上述多边形域Coons面片的定义可知，任一边界曲线$ {{\boldsymbol{C}}_k} $的两端点分别为$ {{\boldsymbol{P}}_k}、{{\boldsymbol{P}}_{k+1}} $，即$ {{\boldsymbol{C}}_k}(0) = {{\boldsymbol{P}}_k}，{{\boldsymbol{C}}_k}(1) = {{\boldsymbol{P}}_{k+1}} $. 由Wachspress坐标的定义及性质可知，在边界曲线$ {{\boldsymbol{C}}_k} $上的点满足:

(15)$ {\lambda _k}+{\lambda _{k+1}} = 1;\;{\lambda _i} = 0,i \ne k,k+1.{\text{ }} $

将式 (15) 代入式 (14)，可得多边形域Coons面片在边界处的情况，代入化简结果如下:

$ \begin{split} &{\boldsymbol{Q}}(0, \cdots ,0,{\lambda _k},{\lambda _{k+1}},0, \cdots ,0) = {\lambda _k}{{\boldsymbol{C}}_{k - 1}}(1)+ \\ &{\text{ }}{{\boldsymbol{C}}_k}({\lambda _{k+1}})+{\lambda _{k+1}}{{\boldsymbol{C}}_{k+1}}(0) - {\lambda _k}{{\boldsymbol{P}}_k} - {\lambda _{k+1}}{{\boldsymbol{P}}_{k+1}} = {\text{ }} \\ &{\text{ }}{{\boldsymbol{C}}_k}({\lambda _{k+1}}). \end{split}$

综上所述，多边形域Coons面片插值任一边界曲线$ {{\boldsymbol{C}}_k} $，即在边界处$ {G^0} $连续.

仅以Wachspress坐标为例，介绍多边形域Coons面片的构造方法，实际上任何能应用于凸多边形域上的广义重心坐标都能应用于多边形域Coons面片，如均值坐标^[26]、调和坐标^[27]、迭代坐标^[28]等.

下面以一些四边网格模型为例，探讨本文所述方法和前人方法在 “补洞” 效果以及计算效率方面的差异. 在四边网格中，度为4的点称为非奇异点，其他点称为奇异点. 以四边网格中$6 \times 6$的非奇异点为控制网格，可以得到双五次Bezier曲面，对所有非奇异点构造双五次Bezier曲面，则奇异点附近会产生 “洞”，且这些 “洞” 的边界曲线都是5次Bezier曲线. 由于提出的多边形域Coons面片不适用于三边的情况，对于三边的 “洞”，采用Gordon-Coons三角面片^[15]填充；对于四边以上的 “洞”，采用本文方法和Salvi^[21]提出的方法填补. 将提出的多边形域Coons面片和Salvi^[21]给出的边界处$ {G^0} $连续的Coons面片分别记为MGC面片和MC面片，利用上述2种方法补全人体、机器人、小鸡和小牛模型中的 “洞”，如图7、8所示. 利用2种方法生成的多边面片与周围面片在边界处都是$ {G^0} $连续，所以二者补洞的效果差异难以用肉眼分辨，但通过图9可以看出，利用2种方法生成的多边面片在曲率分布方面略有不同. 因为2种面片的构造方法不同，计算效率有较大差异. MC面片是利用三边混合方法构建多边面片，对每组相连的3条边界曲线都构造一张矩形域上的双线性Coons面片，其中3条待插值边界曲线是给定的，另一条边界曲线是根据边界信息计算得到的. n边MC面片涉及4n条边和4n个角点的计算; 本文方法只涉及n条边和n个角点的计算，且MC面片需要预先计算有理局部参数，增加了计算量. 理论上MGC面片的计算效率高于MC面片，2种方法的计算时间t_c对比如表1所示，表中数据均在CPU基准速度为3.2 GHz的机器上运行得到，5边和6边面片的采样点数分别为18 602和20 439.

如图9(a)、(b)所示为图7(a)、8(a)中MGC面片和MC面片 “补洞” 的人体模型高斯曲率图. 如图9(c)、(d)所示分别为图9(a)、(b)中方框区域的局部放大. 可以看出，利用2种方法所生成的多边面片曲率分布略有差异. 如表1所示为2种方法在4个模型下的平均计算时间对比. 可以看出，MGC面片在4种模型下的5边和6边面片上的平均计算时间分别为5.04和6.59 s，MC面片的平均计算时间分别为41.27和54.34 s，二者的平均计算时间之比约为3∶25，表明MGC面片的计算效率远高于MC面片.

本文提出基于广义重心坐标的多边形域Coons面片构造方法，该方法具有构造方法简单和计算高效的优点. 该方法的参数域限制于凸多边形，且只插值了边界曲线而没有插值一阶跨界导矢、二阶跨界导矢、角点扭矢等更多的信息，曲面在边界处只能达到${G^0}$连续.如何将参数域推广到一般多边形且插值一阶跨界导矢、二阶跨界导矢、角点扭矢等信息，使它在边界处能够达到${G^1}、{G^2}$连续是进一步的研究方向.