diff --git a/matlab/loadAgilent.m b/matlab/loadAgilent.m new file mode 100644 index 0000000..4067189 --- /dev/null +++ b/matlab/loadAgilent.m @@ -0,0 +1,17 @@ +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 + + diff --git a/stim/envi/agilent_binary.h b/stim/envi/agilent_binary.h index d0d31d2..8c95921 100644 --- a/stim/envi/agilent_binary.h +++ b/stim/envi/agilent_binary.h @@ -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* dst, const agilent_binary* 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 fft(float band_min, float band_max){ + agilent_binary 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."<* cpu_fft = (std::complex*) malloc( R[0] * R[1] * (R[2]/2+1) * sizeof(std::complex) ); 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 result(R[0], R[1], size_i); result.Z[0] = start_i * fft_delta; //set the range for the FFT result diff --git a/stim/envi/envi.h b/stim/envi/envi.h index 5701fa3..994bc64 100644 --- a/stim/envi/envi.h +++ b/stim/envi/envi.h @@ -7,6 +7,7 @@ #include "../envi/bil.h" #include "../math/fd_coefficients.h" #include +#include //#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. diff --git a/stim/envi/envi_header.h b/stim/envi/envi_header.h index 2555a26..98d456f 100644 --- a/stim/envi/envi_header.h +++ b/stim/envi/envi_header.h @@ -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 band_index(double w){ std::vector 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 diff --git a/stim/math/vec3.h b/stim/math/vec3.h index 99e3809..6ddf28c 100644 --- a/stim/math/vec3.h +++ b/stim/math/vec3.h @@ -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; } -- libgit2 0.21.4