#include <math.h>
#include <stdio.h>
#include <stdint.h>
#include <float.h>
#include <limits.h>
#include <stdlib.h> // card > pixar.ppm

struct Vec {
  union {
    struct {
      float x;
      float y;
      float z;
      float w;
    };
    float vec[4];
  };

  Vec(float v = 0) {
    for (int i = 0; i < 4; i++) {
      vec[i] = v;
    }
  }

  Vec(float a, float b, float c = 0) {
    x = a;
    y = b;
    z = c;
    w = 0;
  }

  Vec operator+(Vec r) {
    Vec result;
    for (int i = 0; i < 4; i++) {
      result.vec[i] = vec[i] + r.vec[i];
    }
    return result;
  }

  Vec operator*(Vec r) {
    Vec result;
    for (int i = 0; i < 4; i++) {
      result.vec[i] = vec[i] * r.vec[i];
    }
    return result;
  }

  float operator%(Vec r) { return x * r.x + y * r.y + z * r.z; }

  //inv square root
  Vec operator!() {
    return *this * (1 / sqrtf(*this % *this));
  }
};

float min(float l, float r) { return l < r ? l : r; }

thread_local int seed = 1;
float randomVal() {
    seed = ( 214013 * seed + 2531011);
      return ((seed >> 16) & 0xFFFF) / 66635.0f;
}


// Rectangle CSG equation. Returns minimum signed distance from
// space carved by
// lowerLeft vertex and opposite rectangle vertex upperRight.
float BoxTest(Vec position, Vec lowerLeft, Vec upperRight) {
  lowerLeft = position + lowerLeft * -1;
  upperRight = upperRight + position * -1;
  return -min(
      min(min(lowerLeft.x, upperRight.x), min(lowerLeft.y, upperRight.y)),
      min(lowerLeft.z, upperRight.z));
}

#define HIT_NONE 0
#define HIT_LETTER 1
#define HIT_WALL 2
#define HIT_SUN 3

// Sample the world using Signed Distance Fields.
float QueryDatabase(Vec position, int &hitType) {
  float distance = 1e9;
  Vec f = position; // Flattened position (z=0)
  f.z = 0;
  char letters[15 * 4 + 1] = // 15 two points lines
      "5O5_"
      "5W9W"
      "5_9_" // P (without curve)
      "AOEO"
      "COC_"
      "A_E_" // I
      "IOQ_"
      "I_QO" // X
      "UOY_"
      "Y_]O"
      "WW[W" // A
      "aOa_"
      "aWeW"
      "a_e_"
      "cWiO"; // R (without curve)

  for (int i = 0; i < sizeof(letters); i += 4) {
    Vec begin = Vec(letters[i] - 79, letters[i + 1] - 79) * .5;
    Vec e = Vec(letters[i + 2] - 79, letters[i + 3] - 79) * .5 + begin * -1;
    Vec o =
        f + (begin + e * min(-min((begin + f * -1) % e / (e % e), 0), 1)) * -1;
    distance = min(distance, o % o); // compare squared distance.
  }
  distance = sqrtf(distance); // Get real distance, not square distance.

  // Two curves (for P and R in PixaR) with hard-coded locations.
  Vec curves[] = {Vec(-11, 6), Vec(11, 6)};
  for (int i = 2; i--;) {
    Vec o = f + curves[i] * -1;
    distance = min(distance, o.x > 0 ? fabsf(sqrtf(o % o) - 2)
                                     : (o.y += o.y > 0 ? -2 : 2, sqrtf(o % o)));
  }
  distance = powf(powf(distance, 8) + powf(position.z, 8), .125) - .5;
  hitType = HIT_LETTER;

  float roomDist;
  roomDist = min( // min(A,B) = Union with Constructive solid geometry
                  //-min carves an empty space
      -min(       // Lower room
          BoxTest(position, Vec(-30, -.5, -30), Vec(30, 18, 30)),
          // Upper room
          BoxTest(position, Vec(-25, 17, -25), Vec(25, 20, 25))),
      BoxTest( // Ceiling "planks" spaced 8 units apart.
          Vec(fmodf(fabsf(position.x), 8), position.y, position.z),
          Vec(1.5, 18.5, -25), Vec(6.5, 20, 25)));
  if (roomDist < distance)
    distance = roomDist, hitType = HIT_WALL;

  float sun = 19.9 - position.y; // Everything above 19.9 is light source.
  if (sun < distance)
    distance = sun, hitType = HIT_SUN;

  return distance;
}

// Perform signed sphere marching
// Returns hitType 0, 1, 2, or 3 and update hit position/normal
int RayMarching(Vec origin, Vec direction, Vec &hitPos, Vec &hitNorm) {
  int hitType = HIT_NONE;
  int noHitCount = 0;
  float d; // distance from closest object in world.

  // Signed distance marching
  for (float total_d = 0; total_d < 100; total_d += d)
    if ((d = QueryDatabase(hitPos = origin + direction * total_d, hitType)) <
            .01 ||
        ++noHitCount > 99)
      return hitNorm =
                 !Vec(QueryDatabase(hitPos + Vec(.01, 0), noHitCount) - d,
                      QueryDatabase(hitPos + Vec(0, .01), noHitCount) - d,
                      QueryDatabase(hitPos + Vec(0, 0, .01), noHitCount) - d),
             hitType; // Weird return statement where a variable is also
                      // updated.
  return 0;
}

Vec Trace(Vec origin, Vec direction) {
  Vec sampledPosition, normal, color, attenuation = 1;
  Vec lightDirection(!Vec(.6, .6, 1)); // Directional light

  for (int bounceCount = 3; bounceCount--;) {
    int hitType = RayMarching(origin, direction, sampledPosition, normal);
    if (hitType == HIT_NONE)
      break;                     // No hit. This is over, return color.
    if (hitType == HIT_LETTER) { // Specular bounce on a letter. No color acc.
      direction = direction + normal * (normal % direction * -2);
      origin = sampledPosition + direction * 0.1;
      attenuation = attenuation * 0.2; // Attenuation via distance traveled.
    }
    if (hitType == HIT_WALL) { // Wall hit uses color yellow?
      float incidence = normal % lightDirection;
      float p = 6.283185 * randomVal();
      float c = randomVal();
      float s = sqrtf(1 - c);
      float g = normal.z < 0 ? -1 : 1;
      float u = -1 / (g + normal.z);
      float v = normal.x * normal.y * u;
      direction =
          Vec(v, g + normal.y * normal.y * u, -normal.y) * (cosf(p) * s) +
          Vec(1 + g * normal.x * normal.x * u, g * v, -g * normal.x) *
              (sinf(p) * s) +
          normal * sqrtf(c);
      origin = sampledPosition + direction * .1;
      attenuation = attenuation * 0.2;
      if (incidence > 0 &&
          RayMarching(sampledPosition + normal * .1, lightDirection,
                      sampledPosition, normal) == HIT_SUN)
        color = color + attenuation * Vec(500, 400, 100) * incidence;
    }
    if (hitType == HIT_SUN) { //
      color = color + attenuation * Vec(50, 80, 100);
      break; // Sun Color
    }
  }
  return color;
}

const int w = 960;
const int h = 540;
int samplesCount = 1;
Vec position(-22, 5, 25);
Vec goal = !(Vec(-3, 4, 0) + position * -1);
Vec left = !Vec(goal.z, 0, -goal.x) * (1. / w);

// Cross-product to get the up vector
Vec up(goal.y *left.z - goal.z * left.y, goal.z *left.x - goal.x * left.z,
       goal.x *left.y - goal.y * left.x);

#include <pthread.h>
pthread_t threadIds[h];
Vec img[w * h];

void *plotPixel(void *arg) {
  int y = (int)(size_t)arg;
  seed = y;
  for (int x = w; x--;) {

    Vec color;
    for (int p = samplesCount; p--;)
      color =
          color + Trace(position, !(goal + left * (x - w / 2 + randomVal()) +
                                    up * (y - h / 2 + randomVal())));

    // Reinhard tone mapping
    color = color * (1. / samplesCount) + 14. / 241;
    Vec o = color + 1;
    color = Vec(color.x / o.x, color.y / o.y, color.z / o.z) * 255;

    img[y * w + x] = color;
  }
  return NULL;
}

int main(int argc, char **argv) {

  if (argc > 1) {
    samplesCount = atoi(argv[1]);
  }

  printf("P6 %d %d 255 ", w, h);
  for (int y = h; y--;) {
    pthread_create(&threadIds[y], NULL, plotPixel, (void *)y);
  }

  for (int y = 0; y < h; y++) {
    pthread_join(threadIds[y], NULL);
  }

  for (int y = h - 1; y >= 0; y--) {
    for (int x = w - 1; x >= 0; x--) {
      Vec *p;
      p = &img[w * y + x];
      printf("%c%c%c", (int)p->x, (int)p->y, (int)p->z);
    }
  }
  return 0;
}

// Andrew Kensler