#include <iostream>

#include "perspective.h"
#include "perspectivesandboxview.h"


PerspectiveSandBoxView::PerspectiveSandBoxView():
	persp_p0  (new View::Point( Vector2(-1, -1) )),
	persp_px  (new View::Point( Vector2( 1, -1) )),
	persp_py  (new View::Point( Vector2(-1,  1) )),
	persp_p1  (new View::Point( Vector2( 1,  1) )),
	bounds_p0 (new View::Point( Vector2(-2, -2) )),
	bounds_p1 (new View::Point( Vector2( 2,  2) ))
{
	transform.scale(50, 50);
	points.push_back(persp_p0);
	points.push_back(persp_px);
	points.push_back(persp_py);
	points.push_back(persp_p1);
	points.push_back(bounds_p0);
	points.push_back(bounds_p1);
}


static Vector2 calc_optimal_resolution(
	const Vector2 &ox,
	const Vector2 &oy )
{
	const Real a = ox.x * ox.x;
	const Real b = ox.y * ox.y;
	const Real c = oy.x * oy.x;
	const Real d = oy.y * oy.y;
	Real e = fabs(ox.x*oy.y - ox.y*oy.x);
	if (e < real_precision_sqr)
		return Vector2(); // paranoid check, e must be non-zero when matrix is invertible
	e = 1.0/e;
	
	const Real sum = a*d + b*c - real_precision_sqr;
	Vector2 scale;
	if (2*a*b >= sum && 2*c*d >= sum) {
		if (a*b < c*d) {
			scale.x = sqrt(2*b)*e;
			scale.y = sqrt(2*a)*e;
		} else {
			scale.x = sqrt(2*d)*e;
			scale.y = sqrt(2*c)*e;
		}
	} else
	if (2*a*b >= sum) {
		scale.x = sqrt(2*b)*e;
		scale.y = sqrt(2*a)*e;
	} else
	if (2*c*d >= sum) {
		scale.x = sqrt(2*d)*e;
		scale.y = sqrt(2*c)*e;
	} else {
		const Real dif = a*d - b*c;
		scale.x = sqrt(dif/(a - c))*e;
		scale.y = sqrt(dif/(d - b))*e;
	}
	
	return scale;
}


void PerspectiveSandBoxView::update_background(
	const Matrix3 &matrix,
	int width,
	int height )
{
	DataSurface surface(width, height);
	for(int r = 0; r < height; ++r) {
		for(int c = 0; c < width; ++c) {
			Vector3 v = matrix*Vector3(c, r, 1);
			if (v.z <= real_precision)
				continue;

			const Real k = 1/(v.z*v.z);
			const Vector2 ox = Vector2( v.z*matrix.m00 - matrix.m02*v.x,
										v.z*matrix.m10 - matrix.m12*v.x )*k;
			const Vector2 oy = Vector2( v.z*matrix.m01 - matrix.m02*v.y,
										v.z*matrix.m11 - matrix.m12*v.y )*k;
			const Vector2 resolution = calc_optimal_resolution(ox, oy)*100;
			if (resolution.x <= real_precision || resolution.y <= real_precision)
				continue;
			
			const int lx = (int)round(log2(resolution.x));
			const int ly = (int)round(log2(resolution.y));
			Color color;
			color.r = lx < 0 ? lx%2 + 1 : lx%2;
			color.g = 0;
			color.b = ly < 0 ? ly%2 + 1 : ly%2;
			color.a = 1;
			surface[r][c] = color;
		}
	}
	background = surface.to_cairo_surface();
}


void PerspectiveSandBoxView::update_background_best_perimeter(
	const Matrix3 &matrix,
	int width,
	int height )
{
	const int log_min = -10000;
	const int log_max =  10000;
	const Real log_base = log(1.001);
	
	std::vector<Real> map_pk_max(log_max - log_min + 2, 0);
	std::vector<Real> map_pk_min(log_max - log_min + 2, 1);
	std::vector<Real> map_pk(width*height, 0);
	std::vector<int> map_log_index(width*height, 0);
	
	
	// fill map
	for(int r = 0; r < height; ++r) {
		for(int c = 0; c < width; ++c) {
			Vector3 v = matrix*Vector3(c, r, 1);
			if (v.z <= real_precision)
				continue;
			
			Real l = round(log(v.z)/log_base);
			if (l-1 < log_min || l+1 > log_max)
				continue;
			int log_index = (int)l - log_min + 1;
			assert(log_index > 0 && log_index < (int)map_pk_max.size() && log_index < (int)map_pk_min.size());

			const Real k = 1/(v.z*v.z);
			const Vector2 ox = Vector2( v.z*matrix.m00 - matrix.m02*v.x,
										v.z*matrix.m10 - matrix.m12*v.x )*k;
			const Vector2 oy = Vector2( v.z*matrix.m01 - matrix.m02*v.y,
										v.z*matrix.m11 - matrix.m12*v.y )*k;
			
			const Real area = fabs(ox.x*oy.y - ox.y*oy.x);
			const Real optimal_side = sqrt(area);
			const Real optimal_perimeter = 4*optimal_side;
			const Real perimeter = (ox.length() + oy.length())*2;
			
			const Real pk = optimal_perimeter/perimeter;
			
			map_pk[r*width + c] = pk;
			map_log_index[r*width + c] = log_index;
			if (map_pk_max[log_index] < pk)
				map_pk_max[log_index] = pk;
			if (map_pk_min[log_index] > pk)
				map_pk_min[log_index] = pk;
		}
	}

	// paint
	DataSurface surface(width, height);
	for(int r = 0; r < height; ++r) {
		for(int c = 0; c < width; ++c) {
			const int log_index = map_log_index[r*width + c];
			assert(log_index >= 0 && log_index < (int)map_pk_max.size() && log_index < (int)map_pk_min.size());
			if (log_index) {
				const Real v = map_pk[r*width + c];
				const Real v0 = map_pk_min[log_index];
				const Real v1 = map_pk_max[log_index];
				if (v0 + real_precision < v1) {
					Real k = (v - v0)/(v1 - v0);
					const Real p = 0.5;
					k -= p;
					if (k > 0) {
						k /= 1 - p;
						k = pow(k, 20);
						surface[r][c] = Color(k, k, k, 1);
					}
				}
			}
		}
	}
	
	// paint centers
	Perspective::LayerList layers;
	Perspective::create_layers(layers, matrix.inverted(), IntPair2(IntVector2(), IntVector2(width, height)), 2);
	for(Perspective::LayerList::const_iterator i = layers.begin(); i != layers.end(); ++i)
		Perspective::paint_cross(surface, i->center);
	std::cout << "layers count: " << layers.size() << std::endl;

	background = surface.to_cairo_surface();
}


void
PerspectiveSandBoxView::draw_grid(
	const Cairo::RefPtr<Cairo::Context> &context,
	const Matrix &matrix,
	const Color &color )
{
	const int count = 100;
	const int subcount = 10;

	const Real ps = get_pixel_size();

	for(int i = -count; i < count; ++i) {
		for(int j = -count; j < count; ++j) {
			Vector4 src(i/(Real)subcount, j/(Real)subcount, 0, 1);
			Vector4 dst = matrix*src;
			Real w = dst.w;
			if (w > real_precision) {
				Vector2 pos(dst.x/w, dst.y/w);
				
				Real ka = clamp(1/w, 0, 1);
				
				Real a = color.a*ka;
				Real r = 4*ps;
				if (i) r *= 0.5;
				if (j) r *= 0.5;
				if (i % 10) r *= 0.75;
				if (j % 10) r *= 0.75;
				context->set_source_rgba(color.r, color.g, color.b, a);
				context->arc(pos.x, pos.y, r, 0, 2.0*M_PI);
				context->fill();
			}
		}
	}
}

void
PerspectiveSandBoxView::draw_line(
	const Cairo::RefPtr<Cairo::Context> &context,
	const Pair2 &bounds,
	Real a,
	Real b,
	Real c,
	const Color &color )
{
	// equatation of line is a*x + b*y = c

	int count = 0;
	Vector2 points[4];
	
	if (fabs(a) > real_precision) {
		Real x0 = (c - bounds.p0.y*b)/a;
		if ( x0 + real_precision >= bounds.p0.x
		  && x0 - real_precision <= bounds.p1.x )
			points[count++] = Vector2(x0, bounds.p0.y);
		
		Real x1 = (c - bounds.p1.y*b)/a;
		if ( x1 + real_precision >= bounds.p0.x
		  && x1 - real_precision <= bounds.p1.x )
			points[count++] = Vector2(x1, bounds.p1.y);
	}

	if (fabs(b) > real_precision) {
		Real y0 = (c - bounds.p0.x*a)/b;
		if ( y0 + real_precision >= bounds.p0.y
		  && y0 - real_precision <= bounds.p1.y )
			points[count++] = Vector2(bounds.p0.x, y0);
		
		Real y1 = (c - bounds.p1.x*a)/b;
		if ( y1 + real_precision >= bounds.p0.y
		  && y1 - real_precision <= bounds.p1.y )
			points[count++] = Vector2(bounds.p1.x, y1);
	}

	if (count >= 2) {
		context->set_source_rgba(color.r, color.g, color.b, color.a);
		context->move_to(points[0].x, points[0].y);
		context->line_to(points[1].x, points[1].y);
		context->stroke();
	}
}

void
PerspectiveSandBoxView::draw_subdivisions(
	const Cairo::RefPtr<Cairo::Context> &context,
	const Matrix &matrix,
	const Pair2 &bounds,
	const Color &color )
{
	const Real ps = get_pixel_size();
	
	Vector4 a = matrix.row_x();
	Vector4 b = matrix.row_y();
	Vector4 c = matrix.row_w();

	// calc coefficient for equatation of "horizontal" line: A*x + B*y = C + D/w
	//                     equatation of line of horizon is: A*x + B*y = C
	Real A = a.y*b.w - a.w*b.y;
	Real B = a.w*b.x - a.x*b.w;
	Real C = a.y*b.x - a.x*b.y;
	Real D = a.w*(b.x*c.y - b.y*c.x) + b.w*(a.y*c.x - a.x*c.y) - C*c.w;
	
	if (D < 0) {
		A = -A;
		B = -B;
		C = -C;
		D = -D;
	}

	Real hd = ps*sqrt(A*A + B*B);
	if (hd <= real_precision)
		return; // orthogonal projection - no prespective - no subdiviosions
	hd = D/hd;
	Real horizonw  = hd;
	Real horizonw2 = hd/4;

	Real maxw = -INFINITY, minw = INFINITY;
	Vector2 corners[] = {
		Vector2(bounds.p0.x, bounds.p0.y),
		Vector2(bounds.p1.x, bounds.p0.y),
		Vector2(bounds.p1.x, bounds.p1.y),
		Vector2(bounds.p0.x, bounds.p1.y) };
	for(int i = 0; i < 4; ++i) {
		Real w = A*corners[i].x + B*corners[i].y - C;
		if (fabs(w) > real_precision) {
			w = D/w;
			if (w < 0 || w > horizonw) w = horizonw;
			if (minw > w) minw = w;
			if (maxw < w) maxw = w;
		}
	}

	if (minw >= maxw - real_precision)
		return; // all bounds too thin
	Real maxw2 = maxw < horizonw2 ? maxw : horizonw2;

	// draw limits
	Color limits_color = color;
	limits_color.a *= 0.5;
	draw_line(context, get_bounds(), A, B, C + D/minw, limits_color);
	draw_line(context, get_bounds(), A, B, C + D/maxw, limits_color);

	// split
	int minlog = (int)ceil (log2(minw ) + real_precision);
	int maxlog = (int)floor(log2(maxw2) - real_precision);
	for(int i = minlog; i <= maxlog; ++i)
		draw_line(context, bounds, A, B, C + D/exp2(i), color);
}

void
PerspectiveSandBoxView::on_draw_view(const Cairo::RefPtr<Cairo::Context> &context) {
	context->save();

	const Real ps = get_pixel_size();
	context->set_line_width(ps);

	Vector2 p0 = persp_p0->position;
	Vector2 px = persp_px->position;
	Vector2 py = persp_py->position;
	Vector2 p1 = persp_p1->position;
	
	Pair2 bounds(bounds_p0->position, bounds_p1->position);
	
	
	// perspective matrix
	Matrix matrix;
	{
		Vector2 A = px - p1;
		Vector2 B = py - p1;
		Vector2 C = p0 + p1 - px - py;
		
		Real cw = A.y*B.x - A.x*B.y;
		Real aw = B.x*C.y - B.y*C.x;
		Real bw = A.y*C.x - A.x*C.y;
		
		if (cw < 0) {
			aw = -aw;
			bw = -bw;
			cw = -cw;
		}
		
		Vector2 c = p0*cw;
		Vector2 a = px*(cw + aw) - c;
		Vector2 b = py*(cw + bw) - c;

		matrix.row_x() = Vector4(a, 0, aw);
		matrix.row_y() = Vector4(b, 0, bw);
		matrix.row_w() = Vector4(c, 0, cw);
	}
	
	{ // update background
		Matrix4 m4 = transform_to_pixels()*matrix;
		Matrix3 m3(
			Vector3( m4.m00, m4.m01, m4.m03 ),
			Vector3( m4.m10, m4.m11, m4.m13 ),
			Vector3( m4.m30, m4.m31, m4.m33 ) );
		Matrix3 back_matrix = Perspective::normalize_matrix_by_det( m3.inverted() );
		update_background_best_perimeter(back_matrix, get_width(), get_height());
	}
	
	if (background) {
		context->save();
		context->transform(transform_from_pixels().to_cairo());
		context->set_source(background, 0, 0);
		context->paint();
		context->restore();
	}
	
	draw_grid(context, matrix, Color(0, 1, 0, 1));
	draw_subdivisions(context, matrix, bounds, Color(0, 0, 1, 1));

	// draw frames
	context->set_source_rgba(0, 0, 1, 0.75);
	context->move_to(p0.x, p0.y);
	context->line_to(px.x, px.y);
	context->line_to(p1.x, p1.y);
	context->line_to(py.x, py.y);
	context->close_path();
	context->stroke();

	context->move_to(bounds.p0.x, bounds.p0.y);
	context->line_to(bounds.p1.x, bounds.p0.y);
	context->line_to(bounds.p1.x, bounds.p1.y);
	context->line_to(bounds.p0.x, bounds.p1.y);
	context->close_path();
	context->stroke();
	
	context->restore();
}