Commit 517876d61f3637a07b5e36246be172cc1608ba0f

Authored by heziqi
1 parent d1d8ef91

metrics finished and memory allocation error fixed

Showing 4 changed files with 398 additions and 13 deletions   Show diff stats
envi/bil.h
... ... @@ -60,8 +60,7 @@ public:
60 60 unsigned int S = XY * sizeof(T); //calculate the number of bytes of a band
61 61  
62 62 unsigned page=0; //bands around the wavelength
63   - T * p1;
64   - T * p2;
  63 +
65 64  
66 65 //get the bands numbers around the wavelength
67 66  
... ... @@ -82,6 +81,8 @@ public:
82 81 }
83 82 if ( wavelength < w[page] )
84 83 {
  84 + T * p1;
  85 + T * p2;
85 86 p1=(T*)malloc(S); //memory allocation
86 87 p2=(T*)malloc(S);
87 88 band_index(p1, page - 1);
... ... @@ -90,6 +91,8 @@ public:
90 91 double r = (double) (wavelength - w[page-1]) / (double) (w[page] - w[page-1]);
91 92 p[i] = (p2[i] - p1[i]) * r + p1[i];
92 93 }
  94 + free(p1);
  95 + free(p2);
93 96 }
94 97 else //if the wavelength is equal to a wavelength in header file
95 98 {
... ... @@ -401,8 +404,6 @@ public:
401 404 }
402 405 return true;
403 406 }
404   - //untested
405   - //providing the left and the right bound wavelength, return baseline-corrected band height
406 407 bool height(double lb, double rb, double bandwavelength, T* result){
407 408  
408 409 T* lp;
... ... @@ -417,6 +418,8 @@ public:
417 418  
418 419 baseline_band(lb, rb, lp, rp, bandwavelength, result);
419 420  
  421 + free(lp);
  422 + free(rp);
420 423 return true;
421 424 }
422 425  
... ... @@ -479,6 +482,7 @@ public:
479 482 for(unsigned j = 0; j < XY; j++){
480 483 result[j] += (w[ai] - lab) * (cur[j] + cur2[j]) / 2.0;
481 484 }
  485 +
482 486 //calculate the area
483 487 ai++;
484 488 for(unsigned i = ai; i <= bi ;i++)
... ... @@ -491,6 +495,10 @@ public:
491 495 std::swap(cur,cur2); //swap the band pointers
492 496 }
493 497  
  498 + free(lp);
  499 + free(rp);
  500 + free(cur);
  501 + free(cur2);
494 502 return true;
495 503 }
496 504  
... ... @@ -510,6 +518,9 @@ public:
510 518 else
511 519 result[i] = p1[i] / p2[i];
512 520 }
  521 +
  522 + free(p1);
  523 + free(p2);
513 524 return true;
514 525 }
515 526  
... ... @@ -530,6 +541,9 @@ public:
530 541 else
531 542 result[i] = p1[i] / p2[i];
532 543 }
  544 +
  545 + free(p1);
  546 + free(p2);
533 547 return true;
534 548 }
535 549  
... ... @@ -550,9 +564,108 @@ public:
550 564 else
551 565 result[i] = p1[i] / p2[i];
552 566 }
  567 +
  568 + free(p1);
  569 + free(p2);
553 570 return true;
554 571 }
555 572  
  573 + //x * f(x)
  574 + bool x_area(double lb, double rb, double lab, double rab, T* result){
  575 + T* lp; //left band pointer
  576 + T* rp; //right band pointer
  577 + T* cur; //current band 1
  578 + T* cur2; //current band 2
  579 +
  580 + unsigned XY = R[0] * R[1];
  581 + unsigned S = XY * sizeof(T);
  582 +
  583 + lp = (T*) malloc(S); //memory allocation
  584 + rp = (T*) malloc(S);
  585 + cur = (T*) malloc(S);
  586 + cur2 = (T*) malloc(S);
  587 +
  588 + memset(result, (char)0, S);
  589 +
  590 + //find the wavelenght position in the whole band
  591 + unsigned int n = w.size();
  592 + unsigned int ai = 0; //left bound position
  593 + unsigned int bi = n - 1; //right bound position
  594 +
  595 + //to make sure the left and the right bound are in the bandwidth
  596 + if (lb < w[0] || rb < w[0] || lb > w[n-1] || rb >w[n-1]){
  597 + std::cout<<"ERROR: left bound or right bound out of bandwidth"<<std::endl;
  598 + exit(1);
  599 + }
  600 + //to make sure rigth bound is bigger than left bound
  601 + else if(lb > rb){
  602 + std::cout<<"ERROR: right bound should be bigger than left bound"<<std::endl;
  603 + exit(1);
  604 + }
  605 +
  606 + //get the position of lb and rb
  607 + while (lab >= w[ai]){
  608 + ai++;
  609 + }
  610 + while (rab <= w[bi]){
  611 + bi--;
  612 + }
  613 +
  614 + band(lp, lb);
  615 + band(rp, rb);
  616 +
  617 + //calculate the beginning and the ending part
  618 + baseline_band(lb, rb, lp, rp, rab, cur2); //ending part
  619 + baseline_band(lb, rb, lp, rp, w[bi], cur);
  620 + for(unsigned j = 0; j < XY; j++){
  621 + result[j] += (rab - w[bi]) * (rab + w[bi]) * (cur[j] + cur2[j]) / 4.0;
  622 + }
  623 + baseline_band(lb, rb, lp, rp, lab, cur2); //beginnning part
  624 + baseline_band(lb, rb, lp, rp, w[ai], cur);
  625 + for(unsigned j = 0; j < XY; j++){
  626 + result[j] += (w[ai] - lab) * (w[ai] + lab) * (cur[j] + cur2[j]) / 4.0;
  627 + }
  628 +
  629 + //calculate f(x) times x
  630 + ai++;
  631 + for(unsigned i = ai; i <= bi ;i++)
  632 + {
  633 + baseline_band(lb, rb, lp, rp, w[ai], cur2);
  634 + for(unsigned j = 0; j < XY; j++)
  635 + {
  636 + result[j] += (w[ai] - w[ai-1]) * (w[ai] + w[ai-1]) * (cur[j] + cur2[j]) / 4.0;
  637 + }
  638 + std::swap(cur,cur2); //swap the band pointers
  639 + }
  640 +
  641 + free(lp);
  642 + free(rp);
  643 + free(cur);
  644 + free(cur2);
  645 + return true;
  646 + }
  647 +
  648 + //centroid point
  649 + bool cpoint(double lb, double rb, double lab, double rab, T* result){
  650 + T* p1 = (T*)malloc(R[0] * R[1] * sizeof(T));
  651 + T* p2 = (T*)malloc(R[0] * R[1] * sizeof(T));
  652 +
  653 + //get the area and the peak band
  654 + x_area(lb, rb, lab, rab, p1);
  655 + area(lb, rb, lab, rab, p2);
  656 + //calculate the ratio in result
  657 + for(unsigned i = 0; i < R[0] * R[1]; i++){
  658 + if(p1[i] == 0 && p2[i] ==0)
  659 + result[i] = 1;
  660 + else
  661 + result[i] = p1[i] / p2[i];
  662 + }
  663 +
  664 + free(p1);
  665 + free(p2);
  666 + return true;
  667 + }
  668 +
556 669 //create mask file
557 670 bool mask(unsigned char* p, double mask_band, double threshold){
558 671  
... ... @@ -566,6 +679,7 @@ public:
566 679 p[i] = 255;
567 680 }
568 681  
  682 + free(temp);
569 683 return true;
570 684 }
571 685  
... ...
envi/bip.h
... ... @@ -63,8 +63,7 @@ public:
63 63 unsigned int XY = R[0] * R[1]; //calculate the number of pixels in a band
64 64  
65 65 unsigned page=0; //bands around the wavelength
66   - T * p1;
67   - T * p2;
  66 +
68 67  
69 68 //get the bands numbers around the wavelength
70 69  
... ... @@ -85,6 +84,8 @@ public:
85 84 }
86 85 if ( wavelength < w[page] )
87 86 {
  87 + T * p1;
  88 + T * p2;
88 89 p1=(T*)malloc( XY * sizeof(T)); //memory allocation
89 90 p2=(T*)malloc( XY * sizeof(T));
90 91 band_index(p1, page - 1);
... ... @@ -93,6 +94,8 @@ public:
93 94 double r = (double) (wavelength - w[page-1]) / (double) (w[page] - w[page-1]);
94 95 p[i] = (p2[i] - p1[i]) * r + p1[i];
95 96 }
  97 + free(p1);
  98 + free(p2);
96 99 }
97 100 else //if the wavelength is equal to a wavelength in header file
98 101 {
... ... @@ -108,8 +111,7 @@ public:
108 111 unsigned int B = R[2];
109 112  
110 113 unsigned page=0; //samples around the wavelength
111   - T * p1;
112   - T * p2;
  114 +
113 115  
114 116 //get the bands numbers around the wavelength
115 117  
... ... @@ -136,6 +138,8 @@ public:
136 138 }
137 139 if ( wavelength < w[page] )
138 140 {
  141 + T * p1;
  142 + T * p2;
139 143 p1=(T*)malloc( X * sizeof(T)); //memory allocation
140 144 p2=(T*)malloc( X * sizeof(T));
141 145 //band_index(p1, page - 1);
... ... @@ -153,6 +157,8 @@ public:
153 157 double r = (double) (wavelength - w[page-1]) / (double) (w[page] - w[page-1]);
154 158 p[i] = (p2[i] - p1[i]) * r + p1[i];
155 159 }
  160 + free(p1);
  161 + free(p2);
156 162 }
157 163 else //if the wavelength is equal to a wavelength in header file
158 164 {
... ... @@ -488,8 +494,7 @@ public:
488 494 }
489 495 return true;
490 496 }
491   - //untested
492   - //providing the left and the right bound wavelength, return baseline-corrected band height
  497 +
493 498 bool height(double lb, double rb, double bandwavelength, T* result){
494 499  
495 500 T* lp;
... ... @@ -504,6 +509,8 @@ public:
504 509  
505 510 baseline_band(lb, rb, lp, rp, bandwavelength, result);
506 511  
  512 + free(lp);
  513 + free(rp);
507 514 return true;
508 515 }
509 516  
... ... @@ -566,6 +573,7 @@ public:
566 573 for(unsigned j = 0; j < XY; j++){
567 574 result[j] += (w[ai] - lab) * (cur[j] + cur2[j]) / 2.0;
568 575 }
  576 +
569 577 //calculate the area
570 578 ai++;
571 579 for(unsigned i = ai; i <= bi ;i++)
... ... @@ -578,6 +586,10 @@ public:
578 586 std::swap(cur,cur2); //swap the band pointers
579 587 }
580 588  
  589 + free(lp);
  590 + free(rp);
  591 + free(cur);
  592 + free(cur2);
581 593 return true;
582 594 }
583 595  
... ... @@ -597,6 +609,9 @@ public:
597 609 else
598 610 result[i] = p1[i] / p2[i];
599 611 }
  612 +
  613 + free(p1);
  614 + free(p2);
600 615 return true;
601 616 }
602 617  
... ... @@ -617,6 +632,9 @@ public:
617 632 else
618 633 result[i] = p1[i] / p2[i];
619 634 }
  635 +
  636 + free(p1);
  637 + free(p2);
620 638 return true;
621 639 }
622 640  
... ... @@ -637,8 +655,108 @@ public:
637 655 else
638 656 result[i] = p1[i] / p2[i];
639 657 }
  658 +
  659 + free(p1);
  660 + free(p2);
640 661 return true;
641 662 }
  663 +
  664 + //x * f(x)
  665 + bool x_area(double lb, double rb, double lab, double rab, T* result){
  666 + T* lp; //left band pointer
  667 + T* rp; //right band pointer
  668 + T* cur; //current band 1
  669 + T* cur2; //current band 2
  670 +
  671 + unsigned XY = R[0] * R[1];
  672 + unsigned S = XY * sizeof(T);
  673 +
  674 + lp = (T*) malloc(S); //memory allocation
  675 + rp = (T*) malloc(S);
  676 + cur = (T*) malloc(S);
  677 + cur2 = (T*) malloc(S);
  678 +
  679 + memset(result, (char)0, S);
  680 +
  681 + //find the wavelenght position in the whole band
  682 + unsigned int n = w.size();
  683 + unsigned int ai = 0; //left bound position
  684 + unsigned int bi = n - 1; //right bound position
  685 +
  686 + //to make sure the left and the right bound are in the bandwidth
  687 + if (lb < w[0] || rb < w[0] || lb > w[n-1] || rb >w[n-1]){
  688 + std::cout<<"ERROR: left bound or right bound out of bandwidth"<<std::endl;
  689 + exit(1);
  690 + }
  691 + //to make sure rigth bound is bigger than left bound
  692 + else if(lb > rb){
  693 + std::cout<<"ERROR: right bound should be bigger than left bound"<<std::endl;
  694 + exit(1);
  695 + }
  696 +
  697 + //get the position of lb and rb
  698 + while (lab >= w[ai]){
  699 + ai++;
  700 + }
  701 + while (rab <= w[bi]){
  702 + bi--;
  703 + }
  704 +
  705 + band(lp, lb);
  706 + band(rp, rb);
  707 +
  708 + //calculate the beginning and the ending part
  709 + baseline_band(lb, rb, lp, rp, rab, cur2); //ending part
  710 + baseline_band(lb, rb, lp, rp, w[bi], cur);
  711 + for(unsigned j = 0; j < XY; j++){
  712 + result[j] += (rab - w[bi]) * (rab + w[bi]) * (cur[j] + cur2[j]) / 4.0;
  713 + }
  714 + baseline_band(lb, rb, lp, rp, lab, cur2); //beginnning part
  715 + baseline_band(lb, rb, lp, rp, w[ai], cur);
  716 + for(unsigned j = 0; j < XY; j++){
  717 + result[j] += (w[ai] - lab) * (w[ai] + lab) * (cur[j] + cur2[j]) / 4.0;
  718 + }
  719 +
  720 + //calculate f(x) times x
  721 + ai++;
  722 + for(unsigned i = ai; i <= bi ;i++)
  723 + {
  724 + baseline_band(lb, rb, lp, rp, w[ai], cur2);
  725 + for(unsigned j = 0; j < XY; j++)
  726 + {
  727 + result[j] += (w[ai] - w[ai-1]) * (w[ai] + w[ai-1]) * (cur[j] + cur2[j]) / 4.0;
  728 + }
  729 + std::swap(cur,cur2); //swap the band pointers
  730 + }
  731 +
  732 + free(lp);
  733 + free(rp);
  734 + free(cur);
  735 + free(cur2);
  736 + return true;
  737 + }
  738 +
  739 + //centroid point
  740 + bool cpoint(double lb, double rb, double lab, double rab, T* result){
  741 + T* p1 = (T*)malloc(R[0] * R[1] * sizeof(T));
  742 + T* p2 = (T*)malloc(R[0] * R[1] * sizeof(T));
  743 +
  744 + //get the area and the peak band
  745 + x_area(lb, rb, lab, rab, p1);
  746 + area(lb, rb, lab, rab, p2);
  747 + //calculate the ratio in result
  748 + for(unsigned i = 0; i < R[0] * R[1]; i++){
  749 + if(p1[i] == 0 && p2[i] ==0)
  750 + result[i] = 1;
  751 + else
  752 + result[i] = p1[i] / p2[i];
  753 + }
  754 +
  755 + free(p1);
  756 + free(p2);
  757 + return true;
  758 + }
  759 +
642 760 //create mask file
643 761 bool mask(unsigned char* p, double mask_band, double threshold){
644 762  
... ... @@ -652,6 +770,7 @@ public:
652 770 p[i] = 255;
653 771 }
654 772  
  773 + free(temp);
655 774 return true;
656 775  
657 776 }
... ...
envi/bsq.h
... ... @@ -64,8 +64,7 @@ public:
64 64 unsigned int XY = R[0] * R[1]; //calculate the number of pixels in a band
65 65  
66 66 unsigned page=0; //bands around the wavelength
67   - T * p1;
68   - T * p2;
  67 +
69 68  
70 69 //get the bands numbers around the wavelength
71 70  
... ... @@ -86,6 +85,8 @@ public:
86 85 }
87 86 if ( wavelength < w[page] )
88 87 {
  88 + T * p1;
  89 + T * p2;
89 90 p1=(T*)malloc( XY * sizeof(T)); //memory allocation
90 91 p2=(T*)malloc( XY * sizeof(T));
91 92 band_index(p1, page - 1);
... ... @@ -94,6 +95,8 @@ public:
94 95 double r = (double) (wavelength - w[page-1]) / (double) (w[page] - w[page-1]);
95 96 p[i] = (p2[i] - p1[i]) * r + p1[i];
96 97 }
  98 + free(p1);
  99 + free(p2);
97 100 }
98 101 else //if the wavelength is equal to a wavelength in header file
99 102 {
... ... @@ -327,7 +330,7 @@ public:
327 330 }
328 331 return true;
329 332 }
330   - //untested
  333 +
331 334 //providing the left and the right bound wavelength, return baseline-corrected band height
332 335 bool height(double lb, double rb, double bandwavelength, T* result){
333 336  
... ... @@ -343,6 +346,8 @@ public:
343 346  
344 347 baseline_band(lb, rb, lp, rp, bandwavelength, result);
345 348  
  349 + free(lp);
  350 + free(rp);
346 351 return true;
347 352 }
348 353  
... ... @@ -418,6 +423,10 @@ public:
418 423 std::swap(cur,cur2); //swap the band pointers
419 424 }
420 425  
  426 + free(lp);
  427 + free(rp);
  428 + free(cur);
  429 + free(cur2);
421 430 return true;
422 431 }
423 432  
... ... @@ -437,6 +446,9 @@ public:
437 446 else
438 447 result[i] = p1[i] / p2[i];
439 448 }
  449 +
  450 + free(p1);
  451 + free(p2);
440 452 return true;
441 453 }
442 454  
... ... @@ -457,6 +469,9 @@ public:
457 469 else
458 470 result[i] = p1[i] / p2[i];
459 471 }
  472 +
  473 + free(p1);
  474 + free(p2);
460 475 return true;
461 476 }
462 477  
... ... @@ -477,8 +492,108 @@ public:
477 492 else
478 493 result[i] = p1[i] / p2[i];
479 494 }
  495 +
  496 + free(p1);
  497 + free(p2);
480 498 return true;
481 499 }
  500 +
  501 + //x * f(x)
  502 + bool x_area(double lb, double rb, double lab, double rab, T* result){
  503 + T* lp; //left band pointer
  504 + T* rp; //right band pointer
  505 + T* cur; //current band 1
  506 + T* cur2; //current band 2
  507 +
  508 + unsigned XY = R[0] * R[1];
  509 + unsigned S = XY * sizeof(T);
  510 +
  511 + lp = (T*) malloc(S); //memory allocation
  512 + rp = (T*) malloc(S);
  513 + cur = (T*) malloc(S);
  514 + cur2 = (T*) malloc(S);
  515 +
  516 + memset(result, (char)0, S);
  517 +
  518 + //find the wavelenght position in the whole band
  519 + unsigned int n = w.size();
  520 + unsigned int ai = 0; //left bound position
  521 + unsigned int bi = n - 1; //right bound position
  522 +
  523 + //to make sure the left and the right bound are in the bandwidth
  524 + if (lb < w[0] || rb < w[0] || lb > w[n-1] || rb >w[n-1]){
  525 + std::cout<<"ERROR: left bound or right bound out of bandwidth"<<std::endl;
  526 + exit(1);
  527 + }
  528 + //to make sure rigth bound is bigger than left bound
  529 + else if(lb > rb){
  530 + std::cout<<"ERROR: right bound should be bigger than left bound"<<std::endl;
  531 + exit(1);
  532 + }
  533 +
  534 + //get the position of lb and rb
  535 + while (lab >= w[ai]){
  536 + ai++;
  537 + }
  538 + while (rab <= w[bi]){
  539 + bi--;
  540 + }
  541 +
  542 + band(lp, lb);
  543 + band(rp, rb);
  544 +
  545 + //calculate the beginning and the ending part
  546 + baseline_band(lb, rb, lp, rp, rab, cur2); //ending part
  547 + baseline_band(lb, rb, lp, rp, w[bi], cur);
  548 + for(unsigned j = 0; j < XY; j++){
  549 + result[j] += (rab - w[bi]) * (rab + w[bi]) * (cur[j] + cur2[j]) / 4.0;
  550 + }
  551 + baseline_band(lb, rb, lp, rp, lab, cur2); //beginnning part
  552 + baseline_band(lb, rb, lp, rp, w[ai], cur);
  553 + for(unsigned j = 0; j < XY; j++){
  554 + result[j] += (w[ai] - lab) * (w[ai] + lab) * (cur[j] + cur2[j]) / 4.0;
  555 + }
  556 +
  557 + //calculate f(x) times x
  558 + ai++;
  559 + for(unsigned i = ai; i <= bi ;i++)
  560 + {
  561 + baseline_band(lb, rb, lp, rp, w[ai], cur2);
  562 + for(unsigned j = 0; j < XY; j++)
  563 + {
  564 + result[j] += (w[ai] - w[ai-1]) * (w[ai] + w[ai-1]) * (cur[j] + cur2[j]) / 4.0;
  565 + }
  566 + std::swap(cur,cur2); //swap the band pointers
  567 + }
  568 +
  569 + free(lp);
  570 + free(rp);
  571 + free(cur);
  572 + free(cur2);
  573 + return true;
  574 + }
  575 +
  576 + //centroid point
  577 + bool cpoint(double lb, double rb, double lab, double rab, T* result){
  578 + T* p1 = (T*)malloc(R[0] * R[1] * sizeof(T));
  579 + T* p2 = (T*)malloc(R[0] * R[1] * sizeof(T));
  580 +
  581 + //get the area and the peak band
  582 + x_area(lb, rb, lab, rab, p1);
  583 + area(lb, rb, lab, rab, p2);
  584 + //calculate the ratio in result
  585 + for(unsigned i = 0; i < R[0] * R[1]; i++){
  586 + if(p1[i] == 0 && p2[i] ==0)
  587 + result[i] = 1;
  588 + else
  589 + result[i] = p1[i] / p2[i];
  590 + }
  591 +
  592 + free(p1);
  593 + free(p2);
  594 + return true;
  595 + }
  596 +
482 597 //create mask file
483 598 bool mask(unsigned char* p, double mask_band, double threshold){
484 599  
... ... @@ -492,6 +607,7 @@ public:
492 607 p[i] = 255;
493 608 }
494 609  
  610 + free(temp);
495 611 return true;
496 612  
497 613 }
... ...
envi/envi.h
... ... @@ -419,6 +419,42 @@ public:
419 419 return false;
420 420 }
421 421  
  422 + //center of gravity
  423 + bool cpoint(double lb1, double rb1, double lab1, double rab1, void* result){
  424 + if(header.interleave == envi_header::BSQ){ //if the infile is bsq file
  425 + if(header.data_type ==envi_header::float32)
  426 + return ((bsq<float>*)file)->cpoint(lb1, rb1, lab1, rab1, (float*)result);
  427 + else if(header.data_type == envi_header::float64)
  428 + return ((bsq<double>*)file)->cpoint(lb1, rb1, lab1, rab1, (double*)result);
  429 + else
  430 + std::cout<<"ERROR: unidentified data type"<<std::endl;
  431 + }
  432 +
  433 + else if(header.interleave == envi_header::BIL){ //if the infile is bil file
  434 + if(header.data_type ==envi_header::float32)
  435 + return ((bil<float>*)file)->cpoint(lb1, rb1, lab1, rab1, (float*)result);
  436 + else if(header.data_type == envi_header::float64)
  437 + return ((bil<double>*)file)->cpoint(lb1, rb1, lab1, rab1, (double*)result);
  438 + else
  439 + std::cout<<"ERROR: unidentified data type"<<std::endl;
  440 + }
  441 +
  442 + else if(header.interleave == envi_header::BIP){ //if the infile is bip file
  443 + if(header.data_type ==envi_header::float32)
  444 + return ((bip<float>*)file)->cpoint(lb1, rb1, lab1, rab1,(float*)result);
  445 + else if(header.data_type == envi_header::float64)
  446 + return ((bip<double>*)file)->cpoint(lb1, rb1, lab1, rab1, (double*)result);
  447 + else
  448 + std::cout<<"ERROR: unidentified data type"<<std::endl;
  449 + }
  450 +
  451 + else{
  452 + std::cout<<"ERROR: unidentified file type"<<std::endl;
  453 + exit(1);
  454 + }
  455 + return false;
  456 + }
  457 +
422 458 bool close(){
423 459 if(header.interleave == envi_header::BSQ){
424 460 if(header.data_type ==envi_header::float32)
... ...