release / siproc

#include "stim/parser/filename.h"
#include "stim/envi/envi.h"
#include "stim/visualization/colormap.h"
#include "stim/envi/agilent_binary.h"

#include <thread>
void progress_thread_double(double* e);		//progress bar threaded function



void mosaic_coord(unsigned int &x, unsigned int &y, stim::filename dmd_filename){

	//get the file prefix
	std::string prefix = dmd_filename.get_prefix();

	//these files are in the format:
	//		???????????_xxxxx_yyyyy.dmd

	//first find the x coordinate

	//find the position of the last underscore
	size_t y_start = prefix.find_last_of("_");
	//remove all values between x_start and the end of the string
	std::string y_string = prefix.substr(y_start + 1);

	y = atoi(y_string.c_str());

	//crop out the x coordinate and the last underscore
	prefix = prefix.substr(0, y_start);

	//find the y coordinate using the same method as above
	size_t x_start = prefix.find_last_of("_");
	std::string x_string = prefix.substr(x_start + 1);
	x = atoi(x_string.c_str());
}

//Use the file list to find the number of tiles along each dimension of the mosaic
void mosaic_size(unsigned int &sx, unsigned int &sy, std::vector<stim::filename> dmd_list){
	sx = sy = 0;

	unsigned int xi, yi;

	//for each dmd file
	for(unsigned int i = 0; i < dmd_list.size(); i++){

		//find the coordinate of the tile represented by this file
		mosaic_coord(xi, yi, dmd_list[i]);

		//if the coordinate is larger than the current max, store it as the current max
		if(xi > sx) sx = xi;
		if(yi > sy) sy = yi;
	}

	//add 1 to reflect that the indexing starts at 0
	sx++; sy++;
}

void mosaic_agilent(std::string filemask, std::string outfile, size_t &X, size_t &Y, size_t &B){
	
	double progress = 0;
	std::thread t1(progress_thread_double, &progress);		//start the progress bar thread

	//generate a list of tile files
	stim::filename f = filemask;								//create a filename from the specified file mask
	if(!f.wildcards())											//if the file name doesn't contain wildcards
		f = filemask + "/*.dmd";								//default to files of type DMD (generally used as agilent mosaic tiles)	}
		
	std::vector<stim::filename> tile_list = f.get_list();		//create a list of tile files

	stim::agilent_binary<float> binfile(tile_list[0]);			//load the first agilent binary file
	size_t detector_size = binfile.dim(0);					//get the detector size

	std::cout<<"Building a mosaic with a detector size of "<<detector_size<<std::endl;
	
	//calculate the mosaic size
	unsigned int mx, my;
	mosaic_size(mx, my, tile_list);

	std::cout<<"Mosaic size: "<<mx<<"  "<<my<<std::endl;

	//if(create_header){											//if a header has to be created
	std::ifstream in(tile_list[0].str().c_str(), std::ifstream::ate | std::ifstream::binary);	//open the first dmd
	size_t file_size = in.tellg();							//get the size of the file
	B = (file_size - 1020) / (sizeof(float) * detector_size * detector_size);		//estimate the number of bands by dividing by the frame size
	//}
	//else{	
		//f = f.path() + "/*.hdr";								//create a mask for the header file (the name is unknown)
		//std::vector<stim::filename> hdr_list = f.get_list();	//allocate space for a list of headers that match this pattern
		//if(hdr_list.size() == 0){ 								//throw an error of a header doesn't exist
		//	std::cout<<"No ENVI header file found. A header file created with Resolutions Pro (export ENVI) provides useful wavelength information. If one is unavailable, use the --create-header option to estimate these properties and create a minimal header file."<<std::endl;
		//	exit(1);
		//}
		//std::cout<<"First header name: "<<hdr_list[0].str()<<std::endl;
		//open the representative ENVI header file (this should be the header for the thumbnail produced by Resolutions Pro)
		//header.load(hdr_list[0].str());
	//}

	//change the header dimensions to reflect the detector size
	//	this corresponds to the size of the tile along X and Y
	stim::envi_header header;									//allocate space for an ENVI header
	header.samples = detector_size;
	header.lines = detector_size;
	header.bands = B;

	header.header_offset = 1020;							//Agilent binary files have a 1020 byte header

	//create a list of ENVI file structures
	std::vector<stim::envi> envi_list;
	envi_list.resize(mx * my);				//set the length to the number of tiles

	//open each tile
	for(unsigned int m = 0; m < mx * my; m++){
		envi_list[m].open(tile_list[m].str(), header);
		envi_list[m].close();
	}

	//calculate the size of a single band of the mosaic
	X = mx * detector_size;				//number of samples
	Y = my * detector_size;				//number of lines
	size_t type_size = envi_list[0].type_size();	//byte size of the data type
	size_t tile_width = detector_size * type_size;	//byte size of a tile
	size_t band_width = X * type_size;			//byte width of the band image
	size_t band_size = band_width * Y;			//size of the band in bytes


	//allocate space for a single band of the mosaic
	char* band = (char*)malloc(band_size);

	//allocate space for a single tile of the mosaic
	char* tile_band = (char*)malloc(detector_size * detector_size * type_size);

	//open a binary file for writing
	std::ofstream outstream(outfile.c_str(), std::ios::out|std::ios::binary);

	//***********Combine Tiles into a Band Image**********************

	//for each band
	for(size_t b = 0; b < B; b++){
	//unsigned int b = 0;
		//for each tile along y
		for(size_t yt = 0; yt < my; yt++){
			//for each tile along x
			for(size_t xt = 0; xt < mx; xt++){

				envi_list[xt * my + yt].open();								//open the file for reading (windows has a max # of open files)	
				envi_list[xt * my + yt].band_index((void*)tile_band, b);	//load the tile image
				envi_list[xt * my + yt].close();							//close the file

				size_t x_idx = xt * tile_width;						//calculate the position of this line along x in the output image

				//for each line along y
				for(size_t y = 0; y < detector_size; y++){										
					size_t y_idx = ((my-1) - yt) * detector_size + y;			//calculate the position of this line along y

					//calculate the position of this line in the output image
					size_t dst_i = y_idx * band_width + x_idx;

					size_t src_i = y * tile_width;

					memcpy(band + dst_i, tile_band + src_i, tile_width);
				}

			}

		}
		outstream.write(band, band_size);
		progress = (double)(b + 1) / header.bands * 100;
	}

	progress = 100;

	//close the binary file
	outstream.close();

	//free the memory allocated for the tile and output bands
	free(tile_band);
	free(band);

	//update and save the header
	//header.samples = X;
	//header.lines = Y;
	//header.header_offset = 0;			//the final output doesn't have an offset
	//header.save(outfile + ".hdr");

	t1.join();									//wait for the progress bar thread to finish (it probably already is)
}

void mosaic_agilent_spectrum(){

}

void mosaic_agilent_interferogram(std::string filemask, std::string outfile, double ELWN, int UDR){

	size_t X, Y, B;												//store the size of the mosaic
	mosaic_agilent(filemask, outfile, X, Y, B);					//generate the mosaic

	//create the header
	stim::envi_header header;
	header.samples = X;
	header.lines = Y;
	header.bands = B;
	
	double int_delta = (1.0 / ELWN) * ((double)UDR / 2.0);		//calculate the interferogram spacing

	header.set_wavelengths(0, int_delta);
	header.wavelength_units = "centimeters";
	header.save(outfile + ".hdr");

}