codebase / stimlib

Commit 9d3ba0b10c295c9bd0a6a5b1be4779e482be4515

Authored by David Mayerich
1 parent 798cea97

added stim::hsi as a bridge, fixed warnings for a clean build of HSIproc

Showing 12 changed files with 2096 additions and 1319 deletions   Show diff stats
Inline Side-by-side
stim/envi/bil.h
  Show/Hide comments   View file @9d3ba0b
... ... @@ -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,28 +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   - }
37   -
  36 + }*/
  37 + using hsi<T>::w; //use the wavelength array in stim::hsi
  38 + using hsi<T>::nnz;
38 39 using binary<T>::progress;
39 40  
40 41 /// Call the binary nnz() function for the BIL orientation
41   - unsigned long long nnz(unsigned char* mask){
42   - return binary<T>::nnz(mask, X()*Y());
43   - }
  42 + //unsigned long long nnz(unsigned char* mask){
  43 + // return binary<T>::nnz(mask, X()*Y());
  44 + //}
44 45  
45 46 public:
46 47  
... ... @@ -48,6 +49,8 @@ public:
48 49 using binary<T>::file;
49 50 using binary<T>::R;
50 51  
  52 + bil(){ hsi<T>::init_bil(); }
  53 +
51 54 /// Open a data file for reading using the class interface.
52 55  
53 56 /// @param filename is the name of the binary file on disk
... ... @@ -56,11 +59,16 @@ public:
56 59 /// @param B is the number of samples (bands) along dimension 3
57 60 /// @param header_offset is the number of bytes (if any) in the binary header
58 61 /// @param wavelengths is an optional STL vector of size B specifying a numerical label for each band
59   - 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){
60 68  
61 69 w = wavelengths;
62 70  
63   - return open(filename, vec<unsigned int>(X, B, Y), header_offset);
  71 + return open(filename, vec<unsigned long long>(X, B, Y), header_offset);
64 72  
65 73 }
66 74  
... ... @@ -68,12 +76,12 @@ public:
68 76  
69 77 /// @param p is a pointer to an allocated region of memory at least X * Y * sizeof(T) in size.
70 78 /// @param page <= B is the integer number of the band to be copied.
71   - bool band_index( T * p, unsigned int page, bool PROGRESS = false){
  79 + bool band_index( T * p, unsigned long long page, bool PROGRESS = false){
72 80 //return binary<T>::read_plane_1(p, page);
73 81  
74 82 if(PROGRESS) progress = 0;
75   - unsigned int L = X() * sizeof(T); //caculate the number of bytes in a sample line
76   - unsigned int jump = X() * (Z() - 1) * sizeof(T);
  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);
77 85  
78 86 if (page >= Z()){ //make sure the bank number is right
79 87 std::cout<<"ERROR: page out of range"<<std::endl;
... ... @@ -81,7 +89,7 @@ public:
81 89 }
82 90  
83 91 file.seekg(X() * page * sizeof(T), std::ios::beg);
84   - for (unsigned i = 0; i < Y(); i++)
  92 + for (unsigned long long i = 0; i < Y(); i++)
85 93 {
86 94 file.read((char *)(p + i * X()), L);
87 95 file.seekg( jump, std::ios::cur);
... ... @@ -99,12 +107,12 @@ public:
99 107  
100 108 //if there are no wavelengths in the BSQ file
101 109 if(w.size() == 0)
102   - return band_index(p, (unsigned int)wavelength);
  110 + return band_index(p, (unsigned long long)wavelength);
103 111  
104   - unsigned int XY = X() * Y(); //calculate the number of pixels in a band
105   - 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
106 114  
107   - unsigned page=0; //bands around the wavelength
  115 + unsigned long long page=0; //bands around the wavelength
108 116  
109 117  
110 118 //get the bands numbers around the wavelength
... ... @@ -132,9 +140,9 @@ public:
132 140 p2=(T*)malloc(S);
133 141 band_index(p1, page - 1);
134 142 band_index(p2, page, PROGRESS);
135   - for(unsigned i=0; i < XY; i++){
136   - double r = (double) (wavelength - w[page-1]) / (double) (w[page] - w[page-1]);
137   - p[i] = (p2[i] - p1[i]) * r + p1[i];
  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]);
138 146 }
139 147 free(p1);
140 148 free(p2);
... ... @@ -154,10 +162,10 @@ public:
154 162 /// @param wavelength is the wavelength to retrieve
155 163 bool read_x_from_xz(T* p, T* c, double wavelength)
156 164 {
157   - unsigned int B = Z();
158   - unsigned int L = X() * sizeof(T);
  165 + unsigned long long B = Z();
  166 + unsigned long long L = X() * sizeof(T);
159 167  
160   - unsigned page=0; //samples around the wavelength
  168 + unsigned long long page=0; //samples around the wavelength
161 169 T * p1;
162 170 T * p2;
163 171  
... ... @@ -186,9 +194,9 @@ public:
186 194 memcpy(p1, c + (page - 1) * X(), L);
187 195 memcpy(p2, c + page * X(), L);
188 196  
189   - for(unsigned i=0; i < X(); i++){
  197 + for(unsigned long long i=0; i < X(); i++){
190 198 double r = (double) (wavelength - w[page-1]) / (double) (w[page] - w[page-1]);
191   - p[i] = (p2[i] - p1[i]) * r + p1[i];
  199 + p[i] = (T)(((double)p2[i] - (double)p1[i]) * r + (double)p1[i]);
192 200 }
193 201 }
194 202 else //if the wavelength is equal to a wavelength in header file
... ... @@ -202,7 +210,7 @@ public:
202 210 /// @param p is a pointer to pre-allocated memory at least B * sizeof(T) in size.
203 211 /// @param x is the x-coordinate (dimension 1) of the spectrum.
204 212 /// @param y is the y-coordinate (dimension 2) of the spectrum.
205   - bool spectrum(T * p, unsigned x, unsigned y, bool PROGRESS = false){
  213 + bool spectrum(T * p, unsigned long long x, unsigned long long y, bool PROGRESS = false){
206 214 return binary<T>::read_line_1(p, x, y, PROGRESS);
207 215 }
208 216  
... ... @@ -210,18 +218,18 @@ public:
210 218  
211 219 /// @param p is a pointer to pre-allocated memory at least sizeof(T) in size.
212 220 /// @param n is an integer index to the pixel using linear array indexing.
213   - bool pixel(T * p, unsigned n){
  221 + bool pixel(T * p, unsigned long long n){
214 222  
215 223 //calculate the corresponding x, y
216   - unsigned int x = n % X();
217   - unsigned int y = n / X();
  224 + unsigned long long x = n % X();
  225 + unsigned long long y = n / X();
218 226  
219 227 //get the pixel
220 228 return spectrum(p, x, y);
221 229 }
222 230  
223 231 //given a Y ,return a XZ slice
224   - bool read_plane_y(T * p, unsigned y){
  232 + bool read_plane_y(T * p, unsigned long long y){
225 233 return binary<T>::read_plane_2(p, y);
226 234 }
227 235  
... ... @@ -232,15 +240,15 @@ public:
232 240 /// @param wls is the list of baseline points based on band labels.
233 241 bool baseline(std::string outname, std::vector<double> wls, unsigned char* mask = NULL, bool PROGRESS = false){
234 242  
235   - unsigned N = wls.size(); //get the number of baseline points
  243 + unsigned long long N = wls.size(); //get the number of baseline points
236 244  
237 245 std::ofstream target(outname.c_str(), std::ios::binary); //open the target binary file
238 246 std::string headername = outname + ".hdr"; //the header file name
239 247  
240 248 //simplify image resolution
241   - unsigned int ZX = Z() * X(); //calculate the number of points in a Y slice
242   - unsigned int L = ZX * sizeof(T); //calculate the number of bytes of a Y slice
243   - 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();
244 252  
245 253 T* c; //pointer to the current Y slice
246 254 c = (T*)malloc(L); //memory allocation
... ... @@ -254,7 +262,7 @@ public:
254 262  
255 263 double ai, bi; //stores the two baseline points wavelength surrounding the current band
256 264 double ci; //stores the current band's wavelength
257   - unsigned control;
  265 + unsigned long long control;
258 266  
259 267 if (a == NULL || b == NULL || c == NULL){
260 268 std::cout<<"ERROR: error allocating memory";
... ... @@ -262,7 +270,7 @@ public:
262 270 }
263 271 // loop start correct every y slice
264 272  
265   - for (unsigned k =0; k < Y(); k++)
  273 + for (unsigned long long k =0; k < Y(); k++)
266 274 {
267 275 //get the current y slice
268 276 read_plane_y(c, k);
... ... @@ -287,7 +295,7 @@ public:
287 295  
288 296 //correct every YZ line
289 297  
290   - for(unsigned cii = 0; cii < B; cii++){
  298 + for(unsigned long long cii = 0; cii < B; cii++){
291 299  
292 300 //update baseline points, if necessary
293 301 if( w[cii] >= bi && cii != B - 1) {
... ... @@ -317,7 +325,7 @@ public:
317 325  
318 326 ci = w[cii];
319 327  
320   - unsigned jump = cii * X();
  328 + unsigned long long jump = cii * X();
321 329 //perform the baseline correction
322 330 for(unsigned i=0; i < X(); i++)
323 331 {
... ... @@ -347,11 +355,11 @@ public:
347 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.
348 356 bool normalize(std::string outname, double w, unsigned char* mask = NULL, bool PROGRESS = false)
349 357 {
350   - unsigned int B = Z(); //calculate the number of bands
351   - unsigned int ZX = Z() * X();
352   - unsigned int XY = X() * Y(); //calculate the number of pixels in a band
353   - unsigned int S = XY * sizeof(T); //calculate the number of bytes in a band
354   - 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);
355 363  
356 364 std::ofstream target(outname.c_str(), std::ios::binary); //open the target binary file
357 365 std::string headername = outname + ".hdr"; //the header file name
... ... @@ -364,12 +372,12 @@ public:
364 372  
365 373 band(b, w); //get the certain band into memory
366 374  
367   - for(unsigned j = 0; j < Y(); j++)
  375 + for(unsigned long long j = 0; j < Y(); j++)
368 376 {
369 377 read_plane_y(c, j);
370   - for(unsigned i = 0; i < B; i++)
  378 + for(unsigned long long i = 0; i < B; i++)
371 379 {
372   - for(unsigned m = 0; m < X(); m++)
  380 + for(unsigned long long m = 0; m < X(); m++)
373 381 {
374 382 if( mask != NULL && !mask[m + j * X()] )
375 383 c[m + i * X()] = (T)0.0;
... ... @@ -393,7 +401,7 @@ public:
393 401 /// @param outname is the name of the output BSQ file to be saved to disk.
394 402 bool bsq(std::string outname, bool PROGRESS = false)
395 403 {
396   - 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
397 405  
398 406 std::ofstream target(outname.c_str(), std::ios::binary);
399 407 std::string headername = outname + ".hdr";
... ... @@ -401,7 +409,7 @@ public:
401 409 T * p; //pointer to the current band
402 410 p = (T*)malloc(S);
403 411  
404   - for ( unsigned i = 0; i < Z(); i++)
  412 + for ( unsigned long long i = 0; i < Z(); i++)
405 413 {
406 414 band_index(p, i);
407 415 target.write(reinterpret_cast<const char*>(p), S); //write a band data into target file
... ... @@ -429,13 +437,13 @@ public:
429 437 T * q; //pointer to the current ZX slice for bip file
430 438 q = (T*)malloc(S);
431 439  
432   - for ( unsigned i = 0; i < Y(); i++)
  440 + for ( unsigned long long i = 0; i < Y(); i++)
433 441 {
434 442 read_plane_y(p, i);
435   - for ( unsigned k = 0; k < Z(); k++)
  443 + for ( unsigned long long k = 0; k < Z(); k++)
436 444 {
437   - unsigned ks = k * X();
438   - for ( unsigned j = 0; j < X(); j++)
  445 + unsigned long long ks = k * X();
  446 + for ( unsigned long long j = 0; j < X(); j++)
439 447 q[k + j * Z()] = p[ks + j];
440 448  
441 449 if(PROGRESS) progress = (double)((i+1) * Z() + k+1) / (Z() * Y()) * 100; //store the progress for the current operation
... ... @@ -461,12 +469,12 @@ public:
461 469 /// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size.
462 470 bool baseline_band(double lb, double rb, T* lp, T* rp, double wavelength, T* result){
463 471  
464   - unsigned XY = X() * Y();
  472 + unsigned long long XY = X() * Y();
465 473 band(result, wavelength); //get band
466 474  
467 475 //perform the baseline correction
468 476 double r = (double) (wavelength - lb) / (double) (rb - lb);
469   - for(unsigned i=0; i < XY; i++){
  477 + for(unsigned long long i=0; i < XY; i++){
470 478 result[i] =(T) (result[i] - (rp[i] - lp[i]) * r - lp[i] );
471 479 }
472 480 return true;
... ... @@ -482,8 +490,8 @@ public:
482 490  
483 491 T* lp;
484 492 T* rp;
485   - unsigned XY = X() * Y();
486   - unsigned S = XY * sizeof(T);
  493 + unsigned long long XY = X() * Y();
  494 + unsigned long long S = XY * sizeof(T);
487 495 lp = (T*) malloc(S); //memory allocation
488 496 rp = (T*) malloc(S);
489 497  
... ... @@ -513,8 +521,8 @@ public:
513 521 T* cur; //current band 1
514 522 T* cur2; //current band 2
515 523  
516   - unsigned XY = X() * Y();
517   - unsigned S = XY * sizeof(T);
  524 + unsigned long long XY = X() * Y();
  525 + unsigned long long S = XY * sizeof(T);
518 526  
519 527 lp = (T*) malloc(S); //memory allocation
520 528 rp = (T*) malloc(S);
... ... @@ -524,9 +532,9 @@ public:
524 532 memset(result, (char)0, S);
525 533  
526 534 //find the wavelenght position in the whole band
527   - unsigned int n = w.size();
528   - unsigned int ai = 0; //left bound position
529   - 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
530 538  
531 539  
532 540  
... ... @@ -555,23 +563,23 @@ public:
555 563 //calculate the beginning and the ending part
556 564 baseline_band(lb, rb, lp, rp, rab, cur2); //ending part
557 565 baseline_band(lb, rb, lp, rp, w[bi], cur);
558   - for(unsigned j = 0; j < XY; j++){
559   - 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);
560 568 }
561 569 baseline_band(lb, rb, lp, rp, lab, cur2); //beginnning part
562 570 baseline_band(lb, rb, lp, rp, w[ai], cur);
563   - for(unsigned j = 0; j < XY; j++){
564   - 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);
565 573 }
566 574  
567 575 //calculate the area
568 576 ai++;
569   - for(unsigned i = ai; i <= bi ;i++)
  577 + for(unsigned long long i = ai; i <= bi ;i++)
570 578 {
571 579 baseline_band(lb, rb, lp, rp, w[ai], cur2);
572   - for(unsigned j = 0; j < XY; j++)
  580 + for(unsigned long long j = 0; j < XY; j++)
573 581 {
574   - 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);
575 583 }
576 584 std::swap(cur,cur2); //swap the band pointers
577 585 }
... ... @@ -601,7 +609,7 @@ public:
601 609 height(lb1, rb1, pos1, p1);
602 610 height(lb2, rb2, pos2, p2);
603 611 //calculate the ratio in result
604   - for(unsigned i = 0; i < X() * Y(); i++){
  612 + for(unsigned long long i = 0; i < X() * Y(); i++){
605 613 if(p1[i] == 0 && p2[i] ==0)
606 614 result[i] = 1;
607 615 else
... ... @@ -632,7 +640,7 @@ public:
632 640 area(lb1, rb1, lab1, rab1, p1);
633 641 height(lb2, rb2, pos, p2);
634 642 //calculate the ratio in result
635   - for(unsigned i = 0; i < X() * Y(); i++){
  643 + for(unsigned long long i = 0; i < X() * Y(); i++){
636 644 if(p1[i] == 0 && p2[i] ==0)
637 645 result[i] = 1;
638 646 else
... ... @@ -665,7 +673,7 @@ public:
665 673 area(lb1, rb1, lab1, rab1, p1);
666 674 area(lb2, rb2, lab2, rab2, p2);
667 675 //calculate the ratio in result
668   - for(unsigned i = 0; i < X() * Y(); i++){
  676 + for(unsigned long long i = 0; i < X() * Y(); i++){
669 677 if(p1[i] == 0 && p2[i] ==0)
670 678 result[i] = 1;
671 679 else
... ... @@ -690,8 +698,8 @@ public:
690 698 T* cur; //current band 1
691 699 T* cur2; //current band 2
692 700  
693   - unsigned XY = X() * Y();
694   - unsigned S = XY * sizeof(T);
  701 + unsigned long long XY = X() * Y();
  702 + unsigned long long S = XY * sizeof(T);
695 703  
696 704 lp = (T*) malloc(S); //memory allocation
697 705 rp = (T*) malloc(S);
... ... @@ -701,9 +709,9 @@ public:
701 709 memset(result, (char)0, S);
702 710  
703 711 //find the wavelenght position in the whole band
704   - unsigned int n = w.size();
705   - unsigned int ai = 0; //left bound position
706   - 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
707 715  
708 716 //to make sure the left and the right bound are in the bandwidth
709 717 if (lb < w[0] || rb < w[0] || lb > w[n-1] || rb >w[n-1]){
... ... @@ -730,23 +738,23 @@ public:
730 738 //calculate the beginning and the ending part
731 739 baseline_band(lb, rb, lp, rp, rab, cur2); //ending part
732 740 baseline_band(lb, rb, lp, rp, w[bi], cur);
733   - for(unsigned j = 0; j < XY; j++){
734   - 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);
735 743 }
736 744 baseline_band(lb, rb, lp, rp, lab, cur2); //beginnning part
737 745 baseline_band(lb, rb, lp, rp, w[ai], cur);
738   - for(unsigned j = 0; j < XY; j++){
739   - 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);
740 748 }
741 749  
742 750 //calculate f(x) times x
743 751 ai++;
744   - for(unsigned i = ai; i <= bi ;i++)
  752 + for(unsigned long long i = ai; i <= bi ;i++)
745 753 {
746 754 baseline_band(lb, rb, lp, rp, w[ai], cur2);
747   - for(unsigned j = 0; j < XY; j++)
  755 + for(unsigned long long j = 0; j < XY; j++)
748 756 {
749   - 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);
750 758 }
751 759 std::swap(cur,cur2); //swap the band pointers
752 760 }
... ... @@ -773,7 +781,7 @@ public:
773 781 x_area(lb, rb, lab, rab, p1);
774 782 area(lb, rb, lab, rab, p2);
775 783 //calculate the ratio in result
776   - for(unsigned i = 0; i < X() * Y(); i++){
  784 + for(unsigned long long i = 0; i < X() * Y(); i++){
777 785 if(p1[i] == 0 && p2[i] ==0)
778 786 result[i] = 1;
779 787 else
... ... @@ -797,7 +805,7 @@ public:
797 805 T* temp = (T*)malloc(X() * Y() * sizeof(T)); //allocate memory for the certain band
798 806 band(temp, mask_band, PROGRESS);
799 807  
800   - for (unsigned int i = 0; i < X() * Y(); i++) {
  808 + for (unsigned long long i = 0; i < X() * Y(); i++) {
801 809 if (temp[i] < threshold)
802 810 p[i] = 0;
803 811 else
... ... @@ -817,17 +825,17 @@ public:
817 825 std::ofstream target(outfile.c_str(), std::ios::binary);
818 826  
819 827 //I THINK THIS IS WRONG
820   - unsigned XZ = X() * Z(); //calculate the number of values in a page on disk
821   - 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)
822 830  
823 831 T * temp = (T*)malloc(L); //allocate memory for a temporary page
824 832  
825   - 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)
826 834 {
827 835 read_plane_y(temp, i); //retrieve an ZX slice, stored in temp
828   - 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)
829 837 {
830   - for (unsigned k = 0; k < X(); k++) //for each band
  838 + for (unsigned long long k = 0; k < X(); k++) //for each band
831 839 {
832 840 if(p[i * X() + k] == 0)
833 841 temp[j * X() + k] = 0;
... ... @@ -849,9 +857,9 @@ public:
849 857 std::ofstream target(outfile.c_str(), std::ios::binary);
850 858  
851 859 //for loading pages:
852   - unsigned long XZ = X() * Z(); //calculate the number of values in an XZ page on disk
853   - unsigned long B = Z(); //calculate the number of bands
854   - 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)
855 863  
856 864 //allocate temporary memory for a XZ slice
857 865 T* slice = (T*) malloc(L);
... ... @@ -860,18 +868,18 @@ public:
860 868 T* spec = (T*) malloc(B * sizeof(T));
861 869  
862 870 //for each slice along the y axis
863   - 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()
864 872 {
865 873 read_plane_y(slice, y); //retrieve an ZX page, store in "slice"
866 874  
867 875 //for each sample along X
868   - 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()
869 877 {
870 878 //if the mask != 0 at that xy pixel
871 879 if (p[y * X() + x] != 0) //if the mask != 0 at that XY pixel
872 880 {
873 881 //for each band at that pixel
874   - 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
875 883 {
876 884 spec[b] = slice[b*X() + x]; //Pass the correct spectral value from XZ page into the spectrum to be saved.
877 885 }
... ... @@ -896,20 +904,20 @@ public:
896 904 T* temp = (T*)malloc(sizeof(T) * XZ);
897 905 T* line = (T*)malloc(sizeof(T) * X());
898 906  
899   - for (unsigned i = 0; i < Y(); i++){
  907 + for (unsigned long long i = 0; i < Y(); i++){
900 908 getY(temp, i);
901 909 //initialize x-line
902   - for (unsigned j = 0; j < X(); j++){
  910 + for (unsigned long long j = 0; j < X(); j++){
903 911 line[j] = 0;
904 912 }
905   - unsigned c = 0;
906   - for (unsigned j = 0; j < Z(); j++){
907   - 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++){
908 916 line[k] += temp[c] / (T)Z();
909 917 c++;
910 918 }
911 919 }
912   - for (unsigned j = 0; j < X(); j++){
  920 + for (unsigned long long j = 0; j < X(); j++){
913 921 p[j + i * X()] = line[j];
914 922 }
915 923 }
... ... @@ -925,7 +933,7 @@ public:
925 933 unsigned long long XZ = X() * Z();
926 934 unsigned long long XY = X() * Y();
927 935 T* temp = (T*)malloc(sizeof(T) * XZ);
928   - for (unsigned j = 0; j < Z(); j++){
  936 + for (unsigned long long j = 0; j < Z(); j++){
929 937 p[j] = 0;
930 938 }
931 939 //calculate vaild number in a band
... ... @@ -1038,9 +1046,9 @@ public:
1038 1046 //set the start position for the cropped region
1039 1047 file.seekg((y0 * X() * Z() + b0 * X() + x0) * sizeof(T), std::ios::beg);
1040 1048  
1041   - for (unsigned x = 0; x < lines; x++)
  1049 + for (unsigned long long x = 0; x < lines; x++)
1042 1050 {
1043   - for (unsigned z = b0; z < b1; z++)
  1051 + for (unsigned long long z = b0; z < b1; z++)
1044 1052 {
1045 1053 file.read((char *)(temp + z * samples), sizeof(T) * samples);
1046 1054 file.seekg(jumpb, std::ios::cur); //go to the next band
... ...
stim/envi/binary.h
  Show/Hide comments   View file @9d3ba0b
... ... @@ -25,7 +25,7 @@ protected:
25 25 std::string name; //file name
26 26  
27 27 unsigned long long R[D]; //resolution
28   - unsigned int header; //header size (in bytes)
  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 31 double progress; //stores the progress on the current operation (accessible using a thread)
... ... @@ -33,30 +33,13 @@ protected:
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 40 progress = 0;
41 41 }
42 42  
43   - //calculate the number of non-zero pixels in a mask
44   - unsigned long long nnz(unsigned char* mask, unsigned long long N){
45   -
46   - unsigned long long n = 0; //initialize the counter to 0 (zero)
47   - if(mask == NULL) return N; //if the mask is NULL, assume all pixels are masked
48   -
49   - for(unsigned long long i = 0; i < N; i++){ //for each pixel
50   - if(mask[i] != 0) n++; //increment the counter for every non-zero pixel in the mask
51   - }
52   - return n; //return the number of nonzero pixels
53   - }
54   -
55   - /// Calculate the number of nonzero pixels in a mask over X-Y
56   - unsigned long long nnz(unsigned char* mask){
57   - return nnz(mask, R[0] * R[1]);
58   - }
59   -
60 43 /// Private helper function that returns the size of the file on disk using system functions.
61 44 long long int get_file_size(){
62 45 #ifdef _WIN32
... ... @@ -122,7 +105,7 @@ protected:
122 105  
123 106 public:
124 107  
125   - unsigned int get_progress(){
  108 + double get_progress(){
126 109 return progress;
127 110 }
128 111  
... ... @@ -135,9 +118,9 @@ public:
135 118 /// @param filename is the name of the binary file
136 119 /// @param r is a STIM vector specifying the size of the binary file along each dimension
137 120 /// @param h is the length (in bytes) of any header file (default zero)
138   - 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){
139 122  
140   - 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
141 124 R[i] = r[i];
142 125  
143 126 header = h; //save the header size
... ... @@ -154,17 +137,17 @@ public:
154 137 /// @param filename is the name of the binary file to be created
155 138 /// @param r is a STIM vector specifying the size of the file along each dimension
156 139 /// @offset specifies how many bytes to offset the file (used to leave room for a header)
157   - 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){
158 141  
159 142 std::ofstream target(filename.c_str(), std::ios::binary);
160 143  
161 144 //initialize binary file
162 145 T p = 0;
163   - 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++){
164 147 target.write((char*)(&p), sizeof(T));
165 148 }
166 149  
167   - 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
168 151 R[i] = r[i];
169 152  
170 153 header = offset; //save the header size
... ... @@ -178,7 +161,7 @@ public:
178 161  
179 162 /// @param p is a pointer to the data to be written
180 163 /// @param page is the page number (index of the highest-numbered dimension)
181   - bool write_page( T * p, unsigned int page){
  164 + bool write_page( T * p, unsigned long long page){
182 165  
183 166 if(p == NULL){
184 167 std::cout<<"ERROR: unable to write into file, empty pointer"<<std::endl;
... ... @@ -195,7 +178,7 @@ public:
195 178  
196 179 /// @param p is a pointer to pre-allocated memory equal to the page size
197 180 /// @param page is the index of the page
198   - bool read_page( T * p, unsigned int page, bool PROGRESS = false){
  181 + bool read_page( T * p, unsigned long long page, bool PROGRESS = false){
199 182  
200 183 if(PROGRESS) progress = 0;
201 184  
... ... @@ -217,8 +200,8 @@ public:
217 200 /// @param p is a pointer to pre-allocated memory equal to the line size R[2]
218 201 /// @param x is the x coordinate
219 202 /// @param y is the y coordinate
220   - bool read_line_2( T * p, unsigned int x, unsigned int y, bool PROGRESS = false){
221   - 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;
222 205  
223 206 if(PROGRESS) progress = 0;
224 207  
... ... @@ -244,7 +227,7 @@ public:
244 227 /// @param p is a pointer to pre-allocated memory equal to the line size R[2]
245 228 /// @param x is the y coordinate
246 229 /// @param y is the z coordinate
247   - bool read_line_0(T * p, unsigned int y, unsigned int z, bool PROGRESS = false){
  230 + bool read_line_0(T * p, unsigned long long y, unsigned long long z, bool PROGRESS = false){
248 231 //test to make sure the specified value is within range
249 232 if( y >= R[1] || z >= R[2] ){
250 233 std::cout<<"ERROR: sample or line out of range"<<std::endl;
... ... @@ -262,7 +245,7 @@ public:
262 245 /// @param p is a pointer to pre-allocated memory equal to the line size R[2]
263 246 /// @param x is the y coordinate
264 247 /// @param z is the z coordinate
265   - bool read_line_1(T * p, unsigned int x, unsigned int z, bool PROGRESS = false){
  248 + bool read_line_1(T * p, unsigned long long x, unsigned long long z, bool PROGRESS = false){
266 249 if(PROGRESS) progress = 0;
267 250 //test to make sure the specified value is within range
268 251 if( x >= R[0] || z >= R[2] ){
... ... @@ -271,7 +254,7 @@ public:
271 254 }
272 255  
273 256 file.seekg((z * R[0] * R[1] + x) * sizeof(T), std::ios::beg); //seek to the start of the line
274   - 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
275 258 file.read((char *)(p + i), sizeof(T)); //read the pixel
276 259 file.seekg((R[0] - 1) * sizeof(T), std::ios::cur); //seek to the next pixel in the line
277 260 if(PROGRESS) progress = (double)i / (double)R[1] * 100;
... ... @@ -284,22 +267,22 @@ public:
284 267  
285 268 /// @param p is a pointer to pre-allocated memory of size R[1] * R[2] * sizeof(T)
286 269 /// @param n is the 0-axis coordinate used to retrieve the plane
287   - bool read_plane_0(T* p, unsigned int n, bool PROGRESS = false){
  270 + bool read_plane_0(T* p, unsigned long long n, bool PROGRESS = false){
288 271 if(PROGRESS) progress = 0;
289 272 if (n >= R[0]){ //make sure the number is within the possible range
290 273 std::cout<<"ERROR read_plane_0: page out of range"<<std::endl;
291 274 return false;
292 275 }
293   - 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
294 277  
295 278 //seek to the start of the plane
296 279 file.seekg(n * sizeof(T), std::ios::beg);
297 280  
298   - unsigned int N = R[1] * R[2];
299   - 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++){
300 283 file.read((char*)(p+i), sizeof(T));
301 284 file.seekg(jump, std::ios::cur);
302   - if(PROGRESS) progress = (double)i / N * 100;
  285 + if(PROGRESS) progress = (double)(i+1) / N * 100;
303 286 }
304 287  
305 288 if(PROGRESS) progress = 100;
... ... @@ -312,10 +295,10 @@ public:
312 295  
313 296 /// @param p is a pointer to pre-allocated memory of size R[0] * R[2] * sizeof(T)
314 297 /// @param n is the 1-axis coordinate used to retrieve the plane
315   - bool read_plane_1(T* p, unsigned int n, bool PROGRESS = false){
  298 + bool read_plane_1(T* p, unsigned long long n, bool PROGRESS = false){
316 299 if(PROGRESS) progress = 0;
317   - unsigned int L = R[0] * sizeof(T); //caculate the number of bytes in a sample line
318   - unsigned int jump = R[0] * (R[1] - 1) * sizeof(T);
  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);
319 302  
320 303 if (n >= R[1]){ //make sure the bank number is right
321 304 std::cout<<"ERROR read_plane_1: page out of range"<<std::endl;
... ... @@ -323,7 +306,7 @@ public:
323 306 }
324 307  
325 308 file.seekg(R[0] * n * sizeof(T), std::ios::beg);
326   - for (unsigned i = 0; i < R[2]; i++){
  309 + for (unsigned long long i = 0; i < R[2]; i++){
327 310 if(PROGRESS) progress = (double)i / R[2] * 100;
328 311 file.read((char *)(p + i * R[0]), L);
329 312 file.seekg( jump, std::ios::cur);
... ... @@ -338,7 +321,7 @@ public:
338 321  
339 322 /// @param p is a pointer to pre-allocated memory of size R[0] * R[1] * sizeof(T)
340 323 /// @param n is the 2-axis coordinate used to retrieve the plane
341   - bool read_plane_2(T* p, unsigned int n, bool PROGRESS = false){
  324 + bool read_plane_2(T* p, unsigned long long n, bool PROGRESS = false){
342 325 return read_page(p, n, PROGRESS);
343 326 }
344 327  
... ... @@ -346,7 +329,7 @@ public:
346 329  
347 330 /// @param p is a pointer to pre-allocated memory of size sizeof(T)
348 331 /// @param i is the index to the pixel using linear indexing
349   - bool read_pixel(T* p, unsigned int i){
  332 + bool read_pixel(T* p, unsigned long long i){
350 333 if(i >= R[0] * R[1] * R[2]){
351 334 std::cout<<"ERROR read_pixel: n is out of range"<<std::endl;
352 335 return false;
... ... @@ -363,14 +346,14 @@ public:
363 346 /// @param x is the x (0) axis coordinate
364 347 /// @param y is the y (1) axis coordinate
365 348 /// @param z is the z (2) axis coordinate
366   - 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){
367 350  
368 351 if(x < 0 || x >= R[0] || y < 0 || y >= R[1] || z < 0 || z > R[2]){
369 352 std::cout<<"ERROR read_pixel: (x,y,z) is out of range"<<std::endl;
370 353 return false;
371 354 }
372 355  
373   - 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;
374 357 return read_pixel(p, i);
375 358 }
376 359  
... ...
stim/envi/bip.h
  Show/Hide comments   View file @9d3ba0b
... ... @@ -3,7 +3,7 @@
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  
... ... @@ -23,15 +23,15 @@ namespace stim{
23 23 */
24 24 template <typename T>
25 25  
26   -class bip: public binary<T> {
  26 +class bip: public hsi<T> {
27 27  
28 28 protected:
29 29  
30 30  
31   - std::vector<double> w; //band wavelength
32   - unsigned int offset; //header offset
  31 + //std::vector<double> w; //band wavelength
  32 + unsigned long long offset; //header offset
33 33  
34   - unsigned long long X(){
  34 + /*unsigned long long X(){
35 35 return R[1];
36 36 }
37 37 unsigned long long Y(){
... ... @@ -39,17 +39,18 @@ protected:
39 39 }
40 40 unsigned long long Z(){
41 41 return R[0];
42   - }
43   -
  42 + }*/
  43 + using hsi<T>::w; //use the wavelength array in stim::hsi
  44 + using hsi<T>::nnz;
44 45 using binary<T>::progress;
45 46  
46 47 /// Call the binary nnz() function for the BIP orientation
47   - unsigned long long nnz(unsigned char* mask){
48   - return binary<T>::nnz(mask, X()*Y());
49   - }
  48 + //unsigned long long nnz(unsigned char* mask){
  49 + // return binary<T>::nnz(mask, X()*Y());
  50 + //}
50 51  
51 52 /// Linear interpolation of a spectral value given the bounding spectral samples
52   - T lerp(double w, T low_v, double low_w, T high_v, double high_w){
  53 + /*T lerp(double w, T low_v, double low_w, T high_v, double high_w){
53 54 if(low_w == high_w) return low_v; //if the interval is of zero length, just return one of the bounds
54 55 double alpha = (w - low_w) / (high_w - low_w); //calculate the interpolation factor
55 56 return (1.0 - alpha) * low_v + alpha * high_v; //interpolate
... ... @@ -104,7 +105,7 @@ protected:
104 105 v[n] = interp_spectrum(s, wavelengths[n]); //interpolate the measurement
105 106 }
106 107 return v;
107   - }
  108 + }*/
108 109  
109 110 public:
110 111  
... ... @@ -113,6 +114,8 @@ public:
113 114 using binary<T>::R;
114 115 using binary<T>::read_line_0;
115 116  
  117 + bip(){ hsi<T>::init_bip(); }
  118 +
116 119 /// Open a data file for reading using the class interface.
117 120  
118 121 /// @param filename is the name of the binary file on disk
... ... @@ -121,14 +124,19 @@ public:
121 124 /// @param B is the number of samples (bands) along dimension 3
122 125 /// @param header_offset is the number of bytes (if any) in the binary header
123 126 /// @param wavelengths is an optional STL vector of size B specifying a numerical label for each band
124   - 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){
125 133  
126 134 //copy the wavelengths to the BSQ file structure
127 135 w = wavelengths;
128 136 //copy the offset to the structure
129 137 offset = header_offset;
130 138  
131   - return open(filename, vec<unsigned int>(B, X, Y), header_offset);
  139 + return open(filename, vec<unsigned long long>(B, X, Y), header_offset);
132 140  
133 141 }
134 142  
... ... @@ -136,7 +144,7 @@ public:
136 144  
137 145 /// @param p is a pointer to an allocated region of memory at least X * Y * sizeof(T) in size.
138 146 /// @param page <= B is the integer number of the band to be copied.
139   - bool band_index( T * p, unsigned int page, bool PROGRESS = false){
  147 + bool band_index( T * p, unsigned long long page, bool PROGRESS = false){
140 148 return binary<T>::read_plane_0(p, page, PROGRESS);
141 149 }
142 150  
... ... @@ -148,9 +156,9 @@ public:
148 156  
149 157 //if there are no wavelengths in the BSQ file
150 158 if(w.size() == 0)
151   - return band_index(p, (unsigned int)wavelength, PROGRESS);
  159 + return band_index(p, (unsigned long long)wavelength, PROGRESS);
152 160  
153   - 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
154 162  
155 163 unsigned page=0; //bands around the wavelength
156 164  
... ... @@ -180,9 +188,9 @@ public:
180 188 p2=(T*)malloc( XY * sizeof(T));
181 189 band_index(p1, page - 1);
182 190 band_index(p2, page, PROGRESS);
183   - for(unsigned i=0; i < XY; i++){
  191 + for(unsigned long long i=0; i < XY; i++){
184 192 double r = (double) (wavelength - w[page-1]) / (double) (w[page] - w[page-1]);
185   - p[i] = (p2[i] - p1[i]) * r + p1[i];
  193 + p[i] = (T)(((double)p2[i] - (double)p1[i]) * r + (double)p1[i]);
186 194 }
187 195 free(p1);
188 196 free(p2);
... ... @@ -200,7 +208,7 @@ public:
200 208 /// @param p is a pointer to pre-allocated memory at least B * sizeof(T) in size.
201 209 /// @param x is the x-coordinate (dimension 1) of the spectrum.
202 210 /// @param y is the y-coordinate (dimension 2) of the spectrum.
203   - bool spectrum(T * p, unsigned x, unsigned y, bool PROGRESS = false){
  211 + bool spectrum(T * p, unsigned long long x, unsigned long long y, bool PROGRESS = false){
204 212 return read_line_0(p, x, y, PROGRESS); //read a line in the binary YZ plane (dimension order for BIP is ZXY)
205 213 }
206 214  
... ... @@ -211,16 +219,16 @@ public:
211 219 /// @param wavelength is the wavelength of X values to retrieve
212 220 bool read_x_from_xz(T* p, T* c, double wavelength)
213 221 {
214   - unsigned int B = Z();
  222 + unsigned long long B = Z();
215 223  
216   - unsigned page=0; //samples around the wavelength
  224 + unsigned long long page=0; //samples around the wavelength
217 225  
218 226  
219 227 //get the bands numbers around the wavelength
220 228  
221 229 //if wavelength is smaller than the first one in header file
222 230 if ( w[page] > wavelength ){
223   - for(unsigned j = 0; j < X(); j++)
  231 + for(unsigned long long j = 0; j < X(); j++)
224 232 {
225 233 p[j] = c[j * B];
226 234 }
... ... @@ -232,7 +240,7 @@ public:
232 240 page++;
233 241 //if wavelength is larger than the last wavelength in header file
234 242 if (page == B) {
235   - for(unsigned j = 0; j < X(); j++)
  243 + for(unsigned long long j = 0; j < X(); j++)
236 244 {
237 245 p[j] = c[(j + 1) * B - 1];
238 246 }
... ... @@ -246,17 +254,17 @@ public:
246 254 p1=(T*)malloc( X() * sizeof(T)); //memory allocation
247 255 p2=(T*)malloc( X() * sizeof(T));
248 256 //band_index(p1, page - 1);
249   - for(unsigned j = 0; j < X(); j++)
  257 + for(unsigned long long j = 0; j < X(); j++)
250 258 {
251 259 p1[j] = c[j * B + page - 1];
252 260 }
253 261 //band_index(p2, page );
254   - for(unsigned j = 0; j < X(); j++)
  262 + for(unsigned long long j = 0; j < X(); j++)
255 263 {
256 264 p2[j] = c[j * B + page];
257 265 }
258 266  
259   - for(unsigned i=0; i < X(); i++){
  267 + for(unsigned long long i=0; i < X(); i++){
260 268 double r = (double) (wavelength - w[page-1]) / (double) (w[page] - w[page-1]);
261 269 p[i] = (p2[i] - p1[i]) * r + p1[i];
262 270 }
... ... @@ -266,7 +274,7 @@ public:
266 274 else //if the wavelength is equal to a wavelength in header file
267 275 {
268 276 //band_index(p, page);
269   - for(unsigned j = 0; j < X(); j++)
  277 + for(unsigned long long j = 0; j < X(); j++)
270 278 {
271 279 p[j] = c[j * B + page];
272 280 }
... ... @@ -276,29 +284,30 @@ public:
276 284 }
277 285  
278 286 /// Retrieve a single pixel and store it in a pre-allocated double array.
279   - bool pixeld(double* p, unsigned n){
280   - unsigned bandnum = X() * Y(); //calculate numbers in one band
  287 + bool pixeld(double* p, unsigned long long n){
  288 + unsigned long long bandnum = X() * Y(); //calculate numbers in one band
281 289 if ( n >= bandnum){ //make sure the pixel number is right
282 290 std::cout<<"ERROR: sample or line out of range"<<std::endl;
283 291 return false;
284 292 }
285   - unsigned B = Z();
  293 + unsigned long long B = Z();
286 294  
287 295 T* temp = (T*) malloc(B * sizeof(T)); //allocate space for the raw pixel data
288 296 file.seekg(n * B * sizeof(T), std::ios::beg); //point to the certain pixel
289 297 file.read((char *)temp, sizeof(T) * B); //read the spectrum from disk to the temp pointer
290 298  
291   - for(unsigned int i = 0; i < B; i++) //for each element of the spectrum
  299 + for(unsigned long long i = 0; i < B; i++) //for each element of the spectrum
292 300 p[i] = (double) temp[i]; //cast each element to a double value
  301 + return true;
293 302 }
294 303  
295 304 /// Retrieve a single pixel and stores it in pre-allocated memory.
296 305  
297 306 /// @param p is a pointer to pre-allocated memory at least sizeof(T) in size.
298 307 /// @param n is an integer index to the pixel using linear array indexing.
299   - bool pixel(T * p, unsigned n){
  308 + bool pixel(T * p, unsigned long long n){
300 309  
301   - unsigned N = X() * Y(); //calculate numbers in one band
  310 + unsigned long long N = X() * Y(); //calculate numbers in one band
302 311 if ( n >= N){ //make sure the pixel number is right
303 312 std::cout<<"ERROR: sample or line out of range"<<std::endl;
304 313 return false;
... ... @@ -310,7 +319,7 @@ public:
310 319 }
311 320  
312 321 //given a Y ,return a ZX slice
313   - bool read_plane_y(T * p, unsigned y){
  322 + bool read_plane_y(T * p, unsigned long long y){
314 323 return binary<T>::read_plane_2(p, y);
315 324 }
316 325  
... ... @@ -330,7 +339,6 @@ public:
330 339 std::vector<T> base_vals; //allocate space for the values at each baseline point
331 340 double aw, bw; //surrounding baseline point wavelengths
332 341 T av, bv; //surrounding baseline point values
333   - double alpha;
334 342 unsigned long long ai, bi; //surrounding baseline point band indices
335 343 for(unsigned long long n = 0; n < N; n++){ //for each pixel in the image
336 344 if(mask != NULL && !mask[n]){ //if the pixel isn't masked
... ... @@ -393,14 +401,14 @@ public:
393 401 T* s = (T*)malloc(sizeof(T) * B); //allocate memory to store a pixel
394 402 T nv; //stores the value of the normalized band
395 403 for(unsigned long long n = 0; n < N; n++){ //for each pixel in the image
396   - pixel(s, n); //retrieve the spectrum s
397   - nv = interp_spectrum(s, w); //find the value of the normalization band
398 404 if(mask != NULL && !mask[n]) //if the normalization band is below threshold
399 405 memset(s, 0, sizeof(T) * B); //set the output to zero
400 406 else{
401   - for(unsigned long long b = 0; b < B; b++){ //for each band in the spectrum
  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
402 411 s[b] /= nv; //divide by the normalization value
403   - }
404 412 }
405 413  
406 414 if(PROGRESS) progress = (double)(n+1) / N * 100; //set the current progress
... ... @@ -419,23 +427,23 @@ public:
419 427 /// @param outname is the name of the output BIL file to be saved to disk.
420 428 bool bil(std::string outname, bool PROGRESS = false)
421 429 {
422   - 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
423 431  
424 432 std::ofstream target(outname.c_str(), std::ios::binary);
425   - std::string headername = outname + ".hdr";
  433 + //std::string headername = outname + ".hdr";
426 434  
427 435 T * p; //pointer to the current ZX slice for bip file
428 436 p = (T*)malloc(S);
429 437 T * q; //pointer to the current XZ slice for bil file
430 438 q = (T*)malloc(S);
431 439  
432   - for ( unsigned i = 0; i < Y(); i++)
  440 + for ( unsigned long long i = 0; i < Y(); i++)
433 441 {
434 442 read_plane_y(p, i);
435   - for ( unsigned k = 0; k < Z(); k++)
  443 + for ( unsigned long long k = 0; k < Z(); k++)
436 444 {
437   - unsigned ks = k * X();
438   - for ( unsigned j = 0; j < X(); j++)
  445 + unsigned long long ks = k * X();
  446 + for ( unsigned long long j = 0; j < X(); j++)
439 447 q[ks + j] = p[k + j * Z()];
440 448  
441 449 if(PROGRESS) progress = (double)(i * Z() + k+1) / (Y() * Z()) * 100;
... ... @@ -459,12 +467,12 @@ public:
459 467 /// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size.
460 468 bool baseline_band(double lb, double rb, T* lp, T* rp, double wavelength, T* result){
461 469  
462   - unsigned XY = X() * Y();
  470 + unsigned long long XY = X() * Y();
463 471 band(result, wavelength); //get band
464 472  
465 473 //perform the baseline correction
466 474 double r = (double) (wavelength - lb) / (double) (rb - lb);
467   - for(unsigned i=0; i < XY; i++){
  475 + for(unsigned long long i=0; i < XY; i++){
468 476 result[i] =(T) (result[i] - (rp[i] - lp[i]) * r - lp[i] );
469 477 }
470 478 return true;
... ... @@ -480,8 +488,8 @@ public:
480 488  
481 489 T* lp;
482 490 T* rp;
483   - unsigned XY = X() * Y();
484   - unsigned S = XY * sizeof(T);
  491 + unsigned long long XY = X() * Y();
  492 + unsigned long long S = XY * sizeof(T);
485 493 lp = (T*) malloc(S); //memory allocation
486 494 rp = (T*) malloc(S);
487 495  
... ... @@ -510,8 +518,8 @@ public:
510 518 T* cur; //current band 1
511 519 T* cur2; //current band 2
512 520  
513   - unsigned XY = X() * Y();
514   - unsigned S = XY * sizeof(T);
  521 + unsigned long long XY = X() * Y();
  522 + unsigned long long S = XY * sizeof(T);
515 523  
516 524 lp = (T*) malloc(S); //memory allocation
517 525 rp = (T*) malloc(S);
... ... @@ -521,9 +529,9 @@ public:
521 529 memset(result, (char)0, S);
522 530  
523 531 //find the wavelenght position in the whole band
524   - unsigned int n = w.size();
525   - unsigned int ai = 0; //left bound position
526   - unsigned int bi = n - 1; //right bound position
  532 + unsigned long long n = w.size();
  533 + unsigned long long ai = 0; //left bound position
  534 + unsigned long long bi = n - 1; //right bound position
527 535  
528 536  
529 537  
... ... @@ -552,23 +560,23 @@ public:
552 560 //calculate the beginning and the ending part
553 561 baseline_band(lb, rb, lp, rp, rab, cur2); //ending part
554 562 baseline_band(lb, rb, lp, rp, w[bi], cur);
555   - for(unsigned j = 0; j < XY; j++){
556   - result[j] += (rab - w[bi]) * (cur[j] + cur2[j]) / 2.0;
  563 + for(unsigned long long j = 0; j < XY; j++){
  564 + result[j] += (T)((rab - w[bi]) * ((double)cur[j] + (double)cur2[j]) / 2.0);
557 565 }
558 566 baseline_band(lb, rb, lp, rp, lab, cur2); //beginnning part
559 567 baseline_band(lb, rb, lp, rp, w[ai], cur);
560   - for(unsigned j = 0; j < XY; j++){
561   - result[j] += (w[ai] - lab) * (cur[j] + cur2[j]) / 2.0;
  568 + for(unsigned long long j = 0; j < XY; j++){
  569 + result[j] += (T)((w[ai] - lab) * ((double)cur[j] + (double)cur2[j]) / 2.0);
562 570 }
563 571  
564 572 //calculate the area
565 573 ai++;
566   - for(unsigned i = ai; i <= bi ;i++)
  574 + for(unsigned long long i = ai; i <= bi ;i++)
567 575 {
568 576 baseline_band(lb, rb, lp, rp, w[ai], cur2);
569   - for(unsigned j = 0; j < XY; j++)
  577 + for(unsigned long long j = 0; j < XY; j++)
570 578 {
571   - result[j] += (w[ai] - w[ai-1]) * (cur[j] + cur2[j]) / 2.0;
  579 + result[j] += (T)((w[ai] - w[ai-1]) * ((double)cur[j] + (double)cur2[j]) / 2.0);
572 580 }
573 581 std::swap(cur,cur2); //swap the band pointers
574 582 }
... ... @@ -598,7 +606,7 @@ public:
598 606 height(lb1, rb1, pos1, p1);
599 607 height(lb2, rb2, pos2, p2);
600 608 //calculate the ratio in result
601   - for(unsigned i = 0; i < X() * Y(); i++){
  609 + for(unsigned long long i = 0; i < X() * Y(); i++){
602 610 if(p1[i] == 0 && p2[i] ==0)
603 611 result[i] = 1;
604 612 else
... ... @@ -629,7 +637,7 @@ public:
629 637 area(lb1, rb1, lab1, rab1, p1);
630 638 height(lb2, rb2, pos, p2);
631 639 //calculate the ratio in result
632   - for(unsigned i = 0; i < X() * Y(); i++){
  640 + for(unsigned long long i = 0; i < X() * Y(); i++){
633 641 if(p1[i] == 0 && p2[i] ==0)
634 642 result[i] = 1;
635 643 else
... ... @@ -662,7 +670,7 @@ public:
662 670 area(lb1, rb1, lab1, rab1, p1);
663 671 area(lb2, rb2, lab2, rab2, p2);
664 672 //calculate the ratio in result
665   - for(unsigned i = 0; i < X() * Y(); i++){
  673 + for(unsigned long long i = 0; i < X() * Y(); i++){
666 674 if(p1[i] == 0 && p2[i] ==0)
667 675 result[i] = 1;
668 676 else
... ... @@ -687,8 +695,8 @@ public:
687 695 T* cur; //current band 1
688 696 T* cur2; //current band 2
689 697  
690   - unsigned XY = X() * Y();
691   - unsigned S = XY * sizeof(T);
  698 + unsigned long long XY = X() * Y();
  699 + unsigned long long S = XY * sizeof(T);
692 700  
693 701 lp = (T*) malloc(S); //memory allocation
694 702 rp = (T*) malloc(S);
... ... @@ -698,9 +706,9 @@ public:
698 706 memset(result, (char)0, S);
699 707  
700 708 //find the wavelenght position in the whole band
701   - unsigned int n = w.size();
702   - unsigned int ai = 0; //left bound position
703   - unsigned int bi = n - 1; //right bound position
  709 + unsigned long long n = w.size();
  710 + unsigned long long ai = 0; //left bound position
  711 + unsigned long long bi = n - 1; //right bound position
704 712  
705 713 //to make sure the left and the right bound are in the bandwidth
706 714 if (lb < w[0] || rb < w[0] || lb > w[n-1] || rb >w[n-1]){
... ... @@ -727,23 +735,23 @@ public:
727 735 //calculate the beginning and the ending part
728 736 baseline_band(lb, rb, lp, rp, rab, cur2); //ending part
729 737 baseline_band(lb, rb, lp, rp, w[bi], cur);
730   - for(unsigned j = 0; j < XY; j++){
731   - result[j] += (rab - w[bi]) * (rab + w[bi]) * (cur[j] + cur2[j]) / 4.0;
  738 + for(unsigned long long j = 0; j < XY; j++){
  739 + result[j] += (T)((rab - w[bi]) * (rab + w[bi]) * ((double)cur[j] + (double)cur2[j]) / 4.0);
732 740 }
733 741 baseline_band(lb, rb, lp, rp, lab, cur2); //beginnning part
734 742 baseline_band(lb, rb, lp, rp, w[ai], cur);
735   - for(unsigned j = 0; j < XY; j++){
736   - result[j] += (w[ai] - lab) * (w[ai] + lab) * (cur[j] + cur2[j]) / 4.0;
  743 + for(unsigned long long j = 0; j < XY; j++){
  744 + result[j] += (T)((w[ai] - lab) * (w[ai] + lab) * ((double)cur[j] + (double)cur2[j]) / 4.0);
737 745 }
738 746  
739 747 //calculate f(x) times x
740 748 ai++;
741   - for(unsigned i = ai; i <= bi ;i++)
  749 + for(unsigned long long i = ai; i <= bi ;i++)
742 750 {
743 751 baseline_band(lb, rb, lp, rp, w[ai], cur2);
744   - for(unsigned j = 0; j < XY; j++)
  752 + for(unsigned long long j = 0; j < XY; j++)
745 753 {
746   - result[j] += (w[ai] - w[ai-1]) * (w[ai] + w[ai-1]) * (cur[j] + cur2[j]) / 4.0;
  754 + result[j] += (T)((w[ai] - w[ai-1]) * (w[ai] + w[ai-1]) * ((double)cur[j] + (double)cur2[j]) / 4.0);
747 755 }
748 756 std::swap(cur,cur2); //swap the band pointers
749 757 }
... ... @@ -770,7 +778,7 @@ public:
770 778 x_area(lb, rb, lab, rab, p1);
771 779 area(lb, rb, lab, rab, p2);
772 780 //calculate the ratio in result
773   - for(unsigned i = 0; i < X() * Y(); i++){
  781 + for(unsigned long long i = 0; i < X() * Y(); i++){
774 782 if(p1[i] == 0 && p2[i] ==0)
775 783 result[i] = 1;
776 784 else
... ... @@ -794,7 +802,7 @@ public:
794 802 T* temp = (T*)malloc(X() * Y() * sizeof(T)); //allocate memory for the certain band
795 803 band(temp, mask_band, PROGRESS);
796 804  
797   - for (unsigned int i = 0; i < X() * Y();i++) {
  805 + for (unsigned long long i = 0; i < X() * Y();i++) {
798 806 if (temp[i] < threshold)
799 807 p[i] = 0;
800 808 else
... ... @@ -814,17 +822,17 @@ public:
814 822  
815 823 std::ofstream target(outfile.c_str(), std::ios::binary);
816 824  
817   - unsigned ZX = Z() * X(); //calculate the number of values in a page (XZ in BIP)
818   - unsigned L = ZX * sizeof(T); //calculate the number of bytes in a page
  825 + unsigned long long ZX = Z() * X(); //calculate the number of values in a page (XZ in BIP)
  826 + unsigned long long L = ZX * sizeof(T); //calculate the number of bytes in a page
819 827  
820 828 T * temp = (T*)malloc(L); //allocate space for that page
821 829  
822   - for (unsigned i = 0; i < Y(); i++) //for each page (Y in BIP)
  830 + for (unsigned long long i = 0; i < Y(); i++) //for each page (Y in BIP)
823 831 {
824 832 read_plane_y(temp, i); //load that page (it's pointed to by temp)
825   - for ( unsigned j = 0; j < X(); j++) //for each X value
  833 + for ( unsigned long long j = 0; j < X(); j++) //for each X value
826 834 {
827   - for (unsigned k = 0; k < Z(); k++) //for each B value (band)
  835 + for (unsigned long long k = 0; k < Z(); k++) //for each B value (band)
828 836 {
829 837 if (p[i * X() + j] == 0) //if the mask value is zero
830 838 temp[j * Z() + k] = 0; //set the pixel value to zero
... ... @@ -850,15 +858,15 @@ public:
850 858 std::ofstream target(outfile.c_str(), std::ios::binary);
851 859  
852 860 //allocate space for a single spectrum
853   - unsigned long B = Z();
  861 + unsigned long long B = Z();
854 862 T* spectrum = (T*) malloc(B * sizeof(T));
855 863  
856 864 //calculate the number of pixels in a band
857   - unsigned long XY = X() * Y();
  865 + unsigned long long XY = X() * Y();
858 866  
859 867 //for each pixel
860   - unsigned long skip = 0; //number of spectra to skip
861   - for(unsigned long x = 0; x < XY; x++){
  868 + unsigned long long skip = 0; //number of spectra to skip
  869 + for(unsigned long long x = 0; x < XY; x++){
862 870  
863 871 //if the current pixel isn't masked
864 872 if( mask[x] == 0){
... ... @@ -879,9 +887,8 @@ public:
879 887  
880 888 //write this pixel out
881 889 target.write((char *)spectrum, B * sizeof(T));
882   -
883   - if(PROGRESS) progress = (double) (x+1) / XY * 100;
884 890 }
  891 + if(PROGRESS) progress = (double) (x+1) / XY * 100;
885 892  
886 893 }
887 894  
... ... @@ -892,7 +899,7 @@ public:
892 899 return true;
893 900 }
894 901  
895   - bool unsift(std::string outfile, unsigned char* mask, unsigned int samples, unsigned int lines){
  902 + bool unsift(std::string outfile, unsigned char* mask, unsigned long long samples, unsigned long long lines){
896 903  
897 904 // open an output stream
898 905 std::ofstream target(outfile.c_str(), std::ios::binary);
... ... @@ -901,7 +908,7 @@ public:
901 908 file.seekg(0, std::ios::beg);
902 909  
903 910 //allocate space for a single spectrum
904   - unsigned long B = Z();
  911 + unsigned long long B = Z();
905 912 T* spectrum = (T*) malloc(B * sizeof(T));
906 913  
907 914 //allocate space for a spectrum of zeros
... ... @@ -909,11 +916,11 @@ public:
909 916 memset(zeros, 0, B * sizeof(T));
910 917  
911 918 //calculate the number of pixels in a band
912   - unsigned long XY = samples * lines;
  919 + unsigned long long XY = samples * lines;
913 920  
914 921 //for each pixel
915   - unsigned long skip = 0; //number of spectra to skip
916   - for(unsigned long x = 0; x < XY; x++){
  922 + unsigned long long skip = 0; //number of spectra to skip
  923 + for(unsigned long long x = 0; x < XY; x++){
917 924  
918 925 //if the current pixel isn't masked
919 926 if( mask[x] == 0){
... ... @@ -951,24 +958,23 @@ public:
951 958 /// @param p is a pointer to pre-allocated memory of size [B * sizeof(T)] that stores the mean spectrum
952 959 /// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
953 960 bool avg_band(double* p, unsigned char* mask = NULL, bool PROGRESS = false){
954   - unsigned long long XY = X() * Y();
955   - T* temp = (T*)malloc(sizeof(T) * Z());
956   - //Iinitialize
957   - for (unsigned j = 0; j < Z(); j++){
958   - p[j] = 0;
959   - }
960   - //calculate vaild number in a band
961   - unsigned count = nnz(mask);
962   -
963   - //calculate average number of a band
964   - for (unsigned i = 0; i < XY; i++){
965   - if (mask == NULL || mask[i] != 0){
966   - pixel(temp, i);
967   - for (unsigned j = 0; j < Z(); j++){
968   - p[j] += (double)temp[j] / (double)count;
  961 + unsigned long long XY = X() * Y(); //calculate the total number of pixels in the HSI
  962 + T* temp = (T*)malloc(sizeof(T) * Z()); //allocate space for the current spectrum to be read
  963 + memset(p, 0, sizeof(double) * Z()); //initialize the average spectrum to zero (0)
  964 + //for (unsigned j = 0; j < Z(); j++){
  965 + // p[j] = 0;
  966 + //}
  967 +
  968 + unsigned long long count = nnz(mask); //calculate the number of masked pixels
  969 +
  970 + for (unsigned long long i = 0; i < XY; i++){ //for each pixel in the HSI
  971 + if (mask == NULL || mask[i] != 0){ //if the pixel is masked
  972 + pixel(temp, i); //get the spectrum
  973 + for (unsigned long long j = 0; j < Z(); j++){ //for each spectral component
  974 + p[j] += (double)temp[j] / (double)count; //add the weighted value to the average
969 975 }
970 976 }
971   - if(PROGRESS) progress = (double)(i+1) / XY * 100;
  977 + if(PROGRESS) progress = (double)(i+1) / XY * 100; //increment the progress
972 978 }
973 979  
974 980 free(temp);
... ... @@ -976,6 +982,7 @@ public:
976 982 }
977 983 #ifdef CUDA_FOUND
978 984 /// Calculate the covariance matrix for masked pixels using cuBLAS
  985 + /// Note that cuBLAS only supports integer-sized arrays, so there may be issues with large spectra
979 986 bool co_matrix_cublas(double* co, double* avg, unsigned char *mask, bool PROGRESS = false){
980 987  
981 988 cudaError_t cudaStat;
... ... @@ -994,7 +1001,7 @@ public:
994 1001 cudaStat = cudaMalloc(&A_dev, B * B * sizeof(double)); //allocate space on the CUDA device for the covariance matrix
995 1002 cudaStat = cudaMemset(A_dev, 0, B * B * sizeof(double)); //initialize the covariance matrix to zero (0)
996 1003 cudaStat = cudaMalloc(&avg_dev, B * sizeof(double)); //allocate space on the CUDA device for the average spectrum
997   - stat = cublasSetVector(B, sizeof(double), avg, 1, avg_dev, 1); //copy the average spectrum to the CUDA device
  1004 + stat = cublasSetVector((int)B, sizeof(double), avg, 1, avg_dev, 1); //copy the average spectrum to the CUDA device
998 1005  
999 1006 double ger_alpha = 1.0/(double)XY; //scale the outer product by the inverse of the number of samples (mean outer product)
1000 1007 double axpy_alpha = -1; //multiplication factor for the average spectrum (in order to perform a subtraction)
... ... @@ -1007,22 +1014,22 @@ public:
1007 1014 for (unsigned long long xy = 0; xy < XY; xy++){ //for each pixel
1008 1015 if (mask == NULL || mask[xy] != 0){
1009 1016 pixeld(s, xy); //retreive the spectrum at the current xy pixel location
1010   - stat = cublasSetVector(B, sizeof(double), s, 1, s_dev, 1); //copy the spectrum from the host to the device
1011   - stat = cublasDaxpy(handle, B, &axpy_alpha, avg_dev, 1, s_dev, 1); //subtract the average spectrum
1012   - stat = cublasDsyr(handle, CUBLAS_FILL_MODE_UPPER, B, &ger_alpha, s_dev, 1, A_dev, B); //calculate the covariance matrix (symmetric outer product)
  1017 + stat = cublasSetVector((int)B, sizeof(double), s, 1, s_dev, 1); //copy the spectrum from the host to the device
  1018 + stat = cublasDaxpy(handle, (int)B, &axpy_alpha, avg_dev, 1, s_dev, 1); //subtract the average spectrum
  1019 + stat = cublasDsyr(handle, CUBLAS_FILL_MODE_UPPER, (int)B, &ger_alpha, s_dev, 1, A_dev, (int)B); //calculate the covariance matrix (symmetric outer product)
1013 1020 }
1014 1021 if(PROGRESS) progress = (double)(xy+1) / XY * 100; //record the current progress
1015 1022  
1016 1023 }
1017 1024  
1018   - cublasGetMatrix(B, B, sizeof(double), A_dev, B, co, B); //copy the result from the GPU to the CPU
  1025 + cublasGetMatrix((int)B, (int)B, sizeof(double), A_dev, (int)B, co, (int)B); //copy the result from the GPU to the CPU
1019 1026  
1020 1027 cudaFree(A_dev); //clean up allocated device memory
1021 1028 cudaFree(s_dev);
1022 1029 cudaFree(avg_dev);
1023 1030  
1024   - for(unsigned i = 0; i < B; i++){ //copy the upper triangular portion to the lower triangular portion
1025   - for(unsigned j = i+1; j < B; j++){
  1031 + for(unsigned long long i = 0; i < B; i++){ //copy the upper triangular portion to the lower triangular portion
  1032 + for(unsigned long long j = i+1; j < B; j++){
1026 1033 co[B * i + j] = co[B * j + i];
1027 1034 }
1028 1035 }
... ... @@ -1046,14 +1053,13 @@ public:
1046 1053 unsigned long long XY = X() * Y();
1047 1054 unsigned long long B = Z();
1048 1055 T* temp = (T*)malloc(sizeof(T) * B);
1049   - //count vaild pixels in a band
1050   - unsigned long long count = nnz(mask);
  1056 +
  1057 + unsigned long long count = nnz(mask); //count the number of masked pixels
1051 1058  
1052 1059 //initialize covariance matrix
1053 1060 memset(co, 0, B * B * sizeof(double));
1054 1061  
1055 1062 //calculate covariance matrix
1056   - T i_diff; //stores (temp[i] - avg[i]) to speed math
1057 1063 double* co_half = (double*) malloc(B * B * sizeof(double)); //allocate space for a higher-precision intermediate matrix
1058 1064 double* temp_precise = (double*) malloc(B * sizeof(double));
1059 1065 memset(co_half, 0, B * B * sizeof(double)); //initialize the high-precision matrix with zeros
... ... @@ -1083,10 +1089,10 @@ public:
1083 1089 return true;
1084 1090 }
1085 1091  
1086   - bool project(std::string outfile, double* center, double* basis, unsigned long long M, bool PROGRESS = false){
  1092 + bool project(std::string outfile, double* center, double* basis, unsigned long long M, unsigned char* mask = NULL, bool PROGRESS = false){
1087 1093  
1088 1094 std::ofstream target(outfile.c_str(), std::ios::binary); //open the target binary file
1089   - std::string headername = outfile + ".hdr"; //the header file name
  1095 + //std::string headername = outfile + ".hdr"; //the header file name
1090 1096  
1091 1097 //memory allocation
1092 1098 unsigned long long XY = X() * Y();
... ... @@ -1096,15 +1102,17 @@ public:
1096 1102 T* rs = (T*)malloc(sizeof(T) * M); //allocate space for the projected spectrum
1097 1103 double* bv; //pointer to the current projection vector
1098 1104 for(unsigned long long xy = 0; xy < XY; xy++){ //for each spectrum in the image
1099   - pixel(s, xy); //load the spectrum
1100   - memset(rs, 0, sizeof(T) * M); //initialize the rotated spectrum to zero (0)
1101   - for(unsigned long long m = 0; m < M; m++){ //for each basis vector
1102   - bv = &basis[m * B]; //assign 'bv' to the beginning of the basis vector
1103   - for(unsigned long long b = 0; b < B; b++){ //for each band
1104   - rs[m] += (s[b] - center[b]) * bv[b]; //center the spectrum and perform the projection
  1105 + memset(rs, 0, sizeof(T) * M);
  1106 + if(mask == NULL || mask[xy]){
  1107 + pixel(s, xy); //load the spectrum
  1108 + for(unsigned long long m = 0; m < M; m++){ //for each basis vector
  1109 + bv = &basis[m * B]; //assign 'bv' to the beginning of the basis vector
  1110 + for(unsigned long long b = 0; b < B; b++){ //for each band
  1111 + rs[m] += (T)(((double)s[b] - center[b]) * bv[b]); //center the spectrum and perform the projection
  1112 + }
1105 1113 }
1106 1114 }
1107   -
  1115 +
1108 1116 target.write(reinterpret_cast<const char*>(rs), sizeof(T) * M); //write the projected vector
1109 1117 if(PROGRESS) progress = (double)(xy+1) / XY * 100;
1110 1118 }
... ... @@ -1135,7 +1143,7 @@ public:
1135 1143 for(unsigned long long c = 0; c < C; c++){ //for each basis vector coefficient
1136 1144 bv = &basis[c * B]; //assign 'bv' to the beginning of the basis vector
1137 1145 for(unsigned long long b = 0; b < B; b++){ //for each component of the basis vector
1138   - s[b] += coeff[c] * bv[b] + center[b]; //calculate the contribution of each element of the basis vector in the final spectrum
  1146 + s[b] += (T)((double)coeff[c] * bv[b] + center[b]); //calculate the contribution of each element of the basis vector in the final spectrum
1139 1147 }
1140 1148 }
1141 1149  
... ...
stim/envi/bsq.h
  Show/Hide comments   View file @9d3ba0b
... ... @@ -2,7 +2,7 @@
2 2 #define STIM_BSQ_H
3 3  
4 4 #include "../envi/envi_header.h"
5   -#include "../envi/binary.h"
  5 +#include "../envi/hsi.h"
6 6 #include "../envi/bil.h"
7 7 #include <cstring>
8 8 #include <utility>
... ... @@ -22,27 +22,27 @@ namespace stim{
22 22 */
23 23 template <typename T>
24 24  
25   -class bsq: public binary<T> {
  25 +class bsq: public hsi<T> {
26 26  
27 27  
28 28 protected:
29 29  
30   - std::vector<double> w; //band wavelengths
31   - unsigned int offset;
  30 + //std::vector<double> w; //band wavelengths
  31 + unsigned long long offset;
32 32  
33 33 using binary<T>::R;
34 34  
35   - unsigned long X(){
  35 + /*unsigned long long X(){
36 36 return R[0];
37 37 }
38   - unsigned long Y(){
  38 + unsigned long long Y(){
39 39 return R[1];
40 40 }
41   - unsigned long Z(){
  41 + unsigned long long Z(){
42 42 return R[2];
43   - }
44   -
45   - using binary<T>::nnz;
  43 + }*/
  44 + using hsi<T>::w; //use the wavelength array in stim::hsi
  45 + using hsi<T>::nnz;
46 46 using binary<T>::progress;
47 47  
48 48 public:
... ... @@ -51,8 +51,8 @@ public:
51 51 using binary<T>::file;
52 52 using binary<T>::read_line_2;
53 53 using binary<T>::read_plane_2;
54   - //using binary<T>::getSlice;
55 54  
  55 + bsq(){ hsi<T>::init_bsq(); }
56 56  
57 57 /// Open a data file for reading using the class interface.
58 58  
... ... @@ -62,21 +62,26 @@ public:
62 62 /// @param B is the number of samples (bands) along dimension 3
63 63 /// @param header_offset is the number of bytes (if any) in the binary header
64 64 /// @param wavelengths is an STL vector of size B specifying a numerical label for each band
65   - bool open(std::string filename, unsigned int X, unsigned int Y, unsigned int B, unsigned int header_offset, std::vector<double> wavelengths){
  65 + bool open(std::string filename,
  66 + unsigned long long X,
  67 + unsigned long long Y,
  68 + unsigned long long B,
  69 + unsigned long long header_offset,
  70 + std::vector<double> wavelengths){
66 71  
67 72 //copy the wavelengths to the BSQ file structure
68 73 w = wavelengths;
69 74 //copy the wavelengths to the structure
70 75 offset = header_offset;
71 76  
72   - return open(filename, vec<unsigned int>(X, Y, B), header_offset);
  77 + return open(filename, vec<unsigned long long>(X, Y, B), header_offset);
73 78 }
74 79  
75 80 /// Retrieve a single band (based on index) and stores it in pre-allocated memory.
76 81  
77 82 /// @param p is a pointer to an allocated region of memory at least X * Y * sizeof(T) in size.
78 83 /// @param page <= B is the integer number of the band to be copied.
79   - bool band_index( T * p, unsigned int page){
  84 + bool band_index( T * p, unsigned long long page){
80 85 return read_plane_2(p, page); //call the binary read_plane function (don't let it update the progress)
81 86 }
82 87  
... ... @@ -88,10 +93,10 @@ public:
88 93 if(PROGRESS) progress = 0;
89 94 //if there are no wavelengths in the BSQ file
90 95 if(w.size() == 0)
91   - return band_index(p, (unsigned int)wavelength);
  96 + return band_index(p, (unsigned long long)wavelength);
92 97  
93   - unsigned int XY = X() * Y(); //calculate the number of pixels in a band
94   - unsigned page = 0;
  98 + unsigned long long XY = X() * Y(); //calculate the number of pixels in a band
  99 + unsigned long long page = 0;
95 100  
96 101  
97 102 //get the two neighboring bands (above and below 'wavelength')
... ... @@ -122,9 +127,9 @@ public:
122 127 p2=(T*)malloc( XY * sizeof(T));
123 128 band_index(p1, page - 1);
124 129 band_index(p2, page );
125   - for(unsigned i=0; i < XY; i++){
126   - double r = (double) (wavelength - w[page-1]) / (double) (w[page] - w[page-1]);
127   - p[i] = (p2[i] - p1[i]) * r + p1[i];
  130 + for(unsigned long long i=0; i < XY; i++){
  131 + double r = (wavelength - w[page-1]) / (w[page] - w[page-1]);
  132 + p[i] = (T)(((double)p2[i] - (double)p1[i]) * r + (double)p1[i]);
128 133 }
129 134 free(p1);
130 135 free(p2);
... ... @@ -142,7 +147,7 @@ public:
142 147 /// @param p is a pointer to pre-allocated memory at least B * sizeof(T) in size.
143 148 /// @param x is the x-coordinate (dimension 1) of the spectrum.
144 149 /// @param y is the y-coordinate (dimension 2) of the spectrum.
145   - bool spectrum(T * p, unsigned x, unsigned y, bool PROGRESS = false){
  150 + bool spectrum(T * p, unsigned long long x, unsigned long long y, bool PROGRESS = false){
146 151 return read_line_2(p, x, y, PROGRESS);
147 152 }
148 153  
... ... @@ -150,16 +155,16 @@ public:
150 155  
151 156 /// @param p is a pointer to pre-allocated memory at least sizeof(T) in size.
152 157 /// @param n is an integer index to the pixel using linear array indexing.
153   - bool pixel(T * p, unsigned n){
  158 + bool pixel(T * p, unsigned long long n){
154 159  
155   - unsigned bandnum = X() * Y(); //calculate numbers in one band
  160 + unsigned long long bandnum = X() * Y(); //calculate numbers in one band
156 161 if ( n >= bandnum){ //make sure the pixel number is right
157 162 std::cout<<"ERROR: sample or line out of range"<<std::endl;
158 163 return false;
159 164 }
160 165  
161 166 file.seekg(n * sizeof(T) + binary<T>::header, std::ios::beg); //point to the certain pixel
162   - for (unsigned i = 0; i < Z(); i++)
  167 + for (unsigned long long i = 0; i < Z(); i++)
163 168 {
164 169 file.read((char *)(p + i), sizeof(T));
165 170 file.seekg((bandnum - 1) * sizeof(T), std::ios::cur); //go to the next band
... ... @@ -174,20 +179,20 @@ public:
174 179 /// @param wls is the list of baseline points based on band labels.
175 180 bool baseline(std::string outname, std::vector<double> wls, unsigned char* mask = NULL, bool PROGRESS = false )
176 181 {
177   - unsigned N = wls.size(); //get the number of baseline points
  182 + size_t N = wls.size(); //get the number of baseline points
178 183  
179 184 std::ofstream target(outname.c_str(), std::ios::binary); //open the target binary file
180 185 std::string headername = outname + ".hdr"; //the header file name
181 186  
182 187 //simplify image resolution
183   - unsigned int B = Z(); //calculate the number of bands
184   - unsigned int XY = X() * Y(); //calculate the number of pixels in a band
185   - unsigned int S = XY * sizeof(T); //calculate the number of bytes in a band
  188 + unsigned long long B = Z(); //calculate the number of bands
  189 + unsigned long long XY = X() * Y(); //calculate the number of pixels in a band
  190 + unsigned long long S = XY * sizeof(T); //calculate the number of bytes in a band
186 191  
187 192 double ai, bi; //stores the two baseline points wavelength surrounding the current band
188 193 double ci; //stores the current band's wavelength
189 194  
190   - unsigned control=0;
  195 + unsigned long long control=0;
191 196  
192 197 T * a; //pointers to the high and low band images
193 198 T * b;
... ... @@ -222,7 +227,7 @@ public:
222 227 band(b, bi);
223 228  
224 229 //correct every band
225   - for(unsigned cii = 0; cii < B; cii++){
  230 + for(unsigned long long cii = 0; cii < B; cii++){
226 231  
227 232 //update baseline points, if necessary
228 233 if( w[cii] >= bi && cii != B - 1) {
... ... @@ -255,7 +260,7 @@ public:
255 260 ci = w[cii];
256 261  
257 262 //perform the baseline correction
258   - for(unsigned i=0; i < XY; i++){
  263 + for(unsigned long long i=0; i < XY; i++){
259 264 if(mask != NULL && !mask[i]) //if the pixel is excluded by a mask
260 265 c[i] = 0; //set the value to zero
261 266 else{
... ... @@ -287,9 +292,9 @@ public:
287 292 /// @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.
288 293 bool normalize(std::string outname, double w, unsigned char* mask = NULL, bool PROGRESS = false)
289 294 {
290   - unsigned int B = Z(); //calculate the number of bands
291   - unsigned int XY = X() * Y(); //calculate the number of pixels in a band
292   - unsigned int S = XY * sizeof(T); //calculate the number of bytes in a band
  295 + unsigned long long B = Z(); //calculate the number of bands
  296 + unsigned long long XY = X() * Y(); //calculate the number of pixels in a band
  297 + unsigned long long S = XY * sizeof(T); //calculate the number of bytes in a band
293 298  
294 299 std::ofstream target(outname.c_str(), std::ios::binary); //open the target binary file
295 300 std::string headername = outname + ".hdr"; //the header file name
... ... @@ -302,10 +307,10 @@ public:
302 307  
303 308 band(b, w); //get the certain band into memory
304 309  
305   - for(unsigned j = 0; j < B; j++)
  310 + for(unsigned long long j = 0; j < B; j++)
306 311 {
307 312 band_index(c, j); //get the current band into memory
308   - for(unsigned i = 0; i < XY; i++)
  313 + for(unsigned long long i = 0; i < XY; i++)
309