/* 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. */

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

- LEsize is an unsigned integer type used
  for indexing vertices, triangles, etc.

- The model type uses indices to define triangles
  from vertices and looks something like this:
typedef struct {
	LEsize num_vertex, num_triangle, num_edge;
	float (*vertex)[4];
	LEsize (*indices)[3];
	struct {LEsize edge[3];} *triangle;
	struct {LEsize vertex[2];} *edge;
} Model;

- Could get triangle plane equation at load time,
  since it is invariant to the algorithm
*/

#define _X	[0]
#define _Y	[1]
#define _Z	[2]
#define _W	[3]

#ifndef DOT3
#	define DOT3(u, v) \
		( (u)[0] * (v)[0] \
		+ (u)[1] * (v)[1] \
		+ (u)[2] * (v)[2] )
#endif

#ifndef DOT4
#	define DOT4(u, v) \
		( (u)[0] * (v)[0] \
		+ (u)[1] * (v)[1] \
		+ (u)[2] * (v)[2] \
		+ (u)[3] * (v)[3] )
#endif

#ifndef CROSS3
#	define CROSS3(d, u, v) do { \
		(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]; } while(0)
#endif

#ifndef MOV4
#	define MOV4(d, s) do { \
		(d)[0] = (s)[0]; \
		(d)[1] = (s)[1]; \
		(d)[2] = (s)[2]; \
		(d)[3] = (s)[3]; } while(0)
#endif

/* precalculate all edges in mdl.
   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_triangle * 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_triangle; ++i) {
		/* each edge of that tri */
		for(j = 0; j < 3; ++j) {
			/* verts that make up this edge */
			a = mdl->indices[i][j];
			b = mdl->indices[i][(j + 1) % 3];
			/* edges are order-independent */
			if(a > b) {
				t = a; a = b; b = t;
			}
			/* do we already have this edge? */
			for(k = 0; k < mdl->num_edge; ++k) {
				if(mdl->edge[k].vertex[0] == a
				&& mdl->edge[k].vertex[1] == b) {
					break;
				}
			}
			/* ...no? 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;
			}
		}
	}
	/* assign each triangle's edge to one in the list */
/*	#pragma omp parallel for private(j, k, a, b, t)	*/
	for(i = 0; i < mdl->num_triangle; ++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->indices[i][j];
			b = mdl->indices[i][(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) {
					break;
				}
			}
			mdl->triangle[i].edge[j] = k;
		}
	}
	/* 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;
	}
	return 0;
}

/* get plane equation of a triangle; credit to Eric Lengyel */
static void GetPlane(float plane[4], float tri[3][4]) {
	float 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]);
}

/* extrude vertex 'src' away from 'origin' towards infinity */
static void ExtrudeVertex(float dst[4], float src[4], float origin[4]) {
	/* NOTE this may need to be another
	   value, depending on precision */
	const float ALMOST_ZERO = 0.0f;
	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;
}

/* Build shadow volume from 'mdl' 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. 'mdl' and 'light'
   need to be in the same space, preferably model local. */
static LEsize BuildVolume(float (*out_vertex)[4]
	, LEsize (*out_triangle)[3]
	, signed char *edge_count
	, Model *mdl, float light[4]
	, int caps) {
	float tri[3][4], ext[3][4], v[4];
	LEsize i, p, num = 0;
	unsigned j, a, b;
	if(out_vertex == NULL || out_triangle == NULL
	|| edge_count == NULL
	|| mdl == NULL || light == NULL) {
		return 0;
	}
	/* temp buffer, initialized to zero */
	memset(edge_count, 0, mdl->num_edge);
	/* each triangle may cast a shadow */
	for(i = 0; i < mdl->num_triangle; ++i) {
		/* fetch triangle */
		for(j = 0; j < 3; ++j) {
			p = mdl->indices[i][j];
			MOV4(tri[j], mdl->vertex[p]);
		}
		/* only consider tri when facing away from light */
		GetPlane(v, tri);
		if(DOT4(v, light) >= 0) {
			continue;
		}
		if(caps) {
			for(j = 0; j < 3; ++j) {
				/* cap at model: just the source tri */
				p = (num + 0) * 3 + j;
				MOV4(out_vertex[p], tri[j]);
				/* cap at inf needs to
				   be winded other way */
				p = (num + 1) * 3 + (2 - j);
				ExtrudeVertex(out_vertex[p]
					, tri[j], light);
			}
			num += 2;
		}
		/* update each of the triangle's edges */
		for(j = 0; j < 3; ++j) {
			a = mdl->indices[i][j];
			b = mdl->indices[i][(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;
		}
		/* fill in tri[] and ext[] */
		for(j = 0; j < 2; ++j) {
			p = mdl->edge[i].vertex[j];
			MOV4(tri[j], mdl->vertex[p]);
			ExtrudeVertex(ext[j], mdl->vertex[p], light);
		}
		/* edges facing one way */
		for(; edge_count[i] < 0; ++edge_count[i]) {
			MOV4(out_vertex[num * 3 + 0], tri[1]);
			MOV4(out_vertex[num * 3 + 1], ext[0]);
			MOV4(out_vertex[num * 3 + 2], tri[0]);
			++num;
			MOV4(out_vertex[num * 3 + 0], tri[1]);
			MOV4(out_vertex[num * 3 + 1], ext[1]);
			MOV4(out_vertex[num * 3 + 2], ext[0]);
			++num;
		}
		/* edges facing the other way */
		for(; edge_count[i] > 0; --edge_count[i]) {
			MOV4(out_vertex[num * 3 + 0], tri[0]);
			MOV4(out_vertex[num * 3 + 1], ext[1]);
			MOV4(out_vertex[num * 3 + 2], tri[1]);
			++num;
			MOV4(out_vertex[num * 3 + 0], tri[0]);
			MOV4(out_vertex[num * 3 + 1], ext[0]);
			MOV4(out_vertex[num * 3 + 2], ext[1]);
			++num;
		}
	}
	return num;
}

/* */
static LEsize BuildVolumeIndexed(float (*out_vertex)[4]
	, LEsize (*out_triangle)[3]
	, signed char *edge_count
	, Model *mdl, float light[4]
	, int caps) {
	float tri[3][4], v[4];
	LEsize i, p, q, num = 0;
	unsigned j, a, b;
	if(out_vertex == NULL || out_triangle == NULL
	|| edge_count == NULL
	|| mdl == NULL || light == NULL) {
		return 0;
	}
	/* temp buffer, initialized to zero */
	memset(edge_count, 0, mdl->num_edge);
	/* speculatively extrude each vertex */
	for(i = 0; i < mdl->num_vertex; ++i) {
		MOV4(out_vertex[i], mdl->vertex[i]);
		ExtrudeVertex(out_vertex[i + mdl->num_vertex]
			, mdl->vertex[i], light);
	}
	/* each triangle may cast a shadow */
	for(i = 0; i < mdl->num_triangle; ++i) {
		/* fetch triangle */
		for(j = 0; j < 3; ++j) {
			p = mdl->indices[i][j];
			MOV4(tri[j], mdl->vertex[p]);
		}
		/* only consider triangle when facing away from light */
		GetPlane(v, tri);
		if(DOT4(v, light) >= 0) {
			continue;
		}
		if(caps) {
			for(j = 0; j < 3; ++j) {
				p = mdl->indices[i][j];
				/* cap at model: just the source tri */
				out_triangle[num    ][    j] = p;
				/* cap at inf needs to
				   be winded other way */
				out_triangle[num + 1][2 - j]
					= p + mdl->num_vertex;
			}
			num += 2;
		}
		/* update each of the triangle's edges */
		for(j = 0; j < 3; ++j) {
			a = mdl->indices[i][j];
			b = mdl->indices[i][(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;
		}
		/* vertices that form this edge */
		p = mdl->edge[i].vertex[0];
		q = mdl->edge[i].vertex[1];
		/* edge is facing one way */
		for(; edge_count[i] < 0; ++edge_count[i]) {
			out_triangle[num][0] = q;
			out_triangle[num][1] = p + mdl->num_vertex;
			out_triangle[num][2] = p;
			++num;
			out_triangle[num][0] = q;
			out_triangle[num][1] = q + mdl->num_vertex;
			out_triangle[num][2] = p + mdl->num_vertex;
			++num;
		}
		/* edge is facing the other way */
		for(; edge_count[i] > 0; --edge_count[i]) {
			out_triangle[num][0] = p;
			out_triangle[num][1] = q + mdl->num_vertex;
			out_triangle[num][2] = q;
			++num;
			out_triangle[num][0] = p;
			out_triangle[num][1] = p + mdl->num_vertex;
			out_triangle[num][2] = q + mdl->num_vertex;
			++num;
		}
	}
	return num;
}