better render

Jiaming Guo
1 parent 5bd6c411
Showing 3 changed files with 182 additions and 116 deletions Show diff stats
stim/biomodels/network.h
stim/visualization/gl_network.h
stim/visualization/swc.h
@@ -395,31 +395,40 @@ public:
 		stim::swc<T> S;										// create swc variable
 		S.load(filename);									// load the node information
 		S.create_tree();									// link those node according to their linking relationships as a tree
+		S.resample();
 		//NT.push_back(S.node[0].type);						// set the neuronal_type value to the first vertex in the network
 		std::vector<unsigned> id2vert;						// this list stores the SWC vertex ID associated with each network vertex
 		unsigned i[2];										// temporary, IDs associated with the first and last points
-		for (unsigned int l = 1; l < S.numL(); l++) {		// for every node starts from second one
+
+		for (unsigned int l = 0; l < S.numE(); l++) {		// for every edge
 			//NT.push_back(S.node[l].type);
-			stim::centerline<T> c3(2);						// every fiber contains two vertices
-			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
+
+			std::vector< stim::vec3<T> > c;
+			S.get_points(l, c);
+
+			stim::centerline<T> c3(c.size());				// new fiber
+			
+			for (unsigned j = 0; j < c.size(); j++)
+				c3[j] = c[j];								// copy the points
 			c3.update();									// update the L information
 			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);
+			S.get_radius(l, radius);
+
 			C3.copy_r(radius);
-			edge new_edge(C3);								// contains two points								
+			edge new_edge(C3);								// new edge	
+
+			//create an edge from the given centerline
+			unsigned int I = new_edge.size();				//calculate the number of points on the centerline
 			//get the first and last vertex IDs for the line
-			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
+			i[0] = S.E[l].front();
+			i[1] = S.E[l].back();
 			std::vector<unsigned>::iterator it;				//create an iterator for searching the id2vert array
 			unsigned it_idx;								//create an integer for the id2vert entry index
@@ -441,7 +450,7 @@ public:
 			it = find(id2vert.begin(), id2vert.end(), i[1]);	//look for the second ID
 			if (it == id2vert.end()) {							//if i[1] hasn't already been used
-				vertex new_vertex = new_edge[1];			//create a new vertex, assign it a position
+				vertex new_vertex = new_edge[I - 1];			//create a new vertex, assign it a position
 				new_vertex.e[1].push_back(E.size());			//add the current edge as incoming
 				new_edge.v[1] = V.size();                       //add the new vertex to the edge
 				V.push_back(new_vertex);						//add the new vertex to the vertex list
@@ -620,8 +629,6 @@ public:
 		kdt.create(c, n_data, MaxTreeLevels);
 		for(unsigned e = 0; e < R.E.size(); e++){			//for each edge in A
-			R.E[e].add_mag(0);								//add a new magnitude for the metric
-			size_t errormag_id = R.E[e].nmags() - 1;
 			size_t n = R.E[e].size();						// the number of points in current edge
 			T* query = new T[3 * n];
@@ -634,8 +641,8 @@ public:
 			kdt.cpu_search(queryPt, n, nnIdx, dists);			//find the distance between A and the current network
 			for(unsigned p = 0; p < R.E[e].size(); p++){
-				m1[p] = 1.0f - gaussianFunction((T)dists[p], sigma);		//calculate the metric value based on the distance
-				R.E[e].set_mag(errormag_id, p, m1[p]);					//set the error for the second point in the segment
+				m1[p] = 1.0f - gaussianFunction((T)dists[p], sigma);	//calculate the metric value based on the distance
+				R.E[e].set_r(p, m1[p]);					//set the error for the second point in the segment
 			}
 		}
 #endif
@@ -746,10 +753,8 @@ public:
 			edge_point_num[ee] = B.E[ee].size();
 		for(unsigned e = 0; e < A.E.size(); e++){			//for each edge in A
-			A.E[e].add_mag(0);								//add a new magnitude for the metric
-			size_t errormag_id = A.E[e].nmags() - 1;
-
-			size_t n = A.E[e].size();						// the number of edges in A
+			
+			size_t n = A.E[e].size();						//the number of edges in A
 			T* queryPt = new T[3 * n];
 			T* m1 = new T[n];
@@ -761,7 +766,7 @@ public:
 			for(unsigned p = 0; p < A.E[e].size(); p++){
 				m1[p] = 1.0f - gaussianFunction((T)dists[p], sigma);	//calculate the metric value based on the distance
-				A.E[e].set_mag(errormag_id, p, m1[p]);					//set the error for the second point in the segment
+				A.E[e].set_r(p, m1[p]);									//set the error for the second point in the segment
 				unsigned id = 0;																	//mapping edge's idx
 				size_t num = 0;																		//total number of points before #th edge
@@ -909,9 +914,8 @@ public:
 		for(unsigned e = 0; e < R.E.size(); e++){			// for each edge in A
 			T M;											// the sum of metrics of current edge
-			unsigned errormag_id = R.E[e].nmags() - 1;
 			for(unsigned p = 0; p < R.E[e].size(); p++)
-				M += A.E[e].m(p, errormag_id);
+				M += A.E[e].r(p);
 			M = M / A.E[e].size();
 			if(M > threshold)								
 				C[e] = (unsigned)-1;
@@ -115,6 +115,8 @@ public:
 		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
+			glEnable(GL_COLOR_MATERIAL);
+			glColorMaterial(GL_FRONT, GL_AMBIENT_AND_DIFFUSE);
 			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);
@@ -125,7 +127,7 @@ public:
 					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)
-						glTexCoord1f(E[e].r(p-1));
+						glTexCoord1f(E[e].r(p - 1));
 						glVertex3f(Cp1[i][0], Cp1[i][1], Cp1[i][2]);
 						glTexCoord1f(E[e].r(p));
 						glVertex3f(Cp2[i][0], Cp2[i][1], Cp2[i][2]);
@@ -150,77 +152,67 @@ public:
 		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 glRandColorCylinder1(GLuint &dlist1, std::vector<unsigned> map, std::vector<T> colormap, float sigma, float radius) {
+	///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(dlist1, 1);
-		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(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);
-					}
-				}
-				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(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);
+		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 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], 3*radius, 20);
+			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 {								// 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], 3*radius, 20);
+				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(dlist1);
+		glCallList(dlist+4);
 	}
-	///render the T network cylinder as series of tubes
-	///@param dlist2 is the display list
+	///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 glRandColorCylinder2(GLuint &dlist2, std::vector<unsigned> map, std::vector<T> colormap, float sigma, float radius) {
+	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(dlist2, 1);
-		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 (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
+			glEnable(GL_COLOR_MATERIAL);
+			glColorMaterial(GL_FRONT, GL_AMBIENT_AND_DIFFUSE);
+			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 = 1; p < E[e].size(); p++) {	// for each point on that edge(except last one)
+					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);						// scale the circle to the same
+						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);
@@ -252,52 +244,27 @@ public:
 				}
 			}
 			glEndList();
+			glDisable(GL_COLOR_MATERIAL);
 		}
-		glCallList(dlist2);
+		glCallList(dlist+2);
 	}
-	/// Render the GT network centerline as a series of line strips in random different color
+	///render the network centerline as a series of line strips in random different color
+	///@param I is a indicator: 0 -> GT, 1 -> T
 	///@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 glRandColorCenterline1(GLuint &dlist1, std::vector<unsigned> map, std::vector<T> colormap) {
+	void glRandColorCenterline(int I, 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);
-	}
+				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]);
-	/// 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 glRandColorCenterline2(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]);
@@ -305,7 +272,7 @@ public:
 					glEnd();
 				}
 				else {
-					glColor3f(1.0, 1.0, 1.0);						// white color
+					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]);
@@ -315,7 +282,7 @@ public:
 			}
 			glEndList();
 		}
-		glCallList(dlist2);
+		glCallList(dlist1);
 	}
 	//void glRandColorCenterlineGT(GLuint &dlist1, std::vector<unsigned> map, std::vector<T> colormap){
@@ -86,6 +86,8 @@ namespace stim {
 	class swc {
 	public:
 		std::vector< typename swc_tree::swc_node<T> > node;
+		std::vector< std::vector<int> > E;			// list of nodes
+
 		swc() {};									// default constructor
 		// load the swc as tree nodes
@@ -170,11 +172,104 @@ namespace stim {
 			}
 		}
+		// create a new edge
+		void create_up(std::vector<int>& edge,swc_tree::swc_node<T> cur_node, int target) {
+			while (cur_node.parent_idx != target) {
+				edge.push_back(cur_node.idx);
+				cur_node = node[cur_node.parent_idx - 1];			// move to next node
+			}
+			edge.push_back(cur_node.idx);							// push back the start/end vertex of current edge
+											
+			std::reverse(edge.begin(), edge.end());					// follow the original flow direction
+		}
+		void create_down(std::vector<int>& edge, swc_tree::swc_node<T> cur_node, int j) {
+			while (cur_node.son_idx.size() != 0) {
+				edge.push_back(cur_node.idx);
+				if (cur_node.son_idx.size() > 1)
+					cur_node = node[cur_node.son_idx[j] - 1];		// move to next node
+				else
+					cur_node = node[cur_node.son_idx[0] - 1];
+			}
+			edge.push_back(cur_node.idx);							// push back the start/end vertex of current edge
+		}
+
+		// resample the tree-like SWC
+		void resample() {
+			unsigned n = node.size();
+
+			std::vector<int> joint_node;
+			for (unsigned i = 1; i < n; i++) {			// search all nodes(except the first one) to find joint nodes
+				if (node[i].son_idx.size() > 1) {
+					joint_node.push_back(node[i].idx);	// store the original index
+				}
+			}
+
+			std::vector<int>  new_edge;					// new edge in the network
+			
+			n = joint_node.size();
+			
+			for (unsigned i = 0; i < n; i++) {			// for every joint nodes
+				std::vector<int>  new_edge;									// new edge in the network
+				swc_tree::swc_node<T> cur_node = node[joint_node[i] - 1];	// store current node
+				
+				// go up
+				swc_tree::swc_node<T> tmp_node = node[cur_node.parent_idx - 1];	
+				while (tmp_node.parent_idx != -1 && tmp_node.son_idx.size() == 1) {
+					tmp_node = node[tmp_node.parent_idx - 1];
+				}
+				int target = tmp_node.parent_idx;							// store the go-up target
+				create_up(new_edge, cur_node, target);
+				E.push_back(new_edge);
+				new_edge.clear();
+
+				// go down
+				unsigned t = cur_node.son_idx.size();
+				for (unsigned j = 0; j < t; j++) {							// for every son node
+					tmp_node = node[cur_node.son_idx[j] - 1];				// move down
+					while (tmp_node.son_idx.size() == 1) {
+						tmp_node = node[tmp_node.son_idx[0] - 1];
+					}
+					if (tmp_node.son_idx.size() == 0) {
+						create_down(new_edge, cur_node, j);
+						E.push_back(new_edge);
+						new_edge.clear();
+					}
+				}
+			}
+		}
+
+		// get points in one edge
+		void get_points(int e, std::vector< stim::vec3<T> >& V) {
+			V.resize(E[e].size());
+			for (unsigned i = 0; i < E[e].size(); i++) {
+				unsigned id = E[e][i] - 1;
+
+				V[i][0] = node[id].point[0];
+				V[i][1] = node[id].point[1];
+				V[i][2] = node[id].point[2];
+			}
+		}
+
+		// get radius information in one edge
+		void get_radius(int e, std::vector<T>& radius) {
+			radius.resize(E[e].size());
+			for (unsigned i = 0; i < E[e].size(); i++) {
+				unsigned id = E[e][i] - 1;
+
+				radius[i] = node[id].radius;
+			}
+		}
+
 		// return the number of point in swc
-		unsigned int numL() {
+		unsigned int numP() {
 			return node.size();
 		}
+		// return the number of edges in swc after resample
+		unsigned int numE() {
+			return E.size();
+		}
+
 	};
 }		// end of namespace stim