added binary file framework

David Mayerich
1 parent 0ef519a4
Showing 7 changed files with 550 additions and 394 deletions Show diff stats
envi/bil.h
envi/binary.h
envi/envi.h
math/complexfield.cuh
math/field.cuh
math/function.h
optics/mirst-1d.cuh
+#ifndef RTS_BIL_H
+#define RTS_BIL_H
+
+#include "../envi/envi.h"
+
+namespace rts{
+
+//This class specifies a BIL binary file that can be streamed to and from storage media
+template< typename T >
+class bil{
+
+protected:
+
+	rts::envi header;
+
+
+public:
+
+	bil(){
+		init();
+	}
+
+};
+
+}
+
+#endif
+#ifndef RTS_BINARY_H
+#define RTS_BINARY_H
+
+#include "../envi/envi.h"
+
+#include <fstream>
+#include <sys/stat.h>
+
+namespace rts{
+
+//This class contains a bunch of functions useful for multidimensional binary file access
+template< typename T, unsigned int D = 3 >
+class binary{
+
+protected:
+
+	std::fstream file;		//file stream used for reading and writing
+	std::string name;		//file name
+
+	unsigned int R[D];		//resolution
+	unsigned int header;	//header size (in bytes)
+
+
+	//basic initialization
+	void init(){
+		memset(R, 0, sizeof(unsigned int) * D);		//initialize the resolution to zero
+		header = 0;									//initialize the header size to zero
+	}
+
+	//returns the file size
+	unsigned int get_file_size(){
+
+		struct stat results;
+		if(stat(name.c_str(), &results) == 0)
+			return results.st_size;
+		else return 0;
+	}
+
+	bool test_file_size(){
+		unsigned int sum = header;			//initialize the sum (in bytes) to the header size
+		for(unsigned int i = 0; i<D; i++)	//iterate over each dimension, summing the byte size
+			sum += R[i] * sizeof(T);
+		if(sum == get_file_size()) return true;	//if the byte size matches the file size, we're golden
+		else return false;					//otherwise return an error
+
+	}
+
+	bool open_file(std::string filename){
+		//open the file as binary for reading and writing
+		file.open(filename.c_str(), std::ios::in | std::ios::out | std::ios::binary);
+
+		//if the file is successful
+		if(file){
+			name = filename;		//set the name
+			if(test_file_size())	//test the file size
+				return true;
+		}
+		
+		return false;
+	}
+
+public:
+
+	bool open(std::string filename, vec<unsigned int, D> r, unsigned int h = 0){
+		if(!open_file(filename)) return false;
+
+		for(unsigned int i = 0; i < D; i++)
+			R[i] = r[i];
+
+		header = h;
+
+		return test_file_size();
+	}
+
+	bool open(std::string filename, unsigned int X, unsigned int h = 0){
+		return open(filename, vec<unsigned int, 1>(X), h);
+	}
+
+	bool open(std::string filename, unsigned int X, unsigned int Y, unsigned int h = 0){
+		return open(filename, vec<unsigned int, 2>(X, Y), h);
+
+
+	}
+
+	bool open(std::string filename, unsigned int X, unsigned int Y, unsigned int Z, unsigned int h = 0){
+		return open(filename, vec<unsigned int, 3>(X, Y, Z), h);
+	}
+
+};
+
+}
+
+#endif
-#ifndef RTS_ENVI_H
-#define RTS_ENVI_H
+#ifndef ENVI_HEADER_H
+#define ENVI_HEADER_H
 #include <string>
-#include <iostream>
+#include <iostream>
 #include <fstream>
-#include <sstream>
-
-#include "rts/envi/envi_header.h"
-
-/*	This is a class for reading and writing ENVI binary files. This class will be updated as needed.
-
-What this class CAN currently do:
-	*)  Write band images to a BSQ file.
-	*)  Append band images to a BSQ file.
-	*)  Interpolate and query bands in a float32 BSQ file
-*/
-class EnviFile
+#include <sstream>
+#include <vector>
+#include <stdlib.h>
+
+//information from an ENVI header file
+//A good resource can be found here: http://www.exelisvis.com/docs/enviheaderfiles.html
+
+namespace rts{
+
+struct envi
 {
-	EnviHeader header;
-	FILE* data;
-	std::string mode;
+	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};
-	//a locked file has alread had data written..certain things can't be changed
-	//		accessor functions deal with these issues
-	bool locked;
+	std::string name;
-	void init()
-	{
-		locked = false;
-		data = NULL;
+	std::string description;
+
+	unsigned int samples;	//x-axis
+	unsigned int lines;		//y-axis
+	unsigned int bands;		//spectral axis
+	unsigned int header_offset;		//header offset for binary file (in bytes)
+	std::string file_type;			//should be "ENVI Standard"
+
+	envi::dataType data_type;			//data representation; common value is 4 (32-bit float)
+
+
+	envi::interleaveType interleave;
+
+	std::string sensor_type;		//not really used
+
+	envi::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
+
+	envi(){
+        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";
 	}
-	void readonly()
-	{
-		if(mode == "r")
+	std::string trim(std::string line){
+		//trims whitespace from the beginning and end of line
+		unsigned int 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] != ' ')
+			{
+				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 Error: changes cannot be made to a read-only file."<<std::endl;
+			std::cout<<"ENVI Header error - data not found for token: "<<get_token(line)<<std::endl;
 			exit(1);
 		}
+		return result;
 	}
-	//this function determines if the current software is capable of the requested change
-	bool caps(EnviHeader new_header)
+	std::string get_brace_str(std::string token, std::string line, std::ifstream &file)
 	{
-		readonly();		//if the file is read-only, throw an exception
+		//this function assembles all of the characters between curly braces
+		//this is how strings are defined within an ENVI file
-		//we currently only support BSQ
-		if(new_header.interleave == Envi::BIP || new_header.interleave == Envi::BIL)
+		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<<"This software only supports BSQ.";
-			return false;
+			std::cout<<"ENVI Header error - string token declared without being defined: "<<token<<std::endl;
+			exit(1);
 		}
+		line = line.substr(i+1);
-		//if the number of bands has changed
-		if(header.bands != new_header.bands)
+		//copy character data into the result string until we find a closing brace
+		while(file)
 		{
-			//the new header must be a BSQ file
-			if(new_header.interleave != Envi::BSQ)
+			i = line.find_first_of('}');
+
+
+			if(i != std::string::npos)
 			{
-				std::cout<<"ENVI Error: can only add bands to a BSQ file."<<std::endl;
-				return false;
+				result += line.substr(0, i);
+				break;
 			}
+			else
+				result += line;
+
+			getline(file, line);
 		}
-		//disallow anything that can't be done when the file is locked
-		if(locked)
+		if(i == std::string::npos)
 		{
-			if(header.lines != new_header.lines)
-			{
-				std::cout<<"ENVI Error: cannot change the number of lines from "<<header.lines<<" to "<<new_header.lines<<" in a locked BSQ file."<<std::endl;
-				return false;
-			}
-			if(header.samples != new_header.samples)
-			{
-				std::cout<<"ENVI Error: cannot change the number of samples in a locked BSQ file."<<std::endl;
-				return false;
-			}
-			if(header.data_type != new_header.data_type)
-			{
-				std::cout<<"ENVI Error: cannot change the data type of a locked file."<<std::endl;
-				return false;
-			}
-			if(header.byte_order != new_header.byte_order)
-			{
-				std::cout<<"ENVI Error: cannot change the byte order of a locked file."<<std::endl;
-				return false;
-			}
-		}
-
-		return true;
+			std::cout<<"ENVI Header error - string token declared without a terminating '}': "<<token<<std::endl;
+			exit(1);
+		}
+		return trim(result);
 	}
-	EnviHeader load_header(std::string header_file)
+	std::vector<std::string> get_string_seq(std::string token, std::string sequence)
 	{
-		EnviHeader new_header;
-		new_header.load(header_file);
-		return new_header;
+		//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;
 	}
-	void set_header(EnviHeader new_header)
+	std::vector<double> get_double_seq(std::string token, std::string sequence)
 	{
-		if(caps(new_header))
-			header = new_header;
-		else
-			exit(1);
-	}
-
-	bool test_exist(std::string filename)
-	{
-        //tests for the existance of the specified binary file
-        FILE* f;
-        long sz = 0;
-        f = fopen(filename.c_str(), "rb");
-        if(f == NULL)
-        {
-            //std::cout<<"ENVI Error: error testing file existance."<<std::endl;
-            //exit(1);
-			return false;
-        }
-
-        fseek(f, 0, SEEK_END);
-        sz = ftell(f);
-        fclose(f);
-
-        if(sz>0)
-            return true;
-        else
-            return false;
-
-    }
-
-    void get_surrounding_bands(double wavelength, unsigned int &low, unsigned int &high)
-    {
-        int i;
-        int nBands = header.wavelength.size();
-        low = high = 0;
-        //std::cout<<"Wavelength to search: "<<wavelength<<"     Number of Bands: "<<nBands<<std::endl;
-        for(i=0; i<nBands; i++)
-        {
-            if(header.wavelength[i] < wavelength)
-            {
-                if(header.wavelength[low] > wavelength)
-                    low = i;
-                if(header.wavelength[i] > header.wavelength[low])
-                    low = i;
-            }
-
-            if(header.wavelength[i] > wavelength)
-            {
-                if(header.wavelength[high] < wavelength)
-                    high = i;
-                if(header.wavelength[i] < header.wavelength[high])
-                    high = i;
-            }
-            //std::cout<<"Low: "<<low<<"     High: "<<high<<std::endl;
-        }
-        //exit(1);
-    }
-
-    void interpolate_band(void* result, void* A, void* B, double a)
-    {
-        //interpolate between two bands independent of data type
-        //  CURRENTLY ONLY HANDLES FLOAT
-        int X = header.samples;
-        int Y = header.lines;
-        int N = X * Y;
-
-        if(header.data_type != Envi::float32)
-        {
-            std::cout<<"ENVI Error: this software only handles interpolation of 32-bit floats."<<std::endl;
-            exit(1);
-        }
-
-		float v0, v1;
-        for(int n=0; n<N; n++)
-        {
-			v0 = ((float*)A)[n];
-			v1 = ((float*)B)[n];
-            float r = v0 * a + v1 * (1.0 - a);
-            //((float*)result)[n] = ((float*)A)[n] * a + ((float*)B)[n] * (1.0 - a);
-            ((float*)result)[n] = r;
-        }
-
-    }
-
-public:
-
-	EnviFile()
-	{
-        mode = "w";
-		init();
+		//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;
 	}
-	EnviFile(std::string filename, std::string file_mode = "r")
+	bool load(std::string filename)
 	{
-		init();
-		open(filename, filename + ".hdr", file_mode);
-	}
+		//open the header file
+		std::ifstream file(filename.c_str());
-	void open(std::string file_name, std::string header_name, std::string file_mode = "r")
-	{
-
-        //load the header file
-        //  if the file doesn't exist, that is okay unless it is read-only
-
-        //if the file is being written, create a new header
-        if(file_mode == "w")
-            header.name = header_name;
-        else if(!header.load(header_name))
-        {
-            if(file_mode == "r")
-            {
-                std::cout<<"ENVI Error: header file not found: "<<header_name<<std::endl;
-                exit(1);
-            }
-            //otherwise use the default header and save the name
-            header.name = header_name;
-        }
-
-		mode = file_mode;
-
-		//lock the file if it is read-only
-		if(mode == "r")
-            locked = true;
-		else if(mode == "a" || mode == "r+" || mode == "a+")
-		{
-            //if the file can be edited, lock it if data has already been written
-            if(test_exist(file_name))
-                locked = true;
+		if(!file)
+		{
+            return false;
 		}
-
-
-		//open the data file
-		std::string tmpmode = mode+"b";
-		//std::cout<<"Mode: "<<tmpmode<<std::endl;
-		data = fopen(file_name.c_str(), tmpmode.c_str());
-		if(data == NULL)
+
+		//the first line should just be "ENVI"
+		std::string line;
+		file>>line;
+		if(line != "ENVI")
 		{
-			std::cout<<"ENVI Error: unable to open binary file: "<<file_name<<std::endl;
+			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)
+		{
-	void setDescription(std::string desc)
-	{
-		readonly();		//if the file is read-only, throw an exception
-		header.description = desc;
-	}
-	void setZPlotTitles(std::string spectrum, std::string value)
-	{
-		readonly();		//if the file is read-only, throw an exception
-		header.z_plot_titles[0] = spectrum;
-		header.z_plot_titles[1] = value;
-	}
-	void setPixelSize(double sx, double sy)
-	{
-		readonly();		//if the file is read-only, throw an exception
-		header.pixel_size[0] = sx;
-		header.pixel_size[1] = sy;
-	}
+			//get the token
+			token = get_token(line);
-	void setWavelengthUnits(std::string units)
-	{
-		readonly();		//if the file is read-only, throw an exception
-		header.wavelength_units = units;
-	}
+			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];
+			}
-	void setCoordinates(double x, double y)
-	{
-		readonly();		//if the file is read-only, throw an exception
-		header.x_start = x;
-		header.y_start = y;
+			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
+		unsigned int wavelengths = wavelength.size();
+		if(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);
+        }
+
+		//close the file
+		file.close();
+
+		//set the file name
+		name = filename;
+
+		return true;
 	}
-	//FUNCTIONS THAT MAY BE DISALLOWED IN A LOCKED FILE
-	void setSamples(unsigned int samples)
+	void save(std::string filename)
 	{
-		EnviHeader newHeader = header;
-		newHeader.samples = samples;
-		set_header(newHeader);
+		//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;
+		}
+		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 setLines(unsigned int lines)
-	{
-		EnviHeader newHeader = header;
-		newHeader.lines = lines;
-		set_header(newHeader);
-	}
-
-	unsigned int lines() {return header.lines;}
-	unsigned int samples() {return header.samples;}
-	unsigned int bands() {return header.bands;}
-	double wavelength(unsigned int b) {return header.wavelength[b];}
-
-	void setBands(unsigned int bands)
+	void save()
 	{
-		EnviHeader newHeader = header;
-		newHeader.bands = bands;
-		set_header(newHeader);
-	}
-
-	unsigned int getElementSize()
-	{
-        //return the size of each element in bytes
-        return header.valsize();
-    }
-    unsigned int getBandSize()
-    {
-        //returns the size of a band in bytes
-        return header.valsize() * header.lines * header.samples;
-    }
-
-    Envi::dataType getDataType()
-    {
-        return header.data_type;
-    }
-
-	//DATA MANIPULATION
-	void addBand(void* band, unsigned int samples, unsigned int lines, double wavelength, std::string band_name ="")
-	{
-		//add a band to the file
-
-		EnviHeader newHeader = header;
-		newHeader.bands++;
-		newHeader.samples = samples;
-		newHeader.lines = lines;
-		newHeader.wavelength.push_back(wavelength);
-		if(band_name != "")
-			newHeader.band_names.push_back(band_name);
-		set_header(newHeader);
-
-		//copy the data to the file
-		fwrite(band, header.valsize(), samples * lines, data);
-
-		//since data was written, the file must be locked
-		locked = true;
-	}
-
-	//DATA RETRIEVAL
-	void getBand(void* ptr, unsigned int n)
-	{
-        //reads the n'th band and writes it to the given pointer
-        //  memory has been assumed to be allocated
-
-        if(n > header.wavelength.size())
-        {
-            std::cout<<"ENVI Error: Invalid band number."<<std::endl;
-            return;
-        }
-
-
-        if(header.interleave == Envi::BSQ)
-        {
-
-            //compute the position of the desired band (relative to the beginning of the file)
-            int bandpos = header.header_offset + getBandSize() * n;
-
-            //seek to that position
-            int seek_result = fseek(data, bandpos, SEEK_SET);
-            if(seek_result)
-            {
-                std::cout<<"ENVI Error: Error seeking through data file."<<std::endl;
-                return;
-            }
-
-            //read the band data from the file
-            size_t n_read = fread(ptr, getElementSize(), header.samples * header.lines, data);
-            if(n_read != header.samples * header.lines)
-            {
-                std::cout<<"ENVI Error -- Error reading band data: "<<n_read<<"/"<<samples() * lines()<<" pixels read"<<std::endl;
-                exit(1);
-            }
-
-            return;
-        }
-        else
-        {
-            std::cout<<"ENVI Error: Only BSQ operations are currently supported."<<std::endl;
-            return;
-        }
-    }
-
-    void getWavelength(void* band, double wavelength)
-    {
-        //this function retrieves a band via interpolation of wavelength values
-
-        //get the low and high band numbers
-        unsigned int low, high;
-        get_surrounding_bands(wavelength, low, high);
-        //std::cout<<"High: "<<high<<"   Low: "<<low<<std::endl;
-
-		//calculate the position between the two bands
-		double a;
-		if(low == high)
-			a = 1.0;
-		else
-			a = (wavelength - header.wavelength[low]) / (header.wavelength[high] - header.wavelength[low]);
-
-        //read both bands
-        void* A;
-        A = malloc( getBandSize() );
-        getBand(A, low);
-
-        void* B;
-        B = malloc( getBandSize() );
-        getBand(B, high);
-
-        //interpolate between the bands
-        interpolate_band(band, A, B, a);
-
-        //free the surrounding bands
-        free(A);
-        free(B);
-    }
-
-	//close file
-	void close()
-	{
-        //std::cout<<"ENVI File Mode: "<<mode<<std::endl;
-
-		//close the data stream
-		fclose(data);
-
-		//close the header
-		if(mode != "r")
-			header.save();
+        //std::cout<<"ENVI Header Name: "<<name<<std::endl;
+		save(name);
 	}
-};		//end EnviFile
-
+	//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;
+		}
-#endif
+	}
+};		//end EnviHeader
+}
+#endif
@@ -55,7 +55,7 @@ public:
 		int index;				//result of the max operation
 		rts::complex<T> result;
-		if(sizeof(T) == 4)
+		if(sizeof(T) == 8)
 			stat = cublasIcamax(handle, L, (const cuComplex*)X[n], 1, &index);
 		else
 			stat = cublasIzamax(handle, L, (const cuDoubleComplex*)X[n], 1, &index);
@@ -114,7 +114,7 @@ public:
 	void toImage(std::string filename, attribute type = magnitude, unsigned int n=0){
-		realfield<T, 1, true> rf(R[0], R[1]);
+		field<T, 1> rf(R[0], R[1]);
 		//get cuda parameters
 		dim3 blocks, grids;
@@ -122,7 +122,7 @@ public:
 		if(type == magnitude){
 			gpu_complexfield_mag <<<grids, blocks>>> (rf.ptr(), X[n], R[0], R[1]);
-			rf.toImages(filename, 0);
+			rf.toImage(filename, n, true);
 		}
 	}
@@ -134,4 +134,4 @@ public:
 }	//end namespace rts
-#endif
 \ No newline at end of file
+#endif
@@ -5,10 +5,14 @@
 #include <string>
 #include <sstream>
+#include "cublas_v2.h"
+#include <cuda_runtime.h>
+
 #include "../math/rect.h"
 #include "../cuda/threads.h"
 #include "../cuda/error.h"
 #include "../cuda/devices.h"
+#include "../visualization/colormap.h"
 namespace rts{
@@ -86,6 +90,41 @@ protected:
 		grids = dim3((R[0] + SQRT_BLOCK -1)/SQRT_BLOCK, (R[1] + SQRT_BLOCK - 1)/SQRT_BLOCK);
 	}
+	//find the maximum value of component n
+	T find_max(unsigned int n){
+		cublasStatus_t stat;
+		cublasHandle_t handle;
+
+		//create a CUBLAS handle
+		stat = cublasCreate(&handle);
+		if(stat != CUBLAS_STATUS_SUCCESS){
+			std::cout<<"CUBLAS Error: initialization failed"<<std::endl;
+			exit(1);
+		}
+
+		int L = R[0] * R[1];    //compute the number of discrete points in a slice
+		int index;				//result of the max operation
+		T result;
+
+		if(sizeof(T) == 4)
+			stat = cublasIsamax(handle, L, (const float*)X[n], 1, &index);
+		else
+			stat = cublasIdamax(handle, L, (const double*)X[n], 1, &index);
+
+		index -= 1;        //adjust for 1-based indexing
+
+		//if there was a GPU error, terminate
+		if(stat != CUBLAS_STATUS_SUCCESS){
+			std::cout<<"CUBLAS Error: failure finding maximum value."<<std::endl;
+			exit(1);
+		}
+
+		//retrieve the maximum value for this slice and store it in the maxVal array
+		std::cout<<X[n]<<std::endl;
+		HANDLE_ERROR(cudaMemcpy(&result, X[n] + index, sizeof(T), cudaMemcpyDeviceToHost));
+		return result;
+	}
+
 public:
 	//returns a list of file names given an input string with wild cards
@@ -286,10 +325,20 @@ public:
         	rts::gpu_field_crop <<<dimGrid, dimBlock>>> (result.X[n], X[n], R[0], R[1], width, height);
         return result;
+    }
+
+    //save an image representing component n
+    void toImage(std::string filename, unsigned int n = 0,
+    			 bool positive = false, rts::colormapType cmap = rts::cmBrewer){
+    	T max_val = find_max(n);	//find the maximum value
+    	if(positive)				//if the field is positive, use the range [0 max_val]
+    		rts::gpu2image<T>(X[n], filename, R[0], R[1], 0, max_val, cmap);
+    	else
+    		rts::gpu2image<T>(X[n], filename, R[0], R[1], -max_val, max_val, cmap);
     }
 };
 }		//end namespace rts
-#endif
 \ No newline at end of file
+#endif
@@ -2,7 +2,6 @@
 #define RTS_FUNCTION_H
 #include <string>
-#include <cmath>
 namespace rts{
@@ -59,7 +58,7 @@ protected:
 public:
 	function(){
-		range[0] = range[1] = NAN;
+		range[0] = range[1] = 0;
 		bounding[0] = bounding[1] = 0;
 	}
@@ -182,7 +181,8 @@ public:
 	function<Tx, Ty> & operator= (const Ty & rhs){		
 		X.clear();
 		Y.clear();
-		if(rhs != 0)			//if the RHS is zero, just clear, otherwise add one value of RHS
+		bounding[0] = bounding[1] = rhs;	//set the boundary values to rhs
+		if(rhs != 0)						//if the RHS is zero, just clear, otherwise add one value of RHS
 			insert(0, rhs);
 		return *this;
 #include "../optics/material.h"
 #include "../math/complexfield.cuh"
 #include "../math/constants.h"
+#include "../envi/bil.h"
 #include "cufft.h"
@@ -163,6 +164,8 @@ private:
 	T NA[2];						//numerical aperature (central obscuration and outer diameter)
+	function<T, T> source_profile;	//profile (spectrum) of the source (expressed in inverse centimeters)
+
 	complexfield<T, 1> scratch;		//scratch GPU memory used to build samples, transforms, etc.
 	void fft(int direction = CUFFT_FORWARD){
@@ -251,7 +254,7 @@ private:
 			HANDLE_ERROR(cudaMemcpy( gpuRi + inu, &ri, sizeof(complex<T>), cudaMemcpyHostToDevice ));
 			//save the source profile to the GPU
-			source = 1;
+			source = source_profile(10000 / lambdas[inu]);
 			HANDLE_ERROR(cudaMemcpy( gpuSrc + inu, &source, sizeof(T), cudaMemcpyHostToDevice ));
 		}
@@ -311,6 +314,11 @@ private:
 		scratch = scratch.crop(Z, S);
 	}
+	//save the scratch field as a binary file
+	void to_binary(std::string filename){
+
+	}
+
 public:
@@ -322,6 +330,7 @@ public:
 		NA[0] = 0;
 		NA[1] = 0.8;
 		S = 0;
+		source_profile = 1;
 	}
 	//add a layer, thickness = microns
@@ -334,6 +343,11 @@ public:
 		add_layer(material<T>(filename), thickness);
 	}
+	//adds a profile spectrum for the light source
+	void set_source(std::string filename){
+		source_profile.load(filename);
+	}
+
 	//adds a block of wavenumbers (cm^-1) to the simulation parameters
 	void add_wavenumbers(unsigned int start, unsigned int stop, unsigned int step){
 		unsigned int nu = start;
@@ -368,6 +382,10 @@ public:
 		na(0, out);
 	}
+	rts::function<T, T> get_source(){
+		return source_profile;
+	}
+
 	void save_sample(std::string filename){
 		//create a sample and save the magnitude as an image
 		build_sample();
@@ -375,7 +393,7 @@ public:
 		scratch.toImage(filename, rts::complexfield<T, 1>::magnitude);
 	}
-	void save_mirst(std::string filename){
+	void save_mirst(std::string filename, bool binary = true){
 		//apply the MIRST filter to a sample and save the image
 		//build the sample in the Fourier domain
@@ -390,7 +408,10 @@ public:
 		crop();
 		//save the image
-		scratch.toImage(filename, rts::complexfield<T, 1>::magnitude);
+		if(binary)
+			to_binary(filename);
+		else
+			scratch.toImage(filename, rts::complexfield<T, 1>::magnitude);
 	}
@@ -411,6 +432,9 @@ public:
 		for(unsigned int l=0; l<layers.size(); l++)
 			ss<<"layer "<<l<<": "<<layers[l]<<" um"<<"---------"<<std::endl<<matlist[l].str()<<std::endl;
+		ss<<"source profile-----------"<<std::endl;
+		ss<<get_source().str()<<std::endl;
+
 		return ss.str();