codebase / stimlib

#ifndef STIM_BSQ_H
#define STIM_BSQ_H

#include "../envi/envi_header.h"
#include "../envi/binary.h"
#include "../envi/bil.h"
#include <cstring>
#include <utility>
#include <vector>

namespace stim{

template <typename T>

class bsq: public binary<T> {

protected:
	
	//envi_header header;

	std::vector<double> w;	//band wavelengths
	unsigned int offset;

public:

	using binary<T>::open;
	using binary<T>::file;
	using binary<T>::getSlice;
	using binary<T>::R;

	//open a file, given the file name and dimensions
	bool open(std::string filename, unsigned int X, unsigned int Y, unsigned int B, unsigned int header_offset, std::vector<double> wavelengths){

		//copy the wavelengths to the BSQ file structure
		w = wavelengths;
		//copy the wavelengths to the structure
		offset = header_offset;

		return open(filename, vec<unsigned int>(X, Y, B), header_offset);
		
	}

	//retrieve one band (specified by the band index)
	bool band_index( T * p, unsigned int page){

		if (page >= R[2]){										//make sure the bank number is right
			std::cout<<"ERROR: page out of range"<<std::endl;
			return false;
		}

		file.seekg(R[0] * R[1] * page * sizeof(T));

		file.read((char *)p, sizeof(T) * R[0] * R[1]);

		return true;
	}

	//retrieve one band (specified by the wavelength)
	bool band( T * p, double wavelength){

		//if there are no wavelengths in the BSQ file
		if(w.size() == 0)
			return band_index(p, (unsigned int)wavelength);

		unsigned int XY = R[0] * R[1];	//calculate the number of pixels in a band
		unsigned page = 0;                


		//get the two neighboring bands (above and below 'wavelength')

		//if wavelength is smaller than the first one in header file
		if ( w[page] > wavelength ){
			band_index(p, page);
			return true;
		}

		while( w[page] < wavelength )
		{
			page++;
			//if wavelength is larger than the last wavelength in the header file
			//	(the wavelength is out of bounds)
			if (page == R[2]) {
				band_index(p, R[2]-1);		//return the last band
				return true;
			}
		}
		//when the page counter points to the first band above 'wavelength'
		if ( wavelength < w[page] ){

			//do the interpolation
			T * p1;
			T * p2;
			p1=(T*)malloc( XY * sizeof(T));                     //memory allocation
			p2=(T*)malloc( XY * sizeof(T));
			band_index(p1, page - 1);
			band_index(p2, page );
			for(unsigned i=0; i < XY; i++){
				double r = (double) (wavelength - w[page-1]) / (double) (w[page] - w[page-1]);
				p[i] = (p2[i] - p1[i]) * r + p1[i];
			}
			free(p1);
			free(p2);
		}
		//if the wavelength is equal to a wavelength in header file
		else{
			band_index(p, page);		//return the band
		}

		return true;
	}

	//save one pixel of the file into the memory, and return the pointer
	bool spectrum(T * p, unsigned x, unsigned y){

		unsigned int i;

		if ( x >= R[0] || y >= R[1]){							//make sure the sample and line number is right
			std::cout<<"ERROR: sample or line out of range"<<std::endl;
			return false;
		}

		file.seekg((x + y * R[0]) * sizeof(T), std::ios::beg);           //point to the certain sample and line
		for (i = 0; i < R[2]; i++)
		{
			file.read((char *)(p + i), sizeof(T));
			file.seekg((R[1] * R[0] - 1) * sizeof(T), std::ios::cur);    //go to the next band
		}

		return true;	
	}

	//save one pixel into memory
	bool pixel(T * p, unsigned n){

		unsigned bandnum = R[0] * R[1];		//calculate numbers in one band
		if ( n >= bandnum){							//make sure the pixel number is right
			std::cout<<"ERROR: sample or line out of range"<<std::endl;
			return false;
		}

		file.seekg(n * sizeof(T), std::ios::beg);           //point to the certain pixel
		for (unsigned i = 0; i < R[2]; i++)
		{
			file.read((char *)(p + i), sizeof(T));
			file.seekg((bandnum - 1) * sizeof(T), std::ios::cur);    //go to the next band
		}

		return true;	
	}

	//baseline correction and save it into file
	
	bool baseline(std::string outname, std::vector<double> wls )
	{
		unsigned N = wls.size();			//get the number of baseline points

		std::ofstream target(outname.c_str(), std::ios::binary);	//open the target binary file
		std::string headername = outname + ".hdr";              //the header file name

		//simplify image resolution
		unsigned int B = R[2];		//calculate the number of bands
		unsigned int XY = R[0] * R[1];	//calculate the number of pixels in a band
		unsigned int S = XY * sizeof(T);		//calculate the number of bytes in a band

		double ai, bi;	//stores the two baseline points wavelength surrounding the current band
		double ci;		//stores the current band's wavelength
//	unsigned aii, bii;	//stores the two baseline points number surrounding the current band
		unsigned control=0;

		T * a;					//pointers to the high and low band images
		T * b;
		T * c;				//pointer to the current image

		a = (T*)malloc( S );     //memory allocation
		b = (T*)malloc( S ); 
		c = (T*)malloc( S ); 

		if (a == NULL || b == NULL || c == NULL){
			std::cout<<"ERROR: error allocating memory";
			exit(1);
		}


		//initialize lownum, highnum, low, high		
		ai=w[0];

		//if no baseline point is specified at band 0,
			//set the baseline point at band 0 to 0
		if(wls[0] != w[0]){
			bi = wls[control];			
			memset(a, (char)0, S);
		}
		//else get the low band
		else{
			control += 1;
			band(a, ai);
			bi = wls[control];
		}
		//get the high band
		band(b, bi);

		//correct every band 
		for(unsigned cii = 0; cii < B; cii++){

			//update baseline points, if necessary
			if( w[cii] >= bi && cii != B - 1) {
				//if the high band is now on the last BL point?
				if (control != N-1) {

					control++;		//increment the index

					std::swap(a, b);	//swap the baseline band pointers

					ai = bi;
					bi = wls[control];
					band(b, bi);

				}
				//if the last BL point on the last band of the file?
				else if ( wls[control] < w[B - 1]) {

					std::swap(a, b);	//swap the baseline band pointers

					memset(b, (char)0, S);	//clear the high band

					ai = bi;
					bi = w[B - 1];
				}
			}

			//get the current band
			band_index(c, cii);
			ci = w[cii];
			
			//perform the baseline correction
			for(unsigned i=0; i < XY; i++){
				double r = (double) (ci - ai) / (double) (bi - ai);
				c[i] =(T) ( c[i] - (b[i] - a[i]) * r - a[i] );
			}
			
			target.write(reinterpret_cast<const char*>(c), S);   //write the corrected data into destination
		
		}	

		//header.save(headername);         //save the new header file
		
		free(a);
		free(b);
		free(c);
		target.close();
		return true;
	}

	// normalize the BSQ file
	bool normalize(std::string outname, double w)
	{
		unsigned int B = R[2];		//calculate the number of bands
		unsigned int XY = R[0] * R[1];	//calculate the number of pixels in a band
		unsigned int S = XY * sizeof(T);		//calculate the number of bytes in a band

		std::ofstream target(outname.c_str(), std::ios::binary);	//open the target binary file
		std::string headername = outname + ".hdr";              //the header file name

		T * b;					//pointers to the certain wavelength band
		T * c;				//pointer to the current image

		b = (T*)malloc( S );     //memory allocation
		c = (T*)malloc( S ); 

		band(b, w);             //get the certain band into memory

		for(unsigned j = 0; j < B; j++)
		{
			band_index(c, j);                     //get the current band into memory
			for(unsigned i = 0; i < XY; i++)
			{
				c[i] = c[i] / b[i];
			}
			target.write(reinterpret_cast<const char*>(c), S);   //write normalized data into destination
		}

		//header.save(headername);         //save the new header file
		
		free(b);
		free(c);
		target.close();
		return true;
	}
	
	//convert BSQ file to BIP file and save it
	bool bip(std::string outname)
	{
		std::string temp = outname + "_temp";
		std::string headtemp = temp + ".hdr";
		//first creat a temporary bil file and convert bsq file to bil file
		bil(temp);

		stim::bil<T> n;
		if(n.open(temp, R[0], R[1], R[2], offset, w)==false){        //open infile
			std::cout<<"ERROR: unable to open input file"<<std::endl;
			exit(1);
		}
		//then convert bil file to bip file
		n.bip(outname);
		n.close();
		remove(temp.c_str());
		remove(headtemp.c_str());
		return true;
	}

	//convert BSQ file to BIL file and save it
	bool bil(std::string outname)
	{
		//simplify image resolution
		unsigned int L = R[0] * R[2] * sizeof(T);		//calculate the number of bytes of a ZX slice
		unsigned int jump = (R[1] - 1) * R[0] * sizeof(T);
		
		std::ofstream target(outname.c_str(), std::ios::binary);
		std::string headername = outname + ".hdr";
		
		T * p;			//pointer to the current spectrum
		p = (T*)malloc(L);
		
		for ( unsigned i = 0; i < R[1]; i++)
		{
			file.seekg(R[0] * i * sizeof(T), std::ios::beg);
			for ( unsigned j = 0; j < R[2]; j++ )
			{
				file.read((char *)(p + j * R[0]), sizeof(T) * R[0]);
				file.seekg(jump, std::ios::cur);    //go to the next band	
			}					
			target.write(reinterpret_cast<const char*>(p), L);   //write XZ slice data into target file	
		}
		//header.interleave = rts::envi_header::BIL;  //change the type of file in header file
		//header.save(headername);
		
		free(p);
		target.close();
		return true;
	}

	//providing the left and the right bound data, return baseline-corrected band height
	bool baseline_band(double lb, double rb, T* lp, T* rp, double wavelength, T* result){

		unsigned XY = R[0] * R[1];
		band(result, wavelength);		//get band

		//perform the baseline correction
		double r = (double) (wavelength - lb) / (double) (rb - lb);
		for(unsigned i=0; i < XY; i++){
			result[i] =(T) (result[i] - (rp[i] - lp[i]) * r - lp[i] );
		}
		return true;
	}

	//providing the left and the right bound wavelength, return baseline-corrected band height
	bool height(double lb, double rb, double bandwavelength, T* result){

		T* lp;
		T* rp;
		unsigned XY = R[0] * R[1];
		unsigned S = XY * sizeof(T);
		lp = (T*) malloc(S);			//memory allocation
		rp = (T*) malloc(S);

		band(lp, lb);
		band(rp, rb);		

		baseline_band(lb, rb, lp, rp, bandwavelength, result);

		free(lp);
		free(rp);
		return true;
	}


	//calculate the area between two bound point(including baseline correction)
	bool area(double lb, double rb, double lab, double rab, T* result){

		T* lp;	//left band pointer
		T* rp;	//right band pointer
		T* cur;		//current band 1
		T* cur2;	//current band 2

		unsigned XY = R[0] * R[1];
		unsigned S = XY * sizeof(T);

		lp = (T*) malloc(S);			//memory allocation
		rp = (T*) malloc(S);
		cur = (T*) malloc(S);
		cur2 = (T*) malloc(S);

		memset(result, (char)0, S);

		//find the wavelenght position in the whole band
		unsigned int n = w.size();
		unsigned int ai = 0;		//left bound position
		unsigned int bi = n - 1;		//right bound position



		//to make sure the left and the right bound are in the bandwidth
		if (lb < w[0] || rb < w[0] || lb > w[n-1] || rb >w[n-1]){
			std::cout<<"ERROR: left bound or right bound out of bandwidth"<<std::endl;
			exit(1);
		}
		//to make sure rigth bound is bigger than left bound
		else if(lb > rb){
			std::cout<<"ERROR: right bound should be bigger than left bound"<<std::endl;
			exit(1);
		}

		//get the position of lb and rb
		while (lab >= w[ai]){
			ai++;
		}
		while (rab <= w[bi]){
			bi--;
		}

		band(lp, lb);
		band(rp, rb);

		//calculate the beginning and the ending part
		baseline_band(lb, rb, lp, rp, rab, cur2);		//ending part
		baseline_band(lb, rb, lp, rp, w[bi], cur);
		for(unsigned j = 0; j < XY; j++){
			result[j] += (rab - w[bi]) * (cur[j] + cur2[j]) / 2.0;
		}
		baseline_band(lb, rb, lp, rp, lab, cur2);		//beginnning part
		baseline_band(lb, rb, lp, rp, w[ai], cur);
		for(unsigned j = 0; j < XY; j++){	
			result[j] += (w[ai] - lab) * (cur[j] + cur2[j]) / 2.0;
		}

		//calculate the area
		ai++;
		for(unsigned i = ai; i <= bi ;i++)
		{
			baseline_band(lb, rb, lp, rp, w[ai], cur2);
			for(unsigned j = 0; j < XY; j++)
			{
				result[j] += (w[ai] - w[ai-1]) * (cur[j] + cur2[j]) / 2.0; 
			}
			std::swap(cur,cur2);		//swap the band pointers
		}

		free(lp);
		free(rp);
		free(cur);
		free(cur2);
		return true;
	}

	//peak height ratio
	bool ph_to_ph(double lb1, double rb1, double pos1, double lb2, double rb2, double pos2, T * result){

		T* p1 = (T*)malloc(R[0] * R[1] * sizeof(T));	
		T* p2 = (T*)malloc(R[0] * R[1] * sizeof(T));	
		
		//get the two peak band
		height(lb1, rb1, pos1, p1);
		height(lb2, rb2, pos2, p2);
		//calculate the ratio in result
		for(unsigned i = 0; i < R[0] * R[1]; i++){
			if(p1[i] == 0 && p2[i] ==0)
				result[i] = 1;
			else
				result[i] = p1[i] / p2[i];
		}

		free(p1);
		free(p2);
		return true;	
	}
	
	//peak are to peak height ratio
	bool pa_to_ph(double lb1, double rb1, double lab1, double rab1,
					double lb2, double rb2, double pos, T* result){
		
		T* p1 = (T*)malloc(R[0] * R[1] * sizeof(T));	
		T* p2 = (T*)malloc(R[0] * R[1] * sizeof(T));	
		
		//get the area and the peak band
		area(lb1, rb1, lab1, rab1, p1);
		height(lb2, rb2, pos, p2);
		//calculate the ratio in result
		for(unsigned i = 0; i < R[0] * R[1]; i++){
			if(p1[i] == 0 && p2[i] ==0)
				result[i] = 1;
			else
				result[i] = p1[i] / p2[i];
		}

		free(p1);
		free(p2);
		return true;	
	}		
	
	//peak area to peak area ratio
	bool pa_to_pa(double lb1, double rb1, double lab1, double rab1,
					double lb2, double rb2, double lab2, double rab2, T* result){
		
		T* p1 = (T*)malloc(R[0] * R[1] * sizeof(T));	
		T* p2 = (T*)malloc(R[0] * R[1] * sizeof(T));	
		
		//get the area and the peak band
		area(lb1, rb1, lab1, rab1, p1);
		area(lb2, rb2, lab2, rab2, p2);
		//calculate the ratio in result
		for(unsigned i = 0; i < R[0] * R[1]; i++){
			if(p1[i] == 0 && p2[i] ==0)
				result[i] = 1;
			else
				result[i] = p1[i] / p2[i];
		}

		free(p1);
		free(p2);
		return true;	
	}		

	//x * f(x)
	bool x_area(double lb, double rb, double lab, double rab, T* result){
		T* lp;	//left band pointer
		T* rp;	//right band pointer
		T* cur;		//current band 1
		T* cur2;	//current band 2

		unsigned XY = R[0] * R[1];
		unsigned S = XY * sizeof(T);

		lp = (T*) malloc(S);			//memory allocation
		rp = (T*) malloc(S);
		cur = (T*) malloc(S);
		cur2 = (T*) malloc(S);

		memset(result, (char)0, S);

		//find the wavelenght position in the whole band
		unsigned int n = w.size();
		unsigned int ai = 0;		//left bound position
		unsigned int bi = n - 1;		//right bound position

		//to make sure the left and the right bound are in the bandwidth
		if (lb < w[0] || rb < w[0] || lb > w[n-1] || rb >w[n-1]){
			std::cout<<"ERROR: left bound or right bound out of bandwidth"<<std::endl;
			exit(1);
		}
		//to make sure rigth bound is bigger than left bound
		else if(lb > rb){
			std::cout<<"ERROR: right bound should be bigger than left bound"<<std::endl;
			exit(1);
		}

		//get the position of lb and rb
		while (lab >= w[ai]){
			ai++;
		}
		while (rab <= w[bi]){
			bi--;
		}

		band(lp, lb);
		band(rp, rb);

		//calculate the beginning and the ending part
		baseline_band(lb, rb, lp, rp, rab, cur2);		//ending part
		baseline_band(lb, rb, lp, rp, w[bi], cur);
		for(unsigned j = 0; j < XY; j++){
			result[j] += (rab - w[bi]) * (rab + w[bi]) * (cur[j] + cur2[j]) / 4.0;
		}
		baseline_band(lb, rb, lp, rp, lab, cur2);		//beginnning part
		baseline_band(lb, rb, lp, rp, w[ai], cur);
		for(unsigned j = 0; j < XY; j++){	
			result[j] += (w[ai] - lab) * (w[ai] + lab) * (cur[j] + cur2[j]) / 4.0;
		}

		//calculate f(x) times x
		ai++;
		for(unsigned i = ai; i <= bi ;i++)
		{
			baseline_band(lb, rb, lp, rp, w[ai], cur2);
			for(unsigned j = 0; j < XY; j++)
			{
				result[j] += (w[ai] - w[ai-1]) * (w[ai] + w[ai-1]) * (cur[j] + cur2[j]) / 4.0; 
			}
			std::swap(cur,cur2);		//swap the band pointers
		}

		free(lp);
		free(rp);
		free(cur);
		free(cur2);
		return true;
	}

	//centroid point
	bool cpoint(double lb, double rb, double lab, double rab, T* result){
		T* p1 = (T*)malloc(R[0] * R[1] * sizeof(T));	
		T* p2 = (T*)malloc(R[0] * R[1] * sizeof(T));	
		
		//get the area and the peak band
		x_area(lb, rb, lab, rab, p1);
		area(lb, rb, lab, rab, p2);
		//calculate the ratio in result
		for(unsigned i = 0; i < R[0] * R[1]; i++){
			if(p1[i] == 0 && p2[i] ==0)
				result[i] = 1;
			else
				result[i] = p1[i] / p2[i];
		}

		free(p1);
		free(p2);
		return true;	
	}

	//create mask file
	bool build_mask(double mask_band, double threshold, unsigned char* p = NULL){

		T* temp = (T*)malloc(R[0] * R[1] * sizeof(T));		//allocate memory for the certain band
		band(temp, mask_band);

		for (unsigned int i = 0; i < R[0] * R[1]; i++) {
				if (temp[i] < threshold)
					p[i] = 0;
				else
					p[i] = 255;
		}

		free(temp);
		return true;

	}

	//apply mask
	bool apply_mask(std::string outfile, unsigned char* p){

		std::ofstream target(outfile.c_str(), std::ios::binary);

		unsigned XY = R[0] * R[1];		//calculate number of a band
		unsigned L = XY * sizeof(T);

		T * temp = (T*)malloc(L);

		for (unsigned i = 0; i < R[2]; i++)
		{
			band_index(temp, i);
			for ( unsigned j = 0; j < XY; j++)
			{
				if(p[j] == 0){
					temp[j] = 0;
				}
				else{
					continue;
				}
			}
			target.write(reinterpret_cast<const char*>(temp), L);   //write a band data into target file	
		}
		target.close();
		free(temp);
		return true;
	}

	//calculate the average band
	bool band_avg(T* p){
		unsigned long long XY = R[0] * R[1];
		T* temp = (T*)malloc(sizeof(T) * XY);
		//initialize p
		band_index(p, 0);
		for (unsigned j = 0; j < XY; j++){
			p[j] /= (T)R[2];
		}
		//get every band and add them all
		for (unsigned i = 1; i < R[2]; i++){
			band_index(temp, i);
			for (unsigned j = 0; j < XY; j++){
				p[j] += temp[j]/(T)R[2];
			}
		}
		free(temp);
		return true;
	}

	//calculate the average number of every band
	bool avg_band(T*p, unsigned char* mask){
		unsigned long long XY = R[0] * R[1];
		unsigned count = 0;		//number of vaild pixel in a band
		T* temp = (T*)malloc(sizeof(T) * XY);
		//calculate valid pixel number
		for (unsigned j = 0; j < XY; j++){
			if (mask[j] != 0){
				count++;
			}
		}
		//calculate average of a band
		for (unsigned i = 0; i < R[2]; i++){
			p[i] = 0;
			band_index(temp, i);
			for (unsigned j = 0; j < XY; j++){
				if (mask[j] != 0){
					p[i] += temp[j] / (T)count;
				}				
			}
		}
		free(temp);
		return true;
	}

	//calculate correlated matrix
	bool co_matrix(T* co, T* avg, unsigned char *mask){
		//memory allocation
		unsigned long long xy = R[0] * R[1];
		unsigned int B = R[2];
		T* bandi = (T*)malloc(sizeof(T) * xy);
		T* bandj = (T*)malloc(sizeof(T) * xy);
		
		//count vaild pixels in a band
		unsigned count = 0;
		for (unsigned j = 0; j < xy; j++){
			if (mask[j] != 0){
				count++;
			}
		}
		//calculate correlation coefficient matrix
		for (unsigned i = 0; i < B; i++)
		{
			band_index(bandi, i);
			for (unsigned j = i; j < B; j++){
				band_index(bandj, j);
				T numerator = 0;			//to calculate element in correlation coefficient matrix, numerator part
				//calculate the R(i,j) in correlation coeffient matrix
				for (unsigned k = 0; k < xy; k++){
					if (mask[k] != 0){
						numerator += (bandi[k] - avg[i]) * (bandj[k] - avg[j]) / count;
					}
				}
				co[i*B + j] = numerator;
				co[j*B + i] = numerator;		//because correlated matrix is a symmetric matrix
			}
		}
		free(bandi);
		free(bandj);
		return true;
	}


	//crop specified area the of the original file
	bool crop(std::string outfile, unsigned x0, unsigned y0, unsigned x1, unsigned y1){

		//calculate the new number of samples and lines
		unsigned long long sam = x1 - x0;		//samples
		unsigned long long lin = y1 - y0;		//lines
		unsigned long long L = sam * lin * sizeof(T);
		//get specified band and save
		T* temp = (T*)malloc(L);
		std::ofstream out(outfile.c_str(), std::ios::binary);
		unsigned long long jumpb = R[0] * (R[1] - lin) * sizeof(T);		//jump pointer to the next band
		unsigned long long jumpl = (R[0] - sam) * sizeof(T);	                           	//jump pointer to the next line
		//get start
		file.seekg((y0 * R[0] + x0) * sizeof(T), std::ios::beg);
		for (unsigned i = 0; i < R[2]; i++)
		{
			for (unsigned j = 0; j < lin; j++)
			{
				file.read((char *)(temp + j * sam), sizeof(T) * sam);
				file.seekg(jumpl, std::ios::cur);    //go to the next band	
			}
			out.write(reinterpret_cast<const char*>(temp), L);   //write slice data into target file	
			file.seekg(jumpb, std::ios::cur);
		}
		free(temp);
		return true;
	}


	//close the file
	bool close(){
		file.close();
		return true;
	}

	};
}

#endif