Blame view

src/proc/metric.cpp 5.62 KB
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
  #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");
  
  }