nurbs算法.txt

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

                    &&cff[k][i+1][1]==cff[k][i+2][1]
                    &&cff[k][i+1][2]==cff[k][i+2][2])
                                j++;
            }
           if(j>1) break;
            i++;
        }
        kkk=i;
    if(j>=2)
    {   for(k1=0;k1<=2;k1++) cff[k][kkk][k1]=1;
        deBoor_2D(v2[kkk+3],n,countj+2,v2);
        cc=Nurbs_bzer1(ww,u,v2[kkk+3],v1,v2,1,countj+2);
        nn[kkk][3]=3*(n[kkk][2]/(v2[kkk+3]-v2[kkk])
        -n[kkk+1][2]/(v2[kkk+4]-v2[kkk+1]));
        nn[kkk+1][3]=3*(n[kkk+1][2]/(v2[kkk+4]-v2[kkk+1])
        -n[kkk+2][2]/(v2[kkk+5]-v2[kkk+2]));
        nn[kkk+2][3]=3*(n[kkk+2][2]/(v2[kkk+5]-v2[kkk+2])
        -n[kkk+3][2]/(v2[kkk+6]-v2[kkk+3]));
        ccc=0;
        for(k1=kkk;k1<=kkk+2;k1++) ccc=ccc+ww[k][k1]*nn[k1][3];
        bb[(kkk-2)*3+2]=ww[k][kkk]*(nn[kkk][3]*cc-n[kkk][3]*ccc)/(cc*cc);
        bb[(kkk-2)*3+3]=ww[k][kkk+1]*(nn[kkk+1][3]*cc-n[kkk+1][3]*ccc)/(cc*cc);
        bb[(kkk-2)*3+4]=ww[k][kkk+2]*(nn[kkk+2][3]*cc-n[kkk+1][3]*ccc)/(cc*cc);
    }
    i=i+j;
    if(i<=countj-2) { j=1; goto loop; }
    for(i=0;i<=countj-3;i++)
    {   bx[i]=cff[k][i+1][0];
        by[i]=cff[k][i+1][1];
        bz[i]=cff[k][i+1][2];
    }
    for(i=0;i<=(countj-4)*3+3;i++) bb1[i]=bb[i]; 
    run_run(bb1,countj-2,4+3*(countj-4),bx);
    for(i=0;i<=(countj-4)*3+3;i++) bb1[i]=bb[i]; 
    run_run(bb1,countj-2,4+3*(countj-4),by);
    run_run(bb,countj-2,4+3*(countj-4),bz); 
    for(i=2;i<=countj-1;i++)
    {   coeff[k][i][0]=bx[i-2];
        coeff[k][i][1]=by[i-2];
        coeff[k][i][2]=bz[i-2];
    }
}
CURVE=0;
} 

//NURBS算法，由控制顶点，权因子，节点序列等计算插值点集
void cal_surf(SURF *ps,HDC hdc)
{   GRID *pg;
    float coeffx[15][15],coeffy[15][15],coeffz[15][15];
    int m,n,i,k,j,mark=0;
    float u,w,nn;
    o=0;p=0;
    DEGREE_I=ps->degree_u;
    DEGREE_J=ps->degree_w;
    for(i=0;i<=ps->nums_u-1;i++)
    for(j=0;j<=ps->nums_w-1;j++)
    {   coeffx[i][j]=ps->pro[i][j].x;
        coeffy[i][j]=ps->pro[i][j].y;
        coeffz[i][j]=ps->pro[i][j].z;
    }
    m=ps->sgridn_u;n=ps->sgridn_w;
    if(ROTA_FLG==0) 
    {   ALLOC(ps->pgrid,GRID);
        pg=ps->pgrid;
    } 
        if(CURVE==1) 
    {   if(NC_STYLE==1) { NC_STYLE=0; mark=1; } m=1; }
        for(i=0;i<m+1;i++)
        {   u=(float)(i*1./m);
            p=0;
            for(k=0;k<n+1;k++)
            {    w=(float)(k*1./n); 
                 if(NC_STYLE!=1)
                    {   nn=Nurbs_bzer1(ps->www,u,w,ps->v_u,ps->v_w,ps->nums_u,ps->nums_w);
                        x_[o][p]=L_S_bzer(coeffx,ps->www,u,w,ps->v_u,ps->v_w,ps->nums_u,
                        ps->nums_w)/nn;
                        y_[o][p]=L_S_bzer(coeffy,ps->www,u,w,ps->v_u,ps->v_w,ps->nums_u,
                        ps->nums_w)/nn;
                        z_[o][p]=L_S_bzer(coeffz,ps->www,u,w,ps->v_u,ps->v_w,ps->nums_u,
                        ps->nums_w)/nn;
                    }
                  else
                    {   nn=Nurbs_bzer1(ps->www,w,u,ps->v_u,ps->v_w,ps->nums_u,ps->nums_w);
                        x_[o][p]=L_S_bzer(coeffx,ps->www,w,u,ps->v_u,ps->v_w,ps->nums_u,
                        ps->nums_w)/nn;
                        y_[o][p]=L_S_bzer(coeffy,ps->www,w,u,ps->v_u,ps->v_w,ps->nums_u,
                        ps->nums_w)/nn;
                        z_[o][p]=L_S_bzer(coeffz,ps->www,w,u,ps->v_u,ps->v_w,ps->nums_u,
                        ps->nums_w)/nn;
                    } 
                  if(ROTA_FLG==0) 
                    {     ALLOC(pg->next,GRID);
                          pg=pg->next;
                        // pg->value_u=u; pg->value_w=w;
                          pg->value_g.x=x_[o][p];
                          pg->value_g.y=y_[o][p];
                          pg->value_g.z=z_[o][p];
                    } 
                   if(MAIN_W==1)
                    {    if(p==0) MoveTo(hdc,(int)(x_[o][p]),(int)(y_[o][p]));
                         else LineTo(hdc,(int)(x_[o][p]),(int)(y_[o][p])); 
                                        }
                   if(LEFT_W==1)
                    {      if(p==0) MoveTo(hdc,(int)(z_[o][p]+320),(int)(y_[o][p]));
                           else LineTo(hdc,(int)(z_[o][p]+320),(int)(y_[o][p]));
                    } 
                                      if(DOWN_W==1)
                    {   if(p==0)     MoveTo(hdc,(int)(x_[o][p]),(int)(z_[o][p]+220));
                        else LineTo(hdc,(int)(x_[o][p]),(int)(z_[o][p]+220));
                    } 
                    p++;
            }
            o++;
    }
    if(ROTA_FLG==0) pg->next=NULL;

    if(CURVE!=1) 
    {      if(MAIN_W==1)
            for(i=0;i<=p-1;i++)
            for(j=0;j<=o-1;j++)
            {     if(j==0) MoveTo(hdc,(int)(x_[j][i]),(int)(y_[j][i]));
                  LineTo(hdc,(int)(x_[j][i]),(int)(y_[j][i]));
            }
            if(LEFT_W==1) 
                for(i=0;i<=p-1;i++)
                    for(j=0;j<=o-1;j++)
                    {     if(j==0) MoveTo(hdc,(int)(z_[j][i]+320),(int)(y_[j][i]));
                          LineTo(hdc,(int)(z_[j][i]+320),(int)(y_[j][i]));
                    } 
            if(DOWN_W==1) 
                for(i=0;i<=p-1;i++)
                    for(j=0;j<=o-1;j++)
                        {     if(j==0) MoveTo(hdc,(int)(x_[j][i]),(int)(z_[j][i]+220));
                              LineTo(hdc,(int)(x_[j][i]),(int)(z_[j][i]+220));
                        } 
                }
            if(WATCH_3==1) 
                   { for(i=0;i<4;i++)
                        for(j=0;j<4;j++)
                            {   ff[i][j]=0;
                                if(i==j) ff[i][j]=0.5f;
                            }
                        ff[3][0]=(float)((1-0.25)*320);

                        ff[3][1]=(float)((1-0.25)*220);

                        ff[3][2]=0.25f; ff[3][3]=1;
                        affine3(x_,y_,z_); 
                        for(i=0;i<4;i++)
                                for(j=0;j<4;j++)
                                    {   ff[i][j]=0;
                                        if(i==j) ff[i][j]=1;
                                    }
                                ff[3][0]=-240; ff[3][1]=-140; ff[3][2]=0; ff[3][3]=1;
                        affine3(x_,y_,z_); //XOY
                        carton_draw(x_,y_,hdc); 
                        ff[3][0]=320;ff[3][1]=0;ff[3][2]=480;
                        affine3(x_,y_,z_);
                        carton_draw(z_,y_,hdc);
                        ff[3][0]=-320;ff[3][1]=0;ff[3][2]=-150;
                        affine3(x_,y_,z_);
                        carton_draw(x_,z_,hdc); 

                        ff[3][0]=-190;ff[3][1]=620;ff[3][2]=0;
                        affine3(x_,y_,z_); 
                                                rotate_any(320.,220.,0.,200.,0.,0.,45*PI/180);
                        affine3(x_,y_,z_); 
                                                rotate_any(320.,220.,0.,0.,290.,0.,45*PI/180);
                        affine3(x_,y_,z_);
                        carton_draw(x_,y_,hdc);
                } 
                if(mark==1) NC_STYLE=1; 
    }

//追赶法求解方程组
//b:系数数组。 n:方程个数。 m:系数个数。 d:未知数 

void run_run(b,n,m,d)
int n,m;
float b[],d[];
{   int k,j;
    float s;
    if(m!=(3*n-2)) return;
    for(k=0;k<=n-2;k++)
    {   j=3*k;
        s=b[j];
        if(fabs(s)+1.0==1.0)
            { MessageBox(NULL,"计算错误","系统消息",MB_OK|MB_ICONEXCLAMATION);
        return;
    }
    b[j+1]=b[j+1]/s;
    d[k]=d[k]/s;
    b[j+3]=b[j+3]-b[j+2]*b[j+1];
    d[k+1]=d[k+1]-b[j+2]*d[k];
    }
    s=b[3*n-3];
    if(fabs(s)+1.0==1.0) return;
    d[n-1]=d[n-1]/s;
    for(k=n-2;k>=0;k--)
    d[k]=d[k]-b[3*k+1]*d[k+1];
} 

float Nurbs_bzer1(ww,u,w,v1,v2,counti,countj)
float ww[15][15],v1[20],v2[20];
float u,w;
int counti,countj;
{   float n2[15][4],n1[15][4];
    float c;
    int i,j;
    DEGREE=DEGREE_I;
    deBoor_2D(u,n2,counti,v1);
    DEGREE=DEGREE_J;
    deBoor_2D(w,n1,countj,v2);
    c=0;
    if(CURVE==1) n2[0][DEGREE_I]=1;
    for(i=0;i<=counti-1;i++)
        for(j=0;j<=countj-1;j++)
            c=c+ww[i][j]*n2[i][DEGREE_I]*n1[j][DEGREE_J];
    return(c);
}