added band trimming to ENVI objects

David Mayerich
1 parent 6a46c8ff
Showing 5 changed files with 216 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/envi_header.h
@@ -1098,6 +1098,35 @@ public:
 		return true;
 	}
  
+	/// Remove a list of bands from the ENVI file
+
+	/// @param outfile is the file name for the output hyperspectral image (with trimmed bands)
+	/// @param b is an array of bands to be eliminated
+	void trim(std::string outfile, std::vector<size_t> band_array, bool PROGRESS = false){
+		std::ofstream out(outfile.c_str(), std::ios::binary);	//open the output file for writing
+		file.seekg(0, std::ios::beg);							//move to the beginning of the input file
+
+		size_t Xb = X() * sizeof(T);							//calculate the number of bytes in a line
+		T* line = (T*)malloc(Xb);								//allocate space for a line
+
+		size_t i;												//create an index into the band array
+		for(size_t y = 0; y < Y(); y++){						//for each Y plane
+			i = 0;
+			for(size_t b = 0; b < Z(); b++){					//for each band
+				if(b != band_array[i]){							//if this band isn't trimmed
+					file.read((char*)line, Xb);					//read a line
+					out.write((char*)line, Xb);					//write the line
+				}
+				else{
+					file.seekg(Xb, std::ios::cur);				//if this band is trimmed, skip it
+					i++;
+				}
+			}
+			if(PROGRESS) progress = (double)(y+1) / (double)Y() * 100;
+		}
+		free(line);
+	}
+
  
 	/// Close the file.
 	bool close(){
@@ -1179,6 +1179,43 @@ public:
 		return true;
 	}
  
+	/// Remove a list of bands from the ENVI file
+
+	/// @param outfile is the file name for the output hyperspectral image (with trimmed bands)
+	/// @param b is an array of bands to be eliminated
+	void trim(std::string outfile, std::vector<size_t> band_array, bool PROGRESS = false){
+
+		std::ofstream out(outfile.c_str(), std::ios::binary);	//open the output file for writing
+		file.seekg(0, std::ios::beg);							//move to the beginning of the input file
+
+		size_t B = Z();								//calculate the number of elements in a spectrum
+		size_t Bdst = Z() - band_array.size();		//calculate the number of elements in an output spectrum
+		size_t Bb = B * sizeof(T);					//calculate the number of bytes in a spectrum
+		size_t XY = X() * Y();						//calculate the number of pixels in the image
+		T* src = (T*)malloc(Bb);					//allocate space to store an input spectrum
+		T* dst = (T*)malloc(Bdst * sizeof(T));		//allocate space to store an output spectrum
+
+		size_t i;									//index into the band array
+		size_t bdst;								//index into the output array
+		for(size_t xy = 0; xy < XY; xy++){			//for each pixel
+			i = 0;
+			bdst = 0;
+			file.read((char*)src, Bb);				//read a spectrum
+			for(size_t b = 0; b < B; b++){			//for each band
+				if(b != band_array[i]){				//if the band isn't trimmed
+					dst[bdst] = src[b];				//copy the band value to the output spectrum
+					bdst++;
+				}
+				else i++;							//otherwise increment i
+			}
+			out.write((char*)dst, Bdst * sizeof(T));	//write the output spectrum
+			if(PROGRESS) progress = (double)(xy + 1) / (double) XY * 100;
+		}
+		free(src);
+		free(dst);
+	}
+
+
  
 	/// Close the file.
 	bool close(){
@@ -996,6 +996,35 @@ public:
 		return true;
 	}
  
+	/// Remove a list of bands from the ENVI file
+
+	/// @param outfile is the file name for the output hyperspectral image (with trimmed bands)
+	/// @param b is an array of bands to be eliminated
+	void trim(std::string outfile, std::vector<size_t> band_array, bool PROGRESS = false){
+		
+		std::ofstream out(outfile.c_str(), std::ios::binary);	//open the output file for writing
+		file.seekg(0, std::ios::beg);							//move to the beginning of the input file
+
+		size_t XY = X() * Y();						//calculate the number of elements in a band
+		size_t XYb = XY * sizeof(T);				//calculate the number of bytes in a band
+		T* temp = (T*)malloc(XYb);					//allocate space to store a band
+
+		size_t i = 0;								//store the first index into the band array
+
+		for(size_t b = 0; b < Z(); b++){			//for each band
+			if(b != band_array[i]){					//if this band is not trimmed
+				file.read((char*)temp, XYb);		//read the band
+				out.write((char*)temp, XYb);		//output the band
+			}
+			else{
+				file.seekg(XYb, std::ios::cur);		//otherwise, skip the band
+				i++;
+			}
+			if(PROGRESS) progress = (double)(b+1) / (double) Z() * 100;
+		}
+		free(temp);									//free the scratch space for the band
+	}
+
  
 	/// Close the file.
 	bool close(){
@@ -1235,6 +1235,64 @@ public:
 		return false;
 	}
  
+	/// Remove a list of bands from the ENVI file
+
+	/// @param outfile is the file name for the output hyperspectral image (with trimmed bands)
+	/// @param b is an array of bands to be eliminated
+	void trim(std::string outfile, std::vector<size_t> trimmed, bool PROGRESS = false){
+		
+		envi_header h = header;
+		h.bands = header.bands - trimmed.size();			//calculate the new number of bands
+		if(header.wavelength.size() != 0)
+			h.wavelength.resize(h.bands);
+		if(header.band_names.size() != 0)
+			h.band_names.resize(h.bands);
+		size_t it = 0;									//allocate an index into the trimmed bands array
+		size_t i = 0;
+		for(size_t b = 0; b < header.bands; b++){		//for each band
+			if(b != trimmed[it]){
+				if(h.wavelength.size()) h.wavelength[i] = header.wavelength[b];
+				if(h.band_names.size()) h.band_names[i] = header.band_names[i];
+				i++;
+			}
+			else it++;
+		}
+		h.save(outfile + ".hdr");
+
+		if (header.interleave == envi_header::BSQ){
+			if (header.data_type == envi_header::float32)
+				return ((bsq<float>*)file)->trim(outfile, trimmed, PROGRESS);
+			else if (header.data_type == envi_header::float64)
+				return ((bsq<double>*)file)->trim(outfile, trimmed, 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)->trim(outfile, trimmed, PROGRESS);
+			else if (header.data_type == envi_header::float64)
+				return ((bil<double>*)file)->trim(outfile, trimmed, PROGRESS);
+			else{
+				std::cout << "ERROR: unidentified data type" << std::endl;
+				exit(1);
+			}
+		}
+		else if (header.interleave == envi_header::BIP){
+			if (header.data_type == envi_header::float32)
+				return ((bip<float>*)file)->trim(outfile, trimmed, PROGRESS);
+			else if (header.data_type == envi_header::float64)
+				return ((bip<double>*)file)->trim(outfile, trimmed, PROGRESS);
+			else{
+				std::cout << "ERROR: unidentified data type" << std::endl;
+				exit(1);
+			}
+		}
+		
+
+	}
+
 };
  
 }	//end namespace rts
@@ -6,6 +6,7 @@
 #include <fstream>
 #include <sstream>
 #include <vector>
+#include <algorithm>
 #include <stdlib.h>
  
 //information from an ENVI header file
@@ -431,6 +432,68 @@ struct envi_header
 	std::string interleave_str(){
 		return interleave_str(interleave);
 	}
+
+	/// Convert a wavelength to a band index (or a pair of surrounding band indices)
+	std::vector<size_t> band_index(double w){
+		std::vector<size_t> idx;										//create an empty array of indices
+		if(w < wavelength[0] || w > wavelength[bands-1]) return idx;	//if the wavelength range is outside of the file, return an empty array
+
+		for(size_t b = 0; b < bands; b++){							//for each band in the wavelength vector
+			if(wavelength[b] == w){									//if an exact match is found
+				idx.push_back(b);										//add the band to the array and return it
+				return idx;
+			}
+			if(wavelength[b] >= w){ 								//if the current wavelength exceeds w
+				idx.resize(2);
+				idx[0] = b-1;										//push both the previous and current band index
+				idx[1] = b;
+				return idx;											//return the pair of indices
+			}
+		}
+		return idx;
+	}
+
+	/// Convert a wavelength range to a list of bands
+	std::vector<size_t> band_indices(double w0, double w1){
+
+		size_t idx0 = 0;
+		size_t idx1 = 0;									//declare the interval indices for the band range
+
+		//get the indices for the first wavelength
+		std::vector<size_t> r;
+
+		if(w0 > wavelength[bands-1] || w1 < wavelength[0]){
+			std::cout<<"ERROR in envi_header::band_indices - wavelengths are completely out of range: ["<<w0<<", "<<w1<<"]"<<std::endl;
+			exit(1);
+		}
+		if(w0 < wavelength[0])								//if the lower wavelength is outside the file range
+			idx0 = 0;										//just set it to zero
+		else{
+			r = band_index(w0);
+			if(r.size() == 1)								//if there is an exact match, set the first interval index
+				idx0 = r[0];
+			else
+				idx0 = r[1];								//otherwise save the highest band index (this is the first band with wavelength > w0)
+		}
+		//get the indices for the second wavelength
+		if(w1 > wavelength[bands-1])
+			idx1 = bands-1;
+		else{
+			r = band_index(w1);
+			if(r.size() == 0)
+				idx1 = bands - 1;
+			else
+				idx1 = r[0];									//take the lowest band index (this is the last band < w1)
+		}
+
+		size_t n = idx1 - idx0 + 1;							//calculate the number of bands in this interval
+		r.resize(n);										//resize the band index array
+		for(size_t b = 0; b < n; b++){						//for each band in the interval
+			r[b] = idx0 + b;								//	insert the band in the array
+		}
+
+		return r;											//return the index array
+	}
 };		//end EnviHeader
 }
 #endif