/*
	Utility Routines Used By the Previewer.
	V3Dot, V3Cross, V3Normalize, V3Length,ViewTransform,MatMul taken from GraphicsGems.
-----
Sumant
15 May 1991
*/
#include <stdio.h>
#include <math.h>
#include "preview.h"
extern int view_type;
extern int oneview_flag;
extern float eyepoint_distance;
static Matrix ModelMatrix={
	{1,0,0,0},
	{0,1,0,0},
	{0,0,1,0},
	{0,0,0,1}
};
static double V3SquaredLength(a)
Vector3 *a;
{
        return((a->x * a->x)+(a->y * a->y)+(a->z * a->z));
        }

static double V3Dot(a, b)
Vector3 *a, *b;
{
        return((a->x*b->x)+(a->y*b->y)+(a->z*b->z));
        }
static Vector3 *V3Cross(a, b, c)
Vector3 *a, *b, *c;
{
        c->x = (a->y*b->z) - (a->z*b->y);
        c->y = (a->z*b->x) - (a->x*b->z);
        c->z = (a->x*b->y) - (a->y*b->x);
        return(c);
        }
static double V3Length(a)
Vector3 *a;
{
        return(sqrt(V3SquaredLength(a)));
        }

Vector3 *V3Normalize(v,name)
Vector3 *v;
char *name;
{
double len = V3Length(v);
        if (len != 0.0) { v->x /= len;  v->y /= len; v->z /= len; }
        else fprintf(stderr,"V3Normalize : Zero Length Vector in %s.",name);
        return(v);
        }
static void dump_matrix(T)
Matrix T;
{
	int i,j;
	for (i=0;i<4;i++){
		for(j=0;j<4;j++)
			fprintf(stderr,"%g\t",T[i][j]);
		fprintf(stderr,"\n");
	}
}
static void ViewTransform( VRP, EV, UP, T)
Point3 VRP, EV, UP;
Matrix T;
{
	Vector3	U, V, N;
	float	dot;

	/*
	 *  N = EV  and normalize  N
	 */
	N = EV;
	V3Normalize(&N,"ViewTransform : Normal Calc.");

	/*
	 * V = UP
	 * Make vector  V  orthogonal to  N  and normalize  V
	 */
	V=UP;
	dot = V3Dot(&V,&N);
	V.x -= dot * N.x; V.y -= dot * N.y; V.z -= dot * N.z;
	V3Normalize(&V,"ViewTransform : UP Vector Calc.");


	/*
	 * Compute vector  U = V x N  (cross product)
	 */
	V3Cross(&V,&N,&U);

	/*
	 * Write the vectors U, V, and N as the first three rows of the
	 *       first, second, and third columns of  ViewMatrix, respectively
	 */
	T[0][0] = U.x;		/* column 1 , vector U */
	T[1][0] = U.y;
	T[2][0] = U.z;
	T[0][1] = V.x;		/* column 2 , vector V */
	T[1][1] = V.y;
	T[2][1] = V.z;
	T[0][2] = N.x;		/* column 3 , vector N */
	T[1][2] = N.y;
	T[2][2] = N.z;

	/*
	 * Compute the fourth row of  ViewMatrix  to include the translation of
	 *       VRP  to the origin
	 */
	T[3][0] = - U.x * VRP.x - U.y * VRP.y - U.z * VRP.z;
	T[3][1] = - V.x * VRP.x - V.y * VRP.y - V.z * VRP.z;
	T[3][2] = - N.x * VRP.x - N.y * VRP.y - N.z * VRP.z;

	/* Fourth Column is constant */
	T[0][3] = T[1][3] = T[2][3] = 0.0; T[3][3] = 1.0;
	return;
}
/* multiply together matrices c = ab */
/* note that c must not point to either of the input matrices */
static void MatMul(a, b, c)
Matrix a, b, c;
{
int i, j, k;
        for (i=0; i<4; i++) {
                for (j=0; j<4; j++) {
                        c[i][j] = 0;
                        for (k=0; k<4; k++) c[i][j] +=
                                a[i][k] * b[k][j];
                }
         }
}


static void loadmatrix(mat)
Matrix mat;
{
	int i;
	float *s=(float *)mat;
	float *t=(float *)ModelMatrix;
	for (i = 0; i < 16; i++) *t++ = *s++;
}

void identmatrix(a)
	Matrix 	 a;
{
	int i;
	register float	*p=(float *)a;

	for (i=0; i<16; i++) *p++ = 0.0;

	a[0][0] = a[1][1] = a[2][2] = a[3][3] = 1.0;
}
static void multvector(v, a, b)
Vector	v, a;
Matrix	b;
{
	int i,j;
	for(i=0;i<4;i++) for(j=0,v[i]=0.;j<4;j++) v[i] += a[j]*b[j][i];
}

#define MIN(a,b)        (((a)<(b))?(a):(b))
#define MAX(a,b)        (((a)>(b))?(a):(b))
#define check {\
	multvector(v,a,m);\
	if(v[3] ==0) error("Division By Zero","get_min_max");\
	x = v[0]/v[3]; y = v[1]/v[3];\
	*min_x=MIN(*min_x,x); *max_x=MAX(*max_x,x);\
	*min_y=MIN(*min_y,y); *max_y=MAX(*max_y,y);\
}
static void get_min_max(m,left, right, front, back, bottom,top,min_x,min_y,max_x,max_y)
Matrix m;
float left, right, front, back, bottom,top;
float *min_x,*min_y,*max_x,*max_y;
{
	Vector a,v;
	float x,y;

	*min_x = *min_y = Large; 
	*max_x = *max_y = -Large;

	a[0]=left;a[1]=front;a[2]=bottom;a[3]=1.0; check
	a[0]=right;a[1]=front;a[2]=bottom;a[3]=1.0; check
	a[0]=left;a[1]=back;a[2]=bottom;a[3]=1.0; check
	a[0]=right;a[1]=back;a[2]=bottom;a[3]=1.0; check
	a[0]=left;a[1]=front;a[2]=top;a[3]=1.0; check
	a[0]=right;a[1]=front;a[2]=top;a[3]=1.0; check
	a[0]=left;a[1]=back;a[2]=top;a[3]=1.0; check
	a[0]=right;a[1]=back;a[2]=top;a[3]=1.0; check
}
#undef check

static Matrix ViewMatrix[Maxviews];
static int Vx1[MaxviewsX2],Vy1[MaxviewsX2],Vx2[MaxviewsX2],Vy2[MaxviewsX2];
static float Wx1[Maxviews],Wy1[Maxviews],Wx2[Maxviews],Wy2[Maxviews];

void set_screen_transformation(left, right, front, back, bottom,top)
float left, right, front, back, bottom, top;
{
	extern int window_x_size,window_y_size;
	Point3 VRP,EP,UP;
	Vector box[8];
	float min_x,min_y,max_x,max_y;
	float deltax,deltay,deltaz,maxdelta;
	if ((deltax=(right - left))== 0.0)
		error("Left clipping plane same as right one.",
			"utils.c");

	if ((deltaz=(top - bottom))== 0.0)
		error("Bottom clipping plane same as top one.",
			"utils.c");

	if ((deltay=(back - front))== 0.0)
		error("Front clipping plane same as back one.",
			"utils.c");
	maxdelta = deltax;
	if (deltay > maxdelta) maxdelta = deltay;
	if (deltaz > maxdelta) maxdelta = deltaz;
	left -= (maxdelta-deltax)/2.;
	right += (maxdelta-deltax)/2.;
	bottom -= (maxdelta-deltaz)/2.;
	top += (maxdelta-deltaz)/2.;
	front -= (maxdelta-deltay)/2.;
	back += (maxdelta-deltay)/2.;
	identmatrix(ModelMatrix);
	VRP.x = 0.5*(left+right); VRP.y = 0.5*(front+back); VRP.z = 0.5*(bottom+top);
	/* Elevation */
	UP.x = UP.y = 0.0; UP.z = 1.0;
	EP.y = -1.0; EP.x = EP.z = 0.0;
	ViewTransform(VRP,EP,UP,ViewMatrix[0]);
	get_min_max(ViewMatrix[0],left, right, front, back, bottom, top,
		&min_x,&min_y,&max_x,&max_y);
        Wx1[0] = min_x; Wy1[0] = min_y; Wx2[0] = max_x; Wy2[0] = max_y;
	Vx1[0] = 1; Vy1[0] = window_y_size/2-2; Vx2[0] = window_x_size/2-2; Vy2[0] = 1;
		/* Side */
	EP.x =1.0; EP.z = EP.y = 0.0;
	ViewTransform(VRP,EP,UP,ViewMatrix[1]);
	get_min_max(ViewMatrix[1],left, right, front, back, bottom, top,
		&min_x,&min_y,&max_x,&max_y);
        Wx1[1] = min_x; Wy1[1] = min_y; Wx2[1] = max_x; Wy2[1] = max_y;
	Vx1[1] = window_x_size/2+1; Vy1[1] = window_y_size/2-2;
	Vx2[1] = window_x_size-2; Vy2[1] = 1;
		/* Plan */
	UP.z = UP.x = 0.0; UP.y = 1.0;
	EP.z = 1.0; EP.y = EP.x = 0.0;
	ViewTransform(VRP,EP,UP,ViewMatrix[2]);
	get_min_max(ViewMatrix[2],left, right, front, back, bottom, top,
		&min_x,&min_y,&max_x,&max_y);
        Wx1[2] = min_x; Wy1[2] = min_y; Wx2[2] = max_x; Wy2[2] = max_y;
	Vx1[2] = 1; Vy1[2] = window_y_size-2;
	Vx2[2] = window_x_size/2-2; Vy2[2] = window_y_size/2+1;
		/* Isometry */
	UP.x = UP.y = 0.0; UP.z = 1.0;
	EP.z = 1.0; EP.y = 1.0;EP.x = 1.0;
	ViewTransform(VRP,EP,UP,ViewMatrix[3]);
	ViewTransform(VRP,EP,UP,ViewMatrix[4]);
	get_min_max(ViewMatrix[3],left, right, front, back, bottom, top,
		&min_x,&min_y,&max_x,&max_y);
        Wx1[3] = Wx1[4] = min_x; Wy1[3] = Wy1[4] = min_y;
	Wx2[3] = Wx2[4] = max_x; Wy2[3] = Wy2[4] = max_y;
	Vx1[3] = Vx1[4] = window_x_size/2+1;
		Vy1[3] = Vy1[4] = window_y_size-2;
	Vx2[3] = Vx2[4] = window_x_size-2;
		Vy2[3] = Vy2[4] = window_y_size/2+1;

	Vx1[5] = Vx1[6] = Vx1[7] = Vx1[8] = Vx1[9] = 1;
	Vy1[5] = Vy1[6] = Vy1[7] = Vy1[8] = Vy1[9] = window_y_size-2;
	Vx2[5] = Vx2[6] = Vx2[7] = Vx2[8] = Vx2[9] = window_x_size-2;
	Vy2[5] = Vy2[6] = Vy2[7] = Vy2[8] = Vy2[9] = 1;
}

compute_perspective_matrix(left, right, front, back, bottom, top)
float left, right, front, back, bottom, top;
{
	extern int window_x_size,window_y_size;
	Matrix m,m1;
	Point3 EP,VRP,UP;
	float min_x,min_y,max_x,max_y;

	fprintf(stderr,"Screen Extents are <%g %g %g> to <%g %g %g>.\n",
		bound_min_x,bound_min_y,bound_min_z,
		bound_max_x,bound_max_y,bound_max_z);
	fprintf(stderr,"Give the view params:\n");
	fprintf(stderr,"Eye Point : ");
	scanf("%lf%lf%lf",&(EP.x),&(EP.y),&(EP.z));
	fprintf(stderr,"View Reference Point : ");
	scanf("%lf%lf%lf",&(VRP.x),&(VRP.y),&(VRP.z));
	fprintf(stderr,"Up Point : ");
	scanf("%lf%lf%lf",&(UP.x),&(UP.y),&(UP.z));

	ViewTransform(VRP,EP,UP,ViewMatrix[4]);
	/* Project on a view plane 1 unit away */
	identmatrix(m);
	m[2][2]=m[3][3]=0.; m[2][3] = 1.;
	MatMul(m1,m,ViewMatrix[4]);
	get_min_max(ViewMatrix[4],left, right, front, back, bottom, top,
		&min_x,&min_y,&max_x,&max_y);
       	Wx1[4] = min_x; Wy1[4] = min_y; Wx2[4] = max_x; Wy2[4] = max_y;
	Vx1[4] = window_x_size/2+1; Vy1[4] = window_y_size-2;
	Vx2[4] = window_x_size-2; Vy2[4] = window_y_size/2+1;
}

static int get_screen_xy(n,deltan,x,y,z,sx,sy)
int n,deltan;
float x,y,z;
int *sx,*sy;
{
	Vector v,v1,vv;
	v[0]=x; v[1]=y; v[2]=z; v[3]=1.0;

	/* Model Transform the Point */
	multvector(v1,v,ModelMatrix);
	/* View Transform the point */
	multvector(vv,v1,ViewMatrix[n]);
	if (vv[3] == 0.0){
		warn("get_screen_xy : Homogeneous Coordinate.", "utils.c");
		return(1);
	}
	/* Window to View port Map */
	*sx = Vx1[n+deltan] + (vv[0]/vv[3] - Wx1[n])*
				(Vx2[n+deltan]-Vx1[n+deltan])/(Wx2[n]-Wx1[n]);
	*sy = Vy1[n+deltan] + (vv[1]/vv[3] - Wy1[n])*
				(Vy2[n+deltan]-Vy1[n+deltan])/(Wy2[n]-Wy1[n]);
	return(0);
}
static int previous_screen_x[4]={0,320,0,320},
	   previous_screen_y[4]={319,319,639,639};

#define Code(x,y,c) {\
		c = 0x00;\
		if (x < min_x) c |= left;\
		else if (x > max_x) c |= right;\
		if (y > min_y) c |= bottom;\
		else if (y < max_y) c |= top;\
	}
void clip_and_draw(x1,y1,x2,y2,min_x,min_y,max_x,max_y)
int x1,y1,x2,y2,min_x,max_y,max_x,min_y;
{
	unsigned char left=0x01,right=0x02,bottom=0x04,top=0x08;
	int x,y;
	unsigned char c,c1,c2;
	Code(x1,y1,c1); Code(x2,y2,c2);
	while (c1 || c2){
		if (c1 & c2) return;
		c = c1; if (!c) c = c2;
		if (c & left){
			y = y1 + ((float)(y2 - y1)*(min_x-x1))/(x2-x1);
			x = min_x;
		}
		else if (c & right){
			y = y1 + ((float)(y2 - y1)*(max_x-x1))/(x2-x1);
			x = max_x;
		}	
		else if (c & bottom){
			x = x1 + ((float)(x2 - x1)*(min_y-y1))/(y2-y1);
			y = min_y;
		}	
		else if (c & top){
			x = x1 + ((float)(x2 - x1)*(max_y-y1))/(y2-y1);
			y = max_y;
		}	
		if (c == c1) {
			x1 = x; y1 = y; Code(x,y,c1);
		}
		else{
			x2 = x; y2 = y; Code(x,y,c2);
		}
	}
	Xdraw(x1,y1,x2,y2);
}
#undef Code
void draw(x,y,z)
float x,y,z;
{
	int i;
	int sx,sy;
	if (view_type == ONE){
		if (get_screen_xy(oneview_flag,Maxviews,x,y,z,&sx,&sy))
			return;
		Xdraw(previous_screen_x[0],previous_screen_y[0],sx,sy);
		previous_screen_x[0]=sx; previous_screen_y[0]=sy;
	}
	else for(i=0;i<4;i++){
		if (get_screen_xy(i,0,x,y,z,&sx,&sy))
			continue;
		clip_and_draw(
			previous_screen_x[i],previous_screen_y[i],sx,sy,
			Vx1[i],Vy1[i],Vx2[i],Vy2[i]);
		previous_screen_x[i]=sx; previous_screen_y[i]=sy;
	}
}
void move(x,y,z)
float x,y,z;
{
	int i;
	int sx,sy;
	if (view_type == ONE){
		if (get_screen_xy(oneview_flag,Maxviews,x,y,z,&sx,&sy))
			return;
		previous_screen_x[0]=sx; previous_screen_y[0]=sy;
	}
	else for(i=0;i<4;i++){
		if (get_screen_xy(i,0,x,y,z,&sx,&sy))
			continue;
		previous_screen_x[i]=sx; previous_screen_y[i]=sy;
	}
}
#define MaxStackSize 100
static float stack[MaxStackSize];
static int topofstack=0;
void pushmatrix()
{
	if ((topofstack+16)< MaxStackSize){
		int i;
		float *m=(float *)ModelMatrix;
		for(i=0;i<16;i++) stack[topofstack++] = *m++;
	}
	else error("pushmatrix : Stack overflow","utils.c");
}
void popmatrix()
{
	if (topofstack < 16) error("popmatrix : Empty Stack","utils.c");
	else {
		int i;
		float *m=(float *)ModelMatrix;
		for(i=0,m+=15;i<16;i++) *m--=stack[--topofstack];
	}
}

void translate(x,y,z)
float x,y,z;
{
	Matrix m,m1;
	identmatrix(m);
	m[3][0] = x; m[3][1] = y; m[3][2] = z;
	MatMul(m,ModelMatrix,m1);
	loadmatrix(m1);
}

void rotate(r,axis)
float r;
char axis;
{
	float sintheta=sin((double)(DTOR * r));
	float costheta=cos((double)(DTOR * r));
	float tmp;
	Matrix m,m1;
	identmatrix(m);
	switch(axis){
	case 'x':
	case 'X': 
		  m[1][1]=m[2][2]= costheta;
		  m[1][2]= sintheta;
		  m[2][1]= -sintheta;
		break;
	case 'y':
	case 'Y': 
		  m[0][0] = m[2][2] = costheta;
		  m[0][2] = -sintheta;
		  m[2][0] = sintheta;
		break;
	case 'z':
	case 'Z': 
		  m[0][0] = m[1][1] = costheta;
		  m[0][1] = sintheta;
		  m[1][0] = -sintheta;
		break;
	default: error("rotate : Illegal axis.","utils.c"); break;
	}
	MatMul(m,ModelMatrix,m1);
	loadmatrix(m1);
}

void multmatrix(m)
Matrix m;
{
	Matrix m1;
	MatMul(m,ModelMatrix,m1);
	loadmatrix(m1);
}

void circle(radius)
float radius;
{
	extern int circle_precision;
	int i;
	float   cx, cy, dx, dy;
	float   angle=2.0 * PI /circle_precision;
        float cosine = cos((double)angle);
        float sine = sin((double)angle);
        cx = radius;
        cy = 0.0;
	move(cx,cy,0.0);
        for (i = 0; i < (circle_precision - 1); i++) {
                dx = cx;
                dy = cy;
                cx = dx * cosine - dy * sine;
                cy = dx * sine + dy * cosine;
                draw(cx, cy,0.0);
	}
	draw(radius,0.0,0.0);
}