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), std::ios::beg);
  
  		file.read((char *)p, sizeof(T) * R[0] * R[1]);
  
  //		getSlice(p, 2, page);

  		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;                      //bands around the wavelength

  
  		//get the bands numbers around the 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 header file

  			if (page == R[2]) {
  				getSlice(p, 2, R[2]-1);

  				return true;
  			}
  		}

  		if ( wavelength < w[page] )

  		{

  			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);

  		}
  		else                           //if the wavelength is equal to a wavelength in header file
  		{

  			band_index(p, page);

  		}
  

  		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[3]; 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;	
  	}

  
  	//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;
  	}
  
  

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

  	};
  }

  #endif