#define _DEFAULT_SOURCE // for M_PI
#include <stdio.h>
#include <stdint.h>
#include <string.h>
#include <assert.h>
#include <stdbool.h>
#include <stdlib.h>
#include <math.h>

#define UNUSED(X) (void)(X)

// Defining data structures
typedef struct Vector {
  float data[3];
} Vector;

typedef struct Matrix {
  Vector axis[3];
} Matrix;

enum e_attribute_in {
  AIN_POSITION,
  AIN_COLOR,
  AIN_COUNT,
};

enum e_attribute_out {
  AOUT_POSITION,
  AOUT_NORMAL,
  AOUT_COLOR,
  AOUT_COUNT,
};

typedef struct attribute {
  const Vector* vertex;
  const unsigned (*index)[3];
} Attribute;

typedef struct Geometry {
  Attribute attribute[AIN_COUNT];
  unsigned triangle_count;
} Geometry;

typedef struct Triangle {
  Vector vertex[3];
} Triangle;

typedef struct PolySlice {
  double y, x[2];
  Vector baryzentric[2];
} PolySlice;

typedef struct Uniform {
  Matrix modelview;
  Vector light;
} Uniform;

// Constants

const Matrix indentity_matrix = {{
  {{1,0,0}},
  {{0,1,0}},
  {{0,0,1}},
}};

// Specifying geometry to render
const Geometry box = {
  .attribute = {
    [AIN_POSITION] = {
      .vertex = (const Vector[]){
        {{-1,-1,-1}}, {{ 1,-1,-1}}, {{-1, 1,-1}}, {{ 1, 1,-1}},
        {{-1,-1, 1}}, {{ 1,-1, 1}}, {{-1, 1, 1}}, {{ 1, 1, 1}},
      },
      .index = (const unsigned[][3]){
        {1,0,2},  {1,2,3}, // Front
        {3,2,6},  {3,6,7}, // Bottom
        {7,6,4},  {7,4,5}, // Back
        {5,4,0},  {5,0,1}, // Top
        {0,2,6},  {0,6,4}, // Left
        {3,1,5},  {3,5,7}, // Right
      }
    },
    [AIN_COLOR] = {
      .vertex = (const Vector[]){
        {{1,0,0}}, {{0,1,0}}, {{0,0,1}},
        {{1,1,0}}, {{0,1,1}}, {{1,0,1}},
      },
/*      .index = (const unsigned[][3]){
        {0,0,0},  {0,0,0}, // Front
        {1,1,1},  {1,1,1}, // Bottom
        {2,2,2},  {2,2,2}, // Back
        {3,3,3},  {3,3,3}, // Top
        {4,4,4},  {4,4,4}, // Left
        {5,5,5},  {5,5,5}, // Right
      }*/
      .index = (const unsigned[][3]){
        {3,3,3},  {3,3,3}, // Front
        {3,3,3},  {3,3,3}, // Bottom
        {3,3,3},  {3,3,3}, // Back
        {3,3,3},  {3,3,3}, // Top
        {3,3,3},  {3,3,3}, // Left
        {3,3,3},  {3,3,3}, // Right
      }
    },
  },
  .triangle_count = 12
};

// Utility functions
Vector vmulf (Vector v, double f){ return (Vector){{  v.data[0] * f         , v.data[1] * f         , v.data[2] * f          }}; }
Vector vmul  (Vector a, Vector b){ return (Vector){{  a.data[0] * b.data[0] , a.data[1] * b.data[1] , a.data[2] * b.data[2]  }}; }
Vector vdivf (Vector v, double f){ return (Vector){{  v.data[0] / f         , v.data[1] / f         , v.data[2] / f          }}; }
Vector vdiv  (Vector a, Vector b){ return (Vector){{  a.data[0] / b.data[0] , a.data[1] / b.data[1] , a.data[2] / b.data[2]  }}; }
Vector vaddf (Vector v, double f){ return (Vector){{  v.data[0] + f         , v.data[1] + f         , v.data[2] + f          }}; }
Vector vadd  (Vector a, Vector b){ return (Vector){{  a.data[0] + b.data[0] , a.data[1] + b.data[1] , a.data[2] + b.data[2]  }}; }
Vector vsubf (Vector v, double f){ return (Vector){{  v.data[0] - f         , v.data[1] - f         , v.data[2] - f          }}; }
Vector vsub  (Vector a, Vector b){ return (Vector){{  a.data[0] - b.data[0] , a.data[1] - b.data[1] , a.data[2] - b.data[2]  }}; }
double vdot  (Vector a, Vector b){ return             a.data[0] * b.data[0] + a.data[1] * b.data[1] + a.data[2] * b.data[2]; }
Vector vmaxf (Vector v, double f){ return (Vector){{  fmax(v.data[0],f)     , fmax(v.data[0],f)     , fmax(v.data[0],f)      }}; }
Vector vmaxv (Vector a, Vector b){ return (Vector){{  fmax(a.data[0],b.data[0]), fmax(a.data[1],b.data[1]), fmax(a.data[2],b.data[2]) }}; }
Vector vneg  (Vector v){ return (Vector){{  -v.data[0], -v.data[1], -v.data[2] }}; }
Vector vnormalize(Vector v){   return vdivf(v, sqrt(vdot(v,v)) ); }

Vector mmulv(Matrix m, Vector v){
  return (Vector){{
    m.axis[0].data[0]*v.data[0] + m.axis[1].data[0]*v.data[1] + m.axis[2].data[0]*v.data[2],
    m.axis[0].data[1]*v.data[0] + m.axis[1].data[1]*v.data[1] + m.axis[2].data[1]*v.data[2],
    m.axis[0].data[2]*v.data[0] + m.axis[1].data[2]*v.data[1] + m.axis[2].data[2]*v.data[2],
  }};
}

Vector vcross(Vector a, Vector b){
  return (Vector){{
    a.data[1]*b.data[2] - a.data[2]*b.data[1],
    a.data[2]*b.data[0] - a.data[0]*b.data[2],
    a.data[0]*b.data[1] - a.data[1]*b.data[0],
  }};
}

Vector vinterpolate(Vector a, Vector b, double t){
  return (Vector){{
    a.data[0]*(1.-t) + b.data[0]*t,
    a.data[1]*(1.-t) + b.data[1]*t,
    a.data[2]*(1.-t) + b.data[2]*t,
  }};
}

Vector bcoords_interpolate(Vector v[3], Vector bcoord){
  return (Vector){{
    v[0].data[0]*bcoord.data[0] + v[1].data[0]*bcoord.data[1] + v[2].data[0]*bcoord.data[2],
    v[0].data[1]*bcoord.data[0] + v[1].data[1]*bcoord.data[1] + v[2].data[1]*bcoord.data[2],
    v[0].data[2]*bcoord.data[0] + v[1].data[2]*bcoord.data[1] + v[2].data[2]*bcoord.data[2],
  }};
}

Matrix rotateX(double a){
  a = a / 180 * M_PI;
  return (Matrix){{
    {{ 1,      0,      0 }},
    {{ 0,  cos(a), sin(a) }},
    {{ 0, -sin(a), cos(a) }},
  }};
}

Matrix rotateY(double a){
  a = a / 180 * M_PI;
  return (Matrix){{
    {{  cos(a), 0, -sin(a) }},
    {{       0, 1,      0 }},
    {{  sin(a), 0, cos(a) }},
  }};
}

Matrix rotateZ(double a){
  a = a / 180 * M_PI;
  return (Matrix){{
    {{  cos(a), sin(a), 0}},
    {{ -sin(a), cos(a), 0}},
    {{       0,      0, 1}},
  }};
}

Matrix scale3d(Vector v){
  return (Matrix){{
    {{ v.data[0], 0, 0 }},
    {{ 0, v.data[1], 0 }},
    {{ 0, 0, v.data[2] }},
  }};
}

Matrix scale(double f){
  return (Matrix){{
    {{ f, 0, 0 }},
    {{ 0, f, 0 }},
    {{ 0, 0, f }},
  }};
}

Matrix mmulm(Matrix a, Matrix b){
  return (Matrix){{
    {{
      a.axis[0].data[0]*b.axis[0].data[0] + a.axis[1].data[0]*b.axis[0].data[1] + a.axis[2].data[0]*b.axis[0].data[2],
      a.axis[0].data[1]*b.axis[0].data[0] + a.axis[1].data[1]*b.axis[0].data[1] + a.axis[2].data[1]*b.axis[0].data[2],
      a.axis[0].data[2]*b.axis[0].data[0] + a.axis[1].data[2]*b.axis[0].data[1] + a.axis[2].data[2]*b.axis[0].data[2],
    }},{{
      a.axis[0].data[0]*b.axis[1].data[0] + a.axis[1].data[0]*b.axis[1].data[1] + a.axis[2].data[0]*b.axis[1].data[2],
      a.axis[0].data[1]*b.axis[1].data[0] + a.axis[1].data[1]*b.axis[1].data[1] + a.axis[2].data[1]*b.axis[1].data[2],
      a.axis[0].data[2]*b.axis[1].data[0] + a.axis[1].data[2]*b.axis[1].data[1] + a.axis[2].data[2]*b.axis[1].data[2],
    }},{{
      a.axis[0].data[0]*b.axis[2].data[0] + a.axis[1].data[0]*b.axis[2].data[1] + a.axis[2].data[0]*b.axis[2].data[2],
      a.axis[0].data[1]*b.axis[2].data[0] + a.axis[1].data[1]*b.axis[2].data[1] + a.axis[2].data[1]*b.axis[2].data[2],
      a.axis[0].data[2]*b.axis[2].data[0] + a.axis[1].data[2]*b.axis[2].data[1] + a.axis[2].data[2]*b.axis[2].data[2],
    }},
  }};
}

// Shaders
void triangle_shader(const Uniform*restrict uniform, Triangle out[restrict AOUT_COUNT], const Triangle in[restrict AIN_COUNT]){ // Not something found in regular pipelines, but useful for per-triangle stuff
  UNUSED(uniform);
  Vector normal = vnormalize(vcross(
    vsub(in[AIN_POSITION].vertex[1], in[AIN_POSITION].vertex[0]),
    vsub(in[AIN_POSITION].vertex[2], in[AIN_POSITION].vertex[0])
  ));
  out[AOUT_NORMAL].vertex[0] = normal;
  out[AOUT_NORMAL].vertex[1] = normal;
  out[AOUT_NORMAL].vertex[2] = normal;
}

void vertex_shader(const Uniform*restrict uniform, Vector out[AOUT_COUNT], const Vector in[AIN_COUNT]){
  UNUSED(uniform);
  out[AOUT_POSITION] = mmulv(uniform->modelview, in[AIN_POSITION]);
  out[AOUT_NORMAL] = mmulv(uniform->modelview, out[AOUT_NORMAL]);
  out[AOUT_COLOR] = in[AIN_COLOR];
}

Vector fragment_shader(const Uniform*restrict uniform, double*restrict depth, Vector varying[restrict AOUT_COUNT]){
  UNUSED(depth);
  float ambient_strength = 0.2;
  Vector normal = vnormalize(varying[AOUT_NORMAL]);
  Vector ambient_color = vmulf(varying[AOUT_COLOR], ambient_strength);
  Vector light_direction = vnormalize(vsub(uniform->light, varying[AOUT_POSITION]));
  Vector diffuse_color = vmulf(varying[AOUT_COLOR], fmax(vdot(normal, light_direction), 0.0));
  Vector color = vadd(ambient_color, diffuse_color);
  return color;
}

// Rendering
static bool PolySlice_cut(PolySlice slice[restrict 2], const PolySlice*restrict const next, const double y){
  const double dy = next->y-slice->y;
  if(dy <= 0 || y >= next->y || y <= slice->y)
    return false;
  const double t = (double)(y-slice->y)/dy;
  slice[1] = (PolySlice){
    .y = y,
    .x[0] = slice->x[0]*(1.-t) + next->x[0]*t,
    .x[1] = slice->x[1]*(1.-t) + next->x[1]*t,
    .baryzentric = {
      vinterpolate(slice->baryzentric[0], next->baryzentric[0], t),
      vinterpolate(slice->baryzentric[1], next->baryzentric[1], t),
    },
  };
  return true;
}

void draw_triangle(
  const uint32_t w,
  const uint32_t h,
  uint8_t image[restrict h][w][4],
  double depth_plane[restrict h][w],
  const Uniform*restrict uniform,
  Triangle triangle[AIN_COUNT]
){
  int si = 0;
  PolySlice slice[7];

  int a=0, b=1, c=2;
  { // Split triangle into vertical trapezoid slices, so we can just render what's inside the view port, and don't have to deal with coords outside -1..1 later.
    // FIXME: Currently only handles intersaction of top and bottom plane. It also needs to split at left / right planes, or it will really mess up the geometry & image in those cases.
    if(triangle[AOUT_POSITION].vertex[a].data[1] > triangle[AOUT_POSITION].vertex[b].data[1]){ a=1, b=0; }
    if(triangle[AOUT_POSITION].vertex[a].data[1] > triangle[AOUT_POSITION].vertex[c].data[1]){ c=a, a=2; }
    if(triangle[AOUT_POSITION].vertex[b].data[1] > triangle[AOUT_POSITION].vertex[c].data[1]){ int t=c; c=b, b=t; }
    if(triangle[AOUT_POSITION].vertex[a].data[1] >  1) return;
    if(triangle[AOUT_POSITION].vertex[c].data[1] < -1) return;

    slice[si++] = (PolySlice){
      .y    = triangle[AOUT_POSITION].vertex[a].data[1],
      .x[0] = triangle[AOUT_POSITION].vertex[a].data[0],
      .x[1] = triangle[AOUT_POSITION].vertex[a].data[0],
      .baryzentric = {
        {{1,0,0}},
        {{1,0,0}},
      },
    };
    {
      double bct = (triangle[AOUT_POSITION].vertex[b].data[1] - triangle[AOUT_POSITION].vertex[a].data[1]) / (triangle[AOUT_POSITION].vertex[c].data[1] - triangle[AOUT_POSITION].vertex[a].data[1]);
      double bx2 = bct * (triangle[AOUT_POSITION].vertex[c].data[0] - triangle[AOUT_POSITION].vertex[a].data[0]) + triangle[AOUT_POSITION].vertex[a].data[0];
      PolySlice bslice;
      bslice.y      = triangle[AOUT_POSITION].vertex[b].data[1];
      if(bx2 < triangle[AOUT_POSITION].vertex[b].data[0]){
        bslice.x[1] = triangle[AOUT_POSITION].vertex[b].data[0];
        bslice.x[0] = bx2;
        bslice.baryzentric[1] = (Vector){{0,1,0}};
        bslice.baryzentric[0] = (Vector){{1.-bct,0,bct}};
      }else{
        bslice.x[0] = triangle[AOUT_POSITION].vertex[b].data[0];
        bslice.x[1] = bx2;
        bslice.baryzentric[0] = (Vector){{0,1,0}};
        bslice.baryzentric[1] = (Vector){{1.-bct,0,bct}};
      }
      if(PolySlice_cut(&slice[si-1], &bslice, -1)) si++;
      if(PolySlice_cut(&slice[si-1], &bslice,  1)) si++;
      slice[si++] = bslice;
    }
    {
      PolySlice cslice = {
        .y    = triangle[AOUT_POSITION].vertex[c].data[1],
        .x[0] = triangle[AOUT_POSITION].vertex[c].data[0],
        .x[1] = triangle[AOUT_POSITION].vertex[c].data[0],
        .baryzentric = {
          {{0,0,1}},
          {{0,0,1}},
        },
      };
      if(PolySlice_cut(&slice[si-1], &cslice, -1)) si++;
      if(PolySlice_cut(&slice[si-1], &cslice,  1)) si++;
      slice[si++] = cslice;
    }
  }

  int first = -1, last=-1;
  for(int i=0; i<si; i++){
    const double sx=slice[i].x[0], ex=slice[i].x[1];
    if(ex < -1 || sx > 1 || slice[i].y < -1 || slice[i].y > 1 ){
      if(first != -1)
        break;
      continue;
    }
    if(first == -1)
      first = i;
    last = i;
    Vector sb = slice[i].baryzentric[0];
    Vector eb = slice[i].baryzentric[1];
    if(sx < -1){
      slice[i].x[0] = -1;
      slice[i].baryzentric[0] = vinterpolate(sb, eb, -(sx+1.)/2./(ex-sx));
    }
    if(ex >  1){
      slice[i].x[1] =  1;
      slice[i].baryzentric[1] = vinterpolate(sb, eb, (1.-(sx+1.)/2.)/(ex-sx));
    }
  }
  if(first == last)
    return;

  // Breseham would probably be faster, but this was simpler to figure out & I'm lazy
  for(int i=first; i<=last-1; i++){
    const PolySlice*const restrict s = &slice[i];
    const PolySlice*const restrict e = &slice[i+1];
    const uint32_t sy = (s->y+1.)/2. * (h-1);
    const uint32_t ey = (e->y+1.)/2. * (h-1);
    const uint32_t sxa[2] = { (s->x[0]+1.)/2.*(w-1), (s->x[1]+1.)/2.*(w-1) };
    const uint32_t exa[2] = { (e->x[0]+1.)/2.*(w-1), (e->x[1]+1.)/2.*(w-1) };
    const uint32_t ly = (ey - sy) ?  (ey - sy) : 1;
    for(uint32_t y=sy; y<=ey; y++){
      const uint32_t iy = h-y-1;
      // FIXME: In theory, I'd need 65bit in the absolute worst case
      // And don't use double here, integer arithmetic is used to avoid blank pixels due to non-linear precision errors
      const uint32_t sx = ( (uint64_t)sxa[0]*(ly-(y-sy)) + (uint64_t)exa[0]*(y-sy) )/ly;
      const uint32_t ex = ( (uint64_t)sxa[1]*(ly-(y-sy)) + (uint64_t)exa[1]*(y-sy) )/ly;
      const uint32_t lx = (ex - sx) ? (ex - sx) : 1;
      const double ty = ((double)y-sy)/ly;
      const Vector sb = vinterpolate(s->baryzentric[0], e->baryzentric[0], ty);
      const Vector eb = vinterpolate(s->baryzentric[1], e->baryzentric[1], ty);
      for(uint32_t x=sx; x<=ex; x++){
        const uint32_t ix = w-x-1;
        const double tx = ((double)x-sx)/lx;
        const Vector bcoord = vinterpolate(sb, eb, tx);
        Vector varying[AOUT_COUNT];
        for(unsigned i=0; i<AOUT_COUNT; i++)
          varying[i] = bcoords_interpolate((Vector[]){
            triangle[i].vertex[a],
            triangle[i].vertex[b],
            triangle[i].vertex[c],
          }, bcoord);
        Vector color = {0};
        double depth = varying[AOUT_POSITION].data[2];
        color = fragment_shader(uniform, &depth, varying);
        if(depth != depth || depth > depth_plane[y][x])
          continue;
        depth_plane[y][x] = depth;
        //color.data[0] = -depth*0.5+0.5;
        color = vmulf(color, 0x100);
        if(color.data[0] <= 0x00) color.data[0] = 0x00;
        if(color.data[1] <= 0x00) color.data[1] = 0x00;
        if(color.data[2] <= 0x00) color.data[2] = 0x00;
        if(color.data[0] >= 0xFF) color.data[0] = 0xFF;
        if(color.data[1] >= 0xFF) color.data[1] = 0xFF;
        if(color.data[2] >= 0xFF) color.data[2] = 0xFF;
        // ix / iy are flipped versions of x / y.
        image[iy][ix][2] = color.data[0];
        image[iy][ix][1] = color.data[1];
        image[iy][ix][0] = color.data[2];
        image[iy][ix][3] = 0xFF;
      }
    }
  }
}

// Draw the geometry
void draw(
  const uint32_t w,
  const uint32_t h,
  uint8_t image[h][w][4],
  double depth[h][w],
  const Uniform*restrict uniform,
  const Geometry*const restrict geometry
){
  for(unsigned i=0; i<geometry->triangle_count; i++){
    Triangle triangle_in[AIN_COUNT] = {0};
    for(enum e_attribute_in j=0; j<AIN_COUNT; j++){
      const Attribute attribute = geometry->attribute[j];
      const unsigned*restrict indeces = attribute.index[i];
      for(unsigned k=0; k<3; k++)
        triangle_in[j].vertex[k] = attribute.vertex[indeces[k]];
    }
    Triangle triangle_out[AOUT_COUNT] = {0};
    triangle_shader(uniform, triangle_out, triangle_in);
    for(unsigned k=0; k<3; k++){
      Vector input[AIN_COUNT];
      for(enum e_attribute_in j=0; j<AIN_COUNT; j++)
        input[j] = triangle_in[j].vertex[k];
      Vector output[AOUT_COUNT];
      for(enum e_attribute_out j=0; j<AOUT_COUNT; j++)
        output[j] = triangle_out[j].vertex[k];
      vertex_shader(uniform, output, input);
      for(enum e_attribute_out j=0; j<AOUT_COUNT; j++)
        triangle_out[j].vertex[k] = output[j];
    }
    draw_triangle(w,h,image,depth, uniform, triangle_out);
  }
}

// program logic
bool save_bmp(
  const char*restrict file,
  const uint32_t w,
  const uint32_t h,
  uint8_t image[h][w][4]
){
  FILE* nf = 0;
  FILE* of = 0;
  if(!strcmp(file, "-")){
    of = stdout;
  }else{
    of = nf = fopen(file, "wb");
    if(!nf) return false;
  }

  const size_t ims = sizeof(uint8_t[h][w][4]);
  const uint8_t header[54] = {
    'B','M', (sizeof(header)+ims),(sizeof(header)+ims)>>8,(sizeof(header)+ims)>>16,(sizeof(header)+ims)>>24, 0,0,0,0, sizeof(header),sizeof(header)>>8,sizeof(header)>>16,sizeof(header)>>24,
    40,0,0,0, w,w>>8,w>>16,w>>24, h,h>>8,h>>16,h>>24, 1,0, 32,0, 0,0,0,0, ims,ims>>8,ims>>16,ims>>24, 0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0
  };
  fwrite(header, 1, sizeof(header), of);
  fwrite(image , 1, ims, of);

  if(nf) fclose(nf);
  return true;
}

struct params {
  const char* file;
  uint32_t w;
  uint32_t h;
  int ry, rx;
};

struct params parse_args(int argc, char* argv[]){
  struct params p = {
    .w = 800,
    .h = 600,
    .ry = -20,
    .rx =  25
  };
  for(int i=1; i<argc; i++){
    if(argv[i][0] == '-' && argv[i][1] != '\0'){
      if(argv[i][2] != '\0' || i+1 >= argc)
        goto usage;
      switch(argv[i][1]){
        case 'w': p.w = atoi(argv[++i]); break;
        case 'h': p.h = atoi(argv[++i]); break;
        case 'y': p.ry = atoi(argv[++i]); break;
        case 'x': p.rx = atoi(argv[++i]); break;
        default: goto usage;
      }
    }else{
      if(i+1 != argc)
        goto usage;
      p.file = argv[i];
    }
  }
  if(!p.file)
    goto usage;
  return p;
usage:
  fprintf(stderr, "usage: %s [-w w|-h h|-y ry|-x rx] file.bmp\n", *argv);
  exit(1);
}

int main(int argc, char* argv[]){
  int ret = 0;
  const struct params p = parse_args(argc, argv);

  // Where do we place the light?
  Vector light = {{0.5,-1,1}};

  // We rotate the world
  Matrix m_view = indentity_matrix;
  m_view = mmulm(rotateY(p.ry), m_view);
  m_view = mmulm(rotateX(p.rx), m_view);

  // Initialise screen buffer
  uint8_t (*image)[p.w][4] = calloc(1,sizeof(uint8_t[p.h][p.w][4]));
  if(!image)
    return 1;
  double (*depth)[p.w] = malloc(sizeof(double[p.h][p.w]));
  if(!depth){
    free(image);
    return 1;
  }
  for(uint32_t y=0; y<p.h; y++)
    for(uint32_t x=0; x<p.w; x++)
      depth[y][x] = INFINITY;

  // Draw image
  {
    Matrix m_model = scale(0.5); // We scale down our cube
    draw(p.w,p.h,image,depth, &(Uniform){
      .modelview = mmulm(m_view, m_model),
      .light = mmulv(m_view, light),
    }, &box);
  }

  if(!save_bmp(p.file, p.w,p.h,image))
    ret = 1;

  free(depth);
  free(image);
  return ret;
}