第3章几 何 建 模计算机图形学的早期发展动力源自于CAD/CAE/CAM系统中对几何造型的需求。它们需要通过计算机表示、存储、分析、控制并输出指定形状的几何模型。这里的形状既可以是单个物体的外形，也可以是由一组物体所构成的复杂场景的外形。此外，在各种类型图形的真实感绘制过程中，模型的几何形状也会直接影响其表面光照分布、纹理密度等因素，进而与模型的绘制效果也是密切相关的。因此，长期以来几何建模都是计算机图形学中的重要内容。
本章介绍几何建模涉及的数学基础知识，包括形状表达的数学形式、常用几何性质等。接下来围绕三种重要的几何建模技术: 自由曲线/曲面建模、细分曲面建模和三维重建，介绍相应的概念、模型及其使用方法。此外，还介绍若干其他类型的建模方法，例如基于分形、粒子系统、语言学等的方法。最后，针对计算机内部数据存储特点，介绍面向几何建模的典型数据结构。
3.1数 学 基 础
一般来讲，计算机图形学中的几何建模是指构造物体或场景二维/三维形状的数学表达形式及其在计算机内表示的数据结构。通常意义下，这种数学形式具有直接、明确的函数关系，使得所有满足该关系的形状也具有相应的几何性质。因此，首先需要了解几何形状的数学形式。
3.1.1几何形状的数学形式
数学上，连续几何形状的表达方式主要有三种: 显式表达式、隐式表达式和参数表达式。其中，显式表达式和隐式表达式又称为代数形式。这三种形式的主要区别在于定义函数时，其因变量和自变量的表达方式和相互之间代数制约关系的不同。下面分别介绍这三种表达式。
1. 显式表达式
显式表达式是一种因变量随自变量变化而变化的函数形式。二维平面上几何形状的显式表达式通常为y=f(x)，其中x是自变量，y是因变量。例如，y=2x+1代表平面直线，y=x2+x+1代表平面抛物线。三维空间中几何形状的显式表达式通常为z=f(x,y)，其中x和y是自变量，z是因变量。例如，z=2x+y+1表示了三维空间中的一张平面，z=1-x2-y2表示了半径为1的半球面。从表达形式上看，显式表达式的因变量和自变量分别位于等号两侧，二者具有明显的对应关系。
2. 隐式表达式
隐式表达式是由多个变量共同定义的一种函数形式。这里的变量是没有自变量和因变量之分的，不能单独进行表示。二维平面上几何形状的隐式表达式通常具有f(x,y)=0的形式。例如，x2+y2-1=0表示了平面圆形。三维空间中几何形状的隐式表达式具有f(x,y,z)=0的形式。例如，x2+y2+z2-1=0表示了三维球面。从表达形式上看，隐式表达式的因变量和自变量是混合在一起的，都位于等式的同侧，二者没有明显的对应关系。
〖1〗现代计算机图形学基础第3章几何建模〖3〗〖3〗3. 参数表达式
参数表达式是采用若干独立变量作为自变量的显式表达式组成的集合。这种表达方式下几何形状是由指定集合的参数作为自变量而得到的因变量所表示的。作为自变量的参数也称为参变量。例如，二维平面和三维空间中的几何形状分别具有如下表达形式: x=f(t)
y=g(t)x=f(u,v)
y=g(u,v)
z=h(u,v)(3.1)其中，t、u和v是参变量。它们在给定取值范围内变化，从而产生相应的几何形状。
从代数形式来讲，参数表达式也可以看作为一种显式表达式。但是，参数表达式具有更加直观的几何意义，反映了参变量变化所产生的各个维度的变化情况。例如，曲线可以看作单个参变量在指定区间内变化所形成的空间点的集合；而曲面则可以看作双参变量变化形成的空间点的集合。此外，参数表达式的代数和微积分运算非常方便，适合对几何形状的性质进行分析。
对于一些几何形状，可能存在多种不同形式的表达式。如图3.1所示，三维空间中半径为1的球面，除了上述隐式表达式，也可以采用基于球极坐标的参数表达式，也就是x=sinφcosθ
y=sinφsinθ
z=cosφ(3.2)其中，参变量θ和φ分别对应于方位角和俯仰角。但是对于球面而言，一般不存在直接的显式表达式。
图3.1球面的隐式表达式与参数表达式
例题31平面曲线的数学表达式。
问题: 如右图所示，平面坐标系第一象限内的单位圆弧的显式表达式和参数表达式分别是什么？
解答: 该圆弧的显式表达式可以写为: y=1-x2,0≤x≤1对应的参数表达式可以写为: x(t)=(1-t2)/(1+t2)
y(t)=(2t)/(1+t2),0≤t≤1□
3.1.2几何性质
通过几何形状的数学形式，可以研究模型所具备的几何性质，从而更好地设计和控制建模效果。计算机图形学中几何建模所涉及的形状主要包括两种类型: 曲线和曲面。曲线又包括二维平面曲线和三维空间曲线，而曲面主要是指三维空间中的曲面。这些也是在日常生产生活中，人们接触最多的物体形状。
1. 二维平面曲线
二维平面曲线理论上可以采用显式、隐式或参数形式进行表达。显式形式是具有一个自变量和一个因变量的表达式，也就是y=f(x)。隐式形式是由两个变量构成的方程式，也就是f(x,y)=0。而参数形式则是单参变量表示的二维向量，也就是γ(t)=(x(t),y(t))。例如，平面上椭圆的隐式表达式和参数表达式分别具有如下形式: x2a2+y2b2=1,x=acosθ
y=bsinθ(3.3)其中，a和b是非零常数，θ是参变量。平面曲线常用的几何性质有曲线弧长、曲率等，反映了曲线本身的形状。
以参数形式表达的曲线为例，它的弧长描述了曲线在某一参变量取值区间内的测度，计算公式为l=∫t1t0x·2(t)+y·2(t)dt其中，x·(t)和y·(t)是函数关于参变量t的一阶导数，［t0,t1］是t的取值区间。曲率定义为切线方向角相对于弧长的变化率，计算公式为κ(t)=|x·(t)y··(t)-x··(t)y·(t)|(x·2(t)+y·2(t))32其中，x··(t)和y··(t)是函数关于参变量t的二阶导数。实际上，曲率反映了曲线在某一点的弯曲程度，例如直线的曲率处处为0。此外，平面曲线上某一点的曲率大小等于该点处密切圆半径r的倒数，而密切圆的圆心位于曲线凹向的一侧（如图3.2(a)所示）。 
图3.2平面二维曲线和三维空间曲线的曲率圆
2. 三维空间曲线
三维空间曲线可以采用单个参变量作为参数的三维向量进行表示，也就是γ(t)=(x(t),y(t),z(t))。它直观上反映了质点在三维空间中单自由度的运动轨迹。这里需要注意的是，三维空间曲线一般不存在显式表达式，而其隐式表达式由两张曲面的隐式表达式组成的方程组所构成，可以记为f(x,y,z)=0
g(x,y,z)=0它反映了三维空间曲面可以看作两张曲面相交而产生。弧长、曲率和挠率是空间曲线的三种基本属性，反映了三维曲线的形状。与平面曲线类似，弧长表示了指定区间内的曲线长度，计算公式为l=∫t1t0x·2(t)+y·2(t)+z·2(t)dt曲率反映了曲线在空间中的弯曲程度，计算公式为κ(t)=|γ·(t)×γ··(t)|(|γ·(t)|)32其中，γ·(t)和γ··(t)分别是曲线所对应三维向量的一阶和二阶导数。如图3.2(b)所示，曲率也与该点密切圆半径r成反比关系。与平面曲线不同，空间曲线具有挠率属性，计算公式为τ(t)=(γ·(t)×γ··(t))·γ···(t)|γ·(t)×γ··(t)|2挠率反映了曲线切平面的扭转状况。例如，所有平面曲线的挠率都为0。
3. 三维空间曲面
图3.3环面三维空间曲面理论上可以采用显式、隐式或参数形式进行表达。显式形式是具有两个自变量和一个因变量的表达式，通常记为z=f(x,y)。隐式形式是由三个变量构成的方程式，也就是f(x,y,z)=0。参数形式则是由两个参变量表示的三维向量形式，通常记为π(u,v)=(x(u,v),y(u,v),z(u,v))。例如，图3.3所示环面的隐式表达式和参数表达式分别为: (x2+y2+z2+R2-r2)2-4R2(x2+y2)=0,x(u,v)=(R+rcosv)cosu
y(u,v)=(R+rcosv)sinu
z(u,v)=rsinv(3.4)其中R和r分别表示环截面半径和内环半径。
空间曲面常用的几何性质涉及面积、法曲率、主曲率、高斯曲率、平均曲率等属性。以参数形式表示的曲面为例，其面积是指双参变量在取值范围内的表面测度，计算公式为s=∫u1u0∫v1v0|π·u×π·v|dudv其中［u0,u1］和［v0,v1］分别是参变量u和v的取值范围。曲面上某点处的曲率和经过该点的曲面上曲线的弯曲程度相关，称为法曲率。事实上，曲面在一点处有无穷多个切方向，因此可以定义无穷多个法曲率。其中，法曲率的最大值和最小值称为主曲率，代表了该点处法曲率的极值分布情况，也就是能达到的最大和最小弯曲程度。曲面上某点处主曲率的平均值，称为平均曲率，而主曲率的乘积则称为高斯曲率。事实上，平均曲率和高斯曲率是反映曲面局部形状的两个重要属性。
3.1.3建模工具
根据建模对象的不同，几何建模可以简单地分为自然物体建模和人造物体建模。所谓自然物体建模，是指建模对象来自于自然界中正常存在的物体，例如各种动物、植物、地形、地貌等（如图3.4(a)所示）。这类物体的几何形状往往很难用精确的数学公式直接表示。人造物体建模，是指通过手工交互或自动化的方式，利用相应的工具对基本几何形状进行变换和重组，获得具有新形状的物体（如图3.4(b)所示）。这类物体的几何形状，通常可以借助变换或重组时的数学表达式进行准确描述。
图3.4物体对象几何建模的方式: 插值与拟合
无论自然物体建模，还是人造物体建模，主要有插值和拟合两种方式（如图3.4(c)所示）。这两种方式都是先对物体形状有个大概的主观描述，在此基础上按照精度的要求进行模型的构建。插值，是要求模型能够严格满足给定数据点的几何位置约束。这些数据点也称为型值点，是对建模对象外形的离散采样。拟合则是在一定的几何意义下，建模的形状尽可能逼近给定的型值点。不论对模型进行插值或拟合，常用的数学工具有自由曲线/曲面、细分曲面等，它们为自然物体建模和人造物体建模提供了有效手段。
3.2自由曲线/曲面建模
经典几何形状，如平面、圆柱面、圆锥面等都具有明确的解析表达式。然而，绝大部分的自然物体，例如人脸等，其形状往往无法使用一个解析函数进行显式表示，需要构造新的函数形式来表示其形状。所谓自由曲线/曲面，是指不能用初等解析函数完全准确地表达全部形状的曲线和曲面。因此，自由曲线/曲面的出现，是为了克服使用解析函数进行形状表示的局限性，采用更广泛意义下的数学函数来表示形状。
3.2.1平面三次多项式曲线
多项式曲线是一种最常见的平面曲线表达形式。在给定型值点集合后，计算n次多项式进行插值或拟合。这样就可以转化为线性方程组来求解该多项式。数学上，多项式次数越高，其对应曲线出现的拐点就越多，能表示的形状也就越复杂。但是次数越高，越可能在拟合时出现过拟合，难以控制那些没有采样点处的形状。此外，五次以上多项式的计算比较复杂，不利于快速建模。如果多项式次数较低，又会出现欠拟合情况，导致形状的表现力不足。在实际应用中，三次多项式曲线是使用较广泛的一类多项式曲线。它对应的参数表达式为: p(t)=c0+c1t+c2t2+c3t3(3.5)图3.5由4个型值点定义的
三次多项式曲线
其中,ci=(xi,yi)T称为控制顶点，t∈［0,1］是参变量取值范围。通过公式(3.5)可以看出，{ci}取值不同，就会产生不同的曲线形状。
如图3.5所示，给定4个型值点{p1,p2,p3,p4}，可以计算通过这些型值点的一条三次多项式曲线。这是由于三次多项式总共有4个系数，每个系数同时有x和y两个分量，因此一共有8个未知数。而4个型值点恰好可以提供8个约束条件。如果给定的型值点多于4个，那么可以通过最小二乘法计算一条最优的三次多项式曲线来拟合这些型值点，使得它们距离曲线的总体几何距离最小。
例题32平面三次多项式插值曲线。
问题: 计算图3.5所示的多项式曲线的参数表达式p(t)。
解答: 根据4个型值点{p1,p2,p3,p4}的x和y坐标，分别计算p(t)=(x(t),y(t))控制顶点ci=(xi,yi)T。其中，x和y坐标分别满足如下等式: 0=x0
2=x0+14x1+116x2+164x3
0=x0+34x1+916x2+2764x3
2=x0+x1+x2+x3
0=y0
2=y0+14y1+116y2+164y3
3=y0+34y1+916y2+2764y3
4=y0+y1+y2+y3通过求解上述两个线性方程组，可以得到4个控制顶点的坐标分别是c0=(0,0)、c1=(18,12)、c2=(-48,-18.67)和c3=(32,10.67)。
□
平面三次多项式能够通过插值或拟合的方式对平面曲线建模，但是像公式(3.5)这样的表达式缺少直观的几何解释（例如，其控制顶点无法直接描述曲线形状），这样就不利于工程人员按照主观意图进行曲线建模来设计相应的形状。此外，受单个多项式代数性质的影响，能有效表示的曲线形状有限，并不适合更加丰富多样的几何外形设计。因此，需要更灵活的自由曲线/曲面建模技术，典型的建模技术有Bézier曲线/曲面和B样条曲线/曲面。
3.2.2Bézier曲线/曲面
Bézier曲线是以法国工程师Pierre Bézier命名的。而Bézier曲线的计算方法最早可以追溯到法国雷诺公司工程师Paul de Casteljau。他提出了适用于计算机编程的递归算法，只是没有意识到这种算法最后所生成的曲线形状就是Bézier曲线的形状。1962年，Bézier明确给出了这种曲线的数学公式，并成功地应用于汽车外形设计。
Bézier曲线是第一种在工业界被广泛采用的自由曲线建模工具，其本质也是多项式曲线。但是，这种曲线具有更直观的几何表达形式，能够方便工程师通过操纵控制顶点来进行建模。具体地讲，n次Bézier曲线是一种单参变量的参数表达式，数学上定义为: p(t)=∑ni=0Bni(t)ci=∑ni=0n
i(1-t)n-itici,t∈［0,1］(3.6)其中，{c0,c1,…,cn}是n+1个控制顶点。这些控制顶点依次连接形成控制多边形。单参变量函数Bni(t)也称为Bernstein多项式函数。图3.6分别展示了二次和三次Bézier曲线。
图3.6二次和三次Bézier曲线
数学上，多项式函数空间是可以采用Bernstein函数作为一组基函数。因此，Bézier曲线就是采用这组基函数的线性组合来表示的一类多项式函数。在此基础上，可以将公式(3.5)所表示的任意形式的多项式曲线转化为公式(3.6)表示的Bézier曲线。这个过程实质上就是将多项式基函数{1,t,t2,t3}转化为Bernstein基函数表示。
Bézier曲线之所以在工业设计中广受欢迎，主要是因为公式(3.6)定义的几何形状具有以下良好性质，能够方便设计人员对形状进行控制。
（1） 首末端点插值。曲线经过控制多边形的首末端点，也就是曲线端点满足p(0)=c0和p(1)=cn。此外，曲线在两个端点处的切线方向与控制多边形的边平行，也就是曲线插值首末端点的切向量。
（2） 保凸性。曲线位于控制多边形构成的凸包（凸多边形边界）内。这样在给定控制多边形后，就可以通过凸包来限定生成的Bézier曲线的范围。这个性质方便设计人员通过控制顶点对形状进行调控。
（3） de Casteljau递归算法。在计算公式(3.6)表示的Bézier曲线时，除了直接采用多项式的代数运算，还可以借助de Casteljau递归算法实现更加快速地计算。该递归算法基于以下递推公式:p［r］i(t)=(1-t)·p［r-1］i-1(t)+t·p［r-1］i(t)
p［0］i(t)≡p［0］i=ci(3.7)上述de Casteljau递归算法证明了当r=n时，p［n］n(t)一定是位于Bézier曲线上的点。这种递归算法非常适合计算机编程实现，计算效率很高。
例题33平面三次Bézier曲线。
问题: 写出右图所示的控制顶点定义的三次Bézier曲线，其中4个控制顶点坐标分别是c0(-2,0)，c1(-1,1)，c2(1,1.5)，c3(2,0)。然后，给出采用de Casteljau递归算法计算曲线上任意一点的递归过程。
解答: 将4个控制顶点的坐标代入公式(3.6)，可以得到如下相应的Bézier曲线表达式: p(t)=∑3i=0B3i(t)ci
=-2B30(t)-B31(t)+B32(t)+2B33(t)
B31(t)+1.5B32(t)
=-2(1-t)3-3(1-t)2t+3(1-t)t2+2t3
3(1-t)2t+4.5(1-t)t2其中，t∈［0,1］。进一步，根据公式(3.7)定义的de Casteljau递归算法，对于曲线上任意一点p(t0)，0≤t0≤1，其递归计算过程如下: p(t0)=p［3］3(t0)=(1-t0)p［2］2(t0)+t0p［2］3(t0)
p［2］2(t0)=(1-t0)p［1］1(t0)+t0p［1］2(t0)
p［2］3(t0)=(1-t0)p［1］2(t0)+t0p［1］3(t0)
p［1］1(t0)=(1-t0)p［0］0(t0)+t0p［0］1(t0)
p［1］2(t0)=(1-t0)p［0］1(t0)+t0p［0］2(t0)
p［1］3(t0)=(1-t0)p［0］2(t0)+t0p［0］3(t0)其中，p［0］0(t0)=c0，p［0］1(t0)=c1，p［0］2(t0)=c2以及p［0］3(t0)=c3。
□
进一步，Bézier曲面是由两个参变量的Bernstein混合函数表示的参数曲面。具体来讲，m×n次Bézier曲面是由(m+1)×(n+1)个控制顶点组成的多面体（也称为控制网格）来定义，记为: p(s,t)=∑mi=0∑nj=0cijBmi(s)Bnj(t),s∈［0,1］, t∈［0,1］(3.8)与Bézier曲线类似，Bézier曲面也具有很多良好的性质。例如，Bézier曲面的边界是由4个边界多边形作为控制多边形所定义的4条Bézier曲线；角点插值性，也就是说Bézier曲面在4个角点插值控制网格的顶点；凸包性，也就是说Bézier曲面位于曲面控制网格形成的凸包内。
与一般的平面多项式曲线/曲面表达式相比，Bézier曲线/曲面的主要优点是容易编程实现、具有端点和切向插值等性质，而且方便各种微积分运算的数值求解。但Bézier曲线/曲面的缺点也很明显，就是缺乏对曲线形状的局部可控性。具体来讲，改变其中一个控制顶点的位置，就会改变整个曲线/曲面的形状。这种局限性也限制了Bézier曲线/曲面的几何建模能力。
3.2.3B样条曲线/曲面
为了克服Bézier曲线/曲面的不足，美国Utah大学的William Gordon和Richard Riesenfeld在1974年将B样条参数曲线引入几何建模。样条，源于生产实践，本意是指富有弹性的细长条。样条利用压铁使其通过指定的型值点，并调整样条使它具有满意的形状，然后沿样条画出曲线。B样条曲线则是通过B样条函数表示的曲线，本质上是分段多项式曲线。因此，B样条曲线是由多个在连接处满足一定连续性的多项式曲线所构成的曲线，是一类多项式组合的曲线。在此之前，美国WisconsinMadison大学的Isaac Schoenberg早在1946年就利用B样条进行统计数据的光滑处理，开创了样条逼近的现代理论。
由于分段多项式的性质，B样条曲线的定义除了需要控制顶点{c0,c1,…,cn}，还需要设置节点向量t0≤t1≤…≤tn+k+1。在此基础上，k次（k+1阶）B样条曲线定义为: p(t)=∑ni=0Nki(t)ci,n≥k(3.9)这里，相邻两个节点形成的区间［tk,tk+1)定义了在其每一段上的多项式函数，而这些多项式函数的组合就定义了B样条曲线。图3.7展示了一条三次B样条曲线。
图3.7三次B样条曲线
具体来讲，在公式(3.9)中，Nki(t)是B样条基函数，可以通过如下递推公式来定义: N0i(t)=1,t∈［ti,ti+1)
0,t［ti,ti+1)
Nki(t)=t-titi+k-tiNk-1i(t)+ti+k+1-tti+k+1-ti+1Nk-1i+1(t)(3.10)公式(3.10)称为de BoorCox公式。由公式(3.10)可以看出，基函数Nki(t)只有在位于节点ti和ti+k+1之间的区间内取非负值，而在其他处则取值为零。这也意味着在区间［ti,ti+1)上，总共只有k+1个B样条基函数取非负值，它们依次是Nki-k(t)，Nki-k+1(t)，…，Nki(t)。基于该递归公式，B样条曲线就可以采用类似Bézier曲线的递归方式进行计算，也就是通过逐次迭代的方式计算公式(3.8)表示的B样条曲线上任意点的位置坐标。这种计算方式称为de Boor递归算法。
该递归算法与Bézier曲线的de Casteljau递归算法类似，都是从控制顶点组成的控制多边形cj-kcj-k+1…cj开始，依次执行k次割角操作。其中，第r次割角是用线段p［r］i(t)p［r］i+1(t)割去角p［r-1］i。这里p［r］i(t)=(1-λi,k-r+1)p［r-1］i-1(t)+λi,k-r+1p［r-1］i(t)，而割角时线段端点在控制多边形边上的比例是λi,k-r+1=(t-ti)/(ti+k-r+1-ti)。那么，最后得到的角点p［k］j(t)就是B样条曲线上的点p(t)。整个递归割角过程可以用如下公式表示: 
p(t)=∑ji=j-kNki(t)p［0］i=∑ji=j-k+1Nk-1i(t)p［1］i…=∑ji=j-1N1i(t)p［k-1］i…=p［k］j(3.11)
其中，p［0］i=ci是控制顶点。事实上，该B样条曲线是定义在区间［tj,tj+1)上的一条多项式曲线。该区间也是公式(3.11)中能够使得所有B样条基函数取非零值的区间。
B样条曲线保留了Bézier曲线良好的几何性质，同时具备了对曲线形状的局部可控性。具体来讲，根据B样条基函数的局部支撑性，改动其中一个控制顶点，B样条曲线上仅仅和该控制顶点相关的曲线形状发生变化。此外，在进行多个B样条曲线拼接时，也可以根据节点设置很容易地保持拼接时的几何连续性。因此，B样条曲线具有更加灵活的几何建模能力。
例题34平面三次B样条曲线。
问题: 如右图所示的三次均匀B样条曲线p(t)的控制顶点分别是c0、c1、c2、c3，节点间距相等，设为ti=i。写出t∈［t3,t4)区间上的B样条曲线表达式。然后，给出采用de BoorCox递归算法计算曲线上任意一点的递归过程。
解答: 根据公式(3.10)中B样条基函数的定义，可以得到区间［t3,t4)上4个非零基函数的表达式: N30(t)=(-t3+3t2-3t+1)/6,N31(t)=(3t3-6t2+4)/6
N32(t)=(-3t3+3t2+3t+1)/6,N33(t)=t3/6因此，对应的三次均匀B样条曲线可以写为：p(t)=∑3i=0N3i(t)ci=16(t3,t2,t1,1)M4×4·c0
c1
c2
c3
M4×4=-13-31
3-630
-3030
1410对于曲线上任意一点p(t0)，3≤t0<4，其递归计算过程如下: p(t0)=p［3］3(t0)=(4-t0)p［2］2(t0)+(t0-3)p［2］3(t0)
p［2］2(t0)=4-t02p［1］1(t0)+t0-22p［1］2(t0)
p［2］3(t0)=5-t02p［1］2(t0)+t0-32p［1］3(t0)
p［1］1(t0)=4-t03p［0］0(t0)+t0-13p［0］1(t0)
p［1］2(t0)=5-t03p［0］1(t0)+t0-23p［0］2(t0)
p［1］3(t0)=6-t03p［0］2(t0)+t0-33p［0］3(t0)其中，p［0］0(t0)=c0，p［0］1(t0)=c1，p［0］2(t0)=c2以及p［0］3(t0)=c3。
□
进一步，B样条曲面是一种双参变量B样条函数组成的混合多项式曲面。B样条曲面可以看作B样条曲线在三维空间沿另外一条B样条曲线滑动形成。例如双三次B样条曲面的表达式为p(s,t)=∑3i=0∑3j=0Nki(s)Nkj(t)cij。与B样条曲线类似，B样条曲面也具有很多优良性质，例如局部性，即曲面形状只和最相关的几个控制顶点有关；凸包性，即B样条曲面的每一片都位于定义该片曲面的控制顶点的凸包之内；磨光性，即同一组控制顶点定义的B样条曲面，随着次数的升高越来越光滑。此外，Bézier曲面也可以看作B样条曲面的特例。
然而，B样条曲面无法直接表示圆锥等有理曲面。为此，研究人员又提出了非均匀有理B样条曲面（NURBS）。与普通的B样条曲线/曲面相比，NURBS的特点是采用非均匀节点，同时是一种有理表示形式。与此同时，加入了权重因子，进一步控制曲线/曲面形状（如图3.8所示）。具体来讲，NURBS具有如下表达式: p(s,t)=∑ni=0∑mj=0cijwijNki(s)Nkj(t)∑ni=0∑mj=0wijNki(s)Nkj(t)(3.12)其中cij是控制顶点，wij是权重因子。
图3.8NURBS定义的曲面形状
NURBS兼具B样条曲线/曲面形状局部可调以及连续阶数可调的优点，又能像有理Bézier曲线可精确地表示圆锥曲线的特性。1991年，国际标准化组织（ISO）在其正式发布的工业产品数据交换STEP标准中，把NURBS作为自由曲线/曲面的唯一定义，成为设计工业产品几何形状的唯一数学方法。许多国际著名的CAD软件也把NURBS作为几何造型工具的首选，例如AutoCAD、CATIA等软件。
3.3细分曲面建模
B样条曲线/曲面、NURBS等自由曲线/曲面在CAD中已得到广泛的应用。然而，这类自由曲线/曲面对它们自身的控制多边形或控制网格的拓扑结构有严格要求，通常来说，只能定义在矩形网格上。然而，现实中几何建模对象的拓扑结构往往是复杂多样的，例如计算机动画中要表示人的头和人的手。此外，因为模型是活动的，要在曲面的连接处保持光滑也是需要解决的问题。
针对上述问题，科研人员进一步发明了细分曲面。所谓细分曲面，是指多面体按照指定的细分规则进行无穷细化的极限。细分曲面可以看作自由曲线曲面在任意拓扑定义域上的推广。例如，B样条曲线的递归算法实际上就是对其控制多边形或控制网格的切割磨光过程。因此，人们自然地希望将这一算法推广到任意拓扑的多面体上去，也就是通过几何和拓扑的修改，重新添加边、顶点、面来重定义新的控制网格，以及移动顶点的空间位置来平滑该控制网格，并由此定义细分建模过程。这个过程所产生的极限曲面就是细分曲面。
细分曲面提供了更加灵活、更好光滑性的曲面生成方法，早期主要应用于计算机动画领域。最典型的是Pixar公司的Tony DeRose将细分曲面应用于电影动画Geris Game的制作，并获得了1998年奥斯卡最佳动画短片奖。对于细分曲面，最核心的是如何设计细分规则，也就是顶点、边、面的几何和拓扑生成方式。接下来主要介绍4种典型的细分规则及其生成的极限曲面: CatmullClark、DooSabin、Loop和Butterfly细分曲面。
3.3.1CatmullClark细分曲面
CatmullClark细分曲面将双三次B样条曲面的生成方法推广到任意拓扑的控制网格，由Edwin Catmull和Jim Clark于1978年首次提出。普通的双三次B样条曲面由16个控制顶点{cij}4,4i=1,j=1定义。根据B样条曲线/曲面的递归算法，在每两个节点的中点处嵌入一个新节点，就会产生25个新的控制顶点。这些新的控制顶点就定义了四片新的子曲面。这样，对应于原控制网格中的每一个小四边形（如c11c12c22c21），都有一个新顶点产生，称为面点；对于每一条边，也会有一个新顶点产生，称为边点；对于每一个顶点cij，也产生一个新顶点。这些新的面点、边点和顶点组成一个新的控制网格。该过程可重复进行下去，最终控制网格收敛到双三次B样条曲面。受此启发，Catmull和Clark提出了面向任意拓扑控制网格的细分规则,具体包括如下两方面的几何细分和拓扑细分规则。
（1） 几何规则。如图3.9所示，每个控制网格面在中心处加一个新的面点，如Q11；每条控制网格边加一个新的边点，如Q12；每个顶点用新的位置取代旧的位置，如Q22。
图3.9CatmullClark细分的几何规则和拓扑规则（图片来自［8］）



Q11=(c11+c12+c22+c21)/4
Q12=((Q11+Q13)/2+(c12+c22)/2)/2
Q13=(c12+c13+c22+c23)/4
Q22=Q/4+R/2+c22/4




新顶点Q22的计算公式中有Q和R两个新的变量，其中Q=(Q11+Q13+Q33+Q31)/4，R=(1/4)((c22+c12)/2+(c22+c21)/2+(c22+c32)/2+(c22+c23)/2))。
（2） 拓扑规则。如图3.9中的虚线所示，新的边是由连接每个面点到邻接的边点以及连接每个顶点到邻接的边点所形成。这些新的顶点和边就组成了新的控制网格。
与B样条曲面采用矩形拓扑的控制网格不同，CatmullClark细分曲面不受控制网格的拓扑限制，可作用到任意拓扑的控制网格上去。上述细分规则在作用一次以后，所有的面均变为四边形，而且从此以后度数不为4的顶点（称为奇异点）的个数会保持不变。除了奇异点以外，CatmullClark曲面可以看作由一系列双三次B样条曲面覆盖而成，通过上述的细分规则就能够达到几乎处处连续的曲率。而在奇异点处，也能够保持切平面连续，但需要对细分规则做一些相应调整。
例题35CatmullClark细分曲面。
问题: 右图所示的正立方体的8个顶点坐标分别是A(-1,1,1)、B(1,1,1)、C(1,1,-1)、D(-1,1,-1)、E(-1,-1,1)、F(1,-1,1)、G(1,-1,-1)和H(-1,-1,1)。写出采用CatmullClark细分生成新的控制网格顶点的伪代码。
解答: 假设k-1次细分后顶点集合是Vk-1={vk-11,vk-12,…}，其中V0={A,B,C,D,E,F,G,H}，在k-1次细分后所得的面和边的集合分别记为Fk-1={fk-11,fk-12,…}和Ek-1={ek-1ij}，那么第k次细分生成新的控制网格顶点的伪代码是function CCSubdivision (Vk-1, Fk-1, Ek-1)
for each face fk-1ijmn=(vk-1i,vk-1j,vk-1m,vk-1n) in Fk-1 do/计算新的面点/
v~kijmn← (vk-1i+vk-1j+vk-1m+vk-1n)/4
end for
for each edge ek-1ij=(vk-1i,vk-1j) in Ek-1 do/计算新的边点/
v~kij← ((vk-1i+vk-1j)/2+(vkijmn+vkijpq)/2)/2/vkijmn和vkijpq是面点/
end for
for each vertex vk-1i in Vk-1 do/计算新的顶点/
v~ki← Q/4+R/2+vk-1i/4/Q和R是新变量/
end for
end function□
3.3.2DooSabin细分曲面
DooSabin细分曲面将双二次B样条曲面的生成方法推广到任意拓扑的控制网格，由Daniel Doo和Malcolm Sabin在1978年最先提出的。普通的双二次B样条曲面在递归计算过程中，每一个面上对应于每一个顶点，产生一个新顶点。连接这些新顶点，就会产生对应于原控制网格中的面点、边点和顶点的新面，并由此形成不断加密的控制网格。最终，这个加密过程得到的极限曲面就是双二次B样条曲面。受此启发，DooSabin细分采用如下两方面的几何细分和拓扑细分规则: 
（1） 几何规则。每个控制网格面的K个顶点Q1，Q2，…，QK生成新的对应顶点Q′1，Q′2，…，Q′K，其中Q′k=∑Kj=1αijQj,αij=K+54K,i=j
3+2cos(2(i-j)π/K)4K,i≠j(3.13)这些新的顶点就定义了新的控制网格（如图3.10(a)所示）。
（2） 拓扑规则。如图3.10(b)所示，每个旧控制网格面的新顶点连接形成一个新的面F；每条旧边两侧的4个新顶点连接形成新的面E；每个旧顶点周围的新顶点连接形成新的面V。
图3.10DooSabin细分规则
经过一次DooSabin细分后，每个顶点的度数均变为4。再经过一次细分后，度数不为4的面的个数会保持不变。因此除了在有限个奇异点外，DooSabin细分曲面是由一系列双二次B样条曲面覆盖而成。此外，DooSabin曲面在奇异点处也具有一阶光滑连续性。
3.3.3Loop细分曲面
Loop细分曲面将箱样条推广到三角形组成的控制网格，由Charles Loop在1987年最先提出。相比于CatmullClark细分曲面和DooSabin细分曲面，Loop细分曲面的细分规则较为简单。通过生成新的顶点，并依次将其连接形成新的控制网格（如图3.11所示）。具体来讲，Loop细分采用如下几何和拓扑规则: 
（1） 几何规则。对于控制网格内部顶点Q0，假设其N个相邻顶点为Q1，Q2，…，QN，那么对应的新顶点Q′0=(1-Nβ)Q0+β∑Ni=1Qi，其中β=58-38+14cos2πN2N。对于控制网格边界上顶点Q0，假设与它相连的两个顶点是Q1和Q2，那么对应的新顶点满足Q′0=34Q0+18(Q1+Q2)。假设控制网格内部一条边的两个顶点是Q1和Q2，相对的两个顶点是Q3和Q4，那么新增加的顶点是Q′0=38(Q1+Q2)+18(Q3+Q4)。如果Q1和Q2是控制网格的边界边上的两个顶点，则对应于该边界边所新增加的顶点是Q′0=12(Q1+Q2)。
（2） 拓扑规则。按照三角形控制网格的顶点和边的连接关系，将新顶点连接形成新的控制网格。
Loop细分曲面除了一些特殊点外，几乎处处具有二阶光滑的连续性。
图3.11Loop细分规则
例题36Loop细分曲面。
问题: 右图所示正四面体的4个顶点的坐标分别为A(0,1,0)、B(0,0,1)、C(1,0,0)和D(0,0,0)。写出采用Loop细分生成新的控制网格顶点的伪代码。
解答: 假设k-1次细分后顶点集合是Vk-1={vk-11,vk-12,…}，其中V0={A,B,C,D}。在k-1次细分后所得的边的集合记为Ek-1={ek-1ij}，那么第k次细分生成新的控制网格顶点的伪代码是function LoopSubdivision (Vk-1, Fk-1, Ek-1)
for each vertex vk-1i in Vk-1 do
if vk-1i in vk-1i∈Vint then/计算内部顶点/
v~ki← (1-Nβ)vk-1i+β∑Nj=1vk-1ij/{vk-1ij}是相邻N个顶点/
else/计算边界顶点/
v~ki←34vk-1i+18(vk-1i1+vk-1i2)/vk-1i1和vk-1i2是相邻顶点/
end if
end for
for each edge ek-1ij=(vk-1i,vk-1j) in Ek-1 do
if ek-1ij in Eintthen/计算内部边界顶点/
v~kijpq← 38(vk-1i+vk-1j)+18(vk-1p+vk-1q)/vk-1p和vk-1q是相对顶点/
else/计算边界边顶点/
v~kij← 12(vk-1i+vk-1j)
end if
end for
end function□
3.3.4Butterfly细分曲面
Butterfly细分曲面是由Nira Dyn等人于1990年提出的，主要针对三角形组成的控制网格。具体来讲，对于每个边，使用指定的规则创造一图3.12Butterfly细分的几何
规则和拓扑规则
个新顶点，然后和旧的顶点把一个三角面片转化成四个新的三角面片（如图3.12所示）。Butterfly细分采用如下几何和拓扑规则: 
（1） 几何规则。对于三角形控制网格的每条边，利用可调控的权因子ω定义新的控制网格顶点Q′0=12(Q1+Q2)+ω(Q3+Q4)-ω2(Q5+Q6+Q7+Q8)。这里的权因子ω是可以根据细分曲面的外形要求而设定的。
（2） 拓扑规则。依次将新的顶点和旧控制网格顶点连接，形成新的控制网格。
根据上述几何和拓扑规则得到的Butterfly细分曲面，除了一些特殊点外，几乎处处具有一阶光滑的连续性。
上述细分曲面在细分过程中使用的几何和拓扑规则不同，由此产生不同性质的细分曲面。图3.13展示了同一个控制网格，采用不同细分规则后所生成的细分曲面，可以看出其结果在形状和光滑性上都有着明显的差异。
图3.13由相同控制网格生成的不同细分曲面
3.4三 维 重 建
自由曲线/曲面、细分曲面等几何建模方式是从数学模型出发来设计和构造几何形状，一定程度上要依赖设计人员关于几何形状的直觉和理解。与此相反，三维重建是直接或间接地获取真实世界中物体的形状和表观，使之能够在计算机中存储、处理和显示。因此，三维重建是一种从真实世界物体出发的几何建模方式。最具代表性的是Stanford大学的“米开朗基罗”项目。从1998年到2000年，该项目采用最先进的激光三维扫描仪将文艺复兴时代的著名雕塑作品全部进行几何建模，从而实现了文物的计算机数字化，起到了文物保护和文化传承的作用。
3.4.1被动式与主动式建模
按照数据来源的不同，采用三维重建进行几何建模主要分为两种方式: 被动式和主动式（如图3.14所示）。这两种都属于光学测量的范畴。其中，被动式包括基于图像的三维重建、基于视频的三维重建等；主动式包括基于激光测距的三维重建、基于Kinect的三维重建等。这两种方式各有优缺点，适用于不同的场合。
图3.14典型的三维重建方法
1. 被动式建模
被动式建模不需要与重建对象接触。在进行重建时，输入的数据是图像、视频等视觉信息，然后通过成像方式测量物体表面来推测三维结构。这类方法通常称为Shape from X，X可以指代轮廓、着色、纹理、阴影等能够通过设备获取的视觉信息。
被动式方法的优点是破坏性小、安全、成本较低。在数据采集时不需要和物体接触，所以不会对物体造成破坏，因此也避免了安全隐患。在建模时通常只需要根据拍摄的图像、视频数据作为输入，因而成本较低。但缺点是对物体透明度敏感、不能处理镜面反射和内部折射，因此无法对玻璃材质的物体进行有效的几何建模。
2. 主动式建模
主动式建模是在采集数据时通过机械接触或主动观测等方式获取三维信息，例如传感器标记、结构光、激光、超声波等。按照扫描方式的区别可以分为三维扫描仪、飞时测距、三角测距、结构光等。
接触式三维扫描仪，如坐标测量机等，具有测量精确度高的优点，可达到微米级别。但是这种方式价格昂贵，且需要专业的操作者。飞时测距，是通过发出激光脉冲，并计算这束光返回所需要的时间来测算距离。这种方式的优点是扫描速度快、便携、方便，而且测量范围大；但缺点是精度有限，只能达到毫米级别。三角测距，是通过发射一道激光到待测物上，并利用摄影机记录待测物上的激光光点，然后利用激光光点、摄影机、激光源构成的三角形来测算距离。这种方式的优点是精度较高、适合测量大尺寸物体；但缺点是扫描速度慢，需要花费较长时间来完成三维重建。结构光扫描，如Kinect等，则使用红外线发射器发射红外光线，然后利用红外线传感器接收反射回来的红外光线来获取深度图像。这种方式的优点是设备价格便宜、易于安装，但缺点是精度较差，而且测量范围有限（40cm～3.5m）。
3.4.2基于图像的三维重建
基于图像的三维重建（Imagebased reconstruction，IBR），是指从拍摄的图像对三维物体的外形进行重建。按照输入图像数量的不同，可分为多视角重建和单幅图像重建。多视角重建需要输入多幅图像，然后恢复物体的外形。单幅图像重建仅需要一幅图像作为输入进行重建。一般情况下，图像像素记录了来自不同方向的入射光信息。因此，增加输入图像数量能够提升几何模型重建的精度。
如图3.15所示，多视角重建的算法流程一般包括四个步骤: 摄像机标定、三角测量、从运动恢复结构（稀疏形状估计）以及立体匹配（稠密形状估计）。
图3.15多视角三维重建算法流程
1. 摄像机标定（Calibration）
该步骤从一些已知坐标的三维空间点集，反求相机的内参（如焦距等）和外参（相机的姿态）变量，进而得到与物体三维信息有关的相机运动。摄像机的成像模型通常简化为小孔成像模型，那么从三维空间到二维成像平面的成像可以看作投影变换。因此，摄像机的成像模型表示为: x3×1=P3×4·X4×1x
y
1=p11p12p13p14
p21p22p23p24
p31p32p33p34X
Y
Z
1(3.14)其中，x是二维成像平面上的像素坐标，X是三维空间中的点坐标，P3×4是描述小孔成像过程的投影矩阵。为了方便统一的表示，这里采用齐次坐标形式。摄像机标定的主要任务就变成计算投影矩阵P3×4，进一步可以分解为内参矩阵K3×3和外参矩阵M3×4，即P3×4=K3×3·M3×4。
内参矩阵描述了在摄像机坐标系下，三维空间到二维图像的投影变换。假设摄像机位于原点，朝向与坐标轴重合，而图像坐标系的原点在图像中心，那么内参矩阵K3×3可以表示为一个对角矩阵和上三角矩阵的乘积: K3×3=mx
my
1fpx
fpy
1(3.15)其中，mx和my是由于透镜弯曲产生的畸变系数，f是焦距，(px,py)是图像中心的像素坐标。外参矩阵则描述了摄像机姿态的变化对不同视角成像的影响。通常把摄像机看作刚体，其姿态变化就可以通过旋转和平移变换来实现，也就得到M3×4=［R|t］。其中R3×3是三维空间中的旋转矩阵，t是平移向量。通过该矩阵，就可以描述摄像机在拍摄两幅图像时的空间位置关系。
而在实际中，不同视角图像拍摄时摄像机姿态是不同的，这就需要直接根据多幅图像之间的对应关系计算相应的内参和外参矩阵。这个过程通常包含4个步骤: 角点检测、投影变换计算、参数估计、参数优化。
角点检测是为了获取图像中具有视觉显著性的特征点，然后通过特征点匹配建立不同图像之间的对应关系。为了提高检测和匹配的准确性，通常采用固定模式的物体做参考图像，例如棋盘格。首先进行边缘检测，将各个矩形框边缘拟合成直线；然后，计算直线交点作为角点位置。这样通过角点所在直线的行列排布，就可以直接建立不同图像之间角点的对应关系。进一步，将棋盘格摆放在摄像机前的指定位置，就能获得这些角点在三维空间中的坐标。
投影变换计算则是通过角点的二维图像坐标与三维空间坐标求解变换矩阵P。常用的计算方法有直接线性变换法（Direct linear transform，DLT）、最小二乘优化等，其本质是优化三维空间点在投影变换后与图像角点的几何误差。
参数估计是通过对投影变换的分解，估计内参矩阵和外参矩阵的初始值。根据内参矩阵和外参矩阵的形式，将3×4的投影变换矩阵P表示为3×3的矩阵B和3×1的向量b，也就是写成P=［B|b］，那么，内参矩阵可以表示为K=B·BT。进一步，对B进行QR分解可以得到矩阵R和Q，其中R就对应于外参矩阵的旋转变换，平移变换对应的向量则可以通过t=Q-1·b来计算。
上述参数估计是对单张图像的角点投影误差进行优化求解，而不同视角下的内参矩阵需要保持一致。另一方面，通过矩阵分解计算得到的内外参数，本质上是通过优化代数距离进行的计算，并没有实际的物理意义。这也就是说，图3.16重投影误差
获得的解是有几何偏差的。因此，需要对多视角投影误差做进一步优化，以便进一步提高摄像机标定的准确性。
给定多视角下n组对应角点的二维图像坐标和三维空间点坐标，计算角点检测位置与其通过投影变换模型预测的成像点之间的重投影误差，将其做最小化处理（如图3.16所示），建立关于矩阵元素的优化模型。该过程可以通过非线性最小二乘优化进行计算，也就是求解能量函数E(K,R,t)=minP∑i‖xi-K［Rt］Xi‖2的局部最优解。数学上，这个优化问题可以使用LM算法（LevenbergMarquardt algorithm）进行有效求解。
2. 三角测量（Triangulation）
通过指定棋盘格图像对摄像机进行标定后，就可以根据任意拍摄的图像对物体进行三维重建。首先需要通过标定的内参和外参矩阵，并结合两幅或多幅图像间特征点的对应关系，计算相应的三维空间坐标。三角测量是一种最简洁的方法。它利用已知空间点X在两幅图像上的投射点的坐标和两条投射线与基准线的夹角，估计这些点在三维空间中的位置。基准线是指两个相机坐标系原点的连线，那么形成的三角形的第三个点就是点X的三维坐标。但是由于成像噪声及计算误差的原因，两条直线在三维空间中可能不会产生严格意义的相交。这时，可通过求解距离两条射线最近的点作为近似交点。
3. 从运动恢复结构（Structurefrommotion，SFM）
将摄像机标定和三角测量的结果作为初始值，通过分析物体在不同视角成像的运动，就可以得到更加准确的三维模型的结构信息。这一过程通常是利用同一个三维空间点投影到不同视角图像上二维特征点的投影关系来实现，同时也得到摄像机参数和三维空间坐标。
图3.17光束平差法的重投影误差光速平差法（Bundled adjustment）是常用的方法之一。它是基于三维模型结构和视角参数（如相机位置、朝向、固有标定和径向畸变）的优化问题，用于获得最佳的三维坐标和运动（如相机矩阵）参数估计。这种方法本质上是最小化投影变换所决定的重投影误差。具体来讲，对应于如图3.17所示的一个空间点的重投影误差，可以采用如下公式进行优化求解:
min∑n1∑m1ωij‖Xij-P(Oi,Xj)‖2(3.16)
其中Xij表示三维空间中的点Xj在第i个相机投影后的像素位置，ωij是和相机位置分布有关的权因子，P(·,·)表示对应的相机投影操作。公式(3.16)可以使用最小二乘算法来进行最小误差的计算，最终得到最优的投影矩阵和三维模型。