metric.cpp
5.62 KB
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
#include <iostream>
#include <fstream>
#include <time.h>
#include <thread>
#include <iterator>
#include "stim/envi/envi_header.h"
#include "stim/envi/bip.h"
#include "stim/envi/bil.h"
#include "stim/envi/bsq.h"
#include "stim/envi/envi.h"
#include "stim/envi/binary.h"
void progress_thread_double(double* e);
extern stim::envi ENVI;
extern unsigned char* MASK;
/// Reads a file of numbers separated by lines and white space into a 2D STL array
std::vector<std::string> tokenize(std::string str){
// construct a stream from the string
std::stringstream strstr(str);
// use stream iterators to copy the stream to the vector as whitespace separated strings
std::istream_iterator<std::string> it(strstr);
std::istream_iterator<std::string> end;
std::vector<std::string> results(it, end);
// send the vector to stdout.
//std::ostream_iterator<std::string> oit(std::cout);
//std::copy(results.begin(), results.end(), oit);
return results;
}
/// Convert a metric file into an array of metric parameters
std::vector< std::vector<double> > read_metric_list(std::string filename){
std::ifstream infile(filename);
if(!infile){
std::cout<<"ERROR opening metric file: "<<filename<<std::endl;
exit(1);
}
std::vector< std::vector<double> > metrics;
std::string line; //stores the metric parameters as a series of values in text
while(std::getline(infile, line)){ //for each line in the input file
std::vector<std::string> tokens = tokenize(line); //break the line into tokens
std::vector<double> metric_params; //create an array to store the metric parameters
for(size_t t = 0; t < tokens.size(); t++){ //for each token
metric_params.push_back(atof(tokens[t].c_str()));
}
metrics.push_back(metric_params);
}
return metrics; //return the list of metrics and parameters
}
/// Creates an ENVI file filled with the given metrics, based on the Bhargava format
void create_metric_file(std::vector< std::vector<double> > metrics, std::string outfile){
//warn if the data set is BIP or BIL
if (ENVI.header.interleave == stim::envi_header::BIL ||
ENVI.header.interleave == stim::envi_header::BIP) {
std::cout << "WARNING: metric processing is slow for BIL and BIP files - recommend changing to BSQ" << std::endl;
}
double progress = 0;
std::thread t1(progress_thread_double, &progress); //start the progress bar thread
std::ofstream target(outfile.c_str(), std::ios::binary); //create an output file for writing the metrics
size_t XYbytes = ENVI.header.samples * ENVI.header.lines * ENVI.header.valsize(); //size of a band in bytes
void* M = malloc(XYbytes); //allocate space to store the metric result
stim::envi_header new_header = ENVI.header; //copy the current header
new_header.bands = metrics.size(); //set the number of bands to the number of metrics
new_header.interleave = stim::envi_header::BSQ; //the metric output file will be BSQ
new_header.band_names.resize(metrics.size()); //re-set the band names (names will be filled in below)
new_header.wavelength.clear(); //clear the wavelengths - the output won't have wavelengths
bool success = false;
for(size_t m = 0; m < metrics.size(); m++){ //for each metric
std::stringstream name_str;
switch( (unsigned)metrics[m][0] ){
case 1: //peak height ratio
if (ENVI.ph_to_ph(M,
metrics[m][1],
metrics[m][2],
metrics[m][3],
metrics[m][4],
metrics[m][5],
metrics[m][6],
MASK)){
success = true;
name_str << metrics[m][3] << " : " << metrics[m][6] << "\n";
new_header.band_names[m] = name_str.str();
}
else success = false;
break;
case 2: //peak area to height ratio
if (ENVI.pa_to_ph(M,
metrics[m][1],
metrics[m][2],
metrics[m][3],
metrics[m][4],
metrics[m][5],
metrics[m][6],
metrics[m][7],
MASK)) {
success = true;
name_str << metrics[m][3] << " : [" << metrics[m][6] << ", " << metrics[m][7] << "]\n";
new_header.band_names[m] = name_str.str();
}
else success = false;
break;
case 3: //peak area to area ratio
if (ENVI.pa_to_pa(M,
metrics[m][1],
metrics[m][2],
metrics[m][3],
metrics[m][4],
metrics[m][5],
metrics[m][6],
metrics[m][7],
metrics[m][8],
MASK)) {
success = true;
name_str << "[" << metrics[m][3] << ", " << metrics[m][4] << "] : [" << metrics[m][6] << ", " << metrics[m][7] << "]\n";
new_header.band_names[m] = name_str.str();
}
else success = false;
break;
case 4: //center of gravity
if (ENVI.centroid(M,
metrics[m][1],
metrics[m][2],
metrics[m][3],
metrics[m][4],
MASK)) {
success = true;
name_str << "[" << metrics[m][3] << ", " << metrics[m][4] << "]\n";
new_header.band_names[m] = name_str.str();
}
else success = false;
break;
default:
std::cout << "ERROR loading metric " << m << " from metric file" << std::endl;
exit(1);
}
if (!success) {
progress = 100;
t1.join();
progress = ((double)(m + 1) / (double)metrics.size()) * 100; //save the progress
std::cout << std::endl << "ERROR processing metric " << m + 1 << ":" << std::endl;
std::cout << metrics[m][0] << " "
<< metrics[m][1] << " "
<< metrics[m][2] << " "
<< metrics[m][3] << " "
<< metrics[m][4] << " "
<< metrics[m][5] << " "
<< metrics[m][6] << " "
<< metrics[m][6] << " "
<< metrics[m][6] << std::endl;
exit(1);
}
target.write((const char*)M, XYbytes); //write the metric results to a file
progress = ( (double)(m+1) / (double)metrics.size() ) * 100; //save the progress
}
t1.join(); //end the progress bar thread
new_header.save(outfile + ".hdr");
}