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
|
#include <iostream>
#include <algorithm>
#include <unordered_map>
#include <stim/parser/arguments.h>
#include <stim/envi/agilent_binary.h>
#include <stim/parser/filename.h>
#include <stim/envi/envi.h>
stim::arglist args; //generate a list of accepted arguments
#include <thread>
void progress_thread_double(double* e); //progress bar threaded function
struct speroFile {
int x;
int y;
stim::filename headername;
};
void spero_coord(int &x, int &y, stim::filename filename) {
//get the file prefix
std::string prefix = filename.get_prefix();
//these files are in the format:
// ???????????[xxxxx_yyyyy][??????????].hdr
//first find the x coordinate
//find the position of the first bracket
size_t start = prefix.find_first_of("[") + 1; //find the first element of the y coordinate
prefix = prefix.substr(start, prefix.length() - start);
size_t end = prefix.find_first_of("_"); //find the index of the last element of the y coordinate
std::string x_string = prefix.substr(0, end);
x = atoi(x_string.c_str());
//crop out the x coordinate and the last underscore
prefix = prefix.substr(end + 1, prefix.length() - end);
end = prefix.find_first_of("]");
std::string y_string = prefix.substr(0, end);
y = atoi(y_string.c_str());
//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());*/
}
bool sperosort(speroFile i, speroFile j) {
return i.x < j.x;
}
//This function test to make sure that the given file names are valid for a SPERO reconstruction
bool verify_filenames(std::vector<stim::filename> files) {
size_t idx = files[0].prefix().find('[');
if (idx == std::string::npos) return false;
for (size_t i = 0; i < files.size() - 1; i++) {
if (files[i].prefix().substr(0, idx) != files[i + 1].prefix().substr(0, idx)) return false;
}
return true;
}
void mosaic_spero(std::string directory, std::string outfile) {
//test for wild-cards in the directory
if (directory.find('*') != std::string::npos || directory.find('?') != std::string::npos) {
std::cout << "ERROR: no wild-cards are allowed, just specify the source directory containing SPERO files" << std::endl;
exit(1);
}
std::cout << "Building a mosaic from SPERO frames..." << std::endl;
//generate a list of dmd files
stim::filename f = directory + "/*.hdr";
std::vector<stim::filename> hdr_list = f.get_list();
if (hdr_list.size() == 0) {
std::cout << "ERROR: no files found for assembly" << std::endl;
exit(1);
}
std::cout << hdr_list.size() << " files found" << std::endl;
if (!verify_filenames(hdr_list)) {
std::cout << "ERROR: incorrect file prefix found - all files should have the same prefix" << std::endl;
exit(1);
}
if (hdr_list.size() == 0) {
std::cout << "ERROR: no header files found in directory '" + directory << std::endl;
exit(1);
}
std::vector<speroFile> speroFiles; //create a list of spero files
speroFiles.resize(hdr_list.size()); //pre-allocate to save time (we know how many there are)
for (size_t i = 0; i < hdr_list.size(); i++) { //for each header name
spero_coord(speroFiles[i].x, speroFiles[i].y, hdr_list[i]);
speroFiles[i].headername = hdr_list[i].str();
}
//sort the SPERO files by x coordinate
std::sort(speroFiles.begin(), speroFiles.end(), sperosort);
//calculate the height of the image by looking at how many files have the same x coordinate
size_t height = 0; //stores the width of the array
int first_x; //the first x value encountered
for (size_t i = 0; i < speroFiles.size(); i++) { //for each file
if (i == 0) first_x = speroFiles[i].x; //if this is the first file, store the x coordinate
if (first_x == speroFiles[i].x) height++;
else break;
}
size_t width;
if (height == 0) width = speroFiles.size();
else width = speroFiles.size() / height; //the number of files should be divisible by the width to get the height
if (width * height != speroFiles.size()) {
std::cout << "ERROR: calculated width and height do not match the number of input files: " << width << " x " << height << " = " << width * height << " found, " << speroFiles.size() << " needed" << std::endl;
exit(1);
}
std::vector<int> yvals(height); //stores the height values
for (size_t i = 0; i < height; i++) //go through the first "width" files and store their y values
yvals[i] = speroFiles[i].y;
std::vector<int> xvals(width); //stores the width values
for (size_t i = 0; i < speroFiles.size(); i += height) //go through every "width" file and store the x coordinate
xvals[i / height] = speroFiles[i].x;
std::sort(&yvals[0], &yvals[0] + yvals.size()); //sort both arrays of coordinates
std::sort(&xvals[0], &xvals[0] + xvals.size());
//create a mapping from the x and y coordinates to their indices
std::unordered_map<int, size_t> x2idx;
std::unordered_map<int, size_t> y2idx;
for (size_t i = 0; i < height; i++) y2idx.insert(std::pair<int, size_t>(yvals[i], i)); //unordered maps now give an index for a coordinate
for (size_t i = 0; i < width; i++) x2idx.insert(std::pair<int, size_t>(xvals[i], i));
std::vector<stim::envi> M; //create an array of ENVI files in mosaic order
M.resize(speroFiles.size());
size_t xi, yi;
for (size_t i = 0; i < M.size(); i++) { //for each input file
xi = x2idx[speroFiles[i].x]; //get the x and y indices for this file's position in the mosaic
yi = y2idx[speroFiles[i].y];
size_t idx = yi * width + xi;
M[idx].open(speroFiles[i].headername.str_noext(), speroFiles[i].headername.str()); //open the ENVI file and set parameters
M[idx].close(); //close the file, since we can only open so many
}
std::ofstream target(outfile.c_str(), std::ios::binary); //create a file to store the output
stim::envi_header h = M[0].header; //create a header for the mosaic
h.samples *= width;
h.lines *= height;
size_t tile_width = M[0].header.samples;
size_t tile_height = M[0].header.lines;
double progress = 0;
std::thread t1(progress_thread_double, &progress); //start the progress bar thread
float* band = (float*)malloc(h.samples * h.lines * sizeof(float)); //allocate space for a band image
float* tile_band = (float*)malloc(tile_width * tile_height * sizeof(float));
size_t dst, src;
//reconstruct the array
for (size_t b = 0; b < h.bands; b++) { //for each band
for (size_t yi = 0; yi < height; yi++) { //for each tile in X and Y
for (size_t xi = 0; xi < width; xi++) {
M[yi * width + xi].open(); //open the tile
M[yi * width + xi].band_index(tile_band, b); //fill the tile_band pointer with the tile image at band b
for (size_t y = 0; y < tile_height; y++) { //copy the tile image to the destination band image
//dst = (yi * tile_height + tile_height - y - 1) * h.samples + xi * tile_width;
dst = ((height - yi - 1) * tile_height + y) * h.samples + xi * tile_width;
src = y * tile_width;
memcpy(band + dst, tile_band + src, sizeof(float) * tile_width);
}
M[yi * width + xi].close(); //close the tile
progress = (double)(b * height * width + yi * width + xi + 1) / (width * height * h.bands) * 100; //update the progress bar
}
}
target.write((char*)band, h.samples * h.lines * sizeof(float)); //write the entire band to the output image
}
t1.join(); //wait for the thread to terminate
target.close(); //close the output file
h.save(outfile + ".hdr"); //write the header
free(band);
}
void advertise() {
std::cout << std::endl << std::endl;
std::cout << "=========================================================================" << std::endl;
std::cout << "Thank you for using the SIPROC spectroscopic image processing toolkit!" << std::endl;
std::cout << "Scalable Tissue Imaging and Modeling (STIM) Lab, University of Houston" << std::endl;
std::cout << "Developers: Sam Saki, Ziqi He, David Mayerich, Brad Deutsch" << std::endl;
std::cout << "Source: https://github.com/stimlab/hsiproc.git" << std::endl;
std::cout << "This version has been compiled on " << __DATE__ << " at " << __TIME__ << std::endl;
std::cout << "=========================================================================" << std::endl << std::endl;
std::cout << "Daylight SPERO mosaic assembly usage:" << std::endl;
std::cout << " >> spero directory/to/files mosaic" << std::endl;
}
int main(int argc, char* argv[]) {
args.add("help", "print help");
args.parse(argc, argv); //parse the argument list
if (args["help"]) {
std::cout << args.str() << std::endl;
exit(1);
}
std::string indir = ".";
std::string outfile = "mosaic";
if (args.nargs() == 1) outfile = args.arg(0); //if only one argument is given, assume it is the output file
else if (args.nargs() >= 2) { //if two arguments are given
indir = args.arg(0); //the first argument is the directory
outfile = args.arg(1); //the second argument is the output file
}
std::cout << "SPERO mosaic assembly: " << indir << " --> " << outfile << std::endl;
mosaic_spero(indir, outfile); //call the SPERO mosaic function
}
|
