25d8d20b
David Mayerich
added FFT support...
|
1
2
3
4
5
6
|
//make sure that this header file is only loaded once
#ifndef STIM_AGILENT_BINARY_H
#define STIM_AGILENT_BINARY_H
#include <string>
#include <fstream>
|
b4cea01a
David Mayerich
removed automated...
|
7
|
#include <complex>
|
e38fc0c3
David Mayerich
header file wavel...
|
8
9
|
#include <cstring>
#include <chrono>
|
25d8d20b
David Mayerich
added FFT support...
|
10
11
|
//CUDA
|
2781a48c
David Mayerich
fixed cuBLAS func...
|
12
13
14
15
|
//#ifdef CUDA_FOUND
#include <cuda_runtime.h>
#include "cufft.h"
#include <stim/cuda/cudatools/error.h>
|
9d34229e
David Mayerich
restructured the ...
|
16
|
#include <stim/envi/envi_header.h>
|
2781a48c
David Mayerich
fixed cuBLAS func...
|
17
|
//#endif
|
25d8d20b
David Mayerich
added FFT support...
|
18
19
20
21
22
23
24
25
|
namespace stim{
template<typename T>
class agilent_binary{
protected:
std::string fname;
|
29165fa9
David Mayerich
added ROC curve c...
|
26
27
28
29
|
T* ptr; //pointer to the image data
size_t R[3]; //size of the binary image in X, Y, and Z
static const size_t header = 1020; //header size
double Z[2]; //range of z values (position or wavelength)
|
25d8d20b
David Mayerich
added FFT support...
|
30
31
32
33
34
35
36
37
38
39
|
public:
size_t size(){
return (size_t)R[0] * (size_t)R[1] * (size_t)R[2];
}
size_t bytes(){
return size() * sizeof(T);
}
void alloc(){
|
00a7e339
David Mayerich
fixed agilent bin...
|
40
|
if (ptr != NULL) free(ptr);
|
2097dd82
David Mayerich
fixed agilent_bin...
|
41
|
ptr = NULL;
|
25d8d20b
David Mayerich
added FFT support...
|
42
43
|
ptr = (T*) malloc(bytes());
}
|
00a7e339
David Mayerich
fixed agilent bin...
|
44
|
void alloc(size_t x, size_t y, size_t z){
|
25d8d20b
David Mayerich
added FFT support...
|
45
46
47
48
49
50
|
R[0] = x;
R[1] = y;
R[2] = z;
alloc();
}
|
332925ad
David Mayerich
added the ability...
|
51
52
53
54
|
char* data() {
return (char*)ptr;
}
|
2781a48c
David Mayerich
fixed cuBLAS func...
|
55
56
57
58
|
size_t dim(size_t i){
return R[i];
}
|
48f356fa
David Mayerich
bug fixes for FFT...
|
59
|
/// Create a deep copy of an agileng_binary object
|
25d8d20b
David Mayerich
added FFT support...
|
60
61
62
63
64
65
|
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
}
|
48f356fa
David Mayerich
bug fixes for FFT...
|
66
|
/// Default constructor, sets the resolution to zero and the data pointer to NULL
|
25d8d20b
David Mayerich
added FFT support...
|
67
|
agilent_binary(){
|
25d8d20b
David Mayerich
added FFT support...
|
68
|
ptr = NULL;
|
2097dd82
David Mayerich
fixed agilent_bin...
|
69
70
|
memset(R, 0, sizeof(size_t) * 3); //set the resolution to zero
memset(Z, 0, sizeof(double) * 2);
|
25d8d20b
David Mayerich
added FFT support...
|
71
72
73
|
}
/// Constructor with resolution
|
48f356fa
David Mayerich
bug fixes for FFT...
|
74
|
agilent_binary(size_t x, size_t y, size_t z){
|
2097dd82
David Mayerich
fixed agilent_bin...
|
75
|
ptr = NULL;
|
25d8d20b
David Mayerich
added FFT support...
|
76
|
alloc(x, y, z);
|
e38fc0c3
David Mayerich
header file wavel...
|
77
|
memset(Z, 0, sizeof(double) * 2);
|
25d8d20b
David Mayerich
added FFT support...
|
78
79
80
81
82
|
}
/// Constructor with filename
agilent_binary(std::string filename){
ptr = NULL;
|
e38fc0c3
David Mayerich
header file wavel...
|
83
|
memset(Z, 0, sizeof(double) * 2);
|
25d8d20b
David Mayerich
added FFT support...
|
84
85
86
87
88
|
load(filename);
}
/// Copy constructor
agilent_binary(const agilent_binary<T> &obj){
|
2097dd82
David Mayerich
fixed agilent_bin...
|
89
|
ptr = NULL;
|
25d8d20b
David Mayerich
added FFT support...
|
90
91
92
93
94
95
96
97
98
99
|
deep_copy(this, &obj);
}
agilent_binary<T>& operator=(const agilent_binary<T> rhs){
if(this != &rhs){ //check for self-assignment
deep_copy(this, &rhs); //make a deep copy
}
return *this; //return the result
}
|
332925ad
David Mayerich
added the ability...
|
100
101
102
103
104
|
operator bool() {
if (R[0] == 0 || R[1] == 0 || R[2] == 0) return false;
else return true;
}
|
25d8d20b
David Mayerich
added FFT support...
|
105
|
~agilent_binary(){
|
2097dd82
David Mayerich
fixed agilent_bin...
|
106
107
|
if(ptr != NULL)
free(ptr);
|
25d8d20b
David Mayerich
added FFT support...
|
108
109
110
|
}
void load(std::string filename){
|
2097dd82
David Mayerich
fixed agilent_bin...
|
111
112
|
if(ptr != NULL) free(ptr); //if memory has been allocated, free it
ptr = NULL;
|
25d8d20b
David Mayerich
added FFT support...
|
113
|
|
2097dd82
David Mayerich
fixed agilent_bin...
|
114
|
fname = filename; //save the filename
|
25d8d20b
David Mayerich
added FFT support...
|
115
116
117
118
|
short x, y, z;
std::ifstream infile(fname, std::ios::binary); //open the input file
|
2097dd82
David Mayerich
fixed agilent_bin...
|
119
120
|
if (infile) {
infile.seekg(9, std::ios::beg); //seek past 9 bytes from the beginning of the file
|
25d8d20b
David Mayerich
added FFT support...
|
121
|
|
2097dd82
David Mayerich
fixed agilent_bin...
|
122
|
infile.read((char*)(&z), 2); //read two bytes of data (the number of samples is stored as a 16-bit integer)
|
25d8d20b
David Mayerich
added FFT support...
|
123
|
|
2097dd82
David Mayerich
fixed agilent_bin...
|
124
125
126
|
infile.seekg(13, std::ios::cur); //skip another 13 bytes
infile.read((char*)(&x), 2); //read the X and Y dimensions
infile.read((char*)(&y), 2);
|
25d8d20b
David Mayerich
added FFT support...
|
127
|
|
2097dd82
David Mayerich
fixed agilent_bin...
|
128
|
infile.seekg(header, std::ios::beg); //seek to the start of the data
|
25d8d20b
David Mayerich
added FFT support...
|
129
|
|
a688950c
David Mayerich
fixed a bug intro...
|
130
|
alloc(x, y, z); //allocate the data
|
2097dd82
David Mayerich
fixed agilent_bin...
|
131
|
infile.read((char*)ptr, bytes()); //read the data
|
a688950c
David Mayerich
fixed a bug intro...
|
132
|
infile.close(); //close the file
|
2097dd82
David Mayerich
fixed agilent_bin...
|
133
134
135
|
Z[0] = 1;
Z[1] = (double)R[2];
}
|
25d8d20b
David Mayerich
added FFT support...
|
136
137
138
139
140
141
142
|
}
void save(std::string filename){
std::ofstream outfile(filename, std::ios::binary); //create an output file
char zero = 0;
for(size_t i = 0; i < 9; i++) outfile.write(&zero, 1); //write 9 zeros
|
48f356fa
David Mayerich
bug fixes for FFT...
|
143
|
outfile.write((char*)&R[2], 2);
|
25d8d20b
David Mayerich
added FFT support...
|
144
|
for(size_t i = 0; i < 13; i++) outfile.write(&zero, 1); //write 13 zeros
|
48f356fa
David Mayerich
bug fixes for FFT...
|
145
|
outfile.write((char*)&R[0], 2);
|
25d8d20b
David Mayerich
added FFT support...
|
146
|
outfile.write((char*)&R[1], 2);
|
25d8d20b
David Mayerich
added FFT support...
|
147
|
for(size_t i = 0; i < 992; i++) outfile.write(&zero, 1); //write 992 zeros
|
25d8d20b
David Mayerich
added FFT support...
|
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
|
size_t b = bytes();
outfile.write((char*)ptr, b); //write the data to the output file
outfile.close();
}
stim::envi_header create_header(){
stim::envi_header header;
header.samples = R[0];
header.lines = R[1];
header.bands = R[2];
double z_delta = (double)(Z[1] - Z[0]) / (double)(R[2] - 1);
header.wavelength.resize(R[2]);
for(size_t i = 0; i < R[2]; i++)
header.wavelength[i] = i * z_delta + Z[0];
return header;
}
|
29165fa9
David Mayerich
added ROC curve c...
|
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
|
/// Subtract the mean from each pixel. Generally used for centering an interferogram.
void meancenter(){
size_t Z = R[2]; //store the number of bands
size_t XY = R[0] * R[1]; //store the number of pixels in the image
T sum = (T)0;
T mean;
for(size_t xy = 0; xy < XY; xy++){ //for each pixel
sum = 0;
for(size_t z = 0; z < Z; z++){ //for each band
sum += ptr[ z * XY + xy ]; //add the band value to a running sum
}
mean = sum / (T)Z; //calculate the pixel mean
for(size_t z = 0; z < Z; z++){
ptr[ z * XY + xy ] -= mean; //subtract the mean from each band
}
}
}
/// adds n bands of zero padding to the end of the file
void zeropad(size_t n){
size_t newZ = R[2] + n;
T* temp = (T*) calloc(R[0] * R[1] * newZ, sizeof(T)); //allocate space for the new image
memcpy(temp, ptr, size() * sizeof(T)); //copy the old data to the new image
free(ptr); //free the old data
ptr = temp; //swap in the new data
R[2] = newZ; //set the z-dimension to the new zero value
}
//pads to the nearest power-of-two
void zeropad(){
|
ebd02452
David Mayerich
added cuBLAS erro...
|
199
|
size_t newZ = (size_t)pow(2, ceil(log(R[2])/log(2))); //find the nearest power-of-two
|
29165fa9
David Mayerich
added ROC curve c...
|
200
|
size_t n = newZ - R[2]; //calculate the number of bands to add
|
29165fa9
David Mayerich
added ROC curve c...
|
201
202
203
|
zeropad(n); //add the padding
}
|
25d8d20b
David Mayerich
added FFT support...
|
204
205
206
207
208
209
210
211
212
213
214
|
/// Calculate the absorbance spectrum from the transmission spectrum given a background
void absorbance(stim::agilent_binary<T>* background){
size_t N = size(); //calculate the number of values to be ratioed
if(N != background->size()){
std::cerr<<"ERROR in stim::agilent_binary::absorbance() - transmission image size doesn't match background"<<std::endl;
exit(1);
}
for(size_t i = 0; i < N; i++)
ptr[i] = -log10(ptr[i] / background->ptr[i]);
}
|
332925ad
David Mayerich
added the ability...
|
215
216
217
218
219
220
221
222
223
224
225
|
//crops the image down to a set number of samples
void crop(size_t n) {
if (n < R[2]) { //if the requested size is smaller than the image
R[2] = n; //update the number of bands
T* old_ptr = ptr; //store the old pointer
alloc(); //allocate space for the new image
memcpy(ptr, old_ptr, bytes()); //copy the old data to the new image
free(old_ptr); //free the old data
}
}
|
2781a48c
David Mayerich
fixed cuBLAS func...
|
226
|
//#ifdef CUDA_FOUND
|
25d8d20b
David Mayerich
added FFT support...
|
227
|
/// Perform an FFT and return a binary file with bands in the specified range
|
bb457896
David Mayerich
improved error me...
|
228
|
agilent_binary<T> fft(double band_min, double band_max, double ELWN = 15798, int UDR = 2, int device = 0){
|
25d8d20b
David Mayerich
added FFT support...
|
229
230
231
232
233
234
235
236
237
238
239
240
241
242
|
auto total_start = std::chrono::high_resolution_clock::now();
auto start = std::chrono::high_resolution_clock::now();
T* cpu_data = (T*) malloc( bytes() ); //allocate space for the transposed data
for(size_t b = 0; b < R[2]; b++){
for(size_t x = 0; x < R[0] * R[1]; x++){
cpu_data[x * R[2] + b] = ptr[b * R[0] * R[1] + x];
}
}
auto end = std::chrono::high_resolution_clock::now();
std::chrono::duration<double> diff = end-start;
// std::cout << "Transpose data: " << diff.count() << " s\n";
start = std::chrono::high_resolution_clock::now();
|
bb457896
David Mayerich
improved error me...
|
243
244
245
|
if (device >= 0) { //if a CUDA device is specified
int dev_count;
HANDLE_ERROR(cudaGetDeviceCount(&dev_count)); //get the number of CUDA devices
|
fa8763a6
David Mayerich
removed CUDA devi...
|
246
|
//std::cout << "Number of CUDA devices: " << dev_count << std::endl; //output the number of CUDA devices
|
bb457896
David Mayerich
improved error me...
|
247
|
cudaDeviceProp prop;
|
fa8763a6
David Mayerich
removed CUDA devi...
|
248
|
//std::cout << "CUDA devices----" << std::endl;
|
bb457896
David Mayerich
improved error me...
|
249
250
251
|
for (int d = 0; d < dev_count; d++) { //for each CUDA device
cudaGetDeviceProperties(&prop, d); //get the property of the first device
//float cc = prop.major + prop.minor / 10.0f; //calculate the compute capability
|
fa8763a6
David Mayerich
removed CUDA devi...
|
252
|
//std::cout << d << ": [" << prop.major << "." << prop.minor << "] " << prop.name << std::endl; //display the device information
|
bb457896
David Mayerich
improved error me...
|
253
254
255
256
257
258
259
260
|
//if(cc > best_device_cc){
// best_device_cc = cc; //if this is better than the previous device, use it
// best_device_id = d;
//}
}
if (dev_count > 0 && dev_count > device) { //if the first device is not an emulator
cudaGetDeviceProperties(&prop, device); //get the property of the requested CUDA device
if (prop.major != 9999) {
|
fa8763a6
David Mayerich
removed CUDA devi...
|
261
|
//std::cout << "Using device " << device << std::endl;
|
bb457896
David Mayerich
improved error me...
|
262
263
264
265
|
HANDLE_ERROR(cudaSetDevice(device));
}
}
}
|
25d8d20b
David Mayerich
added FFT support...
|
266
267
268
269
270
271
272
273
274
275
276
277
278
|
cufftHandle plan; //allocate space for a cufft plan
cufftReal* gpu_data; //create a pointer to the data
size_t batch = R[0] * R[1]; //calculate the batch size (X * Y)
HANDLE_ERROR(cudaMalloc((void**)&gpu_data, bytes())); //allocate space on the GPU
HANDLE_ERROR(cudaMemcpy(gpu_data, cpu_data, bytes(), cudaMemcpyHostToDevice)); //copy the data to the GPU
cufftComplex* gpu_fft;
HANDLE_ERROR(cudaMalloc((void**)&gpu_fft, R[0] * R[1] * (R[2]/2 + 1) * sizeof(cufftComplex)));
end = std::chrono::high_resolution_clock::now();
diff = end-start;
//std::cout << "Allocate/copy: " << diff.count() << " s\n";
start = std::chrono::high_resolution_clock::now();
int N[1]; //create an array with the interferogram size (required for cuFFT input)
|
48f356fa
David Mayerich
bug fixes for FFT...
|
279
280
|
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){
|
25d8d20b
David Mayerich
added FFT support...
|
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
|
std::cout<<"cuFFT Error: unable to create 1D plan."<<std::endl;
exit(1);
}
end = std::chrono::high_resolution_clock::now();
diff = end-start;
//std::cout << "Create a plan: " << diff.count() << " s\n";
start = std::chrono::high_resolution_clock::now();
if (cufftExecR2C(plan, gpu_data, gpu_fft) != CUFFT_SUCCESS){ //execute the (implicitly forward) transform
std::cout<<"CUFFT error: ExecR2C Forward failed";
exit(1);
}
end = std::chrono::high_resolution_clock::now();
diff = end-start;
//std::cout << "Perform FFT: " << diff.count() << " s\n";
start = std::chrono::high_resolution_clock::now();
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
|
48f356fa
David Mayerich
bug fixes for FFT...
|
301
|
//double int_delta = 0.00012656; //interferogram sample spacing in centimeters
|
332925ad
David Mayerich
added the ability...
|
302
303
304
305
|
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)
double fft_max = fft_delta * R[2]/2; //get the maximum wavenumber value supported by the specified number of interferogram samples
|
550bb945
David Mayerich
fixed UDR problem...
|
306
|
|
332925ad
David Mayerich
added the ability...
|
307
308
|
if(band_max > fft_max) band_max = fft_max; //the user gave a band outside of the FFT range, reset the band to the maximum available
if (band_min < 0) band_min = 0;
|
25d8d20b
David Mayerich
added FFT support...
|
309
|
|
48f356fa
David Mayerich
bug fixes for FFT...
|
310
311
|
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
|
332925ad
David Mayerich
added the ability...
|
312
|
size_t end_i = start_i + size_i; //last band number
|
25d8d20b
David Mayerich
added FFT support...
|
313
|
agilent_binary<T> result(R[0], R[1], size_i);
|
332925ad
David Mayerich
added the ability...
|
314
|
result.Z[0] = start_i * fft_delta; //set the range for the FFT result
|
25d8d20b
David Mayerich
added FFT support...
|
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
|
result.Z[1] = end_i * fft_delta;
for(size_t b = start_i; b < end_i; b++){
for(size_t x = 0; x < R[0] * R[1]; x++){
result.ptr[(b - start_i) * R[0] * R[1] + x] = abs(cpu_fft[x * (R[2]/2+1) + b]);
}
}
end = std::chrono::high_resolution_clock::now();
diff = end-start;
//std::cout << "Transpose/Crop: " << diff.count() << " s\n";
auto total_end = std::chrono::high_resolution_clock::now();
diff = total_end-total_start;
cufftDestroy(plan);
HANDLE_ERROR(cudaFree(gpu_data));
HANDLE_ERROR(cudaFree(gpu_fft));
free(cpu_data);
free(cpu_fft);
return result;
}
|
9d34229e
David Mayerich
restructured the ...
|
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
|
//saves the binary as an ENVI file with a BIP interleave format
int bip(T* bip_ptr){
//std::ofstream out(outfile.c_str(), std::ios::binary); //create a binary file stream for output
size_t XY = R[0] * R[1];
size_t B = R[2];
size_t b;
for(size_t xy = 0; xy < XY; xy++){
for(b = 0; b < B; b++){
bip_ptr[xy * B + b] = ptr[b * XY + xy];
}
}
return 0;
}
|
2781a48c
David Mayerich
fixed cuBLAS func...
|
352
|
//#endif
|
25d8d20b
David Mayerich
added FFT support...
|
353
354
355
356
357
358
|
};
}
#endif
|