Merge branch 'JACK' into 'master'

add function to render the cylinders perfectly See merge request !22

Merge branch 'JACK' into 'master'
add function to render the cylinders perfectly See merge request !22
David Mayerich
2 parents 83c3121c 6b317c71
Showing 4 changed files with 248 additions and 65 deletions Show diff stats
stim/biomodels/network.h
stim/math/circle.h
stim/visualization/cylinder.h
stim/visualization/gl_network.h
@@ -46,7 +46,6 @@ __global__ void d_findedge(size_t* I, size_t n, unsigned* R, size_t* E, size_t n
  
 //hard-coded factor
 int threshold_fac = 10;
-float metric_fac = 0.6f;							// might be related to the distance field
  
 namespace stim{
 /** This is the a class that interfaces with gl_spider in order to store the currently
@@ -579,7 +578,7 @@ public:
 	/// @param B is the corresponding mapping network
 	/// @param sigma is the user-defined tolerance value - smaller values provide a stricter comparison
 	/// @param device is the device that user want to use
-	void split(stim::network<T> A, stim::network<T> B, float sigma, int device){
+	void split(stim::network<T> A, stim::network<T> B, float sigma, int device, float threshold){
  
 		T *c;						                 	
 		size_t n_data = B.total_points();          				
@@ -728,7 +727,7 @@ public:
 			T G = 0;																					// test to see whether it has its nearest neighbor
 			for(unsigned i = 0; i < E[e].size(); i++)
 				G += E[e].r(i);																			// won't split special edges
-			if(G / E[e].size() > metric_fac)															// should based on the color map
+			if(G / E[e].size() > threshold)															// should based on the color map
 				id = E[e].size() - 1;																	// set split idx to outgoing direction vertex
  
 			std::vector<edge> tmpe;
@@ -761,7 +760,7 @@ public:
 	/// @param B is the network that the current network is going to map to
 	/// @param C is the mapping relationship: C[e1] = _e1 means e1 edge in current network is mapping the _e1 edge in B
 	/// @param device is the device that user want to use
-	void mapping(stim::network<T> B, std::vector<unsigned> &C, int device){
+	void mapping(stim::network<T> B, std::vector<unsigned> &C, int device, float threshold){
 		stim::network<T> A;								//generate a network storing the result of the comparison
 		A = (*this);
  
@@ -803,7 +802,7 @@ public:
 			for(unsigned p = 0; p < A.E[e].size(); p++)
 				M += A.E[e].r(p);
 			M = M / A.E[e].size();
-			if(M > metric_fac)
+			if(M > threshold)
 				C[e] = (unsigned)-1;						//set the nearest edge of impossibly mapping edges to maximum of unsigned
 			else{
 				T* queryPt = new T[3];
@@ -840,7 +839,7 @@ public:
 			for(unsigned p = 0; p < R.E[e].size(); p++)
 				M += A.E[e].m(p, errormag_id);
 			M = M / A.E[e].size();
-			if(M > metric_fac)								// if the sum is larger than the metric_fac, we assume that it doesn't has corresponding edge in B
+			if(M > threshold)								
 				C[e] = (unsigned)-1;
 			else{											// if it should have corresponding edge in B, then...
 				p1 = R.E[e][R.E[e].size()/2];							
@@ -111,6 +111,14 @@ public:
 		R *= factor;
 	}
  
+	///set the radius of circle to a certain value 
+	///@param value: the new radius of the circle
+	CUDA_CALLABLE
+	void set_R(T value) 
+	{
+		R = value;
+	}
+
 	///sets the normal for the cirlce
 	///@param n: x,y,z direction of the normal.
 	CUDA_CALLABLE void
@@ -179,12 +179,12 @@ public:
 		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[0];					//allocate space for the segment length
+		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 - 1] - L[p - 2]);		//calculate the segment length using the L array
+			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
 		}
@@ -43,28 +43,44 @@ 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);					// distance between two points = length/height
-		T ax = 57.2957795*acos(dz / d);						// 180°/pi
-		if (dz < 0.0)
-			ax = -ax;
-		T rx = -dy*dz;
-		T ry = dx*dz;
+	//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);
+	//	glPushMatrix();
+	//	glTranslatef(x1, y1, z1);
+	//	glRotatef(ax, rx, ry, 0.0);
  
-		glutSolidCylinder(radius, d, subdivisions, 1);
-		glPopMatrix();
+	//	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);
@@ -72,7 +88,6 @@ public:
 		glPopMatrix();
 	}
  
-
 	/// Render the network centerline as a series of line strips.
 	/// glCenterline0 is for only one input
 	void glCenterline0(){
@@ -92,7 +107,6 @@ public:
 		glCallList(dlist);					// render the display list
 	}
  
-	/// @param m specifies the magnitude value used as the vertex weight (radius, error, etc.)
 	void glCenterline(){
  
 		if(!glIsList(dlist)){					//if dlist isn't a display list, create it
@@ -109,83 +123,245 @@ public:
 			}
 			glEndList();						//end the display list
 		}		
-		glCallList(dlist);					// render the display list
+		glCallList(dlist);						//render the display list
 	}
  
-	void glRandColorCenterlineGT(GLuint &dlist1, std::vector<unsigned> map, std::vector<T> colormap){
-		if(!glIsList(dlist1)){
+	///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)){
+			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++){
+					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);
+				else {
+					glColor3f(1.0, 1.0, 1.0);						// white color
 					glBegin(GL_LINE_STRIP);
-					for(unsigned p = 0; p < E[e].size(); p++){
+					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)){
+	/// 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)){
+			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++){
+					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);
+				else {
+					glColor3f(1.0, 1.0, 1.0);						// white color
 					glBegin(GL_LINE_STRIP);
-					for(unsigned p = 0; p < E[e].size(); p++){
+					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 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