Commit a2bf1d0853c4497c8a94fec854b7e99892ea5b4f

Authored by David Mayerich
1 parent c8c976a9

general bug fixes and clean compiling and exiting for HSIproc

stim/envi/binary.h
@@ -285,7 +285,6 @@ public: @@ -285,7 +285,6 @@ public:
285 if(PROGRESS) progress = (double)(i+1) / N * 100; 285 if(PROGRESS) progress = (double)(i+1) / N * 100;
286 } 286 }
287 287
288 - if(PROGRESS) progress = 100;  
289 return true; 288 return true;
290 289
291 290
@@ -8,8 +8,10 @@ @@ -8,8 +8,10 @@
8 #include <utility> 8 #include <utility>
9 9
10 //CUDA 10 //CUDA
11 -#include <cuda_runtime.h>  
12 -#include "cublas_v2.h" 11 +#ifdef CUDA_FOUND
  12 + #include <cuda_runtime.h>
  13 + #include "cublas_v2.h"
  14 +#endif
13 15
14 namespace stim{ 16 namespace stim{
15 17
@@ -31,82 +33,10 @@ protected: @@ -31,82 +33,10 @@ protected:
31 //std::vector<double> w; //band wavelength 33 //std::vector<double> w; //band wavelength
32 unsigned long long offset; //header offset 34 unsigned long long offset; //header offset
33 35
34 - /*unsigned long long X(){  
35 - return R[1];  
36 - }  
37 - unsigned long long Y(){  
38 - return R[2];  
39 - }  
40 - unsigned long long Z(){  
41 - return R[0];  
42 - }*/  
43 using hsi<T>::w; //use the wavelength array in stim::hsi 36 using hsi<T>::w; //use the wavelength array in stim::hsi
44 using hsi<T>::nnz; 37 using hsi<T>::nnz;
45 using binary<T>::progress; 38 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 - }  
82 - }  
83 -  
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 - }*/  
109 - 39 +
110 public: 40 public:
111 41
112 using binary<T>::open; 42 using binary<T>::open;
@@ -199,7 +129,6 @@ public: @@ -199,7 +129,6 @@ public:
199 { 129 {
200 band_index(p, page, PROGRESS); 130 band_index(p, page, PROGRESS);
201 } 131 }
202 -  
203 return true; 132 return true;
204 } 133 }
205 134
@@ -899,7 +828,7 @@ public: @@ -899,7 +828,7 @@ public:
899 return true; 828 return true;
900 } 829 }
901 830
902 - bool unsift(std::string outfile, unsigned char* mask, unsigned long long samples, unsigned long long lines){ 831 + bool unsift(std::string outfile, unsigned char* mask, unsigned long long samples, unsigned long long lines, bool PROGRESS = false){
903 832
904 // open an output stream 833 // open an output stream
905 std::ofstream target(outfile.c_str(), std::ios::binary); 834 std::ofstream target(outfile.c_str(), std::ios::binary);
@@ -938,7 +867,7 @@ public: @@ -938,7 +867,7 @@ public:
938 target.write((char *)spectrum, B * sizeof(T)); 867 target.write((char *)spectrum, B * sizeof(T));
939 } 868 }
940 869
941 - progress = (double)x / XY * 100; 870 + if(PROGRESS) progress = (double)(x + 1) / XY * 100;
942 871
943 } 872 }
944 873
@@ -946,7 +875,7 @@ public: @@ -946,7 +875,7 @@ public:
946 target.close(); 875 target.close();
947 free(spectrum); 876 free(spectrum);
948 877
949 - progress = 100; 878 + //progress = 100;
950 879
951 return true; 880 return true;
952 881
@@ -765,13 +765,7 @@ public: @@ -765,13 +765,7 @@ public:
765 unsigned long long L = XY * sizeof(T); //size of XY plane (in bytes) 765 unsigned long long L = XY * sizeof(T); //size of XY plane (in bytes)
766 766
767 //calculate the number of pixels in the mask 767 //calculate the number of pixels in the mask
768 - unsigned long long P = 0;  
769 - if(mask == NULL) P = XY;  
770 - else{  
771 - for(unsigned long long xy = 0; xy < XY; xy++){  
772 - if(mask[xy] != 0) P++;  
773 - }  
774 - } 768 + //unsigned long long P = nnz(mask);
775 769
776 T* band_image = (T*) malloc( XY * sizeof(T)); //allocate space for a single band 770 T* band_image = (T*) malloc( XY * sizeof(T)); //allocate space for a single band
777 771
@@ -784,7 +778,9 @@ public: @@ -784,7 +778,9 @@ public:
784 if(mask == NULL || mask[xy] != 0){ //if the pixel is valid 778 if(mask == NULL || mask[xy] != 0){ //if the pixel is valid
785 matrix[i*Z() + b] = band_image[xy]; //copy it to the appropriate point in the values[] array 779 matrix[i*Z() + b] = band_image[xy]; //copy it to the appropriate point in the values[] array
786 i++; 780 i++;
  781 + //std::cout<<i<<std::endl;
787 } 782 }
  783 +
788 } 784 }
789 } 785 }
790 786
@@ -828,7 +824,7 @@ public: @@ -828,7 +824,7 @@ public:
828 } 824 }
829 825
830 /// Generates a spectral image from a matrix of spectral values in lexicographic order and a mask 826 /// Generates a spectral image from a matrix of spectral values in lexicographic order and a mask
831 - bool unsift(std::string outfile, unsigned char* p, unsigned long long samples, unsigned long long lines){ 827 + bool unsift(std::string outfile, unsigned char* p, unsigned long long samples, unsigned long long lines, bool PROGRESS = false){
832 828
833 //create a binary output stream 829 //create a binary output stream
834 std::ofstream target(outfile.c_str(), std::ios::binary); 830 std::ofstream target(outfile.c_str(), std::ios::binary);
@@ -870,7 +866,7 @@ public: @@ -870,7 +866,7 @@ public:
870 i++; 866 i++;
871 } 867 }
872 //std::cout<<xi<<"/"<<XY<<", "<<b<<"/"<<B<<std::endl; 868 //std::cout<<xi<<"/"<<XY<<", "<<b<<"/"<<B<<std::endl;
873 - progress = (double)(b * XY + xi) / (B * XY) * 100; 869 + if(PROGRESS) progress = (double)((b + 1) * XY + xi + 1) / (B * XY) * 100;
874 } 870 }
875 //std::cout<<b*XY<<"/"<<B*XY<<" "<<B<<"-"<<XY<<std::endl; 871 //std::cout<<b*XY<<"/"<<B*XY<<" "<<B<<"-"<<XY<<std::endl;
876 //write the band image to disk 872 //write the band image to disk
@@ -878,7 +874,7 @@ public: @@ -878,7 +874,7 @@ public:
878 } 874 }
879 875
880 //std::cout<<"unsifted"<<std::endl; 876 //std::cout<<"unsifted"<<std::endl;
881 - progress = 100; 877 + //progress = 100;
882 878
883 return true; 879 return true;
884 } 880 }
@@ -214,7 +214,7 @@ public: @@ -214,7 +214,7 @@ public:
214 /// @param band is the band label to be output 214 /// @param band is the band label to be output
215 /// @param threshold is a threshold value specified such that normalization will only be done to values in the band > threshold (preventing division by small numbers) 215 /// @param threshold is a threshold value specified such that normalization will only be done to values in the band > threshold (preventing division by small numbers)
216 bool normalize(std::string outfile, double band, unsigned char* mask = NULL, bool PROGRESS = false){ 216 bool normalize(std::string outfile, double band, unsigned char* mask = NULL, bool PROGRESS = false){
217 - 217 + header.save(outfile + ".hdr");
218 if(header.interleave == envi_header::BSQ){ //if the infile is bsq file 218 if(header.interleave == envi_header::BSQ){ //if the infile is bsq file
219 if(header.data_type ==envi_header::float32) 219 if(header.data_type ==envi_header::float32)
220 return ((bsq<float>*)file)->normalize(outfile, band, mask, PROGRESS); 220 return ((bsq<float>*)file)->normalize(outfile, band, mask, PROGRESS);
@@ -355,17 +355,17 @@ public: @@ -355,17 +355,17 @@ public:
355 /// @param outfile is the file name for the converted output 355 /// @param outfile is the file name for the converted output
356 /// @param interleave is the interleave format for the destination file 356 /// @param interleave is the interleave format for the destination file
357 bool convert(std::string outfile, stim::envi_header::interleaveType interleave, bool PROGRESS = false){ 357 bool convert(std::string outfile, stim::envi_header::interleaveType interleave, bool PROGRESS = false){
358 - 358 +
359 if(header.interleave == envi_header::BSQ){ //if the infile is bsq file 359 if(header.interleave == envi_header::BSQ){ //if the infile is bsq file
360 360
361 - if(header.data_type ==envi_header::float32){ //if the data type of infile is float 361 + if(header.data_type ==envi_header::float32){ //ERROR
362 if(interleave == envi_header::BSQ){ 362 if(interleave == envi_header::BSQ){
363 std::cout<<"ERROR: is already BSQ file"<<std::endl; 363 std::cout<<"ERROR: is already BSQ file"<<std::endl;
364 exit(1); 364 exit(1);
365 } 365 }
366 - else if(interleave == envi_header::BIL) //if the target file is bil file  
367 - return ((bsq<float>*)file)->bil(outfile, PROGRESS);  
368 - else if(interleave == envi_header::BIP){ //if the target file is bip file 366 + else if(interleave == envi_header::BIL) //convert BSQ -> BIL
  367 + ((bsq<float>*)file)->bil(outfile, PROGRESS);
  368 + else if(interleave == envi_header::BIP){ //ERROR
369 std::cout<<"ERROR: conversion from BSQ to BIP isn't practical; use BSQ->BIL->BIP instead"<<std::endl; 369 std::cout<<"ERROR: conversion from BSQ to BIP isn't practical; use BSQ->BIL->BIP instead"<<std::endl;
370 //return ((bsq<float>*)file)->bip(outfile, PROGRESS); 370 //return ((bsq<float>*)file)->bip(outfile, PROGRESS);
371 exit(1); 371 exit(1);
@@ -373,13 +373,13 @@ public: @@ -373,13 +373,13 @@ public:
373 } 373 }
374 374
375 else if(header.data_type == envi_header::float64){ //if the data type is float 375 else if(header.data_type == envi_header::float64){ //if the data type is float
376 - if(interleave == envi_header::BSQ){ 376 + if(interleave == envi_header::BSQ){ //ERROR
377 std::cout<<"ERROR: is already BSQ file"<<std::endl; 377 std::cout<<"ERROR: is already BSQ file"<<std::endl;
378 exit(1); 378 exit(1);
379 } 379 }
380 - else if(interleave == envi_header::BIL)  
381 - return ((bsq<double>*)file)->bil(outfile, PROGRESS);  
382 - else if(interleave == envi_header::BIP){ //if the target file is bip file 380 + else if(interleave == envi_header::BIL) //convert BSQ -> BIL
  381 + ((bsq<double>*)file)->bil(outfile, PROGRESS);
  382 + else if(interleave == envi_header::BIP){ //ERROR
383 std::cout<<"ERROR: conversion from BSQ to BIP isn't practical; use BSQ->BIL->BIP instead"<<std::endl; 383 std::cout<<"ERROR: conversion from BSQ to BIP isn't practical; use BSQ->BIL->BIP instead"<<std::endl;
384 //return ((bsq<float>*)file)->bip(outfile, PROGRESS); 384 //return ((bsq<float>*)file)->bip(outfile, PROGRESS);
385 exit(1); 385 exit(1);
@@ -394,26 +394,26 @@ public: @@ -394,26 +394,26 @@ public:
394 394
395 else if(header.interleave == envi_header::BIL){ 395 else if(header.interleave == envi_header::BIL){
396 396
397 - if(header.data_type ==envi_header::float32){ //if the data type of infile is float 397 + if(header.data_type ==envi_header::float32){ //ERROR
398 if(interleave == envi_header::BIL){ 398 if(interleave == envi_header::BIL){
399 std::cout<<"ERROR: is already BIL file"<<std::endl; 399 std::cout<<"ERROR: is already BIL file"<<std::endl;
400 exit(1); 400 exit(1);
401 } 401 }
402 - else if(interleave == envi_header::BSQ) //if the target file is bsq file  
403 - return ((bil<float>*)file)->bsq(outfile, PROGRESS);  
404 - else if(interleave == envi_header::BIP) //if the target file is bip file  
405 - return ((bil<float>*)file)->bip(outfile, PROGRESS); 402 + else if(interleave == envi_header::BSQ) //BIL -> BSQ
  403 + ((bil<float>*)file)->bsq(outfile, PROGRESS);
  404 + else if(interleave == envi_header::BIP) //BIL -> BIP
  405 + ((bil<float>*)file)->bip(outfile, PROGRESS);
406 } 406 }
407 407
408 - else if(header.data_type == envi_header::float64){ //if the data type is float  
409 - if(interleave == envi_header::BIL){ 408 + else if(header.data_type == envi_header::float64){
  409 + if(interleave == envi_header::BIL){ //ERROR
410 std::cout<<"ERROR: is already BIL file"<<std::endl; 410 std::cout<<"ERROR: is already BIL file"<<std::endl;
411 exit(1); 411 exit(1);
412 } 412 }
413 - else if(interleave == envi_header::BSQ)  
414 - return ((bil<double>*)file)->bsq(outfile, PROGRESS);  
415 - else if(interleave == envi_header::BIP)  
416 - return ((bil<double>*)file)->bip(outfile, PROGRESS); 413 + else if(interleave == envi_header::BSQ) //BIL -> BSQ
  414 + ((bil<double>*)file)->bsq(outfile, PROGRESS);
  415 + else if(interleave == envi_header::BIP) //BIL -> BIP
  416 + ((bil<double>*)file)->bip(outfile, PROGRESS);
417 } 417 }
418 418
419 else{ 419 else{
@@ -424,28 +424,28 @@ public: @@ -424,28 +424,28 @@ public:
424 424
425 else if(header.interleave == envi_header::BIP){ 425 else if(header.interleave == envi_header::BIP){
426 426
427 - if(header.data_type ==envi_header::float32){ //if the data type of infile is float  
428 - if(interleave == envi_header::BIP){ 427 + if(header.data_type ==envi_header::float32){
  428 + if(interleave == envi_header::BIP){ //ERROR
429 std::cout<<"ERROR: is already BIP file"<<std::endl; 429 std::cout<<"ERROR: is already BIP file"<<std::endl;
430 exit(1); 430 exit(1);
431 } 431 }
432 - else if(interleave == envi_header::BIL) //if the target file is bil file  
433 - return ((bip<float>*)file)->bil(outfile, PROGRESS);  
434 - else if(interleave == envi_header::BSQ){ //if the target file is bip file 432 + else if(interleave == envi_header::BIL) //BIP -> BIL
  433 + ((bip<float>*)file)->bil(outfile, PROGRESS);
  434 + else if(interleave == envi_header::BSQ){ //ERROR
435 std::cout<<"ERROR: conversion from BIP to BSQ isn't practical; use BIP->BIL->BSQ instead"<<std::endl; 435 std::cout<<"ERROR: conversion from BIP to BSQ isn't practical; use BIP->BIL->BSQ instead"<<std::endl;
436 //return ((bsq<float>*)file)->bip(outfile, PROGRESS); 436 //return ((bsq<float>*)file)->bip(outfile, PROGRESS);
437 exit(1); 437 exit(1);
438 } 438 }
439 } 439 }
440 440
441 - else if(header.data_type == envi_header::float64){ //if the data type is float  
442 - if(interleave == envi_header::BIP){ 441 + else if(header.data_type == envi_header::float64){
  442 + if(interleave == envi_header::BIP){ //ERROR
443 std::cout<<"ERROR: is already BIP file"<<std::endl; 443 std::cout<<"ERROR: is already BIP file"<<std::endl;
444 exit(1); 444 exit(1);
445 } 445 }
446 - else if(interleave == envi_header::BIL) //if the target file is bil file  
447 - return ((bip<double>*)file)->bil(outfile, PROGRESS);  
448 - else if(interleave == envi_header::BSQ){ //if the target file is bip file 446 + else if(interleave == envi_header::BIL) //BIP -> BIL
  447 + ((bip<double>*)file)->bil(outfile, PROGRESS);
  448 + else if(interleave == envi_header::BSQ){ //ERROR
449 std::cout<<"ERROR: conversion from BIP to BSQ isn't practical; use BIP->BIL->BSQ instead"<<std::endl; 449 std::cout<<"ERROR: conversion from BIP to BSQ isn't practical; use BIP->BIL->BSQ instead"<<std::endl;
450 //return ((bsq<float>*)file)->bip(outfile, PROGRESS); 450 //return ((bsq<float>*)file)->bip(outfile, PROGRESS);
451 exit(1); 451 exit(1);
@@ -462,7 +462,11 @@ public: @@ -462,7 +462,11 @@ public:
462 std::cout<<"ERROR: unidentified interleave type"<<std::endl; 462 std::cout<<"ERROR: unidentified interleave type"<<std::endl;
463 exit(1); 463 exit(1);
464 } 464 }
465 - return false; 465 + stim::envi_header h = header;
  466 + h.interleave = interleave;
  467 + h.save(outfile + ".hdr");
  468 +
  469 + return true;
466 470
467 } 471 }
468 472
@@ -509,7 +513,7 @@ public: @@ -509,7 +513,7 @@ public:
509 /// @param p is memory of size X*Y containing the mask (0 = false, all other values are true) 513 /// @param p is memory of size X*Y containing the mask (0 = false, all other values are true)
510 bool apply_mask(std::string outfile, unsigned char* p) 514 bool apply_mask(std::string outfile, unsigned char* p)
511 { 515 {
512 - 516 + header.save(outfile + ".hdr");
513 if (header.interleave == envi_header::BSQ){ //if the infile is bsq file 517 if (header.interleave == envi_header::BSQ){ //if the infile is bsq file
514 if (header.data_type == envi_header::float32) 518 if (header.data_type == envi_header::float32)
515 return ((bsq<float>*)file)->apply_mask(outfile, p); 519 return ((bsq<float>*)file)->apply_mask(outfile, p);
@@ -630,7 +634,7 @@ public: @@ -630,7 +634,7 @@ public:
630 return false; 634 return false;
631 } 635 }
632 636
633 - bool unsift(std::string outfile, unsigned char* mask, unsigned long long samples, unsigned long long lines){ 637 + bool unsift(std::string outfile, unsigned char* mask, unsigned long long samples, unsigned long long lines, bool PROGRESS = false){
634 638
635 //create a new header 639 //create a new header
636 envi_header new_header = header; 640 envi_header new_header = header;
@@ -643,9 +647,9 @@ public: @@ -643,9 +647,9 @@ public:
643 647
644 if (header.interleave == envi_header::BSQ){ //if the infile is bsq file 648 if (header.interleave == envi_header::BSQ){ //if the infile is bsq file
645 if (header.data_type == envi_header::float32) 649 if (header.data_type == envi_header::float32)
646 - return ((bsq<float>*)file)->unsift(outfile, mask, samples, lines); 650 + return ((bsq<float>*)file)->unsift(outfile, mask, samples, lines, PROGRESS);
647 else if (header.data_type == envi_header::float64) 651 else if (header.data_type == envi_header::float64)
648 - return ((bsq<double>*)file)->unsift(outfile, mask, samples, lines); 652 + return ((bsq<double>*)file)->unsift(outfile, mask, samples, lines, PROGRESS);
649 else 653 else
650 std::cout << "ERROR: unidentified data type" << std::endl; 654 std::cout << "ERROR: unidentified data type" << std::endl;
651 } 655 }
@@ -658,9 +662,9 @@ public: @@ -658,9 +662,9 @@ public:
658 else if (header.interleave == envi_header::BIP){ //if the infile is bip file 662 else if (header.interleave == envi_header::BIP){ //if the infile is bip file
659 663
660 if (header.data_type == envi_header::float32) 664 if (header.data_type == envi_header::float32)
661 - return ((bip<float>*)file)->unsift(outfile, mask, samples, lines); 665 + return ((bip<float>*)file)->unsift(outfile, mask, samples, lines, PROGRESS);
662 else if (header.data_type == envi_header::float64) 666 else if (header.data_type == envi_header::float64)
663 - return ((bip<double>*)file)->unsift(outfile, mask, samples, lines); 667 + return ((bip<double>*)file)->unsift(outfile, mask, samples, lines, PROGRESS);
664 else 668 else
665 std::cout << "ERROR: unidentified data type" << std::endl; 669 std::cout << "ERROR: unidentified data type" << std::endl;
666 } 670 }
stim/envi/envi_header.h
@@ -411,6 +411,26 @@ struct envi_header @@ -411,6 +411,26 @@ struct envi_header
411 } 411 }
412 412
413 } 413 }
  414 +
  415 + /// Convert an interleave type to a string
  416 + static std::string interleave_str(interleaveType t){
  417 + switch(t){
  418 + case stim::envi_header::BSQ:
  419 + return std::string("BSQ");
  420 + case stim::envi_header::BIL:
  421 + return std::string("BIL");
  422 + case stim::envi_header::BIP:
  423 + return std::string("BIP");
  424 + default:
  425 + std::cout<<"ERROR in stim::envi_header::typestr() - unrecognized type"<<std::endl;
  426 + exit(1);
  427 + }
  428 + }
  429 +
  430 + /// Convert the current interleave type to a string
  431 + std::string interleave_str(){
  432 + return interleave_str(interleave);
  433 + }
414 }; //end EnviHeader 434 }; //end EnviHeader
415 } 435 }
416 #endif 436 #endif
@@ -102,6 +102,17 @@ protected: @@ -102,6 +102,17 @@ protected:
102 } 102 }
103 return v; 103 return v;
104 } 104 }
  105 +
  106 + /// Returns the 1D on-disk index of a specified pixel location and band
  107 + size_t idx(size_t x, size_t y, size_t b){
  108 + size_t c[3]; //generate a coefficient list
  109 +
  110 + c[O[0]] = x; //assign the coordinates based on the coefficient order
  111 + c[O[1]] = y;
  112 + c[O[2]] = b;
  113 +
  114 + return c[2] * R[0] * R[1] + c[1] * R[0] + c[0]; //calculate and return the index (trust me this works)
  115 + }
105 }; 116 };
106 117
107 } //end namespace STIM 118 } //end namespace STIM
stim/image/image.h
@@ -5,6 +5,7 @@ @@ -5,6 +5,7 @@
5 #include <opencv2/highgui/highgui.hpp> 5 #include <opencv2/highgui/highgui.hpp>
6 #include <vector> 6 #include <vector>
7 #include <iostream> 7 #include <iostream>
  8 +#include <limits>
8 9
9 namespace stim{ 10 namespace stim{
10 /// This static class provides the STIM interface for loading, saving, and storing 2D images. 11 /// This static class provides the STIM interface for loading, saving, and storing 2D images.
@@ -19,9 +20,9 @@ class image{ @@ -19,9 +20,9 @@ class image{
19 T* img; //pointer to the image data (assumes RGB for loading/saving) 20 T* img; //pointer to the image data (assumes RGB for loading/saving)
20 size_t R[3]; 21 size_t R[3];
21 22
22 - size_t X(){ return R[1]; }  
23 - size_t Y(){ return R[2]; }  
24 - size_t C(){ return R[0]; } 23 + size_t X() const { return R[1]; }
  24 + size_t Y() const { return R[2]; }
  25 + size_t C() const { return R[0]; }
25 26
26 void init(){ //initializes all variables, assumes no memory is allocated 27 void init(){ //initializes all variables, assumes no memory is allocated
27 memset(R, 0, sizeof(size_t) * 3); //set the resolution and number of channels to zero 28 memset(R, 0, sizeof(size_t) * 3); //set the resolution and number of channels to zero
@@ -67,40 +68,87 @@ class image{ @@ -67,40 +68,87 @@ class image{
67 exit(1); 68 exit(1);
68 } 69 }
69 70
  71 + /// Returns the value for "white" based on the dynamic range (assumes white is 1.0 for floating point images)
  72 + T white(){
  73 + if(std::is_same<T, unsigned char>::value) return UCHAR_MAX;
  74 + if(std::is_same<T, unsigned short>::value) return SHRT_MAX;
  75 + if(std::is_same<T, unsigned>::value) return UINT_MAX;
  76 + if(std::is_same<T, unsigned long>::value) return ULONG_MAX;
  77 + if(std::is_same<T, unsigned long long>::value) return ULLONG_MAX;
  78 + if(std::is_same<T, float>::value) return 1.0f;
  79 + if(std::is_same<T, double>::value) return 1.0;
  80 +
  81 + std::cout<<"ERROR in stim::image::white - no white value known for this data type"<<std::endl;
  82 +
  83 + }
  84 +
70 85
71 public: 86 public:
72 87
73 - //default constructor 88 + /// Default constructor - creates an empty image object
74 image(){ init(); } //initialize all variables to zero, don't allocate any memory 89 image(){ init(); } //initialize all variables to zero, don't allocate any memory
75 90
  91 + /// Constructor with a filename - loads the specified file
76 image(std::string filename){ //constructor initialize the image with an image file 92 image(std::string filename){ //constructor initialize the image with an image file
77 load(filename); 93 load(filename);
78 } 94 }
79 95
80 /// Create a new image from scratch given a number of samples and channels 96 /// Create a new image from scratch given a number of samples and channels
81 image(size_t x, size_t y = 1, size_t c = 1){ 97 image(size_t x, size_t y = 1, size_t c = 1){
82 - img = (T*) malloc( sizeof(T) * x * y * c ); 98 + allocate(x, y, c);
83 } 99 }
84 100
  101 + /// Create a new image with the data given in 'data'
85 image(T* data, size_t x, size_t y, size_t c = 1){ 102 image(T* data, size_t x, size_t y, size_t c = 1){
86 allocate(x, y, c); 103 allocate(x, y, c);
87 memcpy(img, data, bytes()); 104 memcpy(img, data, bytes());
88 } 105 }
89 106
  107 + /// Copy constructor - duplicates an image object
  108 + image(const stim::image<T>& I){
  109 + allocate(I.X(), I.Y(), I.C());
  110 + //allocate(I.R[1], I.R[2], I.R[0]);
  111 + memcpy(img, I.img, bytes());
  112 + }
  113 +
  114 + /// Destructor - clear memory
  115 + ~image(){
  116 + free(img);
  117 + }
  118 +
  119 + stim::image<T> operator=(const stim::image<T>& I){
  120 + if(&I == this) //handle self-assignment
  121 + return *this;
  122 + allocate(I.X(), I.Y(), I.C());
  123 + memcpy(img, I.img, bytes());
  124 + return *this;
  125 + }
  126 +
90 /// Load an image from a file 127 /// Load an image from a file
91 void load(std::string filename){ 128 void load(std::string filename){
92 129
93 cv::Mat cvImage = cv::imread(filename, CV_LOAD_IMAGE_UNCHANGED); //use OpenCV to open the image file 130 cv::Mat cvImage = cv::imread(filename, CV_LOAD_IMAGE_UNCHANGED); //use OpenCV to open the image file
  131 + if(!cvImage.data){
  132 + std::cout<<"ERROR stim::image::load() - unable to find image "<<filename<<std::endl;
  133 + exit(1);
  134 + }
94 allocate(cvImage.cols, cvImage.rows, cvImage.channels()); //allocate space for the image 135 allocate(cvImage.cols, cvImage.rows, cvImage.channels()); //allocate space for the image
95 T* cv_ptr = (T*)cvImage.data; 136 T* cv_ptr = (T*)cvImage.data;
96 - set_interleaved_bgr(cv_ptr, X(), Y()); 137 + if(C() == 1) //if this is a single-color image, just copy the data
  138 + memcpy(img, cv_ptr, bytes());
  139 + if(C() == 3) //if this is a 3-color image, OpenCV uses BGR interleaving
  140 + set_interleaved_bgr(cv_ptr, X(), Y());
97 } 141 }
98 142
99 //save a file 143 //save a file
100 void save(std::string filename){ 144 void save(std::string filename){
101 //OpenCV uses an interleaved format, so convert first and then output 145 //OpenCV uses an interleaved format, so convert first and then output
102 T* buffer = (T*) malloc(bytes()); 146 T* buffer = (T*) malloc(bytes());
103 - data_interleaved_bgr(buffer); 147 +
  148 + if(C() == 1)
  149 + memcpy(buffer, img, bytes());
  150 + else if(C() == 3)
  151 + data_interleaved_bgr(buffer);
104 cv::Mat cvImage((int)Y(), (int)X(), cv_type(), buffer); 152 cv::Mat cvImage((int)Y(), (int)X(), cv_type(), buffer);
105 cv::imwrite(filename, cvImage); 153 cv::imwrite(filename, cvImage);
106 } 154 }
@@ -137,7 +185,7 @@ public: @@ -137,7 +185,7 @@ public:
137 image<T> channel(size_t c){ 185 image<T> channel(size_t c){
138 186
139 //create a new image 187 //create a new image
140 - image<T> r(R[0], R[1], 1); 188 + image<T> r(X(), Y(), 1);
141 189
142 for(size_t x = 0; x < X(); x++){ 190 for(size_t x = 0; x < X(); x++){
143 for(size_t y = 0; y < Y(); y++){ 191 for(size_t y = 0; y < Y(); y++){
@@ -266,6 +314,16 @@ public: @@ -266,6 +314,16 @@ public:
266 return max; 314 return max;
267 } 315 }
268 316
  317 + /// Invert an image by calculating I1 = alpha - I0, where alpha is the maximum image value
  318 + image<T> invert(T white_val){
  319 + size_t N = size(); //calculate the total number of values in the image
  320 + image<T> r(X(), Y(), C()); //allocate space for the resulting image
  321 + for(size_t n = 0; n < N; n++)
  322 + r.img[n] = white_val - img[n]; //perform the inversion
  323 +
  324 + return r; //return the inverted image
  325 + }
  326 +
269 image<T> srgb2lab(){ 327 image<T> srgb2lab(){
270 std::cout<<"ERROR stim::image::srgb2lab - function has been broken, re-implement."<<std::endl; 328 std::cout<<"ERROR stim::image::srgb2lab - function has been broken, re-implement."<<std::endl;
271 exit(1); 329 exit(1);
stim/parser/arguments.h
@@ -82,7 +82,7 @@ namespace stim{ @@ -82,7 +82,7 @@ namespace stim{
82 } 82 }
83 83
84 //return the value of a text option 84 //return the value of a text option
85 - std::string as_string(unsigned int n = 0) 85 + std::string as_string(size_t n = 0)
86 { 86 {
87 if(!flag) 87 if(!flag)
88 { 88 {
@@ -97,7 +97,7 @@ namespace stim{ @@ -97,7 +97,7 @@ namespace stim{
97 } 97 }
98 98
99 //return the value of a floating point option 99 //return the value of a floating point option
100 - float as_float(unsigned int n = 0) 100 + float as_float(size_t n = 0)
101 { 101 {
102 if(!flag) 102 if(!flag)
103 { 103 {
@@ -116,7 +116,7 @@ namespace stim{ @@ -116,7 +116,7 @@ namespace stim{
116 } 116 }
117 117
118 //return the value of an integer option 118 //return the value of an integer option
119 - int as_int(unsigned int n = 0) 119 + int as_int(size_t n = 0)
120 { 120 {
121 if(!flag) 121 if(!flag)
122 { 122 {
stim/parser/filename.h
@@ -299,9 +299,9 @@ public: @@ -299,9 +299,9 @@ public:
299 stim::filename result = *this; //initialize the result filename to the current filename 299 stim::filename result = *this; //initialize the result filename to the current filename
300 size_t loc = result.prefix.find('*'); //find a wild card in the string 300 size_t loc = result.prefix.find('*'); //find a wild card in the string
301 if(loc == std::string::npos) //if a wildcard isn't found 301 if(loc == std::string::npos) //if a wildcard isn't found
302 - return result; //return the original string  
303 -  
304 - result.prefix.replace(loc, 1, str); //replace the wildcard with the string 302 + result.prefix += str; //append the value to the prefix
  303 + else
  304 + result.prefix.replace(loc, 1, str); //replace the wildcard with the string
305 return result; //return the result 305 return result; //return the result
306 } 306 }
307 307
stim/visualization/colormap.h
@@ -3,6 +3,12 @@ @@ -3,6 +3,12 @@
3 3
4 #include <string> 4 #include <string>
5 #include <stdlib.h> 5 #include <stdlib.h>
  6 +#include <cmath>
  7 +
  8 +#ifdef _WIN32
  9 + #include <float.h>
  10 +#endif
  11 +
6 #ifdef __CUDACC__ 12 #ifdef __CUDACC__
7 #include <stim/cuda/cudatools/error.h> 13 #include <stim/cuda/cudatools/error.h>
8 #endif 14 #endif
@@ -226,8 +232,13 @@ static void cpuApplyBrewer(T* cpuSource, unsigned char* cpuDest, size_t N, T min @@ -226,8 +232,13 @@ static void cpuApplyBrewer(T* cpuSource, unsigned char* cpuDest, size_t N, T min
226 a = (cpuSource[i] - minVal) / (maxVal - minVal); 232 a = (cpuSource[i] - minVal) / (maxVal - minVal);
227 else 233 else
228 a = 0.5; 234 a = 0.5;
229 - if(a < 0) a = 0;  
230 - if(a > 1) a = 1; 235 +#ifdef _WIN32
  236 + if(!_finite(a)) a = 1; //deal with infinite and NaN values (return maximum in all cases)
  237 +#else
  238 + if(!std::isfinite(a)) a = 1;
  239 +#endif
  240 + else if(a < 0) a = 0;
  241 + else if(a > 1) a = 1;
231 242
232 float c = a * (float)(BREWER_CTRL_PTS-1); 243 float c = a * (float)(BREWER_CTRL_PTS-1);
233 int ptLow = (int)c; 244 int ptLow = (int)c;