From 724ec3470fd1d9a8570e991a2bfa6e5a97e57487 Mon Sep 17 00:00:00 2001 From: David Mayerich Date: Fri, 12 Jun 2015 14:05:12 -0500 Subject: [PATCH] simplified the ENVI classes, moved more stuff into binary.h --- stim/envi/bil.h | 277 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------------------------------------------------------------------------------------------------------------------------------------------------- stim/envi/binary.h | 143 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---- stim/envi/bip.h | 340 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- stim/envi/bsq.h | 149 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------------------------------------------------------------------------- stim/envi/envi.h | 79 +++++++++++++++++++++++++++++++++++-------------------------------------------- 5 files changed, 495 insertions(+), 493 deletions(-) diff --git a/stim/envi/bil.h b/stim/envi/bil.h index 4ba9dd0..a3fb1ba 100644 --- a/stim/envi/bil.h +++ b/stim/envi/bil.h @@ -25,6 +25,16 @@ protected: std::vector w; //band wavelength + unsigned long X(){ + return R[0]; + } + unsigned long Y(){ + return R[2]; + } + unsigned long Z(){ + return R[1]; + } + public: using binary::open; @@ -43,7 +53,7 @@ public: w = wavelengths; - return open(filename, vec(X, Y, B), header_offset); + return open(filename, vec(X, B, Y), header_offset); } @@ -52,19 +62,20 @@ 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_1(p, page); - unsigned int L = R[0] * sizeof(T); //caculate the number of bytes in a sample line - unsigned int jump = R[0] * (R[2] - 1) * sizeof(T); + unsigned int L = X() * sizeof(T); //caculate the number of bytes in a sample line + unsigned int jump = X() * (Z() - 1) * sizeof(T); - if (page >= R[2]){ //make sure the bank number is right + if (page >= Z()){ //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"<::read_line_02(p, x, y); } /// Retrieve a single pixel and stores it in pre-allocated memory. @@ -207,26 +204,16 @@ public: bool pixel(T * p, unsigned n){ //calculate the corresponding x, y - unsigned int x = n % R[0]; - unsigned int y = n / R[0]; + unsigned int x = n % X(); + unsigned int y = n / X(); //get the pixel return spectrum(p, x, y); } //given a Y ,return a XZ slice - bool getY(T * p, unsigned y) - { - if ( y >= R[1]){ //make sure the line number is right - std::cout<<"ERROR: line out of range"<::read_plane_2(p, y); } @@ -242,19 +229,18 @@ public: std::string headername = outname + ".hdr"; //the header file name //simplify image resolution - unsigned int ZX = R[2] * R[0]; //calculate the number of points in a Y slice + 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 = R[2]; - unsigned int X = R[0]; - + 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)); + 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 @@ -267,10 +253,10 @@ public: } // loop start correct every y slice - for (unsigned k =0; k < R[1]; k++) + for (unsigned k =0; k < Y(); k++) { //get the current y slice - getY(c, k); + read_plane_y(c, k); //initialize lownum, highnum, low, high ai = w[0]; @@ -279,16 +265,16 @@ public: //set the baseline point at band 0 to 0 if(wls[0] != w[0]){ bi = wls[control]; - memset(a, (char)0, X * sizeof(T) ); + memset(a, (char)0, X() * sizeof(T) ); } //else get the low band else{ control++; - getYZ(a, c, ai); + read_x_from_xz(a, c, ai); bi = wls[control]; } //get the high band - getYZ(b, c, bi); + read_x_from_xz(b, c, bi); //correct every YZ line @@ -305,7 +291,7 @@ public: ai = bi; bi = wls[control]; - getYZ(b, c, bi); + read_x_from_xz(b, c, bi); } //if the last BL point on the last band of the file? @@ -313,7 +299,7 @@ public: std::swap(a, b); //swap the baseline band pointers - memset(b, (char)0, X * sizeof(T) ); //clear the high band + memset(b, (char)0, X() * sizeof(T) ); //clear the high band ai = bi; bi = w[B - 1]; @@ -322,9 +308,9 @@ public: ci = w[cii]; - unsigned jump = cii * X; + unsigned jump = cii * X(); //perform the baseline correction - for(unsigned i=0; i < X; i++) + for(unsigned i=0; i < X(); i++) { double r = (double) (ci - ai) / (double) (bi - ai); c[i + jump] =(T) ( c[i + jump] - (b[i] - a[i]) * r - a[i] ); @@ -349,11 +335,9 @@ public: /// @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) { - unsigned int B = R[2]; //calculate the number of bands - unsigned int Y = R[1]; - unsigned int X = R[0]; - unsigned int ZX = R[2] * R[0]; - unsigned int XY = R[0] * R[1]; //calculate the number of pixels in a band + 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); @@ -368,17 +352,17 @@ public: band(b, w); //get the certain band into memory - for(unsigned j = 0; j < Y; j++) + for(unsigned j = 0; j < Y(); j++) { - getY(c, j); + read_plane_y(c, j); for(unsigned i = 0; i < B; i++) { - for(unsigned m = 0; m < X; m++) + for(unsigned m = 0; m < X(); m++) { - if( b[m + j * X] < t ) - c[m + i * X] = (T)0.0; + if( b[m + j * X()] < t ) + c[m + i * X()] = (T)0.0; else - c[m + i * X] = c[m + i * X] / b[m + j * X]; + c[m + i * X()] = c[m + i * X()] / b[m + j * X()]; } } target.write(reinterpret_cast(c), L); //write normalized data into destination @@ -395,7 +379,7 @@ public: /// @param outname is the name of the output BSQ file to be saved to disk. bool bsq(std::string outname) { - unsigned int S = R[0] * R[1] * sizeof(T); //calculate the number of bytes in a band + unsigned int 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"; @@ -403,7 +387,7 @@ public: T * p; //pointer to the current band p = (T*)malloc(S); - for ( unsigned i = 0; i < R[2]; i++) + for ( unsigned i = 0; i < Z(); i++) { band_index(p, i); target.write(reinterpret_cast(p), S); //write a band data into target file @@ -419,7 +403,7 @@ public: /// @param outname is the name of the output BIP file to be saved to disk. bool bip(std::string outname) { - unsigned int S = R[0] * R[2] * sizeof(T); //calculate the number of bytes in a ZX slice + unsigned int 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"; @@ -429,14 +413,14 @@ public: T * q; //pointer to the current ZX slice for bip file q = (T*)malloc(S); - for ( unsigned i = 0; i < R[1]; i++) + for ( unsigned i = 0; i < Y(); i++) { - getY(p, i); - for ( unsigned k = 0; k < R[2]; k++) + read_plane_y(p, i); + for ( unsigned k = 0; k < Z(); k++) { - unsigned ks = k * R[0]; - for ( unsigned j = 0; j < R[0]; j++) - q[k + j * R[2]] = p[ks + j]; + unsigned ks = k * X(); + for ( unsigned j = 0; j < X(); j++) + q[k + j * Z()] = p[ks + j]; } target.write(reinterpret_cast(q), S); //write a band data into target file @@ -460,7 +444,7 @@ public: /// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size. bool baseline_band(double lb, double rb, T* lp, T* rp, double wavelength, T* result){ - unsigned XY = R[0] * R[1]; + unsigned XY = X() * Y(); band(result, wavelength); //get band //perform the baseline correction @@ -481,7 +465,7 @@ public: T* lp; T* rp; - unsigned XY = R[0] * R[1]; + unsigned XY = X() * Y(); unsigned S = XY * sizeof(T); lp = (T*) malloc(S); //memory allocation rp = (T*) malloc(S); @@ -512,7 +496,7 @@ public: T* cur; //current band 1 T* cur2; //current band 2 - unsigned XY = R[0] * R[1]; + unsigned XY = X() * Y(); unsigned S = XY * sizeof(T); lp = (T*) malloc(S); //memory allocation @@ -593,14 +577,14 @@ public: /// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size bool ph_to_ph(double lb1, double rb1, double pos1, double lb2, double rb2, double pos2, T * result){ - T* p1 = (T*)malloc(R[0] * R[1] * sizeof(T)); - T* p2 = (T*)malloc(R[0] * R[1] * 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); //calculate the ratio in result - for(unsigned i = 0; i < R[0] * R[1]; i++){ + for(unsigned i = 0; i < X() * Y(); i++){ if(p1[i] == 0 && p2[i] ==0) result[i] = 1; else @@ -624,14 +608,14 @@ public: bool pa_to_ph(double lb1, double rb1, double lab1, double rab1, double lb2, double rb2, double pos, T* result){ - T* p1 = (T*)malloc(R[0] * R[1] * sizeof(T)); - T* p2 = (T*)malloc(R[0] * R[1] * 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); //calculate the ratio in result - for(unsigned i = 0; i < R[0] * R[1]; i++){ + for(unsigned i = 0; i < X() * Y(); i++){ if(p1[i] == 0 && p2[i] ==0) result[i] = 1; else @@ -657,14 +641,14 @@ public: bool pa_to_pa(double lb1, double rb1, double lab1, double rab1, double lb2, double rb2, double lab2, double rab2, T* result){ - T* p1 = (T*)malloc(R[0] * R[1] * sizeof(T)); - T* p2 = (T*)malloc(R[0] * R[1] * 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); //calculate the ratio in result - for(unsigned i = 0; i < R[0] * R[1]; i++){ + for(unsigned i = 0; i < X() * Y(); i++){ if(p1[i] == 0 && p2[i] ==0) result[i] = 1; else @@ -689,7 +673,7 @@ public: T* cur; //current band 1 T* cur2; //current band 2 - unsigned XY = R[0] * R[1]; + unsigned XY = X() * Y(); unsigned S = XY * sizeof(T); lp = (T*) malloc(S); //memory allocation @@ -765,14 +749,14 @@ 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 cpoint(double lb, double rb, double lab, double rab, T* result){ - T* p1 = (T*)malloc(R[0] * R[1] * sizeof(T)); - T* p2 = (T*)malloc(R[0] * R[1] * 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); //calculate the ratio in result - for(unsigned i = 0; i < R[0] * R[1]; i++){ + for(unsigned i = 0; i < X() * Y(); i++){ if(p1[i] == 0 && p2[i] ==0) result[i] = 1; else @@ -793,10 +777,10 @@ public: /// @param p is a pointer to a pre-allocated array at least X * Y in size bool build_mask(double mask_band, double threshold, unsigned char* p){ - T* temp = (T*)malloc(R[0] * R[1] * sizeof(T)); //allocate memory for the certain band + T* temp = (T*)malloc(X() * Y() * sizeof(T)); //allocate memory for the certain band band(temp, mask_band); - for (unsigned int i = 0; i < R[0] * R[1]; i++) { + for (unsigned int i = 0; i < X() * Y(); i++) { if (temp[i] < threshold) p[i] = 0; else @@ -816,20 +800,20 @@ public: std::ofstream target(outfile.c_str(), std::ios::binary); //I THINK THIS IS WRONG - unsigned XZ = R[0] * R[2]; //calculate the number of values in a page on disk + unsigned XZ = X() * Z(); //calculate the number of values in a page on disk unsigned L = XZ * sizeof(T); //calculate the size of the page (in bytes) T * temp = (T*)malloc(L); //allocate memory for a temporary page - for (unsigned i = 0; i < R[1]; i++) //for each value in R[1] (BIP should be X) + for (unsigned i = 0; i < Y(); i++) //for each value in Y() (BIP should be X) { - getY(temp, i); //retrieve an ZX slice, stored in temp - for ( unsigned j = 0; j < R[2]; j++) //for each R[2] (Y) + read_plane_y(temp, i); //retrieve an ZX slice, stored in temp + for ( unsigned j = 0; j < Z(); j++) //for each Z() (Y) { - for (unsigned k = 0; k < R[0]; k++) //for each band + for (unsigned k = 0; k < X(); k++) //for each band { - if(p[i * R[0] + k] == 0) - temp[j * R[0] + k] = 0; + if(p[i * X() + k] == 0) + temp[j * X() + k] = 0; else continue; } @@ -843,30 +827,29 @@ public: ///Saves to disk only those spectra corresponding to mask values != 0 bool sift_mask(std::string outfile, unsigned char* p){ - // Assume R[0] = X, R[1] = Y, R[2] = Z. + // Assume X() = X, Y() = Y, Z() = Z. std::ofstream target(outfile.c_str(), std::ios::binary); //for loading pages: - unsigned XZ = R[0] * R[2]; //calculate the number of values in an XZ page on disk + unsigned XZ = X() * Z(); //calculate the number of values in an XZ page on disk unsigned L = XZ * sizeof(T); //calculate the size of the page (in bytes) T * temp = (T*)malloc(L); //allocate memory for a temporary page //for saving spectra: - unsigned Z = R[2]; //calculate the number of values in a spectrum - unsigned LZ = Z * sizeof(T); //calculate the size of the spectrum (in bytes) + unsigned LZ = Z() * sizeof(T); //calculate the size of the spectrum (in bytes) T * tempZ = (T*)malloc(LZ); //allocate memory for a temporary spectrum - spectrum(tempZ, R[0] - 1, R[1] - 1); //creates a dummy spectrum by taking the last spectrum in the image + spectrum(tempZ, X() - 1, Y() - 1); //creates a dummy spectrum by taking the last spectrum in the image - for (unsigned i = 0; i < R[1]; i++) //Select a page by choosing Y coordinate, R[1] + for (unsigned i = 0; i < Y(); i++) //Select a page by choosing Y coordinate, Y() { - getY(temp, i); //retrieve an ZX page, store in "temp" - for (unsigned j = 0; j < R[0]; j++) //Select a pixel by choosing X coordinate in the page, R[0] + read_plane_y(temp, i); //retrieve an ZX page, store in "temp" + for (unsigned j = 0; j < X(); j++) //Select a pixel by choosing X coordinate in the page, X() { - if (p[j * R[0] + i] != 0) //if the mask != 0 at that XY pixel + if (p[j * X() + i] != 0) //if the mask != 0 at that XY pixel { - for (unsigned k = 0; k < R[2]; k++) //Select a voxel by choosing Z coordinate at the pixel + for (unsigned k = 0; k < Z(); k++) //Select a voxel by choosing Z coordinate at the pixel { - tempZ[k] = temp[k*R[0] + i]; //Pass the correct spectral value from XZ page into the spectrum to be saved. + tempZ[k] = temp[k*X() + i]; //Pass the correct spectral value from XZ page into the spectrum to be saved. } target.write(reinterpret_cast(tempZ), LZ); //write that spectrum to disk. Size is L2. } @@ -883,25 +866,25 @@ public: /// @param p is a pointer to memory of size X * Y * sizeof(T) that will store the band averages. bool band_avg(T* p){ - unsigned long long XZ = R[0] * R[2]; + unsigned long long XZ = X() * Z(); T* temp = (T*)malloc(sizeof(T) * XZ); - T* line = (T*)malloc(sizeof(T) * R[0]); + T* line = (T*)malloc(sizeof(T) * X()); - for (unsigned i = 0; i < R[1]; i++){ + for (unsigned i = 0; i < Y(); i++){ getY(temp, i); //initialize x-line - for (unsigned j = 0; j < R[0]; j++){ + for (unsigned j = 0; j < X(); j++){ line[j] = 0; } unsigned c = 0; - for (unsigned j = 0; j < R[2]; j++){ - for (unsigned k = 0; k < R[0]; k++){ - line[k] += temp[c] / (T)R[2]; + for (unsigned j = 0; j < Z(); j++){ + for (unsigned k = 0; k < X(); k++){ + line[k] += temp[c] / (T)Z(); c++; } } - for (unsigned j = 0; j < R[0]; j++){ - p[j + i * R[0]] = line[j]; + for (unsigned j = 0; j < X(); j++){ + p[j + i * X()] = line[j]; } } free(temp); @@ -913,10 +896,10 @@ 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(T*p, unsigned char* mask){ - unsigned long long XZ = R[0] * R[2]; - unsigned long long XY = R[0] * R[1]; + unsigned long long XZ = X() * Z(); + unsigned long long XY = X() * Y(); T* temp = (T*)malloc(sizeof(T) * XZ); - for (unsigned j = 0; j < R[2]; j++){ + for (unsigned j = 0; j < Z(); j++){ p[j] = 0; } //calculate vaild number in a band @@ -926,13 +909,13 @@ public: count++; } } - for (unsigned k = 0; k < R[1]; k++){ - getY(temp, k); - unsigned kx = k * R[0]; - for (unsigned i = 0; i < R[0]; i++){ + for (unsigned k = 0; k < Y(); k++){ + read_plane_y(temp, k); + unsigned kx = k * X(); + for (unsigned i = 0; i < X(); i++){ if (mask[kx + i] != 0){ - for (unsigned j = 0; j < R[2]; j++){ - p[j] += temp[j * R[0] + i] / (T)count; + for (unsigned j = 0; j < Z(); j++){ + p[j] += temp[j * X() + i] / (T)count; } } } @@ -948,8 +931,8 @@ public: /// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location bool co_matrix(T* co, T* avg, unsigned char *mask){ //memory allocation - unsigned long long xy = R[0] * R[1]; - unsigned int B = R[2]; + unsigned long long xy = X() * Y(); + unsigned int B = Z(); T* temp = (T*)malloc(sizeof(T) * B); //count vaild pixels in a band unsigned count = 0; @@ -999,16 +982,16 @@ public: //calculate the new number of samples and lines unsigned long long sam = x1 - x0; //samples unsigned long long lin = y1 - y0; //lines - unsigned long long L = sam * R[2] * sizeof(T); + unsigned long long L = sam * Z() * sizeof(T); //get specified band and save T* temp = (T*)malloc(L); std::ofstream out(outfile.c_str(), std::ios::binary); - unsigned long long jumpb = (R[0] - sam) * sizeof(T); //jump pointer to the next band + unsigned long long jumpb = (X() - sam) * sizeof(T); //jump pointer to the next band //get start - file.seekg((y0 * R[0] * R[2] + x0) * sizeof(T), std::ios::beg); + file.seekg((y0 * X() * Z() + x0) * sizeof(T), std::ios::beg); for (unsigned i = 0; i < lin; i++) { - for (unsigned j = 0; j < R[2]; j++) + for (unsigned j = 0; j < Z(); j++) { file.read((char *)(temp + j * sam), sizeof(T) * sam); file.seekg(jumpb, std::ios::cur); //go to the next band diff --git a/stim/envi/binary.h b/stim/envi/binary.h index 308e71d..5925cea 100644 --- a/stim/envi/binary.h +++ b/stim/envi/binary.h @@ -81,6 +81,8 @@ protected: return false; } + + public: /// Open a binary file for streaming. @@ -159,8 +161,141 @@ public: return true; } + + + ///Reads a line Z (slowest dimension) for a given XY value + + /// @param p is a pointer to pre-allocated memory equal to the line size R[2] + /// @param x is the x coordinate + /// @param y is the y coordinate + bool read_line_01( T * p, unsigned int x, unsigned int y){ + unsigned int i; + + if ( x >= 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"<= R[1]){ //make sure the bank number is right + std::cout<<"ERROR read_plane_1: page out of range"<= R[0] * R[1] * R[2]){ + std::cout<<"ERROR read_pixel: n is out of range"<= R[0] || y < 0 || y >= R[1] || z < 0 || z > R[2]){ + std::cout<<"ERROR read_pixel: (x,y,z) is out of range"< w; //band wavelength unsigned int offset; //header offset + unsigned long X(){ + return R[1]; + } + unsigned long Y(){ + return R[2]; + } + unsigned long Z(){ + return R[0]; + } + public: using binary::open; using binary::file; using binary::R; + using binary::read_line_12; /// Open a data file for reading using the class interface. @@ -48,7 +59,7 @@ public: //copy the offset to the structure offset = header_offset; - return open(filename, vec(X, Y, B), header_offset); + return open(filename, vec(B, X, Y), header_offset); } @@ -57,21 +68,7 @@ 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){ - - if (page >= R[2]){ //make sure the bank number is right - std::cout<<"ERROR: page out of range"<::read_plane_0(p, page); } /// Retrieve a single band (by numerical label) and stores it in pre-allocated memory. @@ -84,7 +81,7 @@ public: if(w.size() == 0) return band_index(p, (unsigned int)wavelength); - unsigned int XY = R[0] * R[1]; //calculate the number of pixels in a band + unsigned int XY = X() * Y(); //calculate the number of pixels in a band unsigned page=0; //bands around the wavelength @@ -101,8 +98,8 @@ public: { page++; //if wavelength is larger than the last wavelength in header file - if (page == R[2]) { - band_index(p, R[2]-1); + if (page == Z()) { + band_index(p, Z()-1); return true; } } @@ -128,11 +125,24 @@ public: return true; } - //get YZ line from the a Y slice, Y slice data should be already IN the MEMORY - bool getYZ(T* p, T* c, double wavelength) + + /// Retrieve a single spectrum (Z-axis line) at a given (x, y) location and stores it in pre-allocated memory. + + /// @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) + } + + /// Retrieves a band of x values from a given xz plane. + + /// @param p is a pointer to pre-allocated memory at least X * sizeof(T) in size + /// @param c is a pointer to an existing XZ plane (size X*Z*sizeof(T)) + /// @param wavelength is the wavelength of X values to retrieve + bool read_x_from_xz(T* p, T* c, double wavelength) { - unsigned int X = R[0]; //calculate the number of pixels in a sample - unsigned int B = R[2]; + unsigned int B = Z(); unsigned page=0; //samples around the wavelength @@ -141,7 +151,7 @@ public: //if wavelength is smaller than the first one in header file if ( w[page] > wavelength ){ - for(unsigned j = 0; j < R[0]; j++) + for(unsigned j = 0; j < X(); j++) { p[j] = c[j * B]; } @@ -153,7 +163,7 @@ public: page++; //if wavelength is larger than the last wavelength in header file if (page == B) { - for(unsigned j = 0; j < R[0]; j++) + for(unsigned j = 0; j < X(); j++) { p[j] = c[(j + 1) * B - 1]; } @@ -164,20 +174,20 @@ public: { T * p1; T * p2; - p1=(T*)malloc( X * sizeof(T)); //memory allocation - p2=(T*)malloc( X * sizeof(T)); + p1=(T*)malloc( X() * sizeof(T)); //memory allocation + p2=(T*)malloc( X() * sizeof(T)); //band_index(p1, page - 1); - for(unsigned j = 0; j < X; j++) + for(unsigned j = 0; j < X(); j++) { p1[j] = c[j * B + page - 1]; } //band_index(p2, page ); - for(unsigned j = 0; j < X; j++) + for(unsigned j = 0; j < X(); j++) { p2[j] = c[j * B + page]; } - for(unsigned i=0; i < X; i++){ + for(unsigned 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]; } @@ -187,7 +197,7 @@ public: else //if the wavelength is equal to a wavelength in header file { //band_index(p, page); - for(unsigned j = 0; j < X; j++) + for(unsigned j = 0; j < X(); j++) { p[j] = c[j * B + page]; } @@ -195,98 +205,6 @@ public: return true; } - //given a y and a wavelength, return the y-band data - //I do not use it right now, to accelerate the processing speed, I try to read a YZ whole slice into memory first, and then we can call "getYZ" - bool get_x_wavelength(T * p, unsigned y, double wavelength) - { - unsigned int X = R[0]; //calculate the number of pixels in a sample - - unsigned page=0; //samples around the wavelength - T * p1; - T * p2; - - //get the bands numbers around the wavelength - - //if wavelength is smaller than the first one in header file - if ( w[page] > wavelength ){ - file.seekg( y * R[2] * R[0] * sizeof(T), std::ios::beg); - for(unsigned j = 0; j < R[0]; j++) - { - file.read((char *)(p + j), sizeof(T)); - file.seekg((R[2] - 1) * sizeof(T), std::ios::cur); - } - return true; - } - - while( w[page] < wavelength ) - { - page++; - //if wavelength is larger than the last wavelength in header file - if (page == R[2]) { - file.seekg( (y * R[2] * R[0] + R[2] - 1)* sizeof(T), std::ios::beg); - for(unsigned j = 0; j < R[0]; j++) - { - file.read((char *)(p + j), sizeof(T)); - file.seekg((R[2] - 1) * sizeof(T), std::ios::cur); - } - return true; - } - } - if ( wavelength < w[page] ) - { - p1=(T*)malloc( X * sizeof(T)); //memory allocation - p2=(T*)malloc( X * sizeof(T)); - //band_index(p1, page - 1); - file.seekg( (y * R[2] * R[0] + page - 1)* sizeof(T), std::ios::beg); - for(unsigned j = 0; j < R[0]; j++) - { - file.read((char *)(p1 + j), sizeof(T)); - file.seekg((R[2] - 1) * sizeof(T), std::ios::cur); - } - //band_index(p2, page ); - file.seekg( (y * R[2] * R[0] + page)* sizeof(T), std::ios::beg); - for(unsigned j = 0; j < R[0]; j++) - { - file.read((char *)(p2 + j), sizeof(T)); - file.seekg((R[2] - 1) * sizeof(T), std::ios::cur); - } - - for(unsigned i=0; i < R[0]; i++){ - double r = (double) (wavelength - w[page-1]) / (double) (w[page] - w[page-1]); - p[i] = (p2[i] - p1[i]) * r + p1[i]; - } - } - else //if the wavelength is equal to a wavelength in header file - { - //band_index(p, page); - file.seekg( (y * R[2] * R[0] + page) * sizeof(T), std::ios::beg); - for(unsigned j = 0; j < R[0]; j++) - { - file.read((char *)(p + j), sizeof(T)); - file.seekg((R[2] - 1) * sizeof(T), std::ios::cur); - } - } - return true; - } - - /// Retrieve a single spectrum (B-axis line) at a given (x, y) location and stores it in pre-allocated memory. - - /// @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){ - - if ( x >= R[0] || y >= R[1]){ //make sure the sample and line number is right - std::cout<<"ERROR: sample or line out of range"<= bandnum){ //make sure the pixel number is right std::cout<<"ERROR: sample or line out of range"<= R[1]){ //make sure the line number is right - std::cout<<"ERROR: line out of range"<::read_plane_2(p, y); } /// Perform baseline correction given a list of baseline points and stores the result in a new BSQ file. @@ -332,10 +240,9 @@ public: std::string headername = outname + ".hdr"; //the header file name //simplify image resolution - unsigned int ZX = R[2] * R[0]; //calculate the number of points in a Y slice + 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 = R[2]; - unsigned int X = R[0]; + unsigned int B = Z(); T* c; //pointer to the current Y slice c = (T*)malloc(L); //memory allocation @@ -343,8 +250,8 @@ public: 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)); + 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 @@ -357,10 +264,10 @@ public: } // loop start correct every y slice - for (unsigned k =0; k < R[1]; k++) + for (unsigned k =0; k < Y(); k++) { //get the current y slice - getY(c, k); + read_plane_y(c, k); //initialize lownum, highnum, low, high control=0; @@ -369,16 +276,16 @@ public: //set the baseline point at band 0 to 0 if(wls[0] != w[0]){ bi = wls[control]; - memset(a, (char)0, X * sizeof(T) ); + memset(a, (char)0, X() * sizeof(T) ); } //else get the low band else{ control++; - getYZ(a, c, ai); + read_x_from_xz(a, c, ai); bi = wls[control]; } //get the high band - getYZ(b, c, bi); + read_x_from_xz(b, c, bi); //correct every YZ line @@ -394,7 +301,7 @@ public: ai = bi; bi = wls[control]; - getYZ(b, c, bi); + read_x_from_xz(b, c, bi); } //if the last BL point on the last band of the file? @@ -402,7 +309,7 @@ public: std::swap(a, b); //swap the baseline band pointers - memset(b, (char)0, X * sizeof(T) ); //clear the high band + memset(b, (char)0, X() * sizeof(T) ); //clear the high band ai = bi; bi = w[B - 1]; @@ -412,7 +319,7 @@ public: ci = w[cii]; //perform the baseline correction - for(unsigned i=0; i < X; i++) + 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] ); @@ -439,11 +346,9 @@ public: /// @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) { - unsigned int B = R[2]; //calculate the number of bands - unsigned int Y = R[1]; - unsigned int X = R[0]; - unsigned int ZX = R[2] * R[0]; - unsigned int XY = R[0] * R[1]; //calculate the number of pixels in a band + 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); @@ -458,11 +363,11 @@ public: band(b, w); //get the certain band into memory - for(unsigned j = 0; j < Y; j++) + for(unsigned j = 0; j < Y(); j++) { - getY(c, j); - unsigned jX = j * X; //to avoid calculating it many times - for(unsigned i = 0; i < X; i++) + 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++) @@ -494,7 +399,7 @@ public: bil(temp); stim::bil n; - if(n.open(temp, R[0], R[1], R[2], offset, w)==false){ //open infile + 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 } @@ -550,7 +455,7 @@ public: /// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size. bool baseline_band(double lb, double rb, T* lp, T* rp, double wavelength, T* result){ - unsigned XY = R[0] * R[1]; + unsigned XY = X() * Y(); band(result, wavelength); //get band //perform the baseline correction @@ -571,7 +476,7 @@ public: T* lp; T* rp; - unsigned XY = R[0] * R[1]; + unsigned XY = X() * Y(); unsigned S = XY * sizeof(T); lp = (T*) malloc(S); //memory allocation rp = (T*) malloc(S); @@ -601,7 +506,7 @@ public: T* cur; //current band 1 T* cur2; //current band 2 - unsigned XY = R[0] * R[1]; + unsigned XY = X() * Y(); unsigned S = XY * sizeof(T); lp = (T*) malloc(S); //memory allocation @@ -682,14 +587,14 @@ public: /// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size bool ph_to_ph(double lb1, double rb1, double pos1, double lb2, double rb2, double pos2, T * result){ - T* p1 = (T*)malloc(R[0] * R[1] * sizeof(T)); - T* p2 = (T*)malloc(R[0] * R[1] * 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); //calculate the ratio in result - for(unsigned i = 0; i < R[0] * R[1]; i++){ + for(unsigned i = 0; i < X() * Y(); i++){ if(p1[i] == 0 && p2[i] ==0) result[i] = 1; else @@ -713,14 +618,14 @@ public: bool pa_to_ph(double lb1, double rb1, double lab1, double rab1, double lb2, double rb2, double pos, T* result){ - T* p1 = (T*)malloc(R[0] * R[1] * sizeof(T)); - T* p2 = (T*)malloc(R[0] * R[1] * 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); //calculate the ratio in result - for(unsigned i = 0; i < R[0] * R[1]; i++){ + for(unsigned i = 0; i < X() * Y(); i++){ if(p1[i] == 0 && p2[i] ==0) result[i] = 1; else @@ -746,14 +651,14 @@ public: bool pa_to_pa(double lb1, double rb1, double lab1, double rab1, double lb2, double rb2, double lab2, double rab2, T* result){ - T* p1 = (T*)malloc(R[0] * R[1] * sizeof(T)); - T* p2 = (T*)malloc(R[0] * R[1] * 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); //calculate the ratio in result - for(unsigned i = 0; i < R[0] * R[1]; i++){ + for(unsigned i = 0; i < X() * Y(); i++){ if(p1[i] == 0 && p2[i] ==0) result[i] = 1; else @@ -778,7 +683,7 @@ public: T* cur; //current band 1 T* cur2; //current band 2 - unsigned XY = R[0] * R[1]; + unsigned XY = X() * Y(); unsigned S = XY * sizeof(T); lp = (T*) malloc(S); //memory allocation @@ -854,14 +759,14 @@ 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 cpoint(double lb, double rb, double lab, double rab, T* result){ - T* p1 = (T*)malloc(R[0] * R[1] * sizeof(T)); - T* p2 = (T*)malloc(R[0] * R[1] * 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); //calculate the ratio in result - for(unsigned i = 0; i < R[0] * R[1]; i++){ + for(unsigned i = 0; i < X() * Y(); i++){ if(p1[i] == 0 && p2[i] ==0) result[i] = 1; else @@ -882,10 +787,10 @@ public: /// @param p is a pointer to a pre-allocated array at least X * Y in size bool build_mask(double mask_band, double threshold, unsigned char* p){ - T* temp = (T*)malloc(R[0] * R[1] * sizeof(T)); //allocate memory for the certain band + T* temp = (T*)malloc(X() * Y() * sizeof(T)); //allocate memory for the certain band band(temp, mask_band); - for (unsigned int i = 0; i < R[0] * R[1];i++) { + for (unsigned int i = 0; i < X() * Y();i++) { if (temp[i] < threshold) p[i] = 0; else @@ -905,20 +810,20 @@ public: std::ofstream target(outfile.c_str(), std::ios::binary); - unsigned ZX = R[2] * R[0]; //calculate the number of values in a page (XZ in BIP) + unsigned ZX = Z() * X(); //calculate the number of values in a page (XZ in BIP) unsigned L = ZX * sizeof(T); //calculate the number of bytes in a page T * temp = (T*)malloc(L); //allocate space for that page - for (unsigned i = 0; i < R[1]; i++) //for each page (Y in BIP) + for (unsigned i = 0; i < Y(); i++) //for each page (Y in BIP) { - getY(temp, i); //load that page (it's pointed to by temp) - for ( unsigned j = 0; j < R[0]; j++) //for each X value + read_plane_y(temp, i); //load that page (it's pointed to by temp) + for ( unsigned j = 0; j < X(); j++) //for each X value { - for (unsigned k = 0; k < R[2]; k++) //for each B value (band) + for (unsigned k = 0; k < Z(); k++) //for each B value (band) { - if (p[i * R[0] + j] == 0) //if the mask value is zero - temp[j * R[2] + k] = 0; //set the pixel value to zero + if (p[i * X() + j] == 0) //if the mask value is zero + temp[j * Z() + k] = 0; //set the pixel value to zero else //otherwise just continue continue; } @@ -933,30 +838,29 @@ public: ///Saves to disk only those spectra corresponding to mask values != 0 bool sift_mask(std::string outfile, unsigned char* p){ - // Assume R[0] = X, R[1] = Y, R[2] = Z. + // Assume X() = X, Y() = Y, Z() = Z. std::ofstream target(outfile.c_str(), std::ios::binary); //for loading pages: - unsigned ZX = R[0] * R[2]; //calculate the number of values in an XZ page on disk + unsigned ZX = X() * Z(); //calculate the number of values in an XZ page on disk unsigned L = ZX * sizeof(T); //calculate the size of the page (in bytes) T * temp = (T*)malloc(L); //allocate memory for a temporary page //for saving spectra: - unsigned Z = R[2]; //calculate the number of values in a spectrum - unsigned LZ = Z * sizeof(T); //calculate the size of the spectrum (in bytes) + unsigned LZ = Z() * sizeof(T); //calculate the size of the spectrum (in bytes) T * tempZ = (T*)malloc(LZ); //allocate memory for a temporary spectrum - spectrum(tempZ, R[0] - 1, R[1] - 1); //creates a dummy spectrum by taking the last spectrum in the image + spectrum(tempZ, X() - 1, Y() - 1); //creates a dummy spectrum by taking the last spectrum in the image - for (unsigned i = 0; i < R[1]; i++) //Select a page by choosing Y coordinate, R[1] + for (unsigned i = 0; i < Y(); i++) //Select a page by choosing Y coordinate, Y() { - getY(temp, i); //retrieve an ZX page, store in "temp" - for (unsigned j = 0; j < R[0]; j++) //Select a pixel by choosing X coordinate in the page, R[0] + read_plane_y(temp, i); //retrieve an ZX page, store in "temp" + for (unsigned j = 0; j < X(); j++) //Select a pixel by choosing X coordinate in the page, X() { - if (p[j * R[0] + i] != 0) //if the mask != 0 at that XY pixel + if (p[j * X() + i] != 0) //if the mask != 0 at that XY pixel { - for (unsigned k = 0; k < R[2]; k++) //Select a voxel by choosing Z coordinate at the pixel + for (unsigned k = 0; k < Z(); k++) //Select a voxel by choosing Z coordinate at the pixel { - tempZ[k] = temp[i*R[2] + k]; //Pass the correct spectral value from XZ page into the spectrum to be saved. + tempZ[k] = temp[i*Z() + k]; //Pass the correct spectral value from XZ page into the spectrum to be saved. } target.write(reinterpret_cast(tempZ), LZ); //write that spectrum to disk. Size is L2. } @@ -973,16 +877,16 @@ public: /// @param p is a pointer to memory of size X * Y * sizeof(T) that will store the band averages. bool band_avg(T* p){ - unsigned long long XY = R[0] * R[1]; + unsigned long long XY = X() * Y(); //get every pixel and calculate average value - T* temp = (T*)malloc(sizeof(T) * R[2]); + T* temp = (T*)malloc(sizeof(T) * Z()); T sum; for (unsigned i = 0; i < XY; i++){ pixel(temp, i); //calculate the sum value of every value sum = 0; //initialize sum value - for (unsigned j = 0; j < R[2]; j++){ - sum += temp[j]/(T)R[2]; + for (unsigned j = 0; j < Z(); j++){ + sum += temp[j]/(T)Z(); } p[i] = sum; } @@ -995,10 +899,10 @@ 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(T*p, unsigned char* mask){ - unsigned long long XY = R[0] * R[1]; - T* temp = (T*)malloc(sizeof(T) * R[2]); + unsigned long long XY = X() * Y(); + T* temp = (T*)malloc(sizeof(T) * Z()); //Iinitialize - for (unsigned j = 0; j < R[2]; j++){ + for (unsigned j = 0; j < Z(); j++){ p[j] = 0; } //calculate vaild number in a band @@ -1012,7 +916,7 @@ public: for (unsigned i = 0; i < XY; i++){ if (mask[i] != 0){ pixel(temp, i); - for (unsigned j = 0; j < R[2]; j++){ + for (unsigned j = 0; j < Z(); j++){ p[j] += temp[j] / (T)count; } } @@ -1028,8 +932,8 @@ public: /// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location bool co_matrix(T* co, T* avg, unsigned char *mask){ //memory allocation - unsigned long long xy = R[0] * R[1]; - unsigned int B = R[2]; + unsigned long long xy = X() * Y(); + unsigned int B = Z(); T* temp = (T*)malloc(sizeof(T) * B); //count vaild pixels in a band unsigned count = 0; @@ -1079,17 +983,17 @@ public: //calculate the new number of samples and lines unsigned long long sam = x1 - x0; //samples unsigned long long lin = y1 - y0; //lines - unsigned long long L = R[2] * sizeof(T); + unsigned long long L = Z() * sizeof(T); //get specified band and save T* temp = (T*)malloc(L); std::ofstream out(outfile.c_str(), std::ios::binary); //get start - unsigned long long sp = y0 * R[0] + x0; //start pixel + unsigned long long sp = y0 * X() + x0; //start pixel for (unsigned i = 0; i < lin; i++) { for (unsigned j = 0; j < sam; j++) { - pixel(temp, sp + j + i * R[0]); + pixel(temp, sp + j + i * X()); out.write(reinterpret_cast(temp), L); //write slice data into target file } } diff --git a/stim/envi/bsq.h b/stim/envi/bsq.h index fc82f5e..f13c3e4 100644 --- a/stim/envi/bsq.h +++ b/stim/envi/bsq.h @@ -30,12 +30,26 @@ protected: std::vector w; //band wavelengths unsigned int offset; + using binary::R; + + unsigned long X(){ + return R[0]; + } + unsigned long Y(){ + return R[1]; + } + unsigned long Z(){ + return R[2]; + } + public: using binary::open; using binary::file; + using binary::read_line_01; + using binary::read_plane_2; //using binary::getSlice; - using binary::R; + /// Open a data file for reading using the class interface. @@ -61,17 +75,7 @@ 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){ - - 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"<= bandnum){ //make sure the pixel number is right std::cout<<"ERROR: sample or line out of range"< n; - if(n.open(temp, R[0], R[1], R[2], offset, w)==false){ //open infile + if(n.open(temp, X(), Y(), Z(), offset, w)==false){ //open infile std::cout<<"ERROR: unable to open input file"<(p), L); //write XZ slice data into target file @@ -391,7 +380,7 @@ public: /// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size. bool baseline_band(double lb, double rb, T* lp, T* rp, double wavelength, T* result){ - unsigned XY = R[0] * R[1]; + unsigned XY = X() * Y(); band(result, wavelength); //get band //perform the baseline correction @@ -412,7 +401,7 @@ public: T* lp; T* rp; - unsigned XY = R[0] * R[1]; + unsigned XY = X() * Y(); unsigned S = XY * sizeof(T); lp = (T*) malloc(S); //memory allocation rp = (T*) malloc(S); @@ -442,7 +431,7 @@ public: T* cur; //current band 1 T* cur2; //current band 2 - unsigned XY = R[0] * R[1]; + unsigned XY = X() * Y(); unsigned S = XY * sizeof(T); lp = (T*) malloc(S); //memory allocation @@ -523,14 +512,14 @@ public: /// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size bool ph_to_ph(double lb1, double rb1, double pos1, double lb2, double rb2, double pos2, T * result){ - T* p1 = (T*)malloc(R[0] * R[1] * sizeof(T)); - T* p2 = (T*)malloc(R[0] * R[1] * 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); //calculate the ratio in result - for(unsigned i = 0; i < R[0] * R[1]; i++){ + for(unsigned i = 0; i < X() * Y(); i++){ if(p1[i] == 0 && p2[i] ==0) result[i] = 1; else @@ -554,14 +543,14 @@ public: bool pa_to_ph(double lb1, double rb1, double lab1, double rab1, double lb2, double rb2, double pos, T* result){ - T* p1 = (T*)malloc(R[0] * R[1] * sizeof(T)); - T* p2 = (T*)malloc(R[0] * R[1] * 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); //calculate the ratio in result - for(unsigned i = 0; i < R[0] * R[1]; i++){ + for(unsigned i = 0; i < X() * Y(); i++){ if(p1[i] == 0 && p2[i] ==0) result[i] = 1; else @@ -587,14 +576,14 @@ public: bool pa_to_pa(double lb1, double rb1, double lab1, double rab1, double lb2, double rb2, double lab2, double rab2, T* result){ - T* p1 = (T*)malloc(R[0] * R[1] * sizeof(T)); - T* p2 = (T*)malloc(R[0] * R[1] * 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); //calculate the ratio in result - for(unsigned i = 0; i < R[0] * R[1]; i++){ + for(unsigned i = 0; i < X() * Y(); i++){ if(p1[i] == 0 && p2[i] ==0) result[i] = 1; else @@ -619,7 +608,7 @@ public: T* cur; //current band 1 T* cur2; //current band 2 - unsigned XY = R[0] * R[1]; + unsigned XY = X() * Y(); unsigned S = XY * sizeof(T); lp = (T*) malloc(S); //memory allocation @@ -695,14 +684,14 @@ 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 cpoint(double lb, double rb, double lab, double rab, T* result){ - T* p1 = (T*)malloc(R[0] * R[1] * sizeof(T)); - T* p2 = (T*)malloc(R[0] * R[1] * 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); //calculate the ratio in result - for(unsigned i = 0; i < R[0] * R[1]; i++){ + for(unsigned i = 0; i < X() * Y(); i++){ if(p1[i] == 0 && p2[i] ==0) result[i] = 1; else @@ -723,10 +712,10 @@ public: /// @param p is a pointer to a pre-allocated array at least X * Y in size bool build_mask(double mask_band, double threshold, unsigned char* p = NULL){ - T* temp = (T*)malloc(R[0] * R[1] * sizeof(T)); //allocate memory for the certain band + T* temp = (T*)malloc(X() * Y() * sizeof(T)); //allocate memory for the certain band band(temp, mask_band); - for (unsigned int i = 0; i < R[0] * R[1]; i++) { + for (unsigned int i = 0; i < X() * Y(); i++) { if (temp[i] < threshold) p[i] = 0; else @@ -746,12 +735,12 @@ public: std::ofstream target(outfile.c_str(), std::ios::binary); - unsigned XY = R[0] * R[1]; //calculate number of a band + unsigned XY = X() * Y(); //calculate number of a band unsigned L = XY * sizeof(T); T * temp = (T*)malloc(L); - for (unsigned i = 0; i < R[2]; i++) //for each spectral bin + for (unsigned i = 0; i < Z(); i++) //for each spectral bin { band_index(temp, i); //get the specified band (by index) for ( unsigned j = 0; j < XY; j++) // for each pixel @@ -774,13 +763,13 @@ public: bool sift_mask(std::string outfile, unsigned char* p){ std::ofstream target(outfile.c_str(), std::ios::binary); // open a band (XY plane) - unsigned XY = R[0] * R[1]; //Number of XY pixels + unsigned XY = X() * Y(); //Number of XY pixels unsigned L = XY * sizeof(T); //size of XY pixels T * temp = (T*)malloc(L); //allocate memory for a band T * temp_vox = (T*)malloc(sizeof(T)); //allocate memory for one voxel - for (unsigned i = 0; i < R[2]; i++) //for each spectral bin + for (unsigned i = 0; i < Z(); i++) //for each spectral bin { band_index(temp, i); //get the specified band (XY sheet by index) for (unsigned j = 0; j < XY; j++) // for each pixel @@ -804,18 +793,18 @@ public: /// @param p is a pointer to memory of size X * Y * sizeof(T) that will store the band averages. bool band_avg(T* p){ - unsigned long long XY = R[0] * R[1]; + unsigned long long XY = X() * Y(); T* temp = (T*)malloc(sizeof(T) * XY); //initialize p band_index(p, 0); for (unsigned j = 0; j < XY; j++){ - p[j] /= (T)R[2]; + p[j] /= (T)Z(); } //get every band and add them all - for (unsigned i = 1; i < R[2]; i++){ + for (unsigned i = 1; i < Z(); i++){ band_index(temp, i); for (unsigned j = 0; j < XY; j++){ - p[j] += temp[j]/(T)R[2]; + p[j] += temp[j]/(T)Z(); } } free(temp); @@ -827,7 +816,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(T*p, unsigned char* mask){ - unsigned long long XY = R[0] * R[1]; + unsigned long long XY = X() * Y(); unsigned count = 0; //count will store the number of masked pixels T* temp = (T*)malloc(sizeof(T) * XY); //calculate this loop counts the number of true pixels in the mask @@ -836,9 +825,9 @@ public: count++; } } - //this loops goes through each band in B (R[2]) + //this loops goes through each band in B (Z()) // masked (or valid) pixels from that band are averaged and the average is stored in p - for (unsigned i = 0; i < R[2]; i++){ + for (unsigned i = 0; i < Z(); i++){ p[i] = 0; band_index(temp, i); //get the band image and store it in temp for (unsigned j = 0; j < XY; j++){ //loop through temp, averaging valid pixels @@ -858,8 +847,8 @@ public: /// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location bool co_matrix(T* co, T* avg, unsigned char *mask){ //memory allocation - unsigned long long xy = R[0] * R[1]; - unsigned int B = R[2]; + unsigned long long xy = X() * Y(); + unsigned int B = Z(); T* bandi = (T*)malloc(sizeof(T) * xy); T* bandj = (T*)malloc(sizeof(T) * xy); @@ -909,11 +898,11 @@ public: //get specified band and save T* temp = (T*)malloc(L); std::ofstream out(outfile.c_str(), std::ios::binary); - unsigned long long jumpb = R[0] * (R[1] - lin) * sizeof(T); //jump pointer to the next band - unsigned long long jumpl = (R[0] - sam) * sizeof(T); //jump pointer to the next line + unsigned long long jumpb = X() * (Y() - lin) * sizeof(T); //jump pointer to the next band + unsigned long long jumpl = (X() - sam) * sizeof(T); //jump pointer to the next line //get start - file.seekg((y0 * R[0] + x0) * sizeof(T), std::ios::beg); - for (unsigned i = 0; i < R[2]; i++) + file.seekg((y0 * X() + x0) * sizeof(T), std::ios::beg); + for (unsigned i = 0; i < Z(); i++) { for (unsigned j = 0; j < lin; j++) { diff --git a/stim/envi/envi.h b/stim/envi/envi.h index f1dddaa..71e34c1 100644 --- a/stim/envi/envi.h +++ b/stim/envi/envi.h @@ -710,6 +710,41 @@ public: return false; } + bool spectrum(void* ptr, unsigned int x, unsigned int y){ + + 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); + else if (header.data_type == envi_header::float64) + return ((bsq*)file)->spectrum((double*)ptr, x, y); + 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)->spectrum((float*)ptr, x, y); + else if (header.data_type == envi_header::float64) + return ((bil*)file)->spectrum((double*)ptr, x, y); + else{ + std::cout << "ERROR: unidentified data type" << std::endl; + exit(1); + } + } + else if (header.interleave == envi_header::BIP){ + if (header.data_type == envi_header::float32) + return ((bip*)file)->spectrum((float*)ptr, x, y); + else if (header.data_type == envi_header::float64) + return ((bip*)file)->spectrum((double*)ptr, x, y); + else{ + std::cout << "ERROR: unidentified data type" << std::endl; + exit(1); + } + } + return false; + } + /// Retrieve a single band (based on index) 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. @@ -762,50 +797,6 @@ public: return true; } - //p:start positon; N: number of pixels saved in X; - bool feature_matrix(void * X, unsigned char * mask, unsigned start, unsigned N) - { - //save pixels in X as floating numbers: float * p - float * p = (float*)malloc(header.bands * sizeof(float)); //save pixel information - unsigned pixels = header.samples * header.lines; //calculate pixel numbers in a band - unsigned count = 0; //for counting use - unsigned j = 0; //memory the pointer location in X - - //create two indices into the mask (mask index and pixel index) - unsigned mi = 0; //valid pixel index - unsigned pi = 0; //actual pixel in the mask - - //find the actual pixel index for the mask index "start" - while(mi < start){ - if(mask[pi]) - mi++; - - pi++; - } - - for(unsigned i = pi; i < pixels; i++){ - if(mask[i] != 0){ - pixel(p, i); - //copy p to X - for(unsigned k = 0; k < header.bands; k++){ - ((float *)X)[j] = p[k]; - j++; - } - count++; - if(count == N) - break; - } - else - continue; - } - if(count < N){ - std::cout << "number of valid pixels in the mask : " << count <<"is less than N: "<< N; - exit(1); - } - free(p); - return true; - } - /// Calculate the mean value for all masked (or valid) pixels in a band and returns the average spectrum /// @param p is a pointer to pre-allocated memory of size [B * sizeof(T)] that stores the mean spectrum -- libgit2 0.21.4