// diydsp's threaded version of
// https://github.com/kid0m4n/gorays/blob/940c91f601ef840e6d75ddf272ab6cd3eb8d5531/ref/rays.cpp
// compile with g++ rt02.cpp -o rt02 -O3  and run with ./rt02
// set NUM_THREADS from 1-64, depending on your CPU
/*
As run on my AMD Dual-Core E300 CPU:

Compiled with -g
1 threads ~50 seconds
2 threads ~40 seconds
3 threads ~36 seconds
4 threads ~34 seconds
5 threads ~33 seconds

Compiled with -O
1 thread  ~22 seconds
2 threads ~18 seconds
3 threads ~18 seconds
4 threads ~17 seconds
5 threads ~19 secodns
6 threads ~16 seconds
7 threads ~17 seconds
8 threads ~17 seconds

Compiled with -O3
1 thread  ~17 seconds
2 threads ~12 seconds
3 threads ~12 seconds
4 threads ~11 seconds
5 threads ~12 seconds
6 threads ~12 seconds
*/

#include <stdlib.h>
#include <stdio.h>
#include <pthread.h>
#include <math.h>
#include <list>
#include <cstring>

#define NUM_THREADS 8
typedef struct {  int y_start,y_end,byte_index; } Params;
Params params[NUM_THREADS];

//Define a vector class with constructor and operator: 'v'
struct vector {
  float x,y,z; // Vector has three float attributes.
  vector operator+(vector r){return vector(x+r.x,y+r.y,z+r.z);} //Vector add
  vector operator*(float r){return vector(x*r,y*r,z*r);} //Vector scaling
  float operator%(vector r){return x*r.x+y*r.y+z*r.z;} //Vector dot product
  vector(){} //Empty constructor
  vector operator^(vector r){return vector(y*r.z-z*r.y,z*r.x-x*r.z,x*r.y-y*r.x);} //Cross-product
  vector(float a,float b,float c){x=a;y=b;z=c;} //Constructor
  vector operator!(){return *this*(1 /sqrt(*this%*this));} // Used later for normalizing the vector
};

const char *art[] = {
  " ",
  " 1111 ",
  " 1 1 ",
  " 1 11 ",
  " 1 1 1 ",
  " 1 11 1 1 ",
  " 1 1 1 1 ",
  " 1 1 1 1 ",
  " 11111 11 "
};

struct object {
  int k,j;
  object(int x,int y){k=x;j=y;}
};

std::list<object> objects;

void F() {
  int nr = sizeof(art) / sizeof(char *),
  nc = strlen(art[0]);
  for (int k = nc - 1; k >= 0; k--) {
    for (int j = nr - 1; j >= 0; j--) {
      if(art[j][nc - 1 - k] != ' ') {
        objects.push_back(object(k, nr - 1 - j));
      }
    }
  }
}

unsigned int seed = ~0;

float R() {
  seed += seed;
  seed ^= 1;
  if ((int)seed < 0)
    seed ^= 0x88888eef;
  return (float)(seed % 95) / (float)95;
}

//The intersection test for line [o,v].
// Return 2 if a hit was found (and also return distance t and bouncing ray n).
// Return 0 if no hit was found but ray goes upward
// Return 1 if no hit was found but ray goes downward
int T(vector o,vector d,float& t,vector& n) {
  t=1e9;
  int m=0;
  float p=-o.z/d.z;

  if(.01<p)
    t=p,n=vector(0,0,1),m=1;

  std::list<object>::iterator it;

  for (it = objects.begin(); it != objects.end(); ++it) {
    int k = it->k,
    j = it->j;

    // There is a sphere but does the ray hits it ?
    vector p=o+vector(-k,3,-j-4);
    float b=p%d,c=p%p-1,q=b*b-c;

    // Does the ray hit the sphere ?
    if(q>0) {
      //It does, compute the distance camera-sphere
      float s=-b-sqrt(q);

      if(s<t && s>.01)
      // So far this is the minimum distance, save it. And also
      // compute the bouncing ray vector into 'n'
      t=s, n=!(p+d*t), m=2;
    }
  }

  return m;
}

// (S)ample the world and return the pixel color for
// a ray passing by point o (Origin) and d (Direction)
vector S(vector o,vector d) {
  float t;
  vector n, on;

  //Search for an intersection ray Vs World.
  int m=T(o,d,t,n);
  on = n;

  if(!m) { // m==0
    //No sphere found and the ray goes upward: Generate a sky color
    float p = 1-d.z;
    p = p*p;
    p = p*p;
    return vector(.7,.6,1)*p;
  }

  //A sphere was maybe hit.

  vector h=o+d*t, // h = intersection coordinate
  l=!(vector(9+R(),9+R(),16)+h*-1); // 'l' = direction to light (with random delta for soft-shadows).

  //Calculated the lambertian factor
  float b=l%n;

  //Calculate illumination factor (lambertian coefficient > 0 or in shadow)?
  if(b<0||T(h,l,t,n))
    b=0;

  if(m&1) { //m == 1
    h=h*.2; //No sphere was hit and the ray was going downward: Generate a floor color
    return((int)(ceil(h.x)+ceil(h.y))&1?vector(3,1,1):vector(3,3,3))*(b*.2+.1);
  }

  vector r=d+on*(on%d*-2); // r = The half-vector

  // Calculate the color 'p' with diffuse and specular component
  float p=l%r*(b>0);
  float p33 = p*p;
  p33 = p33*p33;
  p33 = p33*p33;
  p33 = p33*p33;
  p33 = p33*p33;
  p33 = p33*p;
  p = p33*p33*p33;

  //m == 2 A sphere was hit. Cast an ray bouncing from the sphere surface.
  return vector(p,p,p)+S(h,r)*.5; //Attenuate color by 50% since it is bouncing (* .5)
}

  int w = 512,
  h = 512;

  // The '!' are for normalizing each vectors with ! operator.
  vector g=!vector(-5.5,-16,0), // Camera direction
    a=!(vector(0,0,1)^g)*.002, // Camera up vector...Seem Z is pointing up :/ WTF !
    b=!(g^a)*.002, // The right vector, obtained via traditional cross-product
    c=(a+b)*-256+g; // WTF ? See https://news.ycombinator.com/item?id=6425965 for more.

  int s = 3*w*h;
  char *bytes = new char[s];

// The main function. It generates a PPM image to stdout.
// Usage of the program is hence: ./card > erk.ppm
int main(int argc, char **argv) {
  F();


  if (argc > 1) {
    w = atoi(argv[1]);
  }
  if (argc > 2) {
    h = atoi(argv[2]);
  }

  printf("P6 %d %d 255 ", w, h); // The PPM Header is issued




  pthread_t thr[NUM_THREADS];
  //  void *render_thread(void *);
  void *render_thread(void *);

  for(int i=0;i<NUM_THREADS;i++)
  {
    params[i].y_start    = (i+1) * w / NUM_THREADS - 1;
    params[i].y_end      =  i    * w / NUM_THREADS    ;
    params[i].byte_index = params[i].y_end * w * 3;
    pthread_create(&thr[i],NULL,render_thread,&params[i]);
  }

  for(int i=0;i<NUM_THREADS;i++)
  {
    pthread_join(thr[i],NULL);
  }

  fwrite(bytes, 1, s, stdout);
  delete [] bytes;
}

void *render_thread(void *data)
{
  Params *p = (Params *)data;
  int k = p->byte_index;
  volatile int y;

  //printf("\ny_start   : %d ",p->y_start);
  //printf("y_end     : %d ", p->y_end);
  //printf("byte_index: %d\n", p->byte_index);

  for(y = p->y_start ; y >= p->y_end;  y--) //For each column
  {

    //printf(".%d\n",y);

    for(int x=w;x--;)  //For each pixel in a line
    {
      printf("x");
      //Reuse the vector class to store not XYZ but a RGB pixel color
      vector p(13,13,13); // Default pixel color is almost pitch black

      //	Cast 64 rays per pixel (For blur (stochastic sampling) and soft-shadows.
      for(int r=64;r--;)
      {
        // The delta to apply to the origin of the view (For Depth of View blur).
        vector t=a*(R()-.5)*99+b*(R()-.5)*99; // A little bit of delta up/down and left/right
        // Set the camera focal point vector(17,16,8) and Cast the ray
        // Accumulate the color returned in the p variable
        p=S(vector(17,16,8)+t, //Ray Origin
        !(t*-1+(a*(R()+x)+b*(y+R())+c)*16) // Ray Direction with random deltas
                                         // for stochastic sampling
        )*3.5+p; // +p for color accumulation
      }

      bytes[k++] = (char)p.x;
      bytes[k++] = (char)p.y;
      bytes[k++] = (char)p.z;
    }
    //printf(",%d\n",y);

  }
}