From 63fc1df8ebee680b2a16a6709d79f05e2dffe90b Mon Sep 17 00:00:00 2001 From: David Date: Sun, 10 Apr 2016 18:20:05 -0500 Subject: [PATCH] added an inverse PCA transform to the stim::envi classes --- stim/envi/bil.h | 70 +++++++++++++++++++++++++++++++++------------------------------------- stim/envi/binary.h | 73 +++++++++++++++++++++++++++++++++++++++++++++++++------------------------ stim/envi/bip.h | 434 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- stim/envi/bsq.h | 71 +++++++++++++++++++++++++++++++++-------------------------------------- stim/envi/envi.h | 228 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------------------------------------------------------------------------- 5 files changed, 491 insertions(+), 385 deletions(-) diff --git a/stim/envi/bil.h b/stim/envi/bil.h index 6999319..f80f909 100644 --- a/stim/envi/bil.h +++ b/stim/envi/bil.h @@ -35,7 +35,12 @@ protected: return R[1]; } - using binary::thread_data; + using binary::progress; + + /// Call the binary nnz() function for the BIL orientation + unsigned long long nnz(unsigned char* mask){ + return binary::nnz(mask, X()*Y()); + } public: @@ -63,9 +68,10 @@ public: /// @param p is a pointer to an allocated region of memory at least X * Y * sizeof(T) in size. /// @param page <= B is the integer number of the band to be copied. - bool band_index( T * p, unsigned int page){ + bool band_index( T * p, unsigned int page, bool PROGRESS = false){ //return binary::read_plane_1(p, page); + if(PROGRESS) progress = 0; unsigned int L = X() * sizeof(T); //caculate the number of bytes in a sample line unsigned int jump = X() * (Z() - 1) * sizeof(T); @@ -79,6 +85,7 @@ public: { file.read((char *)(p + i * X()), L); file.seekg( jump, std::ios::cur); + if(PROGRESS) progress = (double)(i+1) / Y() * 100; } return true; @@ -88,7 +95,7 @@ public: /// @param p is a pointer to an allocated region of memory at least X * Y * sizeof(T) in size. /// @param wavelength is a floating point value (usually a wavelength in spectral data) used as a label for the band to be copied. - bool band( T * p, double wavelength){ + bool band( T * p, double wavelength, bool PROGRESS = false){ //if there are no wavelengths in the BSQ file if(w.size() == 0) @@ -104,7 +111,7 @@ public: //if wavelength is smaller than the first one in header file if ( w[page] > wavelength ){ - band_index(p, page); + band_index(p, page, PROGRESS); return true; } @@ -124,7 +131,7 @@ public: p1=(T*)malloc(S); //memory allocation p2=(T*)malloc(S); band_index(p1, page - 1); - band_index(p2, page ); + band_index(p2, page, PROGRESS); for(unsigned i=0; i < XY; i++){ double r = (double) (wavelength - w[page-1]) / (double) (w[page] - w[page-1]); p[i] = (p2[i] - p1[i]) * r + p1[i]; @@ -134,7 +141,7 @@ public: } else //if the wavelength is equal to a wavelength in header file { - band_index(p, page); + band_index(p, page, PROGRESS); } return true; @@ -195,8 +202,8 @@ public: /// @param p is a pointer to pre-allocated memory at least B * sizeof(T) in size. /// @param x is the x-coordinate (dimension 1) of the spectrum. /// @param y is the y-coordinate (dimension 2) of the spectrum. - bool spectrum(T * p, unsigned x, unsigned y){ - return binary::read_line_02(p, x, y); + bool spectrum(T * p, unsigned x, unsigned y, bool PROGRESS = false){ + return binary::read_line_1(p, x, y, PROGRESS); } /// Retrieve a single pixel and stores it in pre-allocated memory. @@ -223,7 +230,7 @@ public: /// @param outname is the name of the output file used to store the resulting baseline-corrected data. /// @param wls is the list of baseline points based on band labels. - bool baseline(std::string outname, std::vector wls){ + bool baseline(std::string outname, std::vector wls, bool PROGRESS = false){ unsigned N = wls.size(); //get the number of baseline points @@ -321,7 +328,7 @@ public: }//loop for YZ line end target.write(reinterpret_cast(c), L); //write the corrected data into destination - thread_data = (double)k / Y() * 100; + if(PROGRESS) progress = (double)(k+1) / Y() * 100; }//loop for Y slice end free(a); @@ -329,8 +336,6 @@ public: free(c); target.close(); - thread_data = 100; - return true; } @@ -340,7 +345,7 @@ public: /// @param outname is the name of the output file used to store the resulting baseline-corrected data. /// @param w is the label specifying the band that the hyperspectral image will be normalized to. /// @param t is a threshold specified such that a spectrum with a value at w less than t is set to zero. Setting this threshold allows the user to limit division by extremely small numbers. - bool normalize(std::string outname, double w, double t = 0.0) + bool normalize(std::string outname, double w, double t = 0.0, bool PROGRESS = false) { unsigned int B = Z(); //calculate the number of bands unsigned int ZX = Z() * X(); @@ -374,20 +379,19 @@ public: } target.write(reinterpret_cast(c), L); //write normalized data into destination - thread_data = (double)j / Y() * 100; + if(PROGRESS) progress = (double)(j+1) / Y() * 100; } free(b); free(c); target.close(); - thread_data = 100; return true; } /// Convert the current BIL file to a BSQ file with the specified file name. /// @param outname is the name of the output BSQ file to be saved to disk. - bool bsq(std::string outname) + bool bsq(std::string outname, bool PROGRESS = false) { unsigned int S = X() * Y() * sizeof(T); //calculate the number of bytes in a band @@ -402,11 +406,9 @@ public: band_index(p, i); target.write(reinterpret_cast(p), S); //write a band data into target file - thread_data = (double)i / Z() * 100; //store the progress for the current operation + if(PROGRESS) progress = (double)(i+1) / Z() * 100; //store the progress for the current operation } - thread_data = 100; //set the current progress to 100% - free(p); target.close(); return true; @@ -415,7 +417,7 @@ public: /// Convert the current BIL file to a BIP file with the specified file name. /// @param outname is the name of the output BIP file to be saved to disk. - bool bip(std::string outname) + bool bip(std::string outname, bool PROGRESS = false) { unsigned int S = X() * Z() * sizeof(T); //calculate the number of bytes in a ZX slice @@ -436,15 +438,12 @@ public: for ( unsigned j = 0; j < X(); j++) q[k + j * Z()] = p[ks + j]; - thread_data = (double)(i * Z() + k) / (Z() * Y()) * 100; //store the progress for the current operation + if(PROGRESS) progress = (double)((i+1) * Z() + k+1) / (Z() * Y()) * 100; //store the progress for the current operation } target.write(reinterpret_cast(q), S); //write a band data into target file } - thread_data = 100; - - free(p); free(q); target.close(); @@ -879,14 +878,14 @@ public: target.write((char*)spec, B * sizeof(T)); //write that spectrum to disk. Size is L2. } - thread_data = (double) (y * X() + x) / (Y() * X()) * 100; + progress = (double) (y * X() + x) / (Y() * X()) * 100; } } target.close(); free(slice); free(spec); - thread_data = 100; + progress = 100; return true; } @@ -923,7 +922,7 @@ public: /// @param p is a pointer to pre-allocated memory of size [B * sizeof(T)] that stores the mean spectrum /// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location - bool avg_band(double* p, unsigned char* mask = NULL){ + bool avg_band(double* p, unsigned char* mask = NULL, bool PROGRESS = false){ unsigned long long XZ = X() * Z(); unsigned long long XY = X() * Y(); T* temp = (T*)malloc(sizeof(T) * XZ); @@ -947,6 +946,7 @@ public: } } } + if(PROGRESS) progress = (double)(k+1) / Y() * 100; } free(temp); return true; @@ -957,8 +957,8 @@ public: /// @param co is a pointer to pre-allocated memory of size [B * B] that stores the resulting covariance matrix /// @param avg is a pointer to memory of size B that stores the average spectrum /// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location - bool co_matrix(double* co, double* avg, unsigned char *mask){ - thread_data = 0; + bool co_matrix(double* co, double* avg, unsigned char *mask, bool PROGRESS = false){ + progress = 0; //memory allocation unsigned long long xy = X() * Y(); unsigned long long B = Z(); @@ -986,7 +986,7 @@ public: } } } - thread_data = (double)j / xy * 100; + if(PROGRESS) progress = (double)(j+1) / xy * 100; } //because correlation matrix is symmetric for (unsigned long long i = 0; i < B; i++){ @@ -996,9 +996,6 @@ public: } free(temp); - - thread_data = 100; //processing complete - return true; } @@ -1015,7 +1012,8 @@ public: unsigned long long x1, unsigned long long y1, unsigned long long b0, - unsigned long long b1){ + unsigned long long b1, + bool PROGRESS = false){ //calculate the new image parameters unsigned long long samples = x1 - x0; @@ -1048,7 +1046,7 @@ public: file.read((char *)(temp + z * samples), sizeof(T) * samples); file.seekg(jumpb, std::ios::cur); //go to the next band - thread_data = (double)(x * Z() + z) / (lines * Z()) * 100; + if(PROGRESS) progress = (double)(x * Z() + z+1) / (lines * Z()) * 100; } //write slice data into target file @@ -1061,8 +1059,6 @@ public: //free the temporary frame free(temp); - thread_data = 100; - return true; } diff --git a/stim/envi/binary.h b/stim/envi/binary.h index 3989ed2..e010b50 100644 --- a/stim/envi/binary.h +++ b/stim/envi/binary.h @@ -24,11 +24,11 @@ protected: std::fstream file; //file stream used for reading and writing std::string name; //file name - unsigned long long int R[D]; //resolution + unsigned long long R[D]; //resolution unsigned int header; //header size (in bytes) unsigned char* mask; //pointer to a character array: 0 = background, 1 = foreground (or valid data) - unsigned int thread_data; //unsigned integer used to pass data to threads during processing + double progress; //stores the progress on the current operation (accessible using a thread) /// Private initialization function used to set default parameters in the data structure. @@ -37,7 +37,24 @@ protected: header = 0; //initialize the header size to zero mask = NULL; - thread_data = 0; + progress = 0; + } + + //calculate the number of non-zero pixels in a mask + unsigned long long nnz(unsigned char* mask, unsigned long long N){ + + unsigned long long n = 0; //initialize the counter to 0 (zero) + if(mask == NULL) return N; //if the mask is NULL, assume all pixels are masked + + for(unsigned long long i = 0; i < N; i++){ //for each pixel + if(mask[i] != 0) n++; //increment the counter for every non-zero pixel in the mask + } + return n; //return the number of nonzero pixels + } + + /// Calculate the number of nonzero pixels in a mask over X-Y + unsigned long long nnz(unsigned char* mask){ + return nnz(mask, R[0] * R[1]); } /// Private helper function that returns the size of the file on disk using system functions. @@ -105,12 +122,12 @@ protected: public: - unsigned int get_thread_data(){ - return thread_data; + unsigned int get_progress(){ + return progress; } - void reset_thread_data(){ - thread_data = 0; + void reset_progress(){ + progress = 0; } /// Open a binary file for streaming. @@ -178,7 +195,9 @@ public: /// @param p is a pointer to pre-allocated memory equal to the page size /// @param page is the index of the page - bool read_page( T * p, unsigned int page){ + bool read_page( T * p, unsigned int page, bool PROGRESS = false){ + + if(PROGRESS) progress = 0; if (page >= R[2]){ //make sure the bank number is right std::cout<<"ERROR: page out of range"<= R[0] || y >= R[1]){ //make sure the sample and line number is right std::cout<<"ERROR: sample or line out of range"<= R[1] || z >= R[2] ){ std::cout<<"ERROR: sample or line out of range"<= R[0] || z >= R[2] ){ std::cout<<"ERROR: sample or line out of range"<= R[0]){ //make sure the number is within the possible range std::cout<<"ERROR read_plane_0: page out of range"<::thread_data; + using binary::progress; + + /// Call the binary nnz() function for the BIP orientation + unsigned long long nnz(unsigned char* mask){ + return binary::nnz(mask, X()*Y()); + } + + /// Linear interpolation of a spectral value given the bounding spectral samples + T lerp(double w, T low_v, double low_w, T high_v, double high_w){ + if(low_w == high_w) return low_v; //if the interval is of zero length, just return one of the bounds + double alpha = (w - low_w) / (high_w - low_w); //calculate the interpolation factor + return (1.0 - alpha) * low_v + alpha * high_v; //interpolate + } + + /// Gets the two band indices surrounding a given wavelength + void band_bounds(double wavelength, unsigned long long& low, unsigned long long& high){ + unsigned long long B = Z(); + for(high = 0; high < B; high++){ + if(w[high] > wavelength) break; + } + low = 0; + if(high > 0) + low = high-1; + } + + /// Get the list of band numbers that bound a list of wavelengths + void band_bounds(std::vector wavelengths, + std::vector& low_bands, + std::vector& high_bands){ + + unsigned long long W = w.size(); //get the number of wavelengths in the list + low_bands.resize(W); //pre-allocate space for the band lists + high_bands.resize(W); + + for(unsigned long long wl = 0; wl < W; wl++){ //for each wavelength + band_bounds(wavelengths[wl], low_bands[wl], high_bands[wl]); //find the low and high bands + } + } + + /// Returns the interpolated in the given spectrum based on the given wavelength + + /// @param s is the spectrum in main memory of length Z() + /// @param wavelength is the wavelength value to interpolate out + T interp_spectrum(T* s, double wavelength){ + unsigned long long low, high; //indices for the bands surrounding wavelength + band_bounds(wavelength, low, high); //get the surrounding band indices + + if(high == w.size()) return s[w.size()-1]; //if the high band is above the wavelength range, return the highest wavelength value + + return lerp(wavelength, s[low], w[low], s[high], w[high]); + } + + /// Returns the interpolated value corresponding to the given list of wavelengths + std::vector interp_spectrum(T* s, std::vector wavelengths){ + + unsigned long long N = wavelengths.size(); //get the number of wavelength measurements + + std::vector v; //allocate space for the resulting values + v.resize(wavelengths.size()); + for(unsigned long long n = 0; n < N; n++){ //for each measurement + v[n] = interp_spectrum(s, wavelengths[n]); //interpolate the measurement + } + return v; + } public: using binary::open; using binary::file; using binary::R; - using binary::read_line_12; + using binary::read_line_0; /// Open a data file for reading using the class interface. @@ -73,19 +136,19 @@ public: /// @param p is a pointer to an allocated region of memory at least X * Y * sizeof(T) in size. /// @param page <= B is the integer number of the band to be copied. - bool band_index( T * p, unsigned int page){ - return binary::read_plane_0(p, page); + bool band_index( T * p, unsigned int page, bool PROGRESS = false){ + return binary::read_plane_0(p, page, PROGRESS); } /// Retrieve a single band (by numerical label) and stores it in pre-allocated memory. /// @param p is a pointer to an allocated region of memory at least X * Y * sizeof(T) in size. /// @param wavelength is a floating point value (usually a wavelength in spectral data) used as a label for the band to be copied. - bool band( T * p, double wavelength){ + bool band( T * p, double wavelength, bool PROGRESS = false){ //if there are no wavelengths in the BSQ file if(w.size() == 0) - return band_index(p, (unsigned int)wavelength); + return band_index(p, (unsigned int)wavelength, PROGRESS); unsigned int XY = X() * Y(); //calculate the number of pixels in a band @@ -96,7 +159,7 @@ public: //if wavelength is smaller than the first one in header file if ( w[page] > wavelength ){ - band_index(p, page); + band_index(p, page, PROGRESS); return true; } @@ -105,7 +168,7 @@ public: page++; //if wavelength is larger than the last wavelength in header file if (page == Z()) { - band_index(p, Z()-1); + band_index(p, Z()-1, PROGRESS); return true; } } @@ -116,7 +179,7 @@ public: p1=(T*)malloc( XY * sizeof(T)); //memory allocation p2=(T*)malloc( XY * sizeof(T)); band_index(p1, page - 1); - band_index(p2, page ); + band_index(p2, page, PROGRESS); for(unsigned i=0; i < XY; i++){ double r = (double) (wavelength - w[page-1]) / (double) (w[page] - w[page-1]); p[i] = (p2[i] - p1[i]) * r + p1[i]; @@ -126,7 +189,7 @@ public: } else //if the wavelength is equal to a wavelength in header file { - band_index(p, page); + band_index(p, page, PROGRESS); } return true; @@ -137,8 +200,8 @@ public: /// @param p is a pointer to pre-allocated memory at least B * sizeof(T) in size. /// @param x is the x-coordinate (dimension 1) of the spectrum. /// @param y is the y-coordinate (dimension 2) of the spectrum. - bool spectrum(T * p, unsigned x, unsigned y){ - return read_line_12(p, x, y); //read a line in the binary YZ plane (dimension order for BIP is ZXY) + bool spectrum(T * p, unsigned x, unsigned y, bool PROGRESS = false){ + return read_line_0(p, x, y, PROGRESS); //read a line in the binary YZ plane (dimension order for BIP is ZXY) } /// Retrieves a band of x values from a given xz plane. @@ -235,8 +298,8 @@ public: /// @param n is an integer index to the pixel using linear array indexing. bool pixel(T * p, unsigned n){ - unsigned bandnum = X() * Y(); //calculate numbers in one band - if ( n >= bandnum){ //make sure the pixel number is right + unsigned N = X() * Y(); //calculate numbers in one band + if ( n >= N){ //make sure the pixel number is right std::cout<<"ERROR: sample or line out of range"< wls){ - - unsigned N = wls.size(); //get the number of baseline points + bool baseline(std::string outname, std::vector base_pts, bool PROGRESS = false){ std::ofstream target(outname.c_str(), std::ios::binary); //open the target binary file std::string headername = outname + ".hdr"; //the header file name - //simplify image resolution - unsigned int ZX = Z() * X(); //calculate the number of points in a Y slice - unsigned int L = ZX * sizeof(T); //calculate the number of bytes of a Y slice - unsigned int B = Z(); - - T* c; //pointer to the current Y slice - c = (T*)malloc(L); //memory allocation - - T* a; //pointer to the two YZ lines surrounding the current YZ line - T* b; - - a = (T*)malloc(X() * sizeof(T)); - b = (T*)malloc(X() * sizeof(T)); - - - double ai, bi; //stores the two baseline points wavelength surrounding the current band - double ci; //stores the current band's wavelength - unsigned control; - - if (a == NULL || b == NULL || c == NULL){ - std::cout<<"ERROR: error allocating memory"; - exit(1); - } - // loop start correct every y slice - - for (unsigned k =0; k < Y(); k++) - { - //get the current y slice - read_plane_y(c, k); - - //initialize lownum, highnum, low, high - control=0; - ai = w[0]; - //if no baseline point is specified at band 0, - //set the baseline point at band 0 to 0 - if(wls[0] != w[0]){ - bi = wls[control]; - memset(a, (char)0, X() * sizeof(T) ); - } - //else get the low band - else{ - control++; - read_x_from_xz(a, c, ai); - bi = wls[control]; - } - //get the high band - read_x_from_xz(b, c, bi); - - //correct every YZ line - - for(unsigned cii = 0; cii < B; cii++){ - //update baseline points, if necessary - if( w[cii] >= bi && cii != B - 1) { - //if the high band is now on the last BL point? - if (control != N-1) { - - control++; //increment the index - - std::swap(a, b); //swap the baseline band pointers - - ai = bi; - bi = wls[control]; - read_x_from_xz(b, c, bi); - - } - //if the last BL point on the last band of the file? - else if ( wls[control] < w[B - 1]) { - - std::swap(a, b); //swap the baseline band pointers - - memset(b, (char)0, X() * sizeof(T) ); //clear the high band - - ai = bi; - bi = w[B - 1]; - } + unsigned long long N = X() * Y(); //calculate the total number of pixels to be processed + unsigned long long B = Z(); //get the number of bands + T* s = (T*)malloc(sizeof(T) * B); //allocate memory to store a pixel + T* sbc = (T*)malloc(sizeof(T) * B); //allocate memory to store the baseline corrected spectrum + + std::vector base_vals; //allocate space for the values at each baseline point + double aw, bw; //surrounding baseline point wavelengths + T av, bv; //surrounding baseline point values + double alpha; + unsigned long long ai, bi; //surrounding baseline point band indices + for(unsigned long long n = 0; n < N; n++){ //for each pixel in the image + pixel(s, n); //retrieve the spectrum s + base_vals = interp_spectrum(s, base_pts); //get the values at each baseline point + + ai = bi = 0; + aw = w[0]; //initialize the current baseline points (assume the spectrum starts at 0) + av = s[0]; + bw = base_pts[0]; + for(unsigned long long b = 0; b < B; b++){ //for each band in the spectrum + while(bi < base_pts.size() && base_pts[bi] < w[b]) //while the current wavelength is higher than the second baseline point + bi++; //move to the next baseline point + if(bi < base_pts.size()){ + bw = base_pts[bi]; //set the wavelength for the upper bound baseline point + bv = base_vals[bi]; //set the value for the upper bound baseline point } - - ci = w[cii]; - - //perform the baseline correction - for(unsigned i=0; i < X(); i++) - { - double r = (double) (ci - ai) / (double) (bi - ai); - c[i * B + cii] =(T) ( c[i * B + cii] - (b[i] - a[i]) * r - a[i] ); + if(bi == base_pts.size()){ //if we have passed the last baseline point + bw = w[B-1]; //set the outer bound to the last spectral band + bv = s[B-1]; } - - }//loop for YZ line end - target.write(reinterpret_cast(c), L); //write the corrected data into destination - - thread_data = (double)k / Y() * 100; - }//loop for Y slice end + if(bi != 0){ + ai = bi - 1; //set the lower bound baseline point index + aw = base_pts[ai]; //set the wavelength for the lower bound baseline point + av = base_vals[ai]; //set the value for the lower bound baseline point + } + sbc[b] = s[b] - lerp(w[b], av, aw, bv, bw); //perform the baseline correction and save the new value + } + if(PROGRESS) progress = (double)(n+1) / N * 100; //set the current progress + target.write((char*)sbc, sizeof(T) * B); //write the corrected data into destination + } //end for each pixel - free(a); - free(b); - free(c); + free(s); + free(sbc); target.close(); - - thread_data = 100; + return true; } @@ -371,80 +378,41 @@ public: /// @param outname is the name of the output file used to store the resulting baseline-corrected data. /// @param w is the label specifying the band that the hyperspectral image will be normalized to. /// @param t is a threshold specified such that a spectrum with a value at w less than t is set to zero. Setting this threshold allows the user to limit division by extremely small numbers. - bool normalize(std::string outname, double w, double t = 0.0) + bool normalize(std::string outname, double w, double t = 0.0, bool PROGRESS = false) { - unsigned int B = Z(); //calculate the number of bands - unsigned int ZX = Z() * X(); - unsigned int XY = X() * Y(); //calculate the number of pixels in a band - unsigned int S = XY * sizeof(T); //calculate the number of bytes in a band - unsigned int L = ZX * sizeof(T); - std::ofstream target(outname.c_str(), std::ios::binary); //open the target binary file - std::string headername = outname + ".hdr"; //the header file name - - T * c; //pointer to the current ZX slice - T * b; //pointer to the standard band - - b = (T*)malloc( S ); //memory allocation - c = (T*)malloc( L ); - - band(b, w); //get the certain band into memory - - for(unsigned j = 0; j < Y(); j++) - { - read_plane_y(c, j); - unsigned jX = j * X(); //to avoid calculating it many times - for(unsigned i = 0; i < X(); i++) - { - unsigned iB = i * B; - for(unsigned m = 0; m < B; m++) - { - if( b[i+jX] < t ) - c[m + iB] = (T)0.0; - else - c[m + iB] = c[m + iB] / b[i + jX]; //perform normalization - } + std::string headername = outname + ".hdr"; //the header file name + + unsigned long long N = X() * Y(); //calculate the total number of pixels to be processed + unsigned long long B = Z(); //get the number of bands + T* s = (T*)malloc(sizeof(T) * B); //allocate memory to store a pixel + T nv; //stores the value of the normalized band + for(unsigned long long n = 0; n < N; n++){ //for each pixel in the image + pixel(s, n); //retrieve the spectrum s + nv = interp_spectrum(s, w); //find the value of the normalization band + if(abs(nv) <= t) //if the normalization band is below threshold + memset(s, 0, sizeof(T) * B); //set the output to zero + else{ + for(unsigned long long b = 0; b < B; b++){ //for each band in the spectrum + s[b] /= nv; //divide by the normalization value + } } - target.write(reinterpret_cast(c), L); //write normalized data into destination - - thread_data = (double) j / Y() * 100; - } + + if(PROGRESS) progress = (double)(n+1) / N * 100; //set the current progress + target.write((char*)s, sizeof(T) * B); //write the corrected data into destination + } //end for each pixel - free(b); - free(c); + free(s); target.close(); - thread_data = 100; return true; } - /// Convert the current BIP file to a BSQ file with the specified file name. - - /// @param outname is the name of the output BSQ file to be saved to disk. - bool bsq(std::string outname) - { - std::string temp = outname + "_temp"; - std::string headtemp = temp + ".hdr"; - //first creat a temporary bil file and convert bip file to bil file - bil(temp); - - stim::bil n; - if(n.open(temp, X(), Y(), Z(), offset, w)==false){ //open infile - std::cout<<"ERROR: unable to open input file"<(q), S); //write a band data into target file } - thread_data = 100; - free(p); free(q); target.close(); @@ -909,7 +875,7 @@ public: //write this pixel out target.write((char *)spectrum, B * sizeof(T)); - thread_data = (double) x / XY * 100; + progress = (double) x / XY * 100; } } @@ -918,7 +884,7 @@ public: target.close(); free(spectrum); - thread_data = 100; + progress = 100; return true; } @@ -962,7 +928,7 @@ public: target.write((char *)spectrum, B * sizeof(T)); } - thread_data = (double)x / XY * 100; + progress = (double)x / XY * 100; } @@ -970,7 +936,7 @@ public: target.close(); free(spectrum); - thread_data = 100; + progress = 100; return true; @@ -981,7 +947,7 @@ public: /// @param p is a pointer to pre-allocated memory of size [B * sizeof(T)] that stores the mean spectrum /// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location - bool avg_band(double* p, unsigned char* mask = NULL){ + bool avg_band(double* p, unsigned char* mask = NULL, bool PROGRESS = false){ unsigned long long XY = X() * Y(); T* temp = (T*)malloc(sizeof(T) * Z()); //Iinitialize @@ -989,12 +955,8 @@ public: p[j] = 0; } //calculate vaild number in a band - unsigned count = 0; - for (unsigned j = 0; j < XY; j++){ - if (mask == NULL || mask[j] != 0){ - count++; - } - } + unsigned count = nnz(mask); + //calculate average number of a band for (unsigned i = 0; i < XY; i++){ if (mask == NULL || mask[i] != 0){ @@ -1003,22 +965,21 @@ public: p[j] += (double)temp[j] / (double)count; } } - thread_data = (double)i / XY * 100; + if(PROGRESS) progress = (double)(i+1) / XY * 100; } free(temp); - thread_data = 100; return true; } #ifdef CUDA_FOUND /// Calculate the covariance matrix for masked pixels using cuBLAS - bool cu_co_matrix(double* co, double* avg, unsigned char *mask){ + bool co_matrix_cublas(double* co, double* avg, unsigned char *mask, bool PROGRESS = false){ cudaError_t cudaStat; cublasStatus_t stat; cublasHandle_t handle; - thread_data = 0; //initialize the progress to zero (0) + progress = 0; //initialize the progress to zero (0) unsigned long long XY = X() * Y(); //calculate the number of elements in a band image unsigned long long B = Z(); //calculate the number of spectral elements @@ -1047,7 +1008,7 @@ public: stat = cublasDaxpy(handle, B, &axpy_alpha, avg_dev, 1, s_dev, 1); //subtract the average spectrum stat = cublasDsyr(handle, CUBLAS_FILL_MODE_UPPER, B, &ger_alpha, s_dev, 1, A_dev, B); //calculate the covariance matrix (symmetric outer product) } - thread_data = (double)xy / XY * 100; //record the current progress + if(PROGRESS) progress = (double)(xy+1) / XY * 100; //record the current progress } @@ -1063,7 +1024,6 @@ public: } } - thread_data = 100; //processing complete return true; } #endif @@ -1073,23 +1033,19 @@ public: /// @param co is a pointer to pre-allocated memory of size [B * B] that stores the resulting covariance matrix /// @param avg is a pointer to memory of size B that stores the average spectrum /// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location - bool co_matrix(double* co, double* avg, unsigned char *mask){ + bool co_matrix(double* co, double* avg, unsigned char *mask, bool PROGRESS = false){ #ifdef CUDA_FOUND - return cu_co_matrix(co, avg, mask); //if CUDA is available, use cuBLAS to calculate the covariance matrix + return co_matrix_cublas(co, avg, mask, PROGRESS); //if CUDA is available, use cuBLAS to calculate the covariance matrix #endif - thread_data = 0; + progress = 0; //memory allocation unsigned long long XY = X() * Y(); unsigned long long B = Z(); T* temp = (T*)malloc(sizeof(T) * B); //count vaild pixels in a band - unsigned long long count = 0; - for (unsigned long long j = 0; j < XY; j++){ - if (mask == NULL || mask[j] != 0){ - count++; - } - } + unsigned long long count = nnz(mask); + //initialize covariance matrix memset(co, 0, B * B * sizeof(double)); @@ -1110,7 +1066,7 @@ public: co_half[idx++] += temp_precise[b0] * temp_precise[b1]; } } - thread_data = (double)xy / XY * 100; + if(PROGRESS) progress = (double)(xy+1) / XY * 100; } idx = 0; for (unsigned long long i = 0; i < B; i++){ //copy the precision matrix to both halves of the output matrix @@ -1119,12 +1075,78 @@ public: } } - thread_data = 100; //processing complete free(temp); free(temp_precise); return true; } + bool project(std::string outfile, double* center, double* basis, unsigned long long M, bool PROGRESS = false){ + + std::ofstream target(outfile.c_str(), std::ios::binary); //open the target binary file + std::string headername = outfile + ".hdr"; //the header file name + + //memory allocation + unsigned long long XY = X() * Y(); + unsigned long long B = Z(); + + T* s = (T*)malloc(sizeof(T) * B); //allocate space for the spectrum + T* rs = (T*)malloc(sizeof(T) * M); //allocate space for the projected spectrum + double* bv; //pointer to the current projection vector + for(unsigned long long xy = 0; xy < XY; xy++){ //for each spectrum in the image + pixel(s, xy); //load the spectrum + memset(rs, 0, sizeof(T) * M); //initialize the rotated spectrum to zero (0) + for(unsigned long long m = 0; m < M; m++){ //for each basis vector + bv = &basis[m * B]; //assign 'bv' to the beginning of the basis vector + for(unsigned long long b = 0; b < B; b++){ //for each band + rs[m] += (s[b] - center[b]) * bv[b]; //center the spectrum and perform the projection + } + } + + target.write(reinterpret_cast(rs), sizeof(T) * M); //write the projected vector + if(PROGRESS) progress = (double)(xy+1) / XY * 100; + } + + free(s); //free temporary storage arrays + free(rs); + target.close(); //close the output file + + return true; + } + + bool inverse(std::string outfile, double* center, double* basis, unsigned long long B, unsigned long long C = 0, bool PROGRESS = false){ + + std::ofstream target(outfile.c_str(), std::ios::binary); //open the target binary file + std::string headername = outfile + ".hdr"; //the header file name + + //memory allocation + unsigned long long XY = X() * Y(); + if(C == 0) C = Z(); //if no coefficient number is given, assume all are used + C = std::min(C, Z()); //set the number of coefficients (the user can specify fewer) + + T* coeff = (T*)malloc(sizeof(T) * Z()); //allocate space for the coefficients + T* s = (T*)malloc(sizeof(T) * B); //allocate space for the spectrum + double* bv; //pointer to the current projection vector + for(unsigned long long xy = 0; xy < XY; xy++){ //for each pixel in the image (expressed as a set of coefficients) + pixel(coeff, xy); //load the coefficients + memset(s, 0, sizeof(T) * B); //initialize the spectrum to zero (0) + for(unsigned long long c = 0; c < C; c++){ //for each basis vector coefficient + bv = &basis[c * B]; //assign 'bv' to the beginning of the basis vector + for(unsigned long long b = 0; b < B; b++){ //for each component of the basis vector + s[b] += coeff[c] * bv[b] + center[b]; //calculate the contribution of each element of the basis vector in the final spectrum + } + } + + target.write(reinterpret_cast(s), sizeof(T) * B); //write the projected vector + if(PROGRESS) progress = (double)(xy+1) / XY * 100; + } + + free(coeff); //free temporary storage arrays + free(s); + target.close(); //close the output file + + return true; + } + /// Crop a region of the image and save it to a new file. @@ -1138,7 +1160,8 @@ public: unsigned long long x1, unsigned long long y1, unsigned long long b0, - unsigned long long b1){ + unsigned long long b1, + bool PROGRESS = false){ //calculate the new number of samples, lines, and bands unsigned long long samples = x1 - x0; @@ -1180,14 +1203,13 @@ public: file.seekg(jump_sample, std::ios::cur); - thread_data = (double)(y * samples + x) / (lines * samples) * 100; + if(PROGRESS) progress = (double)((y+1) * samples + x + 1) / (lines * samples) * 100; } file.seekg(jump_line, std::ios::cur); } free(temp); - thread_data = 100; return true; } diff --git a/stim/envi/bsq.h b/stim/envi/bsq.h index 5bc20c9..6c5c180 100644 --- a/stim/envi/bsq.h +++ b/stim/envi/bsq.h @@ -42,13 +42,14 @@ protected: return R[2]; } - using binary::thread_data; + using binary::nnz; + using binary::progress; public: using binary::open; using binary::file; - using binary::read_line_01; + using binary::read_line_2; using binary::read_plane_2; //using binary::getSlice; @@ -76,15 +77,15 @@ public: /// @param p is a pointer to an allocated region of memory at least X * Y * sizeof(T) in size. /// @param page <= B is the integer number of the band to be copied. bool band_index( T * p, unsigned int page){ - return read_plane_2(p, page); + return read_plane_2(p, page); //call the binary read_plane function (don't let it update the progress) } /// Retrieve a single band (by numerical label) and stores it in pre-allocated memory. /// @param p is a pointer to an allocated region of memory at least X * Y * sizeof(T) in size. /// @param wavelength is a floating point value (usually a wavelength in spectral data) used as a label for the band to be copied. - bool band( T * p, double wavelength){ - + bool band( T * p, double wavelength, bool PROGRESS = false){ + if(PROGRESS) progress = 0; //if there are no wavelengths in the BSQ file if(w.size() == 0) return band_index(p, (unsigned int)wavelength); @@ -132,7 +133,7 @@ public: else{ band_index(p, page); //return the band } - + if(PROGRESS) progress = 100; return true; } @@ -141,8 +142,8 @@ public: /// @param p is a pointer to pre-allocated memory at least B * sizeof(T) in size. /// @param x is the x-coordinate (dimension 1) of the spectrum. /// @param y is the y-coordinate (dimension 2) of the spectrum. - bool spectrum(T * p, unsigned x, unsigned y){ - return read_line_01(p, x, y); + bool spectrum(T * p, unsigned x, unsigned y, bool PROGRESS = false){ + return read_line_2(p, x, y, PROGRESS); } /// Retrieve a single pixel and stores it in pre-allocated memory. @@ -171,8 +172,9 @@ public: /// @param outname is the name of the output file used to store the resulting baseline-corrected data. /// @param wls is the list of baseline points based on band labels. - bool baseline(std::string outname, std::vector wls ) + bool baseline(std::string outname, std::vector wls, bool PROGRESS = false ) { + if(PROGRESS) progress = 0; unsigned N = wls.size(); //get the number of baseline points std::ofstream target(outname.c_str(), std::ios::binary); //open the target binary file @@ -261,7 +263,7 @@ public: target.write(reinterpret_cast(c), S); //write the corrected data into destination - thread_data = (double)cii / B * 100; + if(PROGRESS)progress = (double)(cii+1) / B * 100; } @@ -272,7 +274,6 @@ public: free(c); target.close(); - thread_data = 100; return true; } @@ -281,7 +282,7 @@ public: /// @param outname is the name of the output file used to store the resulting baseline-corrected data. /// @param w is the label specifying the band that the hyperspectral image will be normalized to. /// @param t is a threshold specified such that a spectrum with a value at w less than t is set to zero. Setting this threshold allows the user to limit division by extremely small numbers. - bool normalize(std::string outname, double w, double t = 0.0) + bool normalize(std::string outname, double w, double t = 0.0, bool PROGRESS = false) { unsigned int B = Z(); //calculate the number of bands unsigned int XY = X() * Y(); //calculate the number of pixels in a band @@ -310,7 +311,7 @@ public: } target.write(reinterpret_cast(c), S); //write normalized data into destination - thread_data = (double)j / B * 100; + if(PROGRESS) progress = (double)(j+1) / B * 100; } @@ -320,14 +321,13 @@ public: free(b); free(c); target.close(); - thread_data = 100; //make sure that the progress bar is full return true; } /// Convert the current BSQ file to a BIP file with the specified file name. /// @param outname is the name of the output BIP file to be saved to disk. - bool bip(std::string outname) + /*bool bip(std::string outname) { std::string temp = outname + "_temp"; std::string headtemp = temp + ".hdr"; @@ -340,17 +340,18 @@ public: exit(1); } //then convert bil file to bip file + progress = 0; n.bip(outname); n.close(); remove(temp.c_str()); remove(headtemp.c_str()); return true; - } + }*/ /// Convert the current BSQ file to a BIL file with the specified file name. /// @param outname is the name of the output BIL file to be saved to disk. - bool bil(std::string outname) + bool bil(std::string outname, bool PROGRESS = false) { //simplify image resolution unsigned int L = X() * Z() * sizeof(T); //calculate the number of bytes of a ZX slice @@ -370,17 +371,11 @@ public: { file.read((char *)(xz_slice + z * X()), sizeof(T) * X()); //read a line file.seekg(jump, std::ios::cur); //seek to the next band - - - thread_data = (double)(y * Z() + z) / (Z() * Y()) * 100; //update the progress counter + if(PROGRESS) progress = (double)(y * Z() + z + 1) / (Z() * Y()) * 100; //update the progress counter } target.write(reinterpret_cast(xz_slice), xz_bytes); //write the generated XZ slice to the target file } - - - thread_data = 100; - free(xz_slice); target.close(); @@ -869,13 +864,13 @@ public: continue; } - thread_data = (double)(i * XY + j) / (XY * Z()) * 100; + progress = (double)(i * XY + j) / (XY * Z()) * 100; } } target.close(); free(temp); - thread_data = 100; + progress = 100; return true; } @@ -923,7 +918,7 @@ public: i++; } //std::cout<(temp), L); //write slice data into target file file.seekg(jumpb, std::ios::cur); } free(temp); - thread_data = 100; - return true; } diff --git a/stim/envi/envi.h b/stim/envi/envi.h index 1021886..710ec51 100644 --- a/stim/envi/envi.h +++ b/stim/envi/envi.h @@ -37,27 +37,27 @@ public: if(header.interleave == envi_header::BSQ){ //if the infile is bsq file if(header.data_type ==envi_header::float32) - ((bsq*)file)->reset_thread_data(); + ((bsq*)file)->reset_progress(); else if(header.data_type == envi_header::float64) - ((bsq*)file)->reset_thread_data(); + ((bsq*)file)->reset_progress(); else std::cout<<"ERROR: unidentified data type"<*)file)->reset_thread_data(); + ((bil*)file)->reset_progress(); else if(header.data_type == envi_header::float64) - ((bil*)file)->reset_thread_data(); + ((bil*)file)->reset_progress(); else std::cout<<"ERROR: unidentified data type"<*)file)->reset_thread_data(); + ((bip*)file)->reset_progress(); else if(header.data_type == envi_header::float64) - ((bip*)file)->reset_thread_data(); + ((bip*)file)->reset_progress(); else std::cout<<"ERROR: unidentified data type"<*)file)->get_thread_data(); + return ((bsq*)file)->get_progress(); else if(header.data_type == envi_header::float64) - return ((bsq*)file)->get_thread_data(); + return ((bsq*)file)->get_progress(); else std::cout<<"ERROR: unidentified data type"<*)file)->get_thread_data(); + return ((bil*)file)->get_progress(); else if(header.data_type == envi_header::float64) - return ((bil*)file)->get_thread_data(); + return ((bil*)file)->get_progress(); else std::cout<<"ERROR: unidentified data type"<*)file)->get_thread_data(); + return ((bip*)file)->get_progress(); else if(header.data_type == envi_header::float64) - return ((bip*)file)->get_thread_data(); + return ((bip*)file)->get_progress(); else std::cout<<"ERROR: unidentified data type"< threshold (preventing division by small numbers) - bool normalize(std::string outfile, double band, double threshold = 0.0){ + bool normalize(std::string outfile, double band, double threshold = 0.0, bool PROGRESS = false){ if(header.interleave == envi_header::BSQ){ //if the infile is bsq file if(header.data_type ==envi_header::float32) - return ((bsq*)file)->normalize(outfile, band, threshold); + return ((bsq*)file)->normalize(outfile, band, threshold, PROGRESS); else if(header.data_type == envi_header::float64) - return ((bsq*)file)->normalize(outfile,band, threshold); + return ((bsq*)file)->normalize(outfile,band, threshold, PROGRESS); else std::cout<<"ERROR: unidentified data type"<*)file)->normalize(outfile, band); + return ((bil*)file)->normalize(outfile, band, threshold, PROGRESS); else if(header.data_type == envi_header::float64) - return ((bil*)file)->normalize(outfile,band); + return ((bil*)file)->normalize(outfile,band, threshold, PROGRESS); else std::cout<<"ERROR: unidentified data type"<*)file)->normalize(outfile, band); + return ((bip*)file)->normalize(outfile, band, threshold, PROGRESS); else if(header.data_type == envi_header::float64) - return ((bip*)file)->normalize(outfile,band); + return ((bip*)file)->normalize(outfile,band, threshold, PROGRESS); else std::cout<<"ERROR: unidentified data type"< w){ + bool baseline(std::string outfile, std::vector w, bool PROGRESS){ if(header.interleave == envi_header::BSQ){ //if the infile is bsq file if(header.data_type ==envi_header::float32) - return ((bsq*)file)->baseline(outfile, w); + return ((bsq*)file)->baseline(outfile, w, PROGRESS); else if(header.data_type == envi_header::float64) - return ((bsq*)file)->baseline(outfile,w); + return ((bsq*)file)->baseline(outfile,w, PROGRESS); else{ std::cout<<"ERROR: unidentified data type"<*)file)->baseline(outfile, w); + return ((bil*)file)->baseline(outfile, w, PROGRESS); else if(header.data_type == envi_header::float64) - return ((bil*)file)->baseline(outfile, w); + return ((bil*)file)->baseline(outfile, w, PROGRESS); else{ std::cout<<"ERROR: unidentified data type"<*)file)->baseline(outfile, w); + return ((bip*)file)->baseline(outfile, w, PROGRESS); else if(header.data_type == envi_header::float64) - return ((bip*)file)->baseline(outfile, w); + return ((bip*)file)->baseline(outfile, w, PROGRESS); else{ std::cout<<"ERROR: unidentified data type"<*)file)->project(outfile, center, basis, M, PROGRESS); + else if(header.data_type == envi_header::float64) + ((bip*)file)->project(outfile, center, basis, M, PROGRESS); + else{ + std::cout<<"ERROR: unidentified data type"<*)file)->inverse(outfile, center, basis, B, C, PROGRESS); + else if(header.data_type == envi_header::float64) + ((bip*)file)->inverse(outfile, center, basis, B, C, PROGRESS); + else{ + std::cout<<"ERROR: unidentified data type"<*)file)->bil(outfile); - else if(interleave == envi_header::BIP) //if the target file is bip file - return ((bsq*)file)->bip(outfile); + return ((bsq*)file)->bil(outfile, PROGRESS); + else if(interleave == envi_header::BIP){ //if the target file is bip file + std::cout<<"ERROR: conversion from BSQ to BIP isn't practical; use BSQ->BIL->BIP instead"<*)file)->bip(outfile, PROGRESS); + exit(1); + } } else if(header.data_type == envi_header::float64){ //if the data type is float @@ -319,9 +376,12 @@ public: exit(1); } else if(interleave == envi_header::BIL) - return ((bsq*)file)->bil(outfile); - else if(interleave == envi_header::BIP) - return ((bsq*)file)->bip(outfile); + return ((bsq*)file)->bil(outfile, PROGRESS); + else if(interleave == envi_header::BIP){ //if the target file is bip file + std::cout<<"ERROR: conversion from BSQ to BIP isn't practical; use BSQ->BIL->BIP instead"<*)file)->bip(outfile, PROGRESS); + exit(1); + } } else{ @@ -338,9 +398,9 @@ public: exit(1); } else if(interleave == envi_header::BSQ) //if the target file is bsq file - return ((bil*)file)->bsq(outfile); + return ((bil*)file)->bsq(outfile, PROGRESS); else if(interleave == envi_header::BIP) //if the target file is bip file - return ((bil*)file)->bip(outfile); + return ((bil*)file)->bip(outfile, PROGRESS); } else if(header.data_type == envi_header::float64){ //if the data type is float @@ -349,9 +409,9 @@ public: exit(1); } else if(interleave == envi_header::BSQ) - return ((bil*)file)->bsq(outfile); + return ((bil*)file)->bsq(outfile, PROGRESS); else if(interleave == envi_header::BIP) - return ((bil*)file)->bip(outfile); + return ((bil*)file)->bip(outfile, PROGRESS); } else{ @@ -368,9 +428,12 @@ public: exit(1); } else if(interleave == envi_header::BIL) //if the target file is bil file - return ((bip*)file)->bil(outfile); - else if(interleave == envi_header::BSQ) //if the target file is bsq file - return ((bip*)file)->bsq(outfile); + return ((bip*)file)->bil(outfile, PROGRESS); + else if(interleave == envi_header::BSQ){ //if the target file is bip file + std::cout<<"ERROR: conversion from BIP to BSQ isn't practical; use BIP->BIL->BSQ instead"<*)file)->bip(outfile, PROGRESS); + exit(1); + } } else if(header.data_type == envi_header::float64){ //if the data type is float @@ -379,9 +442,12 @@ public: exit(1); } else if(interleave == envi_header::BIL) //if the target file is bil file - return ((bip*)file)->bil(outfile); - else if(interleave == envi_header::BSQ) //if the target file is bsq file - return ((bip*)file)->bsq(outfile); + return ((bip*)file)->bil(outfile, PROGRESS); + else if(interleave == envi_header::BSQ){ //if the target file is bip file + std::cout<<"ERROR: conversion from BIP to BSQ isn't practical; use BIP->BIL->BSQ instead"<*)file)->bip(outfile, PROGRESS); + exit(1); + } } else{ @@ -872,13 +938,13 @@ public: /// @param p is a pointer to an allocated region of memory at least X * Y * sizeof(T) in size. /// @param wavelength is a floating point value (usually a wavelength in spectral data) used as a label for the band to be copied. - bool band(void* ptr, double wavelength){ + bool band(void* ptr, double wavelength, bool PROGRESS = false){ if(header.interleave == envi_header::BSQ){ //if the infile is bsq file if(header.data_type ==envi_header::float32) - return ((bsq*)file)->band((float*)ptr, wavelength); + return ((bsq*)file)->band((float*)ptr, wavelength, PROGRESS); else if (header.data_type == envi_header::float64) - return ((bsq*)file)->band((double*)ptr, wavelength); + return ((bsq*)file)->band((double*)ptr, wavelength, PROGRESS); else{ std::cout << "ERROR: unidentified data type" << std::endl; exit(1); @@ -886,9 +952,9 @@ public: } else if (header.interleave == envi_header::BIL){ if (header.data_type == envi_header::float32) - return ((bil*)file)->band((float*)ptr, wavelength); + return ((bil*)file)->band((float*)ptr, wavelength, PROGRESS); else if (header.data_type == envi_header::float64) - return ((bil*)file)->band((double*)ptr, wavelength); + return ((bil*)file)->band((double*)ptr, wavelength, PROGRESS); else{ std::cout << "ERROR: unidentified data type" << std::endl; exit(1); @@ -896,9 +962,9 @@ public: } else if (header.interleave == envi_header::BIP){ if (header.data_type == envi_header::float32) - return ((bip*)file)->band((float*)ptr, wavelength); + return ((bip*)file)->band((float*)ptr, wavelength, PROGRESS); else if (header.data_type == envi_header::float64) - return ((bip*)file)->band((double*)ptr, wavelength); + return ((bip*)file)->band((double*)ptr, wavelength, PROGRESS); else{ std::cout << "ERROR: unidentified data type" << std::endl; exit(1); @@ -912,13 +978,13 @@ public: /// @param ptr is a pointer to pre-allocated memory of size B*sizeof(T) /// @param x is the x-coordinate of the spectrum /// @param y is the y-coordinate of the spectrum - bool spectrum(void* ptr, unsigned int x, unsigned int y){ + bool spectrum(void* ptr, unsigned int x, unsigned int y, bool PROGRESS = false){ if(header.interleave == envi_header::BSQ){ //if the infile is bsq file if(header.data_type ==envi_header::float32) - return ((bsq*)file)->spectrum((float*)ptr, x, y); + return ((bsq*)file)->spectrum((float*)ptr, x, y, PROGRESS); else if (header.data_type == envi_header::float64) - return ((bsq*)file)->spectrum((double*)ptr, x, y); + return ((bsq*)file)->spectrum((double*)ptr, x, y, PROGRESS); else{ std::cout << "ERROR: unidentified data type" << std::endl; exit(1); @@ -926,9 +992,9 @@ public: } else if (header.interleave == envi_header::BIL){ if (header.data_type == envi_header::float32) - return ((bil*)file)->spectrum((float*)ptr, x, y); + return ((bil*)file)->spectrum((float*)ptr, x, y, PROGRESS); else if (header.data_type == envi_header::float64) - return ((bil*)file)->spectrum((double*)ptr, x, y); + return ((bil*)file)->spectrum((double*)ptr, x, y, PROGRESS); else{ std::cout << "ERROR: unidentified data type" << std::endl; exit(1); @@ -936,9 +1002,9 @@ public: } else if (header.interleave == envi_header::BIP){ if (header.data_type == envi_header::float32) - return ((bip*)file)->spectrum((float*)ptr, x, y); + return ((bip*)file)->spectrum((float*)ptr, x, y, PROGRESS); else if (header.data_type == envi_header::float64) - return ((bip*)file)->spectrum((double*)ptr, x, y); + return ((bip*)file)->spectrum((double*)ptr, x, y, PROGRESS); else{ std::cout << "ERROR: unidentified data type" << std::endl; exit(1); @@ -989,12 +1055,12 @@ public: /// @param p is a pointer to pre-allocated memory of size [B * sizeof(T)] that stores the mean spectrum /// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location - bool avg_band(double * p, unsigned char* mask){ + bool avg_band(double * p, unsigned char* mask, bool PROGRESS = false){ if (header.interleave == envi_header::BSQ){ if (header.data_type == envi_header::float32) - return ((bsq*)file)->avg_band(p, mask); + return ((bsq*)file)->avg_band(p, mask, PROGRESS); else if (header.data_type == envi_header::float64) - return ((bsq*)file)->avg_band(p, mask); + return ((bsq*)file)->avg_band(p, mask, PROGRESS); else{ std::cout << "ERROR: unidentified data type" << std::endl; exit(1); @@ -1002,9 +1068,9 @@ public: } else if (header.interleave == envi_header::BIL){ if (header.data_type == envi_header::float32) - return ((bil*)file)->avg_band(p, mask); + return ((bil*)file)->avg_band(p, mask, PROGRESS); else if (header.data_type == envi_header::float64) - return ((bil*)file)->avg_band(p, mask); + return ((bil*)file)->avg_band(p, mask, PROGRESS); else{ std::cout << "ERROR: unidentified data type" << std::endl; exit(1); @@ -1012,9 +1078,9 @@ public: } else if (header.interleave == envi_header::BIP){ if (header.data_type == envi_header::float32) - return ((bip*)file)->avg_band(p, mask); + return ((bip*)file)->avg_band(p, mask, PROGRESS); else if (header.data_type == envi_header::float64) - return ((bip*)file)->avg_band(p, mask); + return ((bip*)file)->avg_band(p, mask, PROGRESS); else{ std::cout << "ERROR: unidentified data type" << std::endl; exit(1); @@ -1028,22 +1094,24 @@ public: /// @param co is a pointer to pre-allocated memory of size [B * B] that stores the resulting covariance matrix /// @param avg is a pointer to memory of size B that stores the average spectrum /// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location - bool co_matrix(double* co, double* avg, unsigned char* mask){ + bool co_matrix(double* co, double* avg, unsigned char* mask, bool PROGRESS = false){ if (header.interleave == envi_header::BSQ){ - if (header.data_type == envi_header::float32) - return ((bsq*)file)->co_matrix(co, avg, mask); + std::cout<<"ERROR: calculating the covariance matrix for a BSQ file is impractical; convert to BIL or BIP first"<*)file)->co_matrix(co, avg, mask, PROGRESS); else if (header.data_type == envi_header::float64) - return ((bsq*)file)->co_matrix(co, avg, mask); + return ((bsq*)file)->co_matrix(co, avg, mask, PROGRESS); else{ std::cout << "ERROR: unidentified data type" << std::endl; exit(1); - } + }*/ } else if (header.interleave == envi_header::BIL){ if (header.data_type == envi_header::float32) - return ((bil*)file)->co_matrix(co, avg, mask); + return ((bil*)file)->co_matrix(co, avg, mask, PROGRESS); else if (header.data_type == envi_header::float64) - return ((bil*)file)->co_matrix(co, avg, mask); + return ((bil*)file)->co_matrix(co, avg, mask, PROGRESS); else{ std::cout << "ERROR: unidentified data type" << std::endl; exit(1); @@ -1051,9 +1119,9 @@ public: } else if (header.interleave == envi_header::BIP){ if (header.data_type == envi_header::float32) - return ((bip*)file)->co_matrix(co, avg, mask); + return ((bip*)file)->co_matrix(co, avg, mask, PROGRESS); else if (header.data_type == envi_header::float64) - return ((bip*)file)->co_matrix(co, avg, mask); + return ((bip*)file)->co_matrix(co, avg, mask, PROGRESS); else{ std::cout << "ERROR: unidentified data type" << std::endl; exit(1); @@ -1070,7 +1138,7 @@ public: /// @param y0 is the lower-left y pixel coordinate to be included in the cropped image /// @param x1 is the upper-right x pixel coordinate to be included in the cropped image /// @param y1 is the upper-right y pixel coordinate to be included in the cropped image - bool crop(std::string outfile,unsigned x0, unsigned y0, unsigned x1, unsigned y1, unsigned b0, unsigned b1){ + bool crop(std::string outfile,unsigned x0, unsigned y0, unsigned x1, unsigned y1, unsigned b0, unsigned b1, bool PROGRESS = false){ //save the header for the cropped file stim::envi_header new_header = header; @@ -1084,9 +1152,9 @@ public: if (header.interleave == envi_header::BSQ){ if (header.data_type == envi_header::float32) - return ((bsq*)file)->crop(outfile, x0, y0, x1, y1, b0, b1); + return ((bsq*)file)->crop(outfile, x0, y0, x1, y1, b0, b1, PROGRESS); else if (header.data_type == envi_header::float64) - return ((bsq*)file)->crop(outfile, x0, y0, x1, y1, b0, b1); + return ((bsq*)file)->crop(outfile, x0, y0, x1, y1, b0, b1, PROGRESS); else{ std::cout << "ERROR: unidentified data type" << std::endl; exit(1); @@ -1094,9 +1162,9 @@ public: } else if (header.interleave == envi_header::BIL){ if (header.data_type == envi_header::float32) - return ((bil*)file)->crop(outfile, x0, y0, x1, y1, b0, b1); + return ((bil*)file)->crop(outfile, x0, y0, x1, y1, b0, b1, PROGRESS); else if (header.data_type == envi_header::float64) - return ((bil*)file)->crop(outfile, x0, y0, x1, y1, b0, b1); + return ((bil*)file)->crop(outfile, x0, y0, x1, y1, b0, b1, PROGRESS); else{ std::cout << "ERROR: unidentified data type" << std::endl; exit(1); @@ -1104,9 +1172,9 @@ public: } else if (header.interleave == envi_header::BIP){ if (header.data_type == envi_header::float32) - return ((bip*)file)->crop(outfile, x0, y0, x1, y1, b0, b1); + return ((bip*)file)->crop(outfile, x0, y0, x1, y1, b0, b1, PROGRESS); else if (header.data_type == envi_header::float64) - return ((bip*)file)->crop(outfile, x0, y0, x1, y1, b0, b1); + return ((bip*)file)->crop(outfile, x0, y0, x1, y1, b0, b1, PROGRESS); else{ std::cout << "ERROR: unidentified data type" << std::endl; exit(1); -- libgit2 0.21.4