codebase / stimlib

Commit 2aa68315740fc22b7b63d3b1a6629149e72ccd3e

Authored by David Mayerich
1 parent 43dec788

major bug fixes for PCA, covariance, and mean spectrum calculations, added Matla… 

…b code for loading entire ENVI files
Showing 5 changed files with 137 additions and 121 deletions   Show diff stats
Inline Side-by-side
matlab/enviLoadRaw.m 0 → 100644
  Show/Hide comments   View file @2aa6831
  1 +%loads an ENVI file without any manipulation (changing orientation)
  2 +function M = enviLoadRaw(filename, headername)
  3 +
  4 + %if a header isn't provided, assume it's just the filename
  5 + % with '.hdr' added to the end
  6 + if nargin == 1
  7 + headername = [filename '.hdr'];
  8 + end
  9 +
  10 + h = enviLoadHeader(headername);
  11 +
  12 + if strcmp(h.interleave, 'bsq')
  13 + X = h.samples;
  14 + Y = h.lines;
  15 + Z = h.bands;
  16 + elseif strcmp(h.interleave, 'bil')
  17 + X = h.samples;
  18 + Y = h.bands;
  19 + Z = h.lines;
  20 + else
  21 + X = h.bands;
  22 + Y = h.samples;
  23 + Z = h.lines;
  24 + end
  25 +
  26 + fid = fopen(filename);
  27 + M = fread(fid, [X, Y*Z], '*float32');
  28 + M = reshape(M, [X, Y, Z]);
  29 +
  30 + fclose(fid);
0 31 \ No newline at end of file
... ...
stim/envi/bil.h
  Show/Hide comments   View file @2aa6831
... ... @@ -923,7 +923,7 @@ public:
923 923  
924 924 /// @param p is a pointer to pre-allocated memory of size [B * sizeof(T)] that stores the mean spectrum
925 925 /// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
926   - bool avg_band(T*p, unsigned char* mask){
  926 + bool avg_band(double* p, unsigned char* mask = NULL){
927 927 unsigned long long XZ = X() * Z();
928 928 unsigned long long XY = X() * Y();
929 929 T* temp = (T*)malloc(sizeof(T) * XZ);
... ... @@ -931,19 +931,19 @@ public:
931 931 p[j] = 0;
932 932 }
933 933 //calculate vaild number in a band
934   - unsigned count = 0;
935   - for (unsigned j = 0; j < XY; j++){
936   - if (mask[j] != 0){
  934 + unsigned long long count = 0;
  935 + for (unsigned long long j = 0; j < XY; j++){
  936 + if (mask == NULL || mask[j] != 0){
937 937 count++;
938 938 }
939 939 }
940   - for (unsigned k = 0; k < Y(); k++){
  940 + for (unsigned long long k = 0; k < Y(); k++){
941 941 read_plane_y(temp, k);
942   - unsigned kx = k * X();
943   - for (unsigned i = 0; i < X(); i++){
944   - if (mask[kx + i] != 0){
945   - for (unsigned j = 0; j < Z(); j++){
946   - p[j] += temp[j * X() + i] / (T)count;
  942 + unsigned long long kx = k * X();
  943 + for (unsigned long long i = 0; i < X(); i++){
  944 + if (mask == NULL || mask[kx + i] != 0){
  945 + for (unsigned long long j = 0; j < Z(); j++){
  946 + p[j] += temp[j * X() + i] / (double)count;
947 947 }
948 948 }
949 949 }
... ... @@ -957,43 +957,48 @@ public:
957 957 /// @param co is a pointer to pre-allocated memory of size [B * B] that stores the resulting covariance matrix
958 958 /// @param avg is a pointer to memory of size B that stores the average spectrum
959 959 /// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
960   - bool co_matrix(T* co, T* avg, unsigned char *mask){
  960 + bool co_matrix(double* co, double* avg, unsigned char *mask){
  961 + thread_data = 0;
961 962 //memory allocation
962 963 unsigned long long xy = X() * Y();
963   - unsigned int B = Z();
  964 + unsigned long long B = Z();
964 965 T* temp = (T*)malloc(sizeof(T) * B);
965 966 //count vaild pixels in a band
966   - unsigned count = 0;
967   - for (unsigned j = 0; j < xy; j++){
968   - if (mask[j] != 0){
  967 + unsigned long long count = 0;
  968 + for (unsigned long long j = 0; j < xy; j++){
  969 + if (mask == NULL || mask[j] != 0){
969 970 count++;
970 971 }
971 972 }
972 973 //initialize correlation matrix
973   - for (unsigned i = 0; i < B; i++){
974   - for (unsigned k = 0; k < B; k++){
  974 + for (unsigned long long i = 0; i < B; i++){
  975 + for (unsigned long long k = 0; k < B; k++){
975 976 co[i * B + k] = 0;
976 977 }
977 978 }
978 979 //calculate correlation coefficient matrix
979   - for (unsigned j = 0; j < xy; j++){
980   - if (mask[j] != 0){
  980 + for (unsigned long long j = 0; j < xy; j++){
  981 + if (mask == NULL || mask[j] != 0){
981 982 pixel(temp, j);
982   - for (unsigned i = 0; i < B; i++){
983   - for (unsigned k = i; k < B; k++){
984   - co[i * B + k] += (temp[i] - avg[i]) * (temp[k] - avg[k]) / count;
  983 + for (unsigned long long i = 0; i < B; i++){
  984 + for (unsigned long long k = i; k < B; k++){
  985 + co[i * B + k] += ((double)temp[i] - (double)avg[i]) * ((double)temp[k] - (double)avg[k]) / (double)count;
985 986 }
986 987 }
987 988 }
  989 + thread_data = (double)j / xy * 100;
988 990 }
989 991 //because correlation matrix is symmetric
990   - for (unsigned i = 0; i < B; i++){
991   - for (unsigned k = i + 1; k < B; k++){
  992 + for (unsigned long long i = 0; i < B; i++){
  993 + for (unsigned long long k = i + 1; k < B; k++){
992 994 co[k * B + i] = co[i * B + k];
993 995 }
994 996 }
995 997  
996 998 free(temp);
  999 +
  1000 + thread_data = 100; //processing complete
  1001 +
997 1002 return true;
998 1003 }
999 1004  
... ...
stim/envi/bip.h
  Show/Hide comments   View file @2aa6831
... ... @@ -956,32 +956,11 @@ public:
956 956  
957 957 }
958 958  
959   -
960   -
961   - /// @param p is a pointer to memory of size X * Y * sizeof(T) that will store the band averages.
962   - bool band_avg(T* p){
963   - unsigned long long XY = X() * Y();
964   - //get every pixel and calculate average value
965   - T* temp = (T*)malloc(sizeof(T) * Z());
966   - T sum;
967   - for (unsigned i = 0; i < XY; i++){
968   - pixel(temp, i);
969   - //calculate the sum value of every value
970   - sum = 0; //initialize sum value
971   - for (unsigned j = 0; j < Z(); j++){
972   - sum += temp[j]/(T)Z();
973   - }
974   - p[i] = sum;
975   - }
976   - free(temp);
977   - return true;
978   - }
979   -
980 959 /// Calculate the mean value for all masked (or valid) pixels in a band and returns the average spectrum
981 960  
982 961 /// @param p is a pointer to pre-allocated memory of size [B * sizeof(T)] that stores the mean spectrum
983 962 /// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
984   - bool avg_band(T*p, unsigned char* mask){
  963 + bool avg_band(double* p, unsigned char* mask = NULL){
985 964 unsigned long long XY = X() * Y();
986 965 T* temp = (T*)malloc(sizeof(T) * Z());
987 966 //Iinitialize
... ... @@ -991,65 +970,76 @@ public:
991 970 //calculate vaild number in a band
992 971 unsigned count = 0;
993 972 for (unsigned j = 0; j < XY; j++){
994   - if (mask[j] != 0){
  973 + if (mask == NULL || mask[j] != 0){
995 974 count++;
996 975 }
997 976 }
998 977 //calculate average number of a band
999 978 for (unsigned i = 0; i < XY; i++){
1000   - if (mask[i] != 0){
  979 + if (mask == NULL || mask[i] != 0){
1001 980 pixel(temp, i);
1002 981 for (unsigned j = 0; j < Z(); j++){
1003   - p[j] += temp[j] / (T)count;
  982 + p[j] += (double)temp[j] / (double)count;
1004 983 }
1005 984 }
  985 + thread_data = (double)i / XY * 100;
1006 986 }
  987 +
1007 988 free(temp);
  989 + thread_data = 100;
1008 990 return true;
1009 991 }
1010 992  
1011   - /// Calculate the covariance matrix for all masked pixels in the image.
  993 + /// Calculate the covariance matrix for all masked pixels in the image with 64-bit floating point precision.
1012 994  
1013 995 /// @param co is a pointer to pre-allocated memory of size [B * B] that stores the resulting covariance matrix
1014 996 /// @param avg is a pointer to memory of size B that stores the average spectrum
1015 997 /// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
1016   - bool co_matrix(T* co, T* avg, unsigned char *mask){
  998 + bool co_matrix(double* co, double* avg, unsigned char *mask){
  999 + thread_data = 0;
1017 1000 //memory allocation
1018   - unsigned long long xy = X() * Y();
1019   - unsigned int B = Z();
  1001 + unsigned long long XY = X() * Y();
  1002 + unsigned long long B = Z();
1020 1003 T* temp = (T*)malloc(sizeof(T) * B);
1021 1004 //count vaild pixels in a band
1022   - unsigned count = 0;
1023   - for (unsigned j = 0; j < xy; j++){
1024   - if (mask[j] != 0){
  1005 + unsigned long long count = 0;
  1006 + for (unsigned long long j = 0; j < XY; j++){
  1007 + if (mask == NULL || mask[j] != 0){
1025 1008 count++;
1026 1009 }
1027 1010 }
1028   - //initialize correlation matrix
1029   - for (unsigned i = 0; i < B; i++){
1030   - for (unsigned k = 0; k < B; k++){
1031   - co[i * B + k] = 0;
1032   - }
1033   - }
1034   - //calculate correlation coefficient matrix
1035   - for (unsigned j = 0; j < xy; j++){
1036   - if (mask[j] != 0){
1037   - pixel(temp, j);
1038   - for (unsigned i = 0; i < B; i++){
1039   - for (unsigned k = i; k < B; k++){
1040   - co[i * B + k] += (temp[i] - avg[i]) * (temp[k] - avg[k]) / count;
1041   - }
  1011 + //initialize covariance matrix
  1012 + memset(co, 0, B * B * sizeof(double));
  1013 +
  1014 + //calculate covariance matrix
  1015 + T i_diff; //stores (temp[i] - avg[i]) to speed math
  1016 + double* co_half = (double*) malloc(B * B * sizeof(double)); //allocate space for a higher-precision intermediate matrix
  1017 + double* temp_precise = (double*) malloc(B * sizeof(double));
  1018 + memset(co_half, 0, B * B * sizeof(double)); //initialize the high-precision matrix with zeros
  1019 + unsigned long long idx; //stores i*B to speed indexing
  1020 + for (unsigned long long xy = 0; xy < XY; xy++){
  1021 + if (mask == NULL || mask[xy] != 0){
  1022 + pixel(temp, xy); //retreive the spectrum at the current xy pixel location
  1023 + for(unsigned long long b = 0; b < B; b++) //subtract the mean from this spectrum and increase the precision
  1024 + temp_precise[b] = (double)temp[b] - (double)avg[b];
  1025 + idx = 0;
  1026 + for (unsigned long long b0 = 0; b0 < B; b0++){ //for each band
  1027 + for (unsigned long long b1 = b0; b1 < B; b1++)
  1028 + co_half[idx++] += temp_precise[b0] * temp_precise[b1];
1042 1029 }
1043 1030 }
  1031 + thread_data = (double)xy / XY * 100;
1044 1032 }
1045   - //because correlation matrix is symmetric
1046   - for (unsigned i = 0; i < B; i++){
1047   - for (unsigned k = i + 1; k < B; k++){
1048   - co[k * B + i] = co[i * B + k];
  1033 + idx = 0;
  1034 + for (unsigned long long i = 0; i < B; i++){ //copy the precision matrix to both halves of the output matrix
  1035 + for (unsigned long long j = i; j < B; j++){
  1036 + co[j * B + i] = co[i * B + j] = co_half[idx++] / (double) count;
1049 1037 }
1050 1038 }
1051 1039  
  1040 + thread_data = 100; //processing complete
1052 1041 free(temp);
  1042 + free(temp_precise);
1053 1043 return true;
1054 1044 }
1055 1045  
... ...
stim/envi/bsq.h
  Show/Hide comments   View file @2aa6831
... ... @@ -157,7 +157,7 @@ public:
157 157 return false;
158 158 }
159 159  
160   - file.seekg(n * sizeof(T), std::ios::beg); //point to the certain pixel
  160 + file.seekg(n * sizeof(T) + header, std::ios::beg); //point to the certain pixel
161 161 for (unsigned i = 0; i < Z(); i++)
162 162 {
163 163 file.read((char *)(p + i), sizeof(T));
... ... @@ -966,24 +966,24 @@ public:
966 966  
967 967 /// @param p is a pointer to pre-allocated memory of size [B * sizeof(T)] that stores the mean spectrum
968 968 /// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
969   - bool avg_band(T*p, unsigned char* mask){
  969 + bool avg_band(double* p, unsigned char* mask = NULL){
970 970 unsigned long long XY = X() * Y();
971   - unsigned count = 0; //count will store the number of masked pixels
  971 + unsigned long long count = 0; //count will store the number of masked pixels
972 972 T* temp = (T*)malloc(sizeof(T) * XY);
973 973 //calculate this loop counts the number of true pixels in the mask
974 974 for (unsigned j = 0; j < XY; j++){
975   - if (mask[j] != 0){
  975 + if (mask == NULL || mask[j] != 0){
976 976 count++;
977 977 }
978 978 }
979 979 //this loops goes through each band in B (Z())
980 980 // masked (or valid) pixels from that band are averaged and the average is stored in p
981   - for (unsigned i = 0; i < Z(); i++){
  981 + for (unsigned long long i = 0; i < Z(); i++){
982 982 p[i] = 0;
983 983 band_index(temp, i); //get the band image and store it in temp
984   - for (unsigned j = 0; j < XY; j++){ //loop through temp, averaging valid pixels
985   - if (mask[j] != 0){
986   - p[i] += temp[j] / (T)count;
  984 + for (unsigned long long j = 0; j < XY; j++){ //loop through temp, averaging valid pixels
  985 + if (mask == NULL || mask[j] != 0){
  986 + p[i] += (double)temp[j] / (double)count;
987 987 }
988 988 }
989 989 }
... ... @@ -996,39 +996,44 @@ public:
996 996 /// @param co is a pointer to pre-allocated memory of size [B * B] that stores the resulting covariance matrix
997 997 /// @param avg is a pointer to memory of size B that stores the average spectrum
998 998 /// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
999   - bool co_matrix(T* co, T* avg, unsigned char *mask){
  999 + bool co_matrix(double* co, double* avg, unsigned char *mask){
  1000 + thread_data = 0; //initialize thread data for timers and progress bars
1000 1001 //memory allocation
1001 1002 unsigned long long xy = X() * Y();
1002   - unsigned int B = Z();
  1003 + unsigned long long B = Z();
1003 1004 T* bandi = (T*)malloc(sizeof(T) * xy);
1004 1005 T* bandj = (T*)malloc(sizeof(T) * xy);
1005 1006  
1006 1007 //count vaild pixels in a band
1007   - unsigned count = 0;
1008   - for (unsigned j = 0; j < xy; j++){
1009   - if (mask[j] != 0){
  1008 + unsigned long long count = 0;
  1009 + for (unsigned long long j = 0; j < xy; j++){
  1010 + if (mask == NULL || mask[j] != 0){
1010 1011 count++;
1011 1012 }
1012 1013 }
1013 1014 //calculate correlation coefficient matrix
1014   - for (unsigned i = 0; i < B; i++)
  1015 + for (unsigned long long i = 0; i < B; i++)
1015 1016 {
1016 1017 band_index(bandi, i);
1017   - for (unsigned j = i; j < B; j++){
  1018 + for (unsigned long long j = i; j < B; j++){
1018 1019 band_index(bandj, j);
1019   - T numerator = 0; //to calculate element in correlation coefficient matrix, numerator part
  1020 + double numerator = 0; //to calculate element in correlation coefficient matrix, numerator part
1020 1021 //calculate the R(i,j) in correlation coeffient matrix
1021   - for (unsigned k = 0; k < xy; k++){
1022   - if (mask[k] != 0){
1023   - numerator += (bandi[k] - avg[i]) * (bandj[k] - avg[j]) / count;
  1022 + for (unsigned long long k = 0; k < xy; k++){
  1023 + if (mask == NULL || mask[k] != 0){
  1024 + numerator += ((double)bandi[k] - (double)avg[i]) * ((double)bandj[k] - (double)avg[j]) / (double)count;
1024 1025 }
1025 1026 }
1026 1027 co[i*B + j] = numerator;
1027 1028 co[j*B + i] = numerator; //because correlated matrix is a symmetric matrix
1028 1029 }
  1030 + thread_data = (double)i / B * 100;
1029 1031 }
1030 1032 free(bandi);
1031 1033 free(bandj);
  1034 +
  1035 + thread_data = 100; //processing complete
  1036 +
1032 1037 return true;
1033 1038 }
1034 1039  
... ...
stim/envi/envi.h
  Show/Hide comments   View file @2aa6831
... ... @@ -963,30 +963,16 @@ public:
963 963 return false;
964 964 }
965 965  
966   - /// Helper function that loads a mask into memory given a filename.
967   -
968   - /// @param mask is a pointer to pre-allocated memory of size X*Y
969   - /// @param maskname is the file name for the image that will serve as the mask
970   - /*bool load_mask(unsigned char * mask, std::string maskname){
971   - //open the mask file
972   - stim::image<unsigned char> mask_image(maskname);
973   - mask_image.data_noninterleaved(mask);
974   - //save mask file into memory
975   - //memcpy(mask, mask_image.data_noninterleaved(), mask_image.size());
976   - //mask_image.clear();
977   - return true;
978   - }*/
979   -
980 966 /// Calculate the mean value for all masked (or valid) pixels in a band and returns the average spectrum
981 967  
982 968 /// @param p is a pointer to pre-allocated memory of size [B * sizeof(T)] that stores the mean spectrum
983 969 /// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
984   - bool avg_band(void * p, unsigned char* mask){
  970 + bool avg_band(double * p, unsigned char* mask){
985 971 if (header.interleave == envi_header::BSQ){
986 972 if (header.data_type == envi_header::float32)
987   - return ((bsq<float>*)file)->avg_band((float*)p, mask);
  973 + return ((bsq<float>*)file)->avg_band(p, mask);
988 974 else if (header.data_type == envi_header::float64)
989   - return ((bsq<double>*)file)->avg_band((double*)p, mask);
  975 + return ((bsq<double>*)file)->avg_band(p, mask);
990 976 else{
991 977 std::cout << "ERROR: unidentified data type" << std::endl;
992 978 exit(1);
... ... @@ -994,9 +980,9 @@ public:
994 980 }
995 981 else if (header.interleave == envi_header::BIL){
996 982 if (header.data_type == envi_header::float32)
997   - return ((bil<float>*)file)->avg_band((float*)p, mask);
  983 + return ((bil<float>*)file)->avg_band(p, mask);
998 984 else if (header.data_type == envi_header::float64)
999   - return ((bil<double>*)file)->avg_band((double*)p, mask);
  985 + return ((bil<double>*)file)->avg_band(p, mask);
1000 986 else{
1001 987 std::cout << "ERROR: unidentified data type" << std::endl;
1002 988 exit(1);
... ... @@ -1004,9 +990,9 @@ public:
1004 990 }
1005 991 else if (header.interleave == envi_header::BIP){
1006 992 if (header.data_type == envi_header::float32)
1007   - return ((bip<float>*)file)->avg_band((float*)p, mask);
  993 + return ((bip<float>*)file)->avg_band(p, mask);
1008 994 else if (header.data_type == envi_header::float64)
1009   - return ((bip<double>*)file)->avg_band((double*)p, mask);
  995 + return ((bip<double>*)file)->avg_band(p, mask);
1010 996 else{
1011 997 std::cout << "ERROR: unidentified data type" << std::endl;
1012 998 exit(1);
... ... @@ -1020,12 +1006,12 @@ public:
1020 1006 /// @param co is a pointer to pre-allocated memory of size [B * B] that stores the resulting covariance matrix
1021 1007 /// @param avg is a pointer to memory of size B that stores the average spectrum
1022 1008 /// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
1023   - bool co_matrix(void* co, void* avg, unsigned char* mask){
  1009 + bool co_matrix(double* co, double* avg, unsigned char* mask){
1024 1010 if (header.interleave == envi_header::BSQ){
1025 1011 if (header.data_type == envi_header::float32)
1026   - return ((bsq<float>*)file)->co_matrix((float*)co, (float*)avg, mask);
  1012 + return ((bsq<float>*)file)->co_matrix(co, avg, mask);
1027 1013 else if (header.data_type == envi_header::float64)
1028   - return ((bsq<double>*)file)->co_matrix((double*)co, (double*)avg, mask);
  1014 + return ((bsq<double>*)file)->co_matrix(co, avg, mask);
1029 1015 else{
1030 1016 std::cout << "ERROR: unidentified data type" << std::endl;
1031 1017 exit(1);
... ... @@ -1033,9 +1019,9 @@ public:
1033 1019 }
1034 1020 else if (header.interleave == envi_header::BIL){
1035 1021 if (header.data_type == envi_header::float32)
1036   - return ((bil<float>*)file)->co_matrix((float*)co, (float*)avg, mask);
  1022 + return ((bil<float>*)file)->co_matrix(co, avg, mask);
1037 1023 else if (header.data_type == envi_header::float64)
1038   - return ((bil<double>*)file)->co_matrix((double*)co, (double*)avg, mask);
  1024 + return ((bil<double>*)file)->co_matrix(co, avg, mask);
1039 1025 else{
1040 1026 std::cout << "ERROR: unidentified data type" << std::endl;
1041 1027 exit(1);
... ... @@ -1043,9 +1029,9 @@ public:
1043 1029 }
1044 1030 else if (header.interleave == envi_header::BIP){
1045 1031 if (header.data_type == envi_header::float32)
1046   - return ((bip<float>*)file)->co_matrix((float*)co, (float*)avg, mask);
  1032 + return ((bip<float>*)file)->co_matrix(co, avg, mask);
1047 1033 else if (header.data_type == envi_header::float64)
1048   - return ((bip<double>*)file)->co_matrix((double*)co, (double*)avg, mask);
  1034 + return ((bip<double>*)file)->co_matrix(co, avg, mask);
1049 1035 else{
1050 1036 std::cout << "ERROR: unidentified data type" << std::endl;
1051 1037 exit(1);
... ...