Merge branch 'master' of git.stim.ee.uh.edu:codebase/stimlib

cherub
2 parents ea2056c0 2ed641f0
Showing 21 changed files with 800 additions and 146 deletions Show diff stats
stim/biomodels/network.h
stim/cuda/cost.h
stim/cuda/timer.h
stim/envi/bil.h
stim/envi/binary.h
stim/envi/bip.h
stim/envi/bsq.h
stim/envi/envi.h
stim/gl/gl_spider.h
stim/grids/grid.h
stim/math/matrix.h
stim/math/plane.h
stim/math/quad.h
stim/math/rect.h
stim/math/spharmonics.h
stim/math/triangle.h
stim/math/vector.h
stim/ui/progressbar.h
stim/visualization/camera.h
stim/visualization/glObj.h
@@ -18,7 +18,7 @@ namespace stim{
 template<typename T>
 class network{
-	//helper classes
+	//-------------------HELPER CLASSES-----------------------
 	/// Stores information about a geometric point on the network centerline (including point position and radius)
 	//template<typename T>
 	class point : public stim::vec<T>{
@@ -42,10 +42,15 @@ class network{
 	/// Stores information about a single capillary (a length of vessel between two branch or end points)
 	//template<typename T>
-	class fiber{
+	class fiber : public std::list< point >{
+
+		using std::list< point >::begin;
+		using std::list< point >::end;
+		using std::list< point >::size;
+
 	public:
-		std::list< point > P;		//geometric point positions
+		//std::list< point > P;		//geometric point positions
 		typename std::list< t_node >::iterator n[2];				//indices to terminal nodes
 		unsigned int id;
@@ -60,9 +65,9 @@ class network{
 			//for each point
 			typename std::list< point >::iterator i;	//create a point iterator
-			for(i = P.begin(); i != P.end(); i++){		//for each point in the fiber
+			for(i = begin(); i != end(); i++){		//for each point in the fiber
-				if(i == P.begin())						//if this is the first point, just store it
+				if(i == begin())						//if this is the first point, just store it
 					p1 = *i;
 				else{									//if this is any other point
 					p0 = p1;							//shift p1->p0
@@ -84,9 +89,9 @@ class network{
 			//for each point
 			typename std::list< point >::iterator i;	//create a point iterator
-			for(i = P.begin(); i != P.end(); i++){		//for each point in the fiber
+			for(i = begin(); i != end(); i++){		//for each point in the fiber
-				if(i == P.begin()){						//if this is the first point, just store it
+				if(i == begin()){						//if this is the first point, just store it
 					p1 = *i;
 					r1 = i->r;
 				}
@@ -105,7 +110,30 @@ class network{
 			}
 			length = length_sum;						//store the total length
-			return radius_sum / length;					//return the average radius of the fiber
+
+			//if the total length is zero, store a radius of zero
+			if(length == 0)
+				return 0;
+			else
+				return radius_sum / length;					//return the average radius of the fiber
+		}
+
+		std::vector< stim::vec<T> > geometry(){
+
+			std::vector< stim::vec<T> > result;				//create an array to store the fiber geometry
+			result.resize( size() );					//pre-allocate the array
+
+			typename std::list< point >::iterator p;				//create a list iterator
+			unsigned int pi = 0;						//create an index into the result array
+
+			//for each geometric point on the fiber centerline
+			for(p = begin(); p != end(); p++){
+				result[pi] = *p;
+				pi++;
+			}
+
+			return result;			//return the geometry array
+
 		}
 		std::string str(){
@@ -113,7 +141,7 @@ class network{
 			//create an iterator for the point list
 			typename std::list<point>::iterator i;
-			for(i = P.begin(); i != P.end(); i++){
+			for(i = begin(); i != end(); i++){
 				ss<<i->str()<<"  r = "<<i->r<<std::endl;
 			}
@@ -172,7 +200,7 @@ class network{
 	};
-
+//---------------NETWORK CLASS-----------------------------
 protected:
@@ -231,6 +259,8 @@ public:
 	}
 	/// Load a network from an OBJ object
+
+	/// @param object is the object file to be used as the basis for the network
 	void load( stim::obj<T> object){
 		//get the number of vertices in the object
@@ -291,12 +321,98 @@ public:
 			for(unsigned int pi = 0; pi < nP; pi++){
 				point p = (point)L[pi];					//move the geometric coordinates into a point structure
 				p.r = R[pi][0];							//store the radius
-				f->P.push_back(p);						//push the point onto the current fiber
+				f->push_back(p);						//push the point onto the current fiber
 			}
 		}
 	}	//end load()
+	/// Returns an array of node positions
+	std::vector< stim::vec<T> > get_node_positions(){
+
+		std::vector< stim::vec<T> > result;				//create an array to store the result
+		result.resize(N.size());						//set the array size
+
+		t_node_i ni;									//create a terminal node iterator
+		unsigned int vi = 0;							//vertex index into the result array
+
+		//for every terminal node
+		for(ni = N.begin(); ni != N.end(); ni++){
+
+			//create a vector based on the node position
+
+			//if the number of outgoing nodes is nonzero
+			if(ni->out.size() != 0)
+				result[vi] = ni->out.front()->front();
+			else if(ni->in.size() != 0)
+				result[vi] = ni->in.front()->back();
+
+			vi++;										//increment the array index
+		}
+
+		//return the resulting array
+		return result;
+	}
+
+	std::vector< stim::vec<T> > get_fiber_geometry( fiber_i f ){
+		return f->geometry();
+	}
+
+	/// Generate an OBJ file from the network
+
+	stim::obj<T> obj(){
+
+		//create an OBJ object
+		stim::obj<T> object;
+
+		//name the nodes
+		set_names();
+
+		//retrieve a list of terminal node positions
+		std::vector< stim::vec<T> > node_pos = get_node_positions();
+
+		//add the nodes to the obj file
+		object.addV(node_pos);
+
+		//counter for vertex indices in the object class
+		unsigned int nP;
+
+		//for each fiber
+		fiber_i fi;						//create a fiber iterator
+		for(fi = F.begin(); fi != F.end(); fi++){
+
+			//get an array of fiber points
+			std::vector< stim::vec<T> > fiber_p = get_fiber_geometry(fi);
+
+			//create a subset of this array
+			typename std::vector< stim::vec<T> >::iterator start = fiber_p.begin() + 1;
+			typename std::vector< stim::vec<T> >::iterator end = fiber_p.end() - 1;
+			typename std::vector< stim::vec<T> > fiber_subset(start, end);
+
+			//add this subset to the geometry object
+			nP = object.addV(fiber_subset);
+
+			//create an array to hold vertex indices for a line
+			std::vector<unsigned int> line;
+			line.resize(fiber_p.size());
+
+			//add the terminal nodes to the line list (make sure to add 1 to make them compatible with the OBJ)
+			line[0] = fi->n[0]->id + 1;
+			line[line.size() - 1] = fi->n[1]->id + 1;
+
+			//add the intermediate vertex indices to the line array
+			for(unsigned int i = 0; i < fiber_subset.size(); i++){
+				line[1 + i] = nP + i;
+			}
+
+			//add the line list to the object class
+			object.addLine(line);
+
+		}
+
+		return object;
+	}
+
 	/// This function returns the information necessary for a simple graph-based physical (ex. fluid) simulation.
 	/// @param n0 is a array which will contain the list of source nodes
@@ -5,8 +5,9 @@
 #include <stdio.h>
 #include <stim/visualization/colormap.h>
 #include <sstream>
-#include <stim/math/mathvec.h>
-
+#include <stim/math/vector.h>
+#include <stim/cuda/devices.h>
+#include <stim/cuda/threads.h>
 ///Cost function that works with the gl-spider class to find index of the item with min-cost.
 typedef unsigned char uchar;
@@ -38,6 +39,7 @@ float get_sum(float *diff)
 	ret = cublasSetVector(20*10, sizeof(*diff), diff, 1, v_dif, 1);
 	float out;
 	ret = cublasSasum(handle, 20*10, v_dif, 1, &out);
+//	cublasDestroy(ret);
 	cublasDestroy(handle);
 	return out;
 }
@@ -90,7 +92,7 @@ void initArray(cudaGraphicsResource_t src, int DIM_Y)
 		cudaGraphicsMapResources(1, &src)	
 	);
 	HANDLE_ERROR(
-		cudaGraphicsSubResourceGetMappedArray(&srcArray, src,0,0)
+		cudaGraphicsSubResourceGetMappedArray(&srcArray, src, 0, 0)
 		);
 	HANDLE_ERROR(
 		cudaBindTextureToArray(texIn, srcArray)
@@ -109,9 +111,6 @@ void initArray(cudaGraphicsResource_t src, int DIM_Y)
 void cleanUP(cudaGraphicsResource_t src)
 {
 	HANDLE_ERROR(
-		cudaUnbindTexture(texIn)
-	);
-	HANDLE_ERROR(
 		cudaFree(result)
 	);
 	HANDLE_ERROR(
@@ -120,7 +119,13 @@ void cleanUP(cudaGraphicsResource_t src)
 	HANDLE_ERROR(
 		cudaFree(v_dif)
 	);
+	HANDLE_ERROR(
+		cudaUnbindTexture(texIn)
+	);
 }
+
+
+
 ///External access-point to the cuda function
 ///@param src, cudaGraphicsResource that handles the shared OpenGL/Cuda Texture
 ///@param DIM_Y, the number of samples in the template.
@@ -128,17 +133,33 @@ void cleanUP(cudaGraphicsResource_t src)
 extern "C"
 stim::vec<int> get_cost(cudaGraphicsResource_t src, int DIM_Y)
 {
+//	int minGridSize;
+//	int blockSize;
+
+//	cudaOccupancyMaxPotentialBlockSize(&minGridSize, &blockSize, get_diff, 0, 20*DIM_Y*10);
+//	std::cout << blockSize << std::endl;
+//	std::cout << minGridSize << std::endl;
+//	stringstream name;	//for debugging
+//	name << "Test.bmp";
+//	dim3 block(4,4);
+//	dim3 grid(20/4, DIM_Y*10/4);
+//	int gridSize = (DIM_Y*10*20 + 1024 - 1)/1024;
+//	dim3 grid(26, 26);
+//	dim3 grid = GenGrid1D(DIM_Y*10*20);
+//	stim::gpu2image<float>(result, name.str(), 20,DIM_Y*10,0,1);
+//	name.clear();
+//	name << "sample_" << inter << "_" << idx << ".bmp";
+//	stim::gpu2image<float>(v_dif, name.str(), 20,10,0,1);
+
 	float output[DIM_Y];
 	stim::vec<int> ret(0, 0);
 	float mini = 10000000000000000.0;
 	int idx;
-	stringstream name;	//for debugging
-	name << "Test.bmp";
 	initArray(src, DIM_Y*10);
-	dim3 grid(20, DIM_Y*10);
-	dim3 block(1, 1);
+	dim3 grid(20/2, DIM_Y*10/2);
+	dim3 block(2, 2);
+
 	get_diff <<< grid, block >>> (result);
-	stim::gpu2image<float>(result, name.str(), 20,DIM_Y*10,0,1);
 	for (int i = 0; i < DIM_Y; i++){
 		output[i] = get_sum(result+(20*10*i));
 		if(output[i] <= mini){
@@ -147,10 +168,7 @@ stim::vec&lt;int&gt; get_cost(cudaGraphicsResource_t src, int DIM_Y)
 		}
 	}	
-//	name.clear();
-//	name << "sample_" << inter << "_" << idx << ".bmp";
 	output[idx] = get_sum(result+(20*10*idx));
-//	stim::gpu2image<float>(v_dif, name.str(), 20,10,0,1);
 	cleanUP(src);
 	ret[0] = idx; ret[1] = (int) output[idx];
 	return ret;
+#ifndef STIM_CUDA_TIMER
+#define STIM_CUDA_TIMER
+
 static cudaEvent_t tStartEvent;
 static cudaEvent_t tStopEvent;
+namespace stim{
+
+/// These functions calculate the time between GPU functions in milliseconds
 static void gpuStartTimer()
 {
 	//set up timing events
@@ -18,4 +24,8 @@ static float gpuStopTimer()
 	cudaEventDestroy(tStartEvent);
 	cudaEventDestroy(tStopEvent);
 	return elapsedTime;
-}
 \ No newline at end of file
+}
+
+}	//end namespace stim
+
+#endif
 \ No newline at end of file
@@ -35,6 +35,8 @@ protected:
 		return R[1];
 	}
+	using binary<T>::thread_data;
+
 public:
 	using binary<T>::open;
@@ -318,12 +320,17 @@ public:
 			}//loop for YZ line end  
 			target.write(reinterpret_cast<const char*>(c), L);   //write the corrected data into destination	
+
+			thread_data = (double)k / Y() * 100;
 		}//loop for Y slice end		
 		free(a);
 		free(b);
 		free(c);
 		target.close();
+
+		thread_data = 100;
+
 		return true;	
 	}
@@ -366,11 +373,14 @@ public:
 				}								
 			}
 			target.write(reinterpret_cast<const char*>(c), L);   //write normalized data into destination
+
+			thread_data = (double)j / Y() * 100;
 		}
 		free(b);
 		free(c);
 		target.close();
+		thread_data = 100;
 		return true;
 	}
@@ -390,9 +400,13 @@ public:
 		for ( unsigned i = 0; i < Z(); i++)
 		{			
 				band_index(p, i);
-				target.write(reinterpret_cast<const char*>(p), S);   //write a band data into target file	
+				target.write(reinterpret_cast<const char*>(p), S);   //write a band data into target file
+
+				thread_data = (double)i / Z() * 100;	//store the progress for the current operation	
 		}
+		thread_data = 100;							//set the current progress to 100%
+
 		free(p);
 		target.close();
 		return true;
@@ -421,11 +435,15 @@ public:
 				unsigned ks = k * X();
 				for ( unsigned j = 0; j < X(); j++)
 					q[k + j * Z()] = p[ks + j];
+
+				thread_data = (double)(i * Z() + k) / (Z() * Y()) * 100;	//store the progress for the current operation	
 			}
 			target.write(reinterpret_cast<const char*>(q), S);   //write a band data into target file	
 		}
+		thread_data = 100;
+
 		free(p);
 		free(q);
@@ -825,40 +843,50 @@ public:
 		return true;
 	}
-	///Saves to disk only those spectra corresponding to mask values != 0
+	/// Saves to disk only those spectra corresponding to mask values != 0
+	/// Unlike the BIP and BSQ versions of this function, this version saves a different format: the BIL file is saved as a BIP
 	bool sift(std::string outfile, unsigned char* p){
 		// Assume X() = X, Y() = Y, Z() = Z.
 		std::ofstream target(outfile.c_str(), std::ios::binary);
 		//for loading pages:
-		unsigned XZ = X() * Z();		//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
+		unsigned long XZ = X() * Z();		//calculate the number of values in an XZ page on disk
+		unsigned long B = Z();			//calculate the number of bands
+		unsigned long L = XZ * sizeof(T);	//calculate the size of the page (in bytes)
-		//for saving spectra:
-		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, X() - 1, Y() - 1);	//creates a dummy spectrum by taking the last spectrum in the image
+		//allocate temporary memory for a XZ slice
+		T* slice = (T*) malloc(L);
-		for (unsigned i = 0; i < Y(); i++)			//Select a page by choosing Y coordinate, Y()
+		//allocate temporary memory for a spectrum
+		T* spec = (T*) malloc(B * sizeof(T));
+
+		//for each slice along the y axis
+		for (unsigned long y = 0; y < Y(); y++)			//Select a page by choosing Y coordinate, Y()
 		{
-			read_plane_y(temp, i);							//retrieve an ZX page, store in "temp"
-			for (unsigned j = 0; j < X(); j++)		//Select a pixel by choosing X coordinate in the page, X()
+			read_plane_y(slice, y);							//retrieve an ZX page, store in "slice"
+
+			//for each sample along X			
+			for (unsigned long x = 0; x < X(); x++)		//Select a pixel by choosing X coordinate in the page, X()
 			{
-				if (p[j * X() + i] != 0)					//if the mask != 0 at that XY pixel
+				//if the mask != 0 at that xy pixel
+				if (p[y * X() + x] != 0)					//if the mask != 0 at that XY pixel
 				{
-					for (unsigned k = 0; k < Z(); k++)		//Select a voxel by choosing Z coordinate at the pixel
+					//for each band at that pixel
+					for (unsigned long b = 0; b < B; b++)		//Select a voxel by choosing Z coordinate at the pixel
 					{
-						tempZ[k] = temp[k*X() + i];		//Pass the correct spectral value from XZ page into the spectrum to be saved.
+						spec[b] = slice[b*X() + x];		//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.
+					target.write((char*)spec, B * sizeof(T));		//write that spectrum to disk. Size is L2.
 				}
-				else
-					continue;
+
+				thread_data = (double) (y * X() + x) / (Y() * X()) * 100;
 			}
 		}
 		target.close();
-		free(temp);
+		free(slice);
+		free(spec);
+
+		thread_data = 100;
 		return true;
 	}
@@ -995,10 +1023,15 @@ public:
 			{
 				file.read((char *)(temp + j * sam), sizeof(T) * sam);
 				file.seekg(jumpb, std::ios::cur);    //go to the next band	
+
+				thread_data = (double)(i * Z() + j) / (lin * Z()) * 100;
 			}
 			out.write(reinterpret_cast<const char*>(temp), L);   //write slice data into target file	
 		}
 		free(temp);
+
+		thread_data = 100;
+
 		return true;
 	}
@@ -28,7 +28,7 @@ protected:
 	unsigned int header;	//header size (in bytes)
 	unsigned char* mask;	//pointer to a character array: 0 = background, 1 = foreground (or valid data)
-
+	unsigned int thread_data;		//unsigned integer used to pass data to threads during processing
 	/// Private initialization function used to set default parameters in the data structure.
@@ -36,6 +36,8 @@ protected:
 		memset(R, 0, sizeof(unsigned int) * D);		//initialize the resolution to zero
 		header = 0;									//initialize the header size to zero
 		mask = NULL;
+
+		thread_data = 0;
 	}
 	/// Private helper function that returns the size of the file on disk using system functions.
@@ -95,6 +97,14 @@ protected:
 public:
+	unsigned int get_thread_data(){
+		return thread_data;
+	}
+
+	void reset_thread_data(){
+		thread_data = 0;
+	}
+
 	/// Open a binary file for streaming.
 	/// @param filename is the name of the binary file
@@ -37,6 +37,8 @@ protected:
 		return R[0];
 	}
+	using binary<T>::thread_data;
+
 public:
 	using binary<T>::open;
@@ -327,6 +329,8 @@ public:
 			}//loop for YZ line end  
 			target.write(reinterpret_cast<const char*>(c), L);   //write the corrected data into destination	
+
+			thread_data = (double)k / Y() * 100;
 		}//loop for Y slice end
@@ -335,6 +339,8 @@ public:
 		free(b);
 		free(c);
 		target.close();
+
+		thread_data = 100;
 		return true;	
 	}
@@ -379,12 +385,15 @@ public:
 				}								
 			}
 			target.write(reinterpret_cast<const char*>(c), L);   //write normalized data into destination
+
+			thread_data = (double) j / Y() * 100;
 		}
 		free(b);
 		free(c);
 		target.close();
+		thread_data = 100;
 		return true;
 	}
@@ -434,10 +443,13 @@ public:
 				unsigned ks = k * X();
 				for ( unsigned j = 0; j < X(); j++)
 					q[ks + j] = p[k + j * Z()];
+
+				thread_data = (double)(i * Z() + k) / (Y() * Z()) * 100;
 			}
 			target.write(reinterpret_cast<const char*>(q), S);   //write a band data into target file	
 		}
+		thread_data = 100;
 		free(p);
 		free(q);
@@ -837,42 +849,113 @@ public:
 	/// Saves to disk only those spectra corresponding to mask values != 0
-	bool sift(std::string outfile, unsigned char* p){
-		// Assume X() = X, Y() = Y, Z() = Z.
+	bool sift(std::string outfile, unsigned char* mask){
+
+		//reset the file pointer to the beginning of the file
+		file.seekg(0, std::ios::beg);		
+
+		// open an output stream
 		std::ofstream target(outfile.c_str(), std::ios::binary);
-		//for loading pages:
-		unsigned ZX = X() * Z();		//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
+		//allocate space for a single spectrum
+		unsigned long B = Z();
+		T* spectrum = (T*) malloc(B * sizeof(T));
-		//for saving spectra:
-		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, X() - 1, Y() - 1);	//creates a dummy spectrum by taking the last spectrum in the image
+		//calculate the number of pixels in a band
+		unsigned long XY = X() * Y();
-		for (unsigned i = 0; i < Y(); i++)			//Select a page by choosing Y coordinate, Y()
-		{
-			read_plane_y(temp, i);							//retrieve an ZX page, store in "temp"
-			for (unsigned j = 0; j < X(); j++)		//Select a pixel by choosing X coordinate in the page, X()
-			{
-				if (p[j * X() + i] != 0)					//if the mask != 0 at that XY pixel
-				{
-					for (unsigned k = 0; k < Z(); k++)		//Select a voxel by choosing Z coordinate at the pixel
-					{
-						tempZ[k] = temp[i*Z() + 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;
+		//for each pixel
+		unsigned long skip = 0;					//number of spectra to skip
+		for(unsigned long x = 0; x < XY; x++){
+
+			//if the current pixel isn't masked
+			if( mask[x] == 0){
+				//increment the number of skipped pixels
+				skip++;
+			}
+			//if the current pixel is masked
+			else{
+
+				//skip the intermediate pixels
+				file.seekg(skip * B * sizeof(T), std::ios::cur);
+
+				//set the skip value to zero
+				skip = 0;
+
+				//read this pixel into memory
+				file.read((char *)spectrum, B * sizeof(T));
+
+				//write this pixel out
+				target.write((char *)spectrum, B * sizeof(T));
+
+				thread_data = (double) x / XY * 100;
 			}
+
 		}
+
+		//close the output file
 		target.close();
-		free(temp);
+		free(spectrum);
+
+		thread_data = 100;
+
 		return true;
 	}
+	bool unsift(std::string outfile, unsigned char* mask, unsigned int samples, unsigned int lines){
+
+		// open an output stream
+		std::ofstream target(outfile.c_str(), std::ios::binary);
+
+		//reset the file pointer to the beginning of the file
+		file.seekg(0, std::ios::beg);		
+
+		//allocate space for a single spectrum
+		unsigned long B = Z();
+		T* spectrum = (T*) malloc(B * sizeof(T));
+
+		//allocate space for a spectrum of zeros
+		T* zeros = (T*) malloc(B * sizeof(T));
+		memset(zeros, 0, B * sizeof(T));
+
+		//calculate the number of pixels in a band
+		unsigned long XY = samples * lines;
+
+		//for each pixel
+		unsigned long skip = 0;					//number of spectra to skip
+		for(unsigned long x = 0; x < XY; x++){
+
+			//if the current pixel isn't masked
+			if( mask[x] == 0){
+
+				//write a bunch of zeros
+				target.write((char *)zeros, B * sizeof(T));
+			}
+			//if the current pixel is masked
+			else{
+
+				//read a pixel into memory
+				file.read((char *)spectrum, B * sizeof(T));
+
+				//write this pixel out
+				target.write((char *)spectrum, B * sizeof(T));
+			}
+
+			thread_data = (double)x / XY * 100;
+
+		}
+
+		//close the output file
+		target.close();
+		free(spectrum);
+
+		thread_data = 100;
+
+		return true;
+
+
+	}
+
 	/// @param p is a pointer to memory of size X * Y * sizeof(T) that will store the band averages.
@@ -995,9 +1078,13 @@ public:
 			{
 				pixel(temp, sp + j + i * X());
 				out.write(reinterpret_cast<const char*>(temp), L);   //write slice data into target file	
+
+				thread_data = (double)(i * sam + j) / (lin * sam) * 100;
 			}
 		}
 		free(temp);
+
+		thread_data = 100;
 		return true;
 	}
@@ -42,6 +42,8 @@ protected:
 		return R[2];
 	}
+	using binary<T>::thread_data;
+
 public:
 	using binary<T>::open;
@@ -259,6 +261,8 @@ public:
 			}
 			target.write(reinterpret_cast<const char*>(c), S);   //write the corrected data into destination
+
+			thread_data = (double)cii / B * 100;
 		}	
@@ -268,6 +272,8 @@ public:
 		free(b);
 		free(c);
 		target.close();
+
+		thread_data = 100;
 		return true;
 	}
@@ -304,13 +310,18 @@ public:
 					c[i] = c[i] / b[i];
 			}
 			target.write(reinterpret_cast<const char*>(c), S);   //write normalized data into destination
+
+			thread_data = (double)j / B * 100;
 		}
+		
+
 		//header.save(headername);         //save the new header file
 		free(b);
 		free(c);
 		target.close();
+		thread_data = 100;					//make sure that the progress bar is full
 		return true;
 	}
@@ -349,24 +360,31 @@ public:
 		std::ofstream target(outname.c_str(), std::ios::binary);
 		std::string headername = outname + ".hdr";
-		T * p;			//pointer to the current spectrum
-		p = (T*)malloc(L);
-		
-		for ( unsigned i = 0; i < Y(); i++)
+		unsigned int xz_bytes = X() * Z() * sizeof(T);
+		T * xz_slice;						//pointer to the current spectrum
+		xz_slice = (T*)malloc(xz_bytes);
+
+		for ( unsigned y = 0; y < Y(); y++)									//for each y position
 		{
-			file.seekg(X() * i * sizeof(T), std::ios::beg);
-			for ( unsigned j = 0; j < Z(); j++ )
+			file.seekg(y * X() * sizeof(T), std::ios::beg);					//seek to the beginning of the xz slice
+			for ( unsigned z = 0; z < Z(); z++ )							//for each band
 			{
-				file.read((char *)(p + j * X()), sizeof(T) * X());
-				file.seekg(jump, std::ios::cur);    //go to the next band	
+				file.read((char *)(xz_slice + z * X()), sizeof(T) * X());	//read a line
+				file.seekg(jump, std::ios::cur);							//seek to the next band
+
+				
+				thread_data = (double)(y * Z() + z) / (Z() * Y()) * 100;	//update the progress counter
 			}					
-			target.write(reinterpret_cast<const char*>(p), L);   //write XZ slice data into target file	
+			target.write(reinterpret_cast<const char*>(xz_slice), xz_bytes);   //write the generated XZ slice to the target file	
 		}
-		//header.interleave = rts::envi_header::BIL;  //change the type of file in header file
-		//header.save(headername);
+
+
+		
+		thread_data = 100;
-		free(p);
+		free(xz_slice);
 		target.close();
+
 		return true;
 	}
@@ -763,16 +781,16 @@ public:
 	bool sift(std::string outfile, unsigned char* p){
 		std::ofstream target(outfile.c_str(), std::ios::binary);
 		// open a band (XY plane)
-		unsigned XY = X() * Y(); //Number of XY pixels
-		unsigned L = XY * sizeof(T); //size of XY pixels
+		unsigned long XY = X() * Y(); //Number of XY pixels
+		unsigned long 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 < Z(); i++)			//for each spectral bin
+		for (unsigned long i = 0; i < Z(); 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
+			for (unsigned long 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];
@@ -781,10 +799,15 @@ public:
 				else{
 					continue;
 				}
+
+				thread_data = (double)(i * XY + j) / (XY * Z()) * 100;
 			}
 		}
 		target.close();
 		free(temp);
+
+		thread_data = 100;
+
 		return true;
 	}
@@ -803,9 +826,9 @@ public:
 		std::cout<<"started sifting"<<std::endl;
 		//get the number of pixels and bands in the input image
-		unsigned int P = X(); 					//Number of pixels
-		unsigned int B = Z();					//number of bands
-		unsigned int XY = samples * lines;		//total number of pixels in an unsifted image
+		unsigned long P = X(); 					//Number of pixels
+		unsigned long B = Z();					//number of bands
+		unsigned long XY = samples * lines;		//total number of pixels in an unsifted image
 		// allocate memory for a sifted band
 		T * sifted = (T*)malloc(P * sizeof(T)); //allocate memory for a band
@@ -814,28 +837,33 @@ public:
 		T* unsifted = (T*) malloc(XY * sizeof(T));
 		//for each band
-		for(unsigned int b = 0; b < B; b++){
+		for(unsigned long b = 0; b < B; b++){
 			//set the unsifted index value to zero
-			unsigned int i = 0;
+			unsigned long i = 0;
 			//retrieve the sifted band (masked pixels only)
 			band_index(sifted, b);
 			//for each pixel in the final image (treat it as a 1D image)
-			for(unsigned int xi = 0; xi < XY; xi++){
+			for(unsigned long xi = 0; xi < XY; xi++){
 				if( p[xi] == 0 )
 					unsifted[xi] = 0;
 				else{
 					unsifted[xi] = sifted[i];
 					i++;
 				}
+				//std::cout<<xi<<"/"<<XY<<",   "<<b<<"/"<<B<<std::endl;
+				thread_data = (double)(b * XY + xi) / (B * XY) * 100;
 			}
-
+			//std::cout<<b*XY<<"/"<<B*XY<<"      "<<B<<"-"<<XY<<std::endl;
 			//write the band image to disk
 			target.write(reinterpret_cast<const char*>(unsifted), sizeof(T) * XY);
 		}
+		//std::cout<<"unsifted"<<std::endl;
+		thread_data = 100;
+
 		return true;
 	}
@@ -953,17 +981,22 @@ public:
 		unsigned long long jumpl = (X() - sam) * sizeof(T);	                           	//jump pointer to the next line
 		//get start
 		file.seekg((y0 * X() + x0) * sizeof(T), std::ios::beg);
-		for (unsigned i = 0; i < Z(); i++)
+		for (unsigned z = 0; z < Z(); z++)
 		{
 			for (unsigned j = 0; j < lin; j++)
 			{
 				file.read((char *)(temp + j * sam), sizeof(T) * sam);
 				file.seekg(jumpl, std::ios::cur);    //go to the next band	
+
+				thread_data = (double)(z * lin + j) / (Z() * lin) * 100;
 			}
 			out.write(reinterpret_cast<const char*>(temp), L);   //write slice data into target file	
 			file.seekg(jumpb, std::ios::cur);
 		}
 		free(temp);
+
+		thread_data = 100;
+
 		return true;
 	}
@@ -6,7 +6,7 @@
 #include "../envi/bip.h"
 #include "../envi/bil.h"
 #include <iostream>
-#include "../image/image.h"
+//#include "../image/image.h"
 namespace stim{
@@ -23,6 +23,78 @@ public:
 	envi_header header;
+	/// Returns the progress of the current processing operation as a percentage
+	void reset_progress(){
+
+		if(header.interleave == envi_header::BSQ){		//if the infile is bsq file
+			if(header.data_type ==envi_header::float32)
+				((bsq<float>*)file)->reset_thread_data();
+			else if(header.data_type == envi_header::float64)
+				((bsq<double>*)file)->reset_thread_data();
+			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)
+				((bil<float>*)file)->reset_thread_data();
+			else if(header.data_type == envi_header::float64)
+				((bil<double>*)file)->reset_thread_data();
+			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)
+				((bip<float>*)file)->reset_thread_data();
+			else if(header.data_type == envi_header::float64)
+				((bip<double>*)file)->reset_thread_data();
+			else
+				std::cout<<"ERROR: unidentified data type"<<std::endl;
+		}
+
+		else{
+			std::cout<<"ERROR: unidentified file type"<<std::endl;
+			exit(1);
+		}
+	}
+
+	/// Returns the progress of the current processing operation as a percentage
+	unsigned int progress(){
+
+		if(header.interleave == envi_header::BSQ){		//if the infile is bsq file
+			if(header.data_type ==envi_header::float32)
+				return ((bsq<float>*)file)->get_thread_data();
+			else if(header.data_type == envi_header::float64)
+				return ((bsq<double>*)file)->get_thread_data();
+			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)->get_thread_data();
+			else if(header.data_type == envi_header::float64)
+				return ((bil<double>*)file)->get_thread_data();
+			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)->get_thread_data();
+			else if(header.data_type == envi_header::float64)
+				return ((bip<double>*)file)->get_thread_data();
+			else
+				std::cout<<"ERROR: unidentified data type"<<std::endl;
+		}
+
+		else{
+			std::cout<<"ERROR: unidentified file type"<<std::endl;
+		}
+		return 0;
+	}
+
 	/// Allocate memory for a new ENVI file based on the current interleave format (BIP, BIL, BSQ) and data type.
 	bool allocate(){
@@ -385,6 +457,10 @@ public:
 		//create a new header
 		envi_header new_header = header;
+		//if a BIL file is sifted, it's saved as a BIP
+		if(header.interleave == envi_header::BIL)
+			new_header.interleave = envi_header::BIP;		
+
 		//set the number of lines to 1 (this is a matrix with 1 line and N samples)
 		new_header.lines = 1;
 		new_header.samples = nnz;
@@ -455,14 +531,18 @@ public:
 		else if (header.interleave == envi_header::BIP){		//if the infile is bip file
-				std::cout << "ERROR in stim::envi::unsift - BIP files aren't supported yet" << std::endl;
+			if (header.data_type == envi_header::float32)
+				return ((bip<float>*)file)->unsift(outfile, mask, samples, lines);
+			else if (header.data_type == envi_header::float64)
+				return ((bip<double>*)file)->unsift(outfile, mask, samples, lines);
+			else
+				std::cout << "ERROR: unidentified data type" << std::endl;
 		}
 		else{
 			std::cout << "ERROR: unidentified file type" << std::endl;
-			exit(1);
 		}
-
+		return 0;
 	}
 	/// Compute the ratio of two baseline-corrected peaks. The result is stored in a pre-allocated array.
@@ -847,7 +927,7 @@ public:
 	/// @param mask is a pointer to pre-allocated memory of size X*Y
 	/// @param maskname is the file name for the image that will serve as the mask
-	bool load_mask(unsigned char * mask, std::string maskname){
+	/*bool load_mask(unsigned char * mask, std::string maskname){
 		//open the mask file
 		stim::image<unsigned char> mask_image(maskname);
 		mask_image.data_noninterleaved(mask);
@@ -855,7 +935,7 @@ public:
 		//memcpy(mask, mask_image.data_noninterleaved(), mask_image.size());
 		//mask_image.clear();
 		return true;
-	}
+	}*/
 	/// Calculate the mean value for all masked (or valid) pixels in a band and returns the average spectrum
@@ -10,11 +10,12 @@
 #include "stim/gl/gl_texture.h"
 #include "stim/visualization/camera.h"
 #include "stim/gl/error.h"
-#include "stim/math/mathvec.h"
+#include "stim/math/vector.h"
 #include "stim/math/rect.h"
 #include "stim/math/matrix.h"
 #include "stim/cuda/cost.h"
 #include "stim/cuda/glbind.h"
+#include <stim/visualization/obj.h>
 #include <vector>
 #include <iostream>
@@ -55,6 +56,7 @@ class gl_spider
 		GLuint texbufferID;
 		int numSamples;
 		float stepsize = 3.0;
+		int current_cost;
 		/// Method for finding the best scale for the spider.
 		/// changes the x, y, z size of the spider to minimize the cost
@@ -119,7 +121,7 @@ class gl_spider
 			setMatrix();
 			glCallList(dList+3);
-			int best = getCost();
+//			int best = getCost();
 		}
@@ -142,7 +144,7 @@ class gl_spider
 		///Method for populating the vector arrays with sampled vectors.
 		///uses the default d vector <0,0,1>
 		void
-		genDirectionVectors(float solidAngle = M_PI)
+		genDirectionVectors(float solidAngle = 5*M_PI/4)
 		{
 			//ofstream file;
 			//file.open("dvectors.txt");
@@ -492,8 +494,9 @@ class gl_spider
 			createResource();
 			stim::vec<int> cost = 	get_cost(resource, numSamples);
 			destroyResource();
-			if (cost[1] >= 80)
-				exit(0);
+//			if (cost[1] >= 80)
+//				exit(0);
+			current_cost = cost[1];
 			return cost[0];
 		}
@@ -562,7 +565,7 @@ class gl_spider
 			dList = glGenLists(3);
 			glListBase(dList);
 			Bind();
-			genDirectionVectors(M_PI);
+			genDirectionVectors(5*M_PI/4);
 			genPositionVectors();
 			genMagnitudeVectors();
 			DrawCylinder();
@@ -649,7 +652,7 @@ class gl_spider
 		{
 			m[0] = mag;
 			m[1] = mag;
-			m[2] = mag;
+	//		m[2] = mag;
 		}
@@ -719,15 +722,16 @@ class gl_spider
 		}
-		void
+		int
 		Step()
 		{
 			Bind();
 			findOptimalDirection();
 			findOptimalPosition();
 			findOptimalScale();
-//			branchDetection();
+	//		branchDetection();
 			Unbind();
+			return current_cost;
 		}
@@ -7,7 +7,7 @@
 #include <fstream>
 #include <cstdarg>
-#include "../math/mathvec.h"
+#include <stim/math/vector.h>
 namespace stim{
@@ -4,7 +4,7 @@
 //#include "rts/vector.h"
 #include <string.h>
 #include <iostream>
-#include "mathvec.h"
+#include <stim/math/vector.h>
 #include "../cuda/callable.h"
 namespace stim{
@@ -2,7 +2,7 @@
 #define RTS_PLANE_H
 #include <iostream>
-#include "../math/vector.h"
+#include <stim/math/vector.h>
 #include "rts/cuda/callable.h"
@@ -2,10 +2,10 @@
 #define RTS_QUAD_H
 //enable CUDA_CALLABLE macro
-#include "../cuda/callable.h"
-#include "../math/vector.h"
-#include "../math/triangle.h"
-#include "../math/quaternion.h"
+#include <stim/cuda/callable.h>
+#include <stim/math/vector.h>
+#include <stim/math/triangle.h>
+#include <stim/math/quaternion.h>
 #include <iostream>
 #include <iomanip>
 #include <algorithm>
@@ -2,10 +2,10 @@
 #define RTS_RECT_H
 //enable CUDA_CALLABLE macro
-#include "../cuda/callable.h"
-#include "../math/mathvec.h"
-#include "../math/triangle.h"
-#include "../math/quaternion.h"
+#include <stim/cuda/callable.h>
+#include <stim/math/vector.h>
+#include <stim/math/triangle.h>
+#include <stim/math/quaternion.h>
 #include <iostream>
 #include <iomanip>
 #include <algorithm>
 #ifndef STIM_SPH_HARMONICS
 #define STIM_SPH_HARMONICS
-#include <stim/math/mathvec.h>
+#include <stim/math/vector.h>
 #include <boost/math/special_functions/spherical_harmonic.hpp>
 #include <vector>
@@ -2,8 +2,8 @@
 #define RTS_TRIANGLE_H
 //enable CUDA_CALLABLE macro
-#include "../cuda/callable.h"
-#include "../math/mathvec.h"
+#include <stim/cuda/callable.h>
+#include <stim/math/vector.h>
 #include <iostream>
 namespace stim{
@@ -59,7 +59,7 @@ struct vec : public std::vector&lt;T&gt;
 	//copy constructor
 	vec( const vec<T>& other){
 		unsigned int N = other.size();
-		for(int i=0; i<N; i++)
+		for(unsigned int i=0; i<N; i++)
 			push_back(other[i]);
 	}
@@ -112,7 +112,7 @@ struct vec : public std::vector&lt;T&gt;
         //compute and return the vector length
         T sum_sq = (T)0;
-        for(int i=0; i<N; i++)
+        for(unsigned int i=0; i<N; i++)
         {
             sum_sq += pow( at(i), 2 );
         }
@@ -231,7 +231,7 @@ struct vec : public std::vector&lt;T&gt;
         vec<T> result(N);
-        for(int i=0; i<N; i++)
+        for(unsigned int i=0; i<N; i++)
             result[i] = at(i) - rhs[i];
         return result;
@@ -340,7 +340,7 @@ struct vec : public std::vector&lt;T&gt;
 		unsigned int N = size();
 		ss<<"[";
-		for(int i=0; i<N; i++)
+		for(unsigned int i=0; i<N; i++)
 		{
 			ss<<at(i);
 			if(i != N-1)
-#ifndef RTS_PROGRESSBAR_H
-#define RTS_PROGRESSBAR_H
-
+#ifndef RTS_PROGRESSBAR_H
+#define RTS_PROGRESSBAR_H
+
 #include <iostream>
 #include <sstream>
 using namespace std;
-static void rtsProgressBar(int percent)
+static void rtsProgressBar(unsigned int percent)
 {
+	//std::cout<<percent<<std::endl;
 	stringstream bar;
 	static int x = 0;
 	string slash[4];
@@ -25,9 +26,10 @@ static void rtsProgressBar(int percent)
 	bar<<"]";
 	cout << "\r"; // carriage return back to beginning of line
 	cout << bar.str() << " " << slash[x] << " " << percent << " %"; // print the bars and percentage
+	cout.flush();
 	x++; // increment to make the slash appear to rotate
 	if(x == 4)
 	x = 0; // reset slash animation
-}
-
-#endif
+}
+
+#endif
-#include "../math/mathvec.h"
-#include "../math/quaternion.h"
-#include "../math/matrix.h"
+#include <stim/math/vector.h>
+#include <stim/math/quaternion.h>
+#include <stim/math/matrix.h>
 #include <ostream>
+#ifndef STIM_GLOBJ_H
+#define STIM_GLOBJ_H
+
+#include <stim/visualization/obj.h>
+#include <GL/glew.h>
+#include <GL/glut.h>
+#include <stim/math/vector.h>
+
+
+namespace stim
+{
+/** This class uses the loading functionality of the obj class in order to
+ *  Assist with the visualization the segmented vessels.
+ *  Uses
+*/
+
+class objGl : public virtual stim::obj<T>
+{
+private:
+	using stim::obj<T>::L;
+	using stim::obj<T>::P;
+	using stim::obj<T>::F;
+	using vec<T>::size;