5f3cba02
 David Mayerich
initial public co...
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
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
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
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
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
|
//#define _CRTDBG_MAP_ALLOC
//#include <stdlib.h>
//#include <crtdbg.h>
#include <iostream>
#include <stim/parser/arguments.h>
#include <stim/envi/agilent_binary.h>
#include <stim/parser/filename.h>
stim::arglist args; //generate a list of accepted arguments
#include <thread>
void progress_thread_double(double* e); //progress bar threaded function
bool fft = false; //perform an FFT before storing
double wn_low, wn_high;
float elwn;
int udr;
bool absorbance = false;
stim::agilent_binary<float> background; //store the background data
bool crop_samples = false;
size_t n_samples;
void mosaic_coord(unsigned int &x, unsigned int &y, stim::filename dmd_filename) {
//get the file prefix
std::string prefix = dmd_filename.get_prefix();
//these files are in the format:
// ???????????_xxxxx_yyyyy.dmd
//find the position of the last underscore
size_t y_start = prefix.find_last_of("_");
//remove all values between x_start and the end of the string
std::string y_string = prefix.substr(y_start + 1);
y = atoi(y_string.c_str());
//crop out the x coordinate and the last underscore
prefix = prefix.substr(0, y_start);
//find the y coordinate using the same method as above
size_t x_start = prefix.find_last_of("_");
std::string x_string = prefix.substr(x_start + 1);
x = atoi(x_string.c_str());
}
//Use the file list to find the number of tiles along each dimension of the mosaic
void mosaic_size(unsigned int &sx, unsigned int &sy, std::vector<stim::filename> dmd_list) {
sx = sy = 0;
unsigned int xi, yi;
//for each dmd file
for (unsigned int i = 0; i < dmd_list.size(); i++) {
//find the coordinate of the tile represented by this file
mosaic_coord(xi, yi, dmd_list[i]);
//if the coordinate is larger than the current max, store it as the current max
if (xi > sx) sx = xi;
if (yi > sy) sy = yi;
}
//add 1 to reflect that the indexing starts at 0
sx++; sy++;
}
void save_mosaic(std::vector<stim::filename> tile_list, std::string outfile, stim::agilent_binary<float>* background = NULL) {
stim::agilent_binary<float> binfile(tile_list[0]); //load the first agilent binary file
size_t ds = binfile.dim(0); //get the detector size
unsigned int mx, my;
mosaic_size(mx, my, tile_list); //calculate the mosaic size
std::cout << "Building a ["<<mx<<" x "<<my<<"] mosaic with a detector size of (" << ds << " x " << ds <<") pixels"<< std::endl;
std::cout << "Mosaic size: " << mx << " " << my << std::endl; //output the mosaic size
std::ofstream out(outfile, std::ios::binary); //open a binary file for writing
size_t nxo = mx * ds; //number of output pixels along x
size_t nyo = my * ds; //number of output pixels along y
size_t nbo;
size_t bytes_tile_band = ds * ds * sizeof(float); //number of bytes in a tile band
size_t bytes_out_band = bytes_tile_band * mx * my; //number of bytes in an output image band
double progress = 0.0;
std::thread t1(progress_thread_double, &progress); //start the progress bar thread
size_t xo, yo, io;
size_t ii;
stim::envi_header header;
bool header_flag = true;
for (size_t ty = 0; ty < my; ty++) {
for (size_t tx = 0; tx < mx; tx++) {
xo = tx * ds; //calculate the x output position for the current tile
stim::agilent_binary<float> I(tile_list[tx * my + ty].str());
if (crop_samples) I.crop(n_samples); //crop samples if necessary to reduce spectral resolution
if (fft) I = I.fft(wn_low, wn_high, elwn, udr); //if an FFT is requested, perform it
if (background) I.absorbance(background); //if the user provides a background file
if (header_flag) {
header = I.create_header();
header_flag = false;
}
nbo = I.dim(2); //get the number of bands in the output image
for (size_t b = 0; b < nbo; b++) { //for each band in the input
for (size_t yi = 0; yi < ds; yi++) { //for each line in the input image
yo = ty * ds + ds - yi - 1; //calculate the y output position
io = (b * nxo * nyo + yo * nxo + xo); //calculate the output position as a 1d index
ii = (b * ds * ds + yi * ds); //calculate the input position as a 1d index
out.seekp(io * sizeof(float)); //seek to the correct output location
out.write(I.data() + ii * sizeof(float), ds * sizeof(float)); //copy the line from input to output
}
}
progress = (double)(ty * mx + tx + 1) / (double)(mx * my) * 100;
}
}
out.close();
t1.join(); //wait for the progress bar thread to finish (it probably already is)
header.samples = mx * ds;
header.lines = my * ds;
header.save(outfile + ".hdr");
}
void output_files(std::vector<stim::filename> infiles, std::string outfile) {
if (infiles.size() == 1) {
std::cout << "source file: " << infiles[0].str() << std::endl;
}
else {
std::cout << "number of input files: " << infiles.size() << std::endl;
}
std::cout << "destination file: " << outfile << std::endl;
}
int main(int argc, char* argv[]){
args.add("help", "print help");
args.add("spacing", "specify the spectral resolution", "", "positive integer specifying wavenumber spacing");
args.add("fft", "perform a fourier transform of the data before assembling the mosaic", "", "optional wavenumber range (override UDR), ex: 900 3500");
args.add("udr", "udr filter used for imaging", "4", "(4 or 8 supported)");
args.add("elwn", "ELWN tuning number", "15798.0039", "HeNe sampling rate (in cm^(-1))");
args.add("background", "specify a background file for measuring absorbance", "", "filename.seq");
args.add("cuda", "specify a CUDA device for the FFT", "0", "index of CUDA device (-1 forces CPU only)");
args.parse(argc, argv); //parse the argument list
if (args.nargs() < 2 || args["help"]) {
std::cout << "Usage: " << std::endl<<std::endl;
std::cout << "\tGenerate an ENVI mosaic from a set of tiles produced by a Cary FTIR imaging system:" << std::endl;
std::cout << "\t\t>> cary-ftir input_files/*.dmd output_mosaic" << std::endl << std::endl;
std::cout << args.str() << std::endl;
exit(1);
}
stim::filename infile(args.arg(0));
std::vector<stim::filename> in; //create a list of input files
if (args.nargs() == 2) {
if (infile.wildcards()) in = infile.get_list(); //if a wildcard is specified, get the file list
else in.push_back(infile); //if a single file is specified, create a file list of 1
}
//if the user specifies several files on the command line
else {
for (size_t f = 0; f < args.nargs() - 1; f++) { //for each file specified
in.push_back(args.arg(f)); //add the file to the file list
}
}
std::string outfilename = args.arg(args.nargs() - 1); //get the output file
output_files(in, outfilename);
udr = args["udr"].as_int(); //save the UDR filter type
if (udr == 4) {
wn_low = 900;
wn_high = 3900;
}
else {
wn_low = 800;
wn_high = 1900;
}
if (args["fft"].is_set()) {
fft = true; //set the fft flag to true
if (args["fft"].nargs() >= 1) {
if(args["fft"].is_num(0))
wn_low = args["fft"].as_float(0);
else {
std::cout << "ERROR: FFT output range must be a number: " << args["fft"].as_string(0) << std::endl;
return 1;
}
}
if (args["fft"].nargs() >= 2) {
if (args["fft"].is_num(1))
wn_high = args["fft"].as_float(1);
else {
std::cout << "ERROR: FFT output range must be a number: " << args["fft"].as_string(1) << std::endl;
return 1;
}
}
}
elwn = (float)args["elwn"].as_float();
if (args["spacing"]) {
double di = (1.0 / elwn) * ((double)udr / 2.0); //calculate the interferogram spacing
double df = (double)args["spacing"].as_int(); //get the desired spectral spacing
n_samples = (size_t)std::ceil(1.0 / (di * df)); //calculate the number of samples required
std::cout << "Override default spectral spacing with " << df << " cm^(-1) by cropping to " << n_samples << " samples" << std::endl;
crop_samples = true;
}
if (args["background"]) {
if (!args["fft"]) {
std::cout << "ERROR: using a background requires calculating an FFT - use the --fft option as well" << std::endl;
exit(1);
}
stim::agilent_binary<float> background(args["background"].as_string());
if (background) {
if (crop_samples) background.crop(n_samples); //crop the samples if necessary
background = background.fft(wn_low, wn_high, elwn, udr);
save_mosaic(in, outfilename, &background);
//absorbance = true;
}
else {
std::cout << "ERROR loading background file: " << args["background"].as_string() << std::endl;
return 1;
}
}
else {
save_mosaic(in, outfilename);
}
// _CrtDumpMemoryLeaks();
}
|
