/* Subroutine to calculate a Cartesian product B-spline surface using open uniform knot vectors (see Eq. 6.1}). C code for An Introduction to NURBS by David F. Rogers. Copyright (C) 2000 David F. Rogers, All rights reserved. Name: bsplsurf.c Language: C Subroutines called: knot.c, basis.c Book reference: Chapter 6, Ex. 6.1, Alg. p. 303 b[] = array containing the polygon net points b[1] = x-component b[2] = y-component b[3] = z-component Note: Bi,j = b[] has dimensions of n*m*3 with j varying fastest The polygon net is n x m d1nbasis[] = array containing 1st deriv. of the nonrational basis functions for one value of u (see \eq{5--84}) d1mbasis[] = array containing 1st deriv. of the nonrational basis functions for one value of w (see \eq{5--84}) d2nbasis[] = array containing 2nd deriv. of the nonrational basis functions for one value of u (see \eq{5--84}) d2mbasis[] = array containing 2nd deriv. of the nonrational basis functions for one value of w (see \eq{5--84}) k = order in the u direction l = order in the w direction mbasis[] = array containing the nonrational basis functions for one value of w (see \eq{5--84}) mpts = the number of polygon vertices in w direction nbasis[] = array containing the nonrational basis functions for one value of u (see \eq{5--84}) npts = the number of polygon vertices in u direction p1 = number of parametric lines in the u direction p2 = number of parametric lines in the w direction q[] = array containing the resulting surface q[1] = x-component q[2] = y-component q[3] = z-component for a fixed value of u the next m elements contain the values for the curve q[u[sub i],w] q has dimensions of p1*p2*3. The display surface is p1 x p2 qu[] = array containing the u derivative qw[] = array containing the w derivative quu[] = array containing the uu derivative quw[] = array containing the uw derivative qww[] = array containing the ww derivative */ #include dbsurf(b,k,ell,npts,mpts,p1,p2,q,qu,qw,quw,quu,qww) int k,ell; int npts,mpts; int p1,p2; float b[],q[],qu[],qw[],quw[],quu[],qww[]; { /* allows for 20 data points with basis function of order 5 */ int ch; int i,j,j1,jbas; int icount; int uinc,winc; int nplusc,mplusc; int x[30],y[30]; int temp; float nbasis[30],mbasis[30]; float d1nbasis[30],d1mbasis[30]; float d2nbasis[30],d2mbasis[30]; float pbasis, d1nbasismbasis, nbasisd1mbasis, d1nbasisd1mbasis; float d2nbasisd1mbasis, d2nbasismbasis, nbasisd2mbasis; float u,w; float stepu,stepw; /* printf("in bsplsurf.c \n");*/ /* zero and redimension the arrays */ /* printf("k,ell,npts,mpts,p1,p2 = %d %d %d %d %d %d \n",k,ell,npts,mpts,p1,p2);*/ nplusc = npts + k; mplusc = mpts + ell; for (i = 1; i <= nplusc; i++){ x[i] = 0; } for (i = 1; i <= mplusc; i++){ y[i] = 0; } for (i = 1; i <= npts; i++){ nbasis[i] = 0.; d1nbasis[i] = 0.; d2nbasis[i] = 0.; } for (i = 1; i <= mpts; i++){ mbasis[i] = 0.; d1mbasis[i] = 0.; d2mbasis[i] = 0.; } temp = 3*(p1*p2); /* printf("p1*p2*3 = %d \n",temp);*/ for (i = 1; i <= 3*p1*p2; i++){ q[i] = 0.; qu[i] = 0.; qw[i] = 0.; quw[i] = 0.; quu[i] = 0.; qww[i] = 0.; } /* generate the open uniform knot vectors */ knot(npts,k,x); /* calculate u knot vector */ knot(mpts,ell,y); /* calculate w knot vector */ icount = 1; /* calculate the points on the \bsp surface */ stepu = (float)x[nplusc]/(float)(p1-1); stepw = (float)y[mplusc]/(float)(p2-1); u = 0.; for (uinc = 1; uinc <= p1; uinc++){ /* printf("u = %f \n",u); */ if ((float)x[nplusc] - u < 5e-6){ u = (float)x[nplusc]; } dbasis(k,u,npts,x,nbasis,d1nbasis,d2nbasis); /* basis function for this value of u */ w = 0.; for (winc = 1; winc <= p2; winc++){ /* printf("w = %f \n",w);*/ if ((float)y[mplusc] - w < 5e-6){ w = (float)y[mplusc]; } dbasis(ell,w,mpts,y,mbasis,d1mbasis,d2mbasis); /* basis function for this value of w */ for (i = 1; i <= npts; i++){ jbas = 3*mpts*(i-1); for (j = 1; j <= mpts; j++){ j1 = jbas + 3*(j-1) + 1; pbasis = nbasis[i]*mbasis[j]; d1nbasismbasis = d1nbasis[i]*mbasis[j]; nbasisd1mbasis = nbasis[i]*d1mbasis[j]; d1nbasisd1mbasis = d1nbasis[i]*d1mbasis[j]; d2nbasismbasis = d2nbasis[i]*mbasis[j]; nbasisd2mbasis = nbasis[i]*d2mbasis[j]; /* printf("pbasis,d1nbasismbasis,nbasisd1mbasis,d1nbasisd1mbasis,d2nbasisd1mbasis,d2nbasismbasis,nbasisd2mbasis\n"); printf("%f %f %f %f %f %f %f \n",pbasis,d1nbasismbasis,nbasisd1mbasis,d1nbasisd1mbasis,d2nbasisd1mbasis,d2nbasismbasis,nbasisd2mbasis); */ q[icount] = q[icount]+b[j1]*pbasis; /* calculate surface point */ q[icount+1] = q[icount+1]+b[j1+1]*pbasis; q[icount+2] = q[icount+2]+b[j1+2]*pbasis; qu[icount] = qu[icount]+b[j1]*d1nbasismbasis; /* calculate u derivative */ qu[icount+1] = qu[icount+1]+b[j1+1]*d1nbasismbasis; qu[icount+2] = qu[icount+2]+b[j1+2]*d1nbasismbasis; qw[icount] = qw[icount]+b[j1]*nbasisd1mbasis; /* calculate w derivative */ qw[icount+1] = qw[icount+1]+b[j1+1]*nbasisd1mbasis; qw[icount+2] = qw[icount+2]+b[j1+2]*nbasisd1mbasis; quw[icount] = quw[icount]+b[j1]*d1nbasisd1mbasis; /* calculate uw derivative */ quw[icount+1] = quw[icount+1]+b[j1+1]*d1nbasisd1mbasis; quw[icount+2] = quw[icount+2]+b[j1+2]*d1nbasisd1mbasis; quu[icount] = quu[icount]+b[j1]*d2nbasismbasis; /* calculate uu derivative */ quu[icount+1] = quu[icount+1]+b[j1+1]*d2nbasismbasis; quu[icount+2] = quu[icount+2]+b[j1+2]*d2nbasismbasis; qww[icount] = qww[icount]+b[j1]*nbasisd2mbasis; /* calculate ww derivative */ qww[icount+1] = qww[icount+1]+b[j1+1]*nbasisd2mbasis; qww[icount+2] = qww[icount+2]+b[j1+2]*nbasisd2mbasis; /* printf("j1,i,j = %d %d %d \n",j1,i,j); */ } } printf("%f %f %f %f %f \n",u,w,qw[icount],qw[icount+1],qw[icount+2]); icount = icount + 3; w = w + stepw; } u = u + stepu; ch = getchar(); } }