/* 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>

#define PI 3.14159265358979323846

#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;:,\"^`'. "

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 */
};

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
vnorm(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 (det > -FLT_EPSILON && 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;
}

int
main(int argc, char *argv[])
{
	int verbosity = 0;
	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 = raw.n, .a = raw.a, .b = raw.b, .c = raw.c };
	}
	if (verbosity >= 1) fprintf(stderr, "finished reading STL: %d facets\n", len);

	vec3 camera_pos = { -70.0, 60.0, 60.0 };
	float camera_phi = PI * 0.8;
	float camera_theta = PI * 0.5;
	float fov = 80.0 * PI / 180.0;
	vec3 light = { 0.0, 0.0, -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);
	}

	for (int h = height - 1; h >= 0; --h) {
		if (verbosity >= 1) fprintf(stderr, "row %d/%d\r", height - h - 1, height);
		for (uint32_t w = 0; w < width; ++w) {
			float theta = camera_theta - fov/2.0 + fov/height*h;
			float phi = camera_phi - fov/2.0 + fov/width*w;
			vec3 dir = { sinf(theta) * cosf(phi), sinf(theta) * sinf(phi), cosf(theta) };
			dir = vnorm(dir);

			float shortest = INFINITY;
			vec3 shortest_norm = {0};
			for (size_t i = 0; i < len; ++i) {
				float d = ray_triangle_intersect(camera_pos, dir, triangles[i]);
				if (d != NAN && fabs(d) < fabs(shortest)) {
					shortest = d;
					shortest_norm = triangles[i].n;
				}
			}

			float brightness = 0.0;
			if (shortest != INFINITY) {
				float light_angle = acosf(dot(shortest_norm, light) / magnitude(shortest_norm) / magnitude(light));
				brightness = light_angle / PI;
			}
			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;
}