codebase / stimlib

  #ifndef STIM_BSQ_H
  #define STIM_BSQ_H
  
  #include "../envi/envi_header.h"

  #include "../envi/hsi.h"

  #include "../envi/bil.h"

  #include <cstring>
  #include <utility>

  #include <vector>

  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;
  

  	/*unsigned long long X(){

  		return R[0];
  	}

  	unsigned long long Y(){

  		return R[1];
  	}

  	unsigned long long Z(){

  		return R[2];

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

  	using binary<T>::progress;

  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);
  			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

  				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.

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

  		return 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 normalize(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;
  	}

  	/// Convert the current BSQ file to a BIL file with the specified file name.

  
  	/// @param outname is the name of the output BIL file to be saved to disk.

  	bool bil(std::string outname, bool PROGRESS = false)

  	{
  		//simplify image resolution

  		unsigned long long jump = (Y() - 1) * X() * sizeof(T);

  		std::ofstream target(outname.c_str(), std::ios::binary);
  		std::string headername = outname + ".hdr";

  		unsigned long long L = X();

  		T* line = (T*)malloc(sizeof(T) * L);

  		for ( unsigned long long y = 0; y < Y(); y++)									//for each y position

  		{

  			file.seekg(y * X() * sizeof(T), std::ios::beg);					//seek to the beginning of the xz slice

  			for ( unsigned long long z = 0; z < Z(); z++ )							//for each band

  			{

  				file.read((char *)line, sizeof(T) * X());					//read a line
  				target.write((char*)line, sizeof(T) * X());					//write the line to the output file

  				file.seekg(jump, std::ios::cur);							//seek to the next band

  				if(PROGRESS) progress = (double)((y+1) * Z() + z + 1) / (Z() * Y()) * 100;	//update the progress counter

  			}

  		}

  		free(line);

  		target.close();

  		return true;
  	}

  	/// Return a baseline corrected band given two adjacent baseline points and their bands. The result is stored in a pre-allocated array.

  
  	/// @param lb is the label value for the left baseline point
  	/// @param rb is the label value for the right baseline point
  	/// @param lp is a pointer to an array holding the band image for the left baseline point
  	/// @param rp is a pointer to an array holding the band image for the right baseline point
  	/// @param wavelength is the label value for the requested baseline-corrected band
  	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size.

  	bool baseline_band(double lb, double rb, T* lp, T* rp, double wavelength, T* result){
  

  		unsigned long long XY = X() * Y();

  		band(result, wavelength);		//get band
  
  		//perform the baseline correction
  		double r = (double) (wavelength - lb) / (double) (rb - lb);

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

  			result[i] =(T) (result[i] - (rp[i] - lp[i]) * r - lp[i] );
  		}
  		return true;
  	}

  	/// Return a baseline corrected band given two adjacent baseline points. The result is stored in a pre-allocated array.

  
  	/// @param lb is the label value for the left baseline point
  	/// @param rb is the label value for the right baseline point
  	/// @param bandwavelength is the label value for the desired baseline-corrected band
  	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size.

  	bool height(double lb, double rb, double bandwavelength, T* result){

  
  		T* lp;
  		T* rp;

  		unsigned long long XY = X() * Y();
  		unsigned long long S = XY * sizeof(T);

  		lp = (T*) malloc(S);			//memory allocation
  		rp = (T*) malloc(S);

  		band(lp, lb);

  		band(rp, rb);

  
  		baseline_band(lb, rb, lp, rp, bandwavelength, result);
  

  		free(lp);
  		free(rp);

  		return true;
  	}
  

  	/// Calculate the area under the spectrum between two specified points and stores the result in a pre-allocated array.

  
  	/// @param lb is the label value for the left baseline point
  	/// @param rb is the label value for the right baseline point
  	/// @param lab is the label value for the left bound (start of the integration)
  	/// @param rab is the label value for the right bound (end of the integration)
  	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size

  	bool area(double lb, double rb, double lab, double rab, T* result){

  
  		T* lp;	//left band pointer
  		T* rp;	//right band pointer
  		T* cur;		//current band 1
  		T* cur2;	//current band 2

  		unsigned long long XY = X() * Y();
  		unsigned long long S = XY * sizeof(T);

  
  		lp = (T*) malloc(S);			//memory allocation
  		rp = (T*) malloc(S);
  		cur = (T*) malloc(S);
  		cur2 = (T*) malloc(S);

  		//find the wavelength position in the whole band

  		unsigned long long n = w.size();
  		unsigned long long ai = 0;		//left bound position
  		unsigned long long bi = n - 1;		//right bound position

  
  
  
  		//to make sure the left and the right bound are in the bandwidth
  		if (lb < w[0] || rb < w[0] || lb > w[n-1] || rb >w[n-1]){
  			std::cout<<"ERROR: left bound or right bound out of bandwidth"<<std::endl;
  			exit(1);
  		}

  		//to make sure right bound is bigger than left bound

  		else if(lb > rb){
  			std::cout<<"ERROR: right bound should be bigger than left bound"<<std::endl;
  			exit(1);
  		}
  

  		//find the indices of the left and right baseline points
  		while (lab >= w[ai]){ 

  			ai++;
  		}

  		while (rab <= w[bi]){

  			bi--;
  		}
  

  		band(lp, lb);						//get the band images for the left and right baseline points

  		band(rp, rb);
  

  		// calculate the average value of the indexed region
  		memset(result, 0, S);								//initialize the integral to zero (0)
  
  		//integrate the region between the specified bands and the closest indexed band
  		//		this integrates the "tails" of the spectrum that lie outside the main indexed region
  		baseline_band(lb, rb, lp, rp, rab, cur2);		//calculate the image for the right-most band in the integral
  		baseline_band(lb, rb, lp, rp, w[bi], cur);		//calculate the image for the right-most indexed band

  		for(unsigned long long j = 0; j < XY; j++){
  			result[j] += (T)((rab - w[bi]) * ((double)cur[j] + (double)cur2[j]) / 2.0);

  		}

  		baseline_band(lb, rb, lp, rp, lab, cur2);		//beginnning part
  		baseline_band(lb, rb, lp, rp, w[ai], cur);

  		for(unsigned long long j = 0; j < XY; j++){
  			result[j] += (T)((w[ai] - lab) * ((double)cur[j] + (double)cur2[j]) / 2.0);

  		}

  		//integrate the main indexed region

  		ai++;

  		for(unsigned long long i = ai; i <= bi ;i++)	//for each band in the integral

  		{

  			baseline_band(lb, rb, lp, rp, w[ai], cur2);	//get the baselined band
  			for(unsigned long long j = 0; j < XY; j++){	//for each pixel in that band

  				result[j] += (T)((w[ai] - w[ai-1]) * ((double)cur[j] + (double)cur2[j]) / 2.0);

  			}
  			std::swap(cur,cur2);		//swap the band pointers
  		}

  		free(lp);
  		free(rp);
  		free(cur);
  		free(cur2);