/* ISC License */

#include <endian.h>
#include <errno.h>
#include <float.h>
#include <getopt.h>
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

#define die(code, msg) do { fprintf(stderr, "pineapple: %s\n", msg); exit(code); } while (0)
#define diesys(code, msg) do { fprintf(stderr, "pineapple: %s: %s\n", msg, strerror(errno)); exit(code); } while (0)

#define USAGE "pineapple [-v] [-w width] [-h height] [-f | -a]"
#define dieusage() die(100, "usage: " USAGE)

#define ASCII_CHARS "$@B%8&WM#*oahkbdpqwmZO0QLCJUYXzcvunxrjft/\\|()1{}[]?-_+~<>i!lI;:,\"^`'. "

static int verbosity = 0;

typedef struct { float x, y, z; } vec3;

struct triangle {
	vec3 n;
	vec3 a, b, c;
};

struct __attribute__((__packed__)) triangle_raw {
	vec3 n;
	vec3 a, b, c;
	uint16_t _attr_count; /* unused */
};

struct camera {
	vec3 pos;
	float theta;
	float phi;
	float fov;
};

static inline vec3
vadd(vec3 a, vec3 b)
{
	return (vec3){ a.x + b.x, a.y + b.y, a.z + b.z };
}

static inline vec3
vsub(vec3 a, vec3 b)
{
	return (vec3){ a.x - b.x, a.y - b.y, a.z - b.z };
}

static inline vec3
smul(vec3 v, float s)
{
	return (vec3){ v.x * s, v.y * s, v.z * s };
}

static inline float
dot(vec3 a, vec3 b)
{
	return a.x * b.x + a.y * b.y + a.z * b.z;
}

static inline vec3
cross(vec3 a, vec3 b)
{
	return (vec3){ a.y*b.z - b.y*a.z, a.z*b.x - b.z*a.x, a.x*b.y - b.x*a.y };
}

static inline float
magnitude(vec3 v)
{
	return sqrt(dot(v, v));
}

static inline vec3
norm(vec3 v)
{
	float mag = magnitude(v);
	return (vec3){ v.x / mag, v.y / mag, v.z / mag };
}

static float
ray_triangle_intersect(vec3 P, vec3 d, struct triangle tri)
{
	vec3 edge1 = vsub(tri.b, tri.a);
	vec3 edge2 = vsub(tri.c, tri.a);
	vec3 ray_cross_e2 = cross(d, edge2);
	float det = dot(edge1, ray_cross_e2);
	if (-FLT_EPSILON < det && det < FLT_EPSILON) return NAN;
	float inv_det = 1.0 / det;
	vec3 s = vsub(P, tri.a);
	float u = inv_det * dot(s, ray_cross_e2);
	if (u < 0.0 || u > 1.0) return NAN;
	vec3 s_cross_e1 = cross(s, edge1);
	float v = inv_det * dot(d, s_cross_e1);
	if (v < 0.0 || u + v > 1.0) return NAN;
	float t = -inv_det * dot(edge2, s_cross_e1);
	if (t > FLT_EPSILON) return t;
	return NAN;
}

static void
render_frame(struct camera cam, vec3 light,
             const struct triangle *triangles, uint32_t len,
             float *pixels, uint32_t width, uint32_t height)
{
	struct timespec start, now;
	if (verbosity >= 1) clock_gettime(CLOCK_MONOTONIC, &start);

	light = norm(light);
	for (size_t w = 0; w < width; ++w) {
		if (verbosity >= 1) {
			clock_gettime(CLOCK_MONOTONIC, &now);
			int elapsed = now.tv_sec - start.tv_sec;
			int total = elapsed * width / (w+1);
			fprintf(stderr, "%ld/%d %02d:%02d/%02d:%02d\r", w, width, elapsed / 60, elapsed % 60, total / 60, total % 60);
		}
		for (size_t h = 0; h < height; ++h) {
			float theta = cam.theta - cam.fov/2.0 + cam.fov/height*h;
			float phi = cam.phi - cam.fov/2.0 + cam.fov/width*w;
			vec3 dir = { sinf(theta) * cosf(phi), sinf(theta) * sinf(phi), cosf(theta) };
			dir = norm(dir);

			float shortest = INFINITY;
			struct triangle shortest_tri = {0};
			for (size_t i = 0; i < len; ++i) {
				float d = ray_triangle_intersect(cam.pos, dir, triangles[i]);
				if (d != NAN && d < shortest) {
					shortest = d;
					shortest_tri = triangles[i];
				}
			}

			pixels[w*height + h] = 0.0;
			if (shortest != INFINITY) {
				float light_angle = acosf(dot(shortest_tri.n, light));
				pixels[w*height + h] = light_angle / M_PI;
			}
		}
	}
	if (verbosity >= 1) {
		clock_gettime(CLOCK_MONOTONIC, &now);
		int elapsed = now.tv_sec - start.tv_sec;
		fprintf(stderr, "Rendered in %02d:%02d\n", elapsed / 60, elapsed % 60);
	}
}

int
main(int argc, char *argv[])
{
	uint32_t width = 80;
	uint32_t height = 40;
	enum { ASCII, FARBFELD } mode = ASCII;
	int c;
	while ((c = getopt(argc, argv, "afvh:w:")) > 0) {
		errno = 0;
		switch (c) {
			case 'a': mode = ASCII; break;
			case 'f': mode = FARBFELD; break;
			case 'v': verbosity = 1; break;
			case 'h': height = strtoul(optarg, NULL, 10); if (errno) dieusage(); break;
			case 'w': width = strtoul(optarg, NULL, 10); if (errno) dieusage(); break;
			default: dieusage();
		}
	}
	argv += optind; argc -= optind;
	if (argc > 0) dieusage();

	char header[80];
	if (fread(header, 1, 80, stdin) < 80)
		diesys(111, "fatal: unable to read STL header from stdin");
	uint32_t len;
	if (fread(&len, 4, 1, stdin) < 1)
		diesys(111, "fatal: unable to read STL length from stdin");
	len = le32toh(len);
	struct triangle *triangles = calloc(len, sizeof(struct triangle));
	if (!triangles) diesys(111, "fatal: unable to allocate enough memory for STL facets");
	struct triangle_raw raw;
	for (size_t i = 0; i < len; ++i) {
		if (fread(&raw, sizeof(raw), 1, stdin) < 1)
			diesys(111, "fatal: unable to read STL facets");
		triangles[i] = (struct triangle){ .n = norm(raw.n), .a = raw.a, .b = raw.b, .c = raw.c };
	}
	if (verbosity >= 1) fprintf(stderr, "finished reading STL: %d facets\n", len);

	struct camera cam = {
		.pos = { -70.0, 60.0, 60.0 },
		.phi = M_PI * 0.8,
		.theta = M_PI * 0.5,
		.fov = 80.0 * M_PI / 180.0
	};
	vec3 light = { 0.0, 0.5, -1.0 };

	if (mode == FARBFELD) {
		uint32_t wbe = htobe32(width);
		uint32_t hbe = htobe32(height);
		fwrite("farbfeld", 1, 8, stdout);
		fwrite(&wbe, 1, 4, stdout);
		fwrite(&hbe, 1, 4, stdout);
	}

	float *pixels = calloc(sizeof(float), width * height);
	if (!pixels) diesys(111, "fatal: unable to allocate enough memory for pixel buffer");
	render_frame(cam, light, triangles, len, pixels, width, height);

	for (int h = height - 1; h >= 0; --h) {
		for (uint32_t w = 0; w < width; ++w) {
			float brightness = pixels[w*height + h];
			if (mode == ASCII) {
				printf("%c", ASCII_CHARS[strlen(ASCII_CHARS) - 1 - (int)(brightness * strlen(ASCII_CHARS))]);
			} else if (mode == FARBFELD) {
				uint16_t x = htobe16((uint16_t)(brightness * 0xffff));
				uint16_t pixel[4] = { x, x, x, htobe16(0xffff) };
				fwrite(&pixel, 2, 4, stdout);
			}
		}
		if (mode == ASCII) printf("\n");
	}

	return 0;
}