/* object.c */

#include "object.h"

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

#include "config.h"
#include "geo.h"


/* ###########################################################################
   PRIVATE FORWARD DECLARATIONS
   ###########################################################################*/


num object_intersect(object_entry* obj, ray r);
num object_intersect_csg(object_meta_csg* csg, ray r, int *obj_id);


/* ###########################################################################
   GENERAL
   ###########################################################################*/

   
void object_transform(num* matrix, tuple t, num* out) {
	// Assume <matrix> is a 16 element transformation matrix
	num x[4] = {t.a, t.b, t.c, 1};
	
	int i, j;
	for (i = 0; i < 4; i++) for (j = 0; j < 4; j++) {
		out[i] += matrix[(4*i)+j]*x[j];
	}
}


/* ###########################################################################
   NORMALS
   ###########################################################################*/

   
tuple object_normal_polygon(tuple* v) {
	return geo_cross(geo_sub(v[0], v[1]), geo_sub(v[0], v[2]));
}

tuple object_normal(
	object_entry* obj, tuple p
#ifdef ENABLE_CSG_NORMALS
	// Point+Object isn't sufficient to reconstruct the ray for CSG tests.
	// Needs some thought into rearranging this
	, ray r
#endif	
) {
	switch( obj->object_type ) {
	
		case OBJ_TYPE_SPHERE: {
			object_sphere* sphere = obj->data;
			return geo_unit(geo_sub(p, sphere->centre));
		} break;
		
		case OBJ_TYPE_PLANE: {
			object_plane* plane = obj->data;
			return plane->normal;
		} break;
		
		case OBJ_TYPE_TRIANGLE: {
			object_triangle* tri = obj->data;
			return object_normal_polygon(tri->v);
		} break;
		
		case OBJ_TYPE_QUAD: {
			object_quad* quad = obj->data;
			return object_normal_polygon(quad->v);
		} break;
		
		case OBJ_TYPE_BOX: {
			object_box* box = obj->data;

#define BOXTEST(l, op) (fabs( p. l - box->centre. l op box->dimensions. l ) < EPS)			
			if BOXTEST(a,-) return geo_new( 1,0,0);
			if BOXTEST(a,+) return geo_new(-1,0,0);
			if BOXTEST(b,-) return geo_new(0, 1,0);
			if BOXTEST(b,+) return geo_new(0,-1,0);
			if BOXTEST(c,-) return geo_new(0,0, 1);
			if BOXTEST(c,+) return geo_new(0,0,-1);
#undef BOXTEST
			
#ifdef GEO_DEBUG
			fputs("WARNING: Couldn't determine surface normal of box\n", stderr);
#endif
			return geo_new(1,0,0);
		} break;
		
		/* Meta objects */
		
		/*
		case OBJ_TYPE_M_TRANSFORM: {
			
		} break;
		*/
		
		case OBJ_TYPE_CSG: {
#ifndef ENABLE_CSG_NORMALS		
			return geo_new(0,0,0);
#else
			// Super inefficient, recalculates intersection..
			object_meta_csg *csg = obj->data;
			int obj_id = 0;
			/*num t = */object_intersect_csg(csg, r, &obj_id);
			/*
			tuple norm =
				(obj_id & OBJ_HIT_SECOND) ?
				object_normal(csg->obj2, p, r) :
				object_normal(csg->obj1, p, r)
			;
			if (obj_id & ~OBJ_HIT_NEAR) {
				norm = geo_neg(norm);
			}
			return norm;
			*/
			switch (obj_id) {
				case (OBJ_HIT_FIRST | OBJ_HIT_NEAR): {
					return object_normal(csg->obj1, p, r);
				} break;
				case (OBJ_HIT_SECOND | OBJ_HIT_NEAR): {
					return object_normal(csg->obj2, p, r);
				} break;
				case (OBJ_HIT_FIRST | OBJ_HIT_FAR): {
					return geo_neg(object_normal(csg->obj1, p, r));
				} break;
				case (OBJ_HIT_SECOND | OBJ_HIT_FAR): {
					return geo_neg(object_normal(csg->obj2, p, r));
				} break;
				default: {
#ifdef GEO_DEBUG
					fputs("WARNING: Invalid result of CSG intersection\n", stderr);
#endif
					return geo_new(0,0,0);
				} break;
			}
			
#endif
		} break;
		
		/* Fallback */
		
		default: {
#ifdef GEO_DEBUG
			fputs("WARNING: Invalid object type\n", stderr);
#endif
			return geo_new(0, 0, 0);
		} break;
	}
}


/* ###########################################################################
   INTERSECTIONS
   ###########################################################################*/


int object_intersect_box(ray* r, object_box* box, num* out) {

	// Axis-aligned ray-box intersection using the Slab technique
	// Ref: http://www.siggraph.org/education/materials/HyperGraph/raytrace/rtinter3.htm
	// Ref: http://ompf.org/ray/ray_box.html
	// Unrolled version (or "funrolled" as GCC would say :P) Let's annoy the
	// branch predictor a little less where possible
	
	num posc, posf, lmin, lmax;
	
	tuple max = geo_add(box->centre, box->dimensions);
	tuple min = geo_sub(box->centre, box->dimensions);
	
	if (r->start.a >= 0) swapf(&max.a, &min.a);
	if (r->start.b >= 0) swapf(&max.b, &min.b);
	if (r->start.c >= 0) swapf(&max.c, &min.c);
	
	posc = (min.a - r->start.a) / r->dir.a;
	posf = (max.a - r->start.a) / r->dir.a;
	lmin = minf(posc,posf);
	lmax = maxf(posc,posf);

	posc = (min.b - r->start.b) / r->dir.b;
	posf = (max.b - r->start.b) / r->dir.b;
	lmin = maxf(minf(posc,posf), lmin);
	lmax = minf(maxf(posc,posf), lmax);
	
	posc = (min.c - r->start.c) / r->dir.c;
	posf = (max.c - r->start.c) / r->dir.c;
	lmin = maxf(minf(posc,posf), lmin);
	lmax = minf(maxf(posc,posf), lmax);

	*out = lmin;
	
	return ((lmax >= 0.0) && (lmax >= lmin));
}

// Returns 1 if (test) is within the polygon (vertices), or 0 otherwise
int object_intersect_poly(tuple* v, int vertex_count, tuple point) {
		
#ifdef GEO_DEBUG
	if (vertex_count > 4) {
		// We only call this with triangles and quads (don't want to malloc in loop!)
		fputs("WARNING: object_intersect_poly called with vertex count > 4\n", stderr);
	}
#endif
	
	// Flatten vertices to 2D by discarding the coordinate with the smallest
	// range. This works well, but it kills the performance to recalculate it
	// over and over for every ray and reflected ray..
	
	num tx, ty, x[4], y[4];
	int i, j, c;
	tuple range = v[0];
	
	for (i = 1; i < vertex_count; i++) range = geo_add(range, v[i]);
	c = (range.a < range.b) ?
		((range.a < range.c) ? 0 : 2) :
		((range.b < range.c) ? 1 : 2)
	;
				
	switch(c) {
		case 0: {
			for (i = 0; i < vertex_count; i++) {x[i] = v[i].b; y[i] = v[i].c;}
			tx = point.b; ty = point.c;
		}
		case 1: {
			for (i = 0; i < vertex_count; i++) {x[i] = v[i].a; y[i] = v[i].c;}
			tx = point.a; ty = point.c;
		}
		case 2: {
			for (i = 0; i < vertex_count; i++) {x[i] = v[i].a; y[i] = v[i].b;}
			tx = point.a; ty = point.b;
		}
	}

	// Use the even/odd method to test for point inclusion in polygon
	// Ref: http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
	
	c = 0;
	for (i = 0, j = vertex_count-1; i < vertex_count; j = i++) {
		if (
			( (y[i] > ty) != (y[j] > ty) ) &&
			( tx < (x[j]-x[i]) * (ty-y[i]) / (y[j]-y[i]) + x[i])
		) {
			c = !c;
		}
	}
	return c;	
}

num object_intersect_plane(tuple normal, tuple point, ray r) {
	num denom = geo_dot(normal, r.dir);
	if (! denom) return -1;
	num t = geo_dot( geo_sub(point, r.start), normal ) / denom ;
	return (t >= 0) ? t : -1;
}

num object_intersect_sphere(object_sphere* sphere, ray r) {
	tuple q = geo_sub(sphere->centre, r.start);
	num vdq = geo_dot(r.dir, q);
	num sds = geo_dot(q,q) - (sphere->radius * sphere->radius);
	num dsm = (vdq*vdq) - sds;
	if (dsm >= 0) {
		num root = sqrt(dsm);
		num t0   = vdq - root;
		num t1   = vdq + root;
		if (t0 > 0) return t0;
		if (t1 > 0) return t1;
	}
	return -1;
}

num object_intersect_triangle(object_triangle* tri, ray r) {
#ifndef FAST_TRIANGLES

	// Perform generic 3-polygon intersection
	
	num in_plane = object_intersect_plane(
		object_normal_polygon(tri->v), tri->v[0], r
	);
	tuple point = ray_pos(r, in_plane);	
	return object_intersect_poly(tri->v, 3, point) ? in_plane : -1;
	
#else
	// Apply Moeller-Trumbore algorithm for fast ray-triangle intersection
	// Ref: http://www.cs.virginia.edu/~gfx/Courses/2003/ImageSynthesis/papers/
	//      Acceleration/Fast%20MinimumStorage%20RayTriangle%20Intersection.pdf

	tuple e1 = geo_sub(tri->v[2],tri->v[0]);
	tuple e2 = geo_sub(tri->v[1],tri->v[0]);

	tuple p = geo_cross(r.dir, e2);
	num det = geo_dot(e1, p);

	if (det > -EPS && det < EPS) return -1;

	num invdet = 1.0 / det;
	tuple s = geo_sub(r.start,tri->v[0]);
	
	// Calculate intersection point
	num u = invdet * (geo_dot(s, p));
	if (u < 0.0 || u > 1.0) return -1;

	tuple q = geo_cross(s, e1);
	num v = invdet * geo_dot(r.dir, q);

	if (v < 0.0 || u + v > 1.0) return -1;
	
	num t = invdet * geo_dot(e2,q);
	return (t > EPS) ? t : -1;
	
#endif
}

num object_intersect_csg(object_meta_csg* csg, ray r, int *obj_id) {

	/* Get near and far intersections for the two objects */
	
	num n1 = object_intersect(csg->obj1, r);
	num f1 = -1;
	if (n1 != -1) {
		tuple p = ray_pos(r, n1);
		ray _r = NEW_RAY(p, r.dir);
		f1 = object_intersect(csg->obj1, _r);					
	}
	
	num n2 = object_intersect(csg->obj2, r);
	num f2 = -1;
	if (n2 != -1) {
		tuple p = ray_pos(r, n2);
		ray _r = NEW_RAY(p, r.dir);
		f2 = object_intersect(csg->obj2, _r);
	}
	
	/* Perform CSG operation */
	
	switch (csg->csg_operation) {
		case OBJ_CSG_UNION: {
			if (n1 == -1) {
				*obj_id = OBJ_HIT_SECOND | OBJ_HIT_NEAR;
				return n2;
			} else {
				if (n1 > n2) {
					*obj_id = OBJ_HIT_FIRST | OBJ_HIT_NEAR;
					return n1;
				} else {
					*obj_id = OBJ_HIT_SECOND | OBJ_HIT_NEAR;
					return n2;
				}
			}
		} break;
		
		case OBJ_CSG_SUBTRACT: {
			if (n1 == -1) {
				return -1;
			} else {
				if (f2 == -1) {
					*obj_id = OBJ_HIT_FIRST | OBJ_HIT_NEAR;
					return n1;
				} else {
					*obj_id = OBJ_HIT_SECOND | OBJ_HIT_FAR;
					return f2;
				}
			}
		} break;
		
		case OBJ_CSG_INTERSECTION: {
			if (n1 == -1) {
				return -1;
			} else {
				if (n2 < n1) {
					*obj_id = OBJ_HIT_FIRST | OBJ_HIT_NEAR;
					return n2;
				} else {
					*obj_id = OBJ_HIT_SECOND | OBJ_HIT_NEAR;
					return n1;
				}
			}
		} break;
	}
	
	/* Fallback */
#ifdef GEO_DEBUG
	fputs("ERROR: Invalid CSG operation\n", stderr);
#endif
	return n1;
}

num object_intersect(object_entry* obj, ray r ) {
	switch( obj->object_type ) {
	
		case OBJ_TYPE_SPHERE: {
			return object_intersect_sphere(obj->data, r);
		} break;
		
		case OBJ_TYPE_PLANE: {
			object_plane* plane = obj->data;
			num ret = object_intersect_plane(plane->normal, plane->point, r);
			return ret;
		} break;
	
		case OBJ_TYPE_QUAD: {
			object_quad* quad = obj->data;
			num in_plane = object_intersect_plane(
				object_normal_polygon(quad->v), quad->v[0], r
			);
			tuple point  = ray_pos(r, in_plane);
			
			return object_intersect_poly(quad->v, 4, point) ? in_plane : -1;
		} break;
		
		case OBJ_TYPE_TRIANGLE: {
			return object_intersect_triangle(obj->data, r);
		} break;
		
		case OBJ_TYPE_BOX: {
			object_box* box = obj->data;
			num near_intersection;
			if (object_intersect_box(&r, box, &near_intersection)) {
				return near_intersection;
			}
			return -1;
		} break;
		
		/* Meta objects */
		
		case OBJ_TYPE_CSG: {
			object_meta_csg* csg = obj->data;
			int object_id;
			return object_intersect_csg(csg, r, &object_id);
		} break;
		
		default: {
			return -1;
		} break;
	}
}

/* ###########################################################################
   COLOURS
   ###########################################################################*/


tuple object_textured_colour(object_entry* obj, num u, num v) {
	const int w = obj->mat.texture->width, h = obj->mat.texture->height;
	
	// x86 seems to allow negative modulo. That's okay with C89 + C99.
	// This non-branching correction should be faster than "if (x<0) x += w".
	
#ifndef ENABLE_BILINEAR_TEXTURES	
	int x = ((int)(u * w)) % w; x = (x + w) % w;
	int y = ((int)(v * h)) % h; y = (y + h) % h;
	
	return geo_from_u32(obj->mat.texture->pixels[x + (y*w)]);
	
#else

	// New method
	// Ref: http://www.devmaster.net/articles/raytracing_series/part6.php
	// n.b. This source has a bug: fv should be *m_Height, not *m_Width
	
	num fu = u * w;
	num fv = v * h;
	
	int u1 = ((int)fu) % w; u1 = (u1 + w) % w;
	int v1 = ((int)fv) % h;	v1 = (v1 + h) % h;	
	int u2 = (u1 + 1)  % w;
	int v2 = (v1 + 1)  % h;

	// Fetch four colours
	tuple c1 = geo_from_u32(obj->mat.texture->pixels[u1 + v1 * w]);
	tuple c2 = geo_from_u32(obj->mat.texture->pixels[u2 + v1 * w]);
	tuple c3 = geo_from_u32(obj->mat.texture->pixels[u1 + v2 * w]);
	tuple c4 = geo_from_u32(obj->mat.texture->pixels[u2 + v2 * w]);
	
	// Calculate weight factors
	num fracu = fu - floorf( fu );
	num fracv = fv - floorf( fv );
	num w1    = (1 - fracu) * (1 - fracv);
	num w2    = fracu * (1 - fracv);
	num w3    = (1 - fracu) * fracv;
	num w4    = fracu *  fracv;
	
	// Perform bilinear filtering by applying weight factors to colours
	return geo_add(geo_add(geo_add(
		geo_scale(c1, w1) , geo_scale(c2, w2)),
		geo_scale(c3, w3)), geo_scale(c4, w4))
	;
	
#endif
}
   
tuple object_diffuse_colour(object_entry* obj, tuple p) {
#ifdef ENABLE_TEXTURES
	if (obj->mat.texture) {
		//int w = obj->mat.texture->width, h = obj->mat.texture->height;
		switch( obj->object_type ) {
			
			case OBJ_TYPE_SPHERE: {
				object_sphere* sphere = obj->data;
			
				// Ref: http://www.cs.unc.edu/~rademach/xroads-RT/RTarticle.html
			
				tuple Vn = geo_new(0, -1, 0); // Pointing north
				tuple Ve = geo_new(1,  0, 0); // Pointing toward an equator		
				tuple Vp = geo_unit( geo_sub( p, sphere->centre ) );
			
				num phi  = acos( - geo_dot(Vn, Vp) );
				num v    = phi / PI;
				num thta = ( acos( geo_dot(Vp, Ve) / sin(phi) ) ) / (2 * PI);
				num u    = geo_dot( geo_cross( Vn, Ve), Vp) > 0 ? thta : 1-thta;

				return object_textured_colour(obj, u, v);
			
			} break;
			
			case OBJ_TYPE_PLANE: {				
				// New version
				// Ref: http://www.devmaster.net/articles/raytracing_series/part6.php
				
				object_plane* plane = obj->data;
				
				// Get orthogonal vectors in the plane
				tuple u_axis = geo_new(
					plane->normal.b, plane->normal.c, -plane->normal.a
				);
				tuple v_axis = geo_cross(u_axis, plane->normal);
				
				num u = geo_dot(p, u_axis) * 1.0;
				num v = geo_dot(p, v_axis) * 1.0;
				
				return object_textured_colour(obj, u, v);

			} break;
		}
	}
#endif	
	return obj->mat.diffuse_colour;
}

int object_search(
	/* In:  */ object_entry* items, int item_count, ray r,
	/* Out: */ object_intersection* ret
) {
	num ret_t = -1; int closest = 0;
	
	// Naive linear search
	
	int i, found = 0;
	for (i = 0; i < item_count; i++) {
		num t = object_intersect(&items[i], r);
		if (t != -1) {
			found = 1;
			if (ret_t == -1 || t < ret_t) {
				ret_t = t;
				closest = i;
			}
		}
	}
	
	if (found) {
		ret->p    = ray_pos(r,ret_t);
		ret->obj  = &items[closest];
#ifdef ENABLE_CSG_NORMALS
		ret->norm = object_normal( ret->obj, ret->p, r);
#else
		ret->norm = object_normal( ret->obj, ret->p );
#endif
		return 1;
	}
	return 0;
}


/* ###########################################################################
   CONSTRUCTORS
   ###########################################################################*/
   

object_entry* object_new_plane(material_t mat, tuple point, tuple normal) {
	object_plane* dat = malloc(sizeof(object_plane));
	dat->point  = point;
	dat->normal = geo_unit(normal);
	
	object_entry* ret = malloc(sizeof(object_entry));
	ret->object_type = OBJ_TYPE_PLANE;
	ret->data = dat;
	ret->mat = mat;
	
	return ret;
}
object_entry* object_new_triangle(material_t mat, tuple v1, tuple v2, tuple v3) {
	object_triangle* dat = malloc(sizeof(object_triangle));
	dat->v[0] = v1; dat->v[1] = v2; dat->v[2] = v3;
	
	object_entry* ret = malloc(sizeof(object_entry));
	ret->object_type = OBJ_TYPE_TRIANGLE;
	ret->data = dat;
	ret->mat = mat;
	
	return ret;
}
object_entry* object_new_quad(material_t mat, tuple v1, tuple v2, tuple v3, tuple v4) {
	object_quad* dat = malloc(sizeof(object_quad));
	dat->v[0] = v1; dat->v[1] = v2; dat->v[2] = v3; dat->v[3] = v4;
	
	object_entry* ret = malloc(sizeof(object_entry));
	ret->object_type = OBJ_TYPE_QUAD;
	ret->data = dat;
	ret->mat = mat;
	
	return ret;
}
object_entry* object_new_box(material_t mat, tuple centre, tuple dimensions) {
	object_box* dat = malloc(sizeof(object_box));
	dat->centre = centre;
	dat->dimensions = dimensions;
	
	object_entry* ret = malloc(sizeof(object_entry));
	ret->object_type = OBJ_TYPE_BOX;
	ret->data = dat;
	ret->mat = mat;
	
	return ret;
}
object_entry* object_new_sphere(material_t mat, tuple centre, num radius) {
	object_sphere* dat = malloc(sizeof(object_sphere));
	dat->centre = centre;
	dat->radius = radius;
	
	object_entry* ret = malloc(sizeof(object_entry));
	ret->object_type = OBJ_TYPE_SPHERE;
	ret->data = dat;
	ret->mat = mat;
	
	return ret;
}
object_entry* object_new_csg(material_t mat, int operation, object_entry* o1, object_entry* o2) {
	object_meta_csg* dat = malloc(sizeof(object_meta_csg));
	dat->csg_operation = operation;
	dat->obj1 = o1;
	dat->obj2 = o2;
	
	object_entry* ret = malloc(sizeof(object_entry));
	ret->object_type = OBJ_TYPE_CSG;
	ret->data = dat;
	ret->mat = mat;
	
	return ret;
}
void object_free(object_entry* oe) {
	if (oe->object_type == OBJ_TYPE_CSG) {
		object_meta_csg *csg = oe->data;
		free(csg->obj1);
		free(csg->obj2);
	}
	free(oe->data);
	free(oe);
}