/*
	name:	intersct.c

	Intersection-routines
	---------------------

	These routines will calculate the t-parameter of a line,
	which corresponds to the point where the given line intersects with
	the given primitive. The routines return zero (0) if no intersection
	occured.

*/

#include  <stdio.h>

#include  "defs.h"
#include  "extern.h"


/**********************************************************************
 *
 *	Get (possible) intersection with given object
 *
 **********************************************************************/


double Intersect_Object(LINE *Line, OBJECT *Object)
{
	LINE	TransformedLine;
	double	t;
	long	i;
	VECTOR	tempv;

	CopyPoint(&TransformedLine.Origin,&Line->Origin);
	CopyVector(&TransformedLine.Direction,&Line->Direction);

	if(Object->Transform.NumTransforms>0) {
	    for(i=Object->Transform.NumTransforms-1L;i>=0L;i--) {
		switch(Object->Transform.Entry[i].Type) {
		    case TRANSFORM_SCALE:
			TransformedLine.Origin.x=TransformedLine.Origin.x/Object->Transform.Entry[i].Values.x;
			TransformedLine.Origin.y=TransformedLine.Origin.y/Object->Transform.Entry[i].Values.y;
			TransformedLine.Origin.z=TransformedLine.Origin.z/Object->Transform.Entry[i].Values.z;
			TransformedLine.Direction.x=TransformedLine.Direction.x/Object->Transform.Entry[i].Values.x;
			TransformedLine.Direction.y=TransformedLine.Direction.y/Object->Transform.Entry[i].Values.y;
			TransformedLine.Direction.z=TransformedLine.Direction.z/Object->Transform.Entry[i].Values.z;
			break;
		    case TRANSFORM_MOVE:
			TransformedLine.Origin.x=TransformedLine.Origin.x-Object->Transform.Entry[i].Values.x;
			TransformedLine.Origin.y=TransformedLine.Origin.y-Object->Transform.Entry[i].Values.y;
			TransformedLine.Origin.z=TransformedLine.Origin.z-Object->Transform.Entry[i].Values.z;
			break;
		    case TRANSFORM_ROTATE:
			NegVector(&tempv,&Object->Transform.Entry[i].Values);
			RotatePoint(&TransformedLine.Origin,&tempv);
			RotateVector(&TransformedLine.Direction,&tempv);
			break;
		    case TRANSFORM_NONE:
		    default:
			break;
		}
	    }
	}

	switch(Object->ShapeType) {
	    case SHAPE_PLANE:
		t=Intersect_Plane(&TransformedLine, (PLANE *)(Object->Shape));
		break;
	    case SHAPE_SPHERE:
		t=Intersect_Sphere(&TransformedLine, (SPHERE *)(Object->Shape));
		break;
	    case SHAPE_ELLIPSOID:
		t=Intersect_Ellipsoid(&TransformedLine, (ELLIPSOID *)(Object->Shape));
		break;
	    case SHAPE_TRIANGLE:
		t=Intersect_Triangle(&TransformedLine, (TRIANGLE *)(Object->Shape));
		break;
	    case SHAPE_BOX:
		t=Intersect_Box(&TransformedLine, (BOX *)(Object->Shape));
		break;
	    case SHAPE_DISC:
		t=Intersect_Disc(&TransformedLine, (DISC *)(Object->Shape));
		break;
	    case SHAPE_CYLINDER:
		t=Intersect_Cylinder(&TransformedLine, (CYLINDER *)(Object->Shape));
		break;
	    default:
		t=0.0;
	}

	return(t);
}



/**********************************************************************
 *
 *	These are the actual intersection calculation routines
 *
 **********************************************************************/


double Intersect_Plane(LINE *Line, PLANE *Plane)
{
	register double	t,pnx,pny,pnz,a,tmp;
	double	lvx,lvy,lvz,lx0,ly0,lz0;

	PlaneIntersectionAttempts++;

	lvx=Line->Direction.x; lvy=Line->Direction.y; lvz=Line->Direction.z;
	lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;
	pnx=Plane->Normal.x; pny=Plane->Normal.y; pnz=Plane->Normal.z; a=Plane->a;

	tmp=pnx*lvx+pny*lvy+pnz*lvz;
	if(tmp!=0.0) {
		t=-(pnx*lx0+pny*ly0+pnz*lz0-a)/tmp;
		if(t<EPSILON) { t=0.0; }
		else { PlaneIntersections++; }
	}
	else { t=0.0; }

	return(t);
}



double Intersect_Sphere(LINE *Line, SPHERE *Sphere)
{
	double	vx,vy,vz,lx0,ly0,lz0,sx0,sy0,sz0,r;

	register double	s,sroot,vxdx,vydy,vzdz,vxsqr,vysqr,vzsqr;
	double	t;
	double	dx,dy,dz;
	double	dxsqr,dysqr,dzsqr;
	double	vsum,t1,t2;

	SphereIntersectionAttempts++;

	vx=Line->Direction.x; vy=Line->Direction.y; vz=Line->Direction.z;
	lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;
	sx0=Sphere->Centre.x; sy0=Sphere->Centre.y; sz0=Sphere->Centre.z;
	r=Sphere->r;

						/* These are used more than once, */
	dx=lx0-sx0; dy=ly0-sy0; dz=lz0-sz0;	/* thus much time is saved by not */
	vxdx=vx*dx; vydy=vy*dy; vzdz=vz*dz;	/* performing the same operation  */
	dxsqr=dx*dx; dysqr=dy*dy; dzsqr=dz*dz;	/* twice (or even more).          */
	vxsqr=vx*vx; vysqr=vy*vy; vzsqr=vz*vz;

	vsum=(vxsqr+vysqr+vzsqr);

	s=-(vxdx+vydy+vzdz);

	sroot=2*(vxdx*(vydy+vzdz)+vydy*vzdz)+r*r*vsum;
	sroot-=vxsqr*(dysqr+dzsqr)+vysqr*(dxsqr+dzsqr)+vzsqr*(dxsqr+dysqr);

	if(sroot>=0) {
		sroot=sqrt(sroot);

		t1=(s+sroot)/vsum;
		t2=(s-sroot)/vsum;

		if((t1>EPSILON)&&((t1<t2)||(t2<EPSILON))) {
			t=t1;
			SphereIntersections++;
		}
		else {
			if((t2>EPSILON)&&((t2<t1)||(t1<EPSILON))) {
				t=t2;
				SphereIntersections++;
			}
			else { t=(double) 0; }
		}
	}
	else t=(double) 0;

	return(t);
}



double Intersect_Ellipsoid(LINE *Line, ELLIPSOID *Ellipsoid)
{
	double	vx,vy,vz,lx0,ly0,lz0,ex0,ey0,ez0,a,b,c;

	register double	sroot,asqr,bsqr,csqr,vxdx,vydy,vzdz;
	double	t;
	double	dx,dy,dz;
	double	vxsqr,vysqr,vzsqr,dxsqr,dysqr,dzsqr;
	double	r,s,t1,t2;

	EllipsoidIntersectionAttempts++;

	vx=Line->Direction.x; vy=Line->Direction.y; vz=Line->Direction.z;
	lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;
	ex0=Ellipsoid->Centre.x; ey0=Ellipsoid->Centre.y; ez0=Ellipsoid->Centre.z;
	a=Ellipsoid->Radius.x; b=Ellipsoid->Radius.y; c=Ellipsoid->Radius.z;

	asqr=a*a; bsqr=b*b; csqr=c*c;		/* These are used more than once, */
	dx=lx0-ex0; dy=ly0-ey0; dz=lz0-ez0;	/* thus much time is saved by not */
	vxdx=vx*dx; vydy=vy*dy; vzdz=vz*dz;	/* performing the same operation  */
	dxsqr=dx*dx; dysqr=dy*dy; dzsqr=dz*dz;	/* twice (or even more).          */
	vxsqr=vx*vx; vysqr=vy*vy; vzsqr=vz*vz;

	s=-((vxdx/asqr)+(vydy/bsqr)+(vzdz/csqr));

	sroot=csqr*(vxsqr*bsqr+vysqr*asqr)+vzsqr*asqr*bsqr;
	sroot+=2*(vxdx*(csqr*vydy+bsqr*vzdz)+asqr*vydy*vzdz);
	sroot-=vxsqr*(csqr*dysqr+bsqr*dzsqr);
	sroot-=vysqr*(csqr*dxsqr+asqr*dzsqr);
	sroot-=vzsqr*(bsqr*dxsqr+asqr*dysqr);

	if(sroot>=0) {
		sroot=sqrt(sroot)/(a*b*c);

		r=((vxsqr/asqr)+(vysqr/bsqr)+(vzsqr/csqr));

		t1=(s+sroot)/r;
		t2=(s-sroot)/r;

		if((t1>EPSILON)&&((t1<t2)||(t2<EPSILON))) {
			t=t1;
			EllipsoidIntersections++;
		}
		else {
			if((t2>EPSILON)&&((t2<t1)||(t1<EPSILON))) {
				t=t2;
				EllipsoidIntersections++;
			}
			else { t=(double) 0; }
		}
	}
	else t=0.0;

	return(t);
}


double Intersect_Triangle(LINE *Line, TRIANGLE *Triangle)
{
	double	t=0.0,ttmp,pnx,pny,pnz,a;
	double	lvx,lvy,lvz,lx0,ly0,lz0,angle1,angle2,angle3;
	VECTOR	v1,v2,v3;
	POINT	ip;

	TriangleIntersectionAttempts++;

	lvx=Line->Direction.x; lvy=Line->Direction.y; lvz=Line->Direction.z;
	lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;
	pnx=Triangle->Plane.Normal.x; pny=Triangle->Plane.Normal.y; pnz=Triangle->Plane.Normal.z; a=Triangle->Plane.a;

	ttmp=-(pnx*lx0+pny*ly0+pnz*lz0-a)/(pnx*lvx+pny*lvy+pnz*lvz);
	if(ttmp>EPSILON) {
	    ip.x=Line->Origin.x+ttmp*Line->Direction.x;	/* ip = intersection point */
	    ip.y=Line->Origin.y+ttmp*Line->Direction.y;
	    ip.z=Line->Origin.z+ttmp*Line->Direction.z;

	    if((ip.x>Triangle->Min.x)&&(ip.x<Triangle->Max.x)&&(ip.y>Triangle->Min.y)&&(ip.y<Triangle->Max.y)&&(ip.z>Triangle->Min.z)&&(ip.z<Triangle->Max.z)) {
		CreateVector(&v1,Triangle->Corners[1].x-Triangle->Corners[0].x,Triangle->Corners[1].y-Triangle->Corners[0].y,Triangle->Corners[1].z-Triangle->Corners[0].z);
		CreateVector(&v2,Triangle->Corners[2].x-Triangle->Corners[0].x,Triangle->Corners[2].y-Triangle->Corners[0].y,Triangle->Corners[2].z-Triangle->Corners[0].z);
		angle1=VectorsAngle(&v1,&v2);
		CreateVector(&v3,ip.x-Triangle->Corners[0].x,ip.y-Triangle->Corners[0].y,ip.z-Triangle->Corners[0].z);
		angle2=VectorsAngle(&v1,&v3);
		angle3=VectorsAngle(&v2,&v3);
		if((angle2+angle3)<=(angle1+EPSILON)) {
			CreateVector(&v1,Triangle->Corners[1].x-Triangle->Corners[2].x,Triangle->Corners[1].y-Triangle->Corners[2].y,Triangle->Corners[1].z-Triangle->Corners[2].z);
			CreateVector(&v2,Triangle->Corners[0].x-Triangle->Corners[2].x,Triangle->Corners[0].y-Triangle->Corners[2].y,Triangle->Corners[0].z-Triangle->Corners[2].z);
			angle1=VectorsAngle(&v1,&v2);
			CreateVector(&v3,ip.x-Triangle->Corners[2].x,ip.y-Triangle->Corners[2].y,ip.z-Triangle->Corners[2].z);
			angle2=VectorsAngle(&v1,&v3);
			angle3=VectorsAngle(&v2,&v3);
			if((angle2+angle3)<=(angle1+EPSILON)) {
				t=ttmp;
				TriangleIntersections++;
			}
		}
	    }
	}

	return(t);
}


double Intersect_Box(LINE *Line, BOX *Box)
{
	double	t=0.0,told=-1.0,ttmp,pnx,pny,pnz,a;
	double	lvx,lvy,lvz,lx0,ly0,lz0;
	int	i;
	POINT	ip;

	BoxIntersectionAttempts++;

	lvx=Line->Direction.x; lvy=Line->Direction.y; lvz=Line->Direction.z;
	lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;
	for(i=0;i<6;i++) {
		pnx=Box->Planes[i].Normal.x; pny=Box->Planes[i].Normal.y; pnz=Box->Planes[i].Normal.z; a=Box->Planes[i].a;

		ttmp=-(pnx*lx0+pny*ly0+pnz*lz0-a)/(pnx*lvx+pny*lvy+pnz*lvz);
		if((ttmp>EPSILON)&&((told<0.0)||(ttmp<told))) {
			switch(i) {
			    case 0:
			    case 1:
				ip.y=Line->Origin.y+ttmp*Line->Direction.y;
				ip.z=Line->Origin.z+ttmp*Line->Direction.z;
				if((ip.y>=Box->Corners[0].y)&&(ip.y<=Box->Corners[1].y)&&(ip.z>=Box->Corners[0].z)&&(ip.z<=Box->Corners[1].z)) {
					told=ttmp;
				}
				break;
			    case 2:
			    case 3:
				ip.x=Line->Origin.x+ttmp*Line->Direction.x;
				ip.z=Line->Origin.z+ttmp*Line->Direction.z;
				if((ip.x>=Box->Corners[0].x)&&(ip.x<=Box->Corners[1].x)&&(ip.z>=Box->Corners[0].z)&&(ip.z<=Box->Corners[1].z)) {
					told=ttmp;
				}
				break;
			    case 4:
			    case 5:
				ip.x=Line->Origin.x+ttmp*Line->Direction.x;
				ip.y=Line->Origin.y+ttmp*Line->Direction.y;
				if((ip.x>=Box->Corners[0].x)&&(ip.x<=Box->Corners[1].x)&&(ip.y>=Box->Corners[0].y)&&(ip.y<=Box->Corners[1].y)) {
					told=ttmp;
				}
				break;
			}
		}
	}
	if(told>EPSILON) {
		t=told;
		BoxIntersections++;
	}

	return(t);
}


double Intersect_Disc(LINE *Line, DISC *Disc)
{
	register double	t,pnx,pny,pnz,a,tmp;
	double	lvx,lvy,lvz,lx0,ly0,lz0,dx,dy,dz;

	DiscIntersectionAttempts++;

	lvx=Line->Direction.x; lvy=Line->Direction.y; lvz=Line->Direction.z;
	lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;
	pnx=Disc->Plane.Normal.x; pny=Disc->Plane.Normal.y; pnz=Disc->Plane.Normal.z; a=Disc->Plane.a;

	tmp=pnx*lvx+pny*lvy+pnz*lvz;
	if(tmp!=0.0) {
		t=-(pnx*lx0+pny*ly0+pnz*lz0-a)/tmp;
		if(t<EPSILON) { t=0.0; }
		else {
			dx=(Line->Origin.x+Line->Direction.x*t)-Disc->Centre.x;
			dy=(Line->Origin.y+Line->Direction.y*t)-Disc->Centre.y;
			dz=(Line->Origin.z+Line->Direction.z*t)-Disc->Centre.z;
			if((dx*dx+dy*dy+dz*dz)>(Disc->r*Disc->r)) { t=0.0; }
			else { DiscIntersections++; }
		}
	}
	else { t=0.0; }

	return(t);
}


double Intersect_Cylinder(LINE *Line, CYLINDER *Cylinder)
{
	register double	t=0.0,t1,t2,t3,t4,rsqr,lx0lvx,ly0lvy;
	double	lvx,lvy,lvz,lx0,ly0,lz0,tmp,tmp2,tmp3;
	POINT	ip1,ip2;

	CylinderIntersectionAttempts++;

	lvx=Line->Direction.x; lvy=Line->Direction.y; lvz=Line->Direction.z;
	lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;
	rsqr=Cylinder->r*Cylinder->r;
	
	lx0lvx=lx0*lvx; ly0lvy=ly0*lvy;
	tmp=2.0*lx0lvx*ly0lvy-lvx*lvx*ly0*ly0-lvy*lvy*lx0*lx0+lvx*lvx*rsqr+lvy*lvy*rsqr;

	if(tmp>=0.0) {
		tmp=sqrt(tmp);
		tmp2=-(lx0lvx+ly0lvy);
		tmp3=1.0/(lvx*lvx+lvy*lvy);
		t1=(tmp2+tmp)*tmp3;
		t2=(tmp2-tmp)*tmp3;
		ip1.x=lx0+lvx*t1; ip1.y=ly0+lvy*t1; ip1.z=lz0+lvz*t1;
		ip2.x=lx0+lvx*t2; ip2.y=ly0+lvy*t2; ip2.z=lz0+lvz*t2;

		tmp=ip1.x*Cylinder->Discs[0].Plane.Normal.x+ip1.y*Cylinder->Discs[0].Plane.Normal.y+ip1.z*Cylinder->Discs[0].Plane.Normal.z;
		if(tmp>Cylinder->Discs[0].Plane.a) { t1=0.0; }
		else {
			tmp=ip1.x*Cylinder->Discs[1].Plane.Normal.x+ip1.y*Cylinder->Discs[1].Plane.Normal.y+ip1.z*Cylinder->Discs[1].Plane.Normal.z;
			if(tmp>Cylinder->Discs[1].Plane.a) { t1=0.0; }
		}
		tmp=ip2.x*Cylinder->Discs[0].Plane.Normal.x+ip2.y*Cylinder->Discs[0].Plane.Normal.y+ip2.z*Cylinder->Discs[0].Plane.Normal.z;
		if(tmp>Cylinder->Discs[0].Plane.a) { t2=0.0; }
		else {
			tmp=ip2.x*Cylinder->Discs[1].Plane.Normal.x+ip2.y*Cylinder->Discs[1].Plane.Normal.y+ip2.z*Cylinder->Discs[1].Plane.Normal.z;
			if(tmp>Cylinder->Discs[1].Plane.a) { t2=0.0; }
		}

		t=0.0;

		if((t1<EPSILON)||(t2<EPSILON)) {
			t3=Intersect_Disc(Line, &Cylinder->Discs[0]);
			t4=Intersect_Disc(Line, &Cylinder->Discs[1]);
			if(t1>EPSILON) { t=t1; }
			if((t2>EPSILON)&&((t2<t)||(t<EPSILON))) { t=t2; }
			if((t3>EPSILON)&&((t3<t)||(t<EPSILON))) { t=t3; }
			if((t4>EPSILON)&&((t4<t)||(t<EPSILON))) { t=t4; }
		}
		else {
			if(t1>EPSILON) { t=t1; }
			if((t2>EPSILON)&&((t2<t)||(t<EPSILON))) { t=t2; }
		}
		if(t>EPSILON) { CylinderIntersections++; }
	}

	return(t);
}