bug fixes for NetMets

David Mayerich
1 parent e21d1051
Showing 4 changed files with 104 additions and 42 deletions Show diff stats
stim/biomodels/centerline.h
stim/biomodels/network.h
stim/visualization/cylinder.h
stim/visualization/gl_network.h
@@ -190,7 +190,7 @@ public:
 		stim::vec3<T> p;      //- intermediate point to be added
 		stim::vec3<T> p1;   // p1 - starting point of an segment on the fiber,
 		stim::vec3<T> p2;   // p2 - ending point,
-		double sum=0;  //distance summation
+		//double sum=0;  //distance summation
 		size_t N = size();
@@ -33,16 +33,15 @@ class network{
 		public:
 		unsigned v[2];		//unique id's designating the starting and ending
 		// default constructor
-		edge() : cylinder<T>()
-		{
+		edge() : cylinder<T>() {
 			v[1] = (unsigned)(-1); v[0] = (unsigned)(-1);
 		}
 		/// Constructor - creates an edge from a list of points by calling the stim::fiber constructor
 		///@param p is an array of positions in space
-		edge(std::vector< stim::vec3<T> > p) : cylinder<T>(p){}
+		edge(centerline p) : cylinder<T>(p){}
-		/// Copy constructor creates an edge from a fiber
+		/// Copy constructor creates an edge from a cylinder
 		edge(stim::cylinder<T> f) : cylinder<T>(f) {}
 		/// Resamples an edge by calling the fiber resampling function
@@ -291,7 +290,7 @@ public:
 			std::vector< stim::vec<T> > c;						//allocate an array of points for the vessel centerline
 			O.getLine(l, c);							//get the fiber centerline
-			std::vector< stim::vec3<T> > c3(c.size());
+			stim::centerline<T> c3(c.size());
 			for(size_t j = 0; j < c.size(); j++)
 				c3[j] = c[j];
@@ -431,7 +430,7 @@ public:
 	/// @param A is the network to compare to - the field is generated for A
 	/// @param sigma is the user-defined tolerance value - smaller values provide a stricter comparison
-	stim::network<T> compare(stim::network<T> A, float sigma, int device){
+	stim::network<T> compare(stim::network<T> A, float sigma, int device = -1){
 		stim::network<T> R;										//generate a network storing the result of the comparison
 		R = (*this);											//initialize the result with the current network
@@ -471,6 +470,7 @@ public:
 		T* queryPt = new T[3];
 		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
 			for(unsigned p = 0; p < R.E[e].size(); p++){				//for each point in the edge
@@ -479,7 +479,7 @@ public:
 				kdt.search(queryPt, 1, nnIdx, dists);			//find the distance between A and the current network
 				m1 = 1.0f - gaussianFunction((T)dists[0], sigma);		//calculate the metric value based on the distance
-				R.E[e].set_mag(m1, p, 1);					//set the error for the second point in the segment
+				R.E[e].set_mag(errormag_id, p, m1);					//set the error for the second point in the segment
 			}
 		}
@@ -11,19 +11,10 @@ namespace stim
 template<typename T>
 class cylinder : public centerline<T> {
 protected:
-	//stim::circle<T> s;							//an arbitrary circle
-	//std::vector<stim::circle<T> > e;				//an array of circles that store the centerline
-	//std::vector<stim::vec3<T> > norms;
 	std::vector< stim::vec3<T> > U;					//stores the array of U vectors defining the Frenet frame
 	std::vector< std::vector<T> > M;				//stores a list of magnitudes for each point in the centerline (assuming mags[0] is the radius)
-	//std::vector< T > L;								//length of the cylinder at each position (pre-integration)
-
-	//using stim::centerline<T>::c;
-	//using stim::centerline<T>::N;
-	using stim::centerline<T>::size;
-	using stim::centerline<T>::at;
 	using stim::centerline<T>::findIdx;
 	//calculates the U values for each point to initialize the frenet frame
@@ -40,33 +31,18 @@ protected:
 		U[size() - 1] = c.U;							//for the last point, duplicate the final frenet frame vector
 	}
-	///adds a magnitude to each point in the cylinder
-	void add_mag(T val = 0) {
-		if (M.size() == 0) M.resize(size());	//if the magnitude vector isn't initialized, resize it to match the centerline
-		for (size_t i = 0; i < size(); i++)		//for each point
-			M[i].push_back(val);				//add this value to the magnitude vector at each point
-	}
-
-	//adds a magnitude based on a list of magnitudes for each point
-	void add_mag(std::vector<T> val) {
-		if (M.size() == 0) M.resize(size());	//if the magnitude vector isn't initialized, resize it to match the centerline
-		for (size_t i = 0; i < size(); i++)		//for each point
-			M[i].push_back(val[i]);				//add this value to the magnitude vector at each point
-	}
-
 public:
-	std::string str() {
-		std::stringstream ss;
-		size_t N = std::vector< stim::vec3<T> >::size();
-		ss << "---------[" << N << "]---------" << std::endl;
-		for (size_t i = 0; i < N; i++)
-			ss << std::vector< stim::vec3<T> >::at(i) << "   r = "<< M[i][0]<<"   u = "<<U[i]<<std::endl;
-		ss << "--------------------" << std::endl;
-		return ss.str();
+	using stim::centerline<T>::size;
+	using stim::centerline<T>::at;
+
+	cylinder() : centerline(){}
+	
+	cylinder(centerline& c) : centerline(c) {
+		init();
 	}
-	cylinder(centerline& c, T r = 0) : centerline(c) {
+	cylinder(centerline& c, T r) : centerline(c) {
 		init();
 		add_mag(r);
 	}
@@ -81,7 +57,7 @@ public:
 	///Interpolates the position along the line.
 	///@param l: the location of the in the cylinder.
 	///@param idx: integer location of the point closest to l but prior to it.
-	T m(T l, int idx, size_t i = 0) {
+	T m(T l, int idx, size_t i) {
 		T a = (l - L[idx]) / (L[idx + 1] - L[idx]);
 		T v2 = (M[idx][i] + (M[idx + 1][i] - M[idx][i])*a);
@@ -100,6 +76,37 @@ public:
 		return m(l, idx, i);
 	}
+	///	Returns the magnitude at the given index
+	///	@param i is the index of the desired point
+	/// @param m is the index of the magnitude value
+	T m(unsigned i, unsigned m) {
+		return M[i][m];
+	}
+
+	///adds a magnitude to each point in the cylinder
+	void add_mag(T val = 0) {
+		if (M.size() == 0) M.resize(size());	//if the magnitude vector isn't initialized, resize it to match the centerline
+		for (size_t i = 0; i < size(); i++)		//for each point
+			M[i].push_back(val);				//add this value to the magnitude vector at each point
+	}
+
+	//adds a magnitude based on a list of magnitudes for each point
+	void add_mag(std::vector<T> val) {
+		if (M.size() == 0) M.resize(size());	//if the magnitude vector isn't initialized, resize it to match the centerline
+		for (size_t i = 0; i < size(); i++)		//for each point
+			M[i].push_back(val[i]);				//add this value to the magnitude vector at each point
+	}
+
+	//sets the value of magnitude m at point i
+	void set_mag(size_t m, size_t i, T v) {
+		M[i][m] = v;
+	}
+
+	size_t nmags() {
+		if (M.size() == 0) return 0;
+		else return M[0].size();
+	}
+
 	///Returns a circle representing the cylinder cross section at point i
 	stim::circle<T> circ(size_t i, size_t m = 0) {
 		return stim::circle<T>(at(i), M[i][m], d(i), U[i]);
@@ -153,6 +160,61 @@ public:
 		return out;
 	}
+	std::string str() {
+		std::stringstream ss;
+		size_t N = std::vector< stim::vec3<T> >::size();
+		ss << "---------[" << N << "]---------" << std::endl;
+		for (size_t i = 0; i < N; i++)
+			ss << std::vector< stim::vec3<T> >::at(i) << "   r = " << M[i][0] << "   u = " << U[i] << std::endl;
+		ss << "--------------------" << std::endl;
+
+		return ss.str();
+	}
+
+	/// Integrates a magnitude value along the cylinder.
+	/// @param m is the magnitude value to be integrated (this is usually the radius)
+	T integrate(size_t m = 0) {
+
+		T sum = 0;						//initialize the integral to zero
+		T m0, m1;						//allocate space for both magnitudes in a single segment
+		m0 = M[0][m];					//initialize the first point and magnitude to the first point in the cylinder
+		T len = L[0];					//allocate space for the segment length
+
+											
+		for (unsigned p = 1; p < size(); p++) {				//for every consecutive point in the cylinder
+			m1 = M[p][m];
+			if (p > 1) len = (L[p - 1] - L[p - 2]);		//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
+	}
+
+	/// Resamples the cylinder to provide a maximum distance of "spacing" between centerline points. All current
+	///		centerline points are guaranteed to exist in the new cylinder
+	/// @param spacing is the maximum spacing allowed between sample points
+	cylinder<T> resample(T spacing) {
+		cylinder<T> c = stim::centerline<T>::resample(spacing);			//resample the centerline and use it to create a new cylinder
+
+		size_t nm = nmags();											//get the number of magnitude values in the current cylinder
+		if (nm > 0) {													//if there are magnitude values
+			std::vector<T> magvec(nm, 0);							//create a magnitude vector for a single point
+			c.M.resize(c.size(), magvec);									//allocate space for a magnitude vector at each point of the new cylinder
+		}
+
+		T l, t;
+		for (size_t i = 0; i < c.size(); i++) {							//for each point in the new cylinder
+			l = c.L[i];													//get the length along the new cylinder
+			t = l / length();										//calculate the parameter value along the new cylinder
+			for (size_t mag = 0; mag < nm; mag++) {							//for each magnitude value
+				c.M[i][mag] = m(t, mag);								//retrieve the interpolated magnitude from the current cylinder and store it in the new one
+			}
+		}
+
+
+		return c;
+	}
+
 		/*
 		///inits the cylinder from a list of points (std::vector of stim::vec3 --inP) and magnitudes (inM)
@@ -55,7 +55,7 @@ public:
 				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].ri(p, m));						//set the texture coordinate based on the specified magnitude index
+					glTexCoord1f(E[e].m(p, m));						//set the texture coordinate based on the specified magnitude index
 				}
 				glEnd();
 			}