codebase / stimlib

  #ifndef STIM_BIL_H
  #define STIM_BIL_H
  

  #include "../envi/envi_header.h"

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

  #include "../math/fd_coefficients.h"

  #include <stim/cuda/cudatools/error.h>

  #include <cstring>
  #include <utility>

  #include <deque>

  namespace stim{

  /**
  	The BIL class represents a 3-dimensional binary file stored using band interleaved by line (BIL) image encoding. The binary file is stored
  	such that X-Z "frames" are stored sequentially to form an image stack along the y-axis. When accessing the data sequentially on disk,
  	the dimensions read, from fastest to slowest, are X, Z, Y.
  
  	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 bil: public hsi<T> {

  
  protected:

  	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>::R;

  	bil(){ hsi<T>::init_bil(); }
  

  	/// 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 optional 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){

  		w = wavelengths;

  		return open(filename, vec<unsigned long long>(X, B, Y), 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, bool PROGRESS = false){

  		//return binary<T>::read_plane_1(p, page);

  		if(PROGRESS) progress = 0;

  		unsigned long long L = X() * sizeof(T);		//caculate the number of bytes in a sample line
  		unsigned long long jump = X() * (Z() - 1) * sizeof(T);

  		if (page >= Z()){										//make sure the bank number is right

  			std::cout<<"ERROR: page out of range"<<std::endl;
  			return false;
  		}
  

  		file.seekg(X() * page * sizeof(T), std::ios::beg);

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

  		{

  			file.read((char *)(p + i * X()), L);

  			file.seekg( jump, std::ios::cur);

  			if(PROGRESS) progress = (double)(i+1) / Y() * 100;

  		}
  
  		return true;
  	}
  

  	/// 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 there are no wavelengths in the BSQ file
  		if(w.size() == 0)

  			return band_index(p, (unsigned long long)wavelength);

  		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 of a band

  		unsigned long long page=0;                      //bands around the wavelength

  
  		//get the bands numbers around the wavelength
  
  		//if wavelength is smaller than the first one in header file

  		if ( w[page] > wavelength ){

  			band_index(p, page, PROGRESS);

  			return true;
  		}
  

  		while( w[page] < wavelength )

  		{
  			page++;
  			//if wavelength is larger than the last wavelength in header file

  			if (page == Z()) {

  				band_index(p, Z()-1, PROGRESS);

  				return true;
  			}
  		}

  		if ( wavelength < w[page] )

  		{

  			T * p1;
  			T * p2;

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

  			band_index(p2, page, PROGRESS);

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

  		}
  		else                           //if the wavelength is equal to a wavelength in header file
  		{

  			band_index(p, page, PROGRESS);

  		}
  
  		return true;
  	}
  

  	/// Retrieves a band of x values from a given xz plane.
  
  	/// @param p is a pointer to pre-allocated memory at least Z * sizeof(T) in size
  	/// @param c is a pointer to an existing XZ plane (size X*Z*sizeof(T))
  	/// @param wavelength is the wavelength to retrieve
  	bool read_x_from_xz(T* p, T* c, double wavelength)

  	{

  		unsigned long long B = Z();
  		unsigned long long L = X() * sizeof(T);

  		unsigned long long page=0;                      //samples around the wavelength

  		T * p1;
  		T * p2;
  
  		//get the bands numbers around the wavelength
  
  		//if wavelength is smaller than the first one in header file

  		if ( w[page] > wavelength ){

  			memcpy(p, c, L);

  			return true;
  		}
  

  		while( w[page] < wavelength )

  		{
  			page++;
  			//if wavelength is larger than the last wavelength in header file
  			if (page == B) {

  				memcpy(p, c + (B - 1) * X(), L);

  				return true;
  			}
  		}

  		if ( wavelength < w[page] )

  		{
  			p1=(T*)malloc( L );                     //memory allocation
  			p2=(T*)malloc( L );
  

  			memcpy(p1, c + (page - 1) * X(), L);
  			memcpy(p2, c + page * X(), L);

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

  				double r = (double) (wavelength - w[page-1]) / (double) (w[page] - w[page-1]);

  				p[i] = (T)(((double)p2[i] - (double)p1[i]) * r + (double)p1[i]);

  			}
  		}

  		else                           //if the wavelength is equal to a wavelength in header file

  			memcpy(p, c + page * X(), L);

  
  		return true;

  	}
  

  	/// Retrieve a single spectrum (B-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, size_t n, bool PROGRESS = false){
  		size_t y = n / X();
  		size_t x = n - y * X();
  		return binary<T>::read_line_1(p, x, y, PROGRESS);
  	}

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

  		return binary<T>::read_line_1(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){

  
  		//calculate the corresponding x, y

  		unsigned long long x = n % X();
  		unsigned long long y = n / X();

  
  		//get the pixel
  		return spectrum(p, x, y);
  	}

  	//given a Y ,return a XZ slice

  	bool read_plane_xz(T * p, size_t y){

  		return binary<T>::read_plane_2(p, y);

  	}

  	//given a Y, return ZX slice (transposed such that the spectrum is the leading dimension)
  	int read_plane_zx(T* p, size_t y){
  		T* temp = (T*) malloc(X() * Z() * sizeof(T));	//allocate space to store the temporary xz plane
  		binary<T>::read_plane_2(temp, y);					//load the plane from disk
  		size_t z, x;
  		for(z = 0; z < Z(); z++){
  			for(x = 0; x <= z; x++){
  				p[x * Z() + z] = temp[z * X() + x];		//copy to the destination frame
  			}
  		}
  	}
  
  	//load a frame y into a pre-allocated double-precision array
  	int read_plane_xzd(double* f, size_t y){		

  		size_t XB = X() * Z();

  		T* temp = (T*) malloc(XB * sizeof(T));			//create a temporary location to store the plane at current precision
  		if(!read_plane_y(temp, y)) return 1;			//read the plane in its native format, if it fails return a 1
  		for(size_t i = 0; i < XB; i++) f[i] = temp[i];	//convert the plane to a double
  		return 0;
  	}
  
  	//given a Y, return ZX slice (transposed such that the spectrum is the leading dimension)
  	int read_plane_zxd(double* p, size_t y){
  		T* temp = (T*) malloc(X() * Z() * sizeof(T));		//allocate space to store the temporary xz plane
  		binary<T>::read_plane_2(temp, y);					//load the plane from disk
  		size_t z, x;
  		for(z = 0; z < Z(); z++){
  			for(x = 0; x < X(); x++){
  				p[x * Z() + z] = (double)temp[z * X() + x];	//copy to the destination frame
  			}
  		}
  		return 0;
  	}
  

  	/// 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){

  		unsigned long long 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 ZX = Z() * X();		//calculate the number of points in a Y slice
  		unsigned long long L = ZX * sizeof(T);			//calculate the number of bytes of a Y slice
  		unsigned long long B = Z();

  		T* c;			//pointer to the current Y slice
  		c = (T*)malloc(L);  //memory allocation

  		T* a;			//pointer to the two YZ lines surrounding the current YZ line
  		T* b;

  		a = (T*)malloc(X() * sizeof(T));
  		b = (T*)malloc(X() * sizeof(T));

  
  
  		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;

  
  		if (a == NULL || b == NULL || c == NULL){
  			std::cout<<"ERROR: error allocating memory";
  			exit(1);
  		}
  	//	loop start	correct every y slice

  		for (unsigned long long k =0; k < Y(); k++)

  		{
  			//get the current y slice

  			read_plane_xz(c, k);

  
  			//initialize lownum, highnum, low, high

  			ai = w[0];

  			control=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, X() * sizeof(T) );

  			}
  			//else get the low band
  			else{
  				control++;

  				read_x_from_xz(a, c, ai);

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

  			read_x_from_xz(b, c, bi);

  			//correct every YZ line

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

  						read_x_from_xz(b, c, 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, X() * sizeof(T) );	//clear the high band

  						ai = bi;

  						bi = w[B - 1];

  					}
  				}

  				ci = w[cii];

  				unsigned long long jump = cii * X();

  				//perform the baseline correction

  				for(unsigned i=0; i < X(); i++)

  				{
  					double r = (double) (ci - ai) / (double) (bi - ai);

  					c[i + jump] =(T) ( c[i + jump] - (b[i] - a[i]) * r - a[i] );

  				}

  
  			}//loop for YZ line end
  			target.write(reinterpret_cast<const char*>(c), L);   //write the corrected data into destination

  			if(PROGRESS) progress = (double)(k+1) / Y() * 100;

  		}//loop for Y slice end
  

  		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 ZX = Z() * X();
  		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

  		unsigned long long L = ZX * sizeof(T);

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

  		c = (T*)malloc( L );

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

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

  		{

  			read_plane_xz(c, j);

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

  			{

  				for(unsigned long long m = 0; m < X(); m++)

  				{

  					if( mask != NULL && !mask[m + j * X()] )

  						c[m + i * X()] = (T)0.0;

  					else

  						c[m + i * X()] = c[m + i * X()] / b[m + j * X()];

  				}

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

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

  		}

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

  		return true;
  	}

  	void normalize(std::string outfile, unsigned char* mask = NULL, bool PROGRESS = false){
  		std::cout<<"ERROR: algorithm not implemented"<<std::endl;
  		exit(1);
  		//This code is almost done but has to be debugged
  		/*std::ofstream target(outfile.c_str(), std::ios::binary);	//open the target binary file		
  
  		size_t B = Z();									//calculate the number of bands
  		size_t L = X();									//get the number of items in a line
  		size_t Lb = L * sizeof(T);						//calculate the number of bytes in a line
  		size_t XY = X() * Y();							//calculate the number of pixels in an image
  
  		T* line = (T*) malloc(Lb);						//allocate space for a line
  		T* len = (T*) malloc(Lb);						//allocate space for the lengths in a line
  
  		for(size_t y = 0; y < Y(); y++){						//for each line in the image
  			file.seekg(y * Lb, std::ios::beg);			//move the pointer to the current file to the beginning
  			for(size_t b = 0; b < B; b++){				//for each band in this line
  				file.read((char*)line, Lb);				//read a line of pixels
  				for(size_t l = 0; l < L; l++)			//for each element in the line for this band
  					len[l] += line[l] * line[l];		//add the square of the spectral component
  			}
  			for(size_t l = 0; l < L; l++)				//for each length element in the line
  				len[l] = sqrt(len[l]);					//calculate the square root of the sum of squares
  			
  			file.seekg(y * Lb, std::ios::beg);			//move the pointer to the current file to the beginning
  			for(size_t b = 0; b < B; b++){				//for each band in this line
  				file.read((char*)line, Lb);				//read a line of pixels
  				for(size_t l = 0; l < L; l++)			//for each element in the line for this band
  					line[l] /= len[l];					//divide each element by the length
  			}
  			target.write((char*)line, Lb);				//write the normalized line to the target file
  			if(PROGRESS) progress = (double)(y+1) / (double)Y();
  		}*/
  	}

  	/// 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.

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

  	{

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

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

  		T * p;			//pointer to the current band
  		p = (T*)malloc(S);

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

  		{

  				band_index(p, i);

  				target.write(reinterpret_cast<const char*>(p), S);   //write a band data into target file
  

  				if(PROGRESS) progress = (double)(i+1) / Z() * 100;	//store the progress for the current operation

  		}

  		free(p);
  		target.close();
  		return true;
  	}
  

  	/// Convert the current BIL file to a BIP file with the specified file name.
  
  	/// @param outname is the name of the output BIP file to be saved to disk.

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

  	{

  		unsigned long long S = X() * Z() * sizeof(T);		//calculate the number of bytes in a ZX slice