codebase / stimlib

Commit d4f06a7e89a896d8a1d937afdea9c3840237ebf6

Authored by Pavel Govyadinov
2 parents 3e5d3ad3 c8c976a9

Merge branch 'master' of git.stim.ee.uh.edu:codebase/stimlib

Showing 14 changed files with 1325 additions and 56886 deletions   Show diff stats
Inline Side-by-side
matlab/enviSaveRaw.m
  Show/Hide comments   View file @d4f06a7
1   -%loads an ENVI file without any manipulation (changing orientation)
  1 +%saves an ENVI file without any manipulation, assumes (X, Y, S)
2 2 function enviSaveRaw(M, filename, headername)
3 3  
4 4 %if a header isn't provided, assume it's just the filename
... ...
stim/envi/bil.h
  Show/Hide comments   View file @d4f06a7
... ... @@ -2,7 +2,7 @@
2 2 #define STIM_BIL_H
3 3  
4 4 #include "../envi/envi_header.h"
5   -#include "../envi/binary.h"
  5 +#include "../envi/hsi.h"
6 6 #include <cstring>
7 7 #include <utility>
8 8  
... ... @@ -19,23 +19,29 @@ namespace stim{
19 19  
20 20 template <typename T>
21 21  
22   -class bil: public binary<T> {
  22 +class bil: public hsi<T> {
23 23  
24 24 protected:
25 25  
26   - std::vector<double> w; //band wavelength
  26 + //std::vector<double> w; //band wavelengths
27 27  
28   - unsigned long X(){
  28 + /*unsigned long long X(){
29 29 return R[0];
30 30 }
31   - unsigned long Y(){
  31 + unsigned long long Y(){
32 32 return R[2];
33 33 }
34   - unsigned long Z(){
  34 + unsigned long long Z(){
35 35 return R[1];
36   - }
  36 + }*/
  37 + using hsi<T>::w; //use the wavelength array in stim::hsi
  38 + using hsi<T>::nnz;
  39 + using binary<T>::progress;
37 40  
38   - using binary<T>::thread_data;
  41 + /// Call the binary nnz() function for the BIL orientation
  42 + //unsigned long long nnz(unsigned char* mask){
  43 + // return binary<T>::nnz(mask, X()*Y());
  44 + //}
39 45  
40 46 public:
41 47  
... ... @@ -43,6 +49,8 @@ public:
43 49 using binary<T>::file;
44 50 using binary<T>::R;
45 51  
  52 + bil(){ hsi<T>::init_bil(); }
  53 +
46 54 /// Open a data file for reading using the class interface.
47 55  
48 56 /// @param filename is the name of the binary file on disk
... ... @@ -51,11 +59,16 @@ public:
51 59 /// @param B is the number of samples (bands) along dimension 3
52 60 /// @param header_offset is the number of bytes (if any) in the binary header
53 61 /// @param wavelengths is an optional STL vector of size B specifying a numerical label for each band
54   - bool open(std::string filename, unsigned int X, unsigned int Y, unsigned int B, unsigned int header_offset, std::vector<double> wavelengths){
  62 + bool open(std::string filename,
  63 + unsigned long long X,
  64 + unsigned long long Y,
  65 + unsigned long long B,
  66 + unsigned long long header_offset,
  67 + std::vector<double> wavelengths){
55 68  
56 69 w = wavelengths;
57 70  
58   - return open(filename, vec<unsigned int>(X, B, Y), header_offset);
  71 + return open(filename, vec<unsigned long long>(X, B, Y), header_offset);
59 72  
60 73 }
61 74  
... ... @@ -63,11 +76,12 @@ public:
63 76  
64 77 /// @param p is a pointer to an allocated region of memory at least X * Y * sizeof(T) in size.
65 78 /// @param page <= B is the integer number of the band to be copied.
66   - bool band_index( T * p, unsigned int page){
  79 + bool band_index( T * p, unsigned long long page, bool PROGRESS = false){
67 80 //return binary<T>::read_plane_1(p, page);
68 81  
69   - unsigned int L = X() * sizeof(T); //caculate the number of bytes in a sample line
70   - unsigned int jump = X() * (Z() - 1) * sizeof(T);
  82 + if(PROGRESS) progress = 0;
  83 + unsigned long long L = X() * sizeof(T); //caculate the number of bytes in a sample line
  84 + unsigned long long jump = X() * (Z() - 1) * sizeof(T);
71 85  
72 86 if (page >= Z()){ //make sure the bank number is right
73 87 std::cout<<"ERROR: page out of range"<<std::endl;
... ... @@ -75,10 +89,11 @@ public:
75 89 }
76 90  
77 91 file.seekg(X() * page * sizeof(T), std::ios::beg);
78   - for (unsigned i = 0; i < Y(); i++)
  92 + for (unsigned long long i = 0; i < Y(); i++)
79 93 {
80 94 file.read((char *)(p + i * X()), L);
81 95 file.seekg( jump, std::ios::cur);
  96 + if(PROGRESS) progress = (double)(i+1) / Y() * 100;
82 97 }
83 98  
84 99 return true;
... ... @@ -88,23 +103,23 @@ public:
88 103  
89 104 /// @param p is a pointer to an allocated region of memory at least X * Y * sizeof(T) in size.
90 105 /// @param wavelength is a floating point value (usually a wavelength in spectral data) used as a label for the band to be copied.
91   - bool band( T * p, double wavelength){
  106 + bool band( T * p, double wavelength, bool PROGRESS = false){
92 107  
93 108 //if there are no wavelengths in the BSQ file
94 109 if(w.size() == 0)
95   - return band_index(p, (unsigned int)wavelength);
  110 + return band_index(p, (unsigned long long)wavelength);
96 111  
97   - unsigned int XY = X() * Y(); //calculate the number of pixels in a band
98   - unsigned int S = XY * sizeof(T); //calculate the number of bytes of a band
  112 + unsigned long long XY = X() * Y(); //calculate the number of pixels in a band
  113 + unsigned long long S = XY * sizeof(T); //calculate the number of bytes of a band
99 114  
100   - unsigned page=0; //bands around the wavelength
  115 + unsigned long long page=0; //bands around the wavelength
101 116  
102 117  
103 118 //get the bands numbers around the wavelength
104 119  
105 120 //if wavelength is smaller than the first one in header file
106 121 if ( w[page] > wavelength ){
107   - band_index(p, page);
  122 + band_index(p, page, PROGRESS);
108 123 return true;
109 124 }
110 125  
... ... @@ -124,17 +139,17 @@ public:
124 139 p1=(T*)malloc(S); //memory allocation
125 140 p2=(T*)malloc(S);
126 141 band_index(p1, page - 1);
127   - band_index(p2, page );
128   - for(unsigned i=0; i < XY; i++){
129   - double r = (double) (wavelength - w[page-1]) / (double) (w[page] - w[page-1]);
130   - p[i] = (p2[i] - p1[i]) * r + p1[i];
  142 + band_index(p2, page, PROGRESS);
  143 + for(unsigned long long i=0; i < XY; i++){
  144 + double r = (wavelength - w[page-1]) / (w[page] - w[page-1]);
  145 + p[i] = (T)(((double)p2[i] - (double)p1[i]) * r + (double)p1[i]);
131 146 }
132 147 free(p1);
133 148 free(p2);
134 149 }
135 150 else //if the wavelength is equal to a wavelength in header file
136 151 {
137   - band_index(p, page);
  152 + band_index(p, page, PROGRESS);
138 153 }
139 154  
140 155 return true;
... ... @@ -147,10 +162,10 @@ public:
147 162 /// @param wavelength is the wavelength to retrieve
148 163 bool read_x_from_xz(T* p, T* c, double wavelength)
149 164 {
150   - unsigned int B = Z();
151   - unsigned int L = X() * sizeof(T);
  165 + unsigned long long B = Z();
  166 + unsigned long long L = X() * sizeof(T);
152 167  
153   - unsigned page=0; //samples around the wavelength
  168 + unsigned long long page=0; //samples around the wavelength
154 169 T * p1;
155 170 T * p2;
156 171  
... ... @@ -179,9 +194,9 @@ public:
179 194 memcpy(p1, c + (page - 1) * X(), L);
180 195 memcpy(p2, c + page * X(), L);
181 196  
182   - for(unsigned i=0; i < X(); i++){
  197 + for(unsigned long long i=0; i < X(); i++){
183 198 double r = (double) (wavelength - w[page-1]) / (double) (w[page] - w[page-1]);
184   - p[i] = (p2[i] - p1[i]) * r + p1[i];
  199 + p[i] = (T)(((double)p2[i] - (double)p1[i]) * r + (double)p1[i]);
185 200 }
186 201 }
187 202 else //if the wavelength is equal to a wavelength in header file
... ... @@ -195,26 +210,26 @@ public:
195 210 /// @param p is a pointer to pre-allocated memory at least B * sizeof(T) in size.
196 211 /// @param x is the x-coordinate (dimension 1) of the spectrum.
197 212 /// @param y is the y-coordinate (dimension 2) of the spectrum.
198   - bool spectrum(T * p, unsigned x, unsigned y){
199   - return binary<T>::read_line_02(p, x, y);
  213 + bool spectrum(T * p, unsigned long long x, unsigned long long y, bool PROGRESS = false){
  214 + return binary<T>::read_line_1(p, x, y, PROGRESS);
200 215 }
201 216  
202 217 /// Retrieve a single pixel and stores it in pre-allocated memory.
203 218  
204 219 /// @param p is a pointer to pre-allocated memory at least sizeof(T) in size.
205 220 /// @param n is an integer index to the pixel using linear array indexing.
206   - bool pixel(T * p, unsigned n){
  221 + bool pixel(T * p, unsigned long long n){
207 222  
208 223 //calculate the corresponding x, y
209   - unsigned int x = n % X();
210   - unsigned int y = n / X();
  224 + unsigned long long x = n % X();
  225 + unsigned long long y = n / X();
211 226  
212 227 //get the pixel
213 228 return spectrum(p, x, y);
214 229 }
215 230  
216 231 //given a Y ,return a XZ slice
217   - bool read_plane_y(T * p, unsigned y){
  232 + bool read_plane_y(T * p, unsigned long long y){
218 233 return binary<T>::read_plane_2(p, y);
219 234 }
220 235  
... ... @@ -223,17 +238,17 @@ public:
223 238  
224 239 /// @param outname is the name of the output file used to store the resulting baseline-corrected data.
225 240 /// @param wls is the list of baseline points based on band labels.
226   - bool baseline(std::string outname, std::vector<double> wls){
  241 + bool baseline(std::string outname, std::vector<double> wls, unsigned char* mask = NULL, bool PROGRESS = false){
227 242  
228   - unsigned N = wls.size(); //get the number of baseline points
  243 + unsigned long long N = wls.size(); //get the number of baseline points
229 244  
230 245 std::ofstream target(outname.c_str(), std::ios::binary); //open the target binary file
231 246 std::string headername = outname + ".hdr"; //the header file name
232 247  
233 248 //simplify image resolution
234   - unsigned int ZX = Z() * X(); //calculate the number of points in a Y slice
235   - unsigned int L = ZX * sizeof(T); //calculate the number of bytes of a Y slice
236   - unsigned int B = Z();
  249 + unsigned long long ZX = Z() * X(); //calculate the number of points in a Y slice
  250 + unsigned long long L = ZX * sizeof(T); //calculate the number of bytes of a Y slice
  251 + unsigned long long B = Z();
237 252  
238 253 T* c; //pointer to the current Y slice
239 254 c = (T*)malloc(L); //memory allocation
... ... @@ -247,7 +262,7 @@ public:
247 262  
248 263 double ai, bi; //stores the two baseline points wavelength surrounding the current band
249 264 double ci; //stores the current band's wavelength
250   - unsigned control;
  265 + unsigned long long control;
251 266  
252 267 if (a == NULL || b == NULL || c == NULL){
253 268 std::cout<<"ERROR: error allocating memory";
... ... @@ -255,7 +270,7 @@ public:
255 270 }
256 271 // loop start correct every y slice
257 272  
258   - for (unsigned k =0; k < Y(); k++)
  273 + for (unsigned long long k =0; k < Y(); k++)
259 274 {
260 275 //get the current y slice
261 276 read_plane_y(c, k);
... ... @@ -280,7 +295,7 @@ public:
280 295  
281 296 //correct every YZ line
282 297  
283   - for(unsigned cii = 0; cii < B; cii++){
  298 + for(unsigned long long cii = 0; cii < B; cii++){
284 299  
285 300 //update baseline points, if necessary
286 301 if( w[cii] >= bi && cii != B - 1) {
... ... @@ -310,7 +325,7 @@ public:
310 325  
311 326 ci = w[cii];
312 327  
313   - unsigned jump = cii * X();
  328 + unsigned long long jump = cii * X();
314 329 //perform the baseline correction
315 330 for(unsigned i=0; i < X(); i++)
316 331 {
... ... @@ -321,7 +336,7 @@ public:
321 336 }//loop for YZ line end
322 337 target.write(reinterpret_cast<const char*>(c), L); //write the corrected data into destination
323 338  
324   - thread_data = (double)k / Y() * 100;
  339 + if(PROGRESS) progress = (double)(k+1) / Y() * 100;
325 340 }//loop for Y slice end
326 341  
327 342 free(a);
... ... @@ -329,8 +344,6 @@ public:
329 344 free(c);
330 345 target.close();
331 346  
332   - thread_data = 100;
333   -
334 347 return true;
335 348  
336 349 }
... ... @@ -340,13 +353,13 @@ public:
340 353 /// @param outname is the name of the output file used to store the resulting baseline-corrected data.
341 354 /// @param w is the label specifying the band that the hyperspectral image will be normalized to.
342 355 /// @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.
343   - bool normalize(std::string outname, double w, double t = 0.0)
  356 + bool normalize(std::string outname, double w, unsigned char* mask = NULL, bool PROGRESS = false)
344 357 {
345   - unsigned int B = Z(); //calculate the number of bands
346   - unsigned int ZX = Z() * X();
347   - unsigned int XY = X() * Y(); //calculate the number of pixels in a band
348   - unsigned int S = XY * sizeof(T); //calculate the number of bytes in a band
349   - unsigned int L = ZX * sizeof(T);
  358 + unsigned long long B = Z(); //calculate the number of bands
  359 + unsigned long long ZX = Z() * X();
  360 + unsigned long long XY = X() * Y(); //calculate the number of pixels in a band
  361 + unsigned long long S = XY * sizeof(T); //calculate the number of bytes in a band
  362 + unsigned long long L = ZX * sizeof(T);
350 363  
351 364 std::ofstream target(outname.c_str(), std::ios::binary); //open the target binary file
352 365 std::string headername = outname + ".hdr"; //the header file name
... ... @@ -359,14 +372,14 @@ public:
359 372  
360 373 band(b, w); //get the certain band into memory
361 374  
362   - for(unsigned j = 0; j < Y(); j++)
  375 + for(unsigned long long j = 0; j < Y(); j++)
363 376 {
364 377 read_plane_y(c, j);
365   - for(unsigned i = 0; i < B; i++)
  378 + for(unsigned long long i = 0; i < B; i++)
366 379 {
367   - for(unsigned m = 0; m < X(); m++)
  380 + for(unsigned long long m = 0; m < X(); m++)
368 381 {
369   - if( b[m + j * X()] < t )
  382 + if( mask != NULL && !mask[m + j * X()] )
370 383 c[m + i * X()] = (T)0.0;
371 384 else
372 385 c[m + i * X()] = c[m + i * X()] / b[m + j * X()];
... ... @@ -374,22 +387,21 @@ public:
374 387 }
375 388 target.write(reinterpret_cast<const char*>(c), L); //write normalized data into destination
376 389  
377   - thread_data = (double)j / Y() * 100;
  390 + if(PROGRESS) progress = (double)(j+1) / Y() * 100;
378 391 }
379 392  
380 393 free(b);
381 394 free(c);
382 395 target.close();
383   - thread_data = 100;
384 396 return true;
385 397 }
386 398  
387 399 /// Convert the current BIL file to a BSQ file with the specified file name.
388 400  
389 401 /// @param outname is the name of the output BSQ file to be saved to disk.
390   - bool bsq(std::string outname)
  402 + bool bsq(std::string outname, bool PROGRESS = false)
391 403 {
392   - unsigned int S = X() * Y() * sizeof(T); //calculate the number of bytes in a band
  404 + unsigned long long S = X() * Y() * sizeof(T); //calculate the number of bytes in a band
393 405  
394 406 std::ofstream target(outname.c_str(), std::ios::binary);
395 407 std::string headername = outname + ".hdr";
... ... @@ -397,16 +409,14 @@ public:
397 409 T * p; //pointer to the current band
398 410 p = (T*)malloc(S);
399 411  
400   - for ( unsigned i = 0; i < Z(); i++)
  412 + for ( unsigned long long i = 0; i < Z(); i++)
401 413 {
402 414 band_index(p, i);
403 415 target.write(reinterpret_cast<const char*>(p), S); //write a band data into target file
404 416  
405   - thread_data = (double)i / Z() * 100; //store the progress for the current operation
  417 + if(PROGRESS) progress = (double)(i+1) / Z() * 100; //store the progress for the current operation
406 418 }
407 419  
408   - thread_data = 100; //set the current progress to 100%
409   -
410 420 free(p);
411 421 target.close();
412 422 return true;
... ... @@ -415,9 +425,9 @@ public:
415 425 /// Convert the current BIL file to a BIP file with the specified file name.
416 426  
417 427 /// @param outname is the name of the output BIP file to be saved to disk.
418   - bool bip(std::string outname)
  428 + bool bip(std::string outname, bool PROGRESS = false)
419 429 {
420   - unsigned int S = X() * Z() * sizeof(T); //calculate the number of bytes in a ZX slice
  430 + unsigned long long S = X() * Z() * sizeof(T); //calculate the number of bytes in a ZX slice
421 431  
422 432 std::ofstream target(outname.c_str(), std::ios::binary);
423 433 std::string headername = outname + ".hdr";
... ... @@ -427,24 +437,21 @@ public:
427 437 T * q; //pointer to the current ZX slice for bip file
428 438 q = (T*)malloc(S);
429 439  
430   - for ( unsigned i = 0; i < Y(); i++)
  440 + for ( unsigned long long i = 0; i < Y(); i++)
431 441 {
432 442 read_plane_y(p, i);
433   - for ( unsigned k = 0; k < Z(); k++)
  443 + for ( unsigned long long k = 0; k < Z(); k++)
434 444 {
435   - unsigned ks = k * X();
436   - for ( unsigned j = 0; j < X(); j++)
  445 + unsigned long long ks = k * X();
  446 + for ( unsigned long long j = 0; j < X(); j++)
437 447 q[k + j * Z()] = p[ks + j];
438 448  
439   - thread_data = (double)(i * Z() + k) / (Z() * Y()) * 100; //store the progress for the current operation
  449 + if(PROGRESS) progress = (double)((i+1) * Z() + k+1) / (Z() * Y()) * 100; //store the progress for the current operation
440 450 }
441 451  
442 452 target.write(reinterpret_cast<const char*>(q), S); //write a band data into target file
443 453 }
444 454  
445   - thread_data = 100;
446   -
447   -
448 455 free(p);
449 456 free(q);
450 457 target.close();
... ... @@ -462,12 +469,12 @@ public:
462 469 /// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size.
463 470 bool baseline_band(double lb, double rb, T* lp, T* rp, double wavelength, T* result){
464 471  
465   - unsigned XY = X() * Y();
  472 + unsigned long long XY = X() * Y();
466 473 band(result, wavelength); //get band
467 474  
468 475 //perform the baseline correction
469 476 double r = (double) (wavelength - lb) / (double) (rb - lb);
470   - for(unsigned i=0; i < XY; i++){
  477 + for(unsigned long long i=0; i < XY; i++){
471 478 result[i] =(T) (result[i] - (rp[i] - lp[i]) * r - lp[i] );
472 479 }
473 480 return true;
... ... @@ -483,8 +490,8 @@ public:
483 490  
484 491 T* lp;
485 492 T* rp;
486   - unsigned XY = X() * Y();
487   - unsigned S = XY * sizeof(T);
  493 + unsigned long long XY = X() * Y();
  494 + unsigned long long S = XY * sizeof(T);
488 495 lp = (T*) malloc(S); //memory allocation
489 496 rp = (T*) malloc(S);
490 497  
... ... @@ -514,8 +521,8 @@ public:
514 521 T* cur; //current band 1
515 522 T* cur2; //current band 2
516 523  
517   - unsigned XY = X() * Y();
518   - unsigned S = XY * sizeof(T);
  524 + unsigned long long XY = X() * Y();
  525 + unsigned long long S = XY * sizeof(T);
519 526  
520 527 lp = (T*) malloc(S); //memory allocation
521 528 rp = (T*) malloc(S);
... ... @@ -525,9 +532,9 @@ public:
525 532 memset(result, (char)0, S);
526 533  
527 534 //find the wavelenght position in the whole band
528   - unsigned int n = w.size();
529   - unsigned int ai = 0; //left bound position
530   - unsigned int bi = n - 1; //right bound position
  535 + unsigned long long n = w.size();
  536 + unsigned long long ai = 0; //left bound position
  537 + unsigned long long bi = n - 1; //right bound position
531 538  
532 539  
533 540  
... ... @@ -556,23 +563,23 @@ public:
556 563 //calculate the beginning and the ending part
557 564 baseline_band(lb, rb, lp, rp, rab, cur2); //ending part
558 565 baseline_band(lb, rb, lp, rp, w[bi], cur);
559   - for(unsigned j = 0; j < XY; j++){
560   - result[j] += (rab - w[bi]) * (cur[j] + cur2[j]) / 2.0;
  566 + for(unsigned long long j = 0; j < XY; j++){
  567 + result[j] += (T)((rab - w[bi]) * ((double)cur[j] + (double)cur2[j]) / 2.0);
561 568 }
562 569 baseline_band(lb, rb, lp, rp, lab, cur2); //beginnning part
563 570 baseline_band(lb, rb, lp, rp, w[ai], cur);
564   - for(unsigned j = 0; j < XY; j++){
565   - result[j] += (w[ai] - lab) * (cur[j] + cur2[j]) / 2.0;
  571 + for(unsigned long long j = 0; j < XY; j++){
  572 + result[j] += (T)((w[ai] - lab) * ((double)cur[j] + (double)cur2[j]) / 2.0);
566 573 }
567 574  
568 575 //calculate the area
569 576 ai++;
570   - for(unsigned i = ai; i <= bi ;i++)
  577 + for(unsigned long long i = ai; i <= bi ;i++)
571 578 {
572 579 baseline_band(lb, rb, lp, rp, w[ai], cur2);
573   - for(unsigned j = 0; j < XY; j++)
  580 + for(unsigned long long j = 0; j < XY; j++)
574 581 {
575   - result[j] += (w[ai] - w[ai-1]) * (cur[j] + cur2[j]) / 2.0;
  582 + result[j] += (T)((w[ai] - w[ai-1]) * ((double)cur[j] + (double)cur2[j]) / 2.0);
576 583 }
577 584 std::swap(cur,cur2); //swap the band pointers
578 585 }
... ... @@ -602,7 +609,7 @@ public:
602 609 height(lb1, rb1, pos1, p1);
603 610 height(lb2, rb2, pos2, p2);
604 611 //calculate the ratio in result
605   - for(unsigned i = 0; i < X() * Y(); i++){
  612 + for(unsigned long long i = 0; i < X() * Y(); i++){
606 613 if(p1[i] == 0 && p2[i] ==0)
607 614 result[i] = 1;
608 615 else
... ... @@ -633,7 +640,7 @@ public:
633 640 area(lb1, rb1, lab1, rab1, p1);
634 641 height(lb2, rb2, pos, p2);
635 642 //calculate the ratio in result
636   - for(unsigned i = 0; i < X() * Y(); i++){
  643 + for(unsigned long long i = 0; i < X() * Y(); i++){
637 644 if(p1[i] == 0 && p2[i] ==0)
638 645 result[i] = 1;
639 646 else
... ... @@ -666,7 +673,7 @@ public:
666 673 area(lb1, rb1, lab1, rab1, p1);
667 674 area(lb2, rb2, lab2, rab2, p2);
668 675 //calculate the ratio in result
669   - for(unsigned i = 0; i < X() * Y(); i++){
  676 + for(unsigned long long i = 0; i < X() * Y(); i++){
670 677 if(p1[i] == 0 && p2[i] ==0)
671 678 result[i] = 1;
672 679 else
... ... @@ -691,8 +698,8 @@ public:
691 698 T* cur; //current band 1
692 699 T* cur2; //current band 2
693 700  
694   - unsigned XY = X() * Y();
695   - unsigned S = XY * sizeof(T);
  701 + unsigned long long XY = X() * Y();
  702 + unsigned long long S = XY * sizeof(T);
696 703  
697 704 lp = (T*) malloc(S); //memory allocation
698 705 rp = (T*) malloc(S);
... ... @@ -702,9 +709,9 @@ public:
702 709 memset(result, (char)0, S);
703 710  
704 711 //find the wavelenght position in the whole band
705   - unsigned int n = w.size();
706   - unsigned int ai = 0; //left bound position
707   - unsigned int bi = n - 1; //right bound position
  712 + unsigned long long n = w.size();
  713 + unsigned long long ai = 0; //left bound position
  714 + unsigned long long bi = n - 1; //right bound position
708 715  
709 716 //to make sure the left and the right bound are in the bandwidth
710 717 if (lb < w[0] || rb < w[0] || lb > w[n-1] || rb >w[n-1]){
... ... @@ -731,23 +738,23 @@ public:
731 738 //calculate the beginning and the ending part
732 739 baseline_band(lb, rb, lp, rp, rab, cur2); //ending part
733 740 baseline_band(lb, rb, lp, rp, w[bi], cur);
734   - for(unsigned j = 0; j < XY; j++){
735   - result[j] += (rab - w[bi]) * (rab + w[bi]) * (cur[j] + cur2[j]) / 4.0;
  741 + for(unsigned long long j = 0; j < XY; j++){
  742 + result[j] += (T)((rab - w[bi]) * (rab + w[bi]) * ((double)cur[j] + (double)cur2[j]) / 4.0);
736 743 }
737 744 baseline_band(lb, rb, lp, rp, lab, cur2); //beginnning part
738 745 baseline_band(lb, rb, lp, rp, w[ai], cur);
739   - for(unsigned j = 0; j < XY; j++){
740   - result[j] += (w[ai] - lab) * (w[ai] + lab) * (cur[j] + cur2[j]) / 4.0;
  746 + for(unsigned long long j = 0; j < XY; j++){
  747 + result[j] += (T)((w[ai] - lab) * (w[ai] + lab) * ((double)cur[j] + (double)cur2[j]) / 4.0);
741 748 }
742 749  
743 750 //calculate f(x) times x
744 751 ai++;
745   - for(unsigned i = ai; i <= bi ;i++)
  752 + for(unsigned long long i = ai; i <= bi ;i++)
746 753 {
747 754 baseline_band(lb, rb, lp, rp, w[ai], cur2);
748   - for(unsigned j = 0; j < XY; j++)
  755 + for(unsigned long long j = 0; j < XY; j++)
749 756 {
750   - result[j] += (w[ai] - w[ai-1]) * (w[ai] + w[ai-1]) * (cur[j] + cur2[j]) / 4.0;
  757 + result[j] += (T)((w[ai] - w[ai-1]) * (w[ai] + w[ai-1]) * ((double)cur[j] + (double)cur2[j]) / 4.0);
751 758 }
752 759 std::swap(cur,cur2); //swap the band pointers
753 760 }
... ... @@ -774,7 +781,7 @@ public:
774 781 x_area(lb, rb, lab, rab, p1);
775 782 area(lb, rb, lab, rab, p2);
776 783 //calculate the ratio in result
777   - for(unsigned i = 0; i < X() * Y(); i++){
  784 + for(unsigned long long i = 0; i < X() * Y(); i++){
778 785 if(p1[i] == 0 && p2[i] ==0)
779 786 result[i] = 1;
780 787 else
... ... @@ -793,16 +800,16 @@ public:
793 800 /// @param mask_band is the band used to specify the mask
794 801 /// @param threshold is the threshold used to determine if the mask value is true or false
795 802 /// @param p is a pointer to a pre-allocated array at least X * Y in size
796   - bool build_mask(double mask_band, double threshold, unsigned char* p){
  803 + bool build_mask(double mask_band, double threshold, unsigned char* p, bool PROGRESS = false){
797 804  
798 805 T* temp = (T*)malloc(X() * Y() * sizeof(T)); //allocate memory for the certain band
799   - band(temp, mask_band);
  806 + band(temp, mask_band, PROGRESS);
800 807  
801   - for (unsigned int i = 0; i < X() * Y(); i++) {
802   - if (temp[i] < threshold)
803   - p[i] = 0;
804   - else
805   - p[i] = 255;
  808 + for (unsigned long long i = 0; i < X() * Y(); i++) {
  809 + if (temp[i] < threshold)
  810 + p[i] = 0;
  811 + else
  812 + p[i] = 255;
806 813 }
807 814  
808 815 free(temp);
... ... @@ -818,17 +825,17 @@ public:
818 825 std::ofstream target(outfile.c_str(), std::ios::binary);
819 826  
820 827 //I THINK THIS IS WRONG
821   - unsigned XZ = X() * Z(); //calculate the number of values in a page on disk
822   - unsigned L = XZ * sizeof(T); //calculate the size of the page (in bytes)
  828 + unsigned long long XZ = X() * Z(); //calculate the number of values in a page on disk
  829 + unsigned long long L = XZ * sizeof(T); //calculate the size of the page (in bytes)
823 830  
824 831 T * temp = (T*)malloc(L); //allocate memory for a temporary page
825 832  
826   - for (unsigned i = 0; i < Y(); i++) //for each value in Y() (BIP should be X)
  833 + for (unsigned long long i = 0; i < Y(); i++) //for each value in Y() (BIP should be X)
827 834 {
828 835 read_plane_y(temp, i); //retrieve an ZX slice, stored in temp
829   - for ( unsigned j = 0; j < Z(); j++) //for each Z() (Y)
  836 + for ( unsigned long long j = 0; j < Z(); j++) //for each Z() (Y)
830 837 {
831   - for (unsigned k = 0; k < X(); k++) //for each band
  838 + for (unsigned long long k = 0; k < X(); k++) //for each band
832 839 {
833 840 if(p[i * X() + k] == 0)
834 841 temp[j * X() + k] = 0;
... ... @@ -845,14 +852,14 @@ public:
845 852  
846 853 /// Saves to disk only those spectra corresponding to mask values != 0
847 854 /// Unlike the BIP and BSQ versions of this function, this version saves a different format: the BIL file is saved as a BIP
848   - bool sift(std::string outfile, unsigned char* p){
  855 + bool sift(std::string outfile, unsigned char* p, bool PROGRESS = false){
849 856 // Assume X() = X, Y() = Y, Z() = Z.
850 857 std::ofstream target(outfile.c_str(), std::ios::binary);
851 858  
852 859 //for loading pages:
853   - unsigned long XZ = X() * Z(); //calculate the number of values in an XZ page on disk
854   - unsigned long B = Z(); //calculate the number of bands
855   - unsigned long L = XZ * sizeof(T); //calculate the size of the page (in bytes)
  860 + unsigned long long XZ = X() * Z(); //calculate the number of values in an XZ page on disk
  861 + unsigned long long B = Z(); //calculate the number of bands
  862 + unsigned long long L = XZ * sizeof(T); //calculate the size of the page (in bytes)
856 863  
857 864 //allocate temporary memory for a XZ slice
858 865 T* slice = (T*) malloc(L);
... ... @@ -861,32 +868,31 @@ public:
861 868 T* spec = (T*) malloc(B * sizeof(T));
862 869  
863 870 //for each slice along the y axis
864   - for (unsigned long y = 0; y < Y(); y++) //Select a page by choosing Y coordinate, Y()
  871 + for (unsigned long long y = 0; y < Y(); y++) //Select a page by choosing Y coordinate, Y()
865 872 {
866 873 read_plane_y(slice, y); //retrieve an ZX page, store in "slice"
867 874  
868 875 //for each sample along X
869   - for (unsigned long x = 0; x < X(); x++) //Select a pixel by choosing X coordinate in the page, X()
  876 + for (unsigned long long x = 0; x < X(); x++) //Select a pixel by choosing X coordinate in the page, X()
870 877 {
871 878 //if the mask != 0 at that xy pixel
872 879 if (p[y * X() + x] != 0) //if the mask != 0 at that XY pixel
873 880 {
874 881 //for each band at that pixel
875   - for (unsigned long b = 0; b < B; b++) //Select a voxel by choosing Z coordinate at the pixel
  882 + for (unsigned long long b = 0; b < B; b++) //Select a voxel by choosing Z coordinate at the pixel
876 883 {
877 884 spec[b] = slice[b*X() + x]; //Pass the correct spectral value from XZ page into the spectrum to be saved.
878 885 }
879 886 target.write((char*)spec, B * sizeof(T)); //write that spectrum to disk. Size is L2.
880 887 }
881 888  
882   - thread_data = (double) (y * X() + x) / (Y() * X()) * 100;
  889 + if(PROGRESS) progress = (double) ((y+1) * X() + x+1) / (Y() * X()) * 100;
883 890 }
884 891 }
885 892 target.close();
886 893 free(slice);
887 894 free(spec);
888 895  
889   - thread_data = 100;
890 896 return true;
891 897 }
892 898  
... ... @@ -898,20 +904,20 @@ public:
898 904 T* temp = (T*)malloc(sizeof(T) * XZ);
899 905 T* line = (T*)malloc(sizeof(T) * X());
900 906  
901   - for (unsigned i = 0; i < Y(); i++){
  907 + for (unsigned long long i = 0; i < Y(); i++){
902 908 getY(temp, i);
903 909 //initialize x-line
904   - for (unsigned j = 0; j < X(); j++){
  910 + for (unsigned long long j = 0; j < X(); j++){
905 911 line[j] = 0;
906 912 }
907   - unsigned c = 0;
908   - for (unsigned j = 0; j < Z(); j++){
909   - for (unsigned k = 0; k < X(); k++){
  913 + unsigned long long c = 0;
  914 + for (unsigned long long j = 0; j < Z(); j++){
  915 + for (unsigned long long k = 0; k < X(); k++){
910 916 line[k] += temp[c] / (T)Z();
911 917 c++;
912 918 }
913 919 }
914   - for (unsigned j = 0; j < X(); j++){
  920 + for (unsigned long long j = 0; j < X(); j++){
915 921 p[j + i * X()] = line[j];
916 922 }
917 923 }
... ... @@ -923,11 +929,11 @@ public:
923 929  
924 930 /// @param p is a pointer to pre-allocated memory of size [B * sizeof(T)] that stores the mean spectrum
925 931 /// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
926   - bool avg_band(double* p, unsigned char* mask = NULL){
  932 + bool avg_band(double* p, unsigned char* mask = NULL, bool PROGRESS = false){
927 933 unsigned long long XZ = X() * Z();
928 934 unsigned long long XY = X() * Y();
929 935 T* temp = (T*)malloc(sizeof(T) * XZ);
930   - for (unsigned j = 0; j < Z(); j++){
  936 + for (unsigned long long j = 0; j < Z(); j++){
931 937 p[j] = 0;
932 938 }
933 939 //calculate vaild number in a band
... ... @@ -947,6 +953,7 @@ public:
947 953 }
948 954 }
949 955 }
  956 + if(PROGRESS) progress = (double)(k+1) / Y() * 100;
950 957 }
951 958 free(temp);
952 959 return true;
... ... @@ -957,8 +964,8 @@ public:
957 964 /// @param co is a pointer to pre-allocated memory of size [B * B] that stores the resulting covariance matrix
958 965 /// @param avg is a pointer to memory of size B that stores the average spectrum
959 966 /// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
960   - bool co_matrix(double* co, double* avg, unsigned char *mask){
961   - thread_data = 0;
  967 + bool co_matrix(double* co, double* avg, unsigned char *mask, bool PROGRESS = false){
  968 + progress = 0;
962 969 //memory allocation
963 970 unsigned long long xy = X() * Y();
964 971 unsigned long long B = Z();
... ... @@ -986,7 +993,7 @@ public:
986 993 }
987 994 }
988 995 }
989   - thread_data = (double)j / xy * 100;
  996 + if(PROGRESS) progress = (double)(j+1) / xy * 100;
990 997 }
991 998 //because correlation matrix is symmetric
992 999 for (unsigned long long i = 0; i < B; i++){
... ... @@ -996,9 +1003,6 @@ public:
996 1003 }
997 1004  
998 1005 free(temp);
999   -
1000   - thread_data = 100; //processing complete
1001   -
1002 1006 return true;
1003 1007 }
1004 1008  
... ... @@ -1015,7 +1019,8 @@ public:
1015 1019 unsigned long long x1,
1016 1020 unsigned long long y1,
1017 1021 unsigned long long b0,
1018   - unsigned long long b1){
  1022 + unsigned long long b1,
  1023 + bool PROGRESS = false){
1019 1024  
1020 1025 //calculate the new image parameters
1021 1026 unsigned long long samples = x1 - x0;
... ... @@ -1041,14 +1046,14 @@ public:
1041 1046 //set the start position for the cropped region
1042 1047 file.seekg((y0 * X() * Z() + b0 * X() + x0) * sizeof(T), std::ios::beg);
1043 1048  
1044   - for (unsigned x = 0; x < lines; x++)
  1049 + for (unsigned long long x = 0; x < lines; x++)
1045 1050 {
1046   - for (unsigned z = b0; z < b1; z++)
  1051 + for (unsigned long long z = b0; z < b1; z++)
1047 1052 {
1048 1053 file.read((char *)(temp + z * samples), sizeof(T) * samples);
1049 1054 file.seekg(jumpb, std::ios::cur); //go to the next band
1050 1055  
1051   - thread_data = (double)(x * Z() + z) / (lines * Z()) * 100;
  1056 + if(PROGRESS) progress = (double)(x * Z() + z+1) / (lines * Z()) * 100;
1052 1057 }
1053 1058  
1054 1059 //write slice data into target file
... ... @@ -1061,8 +1066,6 @@ public:
1061 1066 //free the temporary frame
1062 1067 free(temp);
1063 1068  
1064   - thread_data = 100;
1065   -
1066 1069 return true;
1067 1070 }
1068 1071  
... ...
stim/envi/binary.h
  Show/Hide comments   View file @d4f06a7
... ... @@ -24,20 +24,20 @@ protected:
24 24 std::fstream file; //file stream used for reading and writing
25 25 std::string name; //file name
26 26  
27   - unsigned long long int R[D]; //resolution
28   - unsigned int header; //header size (in bytes)
  27 + unsigned long long R[D]; //resolution
  28 + unsigned long long header; //header size (in bytes)
29 29 unsigned char* mask; //pointer to a character array: 0 = background, 1 = foreground (or valid data)
30 30  
31   - unsigned int thread_data; //unsigned integer used to pass data to threads during processing
  31 + double progress; //stores the progress on the current operation (accessible using a thread)
32 32  
33 33  
34 34 /// Private initialization function used to set default parameters in the data structure.
35 35 void init(){
36   - memset(R, 0, sizeof(unsigned int) * D); //initialize the resolution to zero
  36 + memset(R, 0, sizeof(unsigned long long) * D); //initialize the resolution to zero
37 37 header = 0; //initialize the header size to zero
38 38 mask = NULL;
39 39  
40   - thread_data = 0;
  40 + progress = 0;
41 41 }
42 42  
43 43 /// Private helper function that returns the size of the file on disk using system functions.
... ... @@ -105,12 +105,12 @@ protected:
105 105  
106 106 public:
107 107  
108   - unsigned int get_thread_data(){
109   - return thread_data;
  108 + double get_progress(){
  109 + return progress;
110 110 }
111 111  
112   - void reset_thread_data(){
113   - thread_data = 0;
  112 + void reset_progress(){
  113 + progress = 0;
114 114 }
115 115  
116 116 /// Open a binary file for streaming.
... ... @@ -118,9 +118,9 @@ public:
118 118 /// @param filename is the name of the binary file
119 119 /// @param r is a STIM vector specifying the size of the binary file along each dimension
120 120 /// @param h is the length (in bytes) of any header file (default zero)
121   - bool open(std::string filename, vec<unsigned int> r, unsigned int h = 0){
  121 + bool open(std::string filename, vec<unsigned long long> r, unsigned long long h = 0){
122 122  
123   - for(unsigned int i = 0; i < D; i++) //set the dimensions of the binary file object
  123 + for(unsigned long long i = 0; i < D; i++) //set the dimensions of the binary file object
124 124 R[i] = r[i];
125 125  
126 126 header = h; //save the header size
... ... @@ -137,17 +137,17 @@ public:
137 137 /// @param filename is the name of the binary file to be created
138 138 /// @param r is a STIM vector specifying the size of the file along each dimension
139 139 /// @offset specifies how many bytes to offset the file (used to leave room for a header)
140   - bool create(std::string filename, vec<unsigned int> r, unsigned int offset = 0){
  140 + bool create(std::string filename, vec<unsigned long long> r, unsigned long long offset = 0){
141 141  
142 142 std::ofstream target(filename.c_str(), std::ios::binary);
143 143  
144 144 //initialize binary file
145 145 T p = 0;
146   - for(unsigned int i =0; i < r[0] * r[1] * r[2]; i++){
  146 + for(unsigned long long i =0; i < r[0] * r[1] * r[2]; i++){
147 147 target.write((char*)(&p), sizeof(T));
148 148 }
149 149  
150   - for(unsigned int i = 0; i < D; i++) //set the dimensions of the binary file object
  150 + for(unsigned long long i = 0; i < D; i++) //set the dimensions of the binary file object
151 151 R[i] = r[i];
152 152  
153 153 header = offset; //save the header size
... ... @@ -161,7 +161,7 @@ public:
161 161  
162 162 /// @param p is a pointer to the data to be written
163 163 /// @param page is the page number (index of the highest-numbered dimension)
164   - bool write_page( T * p, unsigned int page){
  164 + bool write_page( T * p, unsigned long long page){
165 165  
166 166 if(p == NULL){
167 167 std::cout<<"ERROR: unable to write into file, empty pointer"<<std::endl;
... ... @@ -178,7 +178,9 @@ public:
178 178  
179 179 /// @param p is a pointer to pre-allocated memory equal to the page size
180 180 /// @param page is the index of the page
181   - bool read_page( T * p, unsigned int page){
  181 + bool read_page( T * p, unsigned long long page, bool PROGRESS = false){
  182 +
  183 + if(PROGRESS) progress = 0;
182 184  
183 185 if (page >= R[2]){ //make sure the bank number is right
184 186 std::cout<<"ERROR: page out of range"<<std::endl;
... ... @@ -187,7 +189,7 @@ public:
187 189  
188 190 file.seekg(R[1] * R[0] * page * sizeof(T) + header, std::ios::beg); //write into memory from the binary file
189 191 file.read((char *)p, R[0] * R[1] * sizeof(T));
190   -
  192 + if(PROGRESS) progress = 100;
191 193 return true;
192 194 }
193 195  
... ... @@ -198,8 +200,10 @@ public:
198 200 /// @param p is a pointer to pre-allocated memory equal to the line size R[2]
199 201 /// @param x is the x coordinate
200 202 /// @param y is the y coordinate
201   - bool read_line_01( T * p, unsigned int x, unsigned int y){
202   - unsigned int i;
  203 + bool read_line_2( T * p, unsigned long long x, unsigned long long y, bool PROGRESS = false){
  204 + unsigned long long i;
  205 +
  206 + if(PROGRESS) progress = 0;
203 207  
204 208 if ( x >= R[0] || y >= R[1]){ //make sure the sample and line number is right
205 209 std::cout<<"ERROR: sample or line out of range"<<std::endl;
... ... @@ -211,7 +215,9 @@ public:
211 215 {
212 216 file.read((char *)(p + i), sizeof(T));
213 217 file.seekg((R[1] * R[0] - 1) * sizeof(T), std::ios::cur); //go to the next band
  218 + if(PROGRESS) progress = (double)i / (double)R[2] * 100;
214 219 }
  220 + if(PROGRESS) progress = 100;
215 221  
216 222 return true;
217 223 }
... ... @@ -221,7 +227,7 @@ public:
221 227 /// @param p is a pointer to pre-allocated memory equal to the line size R[2]
222 228 /// @param x is the y coordinate
223 229 /// @param y is the z coordinate
224   - bool read_line_12(T * p, unsigned int y, unsigned int z){
  230 + bool read_line_0(T * p, unsigned long long y, unsigned long long z, bool PROGRESS = false){
225 231 //test to make sure the specified value is within range
226 232 if( y >= R[1] || z >= R[2] ){
227 233 std::cout<<"ERROR: sample or line out of range"<<std::endl;
... ... @@ -230,7 +236,7 @@ public:
230 236  
231 237 file.seekg((z * R[0] * R[1] + y * R[0]) * sizeof(T), std::ios::beg); //seek to the start of the line
232 238 file.read((char *)p, sizeof(T) * R[0]); //read the line
233   -
  239 + if(PROGRESS) progress = 100;
234 240 return true;
235 241 }
236 242  
... ... @@ -239,7 +245,8 @@ public:
239 245 /// @param p is a pointer to pre-allocated memory equal to the line size R[2]
240 246 /// @param x is the y coordinate
241 247 /// @param z is the z coordinate
242   - bool read_line_02(T * p, unsigned int x, unsigned int z){
  248 + bool read_line_1(T * p, unsigned long long x, unsigned long long z, bool PROGRESS = false){
  249 + if(PROGRESS) progress = 0;
243 250 //test to make sure the specified value is within range
244 251 if( x >= R[0] || z >= R[2] ){
245 252 std::cout<<"ERROR: sample or line out of range"<<std::endl;
... ... @@ -247,11 +254,12 @@ public:
247 254 }
248 255  
249 256 file.seekg((z * R[0] * R[1] + x) * sizeof(T), std::ios::beg); //seek to the start of the line
250   - for (unsigned int i = 0; i < R[1]; i++){ //for each pixel in the line
  257 + for (unsigned long long i = 0; i < R[1]; i++){ //for each pixel in the line
251 258 file.read((char *)(p + i), sizeof(T)); //read the pixel
252 259 file.seekg((R[0] - 1) * sizeof(T), std::ios::cur); //seek to the next pixel in the line
  260 + if(PROGRESS) progress = (double)i / (double)R[1] * 100;
253 261 }
254   -
  262 + if(PROGRESS) progress = 100;
255 263 return true;
256 264 }
257 265  
... ... @@ -259,23 +267,25 @@ public:
259 267  
260 268 /// @param p is a pointer to pre-allocated memory of size R[1] * R[2] * sizeof(T)
261 269 /// @param n is the 0-axis coordinate used to retrieve the plane
262   - bool read_plane_0(T* p, unsigned int n){
263   -
  270 + bool read_plane_0(T* p, unsigned long long n, bool PROGRESS = false){
  271 + if(PROGRESS) progress = 0;
264 272 if (n >= R[0]){ //make sure the number is within the possible range
265 273 std::cout<<"ERROR read_plane_0: page out of range"<<std::endl;
266 274 return false;
267 275 }
268   - unsigned int jump = (R[0] - 1) * sizeof(T); //number of bytes to skip between samples
  276 + unsigned long long jump = (R[0] - 1) * sizeof(T); //number of bytes to skip between samples
269 277  
270 278 //seek to the start of the plane
271 279 file.seekg(n * sizeof(T), std::ios::beg);
272 280  
273   - unsigned int N = R[1] * R[2];
274   - for(unsigned int i = 0; i<N; i++){
  281 + unsigned long long N = R[1] * R[2];
  282 + for(unsigned long long i = 0; i<N; i++){
275 283 file.read((char*)(p+i), sizeof(T));
276 284 file.seekg(jump, std::ios::cur);
  285 + if(PROGRESS) progress = (double)(i+1) / N * 100;
277 286 }
278 287  
  288 + if(PROGRESS) progress = 100;
279 289 return true;
280 290  
281 291  
... ... @@ -285,10 +295,10 @@ public:
285 295  
286 296 /// @param p is a pointer to pre-allocated memory of size R[0] * R[2] * sizeof(T)
287 297 /// @param n is the 1-axis coordinate used to retrieve the plane
288   - bool read_plane_1(T* p, unsigned int n){
289   -
290   - unsigned int L = R[0] * sizeof(T); //caculate the number of bytes in a sample line
291   - unsigned int jump = R[0] * (R[1] - 1) * sizeof(T);
  298 + bool read_plane_1(T* p, unsigned long long n, bool PROGRESS = false){
  299 + if(PROGRESS) progress = 0;
  300 + unsigned long long L = R[0] * sizeof(T); //caculate the number of bytes in a sample line
  301 + unsigned long long jump = R[0] * (R[1] - 1) * sizeof(T);
292 302  
293 303 if (n >= R[1]){ //make sure the bank number is right
294 304 std::cout<<"ERROR read_plane_1: page out of range"<<std::endl;
... ... @@ -296,13 +306,14 @@ public:
296 306 }
297 307  
298 308 file.seekg(R[0] * n * sizeof(T), std::ios::beg);
299   - for (unsigned i = 0; i < R[2]; i++){
300   -
  309 + for (unsigned long long i = 0; i < R[2]; i++){
  310 + if(PROGRESS) progress = (double)i / R[2] * 100;
301 311 file.read((char *)(p + i * R[0]), L);
302 312 file.seekg( jump, std::ios::cur);
303 313 std::cout<<i<<" ";
304 314 }
305 315  
  316 + if(PROGRESS) progress = 100;
306 317 return true;
307 318 }
308 319  
... ... @@ -310,15 +321,15 @@ public:
310 321  
311 322 /// @param p is a pointer to pre-allocated memory of size R[0] * R[1] * sizeof(T)
312 323 /// @param n is the 2-axis coordinate used to retrieve the plane
313   - bool read_plane_2(T* p, unsigned int n){
314   - return read_page(p, n);
  324 + bool read_plane_2(T* p, unsigned long long n, bool PROGRESS = false){
  325 + return read_page(p, n, PROGRESS);
315 326 }
316 327  
317 328 /// Reads a single pixel, treating the entire data set as a linear array
318 329  
319 330 /// @param p is a pointer to pre-allocated memory of size sizeof(T)
320 331 /// @param i is the index to the pixel using linear indexing
321   - bool read_pixel(T* p, unsigned int i){
  332 + bool read_pixel(T* p, unsigned long long i){
322 333 if(i >= R[0] * R[1] * R[2]){
323 334 std::cout<<"ERROR read_pixel: n is out of range"<<std::endl;
324 335 return false;
... ... @@ -335,14 +346,14 @@ public:
335 346 /// @param x is the x (0) axis coordinate
336 347 /// @param y is the y (1) axis coordinate
337 348 /// @param z is the z (2) axis coordinate
338   - bool read_pixel(T* p, unsigned int x, unsigned int y, unsigned int z){
  349 + bool read_pixel(T* p, unsigned long long x, unsigned long long y, unsigned long long z){
339 350  
340 351 if(x < 0 || x >= R[0] || y < 0 || y >= R[1] || z < 0 || z > R[2]){
341 352 std::cout<<"ERROR read_pixel: (x,y,z) is out of range"<<std::endl;
342 353 return false;
343 354 }
344 355  
345   - unsigned int i = z * R[0] * R[1] + y * R[0] + z;
  356 + unsigned long long i = z * R[0] * R[1] + y * R[0] + z;
346 357 return read_pixel(p, i);
347 358 }
348 359  
... ...
stim/envi/bip.h
  Show/Hide comments   View file @d4f06a7
... ... @@ -3,10 +3,14 @@
3 3  
4 4 #include "../envi/envi_header.h"
5 5 #include "../envi/bil.h"
6   -#include "../envi/binary.h"
  6 +#include "../envi/hsi.h"
7 7 #include <cstring>
8 8 #include <utility>
9 9  
  10 +//CUDA
  11 +#include <cuda_runtime.h>
  12 +#include "cublas_v2.h"
  13 +
10 14 namespace stim{
11 15  
12 16 /**
... ... @@ -19,32 +23,98 @@ namespace stim{
19 23 */
20 24 template <typename T>
21 25  
22   -class bip: public binary<T> {
  26 +class bip: public hsi<T> {
23 27  
24 28 protected:
25 29  
26 30  
27   - std::vector<double> w; //band wavelength
28   - unsigned int offset; //header offset
  31 + //std::vector<double> w; //band wavelength
  32 + unsigned long long offset; //header offset
29 33  
30   - unsigned long X(){
  34 + /*unsigned long long X(){
31 35 return R[1];
32 36 }
33   - unsigned long Y(){
  37 + unsigned long long Y(){
34 38 return R[2];
35 39 }
36   - unsigned long Z(){
  40 + unsigned long long Z(){
37 41 return R[0];
  42 + }*/
  43 + using hsi<T>::w; //use the wavelength array in stim::hsi
  44 + using hsi<T>::nnz;
  45 + using binary<T>::progress;
  46 +
  47 + /// Call the binary nnz() function for the BIP orientation
  48 + //unsigned long long nnz(unsigned char* mask){
  49 + // return binary<T>::nnz(mask, X()*Y());
  50 + //}
  51 +
  52 + /// Linear interpolation of a spectral value given the bounding spectral samples
  53 + /*T lerp(double w, T low_v, double low_w, T high_v, double high_w){
  54 + if(low_w == high_w) return low_v; //if the interval is of zero length, just return one of the bounds
  55 + double alpha = (w - low_w) / (high_w - low_w); //calculate the interpolation factor
  56 + return (1.0 - alpha) * low_v + alpha * high_v; //interpolate
  57 + }
  58 +
  59 + /// Gets the two band indices surrounding a given wavelength
  60 + void band_bounds(double wavelength, unsigned long long& low, unsigned long long& high){
  61 + unsigned long long B = Z();
  62 + for(high = 0; high < B; high++){
  63 + if(w[high] > wavelength) break;
  64 + }
  65 + low = 0;
  66 + if(high > 0)
  67 + low = high-1;
  68 + }
  69 +
  70 + /// Get the list of band numbers that bound a list of wavelengths
  71 + void band_bounds(std::vector<double> wavelengths,
  72 + std::vector<unsigned long long>& low_bands,
  73 + std::vector<unsigned long long>& high_bands){
  74 +
  75 + unsigned long long W = w.size(); //get the number of wavelengths in the list
  76 + low_bands.resize(W); //pre-allocate space for the band lists
  77 + high_bands.resize(W);
  78 +
  79 + for(unsigned long long wl = 0; wl < W; wl++){ //for each wavelength
  80 + band_bounds(wavelengths[wl], low_bands[wl], high_bands[wl]); //find the low and high bands
  81 + }
38 82 }
39 83  
40   - using binary<T>::thread_data;
  84 + /// Returns the interpolated in the given spectrum based on the given wavelength
  85 +
  86 + /// @param s is the spectrum in main memory of length Z()
  87 + /// @param wavelength is the wavelength value to interpolate out
  88 + T interp_spectrum(T* s, double wavelength){
  89 + unsigned long long low, high; //indices for the bands surrounding wavelength
  90 + band_bounds(wavelength, low, high); //get the surrounding band indices
  91 +
  92 + if(high == w.size()) return s[w.size()-1]; //if the high band is above the wavelength range, return the highest wavelength value
  93 +
  94 + return lerp(wavelength, s[low], w[low], s[high], w[high]);
  95 + }
  96 +
  97 + /// Returns the interpolated value corresponding to the given list of wavelengths
  98 + std::vector<T> interp_spectrum(T* s, std::vector<double> wavelengths){
  99 +
  100 + unsigned long long N = wavelengths.size(); //get the number of wavelength measurements
  101 +
  102 + std::vector<T> v; //allocate space for the resulting values
  103 + v.resize(wavelengths.size());
  104 + for(unsigned long long n = 0; n < N; n++){ //for each measurement
  105 + v[n] = interp_spectrum(s, wavelengths[n]); //interpolate the measurement
  106 + }
  107 + return v;
  108 + }*/
41 109  
42 110 public:
43 111  
44 112 using binary<T>::open;
45 113 using binary<T>::file;
46 114 using binary<T>::R;
47   - using binary<T>::read_line_12;
  115 + using binary<T>::read_line_0;
  116 +
  117 + bip(){ hsi<T>::init_bip(); }
48 118  
49 119 /// Open a data file for reading using the class interface.
50 120  
... ... @@ -54,14 +124,19 @@ public:
54 124 /// @param B is the number of samples (bands) along dimension 3
55 125 /// @param header_offset is the number of bytes (if any) in the binary header
56 126 /// @param wavelengths is an optional STL vector of size B specifying a numerical label for each band
57   - bool open(std::string filename, unsigned int X, unsigned int Y, unsigned int B, unsigned int header_offset, std::vector<double> wavelengths){
  127 + bool open(std::string filename,
  128 + unsigned long long X,
  129 + unsigned long long Y,
  130 + unsigned long long B,
  131 + unsigned long long header_offset,
  132 + std::vector<double> wavelengths){
58 133  
59 134 //copy the wavelengths to the BSQ file structure
60 135 w = wavelengths;
61 136 //copy the offset to the structure
62 137 offset = header_offset;
63 138  
64   - return open(filename, vec<unsigned int>(B, X, Y), header_offset);
  139 + return open(filename, vec<unsigned long long>(B, X, Y), header_offset);
65 140  
66 141 }
67 142  
... ... @@ -69,21 +144,21 @@ public:
69 144  
70 145 /// @param p is a pointer to an allocated region of memory at least X * Y * sizeof(T) in size.
71 146 /// @param page <= B is the integer number of the band to be copied.
72   - bool band_index( T * p, unsigned int page){
73   - return binary<T>::read_plane_0(p, page);
  147 + bool band_index( T * p, unsigned long long page, bool PROGRESS = false){
  148 + return binary<T>::read_plane_0(p, page, PROGRESS);
74 149 }
75 150  
76 151 /// Retrieve a single band (by numerical label) and stores it in pre-allocated memory.
77 152  
78 153 /// @param p is a pointer to an allocated region of memory at least X * Y * sizeof(T) in size.
79 154 /// @param wavelength is a floating point value (usually a wavelength in spectral data) used as a label for the band to be copied.
80   - bool band( T * p, double wavelength){
  155 + bool band( T * p, double wavelength, bool PROGRESS = false){
81 156  
82 157 //if there are no wavelengths in the BSQ file
83 158 if(w.size() == 0)
84   - return band_index(p, (unsigned int)wavelength);
  159 + return band_index(p, (unsigned long long)wavelength, PROGRESS);
85 160  
86   - unsigned int XY = X() * Y(); //calculate the number of pixels in a band
  161 + unsigned long long XY = X() * Y(); //calculate the number of pixels in a band
87 162  
88 163 unsigned page=0; //bands around the wavelength
89 164  
... ... @@ -92,7 +167,7 @@ public:
92 167  
93 168 //if wavelength is smaller than the first one in header file
94 169 if ( w[page] > wavelength ){
95   - band_index(p, page);
  170 + band_index(p, page, PROGRESS);
96 171 return true;
97 172 }
98 173  
... ... @@ -101,7 +176,7 @@ public:
101 176 page++;
102 177 //if wavelength is larger than the last wavelength in header file
103 178 if (page == Z()) {
104   - band_index(p, Z()-1);
  179 + band_index(p, Z()-1, PROGRESS);
105 180 return true;
106 181 }
107 182 }
... ... @@ -112,17 +187,17 @@ public:
112 187 p1=(T*)malloc( XY * sizeof(T)); //memory allocation
113 188 p2=(T*)malloc( XY * sizeof(T));
114 189 band_index(p1, page - 1);
115   - band_index(p2, page );
116   - for(unsigned i=0; i < XY; i++){
  190 + band_index(p2, page, PROGRESS);
  191 + for(unsigned long long i=0; i < XY; i++){
117 192 double r = (double) (wavelength - w[page-1]) / (double) (w[page] - w[page-1]);
118   - p[i] = (p2[i] - p1[i]) * r + p1[i];
  193 + p[i] = (T)(((double)p2[i] - (double)p1[i]) * r + (double)p1[i]);
119 194 }
120 195 free(p1);
121 196 free(p2);
122 197 }
123 198 else //if the wavelength is equal to a wavelength in header file
124 199 {
125   - band_index(p, page);
  200 + band_index(p, page, PROGRESS);
126 201 }
127 202  
128 203 return true;
... ... @@ -133,8 +208,8 @@ public:
133 208 /// @param p is a pointer to pre-allocated memory at least B * sizeof(T) in size.
134 209 /// @param x is the x-coordinate (dimension 1) of the spectrum.
135 210 /// @param y is the y-coordinate (dimension 2) of the spectrum.
136   - bool spectrum(T * p, unsigned x, unsigned y){
137   - return read_line_12(p, x, y); //read a line in the binary YZ plane (dimension order for BIP is ZXY)
  211 + bool spectrum(T * p, unsigned long long x, unsigned long long y, bool PROGRESS = false){
  212 + return read_line_0(p, x, y, PROGRESS); //read a line in the binary YZ plane (dimension order for BIP is ZXY)
138 213 }
139 214  
140 215 /// Retrieves a band of x values from a given xz plane.
... ... @@ -144,16 +219,16 @@ public:
144 219 /// @param wavelength is the wavelength of X values to retrieve
145 220 bool read_x_from_xz(T* p, T* c, double wavelength)
146 221 {
147   - unsigned int B = Z();
  222 + unsigned long long B = Z();
148 223  
149   - unsigned page=0; //samples around the wavelength
  224 + unsigned long long page=0; //samples around the wavelength
150 225  
151 226  
152 227 //get the bands numbers around the wavelength
153 228  
154 229 //if wavelength is smaller than the first one in header file
155 230 if ( w[page] > wavelength ){
156   - for(unsigned j = 0; j < X(); j++)
  231 + for(unsigned long long j = 0; j < X(); j++)
157 232 {
158 233 p[j] = c[j * B];
159 234 }
... ... @@ -165,7 +240,7 @@ public:
165 240 page++;
166 241 //if wavelength is larger than the last wavelength in header file
167 242 if (page == B) {
168   - for(unsigned j = 0; j < X(); j++)
  243 + for(unsigned long long j = 0; j < X(); j++)
169 244 {
170 245 p[j] = c[(j + 1) * B - 1];
171 246 }
... ... @@ -179,17 +254,17 @@ public:
179 254 p1=(T*)malloc( X() * sizeof(T)); //memory allocation
180 255 p2=(T*)malloc( X() * sizeof(T));
181 256 //band_index(p1, page - 1);
182   - for(unsigned j = 0; j < X(); j++)
  257 + for(unsigned long long j = 0; j < X(); j++)
183 258 {
184 259 p1[j] = c[j * B + page - 1];
185 260 }
186 261 //band_index(p2, page );
187   - for(unsigned j = 0; j < X(); j++)
  262 + for(unsigned long long j = 0; j < X(); j++)
188 263 {
189 264 p2[j] = c[j * B + page];
190 265 }
191 266  
192   - for(unsigned i=0; i < X(); i++){
  267 + for(unsigned long long i=0; i < X(); i++){
193 268 double r = (double) (wavelength - w[page-1]) / (double) (w[page] - w[page-1]);
194 269 p[i] = (p2[i] - p1[i]) * r + p1[i];
195 270 }
... ... @@ -199,7 +274,7 @@ public:
199 274 else //if the wavelength is equal to a wavelength in header file
200 275 {
201 276 //band_index(p, page);
202   - for(unsigned j = 0; j < X(); j++)
  277 + for(unsigned long long j = 0; j < X(); j++)
203 278 {
204 279 p[j] = c[j * B + page];
205 280 }
... ... @@ -208,14 +283,32 @@ public:
208 283 return true;
209 284 }
210 285  
  286 + /// Retrieve a single pixel and store it in a pre-allocated double array.
  287 + bool pixeld(double* p, unsigned long long n){
  288 + unsigned long long bandnum = X() * Y(); //calculate numbers in one band
  289 + if ( n >= bandnum){ //make sure the pixel number is right
  290 + std::cout<<"ERROR: sample or line out of range"<<std::endl;
  291 + return false;
  292 + }
  293 + unsigned long long B = Z();
  294 +
  295 + T* temp = (T*) malloc(B * sizeof(T)); //allocate space for the raw pixel data
  296 + file.seekg(n * B * sizeof(T), std::ios::beg); //point to the certain pixel
  297 + file.read((char *)temp, sizeof(T) * B); //read the spectrum from disk to the temp pointer
  298 +
  299 + for(unsigned long long i = 0; i < B; i++) //for each element of the spectrum
  300 + p[i] = (double) temp[i]; //cast each element to a double value
  301 + return true;
  302 + }
  303 +
211 304 /// Retrieve a single pixel and stores it in pre-allocated memory.
212 305  
213 306 /// @param p is a pointer to pre-allocated memory at least sizeof(T) in size.
214 307 /// @param n is an integer index to the pixel using linear array indexing.
215   - bool pixel(T * p, unsigned n){
  308 + bool pixel(T * p, unsigned long long n){
216 309  
217   - unsigned bandnum = X() * Y(); //calculate numbers in one band
218   - if ( n >= bandnum){ //make sure the pixel number is right
  310 + unsigned long long N = X() * Y(); //calculate numbers in one band
  311 + if ( n >= N){ //make sure the pixel number is right
219 312 std::cout<<"ERROR: sample or line out of range"<<std::endl;
220 313 return false;
221 314 }
... ... @@ -226,7 +319,7 @@ public:
226 319 }
227 320  
228 321 //given a Y ,return a ZX slice
229   - bool read_plane_y(T * p, unsigned y){
  322 + bool read_plane_y(T * p, unsigned long long y){
230 323 return binary<T>::read_plane_2(p, y);
231 324 }
232 325  
... ... @@ -234,113 +327,61 @@ public:
234 327  
235 328 /// @param outname is the name of the output file used to store the resulting baseline-corrected data.
236 329 /// @param wls is the list of baseline points based on band labels.
237   - bool baseline(std::string outname, std::vector<double> wls){
238   -
239   - unsigned N = wls.size(); //get the number of baseline points
  330 + bool baseline(std::string outname, std::vector<double> base_pts, unsigned char* mask = NULL, bool PROGRESS = false){
240 331  
241 332 std::ofstream target(outname.c_str(), std::ios::binary); //open the target binary file
242   - std::string headername = outname + ".hdr"; //the header file name
243 333  
244   - //simplify image resolution
245   - unsigned int ZX = Z() * X(); //calculate the number of points in a Y slice
246   - unsigned int L = ZX * sizeof(T); //calculate the number of bytes of a Y slice
247   - unsigned int B = Z();
248   -
249   - T* c; //pointer to the current Y slice
250   - c = (T*)malloc(L); //memory allocation
251   -
252   - T* a; //pointer to the two YZ lines surrounding the current YZ line
253   - T* b;
254   -
255   - a = (T*)malloc(X() * sizeof(T));
256   - b = (T*)malloc(X() * sizeof(T));
257   -
258   -
259   - double ai, bi; //stores the two baseline points wavelength surrounding the current band
260   - double ci; //stores the current band's wavelength
261   - unsigned control;
262   -
263   - if (a == NULL || b == NULL || c == NULL){
264   - std::cout<<"ERROR: error allocating memory";
265   - exit(1);
266   - }
267   - // loop start correct every y slice
268   -
269   - for (unsigned k =0; k < Y(); k++)
270   - {
271   - //get the current y slice
272   - read_plane_y(c, k);
273   -
274   - //initialize lownum, highnum, low, high
275   - control=0;
276   - ai = w[0];
277   - //if no baseline point is specified at band 0,
278   - //set the baseline point at band 0 to 0
279   - if(wls[0] != w[0]){
280   - bi = wls[control];
281   - memset(a, (char)0, X() * sizeof(T) );
  334 + unsigned long long N = X() * Y(); //calculate the total number of pixels to be processed
  335 + unsigned long long B = Z(); //get the number of bands
  336 + T* s = (T*)malloc(sizeof(T) * B); //allocate memory to store a pixel
  337 + T* sbc = (T*)malloc(sizeof(T) * B); //allocate memory to store the baseline corrected spectrum
  338 +
  339 + std::vector<T> base_vals; //allocate space for the values at each baseline point
  340 + double aw, bw; //surrounding baseline point wavelengths
  341 + T av, bv; //surrounding baseline point values
  342 + unsigned long long ai, bi; //surrounding baseline point band indices
  343 + for(unsigned long long n = 0; n < N; n++){ //for each pixel in the image
  344 + if(mask != NULL && !mask[n]){ //if the pixel isn't masked
  345 + memset(sbc, 0, sizeof(T) * B); //set the baseline corrected spectrum to zero
282 346 }
283   - //else get the low band
284 347 else{
285   - control++;
286   - read_x_from_xz(a, c, ai);
287   - bi = wls[control];
288   - }
289   - //get the high band
290   - read_x_from_xz(b, c, bi);
291   -
292   - //correct every YZ line
293   -
294   - for(unsigned cii = 0; cii < B; cii++){
295   - //update baseline points, if necessary
296   - if( w[cii] >= bi && cii != B - 1) {
297   - //if the high band is now on the last BL point?
298   - if (control != N-1) {
299   -
300   - control++; //increment the index
301   -
302   - std::swap(a, b); //swap the baseline band pointers
303   -
304   - ai = bi;
305   - bi = wls[control];
306   - read_x_from_xz(b, c, bi);
307   -
  348 +
  349 + pixel(s, n); //retrieve the spectrum s
  350 + base_vals = interp_spectrum(s, base_pts); //get the values at each baseline point
  351 +
  352 + ai = bi = 0;
  353 + aw = w[0]; //initialize the current baseline points (assume the spectrum starts at 0)
  354 + av = s[0];
  355 + bw = base_pts[0];
  356 + for(unsigned long long b = 0; b < B; b++){ //for each band in the spectrum
  357 + while(bi < base_pts.size() && base_pts[bi] < w[b]) //while the current wavelength is higher than the second baseline point
  358 + bi++; //move to the next baseline point
  359 + if(bi < base_pts.size()){
  360 + bw = base_pts[bi]; //set the wavelength for the upper bound baseline point
  361 + bv = base_vals[bi]; //set the value for the upper bound baseline point
308 362 }
309   - //if the last BL point on the last band of the file?
310   - else if ( wls[control] < w[B - 1]) {
311   -
312   - std::swap(a, b); //swap the baseline band pointers
313   -
314   - memset(b, (char)0, X() * sizeof(T) ); //clear the high band
315   -
316   - ai = bi;
317   - bi = w[B - 1];
  363 + if(bi == base_pts.size()){ //if we have passed the last baseline point
  364 + bw = w[B-1]; //set the outer bound to the last spectral band
  365 + bv = s[B-1];
318 366 }
  367 + if(bi != 0){
  368 + ai = bi - 1; //set the lower bound baseline point index
  369 + aw = base_pts[ai]; //set the wavelength for the lower bound baseline point
  370 + av = base_vals[ai]; //set the value for the lower bound baseline point
  371 + }
  372 + sbc[b] = s[b] - lerp(w[b], av, aw, bv, bw); //perform the baseline correction and save the new value
319 373 }
320   -
321   - ci = w[cii];
322   -
323   - //perform the baseline correction
324   - for(unsigned i=0; i < X(); i++)
325   - {
326   - double r = (double) (ci - ai) / (double) (bi - ai);
327   - c[i * B + cii] =(T) ( c[i * B + cii] - (b[i] - a[i]) * r - a[i] );
328   - }
329   -
330   - }//loop for YZ line end
331   - target.write(reinterpret_cast<const char*>(c), L); //write the corrected data into destination
332   -
333   - thread_data = (double)k / Y() * 100;
334   - }//loop for Y slice end
  374 + }
335 375  
  376 + if(PROGRESS) progress = (double)(n+1) / N * 100; //set the current progress
336 377  
  378 + target.write((char*)sbc, sizeof(T) * B); //write the corrected data into destination
  379 + } //end for each pixel
337 380  
338   - free(a);
339   - free(b);
340   - free(c);
  381 + free(s);
  382 + free(sbc);
341 383 target.close();
342   -
343   - thread_data = 100;
  384 +
344 385 return true;
345 386  
346 387 }
... ... @@ -350,107 +391,66 @@ public:
350 391 /// @param outname is the name of the output file used to store the resulting baseline-corrected data.
351 392 /// @param w is the label specifying the band that the hyperspectral image will be normalized to.
352 393 /// @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.
353   - bool normalize(std::string outname, double w, double t = 0.0)
  394 + bool normalize(std::string outname, double w, unsigned char* mask = NULL, bool PROGRESS = false)
354 395 {
355   - unsigned int B = Z(); //calculate the number of bands
356   - unsigned int ZX = Z() * X();
357   - unsigned int XY = X() * Y(); //calculate the number of pixels in a band
358   - unsigned int S = XY * sizeof(T); //calculate the number of bytes in a band
359   - unsigned int L = ZX * sizeof(T);
360   -
361 396 std::ofstream target(outname.c_str(), std::ios::binary); //open the target binary file
362   - std::string headername = outname + ".hdr"; //the header file name
363   -
364   - T * c; //pointer to the current ZX slice
365   - T * b; //pointer to the standard band
366   -
367   - b = (T*)malloc( S ); //memory allocation
368   - c = (T*)malloc( L );
369   -
370   - band(b, w); //get the certain band into memory
371   -
372   - for(unsigned j = 0; j < Y(); j++)
373   - {
374   - read_plane_y(c, j);
375   - unsigned jX = j * X(); //to avoid calculating it many times
376   - for(unsigned i = 0; i < X(); i++)
377   - {
378   - unsigned iB = i * B;
379   - for(unsigned m = 0; m < B; m++)
380   - {
381   - if( b[i+jX] < t )
382   - c[m + iB] = (T)0.0;
383   - else
384   - c[m + iB] = c[m + iB] / b[i + jX]; //perform normalization
385   - }
  397 + std::string headername = outname + ".hdr"; //the header file name
  398 +
  399 + unsigned long long N = X() * Y(); //calculate the total number of pixels to be processed
  400 + unsigned long long B = Z(); //get the number of bands
  401 + T* s = (T*)malloc(sizeof(T) * B); //allocate memory to store a pixel
  402 + T nv; //stores the value of the normalized band
  403 + for(unsigned long long n = 0; n < N; n++){ //for each pixel in the image
  404 + if(mask != NULL && !mask[n]) //if the normalization band is below threshold
  405 + memset(s, 0, sizeof(T) * B); //set the output to zero
  406 + else{
  407 + pixel(s, n); //retrieve the spectrum s
  408 + nv = interp_spectrum(s, w); //find the value of the normalization band
  409 +
  410 + for(unsigned long long b = 0; b < B; b++) //for each band in the spectrum
  411 + s[b] /= nv; //divide by the normalization value
386 412 }
387   - target.write(reinterpret_cast<const char*>(c), L); //write normalized data into destination
388   -
389   - thread_data = (double) j / Y() * 100;
390   - }
  413 +
  414 + if(PROGRESS) progress = (double)(n+1) / N * 100; //set the current progress
391 415  
  416 + target.write((char*)s, sizeof(T) * B); //write the corrected data into destination
  417 + } //end for each pixel
392 418  
393   - free(b);
394   - free(c);
  419 + free(s);
395 420 target.close();
396   - thread_data = 100;
397 421 return true;
398 422 }
399 423  
400   - /// Convert the current BIP file to a BSQ file with the specified file name.
401   -
402   - /// @param outname is the name of the output BSQ file to be saved to disk.
403   - bool bsq(std::string outname)
404   - {
405   - std::string temp = outname + "_temp";
406   - std::string headtemp = temp + ".hdr";
407   - //first creat a temporary bil file and convert bip file to bil file
408   - bil(temp);
409   -
410   - stim::bil<T> n;
411   - if(n.open(temp, X(), Y(), Z(), offset, w)==false){ //open infile
412   - std::cout<<"ERROR: unable to open input file"<<std::endl;
413   - exit(1);
414   - }
415   - //then convert bil file to bsq file
416   - n.bsq(outname);
417   - n.close();
418   - remove(temp.c_str());
419   - remove(headtemp.c_str());
420   - return true;
421   - }
422 424  
423 425 /// Convert the current BIP file to a BIL file with the specified file name.
424 426  
425 427 /// @param outname is the name of the output BIL file to be saved to disk.
426   - bool bil(std::string outname)
  428 + bool bil(std::string outname, bool PROGRESS = false)
427 429 {
428   - unsigned int S = X() * Z() * sizeof(T); //calculate the number of bytes in a ZX slice
  430 + unsigned long long S = X() * Z() * sizeof(T); //calculate the number of bytes in a ZX slice
429 431  
430 432 std::ofstream target(outname.c_str(), std::ios::binary);
431   - std::string headername = outname + ".hdr";
  433 + //std::string headername = outname + ".hdr";
432 434  
433 435 T * p; //pointer to the current ZX slice for bip file
434 436 p = (T*)malloc(S);
435 437 T * q; //pointer to the current XZ slice for bil file
436 438 q = (T*)malloc(S);
437 439  
438   - for ( unsigned i = 0; i < Y(); i++)
  440 + for ( unsigned long long i = 0; i < Y(); i++)
439 441 {
440 442 read_plane_y(p, i);
441   - for ( unsigned k = 0; k < Z(); k++)
  443 + for ( unsigned long long k = 0; k < Z(); k++)
442 444 {
443   - unsigned ks = k * X();
444   - for ( unsigned j = 0; j < X(); j++)
  445 + unsigned long long ks = k * X();
  446 + for ( unsigned long long j = 0; j < X(); j++)
445 447 q[ks + j] = p[k + j * Z()];
446 448  
447   - thread_data = (double)(i * Z() + k) / (Y() * Z()) * 100;
  449 + if(PROGRESS) progress = (double)(i * Z() + k+1) / (Y() * Z()) * 100;
448 450 }
449 451 target.write(reinterpret_cast<const char*>(q), S); //write a band data into target file
450 452 }
451 453  
452   - thread_data = 100;
453   -
454 454 free(p);
455 455 free(q);
456 456 target.close();
... ... @@ -467,12 +467,12 @@ public:
467 467 /// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size.
468 468 bool baseline_band(double lb, double rb, T* lp, T* rp, double wavelength, T* result){
469 469  
470   - unsigned XY = X() * Y();
  470 + unsigned long long XY = X() * Y();
471 471 band(result, wavelength); //get band
472 472  
473 473 //perform the baseline correction
474 474 double r = (double) (wavelength - lb) / (double) (rb - lb);
475   - for(unsigned i=0; i < XY; i++){
  475 + for(unsigned long long i=0; i < XY; i++){
476 476 result[i] =(T) (result[i] - (rp[i] - lp[i]) * r - lp[i] );
477 477 }
478 478 return true;
... ... @@ -488,8 +488,8 @@ public:
488 488  
489 489 T* lp;
490 490 T* rp;
491   - unsigned XY = X() * Y();
492   - unsigned S = XY * sizeof(T);
  491 + unsigned long long XY = X() * Y();
  492 + unsigned long long S = XY * sizeof(T);
493 493 lp = (T*) malloc(S); //memory allocation
494 494 rp = (T*) malloc(S);
495 495  
... ... @@ -518,8 +518,8 @@ public:
518 518 T* cur; //current band 1
519 519 T* cur2; //current band 2
520 520  
521   - unsigned XY = X() * Y();
522   - unsigned S = XY * sizeof(T);
  521 + unsigned long long XY = X() * Y();
  522 + unsigned long long S = XY * sizeof(T);
523 523  
524 524 lp = (T*) malloc(S); //memory allocation
525 525 rp = (T*) malloc(S);
... ... @@ -529,9 +529,9 @@ public:
529 529 memset(result, (char)0, S);
530 530  
531 531 //find the wavelenght position in the whole band
532   - unsigned int n = w.size();
533   - unsigned int ai = 0; //left bound position
534