/** @file sector-layer.cu implementation of host-side sector-layer functions */

#include <float.h>
#include <limits.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "data.h"
#include "sector-layer.h"
#include "vec.h"

// hostdevice__ SectorLayer::SectorLayer
// (int row, int sector, fvec2 p0, fvec2 k, range straws, bool is_skewed) {
// 	//this->row = row;
// 	//this->sector = sector;
// 	this->p0 = p0;
// 	this->k = k;
// 	//this->n = normal_to(k);
// 	this->straws = straws;
// 	this->is_skewed = is_skewed;
// }  // SectorLayer

SectorLayer SectorLayer::from_tubes
(layptr<TlTube, Tube> tubes, range straws) {
	SectorLayer sl;
	memset(&sl, 0, sizeof(SectorLayer));
	sl.straws = straws;
	Tube tube0 = tubes[straws.start];
	//sl.row = tube0.row;
	//sl.sector = tube0.sector;
	sl.p0 = tube0.pos.xy();
	sl.k = normalized(tubes[straws.start + 1].pos.xy() - sl.p0);
	//sl.n = normal_to(sl.k);
	sl.is_skewed = tube0.is_skewed();
	return sl;
}  // from_tubes

SectorLayer SectorLayer::from_skewlet_bins
(const SectorLayer *sls1, int stride, 
 layptr<TlTube, Tube> tubes, range bins) {
	// sector-layers from which to derive the sector-layer for skewed bins
	const SectorLayer* sls[2][2] = {{&sls1[0], &sls1[stride]}, 
															 {&sls1[2 * stride], &sls1[3 * stride]}};
	// initialize known fields
	SectorLayer sl;
	memset(&sl, 0, sizeof(SectorLayer));
	sl.is_skewed = false;  // otherwise, intersections reject immediately
	sl.straws = bins;
	// k: just average over all 4 layers, and then normalize
	sl.k = normalized(sls[0][0]->k + sls[0][1]->k + sls[1][0]->k + sls[1][1]->k);
	// printf("skewed layer k = {%lf, %lf}\n", (double)sls[0][0]->k.x,
	// 			 (double)sls[0][0]->k.y);
	// printf("bin layer k = {%lf, %lf}\n", (double)sl.k.x, (double)sl.k.y);
	// p0: first value comes from average over 4 layers, and it should line on the
	// midline
	fvec2 p0 = 0.25f * (sls[0][0]->p0 + sls[0][1]->p0 + sls[1][0]->p0 +
											sls[1][1]->p0);
	float tmin = FLT_MAX, tmax = -FLT_MAX;
	for(int idouble_layer = 0; idouble_layer < 2; idouble_layer++)
		for(int iskewed_layer = 0; iskewed_layer < 2; iskewed_layer++) {
			const SectorLayer &sl1 = *sls[idouble_layer][iskewed_layer];
			for(int itube = sl1.straws.first(); itube <= sl1.straws.last(); itube++) {
				Tube tube = tubes[itube];
				float tpmin = dot(tube.min_pos().xy() - p0, sl.k);
				float tpmax = dot(tube.max_pos().xy() - p0, sl.k);
				tmin = min(tmin, min(tpmin, tpmax));
				tmax = max(tmax, max(tpmin, tpmax));
			}
		}
	// adjust p0, so that tmin corresponds to position -0.5f * SKEWLET_BIN_CENTER_DIST
	sl.p0 = p0 + sl.k * (tmin + 0.5f * SKEWLET_BIN_CENTER_DIST);
	// printf("tmin = %lf, tmax = %lf\n", (double)tmin, (double)tmax);
	// printf("skewed layer p0 = {%lf, %lf}\n", (double)sls[0][0]->p0.x,
	// 			 (double)sls[0][0]->p0.y);
	// printf("bin layer p0 = {%lf, %lf}\n", (double)sl.p0.x, (double)sl.p0.y);
	return sl;
}  // from_skewlet_bins