// n = (nx, ny, nz): normal of plane
// p = (px, py, pz): point on line
// q = (qx, qy, qz): point on plane
// d = (dx, dy, dz): direction of line
// num= -dot(n, q-p)
// den= dot(n, d)
// t = num/den
C_INLINE void compute_line_vs_plane_numden(
	int nx, int ny, int nz,
	int px, int py, int pz,
	int qx, int qy, int qz,
	int dx, int dy, int dz,
	int *num, int *den)
{
	int qpx= qx-px;
	int qpy= qy-py;
	int qpz= qz-pz;
	*num= -(nx*qpx + ny*qpy + nz*qpz);
	*den= nx*dx + ny*dy + nz*dz;
}

// Robust point-in-triangle test using exact boolean predicates only - no constructions which lose precision.
// AFAIK, it should be impossible for a ray to "fall through a crack" in a triangle mesh using this method.
static boolean hit_triangle(int *num, int *den, int p0x, int p0y, int dx, int dy, int cell_type, int half)
{
	int n= *num, d= *den;

	if(n<0)
	{
		*num= n= -n;
		*den= d= -d;
	}

	if(n<d) // t in [0,1) ?
	{
		int64 i64num= n, i64den= d;
		int64 testx= i64den*p0x + i64num*dx;
		uint64 WORLD_ONE_times_den= i64den<<WORLD_FRACTIONAL_BITS;
		if((uint64)testx < WORLD_ONE_times_den) // x in [0, 1) ?
		{
			int64 testy = i64den*p0y + i64num*dy;
			if((uint64)testy < WORLD_ONE_times_den) // y in [0, 1) ?
			{
				int desired_triangle_index= (2*cell_type)+half;
				// FIXME: possible off-by-one error in coordinate reflection here?
				int hit_triangle_index= cell_type ? ((testy<WORLD_ONE_times_den-testx) ? 2 : 3) : ((testy>testx) ? 0 : 1);
				if(hit_triangle_index == desired_triangle_index) // correct half of the cell ?
				{
					return TRUE;
				}
			}
		}
	}

	return FALSE;
}

// NOTE: x, y must already be bounds checked before this function is called!
static boolean mesh_cell_segment_intersection(
	int x, int y,
	const world_point3d *p0_world, const world_point3d *p1_world,
	world_point3d *out_intersection,
	world_vector3d *out_normal,
	int *out_t)
{
	// (x,y)    (x+1,y)
	//    a      b
	//     +----+
	//     |    |
	//     +----+
	//    c      d
	// (x, y+1) (x+1, y+1)

	struct mesh_cell *a= GET_MESH_CELL_UNSAFE(x, y);
	struct mesh_cell *b= a+1;
	struct mesh_cell *c= a+mesh_width;
	struct mesh_cell *d= b+mesh_width;

	int h0= GET_TRUE_CELL_HEIGHT(a);
	int hb= GET_TRUE_CELL_HEIGHT(b)-h0;
	int hc= GET_TRUE_CELL_HEIGHT(c)-h0;
	int hd= GET_TRUE_CELL_HEIGHT(d)-h0;

	int dx= p1_world->x-p0_world->x;
	int dy= p1_world->y-p0_world->y;
	int dz= p1_world->z-p0_world->z;

	// local coordinate origin for the intersection calculation.
	int world_z= h0;
	int world_x= INTEGER_TO_WORLD(x);
	int world_y= INTEGER_TO_WORLD(y);

	boolean hit_tri0, hit_tri1;
	int cell_type= (x^y)&1;
	int num0, den0, num1, den1;

	world_point3d p0; // line start point in local coordinates
	set_world_point3d(&p0, p0_world->x-world_x, p0_world->y-world_y, p0_world->z-world_z);
	//set_world_point3d(&p1, p1_world->x-world_x, p1_world->y-world_y, p1_world->z-world_z);

	if(cell_type)
	{
		// triangle 0: normal = {-hb, -hc, 1}
		// (0,0,0) (1,0,hb)
		// (0,1,hc)
		compute_line_vs_plane_numden(-hb, -hc, WORLD_ONE, p0.x, p0.y, p0.z, 0, 0, 0, dx, dy, dz, &num0, &den0);
		
		// triangle 1: normal = {hc-hd, hb-hd, 1}
		//          (1,0,hb)
		// (0,1,hc) (1,1,hd)
		compute_line_vs_plane_numden(hc-hd, hb-hd, WORLD_ONE, p0.x, p0.y, p0.z, WORLD_ONE, 0, hb, dx, dy, dz, &num1, &den1);
	}
	else
	{
		// triangle 0: normal = {hd-hc, -hc, 1}
		// (0,0,0)         
		// (0,1,hc) (1,1,hd)
		compute_line_vs_plane_numden(hd-hc, -hc, WORLD_ONE, p0.x, p0.y, p0.z, 0, 0, 0, dx, dy, dz, &num0, &den0);

		// triangle 1: normal = {-hb, hb-hd, 1}
		// (0,0,0) (1,0,hb)
		//         (1,1,hd)
		compute_line_vs_plane_numden(-hb, hb-hd, WORLD_ONE, p0.x, p0.y, p0.z, 0, 0, 0, dx, dy, dz, &num1, &den1);
	}

	hit_tri0= hit_triangle(&num0, &den0, p0.x, p0.y, dx, dy, cell_type, 0);
	hit_tri1= hit_triangle(&num1, &den1, p0.x, p0.y, dx, dy, cell_type, 1);

	// hit both; which one was hit first?
	if(hit_tri0 & hit_tri1)
	{
		boolean t1_less_than_t0= (uint64)num0*den1 < (uint64)num1*den0;
		if(t1_less_than_t0)
		{
			hit_tri0= FALSE;
		}
		else
		{
			hit_tri1= FALSE;
		}
	}

	if(!(hit_tri0 | hit_tri1))
	{
		return FALSE;
	}

	if(out_normal)
	{
		word normal= a->normal;
		int precalculated_normal_index= hit_tri1 ?
			(normal >> PRECALCULATED_NORMAL_BITS) & PRECALCULATED_NORMAL_MASK :
			 normal & PRECALCULATED_NORMAL_MASK;
		*out_normal= precalculated_normal_table[precalculated_normal_index].normal;
	}

	if(out_t || out_intersection)
	{
		int t;
		unsigned num, den;

		if(hit_tri0)
		{
			num= num0; den= den0;
		}
		else
		{
			num= num1; den= den1;
		}

		if(num > (1<<(LONG_BITS-1-WORLD_FRACTIONAL_BITS)))
		{
			num>>= WORLD_FRACTIONAL_BITS;
			den>>= WORLD_FRACTIONAL_BITS;
		}

		t= (num<<WORLD_FRACTIONAL_BITS)/den;

		if(out_t)
		{
			*out_t= t;
		}

		if(out_intersection)
		{
			// Now construct the (inexact) intersection point using finite arithmetic.
			// Idea: Construct a point which is directly above (not on) the triangle which was hit.
			// The exact position over the triangle does not matter all that much, as long as that requirement is satisfied.

			enum { DONT_MIRROR= 0, MIRROR=WORLD_ONE-1 };

			int triangle_index, fx, fy, fz, ub, lb, h0, hx, hy, mirrory, mirrorx;

			fx= p0.x + ((dx*t)>>WORLD_FRACTIONAL_BITS);
			fy= p0.y + ((dy*t)>>WORLD_FRACTIONAL_BITS);

			// the interpolated point might not even be inside the triangle due to precision issues, so pin it in there.
			fx= PIN(fx, 0, WORLD_ONE-1);

			triangle_index= 2*cell_type+hit_tri1;
			switch(triangle_index)
			{
			case 0:
				// y < 1 -x
				lb= 0; ub= WORLD_ONE-1-fx;
				h0= hc; hx= hd; hy= 0;
				mirrorx= DONT_MIRROR; mirrory= MIRROR;
				break;
			case 1:
				// y >= 1-x
				lb= WORLD_ONE-1-fx; ub= WORLD_ONE-1;
				h0= hb; hx= 0; hy= hd;
				mirrorx= MIRROR; mirrory= DONT_MIRROR;
				break;
			case 2:
				// y >= x
				lb= fx; ub= WORLD_ONE-1;
				h0= 0; hx= hb; hy= hc;
				mirrorx= DONT_MIRROR; mirrory= DONT_MIRROR;
				break;
			case 3:
				// y < x
				lb= 0; ub= fx-1;
				h0= hd; hx= hc; hy= hb;
				mirrorx= MIRROR; mirrory= MIRROR;
				break;
			}

			// also pin the y coord inside the triangle
			fy= PIN(fy, lb, ub);

			// calculate the height of the triangle at (fx, fy)
			fz= h0 + (((fx^mirrorx)*(hx-h0) + (fy^mirrory)*(hy-h0))>>WORLD_FRACTIONAL_BITS) + 1; // + 1 to put us above the triangle.

			out_intersection->x= world_x+fx;
			out_intersection->y= world_y+fy;
			out_intersection->z= world_z+fz;
		}
	}

	return TRUE;
}

// Build 411 - DRB - new function, handy for lots of things
boolean mesh_segment_intersection(
	const world_point3d *p0, const world_point3d *p1,
	world_point3d *out_intersection,
	world_vector3d *out_normal,
	int *out_t)
{
	boolean success= FALSE;

	int p0x= WORLD_TO_INTEGER(p0->x), p0y= WORLD_TO_INTEGER(p0->y);
	int p1x= WORLD_TO_INTEGER(p1->x), p1y= WORLD_TO_INTEGER(p1->y);

	if(p0x==p1x && p0y==p1y)
	{
		success= mesh_cell_segment_intersection(p0x, p0y, p0, p1, out_intersection, out_normal, out_t);
	}
	else
	{
		struct grid_iterator_2d cell;
		if(grid_traverse_2d_init_boundschecked(&cell, WORLD_FRACTIONAL_BITS,
			// FIXME: should this be -2 instead of -1? is mesh_width the number of cells, or vertices? (verts=cells+1)
			p0->x, p0->y, p1->x, p1->y, 0, 0, mesh_width-1, mesh_height-1))
		{
			do
			{
				success= mesh_cell_segment_intersection(cell.x, cell.y, p0, p1, out_intersection, out_normal, out_t);
			}
			while(!success && grid_traverse_2d_next_boundschecked(&cell));
		}
	}

	return success;
}