segmentation / netmets

  #ifndef GL_MULTISAMPLE
  #define GL_MULTISAMPLE  0x809D
  #endif
  
  #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

  unsigned swc_ind = 0;					//indicator of rendering swc file as networks

  // 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

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

  //mouse position tracking
  int mouse_x;
  int mouse_y;
  

  // mouse wheel move
  float des = 0.0f;
  
  // 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;
  

  // create display lists
  GLuint dlist1;
  GLuint dlist2;
  

  //OpenGL objects
  GLuint cmap_tex = 0;				//texture name for the color map
  
  //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
  

  	eye[2] += des;									//get camera closer to target by factor des

  	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) {
  

  	//no mapping, just comparing
  	if (ind == 0) {	

  		//-----working with obj files-----

  		if (swc_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)
  			}

  			if (num_nets == 2) {										//if two networks are loaded	

  				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
  
  				glEnable(GL_TEXTURE_1D);										//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_1D, cmap_tex);							//bind the Brewer texture map
  

  				GT.glCylinder();							//render the GT network

  
  				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();								//render the T network

  			}
  		}

  		//-----working with swc files-----

  		else {
  			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_swc();									//render the GT network (the only one loaded)

  			}

  
  			if (num_nets == 2) {										//if two networks are loaded	

  				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

  				glEnable(GL_TEXTURE_1D);										//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_1D, cmap_tex);							//bind the Brewer texture map
  
  				GT.glCenterline();							//render the GT network

  
  				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.glCenterline();							//render the T network
  			}
  		}
  	}
  

  	//do comparing and mapping

  	else {	

  		//-----working with obj files-----

  		if (swc_ind == 0) {
  			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_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
  
  					glEnable(GL_TEXTURE_1D);										//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_1D, cmap_tex);							//bind the Brewer texture map
  
  					_GT.glCylinder();							//render the GT network
  
  					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();							//render the T network
  
  				}
  				else {
  					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.glRandColorCenterlineGT(dlist1, _gt_t, colormap);
  					_GT.glRandColorCylinder1(dlist1, _gt_t, colormap);
  
  					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.glRandColorCenterlineT(dlist2, _t_gt, colormap);
  					_T.glRandColorCylinder2(dlist2, _t_gt, colormap);
  				}
  			}
  		}
  		//-----working with swc files-----
  		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
  				//in compare mode
  				if (compareMode) {
  					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
  
  					glEnable(GL_TEXTURE_1D);										//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_1D, cmap_tex);							//bind the Brewer texture map
  
  					_GT.glCenterline();							//render the GT network
  
  					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.glCenterline();							//render the T network
  
  				}

  				//in mapping mode

  				else {

  					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.glRandColorCenterlineGT(dlist1, _gt_t, colormap);

  					_GT.glRandColorCylinder1_swc(dlist1, _gt_t, colormap);

  
  					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.glRandColorCenterlineT(dlist2, _t_gt, colormap);

  					_T.glRandColorCylinder2_swc(dlist2, _t_gt, colormap);

  				}

  			}
  		}

  	}
  

  	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

  
  		glDisable(GL_TEXTURE_1D);

  		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.0, 1.0, 0.0);											// using green 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();
  	}

  	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) {
  	float cam_move_fac;								// camera move unit length
  	stim::vec3<float> eye = cam.getPosition();		// get the camera position (eye point)
  	stim::vec3<float> focus = cam.getLookAt();		// get the camera focal point
  	cam_move_fac = fabs(focus[2] - eye[2]) / 50;	// divided by 50
  	
  	if (direction > 0)								// if it is button 3(up)
  		des -= cam_move_fac;
  	else											// if it is button 4(down)
  		des += cam_move_fac;
  	glutPostRedisplay();
  }
  
  // define keyboard inputs

  void glut_keyboard(unsigned char key, int x, int y){

  	if(key == 'm')											// if keyboard 'm' is pressed

  	{
  		if(compareMode && !mappingMode){					// if it is in compare mode
  			compareMode = false;
  			mappingMode = true;
  		}
  		else{												// if it is in mapping mode
  			compareMode = true;
  			mappingMode = false;
  		}
  	}

  	if (key == 27) {										// if keyboard "ESC" is pressed
  		exit(0);											// exit
  	}

  	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 | GLUT_MULTISAMPLE);	//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

  	glEnable(GL_MULTISAMPLE);

  	// 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
  		glutKeyboardFunc(glut_keyboard);	//register keyboard callback
  		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
  	}						
  	if (swc_ind == 1)						//only in rendering swc files, mouse wheel will be used
  		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;
  }
  

  void map_swc(float sigma, int device, float threshold) {
  
  	// compare two networks
  	_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
  
  	// mapping two 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
  			swc_ind = 1;				// set the indicator of swc file to 1
  		}			
  		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
  		if (1 == swc_ind) {									//loading swc files

  			int device = args["device"].as_int();				//get the device value from the user
  			num_nets = 2;										//set the number of networks to two

  																//does it need to be resampled??

  			float sigma = args["sigma"].as_float();				//get the sigma value from the user
  			T.load_swc(args.arg(1));                            //load the second (test) network
  			if (args["features"].is_set())						//if the user wants to save features
  				features(args["features"].as_string());
  			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_swc(sigma, device, threshold);
  				//std::cout << "right now networks that are loaded from swc files do not need to be mapped with each other!" << std::endl;
  				//exit(1);

  			}
  			else
  				compare(sigma, device);							//run the comparison algorithm
  		}

  		else {

  			int device = args["device"].as_int();				//get the device value from the user
  			num_nets = 2;										//set the number of networks to two
  			float sigma = args["sigma"].as_float();				//get the sigma value from the user
  			T.load_obj(args.arg(1));                            //load the second (test) network
  			if (args["features"].is_set())						//if the user wants to save features
  				features(args["features"].as_string());
  			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

  		}
  	}

  }