release / siproc

///KNOWN BUGS

///		*) The mouse pointer position gets messed up if any of the band windows are resized
///				This is caused by not knowing which window is clicked when the mouse_band() callback is triggered.
///				If it's possible to get the window, just change the if(n == 0) code to reference the correct window and pull it outside the loop
///		*) Each spectrum should be displayed in a different color
///		*) ENVI files with different numbers of bands aren't supported
///				This can be fixed by changing everything in the HSIview code to deal with wavenumbers rather than bands


#include <iostream>

#include <stim/parser/arguments.h>
#include <stim/envi/envi.h>
#include <stim/image/image.h>
#include <stim/visualization/colormap.h>
#include <stim/parser/filename.h>

#ifdef _WIN32
	#include <GL/glew.h>
#endif

#include <GL/glut.h>
#include <stim/gl/error.h>

#define NN	5

size_t N;										//number of open ENVI files

stim::arglist args;								//program arguments structure

stim::envi ENVI[NN];							//global ENVI file being visualized
size_t SAMPLES;									//dimensions of the loaded ENVI files
size_t LINES;
size_t BANDS[NN];								//number of bands in each spectrum
stim::image<unsigned char> band_image[NN];		//stores an image of the band
float* spectrum[NN];							//array stores the current spectrum
double* wavelengths[NN];						//array stores the independent variable for the spectrum (usually wavelength or frequency)
float sMax, sMin;								//maximum value in the current spectrum
float wMin, wMax;								//minimum and maximum independent spectral values

float colorMin, colorMax;
bool colorManual = false;						//control the color mapping manually

size_t B[NN] = {0, 0};				//current band being visualized
double W = 0.0;									//wavelength current visualized
size_t X, Y;						//current pixel position

double vFactor = 0.2;							//fraction of the spectral window used as space above the spectrum to display numbers

GLuint tex_id[NN];

int winBand[NN];								//Window identifier for the band image
int winSpectrum;								//window identifier for the spectrum

float colorMap[NN][3] = {{1, 1, 1}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}, {0, 1, 1}};

///Output software advertisement and instructions
void advertise(){
	//output advertisement
	std::cout<<std::endl<<std::endl;
	std::cout<<"========================================================================="<<std::endl;
	std::cout<<"Thank you for using the HSIview spectroscopic image visualization software!"<<std::endl;
	std::cout<<"Scalable Tissue Imaging and Modeling (STIM) Lab, University of Houston"<<std::endl;
	std::cout<<"Developers: David Mayerich"<<std::endl;
	std::cout<<"Source: https://github.com/stimlab/hsiproc.git"<<std::endl;
	std::cout << "This version has been compiled on " << __DATE__ << " at " << __TIME__ << std::endl;
	std::cout<<"========================================================================="<<std::endl<<std::endl;
	std::cout<<std::endl<<"usage: hsiview input --option [A B C ...]"<<std::endl;
	std::cout<<std::endl<<std::endl
			<< "examples:"<<std::endl
			<< "	hsiview envifile"<<std::endl
			<< "		Visualize the ENVI file 'envifile'"<<std::endl
			;
	std::cout<<"LEFT click on a pixel to display the associated spectrum"<<std::endl;
	std::cout<<"LEFT click on a position in the spectrum to display the associated band image"<<std::endl;
	std::cout<<"MIDDLE click on the spectrum to turn manual color-mapping on/off"<<std::endl;
	std::cout<<"\tSHIFT + LEFT click in the spectral window to set the minimum spectral value (blue)"<<std::endl;
	std::cout<<"\tSHIFT + RIGHT click in the spectral window to set the maximum spectral value (red)"<<std::endl;
	std::cout<<std::endl<<std::endl;

	std::cout<<args.str();
}

//sets an OpenGL viewport taking up the entire window
void glut_band_projection(){

	glMatrixMode(GL_PROJECTION);									//load the projection matrix for editing
	glLoadIdentity();												//start with the identity matrix
	int X = glutGet(GLUT_WINDOW_WIDTH);								//use the whole screen for rendering
	int Y = glutGet(GLUT_WINDOW_HEIGHT);
	glViewport(0, 0, X, Y);											//specify a viewport for the entire window
	float aspect = (float)X / (float)Y;								//calculate the aspect ratio
	gluOrtho2D((double)0, (double)SAMPLES, (double)LINES, (double)0);		//set up a perspective projection
}

/// Render an image of the current band
void render_band(size_t n){
	glut_band_projection();

	glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();

	glClearColor(0, 0, 0, 0);
	glClear(GL_COLOR_BUFFER_BIT);

	glEnable(GL_TEXTURE_2D);										//enable texture mapping
	glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE);		//texture map will be used as the network color
	glBindTexture(GL_TEXTURE_2D, tex_id[n]);							//bind the Brewer texture map

	glBegin(GL_QUADS);								//draw the band image
		glTexCoord2d(0, 0);
		glVertex2d(0, 0);
		glTexCoord2d(0, 1);
		glVertex2d(0, (double)LINES);
		glTexCoord2d(1, 1);
		glVertex2d((double)SAMPLES, (double)LINES);
		glTexCoord2d(1, 0);
		glVertex2d((double)SAMPLES, 0);
	glEnd();

	//draw crosshairs
	glDisable(GL_TEXTURE_2D);					//disable the texture map so that it isn't rendered in the crosshairs
	glColor3f(0.0f, 0.0f, 0.0f);					//the crosshairs will be black
	glBegin(GL_LINES);							//draw two lines
		glVertex2d(0, (double)Y + 0.5);
		glVertex2d((double)SAMPLES, (double)Y + 0.5);

		glVertex2d((double)X + 0.5, (double)0);
		glVertex2d((double)X + 0.5, (double)LINES);
	glEnd();									//finish drawing the crosshairs

	glutSwapBuffers();							//swap the doubles
}

void render_bands(){
	for(size_t n = 0; n < N; n++)
		render_band(n);
}

//Draw a number at the specified location
void draw_num(double num, double x, double y, void* font){
	glRasterPos2d(x, y);

	std::stringstream ss;											//create a string stream
	ss<<num;														//store the spectral value in the stream
	std::string str = ss.str();										//convert the stream to a string
	for(size_t i = 0; i < str.length(); i++)						//for each character in the string
		glutBitmapCharacter(font, str[i]);	//render that character
}

/// Draw a spectrum
void draw_spectrum(size_t n){

	//render the spectrum
	//glColor3f(1.0, 1.0, 1.0);
	glColor3f(colorMap[n][0], colorMap[n][1], colorMap[n][2]);
	glBegin(GL_LINE_STRIP);
		for(unsigned long long b = 0; b < BANDS[n]; b++)
			glVertex2d((double)wavelengths[n][b], (double)spectrum[n][b]);
	glEnd();

        glColor3f(0.0,1.0,0.0);
	draw_num((double)spectrum[n][B[n]], (double)W, (double)spectrum[n][B[n]], GLUT_BITMAP_TIMES_ROMAN_24);	//render the current spectral value
}

void render_spectra(){
	glutSetWindow(winSpectrum);						//set the current window to the spectrum window
	glClearColor(0, 0, 0, 0);						//set the clear color
	glClear(GL_COLOR_BUFFER_BIT);					//clear the window
	CHECK_OPENGL_ERROR
	glMatrixMode(GL_PROJECTION);					//load the projection matrix for editing
	glLoadIdentity();								//start with the identity matrix
	int X = glutGet(GLUT_WINDOW_WIDTH);				//use the whole screen for rendering
	int Y = glutGet(GLUT_WINDOW_HEIGHT);
	glViewport(0, 0, X, Y);							//specify a viewport for the entire window
	CHECK_OPENGL_ERROR
	if(sMin != sMax)
		gluOrtho2D(wMin, wMax, (double)sMin, (double)sMax + (sMax - sMin) * vFactor);		//set up a perspective projection
	else
		gluOrtho2D(wMin, wMax, (double)sMin - 1.0, (double)sMax + 1.0);
	CHECK_OPENGL_ERROR
	//render the zero line
	glColor3f(0.0, 1.0, 1.0);
	glBegin(GL_LINES);
		glVertex2d(wMin, 0);
		glVertex2d(wMax, 0);
	glEnd();

	for(size_t n = 0; n < N; n++)
		draw_spectrum(n);

	if(colorManual){								//if manual color mapping is applied
		glColor3f(1.0, 0.0, 0.0);					//draw a red max line
		glBegin(GL_LINES);
			glVertex2d(wMin, colorMax);
			glVertex2d(wMax, colorMax);
		glEnd();
		glColor3f(0.0, 0.0, 1.0);					//draw a blue min line
		glBegin(GL_LINES);
			glVertex2d(wMin, colorMin);
			glVertex2d(wMax, colorMin);
		glEnd();
	}
        CHECK_OPENGL_ERROR

	//draw the band line
	glColor3f(0.0, 1.0, 0.0);						//the band line will be green
	glBegin(GL_LINES);								//draw a single band line
		glVertex2d((double)W, (double)sMin);
		glVertex2d((double)W, (double)sMax);
	glEnd();

	double pixelsize = (sMax - sMin) / glutGet(GLUT_WINDOW_HEIGHT);
        //std::cout<<pixelsize<<std::endl;
	draw_num(sMin, wMin, sMin, GLUT_BITMAP_TIMES_ROMAN_24);						//render the minimum value
	draw_num(sMax, wMin, sMax, GLUT_BITMAP_TIMES_ROMAN_24);		//render the maximum value

	glutSwapBuffers();								//swap buffers
}

void update_spectrum_minmax(){
	//update the global min and max values
	wMin = wMax = (float)wavelengths[0][0];							//initialize the min and max independend variables
	sMax = sMin = spectrum[0][0];
	for(size_t i = 0; i < N; i++){								//for each spectrum
		//THIS STRANGE SYNTAX FOR std::min AND std::max are to get around pre-defined windows macros :-(
		wMin = (std::min)(wMin, (float)wavelengths[i][0]);
		wMax = (std::max)(wMax, (float)wavelengths[i][BANDS[i]-1]);

		for(size_t b = 0; b < BANDS[i]; b++){					//find the minimum and maximum spectral values
			if(spectrum[i][b] > sMax) sMax = spectrum[i][b];
			if(spectrum[i][b] < sMin) sMin = spectrum[i][b];
		}
	}
}
void update_spectrum(size_t n){

	if(spectrum[n] == NULL)												//if space hasn't been allocated for the spectrum
		spectrum[n] = (float*) malloc(BANDS[n] * sizeof(float));		//allocate the space
	if(wavelengths[n] == NULL)
		wavelengths[n] = (double*) malloc(BANDS[n] * sizeof(double));
	if(ENVI[n].header.wavelength.size() == 0){							//if there are no wavelength values for the ENVI file
		for(size_t i = 0; i < BANDS[n]; i++)							//for each band
			wavelengths[n][i] = (double)i;								//save the band number as the wavelength
	}
	else
		memcpy(wavelengths[n], &ENVI[n].header.wavelength[0], sizeof(double) * BANDS[n]);	//copy the wavelength values

	ENVI[n].spectrum(spectrum[n], X, Y);						//load the spectrum from the ENVI file

	glutSetWindow(winSpectrum);
	glutPostRedisplay();
}

void update_spectra(){
	for(size_t n = 0; n < N; n++)
		update_spectrum(n);
	update_spectrum_minmax();
}

/// Update the band image displayed for ENVI file n
void update_band(size_t n){
	unsigned long long sx = SAMPLES;				//calculate the image size
	unsigned long long sy = LINES;

	float* buffer = NULL;												//allocate a buffer
	buffer = (float*)malloc( sx * sy * sizeof(float));					//allocate space for the band
	std::vector<size_t> b_idx = ENVI[n].header.band_index(W);			//get the bands closest to the current wavelength
	if (b_idx.size() > 0) {												//if a band exists
		B[n] = b_idx[0];
		if (b_idx.size() > 1) {
			double w0 = ENVI[n].header.wavelength[b_idx[0]];				//get the lower wavelength
			double w1 = ENVI[n].header.wavelength[b_idx[1]];				//get the upper wavelength
			if (abs(w1 - W) < abs(W - w0)) B[n] = b_idx[1];							//if the current wavelength is closer to the higher band, display the higher one
			std::cout << "w0: " << w0 << std::endl;
			std::cout << "w1: " << w1 << std::endl;
		}
		std::cout << "Display band: " << B[n] << std::endl;
	}
	ENVI[n].band_index(buffer, B[n]);								//retrieve the band image
	if(colorManual)
		stim::cpu2cpu<float>(buffer, band_image[n].data(), sx * sy, colorMin, colorMax, stim::cmBrewer);
	else
		stim::cpu2cpu<float>(buffer, band_image[n].data(), sx * sy, stim::cmBrewer);		//convert the band image to a color image
	free(buffer);

	glutSetWindow(winBand[n]);							//switch to the band window
	glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);		//enable linear interpolation
	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP);			//clamp the values at the minimum and maximum
	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);

	//test to see if the texture can be accomodated
	glTexImage2D(GL_PROXY_TEXTURE_2D, 0, GL_RGB, (GLsizei)sx, (GLsizei)sy, 0, GL_RGB, GL_UNSIGNED_BYTE, band_image[n].data());
	int success;
	glGetTexLevelParameteriv(GL_PROXY_TEXTURE_2D, 0, GL_TEXTURE_WIDTH, &success);
	if(success)
		glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, (GLsizei)sx, (GLsizei)sy, 0, GL_RGB, GL_UNSIGNED_BYTE, band_image[n].data());	//upload the texture map to the GPU
	else
		std::cout << "Texture format not supported (likely the image size won't work with OpenGL - which sucks)" << std::endl;
	CHECK_OPENGL_ERROR

	glutPostRedisplay();							//update the visualization
}

/// Updates all band images
void update_bands(){
	for(size_t n = 0; n < N; n++)
		update_band(n);
}

/// Refresh band windows (update overlays and pixel locator positions)
void refresh_band(size_t n){
	glutSetWindow(winBand[n]);
	glutPostRedisplay();
}

void refresh_bands(){
	for(size_t n = 0; n < N; n++)
		refresh_band(n);
}

/// Updates the spectrum when a new pixel position is selected in any band window
void mouse_band(int button, int state, int x, int y){

	for(size_t n = 0; n < N; n++){														//for each band window
		//translate from the window to the image
		glutSetWindow(winBand[n]);													//set the current window to the band window
		if(n == 0){																	//only evaluate the mouse position for the first window
			X = (unsigned long long)(SAMPLES * (double)x / (double)glutGet(GLUT_WINDOW_WIDTH));	//calculate the ENVI x coordinate
			Y = (unsigned long long)(LINES * (double)y / (double)glutGet(GLUT_WINDOW_HEIGHT));	//calculate the ENVI y coordinate
			std::cout<<"("<<X<<", "<<Y<<")"<<std::endl;
		}

		//update_spectrum(n);								//update the current spectrum
		refresh_band(n);
	}
	update_spectra();
}

/// Special keys for the band image
void special_band(int key, int x, int y){
	switch(key){
	case GLUT_KEY_UP:
		Y--;
		if(Y > LINES-1) Y = LINES-1;
		break;
	case GLUT_KEY_DOWN:
		Y++;
		if(Y > LINES-1) Y = 0;
		break;
	case GLUT_KEY_LEFT:
		X--;
		if(X > SAMPLES-1) X = SAMPLES-1;
		break;
	case GLUT_KEY_RIGHT:
		X++;
		if(X > SAMPLES-1) X = 0;
		break;
	}
	refresh_bands();
	update_spectra();
}

/// Special keys for the band image
void key_spectrum(unsigned char key, int x, int y){
	
	std::ofstream outfile("spectrum.csv");
	if(key == 's'){
		for(size_t b = 0; b < BANDS[0]; b++){
			if(b != 0) outfile<<",";
			outfile<<spectrum[0][b];
		}
	}

}

//find the current wavelength to be displayed
void mouse_spectrum(int button, int state, int x, int y){

	//set manual colormap bounds
	if(glutGetModifiers() == GLUT_ACTIVE_SHIFT){
		glutSetWindow(winSpectrum);
		double a = (double)(glutGet(GLUT_WINDOW_HEIGHT) - y) / (double)glutGet(GLUT_WINDOW_HEIGHT);	//calculate the click position in [0 1]
		double aMax = sMax + (sMax - sMin) * vFactor;
		double A = a * (aMax - sMin) + sMin;
		if(button == GLUT_LEFT_BUTTON)
			colorMin = (float)A;
		if(button == GLUT_RIGHT_BUTTON)
			colorMax = (float)A;
	}

	//band selection (left button, unmodified)
	else if(button == GLUT_LEFT_BUTTON){						//left button changes bands
		glutSetWindow(winSpectrum);									//retrieve data from the spectral window
		W = (double)x / (double)glutGet(GLUT_WINDOW_WIDTH) * (wMax - wMin) + wMin;		//calculate the wavelength to be visualized
		std::cout<<"wavelength: "<<W<<std::endl;
	}
	
	//turn manual color mapping on/off
	else if(button == GLUT_MIDDLE_BUTTON && state == GLUT_DOWN)
			colorManual = !colorManual;							//swap the colorManual flag, turning manual color mapping on or off

	update_bands();
        update_spectra();
}

void glut_motion(int x, int y){
}

void texture_initialize(size_t n){
	unsigned long long sx = SAMPLES;				//calculate the image size
	unsigned long long sy = LINES;

	band_image[n] = stim::image<unsigned char>(sx, sy, 3);	//allocate space for the band visualization image

	glGenTextures(1, &tex_id[n]);								//generate a texture map name
	glBindTexture(GL_TEXTURE_2D, tex_id[n]);					//bind the texture map
}

void glut_initialize(std::string filename){

	int myargc = 1;					//GLUT requires arguments, so create some bogus ones
	char* myargv[1];
	myargv[0] = (char*)malloc(1);
	myargv[0][0] = 'h';

	glutInit(&myargc, myargv);									//pass bogus arguments to glutInit()
	glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGBA);				//generate a color buffer, depth buffer, and enable double buffering

	for(size_t n = 0; n < N; n++){
		glutInitWindowPosition(100,100);							//set the initial window position
		glutInitWindowSize(320,320);								//set the initial window size
		winBand[n] = glutCreateWindow((filename + " - HSIview - STIM Lab, UH").c_str());					//set the dialog box title
		glutDisplayFunc(render_bands);								//function executed for rendering - renders networks
		glutMouseFunc(mouse_band);									//executed on a mouse click - sets starting mouse positions for rotations
		glutMotionFunc(glut_motion);								//executed when the mouse is moved while a button is pressed
		glutSpecialFunc(special_band);								//special keys pressed when the band image is selected
		texture_initialize(n);										//set up texture mapping (create texture maps, enable features)
	}

	glutInitWindowPosition(400,500);							//set the initial window position
	glutInitWindowSize(960,320);								//set the initial window size
	winSpectrum = glutCreateWindow("Spectrum");
	glutMouseFunc(mouse_spectrum);
	glutDisplayFunc(render_spectra);
	glutKeyboardFunc(key_spectrum);

#ifdef _WIN32
	GLenum err = glewInit();									//initialize GLEW (necessary for Windows)
	if (GLEW_OK != err){										//eror with GLEW
	  std::cout<<"Error with GLEW: "<<glewGetErrorString(err)<<std::endl;
	  exit(1);
	}
#endif

}

int main(int argc, char* argv[]){

	args.add("help", "prints this help text");
	args.add("agilent", "parse the input file(s) as Agilent binary files");

	args.parse(argc, argv);										//parse the input arguments

	//get the number of input files to display
	N = args.nargs();

	if(N > NN){
		std::cout<<"ERROR - HSIview cannot currently support more than "<<NN<<" simultaneous files."<<std::endl;
		exit(1);
	}

	if(args["help"].is_set() || N == 0){
		advertise();
		exit(1);
	}



	for(size_t n = 0; n < N; n++){
		stim::filename fname = args.arg(n);							//get the file name string and store it as a stim::filename
		std::string envifile;
		std::string enviheader;
		if(args["agilent"].is_set())
			ENVI[n].open_agilent(args.arg(n));
		else if(fname.extension() == "drd" ||						//Agilent binary files: DMD, DRD, SEQ, DAT
				fname.extension() == "dmd" ||
				fname.extension() == "seq" ||
				fname.extension() == "dat")
			ENVI[n].open_agilent(fname.str());
		else if (fname.extension() == "hdr") {
			envifile = fname.extension("").str();
			enviheader = fname.str();
		}
		else {
			envifile = args.arg(n);
			enviheader = args.arg(n) + ".hdr";
		}
		if (!ENVI[n].open(envifile, enviheader, stim::io_in)) {
			std::cout << "ERROR opening ENVI file " << n << ": " << args.arg(n) << std::endl;
			exit(1);
		}
		BANDS[n] = ENVI[n].header.bands;
	}

	SAMPLES = ENVI[0].header.samples;
	LINES = ENVI[0].header.lines;

	

	glut_initialize(args.arg(0));											//create a GLUT window

	//make sure that the image size is supported
	GLint max_tex;
	glGetIntegerv(GL_MAX_TEXTURE_SIZE, &max_tex);
	CHECK_OPENGL_ERROR
	if(SAMPLES > max_tex || LINES > max_tex){
		std::cout<<"ERROR: image too large for OpenGL texture support. Maximum texture size is: "<<max_tex<<std::endl;
		exit(1);
	}


	W = (wMax - wMin)/2.0;
	update_bands();												//set the initial band image to the middle band

	X = SAMPLES/2;												//initialize to the center pixel
	Y = LINES/2;
	update_spectra();											//update the pixel position and spectrum
	glutMainLoop();
	ENVI[0].close();
}