codebase / stimlib

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

namespace rts{

template <typename T>

class bsq: public binary<T> {

protected:
	
	envi header;

	Assistant Professor
	CPRIT Scholar
	Scalable Tissue Imaging + Modeling Laboratory

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 head file"<<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;
		}

		getSlice(p, 2, page);
		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) {
				getSlice(p, 2, 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
		{
			getSlice(p, 2, page);
		}

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

	//save one pixel of the file into the memory, and return the pointer
	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;
			return false;
		}

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

		return true;	
	}
	/*
	// for bands only
	bool baselineold(std::string outname, std::vector<unsigned> bands )
	{
		unsigned N = bands.size();			//get the number of baseline points

		std::ofstream target(outname.c_str(), std::ios::binary);	//open the target binary file

		//simplify image resolution
		unsigned int B = header.bands;		//calculate the number of bands
		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 ai, bi;	//stores the two baseline points 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=0;

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

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

			//update baseline points, if necessary
			if( ci == bi && ci != 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 = bands[control];
					band_index(b, bi);

				}
				//if the last BL point on the last band of the file?
				else if ( bands[control] != B - 1) 
				{

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

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

					ai = bi;
					bi = B - 1;
				}
			}

			//get the current band
			band_index(c, ci);

			//perform the baseline correction
			for(unsigned i=0; i < XY; i++){
				double r = (double) (ci - ai) / (double) (bi - ai);
				c[i] = c[i] - (b[i] - a[i]) * r - a[i];
			}
			
			target.write(reinterpret_cast<const char*>(c), S);   //write the corrected data into destination
		
		}
		
		free(a);
		free(b);
		free(c);
		target.close();
		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

		//simplify image resolution
		unsigned int B = header.bands;		//calculate the number of bands
		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

		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=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, S);
		}
		//else get the low band
		else{
			control += 1;
			getBand(a, ai);
			bi = wls[control];
		}
		//get the high band
		getBand(b, bi);

		//correct every band 
		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];
					getBand(b, 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, S);	//clear the high band

					ai = bi;
					bi = header.wavelength[B - 1];
				}
			}

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

	};
}