/* I apologise in advance for switching terminology several times */
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>

#ifndef M_PI
#define M_PI        3.14159265358979323846      /* pi */
#endif

/* caps to avoid any clash with C++ STL */
typedef struct Vector_struct {
  float	x;
  float y;
  float z;
} Vector;

typedef struct Ray_struct {
  Vector origin;
  Vector direction;
} Ray;

typedef struct Colour_struct {
  float r;
  float g;
  float b;
} Colour;

float intersect_sphere(const Ray *r, Vector *centre, float radius)
{
  /* Return lambda on hit, or zero on miss. */
  /* Because we're only tracing eye rays, */
  /* don't worry about the back surface. */
  /* Ray direction vector is normalised. */
  /* See p.388 Graphics Gems. */
  float v;
  Vector EO;
  float disc;
  float d;
  EO.x = centre->x - r->origin.x;
  EO.y = centre->y - r->origin.y;
  EO.z = centre->z - r->origin.z;
  v =
    EO.x*r->direction.x +
    EO.y*r->direction.y +
    EO.z*r->direction.z;
  disc = radius*radius - (EO.x*EO.x + EO.y*EO.y + EO.z*EO.z - v*v);
  if (disc < 0.0f)
    return 0.0f;
  d = (float) sqrt(disc);
  return v-d;
}

int clamp(int i, int min, int max)
{
  if (i < min) return min;
  if (i > max) return max;
  return i;
}

float fclamp(float f, float min, float max)
{
  if (f < min) return min;
  if (f > max) return max;
  return f;
}

float intersect_ray(const Ray *r)
{
  /* Return lambda on hit, or zero on miss. */
  /* The sphere grid is encoded here. */
  int xincrement; /* Direction encoding for x */
  int xpos, ypos; /* Ray enters sphere layer */
  int xmiss, ymiss; /* Ray leaves sphere layer */
  int xcount, ycount;
  float lambda;
  Vector sphere; /* Sphere centre */
  sphere.z = -0.5; /* Spheres are all 0.5 below the plane */
  /* Miss everything if we point vaguely up */
  /* (otherwise there's a division by zero that can happen) */
  if (r->direction.z >= -0.001)
  {
    return 0.0f;
  }

  /* Horrible business to walk through a subset of the spheres */

  xincrement = r->direction.x >= 0.0 ? 1 : -1;
  /* Find out where we hit the z=0 plane */
  ypos = floor(r->origin.y -
	       (r->origin.z/r->direction.z)*r->direction.y);
  xpos = floor(r->origin.x -
	       (r->origin.z/r->direction.z)*r->direction.x);
  /* Now define the rectangle of possible sphere intersections. */
  /* We don't try to be clever with Bressenham's algorithm or */
  /* anything (which would make this substantially faster) */
  ymiss = floor(r->origin.y -
		((r->origin.z+1.0f)/r->direction.z)*r->direction.y);
  xmiss = floor(r->origin.x -
		((r->origin.z+1.0f)/r->direction.z)*r->direction.x);

  ypos = clamp(ypos, 0, 8);
  ymiss = clamp(ymiss,-1, 7);
  xpos = clamp(xpos, -5, 4);
  xmiss = clamp(xmiss, -4, 3);
  /* fprintf(stderr, "xpos = %d, ypos = %d\n", xpos, ypos); */
  /* r.direction.y likely to be > r.direction.x, so iterate over */
  /* so iterate over y in the inner loop. */
  for (ycount = ypos; ycount <= ymiss; ycount++)
    for (xcount = xpos;
	 xcount * xincrement <= xmiss * xincrement;
	 xcount += xincrement)
    {
      /* Check the current cell's sphere */
      sphere.x = xcount + 0.5;
      sphere.y = ycount + 0.5;
      if (0.0f != (lambda = intersect_sphere(r, &sphere, 0.5f)))
	return lambda;
    }
  return 0.0f;
}

int zsort_helper(const void *p1, const void*p2)
{
  return *(*(float **)p1) > *(*(float **)p2);
}

void cause_confusion(Colour *buff,
		     /* Camera lens width in pixels on the focal plane */
		     const float focalwidth,
		     const float focaldepth,
		     float *zbuff,
		     unsigned int imgwidth,
		     unsigned int imgheight)
{
  /* Overwrite the image with one in which each pixel is mapped onto */
  /* a circle of confusion. */
  /* Note that for large lenses the proper way to do this is with ffts. */
  /* However, this way is easier to program. :-) */
  /* Start by sorting the scene by z */
  float **indirection_array;
  unsigned int xcounter, ycounter;
  unsigned long array_counter;
  Colour *second_buffer;
  float progress_counter = -1.0f;
  float pcinc = 80.0f/((float)imgwidth*imgheight);
  indirection_array = (float **) malloc(imgwidth*imgheight*sizeof(float *));
  /* The following assumes 0.0f == 0 in bitfields - which it is for IEEE */
  second_buffer = (Colour *) calloc(imgwidth*imgheight, sizeof(Colour));
  if (!indirection_array || !second_buffer)
  {
    fprintf(stderr, "Fatal error: insufficient space to sort image by z\n");
    exit(1);
  }
  for (ycounter=0U; ycounter<imgheight; ycounter++)
  {
    for (xcounter=0U; xcounter<imgwidth; xcounter++)
    {
      indirection_array[xcounter+imgwidth*ycounter] =
	zbuff + xcounter + imgwidth*ycounter;
    }
  }
  qsort((void *) indirection_array, imgwidth*imgheight, sizeof(float *),
	&zsort_helper);
  /* Start with the original image (?) */
  /* memcpy((void *)second_buffer, (void *)buff,
	 imgwidth*imgheight*sizeof(Colour)); */
  /* Iterate over the pixels, starting at the back. */
  for (array_counter = (long)imgwidth * imgheight;
       --array_counter && *(indirection_array[array_counter]) != 0.0f;)
  {
    int xcoord;
    int ycoord;
    int xcount, ycount;
    int xmax, ymax;
    float confusionradius; /* (okay, 1/area - a diffuseness measure) */
    float confusionarea;
    ycoord = (indirection_array[array_counter]-zbuff) / imgwidth;
    xcoord = (indirection_array[array_counter]-zbuff) % imgwidth;
    confusionradius = 1.0f + focalwidth * 0.5f *
      fabs(*indirection_array[array_counter] - focaldepth)/focaldepth;
    confusionarea = 1.0f/(M_PI*confusionradius*confusionradius);
    /* I'm forever drawing circles, pretty circles on the screen... */
    ymax = clamp(ycoord + confusionradius+1, 0, imgheight-1);

    for (ycount = clamp(ycoord - confusionradius-1, 0, imgheight-1);
	 ycount <= ymax; ycount++)
    {
      xmax = clamp(xcoord + sqrt((confusionradius+1)*(confusionradius+1)
				 - (ycount-ycoord)*(ycount-ycoord)),
		   0, imgwidth-1);
      for (xcount = clamp(xcoord
			  - sqrt((confusionradius+1)*(confusionradius+1)
				 -(ycount-ycoord)*(ycount-ycoord)),
			  0, imgwidth-1);
	   xcount <= xmax; xcount++)
      {
	float scalefactor;
	scalefactor = (float)
	  sqrt((xcount-xcoord)*(xcount-xcoord) +
	       (ycount-ycoord)*(ycount-ycoord));
	/* Start at radius, and handle blending */
	scalefactor = confusionradius - scalefactor + 1.0f;
	/* Antialias the edge */
	scalefactor = fclamp(scalefactor, 0.0f, 1.0f);
	scalefactor = sqrt(scalefactor);
	/* If we don't go through all this we get stepping effects */
	/* near the focal plane (where pixel jumps are larger) */
	/* Partly occlude what's behind */
	second_buffer[xcount + ycount*imgwidth].r *=
	  (1.0f-confusionarea*scalefactor);
	second_buffer[xcount + ycount*imgwidth].g *=
	  (1.0f-confusionarea*scalefactor);
	second_buffer[xcount + ycount*imgwidth].b *=
	  (1.0f-confusionarea*scalefactor);
	/* Add in the contribution from this pixel */
	second_buffer[xcount + ycount*imgwidth].r += confusionarea *
	  buff[xcoord + ycoord*imgwidth].r * scalefactor;
	second_buffer[xcount + ycount*imgwidth].g += confusionarea *
	  buff[xcoord + ycoord*imgwidth].g * scalefactor;
	second_buffer[xcount + ycount*imgwidth].b += confusionarea *
	  buff[xcoord + ycoord*imgwidth].b * scalefactor;
      }
    }
    progress_counter += pcinc;
    if (progress_counter>0.0f)
    {
      progress_counter -= 1.0f;
      fprintf(stderr, ".");
    }
  }
  fprintf(stderr, "\n");
  /* Done. Copy the second image over the first. */
  memcpy((void *)buff, (void *)second_buffer,
	 imgwidth*imgheight*sizeof(Colour));
  free(second_buffer);
  free(indirection_array);
}

void dump_ppm(FILE *file,
	      Colour *buff,
	      unsigned int width,
	      unsigned int height)
{
  /* Currently do this a very slow and stupid way. */
  /* Might want to consider filling up a buffer */
  /* and using fwrite() at some point. */
  unsigned int wcounter, col=1;
  int r,g,b;
  fprintf(file, "P3\n# CREATOR: renderspheres\n%u %u\n255\n",
	  width, height);
  while (height--)
    for (wcounter=0; wcounter<width; wcounter++)
    {
      r = (int) ((buff+height*width+wcounter)->r*255);
      r = (r < 0 ? 0 : r > 255 ? 255 : r);
      g = (int) ((buff+height*width+wcounter)->g*255);
      g = (g < 0 ? 0 : g > 255 ? 255 : g);
      b = (int) ((buff+height*width+wcounter)->b*255);
      b = (b < 0 ? 0 : b > 255 ? 255 : b);
      fprintf(file, "%c%c%c %c%c%c %c%c%c  %s",
	     '0'+r/100, '0'+(r%100)/10, '0'+r%10,
	     '0'+g/100, '0'+(g%100)/10, '0'+g%10,
	     '0'+b/100, '0'+(b%100)/10, '0'+b%10,
	     col++%5?"":"\n");
    }
}

void parseparams(int argc, char** argv,
		 unsigned int *width,
		 unsigned int *height,
		 float *focaldepth,
		 float *lenswidth,
		 float *hfov,
		 float *vfov,
		 unsigned int *numsamples,
		 int *dozbuffer,
		 char *filename)
{
  while (-- argc > 1)
  {
    argv++;
    if (0 == strcmp(*argv, "width"))
      (void) (sscanf(*++argv, "%u", width) &&
	      fprintf(stderr, "Set width\n"));
    else if (0 == strcmp(*argv, "height"))
      (void) (sscanf(*++argv, "%u", height) &&
	      fprintf(stderr, "Set height\n"));
    else if (0 == strcmp(*argv, "focaldepth"))
      (void) (sscanf(*++argv, "%f", focaldepth) &&
	      fprintf(stderr, "Set focaldepth\n"));
    else if (0 == strcmp(*argv, "lenswidth"))
      (void) (sscanf(*++argv, "%f", lenswidth) &&
	      fprintf(stderr, "Set lenswidth\n"));
    else if (0 == strcmp(*argv, "hfov"))
      (void) (sscanf(*++argv, "%f", hfov) &&
	      fprintf(stderr, "Set hfov\n"));
    else if (0 == strcmp(*argv, "vfov"))
      (void) (sscanf(*++argv, "%f", vfov) &&
	      fprintf(stderr, "Set vfov\n"));
    else if (0 == strcmp(*argv, "numsamples"))
      (void) (sscanf(*++argv, "%u", numsamples) &&
	      fprintf(stderr, "Set numsamples\n"));
    else if (0 == strcmp(*argv, "dozbuffer")) {
      *dozbuffer = 1;
      /* If we're doing the z buffer, we need to take single samples */
      *numsamples = 1;
      fprintf(stderr, "Doing z buffer\n");
      argv++;
    } else if (0 == strcmp(*argv, "filename"))
      (void) (strcpy(filename, *++argv) &&
	      fprintf(stderr, "Set filename to %s\n", filename));
    else 
      {
	fprintf(stderr, "Argument %s not recognised\n", *argv);
	continue;
      }
    argc--;
  }
}

int main(int argc, char ** argv)
{
  char filename[255]="spheres.ppm";
  unsigned int width = 320;
  unsigned int height = 240;
  float focaldepth = 2.0f; /* Focal plane in Z */
  float lenswidth = 0.5f; /* Width of camera lens (at y=0) */
  float hfov = 1.0f; /* Horizontal field of view */
  float vfov = 0.75f; /* Vertical field of view */
  FILE* output;
  Colour *screen;
  float *zbuffer = (float *)0;
  /* Our scene has spheres just below the z=0 plane. */
  /* Our projected screen on the focal plane touches the z=0 plane at */
  /* its topmost extent. */
  /* Our camera touches the z=0 plane at its lowest extent. */

  /* The focal plane is tilted so that it's coplanar with the camera, and */
  /* the centre of the two has a mutually perpendicular vector */
  /* of focaldepth length. */
  /* Otherwise we get lens distortion, and it looks awful. */

  /* Start looking for spheres at the box where the ray hits z=0. */
  /* Go l/r or r/l depending on ray x direction. */
  /* Black if we don't hit anything. */
  /* There are 64 spheres in an 8x8 grid. */
  float screenwidth;
  float screenheight;
  unsigned int numsamples = 16U; /* Number of (ray tracing) samples */
  unsigned int samplecounter;
  /* The screen is 8 units wide, at its focal depth */
  unsigned int imgx, imgy;
  Ray r;
  float lambda=0.0f;
  Vector focalPlaneCentre;
  Vector cameraCentre;
  Vector focalPlaneWidth;
  Vector focalPlaneHeight;
  Vector cameraWidth;
  Vector cameraHeight;
  float fdc; /* Component of focal distance above z=0 plane */

  float previewx=-.5f;
  float previewxinc;
  float previewy=-.5f;
  float previewyinc;

  int doxpreview = 0;
  int dozbuffer = 0;

  parseparams(argc, argv,
	      &width, &height,
	      &focaldepth, &lenswidth,
	      &hfov, &vfov,
	      &numsamples,
	      &dozbuffer,
	      &filename[0]);
  screen = (Colour *) malloc (sizeof(Colour) * width * height);
  if (!screen)
  {
    fprintf(stderr, "Fatal error: Insufficient memory for frame buffer\n");
    exit(1);
  }

  if (dozbuffer)
  {
    zbuffer = (float *) malloc (sizeof(float) * width * height);
    if (!zbuffer)
    {
      fprintf(stderr, "Fatal error: Insufficient memory for z buffer\n");
      exit(1);
    }
  }

  screenwidth = 2.0f*focaldepth*hfov;
  screenheight = 2.0f*focaldepth*vfov;
  previewxinc = (float)72.0f/width;
  previewyinc = (float)22.0f/height;

  /* Correct the vectors for looking downwards. */
  /* This is probably harder than putting the ray through */
  /* some appropriate matrix... */
  /* We're looking straight down the y axis. */
  focalPlaneCentre.x = 0.0f; /* Assumed */
  cameraCentre.x = 0.0f; /* Assumed */
  cameraCentre.y = 0.f; /* Defined */
  focalPlaneWidth.x = screenwidth; /* Assumed */
  focalPlaneWidth.y = focalPlaneWidth.z = 0.0f; /* Assumed */
  cameraWidth.x = lenswidth; /* Assumed */
  cameraWidth.y = cameraWidth.z = 0.0f; /* Assumed */
  /* Okay, now for the off-alignment bits. */
  cameraHeight.x = 0.0f;
  focalPlaneHeight.x = 0.0f;
  /* Draw a little diagram:            c
                                     /\ a
                                    /  \|m
				   /   .X e
                                  / _-^ |\ r
    _____________________________/.-____|_\_a_
    z=0           \           _-Y       |
               f   \       _-^ /        |
                o p \   _-^   /         |
                 c l \-^     /          |
                  a a \     /           |
                   l n \   /            |y=0
                      e \ /             |
                         V              |     */
  /* Screen height is the focal plane diagonal */
  /* Camera height is the camera plane diagonal */
  /* Focal depth is the length of the _-^ line */
  /* Camera is centred on y=0 */
  /* Screen and camera touch z=0 */
  /* Note the many similar triangles */

  /* Component of focal depth in "camera triangle" is */
  /* fdc = focal depth * (lens width / (lens width + screen height) */
  /* Similar triangles within the camera triangle give our starting point */
  /* Projection of camera triangle onto z=0 plane is: */
  /* sqrt(fdc^2 + (lenswidth/2)^2) */
  /* cameraHeight.y = lenswidth*lenswidth*.5/sqrt(fdc^2+(lenswidth/2)^2) */
  fdc = focaldepth*(lenswidth/(lenswidth+screenheight));
  cameraHeight.y = lenswidth *
    (lenswidth*0.5f)/(float)sqrt(fdc*fdc+lenswidth*lenswidth*0.25);
  /* Pythagoras */
  cameraHeight.z = (float) sqrt(lenswidth*lenswidth -
				cameraHeight.y*cameraHeight.y);
  /* Similar triangles... */
  focalPlaneHeight.y = screenheight/lenswidth * cameraHeight.y;
  focalPlaneHeight.z = screenheight/lenswidth * cameraHeight.z;
  /* Now the centres */
  cameraCentre.z = cameraHeight.z*0.5f; /* Touches z=0 */
  focalPlaneCentre.z = - focalPlaneHeight.z*0.5f; /* Also touches z=0 */
  /* Pythagoras and similar triangles */
  focalPlaneCentre.y = (float) sqrt(focaldepth*focaldepth -
				    (cameraCentre.z+focalPlaneCentre.z)*
				    (cameraCentre.z+focalPlaneCentre.z));


  for (imgy = height; imgy-- > 0; )
  {
    if (0.0f < previewy)
    {
      doxpreview = 1;
      previewx = -0.5f;
    }
    else
    {
      doxpreview = 0;
    }
    for (imgx = 0; imgx < width; imgx++)
    {
      screen[imgx + width*imgy].r = 0.0f;
      screen[imgx + width*imgy].g = 0.0f;
      screen[imgx + width*imgy].b = 0.0f;
      for (samplecounter=numsamples; samplecounter--;)
      {
	float rayLength;
	float angle, distance;
	float trigged;
	r.origin.x = cameraCentre.x;
	r.origin.y = cameraCentre.y;
	r.origin.z = cameraCentre.z;
	/* Generate random point within a circle of diameter lenswidth */
	angle = (M_PI * rand())/RAND_MAX;
	if (dozbuffer) distance = 0.0f;
	else distance = (float) sqrt((float)rand()/RAND_MAX) * 0.5;
	trigged = distance * cos(angle);
	r.origin.x += trigged*cameraWidth.x;
	r.origin.y += trigged*cameraWidth.y;
	r.origin.z += trigged*cameraWidth.z;
	trigged = distance * sin(angle);
	r.origin.x += trigged*cameraHeight.x;
	r.origin.y += trigged*cameraHeight.y;
	r.origin.z += trigged*cameraHeight.z;
	/* Ray direction = point on focal plane - ray origin */
	r.direction.x = focalPlaneCentre.x;
	r.direction.y = focalPlaneCentre.y;
	r.direction.z = focalPlaneCentre.z;
	r.direction.x += (imgx * focalPlaneWidth.x)/width - focalPlaneWidth.x*0.5f;
	r.direction.y += (imgx * focalPlaneWidth.y)/width - focalPlaneWidth.y*0.5f;
	r.direction.z += (imgx * focalPlaneWidth.z)/width - focalPlaneWidth.z*0.5f;
	r.direction.x += (imgy * focalPlaneHeight.x)/height - focalPlaneHeight.x * 0.5f;
	r.direction.y += (imgy * focalPlaneHeight.y)/height - focalPlaneHeight.y * 0.5f;
	r.direction.z += (imgy * focalPlaneHeight.z)/height - focalPlaneHeight.z * 0.5f;
	r.direction.x -= r.origin.x;
	r.direction.y -= r.origin.y;
	r.direction.z -= r.origin.z;
	/* Life is slightly easier in the sphere intersection if we */
	/* normalise the ray intersection. */
	rayLength = 1.0f / (float) sqrt(r.direction.x*r.direction.x+
					r.direction.y*r.direction.y+
					r.direction.z*r.direction.z);
	r.direction.x *= rayLength;
	r.direction.y *= rayLength;
	r.direction.z *= rayLength;
	lambda = intersect_ray(&r);
	/* No contribution if the ray missed everything */
	if (lambda != 0.0f)
	{
	  Vector P, normal;
	  float nsize;
	  float intensity;
	  const Vector lightDirection = { -2.0f, -2.0f, 2.0f };
	  const float roughness=0.02f;
	  Vector H; /* For specular calculation */
	  Colour shadingResult;
	  /* Hit a sphere - shade it */
	  /* We have lambda, get sphere normal */
	  /* P is the intersection point */
	  P.x = r.origin.x + r.direction.x * lambda;
	  P.y = r.origin.y + r.direction.y * lambda;
	  P.z = r.origin.z + r.direction.z * lambda;
	  normal.x = P.x - (floor(P.x) + 0.5f);
	  normal.y = P.y - (floor(P.y) + 0.5f);
	  normal.z = P.z - (floor(P.z) + 0.5f);
	  nsize = normal.x*normal.x +
	    normal.y*normal.y +
	    normal.z*normal.z;
	  nsize = 1.0f / sqrt(nsize);
	  normal.x *= nsize;
	  normal.y *= nsize;
	  normal.z *= nsize;
	  nsize = (float) sqrt(lightDirection.x*lightDirection.x +
			       lightDirection.y*lightDirection.y +
			       lightDirection.z*lightDirection.z);
	  /* Hall shading model (for now) */
	  intensity =
	    (normal.x * lightDirection.x +
	     normal.y * lightDirection.y +
	     normal.z * lightDirection.z) / nsize;
	  /* fprintf(stderr,"%f\n", intensity); */
	  shadingResult.r = intensity<0.0f?0.0f:intensity;
	  shadingResult.g = 0.0f;
	  shadingResult.b = 0.0f;
	  /* Iff we can see the light... */
	  if (intensity > 0.0f)
	  {
	    /* Add RenderMan-style specular highlight */
	    H.x = lightDirection.x/nsize-r.direction.x;
	    H.y = lightDirection.y/nsize-r.direction.y;
	    H.z = lightDirection.z/nsize-r.direction.z;
	    /* Normalise */
	    nsize = (float) sqrt(H.x*H.x + H.y*H.y + H.z*H.z);
	    nsize = 1.0f/nsize;
	    H.x *= nsize; H.y *= nsize; H.z *= nsize;
	    intensity = H.x*normal.x + H.y*normal.y + H.z*normal.z;
	    if (intensity > 0.0f)
	    {
	      intensity = pow(intensity, 1.0f/roughness);
	      /* Cyan hightlight (so red doesn't saturate and look weird) */
	      shadingResult.g += intensity;
	      shadingResult.b += intensity;
	    }
	  }
	  screen[imgx + width*imgy].r += shadingResult.r/numsamples;
	  screen[imgx + width*imgy].g += shadingResult.g/numsamples;
	  screen[imgx + width*imgy].b += shadingResult.b/numsamples;
	}
      }
      if (dozbuffer)
      {
	zbuffer[imgx + width*imgy] = lambda;
      }
      if (doxpreview)
      {
	if (0.0f < previewx)
	{
	  float intensity;
	  intensity = (1.0f/3.0f) *
	    (screen[imgx + width*imgy].r +
	     screen[imgx + width*imgy].g +
	     screen[imgx + width*imgy].b);
	  printf("%c", " .-+camHM"[clamp((int)(8.0f*intensity),0,8)]);
	  previewx -= 1.0f;
	}
	previewx += previewxinc;
      }
    }
    if (0.0f < previewy)
    {
      previewy -= 1.0f;
      printf("\n");
    }
    previewy += previewyinc;
  }
  if (dozbuffer)
  {
    fprintf(stderr, "Applying DOF post-process\n");
    cause_confusion(screen,
		    lenswidth * width/screenwidth,
		    focaldepth,
		    zbuffer,
		    width, height);
  }
  output = fopen(filename,"wb");
  if (output)
  {
    fprintf(stderr, "Writing output\n");
    dump_ppm(output, screen, width, height);
    fclose(output);
  }
  return 0;
}