merged Dev into Branch_Detection

Pavel Govyadinov
2 parents 6a13a9c2 b2d10ff1
Showing 8 changed files with 567 additions and 302 deletions Show diff stats
stim/envi/binary.h
stim/envi/envi.h
stim/gl/gl_spider.h
stim/math/matrix.h
stim/parser/arguments.h
stim/visualization/camera.h
stim/visualization/colormap.h
stim/visualization/sph_harmonics.h
@@ -225,6 +225,10 @@ public:
 		return true;
 	}
  
+	/// Reads a plane given a coordinate along the 0-axis (YZ plane)
+
+	/// @param p is a pointer to pre-allocated memory of size R[1] * R[2] * sizeof(T)
+	/// @param n is the 0-axis coordinate used to retrieve the plane
 	bool read_plane_0(T* p, unsigned int n){
  
 		if (n >= R[0]){										//make sure the number is within the possible range
@@ -247,6 +251,10 @@ public:
  
 	}
  
+	/// Reads a plane given a coordinate along the 1-axis (XZ plane)
+
+	/// @param p is a pointer to pre-allocated memory of size R[0] * R[2] * sizeof(T)
+	/// @param n is the 1-axis coordinate used to retrieve the plane
 	bool read_plane_1(T* p, unsigned int n){
  
 		unsigned int L = R[0] * sizeof(T);		//caculate the number of bytes in a sample line
@@ -268,10 +276,18 @@ public:
 		return true;
 	}
  
+	/// Reads a plane given a coordinate along the 2-axis (XY plane)
+
+	/// @param p is a pointer to pre-allocated memory of size R[0] * R[1] * sizeof(T)
+	/// @param n is the 2-axis coordinate used to retrieve the plane
 	bool read_plane_2(T* p, unsigned int n){
 		return read_page(p, n);
 	}
  
+	/// Reads a single pixel, treating the entire data set as a linear array
+
+	/// @param p is a pointer to pre-allocated memory of size sizeof(T)
+	/// @param i is the index to the pixel using linear indexing
 	bool read_pixel(T* p, unsigned int i){
 		if(i >= R[0] * R[1] * R[2]){
 			std::cout<<"ERROR read_pixel: n is out of range"<<std::endl;
@@ -283,6 +299,12 @@ public:
  
 	}
  
+	/// Reads a single pixel, given an x, y, z coordinate
+
+	/// @param p is a pointer to pre-allocated memory of size sizeof(T)
+	/// @param x is the x (0) axis coordinate
+	/// @param y is the y (1) axis coordinate
+	/// @param z is the z (2) axis coordinate
 	bool read_pixel(T* p, unsigned int x, unsigned int y, unsigned int z){
  
 		if(x < 0 || x >= R[0] || y < 0 || y >= R[1] || z < 0 || z > R[2]){
@@ -294,54 +316,6 @@ public:
 		return read_pixel(p, i);
 	}
  
-	//saves a hyperplane orthogonal to dimension d at intersection n
-	/*bool read_plane(T * dest, unsigned int d, unsigned int n){
-
-		//reset the file pointer back to the beginning of the file
-		file.seekg(0, std::ios::beg);
-
-		//compute the contiguous size C for each readable block
-		unsigned int C = 1;
-		for(unsigned int i = 0; i < d; i++)		//for each dimension less than d
-			C *= R[i];							//compute the product
-
-		//compute the non-contiguous size NC for each readable block
-		unsigned int NC = 1;
-		for(unsigned int i = d + 1; i < D; i++)
-			NC *= R[i];
-
-		//for all noncontiguous blocks, read each contiguous block that makes up the hyper-plane
-		for(unsigned int nc = 0; nc < NC; nc++){
-			file.seekg(n * C * sizeof(T), std::ios::cur);	//skip n contiguous blocks
-			file.read( (char*)&dest[nc * C], C * sizeof(T));	//read one contiguous block
-			file.seekg( (R[d] - n - 1) * C * sizeof(T), std::ios::cur);	//skip R[d] - n contiguous blocks
-		}
-
-		return true;
-
-	}*/
-
-	//save one pixel of the file into the memory, and return the pointer
-	/*bool read_spectrum(T * p, unsigned x, unsigned y){
-
-		unsigned int i;
-
-		if ( x >= R[0] || y >= R[1]){							//make sure the sample and line number is right
-			std::cout<<"ERROR: sample or line out of range"<<std::endl;
-			return false;
-		}
-
-		file.seekg((x + y * R[0]) * sizeof(T), std::ios::beg);           //point to the certain sample and line
-		for (i = 0; i < R[2]; i++)
-		{
-			file.read((char *)(p + i), sizeof(T));
-			file.seekg((R[1] * R[0] - 1) * sizeof(T), std::ios::cur);    //go to the next band
-		}
-
-		return true;	
-	}*/
-
-
 };
  
 }
@@ -710,6 +710,11 @@ public:
 		return false;
 	}
  
+	/// Retrieve a spectrum from the specified location
+
+	/// @param ptr is a pointer to pre-allocated memory of size B*sizeof(T)
+	/// @param x is the x-coordinate of the spectrum
+	/// @param y is the y-coordinate of the spectrum
 	bool spectrum(void* ptr, unsigned int x, unsigned int y){
  
 		if(header.interleave == envi_header::BSQ){		//if the infile is bsq file
@@ -29,32 +29,33 @@ namespace stim
 {
  
 template<typename T>
-class gl_spider : public virtual gl_texture<T>
+class gl_spider
 {
 	//doen't use gl_texture really, just needs the GLuint id.
 	//doesn't even need the texture iD really.
 	private:
-		stim::vec<float> position;  	//vector designating the position of the spider.
-		stim::vec<float> direction;	//vector designating the orientation of the spider
+		stim::vec<float> p;  	//vector designating the position of the spider.
+		stim::vec<float> d;	//vector designating the orientation of the spider
 						//always a unit vector.
-		stim::vec<float> magnitude;	//magnitude of the direction vector.
+		stim::vec<float> m;	//magnitude of the spider vector.
 						//mag[0] = length.
 						//mag[1] = width.
-		std::vector<stim::vec<float> > dirVectors;
-		std::vector<stim::vec<float> > posVectors;
-		std::vector<stim::vec<float> > magVectors;
-		stim::matrix<float, 4> currentTransform;
-		using gl_texture<T>::texID;
-		using gl_texture<T>::S;
-		using gl_texture<T>::R;
+		std::vector<stim::vec<float> > dV;
+		std::vector<stim::vec<float> > pV;
+		std::vector<stim::vec<float> > mV;
+		//currentTransform
+		stim::matrix<float, 4> cT;
+		GLuint texID;
+		stim::vec<float> S;
+		stim::vec<float> R;
 		cudaGraphicsResource_t resource;
  
 		GLuint dList;
-		GLubyte list[4];
 		GLuint fboID;
 		GLuint texbufferID;
 		int iter;		//temporary for testing
 		int numSamples;
+		float stepsize = 3.0;
  
 		/// Method for finding the best scale for the spider.
 		/// changes the x, y, z size of the spider to minimize the cost
@@ -62,45 +63,41 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		void
 		findOptimalDirection()
 		{
-			//genTemplate(dirVectors, 0);
-			Bind();
-			positionTemplate();
+			//genTemplate(dV, 0);
+			setMatrix();
 			glCallList(dList);
 			int best = getCost();
 			stim::vec<float, 4> next;
- 			next[0] = dirVectors[best][0]*S[1]*R[1];
-			next[1] = dirVectors[best][1]*S[2]*R[2];
-			next[2] = dirVectors[best][2]*S[3]*R[3];
+ 			next[0] = dV[best][0]*S[0]*R[0];
+			next[1] = dV[best][1]*S[1]*R[1];
+			next[2] = dV[best][2]*S[2]*R[2];
 			next[3] = 1;
-			next = (currentTransform*next).norm();
-			setPosition(	position[0]+next[0]*magnitude[0]/3,
-					position[1]+next[1]*magnitude[0]/3,
-					position[2]+next[2]*magnitude[0]/3);
+			next = (cT*next).norm();
+			setPosition(	p[0]+next[0]*m[0]/stepsize,
+					p[1]+next[1]*m[0]/stepsize,
+					p[2]+next[2]*m[0]/stepsize);
 			setDirection(next[0], next[1], next[2]);
-			Unbind();
 		}
  
-		/// Method for finding the best direction for the spider.
-		/// Not sure if necessary since the next position for the spider
-		/// will be at direction * magnitude.
+		/// Method for finding the best d for the spider.
+		/// Not sure if necessary since the next p for the spider
+		/// will be at d * m.
 		void
 		findOptimalPosition()
 		{
-			Bind();
-			positionTemplate();
+			setMatrix();
 			glCallList(dList+1);
 			int best = getCost();
 			stim::vec<float, 4> next;
- 			next[0] = posVectors[best][0];
-			next[1] = posVectors[best][1];
-			next[2] = posVectors[best][2];
+ 			next[0] = pV[best][0];
+			next[1] = pV[best][1];
+			next[2] = pV[best][2];
 			next[3] = 1;
-			next = currentTransform*next;
-			std::cout << "Optimal Position:"<< next << std::endl;
-			setPosition(	next[0]*S[1]*R[1],
-					next[1]*S[2]*R[2],
-					next[2]*S[3]*R[3]);
-			Unbind();
+			next = cT*next;
+			std::cout << "Optimal p:"<< next << std::endl;
+			setPosition(	next[0]*S[0]*R[0],
+					next[1]*S[1]*R[1],
+					next[2]*S[2]*R[2]);
 		}
  
 		/// Method for finding the best scale for the spider.
@@ -109,19 +106,17 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		void
 		findOptimalScale()
 		{
-			Bind();
-			positionTemplate();
+			setMatrix();
 			glCallList(dList+2);
 			int best = getCost();
 			stim::vec<float, 4> next;
- 			next[0] = magVectors[best][0]*S[1]*R[1];
-			next[1] = magVectors[best][1]*S[2]*R[2];
-			next[2] = magVectors[best][2]*S[3]*R[3];
+ 			next[0] = mV[best][0]*S[0]*R[0];
+			next[1] = mV[best][1]*S[1]*R[1];
+			next[2] = mV[best][2]*S[2]*R[2];
 			next[3] = 1;
-			next = currentTransform*next;
+			next = cT*next;
 			std::cout << "Optimal Scale:"<< next << std::endl;
 			setMagnitude(next[0]);
-			Unbind();
 		}
  
 		void
@@ -140,7 +135,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		void
 		Optimize()
 		{
-			/*find the optimum direction and scale */ 
+			/*find the optimum d and scale */ 
 		}
  
  
@@ -152,35 +147,53 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
  
 		///@param solidAngle, the size of the arc to sample.
 		///Method for populating the vector arrays with sampled vectors.
-		///uses the default direction vector <0,0,1>
+		///uses the default d vector <0,0,1>
 		void
 		genDirectionVectors(float solidAngle = M_PI)
 		{
-			float x = 0.0;
+			//Set up the vectors necessary for Rectangle creation.
 			vec<float> Y(1.0,0.0,0.0);
 			vec<float> pos(0.0,0.0,0.0);
 			vec<float> mag(1.0, 1.0, 1.0);
 			vec<float> dir(0.0, 0.0, 1.0);
  
 			vec<float> d_s = direction.cart2sph().norm();
+			//Set up the variable necessary for vector creation.
+			vec<float> d_s = d.cart2sph().norm();
 			vec<float> temp;
 			int dim = (sqrt(numSamples)-1)/2;
 			float p0  	= -M_PI;
 			float dt  	= solidAngle/(2.0 * ((float)dim + 1.0));
 			float dp  	= p0/(2.0*((float)dim + 1.0));
+			
+			glNewList(dList, GL_COMPILE);
+			//Loop over the space
+			int idx = 0;
 			for(int i = -dim; i <= dim; i++){
 				for(int j = -dim; j <= dim; j++){
  
 					//Create linear index
-					
+					idx = (j+dim)+(i+dim)*((dim*2)+1);	
 					temp[0] = d_s[0]; 			//rotate vector
 					temp[1] = d_s[1]+dp*(float) i;
 					temp[2] = d_s[2]+dt*(float) j;
  
 					temp = (temp.sph2cart()).norm();	//back to cart
-					dirVectors.push_back(temp);
+					dV.push_back(temp);
+				 	if(cos(Y.dot(temp))< 0.087){                                                                             Y[0] = 0.0; Y[1] = 1.0;}
+					else{Y[0] = 1.0; Y[1] = 0.0;}
+
+                                       	hor = stim::rect<float>(mag,
+						 pos, temp,
+                       				((Y.cross(temp)).cross(temp)).norm());
+       				 	ver = stim::rect<float>(mag,
+						 pos, temp,
+                       				hor.n());
+					UpdateBuffer(0.0, 0.0+idx*10.0);
+					CHECK_OPENGL_ERROR
 				}
 			}
+			glEndList();
 		}
  
 		///@param solidAngle, the size of the arc to sample.
@@ -189,38 +202,61 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		void
 		genPositionVectors(float delta = 0.2)
 		{
-			ofstream file;
-			//file.open("dvectors.txt");	
+			//Set up the vectors necessary for Rectangle creation.
+			vec<float> Y(1.0,0.0,0.0);
+			vec<float> pos(0.0,0.0,0.0);
+			vec<float> mag(1.0, 1.0, 1.0);
+			vec<float> dir(0.0, 0.0, 1.0);
+
+			//Set up the variable necessary for vector creation.
 			vec<float> temp;
 			int dim = (sqrt(numSamples)-1)/2;
 			stim::rect<float> samplingPlane =
-				 stim::rect<float>(magnitude[0]*delta, position, direction);
+				 stim::rect<float>(m[0]*delta, p, d);
 			float step = 1.0/(dim);
-			//Loop over the samples, keeping the original position sample
-			//in the center of the resulting texture.
  
+			//Loop over the samples, keeping the original p sample
+			//in the center of the resulting texture.
+			int idx;
+			glNewList(dList+1, GL_COMPILE);
 			for(int i = -dim; i <= dim; i++){
 				for(int j = -dim; j <= dim; j++){
 					//Create linear index
+					idx = (j+dim)+(i+dim)*((dim*2)+1);	
  
 					temp = samplingPlane.p(
 							0.5+step*i,
 								 0.5+step*j
 										);
-					//file << temp[0] << ";" << temp[1]
-					//	 << ";" << temp[2] << endl;
-					posVectors.push_back(temp);
+					pV.push_back(temp);
+                			hor = stim::rect<float>(mag,
+						 temp, dir,
+                      				((Y.cross(d)).cross(d))
+						.norm());
+                			ver = stim::rect<float>(mag,
+						 temp, dir,
+                        			hor.n());
+					UpdateBuffer(0.0, 0.0+idx*10.0);
+				CHECK_OPENGL_ERROR
 				}
 			}
+			glEndList();
 		}
  
 		///@param solidAngle, the size of the arc to sample.
 		///Method for populating the buffer with the sampled texture.
-		///uses the default magnitude <1,1,0>
+		///uses the default m <1,1,0>
 		void
 		genMagnitudeVectors(float delta = 0.5)
 		{
  
+			//Set up the vectors necessary for Rectangle creation.
+			vec<float> Y(1.0,0.0,0.0);
+			vec<float> pos(0.0,0.0,0.0);
+			vec<float> mag(1.0, 1.0, 1.0);
+			vec<float> dir(0.0, 0.0, 1.0);
+
+			//Set up the variable necessary for vector creation.
 			int dim = (sqrt(numSamples)-1)/2;
 			float min 	= 1.0-delta;
 			float max 	= 1.0+delta;
@@ -228,95 +264,23 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 			float factor;
 			vec<float> temp;
  
+			glNewList(dList+2, GL_COMPILE);
 			for(int i = 0; i < numSamples; i++){
 				//Create linear index
-				factor = (min+step*i)*magnitude[0];
+				factor = (min+step*i)*m[0];
 				temp = factor;
-				magVectors.push_back(temp);	
+				mV.push_back(temp);	
+				hor = stim::rect<float>(temp,
+					 pos, dir, 
+       	       				((Y.cross(d)).cross(d))
+					.norm());
+               			ver = stim::rect<float>(temp,
+					 pos, dir,
+                       			hor.n());
+				UpdateBuffer(0.0, 0.0+i*10.0);
+			CHECK_OPENGL_ERROR
 			}
-		}
-		///@param vector of stim::vec in that stores all of the samplable vectors.
-		///@param type, one of three operations, 0 for Direction vectors
-		///	1 for Position, 2 for Magnitude.
-		///Function for filling the buffer up with the data based on the vectors
-		///Each vector represents a two rectangular templates.
-		///Loops through all of the vectors and transfers rect. associated with it
-		///Into buffer-space.
-		void
-		genTemplate(std::vector<stim::vec<float> > in, int type)
-		{
-			float x = 0.0;
-			vec<float> Y(1.0,0.0,0.0);
-			vec<float> pos(0.0,0.0,0.0);
-			vec<float> mag(1.0, 1.0, 1.0);
-			vec<float> dir(0.0, 0.0, 1.0);
-                        switch(type) {
-				case 0:          //Direction
-					Bind();
-					glNewList(dList+type, GL_COMPILE);
-					for(int i = 0; i < in.size(); i++)
-					{
-					 	if(cos(Y.dot(in[i]))< 0.087){                                                                             Y[0] = 0.0; Y[1] = 1.0;}
-						else{Y[0] = 1.0; Y[1] = 0.0;}
-
-                                         	hor = stim::rect<float>(mag,
-							 pos, in[i],
-                        				((Y.cross(in[i])).cross(in[i])).norm());
-              				 	ver = stim::rect<float>(mag,
-							 pos, in[i],
-                        				hor.n());
-						UpdateBuffer(x, x+i*10.0);
-					 }
-					 glEndList();
-					 Unbind();
-               				 break;
-        			case 1:		 //Position
-					Bind();	
-					glNewList(dList+type, GL_COMPILE);
-					if(cos(Y.dot(direction))< 0.087){                                                                             Y[0] = 0.0; Y[1] = 1.0;}
-					else{Y[0] = 1.0; Y[1] = 0.0;}
-
-					for(int i = 0; i < in.size(); i++)
-					{
-                				hor = stim::rect<float>(mag,
-							 in[i], dir,
-                      					((Y.cross(direction)).cross(direction))
-							.norm());
-                				ver = stim::rect<float>(mag,
-							 in[i], dir,
-                        				hor.n());
-						UpdateBuffer(x, x+i*10.0);
-					 }
-					 glEndList();
-					 Unbind();
-                			 break;
-        			case 2:		  //Scale
-					Bind();	
-					glNewList(dList+type, GL_COMPILE);
-					if(cos(Y.dot(direction))< 0.087){                                                                             Y[0] = 0.0; Y[1] = 1.0;}
-					else{Y[0] = 1.0; Y[1] = 0.0;}
-
-                			for(int i = 0; i < in.size(); i++)
-					{
-						hor = stim::rect<float>(in[i],
-							 pos, dir,
-                      					((Y.cross(direction)).cross(direction))
-							.norm());
-                				ver = stim::rect<float>(in[i],
-							 pos, dir,
-                        				hor.n());
-						UpdateBuffer(x, x+i*10.0);
-					}
-					glEndList();
-					Unbind();
-					break;
-        			default:
-                			 std::cout << "unknown case have been passed"
-                         			<< std::endl;
-					 break;
-			}               
-		Unbind();
-		CHECK_OPENGL_ERROR
+			glEndList();
 		}
 		///@param v_x x-coordinate in buffer-space,
 		///@param v_y y-coordinate in buffer-space.
@@ -392,7 +356,6 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 					);
 				glVertex2f(v_x+len, v_y+len);
 			glEnd(); 
-			CHECK_OPENGL_ERROR
 		}
  
  
@@ -468,25 +431,25 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
  
 		///Method for using the gl manipulation to alighn templates from
 		///Template space (-0.5 0.5) to Texture space (0.0, 1.0),
-		///Based on the position of the spider in real space (arbitrary).
-		void positionTemplate()
+		///Based on the p of the spider in real space (arbitrary).
+		void setMatrix()
 		{
-			stim::vec<float, 4> rot = getRotation(direction);
+			stim::vec<float, 4> rot = getRotation(d);
 			glMatrixMode(GL_TEXTURE);
 			glLoadIdentity();
  
 			glRotatef(rot[0], rot[1], rot[2], rot[3]);
-			glScalef(1.0/S[1]/R[1], 1.0/S[2]/R[2], 1.0/S[3]/R[3]);
-			glTranslatef(position[0],
-				     position[1],
-				     position[2]);
-			glScalef(magnitude[0],
-				 magnitude[1],
-				 magnitude[0]);
+			glScalef(1.0/S[0]/R[0], 1.0/S[1]/R[1], 1.0/S[2]/R[2]);
+			glTranslatef(p[0],
+				     p[1],
+				     p[2]);
+			glScalef(m[0],
+				 m[1],
+				 m[0]);
  
 			float curTrans[16];
 			glGetFloatv(GL_TEXTURE_MATRIX, curTrans);
-			fillTransform(curTrans);
+			cT.set(curTrans);
 			printTransform();
  
 			CHECK_OPENGL_ERROR
@@ -545,21 +508,11 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 //--------------------------------------------------------------------------//
  
  
-		///Default Constructor
-		gl_spider
-		()
-		{
-			setPosition(0.0,0.0,0.0);
-			setDirection(0.0,0.0,1.0);
-			setMagnitude(1.0);
-			numSamples = 1089;
-			//numSamples = 9;
-		}
 		///@param samples, the number of samples this spider is going to use.
 		///best results if samples is can create a perfect root.
 		///Default Constructor
 		gl_spider
-		(int samples)
+		(int samples = 1089)
 		{
 			setPosition(0.0,0.0,0.0);
 			setDirection(0.0,0.0,1.0);
@@ -577,9 +530,18 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 			setMagnitude(mag_x);
  
 		}
+	
+		~gl_spider
+		(void)
+		{
+			Unbind();
+			glDeleteTextures(1, &texbufferID);
+			glDeleteBuffers(1, &fboID);
+		}
+
 		///@param GLuint id texture that is going to be sampled.
 		///Attached the spider to the texture with the given GLuint ID.
-		///Samples in the default direction acting as the init method.
+		///Samples in the default d acting as the init method.
 		///Also acts an init.	
 		void
 		attachSpider(GLuint id)
@@ -587,110 +549,116 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 			texID = id;
 			iter = 0; ///for debugging purposes
 			GenerateFBO(20, numSamples*10);
-			genDirectionVectors();
-			genPositionVectors();
-			genMagnitudeVectors();
-
-			gl_texture<T>::setDims(0.6, 0.6, 1.0);
-			setSize(512, 512, 426);
+			setDims(0.6, 0.6, 1.0);
+			setSize(512.0, 512.0, 426.0);
  
-			dList = glGenLists(4);
-			list[0] = 0; list[1] = 1; list[2] = 2; list[3] = 3;
+			dList = glGenLists(3);
 			glListBase(dList);
-			genTemplate(dirVectors, 0);
-			genTemplate(posVectors, 1);
-			genTemplate(magVectors, 2);
+			Bind();
+			genDirectionVectors(5*M_PI/4);
+			genPositionVectors();
+			genMagnitudeVectors();
 			DrawCylinder();
+			Unbind();
 		}
  
 //--------------------------------------------------------------------------//
 //-----------------------------ACCESS METHODS-------------------------------//
 //--------------------------------------------------------------------------//
-		///Returns the position vector.
+		///Returns the p vector.
 		vec<float>
 		getPosition()
 		{
-			return position;
+			return p;
 		}
  
-		///Returns the direction vector.
+		///Returns the d vector.
 		vec<float>
 		getDirection()
 		{
-			return direction;
+			return d;
 		}
  
-		///Returns the magnitude vector.
+		///Returns the m vector.
 		vec<float>
 		getMagnitude()
 		{
-			return magnitude;
+			return m;
 		}
  
-		///@param vector pos, the new position.
-		///Sets the position vector to input vector pos.
+		///@param vector pos, the new p.
+		///Sets the p vector to input vector pos.
 		void
 		setPosition(vec<float> pos)
 		{
-			position = pos;
+			p = pos;
 		}
  
 		///@param x x-coordinate.
 		///@param y y-coordinate.
 		///@param z z-coordinate.
-		///Sets the position vector to the input float coordinates x,y,z.
+		///Sets the p vector to the input float coordinates x,y,z.
 		void
 		setPosition(float x, float y, float z)
 		{
-			position[0] = x;
-			position[1] = y;
-			position[2] = z;
+			p[0] = x;
+			p[1] = y;
+			p[2] = z;
 		}
  
-		///@param vector dir, the new direction.
-		///Sets the direction vector to input vector dir.
+		///@param vector dir, the new d.
+		///Sets the d vector to input vector dir.
 		void
 		setDirection(vec<float> dir)
 		{
-			direction = dir;
+			d = dir;
 		}
  
 		///@param x x-coordinate.
 		///@param y y-coordinate.
 		///@param z z-coordinate.
-		///Sets the direction vector to the input float coordinates x,y,z.
+		///Sets the d vector to the input float coordinates x,y,z.
 		void
 		setDirection(float x, float y, float z)
 		{
-			direction[0] = x;
-			direction[1] = y;
-			direction[2] = z;
+			d[0] = x;
+			d[1] = y;
+			d[2] = z;
 		}
  
-		///@param vector dir, the new direction.
-		///Sets the magnitude vector to the input vector mag.	
+		///@param vector dir, the new d.
+		///Sets the m vector to the input vector mag.	
 		void
 		setMagnitude(vec<float> mag)
 		{
-			magnitude[0] = mag[0];
-			magnitude[1] = mag[0];
+			m[0] = mag[0];
+			m[1] = mag[0];
 		}
  
 		///@param mag size of the sampled region.
-		///Sets the magnitude vector to the input mag for both templates.
+		///Sets the m vector to the input mag for both templates.
 		void
 		setMagnitude(float mag)
 		{
-			magnitude[0] = mag;
-			magnitude[1] = mag;
+			m[0] = mag;
+			m[1] = mag;
 		}
  
+
+		void
+		setDims(float x, float y, float z)
+		{
+			S[0] = x;
+			S[1] = y;
+			S[2] = z;
+		}
+
 		void
-		setSize(int x, int y, int z)
+		setSize(float x, float y, float z)
 		{
-			R[1] = x;
-			R[2] = y;
-			R[3] = z;
+			R[0] = x;
+			R[1] = y;
+			R[2] = z;
 		}
  
 		///@param dir, the vector to which we are rotating
@@ -716,18 +684,6 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 			}
 			return out;
 		}
-		///Function that fills the transform data from GL output to a matrix format.
-		///@param mat, a 16 value matrix representing the transformation in row order.
-		void
-		fillTransform(float mat[16])
-		{
-			for(int r = 0; r < 4; r++){
-				for(int c = 0; c < 4; c++){
-					currentTransform(r,c) = mat[c*4+r];
-				}
-			}
-
-		}
  
 		///Function to get back the framebuffer Object attached to the spider.
 		///For external access.
@@ -746,11 +702,11 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		Update()
 		{
 			vec<float> Y(1.0,0.0,0.0);
-			if(cos(Y.dot(direction))< 0.087){
+			if(cos(Y.dot(d))< 0.087){
 				Y[0] = 0.0; Y[1] = 1.0;}
-			hor = stim::rect<float>(magnitude, position, direction.norm(),
-				((Y.cross(direction)).cross(direction)).norm());
-			ver = stim::rect<float>(magnitude, position, direction.norm(),
+			hor = stim::rect<float>(m, p, d.norm(),
+				((Y.cross(d)).cross(d)).norm());
+			ver = stim::rect<float>(m, p, d.norm(),
 				 hor.n());
 		}
  
@@ -758,17 +714,19 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		void
 		Step()
 		{
-	//		findOptimalDirection();
-	//		findOptimalPosition();
-	//		findOptimalScale();
+			Bind();
+			findOptimalDirection();
+			findOptimalPosition();
+			findOptimalScale();
 			branchDetection();
+			Unbind();
 		}
  
  
 		void
 		printTransform()
 		{
-			std::cout << currentTransform << std::endl;
+			std::cout << cT << std::endl;
 		}
  
 		/* Method for initializing the cuda devices, necessary only
@@ -25,6 +25,17 @@ struct matrix
 					(*this)(r, c) = 0;
 	}
  
+	CUDA_CALLABLE matrix(T rhs[N*N])
+	{
+		memcpy(M,rhs, sizeof(T)*N*N);
+	}
+
+	CUDA_CALLABLE matrix<T,N> set(T rhs[N*N])
+	{
+		memcpy(M, rhs, sizeof(T)*N*N);
+		return *this;
+	}
+
 	CUDA_CALLABLE T& operator()(int row, int col)
 	{
 		return M[col * N + row];
@@ -39,6 +50,7 @@ struct matrix
 		return *this;
 	}
  
+
 	template<typename Y>
 	CUDA_CALLABLE vec<Y, N> operator*(vec<Y, N> rhs)
 	{
@@ -51,18 +63,18 @@ struct matrix
 		return result;
 	}
  
-	CUDA_CALLABLE std::string toStr()
+	std::string toStr()
 	{
 		std::stringstream ss;
  
 		for(int r = 0; r < N; r++)
 		{
-			ss<<"| ";
+			ss << "| ";
 			for(int c=0; c<N; c++)
 			{
-				ss<<(*this)(r, c)<<" ";
+				ss << (*this)(r, c) << " ";
 			}
-			ss<<"|"<<std::endl;
+			ss << "|" << std::endl;
 		}
  
 		return ss.str();
@@ -256,14 +256,16 @@ namespace stim{
             args.add("foo", "foo takes a single integer value", "", "[intval]");
             args.add("bar", "bar takes two floating point values", "", "[value1], [value2]");
  
-        4) You generally want to immediately test for help and output available arguments:
+        4) Parse the command line:
+
+            args.parse(argc, argv);
+
+        5) You generally want to immediately test for help and output available arguments:
  
             if(args["help"].is_set())
                 std::cout<<args.str();
  
-        5) Parse the command line:
-
-            args.parse(argc, argv);
+        
  
         6)  Retrieve values:
  
@@ -436,7 +438,8 @@ namespace stim{
             }
  
             //set the last option
-            set(name, params);
+            if(name != "")
+            	set(name, params);
         }
  
         ///Determines of a parameter has been set and returns true if it has
@@ -167,14 +167,14 @@ public:
 	//output the camera settings
 	const void print(std::ostream& output)
 	{
-		output<<"Position: "<<p<<std::endl;
+		output<<"Position: "<<p.str()<<std::endl;
  
 	}
 	friend std::ostream& operator<<(std::ostream& out, const camera& c)
 	{
-		out<<"Position: "<<c.p<<std::endl;
-		out<<"Direction: "<<c.d<<std::endl;
-		out<<"Up: "<<c.up<<std::endl;
+		out<<"Position: "<<c.p.str()<<std::endl;
+		out<<"Direction: "<<c.d.str()<<std::endl;
+		out<<"Up: "<<c.up.str()<<std::endl;
 		out<<"Focal Distance: "<<c.focus<<std::endl;
 		return out;
 	}
@@ -221,7 +221,11 @@ static void cpuApplyBrewer(T* cpuSource, unsigned char* cpuDest, unsigned int N,
     for(int i=0; i<N; i++)
     {
         //compute the normalized value on [minVal maxVal]
-        float a = (cpuSource[i] - minVal) / (maxVal - minVal);
+		float a;
+		if(minVal != maxVal)
+			a = (cpuSource[i] - minVal) / (maxVal - minVal);
+		else
+			a = 0.5;
         if(a < 0) a = 0;
         if(a > 1) a = 1;
  
@@ -263,11 +267,15 @@ static void cpu2cpu(T* cpuSource, unsigned char* cpuDest, unsigned int nVals, T 
         int i;
         float a;
         float range = valMax - valMin;
+
         for(i = 0; i<nVals; i++)
         {
             //normalize to the range [valMin valMax]
-            a = (cpuSource[i] - valMin) / range;
-
+			if(range != 0)
+				a = (cpuSource[i] - valMin) / range;
+			else
+				a = 0.5;
+	
             if(a < 0) a = 0;
             if(a > 1) a = 1;
  
@@ -279,22 +287,26 @@ static void cpu2cpu(T* cpuSource, unsigned char* cpuDest, unsigned int nVals, T 
 }
  
 template<class T>
-static void cpu2cpu(T* cpuSource, unsigned char* cpuDest, unsigned int nVals, colormapType cm = cmGrayscale, bool positive = false)
+static void cpu2cpu(T* cpuSource, unsigned char* cpuDest, unsigned int nVals, colormapType cm = cmGrayscale)//, bool positive = false)
 {
     //computes the max and min range automatically
  
     //find the largest magnitude value
-    T maxVal = fabs(cpuSource[0]);
-    for(int i=0; i<nVals; i++)
+    T maxVal = cpuSource[0];
+    T minVal = cpuSource[0];
+    for(int i=1; i<nVals; i++)
 	{
-        if(fabs(cpuSource[i]) > maxVal)
-            maxVal = fabs(cpuSource[i]);
+        if(cpuSource[i] > maxVal)
+            maxVal = cpuSource[i];
+        if(cpuSource[i] < minVal)
+            minVal = cpuSource[i];
 	}
  
-    if(positive)
-        cpu2cpu(cpuSource, cpuDest, nVals, (T)0, maxVal, cm);
-    else
-        cpu2cpu(cpuSource, cpuDest, nVals, -maxVal, maxVal, cm);
+    //if(positive)
+    //    cpu2cpu(cpuSource, cpuDest, nVals, (T)0, maxVal, cm);
+    //else
+    //    cpu2cpu(cpuSource, cpuDest, nVals, -maxVal, maxVal, cm);
+    cpu2cpu(cpuSource, cpuDest, nVals, minVal, maxVal, cm);
  
 }
  
+#ifndef STIM_SPH_HARMONICS
+#define STIM_SPH_HARMONICS
+
+#include <GL/gl.h>
+
+#include <stim/gl/error.h>
+#include <stim/visualization/colormap.h>
+#include <vector>
+
+#define PI 3.14159
+#define WIRE_SCALE 1.001
+namespace stim{
+
+	class sph_harmonics{
+
+	private:
+
+		double* func;		//stores the raw function data (samples at each point)
+
+		GLuint color_tex;	//texture map that acts as a colormap for the spherical function
+		
+		unsigned int N;		//resolution of the spherical grid
+
+		std::vector<double> C;	//list of SH coefficients
+
+
+		//evaluates an associated Legendre polynomial (-l <= m <= l)
+		double P(int l,int m,double x)
+		{
+			// evaluate an Associated Legendre Polynomial P(l,m,x) at x
+			double pmm = 1.0;
+			if(m>0) {
+				double somx2 = sqrt((1.0-x)*(1.0+x));
+				double fact = 1.0;
+				for(int i=1; i<=m; i++) {
+					pmm *= (-fact) * somx2;
+					fact += 2.0;
+				}
+			}
+			if(l==m) return pmm;
+			double pmmp1 = x * (2.0*m+1.0) * pmm;
+			if(l==m+1) return pmmp1;
+			double pll = 0.0;
+			for(int ll=m+2; ll<=l; ++ll) {
+				pll = ( (2.0*ll-1.0)*x*pmmp1-(ll+m-1.0)*pmm ) / (ll-m);
+				pmm = pmmp1;
+				pmmp1 = pll;
+			}
+			return pll;
+		}
+
+		//recursively calculate a factorial given a positive integer n
+		unsigned int factorial(unsigned int n) {
+			if (n == 0)
+			   return 1;
+			return n * factorial(n - 1);
+		}
+
+		//calculate the SH scaling constant
+		double K(int l, int m){
+
+			// renormalisation constant for SH function
+			double temp = ((2.0*l+1.0)*factorial(l-m)) / (4.0*PI*factorial(l+m));
+			return sqrt(temp);
+		}
+
+		//calculate the value of the SH basis function (l, m) at (theta, phi)
+			//here, theta = [0, PI], phi = [0, 2*PI]
+		double SH(int l, int m, double theta, double phi){
+			// return a point sample of a Spherical Harmonic basis function
+			// l is the band, range [0..N]
+			// m in the range [-l..l]
+			// theta in the range [0..Pi]
+			// phi in the range [0..2*Pi]
+			const double sqrt2 = sqrt(2.0);
+			if(m==0) return K(l,0)*P(l,m,cos(theta));
+			else if(m>0) return sqrt2*K(l,m)*cos(m*phi)*P(l,m,cos(theta));
+			else return sqrt2*K(l,-m)*sin(-m*phi)*P(l,-m,cos(theta));
+		}
+
+		void gen_function(){
+
+			//initialize the function to zero
+			memset(func, 0, sizeof(double) * N * N);
+
+			double theta, phi;
+			double result;
+			int l, m;
+
+			l = m = 0;
+			for(unsigned int c = 0; c < C.size(); c++){
+				
+
+				for(unsigned int xi = 0; xi < N; xi++)
+					for(unsigned int yi = 0; yi < N; yi++){
+
+						theta = (2 * PI) * ((double)xi / (N-1));
+						phi = PI * ((double)yi / (N-1));
+						result = C[c] * SH(l, m, phi, theta);		//phi and theta are reversed here (damn physicists)
+						func[yi * N + xi] += result;
+					}
+
+				m++;			//increment m
+
+				//if we're in a new tier, increment l and set m = -l
+				if(m > l){		
+					l++;
+					m = -l;
+				}
+			}
+		}
+
+		void gl_prep_draw(){
+
+			//enable depth testing
+				//this has to be used instead of culling because the sphere can have negative values
+			glEnable(GL_DEPTH_TEST);
+			glDepthMask(GL_TRUE);
+			glEnable(GL_TEXTURE_2D);	//enable 2D texture mapping
+		}
+
+		//draw a texture mapped sphere representing the function surface
+		void gl_draw_sphere() {
+
+			//PI is used to convert from spherical to cartesian coordinates
+			//const double PI = 3.14159;
+
+			//bind the 2D texture representing the color map
+			glBindTexture(GL_TEXTURE_2D, color_tex);
+
+			//Draw the Sphere
+			int i, j;
+
+			for(i = 1; i <= N-1; i++) {
+				double phi0 = PI * ((double) (i - 1) / (N-1));
+				double phi1 = PI * ((double) i / (N-1));
+
+				glBegin(GL_QUAD_STRIP);
+				for(j = 0; j <= N; j++) {
+
+					//calculate the indices into the function array
+					int phi0_i = i-1;
+					int phi1_i = i;
+					int theta_i = j;
+					if(theta_i == N)
+						theta_i = 0;
+
+					double v0 = func[phi0_i * N + theta_i];
+					double v1 = func[phi1_i * N + theta_i];
+
+					v0 = fabs(v0);
+					v1 = fabs(v1);
+
+
+					double theta = 2 * PI * (double) (j - 1) / N;
+					double x0 = v0 * cos(theta) * sin(phi0);
+					double y0 = v0 * sin(theta) * sin(phi0);
+					double z0 = v0 * cos(phi0);
+
+					double x1 = v1 * cos(theta) * sin(phi1);
+					double y1 = v1 * sin(theta) * sin(phi1);
+					double z1 = v1 * cos(phi1);
+
+					glTexCoord2f(theta / (2 * PI), phi0 / PI);
+					glVertex3f(x0, y0, z0);
+
+					glTexCoord2f(theta / (2 * PI), phi1 / PI);
+					glVertex3f(x1, y1, z1);
+				}
+				glEnd();
+			}
+		}
+
+		//draw a wire frame sphere representing the function surface
+		void gl_draw_wireframe() {
+
+					//PI is used to convert from spherical to cartesian coordinates
+					//const double PI = 3.14159;
+
+					//bind the 2D texture representing the color map
+					glDisable(GL_TEXTURE_2D);
+					glColor3f(0.0f, 0.0f, 0.0f);
+
+					//Draw the Sphere
+					int i, j;
+
+					for(i = 1; i <= N-1; i++) {
+						double phi0 = PI * ((double) (i - 1) / (N-1));
+						double phi1 = PI * ((double) i / (N-1));
+
+						glBegin(GL_LINE_STRIP);
+						for(j = 0; j <= N; j++) {
+
+							//calculate the indices into the function array
+							int phi0_i = i-1;
+							int phi1_i = i;
+							int theta_i = j;
+							if(theta_i == N)
+								theta_i = 0;
+
+							double v0 = func[phi0_i * N + theta_i];
+							double v1 = func[phi1_i * N + theta_i];
+
+							v0 = fabs(v0);
+							v1 = fabs(v1);
+
+
+							double theta = 2 * PI * (double) (j - 1) / N;
+							double x0 = WIRE_SCALE * v0 * cos(theta) * sin(phi0);
+							double y0 = WIRE_SCALE * v0 * sin(theta) * sin(phi0);
+							double z0 = WIRE_SCALE * v0 * cos(phi0);
+
+							double x1 = WIRE_SCALE * v1 * cos(theta) * sin(phi1);
+							double y1 = WIRE_SCALE * v1 * sin(theta) * sin(phi1);
+							double z1 = WIRE_SCALE * v1 * cos(phi1);
+
+							glTexCoord2f(theta / (2 * PI), phi0 / PI);
+							glVertex3f(x0, y0, z0);
+
+							glTexCoord2f(theta / (2 * PI), phi1 / PI);
+							glVertex3f(x1, y1, z1);
+						}
+						glEnd();
+					}
+				}
+
+		void init(unsigned int n){
+
+			//set the sphere resolution
+			N = n;
+
+			//allocate space for the color map
+			unsigned int bytes = N * N * sizeof(unsigned char) * 3;
+			unsigned char* color_image;
+			color_image = (unsigned char*) malloc(bytes);
+
+			//allocate space for the function
+			func = (double*) malloc(N * N * sizeof(double));
+
+			//generate a function (temporary)
+			gen_function();
+
+			//generate a colormap from the function
+			stim::cpu2cpu<double>(func, color_image, N*N, stim::cmBrewer);
+
+			//prep everything for drawing
+			gl_prep_draw();			
+
+			//generate an OpenGL texture map in the current context
+			glGenTextures(1, &color_tex);
+			//bind the texture
+			glBindTexture(GL_TEXTURE_2D, color_tex);
+
+			//copy the color data from the buffer to the GPU
+			glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, N, N, 0, GL_RGB, GL_UNSIGNED_BYTE, color_image);
+
+			//initialize all of the texture parameters
+			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);
+			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);
+			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
+			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
+			glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE);
+			
+			//free the buffer
+			free(color_image);
+		}
+
+
+	public:
+
+		void glRender(){
+			//set all OpenGL parameters required for drawing
+			gl_prep_draw();
+
+			//draw the sphere
+			gl_draw_sphere();
+			//gl_draw_wireframe();
+
+		}
+
+		void glInit(unsigned int n){
+			init(n);
+		}
+
+		void push(double c){
+			C.push_back(c);
+		}
+
+
+
+
+
+	};		//end class sph_harmonics
+
+
+
+
+}
+
+
+#endif