/** @file sector-layer.cu.h implementation of GPU-side sector-layer methods, to
		be included only once in the entire program */

#include <float.h>
#include <math.h>

#include "sector-layer.h"

/** finds intersection of a line with a circle
		@returns true if there is an intersection and false if there is not
 */
static hostdevice__ bool intersect_line_circle
(fvec2 p0, fvec2 k, fvec2 c, float r, float *t0, float *t1) {
	fvec2 pr = p0 - c;
	// square equation: t^2 + 2 * p * t + q = 0
	float p = dot(k, pr), q = len2(pr) - r * r;
	float discr = p * p - q;
	if(discr < 0)
		return false;
	float sqrt_discr = sqrtf(discr);
	*t0 = -p - sqrt_discr; *t1 = -p + sqrt_discr;
	return true;
}  // intersect_line_circle

/** range of tubes to test, after clamping */
static inline hostdevice__ 
range intersect_range(float inv_center_dist, float t0, float t1) {
	if(t0 <= t1) {
		int start = (int)ceilf(t0 * inv_center_dist);
		int end = (int)floorf(t1 * inv_center_dist) + 1;
		return range::from(start, end);
	} else
		return range::empty();
}  // intersect_range

hostdevice__ range2 SectorLayer::intersect_ring_tubes
(fvec2 center, float radius) const {
	return intersect_ring(TUBE_CENTER_DIST, TUBE_INV_CENTER_DIST, TUBE_RADIUS, 
												center, radius);
}

hostdevice__ range2 SectorLayer::intersect_ring_skewlets
(fvec2 center, float radius) const {
	return intersect_ring(SKEWLET_BIN_CENTER_DIST, SKEWLET_BIN_INV_CENTER_DIST, 
												SKEWLET_BIN_RADIUS, center, radius);
}

hostdevice__ range2 SectorLayer::intersect_ring
(float center_dist, float inv_center_dist, float max_dist, fvec2 center, 
float radius) const {
	// empty range if the sector-layer is skewed
	if(is_skewed)
		return range2::empty();
	// fast test against both circles, using minimum and maximun square
	// distances
	// start with end points
	float tmax = float(straws.len() - 1) * center_dist;
	fvec2 pn = p0 + k * tmax;
	float dist2_p0 = dist2(p0, center), dist2_pn = dist2(pn, center);
	float min_dist2 = min(dist2_p0, dist2_pn);
	float max_dist2 = max(dist2_p0, dist2_pn);
	// add line distance if nearest line point is inside the circle
	float tcenter = dot(center - p0, k);
	fvec2 n = normal_to(k);
	if(0 <= tcenter && tcenter <= tmax) {
		// also add line distance
		float dist_line = dot(center - p0, n);
		float dist2_line = dist_line * dist_line;
		min_dist2 = min(min_dist2, dist2_line);
		max_dist2 = max(max_dist2, dist2_line);
	}
	float r_outer = radius + max_dist, r_inner = radius - max_dist;
	float r_outer2 = r_outer * r_outer, r_inner2 = r_inner * r_inner;

	// check if the range lines completely inside or outside the ring
	if(min_dist2 > r_outer2 || max_dist2 < r_inner2)
		return range2::empty();

	// else do full line-ring intersection
	float t0_outer, t1_outer, t0_inner, t1_inner;
	bool res_inner = intersect_line_circle
		(p0, k, center, radius - max_dist, &t0_inner, &t1_inner);
	bool res_outer = intersect_line_circle
		(p0, k, center, radius + max_dist, &t0_outer, &t1_outer);

	// res_outer is always true, while res-inner can be false; 
	float t00, t01, t10, t11; // coordinates of ranges
	if(res_inner) {
		// 2 ranges
		t00 = t0_outer; t01 = t0_inner;
		t10 = t1_inner; t11 = t1_outer;
	} else {
		// 1 range
		t00 = t0_outer; t01 = t1_outer;
		t10 = 1; t11 = 0;  // second range empty
	}

	// clamp ranges by [0, tmax]
	t00 = max(t00, 0.0f); t01 = min(t01, tmax);
	t10 = max(t10, 0.0f); t11 = min(t11, tmax);

	// get tube ranges
	range r1 = intersect_range(inv_center_dist, t00, t01);
	range r2 = intersect_range(inv_center_dist, t10, t11);
	// exclude common point from the range if it is there
	if(!r1.is_empty() && !r2.is_empty() && r1.last() == r2.first())
		r2.start++;
	// shift ranges to the starting tube
	r1 += straws.start;
	r2 += straws.start;
	return range2(r1, r2);
}  // intersect_ring

hostdevice__ int SectorLayer::point_bin(fvec2 pos) const {
	int bin = straws.first() + 
		(int)floorf(dot(pos - p0, k) * SKEWLET_BIN_INV_CENTER_DIST + 0.5f);
	// if(bin < straws.first() || bin > straws.last()) {
	// 	printf("bin %d outside of allowed range (%d, %d)\n", bin, straws.first(), straws.last());
	// }
	return bin;
}  // point_bin