release / siproc

//#define _CRTDBG_MAP_ALLOC  
//#include <stdlib.h>  
//#include <crtdbg.h>

#include <iostream>

#include <stim/parser/arguments.h>
#include <stim/envi/agilent_binary.h>
#include <stim/parser/filename.h>

stim::arglist args;								//generate a list of accepted arguments

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

bool fft = false;								//perform an FFT before storing
double wn_low, wn_high;
float elwn;
int udr;

bool absorbance = false;
stim::agilent_binary<float> background;						//store the background data

bool crop_samples = false;
size_t n_samples;


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

	//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 save_mosaic(std::vector<stim::filename> tile_list, std::string outfile, stim::agilent_binary<float>* background = NULL) {

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

	unsigned int mx, my;
	mosaic_size(mx, my, tile_list);								//calculate the mosaic size
	std::cout << "Building a ["<<mx<<" x "<<my<<"] mosaic with a detector size of (" << ds << " x " << ds <<") pixels"<< std::endl;
	std::cout << "Mosaic size: " << mx << "  " << my << std::endl;		//output the mosaic size

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

	size_t nxo = mx * ds;										//number of output pixels along x
	size_t nyo = my * ds;										//number of output pixels along y
	size_t nbo;

	size_t bytes_tile_band = ds * ds * sizeof(float);			//number of bytes in a tile band
	size_t bytes_out_band = bytes_tile_band * mx * my;			//number of bytes in an output image band

	double progress = 0.0;
	std::thread t1(progress_thread_double, &progress);		//start the progress bar thread
	size_t xo, yo, io;
	size_t ii;
	stim::envi_header header;
	bool header_flag = true;
	for (size_t ty = 0; ty < my; ty++) {
		for (size_t tx = 0; tx < mx; tx++) {
			xo = tx * ds;															//calculate the x output position for the current tile
			stim::agilent_binary<float> I(tile_list[tx * my + ty].str());
			if (crop_samples) I.crop(n_samples);									//crop samples if necessary to reduce spectral resolution
			if (fft) I = I.fft(wn_low, wn_high, elwn, udr);							//if an FFT is requested, perform it
			if (background)	I.absorbance(background);								//if the user provides a background file
			if (header_flag) {
				header = I.create_header();
				header_flag = false;
			}
			nbo = I.dim(2);															//get the number of bands in the output image
			for (size_t b = 0; b < nbo; b++) {										//for each band in the input
				for (size_t yi = 0; yi < ds; yi++) {								//for each line in the input image
					yo = ty * ds + ds - yi - 1;										//calculate the y output position
					io = (b * nxo * nyo + yo * nxo + xo);							//calculate the output position as a 1d index
					ii = (b * ds * ds + yi * ds);									//calculate the input position as a 1d index
					out.seekp(io * sizeof(float));									//seek to the correct output location
					out.write(I.data() + ii * sizeof(float), ds * sizeof(float));	//copy the line from input to output
				}
			}
			progress = (double)(ty * mx + tx + 1) / (double)(mx * my) * 100;
		}
	}
	out.close();
	t1.join();									//wait for the progress bar thread to finish (it probably already is)

	header.samples = mx * ds;
	header.lines = my * ds;
	header.save(outfile + ".hdr");
}

void output_files(std::vector<stim::filename> infiles, std::string outfile) {
	if (infiles.size() == 1) {
		std::cout << "source file: " << infiles[0].str() << std::endl;
	}
	else {
		std::cout << "number of input files: " << infiles.size() << std::endl;
	}
	std::cout << "destination file: " << outfile << std::endl;
}

int main(int argc, char* argv[]){
	
	args.add("help", "print help");
	args.add("spacing", "specify the spectral resolution", "", "positive integer specifying wavenumber spacing");
	args.add("fft", "perform a fourier transform of the data before assembling the mosaic", "", "optional wavenumber range (override UDR), ex: 900 3500");
	args.add("udr", "udr filter used for imaging", "4", "(4 or 8 supported)");
	args.add("elwn", "ELWN tuning number", "15798.0039", "HeNe sampling rate (in cm^(-1))");
	args.add("background", "specify a background file for measuring absorbance", "", "filename.seq");
	args.add("cuda", "specify a CUDA device for the FFT", "0", "index of CUDA device (-1 forces CPU only)");
	args.parse(argc, argv);								//parse the argument list

	if (args.nargs() < 2 || args["help"]) {
		std::cout << "Usage: " << std::endl<<std::endl;
		std::cout << "\tGenerate an ENVI mosaic from a set of tiles produced by a Cary FTIR imaging system:" << std::endl;
		std::cout << "\t\t>> cary-ftir input_files/*.dmd output_mosaic" << std::endl << std::endl;
		std::cout << args.str() << std::endl;
		exit(1);
	}

	stim::filename infile(args.arg(0));
	std::vector<stim::filename> in;						//create a list of input files
	if (args.nargs() == 2) {
		if (infile.wildcards()) in = infile.get_list();		//if a wildcard is specified, get the file list
		else					in.push_back(infile);		//if a single file is specified, create a file list of 1
	}
	//if the user specifies several files on the command line
	else {
		for (size_t f = 0; f < args.nargs() - 1; f++) {		//for each file specified
			in.push_back(args.arg(f));						//add the file to the file list
		}
	}

	std::string outfilename = args.arg(args.nargs() - 1);	//get the output file

	output_files(in, outfilename);
	
	udr = args["udr"].as_int();							//save the UDR filter type
	if (udr == 4) {
		wn_low = 900;
		wn_high = 3900;
	}
	else {
		wn_low = 800;
		wn_high = 1900;
	}

	if (args["fft"].is_set()) {
		fft = true;										//set the fft flag to true
		
		if (args["fft"].nargs() >= 1) {
			if(args["fft"].is_num(0))
				wn_low = args["fft"].as_float(0);
			else {
				std::cout << "ERROR: FFT output range must be a number: " << args["fft"].as_string(0) << std::endl;
				return 1;
			}
		}
		if (args["fft"].nargs() >= 2) {
			if (args["fft"].is_num(1))
				wn_high = args["fft"].as_float(1);
			else {
				std::cout << "ERROR: FFT output range must be a number: " << args["fft"].as_string(1) << std::endl;
				return 1;
			}
		}
			
	}

	elwn = (float)args["elwn"].as_float();

	if (args["spacing"]) {
		double di = (1.0 / elwn) * ((double)udr / 2.0);						//calculate the interferogram spacing
		double df = (double)args["spacing"].as_int();						//get the desired spectral spacing
		n_samples = (size_t)std::ceil(1.0 / (di * df));						//calculate the number of samples required
		std::cout << "Override default spectral spacing with " << df << " cm^(-1) by cropping to " << n_samples << " samples" << std::endl;
		crop_samples = true;
	}

	if (args["background"]) {
		if (!args["fft"]) {
			std::cout << "ERROR: using a background requires calculating an FFT - use the --fft option as well" << std::endl;
			exit(1);
		}
		stim::agilent_binary<float> background(args["background"].as_string());
		if (background) {
			if (crop_samples) background.crop(n_samples);					//crop the samples if necessary
			background = background.fft(wn_low, wn_high, elwn, udr);
			save_mosaic(in, outfilename, &background);
			//absorbance = true;
		}
		else {
			std::cout << "ERROR loading background file: " << args["background"].as_string() << std::endl;
			return 1;
		}
	}
	else {
		save_mosaic(in, outfilename);
	}
	
//	_CrtDumpMemoryLeaks();
}