136.参数曲线的快速生成算法----毕业设计.doc

毕业论文：参数曲线的快速生成算法 江 南 大 学毕 业 设 计 论 文论文题目：参数曲线的快速生成算法姓 名： 学 院：信息工程学院专 业：计算机科学与技术指导老师： 日 期：2003年6月摘 要本毕业设计主要研究参数曲线的直接快速生成，要直接生成参数曲线就需对参数方程x=f(t),y=g(t),(0t1)的参数t每次增加一个步长，然后计算该点的x和y坐标值并绘制该点。要逐点地生成参数曲线，就要求参数t每次增加的步长要使曲线前进的幅度不得超过一个象素长度，否则有可能跨过一个中间象素而产生断点。为了提高曲线生成算法的速度，本毕业设计针对如何选择最佳的步长进行比较讨论，以使曲线前进的幅度在不超过一个象素的前提下，选择尽量大的步长。为了进一步提高算法的速度，在前面讨论的最佳步长的基础上又采用了双步逐点曲线生成算法，即将上述得到的步长增加一倍，以使算法的循环次数减少一半。由于步长增加一倍，这样当曲线前进一步时，其幅度有时会大于一个象素的长度，这时我们通过插值的方法来确定跨过的那个中间象素。通过上述讨论的算法能够比较快速的逐点生成曲线，为了实现上述算法，本毕业设计使用visual c+6.0为工具并以三次bezier曲线、普通参数曲线x=f(t)=x3t3+x2t2+x1t+x0, y=g(t)=y3t3+y2t2+y1t+y0，以及导师所给的一个特殊的曲线方程为例编程实现上述算法。关键词：参数曲线，逐点，双步，visual c+6. 0 作者： 二零零三年 六月abstract this graduation project main reseach the direct born of the parameter curve x=f(t),y=g(t),0= t loadstandcursor(idc_cross)初始化该变量，该语句的作用是取得windows标准鼠标形状句柄并赋给m_hcross，然后就可以通过classwizard添加鼠标动作的消息处理函数。例如，要为程序添加鼠标左键单击消息处理函数，首先在打开的classwizard对话框中的class name组合框中选择view类，然后在object ids列表框选中第一行，并在message列表框中选中wm_lbuttondown一行，最后单击add function按钮即可。类似的还可以为程序添加wm_rbuttondown（鼠标右键单击）消息处理函数。因为需要逐点地生成曲线，因此在程序中需要使用到mfc（microsoft fundation class）中的cdc（设备上下文）类的成员函数setpixel（）来实现在屏幕上画点，cdc类主要用于在指定设备上下文上（如窗口客户区、打印机）进行绘图、显示文本等操作。第二章 计算机图形学中常用的算法2.1 常用直线的算法 画直线的算法有很多，例如数值微分法，中点画线法，bresenham化线算法等。在这里只介绍一个比较方便常用的中点画线法。为了讨论方便，本小节假定直线斜率在0、1之间。其他情况可参照下述讨论进行处理。如图所示，若直线在x方向增加一个单位，则在y方向上的增量只能在0、1之间。假设x坐标为xp的各象素点中，与直线最近者已确定，为(xp,yp)，用实心小圆表示。那么，下一个与直线最近的象素只能是正右方的p1(xp+1,yp)或右上方的p2(xp+1,yp+1)两者之一，用空心小圆表示。再以m 表示p1、p2的中点，即m=(xp+1,yp+0.5）。又设想q是理想直线与垂直直线x=xp+1的交点。显然，若m在q的下方，则p2离直线近，应取为下一个象素；否则应取p1。这就是中点画线法的基本原理。 p2qmp1p=（xp，yp）中点画线算法每步迭代涉及的象素和中点示意图下面来讨论上述算法的实现。假设直线的起点和终点分别是(x0,y0)，(x1,y1)。则直线的方程为 f(x,y) = ax + by + c=0其中，a = y0-y1，b = x1-x0，c = x0yi-x1y0。对于直线上的点，f(x,y)=0；对于直线上方的点，f(x,y)0；而对于直线下方的点f(x,y)0。因此，欲判断前述q在m的上方还是下方，只要把m代入f(x,y)，并判断它的符号。构成判别式 d=f(m)=f(xp+1,yp+0.5)=a(xp+1) + b(yp+0.5) + c当d0时，则应取正右方的p1。当d=0时，二者一样合适，可以随便取一个。我们约定取正右方的p1。对每一个象素计算判别式d，根据它的符号确定下一个象素。至此可以写出完整的算法。但是注意到d是xp和yp的线形函数，可采用增量计算，提高运算效率。在d0的情况下，取正右方象素p1，欲判断再下一个象素应取哪个，应计算 d1=f(xp+2，yp+0.5) =a(xp+2) + b(yp+0.5)+ c= d + a故d 的增量为a。而若d0，则右上方的象素p2，欲判断再下一个象素，则要计算d2=f(xp+2，yp+1.5)=a(xp+2)+b(yp+1.5)+ c = d + a + b故在第二种情况，d 的增量为a+b。再看d的初值。显然，第一个象素应取左端点(x0,y0)，相应的判别式为 d0=f(x0+1,y0+0.5)=a(x0+1)+ b(y0+0.5) + c = ax0+by0+a+c+0.5b = f(x0，y0)+ a+0.5b但由于(x0,y0)在直线上，故f(x0,y0)=0。因此，d的初值为d0= a+0.5b。 由于使用的只是d 的符号，而且d的增量都是整数，只是其初值包含小数。因此我们可以使用2d代替d，来摆脱小数，写出仅包含整数的算法：midpointline(x0,y0,x1,y1,color)int x0,y0,x1,y1,color; int a,b,delta1,delta2,d,x,y; a= x0- y1; b=x1- x0; delta1=2*a; delta2=2*(a+b); x= x0; y=y0 ; drawpixel(x,y,color); while(x x1) if(d0) x+;y+;d+= delta2;elsex+;d+= delta1;drawpixel(x,y,color); 2.2 常用的参数曲线的绘制算法常用的参数曲线有许多，如bezier，b样条，非均匀有理b样条、圆锥曲线等。本小节将以bezier曲线为代表，讲解其绘制算法。bezier曲线的公式 ：给定空间n+1个点的位置矢量pi，bezier曲线上各点坐标的插值公式是：c(t) =, 0t1pi构成该曲线的特征多边形，是bernstein基函数，也是曲线上各点位置矢量的调和函数。=ti(1-t)n-i=cti(1-t)n-i i=0,1,2,n首先，讨论三次bezie曲线的绘制由上述公式可以推出三次bezier曲线形式：当n=3时,c(t)= = p0(1t)3 + 3p1t (1t)2 + 3p2 t2 (1t) + p3t3 0t1可以将其改写为如下的形式：c3(t) = ( cs p0+ ct p1)s+ct2 p2 )s+ ct3p3其中，s=1-t，则可以通过如下程序绘制曲线：float hornbez(degree,coeff,t)/* input:degree:曲线的次数coeff:保存曲线的系数t: 参数的值 output: 该参数的坐标值*/int degree;float coeff;float t; int i,n_choose_i; float fact,t1,aux; t1=1.0-t; fac=1.0; n_choose_i=1; aux=coeff0*t1; /*开始循环计算坐标值*/ for(i=1idegree;i+) fact = fact*t;n_choose_i= n_choose_i*(degree-i+1)/i;aux=(aux+fact* n_choose_i*coeffi*t1; aux=aux+fact*t*coeffdegree; return aux;下面来讨论n次bezier曲线的绘制。在前面我们介绍了一个程序用于计算相应曲线上的诸点，但其只适用于三次bezie曲线，既不通用而且计算量较大。用如下cas-teljau算法产生曲线上的点列相对要简单许多。cas-teljau算法：给定空间n+1个点pi(i=0,1,2,n)及参数t，则有：(t)=(1-t)(t)+ t(t)r=1,2.n, i=0,1.n-r, 0t1其中(t)即为pi, (t)是曲线上具有参数t的点。 当n=3时，用cas-teljau算法递推出的(t)呈直角三角形，对应如下所示： 该三角形垂直边上的点p0，p1，p2， p3是曲线p(t)在0t1内的控制点，斜边上的点是p(t)在0t1/2内的控制点，水平直角边上的点是p(t)在1/2t1内的控制点。这种用分割bezier曲线控制多边形的方法为离散化的曲线提供了方便。用cas-teljau算法绘制bezier曲线的程序如下：void bez_to_points(degree,npoints,coeff,points)/*input :degree:曲线的次数npoint:控制点的数目coeff:控制点的坐标output: 曲线上点的坐标*/int degree,npointes; float coeff,points; float t,delt; int i; float decas(); delt=1.0/(float)npoints; /*步长*/ t=0.0; for(i=0;i=npoints;i+) pointsi=decas(degree,coeff,t);t=t+delt; float decas(degree,coeff,t)float coeff;float t;int degree; int r,i; float t1; float coeffa10; t1=1.0-t; for(i=1;i=degree;i+) coeffai=coeffi; /*计算(t)地值 for(r=1;r=degree;r+) for(i=1;r=degree-r;i+) coeffai=t1*coeffai+t*coeffat+1; return(coeffa0); 第三章 参数曲线的快速生成算法及实现3.1 普通参数曲线生成算法的介绍本毕业设计所研究的参数曲线生成算法要具备直接、逐点、快速的特点。要直接生成参数曲线就需要对曲线的参数每次增加一个步长，然后计算该点的坐标值x和y并绘制该点。要逐点地生成参数曲线，就要求参数每次增加的步长要使曲线前进的幅度不得超过一个象素的长度，否则有可能跨过一个中间象素而产生断点。要快速生成曲线则要求在不产生断点的情况下尽量减少画点次数以及不重复画一个象素点来达到加快曲线的生成速度。根据参考文献可以得到一个满足上述要求的逐点生成参数曲线的整数算法。设曲线的参数方程为：首先讨论第一个方程x=f( t ),首先把曲线分成n 段（限定t取i/n，其中0in。）根据拉格朗日中值定理：f(x+x)f( x )=f()(x)，当n 满足n时，可以得到如下公式：=1 （3.1）根据公式3.1，可以使得算法中每次循环，即每步前进不超过1个象素，从而保证了曲线的连续性。为了只使用整数运算，可将方程x=f(t)乘一正整数n，使其转换成整数方程： nxi =( i ) + zi ， 其中( i )=nf(),它和xi (x坐标的整数位置)及其余数zi ()均为整数。各个量的含义如下图所示。xi表示与计算出的值f()最近的象素点的x坐标。余数zi 表示xif(i/n)xixi f(i/n)(a) (b) 图1-1 xi与 f()两者的关系 与f()的差距（当然为了整数化，也乘了一个n值）。若是图（a）的情况，则zi为负值，即取曲线左边的象素点；若是图（b）的情况，则 zi为正值，即取曲线右边的象素点。算法逐点计算曲线上的象素。设当前象素点的x坐标是xi，则下一点的x坐标xi+1应该满足 n xi+1 = (i+1) + zi+1，其中 (i+1)=( i )+( i )=nxizi +( i )根据3.1式有：=nn，而，所以n。因此，xi+1的可能取值为xi1，xi 或xi +1。这样就得到了xi+1和zi+1的计算公式：xi+1 zi+1这样，就得到了x坐标和余数z的递推计算公式。y坐标的计算公式可以根据以上的讨论同理得到。步长大小（即1/n）的确定是非常重要的。确定步长的标准是：在使曲线前进的幅度不超过一个象素长度的前提下，选择尽量大的步长（即小的n值），以便算法提高速度。如果步长选得太小，则在绘制曲线时造成取点过密，因而重复绘制了很多象素点，浪费了大量的计算也降低了算法的速度。在参考文献提出的算法中，取n的值为： n = max (nx，ny)其中nx和ny分别对于x坐标和y坐标选择的n的值，为nx= m ny= m 其中m为曲线的次数。而该参考文献又对nx和ny值进行了进一步的优化，减小了n 的值，取： nx= m ny= m 其中x和y是曲线第i个控制点的坐标。可以证明这样选择的n值满足上述的条件n。3.2 最佳的步长值根据参考文献中对上述算法进行编程和分析运行结果，发现该算法在绘制曲线时每次循环使曲线前进的幅度的最大值往往不能接近于1（多数在0.6至0.8之间）。这说明在上一小节介绍的算法中的n的取值并没有达到最佳值。从理论上分析，我们的目标是要找到满足条件n的最小的n值，即寻找的最小上界。上述的m 虽然是的一个上界，但并不是其最小上界。下面将讨论如何求出的最小上界以及最佳的n值。的最小上界显然应该出现在f(t)的极值点处或曲线的端点处，即当参数t = 0或 t =1处，因此应该以在曲线的两端点处和其在两端点之间的极值点处的函数值中的最大者作为其最小上界。由于我们采用了求函数的最大极值作为其上界，因此它是最小上界，因为该上界就在曲线上，如果它再小一点则必有曲线上的一点大于它。下面我们将以三次bezier曲线为例描述最佳n的计算过程。对于三次bezier曲线：f( t ) = x0(1t)3 + 3x1t (1t)2 + 3x2 t2 (1t) + x3t3对其求导,得 f ( t )=3(x1x0)(1t)2 + 2(x2 x1)(1t)t + (x3x2) t2=3x3x22(x2x1)+x1x0 t2+2(x22x1+x0) t + x1x0，要求f(t)的极值点就要将其对t 再次求导并另其导函数为0，得：6x33(x2x1)x0 t + x22x1 + x0= 0 设v= x22x1 + x0，w= x33(x2x1)x0；由上式得，当t=时，f( t )取到极值3x1x0。显然，在端点处( t = 0 和t = 1)的值分别为3和3。所以，nx的取值为当0，1时 nx = 3）当0，1时，nx = 3）同理，根据ny来求得ny 的值。最后取 n = max(nx,ny)，从而得到了最佳的n值。 由上面的讨论可知，n值决定了算法的循环次数，n值愈小，算法的速度愈快。在上面计算得到的n值已经比较理想了，那么，是否还可以再进一步减小n值呢？下面提出一个双步逐点曲线生成算法来进一步减小n值。3.3 双步逐点曲线生成算法从上一小节的描述可知，新算法通过求极值已经找到了满足条件n的最小n值。如果再进一步减小n的值，则在曲线上可能出现断点，即曲线前进的幅度大于了一个象素的长度，从而漏过了一个中间的象素。在前面求得的最小n值的基础上，本小节提出双步逐点曲线生成算法，其基本思想是，将n 值减小二分之一（即将步长增加一倍）。这样当曲线前进一步时，其幅度有时会大于一个象素的长度。为了不使曲线出现漏点，我们通过插值的方法来确定跨过的那个中间象素。用此方法可使上一小节提出的算法的循环次数减少一半，并大大提高了参数曲线的生成速度。下面将具体描述双步逐点曲线生成算法的思想。设n 仍是上节计算出的单步算法的n 值。为了保证算法的双步特性，将3.1式（拉格朗日中值定理）=1乘2，得= 2该式保证了双步算法在将步长增加一倍（即n值减小二分之一）时，曲线每次前进的幅度不会大于两个象素的长度。我们还是先讨论x坐标的情况，y坐标与x坐标同理。设当前点坐标xi已经确定，则下一点的x 坐标xi+1应该满足nxi+1 = (i+1) + zi+1 ,其中 (i+1)=( i ) +( i )= nxizi +( i )由第一小节中的等式=nn可知：= n=（）2n，而，所以 n可以看出，xi+1的可能取值为xi2，xi1，xi，xi +1和xi +2。由上面的讨论就得到了xi+1的计算公式： 当 (k)n (i)zi fillrect(&rect,&brush);long u=3,n,n;long f0,f1,f2,f3,g0,g1,g2,g3;long a1,a2,a3,b1,b2,b3;long v,w,m,t; /计算nx的值long nx=labs(x1-x0)labs(x3-x2)?labs(x1-x0):labs(x3-x2);v=x2-2*x1+x0;w=x3+3*(x1-x2)-x0;t=-v/w;if(t0 & t1)m=x1-x0+v*t;if(mnx) nx=m; /计算ny的值long ny=labs(y1-y0)labs(y3-y2)?labs(y1-y0):labs(y3-y2);v=y2-2*y1+y0;w=y3+3*(y1-y2)-y0;t=-v/w;if(t0 & t1)m=y1-y0+v*t;if(mny) ny=m;nx=u*nx;ny=u*ny;n=nxny?nx:ny;/将n的值减小1/2n=n1;/使f0, f1, f2, f3,g0,g1,g2,g3都为整数n=n*n*n;/计算当t=0，t=1/n，t=2/n，t=3/n时f与g的值f0=n*x0;f1=x0*(n-1)*(n-1)*(n-1)+3*x1*(n-1)*(n-1)+3*x2*(n-1)+x3;f2=x0*(n-2)*(n-2)*(n-2)+6*x1*(n-2)*(n-2)+12*x2*(n-2)+8*x3;f3=x0*(n-3)*(n-3)*(n-3)+9*x1*(n-3)*(n-3)+27*x2*(n-3)+27*x3;g0=n*y0;g1=y0*(n-1)*(n-1)*(n-1)+3*y1*(n-1)*(n-1)+3*y2*(n-1)+y3;g2=y0*(n-2)*(n-2)*(n-2)+6*y1*(n-2)*(n-2)+12*y2*(n-2)+8*y3;g3=y0*(n-3)*(n-3)*(n-3)+9*y1*(n-3)*(n-3)+27*y2*(n-3)+27*y3;/ a1，a2，a3和b1，b2，b3分别表示f和g的一、二、三阶差分a1=f1-f0;a2=f2-2*f1+f0;a3=f3-3*f2+3*f1-f0;b1=g1-g0;b2=g2-2*g1+g0;b3=g3-3*g2+3*g1-g0;long x=x0,y=y0;long z1=0,z2=0; /画曲线的第一个点int fx=0,fy=0;pdc-setpixel(x,y,rgb(10,20,20); /循环画点for(int i=0;i=n;i+)/xi+1=xi-2时的情况 if(2*(a1-z1)-n)if(2*(a1-z1)setpixel(x,y,rgb(10,20,20);fx=0; /xi+1=xi-2时yi+1=yi-2的情况 if(2*(b1-z2)-n) if(2*(b1-z2)setpixel(x,y,rgb(10,20,20); pdc-setpixel(x-1,y-1,rgb(10,20,20); pdc-setpixel(x-2,y-2,rgb(10,20,20); y=y-2; z2=z2-b1-2*n; else /xi+1=xi-2时yi+1=yi-1的情况 if(b1-2*z2-n) /*b1-2*z2-n即b1/2-z2setpixel(x,y,rgb(10,20,20); pdc-setpixel(x-1,y-1,rgb(10,20,20); fx=1; els