NURBS 曲线曲面绘制算法

Question

小弟毕业设计，现在基本完成，就只有NURBS曲线曲面的绘制算法没有搞定，理论部分真的很难懂，其中曲面是在OPENGL的场景中绘制的。也就是绘制曲线曲面不是通过调用现有的函数，而是自己编写的算法。有没有直接可以用的算法，如果有的话请联系[email protected]，小弟在线等啊。希望哪位大侠帮帮忙，非常感谢啊。

Answer 1

在这个论坛上有人给出了算法，我贴出来吧～～～
namespace   sp//样条曲线曲面域名空间 
{
static double eps=1.0e-8;//精度

//设定精度 
void   SetEps(double   Eps) 
{ 
if(Eps >0.1)   return; 
eps=Eps; 
}  //获取精度 
double   GetEps() 
{ 
return   eps; 
}  //判断a的绝对值是否比b大 
inline   int   MBiger(double   a,double   b)
{ 
return   a> b||a <-b; 
}  //   k次规范B样条基函数 
double   N(int   i,int   k,double   su[],double   u) 
{ 
if(k==0) 
{  if(u-su[0] <eps) 
{ 
if(i==0) 
return   1.0; 
else 
return   0.0; 
}  if(u> =su[i]&&u <su[i+1])   
return   1.0; 
else   
return   0.0;  } 
else   if(k> 0) 
{  if(MBiger(su[i+k]-su[i],eps)&&MBiger(su[i+k+1]-su[i+1],eps)) 
{  return   (u-su[i])*N(i,k-1,su,u)/(su[i+k]-su[i]) 
+(su[i+k+1]-u)*N(i+1,k-1,su,u)/(su[i+k+1]-su[i+1]);  } 
else   if(!MBiger(su[i+k]-su[i],eps)&&MBiger(su[i+k+1]-su[i+1],eps)) 
{  return   (su[i+k+1]-u)*N(i+1,k-1,su,u)/(su[i+k+1]-su[i+1]);  } 
else   if(MBiger(su[i+k]-su[i],eps)&&!MBiger(su[i+k+1]-su[i+1],eps)) 
{  return   (u-su[i])*N(i,k-1,su,u)/(su[i+k]-su[i]); 

} 
else   return   0.0; 

} 
else   return   0.0; 
}  //B样条曲线 
double   BSL1(int   k,double   su[],double   d[],int   n,double   u) 
{ 
int   i; 
double   s=0.0; 
for(i=0;i <n;i++)   s=s+d[i]*N(i,k,su,u); 
return   s; 
}  //德布尔算法计算B样条曲线 
inline   double   alpha(int   l,int   j,int   k,double   su[],double   u) 
{ 
if(su[j+k+1-l]-su[j] <eps&&su[j+k+1-l]-su[j]> -eps)   return   0.0; 
return   (u-su[j])/(su[j+k+1-l]-su[j]); 
}  double   d1(int   l,int   j,int   k,double   su[],double   d[],double   u) 
{ 
if(l <=0)   return   d[j]; 
double   a=alpha(l,j,k,su,u); 
return   (1-a)*d1(l-1,j-1,k,su,d,u)+a*d1(l-1,j,k,su,d,u); 
}  double   BSL(int   k,double   su[],double   d[],int   n,double   u) 
{ 
int   i; 
for(i=1;i <n+k;i++) 
{ 
if(u <su[i])   break; 
if(i> k&&u==su[i])   break; 
} 
i--; 
return   d1(k,i,k,su,d,u); 
} 

//德布尔算法计算B样条曲线的导矢 
double   d2(int   l,int   j,int   k,double   su[],double   d[],double   u) 
{ 
if(l <=0)   return   d[j]; 
double   z=su[j+k+1-l]-su[j]; 
if(z <eps&&z> -eps)   return   0.0; 
return   (k+1-l)*(d2(l-1,j,k,su,d,u)-d2(l-1,j-1,k,su,d,u))/z; 
}  double   GetDerValue(int   k,double   su[],double   d[],int   n,int   r,double   u) 
{ 
int   i,j; 
double   s=0.0;  for(i=1;i <n+k;i++) 
{ 
if(u <su[i])   break; 
if(i> k&&u==su[i])   break; 
}  i--;  for(j=i-k+r;j <=i;j++) 
{ 
s+=d2(r,j,k,su,d,u)*N(j,k-r,su,u); 
}  return   s; 
} 
//////////////////////////////////////////////////////////////////////////////////////////// 
//B样条曲面 
double   BSS(int   k,double   su[],int   l,double   sv[],double   **d,int   m,double   u,int   n,double   v) 
{ 
int   i,j; 
double   s=0.0;  for(i=0;i <m;i++) 
for(j=0;j <n;j++)   
s=s+d[i][j]*N(i,k,su,u)*N(j,l,sv,v);  return   s; 
}  //NURBS曲线 
double   NURBSL(int   k,double   su[],double   d[],double   w[],int   n,double   u) 
{ 
int   i,t=k; 
double   s1=0.0,s2=0.0,wN; 
for(i=0;i <n+1;i++) 
{ 
if(u <=su[i]) 
{ 
t=i-1; 
break; 
} 
} 

if(t <k) t=k;  for(i=t-k;i <=t;i++) 
{ 
wN=w[i]*N(i,k,su,u); 
s1=s1+d[i]*wN; 
s2=s2+wN; 
}  if(s2> eps||s2 <-eps)   
return   s1/s2; 
else   
return   0.0; 
}  //NURBS曲面 
double   NURBSS(int   k,double   su[],int   l,double   sv[],double   **d,double   **w,int   m,double   u,int   n,double   v) 
{ 
int   i,j,tu=k,tv=l; 
double   s1=0.0,s2=0.0,wN;  for(i=0;i <m+1;i++) 
{ 
if(u <=su[i]) 
{ 
tu=i-1; 
break; 
} 
}  for(j=0;j <n+1;j++) 
{ 
if(v <=sv[j]) 
{ 
tv=j-1; 
break; 
} 
}  if(tu <k)tu=k; 
if(tv <l)tv=l;  for(i=tu-k;i <=tu;i++) 
{ 
for(j=tv-l;j <=tv;j++) 
{ 
wN=w[i][j]*sp::N(i,k,su,u)*N(j,l,sv,v); 
s1=s1+d[i][j]*wN; 
s2=s2+wN; 
} 
}  if(s2> eps||s2 <-eps)   
return   s1/s2; 
else   
return   0.0; 
}  //均匀B样条曲线 
double   UBSL(int   k,double   d[],int   n,double   u) 
{ 
double   s,*su=new   double[n+k+1]; 
int   i;  for(i=0;i <n+k+1;i++)   su[i]=double(i)/(n+k); 

s=BSL(k,su,d,n,u);  delete   []su; 
return   s; 
}  //准均匀B样条曲线 
double   QUBSL(int   k,double   d[],int   n,double   u) 
{ 
double   s,*su=new   double[n+k+1]; 
int   i;  for(i=0;i <k;i++)   su[i]=0.0; 
for(i=k;i <=n;i++)   su[i]=double(i-k)/(n-k); 
for(i=n+1;i <n+k+1;i++)   su[i]=1.0;  s=BSL(k,su,d,n,u); 

delete   []su; 
return   s; 
}  //均匀B样条曲面 
double   UBSS(int   k,int   l,double   **d,int   m,double   u,int   n,double   v) 
{ 
double   s,*su=new   double[m+k+1],*sv=new   double[n+l+1]; 
int   i;  for(i=0;i <m+k+1;i++)   su[i]=double(i)/(m+k); 
for(i=0;i <n+l+1;i++)   sv[i]=double(i)/(n+l);  s=BSS(k,su,l,sv,d,m,u,n,v); 

delete   []su; 
delete   []sv; 
return   s; 
}  //准均匀B样条曲面 
double   QUBSS(int   k,int   l,double   **d,int   m,double   u,int   n,double   v) 
{ 
double   s,*su=new   double[m+k+1],*sv=new   double[n+l+1]; 
int   i;  for(i=0;i <k;i++)   su[i]=0.0; 
for(i=k;i <=m;i++)   su[i]=double(i-k)/(m-k); 
for(i=m+1;i <m+k+1;i++)   su[i]=1.0;  for(i=0;i <l;i++)   sv[i]=0.0; 
for(i=l;i <=n;i++)   sv[i]=double(i-l)/(n-l); 
for(i=n+1;i <n+l+1;i++)   sv[i]=1.0;  s=BSS(k,su,l,sv,d,m,u,n,v);  delete   []su; 
delete   []sv; 
return   s; 
}  //非均匀三次B样条Riesenfeld方法 
void   RiBSL(double   *dx,double   *dy,int   n,double   u,double   *x,double   *y) 
{ 
int   i; 
double   *su=new   double[n+3]; 
double   *l=new   double[n]; 
double   L=0.0,l1;  l[0]=0.0;  for(i=1;i <n;i++) 
{ 
l[i]=sqrt(pow2(dx[i]-dx[i-1])+pow2(dy[i]-dy[i-1])); 
L=L+l[i]; 
}  su[0]=su[1]=su[2]=su[3]=0; 
l1=l[1];  for(i=4;i <n-1;i++) 
{ 
l1=l1+l[i-2]; 
su[i]=l1/L; 
}  su[n]=su[n+1]=su[n+2]=su[n-1]=1; 
*x=BSL(3,su,dx,n,u); 
*y=BSL(3,su,dy,n,u); 

delete   []su; 
delete   []l; 
} 
}