codebase / stimlib

Commit 11f177d5e4f7a3a6c4938ba837e59a120e16fe21

Authored by heziqi
1 parent dcdf1276

Ziqi added functions in biq.h

Showing 1 changed file with 291 additions and 102 deletions   Show diff stats
Inline Side-by-side
envi/bip.h
  Show/Hide comments   View file @11f177d
... ... @@ -40,7 +40,14 @@ public:
40 40 return false;
41 41 }
42 42  
43   - getSlice(p, 0, page);
  43 + //binary::getSlice(p, 0, page); //I met some problems when I called this function?????
  44 + file.seekg(page * sizeof(T), std::ios::beg);
  45 + for (unsigned i = 0; i < header.samples * header.lines; i++)
  46 + {
  47 + file.read((char *)(p + i), sizeof(T));
  48 + file.seekg( (header.bands - 1) * sizeof(T), std::ios::cur);
  49 + }
  50 +
44 51 return true;
45 52 }
46 53  
... ... @@ -89,15 +96,156 @@ public:
89 96 free(p2);
90 97 return true;
91 98 }
  99 + //get YZ line from the a Y slice, Y slice data should be already IN the MEMORY
  100 + bool getYZ(T* p, T* c, double wavelength)
  101 + {
  102 + unsigned int X = header.samples; //calculate the number of pixels in a sample
  103 + unsigned int B = header.bands;
  104 +
  105 + unsigned page=0; //samples around the wavelength
  106 + T * p1;
  107 + T * p2;
  108 +
  109 + //get the bands numbers around the wavelength
  110 +
  111 + //if wavelength is smaller than the first one in header file
  112 + if ( header.wavelength[page] > wavelength ){
  113 + for(unsigned j = 0; j < header.samples; j++)
  114 + {
  115 + p[j] = c[j * B];
  116 + }
  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 + for(unsigned j = 0; j < header.samples; j++)
  126 + {
  127 + p[j] = c[(j + 1) * B - 1];
  128 + }
  129 + return true;
  130 + }
  131 + }
  132 + if ( wavelength < header.wavelength[page] )
  133 + {
  134 + p1=(T*)malloc( X * sizeof(T)); //memory allocation
  135 + p2=(T*)malloc( X * sizeof(T));
  136 + //band_index(p1, page - 1);
  137 + for(unsigned j = 0; j < X; j++)
  138 + {
  139 + p1[j] = c[j * B + page - 1];
  140 + }
  141 + //band_index(p2, page );
  142 + for(unsigned j = 0; j < X; j++)
  143 + {
  144 + p2[j] = c[j * B + page];
  145 + }
  146 +
  147 + for(unsigned i=0; i < X; i++){
  148 + double r = (double) (wavelength - header.wavelength[page-1]) / (double) (header.wavelength[page] - header.wavelength[page-1]);
  149 + p[i] = (p2[i] - p1[i]) * r + p1[i];
  150 + }
  151 + }
  152 + else //if the wavelength is equal to a wavelength in header file
  153 + {
  154 + //band_index(p, page);
  155 + for(unsigned j = 0; j < X; j++)
  156 + {
  157 + p[j] = c[j * B + page];
  158 + }
  159 + }
  160 +
  161 + return true;
  162 + }
  163 + //given a y and a wavelength, return the y-band data
  164 + //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"
  165 + bool get_x_wavelength(T * p, unsigned y, double wavelength)
  166 + {
  167 + unsigned int X = header.samples; //calculate the number of pixels in a sample
  168 +
  169 + unsigned page=0; //samples around the wavelength
  170 + T * p1;
  171 + T * p2;
  172 +
  173 + //get the bands numbers around the wavelength
  174 +
  175 + //if wavelength is smaller than the first one in header file
  176 + if ( header.wavelength[page] > wavelength ){
  177 + file.seekg( y * header.bands * header.samples * sizeof(T), std::ios::beg);
  178 + for(unsigned j = 0; j < header.samples; j++)
  179 + {
  180 + file.read((char *)(p + j), sizeof(T));
  181 + file.seekg((header.bands - 1) * sizeof(T), std::ios::cur);
  182 + }
  183 + return true;
  184 + }
  185 +
  186 + while( header.wavelength[page] < wavelength )
  187 + {
  188 + page++;
  189 + //if wavelength is larger than the last wavelength in header file
  190 + if (page == header.bands) {
  191 + file.seekg( (y * header.bands * header.samples + header.bands - 1)* sizeof(T), std::ios::beg);
  192 + for(unsigned j = 0; j < header.samples; j++)
  193 + {
  194 + file.read((char *)(p + j), sizeof(T));
  195 + file.seekg((header.bands - 1) * sizeof(T), std::ios::cur);
  196 + }
  197 + return true;
  198 + }
  199 + }
  200 + if ( wavelength < header.wavelength[page] )
  201 + {
  202 + p1=(T*)malloc( X * sizeof(T)); //memory allocation
  203 + p2=(T*)malloc( X * sizeof(T));
  204 + //band_index(p1, page - 1);
  205 + file.seekg( (y * header.bands * header.samples + page - 1)* sizeof(T), std::ios::beg);
  206 + for(unsigned j = 0; j < header.samples; j++)
  207 + {
  208 + file.read((char *)(p1 + j), sizeof(T));
  209 + file.seekg((header.bands - 1) * sizeof(T), std::ios::cur);
  210 + }
  211 + //band_index(p2, page );
  212 + file.seekg( (y * header.bands * header.samples + page)* sizeof(T), std::ios::beg);
  213 + for(unsigned j = 0; j < header.samples; j++)
  214 + {
  215 + file.read((char *)(p2 + j), sizeof(T));
  216 + file.seekg((header.bands - 1) * sizeof(T), std::ios::cur);
  217 + }
  218 +
  219 + for(unsigned i=0; i < header.samples; i++){
  220 + double r = (double) (wavelength - header.wavelength[page-1]) / (double) (header.wavelength[page] - header.wavelength[page-1]);
  221 + p[i] = (p2[i] - p1[i]) * r + p1[i];
  222 + }
  223 + }
  224 + else //if the wavelength is equal to a wavelength in header file
  225 + {
  226 + //band_index(p, page);
  227 + file.seekg( (y * header.bands * header.samples + page) * sizeof(T), std::ios::beg);
  228 + for(unsigned j = 0; j < header.samples; j++)
  229 + {
  230 + file.read((char *)(p + j), sizeof(T));
  231 + file.seekg((header.bands - 1) * sizeof(T), std::ios::cur);
  232 + }
  233 + }
  234 +
  235 + //free(p1);
  236 + //free(p2);
  237 + return true;
  238 + }
92 239  
93 240 //save one pixel of the BIP file into the memory, and return the pointer
  241 + //I have not used it so far
94 242 bool getSpectrum(T * p, unsigned x, unsigned y){
95 243  
96 244 unsigned int i;
97 245  
98 246 if ( x >= header.samples || y >= header.lines){ //make sure the sample and line number is right
99 247 std::cout<<"ERROR: sample or line out of range"<<std::endl;
100   - return false;
  248 + exit(1);
101 249 }
102 250  
103 251 file.seekg((x + y * header.samples) * header.bands * sizeof(T), std::ios::beg); //point to the certain sample and line
... ... @@ -106,158 +254,176 @@ public:
106 254  
107 255 return true;
108 256 }
109   - // not finished yet
110   - //(BIP) baseline correction
111   - bool baseline(std::string outname, std::vector<double> wls )
  257 +
  258 + //given a Y ,return a ZX slice
  259 + bool getY(T * p, unsigned y)
112 260 {
  261 + if ( y >= header.lines){ //make sure the line number is right
  262 + std::cout<<"ERROR: line out of range"<<std::endl;
  263 + exit(1);
  264 + }
  265 +
  266 + file.seekg(y * header.bands * header.samples * sizeof(T), std::ios::beg);
  267 + file.read((char *)p, sizeof(T) * header.bands * header.samples);
  268 +
  269 + return true;
  270 +
  271 + }
  272 +
  273 + //(BIP) baseline correction
  274 + //processing speed is too slow
  275 + bool baseline(std::string outname, std::vector<double> wls){
  276 +
113 277 unsigned N = wls.size(); //get the number of baseline points
114 278  
115 279 std::ofstream target(outname.c_str(), std::ios::binary); //open the target binary file
116 280 std::string headername = outname + ".hdr"; //the header file name
117 281  
118 282 //simplify image resolution
119   - unsigned int B = header.bands; //calculate the number of bands
120   - unsigned int L = B * sizeof(T);
  283 + unsigned int ZX = header.bands * header.samples; //calculate the number of points in a Y slice
  284 + unsigned int L = ZX * sizeof(T); //calculate the number of bytes of a Y slice
  285 + unsigned int B = header.bands;
  286 + unsigned int X = header.samples;
121 287  
122   - T* c; //pointer to the current spectrum
123   - c = (T*)malloc(L);
124   -
125   - }
  288 + T* c; //pointer to the current Y slice
  289 + c = (T*)malloc(L); //memory allocation
126 290  
  291 + T* a; //pointer to the two YZ lines surrounding the current YZ line
  292 + T* b;
127 293  
128   - /*
129   - //(BSQ)baseline correction and save it into file
130   -
131   - bool baseline(std::string outname, std::vector<double> wls )
132   - {
133   - unsigned N = wls.size(); //get the number of baseline points
134   -
135   - std::ofstream target(outname.c_str(), std::ios::binary); //open the target binary file
136   - std::string headername = outname + ".hdr"; //the header file name
  294 + a = (T*)malloc(X * sizeof(T));
  295 + b = (T*)malloc(X * sizeof(T));
137 296  
138   - //simplify image resolution
139   - unsigned int B = header.bands; //calculate the number of bands
140   - unsigned int XY = header.samples * header.lines; //calculate the number of pixels in a band
141   - unsigned int S = XY * sizeof(T); //calculate the number of bytes in a band
142 297  
143 298 double ai, bi; //stores the two baseline points wavelength surrounding the current band
144 299 double ci; //stores the current band's wavelength
145   -// unsigned aii, bii; //stores the two baseline points number surrounding the current band
146 300 unsigned control=0;
147 301  
148   - T * a; //pointers to the high and low band images
149   - T * b;
150   - T * c; //pointer to the current image
151   -
152   - a = (T*)malloc( S ); //memory allocation
153   - b = (T*)malloc( S );
154   - c = (T*)malloc( S );
155   -
156 302 if (a == NULL || b == NULL || c == NULL){
157 303 std::cout<<"ERROR: error allocating memory";
158 304 exit(1);
159 305 }
160   -
161   -
162   - //initialize lownum, highnum, low, high
163   - ai=header.wavelength[0];
164   -
165   - //if no baseline point is specified at band 0,
  306 + // loop start correct every y slice
  307 +
  308 + for (unsigned k =0; k < header.lines; k++)
  309 + {
  310 + //get the current y slice
  311 + getY(c, k);
  312 +
  313 + //initialize lownum, highnum, low, high
  314 + ai = header.wavelength[0];
  315 + //if no baseline point is specified at band 0,
166 316 //set the baseline point at band 0 to 0
167   - if(wls[0] != header.wavelength[0]){
168   - bi = wls[control];
169   - memset(a, (char)0, S);
170   - }
171   - //else get the low band
172   - else{
173   - control += 1;
174   - getBand(a, ai);
175   - bi = wls[control];
176   - }
177   - //get the high band
178   - getBand(b, bi);
179   -
180   - //correct every band
181   - for(unsigned cii = 0; cii < B; cii++){
182   -
183   - //update baseline points, if necessary
184   - if( header.wavelength[cii] >= bi && cii != B - 1) {
185   - //if the high band is now on the last BL point?
186   - if (control != N-1) {
187   -
188   - control++; //increment the index
189   -
190   - std::swap(a, b); //swap the baseline band pointers
191   -
192   - ai = bi;
193   - bi = wls[control];
194   - getBand(b, bi);
  317 + if(wls[0] != header.wavelength[0]){
  318 + bi = wls[control];
  319 + memset(a, (char)0, X * sizeof(T) );
  320 + }
  321 + //else get the low band
  322 + else{
  323 + control++;
  324 + getYZ(a, c, ai);
  325 + bi = wls[control];
  326 + }
  327 + //get the high band
  328 + getYZ(b, c, bi);
  329 +
  330 + //correct every YZ line
  331 +
  332 + for(unsigned cii = 0; cii < B; cii++){
195 333  
  334 + //update baseline points, if necessary
  335 + if( header.wavelength[cii] >= bi && cii != B - 1) {
  336 + //if the high band is now on the last BL point?
  337 + if (control != N-1) {
  338 +
  339 + control++; //increment the index
  340 +
  341 + std::swap(a, b); //swap the baseline band pointers
  342 +
  343 + ai = bi;
  344 + bi = wls[control];
  345 + getYZ(b, c, bi);
  346 +
  347 + }
  348 + //if the last BL point on the last band of the file?
  349 + else if ( wls[control] < header.wavelength[B - 1]) {
  350 +
  351 + std::swap(a, b); //swap the baseline band pointers
  352 +
  353 + memset(b, (char)0, X * sizeof(T) ); //clear the high band
  354 +
  355 + ai = bi;
  356 + bi = header.wavelength[B - 1];
  357 + }
196 358 }
197   - //if the last BL point on the last band of the file?
198   - else if ( wls[control] < header.wavelength[B - 1]) {
199 359  
200   - std::swap(a, b); //swap the baseline band pointers
201   -
202   - memset(b, (char)0, S); //clear the high band
203   -
204   - ai = bi;
205   - bi = header.wavelength[B - 1];
  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);
206 367 }
207   - }
208   -
209   - //get the current band
210   - band_index(c, cii);
211   - ci = header.wavelength[cii];
212 368  
213   - //perform the baseline correction
214   - for(unsigned i=0; i < XY; i++){
215   - double r = (double) (ci - ai) / (double) (bi - ai);
216   - c[i] =(float) ( c[i] - (b[i] - a[i]) * r - a[i] );
217   - }
  369 + ci = header.wavelength[cii];
  370 +
  371 + //perform the baseline correction
  372 + for(unsigned i=0; i < X; i++)
  373 + {
  374 + double r = (double) (ci - ai) / (double) (bi - ai);
  375 + c[i * B + cii] =(float) ( c[i * B + cii] - (b[i] - a[i]) * r - a[i] );
  376 + }
  377 +
  378 + }//loop for YZ line end
218 379  
219   - target.write(reinterpret_cast<const char*>(c), S); //write the corrected data into destination
  380 + target.write(reinterpret_cast<const char*>(c), L); //write the corrected data into destination
  381 + }//loop for Y slice end
220 382  
221   - }
222   -
223 383 header.save(headername); //save the new header file
224 384  
225 385 free(a);
226 386 free(b);
227 387 free(c);
228 388 target.close();
229   - return true;
  389 + return true;
  390 +
230 391 }
231   -*/
232   - // normalize the BSQ file
  392 +
  393 + // normalize the BIP file
233 394 bool normalize(std::string outname, double band)
234 395 {
235 396 unsigned int B = header.bands; //calculate the number of bands
  397 + unsigned int Y = header.lines;
  398 + unsigned int X = header.samples;
  399 + unsigned int ZX = header.bands * header.samples;
236 400 unsigned int XY = header.samples * header.lines; //calculate the number of pixels in a band
237   - unsigned int S = XY * sizeof(T); //calculate the number of bytes in a band
  401 + unsigned int S = XY * sizeof(T); //calculate the number of bytes in a band
  402 + unsigned int L = ZX * sizeof(T);
238 403  
239 404 std::ofstream target(outname.c_str(), std::ios::binary); //open the target binary file
240 405 std::string headername = outname + ".hdr"; //the header file name
241 406  
242   - T * b; //pointers to the certain wavelength band
243   - T * c; //pointer to the current image
  407 + T * c; //pointer to the current ZX slice
  408 + T * b; //pointer to the standard band
244 409  
245 410 b = (T*)malloc( S ); //memory allocation
246   - c = (T*)malloc( S );
  411 + c = (T*)malloc( L );
247 412  
248 413 getBand(b, band); //get the certain band into memory
249 414  
250   - for(unsigned j = 0; j < B; j++)
  415 + for(unsigned j = 0; j < Y; j++)
251 416 {
252   - band_index(c, j); //get the current band into memory
253   - for(unsigned i = 0; i < XY; i++)
  417 + getY(c, j);
  418 + for(unsigned i = 0; i < X; i++)
254 419 {
255   - c[i] = c[i] / b[i];
  420 + for(unsigned m = 0; m < B; m++)
  421 + {
  422 + c[m + i * B] = c[m + i * B] / b[i + j * X];
  423 + }
256 424 }
257   - target.write(reinterpret_cast<const char*>(c), S); //write normalized data into destination
258   - std::cout<<j<<std::endl;
  425 + target.write(reinterpret_cast<const char*>(c), L); //write normalized data into destination
259 426 }
260   -
261 427 header.save(headername); //save the new header file
262 428  
263 429 free(b);
... ... @@ -266,6 +432,29 @@ public:
266 432 return true;
267 433 }
268 434  
  435 + //convert BIP file to BSQ file and save it
  436 + bool bsq(std::string outname)
  437 + {
  438 + unsigned int S = header.samples * header.lines * sizeof(T); //calculate the number of bytes in a band
  439 +
  440 + std::ofstream target(outname.c_str(), std::ios::binary);
  441 + std::string headername = outname + ".hdr";
  442 +
  443 + T * p; //pointer to the current band
  444 + p = (T*)malloc(S);
  445 +
  446 + for ( unsigned i = 0; i < header.bands; i++)
  447 + {
  448 + band_index(p, i);
  449 + target.write(reinterpret_cast<const char*>(p), S); //write a band data into target file
  450 + }
  451 + header.interleave = rts::envi::interleaveType::BSQ; //change the type of file in header file
  452 + header.save(headername);
  453 +
  454 + free(p);
  455 + target.close();
  456 + return true;
  457 + }
269 458  
270 459 };
271 460 }
272 461 \ No newline at end of file
... ...