Previous push did not work correctly. Added function 'sift-mask' which saves onl…

…y those spectral values corresponding to mask != 0. Syntax is 'hsiproc InputFile OutputFile --sift-mask MaskName'

Previous push did not work correctly. Added function 'sift-mask' which saves onl…
…y those spectral values corresponding to mask != 0. Syntax is 'hsiproc InputFile OutputFile --sift-mask MaskName'
Brad Deutsch
1 parent 08ec34fa
Showing 6 changed files with 383 additions and 215 deletions Show diff stats
stim/envi/bil.h
stim/envi/bip.h
stim/envi/bsq.h
stim/envi/envi.h
stim/gl/gl_spider.h
stim/math/rect.h
@@ -228,6 +228,7 @@ public:
 		return true;
 	}
+
 	/// Perform baseline correction given a list of baseline points and stores the result in a new BSQ file.
@@ -823,9 +824,9 @@ public:
 		for (unsigned i = 0; i < R[1]; i++)			//for each value in R[1] (BIP should be X)
 		{
 			getY(temp, i);							//retrieve an ZX slice, stored in temp
-			for ( unsigned j = 0; j < R[2]; j++)	//for each R[2] (BIP should be Y)
+			for ( unsigned j = 0; j < R[2]; j++)	//for each R[2] (Y)
 			{
-				for (unsigned k = 0; k < R[0]; k++)
+				for (unsigned k = 0; k < R[0]; k++) //for each band
 				{
 				if(p[i * R[0] + k] == 0)
 					temp[j * R[0] + k] = 0;
@@ -840,6 +841,44 @@ public:
 		return true;
 	}
+	///Saves to disk only those spectra corresponding to mask values != 0
+	bool sift_mask(std::string outfile, unsigned char* p){
+		// Assume R[0] = X, R[1] = Y, R[2] = Z.
+		std::ofstream target(outfile.c_str(), std::ios::binary);
+
+		//for loading pages:
+		unsigned XZ = R[0] * R[2];		//calculate the number of values in an XZ page on disk
+		unsigned L = XZ * sizeof(T);	//calculate the size of the page (in bytes)
+		T * temp = (T*)malloc(L);		//allocate memory for a temporary page
+
+		//for saving spectra:
+		unsigned Z = R[2];						//calculate the number of values in a spectrum
+		unsigned LZ = Z * sizeof(T);			//calculate the size of the spectrum (in bytes)
+		T * tempZ = (T*)malloc(LZ);				//allocate memory for a temporary spectrum
+		spectrum(tempZ, R[0] - 1, R[1] - 1);	//creates a dummy spectrum by taking the last spectrum in the image
+
+		for (unsigned i = 0; i < R[1]; i++)			//Select a page by choosing Y coordinate, R[1]
+		{
+			getY(temp, i);							//retrieve an ZX page, store in "temp"
+			for (unsigned j = 0; j < R[0]; j++)		//Select a pixel by choosing X coordinate in the page, R[0]
+			{
+				if (p[j * R[0] + i] != 0)					//if the mask != 0 at that XY pixel
+				{
+					for (unsigned k = 0; k < R[2]; k++)		//Select a voxel by choosing Z coordinate at the pixel
+					{
+						tempZ[k] = temp[k*R[0] + i];		//Pass the correct spectral value from XZ page into the spectrum to be saved. 
+					}
+					target.write(reinterpret_cast<const char*>(tempZ), LZ);		//write that spectrum to disk. Size is L2.
+				}
+				else
+					continue;
+			}
+		}
+		target.close();
+		free(temp);
+		return true;
+	}
+
 	/// Calculate the mean band value (average along B) at each pixel location.
 	/// @param p is a pointer to memory of size X * Y * sizeof(T) that will store the band averages.
@@ -917,7 +917,7 @@ public:
 			{
 				for (unsigned k = 0; k < R[2]; k++)	//for each B value (band)
 				{
-				if(p[i * R[0] + k] == 0)			//if the mask value is zero
+					if (p[i * R[0] + j] == 0)	//if the mask value is zero
 					temp[j * R[2] + k] = 0;			//set the pixel value to zero
 				else								//otherwise just continue
 					continue;
@@ -930,7 +930,46 @@ public:
 		return true;						//return true
 	}
-	/// Calculate the mean band value (average along B) at each pixel location.
+
+	///Saves to disk only those spectra corresponding to mask values != 0
+	bool sift_mask(std::string outfile, unsigned char* p){
+		// Assume R[0] = X, R[1] = Y, R[2] = Z.
+		std::ofstream target(outfile.c_str(), std::ios::binary);
+
+		//for loading pages:
+		unsigned ZX = R[0] * R[2];		//calculate the number of values in an XZ page on disk
+		unsigned L = ZX * sizeof(T);	//calculate the size of the page (in bytes)
+		T * temp = (T*)malloc(L);		//allocate memory for a temporary page
+
+		//for saving spectra:
+		unsigned Z = R[2];						//calculate the number of values in a spectrum
+		unsigned LZ = Z * sizeof(T);			//calculate the size of the spectrum (in bytes)
+		T * tempZ = (T*)malloc(LZ);				//allocate memory for a temporary spectrum
+		spectrum(tempZ, R[0] - 1, R[1] - 1);	//creates a dummy spectrum by taking the last spectrum in the image
+
+		for (unsigned i = 0; i < R[1]; i++)			//Select a page by choosing Y coordinate, R[1]
+		{
+			getY(temp, i);							//retrieve an ZX page, store in "temp"
+			for (unsigned j = 0; j < R[0]; j++)		//Select a pixel by choosing X coordinate in the page, R[0]
+			{
+				if (p[j * R[0] + i] != 0)					//if the mask != 0 at that XY pixel
+				{
+					for (unsigned k = 0; k < R[2]; k++)		//Select a voxel by choosing Z coordinate at the pixel
+					{
+						tempZ[k] = temp[i*R[2] + k];		//Pass the correct spectral value from XZ page into the spectrum to be saved. 
+					}
+					target.write(reinterpret_cast<const char*>(tempZ), LZ);		//write that spectrum to disk. Size is L2.
+				}
+				else
+					continue;
+			}
+		}
+		target.close();
+		free(temp);
+		return true;
+	}
+
+
 	/// @param p is a pointer to memory of size X * Y * sizeof(T) that will store the band averages.
 	bool band_avg(T* p){
@@ -751,25 +751,55 @@ public:
 		T * temp = (T*)malloc(L);
-		for (unsigned i = 0; i < R[2]; i++)
+		for (unsigned i = 0; i < R[2]; i++)			//for each spectral bin
 		{
-			band_index(temp, i);
-			for ( unsigned j = 0; j < XY; j++)
+			band_index(temp, i);					//get the specified band (by index)
+			for ( unsigned j = 0; j < XY; j++)		// for each pixel
 			{
-				if(p[j] == 0){
-					temp[j] = 0;
+				if(p[j] == 0){						//if the mask is 0 at that pixel
+					temp[j] = 0;					//set temp to zero
 				}
 				else{
 					continue;
 				}
 			}
-			target.write(reinterpret_cast<const char*>(temp), L);   //write a band data into target file	
+			target.write(reinterpret_cast<const char*>(temp), L);   //write the XY slice at that band to disk
 		}
 		target.close();
 		free(temp);
 		return true;
 	}
+	///Saves to disk only those spectra corresponding to mask values != 0
+	bool sift_mask(std::string outfile, unsigned char* p){
+		std::ofstream target(outfile.c_str(), std::ios::binary);
+		// open a band (XY plane)
+		unsigned XY = R[0] * R[1]; //Number of XY pixels
+		unsigned L = XY * sizeof(T); //size of XY pixels
+
+		T * temp = (T*)malloc(L); //allocate memory for a band
+		T * temp_vox = (T*)malloc(sizeof(T)); //allocate memory for one voxel
+
+		for (unsigned i = 0; i < R[2]; i++)			//for each spectral bin
+		{
+			band_index(temp, i);					//get the specified band (XY sheet by index)
+			for (unsigned j = 0; j < XY; j++)		// for each pixel
+			{
+				if (p[j] != 0){						//if the mask is != 0 at that pixel
+					temp_vox[0] = temp[j];
+					target.write(reinterpret_cast<const char*>(temp_vox), sizeof(T));   //write the XY slice at that band to disk
+				}
+				else{
+					continue;
+				}
+			}
+		}
+		target.close();
+		free(temp);
+		return true;
+	}
+
+
 	/// Calculate the mean band value (average along B) at each pixel location.
 	/// @param p is a pointer to memory of size X * Y * sizeof(T) that will store the band averages.
@@ -335,38 +335,76 @@ public:
 	/// @param outfile is the name of the resulting masked output file
 	/// @param p is memory of size X*Y containing the mask (0 = false, all other values are true)
-	bool apply_mask(std::string outfile, unsigned char* p)	
+	bool apply_mask(std::string outfile, unsigned char* p)
 	{
-		if(header.interleave == envi_header::BSQ){		//if the infile is bsq file
-			if(header.data_type ==envi_header::float32)
+		if (header.interleave == envi_header::BSQ){		//if the infile is bsq file
+			if (header.data_type == envi_header::float32)
 				return ((bsq<float>*)file)->apply_mask(outfile, p);
-			else if(header.data_type == envi_header::float64)
+			else if (header.data_type == envi_header::float64)
 				return ((bsq<double>*)file)->apply_mask(outfile, p);
 			else
-				std::cout<<"ERROR: unidentified data type"<<std::endl;
+				std::cout << "ERROR: unidentified data type" << std::endl;
 		}
-		else if(header.interleave == envi_header::BIL){		//if the infile is bil file
-			if(header.data_type ==envi_header::float32)
+		else if (header.interleave == envi_header::BIL){		//if the infile is bil file
+			if (header.data_type == envi_header::float32)
 				return ((bil<float>*)file)->apply_mask(outfile, p);
-			else if(header.data_type == envi_header::float64)
+			else if (header.data_type == envi_header::float64)
 				return ((bil<double>*)file)->apply_mask(outfile, p);
 			else
-				std::cout<<"ERROR: unidentified data type"<<std::endl;
+				std::cout << "ERROR: unidentified data type" << std::endl;
 		}
-		else if(header.interleave == envi_header::BIP){		//if the infile is bip file
-			if(header.data_type ==envi_header::float32)
+		else if (header.interleave == envi_header::BIP){		//if the infile is bip file
+			if (header.data_type == envi_header::float32)
 				return ((bip<float>*)file)->apply_mask(outfile, p);
-			else if(header.data_type == envi_header::float64)
+			else if (header.data_type == envi_header::float64)
 				return ((bip<double>*)file)->apply_mask(outfile, p);
 			else
-				std::cout<<"ERROR: unidentified data type"<<std::endl;
+				std::cout << "ERROR: unidentified data type" << std::endl;
 		}
 		else{
-			std::cout<<"ERROR: unidentified file type"<<std::endl;
+			std::cout << "ERROR: unidentified file type" << std::endl;
+			exit(1);
+		}
+		return false;
+	}
+
+	/// sift-mask saves in an array only those spectra corresponding to nonzero values of the mask. 
+	bool sift_mask(std::string outfile, unsigned char* p)
+	{
+
+		if (header.interleave == envi_header::BSQ){		//if the infile is bsq file
+			if (header.data_type == envi_header::float32)
+				return ((bsq<float>*)file)->sift_mask(outfile, p);
+			else if (header.data_type == envi_header::float64)
+				return ((bsq<double>*)file)->sift_mask(outfile, p);
+			else
+				std::cout << "ERROR: unidentified data type" << std::endl;
+		}
+
+		else if (header.interleave == envi_header::BIL){		//if the infile is bil file
+			if (header.data_type == envi_header::float32)
+				return ((bil<float>*)file)->sift_mask(outfile, p);
+			else if (header.data_type == envi_header::float64)
+				return ((bil<double>*)file)->sift_mask(outfile, p);
+			else
+				std::cout << "ERROR: unidentified data type" << std::endl;
+		}
+
+		else if (header.interleave == envi_header::BIP){		//if the infile is bip file
+			if (header.data_type == envi_header::float32)
+				return ((bip<float>*)file)->sift_mask(outfile, p);
+			else if (header.data_type == envi_header::float64)
+				return ((bip<double>*)file)->sift_mask(outfile, p);
+			else
+				std::cout << "ERROR: unidentified data type" << std::endl;
+		}
+
+		else{
+			std::cout << "ERROR: unidentified file type" << std::endl;
 			exit(1);
 		}
 		return false;
@@ -41,7 +41,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		cudaGraphicsResource_t resource;
 		GLuint fboID;
 		GLuint texbufferID;
-		int iter = 0;			//temporary for testing
+		int iter = 0;
 		void
 		findOptimalPosition()
@@ -185,45 +185,39 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		}
 		void
-		UpdatePlanes(float v_x, float v_y, vec<float> vctr, int type = 0)
+		Update(float v_x, float v_y, vec<float> dir)
 		{
 			vec<float> Y(1.0,0.0,0.0);
-			if(cos(Y.dot(vctr))< 0.087){                                                                             Y[0] = 0.0; Y[1] = 1.0;}
-                        switch(type) {
-				case 0:
-                                         hor = stim::rect<float>(magnitude, position, vctr.norm(),
-                        			((Y.cross(vctr)).cross(vctr)).norm());
-              				 ver = stim::rect<float>(magnitude, position, vctr.norm(),
-                        			hor.n());
-               				 break;
-        			case 1:
-                			 hor = stim::rect<float>(magnitude, vctr, direction,
-                      				((Y.cross(direction)).cross(direction)).norm());
-                			 ver = stim::rect<float>(magnitude, vctr, direction,
-                        			hor.n());
-                			 break;
-        			case 2:
-                			 hor = stim::rect<float>(vctr, position, direction,
-                      				((Y.cross(direction)).cross(direction)).norm());
-                			 ver = stim::rect<float>(vctr, position, direction,
-                        			hor.n());
-					 break;
-        			default:
-                			 std::cout << "unknown case have been passed"
-                         			<< std::endl;
-					 break;
-			}               
+			if(cos(Y.dot(dir))< 0.087){
+				Y[0] = 0.0; Y[1] = 1.0;}
+			hor = stim::rect<float>(magnitude, position, dir.norm(),
+				((Y.cross(dir)).cross(dir)).norm());
+			ver = stim::rect<float>(magnitude, position, dir.norm(),
+				 hor.n());
 			UpdateBuffer(v_x, v_y);
 		}
-		
-		/*
-		Method for controling the buffer and texture binding in order to properly
-		do the render to texture.
-		*/
+
+
 		void
-		Bind(int numSamples = 1089)
+		Sample(vec<float> in = (0,0,1), int numSamples = 1089, int solidAngle = M_PI/2)
 		{
-			float len = 10.0;
+			GenerateFBO(20, numSamples*10);
+
+			ofstream file;
+			file.open("samples.txt", std::ios_base::app);
+			float samples[numSamples][3];    	//Set up the variables
+							 	//necessary for sample generation
+			vec<float> d_s = in.cart2sph();
+			vec<float> temp;
+			int dim = (sqrt(numSamples)-1)/2;
+			//std::cout << dim << std::endl;
+			float y_0 	= 0.0;
+			float len 	= 10.0;
+			float p0  	= M_PI/3;
+			float dt  	= solidAngle/(1.0 * dim);
+			float dp  	= p0/(1.0*dim);
+	
+			
 			glBindFramebuffer(GL_FRAMEBUFFER, fboID);//set up GL buffer		
 			glFramebufferTexture2D(
 				GL_FRAMEBUFFER,
@@ -245,45 +239,11 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 			gluOrtho2D(0.0,2.0*len,0.0,numSamples*len);
 			glEnable(GL_TEXTURE_3D);
 			glBindTexture(GL_TEXTURE_3D, texID);
-		}
-		
-		void
-		Unbind()
-		{
-			//Finalize GL_buffer
-			glBindTexture(GL_TEXTURE_3D, 0);                      
-			glDisable(GL_TEXTURE_3D);
-			glBindFramebuffer(GL_FRAMEBUFFER,0);
-			glBindTexture(GL_TEXTURE_2D, 0);
-		}
-		/*
-		Method for populating the buffer with the sampled texture.
-		usually uses the default direction vector and then
-		*/
-		void
-		sampleDirection(int numSamples = 1089, float solidAngle = M_PI)
-		{
-			//ofstream file;
-			//file.open("samples.txt", std::ios_base::app);
-			float samples[numSamples][3];    	//Set up the variables
-							 	//necessary for sample generation
-			vec<float> d_s = direction.cart2sph();
-			vec<float> temp;
-			int dim = (sqrt(numSamples)-1)/2;
-			//std::cout << dim << std::endl;
-			float y_0 	= 0.0;
-			float p0  	= M_PI/3;
-			float dt  	= solidAngle/(1.0 * dim);
-			float dp  	= p0/(1.0*dim);
-	
-			Bind();	
-			//file << "Step: Direction" << "\n";
 			//file << "iteration: " << iter << "\n";
 			//file << "starting pos and dir:" << "[" << position[0] << "," <<position[1] << "," << position[2] << "]" << ":" << "["  << direction[0] << "," << direction[1] << "," << direction[2] << "]\n"; 
-			//Loop over the samples such that the original orientation is in the
-			//center of the resulting texture index (544)
-			//file << position[0] << "," <<position[1] << "," << position[2] << "\n"; 
+			//Loop over the samples
+			file << position[0] << "," <<position[1] << "," << position[2] << "\n"; 
 			int idx;
 			for(int i = -dim; i <= dim; i++){
 				for(int j = -dim; j <= dim; j++){
@@ -299,93 +259,28 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 					samples[idx][1] = temp[1];
 					samples[idx][2] = temp[2];
-			//		file << idx << ":"  <<"[" << samples[idx][0] << "," << samples[idx][1] << "," << samples[idx][2] << "]" << "\n";
+					//file << idx << ":"  <<"[" << samples[idx][0] << "," << samples[idx][1] << "," << samples[idx][2] << "]" << "\n";
-					UpdatePlanes(0.0, y_0+(idx)*10, temp); 
+					Update(0.0, y_0+(idx)*10, temp); 
 				}
 			}
-			Unbind();
+			//Finalize GL_buffer
+			glBindTexture(GL_TEXTURE_3D, 0);                      
+			glDisable(GL_TEXTURE_3D);
+			glBindFramebuffer(GL_FRAMEBUFFER,0);
+			glBindTexture(GL_TEXTURE_2D, 0);
 			int nxt = getCost();
-			stim::vec<float> next(samples[nxt][0], samples[nxt][1], samples[nxt][2]);
-			//next[0] = samples[nxt][0];
-			//next[1] = samples[nxt][1];
-			//next[2] = samples[nxt][2];
+			stim::vec<float> next;
+			next[0] = samples[nxt][0];
+			next[1] = samples[nxt][1];
+			next[2] = samples[nxt][2];
 			next.norm();	
 			//file << "next direction" << "[" << next[0] << "," << next[1] << "," << next[2] << "]\n\n"; 
-			setPosition(position[0] + next[0]*magnitude[0]/2,
-					 position[1]+next[1]*magnitude[0]/2,
-						position[2]+next[2]*magnitude[0]/2);
+			setPosition(position[0] + next[0]/500,
+					 position[1]+next[1]/500,
+						position[2]+next[2]/500);
 			setDirection(next[0], next[1], next[2]);
 			//file.close();
-			setDirection(next);
-			//file << "ending pos and dir:" << "[" << position[0] << "," <<position[1] << "," << position[2] << "]" << ":" << "["  << direction[0] << "," << direction[1] << "," << direction[2] << "]\n"; 
-		}
-
-		void
-		samplePosition(int numSamples = 1089, float delta = 0.2)
-		{
-			
-			float samples[numSamples][3];    	//Set up the variables
-			vec<float> temp;
-			int dim = (sqrt(numSamples)-1)/2;
-			float y_0 	= 0.0;			//location of the first sample in the buffer
-			stim::rect<float> samplingPlane =
-				 stim::rect<float>(magnitude[0]*delta, position, direction);
-			float step = 1.0/(dim);
-			Bind();	
-			//Loop over the samples, keeping the original position sample
-			//in the center of the resulting texture.
-
-			int idx;
-			for(int i = -dim; i <= dim; i++){
-				for(int j = -dim; j <= dim; j++){
-					//Create linear index
-					idx = (i+dim)*(dim*2+1) + (j+dim);
-
-					temp = samplingPlane.p(
-							0.5+step*i,
-								 0.5+step*j
-										);
-
-					samples[idx][0] = temp[0];	//save sample vector
-					samples[idx][1] = temp[1];
-					samples[idx][2] = temp[2];
-	
-					UpdatePlanes(0.0, y_0+(idx)*10, temp, 1); 
-				}
-			}
-			Unbind();
-			int nxt = getCost();
-			setPosition(samples[nxt][0], samples[nxt][1], samples[nxt][2]);
-		}
-
-		void
-		sampleMagnitude(int numSamples = 1089, float delta = 0.5)
-		{
-			
-			ofstream file;
-			file.open("samples.txt", std::ios_base::app);
-			float samples[numSamples];    	//Set up the variables
-			int dim = (sqrt(numSamples)-1)/2;
-			float y_0 	= 0.0;			//location of the first sample in the buffer
-			float min 	= 1.0-delta;
-			float max 	= 1.0+delta;
-			float step	= (max-min)/(numSamples-1);
-			float factor;
-			vec<float> temp;
-			Bind();	
-
-			for(int i = 0; i < numSamples; i++){
-					//Create linear index
-					factor = (min+step*i)*magnitude[0];
-					samples[i] = factor;	//save sample vector
-					stim::vec<float> temp(factor);
-					UpdatePlanes(0.0, y_0+(i)*10, temp, 2); 
-			}
-			Unbind();
-			int nxt = getCost();
-			setMagnitude(samples[nxt]);
-			file << position[0] << "," <<position[1] << "," << position[2] << "\n"; 
 		}
 //--------------------------------------------------------------------------//
@@ -433,36 +328,48 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		stim::rect<float> hor;
 		stim::rect<float> ver;	
-		//Default Constructor
+
 		gl_spider
 		()
 		{
 			setPosition(0.0,0.0,0.0);
-			setDirection(0.0,0.0,1.0);
-			setMagnitude(0.03);
+			setDirection(1.0,1.0,1.0);
+			setMagnitude(0.1,0.1);
+			//GenerateFBO(400,200);
+			//Update();
 		}
-		//temporary constructor for convenience, will be removed in further updates.	
+		gl_spider
+		(vec<float> pos, vec<float> dir, vec<float> mag)
+		{
+			position = pos;
+			direction = dir;
+			magnitude = mag;
+			//GenerateFBO(400,200);
+			//Update();
+		}
+		//temporary cost for convenience.	
 		gl_spider
 		(float pos_x, float pos_y, float pos_z, float dir_x, float dir_y, float dir_z,
-			float mag_x)
+			float mag_x, float mag_y)
 		{
 			setPosition(pos_x, pos_y, pos_z);
 			setDirection(dir_x, dir_y, dir_z);
-			setMagnitude(mag_x);
+			setMagnitude(mag_x, mag_y);
+			//GenerateFBO(400,200);
+			//Update();
 		}
-		//Attached the spider to the texture with the given GLuint ID.
-		//Samples in the default direction acting as the init method.	
 		void
-		attachSpider(GLuint id, int numSamples = 1089)
+		attachSpider(GLuint id)
 		{
 			texID = id;
-			GenerateFBO(20, numSamples*10);
-			sampleDirection();
+			Sample(direction);
+			//GenerateFBO(20,10000);
+		//	Update();
+		//	generateVectorField(direction, 4.0);
 		}
-		
-		//temporary Method necessary for visualization and testing.
+
 		void
 		Update()
 		{
@@ -473,37 +380,35 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 				((Y.cross(direction)).cross(direction)).norm());
 			ver = stim::rect<float>(magnitude, position, direction.norm(),
 				 hor.n());
+			//UpdateBuffer();
+			generateVectorField(direction, 4.0);
 		}
-		//Returns the position vector.
+		
 		vec<float>
 		getPosition()
 		{
 			return position;
 		}
-		//Returns the direction vector.
 		vec<float>
 		getDirection()
 		{
 			return direction;
 		}
-		//Returns the magnitude vector.
 		vec<float>
 		getMagnitude()
 		{
 			return magnitude;
 		}
-		//Sets the position vector to input vector pos
 		void
 		setPosition(vec<float> pos)
 		{
 			position = pos;
 		}
-		//Sets the position vector to the input float coordinates x,y,z
 		void
 		setPosition(float x, float y, float z)
 		{
@@ -512,15 +417,12 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 			position[2] = z;
 		}
-		
-		//Sets the direction vector to input vector dir
 		void
 		setDirection(vec<float> dir)
 		{
 			direction = dir;
 		}
-		//Sets the direction vector to the input float coordinates x,y,z
 		void
 		setDirection(float x, float y, float z)
 		{
@@ -529,22 +431,19 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 			direction[2] = z;
 		}
-		//Sets the magnitude vector to the input vector mag.	
 		void
 		setMagnitude(vec<float> mag)
 		{
-			magnitude[0] = mag[0];
-			magnitude[1] = mag[0];
+			magnitude = mag;
 		}
 		void
-		setMagnitude(float mag)
+		setMagnitude(float x, float y)
 		{
-			magnitude[0] = mag;
-			magnitude[1] = mag;
+			magnitude[0] = x;
+			magnitude[1] = y;
 		}
-		
-		//temporary method for visualization.
+	
 		GLuint
 		getFB()
 		{
@@ -554,9 +453,142 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 		void
 		Step()
 		{
+			std::cout << position[0] << "," << position[1] << "," << position[1]
+				<< std::endl;
+			//setPosition(direction*magnitude[1]/2+position);
 			findOptimalDirection();
+			//Update();
+			std::cout << position[0] << "," << position[1] << "," << position[1]
+				<< std::endl;
+			
+		}
+
+		void
+		UpdateBuffer()
+		{	
+			stim::vec<float>p1; 
+        	        stim::vec<float>p2; 
+	                stim::vec<float>p3; 
+                	stim::vec<float>p4;	
+			glBindFramebuffer(GL_FRAMEBUFFER, fboID);
+			glFramebufferTexture2D(
+				GL_FRAMEBUFFER,
+				GL_COLOR_ATTACHMENT0,
+				GL_TEXTURE_2D,
+				texbufferID,
+				0);
+			glBindFramebuffer(GL_FRAMEBUFFER, fboID);
+			GLenum DrawBuffers[1] = {GL_COLOR_ATTACHMENT0};
+			glDrawBuffers(1, DrawBuffers);
+			glBindTexture(GL_TEXTURE_2D, texbufferID);
+			glClearColor(0,0,0,0);
+			glClear(GL_COLOR_BUFFER_BIT);
+			glMatrixMode(GL_PROJECTION);
+			glLoadIdentity();
+			glMatrixMode(GL_MODELVIEW);
+			glLoadIdentity();
+			glViewport(0,0,400,200);
+			gluOrtho2D(0.0,2.0,0.0,2.0);
+			glEnable(GL_TEXTURE_3D);
+			glBindTexture(GL_TEXTURE_3D, texID);
+			p1 = hor.p(1,1);
+			p2 = hor.p(1,0);
+			p3 = hor.p(0,0);
+			p4 = hor.p(0,1);
+			glBegin(GL_QUADS);
+				glTexCoord3f(
+					p1[0],
+					p1[1],
+					p1[2]
+					);
+				glVertex2f(0.0,0.0);
+				glTexCoord3f(
+					p2[0],
+					p2[1],
+					p2[2]
+					);
+				glVertex2f(1.0, 0.0);
+				glTexCoord3f(
+					p3[0],
+					p3[1],
+					p3[2]
+					);
+				glVertex2f(1.0, 2.0);
+				glTexCoord3f(
+					p4[0],
+					p4[1],
+					p4[2]
+					);
+				glVertex2f(0.0, 2.0);
+			 glEnd();
+			 p1 = ver.p(1,1);
+			 p2 = ver.p(1,0);
+			 p3 = ver.p(0,0);
+			 p4 = ver.p(0,1);
+		 	 glBegin(GL_QUADS);
+				glTexCoord3f(
+					p1[0],
+					p1[1],
+					p1[2]
+					);
+				glVertex2f(1.0, 0.0);
+				glTexCoord3f(
+					p2[0],
+					p2[1],
+					p2[2]
+					);
+				glVertex2f(2.0, 0.0);
+				glTexCoord3f(
+					p3[0],
+					p3[1],
+					p3[2]
+					);
+				glVertex2f(2.0, 2.0);
+				glTexCoord3f(
+					p4[0],
+					p4[1],
+					p4[2]
+					);
+				glVertex2f(1.0, 2.0);
+			glEnd();
+			glBindTexture(GL_TEXTURE_3D, 0);                      
+			glDisable(GL_TEXTURE_3D);
+			glBindFramebuffer(GL_FRAMEBUFFER,0);
+			glBindTexture(GL_TEXTURE_2D, 0);
 		}
+		
+		void
+		generateVectorField(stim::vec<float> d, float dim)
+		{
+			vec<float> d_s = d.cart2sph();
+			vec<float> temp;
+			float Dim 	= (float) dim;
+			float y_0 	= 0.0;
+			float x_0 	= 0.0;
+			float len 	= 4.0/(2.0*Dim+1.0);
+			float t0  	= M_PI/2;
+			float p0  	= M_PI/3;
+			float dt  	= t0/Dim;
+			float dp  	= p0/Dim;
+			for(int i = -dim; i <= dim; i++){
+				for(int j = -dim; j <= dim; j++){
+					//field[i+dim][j+dim][0] = d[0];
+					//field[i+dim][j+dim][1] = d[1]+dt*i;
+					//field[i+dim][j+dim][2] = d[2]+dp*j;
+					temp[0] = 1;
+					temp[1] = d_s[1]+dt*i;
+					temp[2] = d_s[2]+dp*j;
+					temp = temp.sph2cart();
+					Update(x_0+2.0*(i+dim)*len, y_0+(j+dim)*len, temp); 
+				}
+			}
+			
+		}
+		                                                                   
+
+
+
 		void
 		findOptimalDirection()
@@ -564,9 +596,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 			/* Method for finding the best direction for the spider.
 			   Uses the camera to rotate. Then Calls Evaluate to find new cost.
 			*/
-			sampleDirection();
-			samplePosition();
-			sampleMagnitude();
+			Sample(direction);
 		}
 		/* Method for initializing the cuda devices, necessary only
@@ -92,20 +92,12 @@ public:
 		Y = directionY;
 	}
-	CUDA_CALLABLE rect(T size, vec<T, N> center, vec<T, N> directionX, vec<T, N> directionY )
+	CUDA_CALLABLE rect(vec<T,N> mag, vec<T, N> center, vec<T, N> directionX, vec<T, N> directionY )
 	{	
 		C = center;
 		X = directionX;
 		Y = directionY;
-		scale(size);
-	}
-
-	CUDA_CALLABLE rect(vec<T,N> size, vec<T, N> center, vec<T, N> directionX, vec<T, N> directionY )
-	{	
-		C = center;
-		X = directionX;
-		Y = directionY;
-		scale(size[0], size[1]);
+		scale(mag[0], mag[1]);
 	}
 	CUDA_CALLABLE void scale(T factor1, T factor2){