segmentation / netmets-old

#include "rtsFiberNetwork.h"
#include "rts/rts_glShaderProgram.h"
#include "GL/glut.h"
#include "rts/rts_glutRenderWindow.h"
#include <time.h>

extern void ComputeNetMets();
extern rtsFiberNetwork* goldNetwork;
extern rtsFiberNetwork* testNetwork;
extern float sigmaG, sigmaC;
float network_span;
CoreGraphList coreGraph;
vector<point3D<float> > sequenceColors; 
int current_sequence = 0;

//shader variables
rts_glShaderProgram Edge_ErrorShader;
rts_glShaderProgram Node_ErrorShader;
rts_glShaderProgram Smooth_Shader;

//display lists
GLuint GT_FibersList=0;
GLuint T_FibersList=0;
GLuint GT_EndCaps=0;
GLuint T_EndCaps=0;
GLuint GT_NodesList=0;
GLuint T_NodesList=0;
GLuint T_PathList=0;
GLuint GT_PathList=0;

//drawing variables
int tube_subdivisions = 20;
float node_radius_factor = 0.7;
float fiber_radius_factor = 0.5;
float cull_test_case_threshold = 1.0;

#define DISPLAY_GT_NETWORK	1
#define DISPLAY_T_NETWORK	2
#define DISPLAY_GT_GRAPH	3
#define DISPLAY_T_GRAPH		4
#define DISPLAY_GT_SELECTED	5
#define DISPLAY_T_SELECTED	6

//menu options
#define NETCOMP_EXIT		0
#define DISPLAY_NETWORK		1
#define DISPLAY_GRAPH		2
#define DISPLAY_CONNECTED	3
#define DISPLAY_SELECTED	4
#define COLORMAP_ISOLUMINANT	5
#define COLORMAP_BLACKBODY		6
#define COLORMAP_BREWER			7
#define COLORMAP_POLAR_CIELAB	8
#define COLORMAP_RAINBOW		9
#define CULL_TEST_CASE			10
#define RECOMPUTE_METRIC		11


//fibers to render in Graph Mode
//list<int> T_DisplayEdges;
//list<int> GT_DisplayEdges;

int DisplayMode = DISPLAY_NETWORK;
float L0_pos[3];
float L1_pos[3];
//GLuint texColorMap=0;
rts_glTextureMap texColorMap;




void makeColormap(int ColorMapType = COLORMAP_BREWER)
{
	//if(texColorMap != 0)
	//	glDeleteTextures(1, &texColorMap);

	point3D<float>* ctrlPts;
	int num_points = 0;
	if(ColorMapType == COLORMAP_ISOLUMINANT)
	{
		//allocate memory for the colormap
		num_points = 2;
		ctrlPts = new point3D<float>[num_points];
		//memset(ctrlPts, 0, num_points*sizeof(point3D<float>));

		ctrlPts[0] = point3D<float>(0.0, 1.0, 0.0);
		ctrlPts[1] = point3D<float>(1.0, 0.0, 0.0);
	}
	else if(ColorMapType == COLORMAP_RAINBOW)
	{
		//allocate memory for the colormap
		num_points = 5;
		ctrlPts = new point3D<float>[num_points];
		//memset(ctrlPts, 0, num_points*sizeof(point3D<float>));

		//ctrlPts[0] = point3D<float>(0.7, 0, 0.7);
		ctrlPts[0] = point3D<float>(0, 0, 1);
		ctrlPts[1] = point3D<float>(0, 0.7, 0.7);
		ctrlPts[2] = point3D<float>(0, 1, 0);
		ctrlPts[3] = point3D<float>(0.7, 0.7, 0);
		ctrlPts[4] = point3D<float>(1, 0, 0);
	}
	else if(ColorMapType == COLORMAP_BLACKBODY)
	{
		//allocate memory for the colormap
		num_points = 4;
		ctrlPts = new point3D<float>[num_points];
		//memset(ctrlPts, 0, num_points*sizeof(point3D<float>));

		ctrlPts[0] = point3D<float>(0.0, 0.0, 0.0);
		ctrlPts[1] = point3D<float>(1.0, 0.0, 0.0);
		ctrlPts[2] = point3D<float>(1.0, 1.0, 0.0);
		ctrlPts[3] = point3D<float>(1.0, 1.0, 1.0);
	}
	else if(ColorMapType == COLORMAP_BREWER)
	{
		//allocate memory for the colormap
		num_points = 11;
		ctrlPts = new point3D<float>[num_points];
		//memset(ctrlPts, 0, num_points*sizeof(point3D<float>));

		ctrlPts[0] = point3D<float>(0.192157, 0.211765, 0.584314);
		ctrlPts[1] = point3D<float>(0.270588, 0.458824, 0.705882);
		ctrlPts[2] = point3D<float>(0.454902, 0.678431, 0.819608);
		ctrlPts[3] = point3D<float>(0.670588, 0.85098, 0.913725);
		ctrlPts[4] = point3D<float>(0.878431, 0.952941, 0.972549);
		ctrlPts[5] = point3D<float>(1, 1, 0.74902);
		ctrlPts[6] = point3D<float>(0.996078, 0.878431, 0.564706);
		ctrlPts[7] = point3D<float>(0.992157, 0.682353, 0.380392);
		ctrlPts[8] = point3D<float>(0.956863, 0.427451, 0.262745);
		ctrlPts[9] = point3D<float>(0.843137, 0.188235, 0.152941);
		ctrlPts[10] = point3D<float>(0.647059, 0, 0.14902);

	}
		else if(ColorMapType == COLORMAP_POLAR_CIELAB)
	{
		//allocate memory for the colormap
		num_points = 33;
		ctrlPts = new point3D<float>[num_points];
		//memset(ctrlPts, 0, num_points*sizeof(point3D<float>));

		ctrlPts[0] = point3D<float>(0.07514311, 0.468049805,1);
		ctrlPts[1] = point3D<float>(0.247872569, 0.498782363,1);
		ctrlPts[2] = point3D<float>(0.339526309, 0.528909511,1);
		ctrlPts[3] = point3D<float>(0.409505078, 0.558608486,1);
		ctrlPts[4] = point3D<float>(0.468487184, 0.588057293,1);
		ctrlPts[5] = point3D<float>(0.520796675, 0.617435078,1);
		ctrlPts[6] = point3D<float>(0.568724526, 0.646924167,1);
		ctrlPts[7] = point3D<float>(0.613686735, 0.676713218,1);
		ctrlPts[8] = point3D<float>(0.656658579, 0.707001303,1);
		ctrlPts[9] = point3D<float>(0.698372844, 0.738002964,1);
		ctrlPts[10] = point3D<float>(0.739424025, 0.769954435,1);
		ctrlPts[11] = point3D<float>(0.780330104, 0.803121429,1);
		ctrlPts[12] = point3D<float>(0.821573924, 0.837809045,1);
		ctrlPts[13] = point3D<float>(0.863634967, 0.874374691,1);
		ctrlPts[14] = point3D<float>(0.907017747, 0.913245283,1);
		ctrlPts[15] = point3D<float>(0.936129275, 0.938743558, 0.983038586);
		ctrlPts[16] = point3D<float>(0.943467973, 0.943498599, 0.943398095);
		ctrlPts[17] = point3D<float>(0.990146732, 0.928791426, 0.917447482);
		ctrlPts[18] = point3D<float>(1, 0.88332677, 0.861943246);
		ctrlPts[19] = point3D<float>(1, 0.833985467, 0.803839606);
		ctrlPts[20] = point3D<float>(1, 0.788626485, 0.750707739);
		ctrlPts[21] = point3D<float>(1, 0.746206642, 0.701389973);
		ctrlPts[22] = point3D<float>(1, 0.70590052, 0.654994046);
		ctrlPts[23] = point3D<float>(1, 0.667019783, 0.610806959);
		ctrlPts[24] = point3D<float>(1, 0.6289553, 0.568237474);
		ctrlPts[25] = point3D<float>(1, 0.591130233, 0.526775617);
		ctrlPts[26] = point3D<float>(1, 0.552955184, 0.485962266);
		ctrlPts[27] = point3D<float>(1, 0.513776083, 0.445364274);
		ctrlPts[28] = point3D<float>(1, 0.472800903, 0.404551679);
		ctrlPts[29] = point3D<float>(1, 0.428977855, 0.363073592);
		ctrlPts[30] = point3D<float>(1, 0.380759558, 0.320428137);
		ctrlPts[31] = point3D<float>(0.961891484, 0.313155629, 0.265499262);
		ctrlPts[32] = point3D<float>(0.916482116, 0.236630659, 0.209939162);

	}

	texColorMap.Init(ctrlPts, GL_TEXTURE_1D, num_points, 0, 0, GL_RGB, GL_RGB, GL_FLOAT);

	//glGenTextures(1, &texColorMap);
	//glBindTexture(GL_TEXTURE_1D, texColorMap);
	//glTexImage1D(GL_TEXTURE_1D, 0, GL_RGB, num_points, 0, GL_RGB, GL_FLOAT, ctrlPts);
}



point3D<float> HSLtoRGB(point3D<float> HSL)
{
	float H = HSL.x;
	float S = HSL.y;
	float L = HSL.z;
	
	float temp2;
	if(L < 0.5)
		temp2 = L*(1.0+S);
	else
		temp2 = L+S - L*S;
		
	float temp1 = 2.0*L - temp2;
	
	point3D<float> temp3(H+1.0/3.0, H, H-1.0/3.0);
	if(temp3.x < 0.0) temp3.x = temp3.x+1.0;
	if(temp3.y < 0.0) temp3.y = temp3.y+1.0;
	if(temp3.z < 0.0) temp3.z = temp3.z+1.0;
	
	if(temp3.x > 1.0) temp3.x = temp3.x - 1.0;
	if(temp3.y > 1.0) temp3.y = temp3.y - 1.0;
	if(temp3.z > 1.0) temp3.z = temp3.z - 1.0;
	
	point3D<float> result;
	if(6.0*temp3.x < 1.0) result.x = temp1 +(temp2 - temp1)*6.0*temp3.x;
	else if(2.0*temp3.x < 1.0) result.x = temp2;
	else if(3.0*temp3.x < 2.0) result.x = temp1+(temp2-temp1)*((2.0/3.0) - temp3.x)*6.0;
	else result.x = temp1;
	
	if(6.0*temp3.y < 1.0) result.y = temp1 +(temp2 - temp1)*6.0*temp3.y;
	else if(2.0*temp3.y < 1.0) result.y = temp2;
	else if(3.0*temp3.y < 2.0) result.y = temp1+(temp2-temp1)*((2.0/3.0) - temp3.y)*6.0;
	else result.y = temp1;
	
	if(6.0*temp3.z < 1.0) result.z = temp1 +(temp2 - temp1)*6.0*temp3.z;
	else if(2.0*temp3.z < 1.0) result.z = temp2;
	else if(3.0*temp3.z < 2.0) result.z = temp1+(temp2-temp1)*((2.0/3.0) - temp3.z)*6.0;
	else result.z = temp1;
	
	//result.a = 0.0;
	return result;
}
void ColorFibers()
{
	//srand(time(NULL));
	sequenceColors.clear();
	//for each fiber
	for(CoreGraphList::iterator i = coreGraph.begin(); i!=coreGraph.end(); i++)
	{
		float random_hue = (double)rand()/(double)RAND_MAX;
		//cout<<"Random Hue: "<<random_hue<<endl;
		float random_saturation = 1.0;//(double)rand()/(double)RAND_MAX;
		point3D<float> rgb = HSLtoRGB(point3D<float>(random_hue, random_saturation, 0.5));
		//point3D<float> rgb((double)rand()/(double)RAND_MAX, (double)rand()/(double)RAND_MAX, (double)rand()/(double)RAND_MAX);
		sequenceColors.push_back(rgb);
	}
	

}
void CenterCameraToSelected()
{
	if(coreGraph.size() == 0)
		return;

	//center the fiber in both networks
	point3D<float> min_pt(9999, 9999, 9999);
	point3D<float> max_pt(-9999, -9999, -9999);

	//iterate through the first edge sequence
	EdgeSequence::iterator i;
	int node;
	point3D<float> test;
	for(i=coreGraph[current_sequence].first.begin(); i!=coreGraph[current_sequence].first.end(); i++)
	{
		node = testNetwork->FiberList[*i].n0;
		test = testNetwork->NodeList[node].p;
		min_pt.x = min(test.x, min_pt.x);
		min_pt.y = min(test.y, min_pt.y);
		min_pt.z = min(test.z, min_pt.z);	
		max_pt.x = max(test.x, max_pt.x);
		max_pt.y = max(test.y, max_pt.y);
		max_pt.z = max(test.z, max_pt.z);

		node = testNetwork->FiberList[*i].n1;
		test = testNetwork->NodeList[node].p;
		min_pt.x = min(test.x, min_pt.x);
		min_pt.y = min(test.y, min_pt.y);
		min_pt.z = min(test.z, min_pt.z);	
		max_pt.x = max(test.x, max_pt.x);
		max_pt.y = max(test.y, max_pt.y);
		max_pt.z = max(test.z, max_pt.z);
	}
	point3D<float> middle = min_pt+0.5*(max_pt - min_pt);

	rts_glut_camera.LookAt(middle);

}
void IncrementSelectedFiber(int i)
{
	//get the currently selected fiber id
	if(coreGraph.size() <= 0)
		return;

	//get the number of fibers
	int end_id = coreGraph.size();

	current_sequence+=i;
	if(current_sequence >= end_id)
		current_sequence = 0;
	if(current_sequence < 0)
		current_sequence = coreGraph.size()-1;

	//print the selected edges
	EdgeSequence::iterator EdgeI;

	for(EdgeI = coreGraph[current_sequence].first.begin(); EdgeI != coreGraph[current_sequence].first.end(); EdgeI++)
		cout<<*EdgeI<<" ";
	cout<<"--->";
	for(EdgeI = coreGraph[current_sequence].second.begin(); EdgeI != coreGraph[current_sequence].second.end(); EdgeI++)
		cout<<*EdgeI<<" ";
	cout<<endl;


	CenterCameraToSelected();

}
void DrawNodeSphere(rtsFiberNetwork* network, int n, float radius)
{
	GLUquadricObj* quadric = gluNewQuadric();
	gluQuadricNormals(quadric, GLU_SMOOTH);

	glMatrixMode(GL_MODELVIEW);

	glPushMatrix();

	point3D<float> p;
	//glColor3f(network->NodeList[n].error, 0.0, 0.0);
	p = network->NodeList[n].p;

	glTranslatef(p.x, p.y, p.z);
	//glutSolidSphere(node_radius*standard_deviation, 20, 20);
	glTexCoord1f(network->NodeList[n].error);
	if(network->NodeList[n].color < 0)
		glColor3f(1.0, 0.0, 0.0);
	else
		glColor3f(1.0, 1.0, 1.0);
	gluSphere(quadric,radius,32,32);

	glPopMatrix();

}
void DrawNodeSpheres(rtsFiberNetwork* network, float radius)
{


	unsigned int n;
	for(n=0; n != network->FiberList.size(); n++)
	{
		if(!network->isCulled(n))
		{
			DrawNodeSphere(network, network->FiberList[n].n0,radius);
			DrawNodeSphere(network, network->FiberList[n].n1,radius);
		}
	}
}

void FrenetFrame(vector3D<float> &x, vector3D<float> &y, vector3D<float> &z)
{
	x = vector3D<float>(0.0, 0.0, 1.0);
	y = x.X(z);
	x = z.X(y);
	x.Normalize();
	y.Normalize();
	z.Normalize();
}

vector3D<float> GetColor(float error)
{
	//This function converts an error value to a color
	//The conversion is done by creating an HSV color from the error value and converting that HSV color to RGB
	float H = (240.0/60.0)*(1.0 - error);
	float S = 1.0;
	float V = 1.0;

	int i = floor(H);
	float f = H - i;
	if(i%2 == 0)
		f = 1-f;
	float m = V*(1 - S);
	float n = V*(1-S*f);
	switch(i)
	{
	case 0:
		return vector3D<float>(V, n, m);
	case 1:
		return vector3D<float>(n, V, m);
	case 2:
		return vector3D<float>(m, V, n);
	case 3:
		return vector3D<float>(m, n, V);
	case 4:
		return vector3D<float>(n, m, V);
	case 5:
		return vector3D<float>(V, m, n);
	default:
		return vector3D<float>(0, 0, 0);
	}

}

void DrawTube(point3D<float> p0, vector3D<float> d0, point3D<float> p1, vector3D<float> d1, float error0, float error1, float radius, int subdiv)
{
	
	//draw the first circle
	vector3D<float> x0, y0, z0, x1, y1, z1;
	z0 = d0;
	FrenetFrame(x0, y0, z0);

	z1 = d1;
	FrenetFrame(x1, y1, z1);

	float t_step = (2*3.14159)/subdiv;

	
	float u, v;

	//get the RGB color
	point3D<float> circle0, circle1;
	vector3D<float> RGB0, RGB1;
	vector3D<float> normal;
	
	//RGB0 = GetColor(color0);	
	//RGB1 = GetColor(color1);

	glBegin(GL_TRIANGLE_STRIP);
	for(int t=0; t<=subdiv; t++)
	{
		u = radius * cos(t*t_step);
		v = radius * sin(t*t_step);
		normal = u*x0 + v*y0;		
		circle0 = p0 + normal;
		normal.Normalize();

		glTexCoord1f(error0);
		//glColor4f(error0, 0.0, 0.0, 1.0);
		glNormal3f(normal.x, normal.y, normal.z);
		glVertex3f(circle0.x, circle0.y, circle0.z);

		normal = u*x1 + v*y1;
		circle1 = p1 + normal;
		normal.Normalize();

		glTexCoord1f(error1);
		//glColor4f(error1, 0.0, 0.0, 1.0);
		glNormal3f(normal.x, normal.y, normal.z);
		glVertex3f(circle1.x, circle1.y, circle1.z);

	}
	glEnd();
	CHECK_OPENGL_ERROR
}

void ExtrudeFiber(rtsFiberNetwork* network, int fiber, float radius)
{
	vector3D<float> x, y, z;
	point3D<float> p0, p1, p2, p3;
	vector3D<float> d1, d2;
	float e1, e2;

	//get the first point
	int node = network->FiberList[fiber].n0;
	p1 = network->NodeList[node].p;
	e1 = network->NodeList[node].error;

	//for each vertex in the fiber
	int num_points = (int)network->FiberList[fiber].pointList.size();
	for(int v=0; v<num_points; v++)
	{
		//get the next point
		p2 = network->FiberList[fiber].pointList[v];
		e2 = network->FiberList[fiber].errorList[v];
		
		if(v<num_points-1)
			p3 = network->FiberList[fiber].pointList[v+1];
		else
		{
			node = network->FiberList[fiber].n1;
			p3 = network->NodeList[node].p;
		}

		d2 = p3-p1;

		//compute the fiber derivatives at p1 and p2
		if(v==0)	//if this is the first fiber
			d1 = p2 - p1;
		else
		{
			d1 = p2 - p0;
		}

		DrawTube(p1, d1, p2, d2, e1, e2, radius, tube_subdivisions);

		//shift
		p0 = p1;
		p1 = p2;
		e1 = e2;
	}
	//make the last tube

	//if there were any points in the pointlist
	if(num_points > 0)
	{
		p2 = p3;
		node = network->FiberList[fiber].n1;
		e2 = network->NodeList[node].error;
		d1 = p2-p0;
		d2 = p2-p1;
		DrawTube(p1, d1, p2, d2, e1, e2, radius, tube_subdivisions);
	}
	//if there are only the two node points
	else
	{
		node = network->FiberList[fiber].n1;
		p2 = network->NodeList[node].p;
		e2 = network->NodeList[node].error;
		d1 = p2 - p1;
		d2 = p2 - p1;
		DrawTube(p1, d1, p2, d2, e1, e2, radius, tube_subdivisions);
	}
}

void DrawLineFiber(rtsFiberNetwork* network, int f)
{
	point3D<float> p;
	int node = network->FiberList[f].n0;
	p = network->NodeList[node].p;

	glBegin(GL_LINE_STRIP);
	glVertex3f(p.x, p.y, p.z);
	for(int v=0; v!=network->FiberList[f].pointList.size(); v++)
	{
		p = network->FiberList[f].pointList[v];
		glVertex3f(p.x, p.y, p.z);		
	}
	node = network->FiberList[f].n1;
	p = network->NodeList[node].p;
	glVertex3f(p.x, p.y, p.z);
	glEnd();

}
void DrawLineNetwork(rtsFiberNetwork* network)
{

	int num_fibers = network->FiberList.size();
	for(int f = 0; f < num_fibers; f++)
	{
		/*if(network->FiberList[f].mapped_to == -1)
			glColor3f(0.0,0.0, 0.0);
		else
			glColor3f(1.0, 0.0, 0.0);
		*/
		DrawLineFiber(network, f);
		CHECK_OPENGL_ERROR
	}


}
void DrawGraphNodes(rtsFiberNetwork* network)
{
	//renders graph nodes, colored based on their node color
	glMatrixMode(GL_MODELVIEW);

	unsigned int n;
	for(n=0; n != network->NodeList.size(); n++)
	{
		glPushMatrix();

		point3D<float> p;
		/*if(network->NodeList[n].color < 0)
			glColor3f(1.0, 0.0, 0.0);
		else
			glColor3f(0.0, 1.0, 0.0);
			*/
		//glColor3f(network->NodeList[n].error, 0.0, 0.0);
		p = network->NodeList[n].p;

		glTranslatef(p.x, p.y, p.z);
		glutSolidSphere(node_radius_factor*sigmaC, 20, 20);

		glPopMatrix();
	}
}
void DrawFiberSequence(rtsFiberNetwork* network, EdgeSequence sequence, float fiber_radius, float node_radius)
{
	//glClear(GL_DEPTH_BUFFER_BIT);
	for(EdgeSequence::iterator i = sequence.begin(); i != sequence.end(); i++)
	{
		ExtrudeFiber(network, *i, fiber_radius);
		glPushAttrib(GL_CURRENT_BIT);
		DrawNodeSphere(network, network->FiberList[*i].n0, node_radius);
		DrawNodeSphere(network, network->FiberList[*i].n1, node_radius);
		glPopAttrib();
	}


}
GLuint CreateFiberDisplayList(rtsFiberNetwork* network, float radius)
{
	GLuint result = glGenLists(1);
	glNewList(result, GL_COMPILE);

	int num_fibers = network->FiberList.size();
	for(int f = 0; f < num_fibers; f++)
	{
		if(!network->isCulled(f))
		{
			ExtrudeFiber(network, f, radius);
			CHECK_OPENGL_ERROR
		}
	}
	glEndList();
	return result;
}

GLuint CreateNodeDisplayList(rtsFiberNetwork* network, float radius)
{
	GLuint result = glGenLists(1);
	glNewList(result, GL_COMPILE);
	DrawNodeSpheres(network, radius);
	glEndList();
	return result;
}


void CreateFiberPathLists(float fiber_radius, float node_radius)
{
	GT_PathList = glGenLists(1);
	glNewList(GT_PathList, GL_COMPILE);

	if(coreGraph.size() > 0)
	{
		for(unsigned int i=0; i<sequenceColors.size(); i++)
		{
			point3D<float> rgb = sequenceColors[i];
			glColor3f(rgb.x, rgb.y, rgb.z);
			DrawFiberSequence(goldNetwork, coreGraph[i].second, fiber_radius, node_radius);
		}
	}
	glEndList();

	T_PathList = glGenLists(1);
	glNewList(T_PathList, GL_COMPILE);

	if(coreGraph.size() > 0)
	{
		for(unsigned int i=0; i<sequenceColors.size(); i++)
		{
			point3D<float> rgb = sequenceColors[i];
			glColor3f(rgb.x, rgb.y, rgb.z);
			DrawFiberSequence(testNetwork, coreGraph[i].first, fiber_radius, node_radius);
		}
	}
	glEndList();

}
void CreateDisplayLists()
{
	if(GT_FibersList != 0)
		glDeleteLists(GT_FibersList, 1);
	if(T_FibersList != 0)
		glDeleteLists(T_FibersList, 1);
	if(GT_NodesList != 0)
		glDeleteLists(GT_NodesList, 1);
	if(T_NodesList != 0)
		glDeleteLists(T_NodesList, 1);
	if(GT_EndCaps != 0)
		glDeleteLists(GT_EndCaps, 1);
	if(T_EndCaps != 0)
		glDeleteLists(T_EndCaps, 1);

	//create the display lists
	GT_FibersList = CreateFiberDisplayList(goldNetwork, sigmaG*fiber_radius_factor);
	T_FibersList = CreateFiberDisplayList(testNetwork, sigmaG*fiber_radius_factor);
	GT_NodesList = CreateNodeDisplayList(goldNetwork, sigmaG*node_radius_factor);
	T_NodesList = CreateNodeDisplayList(testNetwork, sigmaG*node_radius_factor);
	GT_EndCaps = CreateNodeDisplayList(goldNetwork, sigmaG*fiber_radius_factor);
	T_EndCaps = CreateNodeDisplayList(testNetwork, sigmaG*fiber_radius_factor);

	if(GT_PathList != 0)
		glDeleteLists(GT_PathList,1);
	if(T_PathList != 0)
		glDeleteLists(T_PathList,1);
	CreateFiberPathLists(fiber_radius_factor*sigmaG, node_radius_factor*sigmaG);
}
void RenderViewCamera()
{
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	//compute the aspect ratio
	float aspect_ratio = (float)glutGet(GLUT_WINDOW_WIDTH)/2.0/(float)glutGet(GLUT_WINDOW_HEIGHT);
	gluPerspective(rts_glut_camera.getFOV(), aspect_ratio, network_span/10, network_span*10);

	//render the camera
	glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();

	point3D<float> camera_position = rts_glut_camera.getPosition();
	vector3D<float> camera_up = rts_glut_camera.getUp();
	point3D<float> camera_lookat = rts_glut_camera.getLookAt();
	gluLookAt(camera_position.x,
			  camera_position.y,
			  camera_position.z,
			  camera_lookat.x,
			  camera_lookat.y,
			  camera_lookat.z,
			  camera_up.x,
			  camera_up.y,
			  camera_up.z);


	//get the light positions (lights move with the camera)
	vector3D<float> up = rts_glut_camera.getUp();
	vector3D<float> dir = rts_glut_camera.getDirection();
	vector3D<float> side = dir.X(up);
	L0_pos[0] = side.x;
	L0_pos[1] = side.y;
	L0_pos[2] = side.z;

	L1_pos[0] = -dir.x;
	L1_pos[1] = -dir.y;
	L1_pos[2] = -dir.z;

	//scale the viewport to the network
	/*vector3D<float> span = max_point - min_point;
	float scale = 1.0/span.Length();
	//compute center point
	point3D<float> center = min_point + 0.5*span;
	glScalef(scale, scale, scale);
	glTranslatef(-center.x, -center.y, -center.z);*/


}
void RecenterCamera()
{
	point3D<float> min_point0 = goldNetwork->min_pos;
	point3D<float> min_point1 = testNetwork->min_pos;

	point3D<float> max_point0 = goldNetwork->max_pos;
	point3D<float> max_point1 = testNetwork->max_pos;

	point3D<float> min_point(min(min_point0.x, min_point1.x), min(min_point0.y, min_point1.y), min(min_point0.z, min_point1.z));
	point3D<float> max_point(max(max_point0.x, max_point1.x), max(max_point0.y, max_point1.y), max(max_point0.z, max_point1.z));
	point3D<float> center_point = min_point + 0.5*(max_point - min_point);


	network_span = (max_point - point3D<float>(0, 0, 0)).Length();
	rts_glut_camera.setPosition(0, 0, 3*network_span);
	rts_glut_camera.LookAt(center_point, vector3D<float>(0.0, 1.0, 0.0));

}
void DrawNetwork(int Display)
{
	//draw the network fibers
	switch(Display)
	{
	case DISPLAY_GT_NETWORK:
		Edge_ErrorShader.UpdateGlobalUniforms();
		Edge_ErrorShader.BeginProgram();
		glCallList(GT_FibersList);
		Edge_ErrorShader.EndProgram();

		Node_ErrorShader.UpdateGlobalUniforms();
		Node_ErrorShader.BeginProgram();
		glCallList(GT_EndCaps);
		Node_ErrorShader.EndProgram();
		break;
	case DISPLAY_T_NETWORK:
		Edge_ErrorShader.UpdateGlobalUniforms();
		Edge_ErrorShader.BeginProgram();
		glCallList(T_FibersList);	
		Edge_ErrorShader.EndProgram();

		Node_ErrorShader.UpdateGlobalUniforms();
		Node_ErrorShader.BeginProgram();
		glCallList(T_EndCaps);	
		Node_ErrorShader.EndProgram();
		break;
	case DISPLAY_GT_GRAPH:
		glColor3f(1.0, 1.0, 1.0);
		Smooth_Shader.UpdateGlobalUniforms();
		Smooth_Shader.BeginProgram();
		glCallList(GT_PathList);
		glColor3f(1.0, 1.0, 1.0);
		glCallList(GT_FibersList);
		glCallList(GT_NodesList);
		Smooth_Shader.EndProgram();
		break;
	case DISPLAY_T_GRAPH:
		glColor3f(1.0, 1.0, 1.0);
		Smooth_Shader.UpdateGlobalUniforms();
		Smooth_Shader.BeginProgram();
		glCallList(T_PathList);
		glColor3f(1.0, 1.0, 1.0);
		glCallList(T_FibersList);
		glCallList(T_NodesList);
		Smooth_Shader.EndProgram();
		break;
	case DISPLAY_GT_SELECTED:
		glColor3f(1.0, 1.0, 1.0);
		Smooth_Shader.UpdateGlobalUniforms();
		Smooth_Shader.BeginProgram();
		glCallList(GT_FibersList);
		glCallList(GT_NodesList);

		glClear(GL_DEPTH_BUFFER_BIT);
		glColor3f(1.0, 0.0, 1.0);
		DrawFiberSequence(goldNetwork, coreGraph[current_sequence].second, fiber_radius_factor*sigmaG, node_radius_factor*sigmaG);
		Smooth_Shader.EndProgram();
		break;
	case DISPLAY_T_SELECTED:
		glColor3f(1.0, 1.0, 1.0);
		Smooth_Shader.UpdateGlobalUniforms();
		Smooth_Shader.BeginProgram();
		glCallList(T_FibersList);
		glCallList(T_NodesList);

		glClear(GL_DEPTH_BUFFER_BIT);
		glColor3f(1.0, 0.0, 1.0);
		DrawFiberSequence(testNetwork, coreGraph[current_sequence].first, fiber_radius_factor*sigmaG, node_radius_factor*sigmaG);
		Smooth_Shader.EndProgram();
		break;
	default:
		break;
	}

}




void MyDisplayFunction()
{
	//glClearColor(0.95, 1.0, 0.95, 1.0);
	glClearColor(1.0, 1.0, 1.0, 1.0);
	glClear(GL_COLOR_BUFFER_BIT);
	glClear(GL_DEPTH_BUFFER_BIT);

	//set left viewport
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	glViewport(0, 0, glutGet(GLUT_WINDOW_WIDTH)/2, glutGet(GLUT_WINDOW_HEIGHT));

	RenderViewCamera();

	
	if(DisplayMode == DISPLAY_NETWORK)
		DrawNetwork(DISPLAY_GT_NETWORK);
	else if(DisplayMode == DISPLAY_GRAPH)
		DrawNetwork(DISPLAY_GT_GRAPH);
	else if(DisplayMode == DISPLAY_SELECTED)
		DrawNetwork(DISPLAY_GT_SELECTED);

	//set right viewport
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	glViewport(glutGet(GLUT_WINDOW_WIDTH)/2, 0, glutGet(GLUT_WINDOW_WIDTH)/2, glutGet(GLUT_WINDOW_HEIGHT));
	
	RenderViewCamera();
	//glClear(GL_COLOR_BUFFER_BIT);
	//glClear(GL_DEPTH_BUFFER_BIT);

	if(DisplayMode == DISPLAY_NETWORK)
		DrawNetwork(DISPLAY_T_NETWORK);
	else if(DisplayMode == DISPLAY_GRAPH)
		DrawNetwork(DISPLAY_T_GRAPH);
	else if(DisplayMode == DISPLAY_SELECTED)
		DrawNetwork(DISPLAY_T_SELECTED);


	
	//glutSwapBuffers();
}

void GlutMenuCallback(int option)
{
	if(option == NETCOMP_EXIT)
		exit(1);
	if(option >= DISPLAY_NETWORK && option <= DISPLAY_SELECTED)
		DisplayMode = option;
	if(option >= COLORMAP_ISOLUMINANT && option <= COLORMAP_RAINBOW)
		makeColormap(option);
	if(option == CULL_TEST_CASE)
	{
		//get the new threshold
		cout<<"Enter new test case fiber threshold [0 1]: ";
		float cull_value;
		cin>>cull_value;
		testNetwork->setCullValue(cull_value);

		//re-create the display lists
		ComputeNetMets();
		CreateDisplayLists();

	}
	if(option == RECOMPUTE_METRIC)
	{
		cout<<"Please enter a sigma value: ";
		cin>>sigmaG;
		sigmaC = sigmaG;

		ComputeNetMets();
		CreateDisplayLists();
	}

	

	/*
	switch(option)
	{
	case NETCOMP_EXIT:
		exit(1);
		break;
	case NETCOMP_VIEW_GT:
		DisplayMode = DISPLAY_GT_NETWORK;
		break;
	case NETCOMP_VIEW_T:
		DisplayMode = DISPLAY_T_NETWORK;
		break;
	default:
		break;
	}*/

}