Merge branch 'master' of git.stim.ee.uh.edu:codebase/stimlib

David Mayerich
2 parents 721fd9d1 9b766f1f
Showing 12 changed files with 273 additions and 234 deletions Show diff stats
stim/cuda/branch_detection.cuh
stim/cuda/spider_cost.cuh
stim/envi/bil.h
stim/envi/bip.h
stim/envi/bsq.h
stim/envi/hsi.h
stim/grids/image_stack.h
stim/image/image.h
stim/math/circle.h
stim/parser/filename.h
stim/visualization/cylinder.h
stim/visualization/obj.h
@@ -7,7 +7,7 @@
 #include <stim/cuda/cuda_texture.cuh>
 #include <stim/cuda/filter.cuh>
 typedef unsigned int uint;
-typedef unsigned int uchar;
+//typedef unsigned int uchar;
 std::vector< stim::vec<float> >
@@ -3,7 +3,7 @@
 #include <assert.h>
 #include <cuda.h>
-#include <cuda_runtime.h>
+//#include <cuda_runtime.h>
 #include <stdio.h>
 #include <stim/visualization/colormap.h>
 #include <sstream>
@@ -3,8 +3,10 @@
 #include "../envi/envi_header.h"
 #include "../envi/hsi.h"
+#include "../math/fd_coefficients.h"
 #include <cstring>
 #include <utility>
+#include <deque>
 namespace stim{
@@ -22,26 +24,14 @@ template &lt;typename T&gt;
 class bil: public hsi<T> {
 protected:
-	
-	//std::vector<double> w;		//band wavelengths
-	/*unsigned long long X(){
-		return R[0];
-	}
-	unsigned long long Y(){
-		return R[2];
-	}
-	unsigned long long Z(){
-		return R[1];
-	}*/
+
 	using hsi<T>::w;				//use the wavelength array in stim::hsi
 	using hsi<T>::nnz;
 	using binary<T>::progress;
-
-	/// Call the binary nnz() function for the BIL orientation
-	//unsigned long long nnz(unsigned char* mask){
-	//	return binary<T>::nnz(mask, X()*Y());
-	//}
+	using hsi<T>::X;
+	using hsi<T>::Y;
+	using hsi<T>::Z;
 public:
@@ -59,11 +49,11 @@ public:
 	/// @param B is the number of samples (bands) along dimension 3
 	/// @param header_offset is the number of bytes (if any) in the binary header
 	/// @param wavelengths is an optional STL vector of size B specifying a numerical label for each band
-	bool open(std::string filename, 
-			  unsigned long long X, 
-			  unsigned long long Y, 
-			  unsigned long long B, 
-			  unsigned long long header_offset, 
+	bool open(std::string filename,
+			  unsigned long long X,
+			  unsigned long long Y,
+			  unsigned long long B,
+			  unsigned long long header_offset,
 			  std::vector<double> wavelengths){
 		w = wavelengths;
@@ -113,7 +103,7 @@ public:
 		unsigned long long S = XY * sizeof(T);		//calculate the number of bytes of a band
 		unsigned long long page=0;                      //bands around the wavelength
-		
+
 		//get the bands numbers around the wavelength
@@ -173,7 +163,7 @@ public:
 		//if wavelength is smaller than the first one in header file
 		if ( w[page] > wavelength ){
-			memcpy(p, c, L);	
+			memcpy(p, c, L);
 			return true;
 		}
@@ -193,16 +183,16 @@ public:
 			memcpy(p1, c + (page - 1) * X(), L);
 			memcpy(p2, c + page * X(), L);
-			
+
 			for(unsigned long long i=0; i < X(); i++){
 				double r = (double) (wavelength - w[page-1]) / (double) (w[page] - w[page-1]);
 				p[i] = (T)(((double)p2[i] - (double)p1[i]) * r + (double)p1[i]);
 			}
 		}
-		else                           //if the wavelength is equal to a wavelength in header file		
+		else                           //if the wavelength is equal to a wavelength in header file
 			memcpy(p, c + page * X(), L);
-		
-		return true;		
+
+		return true;
 	}
 	/// Retrieve a single spectrum (B-axis line) at a given (x, y) location and stores it in pre-allocated memory.
@@ -227,24 +217,24 @@ public:
 		//get the pixel
 		return spectrum(p, x, y);
 	}
-	
+
 	//given a Y ,return a XZ slice
 	bool read_plane_y(T * p, unsigned long long 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.
 	/// @param outname is the name of the output file used to store the resulting baseline-corrected data.
 	/// @param wls is the list of baseline points based on band labels.
 	bool baseline(std::string outname, std::vector<double> wls, unsigned char* mask = NULL, bool PROGRESS = false){
-	
+
 		unsigned long long 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		
-		
+		std::string headername = outname + ".hdr";              //the header file name
+
 		//simplify image resolution
 		unsigned long long ZX = Z() * X();		//calculate the number of points in a Y slice
 		unsigned long long L = ZX * sizeof(T);			//calculate the number of bytes of a Y slice
@@ -252,10 +242,10 @@ public:
 		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));
@@ -269,19 +259,19 @@ public:
 			exit(1);
 		}
 	//	loop start	correct every y slice
-		
+
 		for (unsigned long long k =0; k < Y(); k++)
 		{
 			//get the current y slice
 			read_plane_y(c, k);
-		
-			//initialize lownum, highnum, low, high		
+
+			//initialize lownum, highnum, low, high
 			ai = w[0];
 			control=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];			
+				bi = wls[control];
 				memset(a, (char)0, X() * sizeof(T) );
 			}
 			//else get the low band
@@ -292,39 +282,39 @@ public:
 			}
 			//get the high band
 			read_x_from_xz(b, c, bi);
-		
+
 			//correct every YZ line
-			
+
 			for(unsigned long long 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];
 						read_x_from_xz(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];
-				
+
 				unsigned long long jump = cii * X();
 				//perform the baseline correction
 				for(unsigned i=0; i < X(); i++)
@@ -332,22 +322,22 @@ public:
 					double r = (double) (ci - ai) / (double) (bi - ai);
 					c[i + jump] =(T) ( c[i + jump] - (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 YZ line end
+			target.write(reinterpret_cast<const char*>(c), L);   //write the corrected data into destination
 			if(PROGRESS) progress = (double)(k+1) / Y() * 100;
-		}//loop for Y slice end		
-		
+		}//loop for Y slice end
+
 		free(a);
 		free(b);
 		free(c);
 		target.close();
-		return true;	
-		
+		return true;
+
 	}
-		
+
 	/// Normalize all spectra based on the value of a single band, storing the result in a new BSQ file.
 	/// @param outname is the name of the output file used to store the resulting baseline-corrected data.
@@ -358,7 +348,7 @@ public:
 		unsigned long long B = Z();		//calculate the number of bands
 		unsigned long long ZX = Z() * X();
 		unsigned long long XY = X() * Y();	//calculate the number of pixels in a band
-		unsigned long long S = XY * sizeof(T);		//calculate the number of bytes in a band		
+		unsigned long long S = XY * sizeof(T);		//calculate the number of bytes in a band
 		unsigned long long L = ZX * sizeof(T);
 		std::ofstream target(outname.c_str(), std::ios::binary);	//open the target binary file
@@ -368,7 +358,7 @@ public:
 		T * b;				//pointer to the standard band
 		b = (T*)malloc( S );     //memory allocation
-		c = (T*)malloc( L ); 
+		c = (T*)malloc( L );
 		band(b, w);             //get the certain band into memory
@@ -383,13 +373,13 @@ public:
 						c[m + i * X()] = (T)0.0;
 					else
 						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
 			if(PROGRESS) progress = (double)(j+1) / Y() * 100;
 		}
-		
+
 		free(b);
 		free(c);
 		target.close();
@@ -397,26 +387,26 @@ public:
 	}
 	void normalize(std::string outfile, unsigned char* mask = NULL, bool PROGRESS = false){}
-	
+
 	/// Convert the current BIL file to a BSQ file with the specified file name.
 	/// @param outname is the name of the output BSQ file to be saved to disk.
 	bool bsq(std::string outname, bool PROGRESS = false)
 	{
 		unsigned long long 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";
-		
+
 		T * p;			//pointer to the current band
 		p = (T*)malloc(S);
-		
+
 		for ( unsigned long long i = 0; i < Z(); i++)
-		{			
+		{
 				band_index(p, i);
 				target.write(reinterpret_cast<const char*>(p), S);   //write a band data into target file
-				if(PROGRESS) progress = (double)(i+1) / Z() * 100;	//store the progress for the current operation	
+				if(PROGRESS) progress = (double)(i+1) / Z() * 100;	//store the progress for the current operation
 		}
 		free(p);
@@ -430,17 +420,17 @@ public:
 	bool bip(std::string outname, bool PROGRESS = false)
 	{
 		unsigned long long 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";
-		
+
 		T * p;			//pointer to the current XZ slice for bil file
 		p = (T*)malloc(S);
 		T * q;			//pointer to the current ZX slice for bip file
 		q = (T*)malloc(S);
-		
+
 		for ( unsigned long long i = 0; i < Y(); i++)
-		{			
+		{
 			read_plane_y(p, i);
 			for ( unsigned long long k = 0; k < Z(); k++)
 			{
@@ -448,10 +438,10 @@ public:
 				for ( unsigned long long j = 0; j < X(); j++)
 					q[k + j * Z()] = p[ks + j];
-				if(PROGRESS) progress = (double)((i+1) * Z() + k+1) / (Z() * Y()) * 100;	//store the progress for the current operation	
+				if(PROGRESS) progress = (double)((i+1) * Z() + k+1) / (Z() * Y()) * 100;	//store the progress for the current operation
 			}
-			
-			target.write(reinterpret_cast<const char*>(q), S);   //write a band data into target file	
+
+			target.write(reinterpret_cast<const char*>(q), S);   //write a band data into target file
 		}
 		free(p);
@@ -498,7 +488,7 @@ public:
 		rp = (T*) malloc(S);
 		band(lp, lb);
-		band(rp, rb);		
+		band(rp, rb);
 		baseline_band(lb, rb, lp, rp, bandwavelength, result);
@@ -570,7 +560,7 @@ public:
 		}
 		baseline_band(lb, rb, lp, rp, lab, cur2);		//beginnning part
 		baseline_band(lb, rb, lp, rp, w[ai], cur);
-		for(unsigned long long j = 0; j < XY; j++){	
+		for(unsigned long long j = 0; j < XY; j++){
 			result[j] += (T)((w[ai] - lab) * ((double)cur[j] + (double)cur2[j]) / 2.0);
 		}
@@ -581,7 +571,7 @@ public:
 			baseline_band(lb, rb, lp, rp, w[ai], cur2);
 			for(unsigned long long j = 0; j < XY; j++)
 			{
-				result[j] += (T)((w[ai] - w[ai-1]) * ((double)cur[j] + (double)cur2[j]) / 2.0); 
+				result[j] += (T)((w[ai] - w[ai-1]) * ((double)cur[j] + (double)cur2[j]) / 2.0);
 			}
 			std::swap(cur,cur2);		//swap the band pointers
 		}
@@ -606,7 +596,7 @@ public:
 		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);
@@ -620,9 +610,9 @@ public:
 		free(p1);
 		free(p2);
-		return true;	
+		return true;
 	}
-	
+
 	/// Compute the ratio between a peak area and peak height.
 	/// @param lb1 is the label value for the left baseline point for the first peak (numerator)
@@ -634,10 +624,10 @@ public:
 	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
 	bool pa_to_ph(T* result, double lb1, double rb1, double lab1, double rab1,
 					double lb2, double rb2, double pos, unsigned char* mask = NULL){
-		
+
 		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);
@@ -651,9 +641,9 @@ public:
 		free(p1);
 		free(p2);
-		return true;	
-	}		
-	
+		return true;
+	}
+
 	/// Compute the ratio between two peak areas.
 	/// @param lb1 is the label value for the left baseline point for the first peak (numerator)
@@ -663,14 +653,14 @@ public:
 	/// @param lb2 is the label value for the left baseline point for the second peak (denominator)
 	/// @param rb2 is the label value for the right baseline point for the second peak (denominator)
 	/// @param lab2 is the label value for the left bound (start of the integration) of the second peak (denominator)
-	/// @param rab2 is the label value for the right bound (end of the integration) of the second peak (denominator)	
+	/// @param rab2 is the label value for the right bound (end of the integration) of the second peak (denominator)
 	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
 	bool pa_to_pa(T* result, double lb1, double rb1, double lab1, double rab1,
 					double lb2, double rb2, double lab2, double rab2, unsigned char* mask = NULL){
-		
+
 		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);
@@ -684,8 +674,8 @@ public:
 		free(p1);
 		free(p2);
-		return true;	
-	}		
+		return true;
+	}
 	/// Compute the definite integral of a baseline corrected peak.
@@ -745,7 +735,7 @@ public:
 		}
 		baseline_band(lb, rb, lp, rp, lab, cur2);		//beginnning part
 		baseline_band(lb, rb, lp, rp, w[ai], cur);
-		for(unsigned long long j = 0; j < XY; j++){	
+		for(unsigned long long j = 0; j < XY; j++){
 			result[j] += (T)((w[ai] - lab) * (w[ai] + lab) * ((double)cur[j] + (double)cur2[j]) / 4.0);
 		}
@@ -756,7 +746,7 @@ public:
 			baseline_band(lb, rb, lp, rp, w[ai], cur2);
 			for(unsigned long long j = 0; j < XY; j++)
 			{
-				result[j] += (T)((w[ai] - w[ai-1]) * (w[ai] + w[ai-1]) * ((double)cur[j] + (double)cur2[j]) / 4.0); 
+				result[j] += (T)((w[ai] - w[ai-1]) * (w[ai] + w[ai-1]) * ((double)cur[j] + (double)cur2[j]) / 4.0);
 			}
 			std::swap(cur,cur2);		//swap the band pointers
 		}
@@ -778,7 +768,7 @@ public:
 	bool centroid(T* result, double lb, double rb, double lab, double rab, unsigned char* mask = NULL){
 		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);
@@ -790,7 +780,7 @@ public:
 		free(p1);
 		free(p2);
-		return true;	
+		return true;
 	}
 	/// Create a mask based on a given band and threshold value.
@@ -843,7 +833,7 @@ public:
 					continue;
 				}
 			}
-			target.write(reinterpret_cast<const char*>(temp), L);   //write a band data into target file	
+			target.write(reinterpret_cast<const char*>(temp), L);   //write a band data into target file
 			if(PROGRESS) progress = (double)(i+1) / (double)Y() * 100;
 		}
 		target.close();
@@ -861,7 +851,7 @@ public:
 		T* line = (T*) malloc( Lbytes );					//allocate space for a line
 		file.seekg(0, std::ios::beg);	//seek to the beginning of the file
-		
+
 		size_t pl;
 		size_t p = 0;										//create counter variables
 		for(size_t y = 0; y < Y(); y++){					//for each line in the data set
@@ -903,7 +893,7 @@ public:
 		{
 			read_plane_y(slice, y);							//retrieve an ZX page, store in "slice"
-			//for each sample along X			
+			//for each sample along X
 			for (unsigned long long x = 0; x < X(); x++)		//Select a pixel by choosing X coordinate in the page, X()
 			{
 				//if the mask != 0 at that xy pixel
@@ -1082,7 +1072,7 @@ public:
 			for (unsigned long long z = b0; z < b1; z++)
 			{
 				file.read((char *)(temp + z * samples), sizeof(T) * samples);
-				file.seekg(jumpb, std::ios::cur);    //go to the next band	
+				file.seekg(jumpb, std::ios::cur);    //go to the next band
 				if(PROGRESS) progress = (double)(x * Z() + z+1) / (lines * Z()) * 100;
 			}
@@ -1193,7 +1183,7 @@ public:
 		size_t N = C.size();										//get the number of bands that the kernel spans
 		std::deque<T*> frame(N, NULL);								//create an array to store pointers to each frame
 		for(size_t f = 0; f < N; f++)
-			frame[f] = (T*)malloc(Xb);								//allocate space for the frame		
+			frame[f] = (T*)malloc(Xb);								//allocate space for the frame
 		T* outline = (T*)malloc(Xb);								//allocate space for the output band
@@ -1219,23 +1209,23 @@ public:
 				frame.pop_front();										//pop the first element
 			}
 			if(PROGRESS) progress = (double)(y+1) / (double)Y() * 100;
-		}		
+		}
 	}
 	/// Approximate the spectral derivative of the image
 	void deriv(std::string outfile, size_t d, size_t order, const std::vector<double> w = std::vector<double>(), unsigned char* mask = NULL, bool PROGRESS = false){
 		std::ofstream out(outfile.c_str(), std::ios::binary);		//open the output file for writing
-				
+
 		size_t Xb = X() * sizeof(T);								//calculate the size of a line (frame) in bytes
 		size_t XBb = Xb * Z();
 		size_t B = Z();
-		file.seekg(0, std::ios::beg);								//seek to the beginning of the file		
+		file.seekg(0, std::ios::beg);								//seek to the beginning of the file
 		size_t N = order + d;										//approximating a derivative requires order + d samples
 		std::deque<T*> frame(N, NULL);								//create an array to store pointers to each frame
 		for(size_t f = 0; f < N; f++)								//for each frame
-			frame[f] = (T*)malloc(Xb);								//allocate space for the frame	
+			frame[f] = (T*)malloc(Xb);								//allocate space for the frame
 		T* outline = (T*)malloc(Xb);								//allocate space for the output band
@@ -1273,7 +1263,7 @@ public:
 					}
 				}
 				out.write((char*)outline, Xb);											//output the band
-				
+
 			}
 			if(PROGRESS) progress = (double)(y+1) / (double)Y() * 100;
 		}
@@ -28,15 +28,18 @@ template &lt;typename T&gt;
 class bip: public hsi<T> {
 protected:
-	
-	
+
+
 	//std::vector<double> w; //band wavelength
 	unsigned long long offset;		//header offset
 	using hsi<T>::w;				//use the wavelength array in stim::hsi
 	using hsi<T>::nnz;
 	using binary<T>::progress;
-	
+	using hsi<T>::X;
+	using hsi<T>::Y;
+	using hsi<T>::Z;
+
 public:
 	using binary<T>::open;
@@ -54,11 +57,11 @@ public:
 	/// @param B is the number of samples (bands) along dimension 3
 	/// @param header_offset is the number of bytes (if any) in the binary header
 	/// @param wavelengths is an optional STL vector of size B specifying a numerical label for each band
-	bool open(std::string filename, 
-			  unsigned long long X, 
-			  unsigned long long Y, 
-			  unsigned long long B, 
-			  unsigned long long header_offset, 
+	bool open(std::string filename,
+			  unsigned long long X,
+			  unsigned long long Y,
+			  unsigned long long B,
+			  unsigned long long header_offset,
 			  std::vector<double> wavelengths){
 		//copy the wavelengths to the BSQ file structure
@@ -67,7 +70,7 @@ public:
 		offset = header_offset;
 		return open(filename, vec<unsigned long long>(B, X, Y), header_offset);
-		
+
 	}
 	/// Retrieve a single band (based on index) and stores it in pre-allocated memory.
@@ -160,7 +163,7 @@ public:
 			for(unsigned long long j = 0; j < X(); j++)
 			{
 				p[j] = c[j * B];
-			}		
+			}
 			return true;
 		}
@@ -192,7 +195,7 @@ public:
 			{
 				p2[j] = c[j * B + page];
 			}
-			
+
 			for(unsigned long long 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];
@@ -209,7 +212,7 @@ public:
 			}
 		}
-		return true;		
+		return true;
 	}
 	/// Retrieve a single pixel and store it in a pre-allocated double array.
@@ -246,20 +249,20 @@ public:
 		file.read((char *)p, sizeof(T) * Z());
 		return true;
 	}
-	
+
 	//given a Y ,return a ZX slice
 	bool read_plane_y(T * p, unsigned long long 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.
 	/// @param outname is the name of the output file used to store the resulting baseline-corrected data.
-	/// @param wls is the list of baseline points based on band labels.	
+	/// @param wls is the list of baseline points based on band labels.
 	bool baseline(std::string outname, std::vector<double> base_pts, unsigned char* mask = NULL, bool PROGRESS = false){
-		
+
 		std::ofstream target(outname.c_str(), std::ios::binary);	//open the target binary file
-		
+
 		unsigned long long N = X() * Y();						//calculate the total number of pixels to be processed
 		unsigned long long B = Z();								//get the number of bands
 		T* s = (T*)malloc(sizeof(T) * B);						//allocate memory to store a pixel
@@ -274,9 +277,9 @@ public:
 				memset(sbc, 0, sizeof(T) * B);					//set the baseline corrected spectrum to zero
 			}
 			else{
-				
+
 				pixel(s, n);										//retrieve the spectrum s
-				base_vals = interp_spectrum(s, base_pts);			//get the values at each baseline point
+				base_vals = hsi<T>::interp_spectrum(s, base_pts);			//get the values at each baseline point
 				ai = bi = 0;
 				aw = w[0];											//initialize the current baseline points (assume the spectrum starts at 0)
@@ -298,23 +301,23 @@ public:
 						aw = base_pts[ai];							//set the wavelength for the lower bound baseline point
 						av = base_vals[ai];							//set the value for the lower bound baseline point
 					}
-					sbc[b] = s[b] - lerp(w[b], av, aw, bv, bw);		//perform the baseline correction and save the new value
+					sbc[b] = s[b] - hsi<T>::lerp(w[b], av, aw, bv, bw);		//perform the baseline correction and save the new value
 				}
 			}
-		
+
 			if(PROGRESS) progress = (double)(n+1) / N * 100;	//set the current progress
 			target.write((char*)sbc, sizeof(T) * B);	//write the corrected data into destination
 		}														//end for each pixel
-		
+
 		free(s);
 		free(sbc);
 		target.close();
-		
-		return true;	
-		
+
+		return true;
+
 	}
-		
+
 	/// Normalize all spectra based on the value of a single band, storing the result in a new BSQ file.
 	/// @param outname is the name of the output file used to store the resulting baseline-corrected data.
@@ -323,8 +326,8 @@ public:
 	bool ratio(std::string outname, double w, unsigned char* mask = NULL, bool PROGRESS = false)
 	{
 		std::ofstream target(outname.c_str(), std::ios::binary);	//open the target binary file
-		std::string headername = outname + ".hdr";              //the header file name		
-		
+		std::string headername = outname + ".hdr";              //the header file name
+
 		unsigned long long N = X() * Y();						//calculate the total number of pixels to be processed
 		unsigned long long B = Z();								//get the number of bands
 		T* s = (T*)malloc(sizeof(T) * B);						//allocate memory to store a pixel
@@ -334,24 +337,24 @@ public:
 				memset(s, 0, sizeof(T) * B);					//set the output to zero
 			else{
 				pixel(s, n);										//retrieve the spectrum s
-				nv = interp_spectrum(s, w);							//find the value of the normalization band
-			
+				nv = hsi<T>::interp_spectrum(s, w);							//find the value of the normalization band
+
 				for(unsigned long long b = 0; b < B; b++)			//for each band in the spectrum
 					s[b] /= nv;										//divide by the normalization value
 			}
-		
+
 			if(PROGRESS) progress = (double)(n+1) / N * 100;	//set the current progress
 			target.write((char*)s, sizeof(T) * B);	//write the corrected data into destination
 		}														//end for each pixel
-		
+
 		free(s);
 		target.close();
 		return true;
 	}
 	void normalize(std::string outfile, unsigned char* mask = NULL, bool PROGRESS = false){}
-	
+
 	/// Convert the current BIP file to a BIL file with the specified file name.
@@ -359,17 +362,17 @@ public:
 	bool bil(std::string outname, bool PROGRESS = false)
 	{
 		unsigned long long 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";
-		
+
 		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 long long i = 0; i < Y(); i++)
-		{			
+		{
 			read_plane_y(p, i);
 			for ( unsigned long long k = 0; k < Z(); k++)
 			{
@@ -379,7 +382,7 @@ public:
 				if(PROGRESS) progress = (double)(i * Z() + k+1) / (Y() * Z()) * 100;
 			}
-			target.write(reinterpret_cast<const char*>(q), S);   //write a band data into target file	
+			target.write(reinterpret_cast<const char*>(q), S);   //write a band data into target file
 		}
 		free(p);
@@ -408,7 +411,7 @@ public:
 		}
 		return true;
 	}
-	
+
 	/// Return a baseline corrected band given two adjacent baseline points. The result is stored in a pre-allocated array.
 	/// @param lb is the label value for the left baseline point
@@ -425,7 +428,7 @@ public:
 		rp = (T*) malloc(S);
 		band(lp, lb);
-		band(rp, rb);		
+		band(rp, rb);
 		baseline_band(lb, rb, lp, rp, bandwavelength, result);
@@ -496,7 +499,7 @@ public:
 		}
 		baseline_band(lb, rb, lp, rp, lab, cur2);		//beginnning part
 		baseline_band(lb, rb, lp, rp, w[ai], cur);
-		for(unsigned long long j = 0; j < XY; j++){	
+		for(unsigned long long j = 0; j < XY; j++){
 			result[j] += (T)((w[ai] - lab) * ((double)cur[j] + (double)cur2[j]) / 2.0);
 		}
@@ -507,7 +510,7 @@ public:
 			baseline_band(lb, rb, lp, rp, w[ai], cur2);
 			for(unsigned long long j = 0; j < XY; j++)
 			{
-				result[j] += (T)((w[ai] - w[ai-1]) * ((double)cur[j] + (double)cur2[j]) / 2.0); 
+				result[j] += (T)((w[ai] - w[ai-1]) * ((double)cur[j] + (double)cur2[j]) / 2.0);
 			}
 			std::swap(cur,cur2);		//swap the band pointers
 		}
@@ -530,9 +533,9 @@ public:
 	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
 	bool ph_to_ph(T* result, double lb1, double rb1, double pos1, double lb2, double rb2, double pos2, unsigned char* mask = NULL){
-		T* p1 = (T*)malloc(X() * Y() * sizeof(T));	
-		T* p2 = (T*)malloc(X() * Y() * 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);
@@ -546,9 +549,9 @@ public:
 		free(p1);
 		free(p2);
-		return true;	
+		return true;
 	}
-	
+
 	/// Compute the ratio between a peak area and peak height.
 	/// @param lb1 is the label value for the left baseline point for the first peak (numerator)
@@ -560,10 +563,10 @@ public:
 	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
 	bool pa_to_ph(T* result, double lb1, double rb1, double lab1, double rab1,
 					double lb2, double rb2, double pos, unsigned char* mask = NULL){
-		
-		T* p1 = (T*)malloc(X() * Y() * sizeof(T));	
-		T* p2 = (T*)malloc(X() * Y() * 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);
@@ -577,9 +580,9 @@ public:
 		free(p1);
 		free(p2);
-		return true;	
-	}		
-	
+		return true;
+	}
+
 	/// Compute the ratio between two peak areas.
 	/// @param lb1 is the label value for the left baseline point for the first peak (numerator)
@@ -589,14 +592,14 @@ public:
 	/// @param lb2 is the label value for the left baseline point for the second peak (denominator)
 	/// @param rb2 is the label value for the right baseline point for the second peak (denominator)
 	/// @param lab2 is the label value for the left bound (start of the integration) of the second peak (denominator)
-	/// @param rab2 is the label value for the right bound (end of the integration) of the second peak (denominator)	
+	/// @param rab2 is the label value for the right bound (end of the integration) of the second peak (denominator)
 	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
 	bool pa_to_pa(T* result, double lb1, double rb1, double lab1, double rab1,
 					double lb2, double rb2, double lab2, double rab2, unsigned char* mask = NULL){
-		
-		T* p1 = (T*)malloc(X() * Y() * sizeof(T));	
-		T* p2 = (T*)malloc(X() * Y() * 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);
@@ -610,8 +613,8 @@ public:
 		free(p1);
 		free(p2);
-		return true;	
-	}		
+		return true;
+	}
 	/// Compute the definite integral of a baseline corrected peak.
@@ -671,7 +674,7 @@ public:
 		}
 		baseline_band(lb, rb, lp, rp, lab, cur2);		//beginnning part
 		baseline_band(lb, rb, lp, rp, w[ai], cur);
-		for(unsigned long long j = 0; j < XY; j++){	
+		for(unsigned long long j = 0; j < XY; j++){
 			result[j] += (T)((w[ai] - lab) * (w[ai] + lab) * ((double)cur[j] + (double)cur2[j]) / 4.0);
 		}
@@ -682,7 +685,7 @@ public:
 			baseline_band(lb, rb, lp, rp, w[ai], cur2);
 			for(unsigned long long j = 0; j < XY; j++)
 			{
-				result[j] += (T)((w[ai] - w[ai-1]) * (w[ai] + w[ai-1]) * ((double)cur[j] + (double)cur2[j]) / 4.0); 
+				result[j] += (T)((w[ai] - w[ai-1]) * (w[ai] + w[ai-1]) * ((double)cur[j] + (double)cur2[j]) / 4.0);
 			}
 			std::swap(cur,cur2);		//swap the band pointers
 		}
@@ -702,9 +705,9 @@ 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 centroid(T* result, double lb, double rb, double lab, double rab, unsigned char* mask = NULL){
-		T* p1 = (T*)malloc(X() * Y() * sizeof(T));	
-		T* p2 = (T*)malloc(X() * Y() * 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);
@@ -716,7 +719,7 @@ public:
 		free(p1);
 		free(p2);
-		return true;	
+		return true;
 	}
 	/// Create a mask based on a given band and threshold value.
@@ -788,7 +791,7 @@ public:
 		T* band = (T*) malloc( Bbytes );					//allocate space for a line
 		file.seekg(0, std::ios::beg);	//seek to the beginning of the file
-		
+
 		size_t p = 0;										//create counter variables
 		for(size_t xy = 0; xy < XY; xy++){					//for each pixel
 			if(mask[xy]){									//if the current pixel is masked
@@ -800,7 +803,7 @@ public:
 				p++;									//increment the pixel pointer
 			}
 			else
-				file.seekg(Bbytes, std::ios::cur);			//otherwise skip this band			
+				file.seekg(Bbytes, std::ios::cur);			//otherwise skip this band
 		}
 		return true;
 	}
@@ -809,7 +812,7 @@ public:
 	bool sift(std::string outfile, unsigned char* mask, bool PROGRESS = false){
 		//reset the file pointer to the beginning of the file
-		file.seekg(0, std::ios::beg);		
+		file.seekg(0, std::ios::beg);
 		// open an output stream
 		std::ofstream target(outfile.c_str(), std::ios::binary);
@@ -862,7 +865,7 @@ public:
 		std::ofstream target(outfile.c_str(), std::ios::binary);
 		//reset the file pointer to the beginning of the file
-		file.seekg(0, std::ios::beg);		
+		file.seekg(0, std::ios::beg);
 		//allocate space for a single spectrum
 		unsigned long long B = Z();
@@ -942,10 +945,10 @@ public:
 	/// Note that cuBLAS only supports integer-sized arrays, so there may be issues with large spectra
 	bool co_matrix_cublas(double* co, double* avg, unsigned char *mask, bool PROGRESS = false){
-		cudaError_t cudaStat;    
+		cudaError_t cudaStat;
 		cublasStatus_t stat;
 		cublasHandle_t handle;
-		
+
 		progress = 0;													//initialize the progress to zero (0)
 		unsigned long long XY = X() * Y();									//calculate the number of elements in a band image
 		unsigned long long B = Z();											//calculate the number of spectral elements
@@ -959,7 +962,7 @@ public:
 		cudaStat = cudaMemset(A_dev, 0, B * B * sizeof(double));			//initialize the covariance matrix to zero (0)
 		cudaStat = cudaMalloc(&avg_dev, B * sizeof(double));				//allocate space on the CUDA device for the average spectrum
 		stat = cublasSetVector((int)B, sizeof(double), avg, 1, avg_dev, 1);		//copy the average spectrum to the CUDA device
-		
+
 		double ger_alpha = 1.0/(double)XY;									//scale the outer product by the inverse of the number of samples (mean outer product)
 		double axpy_alpha = -1;												//multiplication factor for the average spectrum (in order to perform a subtraction)
@@ -994,7 +997,7 @@ public:
 		return true;
 	}
 #endif
-	
+
 	/// Calculate the covariance matrix for all masked pixels in the image with 64-bit floating point precision.
 	/// @param co is a pointer to pre-allocated memory of size [B * B] that stores the resulting covariance matrix
@@ -1074,7 +1077,7 @@ public:
 					}
 				}
 			}
-			
+
 			target.write(reinterpret_cast<const char*>(rs), sizeof(T) * M);					//write the projected vector
 			if(PROGRESS) progress = (double)(xy+1) / XY * 100;
 		}
@@ -1104,7 +1107,7 @@ public:
 			memset(s, 0, sizeof(T) * B);							//initialize the spectrum to zero (0)
 			for(unsigned long long c = 0; c < C; c++){				//for each basis vector coefficient
 				bv = &basis[c * B];									//assign 'bv' to the beginning of the basis vector
-				for(unsigned long long b = 0; b < B; b++){			//for each component of the basis vector					
+				for(unsigned long long b = 0; b < B; b++){			//for each component of the basis vector
 					s[b] += (T)((double)coeff[c] * bv[b] + center[b]);			//calculate the contribution of each element of the basis vector in the final spectrum
 				}
 			}
@@ -1172,7 +1175,7 @@ public:
 				//read the cropped spectral region
 				file.read( (char*) temp, L );
 				//pixel(temp, sp + x + y * X());
-				out.write(reinterpret_cast<const char*>(temp), L);   //write slice data into target file	
+				out.write(reinterpret_cast<const char*>(temp), L);   //write slice data into target file
 				file.seekg(jump_sample, std::ios::cur);
@@ -1253,7 +1256,7 @@ public:
 					memset(spec, 0, spec_bytes);
 				out.write((char*)spec, spec_bytes);					//write to the output file
 			}
-			if(PROGRESS) progress = (double)( (y+1) * S[0] + 1) / (double) (S[0] * S[1]) * 100; 
+			if(PROGRESS) progress = (double)( (y+1) * S[0] + 1) / (double) (S[0] * S[1]) * 100;
 		}
 	}
@@ -1267,7 +1270,7 @@ public:
 	void convolve(std::string outfile, std::vector<double> C, size_t start, size_t end, unsigned char* mask = NULL, bool PROGRESS = false){
 		std::ofstream out(outfile.c_str(), std::ios::binary);		//open the output file for writing
-		
+
 		size_t N = end - start + 1;									//number of bands in the output spectrum
 		size_t Nb = N * sizeof(T);									//size of the output spectrum in bytes
 		size_t B = Z();												//calculate the number of values in a spectrum
@@ -1297,8 +1300,8 @@ public:
 	void deriv(std::string outfile, size_t d, size_t order, const std::vector<double> w, unsigned char* mask = NULL, bool PROGRESS = false){
 		std::ofstream out(outfile.c_str(), std::ios::binary);		//open the output file for writing
-		
-		
+
+
 		size_t B = Z();												//calculate the number of values in a spectrum
 		size_t Bb = B * sizeof(T);									//calculate the size of a spectrum in bytes
@@ -37,6 +37,9 @@ protected:
 	using hsi<T>::w;				//use the wavelength array in stim::hsi
 	using hsi<T>::nnz;
 	using binary<T>::progress;
+	using hsi<T>::X;
+	using hsi<T>::Y;
+	using hsi<T>::Z;
 public:
@@ -55,11 +58,11 @@ public:
 	/// @param B is the number of samples (bands) along dimension 3
 	/// @param header_offset is the number of bytes (if any) in the binary header
 	/// @param wavelengths is an STL vector of size B specifying a numerical label for each band
-	bool open(std::string filename, 
-			  unsigned long long X, 
-			  unsigned long long Y, 
-			  unsigned long long B, 
-			  unsigned long long header_offset, 
+	bool open(std::string filename,
+			  unsigned long long X,
+			  unsigned long long Y,
+			  unsigned long long B,
+			  unsigned long long header_offset,
 			  std::vector<double> wavelengths){
 		//copy the wavelengths to the BSQ file structure
@@ -366,7 +369,7 @@ public:
 			out.write((char*)band, XYb);
 			if(PROGRESS) progress = (double) (b+1) / (double)B * 50 + 50;
 		}
-		
+
 	}
 	/// Convert the current BSQ file to a BIL file with the specified file name.
@@ -489,7 +492,7 @@ public:
 		}
 		//find the indices of the left and right baseline points
-		while (lab >= w[ai]){ 
+		while (lab >= w[ai]){
 			ai++;
 		}
 		while (rab <= w[bi]){
@@ -822,13 +825,13 @@ public:
 					i++;
 					//std::cout<<i<<std::endl;
 				}
-				
+
 			}
 		}
 		return true;
 	}
-	
+
 	/// Saves to disk only those spectra corresponding to mask values != 0
 	/// @param outfile is the name of the sifted ENVI file to be written to disk
 	/// @param unsigned char* p is the mask file used for sifting
@@ -853,7 +856,7 @@ public:
 				}
 				else{
 					continue;
-				}				
+				}
 			}
 			if(PROGRESS) progress = (double)(i+1)/ B * 100;
 		}
@@ -973,7 +976,7 @@ public:
 		free(temp);
 		return true;
 	}
-	
+
 	/// Crop a region of the image and save it to a new file.
 	/// @param outfile is the file name for the new cropped image
@@ -1035,7 +1038,7 @@ public:
 	/// @param outfile is the file name for the output hyperspectral image (with trimmed bands)
 	/// @param b is an array of bands to be eliminated
 	void trim(std::string outfile, std::vector<size_t> band_array, bool PROGRESS = false){
-		
+
 		std::ofstream out(outfile.c_str(), std::ios::binary);	//open the output file for writing
 		file.seekg(0, std::ios::beg);							//move to the beginning of the input file
@@ -1185,7 +1188,7 @@ public:
 	void deriv(std::string outfile, size_t d, size_t order, const std::vector<double> w = std::vector<double>(), unsigned char* mask = NULL, bool PROGRESS = false){
 		std::ofstream out(outfile.c_str(), std::ios::binary);		//open the output file for writing
-				
+
 		size_t XY = X() * Y();										//calculate the number of values in a band
 		size_t XYb = XY * sizeof(T);								//calculate the size of a band (frame) in bytes
 		size_t B = Z();
@@ -1227,7 +1230,7 @@ public:
 					}
 				}
 			}
-			out.write((char*)outband, XYb);							//output the band			
+			out.write((char*)outband, XYb);							//output the band
 			if(PROGRESS) progress = (double)(b+1) / (double)B * 100;
 		}
@@ -32,6 +32,9 @@ protected:
 	std::vector<double> w;		//band wavelengths
+	using binary<T>::R;
+	using binary<T>::progress;
+
 	unsigned long long X(){ return R[O[0]]; }
 	unsigned long long Y(){ return R[O[1]]; }
 	unsigned long long Z(){ return R[O[2]]; }
@@ -70,8 +73,8 @@ protected:
 	}
 	/// Get the list of band numbers that bound a list of wavelengths
-	void band_bounds(std::vector<double> wavelengths, 
-					 std::vector<unsigned long long>& low_bands, 
+	void band_bounds(std::vector<double> wavelengths,
+					 std::vector<unsigned long long>& low_bands,
 					 std::vector<unsigned long long>& high_bands){
 		unsigned long long W = w.size();									//get the number of wavelengths in the list
@@ -92,7 +95,7 @@ protected:
 		band_bounds(wavelength, low, high);							//get the surrounding band indices
 		if(high == w.size()) return s[w.size()-1];					//if the high band is above the wavelength range, return the highest wavelength value
-		
+
 		return lerp(wavelength, s[low], w[low], s[high], w[high]);
 	}
@@ -140,11 +143,11 @@ public:
 		size_t XY = X() * Y();
 		memset(mask, 255, XY * sizeof(unsigned char));				//initialize the mask to zero (0)
 		T* page = (T*)malloc(R[0] * R[1] * sizeof(T));		//allocate space for a page of data
-		
+
 		for(size_t p = 0; p < R[2]; p++){					//for each page
 			binary<T>::read_page(page, p);					//read a page
 			for(size_t i = 0; i < R[0] * R[1]; i++){		//for each pixel in that page
-				
+
 #ifdef _WIN32
 				if(!_finite(page[i])){						//if the value at index i is finite
 #else
@@ -171,7 +174,7 @@ public:
 	}
 	/// Calculates the necessary image size required to combine two images given the specified offset (position) of the second image
-	void calc_combined_size(long long xp, long long yp, long long Sx, long long Sy, 
+	void calc_combined_size(long long xp, long long yp, long long Sx, long long Sy,
 							size_t& Sfx, size_t& Sfy,
 							size_t& p0_x, size_t& p0_y,
 							size_t& p1_x, size_t& p1_y){
@@ -179,7 +182,7 @@ public:
 		calc_combined_size(xp, yp, Sx, Sy, Sfx, Sfy);
 		p0_x = p0_y = p1_x = p1_y = 0;										//initialize all boundary positions to zero
-		
+
 		if(xp < 0) p0_x = -xp;												//set the left positions of the current and source image
 		else p1_x = xp;
 		if(yp < 0) p0_y = -yp;
@@ -191,8 +194,8 @@ public:
 		size_t w = X();										//calculate the number of pixels in a line
 		size_t wb = w * sizeof(T);							//calculate the number of bytes in a line
-		pw = (X() + x0 + x1);								//calculate the width of the padded image
-		pwb = pw * sizeof(T);								//calculate the number of bytes in a line of the padded image
+		size_t pw = (X() + x0 + x1);								//calculate the width of the padded image
+		size_t pwb = pw * sizeof(T);								//calculate the number of bytes in a line of the padded image
 		for(size_t y = 0; y < Y(); y++){								//for each line in the real image
 			memcpy( &padded[ (y + y0) * pw + x0 ], &src[y * w], wb );	//use memcpy to copy the line to the appropriate place in the padded image
@@ -68,7 +68,7 @@ public:
 			stim::image<T> I(file_list[i].str());
 			//retrieve the interlaced data from the image - store it in the grid
-			I.data_interleaved(&ptr[ i * R[0] * R[1] * R[2] ]);
+			I.get_interleaved_rgb(&ptr[ i * R[0] * R[1] * R[2] ]);
 		}
 	}
@@ -32,7 +32,7 @@ class image{
 	void unalloc(){								//frees any resources associated with the image
 		if(img)	free(img);						//if memory has been allocated, free it
 	}
-			
+
 	void clear(){								//clears all image data
 		unalloc();								//unallocate previous memory
@@ -53,7 +53,7 @@ class image{
 	size_t idx(size_t x, size_t y, size_t c = 0){
 		return y * C() * X() + x * C() + c;
 	}
-	
+
 	int cv_type(){
 		if(std::is_same<T, unsigned char>::value)	return CV_MAKETYPE(CV_8U, (int)C());
@@ -63,7 +63,7 @@ class image{
 		if(std::is_same<T, int>::value)				return CV_MAKETYPE(CV_32S, (int)C());
 		if(std::is_same<T, float>::value)			return CV_MAKETYPE(CV_32F, (int)C());
 		if(std::is_same<T, double>::value)			return CV_MAKETYPE(CV_64F, (int)C());
-		
+
 		std::cout<<"ERROR in stim::image::cv_type - no valid data type found"<<std::endl;
 		exit(1);
 	}
@@ -94,7 +94,7 @@ public:
 	}
 	/// Create a new image from scratch given a number of samples and channels
-	image(size_t x, size_t y = 1, size_t c = 1){			
+	image(size_t x, size_t y = 1, size_t c = 1){
 		allocate(x, y, c);
 	}
@@ -148,7 +148,7 @@ public:
 		if(C() == 1)
 			memcpy(buffer, img, bytes());
 		else if(C() == 3)
-			data_interleaved_bgr(buffer);
+			get_interleaved_bgr(buffer);
 		cv::Mat cvImage((int)Y(), (int)X(), cv_type(), buffer);
 		cv::imwrite(filename, cvImage);
 	}
@@ -170,7 +170,7 @@ public:
 		}
 	}
-	void data_interleaved_bgr(T* data){
+	void get_interleaved_bgr(T* data){
 		//for each channel
 		for(size_t y = 0; y < Y(); y++){
@@ -182,6 +182,11 @@ public:
 		}
 	}
+	void get_interleaved_rgb(T* data){
+		memcpy(data, img, bytes());
+	}
+
+
 	image<T> channel(size_t c){
 		//create a new image
@@ -284,7 +289,7 @@ public:
 		return s;			//return the index list
 	}
-	
+
 	/// Returns the maximum pixel value in the image
 	T maxv(){
 		T max_val = img[0];				//initialize the maximum value to the first one
@@ -294,10 +299,10 @@ public:
 			if (img[n] > max_val){			//if the value is higher than the current max
 				max_val = img[n];
-			}	
+			}
 		}
-		return max;	
+		return max_val;
 	}
 	/// Returns the maximum pixel value in the image
@@ -308,10 +313,10 @@ public:
 		for (size_t n=0; n<N; n++){		//for every value
 			if (img[n] < min_val){			//if the value is higher than the current max
 				min_val = img[n];
-			}	
+			}
 		}
-		return max;	
+		return min_val;
 	}
 	/// Invert an image by calculating I1 = alpha - I0, where alpha is the maximum image value
@@ -320,7 +325,7 @@ public:
 		image<T> r(X(), Y(), C());				//allocate space for the resulting image
 		for(size_t n = 0; n < N; n++)
 			r.img[n] = white_val - img[n];		//perform the inversion
-		
+
 		return r;								//return the inverted image
 	}
@@ -330,7 +335,7 @@ public:
 	}
 	image<T> convolve2(image<T> mask){
-		
+
 		std::cout<<"ERROR stim::image::convolve2 - function has been broken, and shouldn't really be in here."<<std::endl;
 		exit(1);
 	}
@@ -347,7 +352,7 @@ public:
 		std::cout<<"ERROR stim::image::set_interleaved3 - stim::image no longer supports 3D images."<<std::endl;
 		exit(1);
 	}
-	
+
 };
 };		//end namespace stim
@@ -175,5 +175,4 @@ public:
 };
 }
->>>>>>> origin/Master_Clone_W_Branch
 #endif
@@ -22,6 +22,9 @@
 #include <iostream>
 #include "../parser/parser.h"
+#ifdef BOOST_PRECOMPILED
+#include <boost/filesystem.hpp>
+#endif
 namespace stim{
 //filename class designed to work with both Windows and Unix
@@ -234,8 +237,40 @@ public:
 		}
 		return file_list;
 #else
-	std::cout<<"File lists aren't currently supported on Unix/Linux systems."<<std::endl;
-	exit(1);
+
+		boost::filesystem::path p(dir());	//create a path from the current filename
+		std::vector<stim::filename> file_list;
+		if(boost::filesystem::exists(p)){
+			if(boost::filesystem::is_directory(p)){
+				typedef std::vector<boost::filesystem::path> vec;             // store paths,
+				vec v;                                // so we can sort them later
+				std::copy(boost::filesystem::directory_iterator(p), boost::filesystem::directory_iterator(), back_inserter(v));
+				std::sort(v.begin(), v.end());             // sort, since directory iteration
+				  // is not ordered on some file systems
+				//compare file names to the current template (look for wild cards)
+				for (vec::const_iterator it(v.begin()), it_end(v.end()); it != it_end; ++it)
+				{
+					//if the filename is a wild card *or* it matches the read file name
+					if( prefix == "*" || prefix == (*it).filename().stem().string()){
+						//if the extension is a wild card *or* it matches the read file extension
+						if( ext == "*" || "." + ext == (*it).filename().extension().string()){
+							file_list.push_back((*it).string());	//include it in the final list
+						}
+
+					}
+
+
+
+				}
+
+			}
+
+		}
+		return file_list;
+
+
+//	std::cout<<"File lists aren't currently supported on Unix/Linux systems."<<std::endl;
+//	exit(1);
 #endif
 	}
 	//**************************************************************************************************
@@ -164,7 +164,7 @@ class cylinder
 		///Returns the number of points on the cylinder centerline
 		unsigned int size(){
-			return pos.size();
+			return e.size();
 		}
@@ -231,7 +231,7 @@ class cylinder
 		/// Adds a new magnitude value to all points
 		/// @param m is the starting value for the new magnitude
 		void add_mag(T m = 0){
-			for(unsigned int p = 0; p < pos.size(); p++)
+			for(unsigned int p = 0; p < e.size(); p++)
 				mags[p].push_back(m);
 		}
@@ -317,7 +317,7 @@ class cylinder
 			T len = L[0];						//allocate space for the segment length
 			//for every consecutive point in the cylinder
-			for(unsigned p = 1; p < pos.size(); p++){
+			for(unsigned p = 1; p < e.size(); p++){
 				//p1 = pos[p];							//get the position and magnitude for the next point
 				m1 = mags[p][m];
@@ -347,7 +347,7 @@ class cylinder
 			std::vector< vec<T> > result;
-			vec<T> p0 = e[i].P;								//initialize p0 to the first point on the centerline
+			vec<T> p0 = e[0].P;								//initialize p0 to the first point on the centerline
 			vec<T> p1;
 			unsigned N = size();							//number of points in the current centerline
@@ -158,6 +158,7 @@ protected:
 		using std::vector<triplet>::size;
 		using std::vector<triplet>::at;
 		using std::vector<triplet>::push_back;
+		using std::vector<triplet>::resize;
 		geometry(){}