added a --select feature that creates a new ENVI file from individual bands

David Mayerich
1 parent 66304f4f
Showing 5 changed files with 109 additions and 0 deletions Show diff stats
stim/envi/bil.h
stim/envi/bip.h
stim/envi/bsq.h
stim/envi/envi.h
stim/envi/hsi.h
@@ -462,6 +462,11 @@ public:
 		}*/
 	}
  
+	bool select(std::string outfile, std::vector<double> bandlist, unsigned char* mask = NULL, bool PROGRESS = NULL) {
+		std::cout << "ERROR: select() not implemented for BIL" << std::endl;
+		exit(1);
+	}
+
 	/// 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.
@@ -392,6 +392,50 @@ public:
 		}
 	}
  
+	/// This function loads a specified set of bands and saves them into a new output file
+	bool select(std::string outfile, std::vector<double> bandlist, unsigned char* mask = NULL, bool PROGRESS = false) {
+		std::ofstream target(outfile.c_str(), std::ios::binary);	//open the target binary file
+		if (!target) {
+			std::cout << "ERROR opening output file: " << outfile << std::endl;
+			return false;
+		}
+		file.seekg(0, std::ios::beg);								//move the pointer to the current file to the beginning
+
+		size_t B = Z();												//number of spectral components
+		size_t XY = X() * Y();										//calculate the number of pixels
+		size_t Bout = bandlist.size();
+		size_t in_bytes = B * sizeof(T);							//number of bytes in a spectrum
+		size_t out_bytes = Bout * sizeof(T);						//number of bytes in an output spectrum
+
+		T* in = (T*)malloc(in_bytes);								//allocate space for the input spectrum
+		T* out = (T*)malloc(out_bytes);								//allocate space for the output spectrum
+
+		double wc;													//register to store the desired wavelength
+		double w0, w1;												//registers to store the wavelengths surrounding the given band
+		size_t b0, b1;												//indices of the bands surrounding the specified wavelength
+		for (size_t xy = 0; xy < XY; xy++) {						//for each pixel
+			//memset(out, 0, out_bytes);								//set the spectrum to zero
+			if (mask == NULL || mask[xy]) {							//if the pixel is masked
+				file.read((char*)in, in_bytes);						//read an input spectrum
+				for (size_t b = 0; b < Bout; b++) {					//for each band
+					wc = bandlist[b];									//set the desired wavelength
+					band_bounds(wc, b0, b1);								//get the surrounding bands
+					w0 = w[b0];											//get the wavelength for the lower band
+					w1 = w[b1];											//get the wavelength for the higher band
+					out[b] = lerp(wc, in[b0], w0, in[b1], w1);			//interpolate the spectral values to get the desired output value
+				}
+			}
+			else
+				file.seekg(Bout, std::ios::cur);						//otherwise skip a spectrum
+			target.write((char*)out, out_bytes);							//output the normalized spectrum
+			if (PROGRESS) progress = (double)(xy + 1) / (double)XY * 100;		//update the progress
+		}
+
+		free(in);
+		free(out);
+		return true;
+	}
+
  
 	/// Convert the current BIP file to a BIL file with the specified file name.
  
@@ -380,6 +380,30 @@ public:
  
 	}
  
+	bool select(std::string outfile, std::vector<double> bandlist, unsigned char* mask = NULL, bool PROGRESS = false) {
+		std::ofstream out(outfile.c_str(), std::ios::binary);		//open the output file
+		if (!out) {
+			std::cout << "ERROR opening output file: " << outfile << std::endl;
+			return false;
+		}
+		file.seekg(0, std::ios::beg);							//move the pointer to the current file to the beginning
+
+		size_t B = Z();								//calculate the number of bands
+		size_t XY = X() * Y();						//calculate the number of pixels in a band
+		size_t in_bytes = XY * sizeof(T);				//calculate the size of a band in bytes
+
+		T* in = (T*)malloc(in_bytes);				//allocate space for the band image
+		size_t nb = bandlist.size();				//get the number of bands in the output image
+		for (size_t b = 0; b < nb; b++) {
+			band(in, bandlist[b]);					//get the band associated with the given wavelength
+			out.write((char*)in, in_bytes);		//write the band to the output file
+			if (PROGRESS) progress = (double)(b + 1) / (double)bandlist.size() * 100;
+		}
+		out.close();
+		free(in);
+		return true;
+	}
+
 	size_t readlines(T* dest, size_t start, size_t n){
 		return hsi<T>::read(dest, 0, start, 0, X(), n, Z());
 	}
@@ -562,6 +562,41 @@ public:
 		}
 	}
  
+	bool select(std::string outfile, std::vector<double> bandlist, unsigned char* MASK = NULL, bool PROGRESS = false) {
+		stim::envi_header new_header = header;					//copy all of the data from the current header file to the new one
+		new_header.bands = bandlist.size();						//the number of bands in the new file is equal to the number of bands provided by the user
+		new_header.wavelength = bandlist;						//the wavelength values in the output file are the same as those specified by the user
+		new_header.band_names.empty();							//no band names will be provided in the output file
+		new_header.save(outfile + ".hdr");						//save the output header file
+
+		if (header.interleave == envi_header::BSQ) {		//if the infile is bip file
+			if (header.data_type == envi_header::float32)
+				return ((bsq<float>*)file)->select(outfile, bandlist, MASK, PROGRESS);
+			else if (header.data_type == envi_header::float64)
+				return ((bsq<double>*)file)->select(outfile, bandlist, MASK, PROGRESS);
+			else
+				std::cout << "ERROR: unidentified data type" << std::endl;
+		}
+		else if (header.interleave == envi_header::BIL) {		//if the infile is bip file
+			if (header.data_type == envi_header::float32)
+				return ((bil<float>*)file)->select(outfile, bandlist, MASK, PROGRESS);
+			else if (header.data_type == envi_header::float64)
+				return ((bil<double>*)file)->select(outfile, bandlist, 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)->select(outfile, bandlist, MASK, PROGRESS);
+			else if (header.data_type == envi_header::float64)
+				return ((bip<double>*)file)->select(outfile, bandlist, MASK, PROGRESS);
+			else
+				std::cout << "ERROR: unidentified data type" << std::endl;
+		}
+
+		return false;
+	}
+
 	/// Performs piecewise linear baseline correction of a hyperspectral file/
  
 	/// @param outfile is the file name for the baseline corrected output
@@ -55,6 +55,7 @@ protected:
 		return n;													//return the number of masked pixels
 	}
  
+	//perform linear interpolation between two bands
 	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