/* This program is free software. It comes without any warranty, to
 * the extent permitted by applicable law. You can redistribute it
 * and/or modify it under the terms of the Do What The Fuck You Want
 * To Public License, Version 2, as published by Sam Hocevar. See
 * http://www.wtfpl.net/ for more details. */

/*
This source code is known to work with a specific renderer but may be
unusable for *your* rendering system. Untested on major renderers!

The model type uses indices to define triangles from vertices
and looks something like this:
typedef struct {
	unsigned long num_edge;
	unsigned num_vert, num_tri;
	float (*vertex)[4];
	struct {
		unsigned vertex[3];
		unsigned long edge[3];
		...normals, colors, etc...
	} *triangle;
	struct {
		unsigned vertex[2];
	} *edge;
} Model;
*/

/* some helper macros */
#define _X	[0]
#define _Y	[1]
#define _Z	[2]
#define _W	[3]
#define DOT3(u, v) \
	 ((u)[0] * (v)[0] + (u)[1] * (v)[1] + (u)[2] * (v)[2])
#define DOT4(u, v) \
	 ((u)[0] * (v)[0] + (u)[1] * (v)[1] \
	+ (u)[2] * (v)[2] + (u)[3] * (v)[3])
#define CROSS3(d, u, v) { \
	(d)[0] = (u)[1] * (v)[2] - (u)[2] * (v)[1]; \
	(d)[1] = (u)[2] * (v)[0] - (u)[0] * (v)[2]; \
	(d)[2] = (u)[0] * (v)[1] - (u)[1] * (v)[0]; }
#define MOV4(d, s) { \
	(d)[0] = (s)[0]; \
	(d)[1] = (s)[1]; \
	(d)[2] = (s)[2]; \
	(d)[3] = (s)[3]; }

static void OutOfMemory(const char *s) {
	/* issue some error somehow */
}

/* precalculate all edges in 'mdl'. assume that mdl->edge is NULL.
   this function has bad time complexity. */
int CacheEdges(Model *mdl) {
	unsigned long k;
	unsigned i, j, a, b, t;
	if(mdl == NULL) {
		return 1;
	}
	/* maximum possible number of edges */
	mdl->edge = realloc(mdl->edge
		, sizeof(*mdl->edge) * mdl->num_tri * 3);
	if(mdl->edge == NULL) {
		OutOfMemory("mesh edges");
		return 1;
	}
	mdl->num_edge = 0;
	/* first, find each unique edge in model */
	for(i = 0; i < mdl->num_tri; ++i) {
		/* each edge of that tri */
		for(j = 0; j < 3; ++j) {
			/* vertices that make up this edge */
			a = mdl->triangle[i].vertex[j];
			b = mdl->triangle[i].vertex[(j + 1) % 3];
			/* edges are order-independent */
			if(a > b) {
				t = a; a = b; b = t;
			}
			/* edge already in list? */
			for(k = 0; k < mdl->num_edge; ++k) {
				if(mdl->edge[k].vertex[0] == a
				&& mdl->edge[k].vertex[1] == b) {
					break;
				}
			}
			/* not in list? push new edge */
			if(k == mdl->num_edge) {
				mdl->edge[mdl->num_edge].vertex[0] = a;
				mdl->edge[mdl->num_edge].vertex[1] = b;
				++mdl->num_edge;
			}
		}
	}
	/* shrink-to-fit edge list */
	mdl->edge = realloc(mdl->edge, sizeof(*mdl->edge) * mdl->num_edge);
	if(mdl->edge == NULL) {
		OutOfMemory("mesh edges/Shrink");
		return 1;
	}
	/* assign each triangle's edge to one in the list */
	for(i = 0; i < mdl->num_tri; ++i) {
		for(j = 0; j < 3; ++j) {
			/* make the compiler happy. we know this will get
			   overwritten with the correct value */
			mdl->triangle[i].edge[j] = 47;
			/* verts that make up this edge */
			a = mdl->triangle[i].vertex[j];
			b = mdl->triangle[i].vertex[(j + 1) % 3];
			/* edges are order-independent */
			if(a > b) {
				t = a; a = b; b = t;
			}
			/* find edge in edge list */
			for(k = 0; k < mdl->num_edge; ++k) {
				if(mdl->edge[k].vertex[0] == a
				&& mdl->edge[k].vertex[1] == b) {
					mdl->triangle[i].edge[j] = k;
					break;
				}
			}
		}
	}
	return 0;
}

/* get plane equation of a triangle */
static void GetPlane(double plane[4], double tri[3][4]) {
	double p[3], q[3];
	int i;
	for(i = 0; i < 3; ++i) {
		p[i] = tri[1][i] - tri[0][i];
		q[i] = tri[2][i] - tri[0][i];
	}
	CROSS3(plane, p, q);
	plane _W = -DOT3(plane, tri[0]);
}

/* multiply 'vector' with 'matrix' */
static void ApplyTransform(double dst[4]
	, const double vector[4]
	, const double matrix[16]) {
	int i;
	for(i = 0; i < 4; ++i) {
		dst[i] = matrix[i  ] * vector[0] + matrix[i+4 ] * vector[1]
		       + matrix[i+8] * vector[2] + matrix[i+12] * vector[3];
	}
}

/* extrude vertex 'src' away from 'origin' towards infinity */
static void ExtrudeVertex(double dst[4], double src[4], double origin[4]) {
	/* NOTE this could/should be another value for float */
	const double ALMOST_ZERO = 0.0;
	dst _X = src _X * origin _W - origin _X * src _W;
	dst _Y = src _Y * origin _W - origin _Y * src _W;
	dst _Z = src _Z * origin _W - origin _Z * src _W;
	dst _W = ALMOST_ZERO;
}

/* pack the vertices 'a', 'b', and 'c' into the triangle 'dst' */
static void CopyTriangle(double dst[3][4]
	, double a[4], double b[4], double c[4]) {
	MOV4(dst[0], a);
	MOV4(dst[1], b);
	MOV4(dst[2], c);
}

/* Transform 'mdl' by 'modelview' and 'projection', then build shadow volume
   from this mesh and 'light'. Store the result in 'dst', return number of
   triangles written to 'dst'. Assume 'dst' has enough space allocated to
   store the whole volume. */
static unsigned long BuildVolume(double (*dst)[3][4], Model *mdl
	, double modelview[16], double projection[16], double light[4]) {
	double tri[3][4], ext[3][4], v[4], w[4];
	signed char *edge_count;
	unsigned long i, p, num = 0;
	unsigned j, a, b;
	if(dst == NULL || mdl == NULL
	|| modelview == NULL || projection == NULL || light == NULL) {
		return 0;
	}
	/* temp buffer, initialized to zero */
	edge_count = calloc(mdl->num_edge, sizeof(*edge_count));
	/* each triangle may cast a shadow */
	for(i = 0; i < mdl->num_tri; ++i) {
		/* put each of the triangle's vertices in modelview */
		for(j = 0; j < 3; ++j) {
			p = mdl->triangle[i].vertex[j];
			MOV4(v, mdl->vertex[p]);
			ApplyTransform(tri[j], v, modelview);
		}
		/* only consider triangle when facing away from light */
		GetPlane(v, tri);
		if(DOT4(v, light) >= 0) {
			continue;
		}
#ifdef NEED_CAP_GEOMETRY
		/* caps */
		for(j = 0; j < 3; ++j) {
			/* already have triangle in modelview */
			MOV4(v, tri[j]);
			/* extruded */
			ExtrudeVertex(w, v, light);
			/* put it in projection */
			ApplyTransform(tri[j], v, projection);
			ApplyTransform(ext[j], w, projection);
			/* cap at model: just the source triangle */
			MOV4(dst[num  ][  j], tri[j]);
			/* cap at infinity needs to be winded other way */
			MOV4(dst[num+1][2-j], ext[j]);
		}
		num += 2;
#endif
		/* update each of the triangle's edges */
		for(j = 0; j < 3; ++j) {
			a = mdl->triangle[i].vertex[j];
			b = mdl->triangle[i].vertex[(j + 1) % 3];
			p = mdl->triangle[i].edge[j];
			if(a < b) {
				++edge_count[p];
			}
			else if(a > b) {
				--edge_count[p];
			}
		}
	}
	/* each possible edge */
	for(i = 0; i < mdl->num_edge; ++i) {
		/* early out */
		if(edge_count[i] == 0) {
			continue;
		}
		/* extrude vertices */
		for(j = 0; j < 2; ++j) {
			/* fetch this vertex */
			p = mdl->edge[i].vertex[j];
			MOV4(w, mdl->vertex[p]);
			/* put it in modelview */
			ApplyTransform(v, w, modelview);
			/* extruded */
			ExtrudeVertex(w, v, light);
			/* put it in projection */
			ApplyTransform(tri[j], v, projection);
			ApplyTransform(ext[j], w, projection);
		}
		/* edges facing one way */
		for(; edge_count[i] < 0; ++edge_count[i]) {
			CopyTriangle(dst[num++], tri[1], ext[0], tri[0]);
			CopyTriangle(dst[num++], tri[1], ext[1], ext[0]);
		}
		/* edges facing the other way */
		for(; edge_count[i] > 0; --edge_count[i]) {
			CopyTriangle(dst[num++], tri[0], ext[1], tri[1]);
			CopyTriangle(dst[num++], tri[0], ext[0], ext[1]);
		}
	}
	free(edge_count);
	return num;
}