计算机图形学中的拉格朗日插值算法

 

大家都知道，通过两个点的一次曲线，有且只有一条，通过三个点的二次曲线，有且仅有一条。概括来讲，通过n+1个点的n次曲线，有且仅有一条。那么给定n+1个点，如何确定这个n次方程的n+1个系数呢？即已知这n+1个点的坐标，求x到y的映射。因为这里有n+1个未知数，n+1个等式，我们可以利用矩阵解方程来求解这些系数，但这计算量很大。拉格朗日提出了一个方法，不需要实现求解任何系数即可直接带入运算来求解。简单来讲就是，

，

则。（我们可以理解为表示在结果中占的权重）

但在实际应用中，点的数量可能非常巨大，最后会导致方程次数非常高，所以通常选定3次方程作为上限，然后每次处理4个点。采用这个思想，我们可以建立y与x之间的映射，然后依次求解。但这时候仍然有个问题，每次面对不同的4个点的时候，必须重新去更新的值。

为了解决这个问题，我们给空间在增加一个维度t（事实上我们可以假定t即表示点的顺序）,然后分别建立x到t和y到t的映射，即和。因为我们每次都是选用的4个点，我们可以每次都给他们指定，这样 的值是不变的，也就是说我们在计算出来后就不必再去更新它了。

 

下面的程序代码是我根据校内教材上的源码修改的，更改了函数接口，并修正了其中的一个边界BUG。

（程序中指定）

void Lugrange(int coordinate[][2],int k, int l_section,CDC *pDC)//coordinate控制点坐标 k点数量 l_section相邻点之间描绘曲线时的中间点数量 { float blend[4][50]; float f_blend[4][50]; float l_blend[4][50]; float xsm[4]; float ysm[4]; float x=0,y=0; for(int i=0;i<l_section;i++) { float u=(float)i/(float)l_section;//根据l_section指定步长 blend[0][i]=u*(u-1)*(u-2)/(-6);//用于中间段的插值 blend[1][i]=(u+1)*(u-1)*(u-2)/2.0; blend[2][i]=(u+1)*u*(u-2)/(-2); blend[3][i]=(u+1)*u*(u-1)/6.0; float v=u-1; f_blend[0][i]=v*(v-1)*(v-2)/(-6);//用于首段的插值 f_blend[1][i]=(v+1)*(v-1)*(v-2)/2.0; f_blend[2][i]=(v+1)*v*(v-2)/(-2); f_blend[3][i]=(v+1)*v*(v-1)/6.0; float w=u+1; l_blend[0][i]=w*(w-1)*(w-2)/(-6.0);//用于末端的插值 l_blend[1][i]=(w+1)*(w-1)*(w-2)/2.0; l_blend[2][i]=(w+1)*w*(w-2)/(-2); l_blend[3][i]=(w+1)*w*(w-1)/6.0; } for(int i=0;i<4;i++) { xsm[i]=coordinate[i][0]; ysm[i]=coordinate[i][1]; } pDC->MoveTo(coordinate[0][0],coordinate[0][1]);//计算首段 for(int j=0;j<l_section;j++) { x=y=0; for(int i=0;i<4;i++) { x+=xsm[i]*f_blend[i][j]; y+=ysm[i]*f_blend[i][j]; } pDC->LineTo(x,y); } for(int m=4;m<=k;m++)//每次往后移一个点，取4个点计算中间段曲线 { for(int j=0;j<=l_section;j++)//4个点计算中间段曲线 { x=y=0; for(int i=0;i<4;i++) { x+=xsm[i]*blend[i][j]; y+=ysm[i]*blend[i][j]; } pDC->LineTo(x,y); } if(m==k) break; for(int i=0;i<3;i++)//往后移一个点 { xsm[i]=xsm[i+1]; ysm[i]=ysm[i+1]; } xsm[3]=coordinate[m][0]; ysm[3]=coordinate[m][1]; } for(int j=0;j<l_section;j++)//计算末端 { x=y=0; for(int i=0;i<4;i++) { x+=xsm[i]*l_blend[i][j]; y+=ysm[i]*l_blend[i][j]; } pDC->LineTo(x,y); } pDC->LineTo(xsm[3],ysm[3]); }