/* nurbsutl.c */

/*
 * Mesa 3-D graphics library
 * Version:  1.2
 * Copyright (C) 1995  Brian Paul  (brianp@ssec.wisc.edu)
 *
 * 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.
 */


/*
$Id: nurbsutl.c,v 1.5 1995/11/03 14:15:13 brianp Exp $

$Log: nurbsutl.c,v $
 * Revision 1.5  1995/11/03  14:15:13  brianp
 * Bogdan's November 3, 1995 updates
 *
 * Revision 1.4  1995/09/20  18:25:57  brianp
 * removed Bogdan's old email address
 *
 * Revision 1.3  1995/08/04  13:09:59  brianp
 * include gluP.h to define NULL, just in case
 *
 * Revision 1.2  1995/07/28  21:37:30  brianp
 * updates from Bogdan on July 28
 *
 * Revision 1.1  1995/07/28  14:45:10  brianp
 * Initial revision
 *
 */


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


#include <math.h>
#include <stdlib.h>
#include "nurbs.h"

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;
}

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;
	GLboolean knot_open_at_begin=GL_FALSE,knot_open_at_end=GL_FALSE;
	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;

	if(the_knot->open_at_begin)
	{
		/* knot open at beggining */
		knot_open_at_begin=GL_TRUE;
		++t_min;
	}
	if(the_knot->open_at_end)
	{
		/* knot open at end */
		knot_open_at_end=GL_TRUE;
		--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;
	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;
	m=n+the_knot->delta_nknots;
	alpha=the_knot->alpha;
	if(the_knot->open_at_begin==GL_TRUE)
	{
		k=0;
	}
	else
	{
		k=the_knot->order-1;
		m-=k;
	}
	/* is knot open at end? */
	if(the_knot->open_at_end==GL_FALSE)
		m-=the_knot->order-1;
	/* 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==3)
	{
		x=(GLdouble)pts[indx];
		y=(GLdouble)pts[indx+1];
		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];
			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;
		}
	}
	else
	{
		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;
		}
	}
	len /= tolerance;
	return (GLint)len;
}

static GLenum
calc_sampling_3D(GLfloat *ctrl, GLint ucnt, GLint vcnt, GLint uorder, GLint vorder,
	GLfloat tolerance, GLint dim, GLint **ufactors, GLint **vfactors)
{
	GLint		ufactor_cnt,vfactor_cnt;
	GLint		tmp_factor,max_factor;
	GLint		offset1,offset2;
	GLint		i,j;

	ufactor_cnt=ucnt/uorder;
	vfactor_cnt=vcnt/vorder;
	if((*ufactors=(GLint *)malloc(sizeof(GLint)*ufactor_cnt))==NULL)
	{
		return GLU_OUT_OF_MEMORY;
	}
	if((*vfactors=(GLint *)malloc(sizeof(GLint)*vfactor_cnt))==NULL)
	{
		free(ufactors);
		return GLU_OUT_OF_MEMORY;
	}
	offset1=vorder*dim;
	offset2=vcnt*dim;
	for(j=0;j<vfactor_cnt;j++)
	{
		for(i=0 , max_factor=1;i<ucnt;i++)
		{
			tmp_factor=calc_factor(ctrl,vorder,j*offset1+i*offset2,
				dim,tolerance,dim);
			if(tmp_factor>max_factor)
				max_factor=tmp_factor;
		}
		(*vfactors)[j]=max_factor;
	}
	offset1=offset2;
	offset2*=uorder;
	for(j=0;j<ufactor_cnt;j++)
	{
		for(i=0 , max_factor=1;i<vcnt;i++)
		{
			tmp_factor=calc_factor(ctrl,uorder,j*offset2+i*dim,
				offset1,tolerance,dim);
			if(tmp_factor>max_factor)
				max_factor=tmp_factor;
		}
		(*ufactors)[j]=max_factor;
	}
	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(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, GLfloat *ctrl, GLint scnt, GLint tcnt,
	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(ctrl,scnt,tcnt,nobj->surface.geom.sorder,
		nobj->surface.geom.torder,nobj->sampling_tolerance,dim,
		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;

	*factors=NULL;
	set_sampling_and_culling(nobj);
	if((err=calc_sampling_2D(ctrl,cnt,order,nobj->sampling_tolerance,dim,
		factors))==GLU_ERROR)
	{
		revert_sampling_and_culling(nobj);
		call_user_error(nobj,err);
		return GLU_ERROR;
	}
	revert_sampling_and_culling(nobj);
	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);
}*/