From ba51ae6a3d666879ecea0dd9fe40a56b7ca47092 Mon Sep 17 00:00:00 2001 From: David Date: Mon, 25 Apr 2016 01:29:38 -0500 Subject: [PATCH] fixed metric calculations, added a finite-value mask, as_float() in stim::arglist returns a double --- stim/envi/bil.h | 52 ++++++++++++++++++++++++++++++++++++++++------------ stim/envi/bip.h | 49 +++++++++++++++++++++++++++++++++++++------------ stim/envi/bsq.h | 108 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++---------------------------------------------------- stim/envi/envi.h | 126 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------------------- stim/envi/hsi.h | 46 ++++++++++++++++++++++++++++++++++++++++++++++ stim/parser/arguments.h | 2 +- 6 files changed, 265 insertions(+), 118 deletions(-) diff --git a/stim/envi/bil.h b/stim/envi/bil.h index 1cc2646..d9874a8 100644 --- a/stim/envi/bil.h +++ b/stim/envi/bil.h @@ -600,7 +600,7 @@ public: /// @param rb2 is the label value for the right baseline point for the second peak (denominator) /// @param pos2 is the label value for the second peak (denominator) position /// @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){ + 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)); @@ -630,8 +630,8 @@ public: /// @param rb2 is the label value for the right baseline point for the second peak (denominator) /// @param pos2 is the label value for the second peak (denominator) position /// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size - bool pa_to_ph(double lb1, double rb1, double lab1, double rab1, - double lb2, double rb2, double pos, T* result){ + 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)); @@ -663,8 +663,8 @@ public: /// @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 result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size - bool pa_to_pa(double lb1, double rb1, double lab1, double rab1, - double lb2, double rb2, double lab2, double rab2, T* result){ + 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)); @@ -773,7 +773,7 @@ public: /// @param lab is the label for the start of the peak /// @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){ + 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)); @@ -782,9 +782,7 @@ public: area(lb, rb, lab, rab, p2); //calculate the ratio in result for(unsigned long long i = 0; i < X() * Y(); i++){ - if(p1[i] == 0 && p2[i] ==0) - result[i] = 1; - else + if(mask == NULL || mask[i]) result[i] = p1[i] / p2[i]; } @@ -800,16 +798,16 @@ public: /// @param mask_band is the band used to specify the mask /// @param threshold is the threshold used to determine if the mask value is true or false /// @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, bool PROGRESS = false){ + bool build_mask(unsigned char* mask, double mask_band, double threshold, bool PROGRESS = false){ T* temp = (T*)malloc(X() * Y() * sizeof(T)); //allocate memory for the certain band band(temp, mask_band, PROGRESS); for (unsigned long long i = 0; i < X() * Y(); i++) { if (temp[i] < threshold) - p[i] = 0; + mask[i] = 0; else - p[i] = 255; + mask[i] = 255; } free(temp); @@ -850,6 +848,36 @@ public: return true; } + /// Copies all spectra corresponding to nonzero values of a mask into a pre-allocated matrix of size (B x P) + /// where P is the number of masked pixels and B is the number of bands. The allocated memory can be accessed + /// using the following indexing: i = p*B + b + /// @param matrix is the destination for the pixel data + /// @param mask is the mask + bool sift(T* matrix, unsigned char* mask){ + size_t Lbytes = sizeof(T) * X(); + 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 + for(size_t b = 0; b < Z(); b++){ //for each band in the data set + pl = 0; //initialize the pixel offset for the current line to zero (0) + file.read( (char*)line, Lbytes ); //read the current line + for(size_t x = 0; x < X(); x++){ + if(mask[y * X() + x]){ //if the current pixel is masked + size_t i = (p + pl) * Z() + b; //calculate the index in the sifted matrix + matrix[i] = line[x]; //store the current value in the line at the correct matrix location + pl++; //increment the pixel pointer + } + } + } + p += pl; //add the line increment to the running pixel index + } + return true; + } + /// Saves to disk only those spectra corresponding to mask values != 0 /// Unlike the BIP and BSQ versions of this function, this version saves a different format: the BIL file is saved as a BIP bool sift(std::string outfile, unsigned char* p, bool PROGRESS = false){ diff --git a/stim/envi/bip.h b/stim/envi/bip.h index 678201d..9807229 100644 --- a/stim/envi/bip.h +++ b/stim/envi/bip.h @@ -526,7 +526,7 @@ public: /// @param rb2 is the label value for the right baseline point for the second peak (denominator) /// @param pos2 is the label value for the second peak (denominator) position /// @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){ + 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)); @@ -556,8 +556,8 @@ public: /// @param rb2 is the label value for the right baseline point for the second peak (denominator) /// @param pos2 is the label value for the second peak (denominator) position /// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size - bool pa_to_ph(double lb1, double rb1, double lab1, double rab1, - double lb2, double rb2, double pos, T* result){ + 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)); @@ -589,8 +589,8 @@ public: /// @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 result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size - bool pa_to_pa(double lb1, double rb1, double lab1, double rab1, - double lb2, double rb2, double lab2, double rab2, T* result){ + 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)); @@ -699,7 +699,7 @@ public: /// @param lab is the label for the start of the peak /// @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){ + 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)); @@ -708,9 +708,7 @@ public: area(lb, rb, lab, rab, p2); //calculate the ratio in result for(unsigned long long i = 0; i < X() * Y(); i++){ - if(p1[i] == 0 && p2[i] ==0) - result[i] = 1; - else + if(mask == NULL || mask[i]) result[i] = p1[i] / p2[i]; } @@ -726,16 +724,16 @@ public: /// @param mask_band is the band used to specify the mask /// @param threshold is the threshold used to determine if the mask value is true or false /// @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, bool PROGRESS = false){ + bool build_mask(unsigned char* mask, double mask_band, double threshold, bool PROGRESS = false){ T* temp = (T*)malloc(X() * Y() * sizeof(T)); //allocate memory for the certain band band(temp, mask_band, PROGRESS); for (unsigned long long i = 0; i < X() * Y();i++) { if (temp[i] < threshold) - p[i] = 0; + mask[i] = 0; else - p[i] = 255; + mask[i] = 255; } free(temp); @@ -776,6 +774,33 @@ public: return true; //return true } + /// Copies all spectra corresponding to nonzero values of a mask into a pre-allocated matrix of size (B x P) + /// where P is the number of masked pixels and B is the number of bands. The allocated memory can be accessed + /// using the following indexing: i = p*B + b + /// @param matrix is the destination for the pixel data + /// @param mask is the mask + bool sift(T* matrix, unsigned char* mask){ + size_t Bbytes = sizeof(T) * Z(); + size_t XY = X() * Y(); + 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 + file.read( (char*)band, Bbytes ); //read the current line + for(size_t b = 0; b < Z(); b++){ //copy each band value to the sifted matrix + size_t i = p * Z() + b; //calculate the index in the sifted matrix + matrix[i] = band[b]; //store the current value in the line at the correct matrix location + } + p++; //increment the pixel pointer + } + else + file.seekg(Bbytes, std::ios::cur); //otherwise skip this band + } + return true; + } /// Saves to disk only those spectra corresponding to mask values != 0 bool sift(std::string outfile, unsigned char* mask, bool PROGRESS = false){ diff --git a/stim/envi/bsq.h b/stim/envi/bsq.h index bd6b529..b7f4e45 100644 --- a/stim/envi/bsq.h +++ b/stim/envi/bsq.h @@ -436,9 +436,7 @@ public: cur = (T*) malloc(S); cur2 = (T*) malloc(S); - memset(result, (char)0, S); - - //find the wavelenght position in the whole band + //find the wavelength position in the whole band unsigned long long n = w.size(); unsigned long long ai = 0; //left bound position unsigned long long bi = n - 1; //right bound position @@ -450,26 +448,30 @@ public: std::cout<<"ERROR: left bound or right bound out of bandwidth"< rb){ std::cout<<"ERROR: right bound should be bigger than left bound"<= w[ai]){ + //find the indices of the left and right baseline points + while (lab >= w[ai]){ ai++; } while (rab <= w[bi]){ bi--; } - band(lp, lb); + band(lp, lb); //get the band images for the left and right baseline points band(rp, rb); - //calculate the beginning and the ending part - baseline_band(lb, rb, lp, rp, rab, cur2); //ending part - baseline_band(lb, rb, lp, rp, w[bi], cur); + // calculate the average value of the indexed region + memset(result, 0, S); //initialize the integral to zero (0) + + //integrate the region between the specified bands and the closest indexed band + // this integrates the "tails" of the spectrum that lie outside the main indexed region + baseline_band(lb, rb, lp, rp, rab, cur2); //calculate the image for the right-most band in the integral + baseline_band(lb, rb, lp, rp, w[bi], cur); //calculate the image for the right-most indexed band for(unsigned long long j = 0; j < XY; j++){ result[j] += (T)((rab - w[bi]) * ((double)cur[j] + (double)cur2[j]) / 2.0); } @@ -479,13 +481,12 @@ public: result[j] += (T)((w[ai] - lab) * ((double)cur[j] + (double)cur2[j]) / 2.0); } - //calculate the area + //integrate the main indexed region ai++; - for(unsigned long long i = ai; i <= bi ;i++) + for(unsigned long long i = ai; i <= bi ;i++) //for each band in the integral { - baseline_band(lb, rb, lp, rp, w[ai], cur2); - for(unsigned long long j = 0; j < XY; j++) - { + baseline_band(lb, rb, lp, rp, w[ai], cur2); //get the baselined band + for(unsigned long long j = 0; j < XY; j++){ //for each pixel in that band result[j] += (T)((w[ai] - w[ai-1]) * ((double)cur[j] + (double)cur2[j]) / 2.0); } std::swap(cur,cur2); //swap the band pointers @@ -507,20 +508,22 @@ public: /// @param rb2 is the label value for the right baseline point for the second peak (denominator) /// @param pos2 is the label value for the second peak (denominator) position /// @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){ + 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)); + size_t XYbytes = X() * Y() * sizeof(T); //calculate the size of the band image (in bytes) + + T* p1 = (T*)malloc(XYbytes); //allocate space for both bands in the ratio + T* p2 = (T*)malloc(XYbytes); + memset(result, 0, XYbytes); //initialize the ratio to zero //get the two peak band height(lb1, rb1, pos1, p1); height(lb2, rb2, pos2, p2); //calculate the ratio in result for(unsigned long long i = 0; i < X() * Y(); i++){ - if(p1[i] == 0 && p2[i] ==0) - result[i] = 1; - else + if(mask == NULL || mask[i]){ result[i] = p1[i] / p2[i]; + } } free(p1); @@ -537,20 +540,20 @@ public: /// @param rb2 is the label value for the right baseline point for the second peak (denominator) /// @param pos2 is the label value for the second peak (denominator) position /// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size - bool pa_to_ph(double lb1, double rb1, double lab1, double rab1, - double lb2, double rb2, double pos, T* result){ + 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)); + size_t bytes = X() * Y() * sizeof(T); + T* p1 = (T*)malloc(bytes); //allocate space for both ratio components + T* p2 = (T*)malloc(bytes); + memset(result, 0, bytes); //initialize the ratio to zero //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 long long i = 0; i < X() * Y(); i++){ - if(p1[i] == 0 && p2[i] ==0) - result[i] = 1; - else + if(mask == NULL || mask[i]) result[i] = p1[i] / p2[i]; } @@ -570,21 +573,21 @@ public: /// @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 result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size - bool pa_to_pa(double lb1, double rb1, double lab1, double rab1, - double lb2, double rb2, double lab2, double rab2, T* result){ + 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)); + size_t bytes = X() * Y() * sizeof(T); + T* p1 = (T*)malloc(bytes); //allocate space for each of the operands + T* p2 = (T*)malloc(bytes); + memset(result, 0, bytes); //initialize the ratio result to zero (0) //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 long long i = 0; i < X() * Y(); i++){ - if(p1[i] == 0 && p2[i] ==0) - result[i] = 1; - else - result[i] = p1[i] / p2[i]; + if(mask == NULL || mask[i]) //if the pixel is masked + result[i] = p1[i] / p2[i]; //calculate the ratio } free(p1); @@ -592,7 +595,7 @@ public: return true; } - /// Compute the definite integral of a baseline corrected peak. + /// Compute the definite integral of a baseline corrected peak weighted by the corresponding wavelength /// @param lb is the label value for the left baseline point /// @param rb is the label value for the right baseline point @@ -613,8 +616,6 @@ public: cur = (T*) malloc(S); cur2 = (T*) malloc(S); - memset(result, (char)0, S); - //find the wavelenght position in the whole band unsigned long long n = w.size(); unsigned long long ai = 0; //left bound position @@ -642,6 +643,8 @@ public: band(lp, lb); band(rp, rb); + memset(result, (char)0, S); //initialize the integral to zero (0) + //calculate the beginning and the ending part baseline_band(lb, rb, lp, rp, rab, cur2); //ending part baseline_band(lb, rb, lp, rp, w[bi], cur); @@ -656,12 +659,11 @@ public: //calculate f(x) times x ai++; - for(unsigned long long i = ai; i <= bi ;i++) - { + for(unsigned long long i = ai; i <= bi ;i++){ 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); + for(unsigned long long j = 0; j < XY; j++){ + T v = (T)((w[ai] - w[ai-1]) * (w[ai] + w[ai-1]) * ((double)cur[j] + (double)cur2[j]) / 4.0); + result[j] += v; } std::swap(cur,cur2); //swap the band pointers } @@ -674,25 +676,27 @@ public: } /// Compute the centroid of a baseline corrected peak. + /// Note that the values for the centroid can be outside of [lab, rab] if the spectrum goes negative. /// @param lb is the label value for the left baseline point /// @param rb is the label value for the right baseline point /// @param lab is the label for the start of the peak /// @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(X() * Y() * sizeof(T)); + bool centroid(T* result, double lb, double rb, double lab, double rab, unsigned char* mask = NULL){ + size_t bytes = X() * Y() * sizeof(T); //calculate the number of bytes in a band image + T* p1 = (T*)malloc(X() * Y() * sizeof(T)); //allocate space for both operands T* p2 = (T*)malloc(X() * Y() * sizeof(T)); + memset(result, 0, bytes); //initialize the ratio result to zero (0) //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 long long i = 0; i < X() * Y(); i++){ - if(p1[i] == 0 && p2[i] ==0) - result[i] = 1; - else + if(mask == NULL || mask[i]){ result[i] = p1[i] / p2[i]; + } } free(p1); @@ -707,16 +711,16 @@ public: /// @param mask_band is the band used to specify the mask /// @param threshold is the threshold used to determine if the mask value is true or false /// @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, bool PROGRESS = false){ + bool build_mask(unsigned char* mask, double mask_band, double threshold, bool PROGRESS = false){ T* temp = (T*)malloc(X() * Y() * sizeof(T)); //allocate memory for the certain band band(temp, mask_band); for (unsigned long long i = 0; i < X() * Y(); i++) { if (temp[i] < threshold) - p[i] = 0; + mask[i] = 0; else - p[i] = 255; + mask[i] = 255; if(PROGRESS) progress = (double) (i+1) / (X() * Y()) * 100; } diff --git a/stim/envi/envi.h b/stim/envi/envi.h index 7664e8d..6d949e1 100644 --- a/stim/envi/envi.h +++ b/stim/envi/envi.h @@ -475,31 +475,31 @@ public: /// @param mask_band is the label for the band that will be used to build the mask /// @param threshold is a value selected such that all band values greater than threshold will have a mask value of 'true' /// @param p is memory of size X*Y that will store the resulting mask - bool build_mask(double mask_band, double threshold, unsigned char* p = NULL, bool PROGRESS = false) { + bool build_mask(unsigned char* mask, double mask_band, double threshold, 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)->build_mask(mask_band, threshold, p, PROGRESS); + return ((bsq*)file)->build_mask(mask, mask_band, threshold, PROGRESS); else if(header.data_type == envi_header::float64) - return ((bsq*)file)->build_mask(mask_band, threshold, p, PROGRESS); + return ((bsq*)file)->build_mask(mask, mask_band, threshold, PROGRESS); else std::cout<<"ERROR: unidentified data type"<*)file)->build_mask(mask_band, threshold, p, PROGRESS); + return ((bil*)file)->build_mask(mask, mask_band, threshold, PROGRESS); else if(header.data_type == envi_header::float64) - return ((bil*)file)->build_mask(mask_band, threshold, p, PROGRESS); + return ((bil*)file)->build_mask(mask, mask_band, threshold, PROGRESS); else std::cout<<"ERROR: unidentified data type"<*)file)->build_mask(mask_band, threshold, p, PROGRESS); + return ((bip*)file)->build_mask(mask, mask_band, threshold, PROGRESS); else if(header.data_type == envi_header::float64) - return ((bip*)file)->build_mask(mask_band, threshold, p, PROGRESS); + return ((bip*)file)->build_mask(mask, mask_band, threshold, PROGRESS); else std::cout<<"ERROR: unidentified data type"<*)file)->mask_finite(mask, PROGRESS); + else if(header.data_type == envi_header::float64) + ((bsq*)file)->mask_finite(mask, PROGRESS); + else + std::cout<<"ERROR: unidentified data type"<*)file)->mask_finite(mask, PROGRESS); + else if(header.data_type == envi_header::float64) + ((bil*)file)->mask_finite(mask, PROGRESS); + else + std::cout<<"ERROR: unidentified data type"<*)file)->mask_finite(mask, PROGRESS); + else if(header.data_type == envi_header::float64) + ((bip*)file)->mask_finite(mask, PROGRESS); + else + std::cout<<"ERROR: unidentified data type"<*)file)->sift((float*)matrix, p); else if (header.data_type == envi_header::float64) return ((bsq*)file)->sift((double*)matrix, p); - else + else{ std::cout << "ERROR: unidentified data type" << std::endl; + exit(1); + } } if (header.interleave == envi_header::BIP){ - std::cout<<"Memory sifting not supported for BIP files"<*)file)->sift((float*)matrix, p); + else if (header.data_type == envi_header::float64) + return ((bip*)file)->sift((double*)matrix, p); + else{ + std::cout << "ERROR: unidentified data type" << std::endl; + exit(1); + } } if (header.interleave == envi_header::BIL){ - std::cout<<"Memory sifting not supported for BIL files"<*)file)->sift((float*)matrix, p); + else if (header.data_type == envi_header::float64) + return ((bil*)file)->sift((double*)matrix, p); + else{ + std::cout << "ERROR: unidentified data type" << std::endl; + exit(1); + } } return false; @@ -684,30 +728,30 @@ public: /// @param rb2 is the label value for the right baseline point for the second peak (denominator) /// @param pos2 is the label value for the second peak (denominator) position /// @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, void * result){ + bool ph_to_ph(void * result, double lb1, double rb1, double pos1, double lb2, double rb2, double pos2, unsigned char* mask){ if(header.interleave == envi_header::BSQ){ //if the infile is bsq file if(header.data_type ==envi_header::float32) - return ((bsq*)file)->ph_to_ph(lb1, rb1, pos1, lb2, rb2, pos2, (float*)result); + return ((bsq*)file)->ph_to_ph((float*)result, lb1, rb1, pos1, lb2, rb2, pos2, mask); else if(header.data_type == envi_header::float64) - return ((bsq*)file)->ph_to_ph(lb1, rb1, pos1, lb2, rb2, pos2, (double*)result); + return ((bsq*)file)->ph_to_ph((double*)result, lb1, rb1, pos1, lb2, rb2, pos2, mask); else std::cout<<"ERROR: unidentified data type"<*)file)->ph_to_ph(lb1, rb1, pos1, lb2, rb2, pos2, (float*)result); + return ((bil*)file)->ph_to_ph((float*)result, lb1, rb1, pos1, lb2, rb2, pos2, mask); else if(header.data_type == envi_header::float64) - return ((bil*)file)->ph_to_ph(lb1, rb1, pos1, lb2, rb2, pos2, (double*)result); + return ((bil*)file)->ph_to_ph((double*)result, lb1, rb1, pos1, lb2, rb2, pos2, mask); else std::cout<<"ERROR: unidentified data type"<*)file)->ph_to_ph(lb1, rb1, pos1, lb2, rb2, pos2, (float*)result); + return ((bip*)file)->ph_to_ph((float*)result, lb1, rb1, pos1, lb2, rb2, pos2, mask); else if(header.data_type == envi_header::float64) - return ((bip*)file)->ph_to_ph(lb1, rb1, pos1, lb2, rb2, pos2, (double*)result); + return ((bip*)file)->ph_to_ph((double*)result, lb1, rb1, pos1, lb2, rb2, pos2, mask); else std::cout<<"ERROR: unidentified data type"<*)file)->pa_to_ph(lb1, rb1, lab1, rab1, lb2, rb2, pos, (float*)result); + return ((bsq*)file)->pa_to_ph((float*)result, lb1, rb1, lab1, rab1, lb2, rb2, pos, mask); else if(header.data_type == envi_header::float64) - return ((bsq*)file)->pa_to_ph(lb1, rb1, lab1, rab1, lb2, rb2, pos, (double*)result); + return ((bsq*)file)->pa_to_ph((double*)result, lb1, rb1, lab1, rab1, lb2, rb2, pos, mask); else std::cout<<"ERROR: unidentified data type"<*)file)->pa_to_ph(lb1, rb1, lab1, rab1, lb2, rb2, pos, (float*)result); + return ((bil*)file)->pa_to_ph((float*)result, lb1, rb1, lab1, rab1, lb2, rb2, pos, mask); else if(header.data_type == envi_header::float64) - return ((bil*)file)->pa_to_ph(lb1, rb1, lab1, rab1, lb2, rb2, pos, (double*)result); + return ((bil*)file)->pa_to_ph((double*)result, lb1, rb1, lab1, rab1, lb2, rb2, pos, mask); else std::cout<<"ERROR: unidentified data type"<*)file)->pa_to_ph(lb1, rb1, lab1, rab1, lb2, rb2, pos, (float*)result); + return ((bip*)file)->pa_to_ph((float*)result, lb1, rb1, lab1, rab1, lb2, rb2, pos, mask); else if(header.data_type == envi_header::float64) - return ((bip*)file)->pa_to_ph(lb1, rb1, lab1, rab1, lb2, rb2, pos, (double*)result); + return ((bip*)file)->pa_to_ph((double*)result, lb1, rb1, lab1, rab1, lb2, rb2, pos, mask); else std::cout<<"ERROR: unidentified data type"<*)file)->pa_to_pa(lb1, rb1, lab1, rab1, lb2, rb2, lab2, rab2, (float*)result); + return ((bsq*)file)->pa_to_pa((float*)result, lb1, rb1, lab1, rab1, lb2, rb2, lab2, rab2, mask); else if(header.data_type == envi_header::float64) - return ((bsq*)file)->pa_to_pa(lb1, rb1, lab1, rab1, lb2, rb2, lab2, rab2, (double*)result); + return ((bsq*)file)->pa_to_pa((double*)result, lb1, rb1, lab1, rab1, lb2, rb2, lab2, rab2, mask); else std::cout<<"ERROR: unidentified data type"<*)file)->pa_to_pa(lb1, rb1, lab1, rab1, lb2, rb2, lab2, rab2, (float*)result); + return ((bil*)file)->pa_to_pa((float*)result, lb1, rb1, lab1, rab1, lb2, rb2, lab2, rab2, mask); else if(header.data_type == envi_header::float64) - return ((bil*)file)->pa_to_pa(lb1, rb1, lab1, rab1, lb2, rb2, lab2, rab2, (double*)result); + return ((bil*)file)->pa_to_pa((double*)result, lb1, rb1, lab1, rab1, lb2, rb2, lab2, rab2, mask); else std::cout<<"ERROR: unidentified data type"<*)file)->pa_to_pa(lb1, rb1, lab1, rab1, lb2, rb2, lab2, rab2, (float*)result); + return ((bip*)file)->pa_to_pa((float*)result, lb1, rb1, lab1, rab1, lb2, rb2, lab2, rab2, mask); else if(header.data_type == envi_header::float64) - return ((bip*)file)->pa_to_pa(lb1, rb1, lab1, rab1, lb2, rb2, lab2, rab2, (double*)result); + return ((bip*)file)->pa_to_pa((double*)result, lb1, rb1, lab1, rab1, lb2, rb2, lab2, rab2, mask); else std::cout<<"ERROR: unidentified data type"<*)file)->cpoint(lb1, rb1, lab1, rab1, (float*)result); + return ((bsq*)file)->centroid((float*)result, lb1, rb1, lab1, rab1, mask); else if(header.data_type == envi_header::float64) - return ((bsq*)file)->cpoint(lb1, rb1, lab1, rab1, (double*)result); + return ((bsq*)file)->centroid((double*)result, lb1, rb1, lab1, rab1, mask); else std::cout<<"ERROR: unidentified data type"<*)file)->cpoint(lb1, rb1, lab1, rab1, (float*)result); + return ((bil*)file)->centroid((float*)result, lb1, rb1, lab1, rab1, mask); else if(header.data_type == envi_header::float64) - return ((bil*)file)->cpoint(lb1, rb1, lab1, rab1, (double*)result); + return ((bil*)file)->centroid((double*)result, lb1, rb1, lab1, rab1, mask); else std::cout<<"ERROR: unidentified data type"<*)file)->cpoint(lb1, rb1, lab1, rab1,(float*)result); + return ((bip*)file)->centroid((float*)result, lb1, rb1, lab1, rab1, mask); else if(header.data_type == envi_header::float64) - return ((bip*)file)->cpoint(lb1, rb1, lab1, rab1, (double*)result); + return ((bip*)file)->centroid((double*)result, lb1, rb1, lab1, rab1, mask); else std::cout<<"ERROR: unidentified data type"< #include +#ifdef _WIN32 + #include +#else + #include +#endif + namespace stim{ /** @@ -113,6 +119,46 @@ protected: return c[2] * R[0] * R[1] + c[1] * R[0] + c[0]; //calculate and return the index (trust me this works) } + + /// Returns the 3D coordinates of a specified index + void xyb(size_t idx, size_t& x, size_t& y, size_t& b){ + + size_t c[3]; + + c[2] = idx / (R[0] * R[1]); + c[1] = (idx - c[2] * R[0] * R[1]) / R[0]; + c[0] = idx - c[2] * R[0] * R[1] - c[1] * R[0]; + + x = c[O[0]]; + y = c[O[1]]; + b = c[O[2]]; + } + +public: + /// Get a mask that has all pixels with inf or NaN values masked out (false) + void mask_finite(unsigned char* mask, bool PROGRESS = false){ + 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 + if(!std::isfinite(page[i])){ //C++11 implementation +#endif + size_t x, y, b; + xyb(p * R[0] * R[1] + i, x, y, b); //find the 3D coordinates of the value + mask[ y * X() + x ] = 0; //mask the pixel (it's not bad) + } + } + if(PROGRESS) progress = (double)(p + 1) / (double)R[2] * 100; + } + } + }; } //end namespace STIM diff --git a/stim/parser/arguments.h b/stim/parser/arguments.h index 3c64e23..41231f6 100644 --- a/stim/parser/arguments.h +++ b/stim/parser/arguments.h @@ -97,7 +97,7 @@ namespace stim{ } //return the value of a floating point option - float as_float(size_t n = 0) + double as_float(size_t n = 0) { if(!flag) { -- libgit2 0.21.4