6aa04ba2
David Mayerich
interleave types
|
1
2
3
|
#ifndef STIM_BIL_H
#define STIM_BIL_H
|
e8eb202f
David Mayerich
added a new ENVI ...
|
4
|
#include "../envi/envi_header.h"
|
c25e7d0d
heziqi
speed of bip.base...
|
5
6
7
8
9
10
11
12
13
14
15
16
|
#include "../envi/binary.h"
#include <cstring>
#include <utility>
namespace rts{
template <typename T>
class bil: public binary<T> {
protected:
|
6708cc25
heziqi
Ziqi added envi c...
|
17
|
std::vector<double> w; //band wavelength
|
c25e7d0d
heziqi
speed of bip.base...
|
18
19
20
21
22
|
public:
using binary<T>::open;
using binary<T>::file;
|
6708cc25
heziqi
Ziqi added envi c...
|
23
|
using binary<T>::R;
|
c25e7d0d
heziqi
speed of bip.base...
|
24
25
|
//open a file, given the file and its header's names
|
6708cc25
heziqi
Ziqi added envi c...
|
26
|
bool open(std::string filename, unsigned int X, unsigned int Y, unsigned int B, unsigned int header_offset, std::vector<double> wavelengths){
|
c25e7d0d
heziqi
speed of bip.base...
|
27
|
|
6708cc25
heziqi
Ziqi added envi c...
|
28
|
w = wavelengths;
|
c25e7d0d
heziqi
speed of bip.base...
|
29
|
|
6708cc25
heziqi
Ziqi added envi c...
|
30
31
32
|
return open(filename, vec<unsigned int>(X, Y, B), header_offset);
return false;
|
c25e7d0d
heziqi
speed of bip.base...
|
33
34
35
36
37
38
|
}
//save one band of the file into the memory, and return the pointer
bool band_index( T * p, unsigned int page){
|
6708cc25
heziqi
Ziqi added envi c...
|
39
40
|
unsigned int L = R[0] * sizeof(T); //caculate the number of bytes in a sample line
unsigned int jump = R[0] * (R[2] - 1) * sizeof(T);
|
c25e7d0d
heziqi
speed of bip.base...
|
41
|
|
6708cc25
heziqi
Ziqi added envi c...
|
42
|
if (page >= R[2]){ //make sure the bank number is right
|
c25e7d0d
heziqi
speed of bip.base...
|
43
44
45
46
|
std::cout<<"ERROR: page out of range"<<std::endl;
return false;
}
|
6708cc25
heziqi
Ziqi added envi c...
|
47
48
|
file.seekg(R[0] * page * sizeof(T), std::ios::beg);
for (unsigned i = 0; i < R[1]; i++)
|
c25e7d0d
heziqi
speed of bip.base...
|
49
|
{
|
6708cc25
heziqi
Ziqi added envi c...
|
50
|
file.read((char *)(p + i * R[0]), L);
|
c25e7d0d
heziqi
speed of bip.base...
|
51
52
53
54
55
56
|
file.seekg( jump, std::ios::cur);
}
return true;
}
|
e933c14e
David Mayerich
cleaned up the code
|
57
|
bool band( T * p, double wavelength){
|
c25e7d0d
heziqi
speed of bip.base...
|
58
|
|
6708cc25
heziqi
Ziqi added envi c...
|
59
|
unsigned int XY = R[0] * R[1]; //calculate the number of pixels in a band
|
c25e7d0d
heziqi
speed of bip.base...
|
60
61
62
63
64
65
66
67
68
|
unsigned int S = XY * sizeof(T); //calculate the number of bytes of a band
unsigned page=0; //bands around the wavelength
T * p1;
T * p2;
//get the bands numbers around the wavelength
//if wavelength is smaller than the first one in header file
|
6708cc25
heziqi
Ziqi added envi c...
|
69
|
if ( w[page] > wavelength ){
|
c25e7d0d
heziqi
speed of bip.base...
|
70
71
72
73
|
band_index(p, page);
return true;
}
|
6708cc25
heziqi
Ziqi added envi c...
|
74
|
while( w[page] < wavelength )
|
c25e7d0d
heziqi
speed of bip.base...
|
75
76
77
|
{
page++;
//if wavelength is larger than the last wavelength in header file
|
6708cc25
heziqi
Ziqi added envi c...
|
78
79
|
if (page == R[2]) {
band_index(p, R[2]-1);
|
c25e7d0d
heziqi
speed of bip.base...
|
80
81
82
|
return true;
}
}
|
6708cc25
heziqi
Ziqi added envi c...
|
83
|
if ( wavelength < w[page] )
|
c25e7d0d
heziqi
speed of bip.base...
|
84
85
86
87
88
89
|
{
p1=(T*)malloc(S); //memory allocation
p2=(T*)malloc(S);
band_index(p1, page - 1);
band_index(p2, page );
for(unsigned i=0; i < XY; i++){
|
6708cc25
heziqi
Ziqi added envi c...
|
90
|
double r = (double) (wavelength - w[page-1]) / (double) (w[page] - w[page-1]);
|
c25e7d0d
heziqi
speed of bip.base...
|
91
92
93
94
95
96
97
98
99
100
101
102
103
104
|
p[i] = (p2[i] - p1[i]) * r + p1[i];
}
}
else //if the wavelength is equal to a wavelength in header file
{
band_index(p, page);
}
return true;
}
//get YZ line from the a Y slice, Y slice data should be already IN the MEMORY
bool getYZ(T* p, T* c, double wavelength)
{
|
6708cc25
heziqi
Ziqi added envi c...
|
105
106
|
unsigned int X = R[0]; //calculate the number of pixels in a sample
unsigned int B = R[2];
|
c25e7d0d
heziqi
speed of bip.base...
|
107
108
109
110
111
112
113
114
115
|
unsigned int L = X * sizeof(T);
unsigned page=0; //samples around the wavelength
T * p1;
T * p2;
//get the bands numbers around the wavelength
//if wavelength is smaller than the first one in header file
|
6708cc25
heziqi
Ziqi added envi c...
|
116
|
if ( w[page] > wavelength ){
|
c25e7d0d
heziqi
speed of bip.base...
|
117
118
119
120
|
memcpy(p, c, L);
return true;
}
|
6708cc25
heziqi
Ziqi added envi c...
|
121
|
while( w[page] < wavelength )
|
c25e7d0d
heziqi
speed of bip.base...
|
122
123
124
125
126
127
128
129
|
{
page++;
//if wavelength is larger than the last wavelength in header file
if (page == B) {
memcpy(p, c + (B - 1) * X, L);
return true;
}
}
|
6708cc25
heziqi
Ziqi added envi c...
|
130
|
if ( wavelength < w[page] )
|
c25e7d0d
heziqi
speed of bip.base...
|
131
132
133
134
135
136
137
138
|
{
p1=(T*)malloc( L ); //memory allocation
p2=(T*)malloc( L );
memcpy(p1, c + (page - 1) * X, L);
memcpy(p2, c + page * X, L);
for(unsigned i=0; i < X; i++){
|
6708cc25
heziqi
Ziqi added envi c...
|
139
|
double r = (double) (wavelength - w[page-1]) / (double) (w[page] - w[page-1]);
|
c25e7d0d
heziqi
speed of bip.base...
|
140
141
142
143
144
145
146
147
148
149
|
p[i] = (p2[i] - p1[i]) * r + p1[i];
}
}
else //if the wavelength is equal to a wavelength in header file
memcpy(p, c + page * X, L);
return true;
}
//save one pixel of the BIP file into the memory, and return the pointer
|
e933c14e
David Mayerich
cleaned up the code
|
150
|
bool spectrum(T * p, unsigned x, unsigned y){
|
c25e7d0d
heziqi
speed of bip.base...
|
151
|
|
6708cc25
heziqi
Ziqi added envi c...
|
152
|
if ( x >= R[0] || y >= R[1]){ //make sure the sample and line number is right
|
c25e7d0d
heziqi
speed of bip.base...
|
153
154
155
156
|
std::cout<<"ERROR: sample or line out of range"<<std::endl;
exit(1);
}
|
6708cc25
heziqi
Ziqi added envi c...
|
157
|
unsigned jump = (R[0] - 1) * sizeof(T);
|
c25e7d0d
heziqi
speed of bip.base...
|
158
|
|
6708cc25
heziqi
Ziqi added envi c...
|
159
|
file.seekg((y * R[0] * R[2] + x) * sizeof(T), std::ios::beg);
|
c25e7d0d
heziqi
speed of bip.base...
|
160
|
|
6708cc25
heziqi
Ziqi added envi c...
|
161
|
for(unsigned i = 0; i < R[2]; i++)
|
c25e7d0d
heziqi
speed of bip.base...
|
162
163
164
165
166
167
168
169
170
171
172
|
{ //point to the certain sample and line
file.read((char *)(p + i), sizeof(T));
file.seekg(jump, std::ios::cur);
}
return true;
}
//given a Y ,return a XZ slice
bool getY(T * p, unsigned y)
{
|
6708cc25
heziqi
Ziqi added envi c...
|
173
|
if ( y >= R[1]){ //make sure the line number is right
|
c25e7d0d
heziqi
speed of bip.base...
|
174
175
176
177
|
std::cout<<"ERROR: line out of range"<<std::endl;
exit(1);
}
|
6708cc25
heziqi
Ziqi added envi c...
|
178
179
|
file.seekg(y * R[2] * R[0] * sizeof(T), std::ios::beg);
file.read((char *)p, sizeof(T) * R[2] * R[0]);
|
c25e7d0d
heziqi
speed of bip.base...
|
180
181
182
183
184
185
186
187
188
189
190
191
192
193
|
return true;
}
//(BIL) baseline correction
bool baseline(std::string outname, std::vector<double> wls){
unsigned N = wls.size(); //get the number of baseline points
std::ofstream target(outname.c_str(), std::ios::binary); //open the target binary file
std::string headername = outname + ".hdr"; //the header file name
//simplify image resolution
|
6708cc25
heziqi
Ziqi added envi c...
|
194
|
unsigned int ZX = R[2] * R[0]; //calculate the number of points in a Y slice
|
c25e7d0d
heziqi
speed of bip.base...
|
195
|
unsigned int L = ZX * sizeof(T); //calculate the number of bytes of a Y slice
|
6708cc25
heziqi
Ziqi added envi c...
|
196
197
|
unsigned int B = R[2];
unsigned int X = R[0];
|
c25e7d0d
heziqi
speed of bip.base...
|
198
199
200
201
202
203
204
205
206
207
208
209
210
|
T* c; //pointer to the current Y slice
c = (T*)malloc(L); //memory allocation
T* a; //pointer to the two YZ lines surrounding the current YZ line
T* b;
a = (T*)malloc(X * sizeof(T));
b = (T*)malloc(X * sizeof(T));
double ai, bi; //stores the two baseline points wavelength surrounding the current band
double ci; //stores the current band's wavelength
|
f6169dea
heziqi
Ziqi completed bi...
|
211
|
unsigned control;
|
c25e7d0d
heziqi
speed of bip.base...
|
212
213
214
215
216
217
218
|
if (a == NULL || b == NULL || c == NULL){
std::cout<<"ERROR: error allocating memory";
exit(1);
}
// loop start correct every y slice
|
6708cc25
heziqi
Ziqi added envi c...
|
219
|
for (unsigned k =0; k < R[1]; k++)
|
c25e7d0d
heziqi
speed of bip.base...
|
220
221
222
223
224
|
{
//get the current y slice
getY(c, k);
//initialize lownum, highnum, low, high
|
6708cc25
heziqi
Ziqi added envi c...
|
225
|
ai = w[0];
|
f6169dea
heziqi
Ziqi completed bi...
|
226
|
control=0;
|
c25e7d0d
heziqi
speed of bip.base...
|
227
228
|
//if no baseline point is specified at band 0,
//set the baseline point at band 0 to 0
|
6708cc25
heziqi
Ziqi added envi c...
|
229
|
if(wls[0] != w[0]){
|
c25e7d0d
heziqi
speed of bip.base...
|
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
|
bi = wls[control];
memset(a, (char)0, X * sizeof(T) );
}
//else get the low band
else{
control++;
getYZ(a, c, ai);
bi = wls[control];
}
//get the high band
getYZ(b, c, bi);
//correct every YZ line
for(unsigned cii = 0; cii < B; cii++){
//update baseline points, if necessary
|
6708cc25
heziqi
Ziqi added envi c...
|
247
|
if( w[cii] >= bi && cii != B - 1) {
|
c25e7d0d
heziqi
speed of bip.base...
|
248
249
250
251
252
253
254
255
256
257
258
259
260
|
//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];
getYZ(b, c, bi);
}
//if the last BL point on the last band of the file?
|
6708cc25
heziqi
Ziqi added envi c...
|
261
|
else if ( wls[control] < w[B - 1]) {
|
c25e7d0d
heziqi
speed of bip.base...
|
262
263
264
265
266
267
|
std::swap(a, b); //swap the baseline band pointers
memset(b, (char)0, X * sizeof(T) ); //clear the high band
ai = bi;
|
6708cc25
heziqi
Ziqi added envi c...
|
268
|
bi = w[B - 1];
|
c25e7d0d
heziqi
speed of bip.base...
|
269
270
|
}
}
|
c25e7d0d
heziqi
speed of bip.base...
|
271
|
|
6708cc25
heziqi
Ziqi added envi c...
|
272
|
ci = w[cii];
|
c25e7d0d
heziqi
speed of bip.base...
|
273
274
275
276
277
278
|
unsigned jump = cii * X;
//perform the baseline correction
for(unsigned i=0; i < X; i++)
{
double r = (double) (ci - ai) / (double) (bi - ai);
|
92e4cf05
heziqi
include file fixed
|
279
|
c[i + jump] =(T) ( c[i + jump] - (b[i] - a[i]) * r - a[i] );
|
c25e7d0d
heziqi
speed of bip.base...
|
280
281
282
|
}
}//loop for YZ line end
|
c25e7d0d
heziqi
speed of bip.base...
|
283
|
target.write(reinterpret_cast<const char*>(c), L); //write the corrected data into destination
|
6708cc25
heziqi
Ziqi added envi c...
|
284
|
}//loop for Y slice end
|
c25e7d0d
heziqi
speed of bip.base...
|
285
286
287
288
289
290
291
292
293
|
free(a);
free(b);
free(c);
target.close();
return true;
}
|
f6169dea
heziqi
Ziqi completed bi...
|
294
|
// normalize the BIL file
|
e933c14e
David Mayerich
cleaned up the code
|
295
|
bool normalize(std::string outname, double w)
|
c25e7d0d
heziqi
speed of bip.base...
|
296
|
{
|
6708cc25
heziqi
Ziqi added envi c...
|
297
298
299
300
301
|
unsigned int B = R[2]; //calculate the number of bands
unsigned int Y = R[1];
unsigned int X = R[0];
unsigned int ZX = R[2] * R[0];
unsigned int XY = R[0] * R[1]; //calculate the number of pixels in a band
|
c25e7d0d
heziqi
speed of bip.base...
|
302
303
304
305
306
307
308
309
310
311
312
313
|
unsigned int S = XY * sizeof(T); //calculate the number of bytes in a band
unsigned int L = ZX * sizeof(T);
std::ofstream target(outname.c_str(), std::ios::binary); //open the target binary file
std::string headername = outname + ".hdr"; //the header file name
T * c; //pointer to the current ZX slice
T * b; //pointer to the standard band
b = (T*)malloc( S ); //memory allocation
c = (T*)malloc( L );
|
e933c14e
David Mayerich
cleaned up the code
|
314
|
band(b, w); //get the certain band into memory
|
c25e7d0d
heziqi
speed of bip.base...
|
315
316
317
318
|
for(unsigned j = 0; j < Y; j++)
{
getY(c, j);
|
f6169dea
heziqi
Ziqi completed bi...
|
319
|
for(unsigned i = 0; i < B; i++)
|
c25e7d0d
heziqi
speed of bip.base...
|
320
|
{
|
f6169dea
heziqi
Ziqi completed bi...
|
321
|
for(unsigned m = 0; m < X; m++)
|
c25e7d0d
heziqi
speed of bip.base...
|
322
|
{
|
f6169dea
heziqi
Ziqi completed bi...
|
323
|
c[m + i * X] = c[m + i * X] / b[m + j * X];
|
c25e7d0d
heziqi
speed of bip.base...
|
324
325
326
327
|
}
}
target.write(reinterpret_cast<const char*>(c), L); //write normalized data into destination
}
|
c25e7d0d
heziqi
speed of bip.base...
|
328
329
330
331
332
333
334
|
free(b);
free(c);
target.close();
return true;
}
|
f6169dea
heziqi
Ziqi completed bi...
|
335
|
//convert BIL file to BSQ file and save it
|
c25e7d0d
heziqi
speed of bip.base...
|
336
337
|
bool bsq(std::string outname)
{
|
6708cc25
heziqi
Ziqi added envi c...
|
338
|
unsigned int S = R[0] * R[1] * sizeof(T); //calculate the number of bytes in a band
|
c25e7d0d
heziqi
speed of bip.base...
|
339
340
341
342
343
344
345
|
std::ofstream target(outname.c_str(), std::ios::binary);
std::string headername = outname + ".hdr";
T * p; //pointer to the current band
p = (T*)malloc(S);
|
6708cc25
heziqi
Ziqi added envi c...
|
346
|
for ( unsigned i = 0; i < R[2]; i++)
|
c25e7d0d
heziqi
speed of bip.base...
|
347
348
349
350
|
{
band_index(p, i);
target.write(reinterpret_cast<const char*>(p), S); //write a band data into target file
}
|
6708cc25
heziqi
Ziqi added envi c...
|
351
|
|
c25e7d0d
heziqi
speed of bip.base...
|
352
353
354
355
356
|
free(p);
target.close();
return true;
}
|
f6169dea
heziqi
Ziqi completed bi...
|
357
358
359
|
//convert bil file to bip file and save it
bool bip(std::string outname)
{
|
6708cc25
heziqi
Ziqi added envi c...
|
360
|
unsigned int S = R[0] * R[2] * sizeof(T); //calculate the number of bytes in a ZX slice
|
f6169dea
heziqi
Ziqi completed bi...
|
361
362
363
364
365
366
367
368
369
|
std::ofstream target(outname.c_str(), std::ios::binary);
std::string headername = outname + ".hdr";
T * p; //pointer to the current XZ slice for bil file
p = (T*)malloc(S);
T * q; //pointer to the current ZX slice for bip file
q = (T*)malloc(S);
|
6708cc25
heziqi
Ziqi added envi c...
|
370
|
for ( unsigned i = 0; i < R[1]; i++)
|
f6169dea
heziqi
Ziqi completed bi...
|
371
372
|
{
getY(p, i);
|
6708cc25
heziqi
Ziqi added envi c...
|
373
|
for ( unsigned k = 0; k < R[2]; k++)
|
f6169dea
heziqi
Ziqi completed bi...
|
374
|
{
|
6708cc25
heziqi
Ziqi added envi c...
|
375
376
377
|
unsigned ks = k * R[0];
for ( unsigned j = 0; j < R[0]; j++)
q[k + j * R[2]] = p[ks + j];
|
f6169dea
heziqi
Ziqi completed bi...
|
378
379
380
381
382
|
}
target.write(reinterpret_cast<const char*>(q), S); //write a band data into target file
}
|
6708cc25
heziqi
Ziqi added envi c...
|
383
|
|
f6169dea
heziqi
Ziqi completed bi...
|
384
385
386
387
388
389
|
free(p);
free(q);
target.close();
return true;
}
|
20c212c0
heziqi
Ziqi added functi...
|
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
|
//providing the left and the right bound data, return baseline-corrected band height
bool baseline_band(double lb, double rb, T* lp, T* rp, double wavelength, T* result){
unsigned XY = R[0] * R[1];
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;
}
//untested
//providing the left and the right bound wavelength, return baseline-corrected band height
bool height(double lb, double rb, double bandwavelength){
T* lp;
T* rp;
T* result;
lp = (T*) malloc(sizeof(R[0] * R[1] * sizeof(T))); //memory allocation
rp = (T*) malloc(sizeof(R[0] * R[1] * sizeof(T)));
|
2b10d15c
David Mayerich
updated framework...
|
414
|
result = (T*) malloc(sizeof(R[0] * R[1] * sizeof(T)));
|
20c212c0
heziqi
Ziqi added functi...
|
415
416
417
418
419
420
421
422
423
|
band(lp, lb); //get the data of the left bounbd and the right bound
band(rp, rb);
baseline_band(lb, rb, lp, rp, bandwavelength, result);
return true;
}
|
20c212c0
heziqi
Ziqi added functi...
|
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
|
//calculate the area between two bound point(including baseline correction)
bool area(double lb, double rb){
T* lp; //left band pointer
T* rp; //right band pointer
T* cur; //current band 1
T* cur2; //current band 2
T* result; //area result
unsigned XY = R[0] * R[1];
unsigned S = XY * sizeof(T);
lp = (T*) malloc(S); //memory allocation
rp = (T*) malloc(S);
cur = (T*) malloc(S);
cur2 = (T*) malloc(S);
result = (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 (lb >= w[ai]){
ai++;
}
while (rb <= w[bi]){
bi--;
}
band(lp, lb);
band(rp, rb);
//calculate the beginning and the ending part
baseline_band(lb, rb, lp, rp, w[bi], cur); //ending part
for(unsigned j = 0; j < XY; j++){
result[j] += (rb - w[bi]) * cur[j] / 2.0;
}
baseline_band(lb, rb, lp, rp, w[ai], cur); //beginnning part
for(unsigned j = 0; j < XY; j++){
result[j] += (w[ai] - lb) * cur[j] / 2.0;
}
//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++)
{
result[j] += (w[ai] - w[ai-1]) * (cur[j] + cur2[j]) / 2.0;
}
std::swap(cur,cur2); //swap the band pointers
}
///////////////for testing use
std::ofstream text("area_result.txt",std::ios::out);
for (unsigned i = 0; i < 640; i++)
{
for (unsigned j = 0; j < 384; j++)
{
text<<result[i*384+j]<<", ";
}
text<<std::endl;
}
///////////////
return true;
}
|
4a6f666c
heziqi
Added mask method
|
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
|
//create mask file
bool mask(unsigned char* p, double mask_band, double threshold){
T* temp = (T*)malloc(R[0] * R[1] * sizeof(T)); //allocate memory for the certain band
band(temp, mask_band);
for (unsigned int i = 0; i < R[0] * R[1]; i++) {
if (temp[i] < threshold)
p[i] = 0;
else
p[i] = 255;
}
return true;
}
|
f6169dea
heziqi
Ziqi completed bi...
|
523
524
525
526
527
528
|
//close the file
bool close(){
file.close();
return true;
}
|
c25e7d0d
heziqi
speed of bip.base...
|
529
|
};
|
cac62fd3
David Mayerich
modified to work ...
|
530
|
}
|
6aa04ba2
David Mayerich
interleave types
|
531
532
|
#endif
|