stim/envi/envi_header.h

#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