added an inverse PCA transform to the stim::envi classes

David Mayerich
1 parent 371bf339
Showing 5 changed files with 491 additions and 385 deletions Show diff stats
stim/envi/bil.h
stim/envi/binary.h
stim/envi/bip.h
stim/envi/bsq.h
stim/envi/envi.h
@@ -35,7 +35,12 @@ protected:
 		return R[1];
 	}
-	using binary<T>::thread_data;
+	using binary<T>::progress;
+
+	/// Call the binary nnz() function for the BIL orientation
+	unsigned long long nnz(unsigned char* mask){
+		return binary<T>::nnz(mask, X()*Y());
+	}
 public:
@@ -63,9 +68,10 @@ public:
 	/// @param p is a pointer to an allocated region of memory at least X * Y * sizeof(T) in size.
 	/// @param page <= B is the integer number of the band to be copied.
-	bool band_index( T * p, unsigned int page){
+	bool band_index( T * p, unsigned int page, bool PROGRESS = false){
 		//return binary<T>::read_plane_1(p, page);
+		if(PROGRESS) progress = 0;
 		unsigned int L = X() * sizeof(T);		//caculate the number of bytes in a sample line
 		unsigned int jump = X() * (Z() - 1) * sizeof(T);
@@ -79,6 +85,7 @@ public:
 		{
 			file.read((char *)(p + i * X()), L);
 			file.seekg( jump, std::ios::cur);
+			if(PROGRESS) progress = (double)(i+1) / Y() * 100;
 		}
 		return true;
@@ -88,7 +95,7 @@ public:
 	/// @param p is a pointer to an allocated region of memory at least X * Y * sizeof(T) in size.
 	/// @param wavelength is a floating point value (usually a wavelength in spectral data) used as a label for the band to be copied.
-	bool band( T * p, double wavelength){
+	bool band( T * p, double wavelength, bool PROGRESS = false){
 		//if there are no wavelengths in the BSQ file
 		if(w.size() == 0)
@@ -104,7 +111,7 @@ public:
 		//if wavelength is smaller than the first one in header file
 		if ( w[page] > wavelength ){
-			band_index(p, page);
+			band_index(p, page, PROGRESS);
 			return true;
 		}
@@ -124,7 +131,7 @@ public:
 			p1=(T*)malloc(S);                     //memory allocation
 			p2=(T*)malloc(S);
 			band_index(p1, page - 1);
-			band_index(p2, page );
+			band_index(p2, page, PROGRESS);
 			for(unsigned i=0; i < XY; i++){
 				double r = (double) (wavelength - w[page-1]) / (double) (w[page] - w[page-1]);
 				p[i] = (p2[i] - p1[i]) * r + p1[i];
@@ -134,7 +141,7 @@ public:
 		}
 		else                           //if the wavelength is equal to a wavelength in header file
 		{
-			band_index(p, page);
+			band_index(p, page, PROGRESS);
 		}
 		return true;
@@ -195,8 +202,8 @@ public:
 	/// @param p is a pointer to pre-allocated memory at least B * sizeof(T) in size.
 	/// @param x is the x-coordinate (dimension 1) of the spectrum.
 	/// @param y is the y-coordinate (dimension 2) of the spectrum.
-	bool spectrum(T * p, unsigned x, unsigned y){
-		return binary<T>::read_line_02(p, x, y);
+	bool spectrum(T * p, unsigned x, unsigned y, bool PROGRESS = false){
+		return binary<T>::read_line_1(p, x, y, PROGRESS);
 	}
 	/// Retrieve a single pixel and stores it in pre-allocated memory.
@@ -223,7 +230,7 @@ public:
 	/// @param outname is the name of the output file used to store the resulting baseline-corrected data.
 	/// @param wls is the list of baseline points based on band labels.
-	bool baseline(std::string outname, std::vector<double> wls){
+	bool baseline(std::string outname, std::vector<double> wls, bool PROGRESS = false){
 		unsigned N = wls.size();			//get the number of baseline points
@@ -321,7 +328,7 @@ 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;
+			if(PROGRESS) progress = (double)(k+1) / Y() * 100;
 		}//loop for Y slice end		
 		free(a);
@@ -329,8 +336,6 @@ public:
 		free(c);
 		target.close();
-		thread_data = 100;
-
 		return true;	
 	}
@@ -340,7 +345,7 @@ public:
 	/// @param outname is the name of the output file used to store the resulting baseline-corrected data.
 	///	@param w is the label specifying the band that the hyperspectral image will be normalized to.
 	///	@param t is a threshold specified such that a spectrum with a value at w less than t is set to zero. Setting this threshold allows the user to limit division by extremely small numbers.
-	bool normalize(std::string outname, double w, double t = 0.0)
+	bool normalize(std::string outname, double w, double t = 0.0, bool PROGRESS = false)
 	{
 		unsigned int B = Z();		//calculate the number of bands
 		unsigned int ZX = Z() * X();
@@ -374,20 +379,19 @@ public:
 			}
 			target.write(reinterpret_cast<const char*>(c), L);   //write normalized data into destination
-			thread_data = (double)j / Y() * 100;
+			if(PROGRESS) progress = (double)(j+1) / Y() * 100;
 		}
 		free(b);
 		free(c);
 		target.close();
-		thread_data = 100;
 		return true;
 	}
 	/// Convert the current BIL file to a BSQ file with the specified file name.
 	/// @param outname is the name of the output BSQ file to be saved to disk.
-	bool bsq(std::string outname)
+	bool bsq(std::string outname, bool PROGRESS = false)
 	{
 		unsigned int S = X() * Y() * sizeof(T);		//calculate the number of bytes in a band
@@ -402,11 +406,9 @@ public:
 				band_index(p, i);
 				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	
+				if(PROGRESS) progress = (double)(i+1) / Z() * 100;	//store the progress for the current operation	
 		}
-		thread_data = 100;							//set the current progress to 100%
-
 		free(p);
 		target.close();
 		return true;
@@ -415,7 +417,7 @@ public:
 	/// Convert the current BIL file to a BIP file with the specified file name.
 	/// @param outname is the name of the output BIP file to be saved to disk.
-	bool bip(std::string outname)
+	bool bip(std::string outname, bool PROGRESS = false)
 	{
 		unsigned int S = X() * Z() * sizeof(T);		//calculate the number of bytes in a ZX slice
@@ -436,15 +438,12 @@ public:
 				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	
+				if(PROGRESS) progress = (double)((i+1) * Z() + k+1) / (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);
 		target.close();
@@ -879,14 +878,14 @@ public:
 					target.write((char*)spec, B * sizeof(T));		//write that spectrum to disk. Size is L2.
 				}
-				thread_data = (double) (y * X() + x) / (Y() * X()) * 100;
+				progress = (double) (y * X() + x) / (Y() * X()) * 100;
 			}
 		}
 		target.close();
 		free(slice);
 		free(spec);
-		thread_data = 100;
+		progress = 100;
 		return true;
 	}
@@ -923,7 +922,7 @@ public:
 	/// @param p is a pointer to pre-allocated memory of size [B * sizeof(T)] that stores the mean spectrum
 	/// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
-	bool avg_band(double* p, unsigned char* mask = NULL){
+	bool avg_band(double* p, unsigned char* mask = NULL, bool PROGRESS = false){
 		unsigned long long XZ = X() * Z();
 		unsigned long long XY = X() * Y();
 		T* temp = (T*)malloc(sizeof(T) * XZ);
@@ -947,6 +946,7 @@ public:
 					}
 				}
 			}
+			if(PROGRESS) progress = (double)(k+1) / Y() * 100;
 		}
 		free(temp);
 		return true;
@@ -957,8 +957,8 @@ public:
 	/// @param co is a pointer to pre-allocated memory of size [B * B] that stores the resulting covariance matrix
 	/// @param avg is a pointer to memory of size B that stores the average spectrum
 	/// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
-	bool co_matrix(double* co, double* avg, unsigned char *mask){
-		thread_data = 0;
+	bool co_matrix(double* co, double* avg, unsigned char *mask, bool PROGRESS = false){
+		progress = 0;
 		//memory allocation
 		unsigned long long xy = X() * Y();
 		unsigned long long B = Z();
@@ -986,7 +986,7 @@ public:
 					}
 				}
 			}
-			thread_data = (double)j / xy * 100;
+			if(PROGRESS) progress = (double)(j+1) / xy * 100;
 		}
 		//because correlation matrix is symmetric
 		for (unsigned long long i = 0; i < B; i++){
@@ -996,9 +996,6 @@ public:
 		}
 		free(temp);
-
-		thread_data = 100;							//processing complete
-
 		return true;
 	}
@@ -1015,7 +1012,8 @@ public:
 								   unsigned long long x1,
 								   unsigned long long y1,
 								   unsigned long long b0,
-								   unsigned long long b1){
+								   unsigned long long b1,
+								   bool PROGRESS = false){
 		//calculate the new image parameters
 		unsigned long long samples = x1 - x0;
@@ -1048,7 +1046,7 @@ public:
 				file.read((char *)(temp + z * samples), sizeof(T) * samples);
 				file.seekg(jumpb, std::ios::cur);    //go to the next band	
-				thread_data = (double)(x * Z() + z) / (lines * Z()) * 100;
+				if(PROGRESS) progress = (double)(x * Z() + z+1) / (lines * Z()) * 100;
 			}
 			//write slice data into target file
@@ -1061,8 +1059,6 @@ public:
 		//free the temporary frame
 		free(temp);
-		thread_data = 100;
-
 		return true;
 	}
@@ -24,11 +24,11 @@ protected:
 	std::fstream file;		//file stream used for reading and writing
 	std::string name;		//file name
-	unsigned long long int R[D];		//resolution
+	unsigned long long R[D];		//resolution
 	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
+	double progress;		//stores the progress on the current operation (accessible using a thread)
 	/// Private initialization function used to set default parameters in the data structure.
@@ -37,7 +37,24 @@ protected:
 		header = 0;									//initialize the header size to zero
 		mask = NULL;
-		thread_data = 0;
+		progress = 0;
+	}
+
+	//calculate the number of non-zero pixels in a mask
+	unsigned long long nnz(unsigned char* mask, unsigned long long N){
+
+		unsigned long long n = 0;						//initialize the counter to 0 (zero)
+		if(mask == NULL) return N;						//if the mask is NULL, assume all pixels are masked
+
+		for(unsigned long long i = 0; i < N; i++){		//for each pixel
+			if(mask[i] != 0) n++;						//increment the counter for every non-zero pixel in the mask
+		}
+		return n;										//return the number of nonzero pixels
+	}
+
+	/// Calculate the number of nonzero pixels in a mask over X-Y
+	unsigned long long nnz(unsigned char* mask){
+		return nnz(mask, R[0] * R[1]);
 	}
 	/// Private helper function that returns the size of the file on disk using system functions.
@@ -105,12 +122,12 @@ protected:
 public:
-	unsigned int get_thread_data(){
-		return thread_data;
+	unsigned int get_progress(){
+		return progress;
 	}
-	void reset_thread_data(){
-		thread_data = 0;
+	void reset_progress(){
+		progress = 0;
 	}
 	/// Open a binary file for streaming.
@@ -178,7 +195,9 @@ public:
 	/// @param p is a pointer to pre-allocated memory equal to the page size
 	/// @param page is the index of the page
-	bool read_page( T * p, unsigned int page){
+	bool read_page( T * p, unsigned int page, bool PROGRESS = false){
+
+		if(PROGRESS) progress = 0;
 		if (page >= R[2]){										//make sure the bank number is right
 			std::cout<<"ERROR: page out of range"<<std::endl;
@@ -187,7 +206,7 @@ public:
 		file.seekg(R[1] * R[0] * page * sizeof(T) + header, std::ios::beg);   //write into memory from the binary file
 		file.read((char *)p, R[0] * R[1] * sizeof(T));
-
+		if(PROGRESS) progress = 100;
 		return true;
 	}
@@ -198,9 +217,11 @@ public:
 	/// @param p is a pointer to pre-allocated memory equal to the line size R[2]
 	/// @param x is the x coordinate
 	/// @param y is the y coordinate
-	bool read_line_01( T * p, unsigned int x, unsigned int y){
+	bool read_line_2( T * p, unsigned int x, unsigned int y, bool PROGRESS = false){
 		unsigned int i;
+		if(PROGRESS) progress = 0;
+
 		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;
@@ -211,7 +232,9 @@ public:
 		{
 			file.read((char *)(p + i), sizeof(T));
 			file.seekg((R[1] * R[0] - 1) * sizeof(T), std::ios::cur);    //go to the next band
+			if(PROGRESS) progress = (double)i / (double)R[2] * 100;
 		}
+		if(PROGRESS) progress = 100;
 		return true;
 	}
@@ -221,7 +244,7 @@ public:
 	/// @param p is a pointer to pre-allocated memory equal to the line size R[2]
 	/// @param x is the y coordinate
 	/// @param y is the z coordinate
-	bool read_line_12(T * p, unsigned int y, unsigned int z){
+	bool read_line_0(T * p, unsigned int y, unsigned int z, bool PROGRESS = false){
 		//test to make sure the specified value is within range
 		if( y >= R[1] || z >= R[2] ){
 			std::cout<<"ERROR: sample or line out of range"<<std::endl;
@@ -230,7 +253,7 @@ public:
 		file.seekg((z * R[0] * R[1] + y * R[0]) * sizeof(T), std::ios::beg);	//seek to the start of the line
 		file.read((char *)p, sizeof(T) * R[0]);									//read the line
-
+		if(PROGRESS) progress = 100;
 		return true;
 	}
@@ -239,7 +262,8 @@ public:
 	/// @param p is a pointer to pre-allocated memory equal to the line size R[2]
 	/// @param x is the y coordinate
 	/// @param z is the z coordinate
-	bool read_line_02(T * p, unsigned int x, unsigned int z){
+	bool read_line_1(T * p, unsigned int x, unsigned int z, bool PROGRESS = false){
+		if(PROGRESS) progress = 0;
 		//test to make sure the specified value is within range
 		if( x >= R[0] || z >= R[2] ){
 			std::cout<<"ERROR: sample or line out of range"<<std::endl;
@@ -250,8 +274,9 @@ public:
 		for (unsigned int i = 0; i < R[1]; i++){									//for each pixel in the line
 			file.read((char *)(p + i), sizeof(T));					//read the pixel
 			file.seekg((R[0] - 1) * sizeof(T), std::ios::cur);		//seek to the next pixel in the line
+			if(PROGRESS) progress = (double)i / (double)R[1] * 100;
 		}
-
+		if(PROGRESS) progress = 100;
 		return true;
 	}
@@ -259,8 +284,8 @@ public:
 	/// @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){
-
+	bool read_plane_0(T* p, unsigned int n, bool PROGRESS = false){
+		if(PROGRESS) progress = 0;
 		if (n >= R[0]){										//make sure the number is within the possible range
 			std::cout<<"ERROR read_plane_0: page out of range"<<std::endl;
 			return false;
@@ -274,10 +299,10 @@ public:
 		for(unsigned int i = 0; i<N; i++){
 			file.read((char*)(p+i), sizeof(T));
 			file.seekg(jump, std::ios::cur);
-			thread_data = (double)i / N * 100;
+			if(PROGRESS) progress = (double)i / N * 100;
 		}
-		thread_data = 100;
+		if(PROGRESS) progress = 100;
 		return true;
@@ -287,8 +312,8 @@ public:
 	/// @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){
-
+	bool read_plane_1(T* p, unsigned int n, bool PROGRESS = false){
+		if(PROGRESS) progress = 0;
 		unsigned int L = R[0] * sizeof(T);		//caculate the number of bytes in a sample line
 		unsigned int jump = R[0] * (R[1] - 1) * sizeof(T);
@@ -299,13 +324,13 @@ public:
 		file.seekg(R[0] * n * sizeof(T), std::ios::beg);
 		for (unsigned i = 0; i < R[2]; i++){
-			thread_data = (double)i / R[2] * 100;
+			if(PROGRESS) progress = (double)i / R[2] * 100;
 			file.read((char *)(p + i * R[0]), L);
 			file.seekg( jump, std::ios::cur);
 			std::cout<<i<<"    ";
 		}
-		thread_data = 100;
+		if(PROGRESS) progress = 100;
 		return true;
 	}
@@ -313,8 +338,8 @@ public:
 	/// @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);
+	bool read_plane_2(T* p, unsigned int n, bool PROGRESS = false){
+		return read_page(p, n, PROGRESS);
 	}
 	/// Reads a single pixel, treating the entire data set as a linear array
@@ -41,14 +41,77 @@ protected:
 		return R[0];
 	}
-	using binary<T>::thread_data;
+	using binary<T>::progress;
+
+	/// Call the binary nnz() function for the BIP orientation
+	unsigned long long nnz(unsigned char* mask){
+		return binary<T>::nnz(mask, X()*Y());
+	}
+
+	/// Linear interpolation of a spectral value given the bounding spectral samples
+	T lerp(double w, T low_v, double low_w, T high_v, double high_w){
+		if(low_w == high_w) return low_v;										//if the interval is of zero length, just return one of the bounds
+		double alpha = (w - low_w) / (high_w - low_w);							//calculate the interpolation factor
+		return (1.0 - alpha) * low_v + alpha * high_v;							//interpolate
+	}
+
+	/// Gets the two band indices surrounding a given wavelength
+	void band_bounds(double wavelength, unsigned long long& low, unsigned long long& high){
+		unsigned long long B = Z();
+		for(high = 0; high < B; high++){
+			if(w[high] > wavelength) break;
+		}
+		low = 0;
+		if(high > 0)
+			low = high-1;
+	}
+
+	/// Get the list of band numbers that bound a list of wavelengths
+	void band_bounds(std::vector<double> wavelengths, 
+					 std::vector<unsigned long long>& low_bands, 
+					 std::vector<unsigned long long>& high_bands){
+
+		unsigned long long W = w.size();							//get the number of wavelengths in the list
+		low_bands.resize(W);											//pre-allocate space for the band lists
+		high_bands.resize(W);
+
+		for(unsigned long long wl = 0; wl < W; wl++){						//for each wavelength
+			band_bounds(wavelengths[wl], low_bands[wl], high_bands[wl]);	//find the low and high bands
+		}
+	}
+
+	/// Returns the interpolated in the given spectrum based on the given wavelength
+
+	/// @param s is the spectrum in main memory of length Z()
+	/// @param wavelength is the wavelength value to interpolate out
+	T interp_spectrum(T* s, double wavelength){
+		unsigned long long low, high;								//indices for the bands surrounding wavelength
+		band_bounds(wavelength, low, high);							//get the surrounding band indices
+
+		if(high == w.size()) return s[w.size()-1];					//if the high band is above the wavelength range, return the highest wavelength value
+		
+		return lerp(wavelength, s[low], w[low], s[high], w[high]);
+	}
+
+	/// Returns the interpolated value corresponding to the given list of wavelengths
+	std::vector<T> interp_spectrum(T* s, std::vector<double> wavelengths){
+
+		unsigned long long N = wavelengths.size();						//get the number of wavelength measurements
+
+		std::vector<T> v;										//allocate space for the resulting values
+		v.resize(wavelengths.size());
+		for(unsigned long long n = 0; n < N; n++){						//for each measurement
+			v[n] = interp_spectrum(s, wavelengths[n]);			//interpolate the measurement
+		}
+		return v;
+	}
 public:
 	using binary<T>::open;
 	using binary<T>::file;
 	using binary<T>::R;
-	using binary<T>::read_line_12;
+	using binary<T>::read_line_0;
 	/// Open a data file for reading using the class interface.
@@ -73,19 +136,19 @@ public:
 	/// @param p is a pointer to an allocated region of memory at least X * Y * sizeof(T) in size.
 	/// @param page <= B is the integer number of the band to be copied.
-	bool band_index( T * p, unsigned int page){
-		return binary<T>::read_plane_0(p, page);
+	bool band_index( T * p, unsigned int page, bool PROGRESS = false){
+		return binary<T>::read_plane_0(p, page, PROGRESS);
 	}
 	/// Retrieve a single band (by numerical label) and stores it in pre-allocated memory.
 	/// @param p is a pointer to an allocated region of memory at least X * Y * sizeof(T) in size.
 	/// @param wavelength is a floating point value (usually a wavelength in spectral data) used as a label for the band to be copied.
-	bool band( T * p, double wavelength){
+	bool band( T * p, double wavelength, bool PROGRESS = false){
 		//if there are no wavelengths in the BSQ file
 		if(w.size() == 0)
-			return band_index(p, (unsigned int)wavelength);
+			return band_index(p, (unsigned int)wavelength, PROGRESS);
 		unsigned int XY = X() * Y();	//calculate the number of pixels in a band
@@ -96,7 +159,7 @@ public:
 		//if wavelength is smaller than the first one in header file
 		if ( w[page] > wavelength ){
-			band_index(p, page);
+			band_index(p, page, PROGRESS);
 			return true;
 		}
@@ -105,7 +168,7 @@ public:
 			page++;
 			//if wavelength is larger than the last wavelength in header file
 			if (page == Z()) {
-				band_index(p, Z()-1);
+				band_index(p, Z()-1, PROGRESS);
 				return true;
 			}
 		}
@@ -116,7 +179,7 @@ public:
 			p1=(T*)malloc( XY * sizeof(T));                     //memory allocation
 			p2=(T*)malloc( XY * sizeof(T));
 			band_index(p1, page - 1);
-			band_index(p2, page );
+			band_index(p2, page, PROGRESS);
 			for(unsigned i=0; i < XY; i++){
 				double r = (double) (wavelength - w[page-1]) / (double) (w[page] - w[page-1]);
 				p[i] = (p2[i] - p1[i]) * r + p1[i];
@@ -126,7 +189,7 @@ public:
 		}
 		else                           //if the wavelength is equal to a wavelength in header file
 		{
-			band_index(p, page);
+			band_index(p, page, PROGRESS);
 		}
 		return true;
@@ -137,8 +200,8 @@ public:
 	/// @param p is a pointer to pre-allocated memory at least B * sizeof(T) in size.
 	/// @param x is the x-coordinate (dimension 1) of the spectrum.
 	/// @param y is the y-coordinate (dimension 2) of the spectrum.
-	bool spectrum(T * p, unsigned x, unsigned y){
-		return read_line_12(p, x, y);				//read a line in the binary YZ plane (dimension order for BIP is ZXY)
+	bool spectrum(T * p, unsigned x, unsigned y, bool PROGRESS = false){
+		return read_line_0(p, x, y, PROGRESS);				//read a line in the binary YZ plane (dimension order for BIP is ZXY)
 	}
 	/// Retrieves a band of x values from a given xz plane.
@@ -235,8 +298,8 @@ public:
 	/// @param n is an integer index to the pixel using linear array indexing.
 	bool pixel(T * p, unsigned n){
-		unsigned bandnum = X() * Y();		//calculate numbers in one band
-		if ( n >= bandnum){							//make sure the pixel number is right
+		unsigned N = X() * Y();					//calculate numbers in one band
+		if ( n >= N){							//make sure the pixel number is right
 			std::cout<<"ERROR: sample or line out of range"<<std::endl;
 			return false;
 		}
@@ -255,113 +318,57 @@ public:
 	/// @param outname is the name of the output file used to store the resulting baseline-corrected data.
 	/// @param wls is the list of baseline points based on band labels.	
-	bool baseline(std::string outname, std::vector<double> wls){
-	
-		unsigned N = wls.size();			//get the number of baseline points
+	bool baseline(std::string outname, std::vector<double> base_pts, bool PROGRESS = false){
 		std::ofstream target(outname.c_str(), std::ios::binary);	//open the target binary file
 		std::string headername = outname + ".hdr";              //the header file name		
-		//simplify image resolution
-		unsigned int ZX = Z() * X();		//calculate the number of points in a Y slice
-		unsigned int L = ZX * sizeof(T);			//calculate the number of bytes of a Y slice
-		unsigned int B = Z();
-		
-		T* c;			//pointer to the current Y slice
-		c = (T*)malloc(L);  //memory allocation
-		
-		T* a;			//pointer to the two YZ lines surrounding the current YZ line
-		T* b;
-		
-		a = (T*)malloc(X() * sizeof(T));
-		b = (T*)malloc(X() * sizeof(T));
-
-
-		double ai, bi;	//stores the two baseline points wavelength surrounding the current band
-		double ci;		//stores the current band's wavelength
-		unsigned control;
-
-		if (a == NULL || b == NULL || c == NULL){
-			std::cout<<"ERROR: error allocating memory";
-			exit(1);
-		}
-	//	loop start	correct every y slice
-		
-		for (unsigned k =0; k < Y(); k++)
-		{
-			//get the current y slice
-			read_plane_y(c, k);
-		
-			//initialize lownum, highnum, low, high		
-			control=0;
-			ai = w[0];
-			//if no baseline point is specified at band 0,
-			//set the baseline point at band 0 to 0
-			if(wls[0] != w[0]){
-				bi = wls[control];			
-				memset(a, (char)0, X() * sizeof(T) );
-			}
-			//else get the low band
-			else{
-				control++;
-				read_x_from_xz(a, c, ai);
-				bi = wls[control];
-			}
-			//get the high band
-			read_x_from_xz(b, c, bi);
-		
-			//correct every YZ line
-			
-			for(unsigned cii = 0; cii < B; cii++){
-				//update baseline points, if necessary
-				if( w[cii] >= bi && cii != B - 1) {
-					//if the high band is now on the last BL point?
-					if (control != N-1) {
-	
-						control++;		//increment the index
-	
-						std::swap(a, b);	//swap the baseline band pointers
-	
-						ai = bi;
-						bi = wls[control];
-						read_x_from_xz(b, c, bi);
-	
-					}
-					//if the last BL point on the last band of the file?
-					else if ( wls[control] < w[B - 1]) {
-	
-						std::swap(a, b);	//swap the baseline band pointers
-	
-						memset(b, (char)0, X() * sizeof(T) );	//clear the high band
-	
-						ai = bi;
-						bi = w[B - 1];
-					}
+		unsigned long long N = X() * Y();						//calculate the total number of pixels to be processed
+		unsigned long long B = Z();								//get the number of bands
+		T* s = (T*)malloc(sizeof(T) * B);						//allocate memory to store a pixel
+		T* sbc = (T*)malloc(sizeof(T) * B);						//allocate memory to store the baseline corrected spectrum
+
+		std::vector<T> base_vals;								//allocate space for the values at each baseline point
+		double aw, bw;											//surrounding baseline point wavelengths
+		T av, bv;												//surrounding baseline point values
+		double alpha;
+		unsigned long long ai, bi;								//surrounding baseline point band indices
+		for(unsigned long long n = 0; n < N; n++){				//for each pixel in the image
+			pixel(s, n);										//retrieve the spectrum s
+			base_vals = interp_spectrum(s, base_pts);			//get the values at each baseline point
+
+			ai = bi = 0;
+			aw = w[0];											//initialize the current baseline points (assume the spectrum starts at 0)
+			av = s[0];
+			bw = base_pts[0];
+			for(unsigned long long b = 0; b < B; b++){			//for each band in the spectrum
+				while(bi < base_pts.size() && base_pts[bi] < w[b])	//while the current wavelength is higher than the second baseline point
+					bi++;											//move to the next baseline point
+				if(bi < base_pts.size()){
+					bw = base_pts[bi];							//set the wavelength for the upper bound baseline point
+					bv = base_vals[bi];							//set the value for the upper bound baseline point
 				}
-
-				ci = w[cii];
-			
-				//perform the baseline correction
-				for(unsigned i=0; i < X(); i++)
-				{
-					double r = (double) (ci - ai) / (double) (bi - ai);
-					c[i * B + cii] =(T) ( c[i * B + cii] - (b[i] - a[i]) * r - a[i] );
+				if(bi == base_pts.size()){						//if we have passed the last baseline point
+					bw = w[B-1];								//set the outer bound to the last spectral band
+					bv = s[B-1];
 				}
-				
-			}//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
+				if(bi != 0){
+					ai = bi - 1;								//set the lower bound baseline point index
+					aw = base_pts[ai];							//set the wavelength for the lower bound baseline point
+					av = base_vals[ai];							//set the value for the lower bound baseline point
+				}
+				sbc[b] = s[b] - lerp(w[b], av, aw, bv, bw);		//perform the baseline correction and save the new value
+			}
+			if(PROGRESS) progress = (double)(n+1) / N * 100;	//set the current progress
+			target.write((char*)sbc, sizeof(T) * B);	//write the corrected data into destination
+		}														//end for each pixel
-		free(a);
-		free(b);
-		free(c);
+		free(s);
+		free(sbc);
 		target.close();
-
-		thread_data = 100;
+		
 		return true;	
 	}
@@ -371,80 +378,41 @@ public:
 	/// @param outname is the name of the output file used to store the resulting baseline-corrected data.
 	///	@param w is the label specifying the band that the hyperspectral image will be normalized to.
 	///	@param t is a threshold specified such that a spectrum with a value at w less than t is set to zero. Setting this threshold allows the user to limit division by extremely small numbers.
-	bool normalize(std::string outname, double w, double t = 0.0)
+	bool normalize(std::string outname, double w, double t = 0.0, bool PROGRESS = false)
 	{
-		unsigned int B = Z();		//calculate the number of bands
-		unsigned int ZX = Z() * X();
-		unsigned int XY = X() * Y();	//calculate the number of pixels in a band
-		unsigned int S = XY * sizeof(T);		//calculate the number of bytes in a band		
-		unsigned int L = ZX * sizeof(T);
-
 		std::ofstream target(outname.c_str(), std::ios::binary);	//open the target binary file
-		std::string headername = outname + ".hdr";              //the header file name
-
-		T * c;				//pointer to the current ZX slice
-		T * b;				//pointer to the standard band
-
-		b = (T*)malloc( S );     //memory allocation
-		c = (T*)malloc( L ); 
-
-		band(b, w);             //get the certain band into memory
-
-		for(unsigned j = 0; j < Y(); j++)
-		{
-			read_plane_y(c, j);
-			unsigned jX = j * X();		//to avoid calculating it many times
-			for(unsigned i = 0; i < X(); i++)
-			{
-				unsigned iB = i * B;		
-				for(unsigned m = 0; m < B; m++)
-				{
-					if( b[i+jX] < t )
-						c[m + iB] = (T)0.0;
-					else
-						c[m + iB] = c[m + iB] / b[i + jX];			//perform normalization
-				}								
+		std::string headername = outname + ".hdr";              //the header file name		
+		
+		unsigned long long N = X() * Y();						//calculate the total number of pixels to be processed
+		unsigned long long B = Z();								//get the number of bands
+		T* s = (T*)malloc(sizeof(T) * B);						//allocate memory to store a pixel
+		T nv;													//stores the value of the normalized band
+		for(unsigned long long n = 0; n < N; n++){				//for each pixel in the image
+			pixel(s, n);										//retrieve the spectrum s
+			nv = interp_spectrum(s, w);							//find the value of the normalization band
+			if(abs(nv) <= t)									//if the normalization band is below threshold
+				memset(s, 0, sizeof(T) * B);					//set the output to zero
+			else{
+				for(unsigned long long b = 0; b < B; b++){			//for each band in the spectrum
+					s[b] /= nv;										//divide by the normalization value
+				}
 			}
-			target.write(reinterpret_cast<const char*>(c), L);   //write normalized data into destination
-
-			thread_data = (double) j / Y() * 100;
-		}
+		
+			if(PROGRESS) progress = (double)(n+1) / N * 100;	//set the current progress
+			target.write((char*)s, sizeof(T) * B);	//write the corrected data into destination
+		}														//end for each pixel
-		free(b);
-		free(c);
+		free(s);
 		target.close();
-		thread_data = 100;
 		return true;
 	}
-	/// Convert the current BIP file to a BSQ file with the specified file name.
-
-	/// @param outname is the name of the output BSQ file to be saved to disk.
-	bool bsq(std::string outname)
-	{	
-		std::string temp = outname + "_temp";
-		std::string headtemp = temp + ".hdr";
-		//first creat a temporary bil file and convert bip file to bil file
-		bil(temp);
-
-		stim::bil<T> n;
-		if(n.open(temp, X(), Y(), Z(), offset, w)==false){        //open infile
-			std::cout<<"ERROR: unable to open input file"<<std::endl;
-			exit(1);
-		}
-		//then convert bil file to bsq file
-		n.bsq(outname);
-		n.close();
-		remove(temp.c_str());
-		remove(headtemp.c_str());
-		return true;
-	}
 	/// Convert the current BIP file to a BIL file with the specified file name.
 	/// @param outname is the name of the output BIL file to be saved to disk.
-	bool bil(std::string outname)
+	bool bil(std::string outname, bool PROGRESS = false)
 	{
 		unsigned int S = X() * Z() * sizeof(T);		//calculate the number of bytes in a ZX slice
@@ -465,13 +433,11 @@ public:
 				for ( unsigned j = 0; j < X(); j++)
 					q[ks + j] = p[k + j * Z()];
-				thread_data = (double)(i * Z() + k) / (Y() * Z()) * 100;
+				if(PROGRESS) progress = (double)(i * Z() + k+1) / (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);
 		target.close();
@@ -909,7 +875,7 @@ public:
 				//write this pixel out
 				target.write((char *)spectrum, B * sizeof(T));
-				thread_data = (double) x / XY * 100;
+				progress = (double) x / XY * 100;
 			}
 		}
@@ -918,7 +884,7 @@ public:
 		target.close();
 		free(spectrum);
-		thread_data = 100;
+		progress = 100;
 		return true;
 	}
@@ -962,7 +928,7 @@ public:
 				target.write((char *)spectrum, B * sizeof(T));
 			}
-			thread_data = (double)x / XY * 100;
+			progress = (double)x / XY * 100;
 		}
@@ -970,7 +936,7 @@ public:
 		target.close();
 		free(spectrum);
-		thread_data = 100;
+		progress = 100;
 		return true;
@@ -981,7 +947,7 @@ public:
 	/// @param p is a pointer to pre-allocated memory of size [B * sizeof(T)] that stores the mean spectrum
 	/// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
-	bool avg_band(double* p, unsigned char* mask = NULL){
+	bool avg_band(double* p, unsigned char* mask = NULL, bool PROGRESS = false){
 		unsigned long long XY = X() * Y();
 		T* temp = (T*)malloc(sizeof(T) * Z());
 		//Iinitialize
@@ -989,12 +955,8 @@ public:
 			p[j] = 0;
 		}
 		//calculate vaild number in a band
-		unsigned count = 0;
-		for (unsigned j = 0; j < XY; j++){
-			if (mask == NULL || mask[j] != 0){
-				count++;
-			}
-		}
+		unsigned count = nnz(mask);
+
 		//calculate average number of a band
 		for (unsigned i = 0; i < XY; i++){
 			if (mask == NULL || mask[i] != 0){			
@@ -1003,22 +965,21 @@ public:
 					p[j] += (double)temp[j] / (double)count;
 				}
 			}
-			thread_data = (double)i / XY * 100;
+			if(PROGRESS) progress = (double)(i+1) / XY * 100;
 		}
 		free(temp);
-		thread_data = 100;
 		return true;
 	}
 #ifdef CUDA_FOUND
 	/// Calculate the covariance matrix for masked pixels using cuBLAS
-	bool cu_co_matrix(double* co, double* avg, unsigned char *mask){
+	bool co_matrix_cublas(double* co, double* avg, unsigned char *mask, bool PROGRESS = false){
 		cudaError_t cudaStat;    
 		cublasStatus_t stat;
 		cublasHandle_t handle;
-		thread_data = 0;													//initialize the progress to zero (0)
+		progress = 0;													//initialize the progress to zero (0)
 		unsigned long long XY = X() * Y();									//calculate the number of elements in a band image
 		unsigned long long B = Z();											//calculate the number of spectral elements
@@ -1047,7 +1008,7 @@ public:
 				stat = cublasDaxpy(handle, B, &axpy_alpha, avg_dev, 1, s_dev, 1);				//subtract the average spectrum
 				stat = cublasDsyr(handle, CUBLAS_FILL_MODE_UPPER, B, &ger_alpha, s_dev, 1, A_dev, B);	//calculate the covariance matrix (symmetric outer product)
 			}
-			thread_data = (double)xy / XY * 100;													//record the current progress
+			if(PROGRESS) progress = (double)(xy+1) / XY * 100;													//record the current progress
 		}
@@ -1063,7 +1024,6 @@ public:
 			}
 		}
-		thread_data = 100;														//processing complete
 		return true;
 	}
 #endif
@@ -1073,23 +1033,19 @@ public:
 	/// @param co is a pointer to pre-allocated memory of size [B * B] that stores the resulting covariance matrix
 	/// @param avg is a pointer to memory of size B that stores the average spectrum
 	/// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
-	bool co_matrix(double* co, double* avg, unsigned char *mask){
+	bool co_matrix(double* co, double* avg, unsigned char *mask, bool PROGRESS = false){
 #ifdef CUDA_FOUND
-		return cu_co_matrix(co, avg, mask);								//if CUDA is available, use cuBLAS to calculate the covariance matrix
+		return co_matrix_cublas(co, avg, mask, PROGRESS);								//if CUDA is available, use cuBLAS to calculate the covariance matrix
 #endif
-		thread_data = 0;
+		progress = 0;
 		//memory allocation
 		unsigned long long XY = X() * Y();
 		unsigned long long B = Z();
 		T* temp = (T*)malloc(sizeof(T) * B);
 		//count vaild pixels in a band
-		unsigned long long count = 0;
-		for (unsigned long long j = 0; j < XY; j++){
-			if (mask == NULL || mask[j] != 0){
-				count++;
-			}
-		}
+		unsigned long long count = nnz(mask);
+
 		//initialize covariance matrix
 		memset(co, 0, B * B * sizeof(double));
@@ -1110,7 +1066,7 @@ public:
 						co_half[idx++] += temp_precise[b0] * temp_precise[b1];
 				}
 			}
-			thread_data = (double)xy / XY * 100;
+			if(PROGRESS) progress = (double)(xy+1) / XY * 100;
 		}
 		idx = 0;
 		for (unsigned long long i = 0; i < B; i++){										//copy the precision matrix to both halves of the output matrix
@@ -1119,12 +1075,78 @@ public:
 			}
 		}
-		thread_data = 100;							//processing complete
 		free(temp);
 		free(temp_precise);
 		return true;
 	}
+	bool project(std::string outfile, double* center, double* basis, unsigned long long M, bool PROGRESS = false){
+
+		std::ofstream target(outfile.c_str(), std::ios::binary);	//open the target binary file
+		std::string headername = outfile + ".hdr";					//the header file name
+
+		//memory allocation
+		unsigned long long XY = X() * Y();
+		unsigned long long B = Z();
+
+		T* s = (T*)malloc(sizeof(T) * B);							//allocate space for the spectrum
+		T* rs = (T*)malloc(sizeof(T) * M);							//allocate space for the projected spectrum
+		double* bv;													//pointer to the current projection vector
+		for(unsigned long long xy = 0; xy < XY; xy++){				//for each spectrum in the image
+			pixel(s, xy);											//load the spectrum
+			memset(rs, 0, sizeof(T) * M);							//initialize the rotated spectrum to zero (0)
+			for(unsigned long long m = 0; m < M; m++){				//for each basis vector
+				bv = &basis[m * B];									//assign 'bv' to the beginning of the basis vector
+				for(unsigned long long b = 0; b < B; b++){			//for each band
+					rs[m] += (s[b] - center[b]) * bv[b];				//center the spectrum and perform the projection
+				}
+			}
+
+			target.write(reinterpret_cast<const char*>(rs), sizeof(T) * M);					//write the projected vector
+			if(PROGRESS) progress = (double)(xy+1) / XY * 100;
+		}
+
+		free(s);													//free temporary storage arrays
+		free(rs);
+		target.close();												//close the output file
+
+		return true;
+	}
+
+	bool inverse(std::string outfile, double* center, double* basis, unsigned long long B, unsigned long long C = 0, bool PROGRESS = false){
+
+		std::ofstream target(outfile.c_str(), std::ios::binary);	//open the target binary file
+		std::string headername = outfile + ".hdr";					//the header file name
+
+		//memory allocation
+		unsigned long long XY = X() * Y();
+		if(C == 0) C = Z();											//if no coefficient number is given, assume all are used
+		C = std::min<unsigned long long>(C, Z());					//set the number of coefficients (the user can specify fewer)
+
+		T* coeff = (T*)malloc(sizeof(T) * Z());						//allocate space for the coefficients
+		T* s = (T*)malloc(sizeof(T) * B);							//allocate space for the spectrum
+		double* bv;													//pointer to the current projection vector
+		for(unsigned long long xy = 0; xy < XY; xy++){				//for each pixel in the image (expressed as a set of coefficients)
+			pixel(coeff, xy);										//load the coefficients
+			memset(s, 0, sizeof(T) * B);							//initialize the spectrum to zero (0)
+			for(unsigned long long c = 0; c < C; c++){				//for each basis vector coefficient
+				bv = &basis[c * B];									//assign 'bv' to the beginning of the basis vector
+				for(unsigned long long b = 0; b < B; b++){			//for each component of the basis vector					
+					s[b] += coeff[c] * bv[b] + center[b];			//calculate the contribution of each element of the basis vector in the final spectrum
+				}
+			}
+
+			target.write(reinterpret_cast<const char*>(s), sizeof(T) * B);					//write the projected vector
+			if(PROGRESS) progress = (double)(xy+1) / XY * 100;
+		}
+
+		free(coeff);												//free temporary storage arrays
+		free(s);
+		target.close();												//close the output file
+
+		return true;
+	}
+
 	/// Crop a region of the image and save it to a new file.
@@ -1138,7 +1160,8 @@ public:
 								   unsigned long long x1,
 								   unsigned long long y1,
 								   unsigned long long b0,
-								   unsigned long long b1){
+								   unsigned long long b1,
+								   bool PROGRESS = false){
 		//calculate the new number of samples, lines, and bands
 		unsigned long long samples = x1 - x0;
@@ -1180,14 +1203,13 @@ public:
 				file.seekg(jump_sample, std::ios::cur);
-				thread_data = (double)(y * samples + x) / (lines * samples) * 100;
+				if(PROGRESS) progress = (double)((y+1) * samples + x + 1) / (lines * samples) * 100;
 			}
 			file.seekg(jump_line, std::ios::cur);
 		}
 		free(temp);
-		thread_data = 100;
 		return true;
 	}
@@ -42,13 +42,14 @@ protected:
 		return R[2];
 	}
-	using binary<T>::thread_data;
+	using binary<T>::nnz;
+	using binary<T>::progress;
 public:
 	using binary<T>::open;
 	using binary<T>::file;
-	using binary<T>::read_line_01;
+	using binary<T>::read_line_2;
 	using binary<T>::read_plane_2;
 	//using binary<T>::getSlice;
@@ -76,15 +77,15 @@ public:
 	/// @param p is a pointer to an allocated region of memory at least X * Y * sizeof(T) in size.
 	/// @param page <= B is the integer number of the band to be copied.
 	bool band_index( T * p, unsigned int page){
-		return read_plane_2(p, page);
+		return read_plane_2(p, page);		//call the binary read_plane function (don't let it update the progress)
 	}
 	/// Retrieve a single band (by numerical label) and stores it in pre-allocated memory.
 	/// @param p is a pointer to an allocated region of memory at least X * Y * sizeof(T) in size.
 	/// @param wavelength is a floating point value (usually a wavelength in spectral data) used as a label for the band to be copied.
-	bool band( T * p, double wavelength){
-
+	bool band( T * p, double wavelength, bool PROGRESS = false){
+		if(PROGRESS) progress = 0;
 		//if there are no wavelengths in the BSQ file
 		if(w.size() == 0)
 			return band_index(p, (unsigned int)wavelength);
@@ -132,7 +133,7 @@ public:
 		else{
 			band_index(p, page);		//return the band
 		}
-
+		if(PROGRESS) progress = 100;
 		return true;
 	}
@@ -141,8 +142,8 @@ public:
 	/// @param p is a pointer to pre-allocated memory at least B * sizeof(T) in size.
 	/// @param x is the x-coordinate (dimension 1) of the spectrum.
 	/// @param y is the y-coordinate (dimension 2) of the spectrum.
-	bool spectrum(T * p, unsigned x, unsigned y){
-		return read_line_01(p, x, y);
+	bool spectrum(T * p, unsigned x, unsigned y, bool PROGRESS = false){
+		return read_line_2(p, x, y, PROGRESS);
 	}
 	/// Retrieve a single pixel and stores it in pre-allocated memory.
@@ -171,8 +172,9 @@ public:
 	/// @param outname is the name of the output file used to store the resulting baseline-corrected data.
 	/// @param wls is the list of baseline points based on band labels.
-	bool baseline(std::string outname, std::vector<double> wls )
+	bool baseline(std::string outname, std::vector<double> wls, bool PROGRESS = false )
 	{
+		if(PROGRESS) progress = 0;
 		unsigned N = wls.size();			//get the number of baseline points
 		std::ofstream target(outname.c_str(), std::ios::binary);	//open the target binary file
@@ -261,7 +263,7 @@ public:
 			target.write(reinterpret_cast<const char*>(c), S);   //write the corrected data into destination
-			thread_data = (double)cii / B * 100;
+			if(PROGRESS)progress = (double)(cii+1) / B * 100;
 		}
@@ -272,7 +274,6 @@ public:
 		free(c);
 		target.close();
-		thread_data = 100;
 		return true;
 	}
@@ -281,7 +282,7 @@ public:
 	/// @param outname is the name of the output file used to store the resulting baseline-corrected data.
 	///	@param w is the label specifying the band that the hyperspectral image will be normalized to.
 	///	@param t is a threshold specified such that a spectrum with a value at w less than t is set to zero. Setting this threshold allows the user to limit division by extremely small numbers.
-	bool normalize(std::string outname, double w, double t = 0.0)
+	bool normalize(std::string outname, double w, double t = 0.0, bool PROGRESS = false)
 	{
 		unsigned int B = Z();		//calculate the number of bands
 		unsigned int XY = X() * Y();	//calculate the number of pixels in a band
@@ -310,7 +311,7 @@ public:
 			}
 			target.write(reinterpret_cast<const char*>(c), S);   //write normalized data into destination
-			thread_data = (double)j / B * 100;
+			if(PROGRESS) progress = (double)(j+1) / B * 100;
 		}
@@ -320,14 +321,13 @@ public:
 		free(b);
 		free(c);
 		target.close();
-		thread_data = 100;					//make sure that the progress bar is full
 		return true;
 	}
 	/// Convert the current BSQ file to a BIP file with the specified file name.
 	/// @param outname is the name of the output BIP file to be saved to disk.
-	bool bip(std::string outname)
+	/*bool bip(std::string outname)
 	{
 		std::string temp = outname + "_temp";
 		std::string headtemp = temp + ".hdr";
@@ -340,17 +340,18 @@ public:
 			exit(1);
 		}
 		//then convert bil file to bip file
+		progress = 0;
 		n.bip(outname);
 		n.close();
 		remove(temp.c_str());
 		remove(headtemp.c_str());
 		return true;
-	}
+	}*/
 	/// Convert the current BSQ file to a BIL file with the specified file name.
 	/// @param outname is the name of the output BIL file to be saved to disk.
-	bool bil(std::string outname)
+	bool bil(std::string outname, bool PROGRESS = false)
 	{
 		//simplify image resolution
 		unsigned int L = X() * Z() * sizeof(T);		//calculate the number of bytes of a ZX slice
@@ -370,17 +371,11 @@ public:
 			{
 				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
+				if(PROGRESS) progress = (double)(y * Z() + z + 1) / (Z() * Y()) * 100;	//update the progress counter
 			}
 			target.write(reinterpret_cast<const char*>(xz_slice), xz_bytes);   //write the generated XZ slice to the target file
 		}
-
-
-		thread_data = 100;
-
 		free(xz_slice);
 		target.close();
@@ -869,13 +864,13 @@ public:
 					continue;
 				}
-				thread_data = (double)(i * XY + j) / (XY * Z()) * 100;
+				progress = (double)(i * XY + j) / (XY * Z()) * 100;
 			}
 		}
 		target.close();
 		free(temp);
-		thread_data = 100;
+		progress = 100;
 		return true;
 	}
@@ -923,7 +918,7 @@ public:
 					i++;
 				}
 				//std::cout<<xi<<"/"<<XY<<",   "<<b<<"/"<<B<<std::endl;
-				thread_data = (double)(b * XY + xi) / (B * XY) * 100;
+				progress = (double)(b * XY + xi) / (B * XY) * 100;
 			}
 			//std::cout<<b*XY<<"/"<<B*XY<<"      "<<B<<"-"<<XY<<std::endl;
 			//write the band image to disk
@@ -931,7 +926,7 @@ public:
 		}
 		//std::cout<<"unsifted"<<std::endl;
-		thread_data = 100;
+		progress = 100;
 		return true;
 	}
@@ -963,7 +958,7 @@ public:
 	/// @param p is a pointer to pre-allocated memory of size [B * sizeof(T)] that stores the mean spectrum
 	/// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
-	bool avg_band(double* p, unsigned char* mask = NULL){
+	bool avg_band(double* p, unsigned char* mask = NULL, bool PROGRESS = false){
 		unsigned long long XY = X() * Y();
 		unsigned long long count = 0;						//count will store the number of masked pixels
 		T* temp = (T*)malloc(sizeof(T) * XY);
@@ -983,6 +978,7 @@ public:
 					p[i] += (double)temp[j] / (double)count;
 				}
 			}
+			if(PROGRESS) progress = (double)(i+1) / Z() * 100;
 		}
 		free(temp);
 		return true;
@@ -993,8 +989,8 @@ public:
 	/// @param co is a pointer to pre-allocated memory of size [B * B] that stores the resulting covariance matrix
 	/// @param avg is a pointer to memory of size B that stores the average spectrum
 	/// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
-	bool co_matrix(double* co, double* avg, unsigned char *mask){
-		thread_data = 0;							//initialize thread data for timers and progress bars
+	/*bool co_matrix(double* co, double* avg, unsigned char *mask){
+		progress = 0;							//initialize thread data for timers and progress bars
 		//memory allocation
 		unsigned long long xy = X() * Y();
 		unsigned long long B = Z();
@@ -1024,15 +1020,15 @@ public:
 				co[i*B + j] = numerator;
 				co[j*B + i] = numerator;		//because correlated matrix is a symmetric matrix
 			}
-			thread_data = (double)i / B * 100;
+			progress = (double)i / B * 100;
 		}
 		free(bandi);
 		free(bandj);
-		thread_data = 100;							//processing complete
+		progress = 100;							//processing complete
 		return true;
-	}
+	}*/
 	/// Crop a region of the image and save it to a new file.
@@ -1047,7 +1043,8 @@ public:
 								   unsigned long long x1,
 								   unsigned long long y1,
 								   unsigned long long b0,
-								   unsigned long long b1){
+								   unsigned long long b1,
+								   bool PROGRESS = false){
 		//calculate the new number of samples, lines, and bands
 		unsigned long long samples = x1 - x0;
@@ -1080,15 +1077,13 @@ public:
 				file.read((char *)(temp + y * samples), sizeof(T) * samples);
 				file.seekg(jumpl, std::ios::cur);    //go to the next band
-				thread_data = (double)(z * lines + y) / (Z() * lines) * 100;
+				if(PROGRESS) progress = (double)((z+1) * lines + y + 1) / ((b1 - b0) * lines) * 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;
 	}
@@ -37,27 +37,27 @@ public:
 		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();
+				((bsq<float>*)file)->reset_progress();
 			else if(header.data_type == envi_header::float64)
-				((bsq<double>*)file)->reset_thread_data();
+				((bsq<double>*)file)->reset_progress();
 			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();
+				((bil<float>*)file)->reset_progress();
 			else if(header.data_type == envi_header::float64)
-				((bil<double>*)file)->reset_thread_data();
+				((bil<double>*)file)->reset_progress();
 			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();
+				((bip<float>*)file)->reset_progress();
 			else if(header.data_type == envi_header::float64)
-				((bip<double>*)file)->reset_thread_data();
+				((bip<double>*)file)->reset_progress();
 			else
 				std::cout<<"ERROR: unidentified data type"<<std::endl;
 		}
@@ -73,27 +73,27 @@ public:
 		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();
+				return ((bsq<float>*)file)->get_progress();
 			else if(header.data_type == envi_header::float64)
-				return ((bsq<double>*)file)->get_thread_data();
+				return ((bsq<double>*)file)->get_progress();
 			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();
+				return ((bil<float>*)file)->get_progress();
 			else if(header.data_type == envi_header::float64)
-				return ((bil<double>*)file)->get_thread_data();
+				return ((bil<double>*)file)->get_progress();
 			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();
+				return ((bip<float>*)file)->get_progress();
 			else if(header.data_type == envi_header::float64)
-				return ((bip<double>*)file)->get_thread_data();
+				return ((bip<double>*)file)->get_progress();
 			else
 				std::cout<<"ERROR: unidentified data type"<<std::endl;
 		}
@@ -213,31 +213,31 @@ public:
 	/// @param outfile is the name of the normalized file to be output
 	/// @param band is the band label to be output
 	/// @param threshold is a threshold value specified such that normalization will only be done to values in the band > threshold (preventing division by small numbers)
-	bool normalize(std::string outfile, double band, double threshold = 0.0){
+	bool normalize(std::string outfile, double band, double threshold = 0.0, bool PROGRESS = false){
 		if(header.interleave == envi_header::BSQ){		//if the infile is bsq file
 			if(header.data_type ==envi_header::float32)
-				return ((bsq<float>*)file)->normalize(outfile, band, threshold);
+				return ((bsq<float>*)file)->normalize(outfile, band, threshold, PROGRESS);
 			else if(header.data_type == envi_header::float64)
-				return ((bsq<double>*)file)->normalize(outfile,band, threshold);
+				return ((bsq<double>*)file)->normalize(outfile,band, threshold, PROGRESS);
 			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)->normalize(outfile, band);
+				return ((bil<float>*)file)->normalize(outfile, band, threshold, PROGRESS);
 			else if(header.data_type == envi_header::float64)
-				return ((bil<double>*)file)->normalize(outfile,band);
+				return ((bil<double>*)file)->normalize(outfile,band, threshold, PROGRESS);
 			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)->normalize(outfile, band);
+				return ((bip<float>*)file)->normalize(outfile, band, threshold, PROGRESS);
 			else if(header.data_type == envi_header::float64)
-				return ((bip<double>*)file)->normalize(outfile,band);
+				return ((bip<double>*)file)->normalize(outfile,band, threshold, PROGRESS);
 			else
 				std::cout<<"ERROR: unidentified data type"<<std::endl;
 		}
@@ -253,13 +253,13 @@ public:
 	/// @param outfile is the file name for the baseline corrected output
 	/// @param w is a list of band labels to serve as baseline points (zero values)
-	bool baseline(std::string outfile, std::vector<double> w){
+	bool baseline(std::string outfile, std::vector<double> w, bool PROGRESS){
 		if(header.interleave == envi_header::BSQ){		//if the infile is bsq file
 			if(header.data_type ==envi_header::float32)
-				return ((bsq<float>*)file)->baseline(outfile, w);
+				return ((bsq<float>*)file)->baseline(outfile, w, PROGRESS);
 			else if(header.data_type == envi_header::float64)
-				return ((bsq<double>*)file)->baseline(outfile,w);
+				return ((bsq<double>*)file)->baseline(outfile,w, PROGRESS);
 			else{
 				std::cout<<"ERROR: unidentified data type"<<std::endl;
 				exit(1);
@@ -268,9 +268,9 @@ public:
 		else if(header.interleave == envi_header::BIL){		//if the infile is bil file
 			if(header.data_type ==envi_header::float32)
-				return ((bil<float>*)file)->baseline(outfile, w);
+				return ((bil<float>*)file)->baseline(outfile, w, PROGRESS);
 			else if(header.data_type == envi_header::float64)
-				return ((bil<double>*)file)->baseline(outfile, w);
+				return ((bil<double>*)file)->baseline(outfile, w, PROGRESS);
 			else{
 				std::cout<<"ERROR: unidentified data type"<<std::endl;
 				exit(1);
@@ -279,9 +279,9 @@ public:
 		else if(header.interleave == envi_header::BIP){		//if the infile is bip file
 			if(header.data_type ==envi_header::float32)
-				return ((bip<float>*)file)->baseline(outfile, w);
+				return ((bip<float>*)file)->baseline(outfile, w, PROGRESS);
 			else if(header.data_type == envi_header::float64)
-				return ((bip<double>*)file)->baseline(outfile, w);
+				return ((bip<double>*)file)->baseline(outfile, w, PROGRESS);
 			else{
 				std::cout<<"ERROR: unidentified data type"<<std::endl;
 				exit(1);
@@ -294,11 +294,65 @@ public:
 		}
 	}
+	void project(std::string outfile, double* center, double* basis, unsigned long long M, bool PROGRESS = false){
+		if(header.interleave == envi_header::BSQ){		//if the infile is bsq file
+			std::cout<<"ERROR: BSQ projection not supported"<<std::endl;
+			exit(1);
+		}
+
+		else if(header.interleave == envi_header::BIL){		//if the infile is bil file
+			std::cout<<"ERROR: BIL projection not supported"<<std::endl;
+			exit(1);
+		}
+
+		else if(header.interleave == envi_header::BIP){		//if the infile is bip file
+			if(header.data_type ==envi_header::float32)
+				((bip<float>*)file)->project(outfile, center, basis, M, PROGRESS);
+			else if(header.data_type == envi_header::float64)
+				((bip<double>*)file)->project(outfile, center, basis, M, PROGRESS);
+			else{
+				std::cout<<"ERROR: unidentified data type"<<std::endl;
+				exit(1);
+			}
+		}
+
+		stim::envi_header out_hdr = header;							
+		out_hdr.bands = M;											//set the number of bands in the output header
+		out_hdr.save(outfile + ".hdr");								//save the output header
+	}
+
+	void inverse(std::string outfile, double* center, double* basis, unsigned long long B, unsigned long long C = 0, bool PROGRESS = false){
+		if(header.interleave == envi_header::BSQ){		//if the infile is bsq file
+			std::cout<<"ERROR: BSQ projection not supported"<<std::endl;
+			exit(1);
+		}
+
+		else if(header.interleave == envi_header::BIL){		//if the infile is bil file
+			std::cout<<"ERROR: BIL projection not supported"<<std::endl;
+			exit(1);
+		}
+
+		else if(header.interleave == envi_header::BIP){		//if the infile is bip file
+			if(header.data_type ==envi_header::float32)
+				((bip<float>*)file)->inverse(outfile, center, basis, B, C, PROGRESS);
+			else if(header.data_type == envi_header::float64)
+				((bip<double>*)file)->inverse(outfile, center, basis, B, C, PROGRESS);
+			else{
+				std::cout<<"ERROR: unidentified data type"<<std::endl;
+				exit(1);
+			}
+		}
+
+		stim::envi_header out_hdr = header;							
+		out_hdr.bands = B;											//set the number of bands in the output header
+		out_hdr.save(outfile + ".hdr");								//save the output header
+	}
+
 	/// Converts ENVI files between interleave types (BSQ, BIL, and BIP)
 	/// @param outfile is the file name for the converted output
 	/// @param interleave is the interleave format for the destination file
-	bool convert(std::string outfile, stim::envi_header::interleaveType interleave){
+	bool convert(std::string outfile, stim::envi_header::interleaveType interleave, bool PROGRESS = false){
 		if(header.interleave == envi_header::BSQ){			//if the infile is bsq file
@@ -308,9 +362,12 @@ public:
 					exit(1);
 				}
 				else if(interleave == envi_header::BIL)			//if the target file is bil file
-					return ((bsq<float>*)file)->bil(outfile);
-				else if(interleave == envi_header::BIP)			//if the target file is bip file
-					return ((bsq<float>*)file)->bip(outfile);
+					return ((bsq<float>*)file)->bil(outfile, PROGRESS);
+				else if(interleave == envi_header::BIP){			//if the target file is bip file
+					std::cout<<"ERROR: conversion from BSQ to BIP isn't practical; use BSQ->BIL->BIP instead"<<std::endl;
+					//return ((bsq<float>*)file)->bip(outfile, PROGRESS);
+					exit(1);
+				}
 			}
 			else if(header.data_type == envi_header::float64){		//if the data type is float
@@ -319,9 +376,12 @@ public:
 					exit(1);
 				}
 				else if(interleave == envi_header::BIL)
-					return ((bsq<double>*)file)->bil(outfile);
-				else if(interleave == envi_header::BIP)
-					return ((bsq<double>*)file)->bip(outfile);
+					return ((bsq<double>*)file)->bil(outfile, PROGRESS);
+				else if(interleave == envi_header::BIP){			//if the target file is bip file
+					std::cout<<"ERROR: conversion from BSQ to BIP isn't practical; use BSQ->BIL->BIP instead"<<std::endl;
+					//return ((bsq<float>*)file)->bip(outfile, PROGRESS);
+					exit(1);
+				}
 			}
 			else{
@@ -338,9 +398,9 @@ public:
 					exit(1);
 				}
 				else if(interleave == envi_header::BSQ)			//if the target file is bsq file
-					return ((bil<float>*)file)->bsq(outfile);
+					return ((bil<float>*)file)->bsq(outfile, PROGRESS);
 				else if(interleave == envi_header::BIP)			//if the target file is bip file
-					return ((bil<float>*)file)->bip(outfile);
+					return ((bil<float>*)file)->bip(outfile, PROGRESS);
 			}
 			else if(header.data_type == envi_header::float64){		//if the data type is float
@@ -349,9 +409,9 @@ public:
 					exit(1);
 				}
 				else if(interleave == envi_header::BSQ)
-					return ((bil<double>*)file)->bsq(outfile);
+					return ((bil<double>*)file)->bsq(outfile, PROGRESS);
 				else if(interleave == envi_header::BIP)
-					return ((bil<double>*)file)->bip(outfile);
+					return ((bil<double>*)file)->bip(outfile, PROGRESS);
 			}
 			else{
@@ -368,9 +428,12 @@ public:
 					exit(1);
 				}
 				else if(interleave == envi_header::BIL)			//if the target file is bil file
-					return ((bip<float>*)file)->bil(outfile);
-				else if(interleave == envi_header::BSQ)			//if the target file is bsq file
-					return ((bip<float>*)file)->bsq(outfile);
+					return ((bip<float>*)file)->bil(outfile, PROGRESS);
+				else if(interleave == envi_header::BSQ){			//if the target file is bip file
+					std::cout<<"ERROR: conversion from BIP to BSQ isn't practical; use BIP->BIL->BSQ instead"<<std::endl;
+					//return ((bsq<float>*)file)->bip(outfile, PROGRESS);
+					exit(1);
+				}
 			}
 			else if(header.data_type == envi_header::float64){		//if the data type is float
@@ -379,9 +442,12 @@ public:
 					exit(1);
 				}
 				else if(interleave == envi_header::BIL)			//if the target file is bil file
-					return ((bip<double>*)file)->bil(outfile);
-				else if(interleave == envi_header::BSQ)			//if the target file is bsq file
-					return ((bip<double>*)file)->bsq(outfile);
+					return ((bip<double>*)file)->bil(outfile, PROGRESS);
+				else if(interleave == envi_header::BSQ){			//if the target file is bip file
+					std::cout<<"ERROR: conversion from BIP to BSQ isn't practical; use BIP->BIL->BSQ instead"<<std::endl;
+					//return ((bsq<float>*)file)->bip(outfile, PROGRESS);
+					exit(1);
+				}
 			}
 			else{
@@ -872,13 +938,13 @@ public:
 	/// @param p is a pointer to an allocated region of memory at least X * Y * sizeof(T) in size.
 	/// @param wavelength is a floating point value (usually a wavelength in spectral data) used as a label for the band to be copied.
-	bool band(void* ptr, double wavelength){
+	bool band(void* ptr, double wavelength, bool PROGRESS = false){
 		if(header.interleave == envi_header::BSQ){		//if the infile is bsq file
 			if(header.data_type ==envi_header::float32)
-				return ((bsq<float>*)file)->band((float*)ptr, wavelength);
+				return ((bsq<float>*)file)->band((float*)ptr, wavelength, PROGRESS);
 			else if (header.data_type == envi_header::float64)
-				return ((bsq<double>*)file)->band((double*)ptr, wavelength);
+				return ((bsq<double>*)file)->band((double*)ptr, wavelength, PROGRESS);
 			else{
 				std::cout << "ERROR: unidentified data type" << std::endl;
 				exit(1);
@@ -886,9 +952,9 @@ public:
 		}
 		else if (header.interleave == envi_header::BIL){
 			if (header.data_type == envi_header::float32)
-				return ((bil<float>*)file)->band((float*)ptr, wavelength);
+				return ((bil<float>*)file)->band((float*)ptr, wavelength, PROGRESS);
 			else if (header.data_type == envi_header::float64)
-				return ((bil<double>*)file)->band((double*)ptr, wavelength);
+				return ((bil<double>*)file)->band((double*)ptr, wavelength, PROGRESS);
 			else{
 				std::cout << "ERROR: unidentified data type" << std::endl;
 				exit(1);
@@ -896,9 +962,9 @@ public:
 		}
 		else if (header.interleave == envi_header::BIP){
 			if (header.data_type == envi_header::float32)
-				return ((bip<float>*)file)->band((float*)ptr, wavelength);
+				return ((bip<float>*)file)->band((float*)ptr, wavelength, PROGRESS);
 			else if (header.data_type == envi_header::float64)
-				return ((bip<double>*)file)->band((double*)ptr, wavelength);
+				return ((bip<double>*)file)->band((double*)ptr, wavelength, PROGRESS);
 			else{
 				std::cout << "ERROR: unidentified data type" << std::endl;
 				exit(1);
@@ -912,13 +978,13 @@ public:
 	/// @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){
+	bool spectrum(void* ptr, unsigned int x, unsigned int y, bool PROGRESS = false){
 		if(header.interleave == envi_header::BSQ){		//if the infile is bsq file
 			if(header.data_type ==envi_header::float32)
-				return ((bsq<float>*)file)->spectrum((float*)ptr, x, y);
+				return ((bsq<float>*)file)->spectrum((float*)ptr, x, y, PROGRESS);
 			else if (header.data_type == envi_header::float64)
-				return ((bsq<double>*)file)->spectrum((double*)ptr, x, y);
+				return ((bsq<double>*)file)->spectrum((double*)ptr, x, y, PROGRESS);
 			else{
 				std::cout << "ERROR: unidentified data type" << std::endl;
 				exit(1);
@@ -926,9 +992,9 @@ public:
 		}
 		else if (header.interleave == envi_header::BIL){
 			if (header.data_type == envi_header::float32)
-				return ((bil<float>*)file)->spectrum((float*)ptr, x, y);
+				return ((bil<float>*)file)->spectrum((float*)ptr, x, y, PROGRESS);
 			else if (header.data_type == envi_header::float64)
-				return ((bil<double>*)file)->spectrum((double*)ptr, x, y);
+				return ((bil<double>*)file)->spectrum((double*)ptr, x, y, PROGRESS);
 			else{
 				std::cout << "ERROR: unidentified data type" << std::endl;
 				exit(1);
@@ -936,9 +1002,9 @@ public:
 		}
 		else if (header.interleave == envi_header::BIP){
 			if (header.data_type == envi_header::float32)
-				return ((bip<float>*)file)->spectrum((float*)ptr, x, y);
+				return ((bip<float>*)file)->spectrum((float*)ptr, x, y, PROGRESS);
 			else if (header.data_type == envi_header::float64)
-				return ((bip<double>*)file)->spectrum((double*)ptr, x, y);
+				return ((bip<double>*)file)->spectrum((double*)ptr, x, y, PROGRESS);
 			else{
 				std::cout << "ERROR: unidentified data type" << std::endl;
 				exit(1);
@@ -989,12 +1055,12 @@ public:
 	/// @param p is a pointer to pre-allocated memory of size [B * sizeof(T)] that stores the mean spectrum
 	/// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
-	bool avg_band(double * p, unsigned char* mask){
+	bool avg_band(double * p, unsigned char* mask, bool PROGRESS = false){
 		if (header.interleave == envi_header::BSQ){
 			if (header.data_type == envi_header::float32)
-				return ((bsq<float>*)file)->avg_band(p, mask);
+				return ((bsq<float>*)file)->avg_band(p, mask, PROGRESS);
 			else if (header.data_type == envi_header::float64)
-				return ((bsq<double>*)file)->avg_band(p,  mask);
+				return ((bsq<double>*)file)->avg_band(p,  mask, PROGRESS);
 			else{
 				std::cout << "ERROR: unidentified data type" << std::endl;
 				exit(1);
@@ -1002,9 +1068,9 @@ public:
 		}
 		else if (header.interleave == envi_header::BIL){
 			if (header.data_type == envi_header::float32)
-				return ((bil<float>*)file)->avg_band(p,  mask);
+				return ((bil<float>*)file)->avg_band(p,  mask, PROGRESS);
 			else if (header.data_type == envi_header::float64)
-				return ((bil<double>*)file)->avg_band(p,  mask);
+				return ((bil<double>*)file)->avg_band(p,  mask, PROGRESS);
 			else{
 				std::cout << "ERROR: unidentified data type" << std::endl;
 				exit(1);
@@ -1012,9 +1078,9 @@ public:
 		}
 		else if (header.interleave == envi_header::BIP){
 			if (header.data_type == envi_header::float32)
-				return ((bip<float>*)file)->avg_band(p, mask);
+				return ((bip<float>*)file)->avg_band(p, mask, PROGRESS);
 			else if (header.data_type == envi_header::float64)
-				return ((bip<double>*)file)->avg_band(p, mask);
+				return ((bip<double>*)file)->avg_band(p, mask, PROGRESS);
 			else{
 				std::cout << "ERROR: unidentified data type" << std::endl;
 				exit(1);
@@ -1028,22 +1094,24 @@ public:
 	/// @param co is a pointer to pre-allocated memory of size [B * B] that stores the resulting covariance matrix
 	/// @param avg is a pointer to memory of size B that stores the average spectrum
 	/// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
-	bool co_matrix(double* co, double* avg, unsigned char* mask){
+	bool co_matrix(double* co, double* avg, unsigned char* mask, bool PROGRESS = false){
 		if (header.interleave == envi_header::BSQ){
-			if (header.data_type == envi_header::float32)
-				return ((bsq<float>*)file)->co_matrix(co, avg, mask);
+			std::cout<<"ERROR: calculating the covariance matrix for a BSQ file is impractical; convert to BIL or BIP first"<<std::endl;
+			exit(1);
+			/*if (header.data_type == envi_header::float32)
+				return ((bsq<float>*)file)->co_matrix(co, avg, mask, PROGRESS);
 			else if (header.data_type == envi_header::float64)
-				return ((bsq<double>*)file)->co_matrix(co, avg, mask);
+				return ((bsq<double>*)file)->co_matrix(co, avg, mask, PROGRESS);
 			else{
 				std::cout << "ERROR: unidentified data type" << std::endl;
 				exit(1);
-			}
+			}*/
 		}
 		else if (header.interleave == envi_header::BIL){
 			if (header.data_type == envi_header::float32)
-				return ((bil<float>*)file)->co_matrix(co, avg, mask);
+				return ((bil<float>*)file)->co_matrix(co, avg, mask, PROGRESS);
 			else if (header.data_type == envi_header::float64)
-				return ((bil<double>*)file)->co_matrix(co, avg, mask);
+				return ((bil<double>*)file)->co_matrix(co, avg, mask, PROGRESS);
 			else{
 				std::cout << "ERROR: unidentified data type" << std::endl;
 				exit(1);
@@ -1051,9 +1119,9 @@ public:
 		}
 		else if (header.interleave == envi_header::BIP){
 			if (header.data_type == envi_header::float32)
-				return ((bip<float>*)file)->co_matrix(co, avg, mask);
+				return ((bip<float>*)file)->co_matrix(co, avg, mask, PROGRESS);
 			else if (header.data_type == envi_header::float64)
-				return ((bip<double>*)file)->co_matrix(co, avg, mask);
+				return ((bip<double>*)file)->co_matrix(co, avg, mask, PROGRESS);
 			else{
 				std::cout << "ERROR: unidentified data type" << std::endl;
 				exit(1);
@@ -1070,7 +1138,7 @@ public:
 	/// @param y0 is the lower-left y pixel coordinate to be included in the cropped image
 	/// @param x1 is the upper-right x pixel coordinate to be included in the cropped image
 	/// @param y1 is the upper-right y pixel coordinate to be included in the cropped image
-	bool crop(std::string outfile,unsigned x0, unsigned y0, unsigned x1, unsigned y1, unsigned b0, unsigned b1){
+	bool crop(std::string outfile,unsigned x0, unsigned y0, unsigned x1, unsigned y1, unsigned b0, unsigned b1, bool PROGRESS = false){
 		//save the header for the cropped file
 		stim::envi_header new_header = header;
@@ -1084,9 +1152,9 @@ public:
 		if (header.interleave == envi_header::BSQ){
 			if (header.data_type == envi_header::float32)
-				return ((bsq<float>*)file)->crop(outfile, x0, y0, x1, y1, b0, b1);
+				return ((bsq<float>*)file)->crop(outfile, x0, y0, x1, y1, b0, b1, PROGRESS);
 			else if (header.data_type == envi_header::float64)
-				return ((bsq<double>*)file)->crop(outfile, x0, y0, x1, y1, b0, b1);
+				return ((bsq<double>*)file)->crop(outfile, x0, y0, x1, y1, b0, b1, PROGRESS);
 			else{
 				std::cout << "ERROR: unidentified data type" << std::endl;
 				exit(1);
@@ -1094,9 +1162,9 @@ public:
 		}
 		else if (header.interleave == envi_header::BIL){
 			if (header.data_type == envi_header::float32)
-				return ((bil<float>*)file)->crop(outfile, x0, y0, x1, y1, b0, b1);
+				return ((bil<float>*)file)->crop(outfile, x0, y0, x1, y1, b0, b1, PROGRESS);
 			else if (header.data_type == envi_header::float64)
-				return ((bil<double>*)file)->crop(outfile, x0, y0, x1, y1, b0, b1);
+				return ((bil<double>*)file)->crop(outfile, x0, y0, x1, y1, b0, b1, PROGRESS);
 			else{
 				std::cout << "ERROR: unidentified data type" << std::endl;
 				exit(1);
@@ -1104,9 +1172,9 @@ public:
 		}
 		else if (header.interleave == envi_header::BIP){
 			if (header.data_type == envi_header::float32)
-				return ((bip<float>*)file)->crop(outfile, x0, y0, x1, y1, b0, b1);
+				return ((bip<float>*)file)->crop(outfile, x0, y0, x1, y1, b0, b1, PROGRESS);
 			else if (header.data_type == envi_header::float64)
-				return ((bip<double>*)file)->crop(outfile, x0, y0, x1, y1, b0, b1);
+				return ((bip<double>*)file)->crop(outfile, x0, y0, x1, y1, b0, b1, PROGRESS);
 			else{
 				std::cout << "ERROR: unidentified data type" << std::endl;
 				exit(1);