codebase / stimlib

  #ifndef STIM_BSQ_H
  #define STIM_BSQ_H
  

  #include <stim/envi/envi_header.h>
  #include <stim/envi/hsi.h>
  #include <stim/envi/bil.h>

  #include <cstring>
  #include <utility>

  #include <vector>

  #include <deque>

  #include <chrono>

  #include <future>

  #include <algorithm>

  namespace stim{

  /**
  	The BIP class represents a 3-dimensional binary file stored using band sequential (BSQ) image encoding. The binary file is stored
  	such that X-Y "frames" are stored sequentially to form an image stack along the z-axis. When accessing the data sequentially on disk,
  	the dimensions read, from fastest to slowest, are X, Y, Z.
  
  	This class is optimized for data streaming, and therefore supports extremely large (terabyte-scale) files. Data is loaded from disk
  	on request. Functions used to access data are written to support efficient reading.
  */

  template <typename T>
  

  class bsq: public hsi<T> {

  protected:

  	//std::vector<double> w;	//band wavelengths
  	unsigned long long offset;

  	using binary<T>::R;
  

  	using hsi<T>::w;				//use the wavelength array in stim::hsi
  	using hsi<T>::nnz;

  	using binary<T>::progress;

  	using hsi<T>::X;
  	using hsi<T>::Y;
  	using hsi<T>::Z;

  public:
  

  	using binary<T>::open;
  	using binary<T>::file;

  	using binary<T>::read_line_2;

  	using binary<T>::read_plane_2;

  	bsq(){ hsi<T>::init_bsq(); }

  	/// Open a data file for reading using the class interface.

  
  	/// @param filename is the name of the binary file on disk
  	/// @param X is the number of samples along dimension 1
  	/// @param Y is the number of samples (lines) along dimension 2
  	/// @param B is the number of samples (bands) along dimension 3
  	/// @param header_offset is the number of bytes (if any) in the binary header

  	/// @param wavelengths is an STL vector of size B specifying a numerical label for each band

  	bool open(std::string filename,
  			  unsigned long long X,
  			  unsigned long long Y,
  			  unsigned long long B,
  			  unsigned long long header_offset,

  			  std::vector<double> wavelengths){

  		//copy the wavelengths to the BSQ file structure
  		w = wavelengths;

  		//copy the wavelengths to the structure
  		offset = header_offset;

  		return open(filename, vec<unsigned long long>(X, Y, B), header_offset);

  	}
  

  	/// Retrieve a single band (based on index) and stores it in pre-allocated memory.

  
  	/// @param p is a pointer to an allocated region of memory at least X * Y * sizeof(T) in size.
  	/// @param page <= B is the integer number of the band to be copied.

  	bool band_index( T * p, unsigned long long page){

  		return read_plane_2(p, page);		//call the binary read_plane function (don't let it update the progress)

  	}
  

  	/// Retrieve a single band (by numerical label) and stores it in pre-allocated memory.

  
  	/// @param p is a pointer to an allocated region of memory at least X * Y * sizeof(T) in size.
  	/// @param wavelength is a floating point value (usually a wavelength in spectral data) used as a label for the band to be copied.

  	bool band( T * p, double wavelength, bool PROGRESS = false){
  		if(PROGRESS) progress = 0;

  		//if there are no wavelengths in the BSQ file

  		if(w.size() == 0){
  			band_index(p, (unsigned long long)wavelength);
  			if(PROGRESS) progress = 100;
  			return true;
  		}

  		unsigned long long XY = X() * Y();	//calculate the number of pixels in a band
  		unsigned long long page = 0;

  		//get the two neighboring bands (above and below 'wavelength')

  
  		//if wavelength is smaller than the first one in header file

  		if ( w[page] > wavelength ){

  			band_index(p, page);

  			if(PROGRESS) progress = 100;

  			return true;
  		}
  

  		while( w[page] < wavelength )

  		{
  			page++;

  			//if wavelength is larger than the last wavelength in the header file
  			//	(the wavelength is out of bounds)

  			if (page == Z()) {
  				band_index(p, Z()-1);		//return the last band

  				if(PROGRESS) progress = 100;

  				return true;
  			}
  		}

  		//when the page counter points to the first band above 'wavelength'
  		if ( wavelength < w[page] ){

  			//do the interpolation

  			T * p1;
  			T * p2;

  			p1=(T*)malloc( XY * sizeof(T));                     //memory allocation
  			p2=(T*)malloc( XY * sizeof(T));
  			band_index(p1, page - 1);
  			band_index(p2, page );

  			for(unsigned long long i=0; i < XY; i++){
  				double r = (wavelength - w[page-1]) / (w[page] - w[page-1]);
  				p[i] = (T)(((double)p2[i] - (double)p1[i]) * r + (double)p1[i]);

  			}

  			free(p1);
  			free(p2);

  		}

  		//if the wavelength is equal to a wavelength in header file
  		else{
  			band_index(p, page);		//return the band

  		}

  		if(PROGRESS) progress = 100;

  		return true;
  	}
  

  	/// Retrieve a single spectrum (Z-axis line) at a given (x, y) location and stores it in pre-allocated memory.

  
  	/// @param p is a pointer to pre-allocated memory at least B * sizeof(T) in size.
  	/// @param x is the x-coordinate (dimension 1) of the spectrum.
  	/// @param y is the y-coordinate (dimension 2) of the spectrum.

  	void spectrum(T* p, size_t n, bool PROGRESS){
  		read_line_2(p, n, PROGRESS);

  	}

  	void spectrum(T * p, unsigned long long x, unsigned long long y, bool PROGRESS = false){
  		read_line_2(p, x, y, PROGRESS);
  	}
  

  	/// Retrieve a single pixel and stores it in pre-allocated memory.

  
  	/// @param p is a pointer to pre-allocated memory at least sizeof(T) in size.
  	/// @param n is an integer index to the pixel using linear array indexing.

  	bool pixel(T * p, unsigned long long n){

  		unsigned long long bandnum = X() * Y();		//calculate numbers in one band

  		if ( n >= bandnum){							//make sure the pixel number is right
  			std::cout<<"ERROR: sample or line out of range"<<std::endl;
  			return false;
  		}
  

  		file.seekg(n * sizeof(T) + binary<T>::header, std::ios::beg);           //point to the certain pixel

  		for (unsigned long long i = 0; i < Z(); i++)

  		{
  			file.read((char *)(p + i), sizeof(T));
  			file.seekg((bandnum - 1) * sizeof(T), std::ios::cur);    //go to the next band
  		}
  

  		return true;

  	}
  

  	/// Perform baseline correction given a list of baseline points and stores the result in a new BSQ file.

  
  	/// @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, unsigned char* mask = NULL, bool PROGRESS = false )

  	{

  		size_t N = wls.size();			//get the number of baseline points

  		std::ofstream target(outname.c_str(), std::ios::binary);	//open the target binary file

  		std::string headername = outname + ".hdr";              //the header file name

  
  		//simplify image resolution

  		unsigned long long B = Z();		//calculate the number of bands
  		unsigned long long XY = X() * Y();	//calculate the number of pixels in a band
  		unsigned long long S = XY * sizeof(T);		//calculate the number of bytes in a band

  
  		double ai, bi;	//stores the two baseline points wavelength surrounding the current band
  		double ci;		//stores the current band's wavelength

  		unsigned long long control=0;

  
  		T * a;					//pointers to the high and low band images
  		T * b;
  		T * c;				//pointer to the current image
  
  		a = (T*)malloc( S );     //memory allocation

  		b = (T*)malloc( S );
  		c = (T*)malloc( S );

  
  		if (a == NULL || b == NULL || c == NULL){
  			std::cout<<"ERROR: error allocating memory";
  			exit(1);
  		}
  
  

  		//initialize lownum, highnum, low, high

  		ai=w[0];

  
  		//if no baseline point is specified at band 0,
  			//set the baseline point at band 0 to 0

  		if(wls[0] != w[0]){

  			bi = wls[control];

  			memset(a, (char)0, S);
  		}
  		//else get the low band
  		else{
  			control += 1;

  			band(a, ai);

  			bi = wls[control];
  		}
  		//get the high band

  		band(b, bi);

  		//correct every band

  		for(unsigned long long cii = 0; cii < B; cii++){

  
  			//update baseline points, if necessary

  			if( w[cii] >= bi && cii != B - 1) {

  				//if the high band is now on the last BL point?
  				if (control != N-1) {
  
  					control++;		//increment the index
  
  					std::swap(a, b);	//swap the baseline band pointers
  
  					ai = bi;
  					bi = wls[control];

  					band(b, bi);

  
  				}
  				//if the last BL point on the last band of the file?

  				else if ( wls[control] < w[B - 1]) {

  
  					std::swap(a, b);	//swap the baseline band pointers
  
  					memset(b, (char)0, S);	//clear the high band
  
  					ai = bi;

  					bi = w[B - 1];

  				}
  			}
  
  			//get the current band
  			band_index(c, cii);

  			ci = w[cii];

  			//perform the baseline correction

  			for(unsigned long long i=0; i < XY; 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

  			if(PROGRESS)progress = (double)(cii+1) / B * 100;

  
  		}

  		//header.save(headername);         //save the new header file

  		free(a);
  		free(b);
  		free(c);
  		target.close();

  		return true;
  	}

  	/// Normalize all spectra based on the value of a single band, storing the result in a new BSQ file.

  
  	/// @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 ratio(std::string outname, double w, unsigned char* mask = NULL, bool PROGRESS = false)

  	{

  		unsigned long long B = Z();		//calculate the number of bands
  		unsigned long long XY = X() * Y();	//calculate the number of pixels in a band
  		unsigned long long S = XY * sizeof(T);		//calculate the number of bytes in a band

  
  		std::ofstream target(outname.c_str(), std::ios::binary);	//open the target binary file
  		std::string headername = outname + ".hdr";              //the header file name
  
  		T * b;					//pointers to the certain wavelength band
  		T * c;				//pointer to the current image
  
  		b = (T*)malloc( S );     //memory allocation

  		c = (T*)malloc( S );

  		band(b, w);             //get the certain band into memory

  		for(unsigned long long j = 0; j < B; j++)

  		{
  			band_index(c, j);                     //get the current band into memory

  			for(unsigned long long i = 0; i < XY; i++)

  			{

  				if(mask != NULL && !mask[i])

  					c[i] = (T)0.0;
  				else
  					c[i] = c[i] / b[i];

  			}
  			target.write(reinterpret_cast<const char*>(c), S);   //write normalized data into destination

  			if(PROGRESS) progress = (double)(j+1) / B * 100;

  		}
  

  		//header.save(headername);         //save the new header file

  		free(b);
  		free(c);
  		target.close();

  		return true;
  	}

  	void normalize(std::string outfile, unsigned char* mask = NULL, bool PROGRESS = false){
  		size_t B = Z();								//calculate the number of bands
  		size_t XY = X() * Y();						//calculate the number of pixels in a band
  		size_t XYb = XY * sizeof(T);				//calculate the size of a band in bytes
  
  		std::ofstream out(outfile.c_str(), std::ios::binary);		//open the output file
  		file.seekg(0, std::ios::beg);							//move the pointer to the current file to the beginning
  
  		T* len = (T*)malloc(XYb);					//allocate space to store the vector length
  		memset(len, 0, XYb);						//initialize the vector length to zero (0)
  
  		T* band = (T*) malloc(XYb);					//allocate space to store a band image
  
  		for(size_t b = 0; b < B; b++){
  			file.read((char*)band, XYb);
  			for(size_t xy = 0; xy < XY; xy++){
  				if(mask == NULL || mask[xy]){
  					len[xy] += pow(band[xy], 2);		//sum the squared value for each pixel value in the band
  				}
  			}
  			if(PROGRESS) progress = (double) (b+1) / (double)B * 50;
  		}
  		for(size_t xy = 0; xy < XY; xy++){				//for each pixel, calculate the square root
  			if(mask == NULL || mask[xy]){
  				len[xy] += pow(band[xy], 2);		//sum the squared value for each pixel value in the band
  			}
  		}
  		file.seekg(0, std::ios::beg);							//move the pointer to the current file to the beginning
  		for(size_t b = 0; b < B; b++){
  			file.read((char*)band, XYb);
  			for(size_t xy = 0; xy < XY; xy++){
  				if(mask == NULL || mask[xy]){
  					band[xy] /= len[xy];						//divide the band by the vector length
  				}
  			}
  			out.write((char*)band, XYb);
  			if(PROGRESS) progress = (double) (b+1) / (double)B * 50 + 50;
  		}

  	}
  

  	size_t readlines(T* dest, size_t start, size_t n){
  		return hsi<T>::read(dest, 0, start, 0, X(), n, Z());
  	}
  
  	size_t writeblock(std::ofstream* f, T* src, size_t n){
  		auto t0 = std::chrono::high_resolution_clock::now();
  		f->write((char*)src, n);
  		auto t1 = std::chrono::high_resolution_clock::now();
  		return std::chrono::duration_cast<std::chrono::milliseconds>(t1-t0).count();

  	}
  

  	/// Convert this BSQ file to a BIL

  	bool bil(std::string outname, bool PROGRESS = false, bool VERBOSE = false, bool OPTIMIZATION = true){

  		const size_t buffers = 4;													//number of buffers required for this algorithm

  		size_t mem_per_batch = binary<T>::buffer_size / buffers;					//calculate the maximum memory available for a batch
  
  		size_t slice_bytes = X() * Z() * sizeof(T);									//number of bytes in an input batch slice (Y-slice in this case)
  		size_t max_slices_per_batch = mem_per_batch / slice_bytes;					//maximum number of slices we can process in one batch given memory constraints

  		std::cout<<"maximum memory available for processing: "<<(double)binary<T>::buffer_size/(double)1000000<<" MB"<<std::endl;
  		std::cout<<"     this supports a batch size of "<<max_slices_per_batch<<" Y-axis slices ("<<X()<<" x "<<Z()<<") = "<<X() * Z() * sizeof(T) * max_slices_per_batch/1000000<<" MB"<<std::endl;

  		if(max_slices_per_batch == 0){														//if there is insufficient memory for a single slice, throw an error

  			std::cout<<"error, insufficient memory for stim::bsq::bil()"<<std::endl;
  			exit(1);
  		}

  		size_t max_batch_bytes = max_slices_per_batch * slice_bytes;				//calculate the amount of memory that will be allocated for all four buffers
  

  		stream_optimizer O(1, max_slices_per_batch);

  		T* src[2];																	//source double-buffer for asynchronous batching
  		src[0] = (T*) malloc(max_batch_bytes);
  		src[1] = (T*) malloc(max_batch_bytes);
  		T* dst[2];																	//destination double-buffer for asynchronous batching
  		dst[0] = (T*) malloc(max_batch_bytes);
  		dst[1] = (T*) malloc(max_batch_bytes);
  
  		size_t N[2];																		//number of slices stored in buffers 0 and 1

  		N[0] = N[1] = std::min<size_t>(Y(), max_slices_per_batch);										//start with the maximum number of slices that can be stored (may be the entire data set)

  
  		std::ofstream target(outname.c_str(), std::ios::binary);					//open an output file for writing
  																		//initialize with buffer 0 (used for double buffering)
  		size_t y_load = 0;
  		size_t y_proc = 0;

  		std::future<size_t> rthread;

  		std::future<std::ostream&> wthread;										//create asynchronous threads for reading and writing

  		std::chrono::high_resolution_clock::time_point t_start, pt_start;						//high-resolution timers
  		std::chrono::high_resolution_clock::time_point t_end, pt_end;

  		size_t t_batch;																//number of milliseconds to process a batch

  		size_t t_total = 0;														//total time for operation
  		size_t pt_total = 0;													//total time spent processing data
  		size_t rt_total = 0;													//total time spent reading data
  		size_t wt_total = 0;

  		size_t dr = 0;

  		
  		rt_total += readlines(src[0], 0, N[0]);					//read the first batch into the 0 source buffer

  		y_load += N[0];											//increment the loaded slice counter
  		int b = 1;												//initialize the double buffer to 0

  		while(y_proc < Y()){													//while there are still slices to be processed

  			t_start = std::chrono::high_resolution_clock::now();					//start the timer for this batch

  			if(y_load < Y()){													//if there are still slices to be loaded, load them

  				//if(y_proc > 0){
  					
  					
  				//}

  				if(y_load + N[b] > Y()) N[b] = Y() - y_load;					//if the next batch would process more than the total slices, adjust the batch size
  				rthread = std::async(std::launch::async, &stim::bsq<T>::readlines, this, src[b], y_load, N[b]);

  				//rt_total += rthread.get();

  				y_load += N[b];													//increment the number of loaded slices

  			}

  			b = !b;																//swap the double-buffer

  			pt_total += binary<T>::permute(dst[b], src[b], X(), N[b], Z(), 0, 2, 1);		//permute the batch to a BIL file

  			wt_total += writeblock(&target, dst[b], N[b] * slice_bytes);			//write the permuted data to the output file

  			y_proc += N[b];														//increment the counter of processed pixels
  			if(PROGRESS) progress = (double)( y_proc + 1 ) / Y() * 100;			//increment the progress counter based on the number of processed pixels

  			if(y_proc < Y()) rt_total += rthread.get();					//if a new batch was set to load, make sure it loads after calculations

  			t_end = std::chrono::high_resolution_clock::now();
  			t_batch = std::chrono::duration_cast<std::chrono::milliseconds>(t_end-t_start).count();
  			t_total += t_batch;

  			if(OPTIMIZATION)
  				N[b] = O.update(N[!b] * slice_bytes, t_batch, binary<T>::data_rate, VERBOSE);					//set the batch size based on optimization
  			//binary<T>::data_rate = dr;

  			//std::cout<<"New N = "<<N[!b]<<" selected with "<<(double)data_rate / 1000000<<" MB/s"<<std::endl;

  		}

  		free(src[0]);															//free buffer resources
  		free(src[1]);
  		free(dst[0]);
  		free(dst[1]);

  		//if(VERBOSE){

  			std::cout<<"total time to execute bsq::bil(): "<<t_total<<" ms"<<std::endl;
  			std::cout<<"     total time spent processing: "<<pt_total<<" ms"<<std::endl;
  			std::cout<<"        total time spent reading: "<<rt_total<<" ms"<<std::endl;
  			std::cout<<"        total time spent writing: "<<wt_total<<" ms"<<std::endl;

  		//}

  		return true;															//return true

  	}
  

  	/// Convert this BSQ file to a BIP

  	bool bip(std::string outname, bool PROGRESS = false, bool VERBOSE = false, bool OPTIMIZATION = true){

  		const size_t buffers = 4;													//number of buffers required for this algorithm

  		size_t mem_per_batch = binary<T>::buffer_size / buffers;					//calculate the maximum memory available for a batch

  		size_t slice_bytes = X() * Z() * sizeof(T);									//number of bytes in an input batch slice (Y-slice in this case)
  		size_t max_slices_per_batch = mem_per_batch / slice_bytes;					//maximum number of slices we can process in one batch given memory constraints