added progress bars to making masks, fixed overflow error in baseline correction

David Mayerich
1 parent c1078e19
Showing 4 changed files with 88 additions and 78 deletions Show diff stats
stim/envi/bil.h
stim/envi/bip.h
stim/envi/bsq.h
stim/envi/envi.h
@@ -230,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 PROGRESS = false){
+	bool baseline(std::string outname, std::vector<double> wls, unsigned char* mask = NULL, bool PROGRESS = false){
 		unsigned N = wls.size();			//get the number of baseline points
@@ -345,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 PROGRESS = false)
+	bool normalize(std::string outname, double w, unsigned char* mask = NULL, bool PROGRESS = false)
 	{
 		unsigned int B = Z();		//calculate the number of bands
 		unsigned int ZX = Z() * X();
@@ -371,7 +371,7 @@ public:
 			{
 				for(unsigned m = 0; m < X(); m++)
 				{
-					if( b[m + j * X()] < t )
+					if( mask != NULL && !mask[m + j * X()] )
 						c[m + i * X()] = (T)0.0;
 					else
 						c[m + i * X()] = c[m + i * X()] / b[m + j * X()];
@@ -792,16 +792,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 build_mask(double mask_band, double threshold, unsigned char* p, bool PROGRESS = false){
 		T* temp = (T*)malloc(X() * Y() * sizeof(T));		//allocate memory for the certain band
-		band(temp, mask_band);
+		band(temp, mask_band, PROGRESS);
 		for (unsigned int i = 0; i < X() * Y(); i++) {
-				if (temp[i] < threshold)
-					p[i] = 0;
-				else
-					p[i] = 255;
+			if (temp[i] < threshold)
+				p[i] = 0;
+			else
+				p[i] = 255;
 		}
 		free(temp);
@@ -31,13 +31,13 @@ protected:
 	std::vector<double> w; //band wavelength
 	unsigned int offset;		//header offset
-	unsigned long X(){
+	unsigned long long X(){
 		return R[1];
 	}
-	unsigned long Y(){
+	unsigned long long Y(){
 		return R[2];
 	}
-	unsigned long Z(){
+	unsigned long long Z(){
 		return R[0];
 	}
@@ -318,10 +318,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> base_pts, bool PROGRESS = false){
+	bool baseline(std::string outname, std::vector<double> base_pts, unsigned char* mask = NULL, 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		
 		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
@@ -334,30 +333,36 @@ public:
 		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
-				}
-				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];
-				}
-				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
+			if(mask != NULL && !mask[n]){						//if the pixel isn't masked
+				memset(sbc, 0, sizeof(T) * B);					//set the baseline corrected spectrum to zero
+			}
+			else{
+				
+				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
+					}
+					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];
+					}
+					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
 				}
-				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
@@ -378,7 +383,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 PROGRESS = false)
+	bool normalize(std::string outname, double w, unsigned char* mask = NULL, 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		
@@ -390,7 +395,7 @@ public:
 		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
+			if(mask != NULL && !mask[n])						//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
@@ -784,16 +789,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 build_mask(double mask_band, double threshold, unsigned char* p, bool PROGRESS = false){
 		T* temp = (T*)malloc(X() * Y() * sizeof(T));		//allocate memory for the certain band
-		band(temp, mask_band);
+		band(temp, mask_band, PROGRESS);
 		for (unsigned int i = 0; i < X() * Y();i++) {
-				if (temp[i] < threshold)
-					p[i] = 0;
-				else
-					p[i] = 255;
+			if (temp[i] < threshold)
+				p[i] = 0;
+			else
+				p[i] = 255;
 		}
 		free(temp);
@@ -172,9 +172,8 @@ 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 PROGRESS = false )
+	bool baseline(std::string outname, std::vector<double> wls, unsigned char* mask = NULL, 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
@@ -257,8 +256,12 @@ public:
 			//perform the baseline correction
 			for(unsigned i=0; i < XY; i++){
-				double r = (double) (ci - ai) / (double) (bi - ai);
-				c[i] =(T) ( c[i] - (b[i] - a[i]) * r - a[i] );
+				if(mask != NULL && !mask[i])								//if the pixel is excluded by a mask
+					c[i] = 0;												//set the value to zero
+				else{
+					double r = (double) (ci - ai) / (double) (bi - ai);
+					c[i] =(T) ( c[i] - (b[i] - a[i]) * r - a[i] );
+				}
 			}
 			target.write(reinterpret_cast<const char*>(c), S);   //write the corrected data into destination
@@ -282,7 +285,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 PROGRESS = false)
+	bool normalize(std::string outname, double w, unsigned char* mask = NULL, 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
@@ -304,7 +307,7 @@ public:
 			band_index(c, j);                     //get the current band into memory
 			for(unsigned i = 0; i < XY; i++)
 			{
-				if(b[i] < t)
+				if(mask != NULL && !mask[i])
 					c[i] = (T)0.0;
 				else
 					c[i] = c[i] / b[i];
@@ -698,16 +701,18 @@ 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 build_mask(double mask_band, double threshold, unsigned char* p = NULL, bool PROGRESS = false){
 		T* temp = (T*)malloc(X() * Y() * sizeof(T));		//allocate memory for the certain band
 		band(temp, mask_band);
 		for (unsigned int i = 0; i < X() * Y(); i++) {
-				if (temp[i] < threshold)
-					p[i] = 0;
-				else
-					p[i] = 255;
+			if (temp[i] < threshold)
+				p[i] = 0;
+			else
+				p[i] = 255;
+
+			if(PROGRESS) progress = (double) (i+1) / (X() * Y());
 		}
 		free(temp);
@@ -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 PROGRESS = false){
+	bool normalize(std::string outfile, double band, unsigned char* mask = NULL, 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, PROGRESS);
+				return ((bsq<float>*)file)->normalize(outfile, band, mask, PROGRESS);
 			else if(header.data_type == envi_header::float64)
-				return ((bsq<double>*)file)->normalize(outfile,band, threshold, PROGRESS);
+				return ((bsq<double>*)file)->normalize(outfile,band, 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)
-				return ((bil<float>*)file)->normalize(outfile, band, threshold, PROGRESS);
+				return ((bil<float>*)file)->normalize(outfile, band, mask, PROGRESS);
 			else if(header.data_type == envi_header::float64)
-				return ((bil<double>*)file)->normalize(outfile,band, threshold, PROGRESS);
+				return ((bil<double>*)file)->normalize(outfile,band, 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)
-				return ((bip<float>*)file)->normalize(outfile, band, threshold, PROGRESS);
+				return ((bip<float>*)file)->normalize(outfile, band, mask, PROGRESS);
 			else if(header.data_type == envi_header::float64)
-				return ((bip<double>*)file)->normalize(outfile,band, threshold, PROGRESS);
+				return ((bip<double>*)file)->normalize(outfile,band, mask, 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 PROGRESS){
-
+	bool baseline(std::string outfile, std::vector<double> w, unsigned char* mask = NULL, bool PROGRESS = false){
+		header.save(outfile + ".hdr");
 		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, PROGRESS);
+				return ((bsq<float>*)file)->baseline(outfile, w, mask, PROGRESS);
 			else if(header.data_type == envi_header::float64)
-				return ((bsq<double>*)file)->baseline(outfile,w, PROGRESS);
+				return ((bsq<double>*)file)->baseline(outfile,w, mask, 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, PROGRESS);
+				return ((bil<float>*)file)->baseline(outfile, w, mask, PROGRESS);
 			else if(header.data_type == envi_header::float64)
-				return ((bil<double>*)file)->baseline(outfile, w, PROGRESS);
+				return ((bil<double>*)file)->baseline(outfile, w, mask, 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, PROGRESS);
+				return ((bip<float>*)file)->baseline(outfile, w, mask, PROGRESS);
 			else if(header.data_type == envi_header::float64)
-				return ((bip<double>*)file)->baseline(outfile, w, PROGRESS);
+				return ((bip<double>*)file)->baseline(outfile, w, mask, PROGRESS);
 			else{
 				std::cout<<"ERROR: unidentified data type"<<std::endl;
 				exit(1);
@@ -469,31 +469,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 build_mask(double mask_band, double threshold, unsigned char* p = NULL, 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);
+				return ((bsq<float>*)file)->build_mask(mask_band, threshold, p, PROGRESS);
 			else if(header.data_type == envi_header::float64)
-				return ((bsq<double>*)file)->build_mask(mask_band, threshold, p);
+				return ((bsq<double>*)file)->build_mask(mask_band, threshold, p, 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);
+				return ((bil<float>*)file)->build_mask(mask_band, threshold, p, PROGRESS);
 			else if(header.data_type == envi_header::float64)
-				return ((bil<double>*)file)->build_mask(mask_band, threshold, p);
+				return ((bil<double>*)file)->build_mask(mask_band, threshold, p, 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);
+				return ((bip<float>*)file)->build_mask(mask_band, threshold, p, PROGRESS);
 			else if(header.data_type == envi_header::float64)
-				return ((bip<double>*)file)->build_mask(mask_band, threshold, p);
+				return ((bip<double>*)file)->build_mask(mask_band, threshold, p, PROGRESS);
 			else
 				std::cout<<"ERROR: unidentified data type"<<std::endl;
 		}