第5章曲线和曲面

第5章 曲线和曲面曲线和曲面可分为由已知方程式表示的规则曲线或曲面和用曲线或曲面拟合方法逼近的不规则曲线或曲面。这些不规则曲线或曲面一般采用分段的多项式参数方程来表示，由此形成一条光滑连续的曲线或曲面称为样条曲线或曲面。本章将主要讨论参数样条曲线和曲面的绘制方法。

5.1曲线和曲面基础从卫星的轨道、导弹的弹道，到汽车和飞机等的外形，直至日常生活中的图案和花样设计，都离不了对曲线的描述和绘制。以至于可以说，几乎没有一张设计图纸上是没有曲线的。在我们遇到的各种各样的曲线中，归纳起来，大概不外乎两类：一类是我们已经比较熟悉的，如圆、椭圆、双曲线、正弦余弦、概率分布、摆线螺线等等。这类曲线均可以用一个曲线方程式来表示，称此类曲线为规则曲线。比如圆的方程可以写成x2＋y2=R2等。而另有一类曲线，我们尚不能确切给出描述整个曲线的方程，它们往往是由一些从实际中测量得到的一系列离散数据点用曲线拟合方法来逼近的，称为不规则曲线。这些曲线一般采用分段的多项式参数方程来表示，由此形成一条光滑连续的曲线称为样条曲线或简称样条。常见的参数样条曲线有抛物样条曲线、Hermite插值样条曲线、Bezier样条曲线和B样条曲线等。当曲线的数学表达方法确定以后，剩下的问题就是如何把这些曲线绘制出来。要绘制一条指定的曲线函数的直接方法是用很多短直线段来逼近曲线。绘出的曲线的光滑度和精确度取决于我们所选择的数据点的精度和数量。点的数量越多，直线段越短，则连成的曲线愈接近于理想曲线。至于点的数量取多少，直线段取多长，则取决于我们对所绘制曲线的精度要求和图形输出设备的精度，但我们对所绘制曲线的精度要求不能逾越图形输出设备所实际具有的精度。5.1.1规则曲线和曲面的3种坐标表示法

一般平面曲线常用直角坐标、极坐标或参数方程表示，如工程上常用的渐开线、摆线以及正弦余弦曲线等，都是用这三种坐标表示的重要曲线。但从计算机图形学和计算几何的角度看，用参数方程绘制曲线比较方便。实际上绘制任何平面曲线．都要将曲线方程用参数方程形式表示，即得到曲线上点坐标x与y的分别计算式。于是计算出点的坐标值，调用画线函数或画点函救绘出曲线上的所有点，便得到一条曲线。

1．直角坐标曲线曲线的直角坐标表示，有显式y=f(x)和隐式f(x，y)

=

0之分。如y=

sin

(x)是显式表示，而x2

+y2

=

1是隐式表示。无论是哪种表示，都要将其转换成参数坐标表示．即x=x(t)

y=y(t)然后可以开始绘制它的图形了。下面对曲线的直角坐标显式和隐式两种表示分别加以讨论。

（1）显式对于显式表示y=f(x)的曲线转换成参数坐标表示，这是非常容易的，即y=k·x+bx=xy=f(x)这里式子右边的x看成参数变量。此两式便是显式表示y=f(x)曲线的参数坐标表示。例如正弦曲线y=

sin(x)

是直角坐标显式表示．其参数坐标表达式为

x=xy=

sin(x)

这时式子右边的x当作参数变量。这样给定一个参数变量x值，就可求得正弦曲线上一个点的坐标x与y值。进而一点一点地绘出正弦曲线。

（2）隐式一般隐式f(x，y)

=

0的曲线转换成参数坐标表示式是很困难的，如下面隐式曲线4x4–3x3

+

2y2–y–x2（a+x）/（a–x）=0

（a>0）

要表示成参数坐标式，至今未能成功，因此无法使用计算机绘制它的图形。不过常用的重要曲线基本上都能用参数坐标表示。例如星形线直角坐标表示式：

x2/3＋y2/3=R2/3

（R正常数）

可写成参数坐标表示式：

x=Rcos3θy=Rsin3θ（0≤θ≤2π）从而可用计算机绘出其曲线图。

F(x,y)=ax+by+c=02．极坐标曲线对任一极坐标曲线ρ=ρ(θ)，可利用极坐标与直角坐标变换关系式

x=ρcosθy=ρsinθ

将此曲线转换成参数坐标表示为

x=ρ(θ)cosθy=ρ(θ)sinθ这里θ成为参数坐标。

例如，重要曲线阿基米德螺线

ρ=aθ

（a正常数）

极坐标与直角坐标变换关系式

x=ρcosθy=ρsinθ

将阿基米德螺线ρ=aθ代入上面两式，便得

x=aθcosθy=aθsinθ这样就将阿基米德螺线极坐标表示转换成了参数坐标表示。

由阿基米德螺线参数坐标表示式

x=aθcosθy=aθsinθ可以计算出其曲线上点的坐标值，然后用这些点的坐标值调用绘图函数就可绘出阿基米德螺线曲线图。

3．参数坐标曲线曲线的参数坐标表示一般为

x=x(t)y=y(t)如弹道曲线

x=V0tcosαy=V0tsinα–gt2/2

（0≤t≤2V0Sinα/g）式中V0，g，α均为常数，t为参数变量。

对于给定一个参数变量t值，就可求得弹道曲线上一个点的x与y坐标值。若给出参数变量t的一系列值，便可求出弹道曲线上一系列点的x和y坐标值。有了x,y坐标值,使用绘图函数把这一系列点绘制出来,就获得一条弹道曲线。

P(t)=P1+(P2-P1)t(0≤t≤1)x(t)=x1+(x2-x1)ty(t)=y1+(y2-y1)t对于某一参数曲线，我们不可能也没有必要去研究参变量t从–∞到+∞的整条曲线，而往往只对其中的某一段感兴趣。通常我们经过对参变量t的规格化变换，使t在[0，1]闭区间内变化，写成t∈[0，1]，对此区间内的参数曲线进行研究。

4．参数曲线的优点在曲线的表示上，参数方程比显式、隐式方程有更多的优越性。

（1）有更大的自由度来控制曲线的形状。如一条二维三次曲线的显式表示为：

y=ax3＋Bx2＋cx＋d其中只有4个系数可用来控制此曲线的形状。而二维三次曲线的参数表达式为：

x=at3＋Bt2＋ct＋dy=et3＋ft2＋gt＋h其中有8个系数可用来控制此曲线的形状。

（2）对非参数方程表示的曲线进行变换，必须对曲线上的每个型值点进行几何变换；而对参数表示的曲线可对其参数方程直接进行几何变换（如平移、比例、旋转），从而节省计算工作量。

（3）便于处理斜率为无限大的问题，不会因此而中断计算。

（4）规格化的参数变量t∈[0，1]，使其相应的几何分量是有界的，而不必用另外的参数去定义其边界。

（5）参数方程中，代数、几何相关和无关的变量是完全分离的，而且对变量个数不限，从而便于用户把低维空间中的曲线扩展到高维空间去。这种变量分离的特点使我们可以用数学公式去处理几何分量，如我们以后用的调和函数就具有此特点。

（6）易于用矢量和矩阵表示几何分量，简化了计算。基于这些优点，我们在以后将用参数表达式来讨论曲线问题。5.1.2参数样条曲线和曲面的常用术语在工程设计中，一般多采用低次的参数样条曲线。这是因为高次参数样条曲线计算费时，其数学模型难于建立且性能不稳定，即任何一点的几何信息的变化都有可能引起曲线形状复杂的变化。

因此，实际工作中常采用二次或三次参数样条曲线，如：二次参数样条曲线：P(t)=A0

+A1t+A2t2三次参数样条曲线：

P(t)=A0

+A1t+A2t2

+A3t3

1．型值点和控制点所谓型值点，是指通过测量或计算得到的曲线上少量描述曲线几何形状的数据点。由于型值点的数量有限，不足以充分描述曲线的形状，因此通常是在求得一些型值点后，采用一定的数学方法，建立曲线的数学模型，从而再根据数学模型去获得曲线上每一点的几何信息。所谓控制点，是指用来控制或调整曲线形状的特殊点，曲线段本身不通过该控制点。

1)切线如图5-1，设空间中的一条参数曲线为P(t)，M、Q为曲线上邻近的两点，分别对应于p(t)与p(t+Δt),当曲线上的点Q趋于M时，割线的极限位置为：2．切线、法线和Frenet标架称为曲线在点M处的切线，而定点M称为切点，切线方向称为切向量或切矢量T。直观上看，切线是通过切点的所有直线当中最贴近曲线的直线。图5-1曲线的切线定义若参数曲线上任一点的坐标为p(t)=[x(t)，y(t)，z(t)]，则该点的切线方程即为参数曲线在该点处的一阶导函数，即p'(t)=[x'(t)，y'(t)，z'(t)]。如果选择弧长s作为参数，我们还可以得到：

由于在极限情况下，|ΔP|=|ΔS|,所以上面式(5-4)和（5-5）中的T是单位长度矢量，其方向为曲线的切线方向，称之为单位切矢量，T的几何意义是：T的方向角度变化反映出曲线的平坦程度,当T的方向恒定不变时，曲线则变为一条直线，参见后面图5-3(b)、(c)。过上述空间曲线P(t)上M点的切线和M点的邻近一点Q还可作一平面，当Q点沿着曲线趋近于P时，平面的极限位置称为曲线在M点的密切平面，如图5-2所示。2)法线法线就是垂直切线方向T且通过切点M的直线，在密切平面上的法线称为主法线，所在的方向称为主法向量N，而垂直于密切平面还存在一条法线，如图5-2所示，该法线称为副法线，所在的方向称为副法向量B。主法向量与副法向量构成法平面，而切向量与副法向量构成从切平面。3)Frenet标架假设N和B均为单位向量，我们把两两相互垂直的单位向量T、N和B称为曲线P(t)上M点的Frenet标架，如图5-2所示。图5-2曲线上任意点的Frenet标架构成图5-1曲线的切线定义yxQMαdαα+dαds(a)曲率的定义(b)曲率的几何意义(c)不同曲率的圆弧曲线比较图5-3曲线的曲率定义与计算3．曲率和挠率1)曲率如图5-3(a)所示，曲线上两点M和Q的切线的夹角δ与弧长MQ之比，当Q趋于M时的极限，即

上式(5-6)的推导表明，曲率既是也是切线的方向角对于弧长的转动率，又是单位切向量对于弧长的旋转速率，参见图5-3(b)。当曲线在一点的弯曲程度越大，切向量对于弧长的旋转速度就越大，因此曲率刻画了曲线的弯曲程度，如图5-3(c)所示。式(5-7)表明，其值实际上为曲线在M处的二阶导数。

一般参数曲线P(t)曲率的计算公式为：

而曲率半径的计算公式则是：ρ=1/k。(5-8)1)曲率对于空间曲线，曲线不仅弯曲，而且还要扭转，挠率便是衡量曲线扭转变形程度的一种度量。具体地说，曲线挠率的大小表现出其离开密切平面的程度大小。空间曲线上任意一点的挠率可以由下列公式表达为：

挠率表达式(5-9)的几何意义是：挠率即为曲线在副法方向上随弧线长度变化的速率，它是一般曲线变形的重要表述参数。不难理解：挠率恒等于零的曲线是平面曲线，而曲率恒等于零的曲线则为直线。另外，切线方向恒定不变的曲线也必定是直线。(5-9)2）挠率而曲线的拟合则是这两种设计方法的统称，是指在曲线的设计过程中，用插值或逼近方法使生成的曲线达到某些设计要求，如在允许的范围内贴近原始的型值点或控制点序列，或曲线看上去很光滑等。

4．参数连续性和几何连续性

为保证分段参数曲线从一段到另一段平滑过渡，我们可以在连接点处要求各种参数连续性条件。

0阶参数连续性，记作C0连续，是指曲线相连，即第一个曲线段的终点与第二个曲线段的起点相同。一阶参数连续性，记作C1连续性，指代表两个相邻曲线段的方程在相交点处有相同的一阶导数（切线）。二阶参数连续性，记作C2连续性，是指两个曲线段在交点处有相同的一阶和二阶导数。

4．参数连续性和几何连续性

1）参数连续性为保证分段参数曲线从一段到另一段平滑过渡，我们可以在连接点处要求各种参数连续性条件。0阶参数连续性，记作C0连续，是指曲线相连，即第一个曲线段的终点与第二个曲线段的起点相同。一阶参数连续性，记作C1连续性，指代表两个相邻曲线段的方程在相交点处有相等的一阶导数(切线)。二阶参数连续性，记作C2连续性，是指两个曲线段在交点处有相等的一阶和二阶导数。参数连续性和几何连续性示意图2）几何连续性连结两个相邻曲线段的另一个方法是指定几何连续性条件。这种情况下，只需两曲线段在相交处的参数导数成比例而不是相等。0阶几何连续性，记为G0连续性，与0阶参数连续性相同，即两个曲线段必在公共点处有相同的坐标。一阶几何连续性，记为G1连续性，指一阶导数在两个相邻段的交点处成比例但不一定相等，方向相同。二阶几何连续性，记为G2连续性，指两个曲线段在相交处其一阶和二阶导数均成比例但不一定相等，

方向相同。在实际的曲线造型应用中，我们要适当地选择曲线段间的连续性，使造型物体既能保证其光滑性的要求，也能保证其美观性的要求。4．参数连续性和几何连续性

4．插值、逼近和拟合及光顺插值与逼近是曲线和曲面设计中的两种不同方法。插值设计方法要求建立的曲线或曲面数学模型，严格通过已知的每一个型值点。常用的插值方法有线性插值和抛物线插值两种。而逼近设计方法，顾名思义，用这种方法建立的曲线或曲面数学模型只是近似地接近已知的型值点。常用的逼近方法有最小二乘法等。而曲线或曲面的拟合则是这两种设计方法的统称，是指在曲线或曲面的设计过程中，用插值或逼近方法使生成的曲线或曲面达到某些设计要求，如在允许的范围内贴近原始的型值点或控制点序列，曲线或曲面看上去很光滑等。光顺就是光滑、顺眼的意思，“光滑”通常指曲线曲面的参数连续性或几何连续性，主要是从数学的角度来考虑，有严格的数学定义；“顺眼”则可以理解为侧重于美学、功能、力学等诸多方面的要求，如曲线的拐点不能太多，拐来拐去就会不顺眼。对于平面曲线而言，一般的光顺准则是：二阶几何连续（G2连续）、没有多余拐点、曲率变化比较均匀。其中“拐点”的意思是指曲线的凸凹发生改变的点，“拐来拐去”的典型曲线如常见的正弦曲线。4．插值、逼近和拟合及光顺光顺示意图5．参数曲线的代数形式、几何形式及调和函数的概念

1）参数形式：P(t)=P1+(P2-P1)t(0≤t≤1)对比x(t)=a1t+a02）几何形式：P(t)=(1-t)P1+tP2(0≤t≤1)3）调和函数f1(t)=1-t;f2(t)=t,f1(t)+f2(t)=15.2二次插值样条曲线

5.2.1二次插值样条曲线的数学表达式在拟合生成样条曲线的众多方法中，我们首先选择较为简单的二次样条曲线即抛物样条曲线的生成方法作为基本方法，来讨论如何用插值方法生成通过给定离散型值点的样条曲线。实际上，二次Bezier曲线和二次B样条曲线也是抛物样条曲线，但它们采用的方法是逼近方法。

由于离散点的要求，我们首先要解决由给定点定义抛物线问题。设有不在同一直线上的三点：P1,P2,P3，现在要求通过该给定的三点定义一条抛物线。如图5.2所示。P1P2P3图5.2过三点的二次曲线

假如我们采用矢量表达式来表示参数化的二次曲线，那么可以把抛物线的表达式写成如下的一般形式为：

P（t）=A1

+A2t+A3t2

（0

≤t≤1）

(5-1)抛物线是一条二次曲线，所以表达式中参数t的最高次数为2，同时让参数t在0—l之间取值。

这就是说，只要确定了式(5–1)中的三个系数：A1,A2和A3，那么就确定了抛物线的表达式，随之抛物线的曲线图形也就可以确定。所以，我们的工作是要通过设定一些已知条件来求出这三个系数。

要确定这三个系数(目前尚为未知数)，必须要有三个独立的条件。我们可以给定这三个独立条件为：

该抛物线过P1，P2，P3三个点，并且：①抛物线段以P1点为始点。即当参变量t=0时，曲线过P1点；②抛物线段以P3点为终点。即当参变量t=1时，曲线过P3点；③当参变量t=0.5时，曲线过P2点，且切矢量等于P3–P1。在这三个设定的条件下，构造的抛物线段如图5.3所示。

P1P2P3QAP2’t=0t=0.5t=1图5.3过3点定义的二次曲线图中的数据是这样的：A点为P1P3的中点，AP2=P2Q，抛物线在P1点处与P1Q相切，在P3点处与QP3相切，曲线在P2点处的切矢P'2与P1P3平行。

根据以上设定的条件，可以列出三个方程：

t=

0：P(0)=A1

=P1t=

1：P(1)=A1＋A2十A3

=P3(5-2)

t=

0.5；P(0.5)=A1＋0.5A2＋0.25A3

=P2解以上三个联立方程：

A1

=P1P3

=A1

+A2

+A3

=P1＋A2

+A3

∴A2

=P3–P1–A3P2

=A1

+

0.5A2

+

0.25A3

亦即：4P2

=

4A1

+

2A2

+A3

=

4P1

+

2(P3–P1–A3)

+A3

=

2P1

+

2P3

–A3∴A3

=

2P1

+

2P3–4P2以上式回代到A2

=P3–P1–A3中，得：

A2

=

4P2–P3–3P1

所以，通过解联立方程，得到的三个系数A1，A2，A3分别为：

A1

=P1A2

=

4P2–P3–3P1(5-3)A3

=

2P1＋2P3–4P2把求出的该三个系数的值，代入到抛物线的表达式(5-1)中，可得：

P(t)=A1

+A2t+A3t2

=P1

+(4P2–P3–3P1)t+

(2P1＋2P3–4P2)t2

=

(2t2–3t+

1)P1

+(4t–4t2)P2

+

(2t2–t)P3(5-4)

(0≤t≤1)

可把式(5-4)改写成矩阵形式为：P(t)=[t2

t1](5-5)以上推导求出的算式，即为我们所要求的过不在一直线上的三点：P1(x1,y1)，P2(x2,y2)和P3(x3,y3)的抛物线方程。这时根据参变量t的取值，我们就可以一一计算出位于曲线上的数据点，然后顺次连线绘出图形。

5.2.2二次插值样条曲线的加权合成设有一离散型值点列Pi(i=1,2,…,n)，我们可以按式(5-5)每经过相邻三点作一段抛物线，由于有n个型值点，所以像这样的抛物线段一共可以作出n–2条。如图5.4所示。

P1P2P3P4P5Pn-2Pn-1Pn图5.4产生n–2段抛物线在这n–2条抛物线段中，第i条抛物线段为经过Pi，Pi+1，Pi+2三点，所以它的表达式应为：

Si(ti)=(2ti2–3ti+1)Pi+(4ti–4ti2)Pi+1+(2ti2–ti)Pi+2

(0≤ti≤1)

(5-7)同理，第i+1条抛物线段为经Pi+1，Pi+2和Pi+3三点，所以它的表达式为：Si+1(ti+1)=(2ti+12–3ti+1+1)Pi+1+(4ti+1–4ti+12)Pi+2+(2ti+12–ti+1)Pi+3

(0≤ti≤1)（5-8）

经过四点所画出的两条抛物线段Si(ti)和Si+1(ti+1)的图形如图5.5所示：

PiPi+1Pi+2Pi+3SiSi+1图5.5Si和Si+1

一般说来，每两段曲线之间的搭接区间，两条抛物线是不可能重合的。例如图5.5中，

Si和Si+1两条抛物线，它们在Pi+1和Pi+2两点之间为搭接区间，在这区间内，Si和si+1不太有可能会自然地重合成一条曲线。

显然，对于拟合曲线来说，整个型值点列必须只能用一条光滑的曲线连接起来。为了做到这一点，在Si和si+1这样两条曲线的共同区间内，必须有一个办法让它们按照一个一定的法则结合成一条曲线，这结合的办法就是加权合成。

在加权合成的过程中，我们首先要选择两个合适的权函数。如果我们在这里选择的两个权函数分别为f(T)和g(T)，加权合成后的曲线为Pi+1(t)，则：

Pi+1(t)

=f(T)·Si(ti)+g(T)·Si+1(ti+1)在抛物样条曲线中，我们选择的权函数f(T)和g(T)是简单的一次函数，且它们之间存在有互补性。它们分别为：

f(T)=l–T

g(T)=T(O≤T≤1)这样，式：

Pi+1(t)=f(T)·Si(ti)＋g(T)·Si+1(ti+1)即可改写为：

Pi+1(t)=(1–T)·Si(ti)＋

T·Si+1(ti+1)(5-9)在表达式（5-9）中，包含了三个参变量，即：T、ti和ti+1。假如这三个参变量不加以统一，则接下去的工作是无法进行的，所以我们首先要统一式中的参变量。

对于曲线段Si(ti)，参变量ti的取值范围为：0≤ti≤1，但曲线段Si(ti)与曲线段Si+1(ti+1)搭接的部分是原曲线段的后半截，即是从点Pi+1到Pi+2之间的区间，在这个区间内，参变量的取值范围应为：0.5≤ti≤1。

同理，对于曲线段Si+1(ti+1)，在点Pi+1到Pi+2之间的区间内，其参变量ti+1的取值范围应为：0≤ti≤0.5。在权函数f(T)和g(T)中，变量T的取值范围定为：0≤T≤1。

为了统一式(5-9)中的三个参变量：T、ti和ti+1，我们选择t作为统一后的参变量，把原有的三个参变量T、ti和ti+1均化成唯一含有t的形式，并要给t规定一个合适的取值范围。假如我们使t的取值范围为：0≤t≤0.5，则上面的三个参变量可统一形式为：

T=2t

ti=0.5＋t0≤t≤0.5

ti+1=t于是，原式(5-9)可根据新的参变量t改写成如下这样的形式：

Pi+1(t)=(1–2t)·Si(t+0.5)＋2t·Si+1(t)(5-10)其中：1–2t=f(T)2t=g(T)

Si(t＋0.5)=(2t2–t)Pi+(1–4t2)Pi+1+(2t2＋t)Pi+2

Si+1(t)=(2t2–3t＋1)Pi+1+(4t–4t2)Pi+2+(2t2–t)Pi+3

把以上四式代入式（5-10），展开、整理后可得：

Pi+1(t)=(–4t3+4t2–t)Pi+(13t3–10t2+1)Pi+1

+(–12t3+8t2＋t)Pi+2+(4t3–2t2)Pi+3(i=1，2，…，

n–3)(0≤t≤0.5)

(5-11)

式(5-11)的实质是：表达了每相邻的四个点可以决定中间的一段抛物样条曲线。见图5.6所示。

PiPi+1Pi+2Pi+3Pi+1(t)图5.64个点决定中间的一段样条曲线

假如一个离散点列Pi具有n个型值点，即i=1，2，…，n。那么根据式(5-11)，可以加权合成后生成n–3段抛物样条曲线。即式(5-11)中的i的取值范围为：i＝1～n–3。

5.2.3二次插值样条曲线的端点条件上面已经说到，在全部点列Pi(i=1，2，…，n)中，我们只能得到n–3段曲线。但n个型值点之间应有n–1个区段。亦即，按照式(5-11)产生的曲线段不足以生成n–1段曲线，因其点列的首、尾两段曲线P1P2和Pn–1Pn段，由于缺乏连续相邻的四点这样的条件而无法产生。

为了要产生首尾两段曲线，一个直接的想法就是在原点列的两端各加一个辅助点P0和Pn+1，如图5.8所示。

P0P1P2Pn-1PnPn+1图5.8两端加点但是，余下的问题是这P0和Pn+1两点是如何加上去的，它必须依据什么原则，这就是所谓的“端点条件”。在这里，我们仅介绍常用的三种方法：

①已知两端的切矢P'1和P'n在前面我们已经说过，在由P1、P2、P3三点所确定的抛物线中，过P2点曲线的切矢P'

2

=P3–P1，即：P1=P3–P'

2这样，在抛物样条曲线中，当条件给出了两端的切矢P'1和P'n之后，根据上面的原理可得：

P'1=P2–P0∴P0=P2–P'1

P'n

=Pn+1–Pn–1∴Pn+1=Pn–1+P'n即可以确定辅助点P0和Pn+1的坐标位置。

这种端点的情况，一般适用于所求的曲线要和已经存在的曲线或直线相连接。

②自由端条件另一种补点的方法的原理比较简单，它让所补之点P0和Pn+1与原两端点P1和Pn分别重合，即：

P0=P1

Pn+1=Pn

这样的补点方法称为自由端条件，这种方法一般适用于对曲线的两端没有什么特殊的要求。

③形成封闭曲线为了在n个型值点之间形成封闭曲线，那末就要生成n段曲线段，而不是原来的n–1段。所以在补点工作中要加三个点，首先让首尾两点重合，然后各向前后延长一点，即：

Pn+1=P1

P0=Pn

Pn+2=P2

上述由加权合成推导的抛物样条曲线是由若干个曲线段组成的，每两个相邻的型值点之间形成一段曲线，每相邻的两段曲线在型值点处相接，并且一阶导数相等，因此抛物样条曲线达到C1连续。

5.2.4二次插值样条曲线的性质

抛物样条曲线的算法程序如下：voidparspl（intp[][2]，intn，intk，inte）｛//P为型值点的坐标数组，n为型值点数，k为插值数，即是把参变量t区间细分的份数。intx，y，i，j，m=n；floatt1，t2，t3，t，a，b，c，d；if（e==1）//自由端{p[0][0]=p[1][0]；p[0][1]=p[1][1]；p[n+1][0]=p[n][0]；p[n+1][1]=p[n][1]；}else//画封闭曲线

{p[0][0]=p[n][0]；p[0][1]=p[n][1]；p[m][0]=p[1][0]；p[m][1]=p[1][1]；p[m+1][0]=p[2][0]；p[m+1][1]=p[2][1]；}t＝0.5／k；Moveto(p[1][0],p[1][1])；for(i=0；i<m–1；i++){for(j＝1；j<k；j++)

｛t1＝j*t；t2＝t1*t1；t3＝t2*t1；a＝4.0*t2–t1–4.0*t3；b＝1.0–10.0*t2＋12.0*t3；c＝t1+8.0*t2–12.0*t3；d＝4.0*t3–2.0*t2；x＝a*p[i][0]+b*p[i+1][0]+c*p[i+2][0]+d*p[i+3][0]；y＝a*p[i][1]+b*p[i+1][1]+c*p[i+2][1]+d*p[i+3][1]；lineto(x,y)；

｝lineto(p[i+2][0],p[i+2][1])；}}5.3三次插值样条曲线

插值样条曲线大多用来建立物体运动路径或提供实体表示方法，有时也可用来设计物体形状。三次插值样条曲线在灵活性和计算速度之间提供了一个合理的折中方案。与更高次样条相比，三次插值样条只需较少的计算和存储，且较稳定。与低次插值样条相比，三次插值样条在模拟任意形状时显得更灵活。

三次插值样条曲线由分段的三次多项式来描述。假如我们设其参变量为t，则分段三次插值样条曲线表达式的一般形式可以写成：

P(t)=B1＋B2t＋B3t2十B4t3

(0≤t≤tm)(5-13)其中的P(ti)=[x(ti)y(ti)z(ti))，可以看作是插值样条曲线上某一点的位置向量，ti是该点相应的参变量，它的三个分量：x(ti),y(ti),z(ti)可以看作是点的坐标值。

式(5-13)中的四个系数B1,B2,B3和B4是待定的，即是说方程中有四个未知数。为了把式(5-13)具体化，必须确定这四个系数，显然，这需要设定四个独立的条件。n+1个控制点共产生n段曲线，每段曲线都有4个系数需要确定。不同的系数确定方法就导致不同的三次插值样条曲线。

5.3.1三次自然样条曲线最先用于图形应用的三次插值样条曲线是三次自然样条曲线。这个插值样条曲线是原始绘图样条的一个数学表达式。用公式表示一个三次自然样条曲线时，需要两个相邻曲线段在公共边界处有相同的一阶和二阶导数，即三次自然样条曲线具有C2连续性。

对于三次自然样条曲线，n+1个型值点插值产生n段曲线，这样共有4n个多项式系数需要确定。对于每一个内型值点(共有n-1个)，有四个边界条件：在该型值点两侧的两个曲线段在该点处有相同的一阶和二阶导数，且两个曲线段都要通过该点。

这给出了由4n个多项式系数组成的4n-4个方程，另外再给出一个从第一个型值点p0(即曲线起点)所得的方程，和从最后一个型值点pn(即曲线终点)所得的方程。我们还需要二个条件才能解出所有系数值。获取这二个额外方程的方法之一是在p0和pn处设二阶导数为0。另一个方法是增加二个“隐含”型值点，各自位于型值点序列的两端，即增加p-1和pn+1两个型值点。这时，所有的原有型值点都成了内点，具有所需的4n个边界条件。

尽管三次自然样条曲线是一条三次插值样条曲线，能够做到曲线通过所有型值点，但它有一个主要缺点，即必须解n次方程组，整个曲线受所有型值点控制,不允许“局部控制”，因此不给出完整的新型值点集，不可能去构造曲线的一部分。所以，在实际应用中很少采用三次自然样条曲线。

5.3.2Hermite样条曲线Hermite样条曲线是以法国数学家CharlesHermite命名的，是一个分段三次多项式并在每个型值点有给定的切线。与三次自然样条曲线不同，Hermite样条曲线可以局部调整，因为每个曲线段仅依赖于端点约束。整个曲线通过所有的型值点，而对于每个曲线段来说，它通过两个相邻的型值点。我们把这两个点作为该段曲线的起点和终点，设为P0点和P1点。并且假定曲线段在两端点处的切矢为已知，分别设为

P'0和

P'1。曲线的参变量t是在两个端点取值0和1之间变化，如图5.11所示。

P0(t0=0)P’0P1(t1)P’1图5.11设定端点条件在这样的前提下，对于每一个三次曲线段，就有了四个独立条件，它们是：两个端点的位置向量及曲线在两端点处的切矢。于是，根据这四个条件，就可以求出分段表达式(5-13)中的四个系数。

P0=B1＋B2t＋B3t2十B4t3=B1(当t=0)

P1=

B1＋B2t＋B3t2十B4t3=B1+B2+B3+B4

(当t=1)

P’0=B2＋2B3t十3B4t2=B2(当t=0)(5-14)

P'1=B2＋2B3t十3B4t2=B2+2B3+3B4

(当t=1)上述公式(5-14)可写成如下矩阵形式：(5-15)

=该方程对多项式系数求解，可得Hermite样条曲线的矩阵表达式：(5-16)

P(t)=[t3t2t1]最后可以将上式展开，得到Hermite混合函数的表达式：P(t)=Pk(2t3-3t2+1)+Pk+1(-2t3+3t2)+Pk’(t3-2t2+t)+Pk+1’(t3-t2)(5-17)Hermite多项式对某些数字化应用有用，这里它不太困难就可以得出或估算出曲线的斜率。但对计算机图形学中的大部分问题而言，除了型值点坐标外，更好的做法是不需要输入曲线斜率值或其他几何信息就能生成样条曲线。因此，就出现了Cardinal样条这种基于Hermite样条的变化形式，它不需要输入控制点上的曲线导数值，而是采用控制点的坐标位置来计算导数。

5.3.3Cardinal样条曲线像Hermite样条曲线一样，Cardinal样条曲线也是分段三次插值曲线，且每曲线段终点处均指定切线。与Hermite样条的区别是，不一定要给出终点的切线值。在Cardinal样条曲线中，一个控制点处斜率值可由两个相邻控制点坐标来计算。

一个Cardinal样条曲线完全由四个连续控制点给出。中间两个控制点是曲线段端点，另二个点用来计算终点斜率。如图5.12所示，设P(t)是两控制点Pk和Pk+1间的参数三次函数式，则从Pk-1到Pk+2间的四个控制点用于建立Cardinal样条曲线段的边界条件：

Pk-1PkPk+2Pk+1P(t)图5.12两控制点Pk和Pk+1间的参数三次函数式P(t)P0=Pk

P1=Pk+1

P’0=1/2(1-ts)(Pk+1-Pk-1)(5-18)

P'1=1/2(1-ts)(Pk+2-Pk)则控制点Pk和Pk+1处的斜率分别与弦Pk-1Pk+1和PkPk+2成正比(图5.13)。参数ts称为张力(tension)参数，因为ts控制Cardinal样条曲线与输入控制点间的松紧程度。Pk-1PkPk+2Pk+1Pk’Pk+1’图5.13控制点Pk和Pk+1处的斜率分别与弦Pk-1Pk+1和PkPk+2成正比图5.14说明了张力ts取很小和很大值时Cardinal曲线的形状。当ts=0时，这类曲线称为Catmull-Rom样条曲线或Overhauser样条曲线。

ts<0(较松曲线)ts>0(较紧曲线)图5.14张力参数ts在Cardinal曲线形状中的作用用类似Hermite样条曲线中的方法，我们可以将边界条件(5-17)转换成矩阵形式：

P(t)=[t3t2t1](5-19)这里s=(1-ts)/2。

将矩阵方程(5-17)展开成多项式形式，有：P(t)=Pk-1(-st3+2st2-st)+Pk[(2-s)t3+(s-3)t2+1]+Pk+1[(s-2)t3+(3-2s)t2+st]+Pk+2(st3-st2)(5-20)以下是Cardinal样条曲线的VC程序：voidCardinalSpLine(CPoint*pp,intn,floattension)//这里数组pp的下标从1开始{intx,y,i,j,k=1000;doublet,t1,t2,t3,a,b,c,d,s;t=1.0/k;s=(1.0-tension)/2.0;pp[0].x=2*pp[1].x-pp[2].x;//端点处理pp[0].y=2*pp[1].y-pp[2].y;pp[n].x=2*pp[n-1].x-pp[n-2].x;pp[n].y=2*pp[n-1].y-pp[n-2].y;MoveTo(pp[1]);for(i=1;i<n-1;i++){ for(j=1;j<=k;j++) {t1=j*t;t2=t1*t1;t3=t2*t1;a=-s*t3+2*s*t2-s*t1;b=(2-s)*t3+(s-3)*t2+1; c=(s-2)*t3+(3-2*s)*t2+s*t1; d=s*t3-s*t2; x=int(a*pp[i-1].x+b*pp[i].x+c*pp[i+1].x+d*pp[i+2].x); y=int(a*pp[i-1].y+b*pp[i].y+c*pp[i+1].y+d*pp[i+2].y); LineTo(x,y);}}}5.4Bézier(贝济埃)曲线

在上面讨论的抛物样条曲线和三次插值样条曲线中，它们的共同特点是：最终生成的曲线通过所有给出的型值点，我们称之为“点点通过”。通过以上的讨论，我们可以看到，这两种曲线的计算和编程都比较简单易行，并且分别使曲线达到了一阶和二阶连续．这已经可以满足一般的工程应用问题了，所以应用较广。

但在外形设计中，比如船舶和汽车的外形设计，初始给出的型值点往往并不精确，并且有些地方完全仅仅是为了外观上的考虑。在这样的一种应用前提下，去一点一点地插值计算并使曲线点点通过不很合算。另外，在外形设计中，最好能使设计的结果灵活直观地反映出来，比如局部修改型值点也能在外形上得到直观的反映，以便于设计者进行交互式设计操作。这些，前面所介绍的两种插值样条曲线都是不够的。

为了改善这种情况，法国Bézier(贝济埃)提出了一种新的参数曲线逼近方法,称作贝济埃曲线。后来Gordon,Riesenfeld和Forrest等人对Bezier曲线进行了修改和发展,提出了B样条曲线。这两种曲线都因能较好地适用于外形设计的特殊要求而获得了广泛的应用。

贝济埃曲线的形状是通过一组多边折线(也称为贝济埃多边形或特征多边形)的各项唯一地定义出来的。

在该多边折线的各顶点中，只有第一点和最后一点是在曲线上的，其余的顶点则用来定义曲线的导数、阶次和形状。第一条边和最后一条边则表示出了曲线在起点处和终点处的切线方向，即第一条边和最后一条边分别和曲线在起点和终点处相切。曲线的形状趋向于多边折线的形状。改变多边折线的顶点位置和曲线形状的变化有着直观的联系。图5.15列举了一些贝济埃多边折线和相应的贝济埃曲线的形状关系。

P0P1P2P3P0P1P2P3P0P1P2P3图5.15贝济埃曲线

Bezier曲线的例子5.4.1Bezier曲线的定义贝济埃曲线是由多项式调和函数推导出来的，通常n＋1个顶点定义一个n次多项式，其参数向量表达式为；

在式（5-20）中，Pi为各顶点的位置向量，Bi,n(t)为伯恩斯坦（Bernstein）基函数，也就是贝济埃多边形的各顶点位置向量之间的调和函数。该函数的表达式为：如果我们规定：00和0！均为1，那么，当t=0时：P(0)=P0B0,n(0)+P1B1,n(0)+P2B2,n(0)+…+PnBn,n(0)Cni的求取方法：杨辉三角形的展开式n=111n=2121n=31331n=414641n=515101051…………例：1、2和3次Bézier曲线1次Bézier曲线2次Bézier曲线3322)1(3PtPtt+-+203)1(3)1()(tPtPttP-+-=3次Bézier曲线1在t=0时除第一项外其余各项均为0，所以当t=1时：P(1)=P0B0,n(1)+P1B1,n(1)+P2B2,n(1)+…+PnBn,n(1)在t=1时除最后一项外其余各项均为0，所以

从以上结果可以得出，曲线通过多边折线的起点和终点。下面，我们通过对基函数求导，来分析一下曲线在两端点处的切矢情况。

于是得：（5-23）在起始点，t=0，公式(5-23)中只有i=0,1两项有效，即：

=n（P1一P0）在终止点，t=1，公式(5-23)中只有i=n-1,n两项有效，同理可得：

P’(1)=n(Pn-Pn–1)

这说明：贝济埃曲线在两端点处的切矢方向是与贝济埃多边折线的第一条边和最后一条边相一致的。我们已经知道，确定m个顶点，可以定义一条m–1次的贝济埃曲线。下面，我们来讨论一下较低次的贝济埃曲线，例如二次贝济埃曲线和三次贝济埃曲线。

5.4.2Bezier曲线的性质如图5.4.2所示，对于n+1个控制点Pi（i=0,1,…,n），构成特征多边形的顶点P0P1…Pn构成曲线P(t)，其中t∈[0,1],我们有Bernstein多项式的性质和Bezier曲线的性质如下：1.Bernstein多项式的性质1）非负性Bi,n(t)≥02）权性

3）对称性Bi,n(t)=Bn-i,n(1－t)

图5-18Bezier曲线定义及性质示意图4)递推性Bi,n(t)=(1－t)Bi，n-1(t)＋Bi-1，n-1(t)，因为5）导函数(5-39)1.Bernstein多项式的性质2.Bezier曲线的性质1）端点的位置矢量由式（5-33）、式（5-34）得到：P(0)=P0，P(1)=Pn2)端点的切矢量由式（5-36）、式（5-37）得到：P’(0)=n(P1－P0)，P’(1)=n(Pn－Pn-1)3）端点的二阶导矢量和曲率参考图5-18，Bezier曲线P(t)在P0、P1的二阶导矢量分别为：P”(0)=n(n-1)(P2-2P1-P0)，P”(1)=n(n-1)(Pn-2Pn-1+Pn-2)(5-40)根据式(5-8)计算曲率公式：，Bezier曲线P(t)在P(0)、P(1)处的曲率分别为：，

其中k(t)为P(t)的曲率，PiPj=Pj－Pi，参考图5-18。3）端点的二阶导矢量和曲率图5-18Bezier曲线定义及性质示意图曲率k(0)、k(1)推导过程如下：根据式(5-8)计算曲率公式：由前面Bezier曲线性质，可以得到：，

4）对称性如图5-18所示，对于特征多边形,我们既可以由P0出发构造Bezier曲线，也可以从Pn开始反向构造Bezier曲线。这两条曲线形状完全相同，但参数化方向相反，即Bezier曲线具有对称性。同时表明，同一特征多边形定义的Bezier曲线是唯一的。图5-18Bezier曲线定义及性质示意图7）变差缩减性若Bezier曲线P(t)的特征多边形是一个平面图形，则平面内任意直线与P(t)的交点个数不多于该直线和特征多边形的交点的个数，这一性质叫做变差缩减性。该性质反映了Bezier曲线比其特征多边形的波动要小，也就是说曲线比其特征多边形所在的折线更为光顺。如图5-19所示，直线L与Bezier曲线P(t)的交点为a、b两点，L与特征多边形的交点数为4，符合变差缩减性。5）凸包性Bezier曲线完全被包容在由特征多边形形成的凸包内，如图5-18所示。6）几何不变性Bezier曲线的形状仅取决于特征多边形的顶点，而与坐标系的选取无关图5-19Bezier曲线变差缩减性质示意图5.4.3二次贝济埃曲线当m=3时，顶点P0，P1，P2可定义一条二次(n=2)贝济埃曲线。此时式（5-21）可以改写成：

P(t)=(1–t)2P0＋2t(1–t)P1＋t2P2

(0≤t≤1)(5-24)写成矩阵形式为：

P(t)=[t2t1]在式(5-24)中，相对应于式(5-21)中的调和函数Bi,n(t)分别为：

B0,2(t)=1－2t＋t2

B1,2(t)=2t－t2

B2,2(t)=t2

根据式(5-23)，当n=2时，二次贝济埃曲线在起点P0处有切向量P'0=P'(t=0)=2(P1–P0)；在终点P2处有切向量P'2=P'(t=1)=2(P2–P1)。同时，当t=1/2时：该式说明，二次贝济埃曲线经过△P0P1P2中的一条中线P1Pm的中点P。综上所述，我们可以看出：二次贝济埃曲线是一条抛物线。见图5.20所示。

P1P0P2PmP图5.20二次贝济埃曲线5.4.4三次贝济埃曲线当m=4时，顶点P0,P1,P2,P3四点可定义一条三次(n=3)贝济埃曲线。此时式(5-30)可以改写为：

P(t)=(1–t)3P0+3t(1–t)2P1+3t2(1–t)P2+t3P3

=(1–3t+3t2-t)3P0+(3t–6t2+3t3)P1+(3t2–3t3)P2+t3P3

(0≤t≤1)(5-44)同样，该式可改写成矩阵表达式为：P(t)=[t3t2t1]在式(5-44)中，相对应于式(5-21)中的调和函数Bi,n(t)分别为：

B0,3(t)=1–3t＋3t2–t3

B1,3(t)=3t–6t2＋3t3

B2,3(t)=3t2–3t3

B3,3(t)=t3

三次贝济埃曲线是二阶连续的。∵P’(t)=3(1–t)2P0+(3–12t+9t2)P1+(6t–9t2)P2+3t2P3P”(t)=6(1–t)P0+(–12+18t)P1+(6–18t)P2+6tP3P”i(1)=6(1–1)P0+(–12+18)P1+(6–18)P2+6P3=6(P1+P3-2P2)P”i+1(0)=61–0)P1+(-12+0)P2+(6–0)P3+0P4=6(P1+P3-2P2)∴P”i(1)=P”i+1(0)对于更高次的贝济埃曲线，在此不再一一列举并讨论。

图5-21和图5-22为三次Bezier曲线的各种不同情况的分析比较。图5-21控制点相同，顺序不同的三次Bezier曲线图5-22移动控制点P2的不同效果比较5.4.5反算Bezier曲线控制点

若给定n+1个型值点Qi(i=0,1,…,n)，要求构造一条Bezier曲线通过这些点。问题是如何求得过Qi的Bezier曲线的控制点Pi(i=0,1,…,n)。通常可取参数t=i/n与点Qi相对应，用以反算Pi。设Qi在曲线C(t)上，且有由此式我们可以得到下面关于Pi(i=0,1,…,n)的n+1个方程组成的线性方程组:

由这组方程可解出，这就是过Qi的Bezier曲线的特征多边形的顶点。5.4.6Bezier曲线的多项式定义生成算法根据式(5-30)，可以编写出绘制任意阶次的Bezier曲线的C语言程序。doublepowi(doublev,intk)//计算vk{doubletemp=1.0;if(k==0||v==0)return1;//00=1else{for(inti=1;i<=k;i++)temp=temp*v;}returntemp;}longfac(intm)//计算m!{inti;longtemp=1;if(m==0)return1;//0!=1else{for(i=2;i<=m;i++)temp=temp*i;}returntemp;}voidBezier(POINT*p,intn){intx,y,i,j,k=100;doublet=1.0/k,t1,u,v;doubletemp,temp1,temp2,bi;moveto(p[0].x,p[0].y);for(j=1;j<k;j++){t1=j*t;u=t1;v=1-u;x=0;y=0;for(i=0;i<=n;i++){temp=double(fac(n)/fac(i)/fac(n–i));

temp1=powi(u,i);temp2=powi(v,n–i);bi=temp*temp1*temp2;

x=x+bi*p[i].x;y=y+bi*p[i].y;}lineto(x,y);}lineto(p[n].x,p[n].y);}5.4.7Bezier曲线的分割递推deCasteljau算法图5-23deCasteljau算法分割递推过程示意图图5-24分割点Pi,r的递推关系图如图5-23结合图5-24所示，对于控制节点数为4，边数为3的控制多边形P0P1P2P3,，依次分别在边P0P1、P1P2、P2P3上找到一点P0,1、P1,1、P2,1，该点分别将所在的边分成2：1，然后，将上述3个分割点构成新的控制多边形P0,1P1,1P2,1,其节点数为3，边数为2；再以同样的方法以及同样的比例对边P0,1P1,1和边P1,1P2,1进行分割得到分割点P0,2与P1,2，最后对边P0,2P1,2进行与前面相同比例的分割，得到点P0,3。P0,3即为由原控制多边形所确定的3次Bezier曲线上的参数为t的点P(t)。如果我们让t在[0,1]间变动，并且让变动Δt取一个较小的增量（例如：Δt=0.1），循环多次就可作出3次Bezier曲线。综合上述我们可以总结出Bezier曲线生成的deCasteljau算法：对于任意控制多边形，在以PiPi+1为端点的第i条边上，找一点Pi,1(t)，把该点分成t：(1-t)比例，则分割点Pi,1(t)=(1-t)Pi+tPi+1,i=0,1,…,n-1,这n个点组成一个新的n-1边形，对该多边形重复上述操作，得到一个新的n-2变形的顶点Pi,2(t)(i=0,1,…n-2)依次类推，连续作n次后，得到一个单点Pi,n(t),该点就是Bezier曲线P(t)上参数为t的点P(t),让t在[0,1]间变动，就得到Bezier曲线。用数学归纳法不难证明，按照式Pi,1(t)=(1-t)Pi+tPi+1，递推n次，最终结果P0,n(t)就是满足式(5-30)的Bezier曲线上的点P(t)，即deCasteljau算法最后的切割点可以表示为：

(5.4.17)图5-25deCasteljau算法和多项式算法绘制Bezier曲线对比5.4.8Bezier曲线的拼接a)G0、C0连续的拼接b)G1连续的拼接图5-26两Bezier曲线的拼接示意图Bezier曲线曲面的拼接技术，在实际复杂曲