segmentation / netmets

#ifndef STIM_GL_NETWORK
#define STIM_GL_NETWORK

#include <GL/glut.h>
#include "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
	}

	///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, stim::vec3<T> P1, stim::vec3<T> P2) {
		glBegin(GL_QUAD_STRIP);
		for (unsigned i = 0; i < C1.size(); i++) {				// for every point on the circle
			stim::vec3<T> n1 = C1[i] - P1;						// set normal vector for every vertex on the quads
			stim::vec3<T> n2 = C2[i] - P2;
			n1 = n1.norm();
			n2 = n2.norm();
			glNormal3f(n1[0], n1[1], n1[2]);
			glVertex3f(C1[i][0], C1[i][1], C1[i][2]);
			glNormal3f(n2[0], n2[1], n2[2]);
			glVertex3f(C2[i][0], C2[i][1], C2[i][2]);											
		}	
		glEnd();
	}

	///render the vertex as sphere based on glut build-in function
	///@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 vertex as sphere based on transformation
	///@param x, y, z are the three coordinates of the center point
	///@param radius is the radius of the sphere
	///@param slice is the number of subdivisions around the z-axis
	///@param stack is the number of subdivisions along the z-axis
	void renderBall(T x, T y, T z, T radius, T slice, T stack) {
		T step_z = stim::PI / slice;					// step angle along z-axis
		T step_xy = 2 * stim::PI / stack;				// step angle in xy-plane
		T xx[4], yy[4], zz[4];							// store coordinates

		T angle_z = 0.0;								// start angle
		T angle_xy = 0.0;

		glBegin(GL_QUADS);
		for (unsigned i = 0; i < slice; i++) {			// around the z-axis
			angle_z = i * step_z;						// step step_z each time

			for (unsigned j = 0; j < stack; j++) {		// along the z-axis
				angle_xy = j * step_xy;					// step step_xy each time, draw floor by floor

				xx[0] = radius * std::sin(angle_z) * std::cos(angle_xy);	// four vertices
				yy[0] = radius * std::sin(angle_z) * std::sin(angle_xy);
				zz[0] = radius * std::cos(angle_z);

				xx[1] = radius * std::sin(angle_z + step_z) * std::cos(angle_xy);
				yy[1] = radius * std::sin(angle_z + step_z) * std::sin(angle_xy);
				zz[1] = radius * std::cos(angle_z + step_z);

				xx[2] = radius * std::sin(angle_z + step_z) * std::cos(angle_xy + step_xy);
				yy[2] = radius * std::sin(angle_z + step_z) * std::sin(angle_xy + step_xy);
				zz[2] = radius * std::cos(angle_z + step_z);

				xx[3] = radius * std::sin(angle_z) * std::cos(angle_xy + step_xy);
				yy[3] = radius * std::sin(angle_z) * std::sin(angle_xy + step_xy);
				zz[3] = radius * std::cos(angle_z);

				for (unsigned k = 0; k < 4; k++) {
					glVertex3f(x + xx[k], y + yy[k], z + zz[k]);			// draw the floor plane
				}
			}
		}
		glEnd();
	}

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

	///render the network centerline as a series of line strips(when loading at least two networks, otherwise using glCenterline0())
	///colors are based on metric values
	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 network cylinder as a series of tubes(when only one network loaded) 
	void glCylinder0(T scale = 1.0f, bool undo = false) {

		float r1, r2;
		if (undo == true)
			glDeleteLists(dlist, 1);								// delete display list 
		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
				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);
					r1 = E[e].r(p - 1);
					r2 = E[e].r(p);
					C1.set_R(scale * r1);								// re-scale the circle to the same
					C2.set_R(scale * r2);
					std::vector< stim::vec3<T> > Cp1 = C1.glpoints(20);	// get 20 points on the circle plane
					std::vector< stim::vec3<T> > Cp2 = C2.glpoints(20);
					glBegin(GL_QUAD_STRIP);
					for (unsigned i = 0; i < Cp1.size(); i++) {
						glVertex3f(Cp1[i][0], Cp1[i][1], Cp1[i][2]);
						glVertex3f(Cp2[i][0], Cp2[i][1], Cp2[i][2]);
					}
					glEnd();
				}													// set the texture coordinate based on the specified magnitude index
			}
			for (unsigned n = 0; n < V.size(); n++) {
				for (unsigned i = 0; i < E.size(); i++) {
					if (E[i].v[0] == n) {
						r1 = E[i].r(0) * scale;
						break;
					}
					else if (E[i].v[1] == n) {
						r1 = E[i].r(E[i].size() - 1) * scale;
						break;
					}
				}
				renderBall(V[n][0], V[n][1], V[n][2], r1, 20);
			}
			glEndList();						// end the display list
		}
		glCallList(dlist);						// render the display list
	}

	///render the network cylinder as a series of tubes
	///colors are based on metric values 
	void glCylinder(float sigma, float radius) {

		if (radius != sigma)					// if render radius was changed by user, create a new display list
			glDeleteLists(dlist, 1);
		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
				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(2*radius);								// re-scale the circle to the same
					C2.set_R(2*radius);
					std::vector< stim::vec3<T> > Cp1 = C1.glpoints(20);// get 20 points on the circle plane
					std::vector< stim::vec3<T> > Cp2 = C2.glpoints(20);
					glBegin(GL_QUAD_STRIP);
					for (unsigned i = 0; i < Cp1.size(); i++) {
						stim::vec3<T> n1 = Cp1[i] - E[e][p - 1];	// set normal vector for every vertex on the quads
						stim::vec3<T> n2 = Cp2[i] - E[e][p];
						n1 = n1.norm();
						n2 = n2.norm();

						glNormal3f(n1[0], n1[1], n1[2]);
						glTexCoord1f(E[e].r(p - 1));
						glVertex3f(Cp1[i][0], Cp1[i][1], Cp1[i][2]);
						glNormal3f(n2[0], n2[1], n2[2]);
						glTexCoord1f(E[e].r(p));
						glVertex3f(Cp2[i][0], Cp2[i][1], Cp2[i][2]);
					}
					glEnd();
				}													// set the texture coordinate based on the specified magnitude index
			}
			for (unsigned n = 0; n < V.size(); n++) {
				size_t num = V[n].e[0].size();					// store the number of outgoing edge of that vertex
				if (num != 0) {									// if it has outgoing edge
					unsigned idx = V[n].e[0][0];				// find the index of first outgoing edge of that vertex
					glTexCoord1f(E[idx].r(0));					// bind the texture as metric of first point on that edge
				}					
				else {
					unsigned idx = V[n].e[1][0];				// find the index of first incoming edge of that vertex
					glTexCoord1f(E[idx].r(E[idx].size() - 1));	// bind the texture as metric of last point on that edge
				}
				renderBall(V[n][0], V[n][1], V[n][2], 2*radius, 20);
			}
			glEndList();						// end the display list
		}
		glCallList(dlist);						// render the display list
	}

	///render a T as a adjoint network of GT in transparancy(darkgreen, overlaid)
	void glAdjointCylinder(float sigma, float radius) {
		
		if (radius != sigma)					// if render radius was changed by user, create a new display list
			glDeleteLists(dlist + 4, 1);
		if (!glIsList(dlist + 4)) {
			glNewList(dlist + 4, GL_COMPILE);
			for (unsigned e = 0; e < E.size(); e++) {				// for each edge in the network
				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(2 * radius);							// scale the circle to the same
					C2.set_R(2 * radius);
					std::vector< stim::vec3<T> >Cp1 = C1.glpoints(20);
					std::vector< stim::vec3<T> >Cp2 = C2.glpoints(20);
					glBegin(GL_QUAD_STRIP);
					for (unsigned i = 0; i < Cp1.size(); i++) {		// for every point on the circle(+1 means closing the circle)
						glVertex3f(Cp1[i][0], Cp1[i][1], Cp1[i][2]);
						glVertex3f(Cp2[i][0], Cp2[i][1], Cp2[i][2]);
					}
					glEnd();
				}								// set the texture coordinate based on the specified magnitude index
			}
			for (unsigned n = 0; n < V.size(); n++) {
				size_t num = V[n].e[0].size();					// store the number of outgoing edge of that vertex
				if (num != 0) {									// if it has outgoing edge
					unsigned idx = V[n].e[0][0];				// find the index of first outgoing edge of that vertex 
				}
				else {
					unsigned idx = V[n].e[1][0];				// find the index of first incoming edge of that vertex
				}
				renderBall(V[n][0], V[n][1], V[n][2], 2 * radius, 20);
			}
			glEndList();
		}
		glCallList(dlist + 4);
	}

	///render the network cylinder as series of tubes
	///@param I is a indicator: 0 -> GT, 1 -> T
	///@param map is the mapping relationship between two networks
	///@param colormap is the random generated color set for render
	void glRandColorCylinder(int I, std::vector<unsigned> map, std::vector<T> colormap, float sigma, float radius) {
		
		if (radius != sigma)									// if render radius was changed by user, create a new display list
			glDeleteLists(dlist + 2, 1);
		if (!glIsList(dlist + 2)) {								// if dlist isn't a display list, create it
			glNewList(dlist + 2, 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)) {
					if (I == 0) {								// if it is to render GT
						glColor3f(colormap[e * 3 + 0], colormap[e * 3 + 1], colormap[e * 3 + 2]);
					}									
					else {										// if it is to render T
						glColor3f(colormap[map[e] * 3 + 0], colormap[map[e] * 3 + 1], colormap[map[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(2*radius);						// re-scale the circle to the same
						C2.set_R(2*radius);
						std::vector< stim::vec3<T> >Cp1 = C1.glpoints(20);
						std::vector< stim::vec3<T> >Cp2 = C2.glpoints(20);
						renderCylinder(Cp1, Cp2, E[e][p - 1], E[e][p]);
					}
				}
				else {
					glColor3f(1.0, 1.0, 1.0);					// 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(2*radius);						// scale the circle to the same
						C2.set_R(2*radius);
						std::vector< stim::vec3<T> >Cp1 = C1.glpoints(20);
						std::vector< stim::vec3<T> >Cp2 = C2.glpoints(20);
						renderCylinder(Cp1, Cp2, E[e][p - 1], E[e][p]);
					}
				}
			}
			for (unsigned n = 0; n < V.size(); n++) {
				size_t num_edge = V[n].e[0].size() + V[n].e[1].size();
				if (num_edge > 1) {					// if it is the joint vertex
					glColor3f(0.3, 0.3, 0.3);		// gray
					renderBall(V[n][0], V[n][1], V[n][2], 3*radius, 20);
				}
				else {								// if it is the terminal vertex
					glColor3f(0.6, 0.6, 0.6);		// more white gray
					renderBall(V[n][0], V[n][1], V[n][2], 3*radius, 20);
				}
			}
			glEndList();
		}
		glCallList(dlist + 2);
	}

	///render the network centerline as a series of line strips in random different color
	///@param I is a indicator: 0 -> GT, 1 -> T
	///@param map is the mapping relationship between two networks
	///@param colormap is the random generated color set for render
	void glRandColorCenterline(int I, std::vector<unsigned> map, std::vector<T> colormap) {
		if (!glIsList(dlist + 2)) {
			glNewList(dlist + 2, GL_COMPILE);
			for (unsigned e = 0; e < E.size(); e++) {
				if (map[e] != unsigned(-1)) {					// if it has corresponding edge in another network
					if (I == 0)									// if it is to render GT
						glColor3f(colormap[e * 3 + 0], colormap[e * 3 + 1], colormap[e * 3 + 2]);
					else										// if it is to render T
						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 for the un-mapping edges
					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(dlist + 2);
	}

	void glAdjointCenterline() {
		if (!glIsList(dlist + 4)) {
			glNewList(dlist + 4, GL_COMPILE);
			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(E[e].r(p));							//set the texture coordinate based on the specified magnitude index
				}
				glEnd();
			}
			glEndList();
		}
		glCallList(dlist + 4);
	}

	// highlight the difference part
	void glDifferenceCylinder(int I, std::vector<unsigned> map, std::vector<T> colormap, float sigma, float radius) {

		if (radius != sigma)									// if render radius was changed by user, create a new display list
			glDeleteLists(dlist + 6, 1);
		if (!glIsList(dlist + 6)) {								// if dlist isn't a display list, create it
			glNewList(dlist + 6, 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)) {
					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

					if (I == 0) {								// if it is to render GT
						glColor4f(colormap[e * 3 + 0], colormap[e * 3 + 1], colormap[e * 3 + 2], 0.1);
					}
					else {										// if it is to render T
						glColor4f(colormap[map[e] * 3 + 0], colormap[map[e] * 3 + 1], colormap[map[e] * 3 + 2], 0.1);
					}

					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(2 * radius);						// re-scale the circle to the same
						C2.set_R(2 * radius);
						std::vector< stim::vec3<T> >Cp1 = C1.glpoints(20);
						std::vector< stim::vec3<T> >Cp2 = C2.glpoints(20);
						renderCylinder(Cp1, Cp2, E[e][p - 1], E[e][p]);
					}
					glDisable(GL_BLEND);
					glEnable(GL_DEPTH_TEST);
				}
				else {
					glColor3f(1.0, 1.0, 1.0);					// 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(2 * radius);						// scale the circle to the same
						C2.set_R(2 * radius);
						std::vector< stim::vec3<T> >Cp1 = C1.glpoints(20);
						std::vector< stim::vec3<T> >Cp2 = C2.glpoints(20);
						renderCylinder(Cp1, Cp2, E[e][p - 1], E[e][p]);
					}
				}
			}
			for (unsigned n = 0; n < V.size(); n++) {

				size_t num_edge = V[n].e[0].size() + V[n].e[1].size();
				if (num_edge > 1) {					// if it is the joint vertex
					glColor4f(0.3, 0.3, 0.3, 0.1);		// gray
					renderBall(V[n][0], V[n][1], V[n][2], 3 * radius, 20);
				}
				else {								// if it is the terminal vertex
					glColor4f(0.6, 0.6, 0.6, 0.1);		// more white gray
					renderBall(V[n][0], V[n][1], V[n][2], 3 * radius, 20);
				}
			}
			glEndList();
		}
		glCallList(dlist + 6);
	}

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


	//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 the network centerline from swc file as a series of strips in different colors based on the neuronal type
	/// glCenterline0_swc is for only one input
	/*void glCenterline0_swc() {
	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++) {
	int type = NT[e];					// get the neuronal type
	switch (type) {
	case 0:
	glColor3f(1.0f, 1.0f, 1.0f);	// white for undefined
	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();
	break;
	case 1:
	glColor3f(1.0f, 0.0f, 0.0f);	// red for soma
	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();
	break;
	case 2:
	glColor3f(1.0f, 0.5f, 0.0f);	// orange for axon
	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();
	break;
	case 3:
	glColor3f(1.0f, 1.0f, 0.0f);	// yellow for undefined
	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();
	break;
	case 4:
	glColor3f(0.0f, 1.0f, 0.0f);	// green for undefined
	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();
	break;
	case 5:
	glColor3f(0.0f, 1.0f, 1.0f);	// verdant for undefined
	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();
	break;
	case 6:
	glColor3f(0.0f, 0.0f, 1.0f);	// blue for undefined
	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();
	break;
	case 7:
	glColor3f(0.5f, 0.0f, 1.0f);	// purple for undefined
	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();
	break;
	}
	}
	glEndList();						//end the display list
	}
	glCallList(dlist);					// render the display list
	}*/

	///render the network cylinder as a series of tubes
	///colors are based on metric values 
	//void glCylinder(float sigma) {
	//	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
	//			for (unsigned p = 1; p < E[e].size() - 1; 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(2.5*sigma);							// scale the circle to the same
	//				C2.set_R(2.5*sigma);
	//				std::vector< stim::vec3<T> >Cp1 = C1.glpoints(20);
	//				std::vector< stim::vec3<T> >Cp2 = C2.glpoints(20);
	//				glBegin(GL_QUAD_STRIP);
	//				for (unsigned i = 0; i < Cp1.size(); i++) {		// for every point on the circle(+1 means closing the circle)
	//					glVertex3f(Cp1[i][0], Cp1[i][1], Cp1[i][2]);
	//					glVertex3f(Cp2[i][0], Cp2[i][1], Cp2[i][2]);
	//					glTexCoord1f(E[e].r(p));
	//				}
	//				glEnd();
	//			}								//set the texture coordinate based on the specified magnitude index
	//		}
	//		for (unsigned n = 0; n < V.size(); n++) {
	//			size_t num = V[n].e[0].size();					//store the number of outgoing edge of that vertex
	//			if (num != 0) {									//if it has outgoing edge
	//				unsigned idx = V[n].e[0][0];				//find the index of first outgoing edge of that vertex 
	//				glTexCoord1f(E[idx].r(0));					//bind the texture as metric of first point on that edge
	//			}
	//			else {
	//				unsigned idx = V[n].e[1][0];				//find the index of first incoming edge of that vertex
	//				glTexCoord1f(E[idx].r(E[idx].size() - 1));	//bind the texture as metric of last point on that edge
	//			}
	//			renderBall(V[n][0], V[n][1], V[n][2], 2.5*sigma, 20);
	//		}
	//		glEndList();						//end the display list
	//	}
	//	glCallList(dlist);						//render the display list
	//}

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



#endif