Commit 70407ea9143be076b3c5e3b24305c87de98c6256
1 parent
f47168a2
Ziqi completed height2heigth area2height area2area ratio
Showing
4 changed files
with
311 additions
and
91 deletions
Show diff stats
envi/bil.h
@@ -403,17 +403,16 @@ public: | @@ -403,17 +403,16 @@ public: | ||
403 | } | 403 | } |
404 | //untested | 404 | //untested |
405 | //providing the left and the right bound wavelength, return baseline-corrected band height | 405 | //providing the left and the right bound wavelength, return baseline-corrected band height |
406 | - bool height(double lb, double rb, double bandwavelength){ | 406 | + bool height(double lb, double rb, double bandwavelength, T* result){ |
407 | 407 | ||
408 | T* lp; | 408 | T* lp; |
409 | T* rp; | 409 | T* rp; |
410 | - T* result; | ||
411 | - | ||
412 | - lp = (T*) malloc(sizeof(R[0] * R[1] * sizeof(T))); //memory allocation | ||
413 | - rp = (T*) malloc(sizeof(R[0] * R[1] * sizeof(T))); | ||
414 | - result = (T*) malloc(sizeof(R[0] * R[1] * sizeof(T))); | 410 | + unsigned XY = R[0] * R[1]; |
411 | + unsigned S = XY * sizeof(T); | ||
412 | + lp = (T*) malloc(S); //memory allocation | ||
413 | + rp = (T*) malloc(S); | ||
415 | 414 | ||
416 | - band(lp, lb); //get the data of the left bounbd and the right bound | 415 | + band(lp, lb); |
417 | band(rp, rb); | 416 | band(rp, rb); |
418 | 417 | ||
419 | baseline_band(lb, rb, lp, rp, bandwavelength, result); | 418 | baseline_band(lb, rb, lp, rp, bandwavelength, result); |
@@ -421,14 +420,14 @@ public: | @@ -421,14 +420,14 @@ public: | ||
421 | return true; | 420 | return true; |
422 | } | 421 | } |
423 | 422 | ||
423 | + | ||
424 | //calculate the area between two bound point(including baseline correction) | 424 | //calculate the area between two bound point(including baseline correction) |
425 | - bool area(double lb, double rb){ | 425 | + bool area(double lb, double rb, double lab, double rab, T* result){ |
426 | 426 | ||
427 | T* lp; //left band pointer | 427 | T* lp; //left band pointer |
428 | T* rp; //right band pointer | 428 | T* rp; //right band pointer |
429 | T* cur; //current band 1 | 429 | T* cur; //current band 1 |
430 | T* cur2; //current band 2 | 430 | T* cur2; //current band 2 |
431 | - T* result; //area result | ||
432 | 431 | ||
433 | unsigned XY = R[0] * R[1]; | 432 | unsigned XY = R[0] * R[1]; |
434 | unsigned S = XY * sizeof(T); | 433 | unsigned S = XY * sizeof(T); |
@@ -437,7 +436,6 @@ public: | @@ -437,7 +436,6 @@ public: | ||
437 | rp = (T*) malloc(S); | 436 | rp = (T*) malloc(S); |
438 | cur = (T*) malloc(S); | 437 | cur = (T*) malloc(S); |
439 | cur2 = (T*) malloc(S); | 438 | cur2 = (T*) malloc(S); |
440 | - result = (T*) malloc(S); | ||
441 | 439 | ||
442 | memset(result, (char)0, S); | 440 | memset(result, (char)0, S); |
443 | 441 | ||
@@ -460,10 +458,10 @@ public: | @@ -460,10 +458,10 @@ public: | ||
460 | } | 458 | } |
461 | 459 | ||
462 | //get the position of lb and rb | 460 | //get the position of lb and rb |
463 | - while (lb >= w[ai]){ | 461 | + while (lab >= w[ai]){ |
464 | ai++; | 462 | ai++; |
465 | } | 463 | } |
466 | - while (rb <= w[bi]){ | 464 | + while (rab <= w[bi]){ |
467 | bi--; | 465 | bi--; |
468 | } | 466 | } |
469 | 467 | ||
@@ -471,13 +469,15 @@ public: | @@ -471,13 +469,15 @@ public: | ||
471 | band(rp, rb); | 469 | band(rp, rb); |
472 | 470 | ||
473 | //calculate the beginning and the ending part | 471 | //calculate the beginning and the ending part |
474 | - baseline_band(lb, rb, lp, rp, w[bi], cur); //ending part | 472 | + baseline_band(lb, rb, lp, rp, rab, cur2); //ending part |
473 | + baseline_band(lb, rb, lp, rp, w[bi], cur); | ||
475 | for(unsigned j = 0; j < XY; j++){ | 474 | for(unsigned j = 0; j < XY; j++){ |
476 | - result[j] += (rb - w[bi]) * cur[j] / 2.0; | 475 | + result[j] += (rab - w[bi]) * (cur[j] + cur2[j]) / 2.0; |
477 | } | 476 | } |
478 | - baseline_band(lb, rb, lp, rp, w[ai], cur); //beginnning part | 477 | + baseline_band(lb, rb, lp, rp, lab, cur2); //beginnning part |
478 | + baseline_band(lb, rb, lp, rp, w[ai], cur); | ||
479 | for(unsigned j = 0; j < XY; j++){ | 479 | for(unsigned j = 0; j < XY; j++){ |
480 | - result[j] += (w[ai] - lb) * cur[j] / 2.0; | 480 | + result[j] += (w[ai] - lab) * (cur[j] + cur2[j]) / 2.0; |
481 | } | 481 | } |
482 | //calculate the area | 482 | //calculate the area |
483 | ai++; | 483 | ai++; |
@@ -490,20 +490,69 @@ public: | @@ -490,20 +490,69 @@ public: | ||
490 | } | 490 | } |
491 | std::swap(cur,cur2); //swap the band pointers | 491 | std::swap(cur,cur2); //swap the band pointers |
492 | } | 492 | } |
493 | - ///////////////for testing use | ||
494 | - std::ofstream text("area_result.txt",std::ios::out); | ||
495 | - for (unsigned i = 0; i < 640; i++) | ||
496 | - { | ||
497 | - for (unsigned j = 0; j < 384; j++) | ||
498 | - { | ||
499 | - text<<result[i*384+j]<<", "; | ||
500 | - } | ||
501 | - text<<std::endl; | ||
502 | - } | ||
503 | - /////////////// | 493 | + |
504 | return true; | 494 | return true; |
505 | } | 495 | } |
506 | 496 | ||
497 | + //peak height ratio | ||
498 | + bool ph_to_ph(double lb1, double rb1, double pos1, double lb2, double rb2, double pos2, T * result){ | ||
499 | + | ||
500 | + T* p1 = (T*)malloc(R[0] * R[1] * sizeof(T)); | ||
501 | + T* p2 = (T*)malloc(R[0] * R[1] * sizeof(T)); | ||
502 | + | ||
503 | + //get the two peak band | ||
504 | + height(lb1, rb1, pos1, p1); | ||
505 | + height(lb2, rb2, pos2, p2); | ||
506 | + //calculate the ratio in result | ||
507 | + for(unsigned i = 0; i < R[0] * R[1]; i++){ | ||
508 | + if(p1[i] == 0 && p2[i] ==0) | ||
509 | + result[i] = 1; | ||
510 | + else | ||
511 | + result[i] = p1[i] / p2[i]; | ||
512 | + } | ||
513 | + return true; | ||
514 | + } | ||
515 | + | ||
516 | + //peak are to peak height ratio | ||
517 | + bool pa_to_ph(double lb1, double rb1, double lab1, double rab1, | ||
518 | + double lb2, double rb2, double pos, T* result){ | ||
519 | + | ||
520 | + T* p1 = (T*)malloc(R[0] * R[1] * sizeof(T)); | ||
521 | + T* p2 = (T*)malloc(R[0] * R[1] * sizeof(T)); | ||
522 | + | ||
523 | + //get the area and the peak band | ||
524 | + area(lb1, rb1, lab1, rab1, p1); | ||
525 | + height(lb2, rb2, pos, p2); | ||
526 | + //calculate the ratio in result | ||
527 | + for(unsigned i = 0; i < R[0] * R[1]; i++){ | ||
528 | + if(p1[i] == 0 && p2[i] ==0) | ||
529 | + result[i] = 1; | ||
530 | + else | ||
531 | + result[i] = p1[i] / p2[i]; | ||
532 | + } | ||
533 | + return true; | ||
534 | + } | ||
535 | + | ||
536 | + //peak area to peak area ratio | ||
537 | + bool pa_to_pa(double lb1, double rb1, double lab1, double rab1, | ||
538 | + double lb2, double rb2, double lab2, double rab2, T* result){ | ||
539 | + | ||
540 | + T* p1 = (T*)malloc(R[0] * R[1] * sizeof(T)); | ||
541 | + T* p2 = (T*)malloc(R[0] * R[1] * sizeof(T)); | ||
542 | + | ||
543 | + //get the area and the peak band | ||
544 | + area(lb1, rb1, lab1, rab1, p1); | ||
545 | + area(lb2, rb2, lab2, rab2, p2); | ||
546 | + //calculate the ratio in result | ||
547 | + for(unsigned i = 0; i < R[0] * R[1]; i++){ | ||
548 | + if(p1[i] == 0 && p2[i] ==0) | ||
549 | + result[i] = 1; | ||
550 | + else | ||
551 | + result[i] = p1[i] / p2[i]; | ||
552 | + } | ||
553 | + return true; | ||
554 | + } | ||
555 | + | ||
507 | //create mask file | 556 | //create mask file |
508 | bool mask(unsigned char* p, double mask_band, double threshold){ | 557 | bool mask(unsigned char* p, double mask_band, double threshold){ |
509 | 558 |
envi/bip.h
@@ -490,17 +490,16 @@ public: | @@ -490,17 +490,16 @@ public: | ||
490 | } | 490 | } |
491 | //untested | 491 | //untested |
492 | //providing the left and the right bound wavelength, return baseline-corrected band height | 492 | //providing the left and the right bound wavelength, return baseline-corrected band height |
493 | - bool height(double lb, double rb, double bandwavelength){ | 493 | + bool height(double lb, double rb, double bandwavelength, T* result){ |
494 | 494 | ||
495 | T* lp; | 495 | T* lp; |
496 | T* rp; | 496 | T* rp; |
497 | - T* result; | ||
498 | - | ||
499 | - lp = (T*) malloc(sizeof(R[0] * R[1] * sizeof(T))); //memory allocation | ||
500 | - rp = (T*) malloc(sizeof(R[0] * R[1] * sizeof(T))); | ||
501 | - result = (T*) malloc(sizeof(R[0] * R[1] * sizeof(T))); | 497 | + unsigned XY = R[0] * R[1]; |
498 | + unsigned S = XY * sizeof(T); | ||
499 | + lp = (T*) malloc(S); //memory allocation | ||
500 | + rp = (T*) malloc(S); | ||
502 | 501 | ||
503 | - band(lp, lb); //get the data of the left bounbd and the right bound | 502 | + band(lp, lb); |
504 | band(rp, rb); | 503 | band(rp, rb); |
505 | 504 | ||
506 | baseline_band(lb, rb, lp, rp, bandwavelength, result); | 505 | baseline_band(lb, rb, lp, rp, bandwavelength, result); |
@@ -510,13 +509,12 @@ public: | @@ -510,13 +509,12 @@ public: | ||
510 | 509 | ||
511 | 510 | ||
512 | //calculate the area between two bound point(including baseline correction) | 511 | //calculate the area between two bound point(including baseline correction) |
513 | - bool area(double lb, double rb){ | 512 | + bool area(double lb, double rb, double lab, double rab, T* result){ |
514 | 513 | ||
515 | T* lp; //left band pointer | 514 | T* lp; //left band pointer |
516 | T* rp; //right band pointer | 515 | T* rp; //right band pointer |
517 | T* cur; //current band 1 | 516 | T* cur; //current band 1 |
518 | T* cur2; //current band 2 | 517 | T* cur2; //current band 2 |
519 | - T* result; //area result | ||
520 | 518 | ||
521 | unsigned XY = R[0] * R[1]; | 519 | unsigned XY = R[0] * R[1]; |
522 | unsigned S = XY * sizeof(T); | 520 | unsigned S = XY * sizeof(T); |
@@ -525,7 +523,6 @@ public: | @@ -525,7 +523,6 @@ public: | ||
525 | rp = (T*) malloc(S); | 523 | rp = (T*) malloc(S); |
526 | cur = (T*) malloc(S); | 524 | cur = (T*) malloc(S); |
527 | cur2 = (T*) malloc(S); | 525 | cur2 = (T*) malloc(S); |
528 | - result = (T*) malloc(S); | ||
529 | 526 | ||
530 | memset(result, (char)0, S); | 527 | memset(result, (char)0, S); |
531 | 528 | ||
@@ -548,10 +545,10 @@ public: | @@ -548,10 +545,10 @@ public: | ||
548 | } | 545 | } |
549 | 546 | ||
550 | //get the position of lb and rb | 547 | //get the position of lb and rb |
551 | - while (lb >= w[ai]){ | 548 | + while (lab >= w[ai]){ |
552 | ai++; | 549 | ai++; |
553 | } | 550 | } |
554 | - while (rb <= w[bi]){ | 551 | + while (rab <= w[bi]){ |
555 | bi--; | 552 | bi--; |
556 | } | 553 | } |
557 | 554 | ||
@@ -559,13 +556,15 @@ public: | @@ -559,13 +556,15 @@ public: | ||
559 | band(rp, rb); | 556 | band(rp, rb); |
560 | 557 | ||
561 | //calculate the beginning and the ending part | 558 | //calculate the beginning and the ending part |
562 | - baseline_band(lb, rb, lp, rp, w[bi], cur); //ending part | 559 | + baseline_band(lb, rb, lp, rp, rab, cur2); //ending part |
560 | + baseline_band(lb, rb, lp, rp, w[bi], cur); | ||
563 | for(unsigned j = 0; j < XY; j++){ | 561 | for(unsigned j = 0; j < XY; j++){ |
564 | - result[j] += (rb - w[bi]) * cur[j] / 2.0; | 562 | + result[j] += (rab - w[bi]) * (cur[j] + cur2[j]) / 2.0; |
565 | } | 563 | } |
566 | - baseline_band(lb, rb, lp, rp, w[ai], cur); //beginnning part | 564 | + baseline_band(lb, rb, lp, rp, lab, cur2); //beginnning part |
565 | + baseline_band(lb, rb, lp, rp, w[ai], cur); | ||
567 | for(unsigned j = 0; j < XY; j++){ | 566 | for(unsigned j = 0; j < XY; j++){ |
568 | - result[j] += (w[ai] - lb) * cur[j] / 2.0; | 567 | + result[j] += (w[ai] - lab) * (cur[j] + cur2[j]) / 2.0; |
569 | } | 568 | } |
570 | //calculate the area | 569 | //calculate the area |
571 | ai++; | 570 | ai++; |
@@ -578,20 +577,68 @@ public: | @@ -578,20 +577,68 @@ public: | ||
578 | } | 577 | } |
579 | std::swap(cur,cur2); //swap the band pointers | 578 | std::swap(cur,cur2); //swap the band pointers |
580 | } | 579 | } |
581 | - ///////////////for testing use | ||
582 | - std::ofstream text("area_result.txt",std::ios::out); | ||
583 | - for (unsigned i = 0; i < 640; i++) | ||
584 | - { | ||
585 | - for (unsigned j = 0; j < 384; j++) | ||
586 | - { | ||
587 | - text<<result[i*384+j]<<", "; | ||
588 | - } | ||
589 | - text<<std::endl; | ||
590 | - } | ||
591 | - /////////////// | 580 | + |
592 | return true; | 581 | return true; |
593 | } | 582 | } |
594 | 583 | ||
584 | + //peak height ratio | ||
585 | + bool ph_to_ph(double lb1, double rb1, double pos1, double lb2, double rb2, double pos2, T * result){ | ||
586 | + | ||
587 | + T* p1 = (T*)malloc(R[0] * R[1] * sizeof(T)); | ||
588 | + T* p2 = (T*)malloc(R[0] * R[1] * sizeof(T)); | ||
589 | + | ||
590 | + //get the two peak band | ||
591 | + height(lb1, rb1, pos1, p1); | ||
592 | + height(lb2, rb2, pos2, p2); | ||
593 | + //calculate the ratio in result | ||
594 | + for(unsigned i = 0; i < R[0] * R[1]; i++){ | ||
595 | + if(p1[i] == 0 && p2[i] ==0) | ||
596 | + result[i] = 1; | ||
597 | + else | ||
598 | + result[i] = p1[i] / p2[i]; | ||
599 | + } | ||
600 | + return true; | ||
601 | + } | ||
602 | + | ||
603 | + //peak are to peak height ratio | ||
604 | + bool pa_to_ph(double lb1, double rb1, double lab1, double rab1, | ||
605 | + double lb2, double rb2, double pos, T* result){ | ||
606 | + | ||
607 | + T* p1 = (T*)malloc(R[0] * R[1] * sizeof(T)); | ||
608 | + T* p2 = (T*)malloc(R[0] * R[1] * sizeof(T)); | ||
609 | + | ||
610 | + //get the area and the peak band | ||
611 | + area(lb1, rb1, lab1, rab1, p1); | ||
612 | + height(lb2, rb2, pos, p2); | ||
613 | + //calculate the ratio in result | ||
614 | + for(unsigned i = 0; i < R[0] * R[1]; i++){ | ||
615 | + if(p1[i] == 0 && p2[i] ==0) | ||
616 | + result[i] = 1; | ||
617 | + else | ||
618 | + result[i] = p1[i] / p2[i]; | ||
619 | + } | ||
620 | + return true; | ||
621 | + } | ||
622 | + | ||
623 | + //peak area to peak area ratio | ||
624 | + bool pa_to_pa(double lb1, double rb1, double lab1, double rab1, | ||
625 | + double lb2, double rb2, double lab2, double rab2, T* result){ | ||
626 | + | ||
627 | + T* p1 = (T*)malloc(R[0] * R[1] * sizeof(T)); | ||
628 | + T* p2 = (T*)malloc(R[0] * R[1] * sizeof(T)); | ||
629 | + | ||
630 | + //get the area and the peak band | ||
631 | + area(lb1, rb1, lab1, rab1, p1); | ||
632 | + area(lb2, rb2, lab2, rab2, p2); | ||
633 | + //calculate the ratio in result | ||
634 | + for(unsigned i = 0; i < R[0] * R[1]; i++){ | ||
635 | + if(p1[i] == 0 && p2[i] ==0) | ||
636 | + result[i] = 1; | ||
637 | + else | ||
638 | + result[i] = p1[i] / p2[i]; | ||
639 | + } | ||
640 | + return true; | ||
641 | + } | ||
595 | //create mask file | 642 | //create mask file |
596 | bool mask(unsigned char* p, double mask_band, double threshold){ | 643 | bool mask(unsigned char* p, double mask_band, double threshold){ |
597 | 644 |
envi/bsq.h
@@ -50,7 +50,11 @@ public: | @@ -50,7 +50,11 @@ public: | ||
50 | return false; | 50 | return false; |
51 | } | 51 | } |
52 | 52 | ||
53 | - getSlice(p, 2, page); | 53 | + file.seekg(R[0] * R[1] * page * sizeof(T), std::ios::beg); |
54 | + | ||
55 | + file.read((char *)p, sizeof(T) * R[0] * R[1]); | ||
56 | + | ||
57 | +// getSlice(p, 2, page); | ||
54 | return true; | 58 | return true; |
55 | } | 59 | } |
56 | 60 | ||
@@ -93,7 +97,7 @@ public: | @@ -93,7 +97,7 @@ public: | ||
93 | } | 97 | } |
94 | else //if the wavelength is equal to a wavelength in header file | 98 | else //if the wavelength is equal to a wavelength in header file |
95 | { | 99 | { |
96 | - getSlice(p, 2, page); | 100 | + band_index(p, page); |
97 | } | 101 | } |
98 | 102 | ||
99 | return true; | 103 | return true; |
@@ -325,17 +329,16 @@ public: | @@ -325,17 +329,16 @@ public: | ||
325 | } | 329 | } |
326 | //untested | 330 | //untested |
327 | //providing the left and the right bound wavelength, return baseline-corrected band height | 331 | //providing the left and the right bound wavelength, return baseline-corrected band height |
328 | - bool height(double lb, double rb, double bandwavelength){ | 332 | + bool height(double lb, double rb, double bandwavelength, T* result){ |
329 | 333 | ||
330 | T* lp; | 334 | T* lp; |
331 | T* rp; | 335 | T* rp; |
332 | - T* result; | ||
333 | - | ||
334 | - lp = (T*) malloc(sizeof(R[0] * R[1] * sizeof(T))); //memory allocation | ||
335 | - rp = (T*) malloc(sizeof(R[0] * R[1] * sizeof(T))); | ||
336 | - result = (T*) malloc(sizeof(R[0] * R[1] * sizeof(T))); | 336 | + unsigned XY = R[0] * R[1]; |
337 | + unsigned S = XY * sizeof(T); | ||
338 | + lp = (T*) malloc(S); //memory allocation | ||
339 | + rp = (T*) malloc(S); | ||
337 | 340 | ||
338 | - band(lp, lb); //get the data of the left bounbd and the right bound | 341 | + band(lp, lb); |
339 | band(rp, rb); | 342 | band(rp, rb); |
340 | 343 | ||
341 | baseline_band(lb, rb, lp, rp, bandwavelength, result); | 344 | baseline_band(lb, rb, lp, rp, bandwavelength, result); |
@@ -345,13 +348,12 @@ public: | @@ -345,13 +348,12 @@ public: | ||
345 | 348 | ||
346 | 349 | ||
347 | //calculate the area between two bound point(including baseline correction) | 350 | //calculate the area between two bound point(including baseline correction) |
348 | - bool area(double lb, double rb){ | 351 | + bool area(double lb, double rb, double lab, double rab, T* result){ |
349 | 352 | ||
350 | T* lp; //left band pointer | 353 | T* lp; //left band pointer |
351 | T* rp; //right band pointer | 354 | T* rp; //right band pointer |
352 | T* cur; //current band 1 | 355 | T* cur; //current band 1 |
353 | T* cur2; //current band 2 | 356 | T* cur2; //current band 2 |
354 | - T* result; //area result | ||
355 | 357 | ||
356 | unsigned XY = R[0] * R[1]; | 358 | unsigned XY = R[0] * R[1]; |
357 | unsigned S = XY * sizeof(T); | 359 | unsigned S = XY * sizeof(T); |
@@ -360,7 +362,6 @@ public: | @@ -360,7 +362,6 @@ public: | ||
360 | rp = (T*) malloc(S); | 362 | rp = (T*) malloc(S); |
361 | cur = (T*) malloc(S); | 363 | cur = (T*) malloc(S); |
362 | cur2 = (T*) malloc(S); | 364 | cur2 = (T*) malloc(S); |
363 | - result = (T*) malloc(S); | ||
364 | 365 | ||
365 | memset(result, (char)0, S); | 366 | memset(result, (char)0, S); |
366 | 367 | ||
@@ -383,10 +384,10 @@ public: | @@ -383,10 +384,10 @@ public: | ||
383 | } | 384 | } |
384 | 385 | ||
385 | //get the position of lb and rb | 386 | //get the position of lb and rb |
386 | - while (lb >= w[ai]){ | 387 | + while (lab >= w[ai]){ |
387 | ai++; | 388 | ai++; |
388 | } | 389 | } |
389 | - while (rb <= w[bi]){ | 390 | + while (rab <= w[bi]){ |
390 | bi--; | 391 | bi--; |
391 | } | 392 | } |
392 | 393 | ||
@@ -394,14 +395,17 @@ public: | @@ -394,14 +395,17 @@ public: | ||
394 | band(rp, rb); | 395 | band(rp, rb); |
395 | 396 | ||
396 | //calculate the beginning and the ending part | 397 | //calculate the beginning and the ending part |
397 | - baseline_band(lb, rb, lp, rp, w[bi], cur); //ending part | 398 | + baseline_band(lb, rb, lp, rp, rab, cur2); //ending part |
399 | + baseline_band(lb, rb, lp, rp, w[bi], cur); | ||
398 | for(unsigned j = 0; j < XY; j++){ | 400 | for(unsigned j = 0; j < XY; j++){ |
399 | - result[j] += (rb - w[bi]) * cur[j] / 2.0; | 401 | + result[j] += (rab - w[bi]) * (cur[j] + cur2[j]) / 2.0; |
400 | } | 402 | } |
401 | - baseline_band(lb, rb, lp, rp, w[ai], cur); //beginnning part | 403 | + baseline_band(lb, rb, lp, rp, lab, cur2); //beginnning part |
404 | + baseline_band(lb, rb, lp, rp, w[ai], cur); | ||
402 | for(unsigned j = 0; j < XY; j++){ | 405 | for(unsigned j = 0; j < XY; j++){ |
403 | - result[j] += (w[ai] - lb) * cur[j] / 2.0; | 406 | + result[j] += (w[ai] - lab) * (cur[j] + cur2[j]) / 2.0; |
404 | } | 407 | } |
408 | + | ||
405 | //calculate the area | 409 | //calculate the area |
406 | ai++; | 410 | ai++; |
407 | for(unsigned i = ai; i <= bi ;i++) | 411 | for(unsigned i = ai; i <= bi ;i++) |
@@ -413,21 +417,68 @@ public: | @@ -413,21 +417,68 @@ public: | ||
413 | } | 417 | } |
414 | std::swap(cur,cur2); //swap the band pointers | 418 | std::swap(cur,cur2); //swap the band pointers |
415 | } | 419 | } |
416 | - ///////////////for testing use | ||
417 | - std::ofstream text("area_result.txt",std::ios::out); | ||
418 | - for (unsigned i = 0; i < 640; i++) | ||
419 | - { | ||
420 | - for (unsigned j = 0; j < 384; j++) | ||
421 | - { | ||
422 | - text<<result[i*384+j]<<", "; | ||
423 | - } | ||
424 | - text<<std::endl; | ||
425 | - } | ||
426 | - /////////////// | 420 | + |
427 | return true; | 421 | return true; |
428 | } | 422 | } |
429 | 423 | ||
424 | + //peak height ratio | ||
425 | + bool ph_to_ph(double lb1, double rb1, double pos1, double lb2, double rb2, double pos2, T * result){ | ||
430 | 426 | ||
427 | + T* p1 = (T*)malloc(R[0] * R[1] * sizeof(T)); | ||
428 | + T* p2 = (T*)malloc(R[0] * R[1] * sizeof(T)); | ||
429 | + | ||
430 | + //get the two peak band | ||
431 | + height(lb1, rb1, pos1, p1); | ||
432 | + height(lb2, rb2, pos2, p2); | ||
433 | + //calculate the ratio in result | ||
434 | + for(unsigned i = 0; i < R[0] * R[1]; i++){ | ||
435 | + if(p1[i] == 0 && p2[i] ==0) | ||
436 | + result[i] = 1; | ||
437 | + else | ||
438 | + result[i] = p1[i] / p2[i]; | ||
439 | + } | ||
440 | + return true; | ||
441 | + } | ||
442 | + | ||
443 | + //peak are to peak height ratio | ||
444 | + bool pa_to_ph(double lb1, double rb1, double lab1, double rab1, | ||
445 | + double lb2, double rb2, double pos, T* result){ | ||
446 | + | ||
447 | + T* p1 = (T*)malloc(R[0] * R[1] * sizeof(T)); | ||
448 | + T* p2 = (T*)malloc(R[0] * R[1] * sizeof(T)); | ||
449 | + | ||
450 | + //get the area and the peak band | ||
451 | + area(lb1, rb1, lab1, rab1, p1); | ||
452 | + height(lb2, rb2, pos, p2); | ||
453 | + //calculate the ratio in result | ||
454 | + for(unsigned i = 0; i < R[0] * R[1]; i++){ | ||
455 | + if(p1[i] == 0 && p2[i] ==0) | ||
456 | + result[i] = 1; | ||
457 | + else | ||
458 | + result[i] = p1[i] / p2[i]; | ||
459 | + } | ||
460 | + return true; | ||
461 | + } | ||
462 | + | ||
463 | + //peak area to peak area ratio | ||
464 | + bool pa_to_pa(double lb1, double rb1, double lab1, double rab1, | ||
465 | + double lb2, double rb2, double lab2, double rab2, T* result){ | ||
466 | + | ||
467 | + T* p1 = (T*)malloc(R[0] * R[1] * sizeof(T)); | ||
468 | + T* p2 = (T*)malloc(R[0] * R[1] * sizeof(T)); | ||
469 | + | ||
470 | + //get the area and the peak band | ||
471 | + area(lb1, rb1, lab1, rab1, p1); | ||
472 | + area(lb2, rb2, lab2, rab2, p2); | ||
473 | + //calculate the ratio in result | ||
474 | + for(unsigned i = 0; i < R[0] * R[1]; i++){ | ||
475 | + if(p1[i] == 0 && p2[i] ==0) | ||
476 | + result[i] = 1; | ||
477 | + else | ||
478 | + result[i] = p1[i] / p2[i]; | ||
479 | + } | ||
480 | + return true; | ||
481 | + } | ||
431 | //create mask file | 482 | //create mask file |
432 | bool mask(unsigned char* p, double mask_band, double threshold){ | 483 | bool mask(unsigned char* p, double mask_band, double threshold){ |
433 | 484 |
envi/envi.h
@@ -310,31 +310,104 @@ public: | @@ -310,31 +310,104 @@ public: | ||
310 | return false; | 310 | return false; |
311 | } | 311 | } |
312 | 312 | ||
313 | - //do baseline correction and caculate the area between two bounds | ||
314 | - bool area(double lb, double rb){ | 313 | + //peak height to peak height ratio |
314 | + bool ph_to_ph(double lb1, double rb1, double pos1, double lb2, double rb2, double pos2, void * result){ | ||
315 | if(header.interleave == envi_header::BSQ){ //if the infile is bsq file | 315 | if(header.interleave == envi_header::BSQ){ //if the infile is bsq file |
316 | if(header.data_type ==envi_header::float32) | 316 | if(header.data_type ==envi_header::float32) |
317 | - return ((bsq<float>*)file)->area(lb, rb); | 317 | + return ((bsq<float>*)file)->ph_to_ph(lb1, rb1, pos1, lb2, rb2, pos2, (float*)result); |
318 | else if(header.data_type == envi_header::float64) | 318 | else if(header.data_type == envi_header::float64) |
319 | - return ((bsq<double>*)file)->area(lb, rb); | 319 | + return ((bsq<double>*)file)->ph_to_ph(lb1, rb1, pos1, lb2, rb2, pos2, (double*)result); |
320 | else | 320 | else |
321 | std::cout<<"ERROR: unidentified data type"<<std::endl; | 321 | std::cout<<"ERROR: unidentified data type"<<std::endl; |
322 | } | 322 | } |
323 | 323 | ||
324 | else if(header.interleave == envi_header::BIL){ //if the infile is bil file | 324 | else if(header.interleave == envi_header::BIL){ //if the infile is bil file |
325 | if(header.data_type ==envi_header::float32) | 325 | if(header.data_type ==envi_header::float32) |
326 | - return ((bil<float>*)file)->area(lb, rb); | 326 | + return ((bil<float>*)file)->ph_to_ph(lb1, rb1, pos1, lb2, rb2, pos2, (float*)result); |
327 | else if(header.data_type == envi_header::float64) | 327 | else if(header.data_type == envi_header::float64) |
328 | - return ((bil<double>*)file)->area(lb, rb); | 328 | + return ((bil<double>*)file)->ph_to_ph(lb1, rb1, pos1, lb2, rb2, pos2, (double*)result); |
329 | else | 329 | else |
330 | std::cout<<"ERROR: unidentified data type"<<std::endl; | 330 | std::cout<<"ERROR: unidentified data type"<<std::endl; |
331 | } | 331 | } |
332 | 332 | ||
333 | else if(header.interleave == envi_header::BIP){ //if the infile is bip file | 333 | else if(header.interleave == envi_header::BIP){ //if the infile is bip file |
334 | if(header.data_type ==envi_header::float32) | 334 | if(header.data_type ==envi_header::float32) |
335 | - return ((bip<float>*)file)->area(lb, rb); | 335 | + return ((bip<float>*)file)->ph_to_ph(lb1, rb1, pos1, lb2, rb2, pos2, (float*)result); |
336 | else if(header.data_type == envi_header::float64) | 336 | else if(header.data_type == envi_header::float64) |
337 | - return ((bip<double>*)file)->area(lb, rb); | 337 | + return ((bip<double>*)file)->ph_to_ph(lb1, rb1, pos1, lb2, rb2, pos2, (double*)result); |
338 | + else | ||
339 | + std::cout<<"ERROR: unidentified data type"<<std::endl; | ||
340 | + } | ||
341 | + | ||
342 | + else{ | ||
343 | + std::cout<<"ERROR: unidentified file type"<<std::endl; | ||
344 | + exit(1); | ||
345 | + } | ||
346 | + return false; | ||
347 | + } | ||
348 | + | ||
349 | + //peak area to peak height ratio | ||
350 | + bool pa_to_ph(double lb1, double rb1, double lab1, double rab1, double lb2, double rb2, double pos, void * result){ | ||
351 | + if(header.interleave == envi_header::BSQ){ //if the infile is bsq file | ||
352 | + if(header.data_type ==envi_header::float32) | ||
353 | + return ((bsq<float>*)file)->pa_to_ph(lb1, rb1, lab1, rab1, lb2, rb2, pos, (float*)result); | ||
354 | + else if(header.data_type == envi_header::float64) | ||
355 | + return ((bsq<double>*)file)->pa_to_ph(lb1, rb1, lab1, rab1, lb2, rb2, pos, (double*)result); | ||
356 | + else | ||
357 | + std::cout<<"ERROR: unidentified data type"<<std::endl; | ||
358 | + } | ||
359 | + | ||
360 | + else if(header.interleave == envi_header::BIL){ //if the infile is bil file | ||
361 | + if(header.data_type ==envi_header::float32) | ||
362 | + return ((bil<float>*)file)->pa_to_ph(lb1, rb1, lab1, rab1, lb2, rb2, pos, (float*)result); | ||
363 | + else if(header.data_type == envi_header::float64) | ||
364 | + return ((bil<double>*)file)->pa_to_ph(lb1, rb1, lab1, rab1, lb2, rb2, pos, (double*)result); | ||
365 | + else | ||
366 | + std::cout<<"ERROR: unidentified data type"<<std::endl; | ||
367 | + } | ||
368 | + | ||
369 | + else if(header.interleave == envi_header::BIP){ //if the infile is bip file | ||
370 | + if(header.data_type ==envi_header::float32) | ||
371 | + return ((bip<float>*)file)->pa_to_ph(lb1, rb1, lab1, rab1, lb2, rb2, pos, (float*)result); | ||
372 | + else if(header.data_type == envi_header::float64) | ||
373 | + return ((bip<double>*)file)->pa_to_ph(lb1, rb1, lab1, rab1, lb2, rb2, pos, (double*)result); | ||
374 | + else | ||
375 | + std::cout<<"ERROR: unidentified data type"<<std::endl; | ||
376 | + } | ||
377 | + | ||
378 | + else{ | ||
379 | + std::cout<<"ERROR: unidentified file type"<<std::endl; | ||
380 | + exit(1); | ||
381 | + } | ||
382 | + return false; | ||
383 | + } | ||
384 | + | ||
385 | + //peak area to peak area ratio | ||
386 | + bool pa_to_pa(double lb1, double rb1, double lab1, double rab1, | ||
387 | + double lb2, double rb2, double lab2, double rab2, void* result){ | ||
388 | + if(header.interleave == envi_header::BSQ){ //if the infile is bsq file | ||
389 | + if(header.data_type ==envi_header::float32) | ||
390 | + return ((bsq<float>*)file)->pa_to_pa(lb1, rb1, lab1, rab1, lb2, rb2, lab2, rab2, (float*)result); | ||
391 | + else if(header.data_type == envi_header::float64) | ||
392 | + return ((bsq<double>*)file)->pa_to_pa(lb1, rb1, lab1, rab1, lb2, rb2, lab2, rab2, (double*)result); | ||
393 | + else | ||
394 | + std::cout<<"ERROR: unidentified data type"<<std::endl; | ||
395 | + } | ||
396 | + | ||
397 | + else if(header.interleave == envi_header::BIL){ //if the infile is bil file | ||
398 | + if(header.data_type ==envi_header::float32) | ||
399 | + return ((bil<float>*)file)->pa_to_pa(lb1, rb1, lab1, rab1, lb2, rb2, lab2, rab2, (float*)result); | ||
400 | + else if(header.data_type == envi_header::float64) | ||
401 | + return ((bil<double>*)file)->pa_to_pa(lb1, rb1, lab1, rab1, lb2, rb2, lab2, rab2, (double*)result); | ||
402 | + else | ||
403 | + std::cout<<"ERROR: unidentified data type"<<std::endl; | ||
404 | + } | ||
405 | + | ||
406 | + else if(header.interleave == envi_header::BIP){ //if the infile is bip file | ||
407 | + if(header.data_type ==envi_header::float32) | ||
408 | + return ((bip<float>*)file)->pa_to_pa(lb1, rb1, lab1, rab1, lb2, rb2, lab2, rab2, (float*)result); | ||
409 | + else if(header.data_type == envi_header::float64) | ||
410 | + return ((bip<double>*)file)->pa_to_pa(lb1, rb1, lab1, rab1, lb2, rb2, lab2, rab2, (double*)result); | ||
338 | else | 411 | else |
339 | std::cout<<"ERROR: unidentified data type"<<std::endl; | 412 | std::cout<<"ERROR: unidentified data type"<<std::endl; |
340 | } | 413 | } |