blake@ccwf.cc.utexas.edu (Blake Freeburg) (03/03/91)
Description:
This file contains various things needed for the generation of nurbs
curves and surfaces. Most of this is a copy of what may be found in the
gl.h header file
History
15 Feb. 1991: Original revision
**************************************************************************/
#define XYZ 0x0001 /* 3D control point */
#define XYZW 0x0002 /* rational 4D control point */
typedef struct fvec2d { double x,y; } FVEC2D;
typedef struct fvec3d { double x,y,z; } FVEC3D;
typedef struct fvec4d { double x,y,z,w; } FVEC4D;
-----------CUT HERE---------------------------CUT HERE---------CUT HERE-----
/************************************************************************
NURBS surface generator
Programmer: Blake Freeburg
Description:
A set of functions that will produce a mesh of points given the
s and t knot vectors, the number of knots, the order of each axis,
byte steps to take in the array, and the control points, in (x,y,z) or
(x,y,z,w) format.
***************************************************************************/
#include<stdio.h>
#include<stdlib.h>
#include <math.h>
#include "nurbs.h"
#define TRIMESH /* turn on triangle generation */
long s_byte_length=0L;/* offset to next control point in s direction */
long t_byte_length=0L;/* offset to next control point in t direction */
long s_knot_count=0L; /* number of knots in knot vector in s directon */
long t_knot_count=0L; /* number of knots in knot vector in t directon */
/*************************************************************************
double N(i,k,knot,t)
long i - the index of the knot vector we are at
long k - the order of the curve
double knot[] - the knot vector or non-decreasing values
double t - the point along the curve from tmin < t < tmax
**************************************************************************/
double N(long i, long k, double knot[], double t)
{
/* temp variables, removed when finished */
double num1, num2, den1, den2, val1, val2;
if(k == 1)
{
if((knot[i]<=t)&&(knot[i+1]>t))
{
return 1.0;
}
else
{
return 0.0;
}
}
else
{
num1 = (t-knot[i])*N(i,k-1,knot,t);
den1 = knot[i+k-1] - knot[i];
num2 = (knot[i+k] - t)*N(i+1,k-1,knot,t);
den2 = knot[i+k] - knot[i+1];
if ((num1==0.0)&&(den1==0.0)) val1 = 0.0;
else val1 = num1/den1;
if ((num2==0.0)&&(den2==0.0)) val2 = 0.0;
else val2 = num2/den2;
return(val1 + val2);
}
}
/**************************************************************************
double S(i,j,s_order,t_order,s_count,t_count,sknot,tknot,ctrlarray,s,t);
long i - the index of the control point in the mesh
long j - the index of the control point in the mesh
long s_order - the order of the curve in s space
long t_order - the order of the curve in t space
long s_count - the number of control points in s space
long t_count - the number of control points in t space
double sknot[] - the knot vector of non-decreasing values in s space
double tknot[] - the knot vector of non-decreasing values in t space
double *ctrlarray - array of data points to be used
double s - the point along the surface from umin < u < umax
double t - the point along the surface from wmin < w < wmax
***************************************************************************/
double S(long i,long j,long s_order,long t_order,long s_count,long t_count,
double sknot[],double tknot[],double *ctrlarray, double s,double t)
{
long x, y;
double den=0.0, num=0.0;
char *position;
/* pointer to find values in array */
position = (char *)ctrlarray;
/* compute the numerator first */
num = (*(FVEC4D *)(position+i*s_byte_length+j*t_byte_length)).w
*N(i,s_order,sknot,s)*M(j,t_order,tknot,t);
/* now compute the denominator */
for(x=0; x<=(s_count-s_order); x++)
for(y=0; y<=(t_count-t_order); y++)
den+= (*(FVEC4D *)(position+i*s_byte_length+j*t_byte_length)).w
*N(i,s_order,sknot,s)*M(j,t_order,tknot,t);
if((num == 0.0) && (den == 0.0)) return(0.0);
else return (num/den);
}
/*************************************************************************
void precomputebasis(order,x,count,knot,basis)
long order the order of the curve in the s direction
double x where we are along the knot vector
long count the length of the knot vector
double knot[] the knot vector itself
double basis[] the array of precomputed knot vectors...
**************************************************************************/
void precomputebasis(long order,double x,long count,double knot[],double basis[])
{
int i;
for(i=0;i<=count;i++)
{
basis[i] = N(i,order,knot,x);
}
}
/*************************************************************************
void sumbasis(position, s_basis, t_basis, s_count, t_count, sum)
char *position where we are in the mesh
double s_basis[] array of s precomputed mesh values
double t_basis[] array of t precomputed mesh values
long s_count[] points in s direction
long t_count[] points in t direction
double *sum the sum of the denominator;
**************************************************************************/
void sumbasis(char *position, double s_basis[], double t_basis[], long s_count, long t_count, double *sum)
{
long i,j;
*sum = 0.0;
for(i=0;i<s_count;i++)
{
if(s_basis[i]!=0.0)
{
for(j=0;j<t_count;j++)
{
if(t_basis[j]!=0.0)
*sum += (*(FVEC4D *)(position+i*s_byte_length+j*t_byte_length)).w * s_basis[i] * t_basis[j];
}
}
}
}
/*************************************************************************
void nurbs_surface(s_count, s_order, sknot, t_count, t_order, tknot,
s_btye_step, t_byte_step, ctrlarray, type)
long s_count; the number of knot points in s direction
long s_order; order of the curve in the s direction
double sknot[]; the knot vector in s direction
long t_count; the number of knot points in t direction
long t_order; order of the curve in the t direction
double tknot[]; the knot vector in t direction
long s_byte_step; step in bytes to next s point
long t_byte_step; step in bytes to next t point
double *ctrlarray; pointer to the array of points
long type; type of curve: XYZ or XYZW
**************************************************************************/
void nthnurbssurface(long s_count, long s_order, double sknot[], long t_count, long t_order, double tknot[],
long s_byte_step, long t_byte_step, double *ctrlarray, long type)
{
int i, j;
int sd, td;
double Sknot;
double s, ds; /* location and step size in the s direction */
double t, dt; /* location and step size in the t direction */
double *s_basis; /* precomputed array of basis functions */
double *t_basis; /* precomputed array of basis functions */
double SumBasis; /* precomputed basis for point in s,t */
long resolution_s = 40;
long resolution_t = 40;
FVEC3D val;
char *position;
double SknotTotal = 0.0;
XPoints *xlist, *current;
long xcross_count = 0;
double leftbound, rightbound;
#ifdef TRIMESH
/* turn on if you want a triangle mesh to be generated */
FVEC3D *lastrow; /* the last row generated. The first row if none present */
FVEC3D *thisrow; /* the current row genrated. moved to last row when finished */
long marker;
lastrow = (FVEC3D *) malloc((resolution_s+1)*sizeof(FVEC3D));
thisrow = (FVEC3D *) malloc((resolution_s+1)*sizeof(FVEC3D));
marker = 0;
#endif
/* allocate the space for the evaluated basis functions */
s_basis = (double *) malloc((s_count+s_order)*sizeof(double));
t_basis = (double *) malloc((t_count+t_order)*sizeof(double));
/* set the position indicator so that we can hop around a bit */
position = (char *)ctrlarray;
/* first set up the global step values */
s_byte_length = s_byte_step;
t_byte_length = t_byte_step;
s_knot_count = s_count;
t_knot_count = t_count;
/* set up the step size and the interval */
ds = (sknot[s_count] - sknot[0])/(double)(resolution_s);
s=sknot[0];
dt = (tknot[t_count] - tknot[0])/(double)(resolution_t);
t=tknot[0];
for(s=tknot[0];s<=sknot[s_count];s+=ds)
{
/* precompute the basis functions */
precomputebasis(s_order,s,(s_count-s_order),sknot,s_basis);
for(t=tknot[0];t<=tknot[t_count];t+=dt)
{
val.x = val.y = val.z = 0.00000;
/* precompute the basis functions */
precomputebasis(t_order,t,(t_count-t_order),tknot,t_basis);
/* check to see if nurbs or bspline curve */
if(type==XYZW) sumbasis(position, s_basis, t_basis, (s_count-s_order), (t_count-t_order), &SumBasis);
for(i=0;i<(s_count-s_order);i++)
{
if(s_basis[i]!=0.0)
{
for(j=0;j<(t_count-t_order);j++)
{
if(t_basis[j]!= 0.0)
{
if(type == XYZW) /* nonuniform rational surface */
{
Sknot = s_basis[i]*t_basis[j]/SumBasis;
val.x += (*(FVEC4D *)(position+i*s_byte_length+j*t_byte_length)).x * Sknot;
val.y += (*(FVEC4D *)(position+i*s_byte_length+j*t_byte_length)).y * Sknot;
val.z += (*(FVEC4D *)(position+i*s_byte_length+j*t_byte_length)).z * Sknot;
}
else /* nonuniform surface - nurbs without the r */
{
Sknot = s_basis[i]*t_basis[j];
val.x += (*(FVEC4D *)(position+i*s_byte_length+j*t_byte_length)).x * Sknot;
val.y += (*(FVEC4D *)(position+i*s_byte_length+j*t_byte_length)).y * Sknot;
val.z += (*(FVEC4D *)(position+i*s_byte_length+j*t_byte_length)).z * Sknot;
}
}
}
}
}
/* At this point we have on point generated */
#ifdef TRIMESH
if(s==0.0) /* if the first in line... */
{
/* this should increment counter from 0 to s_count... */
lastrow[marker].x = val.x;
lastrow[marker].y = val.y;
lastrow[marker].z = val.z;
marker++;
}
else
{
thisrow[marker].x = val.x;
thisrow[marker].y = val.y;
thisrow[marker].z = val.z;
marker++;
}
#endif
}
* At this point we have a row of points that have been generated */
#ifdef TRIMESH
marker = 0;
#endif
#ifdef GL
/* This didn't work, but it was my first try at using GL trimeshes */
if(s!=0.0)
{
bgntmesh();
/* the first triangle */
v3f(&(thisrow[0]));
v3f(&(lastrow[0]));
v3f(&(lastrow[1]));
/* now for the rest of them */
for(i=1;i<resolution_s;i++)
{
v3f(&(thisrow[i]));
v3f(&(lastrow[i+1]));
}
/* end triangle */
v3f(&(thisrow[i+1]));
endtmesh();
free(lastrow);
lastrow = thisrow;
thisrow = (FVEC3D *) malloc((resolution_s+1)*sizeof(FVEC3D));
}
#endif
}
}