#include <cassert>
#include <cmath>
#include <cerrno>
#include <cstdio>
#include <cstring>
#include <random>
#include <stdexcept>
#include <tuple>
#include <vector>

extern "C" {
#include <png.h>
}

template<typename T>
T clamp(const T &x, const T &min, const T &max)
{
	if (x < min)
		return min;
	if (x > max)
		return max;

	return x;
}

struct vector {
	double x, y;

	explicit vector(double x, double y):
		x(x), y(y)
	{
	}

	double normsq() const
	{
		return x * x + y * y;
	}

	double norm() const
	{
		return sqrt(normsq());
	}

	double angle() const
	{
		return atan2(y, x);
	}

	void operator+=(const vector &other)
	{
		x += other.x;
		y += other.y;
	}
};

bool operator==(const vector &a, const vector &b)
{
	return a.x == b.x && a.y == b.y;
}

vector operator-(const vector &a, const vector &b)
{
	return vector(a.x - b.x, a.y - b.y);
}

static double bezier_interpolate(double p0, double p1, double t)
{
	return (1 - t) * p0 + t * p1;
}

static double bezier_interpolate(double p0, double p1, double p2, double t)
{
	double tp = (1 - t);
	return tp * tp * p0 + 2 * tp * t * p1 + t * t * p2;
}

static double bezier_interpolate(double p0, double p1, double p2, double p3, double t)
{
	double tp = 1 - t;
	return tp * tp * tp * p0
		+ 3 * tp * tp * t * p1
		+ 3 * tp * t * t * p2
		+ t * t * t * p3;
}

static const unsigned int width = 320;
static const unsigned int height = 240;

static std::tuple<double, double, double> hsl_to_rgb(double H, double S, double L)
{
	double C = (1 - fabs(2 * L - 1)) * S;
	double X = C * (1 - fabs(fmod(H, 2) - 1));

	double r, g, b;
	if (H >= 0 && H < 1) {
		r = C; g = X; b = 0;
	} else if (H >= 1 && H < 2) {
		r = X; g = C; b = 0;
	} else if (H >= 2 && H < 3) {
		r = 0; g = C; b = X;
	} else if (H >= 3 && H < 4) {
		r = 0; g = X; b = C;
	} else if (H >= 4 && H < 5) {
		r = X; g = 0; b = C;
	} else if (H >= 5 && H <= 6) {
		r = C; g = 0; b = X;
	}

	double m = L - C / 2;

	return std::make_tuple(r + m, g + m, b + m);
}

template<unsigned int nr_corners>
struct blob {
	unsigned int x;
	unsigned int y;
	double size;

	double corners[nr_corners];
	double between[nr_corners];
};

template<typename Random, unsigned int nr_corners>
static void output_frame(Random &rand, const char *filename, const blob<nr_corners> &blue)
{
	png_byte data[width * height * 4];
	for (unsigned int y = 0; y < height; ++y) {
		for (unsigned int x = 0; x < width; ++x) {
			png_bytep rgb = &data[4 * (y * width + x)];
			rgb[0] = 0;
			rgb[1] = 0;
			rgb[2] = 0;
			rgb[3] = 255;
		}
	}

#define LESS

	for (unsigned int y = 0; y < height; ++y) {
		for (unsigned int x = 0; x < width; ++x) {
#ifdef LESS
			if (std::uniform_int_distribution<>(0, 100)(rand) < 50)
				continue;
#endif

			double H = 2 * std::uniform_int_distribution<>(0, 2)(rand);
			double S = std::uniform_real_distribution<>(0, 1)(rand);
#ifdef LESS
			double L = clamp(fabs(std::normal_distribution<>(0, 0.05)(rand)), 0., 1.);
#else
			double L = clamp(0.1 * std::exponential_distribution<>()(rand), 0., 1.);
#endif

			double r, g, b;
			std::tie(r, g, b) = hsl_to_rgb(H, S, L);

			png_bytep rgb = &data[4 * (y * width + x)];
			rgb[0] = 255 * r;
			rgb[1] = 255 * g;
			rgb[2] = 255 * b;
			rgb[3] = 255;
		}
	}

	unsigned int cx = blue.x;
	unsigned int cy = blue.y;
	unsigned int size = blue.size;

	for (unsigned int y = 0; y < 2 * size; ++y) {
		if (cy + y < size)
			continue;
		if (cy + y >= height + size)
			continue;

		for (unsigned int x = 0; x < 2 * size; ++x) {
			if (cx + x < size)
				continue;
			if (cx + x >= width + size)
				continue;

			vector p((double) y - size, (double) x - size);

			double corner = fmod(nr_corners * (M_PI + p.angle()) / 2 / M_PI, nr_corners);
			unsigned int corner_i = floor(corner);

			double t = corner - corner_i;
			assert(t >= 0 && t < 1);

			double norm = size * bezier_interpolate(
				blue.corners[corner_i],
				blue.corners[corner_i] + blue.between[corner_i],
				blue.corners[(corner_i + 1) % nr_corners] - blue.between[(corner_i + 1) % nr_corners],
				blue.corners[(corner_i + 1) % nr_corners],
				t);

			if (p.norm() > norm)
				continue;

			if (size > 100) {
				if (std::uniform_real_distribution<>(0, (p.norm() / norm))(rand) < (size - 100) / 500.)
					continue;
			}

			double H = fmod(std::normal_distribution<>(4.4, .15)(rand), 6);
			double S = clamp(1., 0., 1.);
			double L;
			if (p.norm() > 0.98 * norm)
				L = clamp(std::normal_distribution<>(.5, .1)(rand), 0., 1.);
			else
				L = clamp(std::normal_distribution<>(.3, .1)(rand), 0., 1.);

			double r, g, b;
			std::tie(r, g, b) = hsl_to_rgb(H, S, L);

			png_bytep rgb = &data[4 * ((cy + y - size) * width + (cx + x - size))];
			rgb[0] = 255 * r;
			rgb[1] = 255 * g;
			rgb[2] = 255 * b;
			rgb[3] = 255;
		}
	}

	FILE *fp = fopen(filename, "wb");
	if (!fp)
		throw std::runtime_error(strerror(errno));

	png_structp png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, 0, 0, 0);
	if (!png_ptr)
		throw std::runtime_error("png_create_write_struct");

	png_infop info_ptr = png_create_info_struct(png_ptr);
	if (!info_ptr)
		throw std::runtime_error("png_create_info_struct");

	png_init_io(png_ptr, fp);
	png_set_filter(png_ptr, 0, PNG_FILTER_NONE);
	png_set_IHDR(png_ptr, info_ptr, width, height, 8, PNG_COLOR_TYPE_RGBA,
		PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_DEFAULT,
		PNG_FILTER_TYPE_DEFAULT);

	png_bytep rows[height];
	for (unsigned int i = 0; i < height; ++i)
		rows[i] = &data[i * width * 4];

	png_set_rows(png_ptr, info_ptr, rows);
	png_write_png(png_ptr, info_ptr, PNG_TRANSFORM_IDENTITY, NULL);
	fclose(fp);
}

int main(int argc, char *argv[])
{
	std::random_device rd;
	std::mt19937 rand(rd());

	unsigned int frame_i = 0;

	for (unsigned int i = 0; i < 3; ++i) {
		static const unsigned int nr_corners = 6;
		blob<nr_corners> blue;
		blue.x = std::normal_distribution<>(width / 2, width / 8)(rand);
		blue.y = std::normal_distribution<>(height / 2, height / 8)(rand);
		blue.size = 0;

		for (unsigned int i = 0; i < nr_corners; ++i)
			blue.corners[i] = std::uniform_real_distribution<>(.6, 1)(rand);
		for (unsigned int i = 0; i < nr_corners; ++i)
			blue.between[i] = std::uniform_real_distribution<>(.1, .3)(rand);

		double size_inc = 1;
		for (unsigned int j = 0; blue.size < width; ++j) {
			printf("Frame %u (size %.2f)\n", frame_i, blue.size);

			char png_filename[100];
			sprintf(png_filename, "output/frame-%03u.png", frame_i);

			char gif_filename[100];
			sprintf(gif_filename, "output/frame-%03u.gif", frame_i);

			output_frame(rand, png_filename, blue);
			++frame_i;

			char command[100];
			sprintf(command, "convert %s %s && rm %s", png_filename, gif_filename, png_filename);
			system(command);

			size_inc *= 1.08;
			blue.size += size_inc;
		}
	}

	return 0;
}