nurbs算法.txt

三维图形学NURBS算法实现源代码

应一些网友的要求，drew在这里提供以前在学校写的非均匀有理B样条（NURBS）曲线曲面造型的主要算法函数源码，供感兴趣的朋友参考。

由于这是多年前在学校所做的程序，为方便编程，坐标数据主要采用三维数组形式。当时使用计算机内存很小（好像是4M内存）,所以所取
控制点数有一定的限制(15个），其它未列的函数是一些数据结构，内存分配和空间转换的函数，和算法无关，可以不必考虑。由于这些源
码函数是从drew的论文程序中摘录出来的，所以需要做一些修改才能使用，朋友们可以参考和修改这些算法函数，加入自己的程序。

关于NURBS，曲线曲面造型算法最好的中文资料，drew推荐北航施法中教授编写的《计算机辅助几何设计与非均匀有理B样条(CAGD&NURBS)》.


　

//哈特利(Hartley)-贾德(Judd,1978)的弦长参数化算法求节点矢量U,用于由控制点正算曲线曲面

//coeff:控制顶点。
//v1,v2:U向，W向的节点序列。 count,countj:控制点行数和列数。

void Hartley_Judd(coeff,v1,v2,counti,countj)
float coeff[15][15][3];
float v1[20],v2[20];
int counti,countj;
{
    int i,s,j,k=DEGREE_J,n1=counti-1,n2=countj-1;
    float l1[20],l2[20];
    float ll01,ll11,ll02,ll12;
for(i=1;i<=n1;i++)
    l1[i]=(float)sqrt((coeff[i][0][0]-coeff[i-1][0][0])*
    (coeff[i][0][0]-coeff[i-1][0][0])+
    (coeff[i][0][1]-coeff[i-1][0][1])*
    (coeff[i][0][1]-coeff[i-1][0][1])+
    (coeff[i][0][2]-coeff[i-1][0][2])*
    (coeff[i][0][2]-coeff[i-1][0][2]));
for(i=1;i<=n2;i++)
    l2[i]=(float)sqrt((coeff[0][i][0]-coeff[0][i-1][0])*
    (coeff[0][i][0]-coeff[0][i-1][0])+
    (coeff[0][i][1]-coeff[0][i-1][1])*
    (coeff[0][i][1]-coeff[0][i-1][1])+
    (coeff[0][i][2]-coeff[0][i-1][2])*
    (coeff[0][i][2]-coeff[0][i-1][2]));
ll01=0;
ll11=0;
for(i=0;i<=k;i++) v1[i]=0;
for(i=k+1;i<=n1+1;i++)
{
    for(j=i-k;j<=i-1;j++)
        ll01=ll01+l1[j];
    for(s=k+1;s<=n1+1;s++)
        for(j=s-k;j<=s-1;j++)
            ll11=ll11+l1[j];
    v1[i]=ll01/ll11+v1[i-1];
}
for(i=n1+1;i<=n1+k+1;i++) v1[i]=1;
if(v1[2]==1) v1[2]=0;
ll02=0;
ll12=0;
for(i=0;i<=k;i++) v2[i]=0;
for(i=k+1;i<=n2+1;i++)
{
    for(j=i-k;j<=i-1;j++)
        ll02=ll02+l2[j];
    for(s=k+1;s<=n2+1;s++)
    for(j=s-k;j<=s-1;j++)
        ll12=ll12+l2[j];
    v2[i]=ll02/ll12+v2[i-1];
}
for(i=n2+1;i<=n2+k+1;i++) v2[i]=1;
    if(v2[2]==1) v2[2]=0;
}


// 伯姆(Boehm)节点插入算法：
// coeff:原控制顶点.v:原节点序列.uu:新插入节点.counti:原控制点数.d:新控制顶点


void boehm(coeff,v,uu,counti,d)
float coeff[15][3],v[20],d[15][3],uu;
int counti;
{   int r=0,g=0,i,j;
    float a[15],vd[20];
    for(i=0;i<=counti+DEGREE_J;i++)
        {   if(uu==v[i]) r++;
            else if(uu>=v[i]) g++;
        }
    g--;
for(i=0;i<=g-DEGREE_J;i++) 
    for(j=0;j<=2;j++)
        d[i][j]=coeff[i][j];
for(i=g-DEGREE_J+1;i<=g-r;i++)
    for(j=0;j<=2;j++)
        {   a[i]=(float)((uu-v[i])/(v[i+DEGREE_J]-v[i]));
            d[i][j]=(1-a[i])*coeff[i-1][j]+a[i]*coeff[i][j];
        }
for(i=g-r+1;i<=counti;i++)
    for(j=0;j<=2;j++)
        d[i][j]=coeff[i-1][j];
for(i=0;i<=g;i++)
    vd[i]=v[i];
vd[g+1]=(float)(uu);
for(i=g+1;i<=counti+DEGREE_J+1;i++)
    vd[i+1]=v[i];
for(i=0;i<=counti+DEGREE_J+1;i++) v[i]=vd[i];
} 

//德布尔-考克斯递推公式 for 重节点
void deBoor_2D(float uw,float n[15][4],int count,float v[20])
{   int i,k=0;
    for(i=0;i<=count+1;i++)
    {   if(uw>=v[i]&&uw<=v[i+1]) n[i][0]=1;
        else n[i][0]=0;
}
for(k=1;k<=DEGREE;k++)
    for(i=0;i<=count+1;i++)
    {   if((v[i+k]-v[i])==0)
        {   if((v[i+k+1]-v[i+1])==0) n[i][k]=0;
            else
                n[i][k]=((v[i+k+1]-uw)*n[i+1][k-1])/(v[i+k+1]-v[i+1]);
        }
        else if((v[i+k+1]-v[i+1])==0)
            n[i][k]=((uw-v[i])*n[i][k-1])/(v[i+k]-v[i]);
        else
            n[i][k]=(((uw-v[i])*n[i][k-1])/(v[i+k]-v[i]))
            +(((v[i+k+1]-uw)*n[i+1][k-1])/(v[i+k+1]-v[i+1]));
    }
}

//由给定曲线，曲面型值点反算控制点及节点矢量算法和处理重节点情况算法。
//cff:给定型值点数组. ww:权因子数组。 coeff:反算出的控制顶点。
//v1,v2:U向，W向的节点序列。 count,countj:控制点行数和列数。


void Counter_count(cff,ww,coeff,v1,v2,count,countj)
float cff[15][15][3],coeff[15][15][3],v1[20],v2[20],ww[15][15];
int count,countj;
{   float vu[20],bb[35],bb1[35],n[15][4],nn[15][4];
    float bx[15],by[15],bz[15];
    float u=0,p0_=HEAD_QIE,pn_=END_QIE,cc=0,ccc;
    int i,j,k,k1,kkk;
    DEGREE=3;
    DEGREE_I=3;
    DEGREE_J=3;
    CURVE=1;
    vu[0]=0;
    cc=0;
    for(i=1;i<=count-1;i++)
        cc=cc+(float)sqrt((cff[i][0][0]-cff[i-1][0][0])
        *(cff[i][0][0]-cff[i-1][0][0])
        +(cff[i][0][1]-cff[i-1][0][1])
        *(cff[i][0][1]-cff[i-1][0][1])
        +(cff[i][0][2]-cff[i-1][0][2])
        *(cff[i][0][2]-cff[i-1][0][2]));
    for(i=1;i<=count-1;i++)
        vu[i]=vu[i-1]+(float)sqrt((cff[i][0][0]-cff[i-1][0][0])
        *(cff[i][0][0]-cff[i-1][0][0])
        +(cff[i][0][1]-cff[i-1][0][1])
        *(cff[i][0][1]-cff[i-1][0][1])
        +(cff[i][0][2]-cff[i-1][0][2])
        *(cff[i][0][2]-cff[i-1][0][2]))/cc;
    for(i=0;i<=3;i++) v1[i]=0;
    for(i=4;i<=count-1;i++) v1[i]=vu[i-3];
    for(i=count;i<=count+3;i++) v1[i]=1;
    if(count==2)
    {    for(i=0;i<=3;i++) v1[i]=0;
         for(i=4;i<=7;i++) v1[i]=1.1F;
    }
for(k=0;k<=count-1;k++)
{   vu[0]=0;
    cc=0;
    for(k1=1;k1<=countj-1;k1++)
        cc=cc+(float)sqrt((cff[k][k1][0]-cff[k][k1-1][0])
        *(cff[k][k1][0]-cff[k][k1-1][0])
        +(cff[k][k1][1]-cff[k][k1-1][1])
        *(cff[k][k1][1]-cff[k][k1-1][1])
        +(cff[k][k1][2]-cff[k][k1-1][2])
        *(cff[k][k1][2]-cff[k][k1-1][2]));
    for(k1=1;k1<=countj-1;k1++)
        vu[k1]=vu[k1-1]+(float)sqrt((cff[k][k1][0]-cff[k][k1-1][0])
        *(cff[k][k1][0]-cff[k][k1-1][0])
        +(cff[k][k1][1]-cff[k][k1-1][1])
        *(cff[k][k1][1]-cff[k][k1-1][1])
        +(cff[k][k1][2]-cff[k][k1-1][2])
        *(cff[k][k1][2]-cff[k][k1-1][2]))/cc;
    for(k1=0;k1<=3;k1++) v2[k1]=0;
    for(k1=4;k1<=countj+2;k1++) v2[k1]=vu[k1-3];
    for(k1=countj+2;k1<=countj+5;k1++) v2[k1]=1;
    for(i=0;i<=2;i++)
    {   coeff[k][0][i]=cff[k][0][i];
        coeff[k][countj+1][i]=cff[k][countj-1][i];
    }
    for(i=0;i<=2;i++)
    {   coeff[k][1][i]=p0_*ww[k][0]/(2*ww[k][1])+coeff[k][0][i];
        coeff[k][countj][i]=coeff[k][countj+1][i]
        -pn_*ww[k][countj+1]/(2*ww[k][countj]);
    }
    deBoor_2D(v2[4],n,countj+2,v2);
    cc=Nurbs_bzer1(ww,u,v2[4],v1,v2,1,countj+2);
    bb[0]=ww[k][2]*n[2][3]/cc;
    bb[1]=ww[k][3]*n[3][3]/cc;
    for(i=0;i<=2;i++)
        cff[k][1][i]=cff[k][1][i]-coeff[k][1][i]*ww[k][1]*n[1][3]/cc;
    deBoor_2D(v2[countj+1],n,countj+2,v2);
    cc=Nurbs_bzer1(ww,u,v2[countj+1],v1,v2,1,countj+2);
    bb[(countj-4)*3+2]=ww[k][countj-2]*n[countj-2][3]/cc;
    bb[(countj-4)*3+3]=ww[k][countj-1]*n[countj-1][3]/cc;
    for(i=0;i<=2;i++)
        cff[k][countj-2][i]=cff[k][countj-2][i]-coeff[k][countj][i]
        *ww[k][countj]*n[countj][3]/cc;
    i=0;
    j=2;
    while((j+3)<=countj)
    {   deBoor_2D(v2[j+3],n,countj+2,v2);
        cc=Nurbs_bzer1(ww,u,v2[j+3],v1,v2,1,countj+2);
        bb[i+2]=ww[k][j]*n[j][3]/cc;
        bb[i+3]=ww[k][j+1]*n[j+1][3]/cc;
        bb[i+4]=ww[k][j+2]*n[j+2][3]/cc;
        i=i+3;
        j++;
    }
    j=1;
    i=0;
loop: while(i<=countj-2)
        { if(cff[k][i][0]==cff[k][i+1][0]&&cff[k][i][1]==cff[k][i+1][1]
            &&cff[k][i][2]==cff[k][i+1][2])
            {    j++;
                 if(cff[k][i+1][0]==cff[k][i+2][0]