envi/bip.h

#ifndef STIM_BIP_H
#define STIM_BIP_H
#include "../envi/envi_header.h"
#include "../envi/bil.h"
#include "../envi/binary.h"
#include <cstring>
#include <utility>
namespace stim{
template <typename T>
class bip: public binary<T> {
protected:
	
	
	std::vector<double> w; //band wavelength
	unsigned int offset;		//header offset
public:
	using binary<T>::open;
	using binary<T>::file;
	using binary<T>::R;
	//open a file, given the file and its header's names
	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 offset to the structure
		offset = header_offset;
		return open(filename, vec<unsigned int>(X, Y, B), header_offset);
		
	}
	//save one band of the file into the memory, and return the pointer
	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;
		}
		//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 < R[0] * R[1]; i++)
		{
			file.read((char *)(p + i), sizeof(T));
			file.seekg( (R[2] - 1) * sizeof(T), std::ios::cur);
		}
		return true;
	}
	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]) {
				band_index(p, 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;
	}
	//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 = R[0];	//calculate the number of pixels in a sample
		unsigned int B = R[2];
		unsigned page=0;                      //samples 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 ){
			for(unsigned j = 0; j < R[0]; j++)
			{
				p[j] = c[j * B];
			}		
			return true;
		}
		while( w[page] < wavelength )
		{
			page++;
			//if wavelength is larger than the last wavelength in header file
			if (page == B) {
				for(unsigned j = 0; j < R[0]; j++)
				{
					p[j] = c[(j + 1) * B - 1];
				}
				return true;
			}
		}
		if ( wavelength < w[page] )
		{
			T * p1;
			T * p2;
			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 - 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);
			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 = R[0];	//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 ( w[page] > wavelength ){
			file.seekg( y * R[2] * R[0] * sizeof(T), std::ios::beg);
			for(unsigned j = 0; j < R[0]; j++)
			{
				file.read((char *)(p + j), sizeof(T));
				file.seekg((R[2] - 1) * sizeof(T), std::ios::cur);
			}		
			return true;
		}
		while( w[page] < wavelength )
		{
			page++;
			//if wavelength is larger than the last wavelength in header file
			if (page == R[2]) {
				file.seekg( (y * R[2] * R[0] + R[2] - 1)* sizeof(T), std::ios::beg);
				for(unsigned j = 0; j < R[0]; j++)
				{
					file.read((char *)(p + j), sizeof(T));
					file.seekg((R[2] - 1) * sizeof(T), std::ios::cur);
				}
				return true;
			}
		}
		if ( wavelength < w[page] )
		{
			p1=(T*)malloc( X * sizeof(T));                     //memory allocation
			p2=(T*)malloc( X * sizeof(T));
			//band_index(p1, page - 1);
			file.seekg( (y * R[2] * R[0] + page - 1)* sizeof(T), std::ios::beg);
			for(unsigned j = 0; j < R[0]; j++)
			{
				file.read((char *)(p1 + j), sizeof(T));
				file.seekg((R[2] - 1) * sizeof(T), std::ios::cur);
			}
			//band_index(p2, page );
			file.seekg( (y * R[2] * R[0] + page)* sizeof(T), std::ios::beg);
			for(unsigned j = 0; j < R[0]; j++)
			{
				file.read((char *)(p2 + j), sizeof(T));
				file.seekg((R[2] - 1) * sizeof(T), std::ios::cur);
			}
			
			for(unsigned i=0; i < R[0]; i++){
				double r = (double) (wavelength - w[page-1]) / (double) (w[page] - w[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 * R[2] * R[0] + page) * sizeof(T), std::ios::beg);
			for(unsigned j = 0; j < R[0]; j++)
			{
				file.read((char *)(p + j), sizeof(T));
				file.seekg((R[2] - 1) * sizeof(T), std::ios::cur);
			}
		}
		return true;
	}
	//save one pixel of the BIP file into the memory, and return the pointer
	bool spectrum(T * p, unsigned x, unsigned y){
		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;
			exit(1);
		}
		file.seekg((x + y * R[0]) * R[2] * sizeof(T), std::ios::beg);           //point to the certain sample and line
		
		file.read((char *)p, sizeof(T) * R[2]);
			
		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 * R[2] * sizeof(T), std::ios::beg);           //point to the certain pixel
		file.read((char *)p, sizeof(T) * R[2]);
		return true;	
	}
	
	//given a Y ,return a ZX slice
	bool getY(T * p, unsigned y)
	{
		if ( y >= R[1]){							//make sure the line number is right
			std::cout<<"ERROR: line out of range"<<std::endl;
			exit(1);
		}	
		
		file.seekg(y * R[2] * R[0] * sizeof(T), std::ios::beg);
		file.read((char *)p, sizeof(T) * R[2] * R[0]);
		
		return true;
		
	}
		
	//(BIP) baseline correction
	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 = R[2] * R[0];		//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 = R[2];
		unsigned int X = R[0];
		
		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;
		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 < R[1]; k++)
		{
			//get the current y slice
			getY(c, k);
		
			//initialize lownum, highnum, low, high		
			control=0;
			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, 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( 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];
						getYZ(b, c, 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, X * sizeof(T) );	//clear the high band
	
						ai = bi;
						bi = w[B - 1];
					}
				}
				ci = w[cii];
			
				//perform the baseline correction
				for(unsigned i=0; i < X; i++)
				{
					double r = (double) (ci - ai) / (double) (bi - ai);
					c[i * B + cii] =(T) ( 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
		
		
		free(a);
		free(b);
		free(c);
		target.close();
		return true;	
		
	}
		
	// normalize the BIP file
	bool normalize(std::string outname, double w, double t = 0.0)
	{
		unsigned int B = R[2];		//calculate the number of bands
		unsigned int Y = R[1];
		unsigned int X = R[0];
		unsigned int ZX = R[2] * R[0];
		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		
		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 ); 
		band(b, w);             //get the certain band into memory
		for(unsigned j = 0; j < Y; j++)
		{
			getY(c, j);
			unsigned jX = j * X;		//to avoid calculating it many times
			for(unsigned i = 0; i < X; i++)
			{
				unsigned iB = i * B;		
				for(unsigned m = 0; m < B; m++)
				{
					if( b[i+jX] < t )
						c[m + iB] = (T)0.0;
					else
						c[m + iB] = c[m + iB] / b[i + jX];			//perform normalization
				}								
			}
			target.write(reinterpret_cast<const char*>(c), L);   //write normalized data into destination
		}
		
		free(b);
		free(c);
		target.close();
		return true;
	}
	
	//convert BIP file to BSQ file and save it
	bool bsq(std::string outname)
	{	
		std::string temp = outname + "_temp";
		std::string headtemp = temp + ".hdr";
		//first creat a temporary bil file and convert bip 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 bsq file
		n.bsq(outname);
		n.close();
		remove(temp.c_str());
		remove(headtemp.c_str());
		return true;
	}
	//convert bip file to bil file and save it
	bool bil(std::string outname)
	{
		unsigned int S = R[0] * R[2] * sizeof(T);		//calculate the number of bytes in a ZX slice
		
		std::ofstream target(outname.c_str(), std::ios::binary);
		std::string headername = outname + ".hdr";
		
		T * p;			//pointer to the current ZX slice for bip file
		p = (T*)malloc(S);
		T * q;			//pointer to the current XZ slice for bil file
		q = (T*)malloc(S);
		
		for ( unsigned i = 0; i < R[1]; i++)
		{			
			getY(p, i);
			for ( unsigned k = 0; k < R[2]; k++)
			{
				unsigned ks = k * R[0];
				for ( unsigned j = 0; j < R[0]; j++)
					q[ks + j] = p[k + j * R[2]];
			}
			target.write(reinterpret_cast<const char*>(q), S);   //write a band data into target file	
		}
		
		free(p);
		free(q);
		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;
	}
	
	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){
		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 ZX = R[2] * R[0];		//calculate number of a band
		unsigned L = ZX * sizeof(T);
		T * temp = (T*)malloc(L);
		for (unsigned i = 0; i < R[1]; i++)
		{
			getY(temp, i);
			for ( unsigned j = 0; j < R[0]; j++)
			{
				for (unsigned k = 0; k < R[2]; k++)
				{
				if(p[i * R[0] + k] == 0)
					temp[j * R[2] + k] = 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 of the file
	bool band_avg(T* p){
		unsigned long long XY = R[0] * R[1];
		//get every pixel and calculate average value
		T* temp = (T*)malloc(sizeof(T) * R[2]);
		T sum;
		for (unsigned i = 0; i < XY; i++){
			pixel(temp, i);
			//calculate the sum value of every value
			sum = 0;		//initialize sum value
			for (unsigned j = 0; j < R[2]; j++){
				sum += temp[j]/(T)R[2];
			}
			p[i] = sum;
		}
		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];
		T* temp = (T*)malloc(sizeof(T) * R[2]);
		//Iinitialize
		for (unsigned j = 0; j < R[2]; j++){
			p[j] = 0;
		}
		//calculate vaild number in a band
		unsigned count = 0;
		for (unsigned j = 0; j < XY; j++){
			if (mask[j] != 0){
				count++;
			}
		}
		//calculate average number of a band
		for (unsigned i = 0; i < XY; i++){
			if (mask[i] != 0){			
				pixel(temp, i);
				for (unsigned j = 0; j < R[2]; j++){
					p[j] += temp[j] / (T)count;
				}
			}
		}
		free(temp);
		return true;
	}
	//calculate correlation coefficient 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* temp = (T*)malloc(sizeof(T) * B);
		//count vaild pixels in a band
		unsigned count = 0;
		for (unsigned j = 0; j < xy; j++){
			if (mask[j] != 0){
				count++;
			}
		}
		//initialize correlation matrix
		for (unsigned i = 0; i < B; i++){
			for (unsigned k = 0; k < B; k++){
				co[i * B + k] = 0;
			}
		}
		//calculate correlation coefficient matrix
		for (unsigned j = 0; j < xy; j++){
			if (mask[j] != 0){
				pixel(temp, j);
				for (unsigned i = 0; i < B; i++){
					for (unsigned k = i; k < B; k++){
						co[i * B + k] += (temp[i] - avg[i]) * (temp[k] - avg[k]) / count;
					}
				}
			}
		}
		//because correlation matrix is symmetric
		for (unsigned i = 0; i < B; i++){
			for (unsigned k = i + 1; k < B; k++){
				co[k * B + i] = co[i * B + k];
			}
		}
		free(temp);
		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 = R[2] * sizeof(T);
		//get specified band and save
		T* temp = (T*)malloc(L);
		std::ofstream out(outfile.c_str(), std::ios::binary);
		//get start
		unsigned long long sp = y0 * R[0] + x0;		//start pixel
		for (unsigned i = 0; i < lin; i++)
		{
			for (unsigned j = 0; j < sam; j++)
			{
				pixel(temp, sp + j + i * R[0]);
				out.write(reinterpret_cast<const char*>(temp), L);   //write slice data into target file	
			}
		}
		free(temp);
		return true;
	}
	//close the file
	bool close(){
		file.close();
		return true;
	}
	};
}
#endif