added an axis-aligned bounding box class, as well as classes for rendering networks

David Mayerich
1 parent 58ee2b21
Showing 8 changed files with 348 additions and 89 deletions Show diff stats
stim/biomodels/network.h
stim/parser/arguments.h
stim/visualization/aabb.h
stim/visualization/camera.h
stim/visualization/cylinder.h
stim/visualization/gl_aabb.h
stim/visualization/gl_network.h
stim/visualization/glnetwork.h → stim/visualization/glnetwork-dep.h
@@ -2,7 +2,7 @@
 #define STIM_NETWORK_H
 #include <stdlib.h>
-#include <assert.h> 
+#include <assert.h>
 #include <sstream>
 #include <fstream>
 #include <algorithm>
@@ -50,7 +50,7 @@ class network{
 			return e;						//return the new edge
 		}
-			
+
 		/// Output the edge information as a string
 		std::string str(){
 			std::stringstream ss;
@@ -58,8 +58,8 @@ class network{
 			return ss.str();
 		}
-	};	
-	
+	};
+
 	///Node class that stores the physical position of the node as well as the edges it is connected to (edges that connect to it), As well as any additional data necessary.
 	class vertex : public stim::vec<T>
 	{
@@ -89,15 +89,15 @@ class network{
 				return ss.str();
 			}
-			
+
 	};
-	private:
+protected:
 	std::vector<edge> E;       //list of edges
 	std::vector<vertex> V;	    //list of vertices.
-	
-	public:
+
+public:
 	///Returns the number of edges in the network.
 	unsigned int edges(){
@@ -134,7 +134,7 @@ class network{
 			//get the first and last vertex IDs for the line
 			std::vector< unsigned > id;					//create an array to store the centerline point IDs
-			O.getLinei(l, id);							//get the list of point IDs for the line			
+			O.getLinei(l, id);							//get the list of point IDs for the line
 			i[0] = id.front();							//get the OBJ ID for the first element of the line
 			i[1] = id.back();							//get the OBJ ID for the last element of the line
@@ -169,7 +169,7 @@ class network{
 				V[it_idx].e[1].push_back(E.size());				//add the current edge as incoming
 				new_edge.v[1] = it_idx;
 			}
-			
+
 			E.push_back(new_edge);								//push the edge to the list
 		}
@@ -193,6 +193,7 @@ class network{
 	}
 	/// This function resamples all fibers in a network given a desired minimum spacing
+	/// @param spacing is the minimum distance between two points on the network
 	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
@@ -209,7 +210,7 @@ class network{
-	// total number of points on all fibers after resampling -- function used only on a resampled network
+	/// Calculate the total number of points on all edges.
 	unsigned total_points(){
 		unsigned n = 0;
 		for(unsigned e = 0; e < E.size(); e++)
@@ -227,28 +228,50 @@ class network{
 		a[2] = b[2];
 	}
+	/// Calculate the average magnitude across the entire network.
+	/// @param m is the magnitude value to use. The default is 0 (usually radius).
+	T average(unsigned m = 0){
+
+		T M, L;												//allocate space for the total magnitude and length
+		M = L = 0;											//initialize both the initial magnitude and length to zero
+		for(unsigned e = 0; e < E.size(); e++){				//for each edge in the network
+			M += E[e].integrate(m);							//get the integrated magnitude
+			L += E[e].length();								//get the edge length
+		}
-	/// This function compares two networks and returns a metric
-	float compare(stim::network<T> A, float sigma){
+		return M / L;					//return the average magnitude
+	}
+
+	/// This function compares two networks and returns the percentage of the current network that is missing from A.
+
+	/// @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){
+
+		stim::network<T> R;								//generate a network storing the result of the comparison
+		R = (*this);									//initialize the result with the current network
+
+		//generate a KD-tree for network A
 		float metric = 0.0;                               // initialize metric to be returned after comparing the networks
 		ANNkd_tree* kdt;                                 // initialize a pointer to a kd tree
 		double **c;						                 // centerline (array of double pointers) - points on kdtree must be double
-		unsigned int n_data = total_points();          // set the number of points
+		unsigned int n_data = A.total_points();          // set the number of points
 		c = (double**) malloc(sizeof(double*) * n_data); // allocate the array pointer
 		for(unsigned int i = 0; i < n_data; i++)		 // allocate space for each point of 3 dimensions
 			c[i] = (double*) malloc(sizeof(double) * 3);
-		//std::vector<stim::vec<T> > pointsVector = points_afterResampling(); //get vertices in the network after resampling
+
 		unsigned t = 0;
-		for(unsigned e = 0; e < E.size(); e++){					//for each edge in the network
-			for(unsigned p = 0; p < E[e].size(); p++){			//for each point in the edge
+		for(unsigned e = 0; e < A.E.size(); e++){					//for each edge in the network
+			for(unsigned p = 0; p < A.E[e].size(); p++){			//for each point in the edge
 				for(unsigned d = 0; d < 3; d++){				//for each coordinate
-					c[t][d] = E[e][p][d];
+					c[t][d] = A.E[e][p][d];
 				}
 				t++;
 			}
 		}
+		//compare each point in the current network to the field produced by A
 		ANNpointArray pts = (ANNpointArray)c;           // create an array of data points of type double
 		kdt = new ANNkd_tree(pts, n_data, 3);			// build a KD tree using the annpointarray
 		double eps = 0; // error bound
@@ -261,33 +284,24 @@ class network{
 		float l;												//stores the segment length
 		float L = 0;											//stores the total network length
 		ANNpoint queryPt = annAllocPt(3);
-		for(unsigned e = 0; e < A.E.size(); e++){					//for each edge in A
-			M = 0;												//total metric value for the fiber
-			p0 = A.E[e][0];										//get the first point in the edge
-			stim2ann(queryPt, p0);
-			kdt->annkSearch( queryPt, 1, nnIdx, dists, eps);	//find the distance between A and the current network
-			m0 = 1.0f - gaussianFunction(dists[0], sigma);		//calculate the metric value based on the distance
+		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
-			for(unsigned p = 1; p < A.E[e].size(); p++){			//for each point in the edge
+			for(unsigned p = 0; p < R.E[e].size(); p++){			//for each point in the edge
-				p1 = A.E[e][p];									//get the next point in the edge
+				p1 = R.E[e][p];									//get the next point in the edge
 				stim2ann(queryPt, p1);
 				kdt->annkSearch( queryPt, 1, nnIdx, dists, eps);	//find the distance between A and the current network
 				m1 = 1.0f - gaussianFunction(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
-				//average the metric and scale
-				l = (p1 - p0).len();							//calculate the length of the segment
-				L += l;											//add the length of this segment to the total network length
-				M += (m0 + m1)/2 * l;							//scale the metric based on segment length
-
-				//swap points
-				p0 = p1;
-				m0 = m1;
 			}
 		}
-		return M/L;
+		return R;		//return the resulting network
 	}
+
+
 };		//end stim::network class
 };		//end stim namespace
 #endif
@@ -239,35 +239,29 @@ namespace stim{
     /**The arglist class implements command line arguments.
         Example:
-    
-        1) Create an arglist instance:
-            
-            stim::arglist args;
-        2) Test to see if ANSI output is supported (arguments will use color if it is). Generally, Windows consoles don't support ANSI.
+        1) Create an arglist instance:
-            #ifdef _WIN32
-                args.set_ansi(false);
-            #endif
+            stim::arglist args;
-        3) Add arguments:
+        2) Add arguments:
             args.add("help", "prints this help");
             args.add("foo", "foo takes a single integer value", "", "[intval]");
             args.add("bar", "bar takes two floating point values", "", "[value1], [value2]");
-        4) Parse the command line:
+        3) Parse the command line:
             args.parse(argc, argv);
-        5) You generally want to immediately test for help and output available arguments:
+        4) You generally want to immediately test for help and output available arguments:
             if(args["help"].is_set())
                 std::cout<<args.str();
-        
-        6)  Retrieve values:
+
+        5)  Retrieve values:
             int foo;
             float bar1, bar2;
@@ -300,7 +294,15 @@ namespace stim{
         arglist(){
         	col_width = 0;
-        	ansi = true;
+
+			ansi = true;
+
+			//set ansi to false by default if this is a Windows system
+			//		(console doesn't support ANSI colors)
+			#ifdef _WIN32
+				ansi=false;
+			#endif
+
         }
         ///Sets ANSI support (default is on), which allows colored values in the help output.
+#ifndef STIM_AABB
+#define STIM_AABB
+
+namespace stim{
+
+
+/// This class describes a structure for an axis-aligned bounding box
+template< typename T >
+class aabb{
+
+public:
+	bool set;				//has the bounding box been set to include any points?
+	stim::vec<T> A;			//minimum point in the bounding box
+	stim::vec<T> B;			//maximum point in the bounding box
+
+	aabb(){					//constructor generates an empty bounding box
+		set = false;
+	}
+
+
+	/// Test if a point is inside of the bounding box and returns true if it is.
+
+	/// @param p is the point to be tested
+	bool test(stim::vec<T> p){
+
+		for(unsigned d = 0; d < p.size(); p++){		//for each dimension
+			if(p[d] < A[d]) return false;			//if the point is less than the minimum bound, return false
+			if(p[d] > B[d]) return false;			//if the point is greater than the max bound, return false
+		}
+		return true;
+	}
+
+	/// Expand the bounding box to include the specified point.
+
+	/// @param p is the point to be included
+	void expand(stim::vec<T> p){
+
+		if(!set){							//if the bounding box is empty, fill it with the current point
+			A = B = p;
+			set = true;
+		}
+
+		for(unsigned d = 0; d < p.size(); d++){		//for each dimension
+			if(p[d] < A[d]) A[d] = p[d];			//expand the bounding box as necessary
+			if(p[d] > B[d]) B[d] = p[d];
+		}
+	}
+
+	/// Return the center point of the bounding box as a stim::vec
+	stim::vec<T> center(){
+		return (B + A) * 0.5;
+	}
+
+	/// Return the size of the bounding box as a stim::vec
+	stim::vec<T> size(){
+		return (B - A);
+	}
+
+	/// Generate a string for the bounding box
+	std::string str(){
+		std::stringstream ss;
+		ss<<A.str()<<"----->"<<B.str();
+		return ss.str();
+	}
+
+
+};		//end stim::aabb
+
+
+};		//end namespace stim
+
+#endif
 \ No newline at end of file
@@ -188,6 +188,13 @@ public:
 		focus = 1;
 	}
+
+	/// Outputs the camera information as a string
+	std::string str(){
+		std::stringstream ss;
+		ss<<p.str()<<"----->"<<(p + d * focus).str();
+		return ss.str();
+	}
 };
 }
@@ -15,45 +15,38 @@ class cylinder
 		std::vector< stim::vec<T> > pos;	//positions of the cylinder.
 		std::vector< stim::vec<T> > mags;	//radii at each position
 		std::vector< T > L;			//length of the cylinder at each position.
-	
-		///default init	
+
+		///default init
 		void
-		init()
-		{
+		init(){
 		}
 		///inits the cylinder from a list of points (inP) and radii (inM)
 		void
-		init(std::vector<stim::vec<T> > inP, std::vector<stim::vec<T> > inM)
-		{
+		init(std::vector<stim::vec<T> > inP, std::vector<stim::vec<T> > inM){
 			pos = inP;
 			mags = inM;
 			//calculate each L.
 			L.resize(pos.size()-1);
 			T temp = (T)0;
-			for(int i = 0; i < L.size()-1; i++)
+			for(int i = 0; i < L.size(); i++)
 			{
 				temp += (pos[i] - pos[i+1]).len();
 				L[i] = temp;
 			}
 		}
-		
+
 		///returns the direction vector at point idx.
 		stim::vec<T>
-		d(int idx)
-		{
+		d(int idx){
 			return (pos[idx] - pos[idx+1]).norm();
-			
 		}
-		
-
 		///returns the total length of the line at index j.
 		T
-		getl(int j)
-		{
+		getl(int j){
 			for(int i = 0; i < j-1; ++i)
 			{
 				temp += (pos[i] - pos[i+1]).len();
@@ -64,8 +57,7 @@ class cylinder
 		///finds the index of the point closest to the length l on the lower bound.
 		///binary search.
 		int
-		findIdx(T l)
-		{
+		findIdx(T l){
 			int i = pos.size()/2;
 			while(i > 0 && i < pos.size())
 			{
@@ -85,18 +77,29 @@ class cylinder
 			return i;
 		}
-	public:
-		///default constructor
-		cylinder()
-		{
+		//initializes the length array given the current set of positions
+		void init_length(){
+			vec<T> p0, p1;
+			p0 = pos[0];						//initialize the first point in the segment to the first point in the cylinder
+			T l;								//allocate space for the segment length
+			for(unsigned p = 1; p < pos.size(); p++){		//for each point in the cylinder
+				p1 = pos[p];					//get the second point in the segment
+				l = (p1 - p0).len();			//calculate the length of the segment
+
+				if(p == 1) L[0] = l;			//set the length for the first segment
+				else L[p-1] = L[p-2] + l;		//calculate and set the running length for each additional segment
+			}
 		}
+	public:
+		///default constructor
+		cylinder(){}
+
 		///constructor to create a cylinder from a set of points, radii, and the number of sides for the cylinder.
 		///@param inP:  Vector of stim vecs composing the points of the centerline.
 		///@param inM:  Vector of stim vecs composing the radii of the centerline.
-		cylinder(std::vector<stim::vec<T> > inP, std::vector<stim::vec<T> > inM)
-		{
+		cylinder(std::vector<stim::vec<T> > inP, std::vector<stim::vec<T> > inM){
 			init(inP, inM);
 		}
@@ -104,7 +107,11 @@ class cylinder
 		///@param inP: Vector of stim vecs composing the points of the centerline
 		cylinder(std::vector< stim::vec<T> > inP){
 			std::vector< stim::vec<T> > inM;						//create an array of arbitrary magnitudes
-			inM.resize(inP.size(), 0);								//initialize the magnitude values to zero
+
+			stim::vec<T> zero;
+			zero.push_back(0);
+
+			inM.resize(inP.size(), zero);								//initialize the magnitude values to zero
 			init(inP, inM);
 		}
@@ -120,8 +127,7 @@ class cylinder
 		///interpolates the position along the line.
 		///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1).
 		stim::vec<T>
-		p(T pvalue)
-		{
+		p(T pvalue){
 			if(pvalue < 0.0 || pvalue > 1.0)
 				return;
 			T l = pvalue*L[L.size()-1];
@@ -134,8 +140,7 @@ class cylinder
 		///@param l: the location of the in the cylinder.
 		///@param idx: integer location of the point closest to l but prior to it.
 		stim::vec<T>
-		p(T l, int idx)
-		{
+		p(T l, int idx){
 			return (pos[idx] + (pos[idx+1]-pos[idx])*((l-L[idx])/(L[idx+1]- L[idx])));
 		}
@@ -143,8 +148,7 @@ class cylinder
 		///interpolates the radius along the line.
 		///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1).
 		T
-		r(T pvalue)
-		{
+		r(T pvalue){
 			if(pvalue < 0.0 || pvalue > 1.0)
 				return;
 			T l = pvalue*L[L.size()-1];
@@ -157,11 +161,33 @@ class cylinder
 		///@param l: the location of the in the cylinder.
 		///@param idx: integer location of the point closest to l but prior to it.
 		T
-		r(T l, int idx)
-		{
+		r(T l, int idx){
 			return (mags[idx] + (mags[idx+1]-mags[idx])*((l-L[idx])/(L[idx+1]- L[idx])));
 		}
+		///	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 ri(unsigned i, unsigned m = 0){
+			return mags[i][m];
+		}
+
+		/// Adds a new magnitude value to all points
+		/// @param m is the starting value for the new magnitude
+		void add_mag(T m = 0){
+			for(unsigned int p = 0; p < pos.size(); p++)
+				mags[p].push_back(m);
+		}
+
+		/// Sets a magnitude value
+		/// @param val is the new value for the magnitude
+		/// @param p is the point index for the magnitude to be set
+		/// @param m is the index for the magnitude
+		void set_mag(T val, unsigned p, unsigned m = 0){
+			mags[p][m] = val;
+		}
+
+
 		///returns the position of the point with a given pvalue and theta on the surface
 		///in x, y, z coordinates. Theta is in degrees from 0 to 360.
@@ -174,7 +200,7 @@ class cylinder
 				return;
 			T l = pvalue*L[L.size()-1];
 			int idx = findIdx(l);
-			stim::vec<T> ps = p(l, idx); 
+			stim::vec<T> ps = p(l, idx);
 			T m = r(l, idx);
 			stim::vec<T> dr = d(idx);
 			s = stim::circle<T>(ps, m, dr);
@@ -182,7 +208,7 @@ class cylinder
 		}
 		///returns a vector of points necessary to create a circle at every position in the fiber.
-		///@param sides: the number of sides of each circle.	
+		///@param sides: the number of sides of each circle.
 		std::vector<std::vector<vec<T> > >
 		getPoints(int sides)
 		{
@@ -210,9 +236,7 @@ class cylinder
 		/// Allows a point on the centerline to be accessed using bracket notation
 		vec<T> operator[](unsigned int i){
-
 			return pos[i];
-
 		}
 		/// Returns the total length of the cylinder centerline
@@ -220,13 +244,51 @@ class cylinder
 			return L.back();
 		}
+		/// Integrates a magnitude value along the cylinder.
+		/// @param m is the magnitude value to be integrated (this is usually the radius)
+		T integrate(unsigned m = 0){
+
+			T M = 0;						//initialize the integral to zero
+			T m0, m1;						//allocate space for both magnitudes in a single segment
+
+			//vec<T> p0, p1;					//allocate space for both points in a single segment
+
+			m0 = mags[0][m];				//initialize the first point and magnitude to the first point in the cylinder
+			//p0 = pos[0];
+
+			T len = L[0];						//allocate space for the segment length
+
+			//for every consecutive point in the cylinder
+			for(unsigned p = 1; p < pos.size(); p++){
+
+				//p1 = pos[p];							//get the position and magnitude for the next point
+				m1 = mags[p][m];
+
+				if(p > 1) len = (L[p-1] - L[p-2]);		//calculate the segment length using the L array
+
+				//add the average magnitude, weighted by the segment length
+				M += (m0 + m1)/2.0 * len;
+
+				m0 = m1;								//move to the next segment by shifting points
+			}
+			return M;			//return the integral
+		}
+
+		/// Averages a magnitude value across the cylinder
+		/// @param m is the magnitude value to be averaged (this is usually the radius)
+		T average(unsigned m = 0){			
+
+			//return the average magnitude
+			return integrate(m) / L.back();
+		}
+
 		/// 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){
 			std::vector< vec<T> > result;
-			
+
 			vec<T> p0 = pos[0];								//initialize p0 to the first point on the centerline
 			vec<T> p1;
 			unsigned N = size();							//number of points in the current centerline
@@ -238,7 +300,7 @@ class cylinder
 				vec<T> v = p1 - p0;							//calculate the vector between these two points
 				T d = v.len();								//calculate the distance between these two points (length of the line segment)
-				unsigned nsteps = d / spacing;		//calculate the number of steps to take along the segment to meet the spacing criteria
+				unsigned nsteps = d / spacing+1;		//calculate the number of steps to take along the segment to meet the spacing criteria
 				T stepsize = 1.0 / nsteps;			//calculate the parametric step size between new centerline points
 				//for each step along the line segment
@@ -251,11 +313,11 @@ class cylinder
 			}
 			result.push_back(pos[size() - 1]);			//push the last point in the centerline
-			
+
 			return cylinder<T>(result);
 		}
-		
+
 };
 }
+#ifndef STIM_GL_AABB
+#define STIM_GL_AABB
+
+#include "aabb.h"
+#include "GL/gl.h"
+
+namespace stim{
+
+template <typename T>
+class gl_aabb : public aabb<T>{
+
+public:
+
+	//default constructor
+	gl_aabb() : stim::aabb<T>(){}
+
+	//constructor takes an AABB
+	gl_aabb(stim::aabb<T> b) : stim::aabb<T>(b){}
+
+
+	/// Specifies vertices of the bounding box using CW winding. Use GL_LINE_LOOP for wireframe or GL_QUADS for a solid.
+	void glPointsCW(){
+
+			//front plane (in A[2])
+			glVertex3f(A[0], A[1], A[2]);
+			glVertex3f(A[0], B[1], A[2]);
+			glVertex3f(B[0], B[1], A[2]);
+			glVertex3f(B[0], A[1], A[2]);
+
+			//back plane (in B[2])
+			glVertex3f(B[0], B[1], B[2]);
+			glVertex3f(A[0], B[1], B[2]);
+			glVertex3f(A[0], A[1], B[2]);			
+			glVertex3f(B[0], A[1], B[2]);
+	}
+
+
+};		//end stim::gl_aabb
+
+
+};		//end namespace stim
+
+#endif
 \ No newline at end of file
+#ifndef STIM_GL_NETWORK
+#define STIM_GL_NETWORK
+
+#include <stim/biomodels/network.h>
+#include "aabb.h"
+
+namespace stim{
+
+template <typename T>
+class gl_network : public stim::network<T>{
+
+protected:
+	using stim::network<T>::E;
+	using stim::network<T>::V;
+
+public:
+
+	/// Default constructor
+	gl_network() : stim::network<T>(){}
+
+	/// Constructor creates a gl_network from a stim::network
+	gl_network(stim::network<T> N) : stim::network<T>(N){}
+
+	/// Fills the parameters with the minimum and maximum spatial positions in the network,
+	///     specifying a bounding box for the network geometry
+	aabb<T> boundingbox(){
+
+		aabb<T> bb;								//create a bounding box
+
+		//loop through every edge
+		for(unsigned e = 0; e < E.size(); e++){
+			//loop through every point
+			for(unsigned p = 0; p < E[e].size(); p++)
+				bb.expand(E[e][p]);						//expand the bounding box to include the point
+		}
+
+		return bb;								//return the bounding box
+	}
+
+	/// Render the network centerline as a series of line strips.
+	void glCenterline(){
+
+		for(unsigned e = 0; e < E.size(); e++){				//for each edge in the network
+			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]);
+				glTexCoord1f(E[e].ri(p));
+			}
+			glEnd();
+		}
+
+	}
+
+};		//end stim::gl_network class
+};		//end stim namespace
+
+
+
+#endif
 \ No newline at end of file