Merge branch 'JACK' into 'master'

fix minor error for loading and rendering networks from swc files See merge request !24

Merge branch 'JACK' into 'master'
fix minor error for loading and rendering networks from swc files See merge request !24
David Mayerich
2 parents 8705dfb0 58674036
Showing 4 changed files with 192 additions and 54 deletions Show diff stats
stim/biomodels/network.h
stim/visualization/cylinder.h
stim/visualization/gl_network.h
stim/visualization/swc.h
@@ -390,6 +390,7 @@ public:
 		}
 	}
+	//load a network from an SWC file
 	void load_swc(std::string filename) {
 		stim::swc<T> S;										// create swc variable
 		S.load(filename);									// load the node information
@@ -401,16 +402,23 @@ public:
 		for (unsigned int l = 1; l < S.numL(); l++) {		// for every node starts from second one
 			NT.push_back(S.node[l].type);
 			stim::centerline<T> c3(2);						// every fiber contains two vertices
-			int id = S.node[l].parent_idx - 1;
-			c3[0] = S.node[id].point;						// store the begin vertex
+			int p_idx = S.node[l].parent_idx - 1;			// get the parent node loop idx of current node
+			c3[0] = S.node[p_idx].point;					// store the begin vertex
 			c3[1] = S.node[l].point;						// store the end vertex
 			c3.update();									// update the L information
-
-			edge new_edge(c3);								// contains two points								
-
+			
+			stim::cylinder<T> C3(c3);						// create a new cylinder in order to copy the origin radius information
+			// upadate the R information
+			std::vector<T> radius;
+			radius.push_back(S.node[p_idx].radius);
+			radius.push_back(S.node[l].radius);
+			C3.copy_r(radius);
+
+			edge new_edge(C3);								// contains two points								
+			
 			//get the first and last vertex IDs for the line
-			i[0] = S.node[id].idx;							//get the SWC ID for the start point
+			i[0] = S.node[p_idx].idx;						//get the SWC ID for the start point
 			i[1] = S.node[l].idx;							//get the SWC ID for the end point
 			std::vector<unsigned>::iterator it;				//create an iterator for searching the id2vert array
@@ -471,7 +479,7 @@ public:
 	stim::network<T> resample(T spacing){
 		stim::network<T> n;								//create a new network that will be an exact copy, with resampled fibers
 		n.V = V;									//copy all vertices
-
+		n.NT = NT;										//copy all the neuronal type information
 		n.E.resize(edges());								//allocate space for the edge list
 		//copy all fibers, resampling them in the process
@@ -54,6 +54,12 @@ public:
 		//add_mag(r);
 	}
+	//copy the original radius
+	void copy_r(std::vector<T> radius) {
+		for (unsigned i = 0; i < radius.size(); i++)
+			R[i] = radius[i];
+	}
+
 	///Returns magnitude i at the given length into the fiber (based on the pvalue).
 	///Interpolates the position along the line.
 	///@param l: the location of the in the cylinder.
@@ -175,20 +181,23 @@ public:
 	/// Integrates a magnitude value along the cylinder.
 	/// @param m is the magnitude value to be integrated (this is usually the radius)
 	T integrate() {
+		T sum = 0;							//initialize the integral to zero
+		if (L.size() == 1)
+			return sum;
+		else {					
+			T m0, m1;						//allocate space for both magnitudes in a single segment
+			m0 = R[0];					//initialize the first point and magnitude to the first point in the cylinder
+			T len = L[1];					//allocate space for the segment length
-		T sum = 0;						//initialize the integral to zero
-		T m0, m1;						//allocate space for both magnitudes in a single segment
-		m0 = R[0];					//initialize the first point and magnitude to the first point in the cylinder
-		T len = L[1];					//allocate space for the segment length
-
-											
-		for (unsigned p = 1; p < size(); p++) {				//for every consecutive point in the cylinder
-			m1 = R[p];
-			if (p > 1) len = (L[p] - L[p - 1]);		//calculate the segment length using the L array
-			sum += (m0 + m1) / (T)2.0 * len;				//add the average magnitude, weighted by the segment length
-			m0 = m1;									//move to the next segment by shifting points
+
+			for (unsigned p = 1; p < size(); p++) {				//for every consecutive point in the cylinder
+				m1 = R[p];
+				if (p > 1) len = (L[p] - L[p - 1]);		//calculate the segment length using the L array
+				sum += (m0 + m1) / (T)2.0 * len;				//add the average magnitude, weighted by the segment length
+				m0 = m1;									//move to the next segment by shifting points
+			}
+			return sum;			//return the integral
 		}
-		return sum;			//return the integral
 	}
 	/// Resamples the cylinder to provide a maximum distance of "spacing" between centerline points. All current
@@ -43,28 +43,6 @@ public:
 		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) {
@@ -114,7 +92,7 @@ public:
 			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];
+				int type = NT[e];					// get the neuronal type
 				switch (type) {
 				case 0:
 					glColor3f(1.0f, 1.0f, 1.0f);	// white for undefined
@@ -187,8 +165,8 @@ public:
 		glCallList(dlist);					// render the display list
 	}
-	///render the network centerline as a series of line strips
-	///colors is based on metric values
+	///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
@@ -208,11 +186,39 @@ public:
 		glCallList(dlist);						//render the display list
 	}
+	///render the network cylinder as a series of tubes
+	///colors are based on metric values 
+	void glCylinder() {
+		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(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);
+					glBegin(GL_QUAD_STRIP);
+					for (unsigned i = 0; i < Cp1.size(); i++) {		// for every point on 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
+			}
+			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) {
+	void glRandColorCylinder1(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
@@ -258,11 +264,57 @@ public:
 		glCallList(dlist1);
 	}
+	void glRandColorCylinder1_swc(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(0.5);							// scale the circle to the same
+						C2.set_R(0.5);
+						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(0.5);								// scale the circle to the same
+						C2.set_R(0.5);
+						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], 1, 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], 1, 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) {
+	void glRandColorCylinder2(GLuint &dlist2, std::vector<unsigned> map, std::vector<T> colormap) {
 		if (!glIsList(dlist2)) {
 			dlist2 = glGenLists(1);
 			glNewList(dlist2, GL_COMPILE);
@@ -308,11 +360,57 @@ public:
 		glCallList(dlist2);
 	}
+	void glRandColorCylinder2_swc(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(0.5);								// scale the circle to the same
+						C2.set_R(0.5);
+						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(0.5);								// scale the circle to the same
+						C2.set_R(0.5);
+						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], 1, 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], 1, 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) {
+	void glRandColorCenterline1(GLuint &dlist1, std::vector<unsigned> map, std::vector<T> colormap) {
 		if (!glIsList(dlist1)) {
 			dlist1 = glGenLists(1);
 			glNewList(dlist1, GL_COMPILE);
@@ -343,7 +441,7 @@ public:
 	///@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) {
+	void glRandColorCenterline2(GLuint &dlist2, std::vector<unsigned> map, std::vector<T> colormap) {
 		if (!glIsList(dlist2)) {
 			dlist2 = glGenLists(1);
 			glNewList(dlist2, GL_COMPILE);
@@ -444,6 +542,29 @@ public:
 	//	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();
+	//}
+
 };		//end stim::gl_network class
 };		//end stim namespace
@@ -20,7 +20,7 @@ namespace stim {
 			enum node_information { INTEGER_LABEL, NEURONAL_TYPE, X_COORDINATE, Y_COORDINATE, Z_COORDINATE, RADIUS, PARENT_LABEL }; 
 		public:
-			int idx;							// the index of current node start from 1(ambiguity)
+			int idx;							// the index of current node start from 1(original index from the file)
 			neuronal_type type;					// the type of neuronal segmemt
 			stim::vec3<T> point;				// the point coordinates 
 			T radius;							// the radius at current node
@@ -44,7 +44,7 @@ namespace stim {
 				// for each information contained in the node point we got
 				for (unsigned int i = 0; i < p.size(); i++) {
 					std::stringstream ss(p[i]);	// create a stringstream object for casting
-					
+
 					// store different information
 					switch (i) {
 					case INTEGER_LABEL:
@@ -103,8 +103,8 @@ namespace stim {
 			std::string line;
 			// skip comment
 			while (getline(infile, line)) {
-				if ('#' == line[0])			// if it is comment line
-					continue;
+				if ('#' == line[0] || line.empty())			// if it is comment line or empty line
+					continue;								// keep read
 				else
 					break;
 			}
@@ -144,7 +144,7 @@ namespace stim {
 			std::string line;
 			while (getline(infile, line)) {
-				if ('#' == line[0]) {
+				if ('#' == line[0] || line.empty()) {
 					std::cout << line << std::endl;			// print the comment line by line
 				}
 				else
@@ -165,7 +165,7 @@ namespace stim {
 				if (node[i].parent_idx != node[i - 1].parent_idx)
 					cur_level = node[node[i].parent_idx - 1].level + 1;
 				node[i].level = cur_level;
-				int tmp_parent_idx = node[i].parent_idx - 1;	// get the parent node loop idx of current node
+				int tmp_parent_idx = node[i].parent_idx - 1;		// get the parent node loop idx of current node
 				node[tmp_parent_idx].son_idx.push_back(i + 1);		// son_idx stores the real idx = loop idx + 1
 			}
 		}