/*
 * Copyright (c) 2009 Veselin Georgiev
 *
 * Cutting the bullshit:
 *  all sources here are free software; you can use them for everything you
 *  want. Just don't blame me if something goes wrong; if something happens
 *  (including, but not limited to:
 *     coredumps, formatted harddrives, ping timeouts, winamp sound skips,
 *     fried CPUs, sudden mother-in-law visits, fraglimits, viagra spam,
 *     barking dogs, earthquakes, frog rain, Windows Vista installations,
 *     anrieff playing the bass guitar, dangerous sporks, apocalypse,
 *     concrete donkey, or even Chuck Norris getting your butt kicked...
 *  ), then you know you have been warned. Four times.
 *
 * Btw, Chuck Norris now doesn't code neither Python, nor Ruby. He's hacking
 * with Makefiles.
 */
#include "Recognition/queryengine.h"
#include "Recognition/mtrandom.h"
#include "Recognition/utils.h"
#include <stdio.h>
#include <math.h>
#include <vector>
#include <algorithm>
#include <queue>
#include <cassert>
using namespace std;

MTRandom rnd(42);
const int BRUTEFORCE_THRESH = 850;
int CRUDE_GRID_SIZE = 256;
bool verybig = false;
bool xfinal = false;
int w, h;
float basic_edge;

inline bool inrange(int x, int y) { return (x >= 0 && x < w && y >= 0 && y < h); }
inline int clamp(float f) { int x = (int) round(f); if (x < 0) x = 0; if (x > 255) x = 255; return x; }

void solve_bruteforce(void)
{
	vector<pair<int, int> > points;
	for (int y = 0; y < h; y++)
		for (int x = 0; x < w; x++)
			points.push_back(make_pair(x, y));
	for (int i = 0; i < (int) points.size() - 1; i += 2) {
		q.query(points[i].first, points[i].second, points[i + 1].first, points[i + 1].second);
		if (q.time() > 9.5) break;
	}
	if (points.size() % 2)
		q.query(points[points.size() - 1].first, points[points.size() - 1].second);
	for (int i = 0; i < (w * h); i++) if (q.data[i] == -1) q.data[i] = 127;
	q.submit(q.data);
}

int a[MAXX * MAXY];
int known[MAXX * MAXY];
int untouchable[MAXX * MAXY];
int kw, kh, boxsz;

void resize_bilinear(void)
{
	for (int y = 0; y < h; y++) {
		for (int x = 0; x < w; x++) {
			int xl = x % boxsz;
			int yl = y % boxsz;
			if (q.data[y * w + x] != -1) {
				a[y * w + x] = q.data[y * w + x];
				continue;
			}
			int x0 = x - xl;
			int y0 = y - yl;
			int x1 = x0 + boxsz;
			int y1 = y0 + boxsz;
			int c00 = a[y0 * w + x0];
			int c10, c01, c11;
			if (y1 < h)
				c10 = a[y1 * w + x0];
			else
				c10 = c00;
			if (x1 < w)
				c01 = a[y0 * w + x1];
			else
				c01 = c00;
			if (x1 < w && y1 < h)
				c11 = a[y1 * w + x1];
			else
				c11 = c01;
			float fx = xl / (float) boxsz;
			float fy = yl / (float) boxsz;
			float res = c00 * (1.0f - fx) * (1.0f - fy)
			          + c10 * (1.0f - fx) * (      fy)
			          + c01 * (       fx) * (1.0f - fy)
			          + c11 * fx * fy;
			int pix = (int) round(res);
			if (pix < 0) pix = 0;
			if (pix > 255) pix = 255;
			a[y * w + x] = pix;
		}
	}
}

float get_known(int x, int y)
{
	if (x < 0) x = 0;
	if (x >= kw) x = kw - 1;
	if (y < 0) y = 0;
	if (y >= kh) y = kh - 1;
	return (float)known[y * kw + x];
}

float bicubic(float x)
{
	const float a = -0.5;
	x = fabs(x);
	if (x == 0) return 1.0f;
	if (x >= 2.0f) return 0.0f;
	if (x < 1.0f) {
		return (a + 2.0f) * x * x * x - (a + 3.0f) * x * x + 1.0f;
	} else {
		return a * x * x * x - 5.0f * a * x * x + 8.0f * a * x - 4.0f * a;
	}
}

static inline float sinc(float x)
{
	return sin(x) / x;
}

static float lanczos3(float x)
{
	if (fabs(x) > 3.0) return 0.0f;
	x *= M_PI;
	if (fabs(x) < 1e-16) return 1.0f;
	return sinc(x) * sinc(x/3.0f);
}

void resize_convolution(float (*weight_fn) (float))
{
	for (int y = 0; y < h; y++) {
		for (int x = 0; x < w; x++) {
			if (q.data[y * w + x] != -1) {
				a[y * w + x] = q.data[y * w + x];
				continue;
			}
			int kybase = y / boxsz;
			int kxbase = x / boxsz;
			float wsum = 0.0f;
			float f = 0.0f;
			
			for (int dy = kybase - 3; dy <= kybase + 4; dy++) {
				float wy = weight_fn((y - dy * boxsz) / (float) boxsz);
				if (wy == 0.0f) continue;
				for (int dx = kxbase - 3; dx <= kxbase + 4; dx++) {
					float wx = weight_fn((x - dx * boxsz) / (float) boxsz);
					if (wx == 0.0f) continue;
					float w = wx * wy;
					f += get_known(dx, dy) * w;
					wsum += w;
				}
			}
			int t = (int) round(f / wsum);
			if (t < 0) t = 0;
			if (t > 255) t = 255;
			a[y * w + x] = t;
		}
	}
}

int div_round_up(int a, int b)
{
	if (a % b == 0) return a / b;
	else return a / b + 1;
}

void submit_regular_grid(void)
{
	resize_bilinear();
	q.submit(a);
	resize_convolution(bicubic);
	q.submit(a);
	resize_convolution(lanczos3);
	q.submit(a);
}

void solve_regular_grid(void)
{
	boxsz = 1;
	while ((w / boxsz) * (h / boxsz) > BRUTEFORCE_THRESH) boxsz++;
	kw = div_round_up(w, boxsz);
	kh = div_round_up(h, boxsz);
	vector<pair<int, int> > points;
	for (int y = 0; y < h; y += boxsz) {
		for (int x = 0; x < w; x += boxsz) {
			points.push_back(make_pair(x, y));
		}
	}
	for (int i = 0; i < (int) points.size() - 1; i += 2) {
		q.query(points[i].first, points[i].second, points[i + 1].first, points[i + 1].second);
	}
	if (points.size() % 2)
		q.query(points[points.size() - 1].first, points[points.size() - 1].second);
	for (int i = 0; i < (int) points.size(); i++) {
		known[i] = a[points[i].second * w + points[i].first] = q.data[points[i].second * w + points[i].first];
	}
	submit_regular_grid();
#ifdef PRINT_DEBUGINFO
	printf("First submits after %d queries\n", q.get_queries());
#endif

	int counter = 0;
	while (q.time() < 9.8) {
		int ngen = 0;
		int x[2], y[2];
		int num_tries = 2000;
		while (num_tries > 0 && ngen < 2) {
			--num_tries;
			int tx = rnd.next() % w;
			int ty = rnd.next() % h;
			if (q.data[ty * w + tx] == -1) {
				x[ngen] = tx;
				y[ngen] = ty;
				ngen++;
			}
		}
		if (ngen == 1)
			q.query(x[0], y[0]);
		else if (ngen == 2)
			q.query(x[0], y[0], x[1], y[1]);
		else if (ngen == 0 && q.have_all()) break;
		if (counter / 50 != (counter + ngen)/50) {
			submit_regular_grid();
		}
		counter += ngen;
	}
	submit_regular_grid();
}

const float CENTRICITY_IMPORTANCE = 0.4f;

struct GridPoint {
	int p[2];
	int g[2];
	float gradient;
	float area;
	float centricity;
	float eval() const
	{
		return ((1.0f - CENTRICITY_IMPORTANCE) + CENTRICITY_IMPORTANCE * centricity) * gradient * sqr(area);
	}
};
float eval_centricity(float x, float y)
{
	x /= (w - 1);
	y /= (h - 1);
	float dp = sqrt(sqr(x - 0.5f) + sqr(y - 0.5f));
	return 1.0f - dp / 0.7071067f;
}

bool operator < (const GridPoint& a, const GridPoint& b)
{
	return a.eval() < b.eval();
}

vector<GridPoint> grid;
vector<GridPoint> sortedgrid;
int owned[MAXX * MAXY];
int gvalues[MAXX * MAXY];
int gpoints[MAXX * MAXY][2];
int taboo[MAXX * MAXY];
float xarea[MAXX * MAXY];
int gw, gh;

void get_best_grid_size(int a, int b)
{
	bool swapped = false;
	if (a < b) {
		swapped = true;
		swap(a, b);
	}
	int pgw, pgh;
	for (gw = 1; ; gw++) {
		float xsize = a / (float) gw;
		gh = (int) round(b / xsize);
		gh = max(1, gh);
		int samples_required = (gw + 1) * (gh + 1);
		if (samples_required > 2*CRUDE_GRID_SIZE)
			break;
		pgw = gw;
		pgh = gh;
	}
	gw = pgw;
	gh = pgh;
	float xsize = a / (float) gw;
	float ysize = b / (float) gh;
	basic_edge = sqrt(sqr(xsize) + sqr(ysize));
	if (swapped) swap(gw, gh);
}

GridPoint make_gridpoint(int basex, int basey, int gpx, int gpy, float area)
{
	int gx = basex /*+ rnd.next() % 2*/;
	if (gx >= w) gx = w - 1;
	int gy = basey /*+ rnd.next() % 2*/;
	if (gy >= h) gy = h - 1;
	//
	GridPoint gp;
	gp.p[0] = gx;
	gp.p[1] = gy;
	gp.g[0] = gpx;
	gp.g[1] = gpy;
	q.cached_query(gx, gy);
	gp.area = area;
	xarea[gy * w + gx] = area;
	gp.centricity = eval_centricity(gx, gy);
	return gp;
}

GridPoint make_gridpoint(int gx, int gy)
{
	int basex = (int) round((w - 1) * gx / (float) gw);
	int basey = (int) round((h - 1) * gy / (float) gh);
	return make_gridpoint(basex, basey, gx, gy, 1.0f);
}


bool make_subgridpoint(float fx, float fy, float area, GridPoint& result)
{
	int ix = (int) round(fx);
	int iy = (int) round(fy);
	if (ix >= 0 && ix < w && iy >= 0 && iy < h && q.data[iy * w + ix] == -1) {
		result = make_gridpoint(ix, iy, -1, -1, area);
		return true;
	} else {
		return false;
	}
}

inline float face(float x1, float y1, float x2, float y2, float x3, float y3)
{
	return x1 * y2 + x2 * y3 + x3 * y1 - x3 * y2 - x2 * y1 - x1 * y3;
}

bool inside(float x, float y, float x1, float y1, float x2, float y2, float x3, float y3)
{
	float t[6] = { x1, y1, x2, y2, x3, y3 };
	float area = fabs(face(x1, y1, x2, y2, x3, y3));
	float carea = 0.0f;
	for (int i = 0; i < 3; i++) {
		carea += fabs(face(x, y, t[2 * i], t[2 * i + 1], t[(2 * i + 2) % 6], t[(2 * i + 3) % 6]));
	}
	return (fabs(carea - area) < 1e-4);
}

int triangle_shading(float x, float y, float x1, float y1, float x2, float y2, float x3, float y3, float c1, float c2, float c3)
{
	float hx = x - x1;
	float hy = y - y1;
	float vx = x2 - x1;
	float vy = y2 - y1;
	float ux = x3 - x1;
	float uy = y3 - y1;
	float Dcr = vx * uy - vy * ux;
	if (fabs(Dcr) < 1e-6) {
		// assume points are on a single line
		float d1 = sqrt(sqr(x - x1) + sqr(y - y1));
		float d3 = sqrt(sqr(x - x3) + sqr(y - y3));
		float totald = d1 + d3;
		return clamp((c1 * d3 + c3 * d1) / totald);
	}
	float rDcr = 1.0f / Dcr;
	float p = (hx * uy - hy * ux) * rDcr;
	float q = (vx * hy - vy * hx) * rDcr;
	return clamp(c1 + (c2 - c1) * p + (c3 - c1) * q);
}

float consider(float a, float b, float x)
{
	if (fabs(a-b) < 1e-3) return 0.5f;
	if (a < b) {
		if (x < a) return 0.0f;
		if (x > b) return 1.0f;
		return (x - a) / (b - a);
	} else {
		if (x < b) return 0.0f;
		if (x > a) return 1.0f;
		return (x - b) / (a - b);
	}
}

enum {
	SUBMIT_BRIEF = 0,
	SUBMIT_FINAL
};

bool get_closest_gridpoints(int x, int y, float c[4][2], float colors[], int n)
{
	float closestDist = 1e20;
	int bi = -1;
	for (int i = 0; i < n; i++) {
		float dist = sqrt(sqr(x - c[i][0]) + sqr(y - c[i][1]));
		if (dist < closestDist) {
			closestDist = dist;
			bi = i;
		}
	}
	//
	int gx = (int) c[bi][0];
	int gy = (int) c[bi][1];
	int dx = gx <= x ? +1 : -1;
	int dy = gy <= y ? +1 : -1;
	for (int i = 1; i < 20; i++) {
		int nx = gx + i * dx;
		if (nx < 0 || nx >= w) continue;
		if (q.data[gy * w + nx] == -1) continue;
		for (int j = 1; j < 20; j++) {
			int ny = gy + j * dy;
			if (ny < 0 || ny >= h || q.data[ny * w + gx] == -1) continue;
			if (q.data[ny * w + nx] == -1) continue;
			int x0 = min(gx, nx);
			int x1 = max(gx, nx);
			int y0 = min(gy, ny);
			int y1 = max(gy, ny);
			if (x >= x0 && x <= x1 && y >= y0 && y <= y1) {
				bool good = true;
				for (int k = y0 + 1; k < y1; k++) {
					for (int l = x0 + 1; l < x1; l++) {
						if (q.data[k * w + l] != -1) {
							good = false;
							break;
						}
					}
					if (!good) break;
				}
				if (!good) continue;
				c[0][0] = x0;
				c[0][1] = y0;
				c[1][0] = x1;
				c[1][1] = y0;
				c[2][0] = x0;
				c[2][1] = y1;
				c[3][0] = x1;
				c[3][1] = y1;
				colors[0] = q.data[y0 * w + x0];
				colors[1] = q.data[y0 * w + x1];
				colors[2] = q.data[y1 * w + x0];
				colors[3] = q.data[y1 * w + x1];
				return true;
			}
		}
	}
	return false;
}

void preshade(void)
{
	memset(untouchable, 0, sizeof(untouchable));
	for (int j = 0; j < gh; j++) {
		for (int i = 0; i < gw; i++) {
			int idx = j * (gw + 1) + i;
			int x0 = sortedgrid[idx         ].p[0];
			int y0 = sortedgrid[idx         ].p[1];
			int x1 = sortedgrid[idx + gw + 2].p[0];
			int y1 = sortedgrid[idx + gw + 2].p[1];
			assert(x0 < x1 && y0 < y1);
			int count_known = 0;
			for (int y = y0; y <= y1; y++)
				for (int x = x0; x <= x1; x++)
					if (q.data[y * w + x] != -1) count_known++;
			if (count_known != 4) continue;
			for (int y = y0; y <= y1; y++) {
				for (int x = x0; x <= x1; x++) {
					if (q.data[y * w + x] != -1) continue;
					untouchable[y * w + x] = 1;
					float u = (x - x0) / (float) (x1 - x0);
					float v = (y - y0) / (float) (y1 - y0);
					a[y * w + x] = clamp(
						+ q.data[y0 * w + x0] * (1.0f - u) * (1.0f - v)
						+ q.data[y0 * w + x1] * (       u) * (1.0f - v)
						+ q.data[y1 * w + x0] * (1.0f - u) * (       v)
						+ q.data[y1 * w + x1] * (       u) * (       v)
					);
				}
			}
		}
	}
}

void shade_method1(void)
{
	#define NCELLS 215
	const int nearest_neighs[NCELLS][2] = {
		{0,1},{1,0},{1,1},{0,2},{2,0},{1,2},{2,1},{2,2},{0,3},{3,0},{3,1},{1,3},{3,2},{2,3},{0,4},
		{4,0},{1,4},{4,1},{3,3},{2,4},{4,2},{3,4},{0,5},{5,0},{4,3},{1,5},{5,1},{5,2},{2,5},{4,4},
		{3,5},{5,3},{0,6},{6,0},{6,1},{1,6},{6,2},{2,6},{5,4},{4,5},{3,6},{6,3},{0,7},{7,0},{7,1},
		{1,7},{5,5},{6,4},{4,6},{7,2},{2,7},{7,3},{3,7},{6,5},{5,6},{8,0},{0,8},{7,4},{1,8},{4,7},
		{8,1},{2,8},{8,2},{6,6},{3,8},{8,3},{7,5},{5,7},{4,8},{8,4},{0,9},{9,0},{1,9},{9,1},{7,6},
		{6,7},{2,9},{9,2},{5,8},{8,5},{3,9},{9,3},{4,9},{9,4},{7,7},{6,8},{0,10},{10,0},{8,6},{1,10},
		{10,1},{2,10},{10,2},{5,9},{9,5},{3,10},{10,3},{7,8},{8,7},{10,4},{4,10},{6,9},{9,6},{0,11},{11,0},
		{1,11},{11,1},{5,10},{2,11},{11,2},{10,5},{8,8},{3,11},{7,9},{11,3},{9,7},{6,10},{10,6},{11,4},{4,11},
		{0,12},{12,0},{1,12},{9,8},{8,9},{12,1},{11,5},{5,11},{12,2},{2,12},{7,10},{10,7},{3,12},{12,3},{6,11},
		{11,6},{12,4},{4,12},{9,9},{10,8},{8,10},{0,13},{13,0},{5,12},{12,5},{11,7},{13,1},{7,11},{1,13},{2,13},
		{13,2},{3,13},{13,3},{6,12},{12,6},{10,9},{9,10},{11,8},{4,13},{8,11},{13,4},{7,12},{12,7},{13,5},{5,13},
		{0,14},{14,0},{14,1},{1,14},{10,10},{2,14},{14,2},{9,11},{11,9},{3,14},{14,3},{6,13},{13,6},{8,12},{12,8},
		{4,14},{14,4},{7,13},{13,7},{14,5},{10,11},{11,10},{5,14},{12,9},{0,15},{9,12},{15,0},{15,1},{1,15},{2,15},
		{15,2},{6,14},{14,6},{8,13},{13,8},{3,15},{15,3},{15,4},{4,15},{11,11},{10,12},{12,10},{14,7},{7,14},{5,15},
		{13,9},{15,5},{9,13},{0,16},{16,0},
	};
	for (int y = 0; y < h; y++) {
		for (int x = 0; x < w; x++) {
			if (q.data[y * w + x] != -1) {
				a[y * w + x] = q.data[y * w + x];
				continue;
			}
			if (untouchable[y * w + x]) continue;
			float npoints[4][2];
			int nc = 0;
			float colors[4], dists[4], sumdist = 0.0;
			for (int quadrant = 0; quadrant < 4; quadrant++) {
				int xmul = (quadrant & 1) * 2 - 1;
				int ymul = ((quadrant >> 1) & 1) * 2 - 1;
				int k, nx, ny;
				for (k = 0; k < NCELLS; k++) {
					nx = x + xmul * nearest_neighs[k][0];
					ny = y + ymul * nearest_neighs[k][1];
					if (nx >= 0 && nx < w && ny >= 0 && ny < h && q.data[ny * w + nx] != -1)
						break;
				}
				if (k >= NCELLS) continue;
				npoints[nc][0] = nx;
				npoints[nc][1] = ny;
				float fx = nx, fy = ny;
				colors[nc] = q.data[ny * w + nx];
				dists[nc] = sqrt(sqr(fx - x) + sqr(fy - y));
				nc++;
			}
			if (!get_closest_gridpoints(x, y, npoints, colors, nc)) {
				for (int i = nc - 1; i > 0; i--) {
					for (int j = i - 1; j >= 0; j--) {
						if (npoints[i][0] == npoints[j][0] && npoints[i][1] == npoints[j][1]) {
							for (int k = i; k < nc - 1; k++) {
								npoints[k][0] = npoints[k + 1][0];
								npoints[k][1] = npoints[k + 1][1];
								dists[k] = dists[k + 1];
								colors[k] = colors[k + 1];
							}
							nc--;
							break;
						}
					}
				}
			}
			for (int i = 0; i < nc; i++)
				sumdist += dists[i];
			//
			switch (nc) {
				case 0:
				case 1:
				default:
#ifdef PRINT_DEBUGINFO
					printf("Bad case in submit_adaptive_grid(): %d\n", nc);
#endif
					a[y * w + x] = clamp(colors[0]);
					break;
				case 2:
				{
					a[y * w + x] = clamp(colors[0] * dists[1] / sumdist + colors[1] * dists[0] / sumdist);
					break;
				}
				case 3:
				{
					a[y * w + x] = triangle_shading(
						x, y,
						npoints[0][0], npoints[0][1],
						npoints[1][0], npoints[1][1],
						npoints[2][0], npoints[2][1],
						colors[0], colors[1], colors[2]
					);
					break;
				}
				case 4:
				{
					float u = consider(npoints[0][0], npoints[1][0], x);
					float v = consider(npoints[0][1], npoints[2][1], y);
					a[y * w + x] = clamp(
						+ (1.0f - u) * (1.0f - v) * colors[0]
						+ (       u) * (1.0f - v) * colors[1]
						+ (1.0f - u) * (       v) * colors[2]
						+ (       u) * (       v) * colors[3]);
					break;
				}
			}
		}
	}
}

static inline bool qgood(int x, int y)
{
	return (x >= 0 && x < w && y >= 0 && y < h && q.data[y * w + x] != -1);
}

void shade_method2(void)
{
	float* colors = new float[w*h];
	float* weights = new float[w*h];
	memset(colors, 0, w * h * sizeof(float));
	memset(weights, 0, w * h * sizeof(float));
	
	for (int y = 0; y < h; y++) {
		for (int x = 0; x < w; x++) {
			if (q.data[y * w + x] == -1)
				continue;
			float area = xarea[y * w + x];
			if (area == 0.0f) {
				area = 1.0f;
			}
			area *= basic_edge;
			int xa = (int) ceil(area);
			for (int zy = y - xa; zy <= y + xa; zy++) {
				if (!(zy >= 0 && zy < h)) continue;
				for (int zx = x - xa; zx <= x + xa; zx++) {
					if (!(zx >= 0 && zx < w)) continue;
					float d = sqrt((float) (sqr(x - zx) + sqr(y - zy)));
					if (d > area) continue;
					float f = d/area;
					f = sqr(1.0f - f) / (area * area);
					weights[zy * w + zx] += f;
					colors[zy * w + zx] += f * q.data[y * w + x];
				}
			}
		}
	}

	for (int y = 0; y < h; y++) {
		for (int x = 0; x < w; x++) {
			if (untouchable[y * w + x]) continue;
			if (q.data[y * w + x] != -1) {
				a[y * w + x] = q.data[y * w + x];
			} else {
				a[y * w + x] = 127;
				// final estimate based on neighbour summing
				if (weights[y * w + x] > 0)
					a[y * w + x] = clamp(colors[y * w + x] / weights[y * w + x]);
				
				// see if we can spot nearby grids:
				int grid_sizes[4][2] = { { 2, 2 }, { 2, 3 }, { 3, 2 }, { 3, 3 } };
				bool found = false;
				for (int k = 0; !found && k < 4; k++) {
					int gsx = grid_sizes[k][0];
					int gsy = grid_sizes[k][1];
					for (int offsx = 1; !found && offsx < gsx; offsx++) {
						for (int offsy = 1; !found && offsy < gsy; offsy++) {
							// check xdir:
							if (qgood(x - offsx, y) && qgood(x + gsx - offsx, y)) {
								float u = offsx / (float) gsx;
								a[y * w + x] = clamp(
									+ q.data[y * w + x - offsx] * (1.0f - u)
									+ q.data[y * w + x + gsx - offsx] * u
								);
								found = true; break;
							}
							// check ydir:
							if (qgood(x, y - offsy) && qgood(x, y + gsy - offsy)) {
								float v = offsy / (float) gsy;
								a[y * w + x] = clamp(
									+ q.data[(y - offsy) * w + x] * (1.0f - v)
									+ q.data[(y + gsy - offsy) * w + x] * v
								);
								found = true; break;
							}
							// check grid:
							if (qgood(x - offsx, y - offsy) && qgood(x + gsx - offsx, y - offsy) &&
							    qgood(x - offsx, y + gsy - offsy) && qgood(x + gsx - offsx, y + gsy - offsy))
							{
								float u = offsx / (float) gsx;
								float v = offsy / (float) gsy;
								a[y * w + x] = clamp(
									+ q.data[(y - offsy) * w + x - offsx] * (1.0f - u) * (1.0f - v)
									+ q.data[(y - offsy) * w + x + gsx - offsx] * u * (1.0f - v)
									+ q.data[(y + gsy - offsy) * w + x - offsx] * (1.0f - u) * v
									+ q.data[(y + gsy - offsy) * w + x + gsx - offsx] * u * v
								);
								found = true;
								break;
							}
						}
					}
				}
			}
		}
	}
	delete [] colors;
	delete [] weights;
}

void submit_adaptive_grid(int submit_type)
{
	//preshade();
	//shade_method1();
//	shade_method2();
	int * temp = new int[w*h];
	for (int shade_method = 0; shade_method < 3; shade_method++) {
		if (submit_type == SUBMIT_BRIEF && shade_method != 2) continue;
		switch (shade_method) {
			case 0:
				memset(untouchable, 0, sizeof(untouchable));
				shade_method1();
				break;
			case 1:
				preshade();
				shade_method2();
				break;
			
			case 2:
				memset(untouchable, 0, sizeof(untouchable));
				shade_method2();
				break;
		}
		if (submit_type != SUBMIT_BRIEF) q.submit(a);
		struct BlurKernel {
			int size;
			int kernel[8][8];
		};
		const BlurKernel blur_kernels[3] = {
			{
				3, {
					{ 1, 2, 1 },
					{ 2, 4, 2 },
					{ 1, 2, 1 },
				}
			},
			{
				3, {
					{ 1, 1, 1 },
					{ 1, 1, 1 },
					{ 1, 1, 1 },
				}
			},
			{
				5, {
				// gaussian blur with sigma = 1.0
					{ 0, 1, 2, 1, 0, },
					{ 1, 7, 11, 7, 1, },
					{ 2, 11, 19, 11, 2, },
					{ 1, 7, 11, 7, 1, },
					{ 0, 1, 2, 1, 0, },
				}
			},
		};
		for (int ker = 0; ker < (int)COUNT_OF(blur_kernels); ker++) {
			if (submit_type == SUBMIT_BRIEF && ker > 0) break;
			for (int y = 0; y < h; y++) {
				for (int x = 0; x < w; x++) {
					if (q.data[y * w + x] != -1) {
						temp[y * w + x] = q.data[y * w + x];
						continue;
					}
					int tot_weight = 0, tot_color = 0;
					for (int dy = 0; dy < blur_kernels[ker].size; dy++) {
						for (int dx = 0; dx < blur_kernels[ker].size; dx++) {
							int nx = x + dx - blur_kernels[ker].size / 2;
							int ny = y + dy - blur_kernels[ker].size / 2;
							if (nx < 0 || nx >= w || ny < 0 || ny >= h) continue;
							int weight = blur_kernels[ker].kernel[dy][dx];
							//if (q.data[ny * w + nx] != -1) weight *= 3;
							tot_color += a[ny * w + nx] * weight;
							tot_weight += weight;
						}
					}
					temp[y * w + x] = clamp((float) tot_color / (float) tot_weight);
				}
			}
			q.submit(temp);
		}
	}
	delete [] temp;
}

void calculate_gradient(GridPoint& gp)
{
	// can we apply sobel?
	int x = gp.g[0], y = gp.g[1];
#if 0
	if (x > 0 && x < gw && y > 0 && y < gh) {
		int h_kernel[3][3] = {
			{ -1, 0, +1 },
			{ -2, 0, +2 },
			{ -1, 0, +1 },
		};
		int v_kernel[3][3] = {
			{ +1, +2, +1 },
			{  0,  0,  0 },
			{ -1, -2, -1 },
		};
		float sum_h = 0.0f, sum_v = 0.0f;
		for (int dy = 0; dy < 3; dy++) {
			for (int dx = 0; dx < 3; dx++) {
				float val = gvalues[(y + dy - 1) * (gw + 1) + x + dx - 1];
				sum_h += val * h_kernel[dy][dx];
				sum_v += val * v_kernel[dy][dx];
			}
		}
		gp.gradient = sqrt(sqr(sum_h) + sqr(sum_v)) ; // sobel comes preamplified, so normalize it here.
		return;
	}
#endif
	// shit, we're at the end of something; use laplace-style differences
	float sum_h = 0.0f, sum_v = 0.0f;
	int i = y * (gw + 1) + x;
#if 0
	float h0 = 0.0f, h1 = 0.0f, v0 = 0.0f, v1 = 0.0f;
	if (x > 0)  h0 = gpoints[i][0] - gpoints[i - 1][0];
	if (x < gw) h1 = gpoints[i + 1][0] - gpoints[i][0];
	if (y > 0)  v0 = gpoints[i][1] - gpoints[i - gw - 1][1];
	if (y < gh) v1 = gpoints[i + gw + 1][1] - gpoints[i][1];
#else
	float h0 = 1.0f, h1 = 1.0f, v0 = 1.0f, v1 = 1.0f;
#endif
	if (x > 0 && x < gw) {
		// central difference:
		//sum_h = ((gvalues[i + 1] + gvalues[i]) / h1 - (gvalues[i] + gvalues[i - 1]) / h0) / 2.0f;
		sum_h = (gvalues[i + 1] - gvalues[i - 1]) / (h0 + h1);
	} else {
		if (x > 0)
			sum_h = (gvalues[i - 1] - gvalues[i]) / h0;
		else
			sum_h = (gvalues[i + 1] - gvalues[i]) / h1;
	}
	if (y > 0 && y < gh) {
		// central difference:
		//sum_v = ((gvalues[i + gw + 1] + gvalues[i]) / v1 - (gvalues[i] + gvalues[i - gw - 1]) / v0) / 2.0f;
		sum_v = (gvalues[i + gw + 1] - gvalues[i - gw - 1]) / (v0 + v1);
	} else {
		if (y > 0)
			sum_v = (gvalues[i - gw - 1] - gvalues[i]) / v0;
		else
			sum_v = (gvalues[i + gw + 1] - gvalues[i]) / v1;
	}
	gp.gradient = sqrt(sqr(sum_h) + sqr(sum_v));
	//gp.gradient = abs(sum_h) + abs(sum_v);
}

int search_neighbour(int sx, int sy, int dx, int dy)
{
	for (int i = 1; i < 10; i++) {
		int x = sx + i * dx;
		int y = sy + i * dy;
		if (!inrange(x, y)) break;;
		if (q.data[y * w + x] != -1)
			return q.data[y * w + x];
	}
	for (int i = 1; i < 10; i++) {
		int U = -1, D = -1;
		int fui = 0, fdi = 0;
		if (!inrange(sx + i * dx, sy + i * dy)) break;
		for (int j = 1; j < 7; j++) {
			int x = sx + i * dx - j * dy;
			int y = sy + i * dy + j * dx;
			if (inrange(x, y) && q.data[y * w + x] != -1) {
				U = q.data[y * w + x];
				fui = j;
				break;
			}
		}
		for (int j = 1; j < 7; j++) {
			int x = sx + i * dx + j * dy;
			int y = sy + i * dy - j * dx;
			if (inrange(x, y) && q.data[y * w + x] != -1) {
				D = q.data[y * w + x];
				fdi = j;
				break;
			}
		}
		if (U != -1 && D != -1) {
			int all = fui + fdi;
			float c = (U * fdi + D * fui) / (float) all;
			return clamp(c);
		}
	}
	return -1;
}

void recalculate_gradient(GridPoint& gp)
{
	int x = gp.p[0];
	int y = gp.p[1];
	int L = search_neighbour(x, y, -1,  0);
	int R = search_neighbour(x, y, +1,  0);
	int U = search_neighbour(x, y,  0, -1);
	int D = search_neighbour(x, y,  0, +1);
	int X = q.data[y * w + x];
	float sum_h = 0.0f, sum_v = 0.0f;
	switch ((L != -1) + (R != -1) * 2) {
		case 1:
			sum_h = 2 * (X - L);
			break;
		case 2:
			sum_h = 2 * (X - R);
			break;
		case 3:
			sum_h = (R - L);
			break;
	}
	switch ((U != -1) + (D != -1) * 2) {
		case 1:
			sum_v = 2 * (X - U);
			break;
		case 2:
			sum_v = 2 * (X - D);
			break;
		case 3:
			sum_v = (U - D);
			break;
	}
	gp.gradient = sqrt(sqr(sum_h) + sqr(sum_v));
}

bool cmp_gp(const GridPoint& a, const GridPoint& b)
{
	if (a.g[1] != b.g[1]) return a.g[1] < b.g[1];
	return a.g[0] < b.g[0];
}

void solve_adaptive_grid(void)
{
	int num_pixels = w * h;
	if (num_pixels < 92 * 128)
		CRUDE_GRID_SIZE = 192;
	else
		CRUDE_GRID_SIZE = 192 + 64 * (num_pixels - 92 * 128) / (36 * 128);
#ifdef PRINT_DEBUGINFO
	printf("Selected CRUDE_GRID_SIZE of %d\n", CRUDE_GRID_SIZE);
#endif
	priority_queue<GridPoint> pq, opq;
	memset(owned, 0xff, sizeof(owned));
	get_best_grid_size(w, h);
	float hsize = w / (float) gw;
	float vsize = h / (float) gh;
	int hdi = max(1, (int) (hsize/2));
	int vdi = max(1, (int) (vsize/2));
	
	vector<pair<int, int> > v;
	for (int y = 0; y <= gh; y++) {
		for (int x = 0; x <= gw; x++) {
			v.push_back(make_pair(x, y));
		}
	}
	//
	static bool failsafe_emmited = false;
	random_shuffle(v.begin(), v.end(), stl_rand);
	for (int i = 0; i < (int) v.size(); i++) {
		int x = v[i].first;
		int y = v[i].second;
		GridPoint gp = make_gridpoint(x, y);
		owned[gp.p[1] * w + gp.p[0]] = grid.size();
		grid.push_back(gp);
		if (!failsafe_emmited && q.time() > 8.0) {
			failsafe_emmited = true;
			submit_adaptive_grid(SUBMIT_BRIEF);
		}
	}
	q.flush_cache();
	for (int i = 0; i < (int) grid.size(); i++) {
		GridPoint& gp = grid[i];
		gvalues[gp.g[1] * (gw + 1) + gp.g[0]] = q.data[gp.p[1] * w + gp.p[0]];
		for (int j = 0; j < 2; j++)
			gpoints[gp.g[1] * (gw + 1) + gp.g[0]][j] = gp.p[j];
	}
	float sum_gradients = 0.0;
	for (int i = 0; i < (int) grid.size(); i++) {
		calculate_gradient(grid[i]);
		pq.push(grid[i]);
		sum_gradients += grid[i].gradient;
	}
#ifdef PRINT_DEBUGINFO
	printf("Noise: %.4f\n", sum_gradients / grid.size());
#endif
	sortedgrid = grid;
	sort(sortedgrid.begin(), sortedgrid.end(), cmp_gp);
	/*
	for (int y = 0; y <= gh; y++) {
		for (int x = 0; x <= gw; x++)
			printf("%3d ", gvalues[y * (gw + 1) + x]);
		printf("\n");
	}
	printf("Gradients:\n");
	for (int y = 0; y <= gh; y++) {
		for (int x = 0; x <= gw; x++)
			printf("%6.3f ", sortedgrid[y * (gw + 1) + x].gradient);
		printf("\n");
	}
	*/
	submit_adaptive_grid(SUBMIT_BRIEF);
	int num_iters = 0;
	bool switched_to_full = false;
	while (!pq.empty() && q.time() < 9.4 /* FIXME */) {
		GridPoint gp = pq.top();
		pq.pop();
		if (pq.empty() && !switched_to_full) {
			switched_to_full = true;
			pq = opq;
		}
		if (gp.area < 0.25f) continue;
		/*
		const int dir[4][2] = { { -1, 0 }, { 0, 1 }, { 1, 0 }, { 0, -1 } };
		bool good = true;
		for (int k = 0; k < 4; k++) {
			int ngx = gp.g[0] + dir[k][0];
			int ngy = gp.g[1] + dir[k][1];
			if (ngx >= 0 && ngx <= gw && ngy >= 0 && ngy <= gh && taboo[ngy * (gw + 1) + ngx])
				good = false;
		}
		if (!good) continue;
		taboo[gp.g[1] * (gw + 1) + gp.g[0]] = true;
		*/
		++num_iters;
		q.mark_selected(gp.p[0], gp.p[1]);
		gp.area /= 2;
		vector<GridPoint> t;
		t.push_back(gp);
		for (int dy = -1; dy <= 1; dy++) {
			for (int dx = -1; dx <= 1; dx++) {
				if (dx == 0 && dy == 0) continue;
				GridPoint np;
				float cx, cy;
				if (gp.g[0] != -1) {
					if (verybig) {
						cx = (w - 1) * (gp.g[0] + dx / 2.0f) / (float) gw;
						cy = (h - 1) * (gp.g[1] + dy / 2.0f) / (float) gh;
					} else {
						cx = gp.p[0] + hdi * dx;
						cy = gp.p[1] + vdi * dy;
					}
				} else {
					cx = gp.p[0] + dx * vsize * gp.area;
					cy = gp.p[1] + dy * hsize * gp.area;
				}
				if (make_subgridpoint(cx, cy, gp.area, np)) {
					t.push_back(np);
				}
			}
		}
		
		//q.flush_cache();
		
		for (int i = 0; i < (int) t.size(); i++) {
			recalculate_gradient(t[i]);
			if (switched_to_full)
				pq.push(t[i]);
			else
				opq.push(t[i]);
		}
		
		if (q.time() > 9.0 && num_iters % 4 == 0) {
			submit_adaptive_grid(SUBMIT_BRIEF);
		}
	}
	q.flush_cache();
	xfinal = true;
	submit_adaptive_grid(SUBMIT_FINAL);
}

int main(array<System::String ^> ^args)
{
	q.init(w, h);
	int size = w * h;
	verybig = (size > 128*95);
	//verybig = true;
	if (size < BRUTEFORCE_THRESH) {
		solve_bruteforce();
	} else {
		solve_adaptive_grid();
	}
	return 0;
}
