diff --git a/stim/envi/bil.h b/stim/envi/bil.h index c571d4e..6890885 100644 --- a/stim/envi/bil.h +++ b/stim/envi/bil.h @@ -3,8 +3,10 @@ #include "../envi/envi_header.h" #include "../envi/hsi.h" +#include "../math/fd_coefficients.h" #include #include +#include namespace stim{ @@ -22,26 +24,14 @@ template class bil: public hsi { protected: - - //std::vector w; //band wavelengths - /*unsigned long long X(){ - return R[0]; - } - unsigned long long Y(){ - return R[2]; - } - unsigned long long Z(){ - return R[1]; - }*/ + using hsi::w; //use the wavelength array in stim::hsi using hsi::nnz; 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()); - //} + using hsi::X; + using hsi::Y; + using hsi::Z; public: @@ -59,11 +49,11 @@ public: /// @param B is the number of samples (bands) along dimension 3 /// @param header_offset is the number of bytes (if any) in the binary header /// @param wavelengths is an optional STL vector of size B specifying a numerical label for each band - bool open(std::string filename, - unsigned long long X, - unsigned long long Y, - unsigned long long B, - unsigned long long header_offset, + bool open(std::string filename, + unsigned long long X, + unsigned long long Y, + unsigned long long B, + unsigned long long header_offset, std::vector wavelengths){ w = wavelengths; @@ -113,7 +103,7 @@ public: unsigned long long S = XY * sizeof(T); //calculate the number of bytes of a band unsigned long long page=0; //bands around the wavelength - + //get the bands numbers around the wavelength @@ -173,7 +163,7 @@ public: //if wavelength is smaller than the first one in header file if ( w[page] > wavelength ){ - memcpy(p, c, L); + memcpy(p, c, L); return true; } @@ -193,16 +183,16 @@ public: memcpy(p1, c + (page - 1) * X(), L); memcpy(p2, c + page * X(), L); - + for(unsigned long long i=0; i < X(); i++){ double r = (double) (wavelength - w[page-1]) / (double) (w[page] - w[page-1]); p[i] = (T)(((double)p2[i] - (double)p1[i]) * r + (double)p1[i]); } } - else //if the wavelength is equal to a wavelength in header file + else //if the wavelength is equal to a wavelength in header file memcpy(p, c + page * X(), L); - - return true; + + return true; } /// Retrieve a single spectrum (B-axis line) at a given (x, y) location and stores it in pre-allocated memory. @@ -227,24 +217,24 @@ public: //get the pixel return spectrum(p, x, y); } - + //given a Y ,return a XZ slice bool read_plane_y(T * p, unsigned long long y){ return binary::read_plane_2(p, y); } - + /// Perform baseline correction given a list of baseline points and stores the result in a new BSQ file. /// @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, unsigned char* mask = NULL, bool PROGRESS = false){ - + unsigned long long N = wls.size(); //get the number of baseline points - + std::ofstream target(outname.c_str(), std::ios::binary); //open the target binary file - std::string headername = outname + ".hdr"; //the header file name - + std::string headername = outname + ".hdr"; //the header file name + //simplify image resolution unsigned long long ZX = Z() * X(); //calculate the number of points in a Y slice unsigned long long L = ZX * sizeof(T); //calculate the number of bytes of a Y slice @@ -252,10 +242,10 @@ public: 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)); @@ -269,19 +259,19 @@ public: exit(1); } // loop start correct every y slice - + for (unsigned long long k =0; k < Y(); k++) { //get the current y slice read_plane_y(c, k); - - //initialize lownum, highnum, low, high + + //initialize lownum, highnum, low, high ai = w[0]; control=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]; + bi = wls[control]; memset(a, (char)0, X() * sizeof(T) ); } //else get the low band @@ -292,39 +282,39 @@ public: } //get the high band read_x_from_xz(b, c, bi); - + //correct every YZ line - + for(unsigned long long 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]; } } - + ci = w[cii]; - + unsigned long long jump = cii * X(); //perform the baseline correction for(unsigned i=0; i < X(); i++) @@ -332,22 +322,22 @@ public: double r = (double) (ci - ai) / (double) (bi - ai); c[i + jump] =(T) ( c[i + jump] - (b[i] - a[i]) * r - a[i] ); } - - }//loop for YZ line end - target.write(reinterpret_cast(c), L); //write the corrected data into destination + + }//loop for YZ line end + target.write(reinterpret_cast(c), L); //write the corrected data into destination if(PROGRESS) progress = (double)(k+1) / Y() * 100; - }//loop for Y slice end - + }//loop for Y slice end + free(a); free(b); free(c); target.close(); - return true; - + return true; + } - + /// Normalize all spectra based on the value of a single band, storing the result in a new BSQ file. /// @param outname is the name of the output file used to store the resulting baseline-corrected data. @@ -358,7 +348,7 @@ public: unsigned long long B = Z(); //calculate the number of bands unsigned long long ZX = Z() * X(); unsigned long long XY = X() * Y(); //calculate the number of pixels in a band - unsigned long long S = XY * sizeof(T); //calculate the number of bytes in a band + unsigned long long S = XY * sizeof(T); //calculate the number of bytes in a band unsigned long long L = ZX * sizeof(T); std::ofstream target(outname.c_str(), std::ios::binary); //open the target binary file @@ -368,7 +358,7 @@ public: T * b; //pointer to the standard band b = (T*)malloc( S ); //memory allocation - c = (T*)malloc( L ); + c = (T*)malloc( L ); band(b, w); //get the certain band into memory @@ -383,13 +373,13 @@ public: c[m + i * X()] = (T)0.0; else c[m + i * X()] = c[m + i * X()] / b[m + j * X()]; - } + } } target.write(reinterpret_cast(c), L); //write normalized data into destination if(PROGRESS) progress = (double)(j+1) / Y() * 100; } - + free(b); free(c); target.close(); @@ -397,26 +387,26 @@ public: } void normalize(std::string outfile, unsigned char* mask = NULL, bool PROGRESS = false){} - + /// 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 PROGRESS = false) { unsigned long long S = X() * Y() * sizeof(T); //calculate the number of bytes in a band - + std::ofstream target(outname.c_str(), std::ios::binary); std::string headername = outname + ".hdr"; - + T * p; //pointer to the current band p = (T*)malloc(S); - + for ( unsigned long long i = 0; i < Z(); i++) - { + { band_index(p, i); target.write(reinterpret_cast(p), S); //write a band data into target file - if(PROGRESS) progress = (double)(i+1) / Z() * 100; //store the progress for the current operation + if(PROGRESS) progress = (double)(i+1) / Z() * 100; //store the progress for the current operation } free(p); @@ -430,17 +420,17 @@ public: bool bip(std::string outname, bool PROGRESS = false) { unsigned long long S = X() * Z() * sizeof(T); //calculate the number of bytes in a ZX slice - + std::ofstream target(outname.c_str(), std::ios::binary); std::string headername = outname + ".hdr"; - + T * p; //pointer to the current XZ slice for bil file p = (T*)malloc(S); T * q; //pointer to the current ZX slice for bip file q = (T*)malloc(S); - + for ( unsigned long long i = 0; i < Y(); i++) - { + { read_plane_y(p, i); for ( unsigned long long k = 0; k < Z(); k++) { @@ -448,10 +438,10 @@ public: for ( unsigned long long j = 0; j < X(); j++) q[k + j * Z()] = p[ks + j]; - if(PROGRESS) progress = (double)((i+1) * Z() + k+1) / (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 + + target.write(reinterpret_cast(q), S); //write a band data into target file } free(p); @@ -498,7 +488,7 @@ public: rp = (T*) malloc(S); band(lp, lb); - band(rp, rb); + band(rp, rb); baseline_band(lb, rb, lp, rp, bandwavelength, result); @@ -570,7 +560,7 @@ public: } baseline_band(lb, rb, lp, rp, lab, cur2); //beginnning part baseline_band(lb, rb, lp, rp, w[ai], cur); - for(unsigned long long j = 0; j < XY; j++){ + for(unsigned long long j = 0; j < XY; j++){ result[j] += (T)((w[ai] - lab) * ((double)cur[j] + (double)cur2[j]) / 2.0); } @@ -581,7 +571,7 @@ public: baseline_band(lb, rb, lp, rp, w[ai], cur2); for(unsigned long long j = 0; j < XY; j++) { - result[j] += (T)((w[ai] - w[ai-1]) * ((double)cur[j] + (double)cur2[j]) / 2.0); + result[j] += (T)((w[ai] - w[ai-1]) * ((double)cur[j] + (double)cur2[j]) / 2.0); } std::swap(cur,cur2); //swap the band pointers } @@ -606,7 +596,7 @@ public: T* p1 = (T*)malloc(X() * Y() * sizeof(T)); T* p2 = (T*)malloc(X() * Y() * sizeof(T)); - + //get the two peak band height(lb1, rb1, pos1, p1); height(lb2, rb2, pos2, p2); @@ -620,9 +610,9 @@ public: free(p1); free(p2); - return true; + return true; } - + /// Compute the ratio between a peak area and peak height. /// @param lb1 is the label value for the left baseline point for the first peak (numerator) @@ -634,10 +624,10 @@ public: /// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size bool pa_to_ph(T* result, double lb1, double rb1, double lab1, double rab1, double lb2, double rb2, double pos, unsigned char* mask = NULL){ - + T* p1 = (T*)malloc(X() * Y() * sizeof(T)); T* p2 = (T*)malloc(X() * Y() * sizeof(T)); - + //get the area and the peak band area(lb1, rb1, lab1, rab1, p1); height(lb2, rb2, pos, p2); @@ -651,9 +641,9 @@ public: free(p1); free(p2); - return true; - } - + return true; + } + /// Compute the ratio between two peak areas. /// @param lb1 is the label value for the left baseline point for the first peak (numerator) @@ -663,14 +653,14 @@ public: /// @param lb2 is the label value for the left baseline point for the second peak (denominator) /// @param rb2 is the label value for the right baseline point for the second peak (denominator) /// @param lab2 is the label value for the left bound (start of the integration) of the second peak (denominator) - /// @param rab2 is the label value for the right bound (end of the integration) of the second peak (denominator) + /// @param rab2 is the label value for the right bound (end of the integration) of the second peak (denominator) /// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size bool pa_to_pa(T* result, double lb1, double rb1, double lab1, double rab1, double lb2, double rb2, double lab2, double rab2, unsigned char* mask = NULL){ - + T* p1 = (T*)malloc(X() * Y() * sizeof(T)); T* p2 = (T*)malloc(X() * Y() * sizeof(T)); - + //get the area and the peak band area(lb1, rb1, lab1, rab1, p1); area(lb2, rb2, lab2, rab2, p2); @@ -684,8 +674,8 @@ public: free(p1); free(p2); - return true; - } + return true; + } /// Compute the definite integral of a baseline corrected peak. @@ -745,7 +735,7 @@ public: } baseline_band(lb, rb, lp, rp, lab, cur2); //beginnning part baseline_band(lb, rb, lp, rp, w[ai], cur); - for(unsigned long long j = 0; j < XY; j++){ + for(unsigned long long j = 0; j < XY; j++){ result[j] += (T)((w[ai] - lab) * (w[ai] + lab) * ((double)cur[j] + (double)cur2[j]) / 4.0); } @@ -756,7 +746,7 @@ public: baseline_band(lb, rb, lp, rp, w[ai], cur2); for(unsigned long long j = 0; j < XY; j++) { - result[j] += (T)((w[ai] - w[ai-1]) * (w[ai] + w[ai-1]) * ((double)cur[j] + (double)cur2[j]) / 4.0); + result[j] += (T)((w[ai] - w[ai-1]) * (w[ai] + w[ai-1]) * ((double)cur[j] + (double)cur2[j]) / 4.0); } std::swap(cur,cur2); //swap the band pointers } @@ -778,7 +768,7 @@ public: bool centroid(T* result, double lb, double rb, double lab, double rab, unsigned char* mask = NULL){ T* p1 = (T*)malloc(X() * Y() * sizeof(T)); T* p2 = (T*)malloc(X() * Y() * sizeof(T)); - + //get the area and the peak band x_area(lb, rb, lab, rab, p1); area(lb, rb, lab, rab, p2); @@ -790,7 +780,7 @@ public: free(p1); free(p2); - return true; + return true; } /// Create a mask based on a given band and threshold value. @@ -843,7 +833,7 @@ public: continue; } } - target.write(reinterpret_cast(temp), L); //write a band data into target file + target.write(reinterpret_cast(temp), L); //write a band data into target file if(PROGRESS) progress = (double)(i+1) / (double)Y() * 100; } target.close(); @@ -861,7 +851,7 @@ public: T* line = (T*) malloc( Lbytes ); //allocate space for a line file.seekg(0, std::ios::beg); //seek to the beginning of the file - + size_t pl; size_t p = 0; //create counter variables for(size_t y = 0; y < Y(); y++){ //for each line in the data set @@ -903,7 +893,7 @@ public: { read_plane_y(slice, y); //retrieve an ZX page, store in "slice" - //for each sample along X + //for each sample along X for (unsigned long long x = 0; x < X(); x++) //Select a pixel by choosing X coordinate in the page, X() { //if the mask != 0 at that xy pixel @@ -1082,7 +1072,7 @@ public: for (unsigned long long z = b0; z < b1; z++) { file.read((char *)(temp + z * samples), sizeof(T) * samples); - file.seekg(jumpb, std::ios::cur); //go to the next band + file.seekg(jumpb, std::ios::cur); //go to the next band if(PROGRESS) progress = (double)(x * Z() + z+1) / (lines * Z()) * 100; } @@ -1193,7 +1183,7 @@ public: size_t N = C.size(); //get the number of bands that the kernel spans std::deque frame(N, NULL); //create an array to store pointers to each frame for(size_t f = 0; f < N; f++) - frame[f] = (T*)malloc(Xb); //allocate space for the frame + frame[f] = (T*)malloc(Xb); //allocate space for the frame T* outline = (T*)malloc(Xb); //allocate space for the output band @@ -1219,23 +1209,23 @@ public: frame.pop_front(); //pop the first element } if(PROGRESS) progress = (double)(y+1) / (double)Y() * 100; - } + } } /// Approximate the spectral derivative of the image void deriv(std::string outfile, size_t d, size_t order, const std::vector w = std::vector(), unsigned char* mask = NULL, bool PROGRESS = false){ std::ofstream out(outfile.c_str(), std::ios::binary); //open the output file for writing - + size_t Xb = X() * sizeof(T); //calculate the size of a line (frame) in bytes size_t XBb = Xb * Z(); size_t B = Z(); - file.seekg(0, std::ios::beg); //seek to the beginning of the file + file.seekg(0, std::ios::beg); //seek to the beginning of the file size_t N = order + d; //approximating a derivative requires order + d samples std::deque frame(N, NULL); //create an array to store pointers to each frame for(size_t f = 0; f < N; f++) //for each frame - frame[f] = (T*)malloc(Xb); //allocate space for the frame + frame[f] = (T*)malloc(Xb); //allocate space for the frame T* outline = (T*)malloc(Xb); //allocate space for the output band @@ -1273,7 +1263,7 @@ public: } } out.write((char*)outline, Xb); //output the band - + } if(PROGRESS) progress = (double)(y+1) / (double)Y() * 100; } diff --git a/stim/envi/bip.h b/stim/envi/bip.h index ed6542a..908caa4 100644 --- a/stim/envi/bip.h +++ b/stim/envi/bip.h @@ -28,15 +28,18 @@ template class bip: public hsi { protected: - - + + //std::vector w; //band wavelength unsigned long long offset; //header offset using hsi::w; //use the wavelength array in stim::hsi using hsi::nnz; using binary::progress; - + using hsi::X; + using hsi::Y; + using hsi::Z; + public: using binary::open; @@ -54,11 +57,11 @@ public: /// @param B is the number of samples (bands) along dimension 3 /// @param header_offset is the number of bytes (if any) in the binary header /// @param wavelengths is an optional STL vector of size B specifying a numerical label for each band - bool open(std::string filename, - unsigned long long X, - unsigned long long Y, - unsigned long long B, - unsigned long long header_offset, + bool open(std::string filename, + unsigned long long X, + unsigned long long Y, + unsigned long long B, + unsigned long long header_offset, std::vector wavelengths){ //copy the wavelengths to the BSQ file structure @@ -67,7 +70,7 @@ public: offset = header_offset; return open(filename, vec(B, X, Y), header_offset); - + } /// Retrieve a single band (based on index) and stores it in pre-allocated memory. @@ -160,7 +163,7 @@ public: for(unsigned long long j = 0; j < X(); j++) { p[j] = c[j * B]; - } + } return true; } @@ -192,7 +195,7 @@ public: { p2[j] = c[j * B + page]; } - + for(unsigned long long i=0; i < X(); i++){ double r = (double) (wavelength - w[page-1]) / (double) (w[page] - w[page-1]); p[i] = (p2[i] - p1[i]) * r + p1[i]; @@ -209,7 +212,7 @@ public: } } - return true; + return true; } /// Retrieve a single pixel and store it in a pre-allocated double array. @@ -246,20 +249,20 @@ public: file.read((char *)p, sizeof(T) * Z()); return true; } - + //given a Y ,return a ZX slice bool read_plane_y(T * p, unsigned long long y){ return binary::read_plane_2(p, y); } - + /// Perform baseline correction given a list of baseline points and stores the result in a new BSQ file. /// @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. + /// @param wls is the list of baseline points based on band labels. bool baseline(std::string outname, std::vector base_pts, unsigned char* mask = NULL, bool PROGRESS = false){ - + std::ofstream target(outname.c_str(), std::ios::binary); //open the target binary file - + 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 @@ -274,9 +277,9 @@ public: memset(sbc, 0, sizeof(T) * B); //set the baseline corrected spectrum to zero } else{ - + pixel(s, n); //retrieve the spectrum s - base_vals = interp_spectrum(s, base_pts); //get the values at each baseline point + base_vals = hsi::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) @@ -298,23 +301,23 @@ public: 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 + sbc[b] = s[b] - hsi::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(s); free(sbc); target.close(); - - return true; - + + return true; + } - + /// Normalize all spectra based on the value of a single band, storing the result in a new BSQ file. /// @param outname is the name of the output file used to store the resulting baseline-corrected data. @@ -323,8 +326,8 @@ public: bool ratio(std::string outname, double w, unsigned char* mask = NULL, 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 - + 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 @@ -334,24 +337,24 @@ public: memset(s, 0, sizeof(T) * B); //set the output to zero else{ pixel(s, n); //retrieve the spectrum s - nv = interp_spectrum(s, w); //find the value of the normalization band - + nv = hsi::interp_spectrum(s, w); //find the value of the normalization band + for(unsigned long long b = 0; b < B; b++) //for each band in the spectrum s[b] /= nv; //divide by the normalization value } - + 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(s); target.close(); return true; } void normalize(std::string outfile, unsigned char* mask = NULL, bool PROGRESS = false){} - + /// Convert the current BIP file to a BIL file with the specified file name. @@ -359,17 +362,17 @@ public: bool bil(std::string outname, bool PROGRESS = false) { unsigned long long S = X() * Z() * sizeof(T); //calculate the number of bytes in a ZX slice - + std::ofstream target(outname.c_str(), std::ios::binary); //std::string headername = outname + ".hdr"; - + T * p; //pointer to the current ZX slice for bip file p = (T*)malloc(S); T * q; //pointer to the current XZ slice for bil file q = (T*)malloc(S); - + for ( unsigned long long i = 0; i < Y(); i++) - { + { read_plane_y(p, i); for ( unsigned long long k = 0; k < Z(); k++) { @@ -379,7 +382,7 @@ public: if(PROGRESS) progress = (double)(i * Z() + k+1) / (Y() * Z()) * 100; } - target.write(reinterpret_cast(q), S); //write a band data into target file + target.write(reinterpret_cast(q), S); //write a band data into target file } free(p); @@ -408,7 +411,7 @@ public: } return true; } - + /// Return a baseline corrected band given two adjacent baseline points. The result is stored in a pre-allocated array. /// @param lb is the label value for the left baseline point @@ -425,7 +428,7 @@ public: rp = (T*) malloc(S); band(lp, lb); - band(rp, rb); + band(rp, rb); baseline_band(lb, rb, lp, rp, bandwavelength, result); @@ -496,7 +499,7 @@ public: } baseline_band(lb, rb, lp, rp, lab, cur2); //beginnning part baseline_band(lb, rb, lp, rp, w[ai], cur); - for(unsigned long long j = 0; j < XY; j++){ + for(unsigned long long j = 0; j < XY; j++){ result[j] += (T)((w[ai] - lab) * ((double)cur[j] + (double)cur2[j]) / 2.0); } @@ -507,7 +510,7 @@ public: baseline_band(lb, rb, lp, rp, w[ai], cur2); for(unsigned long long j = 0; j < XY; j++) { - result[j] += (T)((w[ai] - w[ai-1]) * ((double)cur[j] + (double)cur2[j]) / 2.0); + result[j] += (T)((w[ai] - w[ai-1]) * ((double)cur[j] + (double)cur2[j]) / 2.0); } std::swap(cur,cur2); //swap the band pointers } @@ -530,9 +533,9 @@ public: /// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size bool ph_to_ph(T* result, double lb1, double rb1, double pos1, double lb2, double rb2, double pos2, unsigned char* mask = NULL){ - T* p1 = (T*)malloc(X() * Y() * sizeof(T)); - T* p2 = (T*)malloc(X() * Y() * sizeof(T)); - + T* p1 = (T*)malloc(X() * Y() * sizeof(T)); + T* p2 = (T*)malloc(X() * Y() * sizeof(T)); + //get the two peak band height(lb1, rb1, pos1, p1); height(lb2, rb2, pos2, p2); @@ -546,9 +549,9 @@ public: free(p1); free(p2); - return true; + return true; } - + /// Compute the ratio between a peak area and peak height. /// @param lb1 is the label value for the left baseline point for the first peak (numerator) @@ -560,10 +563,10 @@ public: /// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size bool pa_to_ph(T* result, double lb1, double rb1, double lab1, double rab1, double lb2, double rb2, double pos, unsigned char* mask = NULL){ - - T* p1 = (T*)malloc(X() * Y() * sizeof(T)); - T* p2 = (T*)malloc(X() * Y() * sizeof(T)); - + + T* p1 = (T*)malloc(X() * Y() * sizeof(T)); + T* p2 = (T*)malloc(X() * Y() * sizeof(T)); + //get the area and the peak band area(lb1, rb1, lab1, rab1, p1); height(lb2, rb2, pos, p2); @@ -577,9 +580,9 @@ public: free(p1); free(p2); - return true; - } - + return true; + } + /// Compute the ratio between two peak areas. /// @param lb1 is the label value for the left baseline point for the first peak (numerator) @@ -589,14 +592,14 @@ public: /// @param lb2 is the label value for the left baseline point for the second peak (denominator) /// @param rb2 is the label value for the right baseline point for the second peak (denominator) /// @param lab2 is the label value for the left bound (start of the integration) of the second peak (denominator) - /// @param rab2 is the label value for the right bound (end of the integration) of the second peak (denominator) + /// @param rab2 is the label value for the right bound (end of the integration) of the second peak (denominator) /// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size bool pa_to_pa(T* result, double lb1, double rb1, double lab1, double rab1, double lb2, double rb2, double lab2, double rab2, unsigned char* mask = NULL){ - - T* p1 = (T*)malloc(X() * Y() * sizeof(T)); - T* p2 = (T*)malloc(X() * Y() * sizeof(T)); - + + T* p1 = (T*)malloc(X() * Y() * sizeof(T)); + T* p2 = (T*)malloc(X() * Y() * sizeof(T)); + //get the area and the peak band area(lb1, rb1, lab1, rab1, p1); area(lb2, rb2, lab2, rab2, p2); @@ -610,8 +613,8 @@ public: free(p1); free(p2); - return true; - } + return true; + } /// Compute the definite integral of a baseline corrected peak. @@ -671,7 +674,7 @@ public: } baseline_band(lb, rb, lp, rp, lab, cur2); //beginnning part baseline_band(lb, rb, lp, rp, w[ai], cur); - for(unsigned long long j = 0; j < XY; j++){ + for(unsigned long long j = 0; j < XY; j++){ result[j] += (T)((w[ai] - lab) * (w[ai] + lab) * ((double)cur[j] + (double)cur2[j]) / 4.0); } @@ -682,7 +685,7 @@ public: baseline_band(lb, rb, lp, rp, w[ai], cur2); for(unsigned long long j = 0; j < XY; j++) { - result[j] += (T)((w[ai] - w[ai-1]) * (w[ai] + w[ai-1]) * ((double)cur[j] + (double)cur2[j]) / 4.0); + result[j] += (T)((w[ai] - w[ai-1]) * (w[ai] + w[ai-1]) * ((double)cur[j] + (double)cur2[j]) / 4.0); } std::swap(cur,cur2); //swap the band pointers } @@ -702,9 +705,9 @@ public: /// @param rab is the label for the end of the peak /// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size bool centroid(T* result, double lb, double rb, double lab, double rab, unsigned char* mask = NULL){ - T* p1 = (T*)malloc(X() * Y() * sizeof(T)); - T* p2 = (T*)malloc(X() * Y() * sizeof(T)); - + T* p1 = (T*)malloc(X() * Y() * sizeof(T)); + T* p2 = (T*)malloc(X() * Y() * sizeof(T)); + //get the area and the peak band x_area(lb, rb, lab, rab, p1); area(lb, rb, lab, rab, p2); @@ -716,7 +719,7 @@ public: free(p1); free(p2); - return true; + return true; } /// Create a mask based on a given band and threshold value. @@ -788,7 +791,7 @@ public: T* band = (T*) malloc( Bbytes ); //allocate space for a line file.seekg(0, std::ios::beg); //seek to the beginning of the file - + size_t p = 0; //create counter variables for(size_t xy = 0; xy < XY; xy++){ //for each pixel if(mask[xy]){ //if the current pixel is masked @@ -800,7 +803,7 @@ public: p++; //increment the pixel pointer } else - file.seekg(Bbytes, std::ios::cur); //otherwise skip this band + file.seekg(Bbytes, std::ios::cur); //otherwise skip this band } return true; } @@ -809,7 +812,7 @@ public: bool sift(std::string outfile, unsigned char* mask, bool PROGRESS = false){ //reset the file pointer to the beginning of the file - file.seekg(0, std::ios::beg); + file.seekg(0, std::ios::beg); // open an output stream std::ofstream target(outfile.c_str(), std::ios::binary); @@ -862,7 +865,7 @@ public: std::ofstream target(outfile.c_str(), std::ios::binary); //reset the file pointer to the beginning of the file - file.seekg(0, std::ios::beg); + file.seekg(0, std::ios::beg); //allocate space for a single spectrum unsigned long long B = Z(); @@ -942,10 +945,10 @@ public: /// Note that cuBLAS only supports integer-sized arrays, so there may be issues with large spectra bool co_matrix_cublas(double* co, double* avg, unsigned char *mask, bool PROGRESS = false){ - cudaError_t cudaStat; + cudaError_t cudaStat; cublasStatus_t stat; cublasHandle_t handle; - + 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 @@ -959,7 +962,7 @@ public: cudaStat = cudaMemset(A_dev, 0, B * B * sizeof(double)); //initialize the covariance matrix to zero (0) cudaStat = cudaMalloc(&avg_dev, B * sizeof(double)); //allocate space on the CUDA device for the average spectrum stat = cublasSetVector((int)B, sizeof(double), avg, 1, avg_dev, 1); //copy the average spectrum to the CUDA device - + double ger_alpha = 1.0/(double)XY; //scale the outer product by the inverse of the number of samples (mean outer product) double axpy_alpha = -1; //multiplication factor for the average spectrum (in order to perform a subtraction) @@ -994,7 +997,7 @@ public: return true; } #endif - + /// Calculate the covariance matrix for all masked pixels in the image with 64-bit floating point precision. /// @param co is a pointer to pre-allocated memory of size [B * B] that stores the resulting covariance matrix @@ -1074,7 +1077,7 @@ public: } } } - + target.write(reinterpret_cast(rs), sizeof(T) * M); //write the projected vector if(PROGRESS) progress = (double)(xy+1) / XY * 100; } @@ -1104,7 +1107,7 @@ public: 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 + for(unsigned long long b = 0; b < B; b++){ //for each component of the basis vector s[b] += (T)((double)coeff[c] * bv[b] + center[b]); //calculate the contribution of each element of the basis vector in the final spectrum } } @@ -1172,7 +1175,7 @@ public: //read the cropped spectral region file.read( (char*) temp, L ); //pixel(temp, sp + x + y * X()); - out.write(reinterpret_cast(temp), L); //write slice data into target file + out.write(reinterpret_cast(temp), L); //write slice data into target file file.seekg(jump_sample, std::ios::cur); @@ -1253,7 +1256,7 @@ public: memset(spec, 0, spec_bytes); out.write((char*)spec, spec_bytes); //write to the output file } - if(PROGRESS) progress = (double)( (y+1) * S[0] + 1) / (double) (S[0] * S[1]) * 100; + if(PROGRESS) progress = (double)( (y+1) * S[0] + 1) / (double) (S[0] * S[1]) * 100; } } @@ -1267,7 +1270,7 @@ public: void convolve(std::string outfile, std::vector C, size_t start, size_t end, unsigned char* mask = NULL, bool PROGRESS = false){ std::ofstream out(outfile.c_str(), std::ios::binary); //open the output file for writing - + size_t N = end - start + 1; //number of bands in the output spectrum size_t Nb = N * sizeof(T); //size of the output spectrum in bytes size_t B = Z(); //calculate the number of values in a spectrum @@ -1297,8 +1300,8 @@ public: void deriv(std::string outfile, size_t d, size_t order, const std::vector w, unsigned char* mask = NULL, bool PROGRESS = false){ std::ofstream out(outfile.c_str(), std::ios::binary); //open the output file for writing - - + + size_t B = Z(); //calculate the number of values in a spectrum size_t Bb = B * sizeof(T); //calculate the size of a spectrum in bytes diff --git a/stim/envi/bsq.h b/stim/envi/bsq.h index b523dab..b383727 100644 --- a/stim/envi/bsq.h +++ b/stim/envi/bsq.h @@ -37,6 +37,9 @@ protected: using hsi::w; //use the wavelength array in stim::hsi using hsi::nnz; using binary::progress; + using hsi::X; + using hsi::Y; + using hsi::Z; public: @@ -55,11 +58,11 @@ public: /// @param B is the number of samples (bands) along dimension 3 /// @param header_offset is the number of bytes (if any) in the binary header /// @param wavelengths is an STL vector of size B specifying a numerical label for each band - bool open(std::string filename, - unsigned long long X, - unsigned long long Y, - unsigned long long B, - unsigned long long header_offset, + bool open(std::string filename, + unsigned long long X, + unsigned long long Y, + unsigned long long B, + unsigned long long header_offset, std::vector wavelengths){ //copy the wavelengths to the BSQ file structure @@ -366,7 +369,7 @@ public: out.write((char*)band, XYb); if(PROGRESS) progress = (double) (b+1) / (double)B * 50 + 50; } - + } /// Convert the current BSQ file to a BIL file with the specified file name. @@ -489,7 +492,7 @@ public: } //find the indices of the left and right baseline points - while (lab >= w[ai]){ + while (lab >= w[ai]){ ai++; } while (rab <= w[bi]){ @@ -822,13 +825,13 @@ public: i++; //std::cout< band_array, bool PROGRESS = false){ - + std::ofstream out(outfile.c_str(), std::ios::binary); //open the output file for writing file.seekg(0, std::ios::beg); //move to the beginning of the input file @@ -1185,7 +1188,7 @@ public: void deriv(std::string outfile, size_t d, size_t order, const std::vector w = std::vector(), unsigned char* mask = NULL, bool PROGRESS = false){ std::ofstream out(outfile.c_str(), std::ios::binary); //open the output file for writing - + size_t XY = X() * Y(); //calculate the number of values in a band size_t XYb = XY * sizeof(T); //calculate the size of a band (frame) in bytes size_t B = Z(); @@ -1227,7 +1230,7 @@ public: } } } - out.write((char*)outband, XYb); //output the band + out.write((char*)outband, XYb); //output the band if(PROGRESS) progress = (double)(b+1) / (double)B * 100; } diff --git a/stim/envi/hsi.h b/stim/envi/hsi.h index 4a09290..1788cd3 100644 --- a/stim/envi/hsi.h +++ b/stim/envi/hsi.h @@ -32,6 +32,9 @@ protected: std::vector w; //band wavelengths + using binary::R; + using binary::progress; + unsigned long long X(){ return R[O[0]]; } unsigned long long Y(){ return R[O[1]]; } unsigned long long Z(){ return R[O[2]]; } @@ -70,8 +73,8 @@ protected: } /// Get the list of band numbers that bound a list of wavelengths - void band_bounds(std::vector wavelengths, - std::vector& low_bands, + 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 @@ -92,7 +95,7 @@ protected: 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]); } @@ -140,11 +143,11 @@ public: size_t XY = X() * Y(); memset(mask, 255, XY * sizeof(unsigned char)); //initialize the mask to zero (0) T* page = (T*)malloc(R[0] * R[1] * sizeof(T)); //allocate space for a page of data - + for(size_t p = 0; p < R[2]; p++){ //for each page binary::read_page(page, p); //read a page for(size_t i = 0; i < R[0] * R[1]; i++){ //for each pixel in that page - + #ifdef _WIN32 if(!_finite(page[i])){ //if the value at index i is finite #else @@ -171,7 +174,7 @@ public: } /// Calculates the necessary image size required to combine two images given the specified offset (position) of the second image - void calc_combined_size(long long xp, long long yp, long long Sx, long long Sy, + void calc_combined_size(long long xp, long long yp, long long Sx, long long Sy, size_t& Sfx, size_t& Sfy, size_t& p0_x, size_t& p0_y, size_t& p1_x, size_t& p1_y){ @@ -179,7 +182,7 @@ public: calc_combined_size(xp, yp, Sx, Sy, Sfx, Sfy); p0_x = p0_y = p1_x = p1_y = 0; //initialize all boundary positions to zero - + if(xp < 0) p0_x = -xp; //set the left positions of the current and source image else p1_x = xp; if(yp < 0) p0_y = -yp; @@ -191,8 +194,8 @@ public: size_t w = X(); //calculate the number of pixels in a line size_t wb = w * sizeof(T); //calculate the number of bytes in a line - pw = (X() + x0 + x1); //calculate the width of the padded image - pwb = pw * sizeof(T); //calculate the number of bytes in a line of the padded image + size_t pw = (X() + x0 + x1); //calculate the width of the padded image + size_t pwb = pw * sizeof(T); //calculate the number of bytes in a line of the padded image for(size_t y = 0; y < Y(); y++){ //for each line in the real image memcpy( &padded[ (y + y0) * pw + x0 ], &src[y * w], wb ); //use memcpy to copy the line to the appropriate place in the padded image diff --git a/stim/image/image.h b/stim/image/image.h index 898a862..9cb985d 100644 --- a/stim/image/image.h +++ b/stim/image/image.h @@ -32,7 +32,7 @@ class image{ void unalloc(){ //frees any resources associated with the image if(img) free(img); //if memory has been allocated, free it } - + void clear(){ //clears all image data unalloc(); //unallocate previous memory @@ -53,7 +53,7 @@ class image{ size_t idx(size_t x, size_t y, size_t c = 0){ return y * C() * X() + x * C() + c; } - + int cv_type(){ if(std::is_same::value) return CV_MAKETYPE(CV_8U, (int)C()); @@ -63,7 +63,7 @@ class image{ if(std::is_same::value) return CV_MAKETYPE(CV_32S, (int)C()); if(std::is_same::value) return CV_MAKETYPE(CV_32F, (int)C()); if(std::is_same::value) return CV_MAKETYPE(CV_64F, (int)C()); - + std::cout<<"ERROR in stim::image::cv_type - no valid data type found"< max_val){ //if the value is higher than the current max max_val = img[n]; - } + } } - return max; + return max_val; } /// Returns the maximum pixel value in the image @@ -308,10 +308,10 @@ public: for (size_t n=0; n r(X(), Y(), C()); //allocate space for the resulting image for(size_t n = 0; n < N; n++) r.img[n] = white_val - img[n]; //perform the inversion - + return r; //return the inverted image } @@ -330,7 +330,7 @@ public: } image convolve2(image mask){ - + std::cout<<"ERROR stim::image::convolve2 - function has been broken, and shouldn't really be in here."<