simplified the ENVI classes, moved more stuff into binary.h

David Mayerich
1 parent 5f81932b
Showing 5 changed files with 495 additions and 493 deletions Show diff stats
stim/envi/bil.h
stim/envi/binary.h
stim/envi/bip.h
stim/envi/bsq.h
stim/envi/envi.h
@@ -25,6 +25,16 @@ protected:
 	std::vector<double> w;		//band wavelength
+	unsigned long X(){
+		return R[0];
+	}
+	unsigned long Y(){
+		return R[2];
+	}
+	unsigned long Z(){
+		return R[1];
+	}
+
 public:
 	using binary<T>::open;
@@ -43,7 +53,7 @@ public:
 		w = wavelengths;
-		return open(filename, vec<unsigned int>(X, Y, B), header_offset);
+		return open(filename, vec<unsigned int>(X, B, Y), header_offset);
 	}
@@ -52,19 +62,20 @@ public:
 	/// @param p is a pointer to an allocated region of memory at least X * Y * sizeof(T) in size.
 	/// @param page <= B is the integer number of the band to be copied.
 	bool band_index( T * p, unsigned int page){
+		//return binary<T>::read_plane_1(p, page);
-		unsigned int L = R[0] * sizeof(T);		//caculate the number of bytes in a sample line
-		unsigned int jump = R[0] * (R[2] - 1) * sizeof(T);
+		unsigned int L = X() * sizeof(T);		//caculate the number of bytes in a sample line
+		unsigned int jump = X() * (Z() - 1) * sizeof(T);
-		if (page >= R[2]){										//make sure the bank number is right
+		if (page >= Z()){										//make sure the bank number is right
 			std::cout<<"ERROR: page out of range"<<std::endl;
 			return false;
 		}
-		file.seekg(R[0] * page * sizeof(T), std::ios::beg);
-		for (unsigned i = 0; i < R[1]; i++)
+		file.seekg(X() * page * sizeof(T), std::ios::beg);
+		for (unsigned i = 0; i < Y(); i++)
 		{
-			file.read((char *)(p + i * R[0]), L);
+			file.read((char *)(p + i * X()), L);
 			file.seekg( jump, std::ios::cur);
 		}
@@ -81,7 +92,7 @@ public:
 		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 int XY = X() * Y();	//calculate the number of pixels in a band
 		unsigned int S = XY * sizeof(T);		//calculate the number of bytes of a band
 		unsigned page=0;                      //bands around the wavelength
@@ -99,8 +110,8 @@ public:
 		{
 			page++;
 			//if wavelength is larger than the last wavelength in header file
-			if (page == R[2]) {
-				band_index(p, R[2]-1);
+			if (page == Z()) {
+				band_index(p, Z()-1);
 				return true;
 			}
 		}
@@ -127,12 +138,15 @@ public:
 		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)
+	/// Retrieves a band of x values from a given xz plane.
+
+	/// @param p is a pointer to pre-allocated memory at least Z * sizeof(T) in size
+	/// @param c is a pointer to an existing XZ plane (size X*Z*sizeof(T))
+	/// @param wavelength is the wavelength to retrieve
+	bool read_x_from_xz(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 int L = X * sizeof(T);
+		unsigned int B = Z();
+		unsigned int L = X() * sizeof(T);
 		unsigned page=0;                      //samples around the wavelength
 		T * p1;
@@ -151,7 +165,7 @@ public:
 			page++;
 			//if wavelength is larger than the last wavelength in header file
 			if (page == B) {
-				memcpy(p, c + (B - 1) * X, L);
+				memcpy(p, c + (B - 1) * X(), L);
 				return true;
 			}
 		}
@@ -160,16 +174,16 @@ public:
 			p1=(T*)malloc( L );                     //memory allocation
 			p2=(T*)malloc( L );
-			memcpy(p1, c + (page - 1) * X, L);
-			memcpy(p2, c + page * X, L);
+			memcpy(p1, c + (page - 1) * X(), L);
+			memcpy(p2, c + page * X(), L);
-			for(unsigned i=0; i < X; i++){
+			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];
 			}
 		}
 		else                           //if the wavelength is equal to a wavelength in header file		
-			memcpy(p, c + page * X, L);
+			memcpy(p, c + page * X(), L);
 		return true;		
 	}
@@ -180,24 +194,7 @@ public:
 	/// @param x is the x-coordinate (dimension 1) of the spectrum.
 	/// @param y is the y-coordinate (dimension 2) of the spectrum.
 	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);
-		}
-
-		unsigned jump = (R[0] - 1) * sizeof(T);
-
-		file.seekg((y * R[0] * R[2] + x) * sizeof(T), std::ios::beg);
-
-		for(unsigned i = 0; i < R[2]; i++)
-		{       
-			//point to the certain sample and line		
-			file.read((char *)(p + i), sizeof(T));
-			file.seekg(jump, std::ios::cur);  
-		}	
-
-		return true;	
+		return binary<T>::read_line_02(p, x, y);
 	}
 	/// Retrieve a single pixel and stores it in pre-allocated memory.
@@ -207,26 +204,16 @@ public:
 	bool pixel(T * p, unsigned n){
 		//calculate the corresponding x, y
-		unsigned int x = n % R[0];
-		unsigned int y = n / R[0];
+		unsigned int x = n % X();
+		unsigned int y = n / X();
 		//get the pixel
 		return spectrum(p, x, y);
 	}
 	//given a Y ,return a XZ 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;
-		
+	bool read_plane_y(T * p, unsigned y){
+		return binary<T>::read_plane_2(p, y);
 	}
@@ -242,19 +229,18 @@ public:
 		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 ZX = Z() * X();		//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];
-		
+		unsigned int B = Z();
+
 		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));
+		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
@@ -267,10 +253,10 @@ public:
 		}
 	//	loop start	correct every y slice
-		for (unsigned k =0; k < R[1]; k++)
+		for (unsigned k =0; k < Y(); k++)
 		{
 			//get the current y slice
-			getY(c, k);
+			read_plane_y(c, k);
 			//initialize lownum, highnum, low, high		
 			ai = w[0];
@@ -279,16 +265,16 @@ public:
 			//set the baseline point at band 0 to 0
 			if(wls[0] != w[0]){
 				bi = wls[control];			
-				memset(a, (char)0, X * sizeof(T) );
+				memset(a, (char)0, X() * sizeof(T) );
 			}
 			//else get the low band
 			else{
 				control++;
-				getYZ(a, c, ai);
+				read_x_from_xz(a, c, ai);
 				bi = wls[control];
 			}
 			//get the high band
-			getYZ(b, c, bi);
+			read_x_from_xz(b, c, bi);
 			//correct every YZ line
@@ -305,7 +291,7 @@ public:
 						ai = bi;
 						bi = wls[control];
-						getYZ(b, c, bi);
+						read_x_from_xz(b, c, bi);
 					}
 					//if the last BL point on the last band of the file?
@@ -313,7 +299,7 @@ public:
 						std::swap(a, b);	//swap the baseline band pointers
-						memset(b, (char)0, X * sizeof(T) );	//clear the high band
+						memset(b, (char)0, X() * sizeof(T) );	//clear the high band
 						ai = bi;
 						bi = w[B - 1];
@@ -322,9 +308,9 @@ public:
 				ci = w[cii];
-				unsigned jump = cii * X;
+				unsigned jump = cii * X();
 				//perform the baseline correction
-				for(unsigned i=0; i < X; i++)
+				for(unsigned i=0; i < X(); i++)
 				{
 					double r = (double) (ci - ai) / (double) (bi - ai);
 					c[i + jump] =(T) ( c[i + jump] - (b[i] - a[i]) * r - a[i] );
@@ -349,11 +335,9 @@ public:
 	///	@param t is a threshold specified such that a spectrum with a value at w less than t is set to zero. Setting this threshold allows the user to limit division by extremely small numbers.
 	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 B = Z();		//calculate the number of bands
+		unsigned int ZX = Z() * X();
+		unsigned int XY = X() * Y();	//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);
@@ -368,17 +352,17 @@ public:
 		band(b, w);             //get the certain band into memory
-		for(unsigned j = 0; j < Y; j++)
+		for(unsigned j = 0; j < Y(); j++)
 		{
-			getY(c, j);
+			read_plane_y(c, j);
 			for(unsigned i = 0; i < B; i++)
 			{
-				for(unsigned m = 0; m < X; m++)
+				for(unsigned m = 0; m < X(); m++)
 				{
-					if( b[m + j * X] < t )
-						c[m + i * X] = (T)0.0;
+					if( b[m + j * X()] < t )
+						c[m + i * X()] = (T)0.0;
 					else
-						c[m + i * X] = c[m + i * X] / b[m + j * X];
+						c[m + i * X()] = c[m + i * X()] / b[m + j * X()];
 				}								
 			}
 			target.write(reinterpret_cast<const char*>(c), L);   //write normalized data into destination
@@ -395,7 +379,7 @@ public:
 	/// @param outname is the name of the output BSQ file to be saved to disk.
 	bool bsq(std::string outname)
 	{
-		unsigned int S = R[0] * R[1] * sizeof(T);		//calculate the number of bytes in a band
+		unsigned int S = X() * Y() * sizeof(T);		//calculate the number of bytes in a band
 		std::ofstream target(outname.c_str(), std::ios::binary);
 		std::string headername = outname + ".hdr";
@@ -403,7 +387,7 @@ public:
 		T * p;			//pointer to the current band
 		p = (T*)malloc(S);
-		for ( unsigned i = 0; i < R[2]; i++)
+		for ( unsigned i = 0; i < Z(); i++)
 		{			
 				band_index(p, i);
 				target.write(reinterpret_cast<const char*>(p), S);   //write a band data into target file	
@@ -419,7 +403,7 @@ public:
 	/// @param outname is the name of the output BIP file to be saved to disk.
 	bool bip(std::string outname)
 	{
-		unsigned int S = R[0] * R[2] * sizeof(T);		//calculate the number of bytes in a ZX slice
+		unsigned int S = X() * Z() * 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";
@@ -429,14 +413,14 @@ public:
 		T * q;			//pointer to the current ZX slice for bip file
 		q = (T*)malloc(S);
-		for ( unsigned i = 0; i < R[1]; i++)
+		for ( unsigned i = 0; i < Y(); i++)
 		{			
-			getY(p, i);
-			for ( unsigned k = 0; k < R[2]; k++)
+			read_plane_y(p, i);
+			for ( unsigned k = 0; k < Z(); k++)
 			{
-				unsigned ks = k * R[0];
-				for ( unsigned j = 0; j < R[0]; j++)
-					q[k + j * R[2]] = p[ks + j];
+				unsigned ks = k * X();
+				for ( unsigned j = 0; j < X(); j++)
+					q[k + j * Z()] = p[ks + j];
 			}
 			target.write(reinterpret_cast<const char*>(q), S);   //write a band data into target file	
@@ -460,7 +444,7 @@ public:
 	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size.
 	bool baseline_band(double lb, double rb, T* lp, T* rp, double wavelength, T* result){
-		unsigned XY = R[0] * R[1];
+		unsigned XY = X() * Y();
 		band(result, wavelength);		//get band
 		//perform the baseline correction
@@ -481,7 +465,7 @@ public:
 		T* lp;
 		T* rp;
-		unsigned XY = R[0] * R[1];
+		unsigned XY = X() * Y();
 		unsigned S = XY * sizeof(T);
 		lp = (T*) malloc(S);			//memory allocation
 		rp = (T*) malloc(S);
@@ -512,7 +496,7 @@ public:
 		T* cur;		//current band 1
 		T* cur2;	//current band 2
-		unsigned XY = R[0] * R[1];
+		unsigned XY = X() * Y();
 		unsigned S = XY * sizeof(T);
 		lp = (T*) malloc(S);			//memory allocation
@@ -593,14 +577,14 @@ public:
 	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
 	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));	
+		T* p1 = (T*)malloc(X() * Y() * sizeof(T));
+		T* p2 = (T*)malloc(X() * Y() * 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++){
+		for(unsigned i = 0; i < X() * Y(); i++){
 			if(p1[i] == 0 && p2[i] ==0)
 				result[i] = 1;
 			else
@@ -624,14 +608,14 @@ public:
 	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));	
+		T* p1 = (T*)malloc(X() * Y() * sizeof(T));
+		T* p2 = (T*)malloc(X() * Y() * 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++){
+		for(unsigned i = 0; i < X() * Y(); i++){
 			if(p1[i] == 0 && p2[i] ==0)
 				result[i] = 1;
 			else
@@ -657,14 +641,14 @@ public:
 	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));	
+		T* p1 = (T*)malloc(X() * Y() * sizeof(T));
+		T* p2 = (T*)malloc(X() * Y() * 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++){
+		for(unsigned i = 0; i < X() * Y(); i++){
 			if(p1[i] == 0 && p2[i] ==0)
 				result[i] = 1;
 			else
@@ -689,7 +673,7 @@ public:
 		T* cur;		//current band 1
 		T* cur2;	//current band 2
-		unsigned XY = R[0] * R[1];
+		unsigned XY = X() * Y();
 		unsigned S = XY * sizeof(T);
 		lp = (T*) malloc(S);			//memory allocation
@@ -765,14 +749,14 @@ public:
 	/// @param rab is the label for the end of the peak
 	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
 	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));	
+		T* p1 = (T*)malloc(X() * Y() * sizeof(T));
+		T* p2 = (T*)malloc(X() * Y() * 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++){
+		for(unsigned i = 0; i < X() * Y(); i++){
 			if(p1[i] == 0 && p2[i] ==0)
 				result[i] = 1;
 			else
@@ -793,10 +777,10 @@ public:
 	/// @param p is a pointer to a pre-allocated array at least X * Y in size
 	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
+		T* temp = (T*)malloc(X() * Y() * sizeof(T));		//allocate memory for the certain band
 		band(temp, mask_band);
-		for (unsigned int i = 0; i < R[0] * R[1]; i++) {
+		for (unsigned int i = 0; i < X() * Y(); i++) {
 				if (temp[i] < threshold)
 					p[i] = 0;
 				else
@@ -816,20 +800,20 @@ public:
 		std::ofstream target(outfile.c_str(), std::ios::binary);
 		//I THINK THIS IS WRONG
-		unsigned XZ = R[0] * R[2];		//calculate the number of values in a page on disk
+		unsigned XZ = X() * Z();		//calculate the number of values in a page on disk
 		unsigned L = XZ * sizeof(T);	//calculate the size of the page (in bytes)
 		T * temp = (T*)malloc(L);		//allocate memory for a temporary page
-		for (unsigned i = 0; i < R[1]; i++)			//for each value in R[1] (BIP should be X)
+		for (unsigned i = 0; i < Y(); i++)			//for each value in Y() (BIP should be X)
 		{
-			getY(temp, i);							//retrieve an ZX slice, stored in temp
-			for ( unsigned j = 0; j < R[2]; j++)	//for each R[2] (Y)
+			read_plane_y(temp, i);							//retrieve an ZX slice, stored in temp
+			for ( unsigned j = 0; j < Z(); j++)	//for each Z() (Y)
 			{
-				for (unsigned k = 0; k < R[0]; k++) //for each band
+				for (unsigned k = 0; k < X(); k++) //for each band
 				{
-				if(p[i * R[0] + k] == 0)
-					temp[j * R[0] + k] = 0;
+				if(p[i * X() + k] == 0)
+					temp[j * X() + k] = 0;
 				else
 					continue;
 				}
@@ -843,30 +827,29 @@ public:
 	///Saves to disk only those spectra corresponding to mask values != 0
 	bool sift_mask(std::string outfile, unsigned char* p){
-		// Assume R[0] = X, R[1] = Y, R[2] = Z.
+		// Assume X() = X, Y() = Y, Z() = Z.
 		std::ofstream target(outfile.c_str(), std::ios::binary);
 		//for loading pages:
-		unsigned XZ = R[0] * R[2];		//calculate the number of values in an XZ page on disk
+		unsigned XZ = X() * Z();		//calculate the number of values in an XZ page on disk
 		unsigned L = XZ * sizeof(T);	//calculate the size of the page (in bytes)
 		T * temp = (T*)malloc(L);		//allocate memory for a temporary page
 		//for saving spectra:
-		unsigned Z = R[2];						//calculate the number of values in a spectrum
-		unsigned LZ = Z * sizeof(T);			//calculate the size of the spectrum (in bytes)
+		unsigned LZ = Z() * sizeof(T);			//calculate the size of the spectrum (in bytes)
 		T * tempZ = (T*)malloc(LZ);				//allocate memory for a temporary spectrum
-		spectrum(tempZ, R[0] - 1, R[1] - 1);	//creates a dummy spectrum by taking the last spectrum in the image
+		spectrum(tempZ, X() - 1, Y() - 1);	//creates a dummy spectrum by taking the last spectrum in the image
-		for (unsigned i = 0; i < R[1]; i++)			//Select a page by choosing Y coordinate, R[1]
+		for (unsigned i = 0; i < Y(); i++)			//Select a page by choosing Y coordinate, Y()
 		{
-			getY(temp, i);							//retrieve an ZX page, store in "temp"
-			for (unsigned j = 0; j < R[0]; j++)		//Select a pixel by choosing X coordinate in the page, R[0]
+			read_plane_y(temp, i);							//retrieve an ZX page, store in "temp"
+			for (unsigned j = 0; j < X(); j++)		//Select a pixel by choosing X coordinate in the page, X()
 			{
-				if (p[j * R[0] + i] != 0)					//if the mask != 0 at that XY pixel
+				if (p[j * X() + i] != 0)					//if the mask != 0 at that XY pixel
 				{
-					for (unsigned k = 0; k < R[2]; k++)		//Select a voxel by choosing Z coordinate at the pixel
+					for (unsigned k = 0; k < Z(); k++)		//Select a voxel by choosing Z coordinate at the pixel
 					{
-						tempZ[k] = temp[k*R[0] + i];		//Pass the correct spectral value from XZ page into the spectrum to be saved. 
+						tempZ[k] = temp[k*X() + i];		//Pass the correct spectral value from XZ page into the spectrum to be saved.
 					}
 					target.write(reinterpret_cast<const char*>(tempZ), LZ);		//write that spectrum to disk. Size is L2.
 				}
@@ -883,25 +866,25 @@ public:
 	/// @param p is a pointer to memory of size X * Y * sizeof(T) that will store the band averages.
 	bool band_avg(T* p){
-		unsigned long long XZ = R[0] * R[2];
+		unsigned long long XZ = X() * Z();
 		T* temp = (T*)malloc(sizeof(T) * XZ);
-		T* line = (T*)malloc(sizeof(T) * R[0]);
+		T* line = (T*)malloc(sizeof(T) * X());
-		for (unsigned i = 0; i < R[1]; i++){
+		for (unsigned i = 0; i < Y(); i++){
 			getY(temp, i);
 			//initialize x-line
-			for (unsigned j = 0; j < R[0]; j++){
+			for (unsigned j = 0; j < X(); j++){
 				line[j] = 0;
 			}
 			unsigned c = 0;
-			for (unsigned j = 0; j < R[2]; j++){
-				for (unsigned k = 0; k < R[0]; k++){
-					line[k] += temp[c] / (T)R[2];
+			for (unsigned j = 0; j < Z(); j++){
+				for (unsigned k = 0; k < X(); k++){
+					line[k] += temp[c] / (T)Z();
 					c++;
 				}
 			}
-			for (unsigned j = 0; j < R[0]; j++){
-				p[j + i * R[0]] = line[j];
+			for (unsigned j = 0; j < X(); j++){
+				p[j + i * X()] = line[j];
 			}
 		}
 		free(temp);
@@ -913,10 +896,10 @@ public:
 	/// @param p is a pointer to pre-allocated memory of size [B * sizeof(T)] that stores the mean spectrum
 	/// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
 	bool avg_band(T*p, unsigned char* mask){
-		unsigned long long XZ = R[0] * R[2];
-		unsigned long long XY = R[0] * R[1];
+		unsigned long long XZ = X() * Z();
+		unsigned long long XY = X() * Y();
 		T* temp = (T*)malloc(sizeof(T) * XZ);
-		for (unsigned j = 0; j < R[2]; j++){
+		for (unsigned j = 0; j < Z(); j++){
 			p[j] = 0;
 		}
 		//calculate vaild number in a band
@@ -926,13 +909,13 @@ public:
 				count++;
 			}
 		}
-		for (unsigned k = 0; k < R[1]; k++){
-			getY(temp, k);
-			unsigned kx = k * R[0];
-			for (unsigned i = 0; i < R[0]; i++){
+		for (unsigned k = 0; k < Y(); k++){
+			read_plane_y(temp, k);
+			unsigned kx = k * X();
+			for (unsigned i = 0; i < X(); i++){
 				if (mask[kx + i] != 0){
-					for (unsigned j = 0; j < R[2]; j++){
-						p[j] += temp[j * R[0] + i] / (T)count;
+					for (unsigned j = 0; j < Z(); j++){
+						p[j] += temp[j * X() + i] / (T)count;
 					}
 				}
 			}
@@ -948,8 +931,8 @@ public:
 	/// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
 	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];
+		unsigned long long xy = X() * Y();
+		unsigned int B = Z();
 		T* temp = (T*)malloc(sizeof(T) * B);
 		//count vaild pixels in a band
 		unsigned count = 0;
@@ -999,16 +982,16 @@ public:
 		//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 * R[2] * sizeof(T);
+		unsigned long long L = sam * Z() * 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] - sam) * sizeof(T);		//jump pointer to the next band
+		unsigned long long jumpb = (X() - sam) * sizeof(T);		//jump pointer to the next band
 		//get start
-		file.seekg((y0 * R[0] * R[2] + x0) * sizeof(T), std::ios::beg);
+		file.seekg((y0 * X() * Z() + x0) * sizeof(T), std::ios::beg);
 		for (unsigned i = 0; i < lin; i++)
 		{
-			for (unsigned j = 0; j < R[2]; j++)
+			for (unsigned j = 0; j < Z(); j++)
 			{
 				file.read((char *)(temp + j * sam), sizeof(T) * sam);
 				file.seekg(jumpb, std::ios::cur);    //go to the next band	
@@ -81,6 +81,8 @@ protected:
 		return false;
 	}
+
+
 public:
 	/// Open a binary file for streaming.
@@ -159,8 +161,141 @@ public:
 		return true;
 	}
+
+
+	///Reads a line Z (slowest dimension) for a given XY value
+
+	/// @param p is a pointer to pre-allocated memory equal to the line size R[2]
+	/// @param x is the x coordinate
+	/// @param y is the y coordinate
+	bool read_line_01( T * p, unsigned int x, unsigned int 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;
+	}
+
+	///Reads a line X (fastest dimension) for a given YZ value
+
+	/// @param p is a pointer to pre-allocated memory equal to the line size R[2]
+	/// @param x is the y coordinate
+	/// @param y is the z coordinate
+	bool read_line_12(T * p, unsigned int y, unsigned int z){
+		//test to make sure the specified value is within range
+		if( y >= R[1] || z >= R[2] ){
+			std::cout<<"ERROR: sample or line out of range"<<std::endl;
+			return false;
+		}
+
+		file.seekg((z * R[0] * R[1] + y * R[0]) * sizeof(T), std::ios::beg);	//seek to the start of the line
+		file.read((char *)p, sizeof(T) * R[0]);									//read the line
+		
+		return true;
+	}
+
+	///Reads a line Y (second fastest dimension) for a given XZ value
+
+	/// @param p is a pointer to pre-allocated memory equal to the line size R[2]
+	/// @param x is the y coordinate
+	/// @param z is the z coordinate
+	bool read_line_02(T * p, unsigned int x, unsigned int z){
+		//test to make sure the specified value is within range
+		if( x >= R[0] || z >= R[2] ){
+			std::cout<<"ERROR: sample or line out of range"<<std::endl;
+			return false;
+		}
+
+		file.seekg((z * R[0] * R[1] + x) * sizeof(T), std::ios::beg);           //seek to the start of the line
+		for (unsigned int i = 0; i < R[1]; i++){									//for each pixel in the line
+			file.read((char *)(p + i), sizeof(T));					//read the pixel
+			file.seekg((R[0] - 1) * sizeof(T), std::ios::cur);		//seek to the next pixel in the line
+		}
+
+		return true;
+	}
+
+	bool read_plane_0(T* p, unsigned int n){
+
+		if (n >= R[0]){										//make sure the number is within the possible range
+			std::cout<<"ERROR read_plane_0: page out of range"<<std::endl;
+			return false;
+		}
+		unsigned int jump = (R[0] - 1) * sizeof(T);		//number of bytes to skip between samples
+
+		//seek to the start of the plane
+		file.seekg(n * sizeof(T), std::ios::beg);
+
+		unsigned int N = R[1] * R[2];
+		for(unsigned int i = 0; i<N; i++){
+			file.read((char*)(p+i), sizeof(T));
+			file.seekg(jump, std::ios::cur);
+		}
+
+		return true;
+
+
+	}
+
+	bool read_plane_1(T* p, unsigned int n){
+
+		unsigned int L = R[0] * sizeof(T);		//caculate the number of bytes in a sample line
+		unsigned int jump = R[0] * (R[1] - 1) * sizeof(T);
+
+		if (n >= R[1]){										//make sure the bank number is right
+			std::cout<<"ERROR read_plane_1: page out of range"<<std::endl;
+			return false;
+		}
+
+		file.seekg(R[0] * n * sizeof(T), std::ios::beg);
+		for (unsigned i = 0; i < R[2]; i++){
+
+			file.read((char *)(p + i * R[0]), L);
+			file.seekg( jump, std::ios::cur);
+			std::cout<<i<<"    ";
+		}
+
+		return true;
+	}
+
+	bool read_plane_2(T* p, unsigned int n){
+		return read_page(p, n);
+	}
+
+	bool read_pixel(T* p, unsigned int i){
+		if(i >= R[0] * R[1] * R[2]){
+			std::cout<<"ERROR read_pixel: n is out of range"<<std::endl;
+			return false;
+		}
+
+		file.seekg(i * sizeof(T), std::ios::cur);
+		file.read((char*)p, sizeof(T));
+
+	}
+
+	bool read_pixel(T* p, unsigned int x, unsigned int y, unsigned int z){
+
+		if(x < 0 || x >= R[0] || y < 0 || y >= R[1] || z < 0 || z > R[2]){
+			std::cout<<"ERROR read_pixel: (x,y,z) is out of range"<<std::endl;
+			return false;
+		}
+
+		unsigned int i = z * R[0] * R[1] + y * R[0] + z;
+		return read_pixel(p, i);
+	}
+
 	//saves a hyperplane orthogonal to dimension d at intersection n
-	bool read_plane(T * dest, unsigned int d, unsigned int n){
+	/*bool read_plane(T * dest, unsigned int d, unsigned int n){
 		//reset the file pointer back to the beginning of the file
 		file.seekg(0, std::ios::beg);
@@ -184,10 +319,10 @@ public:
 		return true;
-	}
+	}*/
 	//save one pixel of the file into the memory, and return the pointer
-	bool read_spectrum(T * p, unsigned x, unsigned y){
+	/*bool read_spectrum(T * p, unsigned x, unsigned y){
 		unsigned int i;
@@ -204,7 +339,7 @@ public:
 		}
 		return true;	
-	}
+	}*/
 };
@@ -27,11 +27,22 @@ protected:
 	std::vector<double> w; //band wavelength
 	unsigned int offset;		//header offset
+	unsigned long X(){
+		return R[1];
+	}
+	unsigned long Y(){
+		return R[2];
+	}
+	unsigned long Z(){
+		return R[0];
+	}
+
 public:
 	using binary<T>::open;
 	using binary<T>::file;
 	using binary<T>::R;
+	using binary<T>::read_line_12;
 	/// Open a data file for reading using the class interface.
@@ -48,7 +59,7 @@ public:
 		//copy the offset to the structure
 		offset = header_offset;
-		return open(filename, vec<unsigned int>(X, Y, B), header_offset);
+		return open(filename, vec<unsigned int>(B, X, Y), header_offset);
 	}
@@ -57,21 +68,7 @@ public:
 	/// @param p is a pointer to an allocated region of memory at least X * Y * sizeof(T) in size.
 	/// @param page <= B is the integer number of the band to be copied.
 	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;
+		return binary<T>::read_plane_0(p, page);
 	}
 	/// Retrieve a single band (by numerical label) and stores it in pre-allocated memory.
@@ -84,7 +81,7 @@ public:
 		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 int XY = X() * Y();	//calculate the number of pixels in a band
 		unsigned page=0;                      //bands around the wavelength
@@ -101,8 +98,8 @@ public:
 		{
 			page++;
 			//if wavelength is larger than the last wavelength in header file
-			if (page == R[2]) {
-				band_index(p, R[2]-1);
+			if (page == Z()) {
+				band_index(p, Z()-1);
 				return true;
 			}
 		}
@@ -128,11 +125,24 @@ public:
 		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)
+
+	/// Retrieve a single spectrum (Z-axis line) at a given (x, y) location and stores it in pre-allocated memory.
+
+	/// @param p is a pointer to pre-allocated memory at least B * sizeof(T) in size.
+	/// @param x is the x-coordinate (dimension 1) of the spectrum.
+	/// @param y is the y-coordinate (dimension 2) of the spectrum.
+	bool spectrum(T * p, unsigned x, unsigned y){
+		return read_line_12(p, x, y);				//read a line in the binary YZ plane (dimension order for BIP is ZXY)
+	}
+
+	/// Retrieves a band of x values from a given xz plane.
+
+	/// @param p is a pointer to pre-allocated memory at least X * sizeof(T) in size
+	/// @param c is a pointer to an existing XZ plane (size X*Z*sizeof(T))
+	/// @param wavelength is the wavelength of X values to retrieve
+	bool read_x_from_xz(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 int B = Z();
 		unsigned page=0;                      //samples around the wavelength
@@ -141,7 +151,7 @@ public:
 		//if wavelength is smaller than the first one in header file
 		if ( w[page] > wavelength ){
-			for(unsigned j = 0; j < R[0]; j++)
+			for(unsigned j = 0; j < X(); j++)
 			{
 				p[j] = c[j * B];
 			}		
@@ -153,7 +163,7 @@ public:
 			page++;
 			//if wavelength is larger than the last wavelength in header file
 			if (page == B) {
-				for(unsigned j = 0; j < R[0]; j++)
+				for(unsigned j = 0; j < X(); j++)
 				{
 					p[j] = c[(j + 1) * B - 1];
 				}
@@ -164,20 +174,20 @@ public:
 		{
 			T * p1;
 			T * p2;
-			p1=(T*)malloc( X * sizeof(T));                     //memory allocation
-			p2=(T*)malloc( X * sizeof(T));
+			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++)
+			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++)
+			for(unsigned j = 0; j < X(); j++)
 			{
 				p2[j] = c[j * B + page];
 			}
-			for(unsigned i=0; i < X; i++){
+			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];
 			}
@@ -187,7 +197,7 @@ public:
 		else                           //if the wavelength is equal to a wavelength in header file
 		{
 			//band_index(p, page);
-			for(unsigned j = 0; j < X; j++)
+			for(unsigned j = 0; j < X(); j++)
 			{
 				p[j] = c[j * B + page];
 			}
@@ -195,98 +205,6 @@ public:
 		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;
-	}
-
-	/// Retrieve a single spectrum (B-axis line) at a given (x, y) location and stores it in pre-allocated memory.
-
-	/// @param p is a pointer to pre-allocated memory at least B * sizeof(T) in size.
-	/// @param x is the x-coordinate (dimension 1) of the spectrum.
-	/// @param y is the y-coordinate (dimension 2) of the spectrum.
-	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;	
-	}
 	/// Retrieve a single pixel and stores it in pre-allocated memory.
@@ -294,30 +212,20 @@ public:
 	/// @param n is an integer index to the pixel using linear array indexing.
 	bool pixel(T * p, unsigned n){
-		unsigned bandnum = R[0] * R[1];		//calculate numbers in one band
+		unsigned bandnum = X() * Y();		//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;	
+		file.seekg(n * Z() * sizeof(T), std::ios::beg);           //point to the certain pixel
+		file.read((char *)p, sizeof(T) * Z());
+		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;
-		
+	bool read_plane_y(T * p, unsigned y){
+		return binary<T>::read_plane_2(p, y);
 	}
 	/// Perform baseline correction given a list of baseline points and stores the result in a new BSQ file.
@@ -332,10 +240,9 @@ public:
 		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 ZX = Z() * X();		//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];
+		unsigned int B = Z();
 		T* c;			//pointer to the current Y slice
 		c = (T*)malloc(L);  //memory allocation
@@ -343,8 +250,8 @@ public:
 		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));
+		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
@@ -357,10 +264,10 @@ public:
 		}
 	//	loop start	correct every y slice
-		for (unsigned k =0; k < R[1]; k++)
+		for (unsigned k =0; k < Y(); k++)
 		{
 			//get the current y slice
-			getY(c, k);
+			read_plane_y(c, k);
 			//initialize lownum, highnum, low, high		
 			control=0;
@@ -369,16 +276,16 @@ public:
 			//set the baseline point at band 0 to 0
 			if(wls[0] != w[0]){
 				bi = wls[control];			
-				memset(a, (char)0, X * sizeof(T) );
+				memset(a, (char)0, X() * sizeof(T) );
 			}
 			//else get the low band
 			else{
 				control++;
-				getYZ(a, c, ai);
+				read_x_from_xz(a, c, ai);
 				bi = wls[control];
 			}
 			//get the high band
-			getYZ(b, c, bi);
+			read_x_from_xz(b, c, bi);
 			//correct every YZ line
@@ -394,7 +301,7 @@ public:
 						ai = bi;
 						bi = wls[control];
-						getYZ(b, c, bi);
+						read_x_from_xz(b, c, bi);
 					}
 					//if the last BL point on the last band of the file?
@@ -402,7 +309,7 @@ public:
 						std::swap(a, b);	//swap the baseline band pointers
-						memset(b, (char)0, X * sizeof(T) );	//clear the high band
+						memset(b, (char)0, X() * sizeof(T) );	//clear the high band
 						ai = bi;
 						bi = w[B - 1];
@@ -412,7 +319,7 @@ public:
 				ci = w[cii];
 				//perform the baseline correction
-				for(unsigned i=0; i < X; i++)
+				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] );
@@ -439,11 +346,9 @@ public:
 	///	@param t is a threshold specified such that a spectrum with a value at w less than t is set to zero. Setting this threshold allows the user to limit division by extremely small numbers.
 	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 B = Z();		//calculate the number of bands
+		unsigned int ZX = Z() * X();
+		unsigned int XY = X() * Y();	//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);
@@ -458,11 +363,11 @@ public:
 		band(b, w);             //get the certain band into memory
-		for(unsigned j = 0; j < Y; j++)
+		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++)
+			read_plane_y(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++)
@@ -494,7 +399,7 @@ public:
 		bil(temp);
 		stim::bil<T> n;
-		if(n.open(temp, R[0], R[1], R[2], offset, w)==false){        //open infile
+		if(n.open(temp, X(), Y(), Z(), offset, w)==false){        //open infile
 			std::cout<<"ERROR: unable to open input file"<<std::endl;
 			exit(1);
 		}
@@ -511,7 +416,7 @@ public:
 	/// @param outname is the name of the output BIL file to be saved to disk.
 	bool bil(std::string outname)
 	{
-		unsigned int S = R[0] * R[2] * sizeof(T);		//calculate the number of bytes in a ZX slice
+		unsigned int S = X() * Z() * 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";
@@ -521,14 +426,14 @@ public:
 		T * q;			//pointer to the current XZ slice for bil file
 		q = (T*)malloc(S);
-		for ( unsigned i = 0; i < R[1]; i++)
+		for ( unsigned i = 0; i < Y(); i++)
 		{			
-			getY(p, i);
-			for ( unsigned k = 0; k < R[2]; k++)
+			read_plane_y(p, i);
+			for ( unsigned k = 0; k < Z(); k++)
 			{
-				unsigned ks = k * R[0];
-				for ( unsigned j = 0; j < R[0]; j++)
-					q[ks + j] = p[k + j * R[2]];
+				unsigned ks = k * X();
+				for ( unsigned j = 0; j < X(); j++)
+					q[ks + j] = p[k + j * Z()];
 			}
 			target.write(reinterpret_cast<const char*>(q), S);   //write a band data into target file	
 		}
@@ -550,7 +455,7 @@ public:
 	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size.
 	bool baseline_band(double lb, double rb, T* lp, T* rp, double wavelength, T* result){
-		unsigned XY = R[0] * R[1];
+		unsigned XY = X() * Y();
 		band(result, wavelength);		//get band
 		//perform the baseline correction
@@ -571,7 +476,7 @@ public:
 		T* lp;
 		T* rp;
-		unsigned XY = R[0] * R[1];
+		unsigned XY = X() * Y();
 		unsigned S = XY * sizeof(T);
 		lp = (T*) malloc(S);			//memory allocation
 		rp = (T*) malloc(S);
@@ -601,7 +506,7 @@ public:
 		T* cur;		//current band 1
 		T* cur2;	//current band 2
-		unsigned XY = R[0] * R[1];
+		unsigned XY = X() * Y();
 		unsigned S = XY * sizeof(T);
 		lp = (T*) malloc(S);			//memory allocation
@@ -682,14 +587,14 @@ public:
 	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
 	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));	
+		T* p1 = (T*)malloc(X() * Y() * sizeof(T));	
+		T* p2 = (T*)malloc(X() * Y() * 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++){
+		for(unsigned i = 0; i < X() * Y(); i++){
 			if(p1[i] == 0 && p2[i] ==0)
 				result[i] = 1;
 			else
@@ -713,14 +618,14 @@ public:
 	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));	
+		T* p1 = (T*)malloc(X() * Y() * sizeof(T));	
+		T* p2 = (T*)malloc(X() * Y() * 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++){
+		for(unsigned i = 0; i < X() * Y(); i++){
 			if(p1[i] == 0 && p2[i] ==0)
 				result[i] = 1;
 			else
@@ -746,14 +651,14 @@ public:
 	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));	
+		T* p1 = (T*)malloc(X() * Y() * sizeof(T));	
+		T* p2 = (T*)malloc(X() * Y() * 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++){
+		for(unsigned i = 0; i < X() * Y(); i++){
 			if(p1[i] == 0 && p2[i] ==0)
 				result[i] = 1;
 			else
@@ -778,7 +683,7 @@ public:
 		T* cur;		//current band 1
 		T* cur2;	//current band 2
-		unsigned XY = R[0] * R[1];
+		unsigned XY = X() * Y();
 		unsigned S = XY * sizeof(T);
 		lp = (T*) malloc(S);			//memory allocation
@@ -854,14 +759,14 @@ public:
 	/// @param rab is the label for the end of the peak
 	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
 	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));	
+		T* p1 = (T*)malloc(X() * Y() * sizeof(T));	
+		T* p2 = (T*)malloc(X() * Y() * 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++){
+		for(unsigned i = 0; i < X() * Y(); i++){
 			if(p1[i] == 0 && p2[i] ==0)
 				result[i] = 1;
 			else
@@ -882,10 +787,10 @@ public:
 	/// @param p is a pointer to a pre-allocated array at least X * Y in size
 	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
+		T* temp = (T*)malloc(X() * Y() * sizeof(T));		//allocate memory for the certain band
 		band(temp, mask_band);
-		for (unsigned int i = 0; i < R[0] * R[1];i++) {
+		for (unsigned int i = 0; i < X() * Y();i++) {
 				if (temp[i] < threshold)
 					p[i] = 0;
 				else
@@ -905,20 +810,20 @@ public:
 		std::ofstream target(outfile.c_str(), std::ios::binary);
-		unsigned ZX = R[2] * R[0];		//calculate the number of values in a page (XZ in BIP)
+		unsigned ZX = Z() * X();		//calculate the number of values in a page (XZ in BIP)
 		unsigned L = ZX * sizeof(T);	//calculate the number of bytes in a page
 		T * temp = (T*)malloc(L);		//allocate space for that page
-		for (unsigned i = 0; i < R[1]; i++)			//for each page (Y in BIP)
+		for (unsigned i = 0; i < Y(); i++)			//for each page (Y in BIP)
 		{
-			getY(temp, i);							//load that page (it's pointed to by temp)
-			for ( unsigned j = 0; j < R[0]; j++)	//for each X value
+			read_plane_y(temp, i);							//load that page (it's pointed to by temp)
+			for ( unsigned j = 0; j < X(); j++)	//for each X value
 			{
-				for (unsigned k = 0; k < R[2]; k++)	//for each B value (band)
+				for (unsigned k = 0; k < Z(); k++)	//for each B value (band)
 				{
-					if (p[i * R[0] + j] == 0)	//if the mask value is zero
-					temp[j * R[2] + k] = 0;			//set the pixel value to zero
+					if (p[i * X() + j] == 0)	//if the mask value is zero
+					temp[j * Z() + k] = 0;			//set the pixel value to zero
 				else								//otherwise just continue
 					continue;
 				}
@@ -933,30 +838,29 @@ public:
 	///Saves to disk only those spectra corresponding to mask values != 0
 	bool sift_mask(std::string outfile, unsigned char* p){
-		// Assume R[0] = X, R[1] = Y, R[2] = Z.
+		// Assume X() = X, Y() = Y, Z() = Z.
 		std::ofstream target(outfile.c_str(), std::ios::binary);
 		//for loading pages:
-		unsigned ZX = R[0] * R[2];		//calculate the number of values in an XZ page on disk
+		unsigned ZX = X() * Z();		//calculate the number of values in an XZ page on disk
 		unsigned L = ZX * sizeof(T);	//calculate the size of the page (in bytes)
 		T * temp = (T*)malloc(L);		//allocate memory for a temporary page
 		//for saving spectra:
-		unsigned Z = R[2];						//calculate the number of values in a spectrum
-		unsigned LZ = Z * sizeof(T);			//calculate the size of the spectrum (in bytes)
+		unsigned LZ = Z() * sizeof(T);			//calculate the size of the spectrum (in bytes)
 		T * tempZ = (T*)malloc(LZ);				//allocate memory for a temporary spectrum
-		spectrum(tempZ, R[0] - 1, R[1] - 1);	//creates a dummy spectrum by taking the last spectrum in the image
+		spectrum(tempZ, X() - 1, Y() - 1);	//creates a dummy spectrum by taking the last spectrum in the image
-		for (unsigned i = 0; i < R[1]; i++)			//Select a page by choosing Y coordinate, R[1]
+		for (unsigned i = 0; i < Y(); i++)			//Select a page by choosing Y coordinate, Y()
 		{
-			getY(temp, i);							//retrieve an ZX page, store in "temp"
-			for (unsigned j = 0; j < R[0]; j++)		//Select a pixel by choosing X coordinate in the page, R[0]
+			read_plane_y(temp, i);							//retrieve an ZX page, store in "temp"
+			for (unsigned j = 0; j < X(); j++)		//Select a pixel by choosing X coordinate in the page, X()
 			{
-				if (p[j * R[0] + i] != 0)					//if the mask != 0 at that XY pixel
+				if (p[j * X() + i] != 0)					//if the mask != 0 at that XY pixel
 				{
-					for (unsigned k = 0; k < R[2]; k++)		//Select a voxel by choosing Z coordinate at the pixel
+					for (unsigned k = 0; k < Z(); k++)		//Select a voxel by choosing Z coordinate at the pixel
 					{
-						tempZ[k] = temp[i*R[2] + k];		//Pass the correct spectral value from XZ page into the spectrum to be saved. 
+						tempZ[k] = temp[i*Z() + k];		//Pass the correct spectral value from XZ page into the spectrum to be saved. 
 					}
 					target.write(reinterpret_cast<const char*>(tempZ), LZ);		//write that spectrum to disk. Size is L2.
 				}
@@ -973,16 +877,16 @@ public:
 	/// @param p is a pointer to memory of size X * Y * sizeof(T) that will store the band averages.
 	bool band_avg(T* p){
-		unsigned long long XY = R[0] * R[1];
+		unsigned long long XY = X() * Y();
 		//get every pixel and calculate average value
-		T* temp = (T*)malloc(sizeof(T) * R[2]);
+		T* temp = (T*)malloc(sizeof(T) * Z());
 		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];
+			for (unsigned j = 0; j < Z(); j++){
+				sum += temp[j]/(T)Z();
 			}
 			p[i] = sum;
 		}
@@ -995,10 +899,10 @@ public:
 	/// @param p is a pointer to pre-allocated memory of size [B * sizeof(T)] that stores the mean spectrum
 	/// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
 	bool avg_band(T*p, unsigned char* mask){
-		unsigned long long XY = R[0] * R[1];
-		T* temp = (T*)malloc(sizeof(T) * R[2]);
+		unsigned long long XY = X() * Y();
+		T* temp = (T*)malloc(sizeof(T) * Z());
 		//Iinitialize
-		for (unsigned j = 0; j < R[2]; j++){
+		for (unsigned j = 0; j < Z(); j++){
 			p[j] = 0;
 		}
 		//calculate vaild number in a band
@@ -1012,7 +916,7 @@ public:
 		for (unsigned i = 0; i < XY; i++){
 			if (mask[i] != 0){			
 				pixel(temp, i);
-				for (unsigned j = 0; j < R[2]; j++){
+				for (unsigned j = 0; j < Z(); j++){
 					p[j] += temp[j] / (T)count;
 				}
 			}
@@ -1028,8 +932,8 @@ public:
 	/// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
 	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];
+		unsigned long long xy = X() * Y();
+		unsigned int B = Z();
 		T* temp = (T*)malloc(sizeof(T) * B);
 		//count vaild pixels in a band
 		unsigned count = 0;
@@ -1079,17 +983,17 @@ public:
 		//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);
+		unsigned long long L = Z() * 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
+		unsigned long long sp = y0 * X() + x0;		//start pixel
 		for (unsigned i = 0; i < lin; i++)
 		{
 			for (unsigned j = 0; j < sam; j++)
 			{
-				pixel(temp, sp + j + i * R[0]);
+				pixel(temp, sp + j + i * X());
 				out.write(reinterpret_cast<const char*>(temp), L);   //write slice data into target file	
 			}
 		}
@@ -30,12 +30,26 @@ protected:
 	std::vector<double> w;	//band wavelengths
 	unsigned int offset;
+	using binary<T>::R;
+
+	unsigned long X(){
+		return R[0];
+	}
+	unsigned long Y(){
+		return R[1];
+	}
+	unsigned long Z(){
+		return R[2];
+	}
+
 public:
 	using binary<T>::open;
 	using binary<T>::file;
+	using binary<T>::read_line_01;
+	using binary<T>::read_plane_2;
 	//using binary<T>::getSlice;
-	using binary<T>::R;
+	
 	/// Open a data file for reading using the class interface.
@@ -61,17 +75,7 @@ public:
 	/// @param p is a pointer to an allocated region of memory at least X * Y * sizeof(T) in size.
 	/// @param page <= B is the integer number of the band to be copied.
 	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;
-		}
-		//move the file pointer to the start of the page on disk
-		file.seekg(R[0] * R[1] * page * sizeof(T));
-		//read the page from disk
-		file.read((char *)p, sizeof(T) * R[0] * R[1]);
-
-		return true;
+		return read_plane_2(p, page);
 	}
 	/// Retrieve a single band (by numerical label) and stores it in pre-allocated memory.
@@ -84,7 +88,7 @@ public:
 		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 int XY = X() * Y();	//calculate the number of pixels in a band
 		unsigned page = 0;                
@@ -101,8 +105,8 @@ public:
 			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
+			if (page == Z()) {
+				band_index(p, Z()-1);		//return the last band
 				return true;
 			}
 		}
@@ -131,28 +135,13 @@ public:
 		return true;
 	}
-	/// Retrieve a single spectrum (B-axis line) at a given (x, y) location and stores it in pre-allocated memory.
+	/// Retrieve a single spectrum (Z-axis line) at a given (x, y) location and stores it in pre-allocated memory.
 	/// @param p is a pointer to pre-allocated memory at least B * sizeof(T) in size.
 	/// @param x is the x-coordinate (dimension 1) of the spectrum.
 	/// @param y is the y-coordinate (dimension 2) of the spectrum.
 	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;	
+		return read_line_01(p, x, y);
 	}
 	/// Retrieve a single pixel and stores it in pre-allocated memory.
@@ -161,14 +150,14 @@ public:
 	/// @param n is an integer index to the pixel using linear array indexing.
 	bool pixel(T * p, unsigned n){
-		unsigned bandnum = R[0] * R[1];		//calculate numbers in one band
+		unsigned bandnum = X() * Y();		//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++)
+		for (unsigned i = 0; i < Z(); i++)
 		{
 			file.read((char *)(p + i), sizeof(T));
 			file.seekg((bandnum - 1) * sizeof(T), std::ios::cur);    //go to the next band
@@ -189,8 +178,8 @@ public:
 		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 B = Z();		//calculate the number of bands
+		unsigned int XY = X() * Y();	//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
@@ -289,8 +278,8 @@ public:
 	///	@param t is a threshold specified such that a spectrum with a value at w less than t is set to zero. Setting this threshold allows the user to limit division by extremely small numbers.
 	bool normalize(std::string outname, double w, double t = 0.0)
 	{
-		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 B = Z();		//calculate the number of bands
+		unsigned int XY = X() * Y();	//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
@@ -336,7 +325,7 @@ public:
 		bil(temp);
 		stim::bil<T> n;
-		if(n.open(temp, R[0], R[1], R[2], offset, w)==false){        //open infile
+		if(n.open(temp, X(), Y(), Z(), offset, w)==false){        //open infile
 			std::cout<<"ERROR: unable to open input file"<<std::endl;
 			exit(1);
 		}
@@ -354,8 +343,8 @@ public:
 	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);
+		unsigned int L = X() * Z() * sizeof(T);		//calculate the number of bytes of a ZX slice
+		unsigned int jump = (Y() - 1) * X() * sizeof(T);
 		std::ofstream target(outname.c_str(), std::ios::binary);
 		std::string headername = outname + ".hdr";
@@ -363,12 +352,12 @@ public:
 		T * p;			//pointer to the current spectrum
 		p = (T*)malloc(L);
-		for ( unsigned i = 0; i < R[1]; i++)
+		for ( unsigned i = 0; i < Y(); i++)
 		{
-			file.seekg(R[0] * i * sizeof(T), std::ios::beg);
-			for ( unsigned j = 0; j < R[2]; j++ )
+			file.seekg(X() * i * sizeof(T), std::ios::beg);
+			for ( unsigned j = 0; j < Z(); j++ )
 			{
-				file.read((char *)(p + j * R[0]), sizeof(T) * R[0]);
+				file.read((char *)(p + j * X()), sizeof(T) * X());
 				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	
@@ -391,7 +380,7 @@ public:
 	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size.
 	bool baseline_band(double lb, double rb, T* lp, T* rp, double wavelength, T* result){
-		unsigned XY = R[0] * R[1];
+		unsigned XY = X() * Y();
 		band(result, wavelength);		//get band
 		//perform the baseline correction
@@ -412,7 +401,7 @@ public:
 		T* lp;
 		T* rp;
-		unsigned XY = R[0] * R[1];
+		unsigned XY = X() * Y();
 		unsigned S = XY * sizeof(T);
 		lp = (T*) malloc(S);			//memory allocation
 		rp = (T*) malloc(S);
@@ -442,7 +431,7 @@ public:
 		T* cur;		//current band 1
 		T* cur2;	//current band 2
-		unsigned XY = R[0] * R[1];
+		unsigned XY = X() * Y();
 		unsigned S = XY * sizeof(T);
 		lp = (T*) malloc(S);			//memory allocation
@@ -523,14 +512,14 @@ public:
 	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
 	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));	
+		T* p1 = (T*)malloc(X() * Y() * sizeof(T));	
+		T* p2 = (T*)malloc(X() * Y() * 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++){
+		for(unsigned i = 0; i < X() * Y(); i++){
 			if(p1[i] == 0 && p2[i] ==0)
 				result[i] = 1;
 			else
@@ -554,14 +543,14 @@ public:
 	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));	
+		T* p1 = (T*)malloc(X() * Y() * sizeof(T));	
+		T* p2 = (T*)malloc(X() * Y() * 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++){
+		for(unsigned i = 0; i < X() * Y(); i++){
 			if(p1[i] == 0 && p2[i] ==0)
 				result[i] = 1;
 			else
@@ -587,14 +576,14 @@ public:
 	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));	
+		T* p1 = (T*)malloc(X() * Y() * sizeof(T));	
+		T* p2 = (T*)malloc(X() * Y() * 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++){
+		for(unsigned i = 0; i < X() * Y(); i++){
 			if(p1[i] == 0 && p2[i] ==0)
 				result[i] = 1;
 			else
@@ -619,7 +608,7 @@ public:
 		T* cur;		//current band 1
 		T* cur2;	//current band 2
-		unsigned XY = R[0] * R[1];
+		unsigned XY = X() * Y();
 		unsigned S = XY * sizeof(T);
 		lp = (T*) malloc(S);			//memory allocation
@@ -695,14 +684,14 @@ public:
 	/// @param rab is the label for the end of the peak
 	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
 	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));	
+		T* p1 = (T*)malloc(X() * Y() * sizeof(T));	
+		T* p2 = (T*)malloc(X() * Y() * 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++){
+		for(unsigned i = 0; i < X() * Y(); i++){
 			if(p1[i] == 0 && p2[i] ==0)
 				result[i] = 1;
 			else
@@ -723,10 +712,10 @@ public:
 	/// @param p is a pointer to a pre-allocated array at least X * Y in size
 	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
+		T* temp = (T*)malloc(X() * Y() * sizeof(T));		//allocate memory for the certain band
 		band(temp, mask_band);
-		for (unsigned int i = 0; i < R[0] * R[1]; i++) {
+		for (unsigned int i = 0; i < X() * Y(); i++) {
 				if (temp[i] < threshold)
 					p[i] = 0;
 				else
@@ -746,12 +735,12 @@ public:
 		std::ofstream target(outfile.c_str(), std::ios::binary);
-		unsigned XY = R[0] * R[1];		//calculate number of a band
+		unsigned XY = X() * Y();		//calculate number of a band
 		unsigned L = XY * sizeof(T);
 		T * temp = (T*)malloc(L);
-		for (unsigned i = 0; i < R[2]; i++)			//for each spectral bin
+		for (unsigned i = 0; i < Z(); i++)			//for each spectral bin
 		{
 			band_index(temp, i);					//get the specified band (by index)
 			for ( unsigned j = 0; j < XY; j++)		// for each pixel
@@ -774,13 +763,13 @@ public:
 	bool sift_mask(std::string outfile, unsigned char* p){
 		std::ofstream target(outfile.c_str(), std::ios::binary);
 		// open a band (XY plane)
-		unsigned XY = R[0] * R[1]; //Number of XY pixels
+		unsigned XY = X() * Y(); //Number of XY pixels
 		unsigned L = XY * sizeof(T); //size of XY pixels
 		T * temp = (T*)malloc(L); //allocate memory for a band
 		T * temp_vox = (T*)malloc(sizeof(T)); //allocate memory for one voxel
-		for (unsigned i = 0; i < R[2]; i++)			//for each spectral bin
+		for (unsigned i = 0; i < Z(); i++)			//for each spectral bin
 		{
 			band_index(temp, i);					//get the specified band (XY sheet by index)
 			for (unsigned j = 0; j < XY; j++)		// for each pixel
@@ -804,18 +793,18 @@ public:
 	/// @param p is a pointer to memory of size X * Y * sizeof(T) that will store the band averages.
 	bool band_avg(T* p){
-		unsigned long long XY = R[0] * R[1];
+		unsigned long long XY = X() * Y();
 		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];
+			p[j] /= (T)Z();
 		}
 		//get every band and add them all
-		for (unsigned i = 1; i < R[2]; i++){
+		for (unsigned i = 1; i < Z(); i++){
 			band_index(temp, i);
 			for (unsigned j = 0; j < XY; j++){
-				p[j] += temp[j]/(T)R[2];
+				p[j] += temp[j]/(T)Z();
 			}
 		}
 		free(temp);
@@ -827,7 +816,7 @@ public:
 	/// @param p is a pointer to pre-allocated memory of size [B * sizeof(T)] that stores the mean spectrum
 	/// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
 	bool avg_band(T*p, unsigned char* mask){
-		unsigned long long XY = R[0] * R[1];
+		unsigned long long XY = X() * Y();
 		unsigned count = 0;						//count will store the number of masked pixels
 		T* temp = (T*)malloc(sizeof(T) * XY);
 		//calculate this loop counts the number of true pixels in the mask
@@ -836,9 +825,9 @@ public:
 				count++;
 			}
 		}
-		//this loops goes through each band in B (R[2])
+		//this loops goes through each band in B (Z())
 		//	masked (or valid) pixels from that band are averaged and the average is stored in p
-		for (unsigned i = 0; i < R[2]; i++){
+		for (unsigned i = 0; i < Z(); i++){
 			p[i] = 0;
 			band_index(temp, i);				//get the band image and store it in temp
 			for (unsigned j = 0; j < XY; j++){	//loop through temp, averaging valid pixels
@@ -858,8 +847,8 @@ public:
 	/// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
 	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];
+		unsigned long long xy = X() * Y();
+		unsigned int B = Z();
 		T* bandi = (T*)malloc(sizeof(T) * xy);
 		T* bandj = (T*)malloc(sizeof(T) * xy);
@@ -909,11 +898,11 @@ public:
 		//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
+		unsigned long long jumpb = X() * (Y() - lin) * sizeof(T);		//jump pointer to the next band
+		unsigned long long jumpl = (X() - 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++)
+		file.seekg((y0 * X() + x0) * sizeof(T), std::ios::beg);
+		for (unsigned i = 0; i < Z(); i++)
 		{
 			for (unsigned j = 0; j < lin; j++)
 			{
@@ -710,6 +710,41 @@ public:
 		return false;
 	}
+	bool spectrum(void* ptr, unsigned int x, unsigned int y){
+
+		if(header.interleave == envi_header::BSQ){		//if the infile is bsq file
+			if(header.data_type ==envi_header::float32)
+				return ((bsq<float>*)file)->spectrum((float*)ptr, x, y);
+			else if (header.data_type == envi_header::float64)
+				return ((bsq<double>*)file)->spectrum((double*)ptr, x, y);
+			else{
+				std::cout << "ERROR: unidentified data type" << std::endl;
+				exit(1);
+			}
+		}
+		else if (header.interleave == envi_header::BIL){
+			if (header.data_type == envi_header::float32)
+				return ((bil<float>*)file)->spectrum((float*)ptr, x, y);
+			else if (header.data_type == envi_header::float64)
+				return ((bil<double>*)file)->spectrum((double*)ptr, x, y);
+			else{
+				std::cout << "ERROR: unidentified data type" << std::endl;
+				exit(1);
+			}
+		}
+		else if (header.interleave == envi_header::BIP){
+			if (header.data_type == envi_header::float32)
+				return ((bip<float>*)file)->spectrum((float*)ptr, x, y);
+			else if (header.data_type == envi_header::float64)
+				return ((bip<double>*)file)->spectrum((double*)ptr, x, y);
+			else{
+				std::cout << "ERROR: unidentified data type" << std::endl;
+				exit(1);
+			}
+		}
+		return false;
+	}
+
 	/// Retrieve a single band (based on index) and stores it in pre-allocated memory.
 	/// @param p is a pointer to an allocated region of memory at least X * Y * sizeof(T) in size.
@@ -762,50 +797,6 @@ public:
 		return true;
 	}
-	//p:start positon; N: number of pixels saved in X;
-	bool feature_matrix(void * X,  unsigned char * mask, unsigned start, unsigned N)
-	{
-		//save pixels in X as floating numbers: float * p
-		float * p = (float*)malloc(header.bands * sizeof(float));		//save pixel information
-		unsigned pixels = header.samples * header.lines;		//calculate pixel numbers in a band
-		unsigned count = 0;		//for counting use
-		unsigned j = 0;		//memory the pointer location in X
-
-		//create two indices into the mask (mask index and pixel index)
-		unsigned mi = 0;	//valid pixel index
-		unsigned pi = 0;	//actual pixel in the mask
-
-		//find the actual pixel index for the mask index "start"
-		while(mi < start){
-			if(mask[pi])
-				mi++;
-
-			pi++;
-		}
-
-		for(unsigned i = pi; i < pixels; i++){
-			if(mask[i] != 0){
-				pixel(p, i);
-				//copy p to X
-				for(unsigned k = 0; k < header.bands; k++){
-					((float *)X)[j] = p[k];
-					j++;
-				}
-				count++;
-				if(count == N)
-					break;
-			}
-			else
-				continue;
-		}
-		if(count < N){
-			std::cout << "number of valid pixels in the mask : " << count <<"is less than N: "<< N;
-			exit(1);
-		}
-		free(p);
-		return true;
-	}
-
 	/// Calculate the mean value for all masked (or valid) pixels in a band and returns the average spectrum
 	/// @param p is a pointer to pre-allocated memory of size [B * sizeof(T)] that stores the mean spectrum