bug fixes for FFT in stim::agilent_binary

David Mayerich
1 parent e318123b
Showing 5 changed files with 81 additions and 17 deletions Show diff stats
matlab/loadAgilent.m
stim/envi/agilent_binary.h
stim/envi/envi.h
stim/envi/envi_header.h
stim/math/vec3.h
+function S = loadAgilent(filename)
+
+fid = fopen(filename);          %open the file for reading
+fseek(fid, 9, 'bof');           %skip past the first 9 bytes of the header
+
+bands = fread(fid, 1, 'int16'); %read the number of bands in the file
+fseek(fid, 13, 'cof');          %skip the next 13 bytes in the header
+
+samples = fread(fid, 1, 'int16');   %read the number of samples (X)
+lines = fread(fid, 1, 'int16');     %read the number of lines (Y)
+
+fseek(fid, 1020, 'bof');            %skip past the entire header
+S = fread(fid, [samples lines*bands], 'float32');   %read all the data
+S = reshape(S, [samples, lines, bands]);
+fclose(fid);                    %close the file
+
+
@@ -35,26 +35,28 @@ public:
 	void alloc(){
 		ptr = (T*) malloc(bytes());
 	}
-	void alloc(short x, short y, short z){
+	void alloc(size_t x, size_t y, size_t z){
 		R[0] = x;
 		R[1] = y;
 		R[2] = z;
 		alloc();
 	}
  
+	/// Create a deep copy of an agileng_binary object
 	void deep_copy(agilent_binary<T>* dst, const agilent_binary<T>* src){
 		dst->alloc(src->R[0], src->R[1], src->R[2]);			//allocate memory
 		memcpy(dst->ptr, src->ptr, bytes());					//copy the data
 		memcpy(dst->Z, src->Z, sizeof(double) * 2);				//copy the data z range
 	}
  
+	/// Default constructor, sets the resolution to zero and the data pointer to NULL
 	agilent_binary(){
-		memset(R, 0, sizeof(short) * 3);				//set the resolution to zero
+		memset(R, 0, sizeof(size_t) * 3);				//set the resolution to zero
 		ptr = NULL;
 	}
  
 	/// Constructor with resolution
-	agilent_binary(short x, short y, short z){
+	agilent_binary(size_t x, size_t y, size_t z){
 		alloc(x, y, z);
 	}
  
@@ -109,13 +111,11 @@ public:
  
 		char zero = 0;
 		for(size_t i = 0; i < 9; i++) outfile.write(&zero, 1);		//write 9 zeros
-		outfile.write((char*)&R[0], 2);
+		outfile.write((char*)&R[2], 2);
 		for(size_t i = 0; i < 13; i++) outfile.write(&zero, 1);		//write 13 zeros
+		outfile.write((char*)&R[0], 2);
 		outfile.write((char*)&R[1], 2);
-		outfile.write((char*)&R[2], 2);
 		for(size_t i = 0; i < 992; i++) outfile.write(&zero, 1);		//write 992 zeros
-		//char zerovec[1020];
-		//outfile.write((char*)zerovec, 1020);
  
 		size_t b = bytes();
 		outfile.write((char*)ptr, b);							//write the data to the output file
@@ -149,7 +149,7 @@ public:
  
 #ifdef CUDA_FOUND
 	/// Perform an FFT and return a binary file with bands in the specified range
-	agilent_binary<T> fft(float band_min, float band_max){
+	agilent_binary<T> fft(double band_min, double band_max, double ELWN = 15798, int UDR = 2){
 		auto total_start = std::chrono::high_resolution_clock::now();
  
 		auto start = std::chrono::high_resolution_clock::now();
@@ -177,8 +177,8 @@ public:
  
 		start = std::chrono::high_resolution_clock::now();
 		int N[1];					//create an array with the interferogram size (required for cuFFT input)
-		N[0] = R[2];				//set the only array value to the interferogram size
-		if(cufftPlanMany(&plan, 1, N, NULL, 1, R[2], NULL, 1, R[2], CUFFT_R2C, batch) != CUFFT_SUCCESS){
+		N[0] = (int)R[2];				//set the only array value to the interferogram size
+		if(cufftPlanMany(&plan, 1, N, NULL, 1, (int)R[2], NULL, 1, (int)R[2], CUFFT_R2C, (int)batch) != CUFFT_SUCCESS){
 			std::cout<<"cuFFT Error: unable to create 1D plan."<<std::endl;
 			exit(1);
 		}
@@ -199,12 +199,13 @@ public:
 		std::complex<T>* cpu_fft = (std::complex<T>*) malloc( R[0] * R[1] * (R[2]/2+1) * sizeof(std::complex<T>) );
 		HANDLE_ERROR(cudaMemcpy(cpu_fft, gpu_fft, R[0] * R[1] * (R[2]/2+1) * sizeof(cufftComplex), cudaMemcpyDeviceToHost));	//copy data from the host to the device
  
-		double int_delta = 0.00012656;									//interferogram sample spacing in centimeters
+		//double int_delta = 0.00012656;									//interferogram sample spacing in centimeters
+		double int_delta = (1.0 / ELWN) * ((double)UDR / 2.0);			//calculate the interferogram spacing
 		double int_length = int_delta * R[2];							//interferogram length in centimeters
 		double fft_delta = 1/int_length;								//spectrum spacing (in inverse centimeters, wavenumber
  
-		size_t start_i = std::ceil(band_min / fft_delta);				//calculate the first band to store
-		size_t size_i = std::floor(band_max / fft_delta) - start_i;		//calculate the number of bands to store
+		size_t start_i = (size_t)std::ceil(band_min / fft_delta);				//calculate the first band to store
+		size_t size_i = (size_t)std::floor(band_max / fft_delta) - start_i;		//calculate the number of bands to store
 		size_t end_i = start_i + size_i;								//last band number
 		agilent_binary<T> result(R[0], R[1], size_i);
 		result.Z[0] = start_i * fft_delta;								//set the range for the FFT result
@@ -7,6 +7,7 @@
 #include "../envi/bil.h"
 #include "../math/fd_coefficients.h"
 #include <iostream>
+#include <fstream>
 //#include "../image/image.h"
  
 namespace stim{
@@ -224,6 +225,37 @@ public:
  
 	}
  
+	/// Open an Agilent binary file as an ENVI stream
+	bool open_agilent(std::string filename){
+		fname = filename;												//store the file name
+
+		//Open the file temporarily to get the header information
+		FILE* f = fopen(filename.c_str(), "r");							//open the binary file for reading
+		if(f == NULL) return false;										//return false if no file is opened
+
+		fseek(f, 9, SEEK_SET);											//seek to the number of bands
+		short b;														//allocate space for the number of bands
+		fread(&b, sizeof(short), 1, f);									//read the number of bands
+		fseek(f, 13, SEEK_CUR);											//skip the the x and y dimensions
+		short x, y;
+		fread(&x, sizeof(short), 1, f);									//read the image x and y size
+		fread(&y, sizeof(short), 1, f);
+		fclose(f);														//close the file
+
+		//store the information from the Agilent header in the ENVI header
+		header.bands = b;
+		header.samples = x;
+		header.lines = y;
+		header.data_type = envi_header::float32;						//all values are 32-bit floats
+		header.header_offset = 1020;									//number of bytes in an Agilent binary header
+		header.interleave = envi_header::BSQ;							//all Agilent binary files are BSQ
+
+		allocate();														//allocate the streaming file object
+		open();															//open the file for streaming
+
+		return true;
+	}
+
 	/// Open an existing ENVI file given the filename and a header structure
  
 	/// @param filename is the name of the ENVI binary file
@@ -257,7 +289,6 @@ public:
 		//header.load(headername);
  
 		return open(filename, h);
-
 	}
  
 	/// Normalize a hyperspectral ENVI file given a band number and threshold.
@@ -440,9 +440,24 @@ struct envi_header
 	}
  
 	/// Convert a wavelength to a band index (or a pair of surrounding band indices)
+	///		if the file doesn't specify wavelengths, w is assumed to be a band index
 	std::vector<size_t> band_index(double w){
 		std::vector<size_t> idx;										//create an empty array of indices
-		if(w < wavelength[0] || w > wavelength[bands-1]) return idx;	//if the wavelength range is outside of the file, return an empty array
+		if(wavelength.size() == 0){										//if a wavelength vector doesn't exist, assume the passed value is a band
+			if(w < 0 || w > bands-1) return idx;						//if the band is outside the given band range, return an empty vector
+			size_t low, high;											//allocate space for the floor and ceiling
+			low = std::floor(w);										//calculate the floor
+			high = std::ceil(w);										//calculate the ceiling
+			if(low == high)												//if the floor and ceiling are the same
+				idx.push_back(low);										//return a vector with one element (the given w matches a band exactly)
+			else{
+				idx.resize(2);											//otherwise return the floor and ceiling
+				idx[0] = low;
+				idx[1] = high;
+			}
+			return idx;
+		}
+		else if(w < wavelength[0] || w > wavelength[bands-1]) return idx;	//if the wavelength range is outside of the file, return an empty array
  
 		for(size_t b = 0; b < bands; b++){							//for each band in the wavelength vector
 			if(wavelength[b] == w){									//if an exact match is found
@@ -216,7 +216,7 @@ public:
 		return result;
 	}
  
-#ifndef __NVCC__
+//#ifndef __NVCC__
 	/// Outputs the vector as a string
 	std::string str() const{
 		std::stringstream ss;
@@ -234,7 +234,7 @@ public:
  
 		return ss.str();
 	}
-#endif
+//#endif
  
 	size_t size(){ return 3; }