fixed metric calculations, added a finite-value mask, as_float() in stim::arglist returns a double

David Mayerich
1 parent 44f44bb5
Showing 6 changed files with 265 additions and 118 deletions Show diff stats
stim/envi/bil.h
stim/envi/bip.h
stim/envi/bsq.h
stim/envi/envi.h
stim/envi/hsi.h
stim/parser/arguments.h
@@ -600,7 +600,7 @@ public:
 	/// @param rb2 is the label value for the right baseline point for the second peak (denominator)
 	/// @param pos2 is the label value for the second peak (denominator) position
 	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
-	bool ph_to_ph(double lb1, double rb1, double pos1, double lb2, double rb2, double pos2, T * result){
+	bool ph_to_ph(T* result, double lb1, double rb1, double pos1, double lb2, double rb2, double pos2, unsigned char* mask = NULL){
 		T* p1 = (T*)malloc(X() * Y() * sizeof(T));
 		T* p2 = (T*)malloc(X() * Y() * sizeof(T));
@@ -630,8 +630,8 @@ public:
 	/// @param rb2 is the label value for the right baseline point for the second peak (denominator)
 	/// @param pos2 is the label value for the second peak (denominator) position
 	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
-	bool pa_to_ph(double lb1, double rb1, double lab1, double rab1,
-					double lb2, double rb2, double pos, T* result){
+	bool pa_to_ph(T* result, double lb1, double rb1, double lab1, double rab1,
+					double lb2, double rb2, double pos, unsigned char* mask = NULL){
 		T* p1 = (T*)malloc(X() * Y() * sizeof(T));
 		T* p2 = (T*)malloc(X() * Y() * sizeof(T));
@@ -663,8 +663,8 @@ public:
 	/// @param lab2 is the label value for the left bound (start of the integration) of the second peak (denominator)
 	/// @param rab2 is the label value for the right bound (end of the integration) of the second peak (denominator)	
 	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
-	bool pa_to_pa(double lb1, double rb1, double lab1, double rab1,
-					double lb2, double rb2, double lab2, double rab2, T* result){
+	bool pa_to_pa(T* result, double lb1, double rb1, double lab1, double rab1,
+					double lb2, double rb2, double lab2, double rab2, unsigned char* mask = NULL){
 		T* p1 = (T*)malloc(X() * Y() * sizeof(T));
 		T* p2 = (T*)malloc(X() * Y() * sizeof(T));
@@ -773,7 +773,7 @@ public:
 	/// @param lab is the label for the start of the peak
 	/// @param rab is the label for the end of the peak
 	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
-	bool cpoint(double lb, double rb, double lab, double rab, T* result){
+	bool centroid(T* result, double lb, double rb, double lab, double rab, unsigned char* mask = NULL){
 		T* p1 = (T*)malloc(X() * Y() * sizeof(T));
 		T* p2 = (T*)malloc(X() * Y() * sizeof(T));
@@ -782,9 +782,7 @@ public:
 		area(lb, rb, lab, rab, p2);
 		//calculate the ratio in result
 		for(unsigned long long i = 0; i < X() * Y(); i++){
-			if(p1[i] == 0 && p2[i] ==0)
-				result[i] = 1;
-			else
+			if(mask == NULL || mask[i])
 				result[i] = p1[i] / p2[i];
 		}
@@ -800,16 +798,16 @@ public:
 	/// @param mask_band is the band used to specify the mask
 	/// @param threshold is the threshold used to determine if the mask value is true or false
 	/// @param p is a pointer to a pre-allocated array at least X * Y in size
-	bool build_mask(double mask_band, double threshold, unsigned char* p, bool PROGRESS = false){
+	bool build_mask(unsigned char* mask, double mask_band, double threshold, bool PROGRESS = false){
 		T* temp = (T*)malloc(X() * Y() * sizeof(T));		//allocate memory for the certain band
 		band(temp, mask_band, PROGRESS);
 		for (unsigned long long i = 0; i < X() * Y(); i++) {
 			if (temp[i] < threshold)
-				p[i] = 0;
+				mask[i] = 0;
 			else
-				p[i] = 255;
+				mask[i] = 255;
 		}
 		free(temp);
@@ -850,6 +848,36 @@ public:
 		return true;
 	}
+	/// Copies all spectra corresponding to nonzero values of a mask into a pre-allocated matrix of size (B x P)
+	///		where P is the number of masked pixels and B is the number of bands. The allocated memory can be accessed
+	///		using the following indexing: i = p*B + b
+	/// @param matrix is the destination for the pixel data
+	/// @param mask is the mask
+	bool sift(T* matrix, unsigned char* mask){
+		size_t Lbytes = sizeof(T) * X();
+		T* line = (T*) malloc( Lbytes );					//allocate space for a line
+
+		file.seekg(0, std::ios::beg);	//seek to the beginning of the file
+		
+		size_t pl;
+		size_t p = 0;										//create counter variables
+		for(size_t y = 0; y < Y(); y++){					//for each line in the data set
+			for(size_t b = 0; b < Z(); b++){				//for each band in the data set
+				pl = 0;										//initialize the pixel offset for the current line to zero (0)
+				file.read( (char*)line, Lbytes );			//read the current line
+				for(size_t x = 0; x < X(); x++){
+					if(mask[y * X() + x]){					//if the current pixel is masked
+						size_t i = (p + pl) * Z() + b;		//calculate the index in the sifted matrix
+						matrix[i] = line[x];				//store the current value in the line at the correct matrix location
+						pl++;								//increment the pixel pointer
+					}
+				}
+			}
+			p += pl;										//add the line increment to the running pixel index
+		}
+		return true;
+	}
+
 	/// 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, bool PROGRESS = false){
@@ -526,7 +526,7 @@ public:
 	/// @param rb2 is the label value for the right baseline point for the second peak (denominator)
 	/// @param pos2 is the label value for the second peak (denominator) position
 	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
-	bool ph_to_ph(double lb1, double rb1, double pos1, double lb2, double rb2, double pos2, T * result){
+	bool ph_to_ph(T* result, double lb1, double rb1, double pos1, double lb2, double rb2, double pos2, unsigned char* mask = NULL){
 		T* p1 = (T*)malloc(X() * Y() * sizeof(T));	
 		T* p2 = (T*)malloc(X() * Y() * sizeof(T));	
@@ -556,8 +556,8 @@ public:
 	/// @param rb2 is the label value for the right baseline point for the second peak (denominator)
 	/// @param pos2 is the label value for the second peak (denominator) position
 	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
-	bool pa_to_ph(double lb1, double rb1, double lab1, double rab1,
-					double lb2, double rb2, double pos, T* result){
+	bool pa_to_ph(T* result, double lb1, double rb1, double lab1, double rab1,
+					double lb2, double rb2, double pos, unsigned char* mask = NULL){
 		T* p1 = (T*)malloc(X() * Y() * sizeof(T));	
 		T* p2 = (T*)malloc(X() * Y() * sizeof(T));	
@@ -589,8 +589,8 @@ public:
 	/// @param lab2 is the label value for the left bound (start of the integration) of the second peak (denominator)
 	/// @param rab2 is the label value for the right bound (end of the integration) of the second peak (denominator)	
 	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
-	bool pa_to_pa(double lb1, double rb1, double lab1, double rab1,
-					double lb2, double rb2, double lab2, double rab2, T* result){
+	bool pa_to_pa(T* result, double lb1, double rb1, double lab1, double rab1,
+					double lb2, double rb2, double lab2, double rab2, unsigned char* mask = NULL){
 		T* p1 = (T*)malloc(X() * Y() * sizeof(T));	
 		T* p2 = (T*)malloc(X() * Y() * sizeof(T));	
@@ -699,7 +699,7 @@ public:
 	/// @param lab is the label for the start of the peak
 	/// @param rab is the label for the end of the peak
 	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
-	bool cpoint(double lb, double rb, double lab, double rab, T* result){
+	bool centroid(T* result, double lb, double rb, double lab, double rab, unsigned char* mask = NULL){
 		T* p1 = (T*)malloc(X() * Y() * sizeof(T));	
 		T* p2 = (T*)malloc(X() * Y() * sizeof(T));	
@@ -708,9 +708,7 @@ public:
 		area(lb, rb, lab, rab, p2);
 		//calculate the ratio in result
 		for(unsigned long long i = 0; i < X() * Y(); i++){
-			if(p1[i] == 0 && p2[i] ==0)
-				result[i] = 1;
-			else
+			if(mask == NULL || mask[i])
 				result[i] = p1[i] / p2[i];
 		}
@@ -726,16 +724,16 @@ public:
 	/// @param mask_band is the band used to specify the mask
 	/// @param threshold is the threshold used to determine if the mask value is true or false
 	/// @param p is a pointer to a pre-allocated array at least X * Y in size
-	bool build_mask(double mask_band, double threshold, unsigned char* p, bool PROGRESS = false){
+	bool build_mask(unsigned char* mask, double mask_band, double threshold, bool PROGRESS = false){
 		T* temp = (T*)malloc(X() * Y() * sizeof(T));		//allocate memory for the certain band
 		band(temp, mask_band, PROGRESS);
 		for (unsigned long long i = 0; i < X() * Y();i++) {
 			if (temp[i] < threshold)
-				p[i] = 0;
+				mask[i] = 0;
 			else
-				p[i] = 255;
+				mask[i] = 255;
 		}
 		free(temp);
@@ -776,6 +774,33 @@ public:
 		return true;						//return true
 	}
+	/// Copies all spectra corresponding to nonzero values of a mask into a pre-allocated matrix of size (B x P)
+	///		where P is the number of masked pixels and B is the number of bands. The allocated memory can be accessed
+	///		using the following indexing: i = p*B + b
+	/// @param matrix is the destination for the pixel data
+	/// @param mask is the mask
+	bool sift(T* matrix, unsigned char* mask){
+		size_t Bbytes = sizeof(T) * Z();
+		size_t XY = X() * Y();
+		T* band = (T*) malloc( Bbytes );					//allocate space for a line
+
+		file.seekg(0, std::ios::beg);	//seek to the beginning of the file
+		
+		size_t p = 0;										//create counter variables
+		for(size_t xy = 0; xy < XY; xy++){					//for each pixel
+			if(mask[xy]){									//if the current pixel is masked
+				file.read( (char*)band, Bbytes );			//read the current line
+				for(size_t b = 0; b < Z(); b++){			//copy each band value to the sifted matrix
+					size_t i = p * Z() + b;					//calculate the index in the sifted matrix
+					matrix[i] = band[b];					//store the current value in the line at the correct matrix location
+				}
+				p++;									//increment the pixel pointer
+			}
+			else
+				file.seekg(Bbytes, std::ios::cur);			//otherwise skip this band			
+		}
+		return true;
+	}
 	/// Saves to disk only those spectra corresponding to mask values != 0
 	bool sift(std::string outfile, unsigned char* mask, bool PROGRESS = false){
@@ -436,9 +436,7 @@ public:
 		cur = (T*) malloc(S);
 		cur2 = (T*) malloc(S);
-		memset(result, (char)0, S);
-
-		//find the wavelenght position in the whole band
+		//find the wavelength position in the whole band
 		unsigned long long n = w.size();
 		unsigned long long ai = 0;		//left bound position
 		unsigned long long bi = n - 1;		//right bound position
@@ -450,26 +448,30 @@ public:
 			std::cout<<"ERROR: left bound or right bound out of bandwidth"<<std::endl;
 			exit(1);
 		}
-		//to make sure rigth bound is bigger than left bound
+		//to make sure right bound is bigger than left bound
 		else if(lb > rb){
 			std::cout<<"ERROR: right bound should be bigger than left bound"<<std::endl;
 			exit(1);
 		}
-		//get the position of lb and rb
-		while (lab >= w[ai]){
+		//find the indices of the left and right baseline points
+		while (lab >= w[ai]){ 
 			ai++;
 		}
 		while (rab <= w[bi]){
 			bi--;
 		}
-		band(lp, lb);
+		band(lp, lb);						//get the band images for the left and right baseline points
 		band(rp, rb);
-		//calculate the beginning and the ending part
-		baseline_band(lb, rb, lp, rp, rab, cur2);		//ending part
-		baseline_band(lb, rb, lp, rp, w[bi], cur);
+		// calculate the average value of the indexed region
+		memset(result, 0, S);								//initialize the integral to zero (0)
+
+		//integrate the region between the specified bands and the closest indexed band
+		//		this integrates the "tails" of the spectrum that lie outside the main indexed region
+		baseline_band(lb, rb, lp, rp, rab, cur2);		//calculate the image for the right-most band in the integral
+		baseline_band(lb, rb, lp, rp, w[bi], cur);		//calculate the image for the right-most indexed band
 		for(unsigned long long j = 0; j < XY; j++){
 			result[j] += (T)((rab - w[bi]) * ((double)cur[j] + (double)cur2[j]) / 2.0);
 		}
@@ -479,13 +481,12 @@ public:
 			result[j] += (T)((w[ai] - lab) * ((double)cur[j] + (double)cur2[j]) / 2.0);
 		}
-		//calculate the area
+		//integrate the main indexed region
 		ai++;
-		for(unsigned long long i = ai; i <= bi ;i++)
+		for(unsigned long long i = ai; i <= bi ;i++)	//for each band in the integral
 		{
-			baseline_band(lb, rb, lp, rp, w[ai], cur2);
-			for(unsigned long long j = 0; j < XY; j++)
-			{
+			baseline_band(lb, rb, lp, rp, w[ai], cur2);	//get the baselined band
+			for(unsigned long long j = 0; j < XY; j++){	//for each pixel in that band
 				result[j] += (T)((w[ai] - w[ai-1]) * ((double)cur[j] + (double)cur2[j]) / 2.0);
 			}
 			std::swap(cur,cur2);		//swap the band pointers
@@ -507,20 +508,22 @@ public:
 	/// @param rb2 is the label value for the right baseline point for the second peak (denominator)
 	/// @param pos2 is the label value for the second peak (denominator) position
 	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
-	bool ph_to_ph(double lb1, double rb1, double pos1, double lb2, double rb2, double pos2, T * result){
+	bool ph_to_ph(T* result, double lb1, double rb1, double pos1, double lb2, double rb2, double pos2, unsigned char* mask = NULL){
-		T* p1 = (T*)malloc(X() * Y() * sizeof(T));
-		T* p2 = (T*)malloc(X() * Y() * sizeof(T));
+		size_t XYbytes = X() * Y() * sizeof(T);			//calculate the size of the band image (in bytes)
+
+		T* p1 = (T*)malloc(XYbytes);					//allocate space for both bands in the ratio
+		T* p2 = (T*)malloc(XYbytes);
+		memset(result, 0, XYbytes);						//initialize the ratio to zero
 		//get the two peak band
 		height(lb1, rb1, pos1, p1);
 		height(lb2, rb2, pos2, p2);
 		//calculate the ratio in result
 		for(unsigned long long i = 0; i < X() * Y(); i++){
-			if(p1[i] == 0 && p2[i] ==0)
-				result[i] = 1;
-			else
+			if(mask == NULL || mask[i]){
 				result[i] = p1[i] / p2[i];
+			}
 		}
 		free(p1);
@@ -537,20 +540,20 @@ public:
 	/// @param rb2 is the label value for the right baseline point for the second peak (denominator)
 	/// @param pos2 is the label value for the second peak (denominator) position
 	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
-	bool pa_to_ph(double lb1, double rb1, double lab1, double rab1,
-					double lb2, double rb2, double pos, T* result){
+	bool pa_to_ph(T* result, double lb1, double rb1, double lab1, double rab1,
+					double lb2, double rb2, double pos, unsigned char* mask = NULL){
-		T* p1 = (T*)malloc(X() * Y() * sizeof(T));
-		T* p2 = (T*)malloc(X() * Y() * sizeof(T));
+		size_t bytes = X() * Y() * sizeof(T);
+		T* p1 = (T*)malloc(bytes);							//allocate space for both ratio components
+		T* p2 = (T*)malloc(bytes);
+		memset(result, 0, bytes);							//initialize the ratio to zero
 		//get the area and the peak band
 		area(lb1, rb1, lab1, rab1, p1);
 		height(lb2, rb2, pos, p2);
 		//calculate the ratio in result
 		for(unsigned long long i = 0; i < X() * Y(); i++){
-			if(p1[i] == 0 && p2[i] ==0)
-				result[i] = 1;
-			else
+			if(mask == NULL || mask[i])
 				result[i] = p1[i] / p2[i];
 		}
@@ -570,21 +573,21 @@ public:
 	/// @param lab2 is the label value for the left bound (start of the integration) of the second peak (denominator)
 	/// @param rab2 is the label value for the right bound (end of the integration) of the second peak (denominator)
 	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
-	bool pa_to_pa(double lb1, double rb1, double lab1, double rab1,
-					double lb2, double rb2, double lab2, double rab2, T* result){
+	bool pa_to_pa(T* result, double lb1, double rb1, double lab1, double rab1,
+					double lb2, double rb2, double lab2, double rab2, unsigned char* mask = NULL){
-		T* p1 = (T*)malloc(X() * Y() * sizeof(T));
-		T* p2 = (T*)malloc(X() * Y() * sizeof(T));
+		size_t bytes = X() * Y() * sizeof(T);
+		T* p1 = (T*)malloc(bytes);						//allocate space for each of the operands
+		T* p2 = (T*)malloc(bytes);
+		memset(result, 0, bytes);						//initialize the ratio result to zero (0)
 		//get the area and the peak band
 		area(lb1, rb1, lab1, rab1, p1);
 		area(lb2, rb2, lab2, rab2, p2);
 		//calculate the ratio in result
 		for(unsigned long long i = 0; i < X() * Y(); i++){
-			if(p1[i] == 0 && p2[i] ==0)
-				result[i] = 1;
-			else
-				result[i] = p1[i] / p2[i];
+			if(mask == NULL || mask[i])						//if the pixel is masked
+				result[i] = p1[i] / p2[i];					//calculate the ratio
 		}
 		free(p1);
@@ -592,7 +595,7 @@ public:
 		return true;
 	}
-	/// Compute the definite integral of a baseline corrected peak.
+	/// Compute the definite integral of a baseline corrected peak weighted by the corresponding wavelength
 	/// @param lb is the label value for the left baseline point
 	/// @param rb is the label value for the right baseline point
@@ -613,8 +616,6 @@ public:
 		cur = (T*) malloc(S);
 		cur2 = (T*) malloc(S);
-		memset(result, (char)0, S);
-
 		//find the wavelenght position in the whole band
 		unsigned long long n = w.size();
 		unsigned long long ai = 0;		//left bound position
@@ -642,6 +643,8 @@ public:
 		band(lp, lb);
 		band(rp, rb);
+		memset(result, (char)0, S);						//initialize the integral to zero (0)
+
 		//calculate the beginning and the ending part
 		baseline_band(lb, rb, lp, rp, rab, cur2);		//ending part
 		baseline_band(lb, rb, lp, rp, w[bi], cur);
@@ -656,12 +659,11 @@ public:
 		//calculate f(x) times x
 		ai++;
-		for(unsigned long long i = ai; i <= bi ;i++)
-		{
+		for(unsigned long long i = ai; i <= bi ;i++){
 			baseline_band(lb, rb, lp, rp, w[ai], cur2);
-			for(unsigned long long j = 0; j < XY; j++)
-			{
-				result[j] += (T)((w[ai] - w[ai-1]) * (w[ai] + w[ai-1]) * ((double)cur[j] + (double)cur2[j]) / 4.0);
+			for(unsigned long long j = 0; j < XY; j++){
+				T v = (T)((w[ai] - w[ai-1]) * (w[ai] + w[ai-1]) * ((double)cur[j] + (double)cur2[j]) / 4.0);
+				result[j] += v;
 			}
 			std::swap(cur,cur2);		//swap the band pointers
 		}
@@ -674,25 +676,27 @@ public:
 	}
 	/// Compute the centroid of a baseline corrected peak.
+	/// Note that the values for the centroid can be outside of [lab, rab] if the spectrum goes negative.
 	/// @param lb is the label value for the left baseline point
 	/// @param rb is the label value for the right baseline point
 	/// @param lab is the label for the start of the peak
 	/// @param rab is the label for the end of the peak
 	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
-	bool cpoint(double lb, double rb, double lab, double rab, T* result){
-		T* p1 = (T*)malloc(X() * Y() * sizeof(T));
+	bool centroid(T* result, double lb, double rb, double lab, double rab, unsigned char* mask = NULL){
+		size_t bytes = X() * Y() * sizeof(T);		//calculate the number of bytes in a band image
+		T* p1 = (T*)malloc(X() * Y() * sizeof(T));	//allocate space for both operands
 		T* p2 = (T*)malloc(X() * Y() * sizeof(T));
+		memset(result, 0, bytes);					//initialize the ratio result to zero (0)
 		//get the area and the peak band
 		x_area(lb, rb, lab, rab, p1);
 		area(lb, rb, lab, rab, p2);
 		//calculate the ratio in result
 		for(unsigned long long i = 0; i < X() * Y(); i++){
-			if(p1[i] == 0 && p2[i] ==0)
-				result[i] = 1;
-			else
+			if(mask == NULL || mask[i]){
 				result[i] = p1[i] / p2[i];
+			}
 		}
 		free(p1);
@@ -707,16 +711,16 @@ public:
 	/// @param mask_band is the band used to specify the mask
 	/// @param threshold is the threshold used to determine if the mask value is true or false
 	/// @param p is a pointer to a pre-allocated array at least X * Y in size
-	bool build_mask(double mask_band, double threshold, unsigned char* p = NULL, bool PROGRESS = false){
+	bool build_mask(unsigned char* mask, double mask_band, double threshold, bool PROGRESS = false){
 		T* temp = (T*)malloc(X() * Y() * sizeof(T));		//allocate memory for the certain band
 		band(temp, mask_band);
 		for (unsigned long long i = 0; i < X() * Y(); i++) {
 			if (temp[i] < threshold)
-				p[i] = 0;
+				mask[i] = 0;
 			else
-				p[i] = 255;
+				mask[i] = 255;
 			if(PROGRESS) progress = (double) (i+1) / (X() * Y()) * 100;
 		}
@@ -475,31 +475,31 @@ public:
 	/// @param mask_band is the label for the band that will be used to build the mask
 	/// @param threshold is a value selected such that all band values greater than threshold will have a mask value of 'true'
 	/// @param p is memory of size X*Y that will store the resulting mask
-	bool build_mask(double mask_band, double threshold, unsigned char* p = NULL, bool PROGRESS = false)	{
+	bool build_mask(unsigned char* mask, double mask_band, double threshold, 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)->build_mask(mask_band, threshold, p, PROGRESS);
+				return ((bsq<float>*)file)->build_mask(mask, mask_band, threshold, PROGRESS);
 			else if(header.data_type == envi_header::float64)
-				return ((bsq<double>*)file)->build_mask(mask_band, threshold, p, PROGRESS);
+				return ((bsq<double>*)file)->build_mask(mask, mask_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)->build_mask(mask_band, threshold, p, PROGRESS);
+				return ((bil<float>*)file)->build_mask(mask, mask_band, threshold, PROGRESS);
 			else if(header.data_type == envi_header::float64)
-				return ((bil<double>*)file)->build_mask(mask_band, threshold, p, PROGRESS);
+				return ((bil<double>*)file)->build_mask(mask, mask_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)->build_mask(mask_band, threshold, p, PROGRESS);
+				return ((bip<float>*)file)->build_mask(mask, mask_band, threshold, PROGRESS);
 			else if(header.data_type == envi_header::float64)
-				return ((bip<double>*)file)->build_mask(mask_band, threshold, p, PROGRESS);
+				return ((bip<double>*)file)->build_mask(mask, mask_band, threshold, PROGRESS);
 			else
 				std::cout<<"ERROR: unidentified data type"<<std::endl;
 		}
@@ -507,6 +507,36 @@ public:
 		return false;
 	}
+	/// Creates a mask with a true value for all pixels that contain finite values
+	void mask_finite(unsigned char* mask, bool PROGRESS = false){
+		if(header.interleave == envi_header::BSQ){		//if the infile is bsq file
+			if(header.data_type ==envi_header::float32)
+				((bsq<float>*)file)->mask_finite(mask, PROGRESS);
+			else if(header.data_type == envi_header::float64)
+				((bsq<double>*)file)->mask_finite(mask, 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)->mask_finite(mask, PROGRESS);
+			else if(header.data_type == envi_header::float64)
+				((bil<double>*)file)->mask_finite(mask, 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)->mask_finite(mask, PROGRESS);
+			else if(header.data_type == envi_header::float64)
+				((bip<double>*)file)->mask_finite(mask, PROGRESS);
+			else
+				std::cout<<"ERROR: unidentified data type"<<std::endl;
+		}
+	}
+
 	/// Applies a mask to the ENVI file.
 	/// @param outfile is the name of the resulting masked output file
@@ -560,17 +590,31 @@ public:
 				return ((bsq<float>*)file)->sift((float*)matrix, p);
 			else if (header.data_type == envi_header::float64)
 				return ((bsq<double>*)file)->sift((double*)matrix, p);
-			else
+			else{
 				std::cout << "ERROR: unidentified data type" << std::endl;
+				exit(1);
+			}
 		}
 		if (header.interleave == envi_header::BIP){
-			std::cout<<"Memory sifting not supported for BIP files"<<std::endl;
-			exit(1);
+			if (header.data_type == envi_header::float32)
+				return ((bip<float>*)file)->sift((float*)matrix, p);
+			else if (header.data_type == envi_header::float64)
+				return ((bip<double>*)file)->sift((double*)matrix, p);
+			else{
+				std::cout << "ERROR: unidentified data type" << std::endl;
+				exit(1);
+			}
 		}
 		if (header.interleave == envi_header::BIL){
-			std::cout<<"Memory sifting not supported for BIL files"<<std::endl;
-			exit(1);
+			if (header.data_type == envi_header::float32)
+				return ((bil<float>*)file)->sift((float*)matrix, p);
+			else if (header.data_type == envi_header::float64)
+				return ((bil<double>*)file)->sift((double*)matrix, p);
+			else{
+				std::cout << "ERROR: unidentified data type" << std::endl;
+				exit(1);
+			}
 		}
 		return false;
@@ -684,30 +728,30 @@ public:
 	/// @param rb2 is the label value for the right baseline point for the second peak (denominator)
 	/// @param pos2 is the label value for the second peak (denominator) position
 	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
-	bool ph_to_ph(double lb1, double rb1, double pos1, double lb2, double rb2, double pos2, void * result){
+	bool ph_to_ph(void * result, double lb1, double rb1, double pos1, double lb2, double rb2, double pos2, unsigned char* mask){
 		if(header.interleave == envi_header::BSQ){		//if the infile is bsq file
 			if(header.data_type ==envi_header::float32)
-				return ((bsq<float>*)file)->ph_to_ph(lb1, rb1, pos1, lb2, rb2, pos2, (float*)result);
+				return ((bsq<float>*)file)->ph_to_ph((float*)result, lb1, rb1, pos1, lb2, rb2, pos2, mask);
 			else if(header.data_type == envi_header::float64)
-				return ((bsq<double>*)file)->ph_to_ph(lb1, rb1, pos1, lb2, rb2, pos2, (double*)result);
+				return ((bsq<double>*)file)->ph_to_ph((double*)result, lb1, rb1, pos1, lb2, rb2, pos2, mask);
 			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)->ph_to_ph(lb1, rb1, pos1, lb2, rb2, pos2, (float*)result);
+				return ((bil<float>*)file)->ph_to_ph((float*)result, lb1, rb1, pos1, lb2, rb2, pos2, mask);
 			else if(header.data_type == envi_header::float64)
-				return ((bil<double>*)file)->ph_to_ph(lb1, rb1, pos1, lb2, rb2, pos2, (double*)result);
+				return ((bil<double>*)file)->ph_to_ph((double*)result, lb1, rb1, pos1, lb2, rb2, pos2, mask);
 			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)->ph_to_ph(lb1, rb1, pos1, lb2, rb2, pos2, (float*)result);
+				return ((bip<float>*)file)->ph_to_ph((float*)result, lb1, rb1, pos1, lb2, rb2, pos2, mask);
 			else if(header.data_type == envi_header::float64)
-				return ((bip<double>*)file)->ph_to_ph(lb1, rb1, pos1, lb2, rb2, pos2, (double*)result);
+				return ((bip<double>*)file)->ph_to_ph((double*)result, lb1, rb1, pos1, lb2, rb2, pos2, mask);
 			else
 				std::cout<<"ERROR: unidentified data type"<<std::endl;
 		}
@@ -728,30 +772,30 @@ public:
 	/// @param rb2 is the label value for the right baseline point for the second peak (denominator)
 	/// @param pos2 is the label value for the second peak (denominator) position
 	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
-	bool pa_to_ph(double lb1, double rb1, double lab1, double rab1, double lb2, double rb2, double pos, void * result){
+	bool pa_to_ph(void* result, double lb1, double rb1, double lab1, double rab1, double lb2, double rb2, double pos, unsigned char* mask = NULL){
 		if(header.interleave == envi_header::BSQ){		//if the infile is bsq file
 			if(header.data_type ==envi_header::float32)
-				return ((bsq<float>*)file)->pa_to_ph(lb1, rb1, lab1, rab1, lb2, rb2, pos, (float*)result);
+				return ((bsq<float>*)file)->pa_to_ph((float*)result, lb1, rb1, lab1, rab1, lb2, rb2, pos, mask);
 			else if(header.data_type == envi_header::float64)
-				return ((bsq<double>*)file)->pa_to_ph(lb1, rb1, lab1, rab1, lb2, rb2, pos, (double*)result);
+				return ((bsq<double>*)file)->pa_to_ph((double*)result, lb1, rb1, lab1, rab1, lb2, rb2, pos, mask);
 			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)->pa_to_ph(lb1, rb1, lab1, rab1, lb2, rb2, pos, (float*)result);
+				return ((bil<float>*)file)->pa_to_ph((float*)result, lb1, rb1, lab1, rab1, lb2, rb2, pos, mask);
 			else if(header.data_type == envi_header::float64)
-				return ((bil<double>*)file)->pa_to_ph(lb1, rb1, lab1, rab1, lb2, rb2, pos, (double*)result);
+				return ((bil<double>*)file)->pa_to_ph((double*)result, lb1, rb1, lab1, rab1, lb2, rb2, pos, mask);
 			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)->pa_to_ph(lb1, rb1, lab1, rab1, lb2, rb2, pos, (float*)result);
+				return ((bip<float>*)file)->pa_to_ph((float*)result, lb1, rb1, lab1, rab1, lb2, rb2, pos, mask);
 			else if(header.data_type == envi_header::float64)
-				return ((bip<double>*)file)->pa_to_ph(lb1, rb1, lab1, rab1, lb2, rb2, pos, (double*)result);
+				return ((bip<double>*)file)->pa_to_ph((double*)result, lb1, rb1, lab1, rab1, lb2, rb2, pos, mask);
 			else
 				std::cout<<"ERROR: unidentified data type"<<std::endl;
 		}
@@ -774,31 +818,31 @@ public:
 	/// @param lab2 is the label value for the left bound (start of the integration) of the second peak (denominator)
 	/// @param rab2 is the label value for the right bound (end of the integration) of the second peak (denominator)
 	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
-	bool pa_to_pa(double lb1, double rb1, double lab1, double rab1,
-					double lb2, double rb2, double lab2, double rab2, void* result){
+	bool pa_to_pa(void* result, double lb1, double rb1, double lab1, double rab1,
+					double lb2, double rb2, double lab2, double rab2, unsigned char* mask = NULL){
 		if(header.interleave == envi_header::BSQ){		//if the infile is bsq file
 			if(header.data_type ==envi_header::float32)
-				return ((bsq<float>*)file)->pa_to_pa(lb1, rb1, lab1, rab1, lb2, rb2, lab2, rab2, (float*)result);
+				return ((bsq<float>*)file)->pa_to_pa((float*)result, lb1, rb1, lab1, rab1, lb2, rb2, lab2, rab2, mask);
 			else if(header.data_type == envi_header::float64)
-				return ((bsq<double>*)file)->pa_to_pa(lb1, rb1, lab1, rab1, lb2, rb2, lab2, rab2, (double*)result);
+				return ((bsq<double>*)file)->pa_to_pa((double*)result, lb1, rb1, lab1, rab1, lb2, rb2, lab2, rab2, mask);
 			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)->pa_to_pa(lb1, rb1, lab1, rab1, lb2, rb2, lab2, rab2, (float*)result);
+				return ((bil<float>*)file)->pa_to_pa((float*)result, lb1, rb1, lab1, rab1, lb2, rb2, lab2, rab2, mask);
 			else if(header.data_type == envi_header::float64)
-				return ((bil<double>*)file)->pa_to_pa(lb1, rb1, lab1, rab1, lb2, rb2, lab2, rab2, (double*)result);
+				return ((bil<double>*)file)->pa_to_pa((double*)result, lb1, rb1, lab1, rab1, lb2, rb2, lab2, rab2, mask);
 			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)->pa_to_pa(lb1, rb1, lab1, rab1, lb2, rb2, lab2, rab2, (float*)result);
+				return ((bip<float>*)file)->pa_to_pa((float*)result, lb1, rb1, lab1, rab1, lb2, rb2, lab2, rab2, mask);
 			else if(header.data_type == envi_header::float64)
-				return ((bip<double>*)file)->pa_to_pa(lb1, rb1, lab1, rab1, lb2, rb2, lab2, rab2, (double*)result);
+				return ((bip<double>*)file)->pa_to_pa((double*)result, lb1, rb1, lab1, rab1, lb2, rb2, lab2, rab2, mask);
 			else
 				std::cout<<"ERROR: unidentified data type"<<std::endl;
 		}
@@ -817,30 +861,30 @@ public:
 	/// @param lab is the label for the start of the peak
 	/// @param rab is the label for the end of the peak
 	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
-	bool cpoint(double lb1, double rb1, double lab1, double rab1, void* result){
+	bool centroid(void* result, double lb1, double rb1, double lab1, double rab1, unsigned char* mask = NULL){
 		if(header.interleave == envi_header::BSQ){		//if the infile is bsq file
 			if(header.data_type ==envi_header::float32)
-				return ((bsq<float>*)file)->cpoint(lb1, rb1, lab1, rab1, (float*)result);
+				return ((bsq<float>*)file)->centroid((float*)result, lb1, rb1, lab1, rab1, mask);
 			else if(header.data_type == envi_header::float64)
-				return ((bsq<double>*)file)->cpoint(lb1, rb1, lab1, rab1, (double*)result);
+				return ((bsq<double>*)file)->centroid((double*)result, lb1, rb1, lab1, rab1, mask);
 			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)->cpoint(lb1, rb1, lab1, rab1, (float*)result);
+				return ((bil<float>*)file)->centroid((float*)result, lb1, rb1, lab1, rab1, mask);
 			else if(header.data_type == envi_header::float64)
-				return ((bil<double>*)file)->cpoint(lb1, rb1, lab1, rab1, (double*)result);
+				return ((bil<double>*)file)->centroid((double*)result, lb1, rb1, lab1, rab1, mask);
 			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)->cpoint(lb1, rb1, lab1, rab1,(float*)result);
+				return ((bip<float>*)file)->centroid((float*)result, lb1, rb1, lab1, rab1, mask);
 			else if(header.data_type == envi_header::float64)
-				return ((bip<double>*)file)->cpoint(lb1, rb1, lab1, rab1, (double*)result);
+				return ((bip<double>*)file)->centroid((double*)result, lb1, rb1, lab1, rab1, mask);
 			else
 				std::cout<<"ERROR: unidentified data type"<<std::endl;
 		}
@@ -6,6 +6,12 @@
 #include <cstring>
 #include <utility>
+#ifdef _WIN32
+	#include <float.h>
+#else
+	#include<cmath>
+#endif
+
 namespace stim{
 /**
@@ -113,6 +119,46 @@ protected:
 		return c[2] * R[0] * R[1] + c[1] * R[0] + c[0];		//calculate and return the index (trust me this works)
 	}
+
+	/// Returns the 3D coordinates of a specified index
+	void xyb(size_t idx, size_t& x, size_t& y, size_t& b){
+
+		size_t c[3];
+
+		c[2] = idx / (R[0] * R[1]);
+		c[1] = (idx - c[2] * R[0] * R[1]) / R[0];
+		c[0] = idx - c[2] * R[0] * R[1] - c[1] * R[0];
+
+		x = c[O[0]];
+		y = c[O[1]];
+		b = c[O[2]];
+	}
+
+public:
+			/// Get a mask that has all pixels with inf or NaN values masked out (false)
+	void mask_finite(unsigned char* mask, bool PROGRESS = false){
+		size_t XY = X() * Y();
+		memset(mask, 255, XY * sizeof(unsigned char));				//initialize the mask to zero (0)
+		T* page = (T*)malloc(R[0] * R[1] * sizeof(T));		//allocate space for a page of data
+		
+		for(size_t p = 0; p < R[2]; p++){					//for each page
+			binary<T>::read_page(page, p);					//read a page
+			for(size_t i = 0; i < R[0] * R[1]; i++){		//for each pixel in that page
+				
+#ifdef _WIN32
+				if(!_finite(page[i])){						//if the value at index i is finite
+#else
+				if(!std::isfinite(page[i])){					//C++11 implementation
+#endif
+					size_t x, y, b;
+					xyb(p * R[0] * R[1] + i, x, y, b);							//find the 3D coordinates of the value
+					mask[ y * X() + x ] = 0;				//mask the pixel (it's not bad)
+				}
+			}
+			if(PROGRESS) progress = (double)(p + 1) / (double)R[2] * 100;
+		}
+	}
+
 };
 }		//end namespace STIM
@@ -97,7 +97,7 @@ namespace stim{
 		}
         //return the value of a floating point option
-		float as_float(size_t n = 0)
+		double as_float(size_t n = 0)
 		{
             if(!flag)
             {