#include <stdio.h>
#include <string.h>
#include <ctype.h>
#include <math.h>
#include <malloc.h>
#include "GraphicsGems.h"
#include "data_structure.h"
#include "objects.h"

#include "raddecl.h"

/************** Begin object specific intersections **************/
static double intersect_polygon(p,ray_start,ray_direction,pu,pv)
Polygon *p;
Point3 *ray_start;
Vector3 *ray_direction;
double *pu,*pv;
{
	double ray_polygon();
	*pu = *pv = 0.5;
	return(
		ray_polygon(ray_start,ray_direction,
			&(p->plane_normal), p->plane_constant,
			p->dominant_axis_flag,
			p->nvertices,
			p->vertex_list)
	);
}

static double intersect_quadrilateral(q,ray_start,ray_direction,pu,pv)
Quadrilateral *q;
Point3 *ray_start;
Vector3 *ray_direction;
double *pu,*pv;
{
	double ray_quadrilateral();
	return(
		ray_quadrilateral(ray_start,ray_direction,
			&(q->plane_normal), q->plane_constant,
			q->Du0,q->Du1,q->Du2,q->Dux,q->Duy,
			q->Dv0,q->Dv1,q->Dv2,q->Dvx,q->Dvy,
			&(q->Qux),&(q->Quy),&(q->Qvx),&(q->Qvy),
			&(q->Na),&(q->Nb),&(q->Nc),
			pu,pv)
	);
}

static double intersect_sphere(s,ray_start,ray_direction,pu,pv)
Sphere *s;
Point3 *ray_start;
Vector3 *ray_direction;
double *pu,*pv;
{
	double ray_sphere();
	return(
		ray_sphere(ray_start,ray_direction,
			s->radius, &(s->center), 0, pu,pv)
	);
}

static double intersect_bubble(b,ray_start,ray_direction,pu,pv)
Sphere *b;
Point3 *ray_start;
Vector3 *ray_direction;
double *pu,*pv;
{
	double ray_sphere();
	return(
		ray_sphere(ray_start,ray_direction,
			b->radius, &(b->center), 1, pu,pv)
	);
}

static void ray_transform(s_in,d_in,m,s_out,d_out)
Point3 *s_in,*s_out;
Vector3 *d_in,*d_out;
Matrix4 *m;
{
        *d_out = transform_vector(d_in,m);
        *s_out = *s_in;
        V3MulPointByMatrix(s_out,m);
}

static double intersect_cylinder(c,ray_start,ray_direction,pu,pv)
Cylinder *c;
Point3 *ray_start;
Vector3 *ray_direction;
double *pu,*pv;
{
	double ray_cylinder();
	Point3 s;
	Vector3 d;
	ray_transform(ray_start,ray_direction,&(c->world_to_object),&s,&d);
	return(
		ray_cylinder(&s,&d,c->radius,c->height,0,pu,pv)
	);
}

static double intersect_tube(tube,ray_start,ray_direction,pu,pv)
Cylinder *tube;
Point3 *ray_start;
Vector3 *ray_direction;
double *pu,*pv;
{
	double ray_cylinder();
	Point3 s;
	Vector3 d;
	ray_transform(ray_start,ray_direction,&(tube->world_to_object),&s,&d);
	return(
		ray_cylinder(&s,&d,tube->radius,tube->height,1,pu,pv)
	);
}

static double intersect_cone(c,ray_start,ray_direction,pu,pv)
Cone *c;
Point3 *ray_start;
Vector3 *ray_direction;
double *pu,*pv;
{
	double ray_cone();
	Point3 s;
	Vector3 d;
	ray_transform(ray_start,ray_direction,&(c->world_to_object),&s,&d);
	return(
		ray_cone(&s,&d,c->distance,0,pu,pv)
	);
}

static double intersect_cup(c,ray_start,ray_direction,pu,pv)
Cone *c;
Point3 *ray_start;
Vector3 *ray_direction;
double *pu,*pv;
{
	double ray_cone();
	Point3 s;
	Vector3 d;
	ray_transform(ray_start,ray_direction,&(c->world_to_object),&s,&d);
	return(
		ray_cone(&s,&d,c->distance,1,pu,pv)
	);
}

static double intersect_ring(r,ray_start,ray_direction,pu,pv)
Ring *r;
Point3 *ray_start;
Vector3 *ray_direction;
double *pu,*pv;
{
	double ray_ring();
	Point3 s;
	Vector3 d;
	ray_transform(ray_start,ray_direction,&(r->world_to_object),&s,&d);
	return(
		ray_ring(&s,&d,r->inner_radius,r->outer_radius,pu,pv)
	);
}

static double intersect_other_geometry(dummy,ray_start,ray_direction,pu,pv)
void *dummy;
Point3 *ray_start;
Vector3 *ray_direction;
double *pu,*pv;
{
	error("This Geometry is unimplemented.");
}


/************* End object specific intersections **************/
/************* Begin Object Specific Preprocess  **************/
/*
get_view_trans_matrix(U,V,N,m)
Vector3 *U,*V,*N;
Matrix4 *m;
{
        m->element[0][0]=U->x;
        m->element[1][0]=U->y;
        m->element[2][0]=U->z;
        m->element[0][1]=V->x;
        m->element[1][1]=V->y;
        m->element[2][1]=V->z;
        m->element[0][2]=N->x;
        m->element[1][2]=N->y;
        m->element[2][2]=N->z;

	m->element[0][3]= m->element[1][3]= m->element[2][3]=
		m->element[3][0]= m->element[3][1]= m->element[3][2]=0.;
	m->element[3][3]=1.;
}
*/
#define update_scene_bound(bb,b) {\
	b->min.x = MIN(bb->min.x,b->min.x); b->max.x = MAX(bb->max.x,b->max.x);\
	b->min.y = MIN(bb->min.y,b->min.y); b->max.y = MAX(bb->max.y,b->max.y);\
	b->min.z = MIN(bb->min.z,b->min.z); b->max.z = MAX(bb->max.z,b->max.z);\
}

static void null_preprocess(grid,un,vn,spec_struc,bb,b)
Surface_Grid_structure *grid;
int un,vn;
void *spec_struc;
Box3 *bb,*b;
{
	error("No Preprocess Implemented Yet for this Object.");
}
/*
static double area (p1,p2,p3,p4)
Point3 *p1,*p2,*p3,*p4;
{
	Vector3 v1,v2,N;
	double a;
	V3Sub(p2,p1,&v1); V3Sub(p4,p1,&v2); V3Cross(&v1,&v2,&N);
	a = V3Length(&N);
	V3Sub(p4,p3,&v1);V3Sub(p2,p3,&v2); V3Cross(&v1,&v2,&N);
	a += V3Length(&N);
	return(a/2.0);
	}
*/
static int get_dominant_axis(N)
Vector3 *N;
{
	double maxval;
	Point3 anorm;
		/* Computes Dominant axis.
		   0 - for X-axis, 1 - for Y-axis, 2 - for Z axis */
	anorm.x = fabs(N->x);
	anorm.y = fabs(N->y);
	anorm.z = fabs(N->z);
	maxval  = MAX(anorm.x, MAX(anorm.y,anorm.z));
	return((maxval == anorm.x)?0:((maxval == anorm.y)?1:2));
}
static double trapez_area(n,vlist,N,dominant_axis)
/* Polygon area computation by trapezium area summation method */
int n;
Point3 vlist[];
Vector3 *N;
int dominant_axis;
{
	int i;
	double area = 0.;
	Point2 cur,last;
	
	project_point(last,vlist[n-1],dominant_axis);
	for(i=0;i<n;i++,last=cur){
		project_point(cur,vlist[i],dominant_axis);
		area += (cur.x-last.x)*(cur.y+last.y);
	}
	area /= ((dominant_axis==0)?N->x:((dominant_axis==1)?N->y:N->z));
	return(0.5*fabs(area));
}

static void quadrilateral_preprocess(grid,un,vn,q,bb,b)
Surface_Grid_structure *grid;
int un,vn;
Quadrilateral *q;
Box3 *bb,*b;
{
	int i,j,n;
	Point3 Pa, Pb,Pc,Pd;
	double u,v,du=1.0/un,dv=1.0/vn;
	double area();
	/* Compute Bounding Box */
	bb->max.x=MAX(MAX(q->P00.x,q->P10.x),MAX(q->P11.x,q->P01.x));
        bb->max.y=MAX(MAX(q->P00.y,q->P10.y),MAX(q->P11.y,q->P01.y));
	bb->max.z=MAX(MAX(q->P00.z,q->P10.z),MAX(q->P11.z,q->P01.z));

	bb->min.x=MIN(MIN(q->P00.x,q->P10.x),MIN(q->P11.x,q->P01.x));
        bb->min.y=MIN(MIN(q->P00.y,q->P10.y),MIN(q->P11.y,q->P01.y));
	bb->min.z=MIN(MIN(q->P00.z,q->P10.z),MIN(q->P11.z,q->P01.z));
	update_scene_bound(bb,b);
	/* Compute Other parameters. */
       	Pa.x=q->P00.x-q->P10.x+q->P11.x-q->P01.x;
       	Pa.y=q->P00.y-q->P10.y+q->P11.y-q->P01.y;
       	Pa.z=q->P00.z-q->P10.z+q->P11.z-q->P01.z;
       	Pb.x=q->P10.x-q->P00.x;
       	Pb.y=q->P10.y-q->P00.y;
       	Pb.z=q->P10.z-q->P00.z;
       	Pc.x=q->P01.x-q->P00.x;
       	Pc.y=q->P01.y-q->P00.y;
       	Pc.z=q->P01.z-q->P00.z;
       	Pd = q->P00;
       	V3Cross(&Pb,&Pc,&(q->plane_normal));
	V3Normalize(&(q->plane_normal));
	/* Compute grid area. */
	for(n=0,v=0.0,i=0;i<vn;i++,v+=dv){
		Point3 f[4];
		f[0] = get_bilinear(q,(double)0.,v);
		f[3] = get_bilinear(q,(double)0.,v+dv);
		for(u=0.0,j=0;j<un;j++,u+=du,n++){
			f[1] = get_bilinear(q,u+du,v);
			f[2] = get_bilinear(q,u+du,v+dv);
			grid[n].area=trapez_area(4,f,&(q->plane_normal),
				get_dominant_axis(&(q->plane_normal)));
			f[0] = f[1]; f[3] = f[2];
		}
	}
		/* Compute the Coordinate Transformation Matrix */
	/*
	{
		Vector3 U,V,W;
		Matrix4 m;
		V=Pc;
		V3Normalize(&V);
		V3Cross(&V,&(q->plane_normal),&U);
		get_view_trans_matrix(&U,&V,&(q->plane_normal),&m);
		U.x=V.y=W.z=1.0;
		U.y=U.z=V.x=V.z=W.x=W.y=0.0;
		U = transform_vector(&U,&m);
		V = transform_vector(&V,&m);
		W = transform_vector(&W,&m);
		get_view_trans_matrix(&U,&V,&W,&(q->local_trans));
	}
	*/
	align_Z_to_axis(&(q->plane_normal),&(q->local_trans)); 

       	q->plane_constant = -V3Dot(&(q->plane_normal),&(q->P00));

       	V3Cross(&Pa,&(q->plane_normal),&(q->Na));
       	V3Cross(&Pb,&(q->plane_normal),&(q->Nb));
       	V3Cross(&Pc,&(q->plane_normal),&(q->Nc));

       	q->Du0 = V3Dot(&(q->Nc),&Pd);
       	q->Du1 = V3Dot(&(q->Na),&Pd)+V3Dot(&(q->Nc),&Pb);
       	q->Du2 = V3Dot(&(q->Na),&Pb);
       	if (fabs(q->Du2) > EPSILON){
              		double frac = 1.0/q->Du2;
              		q->Quy.x = -q->Nc.x*frac; 
		q->Quy.y = -q->Nc.y*frac; 
		q->Quy.z = -q->Nc.z*frac;
        		q->Duy=q->Du0*frac;
              		frac /= 2.0;
        		q->Qux.x=q->Na.x*frac; 
		q->Qux.y=q->Na.y*frac; 
		q->Qux.z=q->Na.z*frac;
        		q->Dux = -q->Du1*frac;
       	}
       	q->Dv0 = V3Dot(&(q->Nb),&Pd);
       	q->Dv1 = V3Dot(&(q->Na),&Pd)+V3Dot(&(q->Nb),&Pc);
       	q->Dv2 = V3Dot(&(q->Na),&Pc);
       	if (fabs(q->Dv2) > EPSILON){
              		double frac = 1.0/q->Dv2;
              		q->Qvy.x = -q->Nb.x*frac; 
		q->Qvy.y = -q->Nb.y*frac; 
		q->Qvy.z = -q->Nb.z*frac;
              		q->Dvy=q->Dv0*frac;
              		frac /= 2.0;
              		q->Qvx.x=q->Na.x*frac; 
		q->Qvx.y=q->Na.y*frac; 
		q->Qvx.z=q->Na.z*frac;
              		q->Dvx = -q->Dv1*frac;
       	}
}

static void sphere_preprocess(grid,un,vn,s,bb,b)
Surface_Grid_structure *grid;
int un,vn;
Sphere *s;
Box3 *bb,*b;
{
	int i,j,n;
	double area,mult,prevcos=1.0,nextcos;
	double v,du=1.0/un,dv=1.0/vn;
	/* Compute Bounding Box */
	bb->min.x = s->center.x - s->radius;
	bb->min.y = s->center.y - s->radius;
	bb->min.z = s->center.z - s->radius;

	bb->max.x = s->center.x + s->radius;
	bb->max.y = s->center.y + s->radius;
	bb->max.z = s->center.z + s->radius;
	update_scene_bound(bb,b);
	/* Compute Area of the Grid Elements. */
	mult = 2.*PI*s->radius*s->radius*du;
	for(n=0,v=0.0,i=0;i<vn;i++,v+=dv){
		nextcos = cos(PI*(v+dv));
		area=mult*fabs(prevcos-nextcos);
		prevcos=nextcos;
		for(j=0;j<un;j++,n++)
			grid[n].area=area;
	}
}

void translate(dx,dy,dz,m)
double dx,dy,dz;
Matrix4 *m;
{
	m->element[0][0] =
		m->element[1][1] =
			m->element[2][2] =
				m->element[3][3] = 1.;
	m->element[0][1] = m->element[0][2] = m->element[0][3] =
		m->element[1][0] = m->element[1][2] = m->element[1][3] =
			m->element[2][0] = m->element[2][1] = m->element[2][3] = 0.;
	m->element[3][0] = dx;
	m->element[3][1] = dy;
	m->element[3][2] = dz;
}

void scale(sx,sy,sz,m)
double sx,sy,sz;
Matrix4 *m;
{
	m->element[0][0] = sx;
	m->element[1][1] = sy;
	m->element[2][2] = sz;

	m->element[3][3] = 1.;
	m->element[0][1] = m->element[0][2] = m->element[0][3] =
		m->element[1][0] = m->element[1][2] = m->element[1][3] =
			m->element[2][0] = m->element[2][1] = m->element[2][3] =
				m->element[3][0] = m->element[3][1] = m->element[3][2] = 0.;
}

void bound_transform(b,m)
Box3 *b;
Matrix4 *m;
{
	Point3 p;
	Box3 bb;
	bb = *b; 
	b->min.x = b->min.y = b->min.z =  LARGE;
	b->max.x = b->max.y = b->max.z = -LARGE;
	p = bb.min; V3MulPointByMatrix(&p,m);
	b->min.x = MIN(b->min.x,p.x); b->min.y = MIN(b->min.y,p.y); b->min.z = MIN(b->min.z,p.z);
	b->max.x = MAX(b->max.x,p.x); b->max.y = MAX(b->max.y,p.y); b->max.z = MAX(b->max.z,p.z);
	p = bb.min; p.x = bb.max.x; V3MulPointByMatrix(&p,m);
	b->min.x = MIN(b->min.x,p.x); b->min.y = MIN(b->min.y,p.y); b->min.z = MIN(b->min.z,p.z);
	b->max.x = MAX(b->max.x,p.x); b->max.y = MAX(b->max.y,p.y); b->max.z = MAX(b->max.z,p.z);
	p = bb.min; p.y = bb.max.y; V3MulPointByMatrix(&p,m);
	b->min.x = MIN(b->min.x,p.x); b->min.y = MIN(b->min.y,p.y); b->min.z = MIN(b->min.z,p.z);
	b->max.x = MAX(b->max.x,p.x); b->max.y = MAX(b->max.y,p.y); b->max.z = MAX(b->max.z,p.z);
	p = bb.min; p.x = bb.max.x; p.y = bb.max.y; V3MulPointByMatrix(&p,m);
	b->min.x = MIN(b->min.x,p.x); b->min.y = MIN(b->min.y,p.y); b->min.z = MIN(b->min.z,p.z);
	b->max.x = MAX(b->max.x,p.x); b->max.y = MAX(b->max.y,p.y); b->max.z = MAX(b->max.z,p.z);

	p = bb.max; V3MulPointByMatrix(&p,m);
	b->min.x = MIN(b->min.x,p.x); b->min.y = MIN(b->min.y,p.y); b->min.z = MIN(b->min.z,p.z);
	b->max.x = MAX(b->max.x,p.x); b->max.y = MAX(b->max.y,p.y); b->max.z = MAX(b->max.z,p.z);
	p = bb.max; p.x = bb.min.x; V3MulPointByMatrix(&p,m);
	b->min.x = MIN(b->min.x,p.x); b->min.y = MIN(b->min.y,p.y); b->min.z = MIN(b->min.z,p.z);
	b->max.x = MAX(b->max.x,p.x); b->max.y = MAX(b->max.y,p.y); b->max.z = MAX(b->max.z,p.z);
	p = bb.max; p.y = bb.min.y; V3MulPointByMatrix(&p,m);
	b->min.x = MIN(b->min.x,p.x); b->min.y = MIN(b->min.y,p.y); b->min.z = MIN(b->min.z,p.z);
	b->max.x = MAX(b->max.x,p.x); b->max.y = MAX(b->max.y,p.y); b->max.z = MAX(b->max.z,p.z);
	p = bb.max; p.x = bb.min.x; p.y = bb.min.y; V3MulPointByMatrix(&p,m);
	b->min.x = MIN(b->min.x,p.x); b->min.y = MIN(b->min.y,p.y); b->min.z = MIN(b->min.z,p.z);
	b->max.x = MAX(b->max.x,p.x); b->max.y = MAX(b->max.y,p.y); b->max.z = MAX(b->max.z,p.z);
}

static void cylinder_preprocess(grid,un,vn,c,bb,b)
Surface_Grid_structure *grid;
int un,vn;
Cylinder *c;
Box3 *bb,*b;
{
	Vector3 dir;
	Matrix4 m,m1;
	double du=1.0/un,dv=1.0/vn,area;
	int i,j,n;

	/* Compute derived data */
	c->height = V3DistanceBetween2Points(&(c->base_point1),&(c->base_point2));
	/* Compute Transformations */
	V3Sub(&(c->base_point2),&(c->base_point1),&dir); V3Normalize(&dir);
	translate(-c->base_point1.x,-c->base_point1.y,-c->base_point1.z,&m);
	align_axis_to_Z(&dir,&m1);
	V3MatMul(&m,&m1,&(c->world_to_object));
	align_Z_to_axis(&dir,&m);
	translate(c->base_point1.x,c->base_point1.y,c->base_point1.z,&m1);
	V3MatMul(&m,&m1,&(c->object_to_world));

	/* Compute Bounding Box */
	bb->min.x = bb->min.y = -c->radius; bb->min.z = 0.;
	bb->max.x = bb->max.y = c->radius; bb->max.z = c->height;
	bound_transform(bb,&(c->object_to_world));
	update_scene_bound(bb,b);
	/* Compute Area of the Grid Elements */
	area = 2.*PI*du*c->radius*c->height*dv;
	for(n=0,i=0;i<vn;i++)
		for(j=0;j<un;j++,n++)
			grid[n].area=area;
}

static void cone_preprocess(grid,un,vn,c,bb,b)
Surface_Grid_structure *grid;
int un,vn;
Cone *c;
Box3 *bb,*b;
{
	Point3 origin;
	Matrix4 m,m1,m2,m3;
	double H1,H2,tan_theta;
	double v,du=1.0/un,dv=1.0/vn,mult;
	int i,j,n;

	/* Make sure 0 <= apex_radius <= base_radius. */
	if (c->apex_radius > c->base_radius){
		/* interchange */
		Point3 p;
		double r;
		p=c->base_point; c->base_point = c->apex_point; c->apex_point = p;
		r=c->base_radius;c->base_radius=c->apex_radius; c->apex_radius= r;
	}

	/* Compute derived data */
	H1 = V3DistanceBetween2Points(&(c->base_point),&(c->apex_point));
	tan_theta = (c->base_radius-c->apex_radius)/H1;
	H2 = c->apex_radius/tan_theta;
	c->distance = H2/(H1+H2);
	V3Sub(&(c->base_point),&(c->apex_point),&(c->axis)); V3Normalize(&(c->axis));

	/* Compute Transformations */
	origin.x = c->apex_point.x - H2*c->axis.x;
	origin.y = c->apex_point.y - H2*c->axis.y;
	origin.z = c->apex_point.z - H2*c->axis.z;
	translate(-origin.x,-origin.y,-origin.z,&m);
	align_axis_to_Z(&(c->axis),&m1);
	scale(1./c->base_radius,1./c->base_radius,1./(H1+H2),&m2);
	V3MatMul(&m,&m1,&m3);V3MatMul(&m3,&m2,&(c->world_to_object));
	
	scale(c->base_radius,c->base_radius,(H1+H2),&m);
	align_Z_to_axis(&(c->axis),&m1);
	translate(origin.x,origin.y,origin.z,&m2);
	V3MatMul(&m,&m1,&m3);V3MatMul(&m3,&m2,&(c->object_to_world));

	/* Compute Bounding Box */
	bb->min.x = bb->min.y = -1.; bb->min.z = c->distance;
	bb->max.x = bb->max.y = bb->max.z = 1.;
	bound_transform(bb,&(c->object_to_world));
	update_scene_bound(bb,b);

	/* Compute Area of the Grid Elements */
	mult = PI*du*sqrt((c->base_radius-c->apex_radius)*(c->base_radius-c->apex_radius)+H1*H1);
	for(n=0,v=0.0,i=0;i<vn;i++,v+=dv){
		double area=mult*((c->base_radius-c->apex_radius)*((v+dv)*(v+dv)-v*v)+2.0*c->apex_radius*dv);
		for(j=0;j<un;j++,n++)
			grid[n].area=area;
	}
}

static void ring_preprocess(grid,un,vn,r,bb,b)
Surface_Grid_structure *grid;
int un,vn;
Ring *r;
Box3 *bb,*b;
{
	Matrix4 m,m1;
	double extent;
	double v,du=1.0/un,dv=1.0/vn,mult;
	int i,j,n;

	if (r->inner_radius > r->outer_radius){
		double rad = r->outer_radius;
		r->outer_radius = r->inner_radius;
		r->inner_radius = rad;
	}
	V3Normalize(&(r->normal));

	/* Compute Transformations */
	translate(-r->center.x,-r->center.y,-r->center.z,&m);
	align_axis_to_Z(&(r->normal),&m1);
	V3MatMul(&m,&m1,&(r->world_to_object));

	align_Z_to_axis(&(r->normal),&m);
	translate(r->center.x,r->center.y,r->center.z,&m1);
	V3MatMul(&m,&m1,&(r->object_to_world));

	/* Compute Bounding Box */
	extent = r->outer_radius * sqrt(1. - r->normal.x*r->normal.x);
	bb->min.x = r->center.x - extent; 
	bb->max.x = r->center.x + extent;

	extent = r->outer_radius * sqrt(1. - r->normal.y*r->normal.y);
	bb->min.y = r->center.y - extent;
	bb->max.y = r->center.y + extent;

	extent = r->outer_radius * sqrt(1. - r->normal.z*r->normal.z);
	bb->min.z = r->center.z - extent;
	bb->max.z = r->center.z + extent;
	update_scene_bound(bb,b);
	/* Compute Area of the Grid Elements */
	mult = PI*du*(r->outer_radius*r->outer_radius - r->inner_radius*r->inner_radius);
	for(n=0,v=0.0,i=0;i<vn;i++,v+=dv){
		double area=mult*((v+dv)*(v+dv)-v*v);
		for(j=0;j<un;j++,n++)
			grid[n].area=area;
	}
}

static void polygon_preprocess(grid,un,vn,p,bb,b)
Surface_Grid_structure *grid;
int un,vn;
Polygon *p;
Box3 *bb,*b;
{
	int i,j;
	Point3 *point=p->vertex_list;
	/* Compute derived data */
		/* Normal by Newels Method */
	p->plane_normal.x = p->plane_normal.y = p->plane_normal.z = 0.;
	for(i=0,j=1;i<p->nvertices;i++,j++){
		if (j==p->nvertices)j=0;
		p->plane_normal.x += (point[i].y - point[j].y)*(point[i].z+point[j].z);
		p->plane_normal.y += (point[i].z - point[j].z)*(point[i].x+point[j].x);
		p->plane_normal.z += (point[i].x - point[j].x)*(point[i].y+point[j].y);
	}
	V3Normalize(&(p->plane_normal));
	p->plane_constant = -V3Dot(&(p->plane_normal),point);
	p->dominant_axis_flag = get_dominant_axis(&(p->plane_normal));
	align_Z_to_axis(&(p->plane_normal),&(p->local_trans));
	/* Compute Bounding Box */
	bb->min = bb->max = *point;
	for(i=1;i<p->nvertices;i++){
		bb->min.x = MIN(bb->min.x,point[i].x);
		bb->max.x = MAX(bb->max.x,point[i].x);
		bb->min.y = MIN(bb->min.y,point[i].y);
		bb->max.y = MAX(bb->max.y,point[i].y);
		bb->min.z = MIN(bb->min.z,point[i].z);
		bb->max.z = MAX(bb->max.z,point[i].z);
	}
	update_scene_bound(bb,b);
	/* Compute Area of the Grid Elements */
	if (un > 1 || vn > 1)
		error("Does not know any way to grid arbitrary polygon.");
	grid[0].area=
	trapez_area(p->nvertices,p->vertex_list,&(p->plane_normal),p->dominant_axis_flag);
}
#undef update_scene_bound
/************** End Object Specific Preprocess   **************/
/************** Begin Object Specific Input      **************/
/*
		for quadrilateral
			P00 
			P10 
			P11 
			P01 (Vertices Counter Clockwise) 
		for sphere and bubble
			centre
			radius
		for cylinder and tube
			base point1
			base point2
			radius
		for cone and cup
			apex point
			apex radius
			base point
			base radius
		for ring
			center
			normal
			inner radius
			outer radius
		for polygon
			number of vertices
			vertices
			.
			.
		for others
			nill
*/
static void *quadrilateral_input()
{
	Quadrilateral *q=(Quadrilateral *)malloc(sizeof(Quadrilateral));

	if (q==NULL)error("Cannot allocate Quadrilateral.");
	scanf("%lf%lf%lf",&(q->P00.x),&(q->P00.y),&(q->P00.z));
	scanf("%lf%lf%lf",&(q->P10.x),&(q->P10.y),&(q->P10.z));
	scanf("%lf%lf%lf",&(q->P11.x),&(q->P11.y),&(q->P11.z));
	scanf("%lf%lf%lf",&(q->P01.x),&(q->P01.y),&(q->P01.z));

	return((void *)q);
}

static void *polygon_input()
{
	int i;
	Point3 *point;
	Polygon *p=(Polygon *)malloc(sizeof(Polygon));

	if (p==NULL)error("Cannot allocate Polygon.");
	scanf("%d",&(p->nvertices));
	point = p->vertex_list = (Point3 *)malloc(sizeof(Point3) * p->nvertices);
	if (point == NULL) error("cannot allocate Polygon Vertices.");
	for (i=0; i < p->nvertices; i++,point++)
		scanf("%lf%lf%lf",&(point->x),&(point->y),&(point->z));
	return((void *)p);
}

static void *sphere_input()
{
	Sphere *s=(Sphere *)malloc(sizeof(Sphere));

	if (s==NULL)error("Cannot allocate Sphere.");
	scanf("%lf%lf%lf",&(s->center.x),&(s->center.y),&(s->center.z));
	scanf("%lf",&(s->radius));

	return((void *)s);
}

static void *cylinder_input()
{
	Cylinder *c=(Cylinder *)malloc(sizeof(Cylinder));

	if (c==NULL)error("Cannot allocate Cylinder/Tube.");
	scanf("%lf%lf%lf",&(c->base_point1.x),&(c->base_point1.y),&(c->base_point1.z));
	scanf("%lf%lf%lf",&(c->base_point2.x),&(c->base_point2.y),&(c->base_point2.z));
	scanf("%lf",&(c->radius));

	return((void *)c);
}


static void *cone_input()
{
	Cone *c=(Cone *)malloc(sizeof(Cone));

	if (c==NULL)error("Cannot allocate Cone/Cup.");
	scanf("%lf%lf%lf%lf",&(c->apex_point.x),&(c->apex_point.y),&(c->apex_point.z),
			&(c->apex_radius));
	scanf("%lf%lf%lf%lf",&(c->base_point.x),&(c->base_point.y),&(c->base_point.z),
			&(c->base_radius));

	return((void *)c);
}

static void *ring_input()
{
	Ring *r=(Ring *)malloc(sizeof(Ring));

	if (r==NULL)error("Cannot allocate Ring.");
	scanf("%lf%lf%lf",&(r->center.x),&(r->center.y),&(r->center.z));
	scanf("%lf%lf%lf",&(r->normal.x),&(r->normal.y),&(r->normal.z));
	scanf("%lf%lf",&(r->inner_radius),&(r->outer_radius));

	return((void *)r);
}

static void *no_input()
{
	error("This type Not implemented.");
}

/************** End Object Specific Input        **************/
/************** Begin Object Specific Normal     **************/

static Vector3 quadrilateral_normal(q,u,v)
Quadrilateral *q;
double u,v;
{
	return(q->plane_normal);
}

static Vector3 polygon_normal(p,u,v)
Polygon *p;
double u,v;
{
	return(p->plane_normal);
}

static Vector3 sphere_normal(s,u,v)
Sphere *s;
double u,v;
{
	Vector3 N;
	double phi = u * 2. * PI;
	double theta = v * PI;
	double sin_theta = sin(theta);

	N.x = sin_theta*cos(phi);
	N.y = sin_theta*sin(phi);
	N.z = cos(theta);

	return(N);
}

static Vector3 bubble_normal(s,u,v)
Sphere *s;
double u,v;
{
	Vector3 N;

	N=sphere_normal(s,u,v);
	V3Negate(&N);

	return(N);
}

static Vector3 cylinder_normal(c,u,v)
Cylinder *c;
double u,v;
{
	Vector3 N;
	double phi = u * 2. * PI;

	N.x = cos(phi);
	N.y = sin(phi);
	N.z = 0.;
	N = transform_vector(&N,&(c->object_to_world));

	return(N);
}

static Vector3 tube_normal(c,u,v)
Cylinder *c;
double u,v;
{
        Vector3 N;

	N=cylinder_normal(c,u,v);
	V3Negate(&N);

	return(N);
}

static Vector3 cone_normal(c,u,v)
Cone *c;
double u,v;
{
        Vector3 N,v1,v2;
	double phi = u * 2. * PI;

	v1.x = cos(phi); v1.y = sin(phi); v1.z = 1.;
	v1 = transform_vector(&v1,&(c->object_to_world));
	V3Cross(&v1,&(c->axis),&v2);
	V3Cross(&v1,&v2,&N);
	V3Normalize(&N);

        return(N);
}

static Vector3 cup_normal(c,u,v)
Cone *c;
double u,v;
{
        Vector3 N;

	N=cone_normal(c,u,v);
	V3Negate(&N);

        return(N);
}

static Vector3 ring_normal(r,u,v)
Ring *r;
double u,v;
{
        return(r->normal);
}

static Vector3 no_normal(s,u,v)
void *s;
double u,v;
{
	error("No Normal Yet.");
}

/************** End Object Specific Normal       **************/
/************** Begin Object Specific Local Matrix ************/

static Matrix4 quadrilateral_matrix(N,q)
Vector3 *N;
Quadrilateral *q;
{
	return(q->local_trans);
}


static Matrix4 polygon_matrix(N,p)
Vector3 *N;
Polygon *p;
{
	return(p->local_trans);
}

static Matrix4 align_matrix(N,s)
Vector3 *N;
void *s;
{
	Matrix4 m;

	align_Z_to_axis(N,&m);

	return(m);
}

static Matrix4 ring_matrix(N,r)
/*
	CAUTION: This matrix also has a nonzero translation vector.
*/ 
Vector3 *N;
Ring *r;
{
	return(r->object_to_world);
}

static Matrix4 no_matrix(N,s)
Vector3 *N;
void *s;
{
	error("No Matrix Yet.");
}

/************** End Object Specific Local Matrix   ************/
/************** Begin Object Specific Area Sampling ***********/

static Point3 sphere_point(s,u,v)
Sphere *s;
double u,v;
{
	Point3 p;
	double R =s->radius;
	Point3 C;
	C =s->center;
	p.x = R*sin(PI*v)*cos(2.*PI*u);
	p.y = R*sin(PI*v)*sin(2.*PI*u);
	p.z = R*cos(PI*v);
	V3Add(&p,&C,&p);
	return(p);
}

static Point3 cylinder_point(c,u,v)
Cylinder *c;
double u,v;
{
	Point3 p;
	double R = c->radius;
	double H = c->height;
	double phi = 2.*PI*u;
	p.x = R*cos(phi);
	p.y = R*sin(phi);
	p.z = v*H;
	V3MulPointByMatrix(&p,&(c->object_to_world));
	return(p);
}
static Point3 cone_point(c,u,v)
Cone *c;
double u,v;
{
	Point3 posn;
	double phi=u * 2. * PI;
	Point3 p;
	p.z = (1.-c->distance)*v+c->distance;
	p.x = p.z * cos(phi); p.y = p.z * sin(phi);
	V3MulPointByMatrix(&p,&(c->object_to_world));
	return(p);
}
static Point3 ring_point(r,u,v)
Ring *r;
double u,v;
{
	Point3 p;
	double radius=r->inner_radius+
		v*(r->outer_radius-r->inner_radius);
	double phi=2. * u * PI;
	p.x = radius * cos(phi); p.y = radius * sin(phi); p.z = 0.;
	V3MulPointByMatrix(&p,&(r->object_to_world));
	return(p);
}
static Point3 triangle_point(p,u,v)
Point3 *p;
double u,v;
{
	Point3 posn;
	posn.x = p[0].x + (1.-v)*u*(p[1].x-p[0].x)+ v*(p[2].x-p[0].x);
	posn.y = p[0].y + (1.-v)*u*(p[1].y-p[0].y)+ v*(p[2].y-p[0].y);
	posn.z = p[0].z + (1.-v)*u*(p[1].z-p[0].z)+ v*(p[2].z-p[0].z);
	return(posn);
}
static Point3 polygon_point(p,u,v)
Polygon *p;
double u,v;
{
	return(p->vertex_list[0]); /* Arbitrary */
}
static Point3 no_object_point(p,u,v)
void *p;
double u,v;
{
	error("Cannot compute point for this object.");
}
/************** End Object Specific Area Sampling   ***********/

/************** Begin Grid Centre Point             ***********/
#define get_U ((uindex+0.5) * 1./object[n].grid_h_reso)
#define get_V ((vindex+0.5) * 1./object[n].grid_v_reso)
static Point3 quadrilateral_grid_centre_point(n,vindex,uindex)
int n;
int vindex,uindex;
{
	double u = get_U;
	double v = get_V;
	return(get_bilinear(object[n].object_specific_structure,u,v));
}
static Point3 sphere_grid_centre_point(n,vindex,uindex)
int n;
int vindex,uindex;
{
	double u = get_U;
	double v = get_V;
	return(sphere_point(object[n].object_specific_structure,u,v));
}
static Point3 cylinder_grid_centre_point(n,vindex,uindex)
int n;
int vindex,uindex;
{
	double u = get_U;
	double v = get_V;
	return(cylinder_point(object[n].object_specific_structure,u,v));
}
static Point3 cone_grid_centre_point(n,vindex,uindex)
int n;
int vindex,uindex;
{
	double u = get_U;
	double v = get_V;
	return(cone_point(object[n].object_specific_structure,u,v));
}
static Point3 ring_grid_centre_point(n,vindex,uindex)
int n;
int vindex,uindex;
{
	double u = get_U;
	double v = get_V;
	return(ring_point(object[n].object_specific_structure,u,v));
}
static Point3 polygon_grid_centre_point(n,vindex,uindex)
int n;
int vindex,uindex;
{
	Polygon *p=(Polygon *)object[n].object_specific_structure;
	if (p->nvertices==3){/* Triangle */
		double u = get_U;
		double v = get_V;
		return(triangle_point(p->vertex_list,u,v));
	}
	else /* We donot have a parameterisation method for
		arbitrary polygon. So I can send any point 
		on the surface. */
		return(p->vertex_list[0]);
}
static Point3 no_grid_centre_point(n,vindex,uindex)
int n;
int vindex,uindex;
{
	error("Grid center Point Not Implemented.");
}
#undef get_U
#undef get_V
/************** End Grid Centre Point               ***********/

ObjectFunctions ofunc[MAX_GEOMETRY_TYPE]={
/*
	NOTE : Each of the Normal routines must return Normalized Vectors. 
*/
{
	quadrilateral_input,
	quadrilateral_normal,
	intersect_quadrilateral,
	quadrilateral_preprocess,
	quadrilateral_matrix,
	quadrilateral_grid_centre_point,
	get_bilinear
},
{
	sphere_input,
	sphere_normal,
	intersect_sphere,
	sphere_preprocess,
	align_matrix,
	sphere_grid_centre_point,
	sphere_point
},
{
	sphere_input,
	bubble_normal,
	intersect_bubble,
	sphere_preprocess,
	align_matrix,
	sphere_grid_centre_point,
	sphere_point
},
{
	cylinder_input,
	cylinder_normal,
	intersect_cylinder,
	cylinder_preprocess,
	align_matrix,
	cylinder_grid_centre_point,
	cylinder_point
},
{
	cylinder_input,
	tube_normal,
	intersect_tube,
	cylinder_preprocess,
	align_matrix,
	cylinder_grid_centre_point,
	cylinder_point,
},
{
	cone_input,
	cone_normal,
	intersect_cone,
	cone_preprocess,
	align_matrix,
	cone_grid_centre_point,
	cone_point
},
{
	cone_input,
	cup_normal,
	intersect_cup,
	cone_preprocess,
	align_matrix,
	cone_grid_centre_point,
	cone_point
},
{
	ring_input,
	ring_normal,
	intersect_ring,
	ring_preprocess,
	ring_matrix,
	ring_grid_centre_point,
	ring_point
},
{
	polygon_input,
	polygon_normal,
	intersect_polygon,
	polygon_preprocess,
	polygon_matrix,
	polygon_grid_centre_point,
	polygon_point
},
{
	no_input,
	no_normal,
	intersect_other_geometry,
	null_preprocess,
	no_matrix,
	no_grid_centre_point,
	no_object_point
}
};
#undef get_V
#undef get_U
#undef update_scene_bound