今天编写的程序：实现NURBA曲面与曲线的求交２

NURBS曲面与平行于x轴，y轴平面求交的一个简单算例２[@more@]

/*****************************************************************************

* 功能：用得布尔算法求样条曲线上的点

* 输入参数：dim表示维数；kl表示次数加1；l表示相异节点个数；u表示参数值；t[l]表示相异节点数组

* r[l]表示重复度数组；d[dim][kl]所涉及的控制顶点

* 输出参数：p[dim]表示B样条曲线上参数为u值的点

* 调用此函数的主函数中两个数组t，r,d, p都要设置大小t.SetSize(l);r.SetSize(l);

* p.SetSize(dim);

* d.SetSize(dim);

* for(int ii = 0; ii

* d[ii].SetSize(kl);

*********************************************************************************/

double Dbordv(int dim,int kl, int l, double u, CArray& t,

CArray& r, CArray& p,

CArray, CArray&>& d, int m)

{

int k = kl-1;

CArray, CArray& > dt;

dt.SetSize(dim);

for(int ii = 0; ii

dt[ii].SetSize(4);

for(int i= 0; i< dim; i++)

for(int j = 0; j <= k; j++)

{

dt[i][j] = d[i][j];

}

for(int ll =1; ll <= k; ll++)

{

for(int j = 0; j <= k-ll; j++)

{

double denom, alpha;

denom = Knot_Value(l, j+k+1+m, t,r) - Knot_Value(l, j+ll+m, t, r);

if(fabs(denom) <= 1e-5) alpha = 0.0;

else alpha = (u - Knot_Value(l, j+ll+m, t, r))/denom;

for(int i = 0; i < dim; i++)

{

dt[i][j] = (1.0-alpha)*dt[i][j] + alpha*dt[i][j+1];

}

}

}

for(i=0; i

p[i] = dt[i][0];

return p[1];

}

/*****************************************************************************

* 功能：判断所求曲线和矢量区间是不是有交点。

* 输入参数：端点数值，左端点下标-k-1； 所求曲线横坐标或纵坐标

*

* 输出参数：两端点的值减去所求曲线横坐标或纵坐标后的数值

*

*********************************************************************************/

/*double f(int dim,int kl, int l, double u, CArray& t,

CArray& r, CArray& p,

CArray, CArray&>& d,double x,int m)

*/

double f(double& u0,double& u2,int m, double xy)

{

Dbordv( dim, kl, l, u0, t, r, p,d_VCompute,m);

double x1 = p[0]-xy;

Dbordv( dim, kl, l, u2, t, r, p,d_VCompute,m);

double x2 = p[0] -xy;

double result = x1*x2;

return result;

}

int main()

{

t.SetSize(l);

r.SetSize(l);

p.SetSize(dim);

d.SetSize(dim);

result.SetSize(30);

result2.SetSize(30);

d_VCompute.SetSize(dim);

for(int ii=0;ii

{

d[ii].SetSize(kl);

d_VCompute[ii].SetSize(kl);

}

p[0]=0.0;

p[1]=0.0;

t[0]=0.0;

t[1] = 0.5;

t[2] = 1.0;

r[0] =4;

r[1] = 1;

r[2] =4;

int n=0;

/*****************************************************************

*

* 计算x=1或者-1的程序

*

********************************************************************/

for(int y =0; y

{ //本例循环两次

int m1 = y;

for(int x1 = 0; x1 < l-1; x1++) { //本例循环两次

int n1 = x1;

for(n=0; n<30; n++)

{

v =t[n1] +(t[n1+1]-t[n1])*n/30;

for(int i = m1; i

{

for (int k =0; k <3; k++)

{

for(int j= n1; j

{

d[k][j-n1]=ctlpoints[i][j][k];

// if(j==3) cout<

// cout << "d["<";

// if(j==3) cout<

}

// if(k==2) cout <

Dbordv(dim,kl, l, v, t,

r, p,

d, n1);

for( j = 0; j < dim ;j++)

d_VCompute[j][i-m1] = p[j];

}

}//for(int i = m1; i

double a = t[m1];

double b= t[m1+1];

double middle = (a+b)/2;

while(fabs(a-b)>=1e-7)

{

if(f(a,middle,m1,1)<=0)

{

b = middle;

middle =(a+b)/2;

result2[n]=middle;

continue;

}

else if(f(middle,b,m1,1)<=0)

{

a = middle;

middle = (a+b)/2;

result2[n]=middle;

continue;

}

else

{

result2[n]=-1;

break;

}

}

}

for (int i = 0; i< 30; i++)

cout <

} // end x1

}

return 0;

}

来自 “ ITPUB博客 ” ，链接：http://blog.itpub.net/409557/viewspace-893499/，如需转载，请注明出处，否则将追究法律责任。

转载于:http://blog.itpub.net/409557/viewspace-893499/