/*	3D DEM data display with Landsat Overlay
	Program by B.J. Guillot
	Last updated May 5, 2004
*/

#include <stdio.h>
#include <GL/glut.h>

/* hardcoded data for exact pixel sizes of lab 8 data */
#define XDEM 1123
#define YDEM 940
#define NUM_SAMPLES XDEM*YDEM

#define XLANDSAT 1278
#define YLANDSAT 819
#define LANDSAT_SAMPLES XLANDSAT*YLANDSAT
#define NLANDSAT 6

/* global variables saving current state */
int DECIMATION=1;
short int buffer[NUM_SAMPLES] = {0};
unsigned char landsat[NLANDSAT][LANDSAT_SAMPLES] = {0};
float bandstring[6] = {0.47, 0.56, 0.66, 0.83, 1.65, 2.21};
int histLow[6] = {0};
int histHigh[6] = {0};
int maxHeight = 0;
int minHeight = 9999999;
int thetaX = 0;
int thetaY = 0;
int thetaZ = 90;
float translateX = 0;
float translateY = 0;
float translateZ = 0;
float scale = 0.5;
int flatMode = 1;
int bands[3] = {0};
int band = 0;
float stretchB[3];
float stretchM[3];

int loaddata(void) {
	FILE *demfile;
	int status;
	int i, j;

	printf("Loading DEM data...\n");
	demfile = fopen("dem.bil", "rb");
	if (!demfile) {
		printf("ERROR: Cannot open file!");
		return 1;
	}

	status = fread(buffer, 1, sizeof(buffer), demfile);
	fclose(demfile);

	printf("Loading LANDSAT data...\n");
	demfile = fopen("big_pine.bsq", "rb");
	if (!demfile) {
		printf("ERROR: Cannot open Landsat file!");
		return 1;
	}

	for (i = 0; i < NLANDSAT; i++) {
		int histogram[256] = {0};
		int count = 0;
		float percent = 0.0;

		status = fread(landsat[i], 1, LANDSAT_SAMPLES, demfile);
		for (j = 0; j < LANDSAT_SAMPLES; j++) /* generate histogram */
			histogram[landsat[i][j]]++;

		for (j = 0; j < 256; j++) {
			count += histogram[j];
			percent = ((float)count)/((float)LANDSAT_SAMPLES);
			if (percent >= 0.02) {
				printf("2 percent on low end at bin# %d\n", j);
				histLow[i] = j;
				break;
			}
		}

		count = 0;
		for (j = 255; j >= 0; j--) {
			count += histogram[j];
			percent = ((float)count)/((float)LANDSAT_SAMPLES);
			if (percent >= 0.02) {
				printf("2 percent on high end at bin# %d\n", j);
				histHigh[i] = j;
				break;
			}
		}


	}
	fclose(demfile);


	for (i = 0; i < NUM_SAMPLES; i++) {
		char b[2], temp; 	/* swap bytes 1 and 2 */
		memcpy(&b[0], &buffer[i], 2);
		temp = b[0];
		b[0] = b[1];
		b[1] = temp;
		memcpy(&buffer[i], &b[0], 2);
		if (buffer[i] > maxHeight)
			maxHeight = buffer[i];
		if ((buffer[i] > 0) && (buffer[i] < minHeight))
			minHeight = buffer[i];
	}

	for (i = 0; i < NUM_SAMPLES; i++)
		buffer[i] -= minHeight;
	maxHeight -= minHeight;

	return 0;
}

/* calculate point in 2D space */
void calcPoint(int i, int j, float *vertex) { 
	vertex[2] = -0.5;
	vertex[1] = (float)(i-XDEM/2) / XDEM;
	vertex[0] = (float)(j-XDEM/2) / XDEM;
}

/* calculate point in 3D space */
int calcVertex(int i, int j, float *vertex) { 
	int preZ = buffer[j*XDEM+i];
	if (preZ > 0) {
		float x, y, z;
		preZ -= maxHeight/2;
		if (flatMode > 0)
			vertex[2] = 0.0;
		else
			vertex[2] = -((float)preZ / (float)maxHeight) / 5.0; 
		/*bjg 5/4/2004 - add negative sign to vertex[2] to fix upside down orientation */
		vertex[1] = (float)(i-XDEM/2) / XDEM;
		vertex[0] = (float)(j-XDEM/2) / XDEM;
	} else
		return 0; //null vertez
	return 1; //valid vertex
}


/* compute color for point using linear stretch */
void stretch(int i, int j, unsigned char *color) { 
	float col[3];
	col[0] = stretchM[0]*(float)landsat[bands[0]][j*XLANDSAT+i] + stretchB[0]; 
	col[1] = stretchM[1]*(float)landsat[bands[1]][j*XLANDSAT+i] + stretchB[1];
	col[2] = stretchM[2]*(float)landsat[bands[2]][j*XLANDSAT+i] + stretchB[2];
	if (col[0] < 0.0)
		color[0] = 0;
	else if (col[0] > 255.0)
		color[0] = 255;
	else
		color[0] = col[0];

	if (col[1] < 0.0)
		color[1] = 0;
	else if (col[1] > 255.0)
		color[1] = 255;
	else
		color[1] = col[1];

	if (col[2] < 0.0)
		color[2] = 0;
	else if (col[2] > 255.0)
		color[2] = 255;
	else
		color[2] = col[2];

}

void draw() {
	int i, j;
	float v1[3], v2[3], v3[3], v4[3], n[3];
	int ni[XDEM], nj[XDEM]; //north
	unsigned char color[3];

	for (i = 0; i < XDEM; i++)
		ni[i] = -1;
 
	glBegin(GL_QUADS);

		for (i = 0; i < XDEM-DECIMATION; i+=DECIMATION) {
			for (j = 0; j < YDEM-DECIMATION; j+=DECIMATION) {
				if ((calcVertex(i, j, v1)) &&  /* verify vericies are not null */
						(calcVertex(i, j+DECIMATION, v2)) &&
						(calcVertex(i+DECIMATION, j+DECIMATION, v3)) &&
						(calcVertex(i+DECIMATION, j, v4))) {
					 
					if ((i < XLANDSAT) && (j < YLANDSAT)) {
						stretch(i, j, color);
						glColor3ubv(color);
					} else
						glColor4f(0, 0, 0, 0);
				
					glVertex3fv(v1); 
					glVertex3fv(v2);
					glVertex3fv(v3);
					glVertex3fv(v4);

				}
			}
		}
	glEnd();

}



void display(void)
{
	int i;

	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

	glLoadIdentity();    /* clear the matrix */
	glScalef(scale, scale, scale);
	glTranslatef(translateX, translateY, translateZ);
	glRotatef(thetaX, 1, 0, 0);
	glRotatef(thetaY, 0, 1, 0);
	glRotatef(thetaZ, 0, 0, 1);  

	printf("Bands:  R=%f  G=%f  B=%f  Rotation: X=%d  Y=%d  Z=%d\n", 
		bandstring[bands[0]], bandstring[bands[1]], bandstring[bands[2]],
		thetaX, thetaY, thetaZ);

	/* recompute histogram stretch parameters */
	for (i = 0; i < 3; i++) {
		stretchM[i] = 255.0/(float)((float)histHigh[bands[i]]-(float)histLow[bands[i]]);
		stretchB[i] = -histLow[bands[i]]*stretchM[i];
	}
  
	draw();

	glFlush();
	glutSwapBuffers(); /* use double buffers to elim flicker */
}

/* MAIN KEYBOARD EVENT LOOP */
void keyboard(unsigned char key, int x, int y)	
{
	switch (key) {
		case 27: case 'q' : case 'Q' : { exit(0); }
		case 'i' : { thetaX+=5; break; }
		case 'j' : { thetaY+=5; break; }
		case 'k' : { thetaZ+=5; break; }
		case 'I' : { thetaX-=5; break; }
		case 'J' : { thetaY-=5; break; }
		case 'K' : { thetaZ-=5; break; }
		case 'x' : { translateX+=1.0/scale; break; }
		case 'y' : { translateY+=1.0/scale; break; }
		case 'z' : { translateZ+=1.0/scale; break; }
		case 'X' : { translateX-=1.0/scale; break; }
		case 'Y' : { translateY-=1.0/scale; break; }
		case 'Z' : { translateZ-=1.0/scale; break; }
		case 'h' : { translateX=0; translateY=0; translateZ=0; break; } /*home*/
		case 's' : { scale += 1.0; break; }
		case 'S' : { scale -= 1.0; break; }
		case 'f' : { if (flatMode == 0) flatMode = 1; else flatMode = 0; break; }
		case '1' : { bands[band++] = 0; if (band > 2) band = 0; break; }
		case '2' : { bands[band++] = 1; if (band > 2) band = 0;break; }
		case '3' : { bands[band++] = 2; if (band > 2) band = 0; break; }
		case '4' : { bands[band++] = 3; if (band > 2) band = 0; break; }
		case '5' : { bands[band++] = 4; if (band > 2) band = 0; break; }
		case '7' : { bands[band++] = 5; if (band > 2) band = 0; break; }
		case 't' : { bands[0] = 2; bands[1] = 1; bands[2] = 0; break; } /* true color */
		case 'd' : { DECIMATION++; break; }
		case 'D' : { if (DECIMATION > 1) DECIMATION--; break; }
	}

	display();
}

 /* callback when window is re-sized */
void reshape(int w, int h) 
{
	glViewport(0, 0, w, h);
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();

	if (w <= h)
		glOrtho(-2.0, 2.0, -2.0*(GLfloat)h / (GLfloat)w,
			2.0*(GLfloat)h / (GLfloat)w, -20.0, 20.0);
	else
		glOrtho(-2.0*(GLfloat)w / (GLfloat) h,
			2.0*(GLfloat)w / (GLfloat)h, -2.0, 2.0, -20.0, 20.0);
	glMatrixMode(GL_MODELVIEW);
}

void glutinit(int argc, char **argv) {
	glutInit(&argc, argv);
	glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);
	glutInitWindowSize(500, 500);
	glutInitWindowPosition(100, 100);
	glutCreateWindow(argv[0]);
	glutKeyboardFunc(keyboard);
	glutDisplayFunc(display);
	glutReshapeFunc(reshape);

	glClearColor(0, 0, 0, 0);
	glShadeModel(GL_FLAT);
	glEnable(GL_DEPTH_TEST);

	glutMainLoop();
}

int main(int argc, char **argv) {
	loaddata();
	glutinit(argc, argv);
}