//////////////////////////////////////////////////////////////////////
//
//                             Pixie
//
// Copyright © 1999 - 2003, Okan Arikan
//
// Contact: okan@cs.berkeley.edu
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU 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
// General Public License for more details.
// 
// You should have received a copy of the GNU General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//
///////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////
//
//  File				:	subdivision.cpp
//  Classes				:	CSubdivision
//  Description			:	Implements a subdivision surface
//
////////////////////////////////////////////////////////////////////////
#include <math.h>

#include "subdivision.h"
#include "memory.h"
#include "patches.h"
#include "renderer.h"
#include "stats.h"
#include "subdivisionData.h"			// Generated by "precompute"
#include "shading.h"

// This macro is used to fix the degenerate normal vectors
#define	normalFix()	{																				\
	float		*Ng	=	varying[VARIABLE_NG];														\
	int			i;																					\
																									\
	for (i=numVertices;i>0;i--,Ng+=3) {																\
		if (dotvv(Ng,Ng) == 0) {																	\
			const float	*u		=	varying[VARIABLE_U];											\
			const float	*v		=	varying[VARIABLE_V];											\
			const float	*cNg	=	varying[VARIABLE_NG];											\
			const float	cu		=	u[numVertices-i];												\
			const float	cv		=	v[numVertices-i];												\
			float		cd		=	C_INFINITY;														\
			const float	*closest=   Ng;																\
			int			j;																			\
																									\
			for (j=numVertices;j>0;j--,cNg+=3,u++,v++) {											\
				if (dotvv(cNg,cNg) > 0) {															\
					const float	du	=	cu - u[0];													\
					const float	dv	=	cv - v[0];													\
					float		d;																	\
																									\
					d	=	du*du + dv*dv;															\
					if (d < cd) {																	\
						cd		=	d;																\
						closest	=	cNg;															\
					}																				\
				}																					\
			}																						\
																									\
			movvv(Ng,closest);																		\
		}																							\
	}																								\
}

///////////////////////////////////////////////////////////////////////
// Class				:	CSubdivision
// Method				:	CSubdivision
// Description			:	Ctor
// Return Value			:	-
// Comments				:
// Date last edited		:	10/12/2002
CSubdivision::CSubdivision(CAttributes *a,CXform *x,CVertexData *var,CParameter *p,int N,float uOrg,float vOrg,float uMult,float vMult,double *vertex) : CSurface(a,x) {
	const int	K		=	2*N+8;

	stats.numGprims++;
	stats.gprimMemory	+=	sizeof(CSubdivision);

	vertexData			=	var;
	vertexData->attach();
	parameters			=	p;

	assert(vertex != NULL);

	this->N				=	N;
	this->uOrg			=	uOrg;
	this->vOrg			=	vOrg;
	this->uMult			=	uMult;
	this->vMult			=	vMult;

	initv(bmin,C_INFINITY,C_INFINITY,C_INFINITY);
	initv(bmax,-C_INFINITY,-C_INFINITY,-C_INFINITY);

	if (vertexData->moving == FALSE) {
		this->vertex		=	new float[K*vertexData->vertexSize];
		stats.gprimMemory	+=	K*vertexData->vertexSize*sizeof(float);

		projectVertices(this->vertex							,vertex,	0);
	} else {
		this->vertex		=	new float[K*vertexData->vertexSize*2];
		stats.gprimMemory	+=	K*vertexData->vertexSize*2*sizeof(float);

		projectVertices(this->vertex							,vertex,	0);
		projectVertices(this->vertex+K*vertexData->vertexSize	,vertex,	vertexData->vertexSize);
	}

	makeBound(bmin,bmax);
}

///////////////////////////////////////////////////////////////////////
// Class				:	CSubdivision
// Method				:	~CSubdivision
// Description			:	Dtor
// Return Value			:	-
// Comments				:
// Date last edited		:	10/12/2002
CSubdivision::~CSubdivision() {
	const int	K	=	2*N+8;

	delete [] vertex;	

	if (vertexData->moving)	stats.gprimMemory	-=	K*vertexData->vertexSize*2*sizeof(float);
	else					stats.gprimMemory	-=	K*vertexData->vertexSize*sizeof(float);

	if (parameters != NULL)	delete parameters;

	vertexData->detach();

	stats.numGprims--;
	stats.gprimMemory	-=	sizeof(CSubdivision);
}


///////////////////////////////////////////////////////////////////////
// Class				:	CSubdivision
// Method				:	projectVertices
// Description			:	Project a set of vertices into the eigen space of the patch
// Return Value			:	The projected points
// Comments				:
// Date last edited		:	10/12/2002
void	CSubdivision::projectVertices(float *fvertex,double *vertexData,int disp) {
	int			i;
	double		*cVertex;
	const float	*evecs			=	basisData[N].evecs;
	int			K				=	2*N+8;
	const int	vertexSize		=	this->vertexData->vertexSize;
	const int	vs				=	(this->vertexData->moving ? vertexSize*2 : vertexSize);
	double		*vertex			=	(double *) alloca(K*vertexSize*sizeof(double));

	// Clear the new vertex data
	for (i=0;i<K*vertexSize;i++) {
		vertex[i]	=	0;
	}

	// Project the points to the eigen space
	for (cVertex=vertex,i=0;i<K;i++,cVertex+=vertexSize) {
		int		j,k;

		for (j=0;j<K;j++) {	
			double	val	=	evecs[i+j*K];

			for (k=0;k<vertexSize;k++) {
				cVertex[k]	+=	val*vertexData[j*vs+k+disp];
			}
		}
	}

	for (i=0;i<vertexSize*K;i++) {
		fvertex[i]		=	(float) vertex[i];
	}

	for (i=K;i>0;i--,vertexData+=vs) {
		vector	tmp;

		tmp[0]	=	(float) vertexData[0];
		tmp[1]	=	(float) vertexData[1];
		tmp[2]	=	(float) vertexData[2];

		addBox(bmin,bmax,tmp);
	}
}

///////////////////////////////////////////////////////////////////////
// Class				:	CSubdivision
// Method				:	bound
// Description			:	See object.h
// Return Value			:	-
// Comments				:
// Date last edited		:	10/12/2002
void		CSubdivision::bound(float *bmi,float *bma) const {
	movvv(bmi,bmin);
	movvv(bma,bmax);
}


///////////////////////////////////////////////////////////////////////
// Class				:	CSubdivision
// Method				:	tesselate
// Description			:	See object.h
// Return Value			:	-
// Comments				:	-
// Date last edited		:	6/21/2001
void		CSubdivision::tesselate(CShadingContext *context) {
	context->tesselate2D(this);
}

///////////////////////////////////////////////////////////////////////
// Class				:	CSubdivision
// Method				:	sample
// Description			:	See object.h
// Return Value			:	-
// Comments				:
// Date last edited		:	10/12/2002
void		CSubdivision::sample(int start,int numVertices,float **varying,unsigned int &up) const {
	int					i,j,k,t;
	const float			*u						=	varying[VARIABLE_U]+start;
	const float			*v						=	varying[VARIABLE_V]+start;
	const int			vertexSize				=	this->vertexData->vertexSize;
	const CEigenBasis	*cBasis					=	&basisData[N];
	const int			K						=	2*N+8;
	float				*vertexData;
	int					vertexDataStep;

	memBegin();

	if (this->vertexData->moving == 0) {
		vertexData		=	vertex;			// No need for interpolation
		vertexDataStep	=	0;
	} else {									
		if (up & PARAMETER_BEGIN_SAMPLE) {
			vertexData		=	vertex;		// No need for interpolation
			vertexDataStep	=	0;
		} else if (up & PARAMETER_END_SAMPLE) {
			vertexData		=	vertex	+	K*vertexSize;		// No need for interpolation
			vertexDataStep	=	0;
		} else {
												// Interpolate the vertex data in advance
			float			*interpolate;
			const float		*time		=	varying[VARIABLE_TIME] + start;
			const float		*vertex0	=	vertex;
			const float		*vertex1	=	vertex + vertexSize*K;

			vertexData				=	(float *) ralloc(numVertices*K*vertexSize*sizeof(float));
			vertexDataStep			=	K*vertexSize;

			interpolate				=	vertexData;

			for (i=numVertices;i>0;i--) {
				const float ctime		=	*time++;

				for (j=0;j<vertexDataStep;j++) {
					*interpolate++	=	vertex0[j]*(1-ctime) + vertex1[j]*ctime;
				}
			}
		}
	}

	// Note: In order to make sure the normals point the right direction,
	//       I swapped the u and v, which also means swapping dPdu and dPdv
	//       and the order of the cross product for the normal vector.



	{	// Do the vertices
		float	*intr	=	(float *) ralloc(numVertices*vertexSize*sizeof(float));
		float	*tmp	=	intr;
		//float	*dPdu	=	varying[VARIABLE_DPDU] + start*3;
		//float	*dPdv	=	varying[VARIABLE_DPDV] + start*3;
		float	*dPdu	=	varying[VARIABLE_DPDV] + start*3;
		float	*dPdv	=	varying[VARIABLE_DPDU] + start*3;
		float	*N		=	varying[VARIABLE_NG] + start*3;

		// Init the interpolate
		for (i=0;i<numVertices*vertexSize;i++) {
			intr[i]	=	0;
		}

		// FIXME: This definition of real is temporary. After profiling, remove it
#define	real	float

		for (i=0;i<numVertices;i++) {
			real	cu	=	v[i];
			real	cv	=	u[i];
			real	pow2;
			int		n;
			real	u2;
			real	u3;
			real	v2;
			real	v3;
			real	normalScale;

			if ((cu == 0) && (cv == 0)) {
				n		=	10;
			} else {
				n		=		(int)  floor(min(-log(cu),-log(cv))/log(2.0))+1;
				if (n <= 0)	n	=	1;	// Need at least one subdivision
			}

			pow2	=		(real) (1 << (n-1));
			cu		*=		pow2;
			cv		*=		pow2;
			if (cv < 0.5) {
				k	=	0;
				cu	=	2*cu-1;
				cv	=	2*cv;
			} else if (cu < 0.5) {
				k	=	2;
				cu	=	2*cu;
				cv	=	2*cv-1;
			} else {
				k	=	1;
				cu	=	2*cu-1;
				cv	=	2*cv-1;
			}
			initv(dPdu,0,0,0);
			initv(dPdv,0,0,0);
			initv(N,0,0,0);

			u2			=	cu*cu;
			u3			=	u2*cu;
			v2			=	cv*cv;
			v3			=	v2*cv;
			normalScale	=	pow2*2;

			for (j=0;j<K;j++) {
				float	*Psrc;
				real	coef,ducoef,dvcoef;
				real	p;
				real	t0,t1,t2,t3;
				float	*coefs	=	cBasis->basis[k]+j*16;

				t0				=	coefs[0]*v3		+ coefs[1]*v2	+ coefs[2]*cv		+ coefs[3];
				t1				=	coefs[4]*v3		+ coefs[5]*v2	+ coefs[6]*cv		+ coefs[7];
				t2				=	coefs[8]*v3		+ coefs[9]*v2	+ coefs[10]*cv		+ coefs[11];
				t3				=	coefs[12]*v3	+ coefs[13]*v2	+ coefs[14]*cv		+ coefs[15];

				coef			=	u3*t0 + u2*t1 + cu*t2 + t3;
				ducoef			=	3*u2*t0 + 2*cu*t1 + t2;

				t0				=	coefs[0]*3*v2		+ coefs[1]*2*cv		+ coefs[2];
				t1				=	coefs[4]*3*v2		+ coefs[5]*2*cv		+ coefs[6];
				t2				=	coefs[8]*3*v2		+ coefs[9]*2*cv		+ coefs[10];
				t3				=	coefs[12]*3*v2		+ coefs[13]*2*cv	+ coefs[14];

				dvcoef			=	u3*t0 + u2*t1 + cu*t2 + t3;

				p				=	(real) pow(cBasis->evals[j],n-1);

				coef			*=	p;
				ducoef			*=	p*normalScale;
				dvcoef			*=	p*normalScale;
			

				Psrc			=	vertexData + j*vertexSize;

				for (t=0;t<vertexSize;t++) {
					tmp[t]		+=	(float) (Psrc[t]*coef);
				}

				dPdv[0]			+=	(float) (Psrc[0]*dvcoef);
				dPdv[1]			+=	(float) (Psrc[1]*dvcoef);
				dPdv[2]			+=	(float) (Psrc[2]*dvcoef);
				dPdu[0]			+=	(float) (Psrc[0]*ducoef);
				dPdu[1]			+=	(float) (Psrc[1]*ducoef);
				dPdu[2]			+=	(float) (Psrc[2]*ducoef);
			}

			//crossvv(N,dPdu,dPdv);
			crossvv(N,dPdv,dPdu);

			tmp			+=	vertexSize;
			dPdu		+=	3;
			dPdv		+=	3;
			N			+=	3;
			vertexData	+=	vertexDataStep;
		}

#undef real

		this->vertexData->dispatch(intr,start,numVertices,varying);
	}

	// Fix the degenerate normals
	//normalFix();

	up	&=	~(PARAMETER_P | PARAMETER_DPDU | PARAMETER_DPDV | PARAMETER_NG | this->vertexData->parameters);

	memEnd();
}

///////////////////////////////////////////////////////////////////////
// Class				:	CSubdivision
// Method				:	interpolate
// Description			:	See object.h
// Return Value			:	-
// Comments				:
// Date last edited		:	10/12/2002
void		CSubdivision::interpolate(int numVertices,float **varying) const {
	// Correct the parametric range of the primitive
	// do it first so we interpolate varyings on larger patch
	if ((uMult != 1) || (vMult != 1)) {
		float	*u,*v,*du,*dv,*dPdu,*dPdv;
		int		i;

		u		=	varying[VARIABLE_U];
		v		=	varying[VARIABLE_V];
		du		=	varying[VARIABLE_DU];
		dv		=	varying[VARIABLE_DV];
		dPdu	=	varying[VARIABLE_DPDU];
		dPdv	=	varying[VARIABLE_DPDV];

		for (i=numVertices;i>0;i--) {
			*u++	=	(*u) * uMult + uOrg;
			*v++	=	(*v) * vMult + vOrg;
			*du++	*=	uMult;
			*dv++	*=	vMult;
			mulvf(dPdu,uMult);	dPdu	+=	3;
			mulvf(dPdv,vMult);	dPdv	+=	3;
		}
	}
	
	if (parameters != NULL)	parameters->dispatch(numVertices,varying);
}