segmentation / netmets

#include <stdlib.h>
#include <string>
#include <fstream>
#include <algorithm> 

//OpenGL includes
#include <GL/glut.h>
#include <GL/freeglut.h>

//STIM includes
#include <stim/visualization/gl_network.h>
#include <stim/biomodels/network.h>
#include <stim/visualization/gl_aaboundingbox.h>
#include <stim/parser/arguments.h>
#include <stim/visualization/camera.h>

#ifdef __CUDACC__
//CUDA includes
#include <cuda.h>
#endif

//BOOST includes
#include <boost/tuple/tuple.hpp>

//visualization objects
stim::gl_aaboundingbox<float> bb;			//axis-aligned bounding box object
stim::camera cam;					//camera object

// number of networks
unsigned num_nets = 0;

// networks
stim::gl_network<float> GT;			//ground truth network
stim::gl_network<float> T;			//test network
stim::gl_network<float> _GT;		//splitted GT
stim::gl_network<float> _T;			//splitted T

// indicator
unsigned ind = 0;						//indicator of mapping

// relationships
std::vector<unsigned> _gt_t;								// store indices of nearest edge points in _T for _GT
std::vector<unsigned> _t_gt;								// store indices of nearest edge points in _GT for _T

//hard-coded parameters
float resample_rate = 0.5f;			//sample rate for the network (fraction of sigma used as the maximum sample rate)
float camera_factor = 1.2f;			//start point of the camera as a function of X and Y size
float orbit_factor = 0.01f;			//degrees per pixel used to orbit the camera
float zoom_factor = 10.0f;
float radius_factor = 0.5f;

//mouse click
bool LButtonDown = false;			// true when left button down
bool RButtonDown = false;

//mouse position tracking
int mouse_x;
int mouse_y;

// render modes
bool compareMode = true;			// default mode is compare mode
bool mappingMode = false;

// random color set
std::vector<float> colormap;

// special key indicator
int mods;

//OpenGL objects
GLuint cmap_tex = 0;				//texture name for the color map

float delta;
float sigma = 3;					//default sigma
float radius = 3;					//equals to radius
int adjoint_fac = 0;
int light_fac = 0;

//sets an OpenGL viewport taking up the entire window
void glut_render_single_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
	gluPerspective(60, aspect, 0.1, 1000000);		//set up a perspective projection
}

//sets an OpenGL viewport taking up the left half of the window
void glut_render_left_projection(){

	glMatrixMode(GL_PROJECTION);					//load the projection matrix for editing
	glLoadIdentity();								//start with the identity matrix
	int X = glutGet(GLUT_WINDOW_WIDTH) / 2;			//only use half of the screen for the viewport
	int Y = glutGet(GLUT_WINDOW_HEIGHT);
	glViewport(0, 0, X, Y);							//specify the viewport on the left
	float aspect = (float)X / (float)Y;				//calculate the aspect ratio
	gluPerspective(60, aspect, 0.1, 1000000);		//set up a perspective projection
}

//sets an OpenGL viewport taking up the right half of the window
void glut_render_right_projection(){

	glMatrixMode(GL_PROJECTION);					//load the projection matrix for editing
	glLoadIdentity();								//start with the identity matrix
	int X = glutGet(GLUT_WINDOW_WIDTH) / 2;			//only use half of the screen for the viewport
	int Y = glutGet(GLUT_WINDOW_HEIGHT);
	glViewport(X, 0, X, Y);							//specify the viewport on the right
	float aspect = (float)X / (float)Y;				//calculate the aspect ratio
	gluPerspective(60, aspect, 0.1, 1000000);		//set up a perspective projection
}

void glut_render_modelview(){

	glMatrixMode(GL_MODELVIEW);						//load the modelview matrix for editing
	glLoadIdentity();								//start with the identity matrix
	stim::vec3<float> eye = cam.getPosition();		//get the camera position (eye point)
	stim::vec3<float> focus = cam.getLookAt();		//get the camera focal point
	stim::vec3<float> up = cam.getUp();				//get the camera "up" orientation

	gluLookAt(eye[0], eye[1], eye[2], focus[0], focus[1], focus[2], up[0], up[1], up[2]);	//set up the OpenGL camera
}

//draws the network(s)
void glut_render(void) {

	stim::vec3<float> eye = cam.getPosition();		//get the camera position (eye point)
	stim::vec3<float> focus = cam.getLookAt();		//get the camera focal point
	stim::vec3<float> up = cam.getUp();				//get the camera "up" orientation

	stim::vec3<float> v = focus - eye;				//calculate a vector from the camera to the center
	stim::vec3<float> v1 = v.cross(up);				
	stim::vec3<float> v2 = up.cross(v);
	v1 = v1.norm();
	v2 = v2.norm();
	stim::vec3<float> l1 = v1 + up;
	stim::vec3<float> l2 = v2 + -up;
	stim::vec3<float> s = bb.size();
	float r = s.len();

	//stim::vec3<float> p1 = focus + l1*r/5;
	//stim::vec3<float> p2 = focus + l2*r/5;
	stim::vec3<float> p1 = focus + cam.getUp() * 100000;
	stim::vec3<float> p2 = cam.getPosition();

	// light	
	GLfloat global_ambient[] = { 0.5, 0.5, 0.5, 1.0 };
	GLfloat ambient[] = { 0.0, 0.0, 0.0, 1.0 };
	GLfloat diffuse1[] = { 1.0, 1.0, 1.0, 1.0 };
	GLfloat diffuse2[] = { 0.4, 0.4, 0.4, 1.0 };
	GLfloat specular[] = { 1.0, 1.0, 1.0, 1.0 };
	GLfloat position1[] = { p1[0], p1[1], p1[2], 1.0 };		// upper right light source
	GLfloat position2[] = { p2[0], p2[1], p2[2], 1.0 };		// lower left light source

	glClearColor(0.0, 0.0, 0.0, 1.0);
	glShadeModel(GL_SMOOTH);

	glLightModelfv(GL_LIGHT_MODEL_AMBIENT, global_ambient);

	glLightfv(GL_LIGHT0, GL_AMBIENT, ambient);
	glLightfv(GL_LIGHT0, GL_DIFFUSE, diffuse1);
	glLightfv(GL_LIGHT0, GL_SPECULAR, specular);
	glLightfv(GL_LIGHT0, GL_POSITION, position1);

	glLightfv(GL_LIGHT1, GL_AMBIENT, ambient);
	glLightfv(GL_LIGHT1, GL_DIFFUSE, diffuse2);
	glLightfv(GL_LIGHT1, GL_SPECULAR, specular);
	glLightfv(GL_LIGHT1, GL_POSITION, position2);

	//no mapping, just comparing
	if (ind == 0) {	
		if (num_nets == 1) {										//if a single network is loaded
			glEnable(GL_DEPTH_TEST);								//enable depth
			glut_render_single_projection();						//fill the entire viewport
			glut_render_modelview();								//set up the modelview matrix with camera details
			glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);		//clear the screen
			GT.glCenterline0();										//render the GT network (the only one loaded)
			glDisable(GL_DEPTH_TEST);
		}

		if (num_nets == 2) {										//if two networks are loaded	
			glEnable(GL_TEXTURE_1D);										//enable texture mapping
			if (light_fac == 0)
				glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE);	//texture map will be used as the network color
			else
				glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_MODULATE);
			glBindTexture(GL_TEXTURE_1D, cmap_tex);							//bind the Brewer texture map
			
			glEnable(GL_DEPTH_TEST);								//enable depth
			glut_render_left_projection();							//set up a projection for the left half of the window
			glut_render_modelview();								//set up the modelview matrix using camera details
			glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);		//clear the screen

			GT.glCylinder(sigma, radius);							//render the GT network
			if (adjoint_fac == 1) {
				glDisable(GL_TEXTURE_1D);							//disable texture in order to render in other color
				glEnable(GL_BLEND);									//enable color blend
				glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);	//set blend function
				glDisable(GL_DEPTH_TEST);							//should disable depth to render transparancy
				glColor4f(0.0f, 0.3f, 0.0f, 0.2f);
				T.glAdjointCylinder(sigma, radius);
				glDisable(GL_BLEND);
				glEnable(GL_DEPTH_TEST);
				glEnable(GL_TEXTURE_1D);
				glColor4f(1.0f, 1.0f, 1.0f, 1.0f);
			}

			glut_render_right_projection();							//set up a projection for the right half of the window
			glut_render_modelview();								//set up the modelview matrix using camera details

			T.glCylinder(sigma, radius);							//render the T network

			sigma = radius;											// set sigma equal to radius
		}
	}

	//do comparing and mapping
	else {	
		if (num_nets == 1) {											//if a single network is loaded
			std::cout << "You should have at least two networks to do mapping." << std::endl;	//exit program because there isn't enough network
			exit(1);
		}
		if (num_nets == 2) {											//if two networks are loaded
			if (compareMode) {
				glEnable(GL_TEXTURE_1D);										//enable texture mapping
				if (light_fac == 0)
					glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE);	//texture map will be used as the network color
				else
					glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_MODULATE);//map light to texture
				glBindTexture(GL_TEXTURE_1D, cmap_tex);							//bind the Brewer texture map
				
				glEnable(GL_DEPTH_TEST);								//enable depth
				glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);		//clear the screen
				glut_render_left_projection();							//set up a projection for the left half of the window
				glut_render_modelview();								//set up the modelview matrix using camera details

				glPushMatrix();
				glTranslatef(p1[0], p1[1], p1[2]);
				glColor3f(0.0, 1.0, 0.0);
				glutSolidSphere(10, 10, 10);
				glPopMatrix();
				glPushMatrix();
				glTranslatef(p2[0], p2[1], p2[2]);
				glColor3f(0.0, 0.0, 1.0);
				//glutSolidSphere(10, 10, 10);
				glPopMatrix();
				glColor3f(1, 1, 1);

				_GT.glCylinder(sigma, radius);							//render the GT network

				if (adjoint_fac == 1) {
					glDisable(GL_TEXTURE_1D);							//temporarily disable texture
					glEnable(GL_BLEND);									//enable color blend
					glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);	//set blend function
					glDisable(GL_DEPTH_TEST);							//should disable depth
					glColor4f(0.0f, 0.3f, 0.0f, 0.2f);
					_T.glAdjointCylinder(sigma, radius);
					glDisable(GL_BLEND);
					glEnable(GL_DEPTH_TEST);
					glEnable(GL_TEXTURE_1D);							//re-enable texture
					glColor4f(1.0f, 1.0f, 1.0f, 1.0f);
				}

				glut_render_right_projection();							//set up a projection for the right half of the window
				glut_render_modelview();								//set up the modelview matrix using camera details

				_T.glCylinder(sigma, radius);							//render the T network

				sigma = radius;											//set sigma equal to radius
			}
			else {
				glEnable(GL_COLOR_MATERIAL);
				glEnable(GL_DEPTH_TEST);								//enable depth
				glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);		//clear the screen
				glut_render_left_projection();							//set up a projection for the left half of the window
				glut_render_modelview();								//set up the modelview matrix using camera details
						
				_GT.glRandColorCylinder(0, _gt_t, colormap, sigma, radius);

				glut_render_right_projection();							//set up a projection for the right half of the window
				glut_render_modelview();								//set up the modelview matrix using camera details
																		
				_T.glRandColorCylinder(1, _t_gt, colormap, sigma, radius);

				sigma = radius;											//set sigma equal to radius
			}
		}
	}
	glDisable(GL_DEPTH_TEST);
	glDisable(GL_TEXTURE_1D);

	if (num_nets == 2) {											// works only with two networks
		std::ostringstream ss;
		if (mappingMode)											// if it is in mapping mode
			ss << "Mapping Mode";
		else
			ss << "Compare Mode";									// default mode is compare mode

		if (light_fac == 1)
			glDisable(GL_LIGHTING);
		glMatrixMode(GL_PROJECTION);								// set up the 2d viewport for mode text printing
		glPushMatrix();
		glLoadIdentity();
		int X = glutGet(GLUT_WINDOW_WIDTH);							// get the current window width
		int Y = glutGet(GLUT_WINDOW_HEIGHT);						// get the current window height
		glViewport(0, 0, X / 2, Y);									// locate to left bottom corner
		gluOrtho2D(0, X, 0, Y);										// define othogonal aspect
		glColor3f(0.8, 0.0, 0.0);									// using red to show mode

		glMatrixMode(GL_MODELVIEW);
		glPushMatrix();
		glLoadIdentity();

		glRasterPos2f(0, 5);										//print text in the left bottom corner
		glutBitmapString(GLUT_BITMAP_TIMES_ROMAN_24, (const unsigned char*)(ss.str().c_str()));

		glPopMatrix();
		glMatrixMode(GL_PROJECTION);
		glPopMatrix();
		glColor3f(1.0, 1.0, 1.0);									//clear red color
		if (light_fac == 1)
			glEnable(GL_LIGHTING);
	}

	glDisable(GL_COLOR_MATERIAL);

	glutSwapBuffers();
}

// defines camera motion based on mouse dragging
void glut_motion(int x, int y){
	
	if(LButtonDown == true && RButtonDown == false && mods != GLUT_ACTIVE_CTRL){

	float theta = orbit_factor * (mouse_x - x);		//determine the number of degrees along the x-axis to rotate
	float phi = orbit_factor * (y - mouse_y);		//number of degrees along the y-axis to rotate

	cam.OrbitFocus(theta, phi);						//rotate the camera around the focal point

	mouse_x = x;									//update the mouse position
	mouse_y = y;
		
	glutPostRedisplay();							//re-draw the visualization
	}
}

// sets the menu options
void glut_menu(int value) {
	if (value == 1) {								// menu 1 represents comparing mode
		compareMode = true;
		mappingMode = false;
	}
	if (value == 2) {								// menu 2 represents mapping mode
		compareMode = false;
		mappingMode = true;
	}
	if (value == 3) {
		exit(0);
	}
	glutPostRedisplay();
}

// sets the mouse position when clicked
void glut_mouse(int button, int state, int x, int y){
	
	if(button == GLUT_LEFT_BUTTON && state == GLUT_DOWN){
		mouse_x = x;
		mouse_y = y;
		LButtonDown = true;
	}
	else if(button == GLUT_RIGHT_BUTTON && state == GLUT_DOWN){
		mouse_x = x;
		mouse_y = y;
		RButtonDown = true;
	}
	else if(button == GLUT_LEFT_BUTTON && state == GLUT_UP){
		mouse_x = x;
		mouse_y = y;
		LButtonDown = false;
	}
	else if(button == GLUT_RIGHT_BUTTON && state == GLUT_UP){
		mouse_x = x;
		mouse_y = y;
		RButtonDown = false;
	}

	/// implementation of mouse click mapping feedback
	mods = glutGetModifiers();											// get modifier keys
	if (mods == GLUT_ACTIVE_CTRL) 										// if the CTRL key is pressed
		if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN) {
			std::cout << "( " << x << ", " << y << " )" << std::endl;	// if the CTRL key is pressed and LEFT BUTTON is DOWN, print the window coordinates

			GLint    viewport[4];
			GLdouble modelview[16];
			GLdouble projection[16];
			GLdouble winX, winY, winZ;
			GLdouble posX, posY, posZ;

			glGetIntegerv(GL_VIEWPORT, viewport);
			glGetDoublev(GL_MODELVIEW_MATRIX, modelview);
			glGetDoublev(GL_PROJECTION_MATRIX, projection);

			winX = (GLdouble)x;
			winY = viewport[3] - (GLdouble)y;
			glReadPixels((GLint)winX, (GLint)winY, (GLsizei)1, (GLsizei)1, GL_DEPTH_COMPONENT, GL_FLOAT, &winZ);		// need frame buffer FBO
			gluUnProject(winX, winY, winZ, modelview, projection, viewport, &posX, &posY, &posZ);						// not sure why it should add 1 to the winZ
			
			std::cout << "( " << posX << ", " << posY << ", "<< posZ <<" )" << std::endl;
		}
}

// define camera move based on mouse wheel move(actually we can combine this with glut_mouse)
void glut_wheel(int wheel, int direction, int x, int y) {

	if (direction > 0)								// if it is button 3(up), move closer
		delta = zoom_factor;
	else											// if it is button 4(down), leave farther
		delta = -zoom_factor;

	cam.Push(delta);
	glutPostRedisplay();
}

// define keyboard inputs
void glut_keyboard(unsigned char key, int x, int y){
	
	// register different keyboard operation
	switch (key) {
		
		// change render mode
		case 'm':															// if keyboard 'm' is pressed, then change render mode
			if (compareMode && !mappingMode && ind && !adjoint_fac) {		// if current mode is comparing mode
				compareMode = false;
				mappingMode = true;
			}
			else if (!compareMode && mappingMode && ind && !adjoint_fac) {	// if current mode is mapping mode
				compareMode = true;
				mappingMode = false;
			}
			break;

		// zooming
		case 'w':						// if keyboard 'w' is pressed, then move closer
			delta = zoom_factor;
			cam.Push(delta);
			break;
		case 's':						// if keyboard 's' is pressed, then leave farther
			delta = -zoom_factor;
			cam.Push(delta);
			break;

		// resample and re-render the cylinder in different radius
		case 'd':						// if keyboard 'd' is pressed, then increase radius by radius_factor
			radius += radius_factor;
			break;
		case 'a':						// if keyboard 'a' is pressed, then decrease radius by radius_factor
			radius -= radius_factor;
			// get rid of the degenerated case when radius decrease below 0
			if (radius < 0.001f)
				radius = 0.2;
			break;

		// turn on/off the light
		case 'l':						// if keyboard 'l' is pressed, then change the light
			if (!light_fac && !adjoint_fac) {
				light_fac = 1;
				glEnable(GL_LIGHTING);
				glEnable(GL_LIGHT0);
				glEnable(GL_LIGHT1);
			}
			else if (light_fac && !adjoint_fac) {
				light_fac = 0;
				glDisable(GL_LIGHTING);
				glDisable(GL_LIGHT0);
				glDisable(GL_LIGHT1);
			}
			break;

		// render a transparant T very close to GT in compare mode
		case 32:						// if keyboard 'SPACE' is pressed, then change the adjoint_fac
			if (!adjoint_fac && compareMode && !light_fac)
				adjoint_fac = 1;
			else if(adjoint_fac && compareMode && !light_fac)
				adjoint_fac = 0;
			break;

		// close window and exit application
		case 27:						// if keyboard 'ESC' is pressed, then exit
			exit(0);
	}
	glutPostRedisplay();
}

#define BREWER_CTRL_PTS 11							//number of control points in the Brewer map
void texture_initialize(){

	//define the colormap
	static float  brewer_map[BREWER_CTRL_PTS][3] = {			//generate a Brewer color map (blue to red)
		{0.192157f, 0.211765f, 0.584314f},
		{0.270588f, 0.458824f, 0.705882f},
		{0.454902f, 0.678431f, 0.819608f},
		{0.670588f, 0.85098f, 0.913725f},
		{0.878431f, 0.952941f, 0.972549f},
		{1.0f, 1.0f, 0.74902f},
		{0.996078f, 0.878431f, 0.564706f},
		{0.992157f, 0.682353f, 0.380392f},
		{0.956863f, 0.427451f, 0.262745f},
		{0.843137f, 0.188235f, 0.152941f},
		{0.647059f, 0.0f, 0.14902f}
	};

	glGenTextures(1, &cmap_tex);								//generate a texture map name
	glBindTexture(GL_TEXTURE_1D, cmap_tex);						//bind the texture map

	glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);		//enable linear interpolation
	glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
	glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_WRAP_S, GL_CLAMP);			//clamp the values at the minimum and maximum
	glTexImage1D(GL_TEXTURE_1D, 0, 3, BREWER_CTRL_PTS, 0, GL_RGB, GL_FLOAT,	//upload the texture map to the GPU
					brewer_map);
}

//Initialize the OpenGL (GLUT) window, including starting resolution, callbacks, texture maps, and camera
void glut_initialize(){
	
	int myargc = 1;					//GLUT requires arguments, so create some bogus ones
	char* myargv[1];
	myargv [0]=strdup ("netmets");

	glutInit(&myargc, myargv);									//pass bogus arguments to glutInit()
	glutSetOption(GLUT_MULTISAMPLE, 8);
	glutInitDisplayMode(GLUT_DEPTH | GLUT_DOUBLE | GLUT_RGBA);	//generate a color buffer, depth buffer, and enable double buffering
	glutInitWindowPosition(100,100);							//set the initial window position
	glutInitWindowSize(320, 320);								//set the initial window size
	glutCreateWindow("NetMets - STIM Lab, UH");					//set the dialog box title
	
	// register callback functions
	glutDisplayFunc(glut_render);			//function executed for rendering - renders networks
	glutMouseFunc(glut_mouse);				//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
	if (ind == 1) {							//only in mapping mode, keyboard will be used
		glutCreateMenu(glut_menu);			//register menu option callback
		glutAddMenuEntry("Comparing Mode", 1);	//register menu 1 as comparing mode
		glutAddMenuEntry("Mapping Mode", 2);	//register menu 2 as mapping mode
		glutAddMenuEntry("Exit", 3);			//register menu 3 as exiting
		glutAttachMenu(GLUT_RIGHT_BUTTON);		//register right mouse to open menu option
	}		
	glutKeyboardFunc(glut_keyboard);		//register keyboard callback
	glutMouseWheelFunc(glut_wheel);		

	texture_initialize();									//set up texture mapping (create texture maps, enable features)

	stim::vec3<float> c = bb.center();		//get the center of the network bounding box

	//place the camera along the z-axis at a distance determined by the network size along x and y
	cam.setPosition(c + stim::vec<float>(0, 0, camera_factor * std::max(bb.size()[0], bb.size()[1])));
	cam.LookAt(c[0], c[1], c[2]);						//look at the center of the network
}

#ifdef __CUDACC__
// set specific device to work on
void setdevice(int &device){
	int count;
	cudaGetDeviceCount(&count);				// numbers of device that are available
	if(count < device + 1){
	std::cout<<"No such device available, please set another device"<<std::endl;
	exit(1);
	}
}
#else
void setdevice(int &device){
	device = -1;							// set to default -1
}
#endif

//compare both networks and fill the networks with error information
void compare(float sigma, int device){

	GT = GT.compare(T, sigma, device);						//compare the ground truth to the test case - store errors in GT
    T = T.compare(GT, sigma, device);						//compare the test case to the ground truth - store errors in T

	//calculate the metrics
	float FPR = GT.average();						//calculate the metrics
	float FNR = T.average();
	
	std::cout << "FNR: " << FPR << std::endl;		//print false alarms and misses
	std::cout << "FPR: " << FNR << std::endl;
}

//split and map two networks and fill the networks' R with metric information
void map(float sigma, int device, float threshold){

	// compare and split two networks
	_GT.split(GT, T, sigma, device, threshold);
	_T.split(T, GT, sigma, device, threshold);

	// mapping two new splitted networks and get their edge relation
	_GT.mapping(_T, _gt_t, device, threshold);
	_T.mapping(_GT, _t_gt, device, threshold);

	// generate random color set based on the number of edges in GT
	size_t num = _gt_t.size();							// also create random color for unmapping edge, but won't be used though
	colormap.resize(3 * num);							// 3 portions compound RGB
	for(int i = 0; i < 3 * num; i++)
		colormap[i] = rand()/(float)RAND_MAX;			// set to [0, 1]
	
	//calculate the metrics
	float FPR = _GT.average(0);						//calculate the metrics
	float FNR = _T.average(0);
	
	std::cout << "FNR: " << FPR << std::endl;		//print false alarms and misses
	std::cout << "FPR: " << FNR << std::endl;
}

// writes features of the networks i.e average segment length, tortuosity, branching index, contraction, fractal dimension, number of end and branch points to a csv file
// Pranathi wrote this - saves network features to a CSV file
void features(std::string filename){
		double avgL_t, avgL_gt, avgT_t, avgT_gt, avgB_t, avgB_gt, avgC_t, avgC_gt, avgFD_t, avgFD_gt;
		unsigned int e_t, e_gt, b_gt, b_t;
		avgL_gt = GT.Lengths();
		avgT_gt = GT.Tortuosities();
		avgL_t = T.Lengths();
		avgT_t = T.Tortuosities();
		avgB_gt = GT.BranchingIndex();
		avgB_t = T.BranchingIndex();
		avgC_gt = GT.Contractions();
		avgFD_gt = GT.FractalDimensions();
		avgC_t = T.Contractions();
		avgFD_t = T.FractalDimensions();
		e_gt = GT.EndP();
		e_t = T.EndP();
		b_gt = GT.BranchP();
		b_t = T.BranchP();
		std::ofstream myfile;
		myfile.open (filename.c_str());
		myfile << "Length, Tortuosity, Contraction, Fractal Dimension, Branch Points, End points, Branching Index, \n";
		myfile << avgL_gt << "," << avgT_gt << "," << avgC_gt << "," << avgFD_gt << "," << b_gt << "," << e_gt << "," << avgB_gt <<std::endl;
		myfile << avgL_t << "," << avgT_t << "," << avgC_t << "," << avgFD_t << "," << b_t << "," << e_t << "," << avgB_t <<std::endl;
		myfile.close();
}

// Output an advertisement for the lab, authors, and usage information
void advertise(){
	std::cout<<std::endl<<std::endl;
	std::cout<<"========================================================================="<<std::endl;
	std::cout<<"Thank you for using the NetMets network comparison tool!"<<std::endl;
	std::cout<<"Scalable Tissue Imaging and Modeling (STIM) Lab, University of Houston"<<std::endl;
	std::cout<<"Developers: Pranathi Vemuri, David Mayerich, Jiaming Guo"<<std::endl;
	std::cout<<"Source: https://git.stim.ee.uh.edu/segmentation/netmets" <<std::endl;
	std::cout<<"========================================================================="<<std::endl<<std::endl;

	std::cout<<"usage: netmets file1 file2 --sigma 3"<<std::endl;
	std::cout<<"            compare two .obj files with a tolerance of 3 (units defined by the network)"<<std::endl<<std::endl;
	std::cout<<"       netmets file1 --gui"<<std::endl;
	std::cout<<"            load a file and display it using OpenGL"<<std::endl<<std::endl;
	std::cout<<"       netmets file1 file2 --device 0"<<std::endl;
	std::cout<<"            compare two files using device 0 (if there isn't a gpu, use cpu)"<<std::endl<<std::endl;
	std::cout<<"       netmets file1 file2 --mapping value"<<std::endl;
	std::cout<<"            mapping two files in random colors with a threshold of value"<<std::endl<<std::endl;
}

int main(int argc, char* argv[])
{
	stim::arglist args;						//create an instance of arglist

	//add arguments
	args.add("help", "prints this help");
	args.add("sigma", "force a sigma value to specify the tolerance of the network comparison", "3");
	args.add("gui", "display the network or network comparison using OpenGL");
	args.add("device", "choose specific device to run", "0");
	args.add("features", "save features to a CSV file, specify file name");
	args.add("mapping", "mapping input according to similarity");

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

	if(args["help"].is_set() || args.nargs() == 0){			//test for help
		advertise();										//output the advertisement
		std::cout<<args.str();								//output arguments
		exit(1);											//exit
	}
	
	if (args.nargs() >= 1) {			// if at least one network file is specified
		num_nets = 1;					// set the number of networks to one
		std::vector<std::string> tmp = stim::parser::split(args.arg(0), '.');	// split the filename at '.'
		if ("swc" == tmp[1]) 			// loading swc file
			GT.load_swc(args.arg(0));	// load the specified file as the ground truth
		else if ("obj" == tmp[1])		// loading obj file
			GT.load_obj(args.arg(0));	// load the specified file as the ground truth
		else {
			std::cout << "Invalid loading file" << std::endl;
			exit(1);
		}	
	}

	if (args.nargs() == 2) {								//if two files are specified, they will be displayed in neighboring viewports and compared
		int device = args["device"].as_int();				//get the device value from the user
		num_nets = 2;										//set the number of networks to two
		sigma = args["sigma"].as_float();					//get the sigma value from the user
		radius = sigma;
		std::vector<std::string> tmp = stim::parser::split(args.arg(1), '.');	// split the filename at '.'
		if ("swc" == tmp[1]) 								//loading swc files
			T.load_swc(args.arg(1));                        //load the second (test) network
		else if ("obj" == tmp[1])							//loading obj files
			T.load_obj(args.arg(1));
		else {
			std::cout << "Invalid loading file" << std::endl;
			exit(1);
		}
		if (args["features"].is_set())						//if the user wants to save features
			features(args["features"].as_string());
		//does it need to be resampled??
		GT = GT.resample(resample_rate * sigma);			//resample both networks based on the sigma value
		T = T.resample(resample_rate * sigma);
		if (args["mapping"].is_set()) {
			float threshold = args["mapping"].as_float();
			map(sigma, device, threshold);
		}
		else
			compare(sigma, device);							//run the comparison algorithm
		}

	//if a GUI is requested, display the network using OpenGL
	if(args["gui"].is_set()){
		if (args["mapping"].is_set()) {
			ind = 1;								//set indicator of mapping to 1(true)
			bb = _GT.boundingbox();					//generate a bounding volume		
			glut_initialize();						//create the GLUT window and set callback functions		
			glutMainLoop();							//enter GLUT event processing cycle
		}
		else {
			bb = GT.boundingbox();					//generate a bounding volume		
			glut_initialize();						//create the GLUT window and set callback functions		
			glutMainLoop();							//enter GLUT event processing cycle
		}
	}
}