codebase / stimlib

Commit e28fe60df2f39d68dd878fc854d5a69418e4c88c

Authored by heziqi
1 parent 1cba0efe

apply mask in correlation matrix

Showing 4 changed files with 191 additions and 125 deletions   Show diff stats
Inline Side-by-side
envi/bil.h
  Show/Hide comments   View file @e28fe60
... ... @@ -755,19 +755,28 @@ public:
755 755 }
756 756  
757 757 //calculate the average number of every band
758   - bool avg_band(T*p){
  758 + bool avg_band(T*p, unsigned char* mask){
759 759 unsigned long long XZ = R[0] * R[2];
760 760 unsigned long long XY = R[0] * R[1];
761 761 T* temp = (T*)malloc(sizeof(T) * XZ);
762 762 for (unsigned j = 0; j < R[2]; j++){
763 763 p[j] = 0;
764 764 }
  765 + //calculate vaild number in a band
  766 + unsigned count = 0;
  767 + for (unsigned j = 0; j < XY; j++){
  768 + if (mask[j] != 0){
  769 + count++;
  770 + }
  771 + }
765 772 for (unsigned k = 0; k < R[1]; k++){
766 773 getY(temp, k);
767   - for (unsigned j = 0; j < R[2]; j++){
768   - unsigned jx = j * R[0];
769   - for (unsigned i = 0; i < R[0]; i++){
770   - p[j] += temp[jx + i] / (T)XY;
  774 + unsigned kx = k * R[0];
  775 + for (unsigned i = 0; i < R[0]; i++){
  776 + if (mask[kx + i] != 0){
  777 + for (unsigned j = 0; j < R[2]; j++){
  778 + p[j] += temp[j * R[0] + i] / (T)count;
  779 + }
771 780 }
772 781 }
773 782 }
... ... @@ -775,6 +784,46 @@ public:
775 784 return true;
776 785 }
777 786  
  787 + //calculate correlation coefficient matrix
  788 + bool co_matrix(T* co, T* avg, unsigned char *mask){
  789 + //memory allocation
  790 + unsigned long long xy = R[0] * R[1];
  791 + unsigned int B = R[2];
  792 + T* temp = (T*)malloc(sizeof(T) * B);
  793 + //count vaild pixels in a band
  794 + unsigned count = 0;
  795 + for (unsigned j = 0; j < xy; j++){
  796 + if (mask[j] != 0){
  797 + count++;
  798 + }
  799 + }
  800 + //initialize correlation matrix
  801 + for (unsigned i = 0; i < B; i++){
  802 + for (unsigned k = 0; k < B; k++){
  803 + co[i * B + k] = 0;
  804 + }
  805 + }
  806 + //calculate correlation coefficient matrix
  807 + for (unsigned j = 0; j < xy; j++){
  808 + if (mask[j] != 0){
  809 + pixel(temp, j);
  810 + for (unsigned i = 0; i < B; i++){
  811 + for (unsigned k = i; k < B; k++){
  812 + co[i * B + k] += (temp[i] - avg[i]) * (temp[k] - avg[k]) / count;
  813 + }
  814 + }
  815 + }
  816 + }
  817 + //because correlation matrix is symmetric
  818 + for (unsigned i = 0; i < B; i++){
  819 + for (unsigned k = i + 1; k < B; k++){
  820 + co[k * B + i] = co[i * B + k];
  821 + }
  822 + }
  823 +
  824 + free(temp);
  825 + return true;
  826 + }
778 827  
779 828 //close the file
780 829 bool close(){
... ...
envi/bip.h
  Show/Hide comments   View file @e28fe60
... ... @@ -841,24 +841,72 @@ public:
841 841 }
842 842  
843 843 //calculate the average number of every band
844   - bool avg_band(T*p){
  844 + bool avg_band(T*p, unsigned char* mask){
845 845 unsigned long long XY = R[0] * R[1];
846 846 T* temp = (T*)malloc(sizeof(T) * R[2]);
847 847 //Iinitialize
848 848 for (unsigned j = 0; j < R[2]; j++){
849 849 p[j] = 0;
850 850 }
851   -
  851 + //calculate vaild number in a band
  852 + unsigned count = 0;
  853 + for (unsigned j = 0; j < XY; j++){
  854 + if (mask[j] != 0){
  855 + count++;
  856 + }
  857 + }
  858 + //calculate average number of a band
852 859 for (unsigned i = 0; i < XY; i++){
853   - pixel(temp, i);
854   - for (unsigned j = 0; j < R[2]; j++){
855   - p[j] += temp[j] / (T)XY;
  860 + if (mask[i] != 0){
  861 + pixel(temp, i);
  862 + for (unsigned j = 0; j < R[2]; j++){
  863 + p[j] += temp[j] / (T)count;
  864 + }
856 865 }
857 866 }
858 867 free(temp);
859 868 return true;
860 869 }
  870 + //calculate correlation coefficient matrix
  871 + bool co_matrix(T* co, T* avg, unsigned char *mask){
  872 + //memory allocation
  873 + unsigned long long xy = R[0] * R[1];
  874 + unsigned int B = R[2];
  875 + T* temp = (T*)malloc(sizeof(T) * B);
  876 + //count vaild pixels in a band
  877 + unsigned count = 0;
  878 + for (unsigned j = 0; j < xy; j++){
  879 + if (mask[j] != 0){
  880 + count++;
  881 + }
  882 + }
  883 + //initialize correlation matrix
  884 + for (unsigned i = 0; i < B; i++){
  885 + for (unsigned k = 0; k < B; k++){
  886 + co[i * B + k] = 0;
  887 + }
  888 + }
  889 + //calculate correlation coefficient matrix
  890 + for (unsigned j = 0; j < xy; j++){
  891 + if (mask[j] != 0){
  892 + pixel(temp, j);
  893 + for (unsigned i = 0; i < B; i++){
  894 + for (unsigned k = i; k < B; k++){
  895 + co[i * B + k] += (temp[i] - avg[i]) * (temp[k] - avg[k]) / count;
  896 + }
  897 + }
  898 + }
  899 + }
  900 + //because correlation matrix is symmetric
  901 + for (unsigned i = 0; i < B; i++){
  902 + for (unsigned k = i + 1; k < B; k++){
  903 + co[k * B + i] = co[i * B + k];
  904 + }
  905 + }
861 906  
  907 + free(temp);
  908 + return true;
  909 + }
862 910  
863 911 //close the file
864 912 bool close(){
... ...
envi/bsq.h
  Show/Hide comments   View file @e28fe60
... ... @@ -684,20 +684,67 @@ public:
684 684 }
685 685  
686 686 //calculate the average number of every band
687   - bool avg_band(T*p){
  687 + bool avg_band(T*p, unsigned char* mask){
688 688 unsigned long long XY = R[0] * R[1];
  689 + unsigned count = 0; //number of vaild pixel in a band
689 690 T* temp = (T*)malloc(sizeof(T) * XY);
  691 + //calculate valid pixel number
  692 + for (unsigned j = 0; j < XY; j++){
  693 + if (mask[j] != 0){
  694 + count++;
  695 + }
  696 + }
  697 + //calculate average of a band
690 698 for (unsigned i = 0; i < R[2]; i++){
691 699 p[i] = 0;
692 700 band_index(temp, i);
693 701 for (unsigned j = 0; j < XY; j++){
694   - p[i] += temp[j] / (T)XY;
  702 + if (mask[j] != 0){
  703 + p[i] += temp[j] / (T)count;
  704 + }
695 705 }
696 706 }
697 707 free(temp);
698 708 return true;
699 709 }
700 710  
  711 + //calculate correlated matrix
  712 + bool co_matrix(T* co, T* avg, unsigned char *mask){
  713 + //memory allocation
  714 + unsigned long long xy = R[0] * R[1];
  715 + unsigned int B = R[2];
  716 + T* bandi = (T*)malloc(sizeof(T) * xy);
  717 + T* bandj = (T*)malloc(sizeof(T) * xy);
  718 +
  719 + //count vaild pixels in a band
  720 + unsigned count = 0;
  721 + for (unsigned j = 0; j < xy; j++){
  722 + if (mask[j] != 0){
  723 + count++;
  724 + }
  725 + }
  726 + //calculate correlation coefficient matrix
  727 + for (unsigned i = 0; i < B; i++)
  728 + {
  729 + band_index(bandi, i);
  730 + for (unsigned j = i; j < B; j++){
  731 + band_index(bandj, j);
  732 + T numerator = 0; //to calculate element in correlation coefficient matrix, numerator part
  733 + //calculate the R(i,j) in correlation coeffient matrix
  734 + for (unsigned k = 0; k < xy; k++){
  735 + if (mask[k] != 0){
  736 + numerator += (bandi[k] - avg[i]) * (bandj[k] - avg[j]) / count;
  737 + }
  738 + }
  739 + co[i*B + j] = numerator;
  740 + co[j*B + i] = numerator; //because correlated matrix is a symmetric matrix
  741 + }
  742 + }
  743 + free(bandi);
  744 + free(bandj);
  745 + return true;
  746 + }
  747 +
701 748 //close the file
702 749 bool close(){
703 750 file.close();
... ...
envi/envi.h
  Show/Hide comments   View file @e28fe60
... ... @@ -643,12 +643,12 @@ public:
643 643 }
644 644  
645 645 //calculate band average
646   - bool avg_band(void * p){
  646 + bool avg_band(void * p, unsigned char* mask){
647 647 if (header.interleave == envi_header::BSQ){
648 648 if (header.data_type == envi_header::float32)
649   - return ((bsq<float>*)file)->avg_band((float*)p);
  649 + return ((bsq<float>*)file)->avg_band((float*)p, mask);
650 650 else if (header.data_type == envi_header::float64)
651   - return ((bsq<double>*)file)->avg_band((double*)p);
  651 + return ((bsq<double>*)file)->avg_band((double*)p, mask);
652 652 else{
653 653 std::cout << "ERROR: unidentified data type" << std::endl;
654 654 exit(1);
... ... @@ -656,9 +656,9 @@ public:
656 656 }
657 657 else if (header.interleave == envi_header::BIL){
658 658 if (header.data_type == envi_header::float32)
659   - return ((bil<float>*)file)->avg_band((float*)p);
  659 + return ((bil<float>*)file)->avg_band((float*)p, mask);
660 660 else if (header.data_type == envi_header::float64)
661   - return ((bil<double>*)file)->avg_band((double*)p);
  661 + return ((bil<double>*)file)->avg_band((double*)p, mask);
662 662 else{
663 663 std::cout << "ERROR: unidentified data type" << std::endl;
664 664 exit(1);
... ... @@ -666,9 +666,9 @@ public:
666 666 }
667 667 else if (header.interleave == envi_header::BIP){
668 668 if (header.data_type == envi_header::float32)
669   - return ((bip<float>*)file)->avg_band((float*)p);
  669 + return ((bip<float>*)file)->avg_band((float*)p, mask);
670 670 else if (header.data_type == envi_header::float64)
671   - return ((bip<double>*)file)->avg_band((double*)p);
  671 + return ((bip<double>*)file)->avg_band((double*)p, mask);
672 672 else{
673 673 std::cout << "ERROR: unidentified data type" << std::endl;
674 674 exit(1);
... ... @@ -678,116 +678,38 @@ public:
678 678 }
679 679  
680 680 //calculate correlation coefficient matrix with mask
681   - bool co_matrix(void* co, void* avg, unsigned mask){
682   - mask = 1; //delete it later
683   - if (header.data_type == envi_header::float32){
684   - //memory allocation
685   - unsigned long long xy = header.samples * header.lines;
686   - unsigned int b = header.bands;
687   - float* bandi = (float*)malloc(sizeof(float) * xy);
688   - float* bandj = (float*)malloc(sizeof(float) * xy);
689   - float* var = (float*)malloc(sizeof(float) * b); //save variance of each band
690   - //calculate the first line of correlation coefficient matrix and initialize
691   - band_index((void*)bandi, 0);
692   - for (unsigned j = 0; j < header.bands; j++){
693   - band_index((void*)bandj, j);
694   - float numerator = 0; //to calculate element in correlation coefficient matrix, numerator part
695   - float denominator = 0; //denominator part
696   - float denominator1 = 0; //part 1 of the denominator
697   - float denominator2 = 0; //part 2 of the denominator
698   - //calculate the R(i,j) in correlation coefficient matrix
699   - for (unsigned k = 0; k < xy; k++){
700   - if (mask != 0){
701   - numerator += (bandi[k] - ((float*)avg)[0]) * (bandj[k] - ((float*)avg)[j]);
702   - denominator1 += (bandi[k] - ((float*)avg)[0]) * (bandi[k] - ((float*)avg)[0]);
703   - denominator2 += (bandj[k] - ((float*)avg)[j]) * (bandj[k] - ((float*)avg)[j]);
704   - }
705   - }
706   - var[j] = denominator2;
707   - denominator = sqrt(denominator1 * denominator2); //calculate denominator part
708   - ((float*)co)[j] = numerator / denominator;
709   - ((float*)co)[j*header.bands] = ((float*)co)[j]; //because correlated matrix is a symmetric matrix
710   - }
711   -
712   - //calculate the other part of correlated matrix
713   - for (unsigned i = 1; i < header.bands; i++)
714   - {
715   - band_index((void*)bandi, i);
716   - for (unsigned j = i; j < header.bands; j++){
717   - band_index((void*)bandj, j);
718   - float numerator = 0; //to calculate element in correlation coefficient matrix, numerator part
719   - float denominator = 0; //denominator part
720   - //calculate the R(i,j) in correlation coeffient matrix
721   - for (unsigned k = 0; k < xy; k++){
722   - if (mask != 0){
723   - numerator += (bandi[k] - ((float*)avg)[i]) * (bandj[k] - ((float*)avg)[j]);
724   - }
725   - }
726   - denominator = sqrt(var[i] * var[j]); //calculate denominator part
727   - ((float*)co)[i*header.bands + j] = numerator / denominator;
728   - ((float*)co)[j*header.bands + i] = ((float*)co)[i*header.bands + j]; //because correlated matrix is a symmetric matrix
729   - }
730   - }
731   - free(bandi);
732   - free(bandj);
733   - }
734   - //if data type of the binary file is double
735   - else if (header.data_type == envi_header::float64){
736   - //memory allocation
737   - unsigned long long xy = header.samples * header.lines;
738   - unsigned int b = header.bands;
739   - double* bandi = (double*)malloc(sizeof(double) * xy);
740   - double* bandj = (double*)malloc(sizeof(double) * xy);
741   - double* var = (double*)malloc(sizeof(double) * b); //save variance of each band
742   - //calculate the first line of correlation coefficient matrix and initialize
743   - band_index((void*)bandi, 0);
744   - for (unsigned j = 0; j < header.bands; j++){
745   - band_index((void*)bandj, j);
746   - double numerator = 0; //to calculate element in correlation coefficient matrix, numerator part
747   - double denominator = 0; //denominator part
748   - double denominator1 = 0; //part 1 of the denominator
749   - double denominator2 = 0; //part 2 of the denominator
750   - //calculate the R(i,j) in correlation coefficient matrix
751   - for (unsigned k = 0; k < xy; k++){
752   - if (mask != 0){
753   - numerator += (bandi[k] - ((double*)avg)[0]) * (bandj[k] - ((double*)avg)[j]);
754   - denominator1 += (bandi[k] - ((double*)avg)[0]) * (bandi[k] - ((double*)avg)[0]);
755   - denominator2 += (bandj[k] - ((double*)avg)[j]) * (bandj[k] - ((double*)avg)[j]);
756   - }
757   - }
758   - var[j] = denominator2;
759   - denominator = sqrt(denominator1 * denominator2); //calculate denominator part
760   - ((double*)co)[j] = numerator / denominator;
761   - ((double*)co)[j*header.bands] = ((double*)co)[j]; //because correlated matrix is a symmetric matrix
  681 + bool co_matrix(void* co, void* avg, unsigned char* mask){
  682 + if (header.interleave == envi_header::BSQ){
  683 + if (header.data_type == envi_header::float32)
  684 + return ((bsq<float>*)file)->co_matrix((float*)co, (float*)avg, mask);
  685 + else if (header.data_type == envi_header::float64)
  686 + return ((bsq<double>*)file)->co_matrix((double*)co, (double*)avg, mask);
  687 + else{
  688 + std::cout << "ERROR: unidentified data type" << std::endl;
  689 + exit(1);
762 690 }
763   -
764   - //calculate the other part of correlated matrix
765   - for (unsigned i = 1; i < header.bands; i++)
766   - {
767   - band_index((void*)bandi, i);
768   - for (unsigned j = i; j < header.bands; j++){
769   - band_index((void*)bandj, j);
770   - double numerator = 0; //to calculate element in correlation coefficient matrix, numerator part
771   - double denominator = 0; //denominator part
772   - //calculate the R(i,j) in correlation coeffient matrix
773   - for (unsigned k = 0; k < xy; k++){
774   - if (mask != 0){
775   - numerator += (bandi[k] - ((double*)avg)[i]) * (bandj[k] - ((double*)avg)[j]);
776   - }
777   - }
778   - denominator = sqrt(var[i] * var[j]); //calculate denominator part
779   - ((double*)co)[i*header.bands + j] = numerator / denominator;
780   - ((double*)co)[j*header.bands + i] = ((double*)co)[i*header.bands + j]; //because correlated matrix is a symmetric matrix
781   - }
  691 + }
  692 + else if (header.interleave == envi_header::BIL){
  693 + if (header.data_type == envi_header::float32)
  694 + return ((bil<float>*)file)->co_matrix((float*)co, (float*)avg, mask);
  695 + else if (header.data_type == envi_header::float64)
  696 + return ((bil<double>*)file)->co_matrix((double*)co, (double*)avg, mask);
  697 + else{
  698 + std::cout << "ERROR: unidentified data type" << std::endl;
  699 + exit(1);
782 700 }
783   - free(bandi);
784   - free(bandj);
785 701 }
786   - else{
787   - std::cout << "ERROR: unidentified data type" << std::endl;
788   - exit(1);
  702 + else if (header.interleave == envi_header::BIP){
  703 + if (header.data_type == envi_header::float32)
  704 + return ((bip<float>*)file)->co_matrix((float*)co, (float*)avg, mask);
  705 + else if (header.data_type == envi_header::float64)
  706 + return ((bip<double>*)file)->co_matrix((double*)co, (double*)avg, mask);
  707 + else{
  708 + std::cout << "ERROR: unidentified data type" << std::endl;
  709 + exit(1);
  710 + }
789 711 }
790   - return true;
  712 + return false;
791 713 }
792 714 };
793 715  
... ...