/* $Id: nurbsutl.c,v 1.1 1999/08/19 00:55:42 jtg Exp $ */

/*
 * Mesa 3-D graphics library
 * Version:  2.4
 * Copyright (C) 1995-1997  Brian Paul
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Library General Public License for more details.
 *
 * You should have received a copy of the GNU Library General Public
 * License along with this library; if not, write to the Free
 * Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */


/*
 * $Log: nurbsutl.c,v $
 * Revision 1.1  1999/08/19 00:55:42  jtg
 * Initial revision
 *
 * Revision 1.8  1999/06/08 00:44:51  brianp
 * OpenStep updates (pete@ohm.york.ac.uk)
 *
 * Revision 1.7  1998/07/26 02:07:59  brianp
 * updated for Windows compilation per Ted Jump
 *
 * Revision 1.6  1997/10/29 02:02:20  brianp
 * various MS Windows compiler changes (David Bucciarelli, v20 3dfx driver)
 *
 * Revision 1.5  1997/07/24 01:28:44  brianp
 * changed precompiled header symbol from PCH to PC_HEADER
 *
 * Revision 1.4  1997/05/28 02:29:38  brianp
 * added support for precompiled headers (PCH), inserted APIENTRY keyword
 *
 * Revision 1.3  1997/05/27 03:19:54  brianp
 * minor clean-up
 *
 * Revision 1.2  1997/05/27 03:00:16  brianp
 * incorporated Bogdan's new NURBS code
 *
 * Revision 1.1  1996/09/27 01:19:39  brianp
 * Initial revision
 *
 */


/*
 * NURBS implementation written by Bogdan Sikorski (bogdan@cira.it)
 * See README2 for more info.
 */


#ifdef PC_HEADER
#include "all.h"
#else
#include <math.h>
#include <stdlib.h>
#include "gluP.h"
#include "nurbs.h"
#endif


GLenum
test_knot(GLint nknots, GLfloat *knot, GLint order)
{
	GLsizei i;
	GLint knot_mult;
	GLfloat tmp_knot;

	tmp_knot=knot[0];
	knot_mult=1;
	for(i=1;i<nknots;i++)
	{
		if(knot[i] < tmp_knot)
			return GLU_NURBS_ERROR4;
		if(fabs(tmp_knot-knot[i]) > EPSILON)
		{
			if(knot_mult>order)
				return GLU_NURBS_ERROR5;
			knot_mult=1;
			tmp_knot=knot[i];
		}
		else
			++knot_mult;
	}
	return GLU_NO_ERROR;
}

static int
/* qsort function */
#if defined(WIN32) && !defined(OPENSTEP)
__cdecl
#endif
knot_sort(const void *a, const void *b)
{
	GLfloat x,y;

	x=*((GLfloat *)a);
	y=*((GLfloat *)b);
	if(fabs(x-y) < EPSILON)
		return 0;
	if(x > y)
		return 1;
	return -1;
}

/* insert into dest knot all values within the valid range from src knot */
/* that do not appear in dest */
void
collect_unified_knot(knot_str_type *dest, knot_str_type *src,
	GLfloat maximal_min_knot, GLfloat minimal_max_knot)
{
	GLfloat *src_knot,*dest_knot;
	GLint src_t_min,src_t_max,dest_t_min,dest_t_max;
	GLint src_nknots,dest_nknots;
	GLint i,j,k,new_cnt;
	GLboolean	not_found_flag;

	src_knot=src->unified_knot;
	dest_knot=dest->unified_knot;
	src_t_min=src->t_min;
	src_t_max=src->t_max;
	dest_t_min=dest->t_min;
	dest_t_max=dest->t_max;
	src_nknots=src->unified_nknots;
	dest_nknots=dest->unified_nknots;

	k=new_cnt=dest_nknots;
	for(i=src_t_min;i<=src_t_max;i++)
		if(src_knot[i] - maximal_min_knot > -EPSILON &&
		src_knot[i] - minimal_max_knot < EPSILON)
	{
		not_found_flag=GL_TRUE;
		for(j=dest_t_min;j<=dest_t_max;j++)
			if(fabs(dest_knot[j]-src_knot[i]) < EPSILON)
			{
				not_found_flag=GL_FALSE;
				break;
			}
		if(not_found_flag)
		{
			/* knot from src is not in dest - add this knot to dest */
			dest_knot[k++]=src_knot[i];
			++new_cnt;
			++(dest->t_max); /* the valid range widens */
			++(dest->delta_nknots); /* increment the extra knot value counter */
		}
	}
	dest->unified_nknots=new_cnt;
	qsort((void *)dest_knot,(size_t)new_cnt,(size_t)sizeof(GLfloat),
		&knot_sort);
}

/* basing on the new common knot range for all attributes set */
/* t_min and t_max values for each knot - they will be used later on */
/* by explode_knot() and calc_new_ctrl_pts */
static void
set_new_t_min_t_max(knot_str_type *geom_knot, knot_str_type *color_knot,
	knot_str_type *normal_knot, knot_str_type *texture_knot,
	GLfloat maximal_min_knot, GLfloat minimal_max_knot)
{
	GLuint	t_min,t_max,cnt;

	if(minimal_max_knot-maximal_min_knot < EPSILON)
	{
		/* knot common range empty */
		geom_knot->t_min=geom_knot->t_max=0;
		color_knot->t_min=color_knot->t_max=0;
		normal_knot->t_min=normal_knot->t_max=0;
		texture_knot->t_min=texture_knot->t_max=0;
	}
	else
	{
		if(geom_knot->unified_knot!=NULL)
		{
			cnt=geom_knot->unified_nknots;
			for(t_min=0;t_min<cnt;t_min++)
				if(fabs((geom_knot->unified_knot)[t_min] - maximal_min_knot) <
						EPSILON)
					break;
			for(t_max=cnt-1;t_max;t_max--)
				if(fabs((geom_knot->unified_knot)[t_max] - minimal_max_knot) < 
						EPSILON)
					break;
		}
		else
		if(geom_knot->nknots)
		{
			cnt=geom_knot->nknots;
			for(t_min=0;t_min<cnt;t_min++)
				if(fabs((geom_knot->knot)[t_min] - maximal_min_knot) < EPSILON)
					break;
			for(t_max=cnt-1;t_max;t_max--)
				if(fabs((geom_knot->knot)[t_max] - minimal_max_knot) < EPSILON)
					break;
		}
		geom_knot->t_min=t_min;
		geom_knot->t_max=t_max;
		if(color_knot->unified_knot!=NULL)
		{
			cnt=color_knot->unified_nknots;
			for(t_min=0;t_min<cnt;t_min++)
				if(fabs((color_knot->unified_knot)[t_min] - maximal_min_knot) <
						EPSILON)
					break;
			for(t_max=cnt-1;t_max;t_max--)
				if(fabs((color_knot->unified_knot)[t_max] - minimal_max_knot) < 
						EPSILON)
					break;
			color_knot->t_min=t_min;
			color_knot->t_max=t_max;
		}
		if(normal_knot->unified_knot!=NULL)
		{
			cnt=normal_knot->unified_nknots;
			for(t_min=0;t_min<cnt;t_min++)
				if(fabs((normal_knot->unified_knot)[t_min] - maximal_min_knot) <
						EPSILON)
					break;
			for(t_max=cnt-1;t_max;t_max--)
				if(fabs((normal_knot->unified_knot)[t_max] - minimal_max_knot) < 
						EPSILON)
					break;
			normal_knot->t_min=t_min;
			normal_knot->t_max=t_max;
		}
		if(texture_knot->unified_knot!=NULL)
		{
			cnt=texture_knot->unified_nknots;
			for(t_min=0;t_min<cnt;t_min++)
				if(fabs((texture_knot->unified_knot)[t_min] - maximal_min_knot) 
						< EPSILON)
					break;
			for(t_max=cnt-1;t_max;t_max--)
				if(fabs((texture_knot->unified_knot)[t_max] - minimal_max_knot) 
						< EPSILON)
					break;
			texture_knot->t_min=t_min;
			texture_knot->t_max=t_max;
		}
	}
}

/* modify all knot valid ranges in such a way that all have the same */
/* range, common to all knots */
/* do this by knot insertion */
GLenum
select_knot_working_range(GLUnurbsObj *nobj,knot_str_type *geom_knot,
	knot_str_type *color_knot, knot_str_type *normal_knot,
	knot_str_type *texture_knot)
{
	GLint max_nknots;
	GLfloat maximal_min_knot,minimal_max_knot;
	GLint i;

	/* find the maximum modified knot length */
	max_nknots=geom_knot->nknots;
	if(color_knot->unified_knot)
		max_nknots+=color_knot->nknots;
	if(normal_knot->unified_knot)
		max_nknots+=normal_knot->nknots;
	if(texture_knot->unified_knot)
		max_nknots+=texture_knot->nknots;
	maximal_min_knot=(geom_knot->knot)[geom_knot->t_min];
	minimal_max_knot=(geom_knot->knot)[geom_knot->t_max];
	/* any attirb data ? */
	if(max_nknots!=geom_knot->nknots)
	{
		/* allocate space for the unified knots */
		if((geom_knot->unified_knot=
				(GLfloat *)malloc(sizeof(GLfloat)*max_nknots))==NULL)
		{
			call_user_error(nobj,GLU_OUT_OF_MEMORY);
			return GLU_ERROR;
		}
		/* copy the original knot to the unified one */
		geom_knot->unified_nknots=geom_knot->nknots;
		for(i=0;i<geom_knot->nknots;i++)
			(geom_knot->unified_knot)[i]=(geom_knot->knot)[i];
		if(color_knot->unified_knot)
		{
			if((color_knot->knot)[color_knot->t_min] - maximal_min_knot >
					EPSILON)
				maximal_min_knot=(color_knot->knot)[color_knot->t_min];
			if(minimal_max_knot - (color_knot->knot)[color_knot->t_max] >
					EPSILON)
				minimal_max_knot=(color_knot->knot)[color_knot->t_max];
			if((color_knot->unified_knot=
					(GLfloat *)malloc(sizeof(GLfloat)*max_nknots))==NULL)
			{
				free(geom_knot->unified_knot);
				call_user_error(nobj,GLU_OUT_OF_MEMORY);
				return GLU_ERROR;
			}
			/* copy the original knot to the unified one */
			color_knot->unified_nknots=color_knot->nknots;
			for(i=0;i<color_knot->nknots;i++)
				(color_knot->unified_knot)[i]=(color_knot->knot)[i];
		}
		if(normal_knot->unified_knot)
		{
			if((normal_knot->knot)[normal_knot->t_min] - maximal_min_knot >
					EPSILON)
				maximal_min_knot=(normal_knot->knot)[normal_knot->t_min];
			if(minimal_max_knot - (normal_knot->knot)[normal_knot->t_max] >
					EPSILON)
				minimal_max_knot=(normal_knot->knot)[normal_knot->t_max];
			if((normal_knot->unified_knot=
					(GLfloat *)malloc(sizeof(GLfloat)*max_nknots))==NULL)
			{
				free(geom_knot->unified_knot);
				free(color_knot->unified_knot);
				call_user_error(nobj,GLU_OUT_OF_MEMORY);
				return GLU_ERROR;
			}
			/* copy the original knot to the unified one */
			normal_knot->unified_nknots=normal_knot->nknots;
			for(i=0;i<normal_knot->nknots;i++)
				(normal_knot->unified_knot)[i]=(normal_knot->knot)[i];
		}
		if(texture_knot->unified_knot)
		{
			if((texture_knot->knot)[texture_knot->t_min] - maximal_min_knot >
					EPSILON)
				maximal_min_knot=(texture_knot->knot)[texture_knot->t_min];
			if(minimal_max_knot - (texture_knot->knot)[texture_knot->t_max] >
					EPSILON)
				minimal_max_knot=(texture_knot->knot)[texture_knot->t_max];
			if((texture_knot->unified_knot=
					(GLfloat *)malloc(sizeof(GLfloat)*max_nknots))==NULL)
			{
				free(geom_knot->unified_knot);
				free(color_knot->unified_knot);
				free(normal_knot->unified_knot);
				call_user_error(nobj,GLU_OUT_OF_MEMORY);
				return GLU_ERROR;
			}
			/* copy the original knot to the unified one */
			texture_knot->unified_nknots=texture_knot->nknots;
			for(i=0;i<texture_knot->nknots;i++)
				(texture_knot->unified_knot)[i]=(texture_knot->knot)[i];
		}
		/* work on the geometry knot with all additional knot values */
		/* appearing in attirbutive knots */
		if(minimal_max_knot-maximal_min_knot < EPSILON)
		{
			/* empty working range */
			geom_knot->unified_nknots=0;
			color_knot->unified_nknots=0;
			normal_knot->unified_nknots=0;
			texture_knot->unified_nknots=0;
		}
		else
		{
			if(color_knot->unified_knot)
				collect_unified_knot(geom_knot,color_knot,maximal_min_knot,
					minimal_max_knot);
			if(normal_knot->unified_knot)
				collect_unified_knot(geom_knot,normal_knot,maximal_min_knot,
					minimal_max_knot);
			if(texture_knot->unified_knot)
				collect_unified_knot(geom_knot,texture_knot,maximal_min_knot,
					minimal_max_knot);
			/* since we have now built the "unified" geometry knot */
			/* add same knot values to all attributive knots */
			if(color_knot->unified_knot)
				collect_unified_knot(color_knot,geom_knot,maximal_min_knot,
					minimal_max_knot);
			if(normal_knot->unified_knot)
				collect_unified_knot(normal_knot,geom_knot,maximal_min_knot,
					minimal_max_knot);
			if(texture_knot->unified_knot)
				collect_unified_knot(texture_knot,geom_knot,maximal_min_knot,
					minimal_max_knot);
		}
	}
	set_new_t_min_t_max(geom_knot,color_knot,normal_knot,texture_knot,
		maximal_min_knot,minimal_max_knot);
	return GLU_NO_ERROR;
}

void
free_unified_knots(knot_str_type *geom_knot, knot_str_type *color_knot,
	knot_str_type *normal_knot, knot_str_type *texture_knot)
{
	if(geom_knot->unified_knot)
		free(geom_knot->unified_knot);
	if(color_knot->unified_knot)
		free(color_knot->unified_knot);
	if(normal_knot->unified_knot)
		free(normal_knot->unified_knot);
	if(texture_knot->unified_knot)
		free(texture_knot->unified_knot);
}

GLenum
explode_knot(knot_str_type *the_knot)
{
	GLfloat *knot,*new_knot;
	GLint nknots,n_new_knots=0;
	GLint t_min,t_max;
	GLint ord;
	GLsizei i,j,k;
	GLfloat tmp_float;

	if(the_knot->unified_knot)
	{
		knot=the_knot->unified_knot;
		nknots=the_knot->unified_nknots;
	}
	else
	{
		knot=the_knot->knot;
		nknots=the_knot->nknots;
	}
	ord=the_knot->order;
	t_min=the_knot->t_min;
	t_max=the_knot->t_max;

	for(i=t_min;i<=t_max;)
	{
		tmp_float=knot[i];
		for(j=0;j<ord && (i+j)<=t_max;j++)
			if(fabs(tmp_float-knot[i+j])>EPSILON)
				break;
		n_new_knots+=ord-j;
		i+=j;
	}
	/* alloc space for new_knot */
	if((new_knot=(GLfloat *)malloc(sizeof(GLfloat)*(nknots+n_new_knots)))==NULL)
	{
		return GLU_OUT_OF_MEMORY;
	}
	/* fill in new knot */
	for(j=0;j<t_min;j++)
		new_knot[j]=knot[j];
	for(i=j;i<=t_max;i++)
	{
		tmp_float=knot[i];
		for(k=0;k<ord;k++)
		{
			new_knot[j++]=knot[i];
			if(tmp_float==knot[i+1])
				i++;
		}
	}
	for(i=t_max+1;i<(int)nknots;i++)
		new_knot[j++]=knot[i];
	/* fill in the knot structure */
	the_knot->new_knot=new_knot;
	the_knot->delta_nknots+=n_new_knots;
	the_knot->t_max+=n_new_knots;
	return GLU_NO_ERROR;
}

GLenum
calc_alphas(knot_str_type *the_knot)
{
	GLfloat tmp_float;
	int i,j,k,m,n;
	int order;
	GLfloat *alpha,*alpha_new,*tmp_alpha;
	GLfloat denom;
	GLfloat *knot,*new_knot;


	knot=the_knot->knot;
	order=the_knot->order;
	new_knot=the_knot->new_knot;
	n=the_knot->nknots-the_knot->order;
	m=n+the_knot->delta_nknots;
	if((alpha=(GLfloat *)malloc(sizeof(GLfloat)*n*m))==NULL)
	{
		return GLU_OUT_OF_MEMORY;
	}
	if((alpha_new=(GLfloat *)malloc(sizeof(GLfloat)*n*m))==NULL)
	{
		free(alpha);
		return GLU_OUT_OF_MEMORY;
	}
	for(j=0;j<m;j++)
	{
		for(i=0;i<n;i++)
		{
			if((knot[i] <= new_knot[j]) && (new_knot[j] < knot[i+1]))
				tmp_float=1.0;
			else
				tmp_float=0.0;
			alpha[i+j*n]=tmp_float;
		}
	}
	for(k=1;k<order;k++)
	{
		for(j=0;j<m;j++)
			for(i=0;i<n;i++)
			{
				denom=knot[i+k]-knot[i];
				if(fabs(denom)<EPSILON)
					tmp_float=0.0;
				else
					tmp_float=(new_knot[j+k]-knot[i])/denom*
						alpha[i+j*n];
				denom=knot[i+k+1]-knot[i+1];
				if(fabs(denom)>EPSILON)
					tmp_float+=(knot[i+k+1]-new_knot[j+k])/denom*
						alpha[(i+1)+j*n];
				alpha_new[i+j*n]=tmp_float;
			}
		tmp_alpha=alpha_new;
		alpha_new=alpha;
		alpha=tmp_alpha;
	}
	the_knot->alpha=alpha;
	free(alpha_new);
	return GLU_NO_ERROR;
}

GLenum
calc_new_ctrl_pts(GLfloat *ctrl,GLint stride,knot_str_type *the_knot,
	GLint dim,GLfloat **new_ctrl,GLint *ncontrol)
{
	GLsizei i,j,k,l,m,n;
	GLsizei index1,index2;
	GLfloat *alpha;
	GLfloat *new_knot;

	new_knot=the_knot->new_knot;
	n=the_knot->nknots-the_knot->order;
	alpha=the_knot->alpha;

	m=the_knot->t_max+1-the_knot->t_min-the_knot->order;
	k=the_knot->t_min;
	/* allocate space for new control points */
	if((*new_ctrl=(GLfloat *)malloc(sizeof(GLfloat)*dim*m))==NULL)
	{
		return GLU_OUT_OF_MEMORY;
	}
	for(j=0;j<m;j++)
	{
		for(l=0;l<dim;l++)
			(*new_ctrl)[j*dim+l]=0.0;
		for(i=0;i<n;i++)
		{
			index1=i+(j+k)*n;
			index2=i*stride;
			for(l=0;l<dim;l++)
				(*new_ctrl)[j*dim+l]+=alpha[index1]*ctrl[index2+l];
		}
	}
	*ncontrol=(GLint)m;
	return GLU_NO_ERROR;
}

static GLint
calc_factor(GLfloat *pts,GLint order,GLint indx,GLint stride,GLfloat tolerance,
	GLint dim)
{
	GLdouble model[16],proj[16];
	GLint viewport[4];
	GLdouble x,y,z,w,winx1,winy1,winz,winx2,winy2;
	GLint i;
	GLdouble len,dx,dy;

	glGetDoublev(GL_MODELVIEW_MATRIX,model);
	glGetDoublev(GL_PROJECTION_MATRIX,proj);
	glGetIntegerv(GL_VIEWPORT,viewport);
	if(dim==4)
	{
		w=(GLdouble)pts[indx+3];
		x=(GLdouble)pts[indx]/w;
		y=(GLdouble)pts[indx+1]/w;
		z=(GLdouble)pts[indx+2]/w;
		gluProject(x,y,z,model,proj,viewport,&winx1,&winy1,&winz);
		len=0.0;
		for(i=1;i<order;i++)
		{
			w=(GLdouble)pts[indx+i*stride+3];
			x=(GLdouble)pts[indx+i*stride]/w;
			y=(GLdouble)pts[indx+i*stride+1]/w;
			z=(GLdouble)pts[indx+i*stride+2]/w;
			if(gluProject(x,y,z,model,proj,viewport,&winx2,&winy2,&winz))
			{
				dx=winx2-winx1;
				dy=winy2-winy1;
				len+=sqrt(dx*dx+dy*dy);
			}
			winx1=winx2; winy1=winy2;
		}
	}
	else
	{
		x=(GLdouble)pts[indx];
		y=(GLdouble)pts[indx+1];
		if(dim==2)
			z=0.0;
		else
			z=(GLdouble)pts[indx+2];
		gluProject(x,y,z,model,proj,viewport,&winx1,&winy1,&winz);
		len=0.0;
		for(i=1;i<order;i++)
		{
			x=(GLdouble)pts[indx+i*stride];
			y=(GLdouble)pts[indx+i*stride+1];
			if(dim==2)
				z=0.0;
			else
				z=(GLdouble)pts[indx+i*stride+2];
			if(gluProject(x,y,z,model,proj,viewport,&winx2,&winy2,&winz))
			{
				dx=winx2-winx1;
				dy=winy2-winy1;
				len+=sqrt(dx*dx+dy*dy);
			}
			winx1=winx2; winy1=winy2;
		}
	}
	len /= tolerance;
	return ((GLint)len+1);
}

/* we can't use the Mesa evaluators - no way to get the point coords */
/* so we use our own Bezier point calculus routines */
/* because I'm lazy, I reuse the ones from eval.c */

static void
bezier_curve(GLfloat *cp, GLfloat *out, GLfloat t,
	GLuint dim, GLuint order, GLint offset)
{
	GLfloat s, powert;
	GLuint i, k, bincoeff;

	if(order >= 2)
	{ 
		bincoeff = order-1;
		s = 1.0-t;

		for(k=0; k<dim; k++)
			out[k] = s*cp[k] + bincoeff*t*cp[offset+k];

		for(i=2, cp+=2*offset, powert=t*t; i<order; i++, powert*=t, cp +=offset)
		{
			bincoeff *= order-i;
			bincoeff /= i;

			for(k=0; k<dim; k++)
				out[k] = s*out[k] + bincoeff*powert*cp[k];
		}
	}
	else /* order=1 -> constant curve */
	{ 
		for(k=0; k<dim; k++)
			out[k] = cp[k];
	} 
}

static GLint
calc_parametric_factor(GLfloat *pts,GLint order,GLint indx,GLint stride,
	GLfloat tolerance,GLint dim)
{
	GLdouble model[16],proj[16];
	GLint viewport[4];
	GLdouble x,y,z,w,x1,y1,z1,x2,y2,z2,x3,y3,z3;
	GLint i;
	GLint P;
	GLfloat bez_pt[4];
	GLdouble len=0.0,tmp,z_med;

	P = 2*(order+2);
	glGetDoublev(GL_MODELVIEW_MATRIX,model);
	glGetDoublev(GL_PROJECTION_MATRIX,proj);
	glGetIntegerv(GL_VIEWPORT,viewport);
	z_med = (viewport[2] + viewport[3]) * 0.5;
	switch(dim)
	{
		case 4:
			for(i=1;i<P;i++)
			{
				bezier_curve(pts+indx, bez_pt, (GLfloat)i/(GLfloat)P, 4,
					order,stride);
				w = (GLdouble)bez_pt[3];
				x = (GLdouble)bez_pt[0] / w;
				y = (GLdouble)bez_pt[1] / w;
				z = (GLdouble)bez_pt[2] / w;
				gluProject(x,y,z,model,proj,viewport,&x3,&y3,&z3);
				z3 *= z_med;
				bezier_curve(pts+indx, bez_pt, (GLfloat)(i-1)/(GLfloat)P, 4,
					order,stride);
				w = (GLdouble)bez_pt[3];
				x = (GLdouble)bez_pt[0] / w;
				y = (GLdouble)bez_pt[1] / w;
				z = (GLdouble)bez_pt[2] / w;
				gluProject(x,y,z,model,proj,viewport,&x1,&y1,&z1);
				z1 *= z_med;
				bezier_curve(pts+indx, bez_pt, (GLfloat)(i+1)/(GLfloat)P, 4,
					order,stride);
				w = (GLdouble)bez_pt[3];
				x = (GLdouble)bez_pt[0] / w;
				y = (GLdouble)bez_pt[1] / w;
				z = (GLdouble)bez_pt[2] / w;
				gluProject(x,y,z,model,proj,viewport,&x2,&y2,&z2);
				z2 *= z_med;
				/* calc distance between point (x3,y3,z3) and line segment */
				/* <x1,y1,z1><x2,y2,z2> */
				x = x2-x1;
				y = y2-y1;
				z = z2-z1;
				tmp = sqrt(x*x+y*y+z*z);
				x /= tmp;
				y /= tmp;
				z /= tmp;
				tmp = x3*x+y3*y+z3*z-x1*x-y1*y-z1*z;
				x = x1+x*tmp-x3;
				y = y1+y*tmp-y3;
				z = z1+z*tmp-z3;
				tmp = sqrt(x*x+y*y+z*z);
				if(tmp > len)
					len = tmp;
			}
			break;
		case 3:
			for(i=1;i<P;i++)
			{
				bezier_curve(pts+indx, bez_pt, (GLfloat)i/(GLfloat)P, 3,
					order,stride);
				x = (GLdouble)bez_pt[0];
				y = (GLdouble)bez_pt[1];
				z = (GLdouble)bez_pt[2];
				gluProject(x,y,z,model,proj,viewport,&x3,&y3,&z3);
				z3 *= z_med;
				bezier_curve(pts+indx, bez_pt, (GLfloat)(i-1)/(GLfloat)P, 3,
					order,stride);
				x = (GLdouble)bez_pt[0];
				y = (GLdouble)bez_pt[1];
				z = (GLdouble)bez_pt[2];
				gluProject(x,y,z,model,proj,viewport,&x1,&y1,&z1);
				z1 *= z_med;
				bezier_curve(pts+indx, bez_pt, (GLfloat)(i+1)/(GLfloat)P, 3,
					order,stride);
				x = (GLdouble)bez_pt[0];
				y = (GLdouble)bez_pt[1];
				z = (GLdouble)bez_pt[2];
				gluProject(x,y,z,model,proj,viewport,&x2,&y2,&z2);
				z2 *= z_med;
				/* calc distance between point (x3,y3,z3) and line segment */
				/* <x1,y1,z1><x2,y2,z2> */
				x = x2-x1;
				y = y2-y1;
				z = z2-z1;
				tmp = sqrt(x*x+y*y+z*z);
				x /= tmp;
				y /= tmp;
				z /= tmp;
				tmp = x3*x+y3*y+z3*z-x1*x-y1*y-z1*z;
				x = x1+x*tmp-x3;
				y = y1+y*tmp-y3;
				z = z1+z*tmp-z3;
				tmp = sqrt(x*x+y*y+z*z);
				if(tmp > len)
					len = tmp;
			}
			break;
		case 2:
			for(i=1;i<P;i++)
			{
				bezier_curve(pts+indx, bez_pt, (GLfloat)i/(GLfloat)P, 2,
					order,stride);
				x = (GLdouble)bez_pt[0];
				y = (GLdouble)bez_pt[1];
				z = 0.0;
				gluProject(x,y,z,model,proj,viewport,&x3,&y3,&z3);
				z3 *= z_med;
				bezier_curve(pts+indx, bez_pt, (GLfloat)(i-1)/(GLfloat)P, 2,
					order,stride);
				x = (GLdouble)bez_pt[0];
				y = (GLdouble)bez_pt[1];
				z = 0.0;
				gluProject(x,y,z,model,proj,viewport,&x1,&y1,&z1);
				z1 *= z_med;
				bezier_curve(pts+indx, bez_pt, (GLfloat)(i+1)/(GLfloat)P, 2,
					order,stride);
				x = (GLdouble)bez_pt[0];
				y = (GLdouble)bez_pt[1];
				z = 0.0;
				gluProject(x,y,z,model,proj,viewport,&x2,&y2,&z2);
				z2 *= z_med;
				/* calc distance between point (x3,y3,z3) and line segment */
				/* <x1,y1,z1><x2,y2,z2> */
				x = x2-x1;
				y = y2-y1;
				z = z2-z1;
				tmp = sqrt(x*x+y*y+z*z);
				x /= tmp;
				y /= tmp;
				z /= tmp;
				tmp = x3*x+y3*y+z3*z-x1*x-y1*y-z1*z;
				x = x1+x*tmp-x3;
				y = y1+y*tmp-y3;
				z = z1+z*tmp-z3;
				tmp = sqrt(x*x+y*y+z*z);
				if(tmp > len)
					len = tmp;
			}
			break;

	}
	if(len < tolerance)
		return (order);
	else
		return (GLint)(sqrt(len/tolerance)*(order+2)+1);
}

static GLenum
calc_sampling_3D(new_ctrl_type *new_ctrl, GLfloat tolerance, GLint dim,
	GLint uorder, GLint vorder, GLint **ufactors, GLint **vfactors)
{
	GLfloat		*ctrl;
	GLint		tmp_factor1,tmp_factor2;
	GLint		ufactor_cnt,vfactor_cnt;
	GLint		offset1,offset2,offset3;
	GLint		i,j;

	ufactor_cnt=new_ctrl->s_bezier_cnt;
	vfactor_cnt=new_ctrl->t_bezier_cnt;
	if((*ufactors=(GLint *)malloc(sizeof(GLint)*ufactor_cnt*3))
			==NULL)
	{
		return GLU_OUT_OF_MEMORY;
	}
	if((*vfactors=(GLint *)malloc(sizeof(GLint)*vfactor_cnt*3))
			==NULL)
	{
		free(*ufactors);
		return GLU_OUT_OF_MEMORY;
	}
	ctrl=new_ctrl->geom_ctrl;
	offset1=new_ctrl->geom_t_stride*vorder;
	offset2=new_ctrl->geom_s_stride*uorder;
	for(j=0;j<vfactor_cnt;j++)
	{
		*(*vfactors+j*3+1)=tmp_factor1=calc_factor(ctrl,vorder,
			j*offset1,dim,tolerance,dim);
		/* loop ufactor_cnt-1 times */
		for(i=1;i<ufactor_cnt;i++)
		{
			tmp_factor2=calc_factor(ctrl,vorder,
				j*offset1+i*offset2,dim,tolerance,dim);
			if(tmp_factor2>tmp_factor1)
				tmp_factor1=tmp_factor2;
		}
		/* last time for the opposite edge */
		*(*vfactors+j*3+2)=tmp_factor2=calc_factor(ctrl,vorder,
			j*offset1+i*offset2-new_ctrl->geom_s_stride,
			dim,tolerance,dim);
		if(tmp_factor2>tmp_factor1)
			*(*vfactors+j*3)=tmp_factor2;
		else
			*(*vfactors+j*3)=tmp_factor1;
	}
	offset3=new_ctrl->geom_s_stride;
	offset2=new_ctrl->geom_s_stride*uorder;
	for(j=0;j<ufactor_cnt;j++)
	{
		*(*ufactors+j*3+1)=tmp_factor1=calc_factor(ctrl,uorder,
			j*offset2,offset3,tolerance,dim);
		/* loop vfactor_cnt-1 times */
		for(i=1;i<vfactor_cnt;i++)
		{
			tmp_factor2=calc_factor(ctrl,uorder,
				j*offset2+i*offset1,offset3,tolerance,dim);
			if(tmp_factor2>tmp_factor1)
				tmp_factor1=tmp_factor2;
		}
		/* last time for the opposite edge */
		*(*ufactors+j*3+2)=tmp_factor2=calc_factor(ctrl,uorder,
			j*offset2+i*offset1-new_ctrl->geom_t_stride,
			offset3,tolerance,dim);
		if(tmp_factor2>tmp_factor1)
			*(*ufactors+j*3)=tmp_factor2;
		else
			*(*ufactors+j*3)=tmp_factor1;
	}
	return GL_NO_ERROR;
}

static GLenum
calc_sampling_param_3D(new_ctrl_type *new_ctrl, GLfloat tolerance, GLint dim,
	GLint uorder, GLint vorder, GLint **ufactors, GLint **vfactors)
{
	GLfloat		*ctrl;
	GLint		tmp_factor1,tmp_factor2;
	GLint		ufactor_cnt,vfactor_cnt;
	GLint		offset1,offset2,offset3;
	GLint		i,j;

	ufactor_cnt=new_ctrl->s_bezier_cnt;
	vfactor_cnt=new_ctrl->t_bezier_cnt;
	if((*ufactors=(GLint *)malloc(sizeof(GLint)*ufactor_cnt*3))
			==NULL)
	{
		return GLU_OUT_OF_MEMORY;
	}
	if((*vfactors=(GLint *)malloc(sizeof(GLint)*vfactor_cnt*3))
			==NULL)
	{
		free(*ufactors);
		return GLU_OUT_OF_MEMORY;
	}
	ctrl=new_ctrl->geom_ctrl;
	offset1=new_ctrl->geom_t_stride*vorder;
	offset2=new_ctrl->geom_s_stride*uorder;
	for(j=0;j<vfactor_cnt;j++)
	{
		*(*vfactors+j*3+1)=tmp_factor1=calc_parametric_factor(ctrl,vorder,
			j*offset1,dim,tolerance,dim);
		/* loop ufactor_cnt-1 times */
		for(i=1;i<ufactor_cnt;i++)
		{
			tmp_factor2=calc_parametric_factor(ctrl,vorder,
				j*offset1+i*offset2,dim,tolerance,dim);
			if(tmp_factor2>tmp_factor1)
				tmp_factor1=tmp_factor2;
		}
		/* last time for the opposite edge */
		*(*vfactors+j*3+2)=tmp_factor2=calc_parametric_factor(ctrl,vorder,
			j*offset1+i*offset2-new_ctrl->geom_s_stride,
			dim,tolerance,dim);
		if(tmp_factor2>tmp_factor1)
			*(*vfactors+j*3)=tmp_factor2;
		else
			*(*vfactors+j*3)=tmp_factor1;
	}
	offset3=new_ctrl->geom_s_stride;
	offset2=new_ctrl->geom_s_stride*uorder;
	for(j=0;j<ufactor_cnt;j++)
	{
		*(*ufactors+j*3+1)=tmp_factor1=calc_parametric_factor(ctrl,uorder,
			j*offset2,offset3,tolerance,dim);
		/* loop vfactor_cnt-1 times */
		for(i=1;i<vfactor_cnt;i++)
		{
			tmp_factor2=calc_parametric_factor(ctrl,uorder,
				j*offset2+i*offset1,offset3,tolerance,dim);
			if(tmp_factor2>tmp_factor1)
				tmp_factor1=tmp_factor2;
		}
		/* last time for the opposite edge */
		*(*ufactors+j*3+2)=tmp_factor2=calc_parametric_factor(ctrl,uorder,
			j*offset2+i*offset1-new_ctrl->geom_t_stride,
			offset3,tolerance,dim);
		if(tmp_factor2>tmp_factor1)
			*(*ufactors+j*3)=tmp_factor2;
		else
			*(*ufactors+j*3)=tmp_factor1;
	}
	return GL_NO_ERROR;
}

static GLenum
calc_sampling_2D(GLfloat *ctrl, GLint cnt, GLint order,
	GLfloat tolerance, GLint dim, GLint **factors)
{
	GLint		factor_cnt;
	GLint		tmp_factor;
	GLint		offset;
	GLint		i;

	factor_cnt=cnt/order;
	if((*factors=(GLint *)malloc(sizeof(GLint)*factor_cnt))==NULL)
	{
		return GLU_OUT_OF_MEMORY;
	}
	offset=order*dim;
	for(i=0;i<factor_cnt;i++)
	{
		tmp_factor=calc_factor(ctrl,order,i*offset,dim,tolerance,dim);
		if(tmp_factor == 0)
			(*factors)[i]=1;
		else
			(*factors)[i]=tmp_factor;
	}
	return GL_NO_ERROR;
}

static void
set_sampling_and_culling( GLUnurbsObj *nobj )
{
	if(nobj->auto_load_matrix==GL_FALSE)
	{
		GLint i;
		GLfloat m[4];

		glPushAttrib( (GLbitfield) (GL_VIEWPORT_BIT | GL_TRANSFORM_BIT));
		for(i=0;i<4;i++)
			m[i]=nobj->sampling_matrices.viewport[i];
		glViewport(m[0],m[1],m[2],m[3]);
		glMatrixMode(GL_PROJECTION);
		glPushMatrix();
		glLoadMatrixf(nobj->sampling_matrices.proj);
		glMatrixMode(GL_MODELVIEW);
		glPushMatrix();
		glLoadMatrixf(nobj->sampling_matrices.model);
	}
}

static void
revert_sampling_and_culling( GLUnurbsObj *nobj )
{
	if(nobj->auto_load_matrix==GL_FALSE)
	{
		glMatrixMode(GL_MODELVIEW);
		glPopMatrix();
		glMatrixMode(GL_PROJECTION);
		glPopMatrix();
		glPopAttrib();
	}
}

GLenum
glu_do_sampling_3D( GLUnurbsObj *nobj, new_ctrl_type *new_ctrl,
	GLint **sfactors, GLint **tfactors)
{
	GLint	dim;
	GLenum	err;

	*sfactors=NULL;
	*tfactors=NULL;
	dim=nobj->surface.geom.dim;
	set_sampling_and_culling(nobj);
	if((err=calc_sampling_3D(new_ctrl,nobj->sampling_tolerance,dim,
		nobj->surface.geom.sorder,nobj->surface.geom.torder,
		sfactors,tfactors))==GLU_ERROR)
	{
		revert_sampling_and_culling(nobj);
		call_user_error(nobj,err);
		return GLU_ERROR;
	}
	revert_sampling_and_culling(nobj);
	return GLU_NO_ERROR;
}

GLenum
glu_do_sampling_uv( GLUnurbsObj *nobj, new_ctrl_type *new_ctrl,
	GLint **sfactors, GLint **tfactors)
{
	GLint	s_cnt, t_cnt, i;
	GLint	u_steps, v_steps;

	s_cnt = new_ctrl->s_bezier_cnt;
	t_cnt = new_ctrl->t_bezier_cnt;
	*sfactors=NULL;
	*tfactors=NULL;
	if((*sfactors=(GLint *)malloc(sizeof(GLint)*s_cnt*3))
			==NULL)
	{
		return GLU_OUT_OF_MEMORY;
	}
	if((*tfactors=(GLint *)malloc(sizeof(GLint)*t_cnt*3))
			==NULL)
	{
		free(*sfactors);
		return GLU_OUT_OF_MEMORY;
	}
	u_steps = nobj->u_step;
	v_steps = nobj->v_step;
	for(i=0; i<s_cnt; i++)
	{
		*(*sfactors+i*3) = u_steps;
		*(*sfactors+i*3+1) = u_steps;
		*(*sfactors+i*3+2) = u_steps;
	}
	for(i=0; i<t_cnt; i++)
	{
		*(*tfactors+i*3) = v_steps;
		*(*tfactors+i*3+1) = v_steps;
		*(*tfactors+i*3+2) = v_steps;
	}
	return GLU_NO_ERROR;
}

GLenum
glu_do_sampling_param_3D( GLUnurbsObj *nobj, new_ctrl_type *new_ctrl,
	GLint **sfactors, GLint **tfactors)
{
	GLint	dim;
	GLenum	err;

	*sfactors=NULL;
	*tfactors=NULL;
	dim=nobj->surface.geom.dim;
	set_sampling_and_culling(nobj);
	if((err=calc_sampling_param_3D(new_ctrl,nobj->parametric_tolerance,dim,
		nobj->surface.geom.sorder,nobj->surface.geom.torder,
		sfactors,tfactors))==GLU_ERROR)
	{
		revert_sampling_and_culling(nobj);
		call_user_error(nobj,err);
		return GLU_ERROR;
	}
	revert_sampling_and_culling(nobj);
	return GLU_NO_ERROR;
}

GLenum
glu_do_sampling_2D( GLUnurbsObj *nobj, GLfloat *ctrl, GLint cnt, GLint order,
	GLint dim, GLint **factors)
{
	GLenum err;

	set_sampling_and_culling(nobj);
	err=calc_sampling_2D(ctrl,cnt,order,nobj->sampling_tolerance,dim,
		factors);
	revert_sampling_and_culling(nobj);
	return err;
}


GLenum
glu_do_sampling_u( GLUnurbsObj *nobj, GLfloat *ctrl, GLint cnt, GLint order,
	GLint dim, GLint **factors)
{
	GLint	i;
	GLint	u_steps;

	cnt /= order;
	if((*factors=(GLint *)malloc(sizeof(GLint)*cnt))
			==NULL)
	{
		return GLU_OUT_OF_MEMORY;
	}
	u_steps = nobj->u_step;
	for(i=0; i<cnt; i++)
		(*factors)[i] = u_steps;
	return GLU_NO_ERROR;
}

GLenum
glu_do_sampling_param_2D( GLUnurbsObj *nobj, GLfloat *ctrl, GLint cnt,
	GLint order, GLint dim, GLint **factors)
{
	GLint	i;
	GLint	u_steps;
	GLfloat	tolerance;

	set_sampling_and_culling(nobj);
	tolerance = nobj->parametric_tolerance;
	cnt /= order;
	if((*factors=(GLint *)malloc(sizeof(GLint)*cnt))
			==NULL)
	{
		revert_sampling_and_culling(nobj);
		return GLU_OUT_OF_MEMORY;
	}
	u_steps = nobj->u_step;
	for(i=0; i<cnt; i++)
	{
		(*factors)[i] = calc_parametric_factor(ctrl,order,0,
			dim,tolerance,dim);

	}
	revert_sampling_and_culling(nobj);
	return GLU_NO_ERROR;
}

GLenum
glu_do_sampling_crv( GLUnurbsObj *nobj, GLfloat *ctrl, GLint cnt, GLint order,
	GLint dim, GLint **factors)
{
	GLenum err;

	*factors=NULL;
	switch(nobj->sampling_method)
	{
		case GLU_PATH_LENGTH:
			if((err=glu_do_sampling_2D(nobj,ctrl,cnt,order,dim,factors))!=
				GLU_NO_ERROR)
			{
				call_user_error(nobj,err);
				return GLU_ERROR;
			}
			break;
		case GLU_DOMAIN_DISTANCE:
			if((err=glu_do_sampling_u(nobj,ctrl,cnt,order,dim,factors))!=
				GLU_NO_ERROR)
			{
				call_user_error(nobj,err);
				return GLU_ERROR;
			}
			break;
		case GLU_PARAMETRIC_ERROR:
			if((err=glu_do_sampling_param_2D(nobj,ctrl,cnt,order,dim,factors))!=
				GLU_NO_ERROR)
			{
				call_user_error(nobj,err);
				return GLU_ERROR;
			}
			break;
		default:
			abort();
	}

	return GLU_NO_ERROR;
}

/* TODO - i don't like this culling - this one just tests if at least one */
/* ctrl point lies within the viewport . Also the point_in_viewport() */
/* should be included in the fnctions for efficiency reasons */

static GLboolean
point_in_viewport(GLfloat *pt, GLint dim)
{
	GLdouble model[16],proj[16];
	GLint viewport[4];
	GLdouble x,y,z,w,winx,winy,winz;

	glGetDoublev(GL_MODELVIEW_MATRIX,model);
	glGetDoublev(GL_PROJECTION_MATRIX,proj);
	glGetIntegerv(GL_VIEWPORT,viewport);
	if(dim==3)
	{
		x=(GLdouble)pt[0];
		y=(GLdouble)pt[1];
		z=(GLdouble)pt[2];
		gluProject(x,y,z,model,proj,viewport,&winx,&winy,&winz);
	}
	else
	{
		w=(GLdouble)pt[3];
		x=(GLdouble)pt[0]/w;
		y=(GLdouble)pt[1]/w;
		z=(GLdouble)pt[2]/w;
		gluProject(x,y,z,model,proj,viewport,&winx,&winy,&winz);
	}
	if((GLint)winx >= viewport[0] && (GLint)winx < viewport[2] &&
			(GLint)winy >= viewport[1] && (GLint)winy < viewport[3])
		return GL_TRUE;
	return GL_FALSE;
}

GLboolean
fine_culling_test_3D(GLUnurbsObj *nobj,GLfloat *pts,GLint s_cnt,GLint t_cnt,
	GLint s_stride,GLint t_stride, GLint dim)
{
	GLint 	i,j;

	if(nobj->culling==GL_FALSE)
		return GL_FALSE;
	set_sampling_and_culling(nobj);
	
	if(dim==3)
	{
		for(i=0;i<s_cnt;i++)
			for(j=0;j<t_cnt;j++)
				if(point_in_viewport(pts+i*s_stride+j*t_stride,dim))
				{
					revert_sampling_and_culling(nobj);
					return GL_FALSE;
				}
	}
	else
	{
		for(i=0;i<s_cnt;i++)
			for(j=0;j<t_cnt;j++)
				if(point_in_viewport(pts+i*s_stride+j*t_stride,dim))
				{
					revert_sampling_and_culling(nobj);
					return GL_FALSE;
				}
	}
	revert_sampling_and_culling(nobj);
	return GL_TRUE;
}

/*GLboolean
fine_culling_test_3D(GLUnurbsObj *nobj,GLfloat *pts,GLint s_cnt,GLint t_cnt,
	GLint s_stride,GLint t_stride, GLint dim)
{
	GLint		visible_cnt;
	GLfloat		feedback_buffer[5];
	GLsizei		buffer_size;
	GLint 		i,j;

	if(nobj->culling==GL_FALSE)
		return GL_FALSE;
	buffer_size=5;
	set_sampling_and_culling(nobj);
	
	glFeedbackBuffer(buffer_size,GL_2D,feedback_buffer);
	glRenderMode(GL_FEEDBACK);
	if(dim==3)
	{
		for(i=0;i<s_cnt;i++)
		{
			glBegin(GL_LINE_LOOP);
			for(j=0;j<t_cnt;j++)
				glVertex3fv(pts+i*s_stride+j*t_stride);
			glEnd();
		}
		for(j=0;j<t_cnt;j++)
		{
			glBegin(GL_LINE_LOOP);
			for(i=0;i<s_cnt;i++)
				glVertex3fv(pts+i*s_stride+j*t_stride);
			glEnd();
		}
	}
	else
	{
		for(i=0;i<s_cnt;i++)
		{
			glBegin(GL_LINE_LOOP);
			for(j=0;j<t_cnt;j++)
				glVertex4fv(pts+i*s_stride+j*t_stride);
			glEnd();
		}
		for(j=0;j<t_cnt;j++)
		{
			glBegin(GL_LINE_LOOP);
			for(i=0;i<s_cnt;i++)
				glVertex4fv(pts+i*s_stride+j*t_stride);
			glEnd();
		}
	}
	visible_cnt=glRenderMode(GL_RENDER);

	revert_sampling_and_culling(nobj);
	return (GLboolean)(visible_cnt==0);
}*/

GLboolean
fine_culling_test_2D(GLUnurbsObj *nobj,GLfloat *pts,GLint cnt,
	GLint stride, GLint dim)
{
	GLint 	i;

	if(nobj->culling==GL_FALSE)
		return GL_FALSE;
	set_sampling_and_culling(nobj);
	
	if(dim==3)
	{
		for(i=0;i<cnt;i++)
			if(point_in_viewport(pts+i*stride,dim))
			{
				revert_sampling_and_culling(nobj);
				return GL_FALSE;
			}
	}
	else
	{
		for(i=0;i<cnt;i++)
			if(point_in_viewport(pts+i*stride,dim))
			{
				revert_sampling_and_culling(nobj);
				return GL_FALSE;
			}
	}
	revert_sampling_and_culling(nobj);
	return GL_TRUE;
}

/*GLboolean
fine_culling_test_2D(GLUnurbsObj *nobj,GLfloat *pts,GLint cnt,
	GLint stride, GLint dim)
{
	GLint		visible_cnt;
	GLfloat		feedback_buffer[5];
	GLsizei		buffer_size;
	GLint 		i;

	if(nobj->culling==GL_FALSE)
		return GL_FALSE;
	buffer_size=5;
	set_sampling_and_culling(nobj);
	
	glFeedbackBuffer(buffer_size,GL_2D,feedback_buffer);
	glRenderMode(GL_FEEDBACK);
	glBegin(GL_LINE_LOOP);
	if(dim==3)
	{
		for(i=0;i<cnt;i++)
			glVertex3fv(pts+i*stride);
	}
	else
	{
		for(i=0;i<cnt;i++)
			glVertex4fv(pts+i*stride);
	}
	glEnd();
	visible_cnt=glRenderMode(GL_RENDER);

	revert_sampling_and_culling(nobj);
	return (GLboolean)(visible_cnt==0);
}*/