Blame view

stim/envi/bsq.h 35.1 KB
e8eb202f   David Mayerich   added a new ENVI ...
1
2
3
4
  #ifndef STIM_BSQ_H
  #define STIM_BSQ_H
  
  #include "../envi/envi_header.h"
3fc1d461   heziqi   Ziqi added saveBa...
5
  #include "../envi/binary.h"
e8eb202f   David Mayerich   added a new ENVI ...
6
  #include "../envi/bil.h"
fc712213   David Mayerich   code simplifications
7
8
  #include <cstring>
  #include <utility>
e8eb202f   David Mayerich   added a new ENVI ...
9
  #include <vector>
3fc1d461   heziqi   Ziqi added saveBa...
10
  
77060b79   David Mayerich   started DOxygen d...
11
12
  
  
8a86bd56   David Mayerich   changed rts names...
13
  namespace stim{
3fc1d461   heziqi   Ziqi added saveBa...
14
  
77060b79   David Mayerich   started DOxygen d...
15
16
17
18
19
20
21
22
  /**
  	The BIP class represents a 3-dimensional binary file stored using band sequential (BSQ) image encoding. The binary file is stored
  	such that X-Y "frames" are stored sequentially to form an image stack along the z-axis. When accessing the data sequentially on disk,
  	the dimensions read, from fastest to slowest, are X, Y, Z.
  
  	This class is optimized for data streaming, and therefore supports extremely large (terabyte-scale) files. Data is loaded from disk
  	on request. Functions used to access data are written to support efficient reading.
  */
3fc1d461   heziqi   Ziqi added saveBa...
23
24
25
26
  template <typename T>
  
  class bsq: public binary<T> {
  
77060b79   David Mayerich   started DOxygen d...
27
  
3fc1d461   heziqi   Ziqi added saveBa...
28
  protected:
43dec788   David Mayerich   added code for si...
29
  
e8eb202f   David Mayerich   added a new ENVI ...
30
  	std::vector<double> w;	//band wavelengths
6708cc25   heziqi   Ziqi added envi c...
31
  	unsigned int offset;
3fc1d461   heziqi   Ziqi added saveBa...
32
  
724ec347   David Mayerich   simplified the EN...
33
34
35
36
37
38
39
40
41
42
43
44
  	using binary<T>::R;
  
  	unsigned long X(){
  		return R[0];
  	}
  	unsigned long Y(){
  		return R[1];
  	}
  	unsigned long Z(){
  		return R[2];
  	}
  
cbce396e   David Mayerich   added support for...
45
46
  	using binary<T>::thread_data;
  
3fc1d461   heziqi   Ziqi added saveBa...
47
48
  public:
  
fc712213   David Mayerich   code simplifications
49
50
  	using binary<T>::open;
  	using binary<T>::file;
724ec347   David Mayerich   simplified the EN...
51
52
  	using binary<T>::read_line_01;
  	using binary<T>::read_plane_2;
32504c2c   David Mayerich   fixed some functi...
53
  	//using binary<T>::getSlice;
43dec788   David Mayerich   added code for si...
54
  
fc712213   David Mayerich   code simplifications
55
  
a08c98c3   David Mayerich   finished document...
56
  	/// Open a data file for reading using the class interface.
77060b79   David Mayerich   started DOxygen d...
57
58
59
60
61
62
  
  	/// @param filename is the name of the binary file on disk
  	/// @param X is the number of samples along dimension 1
  	/// @param Y is the number of samples (lines) along dimension 2
  	/// @param B is the number of samples (bands) along dimension 3
  	/// @param header_offset is the number of bytes (if any) in the binary header
43dec788   David Mayerich   added code for si...
63
  	/// @param wavelengths is an STL vector of size B specifying a numerical label for each band
e8eb202f   David Mayerich   added a new ENVI ...
64
  	bool open(std::string filename, unsigned int X, unsigned int Y, unsigned int B, unsigned int header_offset, std::vector<double> wavelengths){
3fc1d461   heziqi   Ziqi added saveBa...
65
  
e8eb202f   David Mayerich   added a new ENVI ...
66
67
  		//copy the wavelengths to the BSQ file structure
  		w = wavelengths;
6708cc25   heziqi   Ziqi added envi c...
68
69
  		//copy the wavelengths to the structure
  		offset = header_offset;
e8eb202f   David Mayerich   added a new ENVI ...
70
71
  
  		return open(filename, vec<unsigned int>(X, Y, B), header_offset);
3fc1d461   heziqi   Ziqi added saveBa...
72
73
  	}
  
a08c98c3   David Mayerich   finished document...
74
  	/// Retrieve a single band (based on index) and stores it in pre-allocated memory.
77060b79   David Mayerich   started DOxygen d...
75
76
77
  
  	/// @param p is a pointer to an allocated region of memory at least X * Y * sizeof(T) in size.
  	/// @param page <= B is the integer number of the band to be copied.
0e48cc9c   heziqi   Ziqi added getBan...
78
  	bool band_index( T * p, unsigned int page){
724ec347   David Mayerich   simplified the EN...
79
  		return read_plane_2(p, page);
3fc1d461   heziqi   Ziqi added saveBa...
80
81
  	}
  
a08c98c3   David Mayerich   finished document...
82
  	/// Retrieve a single band (by numerical label) and stores it in pre-allocated memory.
77060b79   David Mayerich   started DOxygen d...
83
84
85
  
  	/// @param p is a pointer to an allocated region of memory at least X * Y * sizeof(T) in size.
  	/// @param wavelength is a floating point value (usually a wavelength in spectral data) used as a label for the band to be copied.
e933c14e   David Mayerich   cleaned up the code
86
  	bool band( T * p, double wavelength){
0e48cc9c   heziqi   Ziqi added getBan...
87
  
ee4dea28   David Mayerich   fixed errors in c...
88
89
90
  		//if there are no wavelengths in the BSQ file
  		if(w.size() == 0)
  			return band_index(p, (unsigned int)wavelength);
0e48cc9c   heziqi   Ziqi added getBan...
91
  
724ec347   David Mayerich   simplified the EN...
92
  		unsigned int XY = X() * Y();	//calculate the number of pixels in a band
43dec788   David Mayerich   added code for si...
93
  		unsigned page = 0;
517876d6   heziqi   metrics finished ...
94
  
0e48cc9c   heziqi   Ziqi added getBan...
95
  
ebe94c9d   David Mayerich   modified comments
96
  		//get the two neighboring bands (above and below 'wavelength')
0e48cc9c   heziqi   Ziqi added getBan...
97
98
  
  		//if wavelength is smaller than the first one in header file
e8eb202f   David Mayerich   added a new ENVI ...
99
  		if ( w[page] > wavelength ){
0e48cc9c   heziqi   Ziqi added getBan...
100
101
102
103
  			band_index(p, page);
  			return true;
  		}
  
e8eb202f   David Mayerich   added a new ENVI ...
104
  		while( w[page] < wavelength )
0e48cc9c   heziqi   Ziqi added getBan...
105
106
  		{
  			page++;
ebe94c9d   David Mayerich   modified comments
107
108
  			//if wavelength is larger than the last wavelength in the header file
  			//	(the wavelength is out of bounds)
724ec347   David Mayerich   simplified the EN...
109
110
  			if (page == Z()) {
  				band_index(p, Z()-1);		//return the last band
0e48cc9c   heziqi   Ziqi added getBan...
111
112
113
  				return true;
  			}
  		}
ebe94c9d   David Mayerich   modified comments
114
115
  		//when the page counter points to the first band above 'wavelength'
  		if ( wavelength < w[page] ){
193ff00e   David Mayerich   dealth with preci...
116
  
ebe94c9d   David Mayerich   modified comments
117
  			//do the interpolation
517876d6   heziqi   metrics finished ...
118
119
  			T * p1;
  			T * p2;
0e48cc9c   heziqi   Ziqi added getBan...
120
121
122
123
124
  			p1=(T*)malloc( XY * sizeof(T));                     //memory allocation
  			p2=(T*)malloc( XY * sizeof(T));
  			band_index(p1, page - 1);
  			band_index(p2, page );
  			for(unsigned i=0; i < XY; i++){
e8eb202f   David Mayerich   added a new ENVI ...
125
  				double r = (double) (wavelength - w[page-1]) / (double) (w[page] - w[page-1]);
0e48cc9c   heziqi   Ziqi added getBan...
126
127
  				p[i] = (p2[i] - p1[i]) * r + p1[i];
  			}
517876d6   heziqi   metrics finished ...
128
129
  			free(p1);
  			free(p2);
0e48cc9c   heziqi   Ziqi added getBan...
130
  		}
ebe94c9d   David Mayerich   modified comments
131
132
133
  		//if the wavelength is equal to a wavelength in header file
  		else{
  			band_index(p, page);		//return the band
0e48cc9c   heziqi   Ziqi added getBan...
134
135
  		}
  
0e48cc9c   heziqi   Ziqi added getBan...
136
137
138
  		return true;
  	}
  
724ec347   David Mayerich   simplified the EN...
139
  	/// Retrieve a single spectrum (Z-axis line) at a given (x, y) location and stores it in pre-allocated memory.
77060b79   David Mayerich   started DOxygen d...
140
141
142
143
  
  	/// @param p is a pointer to pre-allocated memory at least B * sizeof(T) in size.
  	/// @param x is the x-coordinate (dimension 1) of the spectrum.
  	/// @param y is the y-coordinate (dimension 2) of the spectrum.
e933c14e   David Mayerich   cleaned up the code
144
  	bool spectrum(T * p, unsigned x, unsigned y){
724ec347   David Mayerich   simplified the EN...
145
  		return read_line_01(p, x, y);
3fc1d461   heziqi   Ziqi added saveBa...
146
  	}
f0d5a769   heziqi   Ziqi changed the ...
147
  
a08c98c3   David Mayerich   finished document...
148
  	/// Retrieve a single pixel and stores it in pre-allocated memory.
77060b79   David Mayerich   started DOxygen d...
149
150
151
  
  	/// @param p is a pointer to pre-allocated memory at least sizeof(T) in size.
  	/// @param n is an integer index to the pixel using linear array indexing.
bfe9c6db   heziqi   added feature_mat...
152
153
  	bool pixel(T * p, unsigned n){
  
724ec347   David Mayerich   simplified the EN...
154
  		unsigned bandnum = X() * Y();		//calculate numbers in one band
bfe9c6db   heziqi   added feature_mat...
155
156
157
158
159
  		if ( n >= bandnum){							//make sure the pixel number is right
  			std::cout<<"ERROR: sample or line out of range"<<std::endl;
  			return false;
  		}
  
eb5dfb2b   David Mayerich   fixed linux compa...
160
  		file.seekg(n * sizeof(T) + binary<T>::header, std::ios::beg);           //point to the certain pixel
724ec347   David Mayerich   simplified the EN...
161
  		for (unsigned i = 0; i < Z(); i++)
bfe9c6db   heziqi   added feature_mat...
162
163
164
165
166
  		{
  			file.read((char *)(p + i), sizeof(T));
  			file.seekg((bandnum - 1) * sizeof(T), std::ios::cur);    //go to the next band
  		}
  
43dec788   David Mayerich   added code for si...
167
  		return true;
bfe9c6db   heziqi   added feature_mat...
168
169
  	}
  
a08c98c3   David Mayerich   finished document...
170
  	/// Perform baseline correction given a list of baseline points and stores the result in a new BSQ file.
77060b79   David Mayerich   started DOxygen d...
171
172
  
  	/// @param outname is the name of the output file used to store the resulting baseline-corrected data.
43dec788   David Mayerich   added code for si...
173
  	/// @param wls is the list of baseline points based on band labels.
f0d5a769   heziqi   Ziqi changed the ...
174
175
176
  	bool baseline(std::string outname, std::vector<double> wls )
  	{
  		unsigned N = wls.size();			//get the number of baseline points
31394e2a   heziqi   Ziqi added method...
177
  
f0d5a769   heziqi   Ziqi changed the ...
178
  		std::ofstream target(outname.c_str(), std::ios::binary);	//open the target binary file
0c748c81   heziqi   Ziqi added save h...
179
  		std::string headername = outname + ".hdr";              //the header file name
f0d5a769   heziqi   Ziqi changed the ...
180
181
  
  		//simplify image resolution
724ec347   David Mayerich   simplified the EN...
182
183
  		unsigned int B = Z();		//calculate the number of bands
  		unsigned int XY = X() * Y();	//calculate the number of pixels in a band
f0d5a769   heziqi   Ziqi changed the ...
184
185
186
187
  		unsigned int S = XY * sizeof(T);		//calculate the number of bytes in a band
  
  		double ai, bi;	//stores the two baseline points wavelength surrounding the current band
  		double ci;		//stores the current band's wavelength
77060b79   David Mayerich   started DOxygen d...
188
  
f0d5a769   heziqi   Ziqi changed the ...
189
190
191
192
193
194
195
  		unsigned control=0;
  
  		T * a;					//pointers to the high and low band images
  		T * b;
  		T * c;				//pointer to the current image
  
  		a = (T*)malloc( S );     //memory allocation
43dec788   David Mayerich   added code for si...
196
197
  		b = (T*)malloc( S );
  		c = (T*)malloc( S );
f0d5a769   heziqi   Ziqi changed the ...
198
199
200
201
202
203
204
  
  		if (a == NULL || b == NULL || c == NULL){
  			std::cout<<"ERROR: error allocating memory";
  			exit(1);
  		}
  
  
43dec788   David Mayerich   added code for si...
205
  		//initialize lownum, highnum, low, high
e8eb202f   David Mayerich   added a new ENVI ...
206
  		ai=w[0];
f0d5a769   heziqi   Ziqi changed the ...
207
208
209
  
  		//if no baseline point is specified at band 0,
  			//set the baseline point at band 0 to 0
e8eb202f   David Mayerich   added a new ENVI ...
210
  		if(wls[0] != w[0]){
43dec788   David Mayerich   added code for si...
211
  			bi = wls[control];
f0d5a769   heziqi   Ziqi changed the ...
212
213
214
215
216
  			memset(a, (char)0, S);
  		}
  		//else get the low band
  		else{
  			control += 1;
e933c14e   David Mayerich   cleaned up the code
217
  			band(a, ai);
f0d5a769   heziqi   Ziqi changed the ...
218
219
220
  			bi = wls[control];
  		}
  		//get the high band
e933c14e   David Mayerich   cleaned up the code
221
  		band(b, bi);
f0d5a769   heziqi   Ziqi changed the ...
222
  
43dec788   David Mayerich   added code for si...
223
  		//correct every band
f0d5a769   heziqi   Ziqi changed the ...
224
225
226
  		for(unsigned cii = 0; cii < B; cii++){
  
  			//update baseline points, if necessary
e8eb202f   David Mayerich   added a new ENVI ...
227
  			if( w[cii] >= bi && cii != B - 1) {
f0d5a769   heziqi   Ziqi changed the ...
228
229
230
231
232
233
234
235
236
  				//if the high band is now on the last BL point?
  				if (control != N-1) {
  
  					control++;		//increment the index
  
  					std::swap(a, b);	//swap the baseline band pointers
  
  					ai = bi;
  					bi = wls[control];
e933c14e   David Mayerich   cleaned up the code
237
  					band(b, bi);
f0d5a769   heziqi   Ziqi changed the ...
238
239
240
  
  				}
  				//if the last BL point on the last band of the file?
e8eb202f   David Mayerich   added a new ENVI ...
241
  				else if ( wls[control] < w[B - 1]) {
f0d5a769   heziqi   Ziqi changed the ...
242
243
244
245
246
247
  
  					std::swap(a, b);	//swap the baseline band pointers
  
  					memset(b, (char)0, S);	//clear the high band
  
  					ai = bi;
e8eb202f   David Mayerich   added a new ENVI ...
248
  					bi = w[B - 1];
f0d5a769   heziqi   Ziqi changed the ...
249
250
251
252
253
  				}
  			}
  
  			//get the current band
  			band_index(c, cii);
e8eb202f   David Mayerich   added a new ENVI ...
254
  			ci = w[cii];
43dec788   David Mayerich   added code for si...
255
  
f0d5a769   heziqi   Ziqi changed the ...
256
257
258
  			//perform the baseline correction
  			for(unsigned i=0; i < XY; i++){
  				double r = (double) (ci - ai) / (double) (bi - ai);
92e4cf05   heziqi   include file fixed
259
  				c[i] =(T) ( c[i] - (b[i] - a[i]) * r - a[i] );
f0d5a769   heziqi   Ziqi changed the ...
260
  			}
43dec788   David Mayerich   added code for si...
261
  
f0d5a769   heziqi   Ziqi changed the ...
262
  			target.write(reinterpret_cast<const char*>(c), S);   //write the corrected data into destination
cbce396e   David Mayerich   added support for...
263
264
  
  			thread_data = (double)cii / B * 100;
43dec788   David Mayerich   added code for si...
265
266
  
  		}
0c748c81   heziqi   Ziqi added save h...
267
  
e8eb202f   David Mayerich   added a new ENVI ...
268
  		//header.save(headername);         //save the new header file
43dec788   David Mayerich   added code for si...
269
  
f0d5a769   heziqi   Ziqi changed the ...
270
271
272
273
  		free(a);
  		free(b);
  		free(c);
  		target.close();
cbce396e   David Mayerich   added support for...
274
275
  
  		thread_data = 100;
f0d5a769   heziqi   Ziqi changed the ...
276
277
  		return true;
  	}
089d9c81   heziqi   Ziqi added normal...
278
  
a08c98c3   David Mayerich   finished document...
279
  	/// Normalize all spectra based on the value of a single band, storing the result in a new BSQ file.
77060b79   David Mayerich   started DOxygen d...
280
281
282
283
  
  	/// @param outname is the name of the output file used to store the resulting baseline-corrected data.
  	///	@param w is the label specifying the band that the hyperspectral image will be normalized to.
  	///	@param t is a threshold specified such that a spectrum with a value at w less than t is set to zero. Setting this threshold allows the user to limit division by extremely small numbers.
93503ec6   David Mayerich   allowed masking d...
284
  	bool normalize(std::string outname, double w, double t = 0.0)
089d9c81   heziqi   Ziqi added normal...
285
  	{
724ec347   David Mayerich   simplified the EN...
286
287
  		unsigned int B = Z();		//calculate the number of bands
  		unsigned int XY = X() * Y();	//calculate the number of pixels in a band
089d9c81   heziqi   Ziqi added normal...
288
289
290
291
292
293
294
295
296
  		unsigned int S = XY * sizeof(T);		//calculate the number of bytes in a band
  
  		std::ofstream target(outname.c_str(), std::ios::binary);	//open the target binary file
  		std::string headername = outname + ".hdr";              //the header file name
  
  		T * b;					//pointers to the certain wavelength band
  		T * c;				//pointer to the current image
  
  		b = (T*)malloc( S );     //memory allocation
43dec788   David Mayerich   added code for si...
297
  		c = (T*)malloc( S );
089d9c81   heziqi   Ziqi added normal...
298
  
e933c14e   David Mayerich   cleaned up the code
299
  		band(b, w);             //get the certain band into memory
089d9c81   heziqi   Ziqi added normal...
300
301
302
303
304
305
  
  		for(unsigned j = 0; j < B; j++)
  		{
  			band_index(c, j);                     //get the current band into memory
  			for(unsigned i = 0; i < XY; i++)
  			{
93503ec6   David Mayerich   allowed masking d...
306
307
308
309
  				if(b[i] < t)
  					c[i] = (T)0.0;
  				else
  					c[i] = c[i] / b[i];
089d9c81   heziqi   Ziqi added normal...
310
311
  			}
  			target.write(reinterpret_cast<const char*>(c), S);   //write normalized data into destination
cbce396e   David Mayerich   added support for...
312
313
  
  			thread_data = (double)j / B * 100;
089d9c81   heziqi   Ziqi added normal...
314
315
  		}
  
43dec788   David Mayerich   added code for si...
316
  
cbce396e   David Mayerich   added support for...
317
  
e8eb202f   David Mayerich   added a new ENVI ...
318
  		//header.save(headername);         //save the new header file
43dec788   David Mayerich   added code for si...
319
  
089d9c81   heziqi   Ziqi added normal...
320
321
322
  		free(b);
  		free(c);
  		target.close();
cbce396e   David Mayerich   added support for...
323
  		thread_data = 100;					//make sure that the progress bar is full
089d9c81   heziqi   Ziqi added normal...
324
325
  		return true;
  	}
43dec788   David Mayerich   added code for si...
326
  
a08c98c3   David Mayerich   finished document...
327
  	/// Convert the current BSQ file to a BIP file with the specified file name.
77060b79   David Mayerich   started DOxygen d...
328
329
  
  	/// @param outname is the name of the output BIP file to be saved to disk.
421219ba   heziqi   Ziqi added bip.h ...
330
331
  	bool bip(std::string outname)
  	{
f6169dea   heziqi   Ziqi completed bi...
332
333
334
335
336
  		std::string temp = outname + "_temp";
  		std::string headtemp = temp + ".hdr";
  		//first creat a temporary bil file and convert bsq file to bil file
  		bil(temp);
  
8a86bd56   David Mayerich   changed rts names...
337
  		stim::bil<T> n;
724ec347   David Mayerich   simplified the EN...
338
  		if(n.open(temp, X(), Y(), Z(), offset, w)==false){        //open infile
f6169dea   heziqi   Ziqi completed bi...
339
340
  			std::cout<<"ERROR: unable to open input file"<<std::endl;
  			exit(1);
421219ba   heziqi   Ziqi added bip.h ...
341
  		}
f6169dea   heziqi   Ziqi completed bi...
342
343
344
345
346
  		//then convert bil file to bip file
  		n.bip(outname);
  		n.close();
  		remove(temp.c_str());
  		remove(headtemp.c_str());
421219ba   heziqi   Ziqi added bip.h ...
347
348
  		return true;
  	}
c25e7d0d   heziqi   speed of bip.base...
349
  
a08c98c3   David Mayerich   finished document...
350
  	/// Convert the current BSQ file to a BIL file with the specified file name.
77060b79   David Mayerich   started DOxygen d...
351
352
  
  	/// @param outname is the name of the output BIL file to be saved to disk.
c25e7d0d   heziqi   speed of bip.base...
353
354
355
  	bool bil(std::string outname)
  	{
  		//simplify image resolution
724ec347   David Mayerich   simplified the EN...
356
357
  		unsigned int L = X() * Z() * sizeof(T);		//calculate the number of bytes of a ZX slice
  		unsigned int jump = (Y() - 1) * X() * sizeof(T);
43dec788   David Mayerich   added code for si...
358
  
c25e7d0d   heziqi   speed of bip.base...
359
360
  		std::ofstream target(outname.c_str(), std::ios::binary);
  		std::string headername = outname + ".hdr";
43dec788   David Mayerich   added code for si...
361
  
cbce396e   David Mayerich   added support for...
362
363
364
365
366
  		unsigned int xz_bytes = X() * Z() * sizeof(T);
  		T * xz_slice;						//pointer to the current spectrum
  		xz_slice = (T*)malloc(xz_bytes);
  
  		for ( unsigned y = 0; y < Y(); y++)									//for each y position
c25e7d0d   heziqi   speed of bip.base...
367
  		{
cbce396e   David Mayerich   added support for...
368
369
  			file.seekg(y * X() * sizeof(T), std::ios::beg);					//seek to the beginning of the xz slice
  			for ( unsigned z = 0; z < Z(); z++ )							//for each band
c25e7d0d   heziqi   speed of bip.base...
370
  			{
cbce396e   David Mayerich   added support for...
371
372
373
  				file.read((char *)(xz_slice + z * X()), sizeof(T) * X());	//read a line
  				file.seekg(jump, std::ios::cur);							//seek to the next band
  
43dec788   David Mayerich   added code for si...
374
  
3150c537   David Mayerich   added further pro...
375
  				thread_data = (double)(y * Z() + z) / (Z() * Y()) * 100;	//update the progress counter
43dec788   David Mayerich   added code for si...
376
377
  			}
  			target.write(reinterpret_cast<const char*>(xz_slice), xz_bytes);   //write the generated XZ slice to the target file
c25e7d0d   heziqi   speed of bip.base...
378
  		}
3150c537   David Mayerich   added further pro...
379
380
  
  
43dec788   David Mayerich   added code for si...
381
  
cbce396e   David Mayerich   added support for...
382
  		thread_data = 100;
43dec788   David Mayerich   added code for si...
383
  
cbce396e   David Mayerich   added support for...
384
  		free(xz_slice);
c25e7d0d   heziqi   speed of bip.base...
385
  		target.close();
3150c537   David Mayerich   added further pro...
386
  
c25e7d0d   heziqi   speed of bip.base...
387
388
  		return true;
  	}
20c212c0   heziqi   Ziqi added functi...
389
  
a08c98c3   David Mayerich   finished document...
390
  	/// Return a baseline corrected band given two adjacent baseline points and their bands. The result is stored in a pre-allocated array.
77060b79   David Mayerich   started DOxygen d...
391
392
393
394
395
396
397
  
  	/// @param lb is the label value for the left baseline point
  	/// @param rb is the label value for the right baseline point
  	/// @param lp is a pointer to an array holding the band image for the left baseline point
  	/// @param rp is a pointer to an array holding the band image for the right baseline point
  	/// @param wavelength is the label value for the requested baseline-corrected band
  	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size.
20c212c0   heziqi   Ziqi added functi...
398
399
  	bool baseline_band(double lb, double rb, T* lp, T* rp, double wavelength, T* result){
  
724ec347   David Mayerich   simplified the EN...
400
  		unsigned XY = X() * Y();
20c212c0   heziqi   Ziqi added functi...
401
402
403
404
405
406
407
408
409
  		band(result, wavelength);		//get band
  
  		//perform the baseline correction
  		double r = (double) (wavelength - lb) / (double) (rb - lb);
  		for(unsigned i=0; i < XY; i++){
  			result[i] =(T) (result[i] - (rp[i] - lp[i]) * r - lp[i] );
  		}
  		return true;
  	}
517876d6   heziqi   metrics finished ...
410
  
a08c98c3   David Mayerich   finished document...
411
  	/// Return a baseline corrected band given two adjacent baseline points. The result is stored in a pre-allocated array.
77060b79   David Mayerich   started DOxygen d...
412
413
414
415
416
  
  	/// @param lb is the label value for the left baseline point
  	/// @param rb is the label value for the right baseline point
  	/// @param bandwavelength is the label value for the desired baseline-corrected band
  	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size.
70407ea9   heziqi   Ziqi completed he...
417
  	bool height(double lb, double rb, double bandwavelength, T* result){
20c212c0   heziqi   Ziqi added functi...
418
419
420
  
  		T* lp;
  		T* rp;
724ec347   David Mayerich   simplified the EN...
421
  		unsigned XY = X() * Y();
70407ea9   heziqi   Ziqi completed he...
422
423
424
  		unsigned S = XY * sizeof(T);
  		lp = (T*) malloc(S);			//memory allocation
  		rp = (T*) malloc(S);
20c212c0   heziqi   Ziqi added functi...
425
  
70407ea9   heziqi   Ziqi completed he...
426
  		band(lp, lb);
43dec788   David Mayerich   added code for si...
427
  		band(rp, rb);
20c212c0   heziqi   Ziqi added functi...
428
429
430
  
  		baseline_band(lb, rb, lp, rp, bandwavelength, result);
  
517876d6   heziqi   metrics finished ...
431
432
  		free(lp);
  		free(rp);
20c212c0   heziqi   Ziqi added functi...
433
434
435
  		return true;
  	}
  
c359422d   heziqi   Ziqi added functi...
436
  
a08c98c3   David Mayerich   finished document...
437
  	/// Calculate the area under the spectrum between two specified points and stores the result in a pre-allocated array.
77060b79   David Mayerich   started DOxygen d...
438
439
440
441
442
443
  
  	/// @param lb is the label value for the left baseline point
  	/// @param rb is the label value for the right baseline point
  	/// @param lab is the label value for the left bound (start of the integration)
  	/// @param rab is the label value for the right bound (end of the integration)
  	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
70407ea9   heziqi   Ziqi completed he...
444
  	bool area(double lb, double rb, double lab, double rab, T* result){
20c212c0   heziqi   Ziqi added functi...
445
446
447
448
449
  
  		T* lp;	//left band pointer
  		T* rp;	//right band pointer
  		T* cur;		//current band 1
  		T* cur2;	//current band 2
20c212c0   heziqi   Ziqi added functi...
450
  
724ec347   David Mayerich   simplified the EN...
451
  		unsigned XY = X() * Y();
20c212c0   heziqi   Ziqi added functi...
452
453
454
455
456
457
  		unsigned S = XY * sizeof(T);
  
  		lp = (T*) malloc(S);			//memory allocation
  		rp = (T*) malloc(S);
  		cur = (T*) malloc(S);
  		cur2 = (T*) malloc(S);
20c212c0   heziqi   Ziqi added functi...
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
  
  		memset(result, (char)0, S);
  
  		//find the wavelenght position in the whole band
  		unsigned int n = w.size();
  		unsigned int ai = 0;		//left bound position
  		unsigned int bi = n - 1;		//right bound position
  
  
  
  		//to make sure the left and the right bound are in the bandwidth
  		if (lb < w[0] || rb < w[0] || lb > w[n-1] || rb >w[n-1]){
  			std::cout<<"ERROR: left bound or right bound out of bandwidth"<<std::endl;
  			exit(1);
  		}
  		//to make sure rigth bound is bigger than left bound
  		else if(lb > rb){
  			std::cout<<"ERROR: right bound should be bigger than left bound"<<std::endl;
  			exit(1);
  		}
  
  		//get the position of lb and rb
70407ea9   heziqi   Ziqi completed he...
480
  		while (lab >= w[ai]){
20c212c0   heziqi   Ziqi added functi...
481
482
  			ai++;
  		}
70407ea9   heziqi   Ziqi completed he...
483
  		while (rab <= w[bi]){
20c212c0   heziqi   Ziqi added functi...
484
485
486
487
488
489
490
  			bi--;
  		}
  
  		band(lp, lb);
  		band(rp, rb);
  
  		//calculate the beginning and the ending part
70407ea9   heziqi   Ziqi completed he...
491
492
  		baseline_band(lb, rb, lp, rp, rab, cur2);		//ending part
  		baseline_band(lb, rb, lp, rp, w[bi], cur);
20c212c0   heziqi   Ziqi added functi...
493
  		for(unsigned j = 0; j < XY; j++){
70407ea9   heziqi   Ziqi completed he...
494
  			result[j] += (rab - w[bi]) * (cur[j] + cur2[j]) / 2.0;
20c212c0   heziqi   Ziqi added functi...
495
  		}
70407ea9   heziqi   Ziqi completed he...
496
497
  		baseline_band(lb, rb, lp, rp, lab, cur2);		//beginnning part
  		baseline_band(lb, rb, lp, rp, w[ai], cur);
43dec788   David Mayerich   added code for si...
498
  		for(unsigned j = 0; j < XY; j++){
70407ea9   heziqi   Ziqi completed he...
499
  			result[j] += (w[ai] - lab) * (cur[j] + cur2[j]) / 2.0;
20c212c0   heziqi   Ziqi added functi...
500
  		}
70407ea9   heziqi   Ziqi completed he...
501
  
20c212c0   heziqi   Ziqi added functi...
502
503
504
505
506
507
508
  		//calculate the area
  		ai++;
  		for(unsigned i = ai; i <= bi ;i++)
  		{
  			baseline_band(lb, rb, lp, rp, w[ai], cur2);
  			for(unsigned j = 0; j < XY; j++)
  			{
43dec788   David Mayerich   added code for si...
509
  				result[j] += (w[ai] - w[ai-1]) * (cur[j] + cur2[j]) / 2.0;
20c212c0   heziqi   Ziqi added functi...
510
511
512
  			}
  			std::swap(cur,cur2);		//swap the band pointers
  		}
70407ea9   heziqi   Ziqi completed he...
513
  
517876d6   heziqi   metrics finished ...
514
515
516
517
  		free(lp);
  		free(rp);
  		free(cur);
  		free(cur2);
20c212c0   heziqi   Ziqi added functi...
518
519
520
  		return true;
  	}
  
a08c98c3   David Mayerich   finished document...
521
  	/// Compute the ratio of two baseline-corrected peaks. The result is stored in a pre-allocated array.
77060b79   David Mayerich   started DOxygen d...
522
523
524
525
526
527
528
529
  
  	/// @param lb1 is the label value for the left baseline point for the first peak (numerator)
  	/// @param rb1 is the label value for the right baseline point for the first peak (numerator)
  	/// @param pos1 is the label value for the first peak (numerator) position
  	/// @param lb2 is the label value for the left baseline point for the second peak (denominator)
  	/// @param rb2 is the label value for the right baseline point for the second peak (denominator)
  	/// @param pos2 is the label value for the second peak (denominator) position
  	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
70407ea9   heziqi   Ziqi completed he...
530
  	bool ph_to_ph(double lb1, double rb1, double pos1, double lb2, double rb2, double pos2, T * result){
20c212c0   heziqi   Ziqi added functi...
531
  
43dec788   David Mayerich   added code for si...
532
533
534
  		T* p1 = (T*)malloc(X() * Y() * sizeof(T));
  		T* p2 = (T*)malloc(X() * Y() * sizeof(T));
  
70407ea9   heziqi   Ziqi completed he...
535
536
537
538
  		//get the two peak band
  		height(lb1, rb1, pos1, p1);
  		height(lb2, rb2, pos2, p2);
  		//calculate the ratio in result
724ec347   David Mayerich   simplified the EN...
539
  		for(unsigned i = 0; i < X() * Y(); i++){
70407ea9   heziqi   Ziqi completed he...
540
541
542
543
544
  			if(p1[i] == 0 && p2[i] ==0)
  				result[i] = 1;
  			else
  				result[i] = p1[i] / p2[i];
  		}
517876d6   heziqi   metrics finished ...
545
546
547
  
  		free(p1);
  		free(p2);
43dec788   David Mayerich   added code for si...
548
  		return true;
70407ea9   heziqi   Ziqi completed he...
549
  	}
43dec788   David Mayerich   added code for si...
550
  
a08c98c3   David Mayerich   finished document...
551
  	/// Compute the ratio between a peak area and peak height.
77060b79   David Mayerich   started DOxygen d...
552
553
554
555
556
557
558
559
  
  	/// @param lb1 is the label value for the left baseline point for the first peak (numerator)
  	/// @param rb1 is the label value for the right baseline point for the first peak (numerator)
  	/// @param pos1 is the label value for the first peak (numerator) position
  	/// @param lb2 is the label value for the left baseline point for the second peak (denominator)
  	/// @param rb2 is the label value for the right baseline point for the second peak (denominator)
  	/// @param pos2 is the label value for the second peak (denominator) position
  	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
70407ea9   heziqi   Ziqi completed he...
560
561
  	bool pa_to_ph(double lb1, double rb1, double lab1, double rab1,
  					double lb2, double rb2, double pos, T* result){
43dec788   David Mayerich   added code for si...
562
563
564
565
  
  		T* p1 = (T*)malloc(X() * Y() * sizeof(T));
  		T* p2 = (T*)malloc(X() * Y() * sizeof(T));
  
70407ea9   heziqi   Ziqi completed he...
566
567
568
569
  		//get the area and the peak band
  		area(lb1, rb1, lab1, rab1, p1);
  		height(lb2, rb2, pos, p2);
  		//calculate the ratio in result
724ec347   David Mayerich   simplified the EN...
570
  		for(unsigned i = 0; i < X() * Y(); i++){
70407ea9   heziqi   Ziqi completed he...
571
572
573
574
575
  			if(p1[i] == 0 && p2[i] ==0)
  				result[i] = 1;
  			else
  				result[i] = p1[i] / p2[i];
  		}
517876d6   heziqi   metrics finished ...
576
577
578
  
  		free(p1);
  		free(p2);
43dec788   David Mayerich   added code for si...
579
580
581
  		return true;
  	}
  
a08c98c3   David Mayerich   finished document...
582
  	/// Compute the ratio between two peak areas.
0722eb39   David Mayerich   further updates t...
583
584
585
586
587
588
589
590
  
  	/// @param lb1 is the label value for the left baseline point for the first peak (numerator)
  	/// @param rb1 is the label value for the right baseline point for the first peak (numerator)
  	/// @param lab1 is the label value for the left bound (start of the integration) of the first peak (numerator)
  	/// @param rab1 is the label value for the right bound (end of the integration) of the first peak (numerator)
  	/// @param lb2 is the label value for the left baseline point for the second peak (denominator)
  	/// @param rb2 is the label value for the right baseline point for the second peak (denominator)
  	/// @param lab2 is the label value for the left bound (start of the integration) of the second peak (denominator)
43dec788   David Mayerich   added code for si...
591
  	/// @param rab2 is the label value for the right bound (end of the integration) of the second peak (denominator)
0722eb39   David Mayerich   further updates t...
592
  	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
70407ea9   heziqi   Ziqi completed he...
593
594
  	bool pa_to_pa(double lb1, double rb1, double lab1, double rab1,
  					double lb2, double rb2, double lab2, double rab2, T* result){
43dec788   David Mayerich   added code for si...
595
596
597
598
  
  		T* p1 = (T*)malloc(X() * Y() * sizeof(T));
  		T* p2 = (T*)malloc(X() * Y() * sizeof(T));
  
70407ea9   heziqi   Ziqi completed he...
599
600
601
602
  		//get the area and the peak band
  		area(lb1, rb1, lab1, rab1, p1);
  		area(lb2, rb2, lab2, rab2, p2);
  		//calculate the ratio in result
724ec347   David Mayerich   simplified the EN...
603
  		for(unsigned i = 0; i < X() * Y(); i++){
70407ea9   heziqi   Ziqi completed he...
604
605
606
607
608
  			if(p1[i] == 0 && p2[i] ==0)
  				result[i] = 1;
  			else
  				result[i] = p1[i] / p2[i];
  		}
517876d6   heziqi   metrics finished ...
609
610
611
  
  		free(p1);
  		free(p2);
43dec788   David Mayerich   added code for si...
612
613
  		return true;
  	}
517876d6   heziqi   metrics finished ...
614
  
a08c98c3   David Mayerich   finished document...
615
  	/// Compute the definite integral of a baseline corrected peak.
0722eb39   David Mayerich   further updates t...
616
617
618
619
620
621
  
  	/// @param lb is the label value for the left baseline point
  	/// @param rb is the label value for the right baseline point
  	/// @param lab is the label for the start of the definite integral
  	/// @param rab is the label for the end of the definite integral
  	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
517876d6   heziqi   metrics finished ...
622
623
624
625
626
627
  	bool x_area(double lb, double rb, double lab, double rab, T* result){
  		T* lp;	//left band pointer
  		T* rp;	//right band pointer
  		T* cur;		//current band 1
  		T* cur2;	//current band 2
  
724ec347   David Mayerich   simplified the EN...
628
  		unsigned XY = X() * Y();
517876d6   heziqi   metrics finished ...
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
  		unsigned S = XY * sizeof(T);
  
  		lp = (T*) malloc(S);			//memory allocation
  		rp = (T*) malloc(S);
  		cur = (T*) malloc(S);
  		cur2 = (T*) malloc(S);
  
  		memset(result, (char)0, S);
  
  		//find the wavelenght position in the whole band
  		unsigned int n = w.size();
  		unsigned int ai = 0;		//left bound position
  		unsigned int bi = n - 1;		//right bound position
  
  		//to make sure the left and the right bound are in the bandwidth
  		if (lb < w[0] || rb < w[0] || lb > w[n-1] || rb >w[n-1]){
  			std::cout<<"ERROR: left bound or right bound out of bandwidth"<<std::endl;
  			exit(1);
  		}
  		//to make sure rigth bound is bigger than left bound
  		else if(lb > rb){
  			std::cout<<"ERROR: right bound should be bigger than left bound"<<std::endl;
  			exit(1);
  		}
  
  		//get the position of lb and rb
  		while (lab >= w[ai]){
  			ai++;
  		}
  		while (rab <= w[bi]){
  			bi--;
  		}
  
  		band(lp, lb);
  		band(rp, rb);
  
  		//calculate the beginning and the ending part
  		baseline_band(lb, rb, lp, rp, rab, cur2);		//ending part
  		baseline_band(lb, rb, lp, rp, w[bi], cur);
  		for(unsigned j = 0; j < XY; j++){
  			result[j] += (rab - w[bi]) * (rab + w[bi]) * (cur[j] + cur2[j]) / 4.0;
  		}
  		baseline_band(lb, rb, lp, rp, lab, cur2);		//beginnning part
  		baseline_band(lb, rb, lp, rp, w[ai], cur);
43dec788   David Mayerich   added code for si...
673
  		for(unsigned j = 0; j < XY; j++){
517876d6   heziqi   metrics finished ...
674
675
676
677
678
679
680
681
682
683
  			result[j] += (w[ai] - lab) * (w[ai] + lab) * (cur[j] + cur2[j]) / 4.0;
  		}
  
  		//calculate f(x) times x
  		ai++;
  		for(unsigned i = ai; i <= bi ;i++)
  		{
  			baseline_band(lb, rb, lp, rp, w[ai], cur2);
  			for(unsigned j = 0; j < XY; j++)
  			{
43dec788   David Mayerich   added code for si...
684
  				result[j] += (w[ai] - w[ai-1]) * (w[ai] + w[ai-1]) * (cur[j] + cur2[j]) / 4.0;
517876d6   heziqi   metrics finished ...
685
686
687
688
689
690
691
692
693
694
695
  			}
  			std::swap(cur,cur2);		//swap the band pointers
  		}
  
  		free(lp);
  		free(rp);
  		free(cur);
  		free(cur2);
  		return true;
  	}
  
a08c98c3   David Mayerich   finished document...
696
  	/// Compute the centroid of a baseline corrected peak.
0722eb39   David Mayerich   further updates t...
697
698
699
700
701
702
  
  	/// @param lb is the label value for the left baseline point
  	/// @param rb is the label value for the right baseline point
  	/// @param lab is the label for the start of the peak
  	/// @param rab is the label for the end of the peak
  	/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
517876d6   heziqi   metrics finished ...
703
  	bool cpoint(double lb, double rb, double lab, double rab, T* result){
43dec788   David Mayerich   added code for si...
704
705
706
  		T* p1 = (T*)malloc(X() * Y() * sizeof(T));
  		T* p2 = (T*)malloc(X() * Y() * sizeof(T));
  
517876d6   heziqi   metrics finished ...
707
708
709
710
  		//get the area and the peak band
  		x_area(lb, rb, lab, rab, p1);
  		area(lb, rb, lab, rab, p2);
  		//calculate the ratio in result
724ec347   David Mayerich   simplified the EN...
711
  		for(unsigned i = 0; i < X() * Y(); i++){
517876d6   heziqi   metrics finished ...
712
713
714
715
716
717
718
719
  			if(p1[i] == 0 && p2[i] ==0)
  				result[i] = 1;
  			else
  				result[i] = p1[i] / p2[i];
  		}
  
  		free(p1);
  		free(p2);
43dec788   David Mayerich   added code for si...
720
  		return true;
517876d6   heziqi   metrics finished ...
721
722
  	}
  
a08c98c3   David Mayerich   finished document...
723
  	/// Create a mask based on a given band and threshold value.
0722eb39   David Mayerich   further updates t...
724
725
726
727
728
729
  
  	/// All pixels in the
  	/// specified band greater than the threshold are true and all pixels less than the threshold are false.
  	/// @param mask_band is the band used to specify the mask
  	/// @param threshold is the threshold used to determine if the mask value is true or false
  	/// @param p is a pointer to a pre-allocated array at least X * Y in size
f4c5d71b   David Mayerich   started working w...
730
  	bool build_mask(double mask_band, double threshold, unsigned char* p = NULL){
4a6f666c   heziqi   Added mask method
731
  
724ec347   David Mayerich   simplified the EN...
732
  		T* temp = (T*)malloc(X() * Y() * sizeof(T));		//allocate memory for the certain band
4a6f666c   heziqi   Added mask method
733
734
  		band(temp, mask_band);
  
724ec347   David Mayerich   simplified the EN...
735
  		for (unsigned int i = 0; i < X() * Y(); i++) {
4a6f666c   heziqi   Added mask method
736
737
738
739
740
741
  				if (temp[i] < threshold)
  					p[i] = 0;
  				else
  					p[i] = 255;
  		}
  
517876d6   heziqi   metrics finished ...
742
  		free(temp);
4a6f666c   heziqi   Added mask method
743
744
745
  		return true;
  
  	}
740f8cd2   heziqi   added apply_mask
746
  
a08c98c3   David Mayerich   finished document...
747
  	/// Apply a mask file to the BSQ image, setting all values outside the mask to zero.
0722eb39   David Mayerich   further updates t...
748
749
750
  
  	/// @param outfile is the name of the masked output file
  	/// @param p is a pointer to memory of size X * Y, where p(i) = 0 for pixels that will be set to zero.
740f8cd2   heziqi   added apply_mask
751
752
753
754
  	bool apply_mask(std::string outfile, unsigned char* p){
  
  		std::ofstream target(outfile.c_str(), std::ios::binary);
  
724ec347   David Mayerich   simplified the EN...
755
  		unsigned XY = X() * Y();		//calculate number of a band
740f8cd2   heziqi   added apply_mask
756
757
758
759
  		unsigned L = XY * sizeof(T);
  
  		T * temp = (T*)malloc(L);
  
724ec347   David Mayerich   simplified the EN...
760
  		for (unsigned i = 0; i < Z(); i++)			//for each spectral bin
740f8cd2   heziqi   added apply_mask
761
  		{
e843658b   Brad Deutsch   Previous push did...
762
763
  			band_index(temp, i);					//get the specified band (by index)
  			for ( unsigned j = 0; j < XY; j++)		// for each pixel
740f8cd2   heziqi   added apply_mask
764
  			{
e843658b   Brad Deutsch   Previous push did...
765
766
  				if(p[j] == 0){						//if the mask is 0 at that pixel
  					temp[j] = 0;					//set temp to zero
740f8cd2   heziqi   added apply_mask
767
768
769
770
771
  				}
  				else{
  					continue;
  				}
  			}
e843658b   Brad Deutsch   Previous push did...
772
  			target.write(reinterpret_cast<const char*>(temp), L);   //write the XY slice at that band to disk
740f8cd2   heziqi   added apply_mask
773
774
775
776
777
778
  		}
  		target.close();
  		free(temp);
  		return true;
  	}
  
3d0b6243   David Mayerich   fixed hsiproc bug...
779
780
781
782
783
784
  	/// Copies all spectra corresponding to nonzero values of a mask into a pre-allocated matrix of size (B x P)
  	///		where P is the number of masked pixels and B is the number of bands. The allocated memory can be accessed
  	///		using the following indexing: i = p*B + b
  	/// @param matrix is the destination for the pixel data
  	/// @param mask is the mask
  	bool sift(T* matrix, unsigned char* mask){
43dec788   David Mayerich   added code for si...
785
786
  		unsigned long long XY = X() * Y(); 					//Number of XY pixels
  		unsigned long long L = XY * sizeof(T); 				//size of XY plane (in bytes)
3d0b6243   David Mayerich   fixed hsiproc bug...
787
788
789
790
791
792
793
  
  		//calculate the number of pixels in the mask
  		unsigned long long P = 0;
  		if(mask == NULL) P = XY;
  		else{
  			for(unsigned long long xy = 0; xy < XY; xy++){
  				if(mask[xy] != 0) P++;
43dec788   David Mayerich   added code for si...
794
  			}
3d0b6243   David Mayerich   fixed hsiproc bug...
795
  		}
eb5dfb2b   David Mayerich   fixed linux compa...
796
  
3d0b6243   David Mayerich   fixed hsiproc bug...
797
798
799
800
801
802
803
804
805
806
807
  		T* band_image = (T*) malloc( XY * sizeof(T));		//allocate space for a single band
  
  		unsigned long long i;								//pixel index into the sifted array
  		for(unsigned long long b = 0; b < Z(); b++){		//for each band in the data set
  			band_index(band_image, b);						//retrieve an image of that band
  
  			i = 0;
  			for(unsigned long long xy = 0; xy < XY; xy++){
  				if(mask == NULL || mask[xy] != 0){				//if the pixel is valid
  					matrix[i*Z() + b] = band_image[xy];			//copy it to the appropriate point in the values[] array
  					i++;
43dec788   David Mayerich   added code for si...
808
809
810
  				}
  			}
  		}
eb5dfb2b   David Mayerich   fixed linux compa...
811
  
43dec788   David Mayerich   added code for si...
812
813
814
815
816
817
  		return true;
  	}
  
  	/// Saves only those spectra corresponding to mask value != 0 to pre-allocated memory
  	/// @param matrix is the destination for the sifted pixels
  	/// @param p is the mask file used for sifting
3d0b6243   David Mayerich   fixed hsiproc bug...
818
  	/*bool sift(T* matrix, unsigned char* p){
43dec788   David Mayerich   added code for si...
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
  
  		// open a band (XY plane)
  		unsigned long XY = X() * Y(); 					//Number of XY pixels
  		unsigned long L = XY * sizeof(T); 				//size of XY pixels
  
  		//count the number of masked pixels
  		unsigned long P = 0;												//allocate space for the number of pixels
  		for(unsigned long xy = 0; xy < XY; xy++){ 	//count the number of pixels
  			if(p[xy] != 0) P++;
  		}
  
  		T* bandvals = (T*) malloc(sizeof(T) * P);			//allocate a temporary array to store the pixels for a single band
  		memset(matrix, 0, sizeof(T) * P * Z());
  		for (unsigned long z = 0; z < Z(); z++){		//for each band
  
  			if(!sift_band(bandvals, z, p)){
  				std::cout<<"ERROR sifting band index "<<z<<std::endl;
  				exit(1);
  			}
  
  			for(unsigned long p = 0; p < P; p++)		//for each pixel in the array representing a single band
  				matrix[p * Z() + z] = bandvals[p];			//copy it to the appropriate location in the matrix
  		}
  
  		free(bandvals);				//clean up temporary memory
  		return true;
3d0b6243   David Mayerich   fixed hsiproc bug...
845
  	}*/
43dec788   David Mayerich   added code for si...
846
  
2fcab4f0   David Mayerich   added unsifting f...
847
  	/// Saves to disk only those spectra corresponding to mask values != 0
43dec788   David Mayerich   added code for si...
848
849
  	/// @param outfile is the name of the sifted ENVI file to be written to disk
  	/// @param unsigned char* p is the mask file used for sifting
2fcab4f0   David Mayerich   added unsifting f...
850
  	bool sift(std::string outfile, unsigned char* p){
e843658b   Brad Deutsch   Previous push did...
851
852
  		std::ofstream target(outfile.c_str(), std::ios::binary);
  		// open a band (XY plane)
3150c537   David Mayerich   added further pro...
853
854
  		unsigned long XY = X() * Y(); //Number of XY pixels
  		unsigned long L = XY * sizeof(T); //size of XY pixels
e843658b   Brad Deutsch   Previous push did...
855
856
857
858
  
  		T * temp = (T*)malloc(L); //allocate memory for a band
  		T * temp_vox = (T*)malloc(sizeof(T)); //allocate memory for one voxel
  
3150c537   David Mayerich   added further pro...
859
  		for (unsigned long i = 0; i < Z(); i++)			//for each spectral bin
e843658b   Brad Deutsch   Previous push did...
860
861
  		{
  			band_index(temp, i);					//get the specified band (XY sheet by index)
3150c537   David Mayerich   added further pro...
862
  			for (unsigned long j = 0; j < XY; j++)		// for each pixel
e843658b   Brad Deutsch   Previous push did...
863
864
865
866
867
868
869
870
  			{
  				if (p[j] != 0){						//if the mask is != 0 at that pixel
  					temp_vox[0] = temp[j];
  					target.write(reinterpret_cast<const char*>(temp_vox), sizeof(T));   //write the XY slice at that band to disk
  				}
  				else{
  					continue;
  				}
3150c537   David Mayerich   added further pro...
871
872
  
  				thread_data = (double)(i * XY + j) / (XY * Z()) * 100;
e843658b   Brad Deutsch   Previous push did...
873
874
875
876
  			}
  		}
  		target.close();
  		free(temp);
3150c537   David Mayerich   added further pro...
877
878
879
  
  		thread_data = 100;
  
e843658b   Brad Deutsch   Previous push did...
880
881
882
  		return true;
  	}
  
2fcab4f0   David Mayerich   added unsifting f...
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
  	/// Generates a spectral image from a matrix of spectral values in lexicographic order and a mask
  	bool unsift(std::string outfile, unsigned char* p, unsigned int samples, unsigned int lines){
  
  		//create a binary output stream
  		std::ofstream target(outfile.c_str(), std::ios::binary);
  
  		//make sure that there's only one line
  		if(Y() != 1){
  			std::cout<<"ERROR in stim::bsq::sift() - number of lines does not equal 1"<<std::endl;
  			return false;
  		}
  
  		std::cout<<"started sifting"<<std::endl;
  
  		//get the number of pixels and bands in the input image
3150c537   David Mayerich   added further pro...
898
899
900
  		unsigned long P = X(); 					//Number of pixels
  		unsigned long B = Z();					//number of bands
  		unsigned long XY = samples * lines;		//total number of pixels in an unsifted image
2fcab4f0   David Mayerich   added unsifting f...
901
902
903
904
905
906
907
908
  
  		// allocate memory for a sifted band
  		T * sifted = (T*)malloc(P * sizeof(T)); //allocate memory for a band
  
  		//allocate memory for an unsifted band image
  		T* unsifted = (T*) malloc(XY * sizeof(T));
  
  		//for each band
3150c537   David Mayerich   added further pro...
909
  		for(unsigned long b = 0; b < B; b++){
2fcab4f0   David Mayerich   added unsifting f...
910
911
  
  			//set the unsifted index value to zero
3150c537   David Mayerich   added further pro...
912
  			unsigned long i = 0;
2fcab4f0   David Mayerich   added unsifting f...
913
914
915
916
917
  
  			//retrieve the sifted band (masked pixels only)
  			band_index(sifted, b);
  
  			//for each pixel in the final image (treat it as a 1D image)
3150c537   David Mayerich   added further pro...
918
  			for(unsigned long xi = 0; xi < XY; xi++){
2fcab4f0   David Mayerich   added unsifting f...
919
920
921
922
923
924
  				if( p[xi] == 0 )
  					unsifted[xi] = 0;
  				else{
  					unsifted[xi] = sifted[i];
  					i++;
  				}
3150c537   David Mayerich   added further pro...
925
926
  				//std::cout<<xi<<"/"<<XY<<",   "<<b<<"/"<<B<<std::endl;
  				thread_data = (double)(b * XY + xi) / (B * XY) * 100;
2fcab4f0   David Mayerich   added unsifting f...
927
  			}
3150c537   David Mayerich   added further pro...
928
  			//std::cout<<b*XY<<"/"<<B*XY<<"      "<<B<<"-"<<XY<<std::endl;
2fcab4f0   David Mayerich   added unsifting f...
929
930
931
932
  			//write the band image to disk
  			target.write(reinterpret_cast<const char*>(unsifted), sizeof(T) * XY);
  		}
  
3150c537   David Mayerich   added further pro...
933
934
935
  		//std::cout<<"unsifted"<<std::endl;
  		thread_data = 100;
  
2fcab4f0   David Mayerich   added unsifting f...
936
937
938
  		return true;
  	}
  
e843658b   Brad Deutsch   Previous push did...
939
  
a08c98c3   David Mayerich   finished document...
940
  	/// Calculate the mean band value (average along B) at each pixel location.
0722eb39   David Mayerich   further updates t...
941
942
  
  	/// @param p is a pointer to memory of size X * Y * sizeof(T) that will store the band averages.
a43c4fe1   heziqi   Added crop in env...
943
  	bool band_avg(T* p){
724ec347   David Mayerich   simplified the EN...
944
  		unsigned long long XY = X() * Y();
a43c4fe1   heziqi   Added crop in env...
945
946
947
948
  		T* temp = (T*)malloc(sizeof(T) * XY);
  		//initialize p
  		band_index(p, 0);
  		for (unsigned j = 0; j < XY; j++){
724ec347   David Mayerich   simplified the EN...
949
  			p[j] /= (T)Z();
a43c4fe1   heziqi   Added crop in env...
950
951
  		}
  		//get every band and add them all
724ec347   David Mayerich   simplified the EN...
952
  		for (unsigned i = 1; i < Z(); i++){
a43c4fe1   heziqi   Added crop in env...
953
954
  			band_index(temp, i);
  			for (unsigned j = 0; j < XY; j++){
724ec347   David Mayerich   simplified the EN...
955
  				p[j] += temp[j]/(T)Z();
a43c4fe1   heziqi   Added crop in env...
956
957
958
959
960
961
  			}
  		}
  		free(temp);
  		return true;
  	}
  
a08c98c3   David Mayerich   finished document...
962
963
964
965
  	/// Calculate the mean value for all masked (or valid) pixels in a band and returns the average spectrum
  
  	/// @param p is a pointer to pre-allocated memory of size [B * sizeof(T)] that stores the mean spectrum
  	/// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
2aa68315   David Mayerich   major bug fixes f...
966
  	bool avg_band(double* p, unsigned char* mask = NULL){
724ec347   David Mayerich   simplified the EN...
967
  		unsigned long long XY = X() * Y();
2aa68315   David Mayerich   major bug fixes f...
968
  		unsigned long long count = 0;						//count will store the number of masked pixels
a43c4fe1   heziqi   Added crop in env...
969
  		T* temp = (T*)malloc(sizeof(T) * XY);
a08c98c3   David Mayerich   finished document...
970
  		//calculate this loop counts the number of true pixels in the mask
a43c4fe1   heziqi   Added crop in env...
971
  		for (unsigned j = 0; j < XY; j++){
2aa68315   David Mayerich   major bug fixes f...
972
  			if (mask == NULL || mask[j] != 0){
a43c4fe1   heziqi   Added crop in env...
973
974
975
  				count++;
  			}
  		}
724ec347   David Mayerich   simplified the EN...
976
  		//this loops goes through each band in B (Z())
a08c98c3   David Mayerich   finished document...
977
  		//	masked (or valid) pixels from that band are averaged and the average is stored in p
2aa68315   David Mayerich   major bug fixes f...
978
  		for (unsigned long long i = 0; i < Z(); i++){
a43c4fe1   heziqi   Added crop in env...
979
  			p[i] = 0;
a08c98c3   David Mayerich   finished document...
980
  			band_index(temp, i);				//get the band image and store it in temp
2aa68315   David Mayerich   major bug fixes f...
981
982
983
  			for (unsigned long long j = 0; j < XY; j++){	//loop through temp, averaging valid pixels
  				if (mask == NULL || mask[j] != 0){
  					p[i] += (double)temp[j] / (double)count;
43dec788   David Mayerich   added code for si...
984
  				}
a43c4fe1   heziqi   Added crop in env...
985
986
987
988
989
990
  			}
  		}
  		free(temp);
  		return true;
  	}
  
a08c98c3   David Mayerich   finished document...
991
992
993
994
995
  	/// Calculate the covariance matrix for all masked pixels in the image.
  
  	/// @param co is a pointer to pre-allocated memory of size [B * B] that stores the resulting covariance matrix
  	/// @param avg is a pointer to memory of size B that stores the average spectrum
  	/// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
2aa68315   David Mayerich   major bug fixes f...
996
997
  	bool co_matrix(double* co, double* avg, unsigned char *mask){
  		thread_data = 0;							//initialize thread data for timers and progress bars
a43c4fe1   heziqi   Added crop in env...
998
  		//memory allocation
724ec347   David Mayerich   simplified the EN...
999
  		unsigned long long xy = X() * Y();
2aa68315   David Mayerich   major bug fixes f...
1000
  		unsigned long long B = Z();
a43c4fe1   heziqi   Added crop in env...
1001
1002
  		T* bandi = (T*)malloc(sizeof(T) * xy);
  		T* bandj = (T*)malloc(sizeof(T) * xy);
43dec788   David Mayerich   added code for si...
1003
  
a43c4fe1   heziqi   Added crop in env...
1004
  		//count vaild pixels in a band
2aa68315   David Mayerich   major bug fixes f...
1005
1006
1007
  		unsigned long long count = 0;
  		for (unsigned long long j = 0; j < xy; j++){
  			if (mask == NULL || mask[j] != 0){
a43c4fe1   heziqi   Added crop in env...
1008
1009
1010
1011
  				count++;
  			}
  		}
  		//calculate correlation coefficient matrix
2aa68315   David Mayerich   major bug fixes f...
1012
  		for (unsigned long long i = 0; i < B; i++)
a43c4fe1   heziqi   Added crop in env...
1013
1014
  		{
  			band_index(bandi, i);
2aa68315   David Mayerich   major bug fixes f...
1015
  			for (unsigned long long j = i; j < B; j++){
a43c4fe1   heziqi   Added crop in env...
1016
  				band_index(bandj, j);
2aa68315   David Mayerich   major bug fixes f...
1017
  				double numerator = 0;			//to calculate element in correlation coefficient matrix, numerator part
a43c4fe1   heziqi   Added crop in env...
1018
  				//calculate the R(i,j) in correlation coeffient matrix
2aa68315   David Mayerich   major bug fixes f...
1019
1020
1021
  				for (unsigned long long k = 0; k < xy; k++){
  					if (mask == NULL || mask[k] != 0){
  						numerator += ((double)bandi[k] - (double)avg[i]) * ((double)bandj[k] - (double)avg[j]) / (double)count;
a43c4fe1   heziqi   Added crop in env...
1022
1023
1024
1025
1026
  					}
  				}
  				co[i*B + j] = numerator;
  				co[j*B + i] = numerator;		//because correlated matrix is a symmetric matrix
  			}
2aa68315   David Mayerich   major bug fixes f...
1027
  			thread_data = (double)i / B * 100;
a43c4fe1   heziqi   Added crop in env...
1028
1029
1030
  		}
  		free(bandi);
  		free(bandj);
2aa68315   David Mayerich   major bug fixes f...
1031
1032
1033
  
  		thread_data = 100;							//processing complete
  
a43c4fe1   heziqi   Added crop in env...
1034
1035
  		return true;
  	}
0df38ff3   heziqi   fixed head detached
1036
  
a43c4fe1   heziqi   Added crop in env...
1037
  
a08c98c3   David Mayerich   finished document...
1038
1039
1040
1041
1042
1043
1044
  	/// Crop a region of the image and save it to a new file.
  
  	/// @param outfile is the file name for the new cropped image
  	/// @param x0 is the lower-left x pixel coordinate to be included in the cropped image
  	/// @param y0 is the lower-left y pixel coordinate to be included in the cropped image
  	/// @param x1 is the upper-right x pixel coordinate to be included in the cropped image
  	/// @param y1 is the upper-right y pixel coordinate to be included in the cropped image
e9bddc57   David Mayerich   added band croppi...
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
  	bool crop(std::string outfile, unsigned long long x0,
  								   unsigned long long y0,
  								   unsigned long long x1,
  								   unsigned long long y1,
  								   unsigned long long b0,
  								   unsigned long long b1){
  
  		//calculate the new number of samples, lines, and bands
  		unsigned long long samples = x1 - x0;
  		unsigned long long lines = y1 - y0;
  		unsigned long long bands = b1 - b0;
  
  		//calculate the size of a single band
  		unsigned long long L = samples * lines * sizeof(T);
  
  		//allocate space for a single band
a43c4fe1   heziqi   Added crop in env...
1061
  		T* temp = (T*)malloc(L);
e9bddc57   David Mayerich   added band croppi...
1062
1063
  
  		//create an output stream to store the output file
a43c4fe1   heziqi   Added crop in env...
1064
  		std::ofstream out(outfile.c_str(), std::ios::binary);
e9bddc57   David Mayerich   added band croppi...
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
  
  		//calculate the distance required to jump from the end of one band to the beginning of another
  		unsigned long long jumpb = X() * (Y() - lines) * sizeof(T);
  
  		//calculate the distance required to jump from the end of one line to the beginning of another
  		unsigned long long jumpl = (X() - samples) * sizeof(T);
  
  		//seek to the start of the cropped region in the input file
  		file.seekg( (b0 * X() * Y() + y0 * X() + x0) * sizeof(T), std::ios::beg);
  
  		//for each band
  		for (unsigned long long z = b0; z < b1; z++)
a43c4fe1   heziqi   Added crop in env...
1077
  		{
e9bddc57   David Mayerich   added band croppi...
1078
  			for (unsigned y = 0; y < lines; y++)
a43c4fe1   heziqi   Added crop in env...
1079
  			{
e9bddc57   David Mayerich   added band croppi...
1080
  				file.read((char *)(temp + y * samples), sizeof(T) * samples);
43dec788   David Mayerich   added code for si...
1081
  				file.seekg(jumpl, std::ios::cur);    //go to the next band
3150c537   David Mayerich   added further pro...
1082
  
e9bddc57   David Mayerich   added band croppi...
1083
  				thread_data = (double)(z * lines + y) / (Z() * lines) * 100;
a43c4fe1   heziqi   Added crop in env...
1084
  			}
43dec788   David Mayerich   added code for si...
1085
  			out.write(reinterpret_cast<const char*>(temp), L);   //write slice data into target file
a43c4fe1   heziqi   Added crop in env...
1086
1087
1088
  			file.seekg(jumpb, std::ios::cur);
  		}
  		free(temp);
3150c537   David Mayerich   added further pro...
1089
1090
1091
  
  		thread_data = 100;
  
a43c4fe1   heziqi   Added crop in env...
1092
1093
  		return true;
  	}
0df38ff3   heziqi   fixed head detached
1094
  
740f8cd2   heziqi   added apply_mask
1095
  
a08c98c3   David Mayerich   finished document...
1096
  	/// Close the file.
f6169dea   heziqi   Ziqi completed bi...
1097
1098
1099
1100
  	bool close(){
  		file.close();
  		return true;
  	}
31394e2a   heziqi   Ziqi added method...
1101
  
3fc1d461   heziqi   Ziqi added saveBa...
1102
1103
  	};
  }
31394e2a   heziqi   Ziqi added method...
1104
  
e8eb202f   David Mayerich   added a new ENVI ...
1105
  #endif
43dec788   David Mayerich   added code for si...