/*
 * A 5-Chemical Reaction-Diffusion simulator for PLY objects. (rdstripes)
 *
 *		for usage info, run program without any parameters
 *
 */

#include <GL/glut.h>
#include <stdlib.h>
#include <stdio.h>
#include "rply.h"

/***** extra stuff for mouse picking */
/* size of mouse picking selection buffer */
#define SELECT_MAX 250

/* window size (its square) */
GLint winx = 250, winy = 250;

/* the buffer itself */
GLuint selectBuf[SELECT_MAX];

/* viewport array */
GLint vp[4];

/* array of initiator cells */
long initiators[10];
/* index to that array and number of initiators */
int initiatorindex = 0, maxinitiators = 10;

int runningRD = 0;

/***** end of mouse picking stuff */

/* R-D parameters and arrays */
GLfloat p1 = 0.04;
GLfloat p2 = 0.06;
GLfloat p3 = 0.04;

GLfloat Dg = 0.05; //0.1; //0.009;
GLfloat Ds = 0.2; //0.2;
GLfloat speed = 1.0;

GLfloat kg,kr,ks;

GLfloat g1rand = 0.02;

GLfloat *BETA;
GLfloat *g1;
GLfloat *g2;
GLfloat *r;
GLfloat *s1;
GLfloat *s2;
GLfloat *dg1;
GLfloat *dg2;
GLfloat *dr;
GLfloat *ds1;
GLfloat *ds2;
GLfloat *activator;

/* 3d object data arrays */
GLfloat **vertices;
GLfloat **normals;
GLfloat **diffusion;
long **faces;
long **neighbours;

/* just some vars for extracting data from PLY (could be done better) */
long vertex=0, axis=0, face=0, fvertex=0, home=0, neighb=0;
GLfloat maxv = -100;
GLfloat minv = 100;

/* angles and stuff for rotating and zooming the object */
GLfloat xangle = 0.0;
GLfloat yangle = 0.0;
GLfloat zangle = 0.0;
GLfloat zoom = 1.5;
GLfloat atX = 0.0;
GLfloat atY = 0.0;

/* global for number of vertices and faces */
long nvertices,nfaces;
e_ply_type vertex_t, face_length_t, face_value_t;

/* GLUT window identifier and viewport size (so we see whole object)*/
static int win;
GLfloat view;

/* input PLY filename */
char *filename;
char *nbsfile;

/* just some global variables for stuff that should be implemented better */
int plyHasNeighbs=0, plyHasNormals=0, plyHasDiff=0;


/*
 *
 * Display usage info
 *
 */
void printUsage()
{
	printf("Reaction-Diffusion Simulator (Stripes)\n\nusage:\trdstripes <plyfile> [-p <paramfile>]\n\n");
	printf("\tOptions:\n\t-p\t\t- Use R-D parameters in <paramfile>\n\t\t\t  [If not specified defaults are used]\n\n");
	printf("\tControls:\n\t\tRotate on X axis: 'a' and 'd'\n\t\tRotate on Y axis: 's' and 'w'\n\t\tRotate on Z axis: 'z' and 'c'\n\t\tZoom in and out: '+', '-'\n\t\tMove object: Left, Right, Up and Down arrows\n\t\tReset view: 'r'\n\t\tExport concentration levels: 'e'\n\t\tSelect stripe initiator cells: Mouse\n\t\tPause/resume simulation: 'p'\n\t\tQuit: 'q'\n\n");
	printf("\tParameter file format:\n\t\t<p1 value>\t\t- ~0.04\n\t\t<p2 value>\t\t- ~0.06\n\t\t<p3 value>\t\t- ~0.04\n\t\t<Dg value>\t\t- ~0.05\n\t\t<Ds value>\t\t- ~0.2\n\n");
	printf("If the PLY input file does not contain additional data about vertex normals and neighbour cells, this will be generated before R-D simulation begins. In this case a PLY file with the extra data is stored at <plyfile>_nbs.ply, and this file should be used in future to save on computation time.\n\n");
}


/*
 *
 * The Reaction-Diffusion Jazz
 *
 */
void updateConc()
/*
 * updateConc(): calculate the new values for each of the chemicals involved using
 * Meinhardt's R-D equations.
 */
{
	long i;
	int j, frozen=0;
	GLfloat ddg1,ddg2,dds1,dds2;
	GLfloat temp1, temp2, c=0.01;
	
	/* run through each triangular face (cell) */
	for (i=0;i<nfaces;i++) {
		/* calculate the hell out of this stuff (differentials) */
		ddg1 = diffusion[i][0] * g1[neighbours[i][0]] + diffusion[i][1] * g1[neighbours[i][1]] + diffusion[i][2] * g1[neighbours[i][2]];
		ddg2 = diffusion[i][0] * g2[neighbours[i][0]] + diffusion[i][1] * g2[neighbours[i][1]] + diffusion[i][2] * g2[neighbours[i][2]];
		dds1 = diffusion[i][0] * s1[neighbours[i][0]] + diffusion[i][1] * s1[neighbours[i][1]] + diffusion[i][2] * s1[neighbours[i][2]];
		dds2 = diffusion[i][0] * s2[neighbours[i][0]] + diffusion[i][1] * s2[neighbours[i][1]] + diffusion[i][2] * s2[neighbours[i][2]];
		
		temp1 = c * g1[i] * g1[i] * s2[i];
		/* if no stripe initiators are set we introduce some randomness to cause stripes to form */
		if (maxinitiators == 0) temp1 *= BETA[i];
		temp2 = c * g2[i] * g2[i] * s1[i];
		
		/* put it all together */
		dg1[i] = g1[i] * kg + Dg * ddg1 + temp1/r[i];
		dg2[i] = g2[i] * kg + Dg * ddg2 + temp2/r[i];
		dr[i] = r[i] * kr + temp1 + temp2;
		ds1[i] = s1[i] * ks + Ds * dds1 + p3 * g1[i];
		ds2[i] = s2[i] * ks + Ds * dds2 + p3 * g2[i];
	}
	
	/* put the diffs into play */
	for (i=0;i<nfaces;i++) {
		frozen = 0;
		/* cell initators are frozen */
		for (j=0;j<maxinitiators;j++)
			if (initiators[j] == i) frozen = 1;
			
		if (!frozen) {
			g1[i] += (speed * dg1[i]);
			g2[i] += (speed * dg2[i]);
			r[i] += (speed * dr[i]);
			s1[i] += (speed * ds1[i]);
			s2[i] += (speed * ds2[i]);
		}
		
		activator[i] = g1[i];
	}
}

void setStart()
/*
 * setStart(): setup the chemical concentration levels to the initial state
 */
{
	long i;
	for(i=0;i<nfaces;i++) {
		g1[i] = p2/(2*p1);
		g2[i] = g1[i];
		r[i] = 0.02 * g1[i] * g1[i] * g1[i] / p2;
		s1[i] = g1[i];
		s2[i] = g1[i];
		BETA[i] = (g1rand * rand()/(RAND_MAX+1.0))+(1 - g1rand/2);
		
		activator[i] = g1[i];
	}
	
	/* set the k constants */
	kg = -p1 - 3 * Dg;
	kr = -p2;
	ks = -p3 - 3 * Ds;
}


/*
 *
 * Output the concentration levels to a file. This is for use in rdviewer.
 *
 */
void outputConc(void)
{
	long i;
	char outfilename[20];
	FILE *fp;
	
	printf("Enter filename for pattern export:\n\t");
	scanf("%s",outfilename);
	
	if (fp = fopen(outfilename,"w")) {
		printf("Exporting concentration levels to %s... ",outfilename);
		if ((!plyHasNeighbs) || (!plyHasNormals))
			fprintf(fp, "%s\n", nbsfile);
		else
			fprintf(fp, "%s\n", filename);
		
		for (i=0;i<nfaces;i++)
			fprintf(fp, "%f\n", activator[i]);
		
		fclose(fp);
		printf("Done\n");
	} else
		printf("Error writing to file pattern.cnc.\n");
}


/*
 *
 * GLUT routines, callbacks and whatnot
 *
 */

void update(void)
/*
 * update(): GLUT callback function. called every refresh. We use it to run R-D
 * simulation then call for a redisplay.
 */
{
	if (runningRD) updateConc();
		
	glutPostRedisplay();
}


void face3(long face)
/*
 * face3(face): Render the triangular polygon indexed to in the faces array by
 * input parameter face. Colour is set according to activator[] value for the face.
 * The glMaterial and glNormal calls are required for lighting.
 */
{
	GLfloat colour[] = {activator[face], activator[face], activator[face]};
	GLfloat material[] = {activator[face], activator[face], activator[face], 1.0};
	
	glColor3fv(colour);
	glMaterialfv(GL_FRONT, GL_SPECULAR, material);
	glMaterialfv(GL_FRONT, GL_AMBIENT, material);
	glMaterialfv(GL_FRONT, GL_DIFFUSE, material);
	glBegin(GL_POLYGON);
		glNormal3fv(normals[faces[face][0]]);
		glVertex3fv(vertices[faces[face][0]]);
		glNormal3fv(normals[faces[face][1]]);
		glVertex3fv(vertices[faces[face][1]]);
		glNormal3fv(normals[faces[face][2]]);
		glVertex3fv(vertices[faces[face][2]]);
	glEnd();
}

void display(void)
/*
 * display(): GLUT display callback, it sets up the view, and rotation of the scene,
 * then it renders all the polygons, and swaps buffers (we're using double-buffering).
 */
{
	long i;
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	
	/* Set up the modelview transformation. */
	glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();
	//gluLookAt(6.0, 4.0, 14.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0);
	gluLookAt(0, 0, zoom * view, atX, atY, 0.0, 0.0, 1.0, 0.0);
	glRotatef(xangle, 1.0, 0.0, 0.0);
	glRotatef(yangle, 0.0, 1.0, 0.0);
	glRotatef(zangle, 0.0, 0.0, 1.0);
	
	for (i=0;i<nfaces;i++)
		face3(i);
	
	glutSwapBuffers();
}

void reshape(int w, int h)
/*
 * reshape(w, h): this is needed for when the window is resized, it keeps the
 * shape and perspective of the viewport correct.
 */
{
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	gluPerspective(40.0, (GLdouble)w/(GLdouble)h, 0.01, 2000.0);
	glViewport(0, 0, (GLsizei)w, (GLsizei)h);
	
	winx = w;
	winy = h;
	glGetIntegerv(GL_VIEWPORT, vp);
}

/* cell selection using the mouse (experimental) */
GLint selectcell(GLint x, GLint y)
{
  GLint hits;
	long i,cell;
	int depth;

  glSelectBuffer(SELECT_MAX, selectBuf);
  glRenderMode(GL_SELECT);
  glInitNames();
  glPushName(~0);
	
	glMatrixMode(GL_PROJECTION);
	glPushMatrix();
	glLoadIdentity();
	gluPickMatrix(x, winy - y, 3, 3, vp);
	gluPerspective(40.0, (GLdouble)winx/(GLdouble)winy, 0.01, 2000.0);
	glViewport(0, 0, (GLsizei)winx, (GLsizei)winy);
  
  glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();
	gluLookAt(0, 0, zoom * view, atX, atY, 0.0, 0.0, 1.0, 0.0);
	glRotatef(xangle, 1.0, 0.0, 0.0);
	glRotatef(yangle, 0.0, 1.0, 0.0);
  glRotatef(zangle, 0, 0, 1);

	for (i=0;i<nfaces;i++) {
		glLoadName(i);
		face3(i);
	}
	
	glMatrixMode(GL_PROJECTION);
	glPopMatrix();

  hits = glRenderMode(GL_RENDER);
  if (hits <= 0) {
		printf("\tNo cell found there. Try again (maybe zoom in first).\n");
    return -1;
  }
	
	depth = selectBuf[1];
	cell = selectBuf[3];
	for (i=0;i<hits;i++) {
		if (selectBuf[i*4+1] < depth) {
			cell = selectBuf[i*4+3];			// Select The Closer Object
			depth = selectBuf[i*4+1];			// Store How Far Away It Is
		}
	}
	return cell;
	
	//printf("hits: %d = %d\n",hits,selectBuf[(hits - 1) * 4 + 3]);
  //return selectBuf[(hits - 1) * 4 + 3];
}


void keyb(unsigned char key, int x, int y)
/*
 * keyb(key, x, y): GLUT keyboard callback, we use it to rotate the scene, zoom in
 * and out, perform other user functions, and finally to free all the malloced memory.
 */
{
	long i;
	
	switch(key) {
		case 'a':
			yangle -= 10;
			break;
		case 'd':
			yangle += 10;
			break;
		case 'w':
			xangle -= 10;
			break;
		case 's':
			xangle += 10;
			break;
		case 'c':
			zangle -= 10;
			break;
		case 'z':
			zangle += 10;
			break;
		case '+':
		case '=':
			zoom -= 0.2;
			break;
		case '-':
			zoom += 0.2;
			break;
		case 'r':
			xangle = yangle = zangle = atX = atY = 0.0;
			zoom = 1.5;
			break;
		case 'e':
			outputConc();
			break;
		case ' ':
		case 'p':
			if (initiatorindex >= maxinitiators) {
				if (!runningRD) {
					runningRD = 1;
					printf("R-D simulation running...\n");
				} else {
					runningRD = 0;
					printf("R-D simulation paused. Press 'p' to resume.\n");
				}
			}
			break;
		case 'q':
		case 'Q':	
			printf("Freeing memory\n");
			for(i=0;i<nvertices;i++) {
				free(vertices[i]);
				free(normals[i]);
			}
			free(vertices);
			free(normals);
			for(i=0;i<nfaces;i++) {
				free(faces[i]);
				free(neighbours[i]);
				free(diffusion[i]);
			}
			free(faces);
			free(neighbours);
			free(diffusion);
			
			free(BETA);
			free(g1);
			free(g2);
			free(r);
			free(s1);
			free(s2);
			free(dg1);
			free(dg2);
			free(dr);
			free(ds1);
			free(ds2);
			free(activator);
			
			printf("All done.\n");
			glutDestroyWindow(win);
			exit(0);
	}
}

void special(int key, int x, int y)
/*
 * special(key, x, y): GLUT special keyboard callback, used for arrow keys to move scene.
 */
{
	switch(key) {
		case GLUT_KEY_RIGHT:
			atX=atX-0.02*view;
			break;
		case GLUT_KEY_LEFT:
			atX=atX+0.02*view;
			break;
		case GLUT_KEY_UP:
			atY=atY-0.02*view;
			break;
		case GLUT_KEY_DOWN:
			atY=atY+0.02*view;
			break;
		default:
			return;
	}
}

void mouse(int button, int state, int mouseX, int mouseY)
/*
 * mouse(button, state, mouseX, mouseY): GLUT mouse callback, used to select stripe
 * initiator cells
 */
{
  GLint hit;

	if (!runningRD) {
		if (state == GLUT_DOWN) {
			hit = selectcell((GLint) mouseX, (GLint) mouseY);
			if ((hit != -1) && (initiatorindex < maxinitiators)) {
				if (button == GLUT_LEFT_BUTTON) {
					printf("\tSetting cell %d as white stripe initiator cell. Choose %d more...\n",hit,maxinitiators-initiatorindex-1);
					initiators[initiatorindex++] = hit;
					g1[hit] = activator[hit] = 255;
				} else if (button == GLUT_RIGHT_BUTTON) {
					printf("\tSetting cell %d as black stripe initiator cell. Choose %d more...\n",hit,maxinitiators-initiatorindex-1);
					initiators[initiatorindex++] = hit;
					g1[hit] = activator[hit] = 0;
				}
				glutPostRedisplay();
			}
			if (initiatorindex >= maxinitiators)
				printf("Stripe initiators set. Press the spacebar to begin R-D simulation...\n");
		}
	}
}
	
/*
 *
 * All the PLY model handling routines
 *
 */
static int vertex_cb(p_ply_argument argument)
/*
 * vertex_cb(argument): RPly vertex callback. This is called everytime RPly encounters
 * a vertex as it parses the PLY file. We store the vertex coordinates and (rather crudely)
 * find the max and min coordinates, for use in setting the viewport correctly.
 */
{
	long eol;
	p_ply_property property;
	
	ply_get_argument_user_data(argument, NULL, &eol);
	ply_get_argument_property(argument, &property, NULL, NULL);
	ply_get_property_info(property, NULL, &vertex_t, NULL, NULL);
	
	vertices[vertex][axis] = (GLfloat)ply_get_argument_value(argument);

	if (vertices[vertex][axis] > maxv)
		maxv = vertices[vertex][axis];
	if (vertices[vertex][axis] < minv)
		minv = vertices[vertex][axis];
	
	axis++;
	if (eol) {
		vertex++;
		axis=0;
	}
	return 1;
}

static int normal_cb(p_ply_argument argument)
/*
 * normal_cb(argument): RPly normal callback. Store the normal for each vertex.
 */
{
	long eol;
	
	ply_get_argument_user_data(argument, NULL, &eol);
	normals[vertex-1][axis] = (GLfloat)ply_get_argument_value(argument);
	
	axis++;
	if (eol) {
		//vertex++;
		axis=0;
	}
	return 1;
}

static int face_cb(p_ply_argument argument)
/*
 * face_cb(argument): RPly face callback. Store vertex indices for each face.
 */
{
	long length, value_index;
	p_ply_property property;
	
	ply_get_argument_property(argument, &property, &length, &value_index);
	ply_get_property_info(property, NULL, NULL, &face_length_t, &face_value_t);
		
  switch (value_index) {
		case 0:
		case 1: 
			faces[face][fvertex++] = (long)ply_get_argument_value(argument);
			break;
		case 2:
			faces[face][fvertex++] = (long)ply_get_argument_value(argument);
			face++;
			fvertex=0;
			break;
		default: 
			break;
	}
	return 1;
}

static int neighbours_cb(p_ply_argument argument)
/*
 * neighbours_cb(argument): RPly neighbour callback. Retrieves the neighbour indices.
 * home variable indexs current cell, neighb indexs its neighbours. Bit of a cludge to
 * cover the case where neighbours ar found but not diffusion weightings.
 */
{
	long value_index;
	
	ply_get_argument_property(argument, NULL, NULL, &value_index);
  
  switch (value_index) {
		case 0:
		case 1: 
			neighbours[home][neighb++] = (long)ply_get_argument_value(argument);
			break;
		case 2:
			neighbours[home][neighb++] = (long)ply_get_argument_value(argument);
			if (plyHasDiff == 0) home++;
			neighb=0;
			break;
		default:
			break;
	}
	return 1;
}

static int diffusion_cb(p_ply_argument argument)
/*
 * diffusion_cb(argument): RPly diffusion weightings callback.
 */
{
	long value_index;
	
	ply_get_argument_property(argument, NULL, NULL, &value_index);
  
  switch (value_index) {
		case 0:
		case 1: 
			diffusion[home][neighb++] = (GLfloat)ply_get_argument_value(argument);
			break;
		case 2:
			diffusion[home][neighb++] = (GLfloat)ply_get_argument_value(argument);
			home++;
			neighb=0;
			break;
		default:
			break;
	}
	return 1;
}

int outputPLY(char *filename)
/*
 * outputPLY(filename): calls all the RPly output functions in correct order to write
 * a PLY file with all additional data (neighbours, normals, diffusion weightings).
 */
{
	int i,j;
	p_ply ply;

	ply = ply_create(filename, PLY_ASCII, NULL);
	
	/* add vertex element header info */
	if (!ply_add_element(ply, "vertex", nvertices)) printf("cant add vertex element\n");
	ply_add_scalar_property(ply, "x", vertex_t);
	ply_add_scalar_property(ply, "y", vertex_t);
	ply_add_scalar_property(ply, "z", vertex_t);
	ply_add_scalar_property(ply, "xnorm", vertex_t);
	ply_add_scalar_property(ply, "ynorm", vertex_t);
	ply_add_scalar_property(ply, "znorm", vertex_t);

	/* add faces element info */
	if (!ply_add_element(ply, "face", nfaces)) printf("cant add face element\n");
	ply_add_list_property(ply, "vertex_indices", face_length_t, face_value_t);
	
	/* add neighbours element info */
	if (!ply_add_element(ply, "neighbours", nfaces)) printf("cant add neighbours element\n");
	ply_add_list_property(ply, "vertex_indices", face_length_t, face_value_t);
	ply_add_list_property(ply, "diffusion_weighting", face_length_t, vertex_t);
	
	/* add comment saying created by our prog and added neighbours */
	ply_add_comment(ply, "Extra neighbour and vertex normal data generated for Reaction-Diffusion simulations");
	
	/* write the header */
	if (!ply_write_header(ply)) printf("cant write header\n");
	
	/* write the properties one by one */
	for (i=0;i<nvertices;i++) {
		ply_write(ply, vertices[i][0]);
		ply_write(ply, vertices[i][1]);
		ply_write(ply, vertices[i][2]);
		ply_write(ply, normals[i][0]);
		ply_write(ply, normals[i][1]);
		ply_write(ply, normals[i][2]);
	}
	
	for (i=0;i<nfaces;i++) {
		ply_write(ply, 3);
		ply_write(ply, faces[i][0]);
		ply_write(ply, faces[i][1]);
		ply_write(ply, faces[i][2]);
	}
	
	for (i=0;i<nfaces;i++) {
		ply_write(ply, 3);
		ply_write(ply, neighbours[i][0]);
		ply_write(ply, neighbours[i][1]);
		ply_write(ply, neighbours[i][2]);
		ply_write(ply, 3);
		ply_write(ply, diffusion[i][0]);
		ply_write(ply, diffusion[i][1]);
		ply_write(ply, diffusion[i][2]);
	}
	
	ply_close(ply);
	printf("PLY file with neighbour and normal data written to %s\n",filename);
}

/*
 *
 * Neighbour finding and normal vector algorithms. Damn slow for big models.
 *
 */
void findNeighbours(void)
{
	long i,j,k,p,v,v1,v2,v3,v1st,v2nd,matches;
	GLfloat diff, norm;
	
	printf("Finding neighbours (may take some time)...\n");
	for(i=0;i<nfaces;i++) {
		v1 = faces[i][0];
		v2 = faces[i][1];
		v3 = faces[i][2];
		p = 0;
		for(j=0;j<nfaces;j++) {
			matches = 0;
			for(k=0;k<3;k++) {
				v = faces[j][k];
				if ((v == v1) || (v == v2) || (v == v3)) {
					if (matches == 0) v1st = v;
					if (matches == 1) v2nd = v;
					matches++;
				}
			}
			if (matches == 2) {
				neighbours[i][p] = j;
				diff = sqrt(((vertices[v1st][0] - vertices[v2nd][0]) * (vertices[v1st][0] - vertices[v2nd][0])) + ((vertices[v1st][1] - vertices[v2nd][1]) * (vertices[v1st][1] - vertices[v2nd][1])) + ((vertices[v1st][2] - vertices[v2nd][2]) * (vertices[v1st][2] - vertices[v2nd][2])));
				diffusion[i][p++] = diff;
			}
		}
		norm = 3/(diffusion[i][0] + diffusion[i][1] + diffusion[i][2]);
		diffusion[i][0] = norm * diffusion[i][0];
		diffusion[i][1] = norm * diffusion[i][1];
		diffusion[i][2] = norm * diffusion[i][2];
	}
	printf("Done.\n");
}


GLfloat faceNorm(long face, int axis)
{
	/* [v1 - v2] x [v2 - v3] */
	
	int i;
	GLfloat lhs[3];
	GLfloat rhs[3];
	
	for(i=0;i<3;i++) {
		lhs[i] = vertices[faces[face][0]][i] - vertices[faces[face][1]][i];
		rhs[i] = vertices[faces[face][1]][i] - vertices[faces[face][2]][i];
	}
	
	switch (axis) {
		case 0:
			return (lhs[1] * rhs[2] - lhs[2] * rhs[1]);
		case 1:
			return (lhs[2] * rhs[0] - lhs[0] * rhs[2]);
		case 2:
			return (lhs[0] * rhs[1] - lhs[1] * rhs[0]);
	}
}


void findNormals(void)
{
	long i,j;
	int count=0;
	GLfloat normal[3];
	GLfloat magnitude;
	
	printf("Calculating normals (may take some time)...\n");
	
	for(i=0;i<nvertices;i++) {
		count=0;
		normal[0] = normal[1] = normal[2] = 0.0;
		for(j=0;j<nfaces;j++) {
			if ((faces[j][0] == i) || (faces[j][1] == i) || (faces[j][2] == i)) {
				count++;
				normal[0] += faceNorm(j, 0);
				normal[1] += faceNorm(j, 1);
				normal[2] += faceNorm(j, 2);
			}
		}
		magnitude = sqrt((normal[0] * normal[0]) + (normal[1] * normal[1]) + (normal[2] * normal[2]));
		if (magnitude == 0) {
			normals[i][0] = normal[0];
			normals[i][1] = normal[1];
			normals[i][2] = normal[2];
		} else { 
			normals[i][0] = normal[0]/magnitude;
			normals[i][1] = normal[1]/magnitude;
			normals[i][2] = normal[2]/magnitude;
		}
	}
	printf("Done.\n");
}


/*
 *
 * Lots of initialisation stuff. Loading the object, etc.
 *
 */
void init(void)
{
	long i,j;
	p_ply ply;
	GLfloat mat_specular[] = { 0.8, 0.8, 0.8, 1.0 };
	GLfloat mat_shininess[] = { 10.0 };
	GLfloat light_ambient[] = {0.1, 0.1, 0.1, 0.5};
  GLfloat light_position[] = { 1.0, 1.0, 1.0, 0.0 };
	 
	/* some GL initialisation */
	glClearColor(0.0, 0.0, 0.5, 1.0);
	glEnable(GL_DEPTH_TEST);
	
	/* load PLY file */
	ply = ply_open(filename, NULL);
	if ((!ply) || (!ply_read_header(ply))) {
		printf("Error opening PLY file or reading header.\n");
   	exit;
	}
  
	/* set callbacks for PLY reading and get number of vertices*/
  nvertices = ply_set_read_cb(ply, "vertex", "x", vertex_cb, NULL, 0);
  ply_set_read_cb(ply, "vertex", "y", vertex_cb, NULL, 0);
  ply_set_read_cb(ply, "vertex", "z", vertex_cb, NULL, 1);
	
	/* test if PLY file contains vertex normals, if so set callbacks */
	if (ply_set_read_cb(ply, "vertex", "xnorm", normal_cb, NULL, 0)) {
		printf("This PLY file defines normals to vertices. Good work.\n");
		plyHasNormals = 1;
	}
	ply_set_read_cb(ply, "vertex", "ynorm", normal_cb, NULL, 0);
	ply_set_read_cb(ply, "vertex", "znorm", normal_cb, NULL, 1);
	
	/* set face callback and get number of faces */
  nfaces = ply_set_read_cb(ply, "face", "vertex_indices", face_cb, NULL, 0);
  /* try again with element label "vertex_index" if it fails */
	if (nfaces == 0) nfaces = ply_set_read_cb(ply, "face", "vertex_index", face_cb, NULL, 0);
	
	/* test for neighbour and diffusion weighting data and register callbacks */
	if (ply_set_read_cb(ply, "neighbours", "vertex_indices", neighbours_cb, NULL, 0)) {
		printf("The PLY file contains neighbour data. Hoorah!\n");
		plyHasNeighbs = 1;
	}
	if (ply_set_read_cb(ply, "neighbours", "diffusion_weighting", diffusion_cb, NULL, 0))
		plyHasDiff = 1;
	printf("Vertices: %ld\nFaces: %ld\n", nvertices, nfaces);
   
  /* mallocate memory for arrays */
  /* 3d structure data */
  vertices = malloc(nvertices * sizeof(GLfloat *));
	normals = malloc(nvertices * sizeof(GLfloat *));
  for(i=0;i<nvertices;i++) {
  	vertices[i] = malloc(3 * sizeof(GLfloat));
		normals[i] = malloc(3 * sizeof(GLfloat));
	}
  faces = malloc(nfaces * sizeof(long *));
  neighbours = malloc(nfaces * sizeof(long *));
	diffusion = malloc(nfaces * sizeof(GLfloat *));
  for(i=0;i<nfaces;i++) {
		faces[i] = malloc(3 * sizeof(long));
  	neighbours[i] = malloc(3 * sizeof(long));
		diffusion[i] = malloc(3 * sizeof(GLfloat));
  }
   
  /* r-d arrays */
  BETA = malloc(nfaces * sizeof(GLfloat));
  g1 = malloc(nfaces * sizeof(GLfloat));
  g2 = malloc(nfaces * sizeof(GLfloat));
	r = malloc(nfaces * sizeof(GLfloat));
	s1 = malloc(nfaces * sizeof(GLfloat));
	s2 = malloc(nfaces * sizeof(GLfloat));
  dg1 = malloc(nfaces * sizeof(GLfloat));
  dg2 = malloc(nfaces * sizeof(GLfloat));
	dr = malloc(nfaces * sizeof(GLfloat));
	ds1 = malloc(nfaces * sizeof(GLfloat));
	ds2 = malloc(nfaces * sizeof(GLfloat));
  activator = malloc(nfaces * sizeof(GLfloat));
  
  printf("Done mallocating memory.\n");
  
	/* a little bit of error checking */
  if (!ply_read(ply)) printf("Error extracting vertices.\n");
  ply_close(ply);
  
  /* making sure we get everything in view (very basic may need to be improved) */
  view = maxv - minv; 
  //printf("max: %.2f %.2f\nView: %.2f\n",maxv,minv,view);
	
	/* setup for lighting */
	glShadeModel (GL_SMOOTH);
	glMaterialfv(GL_FRONT, GL_SPECULAR, mat_specular);
  glMaterialfv(GL_FRONT, GL_SHININESS, mat_shininess);
	light_position[0] = light_position[1] = light_position[2] = -view;
  glLightfv(GL_LIGHT0, GL_POSITION, light_position);
	glLightfv(GL_LIGHT0, GL_AMBIENT, light_ambient);
  glEnable(GL_LIGHTING);
  glEnable(GL_LIGHT0);
	
	/* if we dont have neighbours or normals find them then save the ply file */
	if ((!plyHasNeighbs) || (!plyHasNormals)) {
		printf("PLY file does not contain neighbours and/or normals data.\n");
		if (!plyHasNeighbs) findNeighbours();
		if (!plyHasNormals) findNormals();
		
		/* set filename for outputting with nbs data (this is stupidly complicated) */
		nbsfile = (char *)calloc(strlen(filename)+5, sizeof(char));
		strncpy(nbsfile,filename,strlen(filename)-4);
		strcat(nbsfile, "_nbs.ply");
		
		/* output a PLY file */
		outputPLY(nbsfile);
	}
	
	/* set the initial chemical level values */
	setStart();
	
	printf("Initialisation done.\n");
	printf("Controls:\n\tRotate on X axis: 'a' and 'd'\n\tRotate on Y axis: 's' and 'w'\n\tRotate on Z axis: 'z' and 'c'\n\tZoom in and out: '+', '-'\n\tMove object: Left, Right, Up and Down arrows\n\tReset view: 'r'\n\tExport concentration levels: 'e'\n\tSelect stripe initiator cells: Mouse\n\tPause/resume simulation: 'p'\n\tQuit: 'q'\n");
	
	/* determine how many stripe initiator cells are required, if none start R-D going */
	printf("\nHow many stripe initiator cells? [0-10] [0 = random stripes] ");
	scanf("%d", &maxinitiators);
	if (maxinitiators > 0) {
		printf("\nSelect %d stripe initator cells using the mouse pointer...\n\t[left-click to set white; right-click to set black]\n", maxinitiators);
	} else {
		printf("\nRandom stripes will be generated.\nR-D simulation running...\n");
		runningRD = 1;
	}
}


int main(int argc, char** argv)
{
	FILE *fp;
	
	/* print out usage information */
	if ((argc == 1) || (argv[1][0] == '-')) {
		printUsage();
		return 0;
	}
	
	/* get PLY filename from command-line */
	filename = argv[1];
	printf("%s\n",filename);
	
	/* test if a parameter file is given, otherwise use defaults */
	if (argv[2] != NULL) {
		switch (argv[2][1]) {
			case 'p':
				if (fp = fopen(argv[3],"r")) {
					fscanf(fp, "%f\n%f\n%f\n%f\n%f", &p1, &p2, &p3, &Dg, &Ds);
					printf("Taking parameters from %s:\n", argv[3]);
					fclose(fp);
					break;
				} else
					printf("Parameter file %s not found. ",argv[3]);
			default:
				printf("Using default parameters:\n");
		}
	} else {
		printf("Using default parameters:\n");
	}
	printf("\tp1: %f\n\tp2: %f\n\tp3: %f\n\tDg: %f\n\tDs: %f\n",p1,p2,p3,Dg,Ds);
	
	/* intialise GLUT and register callbacks, then enter main loop */
	glutInit(&argc, argv);
	glutInitDisplayMode(GLUT_DOUBLE | GLUT_DEPTH | GLUT_RGB);
	glutInitWindowSize(winx, winy);
	glutInitWindowPosition(300,200);
	win = glutCreateWindow("Stripey R-D Objects");
	init();
	glutDisplayFunc(display);
	glutReshapeFunc(reshape);
	glutIdleFunc(update);
	glutKeyboardFunc(keyb);
	glutSpecialFunc(special);
	glutMouseFunc(mouse);
	glutMainLoop();
	
	return 0;
}