codebase / stimlib

Commit c25e7d0d7c5ce33f99b99126de0fd426be193ec1

Authored by heziqi
1 parent 6dbb23ce

speed of bip.baseline fixed

Showing 3 changed files with 403 additions and 13 deletions   Show diff stats
Inline Side-by-side
envi/bil.h
  Show/Hide comments   View file @c25e7d0
  1 +#include "../envi/envi.h"
  2 +#include "../envi/binary.h"
  3 +#include <cstring>
  4 +#include <utility>
  5 +
  6 +namespace rts{
  7 +
  8 +template <typename T>
  9 +
  10 +class bil: public binary<T> {
  11 +
  12 +protected:
  13 +
  14 + envi header;
  15 +
  16 +public:
  17 +
  18 + using binary<T>::open;
  19 + using binary<T>::file;
  20 + using binary<T>::getSlice;
  21 +
  22 + //open a file, given the file and its header's names
  23 + bool open(std::string filename, std::string headername){
  24 +
  25 + if (header.load(headername)==false){
  26 + std::cout<<"ERROR: unable to load header file: "<<headername<<std::endl;
  27 + return false;
  28 + }
  29 +
  30 + open(filename, vec<unsigned int>(header.samples, header.lines, header.bands), header.header_offset);
  31 + return true;
  32 +
  33 + }
  34 +
  35 + //save one band of the file into the memory, and return the pointer
  36 + bool band_index( T * p, unsigned int page){
  37 +
  38 + unsigned int L = header.samples * sizeof(T); //caculate the number of bytes in a sample line
  39 + unsigned int jump = header.samples * (header.bands - 1) * sizeof(T);
  40 +
  41 + if (page >= header.bands){ //make sure the bank number is right
  42 + std::cout<<"ERROR: page out of range"<<std::endl;
  43 + return false;
  44 + }
  45 +
  46 + file.seekg(header.samples * page * sizeof(T), std::ios::beg);
  47 + for (unsigned i = 0; i < header.lines; i++)
  48 + {
  49 + file.read((char *)(p + i * header.samples), L);
  50 + file.seekg( jump, std::ios::cur);
  51 + }
  52 +
  53 + return true;
  54 + }
  55 +
  56 + bool getBand( T * p, double wavelength){
  57 +
  58 + unsigned int XY = header.samples * header.lines; //calculate the number of pixels in a band
  59 + unsigned int S = XY * sizeof(T); //calculate the number of bytes of a band
  60 +
  61 + unsigned page=0; //bands around the wavelength
  62 + T * p1;
  63 + T * p2;
  64 +
  65 + //get the bands numbers around the wavelength
  66 +
  67 + //if wavelength is smaller than the first one in header file
  68 + if ( header.wavelength[page] > wavelength ){
  69 + band_index(p, page);
  70 + return true;
  71 + }
  72 +
  73 + while( header.wavelength[page] < wavelength )
  74 + {
  75 + page++;
  76 + //if wavelength is larger than the last wavelength in header file
  77 + if (page == header.bands) {
  78 + band_index(p, header.bands-1);
  79 + return true;
  80 + }
  81 + }
  82 + if ( wavelength < header.wavelength[page] )
  83 + {
  84 + p1=(T*)malloc(S); //memory allocation
  85 + p2=(T*)malloc(S);
  86 + band_index(p1, page - 1);
  87 + band_index(p2, page );
  88 + for(unsigned i=0; i < XY; i++){
  89 + double r = (double) (wavelength - header.wavelength[page-1]) / (double) (header.wavelength[page] - header.wavelength[page-1]);
  90 + p[i] = (p2[i] - p1[i]) * r + p1[i];
  91 + }
  92 + }
  93 + else //if the wavelength is equal to a wavelength in header file
  94 + {
  95 + band_index(p, page);
  96 + }
  97 +
  98 + return true;
  99 + }
  100 +
  101 + //get YZ line from the a Y slice, Y slice data should be already IN the MEMORY
  102 + bool getYZ(T* p, T* c, double wavelength)
  103 + {
  104 + unsigned int X = header.samples; //calculate the number of pixels in a sample
  105 + unsigned int B = header.bands;
  106 + unsigned int L = X * sizeof(T);
  107 +
  108 + unsigned page=0; //samples around the wavelength
  109 + T * p1;
  110 + T * p2;
  111 +
  112 + //get the bands numbers around the wavelength
  113 +
  114 + //if wavelength is smaller than the first one in header file
  115 + if ( header.wavelength[page] > wavelength ){
  116 + memcpy(p, c, L);
  117 + return true;
  118 + }
  119 +
  120 + while( header.wavelength[page] < wavelength )
  121 + {
  122 + page++;
  123 + //if wavelength is larger than the last wavelength in header file
  124 + if (page == B) {
  125 + memcpy(p, c + (B - 1) * X, L);
  126 + return true;
  127 + }
  128 + }
  129 + if ( wavelength < header.wavelength[page] )
  130 + {
  131 + p1=(T*)malloc( L ); //memory allocation
  132 + p2=(T*)malloc( L );
  133 +
  134 + memcpy(p1, c + (page - 1) * X, L);
  135 + memcpy(p2, c + page * X, L);
  136 +
  137 + for(unsigned i=0; i < X; i++){
  138 + double r = (double) (wavelength - header.wavelength[page-1]) / (double) (header.wavelength[page] - header.wavelength[page-1]);
  139 + p[i] = (p2[i] - p1[i]) * r + p1[i];
  140 + }
  141 + }
  142 + else //if the wavelength is equal to a wavelength in header file
  143 + memcpy(p, c + page * X, L);
  144 +
  145 + return true;
  146 + }
  147 +
  148 + //save one pixel of the BIP file into the memory, and return the pointer
  149 + bool getSpectrum(T * p, unsigned x, unsigned y){
  150 +
  151 + if ( x >= header.samples || y >= header.lines){ //make sure the sample and line number is right
  152 + std::cout<<"ERROR: sample or line out of range"<<std::endl;
  153 + exit(1);
  154 + }
  155 +
  156 + unsigned jump = (header.samples - 1) * sizeof(T);
  157 +
  158 + file.seekg((y * header.samples * header.bands + x) * sizeof(T), std::ios::beg);
  159 +
  160 + for(unsigned i = 0; i < header.bands; i++)
  161 + { //point to the certain sample and line
  162 + file.read((char *)(p + i), sizeof(T));
  163 + file.seekg(jump, std::ios::cur);
  164 + }
  165 +
  166 + return true;
  167 + }
  168 +
  169 + //given a Y ,return a XZ slice
  170 + bool getY(T * p, unsigned y)
  171 + {
  172 + if ( y >= header.lines){ //make sure the line number is right
  173 + std::cout<<"ERROR: line out of range"<<std::endl;
  174 + exit(1);
  175 + }
  176 +
  177 + file.seekg(y * header.bands * header.samples * sizeof(T), std::ios::beg);
  178 + file.read((char *)p, sizeof(T) * header.bands * header.samples);
  179 +
  180 + return true;
  181 +
  182 + }
  183 +
  184 + //(BIL) baseline correction
  185 + bool baseline(std::string outname, std::vector<double> wls){
  186 +
  187 + unsigned N = wls.size(); //get the number of baseline points
  188 +
  189 + std::ofstream target(outname.c_str(), std::ios::binary); //open the target binary file
  190 + std::string headername = outname + ".hdr"; //the header file name
  191 +
  192 + //simplify image resolution
  193 + unsigned int ZX = header.bands * header.samples; //calculate the number of points in a Y slice
  194 + unsigned int L = ZX * sizeof(T); //calculate the number of bytes of a Y slice
  195 + unsigned int B = header.bands;
  196 + unsigned int X = header.samples;
  197 +
  198 + T* c; //pointer to the current Y slice
  199 + c = (T*)malloc(L); //memory allocation
  200 +
  201 + T* a; //pointer to the two YZ lines surrounding the current YZ line
  202 + T* b;
  203 +
  204 + a = (T*)malloc(X * sizeof(T));
  205 + b = (T*)malloc(X * sizeof(T));
  206 +
  207 +
  208 + double ai, bi; //stores the two baseline points wavelength surrounding the current band
  209 + double ci; //stores the current band's wavelength
  210 + unsigned control=0;
  211 +
  212 + if (a == NULL || b == NULL || c == NULL){
  213 + std::cout<<"ERROR: error allocating memory";
  214 + exit(1);
  215 + }
  216 + // loop start correct every y slice
  217 +
  218 + for (unsigned k =0; k < header.lines; k++)
  219 + {
  220 + //get the current y slice
  221 + getY(c, k);
  222 +
  223 + //initialize lownum, highnum, low, high
  224 + ai = header.wavelength[0];
  225 + //if no baseline point is specified at band 0,
  226 + //set the baseline point at band 0 to 0
  227 + if(wls[0] != header.wavelength[0]){
  228 + bi = wls[control];
  229 + memset(a, (char)0, X * sizeof(T) );
  230 + }
  231 + //else get the low band
  232 + else{
  233 + control++;
  234 + getYZ(a, c, ai);
  235 + bi = wls[control];
  236 + }
  237 + //get the high band
  238 + getYZ(b, c, bi);
  239 +
  240 + //correct every YZ line
  241 +
  242 + for(unsigned cii = 0; cii < B; cii++){
  243 +
  244 + //update baseline points, if necessary
  245 + if( header.wavelength[cii] >= bi && cii != B - 1) {
  246 + //if the high band is now on the last BL point
  247 + if (control != N-1) {
  248 +
  249 + control++; //increment the index
  250 +
  251 + std::swap(a, b); //swap the baseline band pointers
  252 +
  253 + ai = bi;
  254 + bi = wls[control];
  255 + getYZ(b, c, bi);
  256 +
  257 + }
  258 + //if the last BL point on the last band of the file?
  259 + else if ( wls[control] < header.wavelength[B - 1]) {
  260 +
  261 + std::swap(a, b); //swap the baseline band pointers
  262 +
  263 + memset(b, (char)0, X * sizeof(T) ); //clear the high band
  264 +
  265 + ai = bi;
  266 + bi = header.wavelength[B - 1];
  267 + }
  268 + }
  269 +
  270 + //get the current YZ line
  271 + //band_index(c, cii);???????????????????????????????????
  272 + /*
  273 + file.seekg( (k * header.bands * X + cii)* sizeof(T), std::ios::beg);
  274 + for(unsigned j = 0; j < X; j++)
  275 + {
  276 + file.read((char *)(c + j), sizeof(T));
  277 + file.seekg((B - 1) * sizeof(T), std::ios::cur);
  278 + }*/
  279 +
  280 + ci = header.wavelength[cii];
  281 +
  282 + unsigned jump = cii * X;
  283 + //perform the baseline correction
  284 + for(unsigned i=0; i < X; i++)
  285 + {
  286 + double r = (double) (ci - ai) / (double) (bi - ai);
  287 + c[i + jump] =(float) ( c[i + jump] - (b[i] - a[i]) * r - a[i] );
  288 + }
  289 +
  290 + }//loop for YZ line end
  291 + std::cout<<k<<std::endl;
  292 + target.write(reinterpret_cast<const char*>(c), L); //write the corrected data into destination
  293 + }//loop for Y slice end
  294 +
  295 + header.save(headername); //save the new header file
  296 +
  297 + free(a);
  298 + free(b);
  299 + free(c);
  300 + target.close();
  301 + return true;
  302 +
  303 + }
  304 +
  305 + // normalize the BIP file
  306 + bool normalize(std::string outname, double band)
  307 + {
  308 + unsigned int B = header.bands; //calculate the number of bands
  309 + unsigned int Y = header.lines;
  310 + unsigned int X = header.samples;
  311 + unsigned int ZX = header.bands * header.samples;
  312 + unsigned int XY = header.samples * header.lines; //calculate the number of pixels in a band
  313 + unsigned int S = XY * sizeof(T); //calculate the number of bytes in a band
  314 + unsigned int L = ZX * sizeof(T);
  315 +
  316 + std::ofstream target(outname.c_str(), std::ios::binary); //open the target binary file
  317 + std::string headername = outname + ".hdr"; //the header file name
  318 +
  319 + T * c; //pointer to the current ZX slice
  320 + T * b; //pointer to the standard band
  321 +
  322 + b = (T*)malloc( S ); //memory allocation
  323 + c = (T*)malloc( L );
  324 +
  325 + getBand(b, band); //get the certain band into memory
  326 +
  327 + for(unsigned j = 0; j < Y; j++)
  328 + {
  329 + getY(c, j);
  330 + for(unsigned i = 0; i < X; i++)
  331 + {
  332 + for(unsigned m = 0; m < B; m++)
  333 + {
  334 + c[m + i * B] = c[m + i * B] / b[i + j * X];
  335 + }
  336 + }
  337 + target.write(reinterpret_cast<const char*>(c), L); //write normalized data into destination
  338 + }
  339 + header.save(headername); //save the new header file
  340 +
  341 + free(b);
  342 + free(c);
  343 + target.close();
  344 + return true;
  345 + }
  346 +
  347 + //convert BIP file to BSQ file and save it
  348 + bool bsq(std::string outname)
  349 + {
  350 + unsigned int S = header.samples * header.lines * sizeof(T); //calculate the number of bytes in a band
  351 +
  352 + std::ofstream target(outname.c_str(), std::ios::binary);
  353 + std::string headername = outname + ".hdr";
  354 +
  355 + T * p; //pointer to the current band
  356 + p = (T*)malloc(S);
  357 +
  358 + for ( unsigned i = 0; i < header.bands; i++)
  359 + {
  360 + band_index(p, i);
  361 + target.write(reinterpret_cast<const char*>(p), S); //write a band data into target file
  362 + }
  363 + header.interleave = rts::envi::interleaveType::BSQ; //change the type of file in header file
  364 + header.save(headername);
  365 +
  366 + free(p);
  367 + target.close();
  368 + return true;
  369 + }
  370 +
  371 + };
  372 +}
0 373 \ No newline at end of file
... ...
envi/bip.h
  Show/Hide comments   View file @c25e7d0
... ... @@ -241,8 +241,6 @@ public:
241 241 //I have not used it so far
242 242 bool getSpectrum(T * p, unsigned x, unsigned y){
243 243  
244   - unsigned int i;
245   -
246 244 if ( x >= header.samples || y >= header.lines){ //make sure the sample and line number is right
247 245 std::cout<<"ERROR: sample or line out of range"<<std::endl;
248 246 exit(1);
... ... @@ -271,7 +269,6 @@ public:
271 269 }
272 270  
273 271 //(BIP) baseline correction
274   - //processing speed is too slow
275 272 bool baseline(std::string outname, std::vector<double> wls){
276 273  
277 274 unsigned N = wls.size(); //get the number of baseline points
... ... @@ -357,15 +354,6 @@ public:
357 354 }
358 355 }
359 356  
360   - //get the current YZ line
361   - //band_index(c, cii);
362   - file.seekg( (k * header.bands * X + cii)* sizeof(T), std::ios::beg);
363   - for(unsigned j = 0; j < X; j++)
364   - {
365   - file.read((char *)(c + j), sizeof(T));
366   - file.seekg((B - 1) * sizeof(T), std::ios::cur);
367   - }
368   -
369 357 ci = header.wavelength[cii];
370 358  
371 359 //perform the baseline correction
... ... @@ -376,7 +364,6 @@ public:
376 364 }
377 365  
378 366 }//loop for YZ line end
379   -
380 367 target.write(reinterpret_cast<const char*>(c), L); //write the corrected data into destination
381 368 }//loop for Y slice end
382 369  
... ...
envi/bsq.h
  Show/Hide comments   View file @c25e7d0
... ... @@ -275,6 +275,37 @@ public:
275 275 target.close();
276 276 return true;
277 277 }
  278 +
  279 + //convert BSQ file to BIL file and save it
  280 + bool bil(std::string outname)
  281 + {
  282 + //simplify image resolution
  283 + unsigned int L = header.samples * header.bands * sizeof(T); //calculate the number of bytes of a ZX slice
  284 + unsigned int jump = (header.lines - 1) * header.samples * sizeof(T);
  285 +
  286 + std::ofstream target(outname.c_str(), std::ios::binary);
  287 + std::string headername = outname + ".hdr";
  288 +
  289 + T * p; //pointer to the current spectrum
  290 + p = (T*)malloc(L);
  291 +
  292 + for ( unsigned i = 0; i < header.lines; i++)
  293 + {
  294 + file.seekg(header.samples * i * sizeof(T), std::ios::beg);
  295 + for ( unsigned j = 0; j < header.bands; j++ )
  296 + {
  297 + file.read((char *)(p + j * header.samples), sizeof(T) * header.samples);
  298 + file.seekg(jump, std::ios::cur); //go to the next band
  299 + }
  300 + target.write(reinterpret_cast<const char*>(p), L); //write XZ slice data into target file
  301 + }
  302 + header.interleave = rts::envi::interleaveType::BIL; //change the type of file in header file
  303 + header.save(headername);
  304 +
  305 + free(p);
  306 + target.close();
  307 + return true;
  308 + }
278 309  
279 310  
280 311 };
... ...