Path: utzoo!attcan!uunet!cs.utexas.edu!ut-emx!ccwf.cc.utexas.edu
From: blake@ccwf.cc.utexas.edu (Blake Freeburg)
Newsgroups: comp.graphics
Subject: C code for NURBS surface - Here it is, have fun
Keywords: NURBS source C 
 This is the source code that I used to generate nurbs surfaces.  Since I
	Project:  Nurbs curves and surfaces
	File:  NURBS.H
	Programmer:  Blake Freeburg
	Date:  Feb 15, 1991
Message-ID: <44983@ut-emx.uucp>
Date: 2 Mar 91 17:31:35 GMT
Sender: news@ut-emx.uucp
Reply-To: blake@ccwf.cc.utexas.edu (Blake Freeburg)
Distribution: comp.graphics
Organization: The University of Texas at Austin
Lines: 317


	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
	}		

}