envi/bip.h

#include "../envi/envi.h"
#include "../envi/binary.h"
#include <cstring>
#include <utility>
namespace rts{
template <typename T>
class bip: public binary<T> {
protected:
	
	envi header;
public:
	using binary<T>::open;
	using binary<T>::file;
	using binary<T>::getSlice;
	//open a file, given the file and its header's names
	bool open(std::string filename, std::string headername){
		if (header.load(headername)==false){
			std::cout<<"ERROR: unable to load header file: "<<headername<<std::endl;
			return false;
		}
		open(filename, vec<unsigned int>(header.samples, header.lines, header.bands), header.header_offset);
		return true;
		
	}
	//save one band of the file into the memory, and return the pointer
	bool band_index( T * p, unsigned int page){
		if (page >= header.bands){										//make sure the bank number is right
			std::cout<<"ERROR: page out of range"<<std::endl;
			return false;
		}
		//binary::getSlice(p, 0, page);         //I met some problems when I called this function?????
		file.seekg(page * sizeof(T), std::ios::beg);
		for (unsigned i = 0; i < header.samples * header.lines; i++)
		{
			file.read((char *)(p + i), sizeof(T));
			file.seekg( (header.bands - 1) * sizeof(T), std::ios::cur);
		}
		return true;
	}
	bool getBand( T * p, double wavelength){
		unsigned int XY = header.samples * header.lines;	//calculate the number of pixels in a band
		unsigned page=0;                      //bands around the wavelength
		T * p1;
		T * p2;
		//get the bands numbers around the wavelength
		//if wavelength is smaller than the first one in header file
		if ( header.wavelength[page] > wavelength ){
			band_index(p, page);
			return true;
		}
		while( header.wavelength[page] < wavelength )
		{
			page++;
			//if wavelength is larger than the last wavelength in header file
			if (page == header.bands) {
				band_index(p, header.bands-1);
				return true;
			}
		}
		if ( wavelength < header.wavelength[page] )
		{
			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 - header.wavelength[page-1]) / (double) (header.wavelength[page] - header.wavelength[page-1]);
				p[i] = (p2[i] - p1[i]) * r + p1[i];
			}
		}
		else                           //if the wavelength is equal to a wavelength in header file
		{
			band_index(p, page);
		}
		free(p1);
		free(p2);
		return true;
	}
	//get YZ line from the a Y slice, Y slice data should be already IN the MEMORY
	bool getYZ(T* p, T* c, double wavelength)
	{
		unsigned int X = header.samples;	//calculate the number of pixels in a sample
		unsigned int B = header.bands;
		unsigned page=0;                      //samples around the wavelength
		T * p1;
		T * p2;
		//get the bands numbers around the wavelength
		//if wavelength is smaller than the first one in header file
		if ( header.wavelength[page] > wavelength ){
			for(unsigned j = 0; j < header.samples; j++)
			{
				p[j] = c[j * B];
			}		
			return true;
		}
		while( header.wavelength[page] < wavelength )
		{
			page++;
			//if wavelength is larger than the last wavelength in header file
			if (page == B) {
				for(unsigned j = 0; j < header.samples; j++)
				{
					p[j] = c[(j + 1) * B - 1];
				}
				return true;
			}
		}
		if ( wavelength < header.wavelength[page] )
		{
			p1=(T*)malloc( X * sizeof(T));                     //memory allocation
			p2=(T*)malloc( X * sizeof(T));
			//band_index(p1, page - 1);
			for(unsigned j = 0; j < X; j++)
			{
				p1[j] = c[j * B + page - 1];
			}
			//band_index(p2, page );
			for(unsigned j = 0; j < X; j++)
			{
				p2[j] = c[j * B + page];
			}
			
			for(unsigned i=0; i < X; i++){
				double r = (double) (wavelength - header.wavelength[page-1]) / (double) (header.wavelength[page] - header.wavelength[page-1]);
				p[i] = (p2[i] - p1[i]) * r + p1[i];
			}
		}
		else                           //if the wavelength is equal to a wavelength in header file
		{
			//band_index(p, page);
			for(unsigned j = 0; j < X; j++)
			{
				p[j] = c[j * B + page];
			}
		}
		return true;		
	}
	//given a y and a wavelength, return the y-band data
	//I do not use it right now, to accelerate the processing speed, I try to read a YZ whole slice into memory first, and then we can call "getYZ"
	bool get_x_wavelength(T * p, unsigned y, double wavelength)
	{
		unsigned int X = header.samples;	//calculate the number of pixels in a sample
		unsigned page=0;                      //samples around the wavelength
		T * p1;
		T * p2;
		//get the bands numbers around the wavelength
		//if wavelength is smaller than the first one in header file
		if ( header.wavelength[page] > wavelength ){
			file.seekg( y * header.bands * header.samples * sizeof(T), std::ios::beg);
			for(unsigned j = 0; j < header.samples; j++)
			{
				file.read((char *)(p + j), sizeof(T));
				file.seekg((header.bands - 1) * sizeof(T), std::ios::cur);
			}		
			return true;
		}
		while( header.wavelength[page] < wavelength )
		{
			page++;
			//if wavelength is larger than the last wavelength in header file
			if (page == header.bands) {
				file.seekg( (y * header.bands * header.samples + header.bands - 1)* sizeof(T), std::ios::beg);
				for(unsigned j = 0; j < header.samples; j++)
				{
					file.read((char *)(p + j), sizeof(T));
					file.seekg((header.bands - 1) * sizeof(T), std::ios::cur);
				}
				return true;
			}
		}
		if ( wavelength < header.wavelength[page] )
		{
			p1=(T*)malloc( X * sizeof(T));                     //memory allocation
			p2=(T*)malloc( X * sizeof(T));
			//band_index(p1, page - 1);
			file.seekg( (y * header.bands * header.samples + page - 1)* sizeof(T), std::ios::beg);
			for(unsigned j = 0; j < header.samples; j++)
			{
				file.read((char *)(p1 + j), sizeof(T));
				file.seekg((header.bands - 1) * sizeof(T), std::ios::cur);
			}
			//band_index(p2, page );
			file.seekg( (y * header.bands * header.samples + page)* sizeof(T), std::ios::beg);
			for(unsigned j = 0; j < header.samples; j++)
			{
				file.read((char *)(p2 + j), sizeof(T));
				file.seekg((header.bands - 1) * sizeof(T), std::ios::cur);
			}
			
			for(unsigned i=0; i < header.samples; i++){
				double r = (double) (wavelength - header.wavelength[page-1]) / (double) (header.wavelength[page] - header.wavelength[page-1]);
				p[i] = (p2[i] - p1[i]) * r + p1[i];
			}
		}
		else                           //if the wavelength is equal to a wavelength in header file
		{
			//band_index(p, page);
			file.seekg( (y * header.bands * header.samples + page) * sizeof(T), std::ios::beg);
			for(unsigned j = 0; j < header.samples; j++)
			{
				file.read((char *)(p + j), sizeof(T));
				file.seekg((header.bands - 1) * sizeof(T), std::ios::cur);
			}
		}
		//free(p1);
		//free(p2);
		return true;
	}
	//save one pixel of the BIP file into the memory, and return the pointer
	//I have not used it so far
	bool getSpectrum(T * p, unsigned x, unsigned y){
		unsigned int i;
		if ( x >= header.samples || y >= header.lines){							//make sure the sample and line number is right
			std::cout<<"ERROR: sample or line out of range"<<std::endl;
			exit(1);
		}
		file.seekg((x + y * header.samples) * header.bands * sizeof(T), std::ios::beg);           //point to the certain sample and line
		
		file.read((char *)p, sizeof(T) * header.bands);
			
		return true;	
	}
	
	//given a Y ,return a ZX slice
	bool getY(T * p, unsigned y)
	{
		if ( y >= header.lines){							//make sure the line number is right
			std::cout<<"ERROR: line out of range"<<std::endl;
			exit(1);
		}	
		
		file.seekg(y * header.bands * header.samples * sizeof(T), std::ios::beg);
		file.read((char *)p, sizeof(T) * header.bands * header.samples);
		
		return true;
		
	}
		
	//(BIP) baseline correction
	//processing speed is too slow
	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 ZX = header.bands * header.samples;		//calculate the number of points in a Y slice
		unsigned int L = ZX * sizeof(T);			//calculate the number of bytes of a Y slice
		unsigned int B = header.bands;
		unsigned int X = header.samples;
		
		T* c;			//pointer to the current Y slice
		c = (T*)malloc(L);  //memory allocation
		
		T* a;			//pointer to the two YZ lines surrounding the current YZ line
		T* b;
		
		a = (T*)malloc(X * sizeof(T));
		b = (T*)malloc(X * sizeof(T));
		double ai, bi;	//stores the two baseline points wavelength surrounding the current band
		double ci;		//stores the current band's wavelength
		unsigned control=0;
		if (a == NULL || b == NULL || c == NULL){
			std::cout<<"ERROR: error allocating memory";
			exit(1);
		}
	//	loop start	correct every y slice
		
		for (unsigned k =0; k < header.lines; k++)
		{
			//get the current y slice
			getY(c, k);
		
			//initialize lownum, highnum, low, high		
			ai = header.wavelength[0];
			//if no baseline point is specified at band 0,
			//set the baseline point at band 0 to 0
			if(wls[0] != header.wavelength[0]){
				bi = wls[control];			
				memset(a, (char)0, X * sizeof(T) );
			}
			//else get the low band
			else{
				control++;
				getYZ(a, c, ai);
				bi = wls[control];
			}
			//get the high band
			getYZ(b, c, bi);
		
			//correct every YZ line
			
			for(unsigned cii = 0; cii < B; cii++){
				//update baseline points, if necessary
				if( header.wavelength[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];
						getYZ(b, c, bi);
	
					}
					//if the last BL point on the last band of the file?
					else if ( wls[control] < header.wavelength[B - 1]) {
	
						std::swap(a, b);	//swap the baseline band pointers
	
						memset(b, (char)0, X * sizeof(T) );	//clear the high band
	
						ai = bi;
						bi = header.wavelength[B - 1];
					}
				}
				//get the current YZ line
				//band_index(c, cii);
				file.seekg( (k * header.bands * X + cii)* sizeof(T), std::ios::beg);
				for(unsigned j = 0; j < X; j++)
				{
					file.read((char *)(c + j), sizeof(T));
					file.seekg((B - 1) * sizeof(T), std::ios::cur);
				}
			
				ci = header.wavelength[cii];
			
				//perform the baseline correction
				for(unsigned i=0; i < X; i++)
				{
					double r = (double) (ci - ai) / (double) (bi - ai);
					c[i * B + cii] =(float) ( c[i * B + cii] - (b[i] - a[i]) * r - a[i] );
				}
				
			}//loop for YZ line end  
			
			target.write(reinterpret_cast<const char*>(c), L);   //write the corrected data into destination	
		}//loop for Y slice end
		
		header.save(headername);         //save the new header file
		
		free(a);
		free(b);
		free(c);
		target.close();
		return true;	
		
	}
		
	// normalize the BIP file
	bool normalize(std::string outname, double band)
	{
		unsigned int B = header.bands;		//calculate the number of bands
		unsigned int Y = header.lines;
		unsigned int X = header.samples;
		unsigned int ZX = header.bands * header.samples;
		unsigned int XY = header.samples * header.lines;	//calculate the number of pixels in a band
		unsigned int S = XY * sizeof(T);		//calculate the number of bytes in a band		
		unsigned int L = ZX * sizeof(T);
		std::ofstream target(outname.c_str(), std::ios::binary);	//open the target binary file
		std::string headername = outname + ".hdr";              //the header file name
		T * c;				//pointer to the current ZX slice
		T * b;				//pointer to the standard band
		b = (T*)malloc( S );     //memory allocation
		c = (T*)malloc( L ); 
		getBand(b, band);             //get the certain band into memory
		for(unsigned j = 0; j < Y; j++)
		{
			getY(c, j);
			for(unsigned i = 0; i < X; i++)
			{
				for(unsigned m = 0; m < B; m++)
				{
					c[m + i * B] = c[m + i * B] / b[i + j * X];
				}								
			}
			target.write(reinterpret_cast<const char*>(c), L);   //write normalized data into destination
		}
		header.save(headername);         //save the new header file
		
		free(b);
		free(c);
		target.close();
		return true;
	}
	
	//convert BIP file to BSQ file and save it
	bool bsq(std::string outname)
	{
		unsigned int S = header.samples * header.lines * sizeof(T);		//calculate the number of bytes in a band
		
		std::ofstream target(outname.c_str(), std::ios::binary);
		std::string headername = outname + ".hdr";
		
		T * p;			//pointer to the current band
		p = (T*)malloc(S);
		
		for ( unsigned i = 0; i < header.bands; i++)
		{			
				band_index(p, i);
				target.write(reinterpret_cast<const char*>(p), S);   //write a band data into target file	
		}
		header.interleave = rts::envi::interleaveType::BSQ;  //change the type of file in header file
		header.save(headername);
		
		free(p);
		target.close();
		return true;
	}
	};
}