Adapted classes to work with the new stim::cylinder class.

David Mayerich
1 parent 32ba91f2
Showing 5 changed files with 66 additions and 28 deletions Show diff stats
stim/biomodels/centerline.h
stim/biomodels/network.h
stim/visualization/camera.h
stim/visualization/cylinder.h
stim/visualization/gl_network.h
@@ -57,10 +57,14 @@ protected:
 	///finds the index of the point closest to the length l on the lower bound.
 	///binary search.
 	size_t findIdx(T l) {
-		size_t i = L.size() / 2;
+		for (size_t i = 1; i < L.size(); i++) {				//for each point in the centerline
+			if (L[i] > l) return i - 1;						//if we have passed the desired length value, return i-1
+		}
+		return L.size() - 1;
+		/*size_t i = L.size() / 2;
 		size_t max = L.size() - 1;
 		size_t min = 0;
-		while (i > 0 && i < L.size() - 1){
+		while (i < L.size() - 1){
 			if (l < L[i]) {
 				max = i;
 				i = min + (max - min) / 2;
@@ -73,7 +77,7 @@ protected:
 				i = min + (max - min) / 2;
 			}
 		}
-		return i;
+		return i;*/
 	}
 	///Returns a position vector at the given length into the fiber (based on the pvalue).
@@ -118,6 +122,10 @@ public:
 		return p(l, idx);
 	}
+	///Update centerline internal parameters (currently the L vector)
+	void update() {
+		init();
+	}
 	///Return the length of the entire centerline
 	T length() {
 		return L.back();
@@ -338,6 +338,7 @@ public:
 			stim::centerline<T> c3(c.size());
 			for(size_t j = 0; j < c.size(); j++)
 				c3[j] = c[j];
+			c3.update();
 	//		edge new_edge = c3;		///This is dangerous.
 			edge new_edge(c3);
@@ -512,8 +513,8 @@ public:
 		kdt.create(c, n_data, MaxTreeLevels);				// build a KD tree
 		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;		//get the id for the new magnitude
+			//R.E[e].add_mag(0);							//add a new magnitude for the metric
+			//size_t errormag_id = R.E[e].nmags() - 1;		//get the id for the new magnitude
 			size_t n = R.E[e].size();						// the number of points in current edge
 			T* queryPt = new T[3 * n];
@@ -540,7 +541,7 @@ public:
 			cudaMemcpy(m1, d_m1, n * sizeof(T), cudaMemcpyDeviceToHost);
 			for(unsigned p = 0; p < n; p++){
-				R.E[e].set_mag(errormag_id, p, m1[p]);
+				R.E[e].set_r(p, m1[p]);
 			}
 			//d_set_mag<<<blocks, threads>>>(R.E[e].M, errormag_id, n, m1);
@@ -608,8 +609,8 @@ public:
 		relation.resize(A.E.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;
+			//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
@@ -643,7 +644,7 @@ public:
 			cudaMemcpy(m1, d_m1, n * sizeof(T), cudaMemcpyDeviceToHost);
 			for(unsigned p = 0; p < n; p++){
-				A.E[e].set_mag(errormag_id, p, m1[p]);											// set the error(radius) value to every point in current edge
+				A.E[e].set_r(p, m1[p]);											// set the error(radius) value to every point in current edge
 			}
 			relation[e].resize(n);																// resize every edge relation size
@@ -723,10 +724,10 @@ public:
 				if(p == E[e].size() - 2)																// if there is no splitting index, set the id to the last point index of current edge
 					id = p + 1;
 			}
-			unsigned errormag_id = E[e].nmags() - 1;
+			//unsigned errormag_id = E[e].nmags() - 1;
 			T G = 0;																					// test to see whether it has its nearest neighbor
 			for(unsigned i = 0; i < E[e].size(); i++)
-				G += E[e].m(i, errormag_id);															// won't split special edges
+				G += E[e].r(i);																			// won't split special edges
 			if(G / E[e].size() > metric_fac)															// should based on the color map
 				id = E[e].size() - 1;																	// set split idx to outgoing direction vertex
@@ -795,12 +796,12 @@ public:
 		kdt.create(c, n_data, MaxTreeLevels);				// build a KD tree
 		for(unsigned e = 0; e < n; e++){					//for each edge in A
-			size_t errormag_id = A.E[e].nmags() - 1;		//get the id for the new magnitude
+			//size_t errormag_id = A.E[e].nmags() - 1;		//get the id for the new magnitude
 			//pre-judge to get rid of impossibly mapping edges
 			T M = 0;
 			for(unsigned p = 0; p < A.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 > metric_fac)
 				C[e] = (unsigned)-1;						//set the nearest edge of impossibly mapping edges to maximum of unsigned
@@ -862,9 +863,9 @@ public:
 	}
 	/// Returns the number of magnitude values stored in each edge. This should be uniform across the network.
-	unsigned nmags(){
-		return E[0].nmags();
-	}
+	//unsigned nmags(){
+	//	return E[0].nmags();
+	//}
 	// split a string in text by the character sep
 	stim::vec<T> split(std::string &text, char sep) 
 	{
 #include <stim/math/vector.h>
 #include <stim/math/quaternion.h>
-#include <stim/math/matrix.h>
+#include <stim/math/matrix_sq.h>
 #include <ostream>
@@ -92,7 +92,7 @@ public:
 		quaternion<float> final = qz*qy*qx;
 		//get the rotation matrix
-		matrix<float, 3> rot_matrix = final.toMatrix3();
+		matrix_sq<float, 3> rot_matrix = final.toMatrix3();
 		//apply the rotation
 		d = rot_matrix*d;
@@ -21,6 +21,7 @@ protected:
 	//	this function assumes that the centerline has already been set
 	void init() {
 		U.resize(size());								//allocate space for the frenet frame vectors
+		R.resize(size());
 		stim::circle<T> c;								//create a circle
 		stim::vec3<T> d0, d1;
@@ -44,13 +45,13 @@ public:
 	cylinder(centerline& c, T r) : centerline(c) {
 		init();
-		add_mag(r);
+		//add_mag(r);
 	}
 	//initialize a cylinder with a list of points and magnitude values
 	cylinder(centerline& c, std::vector<T> r) : centerline(c){
 		init();
-		add_mag(r);
+		//add_mag(r);
 	}
 	///Returns magnitude i at the given length into the fiber (based on the pvalue).
@@ -73,15 +74,15 @@ public:
 		T l = pvalue*L[L.size() - 1];
 		int idx = findIdx(l);
-		return r(l, idx, i);
+		return r(l, idx);
 	}
 	///	Returns the magnitude at the given index
 	///	@param i is the index of the desired point
 	/// @param r is the index of the magnitude value
-	T r(unsigned i) {
-		return R[i];
-	}
+	//T r(unsigned i) {
+	//	return R[i];
+	//}
 	///adds a magnitude to each point in the cylinder
 	/*void add_mag(V val = 0) {
@@ -210,9 +211,37 @@ public:
 			c.R[i] = r(t);								//retrieve the interpolated magnitude from the current cylinder and store it in the new one
 			//}
 		}
+		return c;
+	}
+	std::vector< stim::cylinder<T> > split(unsigned int idx) {
+
+		unsigned N = size();
+		std::vector< stim::centerline<T> > LL;
+		LL.resize(2);
+		LL = (*this).centerline<T>::split(idx);
+		std::vector< stim::cylinder<T> > C(LL.size());
+		unsigned i = 0;
+		C[0] = LL[0];
+		//C[0].R.resize(idx);
+		for (; i < idx; i++) {
+			//for(unsigned d = 0; d < 3; d++)
+			//C[0][i][d] = LL[0].c[i][d];
+			C[0].R[i] = R[i];
+			//C[0].R[i].resize(1);
+		}
+		if (C.size() == 2) {
+			C[1] = LL[1];
+			//C[1].M.resize(N - idx);
+			for (; i < N; i++) {
+				//for(unsigned d = 0; d < 3; d++)
+				//C[1][i - idx][d] = LL[1].c[i - idx][d];
+				C[1].R[i - idx] = R[i];
+				//C[1].M[i - idx].resize(1);
+			}
+		}
-		return c;
+		return C;
 	}
@@ -63,17 +63,17 @@ public:
 	}
 	/// @param m specifies the magnitude value used as the vertex weight (radius, error, etc.)
-	void glCenterline(unsigned m = 0){
+	void glCenterline(){
 		if(!glIsList(dlist)){					//if dlist isn't a display list, create it
 			dlist = glGenLists(3);				//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;
+				//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].m(p, errormag_id));						//set the texture coordinate based on the specified magnitude index
+					glTexCoord1f(E[e].r(p));						//set the texture coordinate based on the specified magnitude index
 				}
 				glEnd();
 			}