#ifndef STIM_GL_NETWORK
#define STIM_GL_NETWORK

#include <stim/biomodels/network.h>
#include <stim/visualization/aaboundingbox.h>

namespace stim{

template <typename T>
class gl_network : public stim::network<T>{

protected:
	using stim::network<T>::E;
	using stim::network<T>::V;

	GLuint dlist;

public:

	/// Default constructor
	gl_network() : stim::network<T>(){
		dlist = 0;
	}

	/// Constructor creates a gl_network from a stim::network
	gl_network(stim::network<T> N) : stim::network<T>(N){
		dlist = 0;
	}

	/// Fills the parameters with the minimum and maximum spatial positions in the network,
	///     specifying a bounding box for the network geometry
	aaboundingbox<T> boundingbox(){

		aaboundingbox<T> bb;								//create a bounding box

		//loop through every edge
		for(unsigned e = 0; e < E.size(); e++){
			//loop through every point
			for(unsigned p = 0; p < E[e].size(); p++)
				bb.expand(E[e][p]);						//expand the bounding box to include the point
		}

		return bb;								//return the bounding box
	}

	//void renderCylinder(T x1, T y1, T z1, T x2, T y2, T z2, T radius, int subdivisions) {
	//	T dx = x2 - x1;
	//	T dy = y2 - y1;
	//	T dz = z2 - z1;
	//	/// handle the degenerate case with an approximation
	//	if (dz == 0)
	//		dz = .00000001;
	//	T d = sqrt(dx*dx + dy*dy + dz*dz);					
	//	T ax = 57.2957795*acos(dz / d);						// 180°/pi
	//	if (dz < 0.0)
	//		ax = -ax;
	//	T rx = -dy*dz;
	//	T ry = dx*dz;

	//	glPushMatrix();
	//	glTranslatef(x1, y1, z1);
	//	glRotatef(ax, rx, ry, 0.0);

	//	glutSolidCylinder(radius, d, subdivisions, 1);
	//	glPopMatrix();
	//}

	///render cylinder based on points from the top/bottom hat
	///@param C1 set of points from one of the hat
	void renderCylinder(std::vector< stim::vec3<T> > C1, std::vector< stim::vec3<T> > C2) {
		glBegin(GL_QUAD_STRIP);
		for (unsigned i = 0; i < C1.size(); i++) {			// for every point on the circle
			glVertex3f(C1[i][0], C1[i][1], C1[i][2]);
			glVertex3f(C2[i][0], C2[i][1], C2[i][2]);
		}	
		glEnd();
		//glFlush();
	}

	///render the vertex as sphere
	///@param x, y, z are the three coordinates of the center point
	///@param radius is the radius of the sphere
	///@param subdivisions is the slice/stride along/around z-axis
	void renderBall(T x, T y, T z, T radius, int subdivisions) {
		glPushMatrix();
		glTranslatef(x, y, z);
		glutSolidSphere(radius, subdivisions, subdivisions);
		glPopMatrix();
	}

	/// Render the network centerline as a series of line strips.
	/// glCenterline0 is for only one input
	void glCenterline0(){
		if (!glIsList(dlist)) {					//if dlist isn't a display list, create it
			dlist = glGenLists(1);				//generate a display list
			glNewList(dlist, GL_COMPILE);		//start a new display list
			for (unsigned e = 0; e < E.size(); e++) {				//for each edge in the network
				glBegin(GL_LINE_STRIP);
				for (unsigned p = 0; p < E[e].size(); p++) {			//for each point on that edge
					glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);		//set the vertex position based on the current point
					glTexCoord1f(0);									//set white color
				}
				glEnd();
			}
			glEndList();						//end the display list
		}
		glCallList(dlist);					// render the display list
	}

	void glCenterline(){

		if(!glIsList(dlist)){					//if dlist isn't a display list, create it
			dlist = glGenLists(1);				//generate a display list
			glNewList(dlist, GL_COMPILE);		//start a new display list
			for(unsigned e = 0; e < E.size(); e++){				//for each edge in the network
				//unsigned errormag_id = E[e].nmags() - 1;
				glBegin(GL_LINE_STRIP);
				for(unsigned p = 0; p < E[e].size(); p++){				//for each point on that edge
					glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);		//set the vertex position based on the current point
					glTexCoord1f(E[e].r(p));						//set the texture coordinate based on the specified magnitude index
				}
				glEnd();
			}
			glEndList();						//end the display list
		}		
		glCallList(dlist);						//render the display list
	}

	///render the GT network cylinder as series of tubes
	///@param dlist1 is the display list
	///@param map is the mapping relationship between two networks
	///@param colormap is the random generated color set for render
	void glCylinderGT(GLuint &dlist1, std::vector<unsigned> map, std::vector<T> colormap) {
		if (!glIsList(dlist1)) {								// if dlist1 isn't a display list, create it
			dlist1 = glGenLists(1);								// generate a display list
			glNewList(dlist1, GL_COMPILE);						// start a new display list
			for (unsigned e = 0; e < E.size(); e++) {			// for each edge in the network
				if (map[e] != unsigned(-1)) {
					glColor3f(colormap[e * 3 + 0], colormap[e * 3 + 1], colormap[e * 3 + 2]);
					for (unsigned p = 1; p < E[e].size(); p++) {// for each point on that edge
						stim::circle<T> C1 = E[e].circ(p - 1);
						stim::circle<T> C2 = E[e].circ(p);
						C1.set_R(10);								// scale the circle to the same
						C2.set_R(10);
						std::vector< stim::vec3<T> >Cp1 = C1.points(20);
						std::vector< stim::vec3<T> >Cp2 = C2.points(20);
						renderCylinder(Cp1, Cp2);
					}
				}
				else {
					glColor3f(1.0f, 1.0f, 1.0f);					// white color for the un-mapping edges
					for (unsigned p = 1; p < E[e].size(); p++) {// for each point on that edge
						stim::circle<T> C1 = E[e].circ(p - 1);
						stim::circle<T> C2 = E[e].circ(p);
						C1.set_R(10);								// scale the circle to the same
						C2.set_R(10);
						std::vector< stim::vec3<T> >Cp1 = C1.points(20);
						std::vector< stim::vec3<T> >Cp2 = C2.points(20);
						renderCylinder(Cp1, Cp2);
					}
				}
			}
			for (unsigned v = 0; v < V.size(); v++) {
				size_t num_edge = V[v].e[0].size() + V[v].e[1].size();
				if (num_edge > 1) {					// if it is the joint vertex
					glColor3f(0.3, 0.3, 0.3);		// gray color 
					renderBall(V[v][0], V[v][1], V[v][2], 20, 20);
				}
				else {								// if it is the terminal vertex
					glColor3f(0.6, 0.6, 0.6);		// more white gray
					renderBall(V[v][0], V[v][1], V[v][2], 20, 20);
				}
			}
			glEndList();
		}
		glCallList(dlist1);
	}

	///render the T network cylinder as series of tubes
	///@param dlist2 is the display list
	///@param map is the mapping relationship between two networks
	///@param colormap is the random generated color set for render
	void glCylinderT(GLuint &dlist2, std::vector<unsigned> map, std::vector<T> colormap) {
		if (!glIsList(dlist2)) {
			dlist2 = glGenLists(1);
			glNewList(dlist2, GL_COMPILE);
			for (unsigned e = 0; e < E.size(); e++) {				// for each edge in the network
				if (map[e] != unsigned(-1)) {
					glColor3f(colormap[map[e] * 3 + 0], colormap[map[e] * 3 + 1], colormap[map[e] * 3 + 2]);
					for (unsigned p = 0; p < E[e].size() - 1; p++) {// for each point on that edge
						stim::circle<T> C1 = E[e].circ(p);
						stim::circle<T> C2 = E[e].circ(p + 1);
						C1.set_R(10);								// scale the circle to the same
						C2.set_R(10);
						std::vector< stim::vec3<T> >Cp1 = C1.points(20);
						std::vector< stim::vec3<T> >Cp2 = C2.points(20);
						renderCylinder(Cp1, Cp2);
					}
				}
				else {
					glColor3f(1.0f, 1.0f, 1.0f);					// white color for the un-mapping edges
					for (unsigned p = 0; p < E[e].size() - 1; p++) {// for each point on that edge
						stim::circle<T> C1 = E[e].circ(p);
						stim::circle<T> C2 = E[e].circ(p + 1);
						C1.set_R(10);								// scale the circle to the same
						C2.set_R(10);
						std::vector< stim::vec3<T> >Cp1 = C1.points(20);
						std::vector< stim::vec3<T> >Cp2 = C2.points(20);
						renderCylinder(Cp1, Cp2);
					}
				}
			}
			for (unsigned v = 0; v < V.size(); v++) {
				size_t num_edge = V[v].e[0].size() + V[v].e[1].size();
				if (num_edge > 1) {					// if it is the joint vertex
					glColor3f(0.3, 0.3, 0.3);		// gray color 
					renderBall(V[v][0], V[v][1], V[v][2], 20, 20);
				}
				else {								// if it is the terminal vertex
					glColor3f(0.6, 0.6, 0.6);		// more white gray
					renderBall(V[v][0], V[v][1], V[v][2], 20, 20);
				}
			}
			glEndList();
		}
		glCallList(dlist2);
	}

	/// Render the GT network centerline as a series of line strips in random different color
	///@param dlist1 is the display list
	///@param map is the mapping relationship between two networks
	///@param colormap is the random generated color set for render
	void glRandColorCenterlineGT(GLuint &dlist1, std::vector<unsigned> map, std::vector<T> colormap) {
		if (!glIsList(dlist1)) {
			dlist1 = glGenLists(1);
			glNewList(dlist1, GL_COMPILE);
			for (unsigned e = 0; e < E.size(); e++) {
				if (map[e] != unsigned(-1)) {						// if it has corresponding edge in another network
					glColor3f(colormap[e * 3 + 0], colormap[e * 3 + 1], colormap[e * 3 + 2]);
					glBegin(GL_LINE_STRIP);
					for (unsigned p = 0; p < E[e].size(); p++) {
						glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
					}
					glEnd();
				}
				else {
					glColor3f(1.0, 1.0, 1.0);						// white color
					glBegin(GL_LINE_STRIP);
					for (unsigned p = 0; p < E[e].size(); p++) {
						glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
					}
					glEnd();
				}
			}
			glEndList();
		}
		glCallList(dlist1);
	}

	/// Render the T network centerline as a series of line strips in random different color
	///@param dlist2 is the display list
	///@param map is the mapping relationship between two networks
	///@param colormap is the random generated color set for render
	void glRandColorCenterlineT(GLuint &dlist2, std::vector<unsigned> map, std::vector<T> colormap) {
		if (!glIsList(dlist2)) {
			dlist2 = glGenLists(1);
			glNewList(dlist2, GL_COMPILE);
			for (unsigned e = 0; e < E.size(); e++) {
				if (map[e] != unsigned(-1)) {						// if it has corresponding edge in another network
					glColor3f(colormap[map[e] * 3 + 0], colormap[map[e] * 3 + 1], colormap[map[e] * 3 + 2]);
					glBegin(GL_LINE_STRIP);
					for (unsigned p = 0; p < E[e].size(); p++) {
						glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
					}
					glEnd();
				}
				else {
					glColor3f(1.0, 1.0, 1.0);						// white color
					glBegin(GL_LINE_STRIP);
					for (unsigned p = 0; p < E[e].size(); p++) {
						glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
					}
					glEnd();
				}
			}
			glEndList();
		}
		glCallList(dlist2);
	}

	//void glRandColorCenterlineGT(GLuint &dlist1, std::vector<unsigned> map, std::vector<T> colormap){
	//	if(!glIsList(dlist1)){
	//		dlist1 = glGenLists(1);
	//		glNewList(dlist1, GL_COMPILE);
	//		for(unsigned e = 0; e < E.size(); e++){
	//			if(map[e] != unsigned(-1)){
	//				glColor3f(colormap[e * 3 + 0], colormap[e * 3 + 1], colormap[e * 3 + 2]);
	//				glBegin(GL_LINE_STRIP);
	//				for(unsigned p = 0; p < E[e].size(); p++){
	//					glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
	//				}
	//				glEnd();
	//				for (unsigned p = 0; p < E[e].size() - 1; p++) {
	//					renderCylinder(E[e][p][0], E[e][p][1], E[e][p][2], E[e][p + 1][0], E[e][p + 1][1], E[e][p + 1][2], 10, 20);
	//				}
	//			}
	//			else{
	//				glColor3f(1.0, 1.0, 1.0);
	//				glBegin(GL_LINE_STRIP);
	//				for(unsigned p = 0; p < E[e].size(); p++){
	//					glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
	//				}
	//				glEnd();
	//			}
	//		}
	//		for (unsigned v = 0; v < V.size(); v++) {
	//			size_t num_edge = V[v].e[0].size() + V[v].e[1].size();
	//			if (num_edge > 1) {
	//				glColor3f(0.3, 0.3, 0.3);		// gray color for vertex
	//				renderBall(V[v][0], V[v][1], V[v][2], 20, 20);
	//			}
	//		}
	//		glEndList();
	//	}
	//	glCallList(dlist1);
	//}

	//void glRandColorCenterlineT(GLuint &dlist2, std::vector<unsigned> map, std::vector<T> colormap){
	//	if(!glIsList(dlist2)){
	//		dlist2 = glGenLists(1);
	//		glNewList(dlist2, GL_COMPILE);
	//		for(unsigned e = 0; e < E.size(); e++){
	//			if(map[e] != unsigned(-1)){
	//				glColor3f(colormap[map[e] * 3 + 0], colormap[map[e] * 3 + 1], colormap[map[e] * 3 + 2]);
	//				glBegin(GL_LINE_STRIP);
	//				for(unsigned p = 0; p < E[e].size(); p++){
	//					glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
	//				}
	//				glEnd();
	//				for (unsigned p = 0; p < E[e].size() - 1; p++) {
	//					renderCylinder(E[e][p][0], E[e][p][1], E[e][p][2], E[e][p + 1][0], E[e][p + 1][1], E[e][p + 1][2], 10, 20);
	//				}
	//			}
	//			else{
	//				glColor3f(1.0, 1.0, 1.0);
	//				glBegin(GL_LINE_STRIP);
	//				for(unsigned p = 0; p < E[e].size(); p++){
	//					glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
	//				}
	//				glEnd();
	//			}
	//		}
	//		for (unsigned v = 0; v < V.size(); v++) {
	//			size_t num_edge = V[v].e[0].size() + V[v].e[1].size();
	//			if (num_edge > 1) {
	//				glColor3f(0.3, 0.3, 0.3);		// gray color for vertex
	//				renderBall(V[v][0], V[v][1], V[v][2], 20, 20);
	//			}
	//		}
	//		glEndList();
	//	}
	//	glCallList(dlist2);
	//}

};		//end stim::gl_network class
};		//end stim namespace



#endif