added sifting and unsifting support for BIP files

David Mayerich
1 parent 2fcab4f0
Showing 2 changed files with 93 additions and 27 deletions Show diff stats
stim/envi/bip.h
stim/envi/envi.h
@@ -837,40 +837,101 @@ public:
  
  
 	/// Saves to disk only those spectra corresponding to mask values != 0
-	bool sift(std::string outfile, unsigned char* p){
-		// Assume X() = X, Y() = Y, Z() = Z.
+	bool sift(std::string outfile, unsigned char* mask){
+
+		//reset the file pointer to the beginning of the file
+		file.seekg(0, std::ios::beg);		
+
+		// open an output stream
 		std::ofstream target(outfile.c_str(), std::ios::binary);
  
-		//for loading pages:
-		unsigned ZX = X() * Z();		//calculate the number of values in an XZ page on disk
-		unsigned L = ZX * sizeof(T);	//calculate the size of the page (in bytes)
-		T * temp = (T*)malloc(L);		//allocate memory for a temporary page
+		//allocate space for a single spectrum
+		unsigned int B = Z();
+		T* spectrum = (T*) malloc(B * sizeof(T));
  
-		//for saving spectra:
-		unsigned LZ = Z() * sizeof(T);			//calculate the size of the spectrum (in bytes)
-		T * tempZ = (T*)malloc(LZ);				//allocate memory for a temporary spectrum
-		spectrum(tempZ, X() - 1, Y() - 1);	//creates a dummy spectrum by taking the last spectrum in the image
+		//calculate the number of pixels in a band
+		unsigned int XY = X() * Y();
  
-		for (unsigned i = 0; i < Y(); i++)			//Select a page by choosing Y coordinate, Y()
-		{
-			read_plane_y(temp, i);							//retrieve an ZX page, store in "temp"
-			for (unsigned j = 0; j < X(); j++)		//Select a pixel by choosing X coordinate in the page, X()
-			{
-				if (p[j * X() + i] != 0)					//if the mask != 0 at that XY pixel
-				{
-					for (unsigned k = 0; k < Z(); k++)		//Select a voxel by choosing Z coordinate at the pixel
-					{
-						tempZ[k] = temp[i*Z() + k];		//Pass the correct spectral value from XZ page into the spectrum to be saved. 
-					}
-					target.write(reinterpret_cast<const char*>(tempZ), LZ);		//write that spectrum to disk. Size is L2.
-				}
-				else
-					continue;
+		//for each pixel
+		unsigned int skip = 0;					//number of spectra to skip
+		for(unsigned int x = 0; x < XY; x++){
+
+			//if the current pixel isn't masked
+			if( mask[x] == 0){
+				//increment the number of skipped pixels
+				skip++;
 			}
+			//if the current pixel is masked
+			else{
+
+				//skip the intermediate pixels
+				file.seekg(skip * B * sizeof(T), std::ios::cur);
+
+				//set the skip value to zero
+				skip = 0;
+
+				//read this pixel into memory
+				file.read((char *)spectrum, B * sizeof(T));
+
+				//write this pixel out
+				target.write((char *)spectrum, B * sizeof(T));
+			}
+
 		}
+
+		//close the output file
 		target.close();
-		free(temp);
+		free(spectrum);
+	}
+
+	bool unsift(std::string outfile, unsigned char* mask, unsigned int samples, unsigned int lines){
+
+		// open an output stream
+		std::ofstream target(outfile.c_str(), std::ios::binary);
+
+		//reset the file pointer to the beginning of the file
+		file.seekg(0, std::ios::beg);		
+
+		//allocate space for a single spectrum
+		unsigned int B = Z();
+		T* spectrum = (T*) malloc(B * sizeof(T));
+
+		//allocate space for a spectrum of zeros
+		T* zeros = (T*) malloc(B * sizeof(T));
+		memset(zeros, 0, B * sizeof(T));
+
+		//calculate the number of pixels in a band
+		unsigned int XY = samples * lines;
+
+		//for each pixel
+		unsigned int skip = 0;					//number of spectra to skip
+		for(unsigned int x = 0; x < XY; x++){
+
+			//if the current pixel isn't masked
+			if( mask[x] == 0){
+
+				//write a bunch of zeros
+				target.write((char *)zeros, B * sizeof(T));
+			}
+			//if the current pixel is masked
+			else{
+
+				//read a pixel into memory
+				file.read((char *)spectrum, B * sizeof(T));
+
+				//write this pixel out
+				target.write((char *)spectrum, B * sizeof(T));
+			}
+
+		}
+
+		//close the output file
+		target.close();
+		free(spectrum);
+
 		return true;
+
+
 	}
  
  
@@ -455,7 +455,12 @@ public:
  
 		else if (header.interleave == envi_header::BIP){		//if the infile is bip file
  
-				std::cout << "ERROR in stim::envi::unsift - BIP files aren't supported yet" << std::endl;
+			if (header.data_type == envi_header::float32)
+				return ((bip<float>*)file)->unsift(outfile, mask, samples, lines);
+			else if (header.data_type == envi_header::float64)
+				return ((bip<double>*)file)->unsift(outfile, mask, samples, lines);
+			else
+				std::cout << "ERROR: unidentified data type" << std::endl;
 		}
  
 		else{