第3章 离散点绘制平面曲线

第三章离散点绘制平面曲线§1概述§2贝塞尔(Bezier)曲线§3B样条曲线§4抛物线调配曲线§5三次参数样条曲线

本章小结一、规则曲线与不规则曲线平面曲线一般分规则和不规则曲线两类。

规则曲线是指可以用一个方程或描述的曲线，即曲线的方程已知。绘制方法：以足够小的步长取曲线上足够多的点，然后利用曲线方程求出这些点的坐标，最后用直线连接相邻的点即可绘出曲线。曲线的绘制精度取决于曲线上取点的密度，密度越大，精度就越高，曲线就越光滑。§1概述不规则曲线（拟合曲线）：指已知平面一些离散点的坐标，但曲线方程未知，需要人为设计曲线方程对这些点进行拟合形成的曲线。对于这类曲线的绘制，首先要找出一种合理的拟合方法来设计曲线方程。拟合方法包括：贝塞尔曲线法、B样条曲线法、抛物线调配曲线法、三次参数样条曲线法，最小二乘法等。拟合方法不同→曲线拟合方程不同→绘制的曲线形状也不同。但是，一旦拟合方法确定并得到相应的曲线拟合方程，不规则曲线也就变成了规则曲线。本章的内容就是介绍：如何根据离散点的坐标，利用拟合方法建立曲线拟合方程，绘制不规则曲线。在用拟合方法建立曲线拟合方程时，通常把不规则曲线分为两类：1.“点点通过”式当已知离散点的位置较精确时，拟合的曲线通过所有的已知点。曲线方程确定后→规则曲线→以足够小的步长获取相邻离散点之间若干个数据点（插值点）的坐标，并用直线连接它们。当步长极小时，所绘出的直线连线在视觉上便是一条光滑的曲线。求相邻离散点之间若干数据点的问题称为插值问题。二、不规则曲线（拟合曲线）的分类2.“平均通过”式当已知数据点有一定误差时，所拟合的曲线不通过所有已知点，曲线代表的是这些数据点的变化趋势。1000-2000-3000-4000-

102030402.02.22.4深度/m总孔隙度（%）密度

要求：设计的曲线方程与所有已知点的“距离”总和最小。即：曲线方程是对所有已知点的“逼近”。典型例子：地层孔隙度－深度变化曲线。

“点点通过”式也称为插值曲线，“平均通过”式也称为逼近曲线。一般可分为直角坐标方程和参数方程两类。例如：对于圆心坐标为(x0,y0)，半径为r的圆，其直角坐标方程为：(x-x0)2+(y-y0)2=r2参数方程为：在计算机绘图时，使用参数方程要比直角坐标方程方便。

参数方程还有另外一种形式。三、曲线的方程分类

x=x0+rcos(t)

y=y0+rsin(t)

（0≤t≤2π）例如：对于二次抛物线曲线，其参数方程可表示为设上述参数方程可写为：方程②称为曲线的参数矢量方程。

在构建曲线方程时，通常采用方程②的矢量形式；在绘制曲线时，通常采用方程①的分量形式。①②当平面上已知数据点较多时，针对所有点拟合曲线方程有时非常困难，或者得到的曲线方程非常复杂、不实用。可分别针对部分点进行分段拟合，分段拟合得到的曲线段会涉及到曲线段的连接问题。例如：若已知3个平面离散点，用1条抛物线段对其拟合；若有5个点，可用2条抛物线段进行分段拟合。两条曲线段在连接点P3处，并非光滑，需对该点进行光滑处理。在光滑处理时，达到什么标准为“光滑”呢？给出2个一般标准：四、曲线段的光滑连接1.C

1连续在连接点pj

处，若两曲线段的切线斜率相等（相同的切线），或者说一阶导数连续，即称两曲线段在连接点pj

处的光滑连接达到C1连续。2.C2连续在连接点

pj

处，不仅两曲线段的切线斜率相同，而且切线斜率的变化率也相同，即称两曲线段在连接点pj

处的光滑连接达到C2连续。显然C2连续比C

1连续要求更高，曲线的连接更光滑。另外还有更高的连续标准，但对一般绘图，曲线段的连接满足C

1或C

2连续，其光滑已足够。。pj及可以设计一个光滑的曲线段去逼近这个特征多边形。

贝塞尔曲线就属于这类曲线。该曲线由法国汽车工程师Bezier首先提出，最初用于汽车零件外形的设计。目前广泛应用于与计算机绘图相关的各个领域。§2贝塞尔曲线一、Bezier曲线1.特征多边形特征多边形是用直线段依次连接平面上离散点所形成的折线多边形。它反映了所要设计曲线的大致轮廓。P0P1P2P4P3P52.Bezier曲线的参数方程已知三个平面离散点P0、P1、P2，那么由这三点可以定义二次抛物线段，其参数矢量方程可表示为：其中P0P1P2经过推导（详见教材P58），可求出A0、A1、A2：

⑴

t=0时，抛物线经过P0，并与P0P1相切。

⑵

t=1时，抛物线经过P2，并与P1P2相切。设上述抛物曲线满足下列条件：因此，该抛物线段的参数矢量方程为：或用矩阵形式表示：该方程描述的曲线称为二次贝塞尔曲线。其分量形式为：

绘制方法：将参数t的区间[0,1]划分为n等份，依次取t=1/n,2/n,3/n,…，利用曲线参数方程计算对应的各点坐标，并用直线段依次连接各点。在已知三个平面离散点时，可绘制二次Bezier曲线；若有n+1个平面离散点时，则可绘制n次Bezier曲线。其中

n次Bezier曲线参数矢量方程的一般形式：其中：Pi（i=0,1,…,n）为n+1个平面离散点，;，其分量形式为：实际应用中，对多个离散点（如：7个离散点）绘制Bezier曲线时：可以用高次Bezier曲线进行绘制（如：7个离散点使用六次Bezier

曲线）；可以用低次Bezier曲线进行分段拟合（常使用二次和三次Bezier

曲线）；需要对连接点进行光滑处理（C

1连续、C

2连续），已有计算机图形学方面的论文研究了此问题。地质绘图中常采用低次Bezier

曲线的分段拟合。二、二次Bezier曲线的程序设计

已知：平面上3个离散点P0、P1、P2的坐标(x0,y0)、(x1,y1)、(x2,y2)。二次Bezier曲线参数方程的分量形式：其中

编程思路：将参数t的区间[0,1]划分为n等份，依次取t=1/n,2/n,3/n,…，利用方程计算对应各点的坐标，并用直线段依次连接各点。1.二次Bezier曲线子程序subroutinebezier2

(x0,y0

,x1,y1,x2,y2,n)

dt=1.0/n

a0=x0

a1=2*(x1-x0)

a2=x0-2*x1+x2n等份（将参数t的区间[0,1]n等份）b0=y0

b1=2*(y1-y0)

b2=y0-2*y1+y2

callmovea(x0,y0)

do100k=1,n

t=k*dt

t2=t*t

x=a0+a1*t+a2*t2y=b0+b1*t+b2*t2calllinea(x,y)

100continueend绘制Bezier曲线2.主程序：

callin

callfact(10.0)

write(*,*)'请输入3个点的x、y坐标：'

read(*,*)x0,y0,x1,y1,x2,y2

callbezier2(x0,y0,x1,y1,x2,y2,100)

end①曲线通过始点和终点，并与特征多边形首末两边相切于始点和终点，中间点将曲线拉向自己。②平面离散点控制曲线的形状，改变一个离散点的坐标，曲线的形状将随之改变（点对曲线具有整体控制性）。③曲线落在特征多边形的凸包之内，它比特征多边形更趋于光滑。三、Bezier曲线的特点P0P1P2④贝塞尔曲线属于“平均通过”式曲线。§3B样条曲线B样条曲线是在Bezier曲线基础上发展起来的一类曲线，它克服了Bezier曲线整体控制性所带来的不便，最常用的是二次和三次B样条曲线。一、二次B样条曲线

1、二次B样条曲线的特点①起点为P0、P1点的中点，并与线段P0P1相切；②终点为P1、P2点的中点，并与线段P1P2相切；P0P1P2③除起点、终点外，中间点将曲线拉向自己。④二次B样条曲线是“平均通过式”曲线。二次Bezier曲线参数方程：2、与贝赛尔曲线有什么异同？P0P1P2其中P0P1P2其中二次B样条曲线参数方程的分量形式：3.参数方程已知三个平面离散点P0、P1、P2，由这三点可以定义二次抛物线段，其参数矢量方程形式为：

二次B样条曲线参数方程的矩阵形式：4.多点时二次B样条曲线的应用二次B样条曲线是针对3个离散点设计的，如果多于3个离散点，则采用二次B样条曲线进行分段拟合。假设有n+1个离散点，记为Pi（i=0,1,…,n），采用二次B样条曲线进行分段拟合的过程是：以P0、P1、P2为控制点绘制第1条二次B样条曲线；以P1、P2、P3

…………2………；┇…┇……；以Pn-2、Pn-1、Pn

……n-1………；6个离散点绘制二次B样条曲线的情况如下图：P4P3P2P1P0P5

问题：如何使曲线始于P0而止于P5？

由于二次B样条曲线在连接点处具有相同的切线，故多点的二次B样条曲线具有C1级连续。5.二次B样条曲线的边界处理（边界问题）思路：增加两点Ps、Pe；第1条曲线用Ps代替P0，最后一条曲线时用Pe代替最后一个点。

按照上述方法处理后，可绘制始于起点而止于终点的二次B样条曲线。5.多点二次B样条曲线程序设计设计思想：在、的延长线上，且在、的延长线上，且PeP0Pn边界处理示意图……PsP1Pn-1

①参数存放点坐标的数组；点的个数；[0,1]区间等分数。②边界处理③计算方程系数

④绘制B样条曲线

移动到第一个点边界处理段内按步长画线计算第K段分量式参数方程的系数nt步完否?n-2段完否?NNYY

subroutinebspline2(x,y,n,nt)real

x(n),y(n)

dt=1.0/nt

callmovea(x(1),y(1))x(1)=2*

x(1)-x(2)y(1)=2*

y(1)-y(2)x(n)=2*

x(n)-x(n-1)y(n)=2*

y(n)-y(n-1)边界处理do10k=1,n-2a0=0.5*(x(k)+x(k+1))a1=x(k+1)-x(k)a2=0.5*(x(k)-2.0*x(k+1)+x(k+2))b0=0.5*(y(k)+y(k+1))b1=y(k+1)-y(k)b2=0.5*(y(k)-2.0*

y(k+1)+y(k+2))do20j=1,ntt=j*dtt2=t*txx=a0+a1*

t+a2*t2yy=b0+b1*

t+b2*t2calllinea(xx,yy)20continue10continueend

参数方程计算式

主程序：

realx(1000),y(1000)

character*30filenamewrite(*,*)'请输入数据文件名：'read(*,*)filenameopen(1,file=filename)read(1,*)ndo10i=1,nread(1,*)x(i),y(i)10continue

callin

callfact(10.0)callbspline2(x,y,n,100)end假设有5个已知离散点：（2,3）、（5,6）、（7,1）、（8,4）、（9,2），试绘制二次B样条曲线。二次贝塞尔曲线由3个平面离散点确定曲线特点多点时可进行分段拟合，连接点一般不光滑二次B样条曲线由3个平面离散点确定曲线特点多点时可进行分段拟合，连接点为C1级光滑边界问题

二、三次B样条曲线三个平面离散点确定一条二次B样条曲线，四个平面离散点则确定一条三次B样条曲线。1.三次B样条曲线的参数表达式矩阵形式：上式中分量形式：设中点为：

①曲线起点S位于的中线

上，距

点处；曲线终点E位于的中线

上,距

点处；②曲线起点切线平行于

，终点切线平行于

；P0P1P2P3M1M2始、终点切线设四个离散点为P0、P1、P2、P3；③为“平均通过式”曲线。2.三次B样条曲线的特点3.多点三次B样条曲线的分段拟合设已知n个平面离散点，记为Pi（i=1,2,…,n）。以P1、P2、P3、P4绘制第1条三次B样条曲线；以P2、P3、P4、P5绘制第2条三次B样条曲线；┇┇┇以Pn-3、Pn-2、Pn-1、Pn绘制第n-3条三次B样条曲线。用上述方法绘出的各段曲线自然衔接，具C

2级连续。P0P1P2P3M1M2PsPe

4.三次B样条曲线的边界处理要使曲线起于始点，止于终点，必须对进行边界处理。思路：与二次B样条曲线类似，在两端各增加一个点；以新的点序列Ps、P0、P1、P2、P3、Pe为控制点分别绘制三次B样条曲线。5.多点三次B样条曲线程序设计自学。§4抛物线调配曲线Bezier曲线和B样条曲线并不通过平面上所有已知的离散点，也就是说，它们是一种“平均通过”式曲线。

如果要求曲线通过所有的已知点，上述两种方法就不能使用，为此再介绍抛物线调配曲线，它是一条“点点通过”式曲线。一、抛物线调配曲线的绘图原理设已知n个平面离散点Pi（i=1,2,…,n），利用二次抛物线进行分段拟合:以P1、P2、P3为控制点绘制第一条过三点的抛物线段；以P2、P3、P4为控制点绘制第二条过三点的抛物线段；┇┇┇最后，以Pn-2、Pn-1、Pn为控制点绘制第n-2条过三点的抛物线段。每相邻两条抛物线段的相交部分会出现2条曲线段。Pi+1Pi+2Pi+3Pi合并

调配：对相交部分的2条曲线段进行合并（调配）处理。抛物线调配曲线需要解决的问题：①过三点的二次抛物线参数方程；②相邻抛物线相交部分2条曲线段的调配（合并）方法；③合并后连接点的光滑问题。Pi+1Pi+2Pi+3Pi合并二、过三点的抛物线参数方程已知三个平面离散点P0、P1、P2，那么由这三点定义的二次抛物线参数矢量方程为：条件为：⑴t=0时，抛物线经过P0点；⑵t=½时，抛物线经过P1点，且切线矢量=P2-P0；⑶t=1时，抛物线经过P2点。二次抛物线参数方程的矩阵形式：展开：

三、调配曲线的计算公式设已知n个平面离散点Pi（i=1,2,…,n），依次取P1、P2、P3、…、Pn-2作为起点绘制过三点的抛物线，共绘出n-2条抛物线段。Pj+2

Pj+3

Pj+1

PjSj+1(tj+1)Sj(tj)Pj+2

Pj+3

Pj+1

PjSj+1(tj+1)Sj(tj)第j条抛物线参数方程为：第j+1条抛物线的参数方程为：

合并(调配)曲线取为第j条曲线与第j+1条曲线的加权平均。设权函数为：则调配(合并)曲线（起点为Pj+1、终点为Pj+2）的方程为：方程右边三个参数：T、tj、tj+1的取值范围为：Pj+2

Pj+3

Pj+1

PjSj+1(tj+1)Sj(tj)Pj+2

Pj+3

Pj+1

PjSj+1(tj+1)Sj(tj)三个参数T、tj、tj+1的取值范围：①T的范围：0≤T≤1

②tj的范围：调配曲线合并的是第j条曲线的后半段，故0.5≤tj≤1

③tj+1的范围：调配曲线合并的是第j+1条曲线的前半段，故0≤tj+1≤0.5设参数t的范围为0≤t≤0.5，有：调配曲线方程变为：再将前面Sj和Sj+1的方程代入，整理后得到：上式称为抛物线调配曲线的参数方程。当j=1时，调配曲线位于P2、P3之间；当j=2时，调配曲线位于P3、P4之间；……当j=n-3时，调配曲线位于Pn-2、Pn-1之间。因此，n个点可得到n-3条调配曲线，首尾两条没有调配。…如何利用调配方程对首尾两条进行调配？分别在首尾补上两点P0、Pn+1。依据补点后新的点序列P0、P1、…、Pn、Pn+1，可绘制n-1调配曲线。……如何确定P0、Pn+1的坐标？对于P0点，满足：对于Pn+1点，满足：常设：因此：P0=P2，Pn+1=Pn-1同时，绘出的n-1条调配曲线在连接点处具有C1级连续（证明略）。四、抛物线调配曲线子程序设计

subroutineparspl(x,y,n,nt)realx(n),y(n),p1(102),p2(102)

do10i=1,np1(i+1)=x(i)10p2(i+1)=y(i)p1(1)=x(2)p2(1)=y(2)p1(n+2)=x(n-1)p2(n+2)=y(n-1)

dt=0.5/ntcallmovea(x(1),y(1))补上首末2点的坐标x、y是存放点坐标的数组；n是点个数；nt是绘制每条调配曲线的等分数 do30j

=1,

n-1 do20i

=1,nt t=i*dt f1=4.0*t*t-t-4.0*t*t*t f2=1.0-10.0*t*t+12.0*t*t*t f3=t+8.0*t*t-12.0*t*t*t f4=4.0*t*t*t-2.0*t*t xx=f1*p1(j)+f2*p1(j+1)+f3*p1(j+2)+f4*p1(j+3) yy=f1*p2(j)+f2*p2(j+1)+f3*p2(j+2)+f4*p2(j+3)calllinea(xx,yy)20continue30continueend主程序：realx(100),y(100)

character*30filenamewrite(*,*)'请输入数据文件名：'read(*,*)filenameopen(1,file=filename)read(1,*)ndo10i=1,n10read(1,*)x(i),y(i)callincallfact(10.0)callparspl(x,y,n,50)end以教材p78表3-3数据点绘制正弦曲线。来源于早期的手工绘图：绘图员把一根富有弹性的细木条放在数据点的位置上，并施以力量，从而使木条弯曲通过每一个数据点，木条变形后的形态就是三次参数样条曲线。§5三次参数样条曲线由弹性力学理论可以证明：曲线近似为分段三次多项式，且在连接点处具有C2连续。绘制三次参数样条曲线的基本思路：首先，构建2点间曲线的三次多项式参数方程；然后，进行分段拟合，依次绘制相邻两点间的曲线段；最后，对各曲线段进行光滑连接处理（C

2连续）。

主要问题：构建曲线的三次多项式参数方程。一、过两点的三次多项式参数方程设已知两点P1、P2，则过两点的曲线三次参数矢量方程为：t1～P1、P2点之间距离（弦长）t～P1、P点之间距离（弦长）（曲线上的任一点）条件为：⑴t=0时通过P1点，且曲线段在该点的切矢量为P1′；⑵t=t1时通过P2点，且曲线段在该点的切矢量为P2′。P1′P2′条件为：⑴t=0时通过P1点，且曲线段在该点的切矢量为P1′；⑵t=t1时通过P2点，且曲线段在该点的切矢量为P2′。当t=0时，有

①

②当t=t1时，有

③

④解出B0、B1、B2、B3：将系数代入方程中整理，可得到：上式即为过两点的三次多项式参数方程的矢量形式。其分量形式为：若有n个平面离散点Pi（i=1,2,…,n），进行分段拟合，每相邻两点之间使用一个三次多项式参数方程进行曲线绘制，共有n-1个曲线段。一般情况下，对于其中由Pi、Pi+1两点构成的第i条曲线段的三次多项式参数方程可表示为：

在这n-1个曲线方程中，涉及到P1′、P2′、…、Pn′共n个切矢量。二、n个切矢量的求解

求解思路：利用相邻两条曲线段之间的连接具有C2级光滑。设第i-1条曲线与第i条曲线相连：两条曲线在连接点Pi具有C2级连续，即将两条曲线的参数方程代入，整理后可得到关系式：Pi-1PiPi+1设第i-1条曲线与第i条曲线相连，可得到关系式：

n个点有n-1条曲线，每两条曲线可以得到一个关系方程，共有n-2个方程。需要对P1′、P2′、…、Pn′共n个切矢量求解，尚缺少两个方程。补充首尾两点P1、Pn的边界条件。…补充首尾两点P1、Pn的边界条件有三种方法：⑴固定端边界条件

直接给定两端的切矢量P1′、Pn′。⑵抛物端边界条件两个端点的三阶导数为零，即⑶自由端边界条件两个端点的二阶导数为零，即由上页的n-2个方程和2个边界条件构成一个方程组，称之为切矢量方程组，对其求解可得到n个切矢量P1′、P2′、…、Pn′的值。切矢量方程组的求解可采用一般线性方程组的求解方法。教材采用追赶法，抛物线边界条件，详细求解过程略。最后可得到一个递推公式：初始条件：，，，由Pi、Pi+1两点构成的三次多项式参数方程：

绘制思路：将参数t的区间[0,ti]划分为k等分，每一等分dt=ti/k，依次取t=dt,2dt,3dt,…j*dt，…，利用方程计算对应各点的坐标，并用直线段依次连接各点。对于n个点Pi（i=1,2,…,n），进行分段拟合，每相邻两点之间使用一个参数方程进行绘制。三、三次参数样条曲线子程序设计以直线连接所有等分点把全部等分点坐标存入数组xx、yy计算n个切矢量Pi（i=1,2,…,n）计算每相邻两点之间的弦长ti求相邻两点间的k个等分点坐标n-1条曲线subroutinespline3(x,y,n,k)dimensionx(50),y(50),t(50),a(50),c(50),d(2,50),

*q(250),r(2,50),v(2,50),xx(1000),yy(1000)do10i=1,n-1t(i)=sqrt((x(i+1)-x(i))**2+(y(i+1)-y(i))**2)10continuedo20i=2,n-1a(i)=t(i)/(t(i-1)+t(i))c(i)=1.0-a(i)d(1,i)=3.*a(i)*(x(i)-x(i-1))/t(i-1)+3.*c(i)*(x(i+1)-x(i))/t(i)d(2,i)=3.*a(i)*(y(i)-y(i-1))/t(i-1)+3.*c(i)*(y(i+1)-y(i))/t(i)20continuex、y点坐标数组，n是点个数k是绘制每条曲线的等分数， q(1)=1. q(n)=1. r(1,1)=2.*(x(2)-x(1))/t(1) r(2,1)=2.*(y(2)-y(1))/t(1) r(1,n)=2.*(x(n)-x(n-1))/t(n-1) r(2,n)=2.*(y(n)-y(n-1))/t(n-1) do70i=2,n-1 w=2.0-a(i)*q(i-1) q(i)=c(i)/w do75j=1,2 r(j,i)=(d(j,i)-a(i)*r(j,i-1))/w75 continue70 continue v(1,n)=(r(1,n)-r(1,n-1))/(q(n)-q(n-1