codebase / stimlib

  #ifndef ENVI_HEADER_H
  #define ENVI_HEADER_H
  
  #include <string>
  #include <iostream>
  #include <fstream>
  #include <sstream>
  #include <vector>
  #include <algorithm>
  #include <stdlib.h>
  #include <cmath>
  
  //information from an ENVI header file
  //A good resource can be found here: http://www.exelisvis.com/docs/enviheaderfiles.html
  
  namespace stim{
  
  struct envi_header
  {
  	enum dataType {dummy0, int8, int16, int32, float32, float64, complex32, dummy7, dummy8, complex64, dummy10, dummy11, uint16, uint32, int64, uint64};
  	enum interleaveType {BIP, BIL, BSQ};	//bip = Z,X,Y; bil = X,Z,Y; bsq = X,Y,Z
  	enum endianType {endianLittle, endianBig};
  
  	std::string name;
  
  	std::string description;
  
  	size_t samples;	//x-axis
  	size_t lines;		//y-axis
  	size_t bands;		//spectral axis
  	size_t header_offset;		//header offset for binary file (in bytes)
  	std::string file_type;			//should be "ENVI Standard"
  
  	envi_header::dataType data_type;			//data representation; common value is 4 (32-bit float)
  
  
  	envi_header::interleaveType interleave;
  
  	std::string sensor_type;		//not really used
  
  	envi_header::endianType byte_order;			//little = least significant bit first (most systems)
  
  	double x_start, y_start;		//coordinates of the upper-left corner of the image
  	std::string wavelength_units;	//stored wavelength units
  	std::string z_plot_titles[2];
  
  	double pixel_size[2];			//pixel size along X and Y
  
  	std::vector<std::string> band_names;	//name for each band in the image
  	std::vector<double> wavelength;		//wavelength for each band
  
  	void init(){
          name = "";
  
  		//specify default values for a new or empty ENVI file
  		samples = 0;
  		lines = 0;
  		bands = 0;
  		header_offset = 0;
  		data_type = float32;
  		interleave = BSQ;
  		byte_order = endianLittle;
  		x_start = y_start = 0;
  		pixel_size[0] = pixel_size[1] = 1;
  
  		//strings
  		file_type = "ENVI Standard";
  		sensor_type = "Unknown";
  		wavelength_units = "Unknown";
  		z_plot_titles[0] = z_plot_titles[1] = "Unknown";
  	}
  
  	envi_header(){
  		init();
  	}
  	envi_header(std::string name){
  		init();
  		load(name);
  	}
  
  	//sets the wavelength vector given a starting value and uniform step size
  	void set_wavelengths(double start, double step){
  		size_t B = bands;						//get the number of bands
  		wavelength.resize(B);
  		for(size_t b = 0; b < B; b++)
  			wavelength[b] = start + b * step;
  	}
  
  	std::string trim(std::string line){
  
  		if(line.length() == 0)
  			return line;
  		//trims whitespace from the beginning and end of line
  		size_t start_i, end_i;
  		for(start_i=0; start_i < line.length(); start_i++)
  			if(line[start_i] != 32)
  			{
  				break;
  			}
  
  		for(end_i = line.length()-1; end_i >= start_i; end_i--)
  			if(line[end_i] != ' ' && line[end_i] != '\r')
  			{
  				break;
  			}
  
  		return line.substr(start_i, end_i - start_i+1);
  	}
  
  
  	std::string get_token(std::string line){
  		//returns a variable name; in this case we look for the '=' sign
  		size_t i = line.find_first_of('=');
  
  		std::string result;
  		if(i != std::string::npos)
  			result = trim(line.substr(0, i-1));
  
  		return result;
  	}
  
  	std::string get_data_str(std::string line){
  		size_t i = line.find_first_of('=');
  
  		std::string result;
  		if(i != std::string::npos)
  			result = trim(line.substr(i+1));
  		else
  		{
  			std::cout<<"ENVI Header error - data not found for token: "<<get_token(line)<<std::endl;
  			exit(1);
  		}
  		return result;
  	}
  
  	std::string get_brace_str(std::string token, std::string line, std::ifstream &file)
  	{
  		//this function assembles all of the characters between curly braces
  		//this is how strings are defined within an ENVI file
  
  		std::string result;
  
  		//first, find the starting brace
  		size_t i;
  		do
  		{
  			i = line.find_first_of('{');
  			if(i != std::string::npos)
  				break;
  		}while(file);
  
  		//if i is still npos, we have reached the end of the file without a brace...something is wrong
  		if(i == std::string::npos)
  		{
  			std::cout<<"ENVI Header error - string token declared without being defined: "<<token<<std::endl;
  			exit(1);
  		}
  		line = line.substr(i+1);
  
  		//copy character data into the result string until we find a closing brace
  		while(file)
  		{
  			i = line.find_first_of('}');
  
  
  			if(i != std::string::npos)
  			{
  				result += line.substr(0, i);
  				break;
  			}
  			else
  				result += line;
  
  			getline(file, line);
  		}
  
  		if(i == std::string::npos)
  		{
  			std::cout<<"ENVI Header error - string token declared without a terminating '}': "<<token<<std::endl;
  			exit(1);
  		}
  
  		return trim(result);
  	}
  
  	std::vector<std::string> get_string_seq(std::string token, std::string sequence)
  	{
  		//this function returns a sequence of comma-delimited strings
  		std::vector<std::string> result;
  
  		std::string entry;
  		size_t i;
  		do
  		{
  			i = sequence.find_first_of(',');
  			entry = sequence.substr(0, i);
  			sequence = sequence.substr(i+1);
  			result.push_back(trim(entry));
  		}while(i != std::string::npos);
  
  		return result;
  	}
  
  	std::vector<double> get_double_seq(std::string token, std::string sequence)
  	{
  		//this function returns a sequence of comma-delimited strings
  		std::vector<double> result;
  		std::string entry;
  		size_t i;
  		do
  		{
  			i = sequence.find_first_of(',');
  			entry = sequence.substr(0, i);
  			sequence = sequence.substr(i+1);
  			result.push_back(atof(entry.c_str()));
  			//std::cout<<entry<<"   ";
  		}while(i != std::string::npos);
  
  		return result;
  	}
  
  	bool load(std::string filename)
  	{
  		//open the header file
  		std::ifstream file(filename.c_str());
  
  		if(!file)
  		{
  			std::cout<<"ERROR: unable to open header file:   "<<filename<<std::endl;
              return false;
  		}
  
  		//the first line should just be "ENVI"
  		std::string line;
  		file>>line;
  		if(line != "ENVI")
  		{
  			std::cout<<"ENVI Header Error: The header doesn't appear to be an ENVI file. The first line should be 'ENVI'."<<std::endl;
  			exit(1);
  		}
  
  		//for each line in the file, get the token
  		std::string token;
  
          //get a line
          getline(file, line);
  		while(file)
  		{
  
  
  
  			//get the token
  			token = get_token(line);
  
  			if(token == "description")
  				description = get_brace_str(token, line, file);
  			else if(token == "band names")
  			{
  				std::string string_sequence = get_brace_str(token, line, file);
  				band_names = get_string_seq(token, string_sequence);
  			}
  			else if(token == "wavelength")
  			{
  				std::string string_sequence = get_brace_str(token, line, file);
  				wavelength = get_double_seq(token, string_sequence);
  			}
  			else if(token == "pixel size")
  			{
  				std::string string_sequence = get_brace_str(token, line, file);
  				std::vector<double> pxsize = get_double_seq(token, string_sequence);
  				pixel_size[0] = pxsize[0];
  				pixel_size[1] = pxsize[1];
  			}
  			else if(token == "z plot titles")
  			{
  				std::string string_sequence = get_brace_str(token, line, file);
  				std::vector<std::string> titles = get_string_seq(token, string_sequence);
  				z_plot_titles[0] = titles[0];
  				z_plot_titles[1] = titles[1];
  			}
  
  			else if(token == "samples")
  				samples = atoi(get_data_str(line).c_str());
  			else if(token == "lines")
  				lines = atoi(get_data_str(line).c_str());
  			else if(token == "bands")
  				bands = atoi(get_data_str(line).c_str());
  			else if(token == "header offset")
  				header_offset = atoi(get_data_str(line).c_str());
  			else if(token == "file type")
  				file_type = get_data_str(line);
  			else if(token == "data type")
  				data_type = (dataType)atoi(get_data_str(line).c_str());
  			else if(token == "interleave")
  			{
  				std::string interleave_str = get_data_str(line);
  				if(interleave_str == "bip")
  					interleave = BIP;
  				else if(interleave_str == "bil")
  					interleave = BIL;
  				else if(interleave_str == "bsq")
  					interleave = BSQ;
  			}
  			else if(token == "sensor type")
  				sensor_type = get_data_str(line);
  			else if(token == "byte order")
  				byte_order = (endianType)atoi(get_data_str(line).c_str());
  			else if(token == "x start")
  				x_start = atof(get_data_str(line).c_str());
  			else if(token == "y start")
  				y_start = atof(get_data_str(line).c_str());
  			else if(token == "wavelength units")
  				wavelength_units = get_data_str(line);
  
              //get the next line
              getline(file, line);
  		}
  
  		//make sure the number of bands matches the number of wavelengths
  		size_t wavelengths = wavelength.size();
  		if(wavelengths && bands != wavelengths)
  		{
              std::cout<<"ENVI Header Error -- Number of wavelengths doesn't match the number of bands.  Bands = "<<bands<<", Wavelengths = "<<wavelength.size()<<std::endl;
              exit(1);
          }
  
  		if(wavelength.size() == 0){											//if no list of wavelengths is given, create one based on band numbers
  			wavelength.resize(bands);
  			for(size_t b = 0; b < bands; b++)
  				wavelength[b] = (double)b;
  		}
  
  		//close the file
  		file.close();
  
  		//set the file name
  		name = filename;
  
  		return true;
  	}
  
  	void save(std::string filename)
  	{
  		//open a file
  		std::ofstream outfile(filename.c_str());
  
  		//write the ENVI type identifier
  		outfile<<"ENVI"<<std::endl;
  
  		//output all of the data
  		outfile<<"description = {"<<std::endl;
  		outfile<<"  "<<description<<"}"<<std::endl;
  
  		outfile<<"samples = "<<samples<<std::endl;
  		outfile<<"lines = "<<lines<<std::endl;
  		outfile<<"bands = "<<bands<<std::endl;
  		outfile<<"header offset = "<<header_offset<<std::endl;
  		outfile<<"file type = "<<file_type<<std::endl;
  		outfile<<"data type = "<<data_type<<std::endl;
  		outfile<<"interleave = ";
  		if(interleave == BIP)
  			outfile<<"bip";
  		if(interleave == BIL)
  			outfile<<"bil";
  		if(interleave == BSQ)
  			outfile<<"bsq";
  			outfile<<std::endl;
  		outfile<<"sensor type = "<<sensor_type<<std::endl;
  		outfile<<"byte order = "<<byte_order<<std::endl;
  		outfile<<"x start = "<<x_start<<std::endl;
  		outfile<<"y start = "<<y_start<<std::endl;
  		outfile<<"wavelength units = "<<wavelength_units<<std::endl;
  		outfile<<"z plot titles = {";
  			outfile<<z_plot_titles[0]<<", "<<z_plot_titles[1]<<"}"<<std::endl;
  		outfile<<"pixel size = {"<<pixel_size[0]<<", "<<pixel_size[1]<<", units=Meters}"<<std::endl;
  		if(band_names.size() > 0)
  		{
  			outfile<<"band names = {"<<std::endl;
  			for(unsigned int i=0; i<band_names.size(); i++)
  			{
  				outfile<<band_names[i];
  				if(i < band_names.size() - 1)
  					outfile<<", ";
  			}
  			outfile<<"}"<<std::endl;
  		}
  		//if there is a list of wavelengths, output it
  		if(wavelength.size()){
  			outfile<<"wavelength = {"<<std::endl;
  				for(unsigned int i=0; i<wavelength.size()-1; i++)
  					outfile<<wavelength[i]<<", ";
  				outfile<<wavelength.back()<<"}"<<std::endl;
  		}
  
  		outfile.close();
  	}
  
  	void save()
  	{
          //std::cout<<"ENVI Header Name: "<<name<<std::endl;
  		save(name);
  	}
  
  	//returns the size of a single value (in bytes)
  	unsigned int valsize()
  	{
  		switch(data_type)
  		{
  		case int8:			//1 = 8-bit byte
  			return 1;
  		case int16:			//16-bit signed integer
  		case uint16:		//16-bit unsigned integer
  			return 2;
  		case int32:			//32-bit signed long integer
  		case float32:		//32-bit floating point
  		case uint32:		//32-bit unsigned long integer
  			return 4;
  		case float64:		//64-bit double-precision floating point
  		case complex32:		//32-bit complex value
  		case int64:		    //64-bit signed long integer
  		case uint64:		//64-bit unsigned long integer
  			return 8;
  		case complex64:		//64-bit complex value
  			return 16;
  		default:
  			return 0;
  		}
  	}
  
  	//return the number of bytes that SHOULD be in the data file
  	size_t data_bytes(){
  		return samples * lines * bands * valsize() + header_offset;
  	}
  	
  
  	/// Convert an interleave type to a string
  	static std::string interleave_str(interleaveType t){
  		switch(t){
  		case stim::envi_header::BSQ:
  			return std::string("BSQ");
  		case stim::envi_header::BIL:
  			return std::string("BIL");
  		case stim::envi_header::BIP:
  			return std::string("BIP");
  		default:
  			std::cout<<"ERROR in stim::envi_header::typestr() - unrecognized type"<<std::endl;
  			exit(1);
  		}
  	}
  
  	/// Convert the current interleave type to a string
  	std::string interleave_str(){
  		return interleave_str(interleave);
  	}
  
  	/// Convert a wavelength to a band index (or a pair of surrounding band indices)
  	///		if the file doesn't specify wavelengths, w is assumed to be a band index
  	std::vector<size_t> band_index(double w){
  		std::vector<size_t> idx;										//create an empty array of indices
  		if(wavelength.size() == 0){										//if a wavelength vector doesn't exist, assume the passed value is a band
  			if(w < 0 || w > bands-1) return idx;						//if the band is outside the given band range, return an empty vector
  			size_t low, high;											//allocate space for the floor and ceiling
  			low = (size_t)std::floor(w);										//calculate the floor
  			high = (size_t)std::ceil(w);										//calculate the ceiling
  			if(low == high)												//if the floor and ceiling are the same
  				idx.push_back(low);										//return a vector with one element (the given w matches a band exactly)
  			else{
  				idx.resize(2);											//otherwise return the floor and ceiling
  				idx[0] = low;
  				idx[1] = high;
  			}
  			return idx;
  		}
  		else if(w < wavelength[0] || w > wavelength[bands-1]) return idx;	//if the wavelength range is outside of the file, return an empty array
  
  		for(size_t b = 0; b < bands; b++){							//for each band in the wavelength vector
  			if(wavelength[b] == w){									//if an exact match is found
  				idx.push_back(b);										//add the band to the array and return it
  				return idx;
  			}
  			if(wavelength[b] >= w){ 								//if the current wavelength exceeds w
  				idx.resize(2);
  				idx[0] = b-1;										//push both the previous and current band index
  				idx[1] = b;
  				return idx;											//return the pair of indices
  			}
  		}
  		return idx;
  	}
  
  	/// Convert a wavelength range to a list of bands
  	std::vector<size_t> band_indices(double w0, double w1){
  
  		size_t idx0 = 0;
  		size_t idx1 = 0;									//declare the interval indices for the band range
  
  		//get the indices for the first wavelength
  		std::vector<size_t> r;
  
  		if(w0 > wavelength[bands-1] || w1 < wavelength[0]){
  			std::cout<<"ERROR in envi_header::band_indices - wavelengths are completely out of range: ["<<w0<<", "<<w1<<"]"<<std::endl;
  			exit(1);
  		}
  		if(w0 < wavelength[0])								//if the lower wavelength is outside the file range
  			idx0 = 0;										//just set it to zero
  		else{
  			r = band_index(w0);
  			if(r.size() == 1)								//if there is an exact match, set the first interval index
  				idx0 = r[0];
  			else
  				idx0 = r[1];								//otherwise save the highest band index (this is the first band with wavelength > w0)
  		}
  		//get the indices for the second wavelength
  		if(w1 > wavelength[bands-1])
  			idx1 = bands-1;
  		else{
  			r = band_index(w1);
  			if(r.size() == 0)
  				idx1 = bands - 1;
  			else
  				idx1 = r[0];									//take the lowest band index (this is the last band < w1)
  		}
  
  		size_t n = idx1 - idx0 + 1;							//calculate the number of bands in this interval
  		r.resize(n);										//resize the band index array
  		for(size_t b = 0; b < n; b++){						//for each band in the interval
  			r[b] = idx0 + b;								//	insert the band in the array
  		}
  
  		return r;											//return the index array
  	}
  };		//end EnviHeader
  }
  #endif