merged changes in vec3

David Mayerich
2 parents 48f356fa 9f5c0d4a
Showing 10 changed files with 382 additions and 125 deletions Show diff stats
stim/biomodels/network.h
stim/cuda/cuda_texture.cuh
stim/cuda/filter.cuh
stim/cuda/ivote/down_sample.cuh
stim/gl/gl_spider.h
stim/grids/image_stack.h
stim/math/matrix.h
stim/math/vec3.h
stim/parser/filename.h
stim/visualization/gl_aaboundingbox.h
@@ -33,7 +33,10 @@ class network{
 		public:
 		unsigned v[2];		//unique id's designating the starting and ending
 		// default constructor
-		 edge() : cylinder<T>(){v[1] = -1; v[0] = -1;}
+		edge() : cylinder<T>()
+		{
+			v[1] = -1; v[0] = -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
@@ -54,7 +57,7 @@ class network{
 		/// Output the edge information as a string
 		std::string str(){
 			std::stringstream ss;
-			ss<<"("<<cylinder<T>::size()<<")\tl = "<<length()<<"\t"<<v[0]<<"----"<<v[1];
+			ss<<"("<<cylinder<T>::size()<<")\tl = "<<this.length()<<"\t"<<v[0]<<"----"<<v[1];
 			return ss.str();
 		}
@@ -64,9 +67,9 @@ class network{
 	class vertex : public stim::vec3<T>
 	{
 		public:
-			//std::vector<unsigned int> edges;		//indices of edges connected to this node.
-			std::vector<unsigned int> e[2];			//indices of edges going out (e[0]) and coming in (e[1])
-			//stim::vec3<T> p;						//position of this node in physical space.
+			//std::vector<unsigned int> edges;					//indices of edges connected to this node.
+			std::vector<unsigned int> e[2];						//indices of edges going out (e[0]) and coming in (e[1])
+			//stim::vec3<T> p;							//position of this node in physical space.
 			//constructor takes a stim::vec
 			vertex(stim::vec3<T> p) : stim::vec3<T>(p){}
@@ -90,6 +93,7 @@ class network{
 				return ss.str();
 			}
+
 	};
 protected:
@@ -99,6 +103,18 @@ protected:
 public:
+	///default constructor
+	network()
+	{
+		
+	}
+
+	///constructor with a file to load.
+	network(std::string fileLocation)
+	{
+		load_obj(fileLocation);
+	}
+
 	///Returns the number of edges in the network.
 	unsigned int edges(){
 		return E.size();
@@ -109,72 +125,219 @@ public:
 		return V.size();
 	}
+	std::vector<vertex> operator*(T s){
+		for (unsigned i=0; i< vertices; i ++ ){
+			V[i] = V[i] * s;
+		}
+		return V;
+	}
+
+	std::vector<vertex> operator*(vec<T> s){
+		for (unsigned i=0; i< vertices; i ++ ){
+			for (unsigned dim = 0 ; dim< 3; dim ++){
+				V[i][dim] = V[i][dim] * s[dim];
+			}
+		}
+		return V;
+	}
+
+	// Returns an average of branching index in the network
+
+	double BranchingIndex(){
+		double B=0;
+		for(unsigned v=0; v < V.size(); v ++){
+			B += ((V[v].e[0].size()) + (V[v].e[1].size()));
+		}
+		B = B / V.size();
+		return B;
+
+	}
+
+	// Returns number of branch points in thenetwork
+
+	unsigned int BranchP(){
+		unsigned int B=0;
+		unsigned int c;
+		for(unsigned v=0; v < V.size(); v ++){
+			c = ((V[v].e[0].size()) + (V[v].e[1].size()));
+			if (c > 2){
+			B += 1;}
+		}		
+		return B;
+
+	}
+
+	// Returns number of end points (tips) in thenetwork
+
+	unsigned int EndP(){
+		unsigned int B=0;
+		unsigned int c;
+		for(unsigned v=0; v < V.size(); v ++){
+			c = ((V[v].e[0].size()) + (V[v].e[1].size()));
+			if (c == 1){
+			B += 1;}
+		}		
+		return B;
+
+	}
+
+	//// Returns a dictionary with the key as the vertex
+	//std::map<std::vector<vertex>,unsigned int> DegreeDict(){
+	//	std::map<std::vector<vertex>,unsigned int> dd;
+	//	unsigned int c = 0;
+	//	for(unsigned v=0; v < V.size(); v ++){
+	//		c = ((V[v].e[0].size()) + (V[v].e[1].size()));
+	//		dd[V[v]] = c;
+	//	}
+	//	return dd;
+	//}
+
+	//// Return number of branching stems
+	//unsigned int Stems(){
+	//	unsigned int s = 0;
+	//	std::map<std::vector<vertex>,unsigned int> dd;
+	//	dd = DegreeDict();
+	//	//for(unsigned v=0; v < V.size(); v ++){
+	//	//	V[v].e[0].
+	//	return s;
+	//}
+
+
+	// Returns an average of fiber/edge lengths in the network
+	double Lengths(){
+		stim::vec<T> L;double sumLength = 0;
+		for(unsigned e = 0; e < E.size(); e++){				//for each edge in the network
+			L.push_back(E[e].length());						//append the edge length
+			sumLength = sumLength + E[e].length();
+		}
+		double avg = sumLength / E.size();
+		return avg;
+	}
+
+
+	// Returns an average of tortuosities in the network
+	double Tortuosities(){
+		stim::vec<T> t;
+		stim::vec<T> id1, id2;                        // starting and ending vertices of the edge
+		double distance;double tortuosity;double sumTortuosity = 0;
+		for(unsigned e = 0; e < E.size(); e++){				//for each edge in the network
+			id1 = E[e][0];									//get the edge starting point
+			id2 = E[e][E[e].size() - 1];					//get the edge ending point
+			distance = (id1 - id2).len();                   //displacement between the starting and ending points
+			if(distance > 0){
+				tortuosity = E[e].length()/	distance	;		// tortuoisty = edge length / edge displacement
+			}
+			else{
+				tortuosity = 0;}
+			t.push_back(tortuosity);
+			sumTortuosity += tortuosity;
+		}
+		double avg = sumTortuosity / E.size();
+		return avg;
+	}
+
+	// Returns average contraction of the network
+	double Contractions(){
+	stim::vec<T> t;
+	stim::vec<T> id1, id2;                        // starting and ending vertices of the edge
+	double distance;double contraction;double sumContraction = 0;
+	for(unsigned e = 0; e < E.size(); e++){				//for each edge in the network
+		id1 = E[e][0];									//get the edge starting point
+		id2 = E[e][E[e].size() - 1];					//get the edge ending point
+		distance = (id1 - id2).len();                   //displacement between the starting and ending points
+		contraction = distance / E[e].length();		// tortuoisty = edge length / edge displacement
+		t.push_back(contraction);
+		sumContraction += contraction;
+	}
+	double avg = sumContraction / E.size();
+	return avg;
+	}
+
+	// returns average fractal dimension of the branches of the network
+	double FractalDimensions(){
+	stim::vec<T> t;
+	stim::vec<T> id1, id2;                        // starting and ending vertices of the edge
+	double distance;double fract;double sumFractDim = 0;
+	for(unsigned e = 0; e < E.size(); e++){				//for each edge in the network
+		id1 = E[e][0];									//get the edge starting point
+		id2 = E[e][E[e].size() - 1];					//get the edge ending point
+		distance = (id1 - id2).len();                   //displacement between the starting and ending points
+		fract = std::log(distance) / std::log(E[e].length());		// tortuoisty = edge length / edge displacement
+		t.push_back(sumFractDim);
+		sumFractDim += fract;
+	}
+	double avg = sumFractDim / E.size();
+	return avg;
+	}
 	stim::cylinder<T> get_cylinder(unsigned f){
-		return E[f];					//return the specified edge (casting it to a fiber)
+		return E[f];									//return the specified edge (casting it to a fiber)
 	}
 	//load a network from an OBJ file
 	void load_obj(std::string filename){
-		stim::obj<T> O;			//create an OBJ object
-		O.load(filename);		//load the OBJ file as an object
+		stim::obj<T> O;									//create an OBJ object
+		O.load(filename);								//load the OBJ file as an object
-		std::vector<unsigned> id2vert;	//this list stores the OBJ vertex ID associated with each network vertex
+		std::vector<unsigned> id2vert;							//this list stores the OBJ vertex ID associated with each network vertex
-		unsigned i[2];					//temporary, IDs associated with the first and last points in an OBJ line
+		unsigned i[2];									//temporary, IDs associated with the first and last points in an OBJ line
 		//for each line in the OBJ object
 		for(unsigned int l = 1; l <= O.numL(); l++){
-			std::vector< stim::vec<T> > c;				//allocate an array of points for the vessel centerline
+			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());
 			for(size_t j = 0; j < c.size(); j++)
 				c3[j] = c[j];
-			edge new_edge = c3;							//create an edge from the given centerline
-			unsigned int I = new_edge.size();			//calculate the number of points on the centerline
+	//		edge new_edge = c3;		///This is dangerous.
+			edge new_edge(c3);
+					
+			//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
-			std::vector< unsigned > id;					//create an array to store the centerline point IDs
+			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
 			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
-			std::vector<unsigned>::iterator it;			//create an iterator for searching the id2vert array
+			std::vector<unsigned>::iterator it;					//create an iterator for searching the id2vert array
 			unsigned it_idx;							//create an integer for the id2vert entry index
 			//find out if the nodes for this fiber have already been created
 			it = find(id2vert.begin(), id2vert.end(), i[0]);	//look for the first node
-			it_idx = std::distance(id2vert.begin(), it);
 			if(it == id2vert.end()){							//if i[0] hasn't already been used
 				vertex new_vertex = new_edge[0];				//create a new vertex, assign it a position
-				new_vertex.e[0].push_back(E.size());			//add the current edge as outgoing
-				new_edge.v[0] = V.size();						//add the new vertex to the edge
-				V.push_back(new_vertex);						//add the new vertex to the vertex list
-				id2vert.push_back(i[0]);						//add the ID to the ID->vertex conversion list
+				new_vertex.e[0].push_back(E.size());				//add the current edge as outgoing
+				new_edge.v[0] = V.size();					//add the new edge to the edge
+				V.push_back(new_vertex);					//add the new vertex to the vertex list
+				id2vert.push_back(i[0]);					//add the ID to the ID->vertex conversion list
 			}
-			else{												//if the vertex already exists
+			else{									//if the vertex already exists
+				it_idx = std::distance(id2vert.begin(), it);
 				V[it_idx].e[0].push_back(E.size());				//add the current edge as outgoing
 				new_edge.v[0] = it_idx;
 			}
-			it = find(id2vert.begin(), id2vert.end(), i[1]);	//look for the second ID
-			it_idx = std::distance(id2vert.begin(), it);
-			if(it == id2vert.end()){							//if i[1] hasn't already been used
-				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();
-				V.push_back(new_vertex);									//add the new vertex to the vertex list
-				id2vert.push_back(i[1]);						//add the ID to the ID->vertex conversion list
+			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[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
+				id2vert.push_back(i[1]);					//add the ID to the ID->vertex conversion list
 			}
-			else{												//if the vertex already exists
+			else{									//if the vertex already exists
+				it_idx = std::distance(id2vert.begin(), it);
 				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
+			E.push_back(new_edge);							//push the edge to the list
 		}
 	}
@@ -199,17 +362,17 @@ public:
 	/// 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
+		stim::network<T> n;								//create a new network that will be an exact copy, with resampled fibers
+		n.V = V;									//copy all vertices
-		n.E.resize(edges());				//allocate space for the edge list
+		n.E.resize(edges());								//allocate space for the edge list
 		//copy all fibers, resampling them in the process
-		for(unsigned e = 0; e < edges(); e++){				//for each edge in the edge list
-			n.E[e] = E[e].resample(spacing);				//resample the edge and copy it to the new network
+		for(unsigned e = 0; e < edges(); e++){						//for each edge in the edge list
+			n.E[e] = E[e].resample(spacing);					//resample the edge and copy it to the new network
 		}
-		return n;							              //return the resampled network
+		return n;							              	//return the resampled network
 	}
@@ -236,14 +399,14 @@ public:
 	/// @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
+		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
+			L += E[e].length();							//get the edge length
 		}
-		return M / L;					//return the average magnitude
+		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.
@@ -256,17 +419,17 @@ public:
 		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 = 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
+		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 = 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);
 		unsigned t = 0;
 		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 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] = A.E[e][p][d];
@@ -276,27 +439,27 @@ public:
 		}
 		//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
+		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
-		ANNdistArray dists = new ANNdist[1];     // near neighbor distances
-		ANNidxArray nnIdx = new ANNidx[1];				// near neighbor indices // allocate near neigh indices
+		ANNdistArray dists = new ANNdist[1];     					// near neighbor distances
+		ANNidxArray nnIdx = new ANNidx[1];						// near neighbor indices // allocate near neigh indices
 		stim::vec3<T> p0, p1;
 		float m1;
-		float M = 0;											//stores the total metric value
-		float L = 0;											//stores the total network length
+		float M = 0;									//stores the total metric value
+		float L = 0;									//stores the total network length
 		ANNpoint queryPt = annAllocPt(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
+			R.E[e].add_mag(0);							//add a new magnitude for the metric
-			for(unsigned p = 0; p < R.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 = R.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
+				kdt->annkSearch( queryPt, 1, nnIdx, dists, eps);		//find the distance between A and the current network
 				m1 = 1.0f - gaussianFunction((float)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(m1, p, 1);					//set the error for the second point in the segment
 			}
 		}
@@ -308,8 +471,106 @@ public:
 	unsigned nmags(){
 		return E[0].nmags();
 	}
+	// split a string in text by the character sep
+	stim::vec<T> split(std::string &text, char sep) 
+	{
+		stim::vec<T> tokens;
+		std::size_t start = 0, end = 0;
+		while ((end = text.find(sep, start)) != std::string::npos) {
+		tokens.push_back(atof(text.substr(start, end - start).c_str()));
+		start = end + 1;
+		}
+		tokens.push_back(atof(text.substr(start).c_str()));
+		return tokens;
+	}
+	// load a network in text file to a network class
+	void load_txt(std::string filename)
+	{
+		std::vector <std::string> file_contents;
+		std::ifstream file(filename);
+		std::string line;
+		std::vector<unsigned> id2vert;	//this list stores the vertex ID associated with each network vertex
+		//for each line in the text file, store them as strings in file_contents
+		while (std::getline(file, line))
+		{
+			std::stringstream ss(line);
+			file_contents.push_back(ss.str());
+		}
+		unsigned int numEdges = atoi(file_contents[0].c_str()); //number of edges in the network
+		unsigned int I = atoi(file_contents[1].c_str()) ;		//calculate the number of points3d on the first edge
+		unsigned int count = 1; unsigned int k = 2; // count is global counter through the file contents, k is for the vertices on the edges
+		// for each edge in the network.
+		for (unsigned int i = 0; i < numEdges; i ++ )
+		{
+			// pre allocate a position vector p with number of points3d on the edge p
+			std::vector< stim::vec<T> > p(0, I);
+			// for each point on the nth edge
+		  for (unsigned int j = k; j < I + k; j++)
+		  {
+			 // split the points3d of floats with separator space and form a float3 position vector out of them
+			  p.push_back(split(file_contents[j], ' '));
+		  }
+			count += p.size() + 1; // increment count to point at the next edge in the network
+			I = atoi(file_contents[count].c_str()); // read in the points3d at the next edge and convert it to an integer
+			k = count + 1;
+			edge new_edge = p; // create an edge with a vector of points3d  on the edge
+			E.push_back(new_edge); // push the edge into the network
+		}
+		unsigned int numVertices = atoi(file_contents[count].c_str()); // this line in the text file gives the number of distinct vertices
+		count = count + 1; // this line of text file gives the first verrtex
+		// push each vertex into V
+		for (unsigned int i = 0; i < numVertices; i ++)
+		{
+			vertex new_vertex = split(file_contents[count], ' ');
+			V.push_back(new_vertex);
+			count += atoi(file_contents[count + 1].c_str()) + 2; // Skip number of edge ids + 2 to point to the next vertex
+		}
+	} // end load_txt function
-
+	// strTxt returns a string of edges
+	std::string
+	strTxt(std::vector< stim::vec<T> > p)
+	{
+		std::stringstream ss;
+		std::stringstream oss;
+		for(unsigned int i = 0; i < p.size(); i++){
+			ss.str(std::string());
+			for(unsigned int d = 0; d < 3; d++){
+				ss<<p[i][d];
+			}
+			ss < "\n";
+		}
+		return ss.str();
+	}
+	// removes specified character from string
+	void removeCharsFromString(std::string &str, char* charsToRemove ) {
+	   for ( unsigned int i = 0; i < strlen(charsToRemove); ++i ) {
+		  str.erase( remove(str.begin(), str.end(), charsToRemove[i]), str.end() );
+	   }
+	}
+	//exports network to txt file
+	void
+	to_txt(std::string filename)
+	{
+		std::ofstream ofs(filename, std::ofstream::out | std::ofstream::app);
+		int num;
+		ofs << (E.size()).str() << "\n";
+		for(unsigned int i = 0; i < E.size(); i++)
+		{
+			 std::string str;
+			 ofs << (E[i].size()).str() << "\n";
+			 str = E[i].strTxt();
+             ofs << str << "\n";
+		} 
+		for(int i = 0; i < V.size(); i++)
+		{
+			std::string str;
+			str = V[i].str();
+			removeCharsFromString(str, "[],");
+			ofs << str << "\n";
+		}
+		ofs.close();
+	}
 };		//end stim::network class
 };		//end stim namespace
 #endif
@@ -134,15 +134,15 @@ namespace stim
 		void
 		UnmapCudaTexture()
 		{
-			HANDLE_ERROR(
-				cudaGraphicsUnmapResources(1, &resource)
-			);
-			HANDLE_ERROR(
-				cudaGraphicsUnregisterResource(resource)
-			);	
-			HANDLE_ERROR(
-				cudaDestroyTextureObject(tObj)
-			);
+		//	HANDLE_ERROR(
+		//		cudaGraphicsUnmapResources(1, &resource)
+		//	);
+		//	HANDLE_ERROR(
+		//		cudaGraphicsUnregisterResource(resource)
+		//	);	
+		//	HANDLE_ERROR(
+		//		cudaDestroyTextureObject(tObj)
+		//	);
 //			HANDLE_ERROR(
 //				cudaFreeArray(srcArray)
 //			);
@@ -7,6 +7,7 @@
 #include <stdio.h>
 #include <stim/visualization/colormap.h>
 #include <sstream>
+#include <stim/math/constants.h>
 #include <stim/cuda/cudatools/devices.h>
 #include <stim/cuda/cudatools/threads.h>
 #include <stim/cuda/cuda_texture.cuh>
@@ -14,7 +15,6 @@
 #include <stim/cuda/arraymath.cuh>
 #define IMAD(a,b,c) ( __mul24((a), (b)) + (c) )
-#define M_PI 3.141592654f
 namespace stim
@@ -73,7 +73,7 @@ namespace stim
 				idx = j*kl+i;
 				x = i - kr - 0.5;
 				y = j - kr - 0.5;
-				LoG[idx] = (-1.0/M_PI/powf(sigma, 4))* (1 - (powf(x,2)+powf(y,2))/2.0/powf(sigma, 2))
+				LoG[idx] = (-1.0/PI/powf(sigma, 4))* (1 - (powf(x,2)+powf(y,2))/2.0/powf(sigma, 2))
 						*expf(-(powf(x,2)+powf(y,2))/2/powf(sigma,2));	
 				t +=LoG[idx];
 			}
@@ -98,7 +98,7 @@ namespace stim
 		int y = blockIdx.y;
 		int xi = threadIdx.x;
 		int yi = threadIdx.y;
-		float val = 0;
+	//	float val = 0;
 		float tu = (x-kr+xi)/(float)DIM_X;
 		float tv = (y-kr+yi)/(float)DIM_Y;
 		shared[xi][yi] = gpuLoG[yi*kl+xi]*(255.0-(float)tex2D<unsigned char>(texIn, tu, tv));
@@ -111,7 +111,7 @@ namespace stim
 		//y = min(y, height - 1);
 		int idx = y*DIM_X+x;
-		int k_idx;
+	//	int k_idx;
                 for(unsigned int step = blockDim.x/2; step >= 1; step >>= 1)
                 {
                         __syncthreads();
@@ -44,7 +44,7 @@ namespace stim{
 			unsigned int y_ds = (y/sigma_ds + (y %sigma_ds == 0 ? 0:1));
 			//get the number of pixels in the image
-			unsigned int pixels_ds = x_ds * y_ds;
+//			unsigned int pixels_ds = x_ds * y_ds;
 			unsigned int max_threads = stim::maxThreadsPerBlock();
 			dim3 threads(max_threads, 1);
@@ -97,4 +97,4 @@ namespace stim{
 	}
 }
-#endif
 \ No newline at end of file
+#endif
@@ -14,6 +14,7 @@
 #include <stim/math/vec3.h>
 #include <stim/math/rect.h>
 #include <stim/math/matrix.h>
+#include <stim/math/constants.h>
 #include <stim/cuda/spider_cost.cuh>
 #include <stim/cuda/cudatools/glbind.h>
 #include <stim/cuda/arraymath.cuh>
@@ -369,7 +370,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		///Stored in a display list.
 		///uses the default d vector <0,0,1>
 		void
-		genDirectionVectors(float solidAngle = 5/M_PI*4)
+		genDirectionVectors(float solidAngle = M_PI/2)
 		{
 			//Set up the vectors necessary for Rectangle creation.
@@ -584,7 +585,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 			glDeleteFramebuffers(1, &framebufferID);
 			glGenFramebuffers(1, &framebufferID);
 			glBindFramebuffer(GL_FRAMEBUFFER, framebufferID);
-			int numChannels = 1;
+//			int numChannels = 1;
 //			unsigned char* texels = new unsigned char[width * height * numChannels];
 			glGenTextures(1, &textureID);
 			glBindTexture(GL_TEXTURE_2D, textureID);
@@ -609,7 +610,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		{
 			glGenFramebuffers(1, &fboID);
 			glBindFramebuffer(GL_FRAMEBUFFER, fboID);
-			int numChannels = 1;
+//			int numChannels = 1;
 //			unsigned char* texels = new unsigned char[width * height * numChannels];
 			glGenTextures(1, &texbufferID);
 			glBindTexture(GL_TEXTURE_2D, texbufferID);
@@ -938,7 +939,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 			network_time = 0;
 			hit_time = 0;
 #endif
-			stepsize = 3.0;
+			stepsize = 2.5;
 			t_length = 16.0;
 			srand(100);	
@@ -1392,13 +1393,9 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 			Z[0] = cos(PHI[0]);
 			Z[1] = cos(PHI[1]);
-
 			range = Z[0] - Z[1];
-
-			float z, theta, phi;
-			
 			std::vector<stim::vec3<float> > vecsUni;
 			for(int i = 0; i < numSamplesPos; i++)
 			{
@@ -1490,8 +1487,6 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		void
 		trace(int min_cost)
 		{	
-			bool sEmpty = true;
-			float lastmag = 16.0;;
 			stim::vec3<float> curSeed; 
 			stim::vec3<float> curSeedVec;			
 			float curSeedMag;
@@ -1527,7 +1522,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 				gpuStartTimer();
 			#endif
-			float s = 3.0;
+//			float s = 3.0;
 			GLuint selectBuf[2048];
 			GLint hits;
 			glSelectBuffer(2048, selectBuf);
@@ -1541,17 +1536,17 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 			CHECK_OPENGL_ERROR
 		//What would that vessel see in front of it.
 				camSel.setPosition(loc);
-				camSel.setFocalDistance(mag/s);
-				camSel.LookAt((loc[0]+dir[0]*mag/s),
-					 (loc[1]+dir[1]*mag/s),
-					 (loc[2]+dir[2]*mag/s));
+				camSel.setFocalDistance(mag/stepsize);
+				camSel.LookAt((loc[0]+dir[0]*mag/stepsize),
+					 (loc[1]+dir[1]*mag/stepsize),
+					 (loc[2]+dir[2]*mag/stepsize));
 				ps  = camSel.getPosition();
 				ups = camSel.getUp();
 				ds  = camSel.getLookAt();
 				glMatrixMode(GL_PROJECTION);
 				glPushMatrix();
 				glLoadIdentity();
-				glOrtho(-mag/s/2.0, mag/s/2.0, -mag/s/2.0, mag/s/2.0, 0.0, mag/s/2.0);
+				glOrtho(-mag/stepsize/2.0, mag/stepsize/2.0, -mag/stepsize/2.0, mag/stepsize/2.0, 0.0, mag/stepsize/2.0);
 				glMatrixMode(GL_MODELVIEW);
 				glPushMatrix();
 				glLoadIdentity();
@@ -1597,21 +1592,11 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		int
 		processHits(GLint hits, GLuint buffer[])
 		{
-			GLuint names, *ptr;
-			//printf("hits = %u\n", hits);
+			GLuint *ptr;
 			ptr = (GLuint *) buffer;
-		//	for (int i = 0; i < hits; i++) { /*  for each hit  */
-			names = *ptr;
-		//		printf (" number of names for hit = %u\n", names);
 			ptr++;
 			ptr++; //Skip the minimum depth value.  
 			ptr++; //Skip the maximum depth value.
-		//		printf ("   the name is ");                                                        
-		//		for (int j = 0; j < names; j++) {     /*  for each name */
-		//			printf ("%u ", *ptr); ptr++;
-		//		}
-		//		printf ("\n");
-	//		}
 			if(hits == 0)
@@ -1675,11 +1660,18 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 			}		
 			#ifdef TIMING
-				double network_time = (std::clock() - s) / (double) CLOCKS_PER_SEC;
-				network_time += network_time * 1000.0;
+				double nt = (std::clock() - s) / (double) CLOCKS_PER_SEC;
+				network_time += nt * 1000.0;
 			#endif
 		}
+//		void
+//		addToNetwork(pair<stim::fiber<float>, int> in, stim::vec3<float> spos, 
+//				stim::vec<float> smag, stim::vec3<float> sdir)
+//		{
+//
+//		}
+
 		void
 		printSizes()
@@ -1783,9 +1775,6 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
                  		}
          		}
 		}
-
-			
-
 };
 }
 #endif
@@ -33,10 +33,10 @@ public:
 		stim::filename file_path(file_mask);
 		//if the file path is relative, update it with the current working directory
-		if(file_path.is_relative()){
-			stim::filename wd = stim::filename::cwd();
-			file_path = wd.get_relative(file_mask);
-		}
+//		if(file_path.is_relative()){
+//			stim::filename wd = stim::filename::cwd();
+//			file_path = wd.get_relative(file_mask);
+//		}
 		//get the list of files
 		std::vector<stim::filename> file_list = file_path.get_list();
@@ -134,10 +134,10 @@ public:
 		stim::filename file_path(file_mask);
 		//if the file path is relative, update it with the current working directory
-		if(file_path.is_relative()){
-			stim::filename wd = stim::filename::cwd();
-			file_path = wd.get_relative(file_mask);
-		}
+//		if(file_path.is_relative()){
+//			stim::filename wd = stim::filename::cwd();
+//			file_path = wd.get_relative(file_mask);
+//		}
 		//create a list of file names
 		std::vector<std::string> file_list = stim::wildcards::increment(file_path.str(), 0, R[3]-1, 1);
@@ -53,13 +53,13 @@ struct matrix
 	template<typename Y>
 	vec<Y> operator*(vec<Y> rhs){
-		unsigned int N = rhs.size();
+		unsigned int M = rhs.size();
 		vec<Y> result;
-		result.resize(N);
+		result.resize(M);
-		for(int r=0; r<N; r++)
-			for(int c=0; c<N; c++)
+		for(int r=0; r<M; r++)
+			for(int c=0; c<M; c++)
 				result[r] += (*this)(r, c) * rhs[c];
 		return result;
@@ -3,6 +3,7 @@
 #include <stim/cuda/cudatools/callable.h>
+#include <cmath>
 namespace stim{
@@ -203,6 +203,9 @@ public:
 		operator=(s);
 	}
+	bool is_relative(){
+		return false;
+	}
 	std::string str(){
@@ -331,7 +334,7 @@ public:
 #elif BOOST_PRECOMPILED
-		boost::filesystem::path p(dir());	//create a path from the current filename
+		boost::filesystem::path p(path());	//create a path from the current filename
 		std::vector<stim::filename> file_list;
 		if(boost::filesystem::exists(p)){
 			if(boost::filesystem::is_directory(p)){
@@ -344,9 +347,9 @@ public:
 				for (vec::const_iterator it(v.begin()), it_end(v.end()); it != it_end; ++it)
 				{
 					//if the filename is a wild card *or* it matches the read file name
-					if( prefix == "*" || prefix == (*it).filename().stem().string()){
+					if( _prefix == "*" || _prefix == (*it).filename().stem().string()){
 						//if the extension is a wild card *or* it matches the read file extension
-						if( ext == "*" || "." + ext == (*it).filename().extension().string()){
+						if( _extension == "*" || "." + _extension == (*it).filename().extension().string()){
 							file_list.push_back((*it).string());	//include it in the final list
 						}
@@ -368,4 +371,4 @@ public:
 	}
 };				//end filename
 }				//end namespace stim
-#endif
 \ No newline at end of file
+#endif
@@ -11,6 +11,9 @@ class gl_aaboundingbox : public aaboundingbox&lt;T&gt;{
 public:
+	using stim::aaboundingbox<T>::A;
+	using stim::aaboundingbox<T>::B;
+
 	//default constructor
 	gl_aaboundingbox() : stim::aaboundingbox<T>(){}
@@ -57,4 +60,4 @@ public:
 };		//end namespace stim
-#endif
 \ No newline at end of file
+#endif