//#include "lapacke.h"
//#include "cblas.h"

//#include <string>
#include "linalg.h"

#include <iostream>
#include <fstream>
#include <random>
#include <algorithm>

#include "stim/parser/arguments.h"
#include "stim/envi/envi_header.h"
#include "stim/envi/envi.h"
#include "stim/visualization/colormap.h"
#include <stim/image/image.h>
#include <stim/parser/table.h>
#include <stim/parser/filename.h>
#include <stim/envi/agilent_binary.h>
#include <stim/math/matrix.h>
#include <time.h>
#include <thread>

//LAPACKE support for Visual Studio
//#include <complex>
//#ifndef LAPACK_COMPLEX_CUSTOM
//#define LAPACK_COMPLEX_CUSTOM
//#define lapack_complex_float std::complex<float>
//#define lapack_complex_double std::complex<double>
//#endif


void baseline(std::string infile, std::string outfile, std::string headerfile, std::vector<double> points, unsigned char* mask);
void normalize(std::string infile, std::string outfile, std::string headerfile, double band, unsigned char* mask);
void convert(std::string infile, std::string outfile, std::string headerfile, stim::envi_header::interleaveType type);
std::vector< std::vector<double> > read_metric_list(std::string filename);
void create_metric_file(std::vector< std::vector<double> > metrics, std::string outfile);
stim::colormapType str2cmap(std::string str);
void bandimage(std::string infile, std::string headername, double band, std::string filename, stim::colormapType cmap, unsigned char* MASK);
void bandimage(std::string infile, std::string headername, double band, std::string filename, double minVal, double maxVal, stim::colormapType cmap, unsigned char* MASK);
//void mosaic_agilent(std::string directory, std::string outfile, bool create_header);
void mosaic_agilent_interferogram(std::string filemask, std::string outfile, double ELWN, int UDR);
void mosaic_spero(std::string directory, std::string outfile);
void mnf(std::string outfile, int keptComponents, std::string NoiseFractions, int cuda_device);

void progress_thread_envi(stim::envi* e);		//progress bar threaded function
void progress_thread_double(double* e);		//progress bar threaded function

//CUDA externs
//void gpu_bsq2bip();


# define MAX_CLASSES 20
# define per 20000

unsigned char* MASK = NULL;					//pointer to an input file mask if one is provided
stim::envi ENVI;								//input ENVI file to be processed
unsigned long long X, Y, B;						//registers to quickly access the image size

bool progressbar = true;
bool verbose = false;
bool optimization = true;

void advertisement(){
	std::cout<<std::endl<<std::endl;
	std::cout<<"========================================================================="<<std::endl;
	std::cout<<"Thank you for using the HSIPROC spectroscopic image processing toolkit!"<<std::endl;
	std::cout<<"Scalable Tissue Imaging and Modeling (STIM) Lab, University of Houston"<<std::endl;
	std::cout<<"Developers: Sam Saki, Ziqi He, David Mayerich, Brad Deutsch"<<std::endl;
	std::cout<<"Source: https://github.com/stimlab/hsiproc.git"<<std::endl;
	std::cout << "This version has been compiled on " << __DATE__ << " at " << __TIME__ << std::endl;
	std::cout<<"========================================================================="<<std::endl<<std::endl;
}

int test(int i){
	std::cout<<"Thread started..."<<std::endl;
	return i * 10;
}

void set_mask(std::string filename, size_t N = 0){
	size_t bytes = X * Y * sizeof(unsigned char);					//calculate the number of bytes in the mask
	MASK = (unsigned char*) malloc(bytes);							//allocate space for the mask
	memset(MASK, 0, bytes);											//initialize the mask to 0 (zeros)
	std::string maskfile = filename;				//mask file name
	stim::image<unsigned char> image(maskfile);						//create an image from the file
	stim::image<unsigned char> mask = image.channel(0);				//grab the first channel of the image
	if (mask.width() != X || mask.height() != Y) {
		std::cout << "ERROR - mask size doesn't match image size: mask = [" << mask.width() << " x " << mask.height() << "], input = [" << X << " x " << Y << "]" << std::endl;
		exit(1);
	}
	if(N == 0){														//if only the mask name is given, use all pixels in the mask
		for(size_t xy = 0; xy < X*Y; xy++)								//for each pixel in the image
			if(mask.data()[xy]) MASK[xy] = mask.data()[xy];				//copy it to the MASK pointer
	}
	else{																//if a number of values is also specified
		N = (std::min)(N, mask.nnz());									//calculate the minimum (this odd notation is used to defeat a macro implemented in Visual Studio)
		std::vector<size_t> idx = image.sparse_idx();					//get the indices for the nonzero values in the image
		std::random_device rd;											//create a random number generation function
		std::mt19937 g(rd());
		std::shuffle(idx.begin(), idx.end(), g);						//shuffle the index values
		for(size_t i = 0; i < N; i++){									//for each of the first n random pixels
			MASK[idx[i]] = mask.data()[idx[i]];									//store the pixel value from the image into the mask
		}
	}

}

int main(int argc, char** argv){

	//create an argument list
	stim::arglist args;

#ifdef _WIN32
	args.set_ansi(false);
#endif

	//basic arguments
	args.add("help", "prints this help");

	args.section("Binary File Manipulation");
		args.add("convert", "convert file to bip, bsq, or bil", "", "[bsq], [bip], or [bil]");
		args.add("sift", "create a matrix of masked spectra in lexicographic order", "", "[mask_filename]");
		args.add("unsift", "creates a 2D image from a matrix of pixels in lexicographic order", "", "[mask_filename]");
		args.add("crop", "crop part of original file, use '-' to specify the limit, spectral coordinates are given in wavelength by default", "", "[x nx y ny w nw]");
		args.add("subimages", "extract a subimage at each mask point", "", "[maskfile width height], if no height is given, height = width");
		args.add("trim", "remove bands", "", "[b0 b1 b2-b3], removes b0, b1, and all bands between b2 and b3");
		args.add("select", "selects a set of bands from the input image and uses them to create a new output image", "", "list of wavelengths or bands, ex: 1250 1650 3200 ...");
		args.add("combine", "combines two images by placing the second at a specified position relative to the first", "", "[filename px py], if (px, py) isn't provided, combining is done along Y");
		args.add("append", "appends the bands from an image to the input image", "", "[filename] - the input and append file must be the same size in X and Y");

	args.section("Arithmetic");
		args.add("multiply", "multiply values in the image by some number n", "", "[n] where n is any real value");
		args.add("add", "add values in the image by some number n", "", "any real value");

	args.section("Data Retrieval / Analysis");
		args.add("image", "wavelength and filename for a band image (.raw files are X x Y float32)", "", "[wavelength] or [wavelength min max]");
		args.add("colormap", "specify the colormap for an image", "brewer", "brewer, grey");
		args.add("spectrum", "saves a CSV spectrum", "", "[x y]");
		args.add("mean", "saves the mean spectrum and variance (sigma^2) as a CSV file - mask can (and should) be applied");
		args.add("median", "saves the median spectrum as a CSV file - mask can be applied, only works for BSQ files");
		args.add("covariance", "calculate the covariance matrix of the spectrum and saves it as a CSV file (BIP recommended)", "", "[mean.csv], optional CSV file name storing the mean spectrum");
		args.add("deriv", "approximate the d-th derivative at order n (default n=2) using finite differences", "", "[d n]");

	args.section("Convolution (assumes equally spaced points)");
		args.add("convolve", "convolve the spectra with a given set of coefficients", "", "[c1 c2 c3...cn]");
		args.add("sg", "applies a Savitzky-Golay filter of width w and order n, where w is odd and n < w/2", "", "[w n]");

	args.section("Hyperspectral Image Correction");
		args.add("baseline", "piecewise linear baseline correction", "", "[baseline-file.txt] or [p0 p1 p2 p3 ...]");
		args.add("mnf", "Generate a basis for MNF noise removal. Use --project to apply the basis projection to a noisy file", "", "[keptComponents noisefile] , ex. --mnf 7 , in most cases: 3 < keptComponents < 23");
		args.add("matcond", "Specify the threshold for the MNF transform condition number c - a warning will be given if c is below this threshold.", "0.01");
		args.add("normalize", "normalize spectra using vector normalization (or a band ratio of a band is specified)", "", "[band]");
	
	args.section("Masking");
		args.add("threshold", "create a mask that includes all pixel in band such that min < val < max", "", "[band min max]");
		args.add("mask-finite", "create a mask - all pixels without inf or NaN values are masked");
		args.add("apply-mask","apply a mask (or set of masks) to an image (any false value in a mask is set to 0)", "", "[image_1 image_2 image_3]");
		args.add("mask", "limit clustering to a masked region specified by an image file", "", "[filename n], where n (optional) limits the mask to n random samples of the image");

	args.section("Dimension Reduction");
		args.add("pca", "calculate and save the PCA basis transform (BIP recommended)");
		args.add("project", "perform a forward projection given a translation and set of basis vectors", "", "[stats_file #vectors], ex. PCA: [pca.sta 30]");
		args.add("inverse", "perform an inverse projection given a translation and set of basis vectors", "", "[stats_file #vectors]");
		args.add("metrics", "calculate metrics based on Bhargava et al. metric file", "", "[metric_file]\n\t\t\t"
															"metric format:\n\t\t\t"
															"peak height ratio: \t1 LB1 RB1 P1 LB2 RB2 P2 0 0\n\t\t\t"
															"peak area ratio: \t2 LB1 RB1 LP1 RP1 LB2 RB2 P2 0\n\t\t\t"
															"area to area ratio: \t3 LB1 RB1 LP1 RP1 LB2 RB2 LP2 RP2\n\t\t\t"
															"center of gravity: \t4 LB1 RB1 LP1 RP1 0 0 0 0");
		
	args.section("Hardware-Specific Processing");
		args.add("create-header", "create a basic header file to represent the mosaic");
	args.section("Debugging Parameters");
		args.add("verbose", "provide verbose debug messages and output");
		args.add("noprogress", "removes the progress bar from the output");
		args.add("nooptimization", "turns off batch optimization, max batch size is used");
		args.add("mempct", "percentage of available memory to use for processing", "40", "values over 90% are not recommended");
		args.add("memraw", "specify a raw amount of memory (in MB) to use for processing (caution)");
		args.add("cuda", "selects the default CUDA device used for calculations", "0", "integer ID, (-1 prevents CUDA use even if available)");

	//parse the command line arguments
	args.parse(argc, argv);
	if(args["help"].is_set()){				//display the help text if requested
		advertisement();
		std::cout<<std::endl<<"usage: hsiproc input output --option [A B C ...]"<<std::endl;
		std::cout<<std::endl<<std::endl
				  << "examples: siproc bsqfile bipfile --convert bip"<<std::endl;
		std::cout << "                      convert the input file bsqfile to the output bipfile, changing the interleave from BSQ to BIP" << std::endl;
		std::cout << "          siproc bsqfile bandfile.bmp --image 1650" << std::endl;
		std::cout << "                      create a bitmap image of bsqfile at wavelength 1650" << std::endl;
		std::cout << "          siproc mosaic_in mosaic_cropped --crop 1000 2000 - 1000 1500 120" << std::endl;
		std::cout << "                      crop the input file mosaic_in to generate a 1000 x 1000 image spanning wavelengths 1500 to 1620" << std::endl;
		std::cout<<std::endl<<std::endl;
		std::cout<<args.str()<<std::endl;
		exit(1);
	}

	////Address Input and Output file issues
	//In some cases, the input and output can be inferred which means that the user doesn't HAVE
	//	to specify 2 arguments (input and output file). Check for these cases here
	std::string infile;
	std::string outfile;

	if(args.nargs() == 0){									//if no arguments are provided
		advertisement();
		std::cout<<"ERROR: No input or output specified, enter hsiproc --help for options."<<std::endl<<std::endl;
		return 1;
	}
	else if(args.nargs() == 1){									//if only one argument is available
		std::cout<<"ERROR: missing input/output file parameter (2 are required for this algorithm)."<<std::endl;
		return 1;
	}
	else{													//two (or more) arguments are available (ignore more than two)
		infile = args.arg(0);
		outfile = args.arg(1);
	}

	
	
	//////////////// Debugging Parameters /////////////////////////////////////
	if(args["verbose"].is_set()){
		verbose = true;
		progressbar = false;
	}
	if(args["nooptimization"].is_set()){
		optimization = false;
	}

	std::string headername = infile + ".hdr";              	//guess the header file name
	
	size_t optcount = args.nopts();

	
	//open the input ENVI file
	if(stim::envi::is_envi(infile)){						//if the input file is an ENVI file, open it
		ENVI.open(infile, infile + ".hdr");
		if(!ENVI){
			std::cout<<"ERROR: failed to open input ENVI file: "<<infile<<std::endl;
			return 1;
		}
		X = ENVI.header.samples;
		Y = ENVI.header.lines;
		B = ENVI.header.bands;
		if(args["memraw"].is_set()) ENVI.set_buffer_raw((size_t)args["memraw"].as_int(0) * 1000000);
		else                         ENVI.set_buffer_frac(args["mempct"].as_float(0) / 100.0);				//set the buffer size for processing the ENVI file
	}
	else {
		std::cout << "ERROR: the input file does not appear to be an ENVI hyperspectral file: " << infile << std::endl;
		return 1;
	}

	//display the input and output files
	std::cout<<"input:  "<<infile<<std::endl;
	std::cout<<"output: "<<outfile<<std::endl;

	time_t t1=time(NULL);			//start time

	//open the mask, if set
	if(args["mask"].is_set()){											//if a mask file is specified

		if(args["mask"].nargs() == 1)
			set_mask(args["mask"].as_string(0));
		else if(args["mask"].nargs() == 2)
			set_mask(args["mask"].as_string(0), args["mask"].as_int(1));
	}
	
	if(args["threshold"].is_set()){

		double mask_band = args["threshold"].as_float(0);
		//double threshold = args["threshold"].as_float(1);
		double lower, upper;
		lower = args["threshold"].as_float(1);
		upper = std::numeric_limits<double>::max();				//set the default maximum value to the maximum possible double value
		if(args["threshold"].nargs() == 3){
			upper = args["threshold"].as_float(2);
		}

		unsigned long long N = X * Y;

		//get the mask size (samples * lines)
		unsigned char* mask = (unsigned char*)malloc(N * sizeof(unsigned char));			//allocate memory for the mask
		std::thread t1(progress_thread_envi, &ENVI);		//start the progress bar thread
		std::cout<<"Creating a mask at band "<<mask_band<<" with threshold ["<<lower<<", "<<upper<<"]..."<<std::endl;
		ENVI.build_mask(mask, mask_band, lower, upper, MASK, true);	//get the mask
		t1.join();

		
		stim::image<unsigned char> mask_image(mask, X, Y);
		mask_image.save(outfile.c_str());								//save the mask image
		
		//count the number of pixels that are masked
		unsigned long long P = 0;
		for(unsigned long long i = 0; i < N; i++)
			if(mask[i] != 0) P++;

		std::cout<<P<<" out of "<<N<<" pixels masked ("<<(double)P/N * 100<<"%)"<<std::endl;
	}
	else if( args["mask-finite"].is_set() ){
		unsigned long long N = X * Y;

		std::cout<<"Creating a mask for finite values..."<<std::endl;
		std::thread t1(progress_thread_envi, &ENVI);		//start the progress bar thread
		unsigned char* mask = (unsigned char*)malloc(N * sizeof(unsigned char));
		ENVI.mask_finite(mask, MASK, true);
		t1.join();

		stim::image<unsigned char> mask_image(mask, X, Y);
		mask_image.save(outfile.c_str());								//save the mask image
		
		//count the number of pixels that are masked
		unsigned long long P = 0;
		for(unsigned long long i = 0; i < N; i++)
			if(mask[i] != 0) P++;

		std::cout<<N-P<<" pixels have been marked as not finite ("<<(double)(N-P)/(double)N * 100<<"%)"<<std::endl;
		//std::cout<<P<<" out of "<<N<<" pixels masked ("<<(double)P/N * 100<<"%)"<<std::endl;
	}

	//if baseline input is provided
	else if(args["baseline"].is_set()){

		//if only a single argument is provided, assume it is a baseline file
		if(args["baseline"].nargs() == 1){
			std::vector<double> blpts;		//create a vector for the baseline points

			std::ifstream blpt_file;			//open the file containing the baseline points
			blpt_file.open(args["baseline"].as_string().c_str());

			//get each baseline point and push it into the vector
			double pt;
			while(blpt_file>>pt){
				blpts.push_back(pt);
			}

			baseline(infile, outfile, headername, blpts, MASK);


		}
		//otherwise assume the points are provided individually
		else{

			std::vector<double> blpts;		//create a vector for the baseline points

			size_t npts = args["baseline"].nargs();	//get the number of baseline points
			for(unsigned p = 0; p < npts; p++){			//insert each point into the baseline point vector
				blpts.push_back(args["baseline"].as_float(p));
			}

			baseline(infile, outfile, headername, blpts, MASK);
		}

	}

	else if(args["normalize"].is_set()){

		//if two parameters are provided for normalization, the second is a mask value
		//	in that case, normalize to a temporary file and then apply a mask

		if(args["normalize"].nargs() == 0){
			std::cout<<"Calculating vector norm..."<<std::endl;
			std::thread t1(progress_thread_envi, &ENVI);		//start the progress bar thread
			ENVI.normalize(outfile, MASK, true);		//perform normalization
			t1.join();									//wait for the progress bar thread to finish (it probably already is)
		}
		else{
			double normband = args["normalize"].as_float(0);
			std::cout<<"Normalizing to band "<<normband<<"..."<<std::endl;
			std::thread t1(progress_thread_envi, &ENVI);		//start the progress bar thread
			ENVI.ratio(outfile, normband, MASK, true);		//perform normalization
			t1.join();									//wait for the progress bar thread to finish (it probably already is)
		}

	}
	else if (args["select"]) {										//if the user specifies the --select option (use bands to create a new file)
		size_t nb = args["select"].nargs();							//get the number of bands that the user wants to select
		std::vector<double> bandlist(nb);							//allocate an array to store the bands selected from the input file
		for (size_t b = 0; b < nb; b++)								//for each band given by the user
			bandlist[b] = args["select"].as_float(b);				//store that band in the bandlist array
		std::thread t1(progress_thread_envi, &ENVI);		//start the progress bar thread
		ENVI.select(outfile, bandlist, MASK, true);					//call an ENVI function to
		t1.join();									//wait for the progress bar thread to finish (it probably already is)

	}

	else if(args["apply-mask"].is_set()){//can mostly copy/paste this code
		
		MASK = (unsigned char*) malloc(X * Y * sizeof(unsigned char));	//allocate space for the mask
		memset(MASK, 255, X*Y*sizeof(unsigned char));
		for(size_t m = 0; m < args["apply-mask"].nargs(); m++){
			std::string maskfile = args["apply-mask"].as_string(m);			//mask file name
			stim::image<unsigned char> mask_image(maskfile);
			for(size_t xy = 0; xy < X*Y; xy++)
				if(!mask_image.data()[xy]) MASK[xy] = 0;
		}
		
		//run the function to apply the mask
		//	this outputs a new binary file called 'outfile' with the mask applied
		std::cout<<"Applying mask(s) to image..."<<std::endl;
		std::thread t1(progress_thread_envi, &ENVI);		//start the progress bar thread
		ENVI.apply_mask(outfile, MASK, true);
		t1.join();									//wait for the progress bar thread to finish (it probably already is)
	}

	else if (args["sift"].is_set()){

		//get the name of the mask file (an image file) as a string
		std::string maskname = args["sift"].as_string();

		stim::image<unsigned char> mask(maskname);
		
		//run the function to apply the mask and save sifted spectra
		//	this outputs a new binary file called 'outfile' (and associated header) with the mask applied.

		std::cout<<"Sifting to a [ P x B ] = [ "<<mask.nnz()<<" x "<<B<<"] matrix..."<<std::endl;
		std::thread t1(progress_thread_envi, &ENVI);		//start the progress bar thread
		ENVI.sift(outfile, mask.data(), true);
		t1.join();									//wait for the progress bar thread to finish (it probably already is)

	}

	else if (args["unsift"].is_set()){

		//get the name of the mask file (an image file) as a string
		std::string maskname = args["unsift"].as_string();

		//load the mask
		stim::image<unsigned char> mask(maskname);

		std::thread t1(progress_thread_envi, &ENVI);		//start the progress bar thread

		//unsift the mask, creating the image outfile
		ENVI.unsift(outfile, mask.data(), mask.width(), mask.height(), true);

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


	}

	else if(args["convert"].is_set()){

		//get the orientation of the destination file
		std::string interleave_name = args["convert"].as_string();
		stim::envi_header::interleaveType interleave_type;
		if(interleave_name == "bip")
			interleave_type = stim::envi_header::BIP;
		else if(interleave_name == "bil")
			interleave_type = stim::envi_header::BIL;
		else if(interleave_name == "bsq")
			interleave_type = stim::envi_header::BSQ;
		else{
			std::cout<<"ERROR - unrecognized interleave format: "<<interleave_name<<std::endl;
			exit(1);
		}


		//convert(infile, outfile, headername, interleave_type);
		std::cout<<"Converting..."<<std::endl;
		std::thread t1(progress_thread_envi, &ENVI);		//start the progress bar thread
		ENVI.convert(outfile, interleave_type, progressbar, verbose, optimization);		//perform the conversion
		ENVI.close();
		t1.detach();									//wait for the progress bar thread to finish (it probably already is)

	}

	//given a text file, output all metrics into a binary file
	else if(args["metrics"]){
		std::vector< std::vector<double> > metrics;					//create a vector to store metrics
		if(args["metrics"].is_num(0)){								//if the first argument of the metrics option is a number, a metric is specified
			metrics.resize(1);										//create one metric
			metrics[0].resize(9, 0);								//the allocate space for the entire metric
			for(size_t i = 0; i < 9; i++) metrics[0][i] = args["metrics"].as_float(i);		//store the metric in the vector
		}
		else{
			std::string metric_file = args["metrics"].as_string(0);
			metrics = read_metric_list(metric_file);
			std::cout<<"Computing "<<metrics.size()<<" metrics..."<<std::endl;
		}
		create_metric_file(metrics, outfile);
	}

	//output an image
	else if(args["image"].is_set()){

		if (args["image"].nargs() == 0){
			std::cout << "ERROR: no image parameters specified" << std::endl;
			exit(1);
		}
		double band = args["image"].as_float(0);
		stim::colormapType cmap = str2cmap(args["colormap"].as_string(0));
		
		if(args["image"].nargs() < 3)
			bandimage(infile, headername, band, outfile, cmap, MASK);
		else{
			double minVal = args["image"].as_float(1);
			double maxVal = args["image"].as_float(2);
			bandimage(infile, headername, band, outfile, minVal, maxVal, cmap, MASK);
		}

	}

	//Output a single spectrum to a CSV file
	else if(args["spectrum"].is_set()){

		if (args["spectrum"].nargs() != 2){
			std::cout << "ERROR: spectrum requires 2 parameters" << std::endl;
			exit(1);
		}
		unsigned int x = args["spectrum"].as_int(0);							//get the x and y positions of the desired spectrum
		unsigned int y = args["spectrum"].as_int(1);

		float* spectrum = (float*)malloc(ENVI.header.bands * sizeof(float));	//allocate space to store the spectrum
		std::thread t1(progress_thread_envi, &ENVI);							//start the progress bar thread
		ENVI.spectrum(spectrum, x, y, true);									//get the spectrum
		t1.join();																//end the progress bar thread

		//output the average spectrum
		std::ofstream csv(outfile.c_str());										//open a CSV file to write the mean

		for (unsigned long long b = 0; b < B; b++) {							//for each band
			if (ENVI.header.wavelength.size() == B)								//output the wavelength, if available
				csv << ENVI.header.wavelength[b] << ',';
			else
				csv << b << ",";
			csv << spectrum[b] << std::endl;									//output the next variable
		}
		csv.close();															//close the file


	}
	//calculate the average spectrum and save to a CSV file
	else if(args["mean"].is_set()){

		unsigned long long XY = X * Y;									//number of pixels in the image		

		double* m = (double*)malloc(sizeof(double) * B);				//allocate space to store the mean spectrum
		double* std = (double*) malloc(sizeof(double) * B);				//allocate space to store the standard deviation

		std::cout<<"Calculating mean spectrum..."<<std::endl;
		std::thread t1(progress_thread_envi, &ENVI);					//start the progress bar thread
		ENVI.mean_spectrum(m, std, MASK, true);							//calculate average spectrum
		t1.join();														//end the progress bar thread

		//output the spectrum
		std::ofstream csv(outfile.c_str());								//open a CSV file to write the mean
		for (unsigned long long b = 1; b < B; b++) {					//for each band
			if (ENVI.header.wavelength.size() == B)						//output the wavelength, if available
				csv << ENVI.header.wavelength[b] << ',';
			else
				csv << b << ",";										//otherwise output a band number
			csv << m[b] << "," << std[b] << std::endl;
		}

		csv.close();													//close the output file

		free(m);														//free memory storing the mean spectrum
		free(std);
	}

	//calculate the median spectrum and save to a CSV file
	else if(args["median"].is_set()){
		double* m = (double*)malloc(B * sizeof(double));				//allocate space for the median spectrum

		std::cout<<"Calculating mean spectrum..."<<std::endl;
		std::thread t1(progress_thread_envi, &ENVI);					//start the progress bar thread
		ENVI.median_spectrum(m, MASK, true);							//calculate average spectrum
		t1.join();														//end the progress bar thread

		//output the median spectrum
		std::ofstream csv(outfile.c_str());								//open a CSV file to write the mean
		csv<<m[0];														//output the first variable
		for(unsigned long long b = 1; b < B; b++)						//for each band
			csv<<","<<m[b];												//output the next variable
		free(m);
	}

	//calculate the covariance matrix
	else if(args["covariance"].is_set()){
		if(ENVI.header.interleave == stim::envi_header::BSQ){
			std::cout<<"ERROR: covariance matrix calculation is not implemented for BSQ files"<<std::endl;
			exit(1);
		}
		//calculate the mean
		unsigned long long XY = X * Y;									//number of pixels in the image

		double* mu = (double*)malloc(sizeof(double) * B);
		double* std = (double*)malloc(sizeof(double) * B);

		std::thread t1;													//create a thread to manage the progress bar
		if(args["covariance"].nargs() == 1){							//if a mean spectrum is given
			std::cout<<"Loading mean spectrum from "<<args["covariance"].as_string(0)<<"..."<<std::endl;
			stim::table mu_table;										//create a CSV table
			mu_table.read_ascii(args["covariance"].as_string(0));		//read the table
			std::vector< std::vector<double> > mu_vector = mu_table.get_vector<double>();	//get the table contents as an STL vector
			for(size_t b = 0; b < B; b++){
				mu[b] = mu_vector[0][b];
			}
		}
		else{
			t1 = std::thread(progress_thread_envi, &ENVI);					//start the progress bar thread
			std::cout<<"Calculating mean spectrum..."<<std::endl;
			ENVI.mean_spectrum(mu, std, MASK, true);									//calculate average spectrum
			t1.join();														//wait for the progress bar thread to finish (it probably already is)
		}

		double *co = (double*)malloc(sizeof(double)* B * B);				//allocate space for the covariance matrix
		
		std::cout<<"Calculating covariance matrix..."<<std::endl;
		t1 = std::thread(progress_thread_envi, &ENVI);							//start the progress bar thread
		int cuda_device = args["cuda"].as_int();
		ENVI.co_matrix(co, mu, MASK, cuda_device, true);						//calculate covariance matrix
		t1.join();														//wait for the progress bar thread to finish (it probably already is)

		//output the covariance matrix
		std::ofstream csv(outfile.c_str());								//open a CSV file to write the mean
		for(unsigned long long b0 = 0; b0 < B; b0++){
			csv<<co[b0 * B];
			for(unsigned long long b1 = 1; b1 < B; b1++)				//for each band
				csv<<","<<co[b0 * B + b1];								//output the next variable
			csv<<std::endl;
		}	
			
		csv.close();													//close the output file	
		free(mu);				//free memory for the mean spectrum
		free(co);				//free memory for the covariance matrix
	}

	else if(args["convolve"].is_set()){
		std::vector<double> coef(args["convolve"].nargs());				//allocate an array of values to store the coefficients
		for(size_t c = 0; c < coef.size(); c++){						//for each coefficient
			std::string s = args["convolve"].as_string(c);
			if(s.find('/') != std::string::npos){
				size_t slash = s.find('/');
				double a = atof(s.substr(0, slash).c_str());
				double b = atof(s.substr(slash+1, s.length() - slash+1).c_str());
				std::cout<<s.substr(slash+1, s.length() - slash+1)<<"-----"<<std::endl;
				coef[c] = a/b;
			}
			else coef[c] = args["convolve"].as_float(c);				//get the floating point value from the parser
			std::cout<<coef[c]<<std::endl;
		}
		std::thread t1(progress_thread_envi, &ENVI);					//start the progress bar thread
		ENVI.convolve(outfile, coef, 0, ENVI.header.bands - coef.size(), 0, MASK, true);
		t1.join();														//wait for the progress bar thread to finish (it probably already is)
	}

	else if(args["sg"].is_set()){
		int N = 5;
		if(args["sg"].nargs() > 0)
			N = args["sg"].as_int(0);
		int D = 3;
		if(args["sg"].nargs() > 1)
			D = args["sg"].as_int(1);

		if(N % 2 == 0 || N <= 0){
			std::cout<<"ERROR: number of points in a Savitzky-Golay filter must be odd and positive"<<std::endl;
			exit(1);
		}
		if(D <= 0){
			std::cout<<"ERROR: the order of a Savitzky-Golay filter must be positive"<<std::endl;
			exit(1);
		}
		if(N < D){
			std::cout<<"ERROR: the number of points for a Savitzky-Golay filter must be greater than the order"<<std::endl;
			exit(1);
		}

		

				
		// calculate the coefficients using the method described in:
		//	https://en.wikipedia.org/wiki/Savitzky%E2%80%93Golay_filter#Derivation_of_convolution_coefficients
		//	C = (J^T * J)^(-1) * Jt
		stim::matrix<double> J(N, D+1);								//create the J matrix
		int r = -N / 2;
		for(int n = 0; n < N; n++){
			for(int d = 0; d < D + 1; d++){
				J(n, d) = pow(n + r, d);
			}
		}
		
		stim::matrix<double> Jt = J.transpose();					//calculate the matrix transpose Jt


		stim::matrix<double> JtJ = Jt * J;							//multiply the transpose by J

		stim::matrix<double> JtJi = JtJ;							//allocate space for the matrix inverse
		int* piv = (int*) malloc(JtJi.rows() * sizeof(int));		//allocate space to store the LU decomposition pivot indices
		//LAPACKE_dgetrf(LAPACK_COL_MAJOR, (int)JtJi.rows(), (int)JtJi.cols(), JtJi.data(), (int)JtJi.rows(), piv);  //use LAPACK for LU decomposition
		//LAPACKE_dgetri(LAPACK_COL_MAJOR, (int)JtJi.rows(), JtJi.data(), (int)JtJi.rows(), piv);					//use LAPACK to solve the inverse
		
		// REPLACED THE ABOVE LAPACKE functions with new CLAPACK functions in linalg.cpp
		LINALG_dgetrf((int)JtJi.rows(), (int)JtJi.cols(), JtJi.data(), (int)JtJi.rows(), piv);  //use LAPACK for LU decomposition
		LINALG_dgetri((int)JtJi.rows(), JtJi.data(), (int)JtJi.rows(), piv);					//use LAPACK to solve the inverse

	
		stim::matrix<double> C = JtJi * Jt;							//calculate C

		std::vector<double> coef(N);								//initialize the coefficient list
		for(int c = 0; c < N; c++)
			coef[c] = C(0, c);

		std::thread t1(progress_thread_envi, &ENVI);					//start the progress bar thread
		ENVI.convolve(outfile, coef, 0, ENVI.header.bands - coef.size(), coef.size()/2, MASK, true);
		t1.join();
				
	}

	else if(args["deriv"].is_set()){
		size_t d = 1;
		if(args["deriv"].nargs() > 0)
			d = args["deriv"].as_int(0);
		size_t order = 2;
		if(args["deriv"].nargs() > 1)
			order = args["deriv"].as_int(1);
		
		std::thread t1(progress_thread_envi, &ENVI);							//start the progress bar thread
		ENVI.deriv(outfile, d, order, MASK, true);
		t1.join();														//wait for the progress bar thread to finish (it probably already is)
	}

	//do PCA and save PCs into a file
	else if (args["pca"].is_set()){

		if(ENVI.header.interleave != stim::envi_header::BIP){
			std::cout<<"Error: Calculating PCA using BSQ or BIL files is impractical. Convert to a BIP format using --convert bip."<<std::endl;
			exit(1);
		}

		unsigned long long XY = ENVI.header.lines * ENVI.header.samples;
		unsigned long long B = ENVI.header.bands;

		//calculate correlation matrix
		double * mu = (double*)malloc(sizeof(double) * B);
		double * std = (double*)malloc(sizeof(double) * B);
		double *co = (double*)malloc(sizeof(double)* B * B);
		std::cout<<"Averaging..."<<std::endl;
		std::thread t1(progress_thread_envi, &ENVI);		//start the progress bar thread
		ENVI.mean_spectrum(mu, std, MASK, true);						//calculate average spectrum

		for(size_t b = 0; b < B; b++){						//for each band, test for non-finite values, exit if found
#ifdef _WIN32
			if(!_finite(mu[b])){							//if the value at index i is not finite
#else
			if(!std::isfinite(mu[b])){					//C++11 implementation
#endif
				std::cout<<"error calculating PCA: the data set contains non-finite values in band "<<b<<". Use --mask-finite to create a mask of finite values."<<std::endl;
				exit(1);
			}
		}
		t1.join();											//wait for the progress bar thread to finish (it probably already is)


		std::cout<<"Calculating the covariance matrix..."<<std::endl;
		t1 =std::thread(progress_thread_envi, &ENVI);		//start the progress bar thread
		int cuda_device = args["cuda"].as_int();
		ENVI.co_matrix(co, mu, MASK, cuda_device, true);		//calculate correlation coefficient
		t1.join();									//wait for the progress bar thread to finish (it probably already is)


		//calculate eigen values and eigen vectors
		//cv::Mat correlation((int)B, (int)B, CV_64FC1, co);
		//cv::Mat eigenvalues((int)B, 1, CV_64FC1);
		//cv::Mat eigenvectors((int)B, (int)B, CV_64FC1);
		std::cout<<"Calculating eigenvectors (LAPACK)..."<<std::endl;
		//cv::eigen(correlation, eigenvalues, eigenvectors);		//calculate eigen values and eigen vectors

		double* lambda_real = (double*) malloc(B * sizeof(double));
		double* lambda_imag = (double*) malloc(B * sizeof(double));
		double* evec = (double*) malloc(B * B * sizeof(double));
		//LAPACKE_dgeev(LAPACK_COL_MAJOR, 'N', 'V', (int)B, co, (int)B, lambda_real, lambda_imag, NULL, (int)B, evec, (int)B);

		// REPLACED THE ABOVE LAPACKE functions with new CLAPACK functions in linalg.cpp
		LINALG_dgeev('N', 'V', (int)B, co, (int)B, lambda_real, lambda_imag, NULL, (int)B, evec, (int)B);


		std::ofstream csv(outfile.c_str(), std::ios::out);		//create a text file to store the PCA stats (mean and covariance matrix)

		csv<<mu[0];												//output the mean spectrum
		for(unsigned long long i = 1; i < B; i++)
			csv<<","<<mu[i];
		csv<<std::endl;

		for(size_t j = 0; j < B; j++){
			csv<<evec[j * B + 0];							//output the first element of the eigenvector
			for(size_t i = 1; i < B; i++){
				csv<<","<<evec[j * B + i];					//output all consecutive elements
			}
			csv<<std::endl;
		}
		csv.close();		//close the text file

		free(mu);
		free(co);
	}

	else if (args["project"].is_set()){
		std::string pcfile = args["project"].as_string(0);

		//recommend user to convert file to BIP file before do rotation
		if (ENVI.header.interleave == stim::envi_header::BSQ || ENVI.header.interleave == stim::envi_header::BIL){
			std::cout << "ERROR: Rotation is only practical for a BIP image; convert the image to BIP" << std::endl;
			exit(1);
		}
		unsigned long long XY = ENVI.header.lines * ENVI.header.samples;					//calculate the size of the image
		unsigned long long B = ENVI.header.bands;										//calculate the number of bands
		unsigned long long N = B;													//set the default number of PCs to the maximum
		if(args["project"].nargs() > 1)											//if the number of PCs is specified
			N = args["project"].as_int(1);										//retrieve the specified number of PCs

		std::ifstream csv(pcfile.c_str(), std::ios::in);							//open the statistics file
		if(!csv){ std::cout<<"ERROR reading statistics file: "<<pcfile<<std::endl; exit(1); }		//make sure the stats file is valid
		double* basis = (double*)malloc(sizeof(double) * B * N);					//allocate space for the basis matrix

		double *mu = (double*)malloc(sizeof(double) * B);								//allocate space for the mean feature vector

		std::cout<<"Loading basis statistics..."<<std::endl;
		
		std::string line, token;
		std::getline(csv, line);
		std::stringstream ss(line);
		for(unsigned i = 0; i < B; i++){											//load the mean feature vector
			std::getline(ss, token, ',');
			mu[i] = atof(token.c_str());						
		}

		for (unsigned long long n = 0; n < N; n++){										//for each component
			std::getline(csv, line);
			//std::cout<<N<<"     "<<B<<"     "<<prog<<std::endl;
			
			std::stringstream ss1(line);
 			//ss = std::stringstream(line);
			for (unsigned long long b = 0; b < B; b++){									//for each feature
				std::getline(ss1, token, ',');
				basis[n * B + b] = atof(token.c_str());
				//prog = (double)((n+1) * B + b+1) / (double)(B * N) * 100;
				//std::cout<<b<<", "<<n<<", "<<B<<", "<<N<<", "<<prog<<std::endl;
			}
			//exit(1);
		}
		//t1.join();

		std::cout<<"Projecting onto the new basis..."<<std::endl;
		//double prog = 0;
		std::thread t1(progress_thread_envi, &ENVI);		//start the progress bar thread
		std::vector<double> pcs(N);
		if (N == B) pcs = ENVI.header.wavelength;
		else {
			for (size_t n = 0; n < N; n++)
				pcs[n] = (double)(n + 1);
		}
		int cuda_device = args["cuda"].as_int();
		ENVI.project(outfile, mu, basis, N, pcs, MASK, cuda_device, true);
		t1.join();
	}

	else if (args["inverse"].is_set()){
		std::string pcfile = args["inverse"].as_string(0);

		//recommend user to convert file to BIP file before do rotation
		if (ENVI.header.interleave == stim::envi_header::BSQ || ENVI.header.interleave == stim::envi_header::BIL){
			std::cout << "ERROR: Rotation is only practical for a BIP image; convert the image to BIP" << std::endl;
			exit(1);
		}
		unsigned long long XY = ENVI.header.lines * ENVI.header.samples;	//calculate the size of the image
		unsigned long long M = ENVI.header.bands;							//calculate the number of bands
		unsigned long long C = M;											//set the default number of coefficients to the maximum
		if(args["inverse"].nargs() > 1)										//if the number of coefficients to be used is specified
			C = args["inverse"].as_int(1);									//retrieve the specified number of coefficients

		std::ifstream csv(pcfile.c_str(), std::ios::in);					//open the statistics file
		std::string line, token;
		std::getline(csv, line);
		unsigned long long B = 1;											//count the number of spectral components that will be produced
		for(unsigned long long c = 0; c < line.size(); c++){				//for each character in the first line
			if(line[c] == ',') B++;												//count the number of commas
		}

		double* basis = (double*)malloc(sizeof(double) * B * C);			//allocate space for the basis matrix
		double *mu = (double*)malloc(sizeof(double) * B);					//allocate space for the mean feature vector

		std::stringstream ss(line);
		for(unsigned long long i = 0; i < B; i++){									//load the mean feature vector
			std::getline(ss, token, ',');
			mu[i] = atof(token.c_str());						
		}

		for (unsigned long long n = 0; n < C; n++){								//for each component
			std::getline(csv, line);
			//ss = std::stringstream(line);
			
			std::stringstream ss1(line);
			for (unsigned long long b = 0; b < B; b++){							//for each feature
				std::getline(ss1, token, ',');
				basis[n * B + b] = atof(token.c_str());
			}
		}

		std::thread t1(progress_thread_envi, &ENVI);		//start the progress bar thread
		ENVI.inverse(outfile, mu, basis, B, C, true);
		t1.join();
	}

	else if (args["crop"].is_set()){
		//initialize the values to the maximum extents
		unsigned long long x = 0;
		unsigned long long y = 0;
		unsigned long long b = 0;
		unsigned long long nx = X;
		unsigned long long ny = Y;
		unsigned long long nb = B;

		// Define the default behavior for the parameters
		//		(I'm trying to be intuitive)

		if(args["crop"].as_string(0)[0] != '-'){
			x = args["crop"].as_int(0);
			nx = X - x;
		}
		if(args["crop"].nargs() > 1 && args["crop"].as_string(1)[0] != '-')
			nx = args["crop"].as_int(1);

		if(args["crop"].nargs() > 2 && args["crop"].as_string(2)[0] != '-'){
			y = args["crop"].as_int(2);
			ny = Y - y;
		}
		if(args["crop"].nargs() > 3 && args["crop"].as_string(3)[0] != '-')
			ny = args["crop"].as_int(3);

		if(args["crop"].nargs() > 4 && args["crop"].as_string(4)[0] != '-'){
			if (ENVI.header.wavelength.size()) {									//if wavelength values are specified in the file
				double wavelength = args["crop"].as_float(4);							//assume that the coordinates are given in wavelength
				size_t low, high;
				ENVI.band_bounds(wavelength, low, high);
				b = low;
			}
			else {																	//if wavelength values aren't given, assume the user parameter is in bands
				b = args["crop"].as_int(4);
			}
			nb = B - b;
		}
		if (args["crop"].nargs() > 5 && args["crop"].as_string(5)[0] != '-') {
			if (ENVI.header.wavelength.size()) {									//if wavelength values are specified in the file
				double nw = args["crop"].as_float(5);							//assume that the coordinates are given in wavelength
				double wavelength = args["crop"].as_float(4) + nw;
				size_t low, high;
				ENVI.band_bounds(wavelength, low, high);
				nb = high - b;
			}
			else {																	//if wavelength values aren't given, assume the user parameter is in bands
				b = args["crop"].as_int(5);
			}
		}

		if(x + nx > X || y + ny > Y || b + nb > B){
			std::cout<<"ERROR: specified image is out of bounds."<<std::endl;
			exit(1);
		}
		std::cout<<"Cropping image to [ "<<nx<<" x "<<ny<<" x "<<nb<<" ]..."<<std::endl;
		std::thread t1(progress_thread_envi, &ENVI);		//start the progress bar thread
		ENVI.crop(outfile, x, y, x + nx - 1, y + ny - 1, b, b + nb - 1, true);	//call the ENVI cropping function

		t1.join();									//wait for the progress bar thread to finish (it probably already is)
	}
	//generate a set of subimages from a mask file
	else if(args["subimages"].is_set()){
		std::thread t1(progress_thread_envi, &ENVI);			//start the progress bar thread

		set_mask(args["subimages"].as_string(0));
		
		size_t nx = 11;
		if(args["subimages"].nargs() >= 2)
			nx = args["subimages"].as_int(1);				//get the width of the subimages
		size_t ny = nx;											//assume that the subimages are square
		if(args["subimages"].nargs() >= 3)						//if the image height is specified
			ny = args["subimages"].as_int(2);					//set the height

		ENVI.subimages(outfile, nx, ny, MASK, true);
		t1.join();									//wait for the progress bar thread to finish (it probably already is)
	}
	else if(args["trim"].is_set()){
		std::vector<size_t> bands;								//create a list of band indices to trim
		std::vector<size_t> idx;
		//std::vector< std::pair<double, double> > trim_ranges;				//create an array of pairs
		double w0, w1;														//create a pair of wavelength values
		size_t i;
		for(size_t p = 0; p < args["trim"].nargs(); p++){
			std::string s = args["trim"].as_string(p);						//get the first value as a string
			i = s.find('-');												//get the index of the dash (if it exists)
			if(i != std::string::npos && i != 0){							//if the string contains a dash (-) and that dash isn't the first value (negative number)
				s[i] = ' ';													//replace the dash with a space
				std::stringstream ss(s);									//create a string stream
				ss>>w0;														//store the first and second values as a range
				ss>>w1;
				idx = ENVI.header.band_indices(w0, w1);						//get the band indices in this range
				bands.insert(bands.end(), idx.begin(), idx.end());			//insert them into the band array
			}
			else{													//otherwise
				w0 = args["trim"].as_float(p);					//otherwise store the single wavelength
				idx = ENVI.header.band_index(w0);				//get the band(s) associated with that wavelength
				if(idx.size() == 1)								//if only one band exists
					bands.push_back(idx[0]);					//add it to the band array (if there are two indices, there isn't a single band associated with the wavelength)
			}
		}
		std::cout<<"Trimming "<<bands.size()<<" / "<<ENVI.header.bands<<" bands"<<std::endl;
		std::thread t1(progress_thread_envi, &ENVI);		//start the progress bar thread
		ENVI.trim(outfile, bands, true);					//trim the bands and store the result
		t1.join();									//wait for the progress bar thread to finish (it probably already is)
	}
	else if(args["combine"].is_set()){
		std::string otherfile(args["combine"].as_string(0));			//get the name for the file to combine with this one
		stim::envi C;
		C.open(otherfile, otherfile + ".hdr");							//open the file to combine

		long long px = 0;												//get the position of the second image
		long long py = 0;
		if(args["combine"].nargs() == 1){								//if no arguments are given, the files are concatenated along Y
			px = 0;
			py = ENVI.header.lines;
		}
		if(args["combine"].nargs() > 1){								//if one argument is provided, assume it is X
			px = args["combine"].as_int(1);
			if(px < (long long)ENVI.header.samples) py = ENVI.header.samples;		//if no Y is given, and x is less than the maximum X, concatenate along Y
		}
		if(args["combine"].nargs() > 2)
			py = args["combine"].as_int(2);
		if( (-px < (long long)C.header.samples)		&&
			( px < (long long)ENVI.header.samples)	&&
			(-py < (long long)C.header.lines)		&&
			( py < (long long)ENVI.header.lines)){
			std::cout<<"ERROR: The files to be combined are overlapping. The extent for the input image is: ["<<ENVI.header.samples<<", "<<ENVI.header.lines<<"]";
			exit(1);
		}
		std::cout<<"Combining images..."<<std::endl;
		std::thread t1(progress_thread_envi, &ENVI);		//start the progress bar thread
		ENVI.combine(outfile, C, px, py, true);
		t1.join();											//wait for the progress bar thread to finish (it probably already is)
	}
	else if (args["append"]) {								//append a file to the input in the band dimension
		std::string otherfile(args["append"].as_string(0));	//get the name of the file to append to this one
		stim::envi C;										//create an ENVI file object
		C.open(otherfile, otherfile + ".hdr");				//open the appending file

		std::thread t1(progress_thread_envi, &ENVI);		//start the progress bar thread
		ENVI.append(outfile, C, true);						//append the files and write the output
		t1.join();											//wait for the progress bar thread to finish (it probably already is)
		C.close();											//close the appended file
	}
	
	else if(args["mnf"].is_set()){
		if (args["mnf"].nargs() == 0){
			std::cout << "ERROR: no mnf parameters specified" << std::endl;
			exit(1);
		}
		int keptComponents = args["mnf"].as_int(0);						//get the user-specified number of components to keep
		std::string component_file;
		if(args["mnf"].nargs() == 2)
			component_file = args["mnf"].as_string(1);
		int device = args["cuda"].as_int();
		mnf(outfile, keptComponents, component_file, device);
	}

	else if(args["multiply"].is_set()){
		std::thread t1(progress_thread_envi, &ENVI);					//start the progress bar thread
		double val = args["multiply"].as_float(0);						//get the value to multiply by
		ENVI.multiply(outfile, val, MASK, true);						//multiply by val
		t1.join();														//wait for the progress bar thread to finish
	}
	else if(args["add"].is_set()){
		std::thread t1(progress_thread_envi, &ENVI);					//start the progress bar thread
		double val = args["add"].as_float(0);						//get the value to multiply by
		ENVI.add(outfile, val, MASK, true);						//multiply by val
		t1.join();														//wait for the progress bar thread to finish
	}
	else{
		std::cout<<"ERROR: no algorithm is specified"<<std::endl;
		exit(1);
	}

	time_t t2=time(NULL);			//end time
	size_t processing_time = t2 - t1;
	std::cout<<"    processing time : "<<processing_time<<" s"<<std::endl<<std::endl<<std::endl;
	std::cout<<"    data rate:        "<<(double)ENVI.bytes() / 1000000 / processing_time<<" MB/s"<<std::endl;

	//close the input file
	ENVI.close();
}