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);
|
25d8d20b
David Mayerich
added FFT support...
|
41
42
|
ptr = (T*) malloc(bytes());
}
|
00a7e339
David Mayerich
fixed agilent bin...
|
43
|
void alloc(size_t x, size_t y, size_t z){
|
25d8d20b
David Mayerich
added FFT support...
|
44
45
46
47
48
49
|
R[0] = x;
R[1] = y;
R[2] = z;
alloc();
}
|
332925ad
David Mayerich
added the ability...
|
50
51
52
53
|
char* data() {
return (char*)ptr;
}
|
2781a48c
David Mayerich
fixed cuBLAS func...
|
54
55
56
57
|
size_t dim(size_t i){
return R[i];
}
|
48f356fa
David Mayerich
bug fixes for FFT...
|
58
|
/// Create a deep copy of an agileng_binary object
|
25d8d20b
David Mayerich
added FFT support...
|
59
60
61
62
63
64
|
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...
|
65
|
/// Default constructor, sets the resolution to zero and the data pointer to NULL
|
25d8d20b
David Mayerich
added FFT support...
|
66
|
agilent_binary(){
|
48f356fa
David Mayerich
bug fixes for FFT...
|
67
|
memset(R, 0, sizeof(size_t) * 3); //set the resolution to zero
|
e38fc0c3
David Mayerich
header file wavel...
|
68
|
memset(Z, 0, sizeof(double) * 2);
|
25d8d20b
David Mayerich
added FFT support...
|
69
70
71
72
|
ptr = NULL;
}
/// Constructor with resolution
|
48f356fa
David Mayerich
bug fixes for FFT...
|
73
|
agilent_binary(size_t x, size_t y, size_t z){
|
25d8d20b
David Mayerich
added FFT support...
|
74
|
alloc(x, y, z);
|
e38fc0c3
David Mayerich
header file wavel...
|
75
|
memset(Z, 0, sizeof(double) * 2);
|
25d8d20b
David Mayerich
added FFT support...
|
76
77
78
79
80
|
}
/// Constructor with filename
agilent_binary(std::string filename){
ptr = NULL;
|
e38fc0c3
David Mayerich
header file wavel...
|
81
|
memset(Z, 0, sizeof(double) * 2);
|
25d8d20b
David Mayerich
added FFT support...
|
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
|
load(filename);
}
/// Copy constructor
agilent_binary(const agilent_binary<T> &obj){
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...
|
97
98
99
100
101
|
operator bool() {
if (R[0] == 0 || R[1] == 0 || R[2] == 0) return false;
else return true;
}
|
25d8d20b
David Mayerich
added FFT support...
|
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
|
~agilent_binary(){
free(ptr);
}
void load(std::string filename){
if(ptr != NULL) free(ptr); //if memory has been allocated, free it
fname = filename; //save the filename
short x, y, z;
std::ifstream infile(fname, std::ios::binary); //open the input file
infile.seekg(9, std::ios::beg); //seek past 9 bytes from the beginning of the file
infile.read((char*)(&z), 2); //read two bytes of data (the number of samples is stored as a 16-bit integer)
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);
infile.seekg(header, std::ios::beg); //seek to the start of the data
alloc(x, y, z);
ptr = (T*) malloc(bytes()); //allocate space for the data
infile.read((char*)ptr, bytes()); //read the data
infile.close();
|
e38fc0c3
David Mayerich
header file wavel...
|
128
|
Z[0] = 1;
|
3c17c706
David Mayerich
added a crop func...
|
129
|
Z[1] = (double)R[2];
|
25d8d20b
David Mayerich
added FFT support...
|
130
131
132
133
134
135
136
|
}
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...
|
137
|
outfile.write((char*)&R[2], 2);
|
25d8d20b
David Mayerich
added FFT support...
|
138
|
for(size_t i = 0; i < 13; i++) outfile.write(&zero, 1); //write 13 zeros
|
48f356fa
David Mayerich
bug fixes for FFT...
|
139
|
outfile.write((char*)&R[0], 2);
|
25d8d20b
David Mayerich
added FFT support...
|
140
|
outfile.write((char*)&R[1], 2);
|
25d8d20b
David Mayerich
added FFT support...
|
141
|
for(size_t i = 0; i < 992; i++) outfile.write(&zero, 1); //write 992 zeros
|
25d8d20b
David Mayerich
added FFT support...
|
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
|
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...
|
162
163
164
165
166
167
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
|
/// 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...
|
193
|
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...
|
194
|
size_t n = newZ - R[2]; //calculate the number of bands to add
|
29165fa9
David Mayerich
added ROC curve c...
|
195
196
197
|
zeropad(n); //add the padding
}
|
25d8d20b
David Mayerich
added FFT support...
|
198
199
200
201
202
203
204
205
206
207
208
|
/// 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...
|
209
210
211
212
213
214
215
216
217
218
219
|
//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...
|
220
|
//#ifdef CUDA_FOUND
|
25d8d20b
David Mayerich
added FFT support...
|
221
|
/// Perform an FFT and return a binary file with bands in the specified range
|
48f356fa
David Mayerich
bug fixes for FFT...
|
222
|
agilent_binary<T> fft(double band_min, double band_max, double ELWN = 15798, int UDR = 2){
|
25d8d20b
David Mayerich
added FFT support...
|
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
|
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();
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...
|
250
251
|
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...
|
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
|
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...
|
272
|
//double int_delta = 0.00012656; //interferogram sample spacing in centimeters
|
332925ad
David Mayerich
added the ability...
|
273
274
275
276
|
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...
|
277
|
|
332925ad
David Mayerich
added the ability...
|
278
279
|
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...
|
280
|
|
48f356fa
David Mayerich
bug fixes for FFT...
|
281
282
|
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...
|
283
|
size_t end_i = start_i + size_i; //last band number
|
25d8d20b
David Mayerich
added FFT support...
|
284
|
agilent_binary<T> result(R[0], R[1], size_i);
|
332925ad
David Mayerich
added the ability...
|
285
|
result.Z[0] = start_i * fft_delta; //set the range for the FFT result
|
25d8d20b
David Mayerich
added FFT support...
|
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
|
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 ...
|
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
|
//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...
|
323
|
//#endif
|
25d8d20b
David Mayerich
added FFT support...
|
324
325
326
327
328
329
|
};
}
#endif
|