#include <arpa/inet.h>
#include <fcntl.h>
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <sys/mman.h>
#include <sys/stat.h>
#include <sys/types.h>
#include <unistd.h>

// https://en.wikipedia.org/wiki/Hilbert_curve
// rotate/flip a quadrant appropriately
void rot(int n, int * x, int * y, int rx, int ry) {
	if(ry == 0) {
		if(rx == 1) {
			*x = n - 1 - *x;
			*y = n - 1 - *y;
		}

		// Swap x and y
		int t = *x;
		*x    = *y;
		*y    = t;
	}
}

// convert d to (x,y)
void d2xy(int n, int d, int * x, int * y) {
	int rx, ry, s, t = d;
	*x = *y = 0;
	for(s = 1; s < n; s *= 2) {
		rx = 1 & (t / 2);
		ry = 1 & (t ^ rx);
		rot(s, x, y, rx, ry);
		*x += s * rx;
		*y += s * ry;
		t /= 4;
	}
}

typedef float rgb_t[3];

int main(int argc, char ** argv) {
	int input = open(argv[1], O_RDONLY);
	if(input < 0) {
		printf("open(%s): %m\n", argv[1] ?: "NULL");
		return 1;
	}

	struct stat input_stat;
	if(fstat(input, &input_stat) < 0) {
		printf("fstat(input): %m\n");
		return 2;
	}

	off_t samples = input_stat.st_size / sizeof(float);
	uint32_t wh   = pow(2, ceil(log2(sqrt(samples))));
	printf("samples: %ld\n", samples);
	printf("size:    %dx\n", wh);

	float * track = mmap(NULL, input_stat.st_size, PROT_READ, MAP_SHARED, input, 0);
	if(track == MAP_FAILED) {
		printf("mmap(input): %m\n");
		return 5;
	}

	int output = open(argv[2], O_RDWR | O_CREAT | O_TRUNC, 0666);
	if(output < 0) {
		printf("open(%s): %m\n", argv[2] ?: "NULL");
		return 3;
	}

	size_t len = ((size_t)wh) * ((size_t)wh * 4 /* channels */ * sizeof(uint16_t)) + 8 + sizeof(uint32_t) * 2;
	// if(lseek(output, len - 1, SEEK_SET) == (off_t)-1) {
	if(ftruncate(output, len) == -1) {
		printf("lseek(output): %m\n");
		return 3;
	}

	printf("len:     %ld\n", len);
	uint8_t * image_raw = mmap(NULL, len, PROT_READ | PROT_WRITE, MAP_SHARED, output, 0);
	if(image_raw == MAP_FAILED) {
		printf("mmap(output): %m\n");
		return 4;
	}

	*(image_raw++) = 'f';
	*(image_raw++) = 'a';
	*(image_raw++) = 'r';
	*(image_raw++) = 'b';
	*(image_raw++) = 'f';
	*(image_raw++) = 'e';
	*(image_raw++) = 'l';
	*(image_raw++) = 'd';

	{
		uint32_t * image = (void *)image_raw;
		image[0] = image[1] = htonl(wh);
	}
	image_raw += 8;

	{
		uint16_t * image = (void *)image_raw;
		for(uint32_t x = 0; x < wh; ++x)
			for(uint32_t y = 0; y < wh; ++y)
				image[((y * wh) + x) * 4 /* channels */ + 3] = htons(0xFFFF);
	}

	int cmyk_chan = atoi(argv[3]);
	for(off_t i = 0; i < samples; ++i) {
		// nabijaczleweli@tarta:~/uwu/pix$ ls -lb cinque*
		// -rwxrw-r-- 1 nabijaczleweli users 409836400 Jul  6 01:43 cinque
		// -rwxrw-r-- 1 nabijaczleweli users  27552000 Jul 16 03:49 cinque2
		// -rwxrw-r-- 1 nabijaczleweli users 269213804 Jul 17 20:28 cinque3
		// -rwxrw-r-- 1 nabijaczleweli users  80406824 Jul 18 01:20 cinque4
		// -rw-r--r-- 1 nabijaczleweli users 787009028 Jul 18 01:50 cinque_all

		uint16_t * image = (void *)image_raw;
		int x, y;
		d2xy(wh, i, &x, &y);
		// printf("%ld: (%d, %d)\n", i, x, y);
		rgb_t rgb;
		if(cmyk_chan == 3) {
			rgb[0] = fabs(track[i]);
			rgb[1] = rgb[0] * (27552000.f / 409836400.f);
			rgb[2] = rgb[0] * (269213804.f / 409836400.f);
		} else
			rgb[cmyk_chan] = fabs(track[i]);
		image[((y * wh) + x) * 4 /* channels */ + 0] = htons(rgb[0] * 0xFFFF);
		image[((y * wh) + x) * 4 /* channels */ + 1] = htons(rgb[1] * 0xFFFF);
		image[((y * wh) + x) * 4 /* channels */ + 2] = htons(rgb[2] * 0xFFFF);
		image[((y * wh) + x) * 4 /* channels */ + 3] = htons(0xFFFF);
	}
}