codebase / stimlib

  #ifndef STIM_ENVI_H
  #define STIM_ENVI_H
  
  #include "../envi/envi_header.h"
  #include "../envi/bsq.h"
  #include "../envi/bip.h"
  #include "../envi/bil.h"

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

  #include <stim/parser/filename.h>

  #include <stim/util/filesize.h>

  #include <iostream>

  #include <fstream>

  //#include "../image/image.h"

  namespace stim{

  /** This class implements reading of ENVI hyperspectral files. These files can be stored in multiple orientations
  	(including BSQ, BIP, and BIL) in order to optimize streaming speed depending on applications. Basic ENVI
  	files are stored on disk as a large binary file with a corresponding header. Code for reading and processing
  	ENVI header files is in the envi_header class.
  */

  class envi{

  
  	void* file;		//void pointer to the relevant file reader (bip, bsq, or bil - with appropriate data type)

  	std::string fname;					//file name used for repeated opening and closing

  	//allocate sufficient space for a spectrum based on the data type and number of bands
  	void* alloc_array(size_t len){
  		switch(header.data_type){
  		case envi_header::int8:
  			return malloc(len);
  		case envi_header::int16:
  		case envi_header::uint16:
  			return malloc(2 * len);
  		case envi_header::int32:
  		case envi_header::uint32:
  		case envi_header::float32:
  			return malloc(4 * len);
  		case envi_header::int64:
  		case envi_header::uint64:
  		case envi_header::float64:
  		case envi_header::complex32:
  			return malloc(8 * len);
  		case envi_header::complex64:
  			return malloc(16 * len);
  		default:
  			std::cout<<"ERROR stim::envi data type not recognized for spectral allocation"<<std::endl;
  			exit(1);
  		}
  	}
  
  	//cast a value to DEST type from SRC type
  	template<typename DST, typename SRC>
  	inline void cast(DST* dst, SRC* src){
  		(*dst) = (DST)(*src);
  	}
  
  	//cast an array from type SRC to type DEST
  	template<typename DST, typename SRC>
  	inline void cast(DST* dst, SRC* src, size_t len){
  		for(size_t i = 0; i < len; i++)
  			cast(&dst[i], &src[i]);
  	}

  public:

  	envi_header header;

  	/// Default constructor
  	envi(){
  		file = NULL;				//set the file pointer to NULL
  	}
  

  	envi(std::string filename, std::string headername){
  		header.load(headername);
  
  		fname = filename;					//save the filename
  
  		allocate();
  	}

  	//used to test if the current ENVI file is valid
  	operator bool(){
  		if(file == NULL) return false;
  		return true;
  	}

  	//test to determine if the specified file is an ENVI file
  	static bool is_envi(std::string fname, std::string hname = ""){

  		stim::filename data_file(fname);
  		stim::filename header_file;
  		if(hname == ""){								//if the header isn't provided
  			header_file = data_file;					//assume that it's the same name as the data file, with a .hdr extension
  			header_file = header_file.extension("hdr");
  		}
  		else header_file = hname;						//otherwise load the passed header
  
  		stim::envi_header H;
  		if(H.load(header_file) == false)				//load the header file, if it doesn't load return false
  			return false;
  		size_t targetBytes = H.data_bytes();			//get the number of bytes that SHOULD be in the data file

  		size_t bytes = stim::file_size(fname);

  		if(bytes != targetBytes) return false;			//if the data doesn't match the header, return false
  		return true;									//otherwise everything looks fine

  	}

  	void* malloc_spectrum(){
  		return alloc_array(header.bands);
  	}
  
  	void* malloc_band(){
  		return alloc_array(header.samples * header.lines);
  	}
  

  	void set_buffer_frac(double memfrac = 0.5){

  		if(header.interleave == envi_header::BSQ){		//if the infile is bsq file
  			if(header.data_type ==envi_header::float32)

  				((bsq<float>*)file)->set_buffer_frac(memfrac);

  			else if(header.data_type == envi_header::float64)

  				((bsq<double>*)file)->set_buffer_frac(memfrac);

  			else
  				std::cout<<"ERROR: unidentified data type"<<std::endl;
  		}
  
  		else if(header.interleave == envi_header::BIL){		//if the infile is bil file
  			if(header.data_type ==envi_header::float32)

  				((bil<float>*)file)->set_buffer_frac(memfrac);

  			else if(header.data_type == envi_header::float64)

  				((bil<double>*)file)->set_buffer_frac(memfrac);

  			else
  				std::cout<<"ERROR: unidentified data type"<<std::endl;
  		}
  
  		else if(header.interleave == envi_header::BIP){		//if the infile is bip file
  			if(header.data_type ==envi_header::float32)

  				((bip<float>*)file)->set_buffer_frac(memfrac);

  			else if(header.data_type == envi_header::float64)

  				((bip<double>*)file)->set_buffer_frac(memfrac);
  			else
  				std::cout<<"ERROR: unidentified data type"<<std::endl;
  		}
  
  		else{
  			std::cout<<"ERROR: unidentified file type"<<std::endl;
  			exit(1);
  		}
  	}
  
  	void set_buffer_raw(size_t bytes){
  		if(header.interleave == envi_header::BSQ){		//if the infile is bsq file
  			if(header.data_type ==envi_header::float32)
  				((bsq<float>*)file)->set_buffer_raw(bytes);
  			else if(header.data_type == envi_header::float64)
  				((bsq<double>*)file)->set_buffer_raw(bytes);
  			else
  				std::cout<<"ERROR: unidentified data type"<<std::endl;
  		}
  
  		else if(header.interleave == envi_header::BIL){		//if the infile is bil file
  			if(header.data_type ==envi_header::float32)
  				((bil<float>*)file)->set_buffer_raw(bytes);
  			else if(header.data_type == envi_header::float64)
  				((bil<double>*)file)->set_buffer_raw(bytes);
  			else
  				std::cout<<"ERROR: unidentified data type"<<std::endl;
  		}
  
  		else if(header.interleave == envi_header::BIP){		//if the infile is bip file
  			if(header.data_type ==envi_header::float32)
  				((bip<float>*)file)->set_buffer_raw(bytes);
  			else if(header.data_type == envi_header::float64)
  				((bip<double>*)file)->set_buffer_raw(bytes);

  			else
  				std::cout<<"ERROR: unidentified data type"<<std::endl;
  		}
  
  		else{
  			std::cout<<"ERROR: unidentified file type"<<std::endl;
  			exit(1);
  		}
  	}
  

  	/// Returns the size of the data type in bytes
  	unsigned int type_size(){
  		if(header.data_type == envi_header::float32) return 4;
  		if(header.data_type == envi_header::float64) return 8;
  
  		exit(1);
  	}
  

  	size_t X(){ return header.samples; }
  	size_t Y(){ return header.lines; }
  	size_t Z(){ return header.bands; }
  	size_t B(){ return Z();	}
  
  	/// Return the size of the data set in bytes
  	size_t bytes(){
  		return X() * Y() * Z() * type_size();
  	}
  

  	/// Returns the progress of the current processing operation as a percentage
  	void reset_progress(){
  
  		if(header.interleave == envi_header::BSQ){		//if the infile is bsq file
  			if(header.data_type ==envi_header::float32)

  				((bsq<float>*)file)->reset_progress();

  			else if(header.data_type == envi_header::float64)

  				((bsq<double>*)file)->reset_progress();

  			else
  				std::cout<<"ERROR: unidentified data type"<<std::endl;
  		}
  
  		else if(header.interleave == envi_header::BIL){		//if the infile is bil file
  			if(header.data_type ==envi_header::float32)

  				((bil<float>*)file)->reset_progress();

  			else if(header.data_type == envi_header::float64)

  				((bil<double>*)file)->reset_progress();

  			else
  				std::cout<<"ERROR: unidentified data type"<<std::endl;
  		}
  
  		else if(header.interleave == envi_header::BIP){		//if the infile is bip file
  			if(header.data_type ==envi_header::float32)

  				((bip<float>*)file)->reset_progress();

  			else if(header.data_type == envi_header::float64)

  				((bip<double>*)file)->reset_progress();

  			else
  				std::cout<<"ERROR: unidentified data type"<<std::endl;
  		}
  
  		else{
  			std::cout<<"ERROR: unidentified file type"<<std::endl;
  			exit(1);
  		}
  	}
  
  	/// Returns the progress of the current processing operation as a percentage

  	double progress(){

  
  		if(header.interleave == envi_header::BSQ){		//if the infile is bsq file
  			if(header.data_type ==envi_header::float32)

  				return ((bsq<float>*)file)->get_progress();

  			else if(header.data_type == envi_header::float64)

  				return ((bsq<double>*)file)->get_progress();

  			else
  				std::cout<<"ERROR: unidentified data type"<<std::endl;
  		}
  
  		else if(header.interleave == envi_header::BIL){		//if the infile is bil file
  			if(header.data_type ==envi_header::float32)

  				return ((bil<float>*)file)->get_progress();

  			else if(header.data_type == envi_header::float64)

  				return ((bil<double>*)file)->get_progress();

  			else
  				std::cout<<"ERROR: unidentified data type"<<std::endl;
  		}
  
  		else if(header.interleave == envi_header::BIP){		//if the infile is bip file
  			if(header.data_type ==envi_header::float32)

  				return ((bip<float>*)file)->get_progress();

  			else if(header.data_type == envi_header::float64)

  				return ((bip<double>*)file)->get_progress();

  			else
  				std::cout<<"ERROR: unidentified data type"<<std::endl;
  		}
  
  		else{
  			std::cout<<"ERROR: unidentified file type"<<std::endl;
  		}
  		return 0;
  	}
  

  	/// Returns the progress of the current processing operation as a percentage
  	size_t data_rate(){
  
  		if(header.interleave == envi_header::BSQ){		//if the infile is bsq file
  			if(header.data_type ==envi_header::float32)
  				return ((bsq<float>*)file)->get_data_rate();
  			else if(header.data_type == envi_header::float64)
  				return ((bsq<double>*)file)->get_data_rate();
  			else
  				std::cout<<"ERROR: unidentified data type"<<std::endl;
  		}
  
  		else if(header.interleave == envi_header::BIL){		//if the infile is bil file
  			if(header.data_type ==envi_header::float32)
  				return ((bil<float>*)file)->get_data_rate();
  			else if(header.data_type == envi_header::float64)
  				return ((bil<double>*)file)->get_data_rate();
  			else
  				std::cout<<"ERROR: unidentified data type"<<std::endl;
  		}
  
  		else if(header.interleave == envi_header::BIP){		//if the infile is bip file
  			if(header.data_type ==envi_header::float32)
  				return ((bip<float>*)file)->get_data_rate();
  			else if(header.data_type == envi_header::float64)
  				return ((bip<double>*)file)->get_data_rate();
  			else
  				std::cout<<"ERROR: unidentified data type"<<std::endl;
  		}
  
  		else{
  			std::cout<<"ERROR: unidentified file type"<<std::endl;
  		}
  		return 0;
  	}
  

  	/// Allocate memory for a new ENVI file based on the current interleave format (BIP, BIL, BSQ) and data type.

  	void allocate(){

  
  		file = NULL;	//set file to a NULL pointer
  

  		if(header.interleave == envi_header::BSQ){

  			if(header.data_type ==envi_header::float32)

  				file = new bsq<float>();

  			else if(header.data_type == envi_header::float64)

  				file = new bsq<double>();

  		}
  		else if(header.interleave == envi_header::BIP){
  			if(header.data_type ==envi_header::float32)

  				file = new bip<float>();

  			else if(header.data_type == envi_header::float64)

  				file = new bip<double>();

  		}
  		else if(header.interleave == envi_header::BIL){
  			if(header.data_type ==envi_header::float32)

  				file = new bil<float>();

  			else if(header.data_type == envi_header::float64)

  				file = new bil<double>();

  		}

  	}
  

  	/// Open a previously opened ENVI file
  	bool open(){

  
  		//load the file

  		if(header.interleave == envi_header::BSQ) {		//if the infile is bsq file
  			if(header.data_type == envi_header::float32) {

  				return ((bsq<float>*)file)->open(fname, header.samples, header.lines, header.bands, header.header_offset, header.wavelength);

  			}
  			else if(header.data_type == envi_header::float64) {

  				return ((bsq<double>*)file)->open(fname, header.samples, header.lines, header.bands, header.header_offset, header.wavelength);

  			}
  			else
  				return false;
  		}

  		else if(header.interleave == envi_header::BIL) {		//if the infile is bil file
  			if(header.data_type == envi_header::float32) {

  				return ((bil<float>*)file)->open(fname, header.samples, header.lines, header.bands, header.header_offset, header.wavelength);

  			}
  			else if(header.data_type == envi_header::float64) {

  				return ((bil<double>*)file)->open(fname, header.samples, header.lines, header.bands, header.header_offset, header.wavelength);

  			}
  			else
  				return false;
  		}

  		else if(header.interleave == envi_header::BIP) {		//if the infile is bip file
  			if(header.data_type == envi_header::float32) {

  				return ((bip<float>*)file)->open(fname, header.samples, header.lines, header.bands, header.header_offset, header.wavelength);

  			}
  			else if(header.data_type == envi_header::float64) {

  				return ((bip<double>*)file)->open(fname, header.samples, header.lines, header.bands, header.header_offset, header.wavelength);

  			}
  			else

  				return false;

  		}
  

  		return true;
  

  
  	}
  

  	/// Open an Agilent binary file as an ENVI stream
  	bool open_agilent(std::string filename){
  		fname = filename;												//store the file name
  
  		//Open the file temporarily to get the header information
  		FILE* f = fopen(filename.c_str(), "r");							//open the binary file for reading
  		if(f == NULL) return false;										//return false if no file is opened
  
  		fseek(f, 9, SEEK_SET);											//seek to the number of bands
  		short b;														//allocate space for the number of bands

  		size_t nread = fread(&b, sizeof(short), 1, f);					//read the number of bands
  		if(nread != 1){
  			std::cout<<"Error reading band number from Agilent file."<<std::endl;
  			exit(1);
  		}

  		fseek(f, 13, SEEK_CUR);											//skip the the x and y dimensions
  		short x, y;

  		nread = fread(&x, sizeof(short), 1, f);									//read the image x and y size
  		if(nread != 1){
  			std::cout<<"Error reading X dimension from Agilent file."<<std::endl;
  			exit(1);
  		}
  		nread = fread(&y, sizeof(short), 1, f);
  		if(nread != 1){
  			std::cout<<"Error reading Y dimension from Agilent file."<<std::endl;
  			exit(1);
  		}

  		fclose(f);														//close the file
  
  		//store the information from the Agilent header in the ENVI header
  		header.bands = b;
  		header.samples = x;
  		header.lines = y;
  		header.data_type = envi_header::float32;						//all values are 32-bit floats
  		header.header_offset = 1020;									//number of bytes in an Agilent binary header
  		header.interleave = envi_header::BSQ;							//all Agilent binary files are BSQ
  
  		allocate();														//allocate the streaming file object
  		open();															//open the file for streaming
  
  		return true;
  	}
  

  	/// Open an existing ENVI file given the filename and a header structure
  
  	/// @param filename is the name of the ENVI binary file
  	/// @param header is an ENVI header structure
  	bool open(std::string filename, stim::envi_header h){

  
  		header = h;							//store the header
  		fname = filename;					//save the filename
  

  		allocate();
  

  		return open();						//open the ENVI file;
  		
  

  	}
  
  	/// Open an existing ENVI file given the file and header names.
  
  	/// @param filename is the name of the ENVI binary file
  	/// @param headername is the name of the ENVI header file
  	bool open(std::string filename, std::string headername){
  
  		//allocate memory
  		//allocate();
  
  		stim::envi_header h;
  		h.load(headername);
  
  		//load the header
  		//header.load(headername);
  
  		return open(filename, h);

  	}
  

  	/// Normalize a hyperspectral ENVI file given a band number and threshold.
  
  	/// @param outfile is the name of the normalized file to be output
  	/// @param band is the band label to be output
  	/// @param threshold is a threshold value specified such that normalization will only be done to values in the band > threshold (preventing division by small numbers)

  	bool ratio(std::string outfile, double band, unsigned char* mask = NULL, bool PROGRESS = false){

  		header.save(outfile + ".hdr");

  		if(header.interleave == envi_header::BSQ){		//if the infile is bsq file

  			if(header.data_type ==envi_header::float32)

  				return ((bsq<float>*)file)->ratio(outfile, band, mask, PROGRESS);

  			else if(header.data_type == envi_header::float64)

  				return ((bsq<double>*)file)->ratio(outfile,band, mask, PROGRESS);

  			else
  				std::cout<<"ERROR: unidentified data type"<<std::endl;
  		}

  		else if(header.interleave == envi_header::BIL){		//if the infile is bil file
  			if(header.data_type ==envi_header::float32)

  				return ((bil<float>*)file)->ratio(outfile, band, mask, PROGRESS);

  			else if(header.data_type == envi_header::float64)

  				return ((bil<double>*)file)->ratio(outfile,band, mask, PROGRESS);

  			else
  				std::cout<<"ERROR: unidentified data type"<<std::endl;
  		}
  
  		else if(header.interleave == envi_header::BIP){		//if the infile is bip file
  			if(header.data_type ==envi_header::float32)

  				return ((bip<float>*)file)->ratio(outfile, band, mask, PROGRESS);

  			else if(header.data_type == envi_header::float64)

  				return ((bip<double>*)file)->ratio(outfile,band, mask, PROGRESS);

  			else
  				std::cout<<"ERROR: unidentified data type"<<std::endl;
  		}
  
  		else{
  			std::cout<<"ERROR: unidentified file type"<<std::endl;
  			exit(1);
  		}

  		return false;

  	}
  

  	/// Perform vector normalization on the ENVI image
  	void normalize(std::string outfile, unsigned char* mask = NULL, bool PROGRESS = false){
  		header.save(outfile + ".hdr");
  		if(header.interleave == envi_header::BSQ){		//if the infile is bsq file
  			if(header.data_type ==envi_header::float32)
  				((bsq<float>*)file)->normalize(outfile, mask, PROGRESS);
  			else if(header.data_type == envi_header::float64)
  				((bsq<double>*)file)->normalize(outfile, mask, PROGRESS);
  			else
  				std::cout<<"ERROR: unidentified data type"<<std::endl;
  		}
  
  		else if(header.interleave == envi_header::BIL){		//if the infile is bil file
  			if(header.data_type ==envi_header::float32)
  				((bil<float>*)file)->normalize(outfile, mask, PROGRESS);
  			else if(header.data_type == envi_header::float64)
  				((bil<double>*)file)->normalize(outfile, mask, PROGRESS);
  			else
  				std::cout<<"ERROR: unidentified data type"<<std::endl;
  		}
  
  		else if(header.interleave == envi_header::BIP){		//if the infile is bip file
  			if(header.data_type ==envi_header::float32)
  				((bip<float>*)file)->normalize(outfile, mask, PROGRESS);
  			else if(header.data_type == envi_header::float64)
  				((bip<double>*)file)->normalize(outfile, mask, PROGRESS);
  			else
  				std::cout<<"ERROR: unidentified data type"<<std::endl;
  		}
  
  		else{
  			std::cout<<"ERROR: unidentified file type"<<std::endl;
  			exit(1);
  		}
  	}
  

  	/// Performs piecewise linear baseline correction of a hyperspectral file/
  
  	/// @param outfile is the file name for the baseline corrected output
  	/// @param w is a list of band labels to serve as baseline points (zero values)

  	bool baseline(std::string outfile, std::vector<double> w, unsigned char* mask = NULL, bool PROGRESS = false){
  		header.save(outfile + ".hdr");

  		if(header.interleave == envi_header::BSQ){		//if the infile is bsq file
  			if(header.data_type ==envi_header::float32)

  				return ((bsq<float>*)file)->baseline(outfile, w, mask, PROGRESS);

  			else if(header.data_type == envi_header::float64)

  				return ((bsq<double>*)file)->baseline(outfile,w, mask, PROGRESS);

  			else{
  				std::cout<<"ERROR: unidentified data type"<<std::endl;
  				exit(1);
  			}
  		}
  
  		else if(header.interleave == envi_header::BIL){		//if the infile is bil file
  			if(header.data_type ==envi_header::float32)

  				return ((bil<float>*)file)->baseline(outfile, w, mask, PROGRESS);

  			else if(header.data_type == envi_header::float64)

  				return ((bil<double>*)file)->baseline(outfile, w, mask, PROGRESS);

  			else{
  				std::cout<<"ERROR: unidentified data type"<<std::endl;
  				exit(1);
  			}
  		}
  

  		else if(header.interleave == envi_header::BIP){		//if the infile is bip file

  			if(header.data_type ==envi_header::float32)

  				return ((bip<float>*)file)->baseline(outfile, w, mask, PROGRESS);

  			else if(header.data_type == envi_header::float64)

  				return ((bip<double>*)file)->baseline(outfile, w, mask, PROGRESS);

  			else{
  				std::cout<<"ERROR: unidentified data type"<<std::endl;
  				exit(1);
  			}
  		}
  
  		else{
  			std::cout<<"ERROR: unidentified file type"<<std::endl;
  			exit(1);
  		}
  	}
  

  	/// Project an array of coefficients onto a basis matrix
  
  	void project(std::string outfile, double* center, double* basis, size_t M, std::vector<double> bands, unsigned char* mask, int cuda_device = 0, bool PROGRESS = false) {

  		if(header.interleave == envi_header::BSQ){		//if the infile is bsq file
  			std::cout<<"ERROR: BSQ projection not supported"<<std::endl;
  			exit(1);
  		}
  
  		else if(header.interleave == envi_header::BIL){		//if the infile is bil file
  			std::cout<<"ERROR: BIL projection not supported"<<std::endl;
  			exit(1);
  		}
  
  		else if(header.interleave == envi_header::BIP){		//if the infile is bip file
  			if(header.data_type ==envi_header::float32)

  				((bip<float>*)file)->project(outfile, center, basis, M, mask, cuda_device, PROGRESS);

  			else if(header.data_type == envi_header::float64)

  				((bip<double>*)file)->project(outfile, center, basis, M, mask, cuda_device, PROGRESS);

  			else{
  				std::cout<<"ERROR: unidentified data type"<<std::endl;
  				exit(1);
  			}
  		}
  
  		stim::envi_header out_hdr = header;							
  		out_hdr.bands = M;											//set the number of bands in the output header

  		out_hdr.wavelength = bands;

  		out_hdr.band_names.clear();

  		out_hdr.save(outfile + ".hdr");								//save the output header
  	}
  
  	void inverse(std::string outfile, double* center, double* basis, unsigned long long B, unsigned long long C = 0, bool PROGRESS = false){
  		if(header.interleave == envi_header::BSQ){		//if the infile is bsq file
  			std::cout<<"ERROR: BSQ projection not supported"<<std::endl;
  			exit(1);
  		}
  
  		else if(header.interleave == envi_header::BIL){		//if the infile is bil file
  			std::cout<<"ERROR: BIL projection not supported"<<std::endl;
  			exit(1);
  		}
  
  		else if(header.interleave == envi_header::BIP){		//if the infile is bip file
  			if(header.data_type ==envi_header::float32)
  				((bip<float>*)file)->inverse(outfile, center, basis, B, C, PROGRESS);
  			else if(header.data_type == envi_header::float64)
  				((bip<double>*)file)->inverse(outfile, center, basis, B, C, PROGRESS);
  			else{
  				std::cout<<"ERROR: unidentified data type"<<std::endl;
  				exit(1);
  			}
  		}
  
  		stim::envi_header out_hdr = header;							
  		out_hdr.bands = B;											//set the number of bands in the output header
  		out_hdr.save(outfile + ".hdr");								//save the output header
  	}
  

  	/// Converts ENVI files between interleave types (BSQ, BIL, and BIP)
  
  	/// @param outfile is the file name for the converted output
  	/// @param interleave is the interleave format for the destination file

  	bool convert(std::string outfile, stim::envi_header::interleaveType interleave, bool PROGRESS = false, bool VERBOSE = false, bool OPTIMIZATION = true){

  		if(header.interleave == envi_header::BSQ){			//if the infile is bsq file
  

  			if(header.data_type ==envi_header::float32){		//ERROR

  				if(interleave == envi_header::BSQ){
  					std::cout<<"ERROR:  is already BSQ file"<<std::endl;
  					exit(1);
  				}

  				else if(interleave == envi_header::BIL)			//convert BSQ -> BIL

  					((bsq<float>*)file)->bil(outfile, PROGRESS, VERBOSE, OPTIMIZATION);

  				else if(interleave == envi_header::BIP){			//ERROR

  					//std::cout<<"ERROR: conversion from BSQ to BIP isn't practical; use BSQ->BIL->BIP instead"<<std::endl;