#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#define GL_GLEXT_PROTOTYPES     //allow GL extensions, shaders, FBOs etc
#include <GL/gl.h>
#include <GL/glu.h>
#include <GL/glext.h>

#ifdef OSMESA
#include <GL/osmesa.h>
#else
#include <GL/glut.h>
#endif

#define frand (rand() / (float) RAND_MAX)

void reshape(int w, int h);
void draw(void);
void particleTex(void);
void snapshot(int width, int height, const char* path); 
void drawPoints(int dir);

GLuint tex;

/* view params */
float modelView[16];

#define FOV 45.0 //60.0        /* Field of view */
float window_h = 600; /* Window height */
float window_w = 800; /* Window width */
float particleRadius = 1.0f;  /* Particle radius */
int numParticles = 1000;  
//int numParticles = 1000000;  


/* Particle, vertex xyz, distance from eye */
typedef struct {
   float x, y, z;
   float distance;
} Particle;

/* Particle data */
Particle *particles;

int compare_dist(const void *a, const void *b)
{
   Particle *p1 = (Particle*)a;
   Particle *p2 = (Particle*)b;
   if (p1->distance == p2->distance) return 0;
   return p1->distance < p2->distance ? -1 : 1;
}

int main(int argc, char **argv) {
#ifdef OSMESA
	// initialize OSMesa
   void * pixelBuffer = malloc(4 * window_w * window_h);
	OSMesaContext osContext = OSMesaCreateContextExt( OSMESA_RGBA, 24, 0, 0, NULL);
	if (osContext == 0)
	{
		fprintf(stderr, "Error creating OSMesa context!\n");
		exit(1);
	}
	if (OSMesaMakeCurrent(osContext, pixelBuffer, GL_UNSIGNED_BYTE, window_w, window_h ) != GL_TRUE)
	{
		fprintf(stderr, "Error binding OSMesa context!\n");
		exit(1);
	}
#else
	// initialize glut
	glutInit(&argc, argv);
	glutInitDisplayMode(GLUT_RGBA | GLUT_ALPHA | GLUT_DEPTH | GLUT_DOUBLE);
	glutInitWindowSize(window_w, window_h);
	glutCreateWindow("Particle Test");
   //fprintf(stderr, "Window created\n");
#endif

   // GL Init stuff
   glEnable(GL_DEPTH_TEST);
   glClearColor(1, 1, 1, 0);

   //Transparency
   glEnable(GL_BLEND);
   glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);

   // Particle positions
   particles = malloc(sizeof(Particle) * numParticles);
   int i;
   for(i=0; i < numParticles; i++) 
   {
      //particles[i].x = (3 * (frand - 0.5f));
      //particles[i].y = (1 * (frand - 0.5f));
      //particles[i].z = (3 * (frand - 0.5f));
      particles[i].x = 0; 
      particles[i].y = 1.0;
      particles[i].z = 3 * (i / (float)numParticles) - 0.5;
   }

   particleTex();

   reshape(window_w, window_h);
	draw();
   /* Save image */
   snapshot(window_w, window_h, "snapshot.tga");
   if (particles) free(particles);

#ifdef OSMESA
   if (osContext)
   {
		OSMesaDestroyContext( osContext );
		free(pixelBuffer);
   }
#endif
   return 0;
}

void reshape(int w, int h)
{
   glMatrixMode(GL_PROJECTION);
   glLoadIdentity();
   gluPerspective((float)FOV, (float) w / (float) h, 0.5, 10.0);
   glMatrixMode(GL_MODELVIEW);
   glViewport(0, 0, w, h);
   window_h = h;
   window_w = w;
}

// Renderer
void draw(void) {
   glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

   // View transform
   glMatrixMode(GL_MODELVIEW);
   glLoadIdentity();

   glTranslatef(0, -0.1, -3.6);
   glRotatef(25, 1.0, 0.0, 0.0);
   //glRotatef(45, 0.0, 1.0, 0.0);

   // Update and depth sort particles..
   fprintf(stderr, "Calculating distances\n");
   glGetFloatv(GL_MODELVIEW_MATRIX, modelView);
   float eye[3];
   int i;
   for (i = 0; i < numParticles; i++) {
      //Get eye pos vector by multiplying particle pos by model view matrix
//      eye[0] = modelView[0] * particles[i].x  + modelView[4] * particles[i].y
//                + modelView[8] * particles[i].z  + modelView[12]; // * w = 1.0
//      eye[1] = modelView[1] * particles[i].x  + modelView[5] * particles[i].y 
//                + modelView[9] * particles[i].z  + modelView[13]; // * w
      eye[2] = modelView[2] * particles[i].x  + modelView[6] * particles[i].y 
                + modelView[10] * particles[i].z  + modelView[14]; // * w

      //Calculate distance from viewing plane (negative z-component of vector: Z.v where Z is Z-axis vector 0,0,-1)
      particles[i].distance = -eye[2]; //Z.v = Z.x*v.x + Z.y*v.y + Z.z*v.z = 0 + 0 + -1*v.z
   }

   fprintf(stderr, "Depth sorting\n");
   qsort(particles, numParticles, sizeof(Particle), compare_dist);

   // Normal blending for rgb data, accumulate opacity in alpha channel with additive blending
   glBlendFuncSeparate(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA, GL_ONE, GL_ONE_MINUS_SRC_ALPHA); 

   glAlphaFunc(GL_GREATER, 0.1); 
   glEnable(GL_ALPHA_TEST);
 
   fprintf(stderr, "Rendering\n");
   // Textured point sprites 
   glEnable(GL_POINT_SPRITE);
   //Diffuse texture
   glActiveTexture(GL_TEXTURE0);
   glEnable(GL_TEXTURE_2D);
   glBindTexture(GL_TEXTURE_2D, tex);


   glTexEnvi(GL_POINT_SPRITE, GL_COORD_REPLACE, GL_TRUE);
    
   //Point size attenuation
   glPointSize(particleRadius);   
   static GLfloat attenuate[3] = { 0.0, 0.0, 0.05 };  //Const, linear, quadratic
   glPointParameterfv(GL_POINT_DISTANCE_ATTENUATION, attenuate);
   //glPointParameterf(GL_POINT_FADE_THRESHOLD_SIZE, 60.0);
 	//glPointParameterf( GL_POINT_SIZE_MIN, 1.0f );
   //glPointParameterf( GL_POINT_SIZE_MAX, 64.0f );

   //bp   glColor4f(1.0, 0.0, 0.0, 0.5);
   glColor4f(1.0, 0.0, 0.0, 1.0);
   glBegin(GL_POINTS);
   for(i=numParticles-1; i>=0; i--) //Reverse order farthest to nearest
      glVertex3f(particles[i].x, particles[i].y, particles[i].z);
   glEnd();

   //Restore state
   glDisable(GL_POINT_SPRITE);
   glDisable(GL_ALPHA_TEST);

}

void particleTex()
{
   //Pre-calculate particle lighting model, save in texture
   int x, y;
   float data[64][64][2];
   for (x=0; x<64; x++)
   {
      for (y=0; y<64; y++)
      {
         float Nx, Ny, Nz, mag, alpha, diffuse;

         Nx = (x/64.0) * 2.0 - 1.0;
         Ny = (y/64.0) * -2.0 + 1.0;
         mag = Nx * Nx + Ny * Ny;
         //printf("Normal at (%d,%d) = (%f,%f) |N| = %f\n", x, y, Nx, Ny, mag);
         diffuse = 0;
         alpha = 0.0;
         if (mag <= 1.0)
         {
            Nz = sqrt(1.0 - mag);
            diffuse = 0.577 * (Nx + Ny + Nz);   //Light pos (1,1,1) normalised = (0.577,0577,0.577)
            if (diffuse < 0.0) 
               diffuse = 0;
            alpha = 1.0;
         }
         data[y][x][0] = diffuse; 
         data[y][x][1] = alpha;
      }    
   }
   glGenTextures(1, &tex);
   glActiveTexture(GL_TEXTURE0);
   glBindTexture(GL_TEXTURE_2D, tex);
   glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
   glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
//    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR_MIPMAP_LINEAR);
//  glTexParameteri(GL_TEXTURE_2D, GL_GENERATE_MIPMAP, GL_TRUE);    //set so texImage2d will gen mipmaps
   glTexImage2D(GL_TEXTURE_2D, 0, GL_LUMINANCE_ALPHA, 64, 64, 0, GL_LUMINANCE_ALPHA, GL_FLOAT, data);
}

void snapshot(int width, int height, const char* path) 
{
   /* WRITE TGA */
   GLubyte *image = (GLubyte *) malloc(width * height * sizeof(GLubyte) * 4);
	glPixelStorei(GL_PACK_ALIGNMENT, 1);
	glReadPixels(0, 0, width, height, GL_BGRA, GL_UNSIGNED_BYTE, image);
	FILE *fp = fopen(path, "wb");
	if (!fp) {printf("[write_tga_file] File could not be opened for writing"); return;}
   //write header
   char tga_header[18] = {0,0,2,0,0,0,0,0,0,0,0,0,
                         width & 255, (width >> 8) & 255,
                         height & 255, (height >> 8) & 255,
                         32,0};
   fwrite(tga_header, 1, 18, fp);
   //write image data
   fwrite(image, sizeof(GLubyte) * 4, width * height, fp);
   fclose(fp);
   free(image);
}

