codebase / stimlib

Commit 9d77bbd9924fb66c38e7af7a3702bc286c909216

Authored by David Mayerich
1 parent b37665a9

updates for HSIview to support multiple files, enabled masking for convolution a… 

…nd derivative approximations
Showing 4 changed files with 70 additions and 57 deletions   Show diff stats
Inline Side-by-side
stim/envi/bil.h
  Show/Hide comments   View file @9d77bbd
... ... @@ -1182,7 +1182,7 @@ public:
1182 1182 /// @param nbands is the number of bands to process
1183 1183 /// @param center is the index for the center coefficient for the kernel (used to set the wavelengths in the output file)
1184 1184  
1185   - void convolve(std::string outfile, std::vector<double> C, size_t start, size_t end, bool PROGRESS = false){
  1185 + void convolve(std::string outfile, std::vector<double> C, size_t start, size_t end, unsigned char* mask = NULL, bool PROGRESS = false){
1186 1186 std::ofstream out(outfile.c_str(), std::ios::binary); //open the output file for writing
1187 1187 size_t B = end - start + 1;
1188 1188 size_t Xb = X() * sizeof(T); //calculate the size of a band (frame) in bytes
... ... @@ -1206,7 +1206,9 @@ public:
1206 1206 memset(outline, 0, Xb); //set the output band to zero (0)
1207 1207 for(size_t c = 0; c < N; c++){ //for each frame (corresponding to each coefficient)
1208 1208 for(size_t x = 0; x < X(); x++){ //for each pixel
1209   - outline[x] += (T)(C[c] * frame[c][x]); //calculate the contribution of the current frame (scaled by the corresponding coefficient)
  1209 + if(mask == NULL || mask[y * X() + x]){
  1210 + outline[x] += (T)(C[c] * frame[c][x]); //calculate the contribution of the current frame (scaled by the corresponding coefficient)
  1211 + }
1210 1212 }
1211 1213 }
1212 1214 out.write((char*)outline, Xb); //output the band
... ... @@ -1219,8 +1221,7 @@ public:
1219 1221 }
1220 1222  
1221 1223 /// Approximate the spectral derivative of the image
1222   -
1223   - void deriv(std::string outfile, size_t d, size_t order, const std::vector<double> w = std::vector<double>(), bool PROGRESS = false){
  1224 + void deriv(std::string outfile, size_t d, size_t order, const std::vector<double> w = std::vector<double>(), unsigned char* mask = NULL, bool PROGRESS = false){
1224 1225 std::ofstream out(outfile.c_str(), std::ios::binary); //open the output file for writing
1225 1226  
1226 1227 size_t Xb = X() * sizeof(T); //calculate the size of a line (frame) in bytes
... ... @@ -1264,7 +1265,9 @@ public:
1264 1265 memset(outline, 0, Xb); //set the output band to zero (0)
1265 1266 for(size_t c = 0; c < N; c++){ //for each frame (corresponding to each coefficient)
1266 1267 for(size_t x = 0; x < X(); x++){ //for each pixel
1267   - outline[x] += (T)(C[c] * frame[c][x]); //calculate the contribution of the current frame (scaled by the corresponding coefficient)
  1268 + if(mask == NULL || mask[y * X() + x]){
  1269 + outline[x] += (T)(C[c] * frame[c][x]); //calculate the contribution of the current frame (scaled by the corresponding coefficient)
  1270 + }
1268 1271 }
1269 1272 }
1270 1273 out.write((char*)outline, Xb); //output the band
... ...
stim/envi/bip.h
  Show/Hide comments   View file @9d77bbd
... ... @@ -1263,7 +1263,7 @@ public:
1263 1263 /// @param nbands is the number of bands to process
1264 1264 /// @param center is the index for the center coefficient for the kernel (used to set the wavelengths in the output file)
1265 1265  
1266   - void convolve(std::string outfile, std::vector<double> C, size_t start, size_t end, bool PROGRESS = false){
  1266 + void convolve(std::string outfile, std::vector<double> C, size_t start, size_t end, unsigned char* mask = NULL, bool PROGRESS = false){
1267 1267 std::ofstream out(outfile.c_str(), std::ios::binary); //open the output file for writing
1268 1268  
1269 1269 size_t N = end - start + 1; //number of bands in the output spectrum
... ... @@ -1281,9 +1281,11 @@ public:
1281 1281 for(size_t xy = 0; xy < XY; xy++){ //for each pixel
1282 1282 file.read((char*)inspec, Bb); //read an input spectrum
1283 1283 memset(outspec, 0, Nb); //set the output spectrum to zero (0)
1284   - for(size_t b = 0; b < N; b++){ //for each component of the spectrum
1285   - for(size_t c = 0; c < nC; c++){ //for each coefficient in the kernel
1286   - outspec[b] += (T)(inspec[b + start + c] * C[c]); //perform the sum/multiply part of the convolution
  1284 + if(mask == NULL || mask[xy]){
  1285 + for(size_t b = 0; b < N; b++){ //for each component of the spectrum
  1286 + for(size_t c = 0; c < nC; c++){ //for each coefficient in the kernel
  1287 + outspec[b] += (T)(inspec[b + start + c] * C[c]); //perform the sum/multiply part of the convolution
  1288 + }
1287 1289 }
1288 1290 }
1289 1291 out.write((char*)outspec, Nb); //output the band
... ... @@ -1291,7 +1293,7 @@ public:
1291 1293 }
1292 1294 }
1293 1295  
1294   - void deriv(std::string outfile, size_t d, size_t order, const std::vector<double> w, bool PROGRESS = false){
  1296 + void deriv(std::string outfile, size_t d, size_t order, const std::vector<double> w, unsigned char* mask = NULL, bool PROGRESS = false){
1295 1297 std::ofstream out(outfile.c_str(), std::ios::binary); //open the output file for writing
1296 1298  
1297 1299  
... ... @@ -1316,33 +1318,33 @@ public:
1316 1318 for(size_t xy = 0; xy < XY; xy++){ //for each pixel
1317 1319 file.read((char*)inspec, Bb); //read an input spectrum
1318 1320 memset(outspec, 0, Bb); //set the output spectrum to zero (0)
1319   -
1320   - iw = 0;
1321   - for(size_t b = 0; b < mid; b++){ //for each component of the spectrum
1322   - std::vector<double> w_pts(w.begin() + iw, w.begin() + iw + nC); //get the wavelengths corresponding to each sample
1323   - std::vector<double> C = diff_coefficients(w[b], w_pts, d); //get the optimal sample weights
1324   - for(size_t c = 0; c < nC; c++) //for each coefficient in the kernel
1325   - outspec[b] += (T)(inspec[iw + c] * C[c]); //perform the sum/multiply part of the convolution
1326   - }
1327   - std::vector<double> w_pts(w.begin(), w.begin() + nC); //get the wavelengths corresponding to each sample
1328   - std::vector<double> C = diff_coefficients(w[0], w_pts, d); //get the optimal sample weights
1329   - for(size_t b = mid; b <= B - (nC - mid); b++){
1330   - iw = b - mid;
1331   - if(!UNIFORM){ //if the spacing is non-uniform, we have to re-calculate these points every iteration
  1321 + if(mask == NULL || mask[xy]){
  1322 + iw = 0;
  1323 + for(size_t b = 0; b < mid; b++){ //for each component of the spectrum
1332 1324 std::vector<double> w_pts(w.begin() + iw, w.begin() + iw + nC); //get the wavelengths corresponding to each sample
1333 1325 std::vector<double> C = diff_coefficients(w[b], w_pts, d); //get the optimal sample weights
  1326 + for(size_t c = 0; c < nC; c++) //for each coefficient in the kernel
  1327 + outspec[b] += (T)(inspec[iw + c] * C[c]); //perform the sum/multiply part of the convolution
  1328 + }
  1329 + std::vector<double> w_pts(w.begin(), w.begin() + nC); //get the wavelengths corresponding to each sample
  1330 + std::vector<double> C = diff_coefficients(w[0], w_pts, d); //get the optimal sample weights
  1331 + for(size_t b = mid; b <= B - (nC - mid); b++){
  1332 + iw = b - mid;
  1333 + if(!UNIFORM){ //if the spacing is non-uniform, we have to re-calculate these points every iteration
  1334 + std::vector<double> w_pts(w.begin() + iw, w.begin() + iw + nC); //get the wavelengths corresponding to each sample
  1335 + std::vector<double> C = diff_coefficients(w[b], w_pts, d); //get the optimal sample weights
  1336 + }
  1337 + for(size_t c = 0; c < nC; c++) //for each coefficient in the kernel
  1338 + outspec[b] += (T)(inspec[iw + c] * C[c]); //perform the sum/multiply part of the convolution
  1339 + }
  1340 + iw = B - nC;
  1341 + for(size_t b = B - (nC - mid) + 1; b < B; b++){
  1342 + std::vector<double> w_pts(w.begin() + iw, w.begin() + iw + nC); //get the wavelengths corresponding to each sample
  1343 + std::vector<double> C = diff_coefficients(w[b], w_pts, d); //get the optimal sample weights
  1344 + for(size_t c = 0; c < nC; c++) //for each coefficient in the kernel
  1345 + outspec[b] += (T)(inspec[iw + c] * C[c]); //perform the sum/multiply part of the convolution
1334 1346 }
1335   - for(size_t c = 0; c < nC; c++) //for each coefficient in the kernel
1336   - outspec[b] += (T)(inspec[iw + c] * C[c]); //perform the sum/multiply part of the convolution
1337   - }
1338   - iw = B - nC;
1339   - for(size_t b = B - (nC - mid) + 1; b < B; b++){
1340   - std::vector<double> w_pts(w.begin() + iw, w.begin() + iw + nC); //get the wavelengths corresponding to each sample
1341   - std::vector<double> C = diff_coefficients(w[b], w_pts, d); //get the optimal sample weights
1342   - for(size_t c = 0; c < nC; c++) //for each coefficient in the kernel
1343   - outspec[b] += (T)(inspec[iw + c] * C[c]); //perform the sum/multiply part of the convolution
1344 1347 }
1345   -
1346 1348 out.write((char*)outspec, Bb); //output the band
1347 1349 if(PROGRESS) progress = (double)(xy+1) / (double)XY * 100;
1348 1350 }
... ...
stim/envi/bsq.h
  Show/Hide comments   View file @9d77bbd
... ... @@ -8,6 +8,7 @@
8 8 #include <utility>
9 9 #include <vector>
10 10 #include <deque>
  11 +#include <chrono>
11 12  
12 13  
13 14  
... ... @@ -1070,7 +1071,7 @@ public:
1070 1071 /// @param nbands is the number of bands to process
1071 1072 /// @param center is the index for the center coefficient for the kernel (used to set the wavelengths in the output file)
1072 1073  
1073   - void convolve(std::ofstream& out, std::vector<double> C, size_t start, size_t end, bool PROGRESS = false){
  1074 + void convolve(std::ofstream& out, std::vector<double> C, size_t start, size_t end, unsigned char* mask = NULL, bool PROGRESS = false){
1074 1075 size_t nbands = end - start + 1;
1075 1076 size_t XY = X() * Y(); //calculate the number of values in a band
1076 1077 size_t XYb = XY * sizeof(T); //calculate the size of a band (frame) in bytes
... ... @@ -1091,9 +1092,14 @@ public:
1091 1092 // re-inserted into the deque.
1092 1093 for(size_t b = 0; b < nbands; b++){ //for each band
1093 1094 memset(outband, 0, XYb); //set the output band to zero (0)
1094   - for(size_t c = 0; c < nframes; c++){ //for each frame (corresponding to each coefficient)
1095   - for(size_t xy = 0; xy < XY; xy++){ //for each pixel
1096   - outband[xy] += (T)(C[c] * frame[c][xy]); //calculate the contribution of the current frame (scaled by the corresponding coefficient)
  1095 + size_t c, xy;
  1096 + double coeff;
  1097 + for(c = 0; c < nframes; c++){ //for each frame (corresponding to each coefficient)
  1098 + coeff = C[c];
  1099 + for(xy = 0; xy < XY; xy++){ //for each pixel
  1100 + if(mask == NULL || mask[xy]){
  1101 + outband[xy] += (T)(coeff * frame[c][xy]); //calculate the contribution of the current frame (scaled by the corresponding coefficient)
  1102 + }
1097 1103 }
1098 1104 }
1099 1105 out.write((char*)outband, XYb); //output the band
... ... @@ -1110,9 +1116,9 @@ public:
1110 1116 /// @param C is an array of coefficients
1111 1117 /// @param start is the band to start processing (the first coefficient starts here)
1112 1118 /// @param nbands is the number of bands to process
1113   - void convolve(std::string outfile, std::vector<double> C, size_t start, size_t end, bool PROGRESS = false){
  1119 + void convolve(std::string outfile, std::vector<double> C, size_t start, size_t end, unsigned char* mask = NULL, bool PROGRESS = false){
1114 1120 std::ofstream out(outfile.c_str(), std::ios::binary); //open the output file for writing
1115   - convolve(out, C, start, end, PROGRESS); //start the convolution
  1121 + convolve(out, C, start, end, mask, PROGRESS); //start the convolution
1116 1122 out.close();
1117 1123 }
1118 1124  
... ... @@ -1122,21 +1128,21 @@ public:
1122 1128 /// @param C is an array containing an array of coefficients for each kernel
1123 1129 /// @param start is the list of start bands for each kernel
1124 1130 /// @param end is the list of end bands for each kernel
1125   - void convolve(std::string outfile, std::vector< std::vector<double> > C, std::vector<size_t> start, std::vector<size_t> end, bool PROGRESS = false){
  1131 + void convolve(std::string outfile, std::vector< std::vector<double> > C, std::vector<size_t> start, std::vector<size_t> end, unsigned char* mask = NULL, bool PROGRESS = false){
1126 1132 std::ofstream out(outfile.c_str(), std::ios::binary); //open the output file for writing
1127 1133  
1128 1134 size_t K = C.size(); //get the number of kernels
1129 1135 for(size_t k = 0; k < K; k++){
1130 1136 size_t b0 = start[k]; //calculate the range of the convolution
1131 1137 size_t b1 = end[k];
1132   - convolve(out, C[k], b0, b1, PROGRESS); //perform the convolution with the current kernel in the given range
  1138 + convolve(out, C[k], b0, b1, mask, PROGRESS); //perform the convolution with the current kernel in the given range
1133 1139 }
1134 1140 out.close();
1135 1141 }
1136 1142  
1137 1143 /// Approximate the spectral derivative of the image
1138 1144  
1139   - void deriv(std::string outfile, size_t d, size_t order, const std::vector<double> w = std::vector<double>(), bool PROGRESS = false){
  1145 + void deriv(std::string outfile, size_t d, size_t order, const std::vector<double> w = std::vector<double>(), unsigned char* mask = NULL, bool PROGRESS = false){
1140 1146 std::ofstream out(outfile.c_str(), std::ios::binary); //open the output file for writing
1141 1147  
1142 1148 size_t XY = X() * Y(); //calculate the number of values in a band
... ... @@ -1175,7 +1181,9 @@ public:
1175 1181 memset(outband, 0, XYb); //set the output band to zero (0)
1176 1182 for(size_t c = 0; c < N; c++){ //for each frame (corresponding to each coefficient)
1177 1183 for(size_t xy = 0; xy < XY; xy++){ //for each pixel
1178   - outband[xy] += (T)(C[c] * frame[c][xy]); //calculate the contribution of the current frame (scaled by the corresponding coefficient)
  1184 + if(mask == NULL || mask[xy]){
  1185 + outband[xy] += (T)(C[c] * frame[c][xy]); //calculate the contribution of the current frame (scaled by the corresponding coefficient)
  1186 + }
1179 1187 }
1180 1188 }
1181 1189 out.write((char*)outband, XYb); //output the band
... ...
stim/envi/envi.h
  Show/Hide comments   View file @9d77bbd
... ... @@ -1350,7 +1350,7 @@ public:
1350 1350 /// @param start is the band to start processing (the first coefficient starts here)
1351 1351 /// @param nbands is the number of bands to process
1352 1352 /// @param center is the index for the center coefficient for the kernel (used to set the wavelengths in the output file)
1353   - void convolve(std::string outfile, std::vector<double> C, size_t start, size_t end, size_t center = 0, bool PROGRESS = false){
  1353 + void convolve(std::string outfile, std::vector<double> C, size_t start, size_t end, size_t center = 0, unsigned char* mask = NULL, bool PROGRESS = false){
1354 1354 size_t nbands = end - start + 1;
1355 1355 envi_header h = header; //copy the current header
1356 1356 h.bands = nbands; //set the number of new bands
... ... @@ -1368,9 +1368,9 @@ public:
1368 1368  
1369 1369 if (header.interleave == envi_header::BSQ){
1370 1370 if (header.data_type == envi_header::float32)
1371   - ((bsq<float>*)file)->convolve(outfile, C, start, end, PROGRESS);
  1371 + ((bsq<float>*)file)->convolve(outfile, C, start, end, mask, PROGRESS);
1372 1372 else if (header.data_type == envi_header::float64)
1373   - ((bsq<double>*)file)->convolve(outfile, C, start, end, PROGRESS);
  1373 + ((bsq<double>*)file)->convolve(outfile, C, start, end, mask, PROGRESS);
1374 1374 else{
1375 1375 std::cout << "ERROR: unidentified data type" << std::endl;
1376 1376 exit(1);
... ... @@ -1378,9 +1378,9 @@ public:
1378 1378 }
1379 1379 else if (header.interleave == envi_header::BIL){
1380 1380 if (header.data_type == envi_header::float32)
1381   - ((bil<float>*)file)->convolve(outfile, C, start, end, PROGRESS);
  1381 + ((bil<float>*)file)->convolve(outfile, C, start, end, mask, PROGRESS);
1382 1382 else if (header.data_type == envi_header::float64)
1383   - ((bil<double>*)file)->convolve(outfile, C, start, end, PROGRESS);
  1383 + ((bil<double>*)file)->convolve(outfile, C, start, end, mask, PROGRESS);
1384 1384 else{
1385 1385 std::cout << "ERROR: unidentified data type" << std::endl;
1386 1386 exit(1);
... ... @@ -1388,9 +1388,9 @@ public:
1388 1388 }
1389 1389 else if (header.interleave == envi_header::BIP){
1390 1390 if (header.data_type == envi_header::float32)
1391   - ((bip<float>*)file)->convolve(outfile, C, start, end, PROGRESS);
  1391 + ((bip<float>*)file)->convolve(outfile, C, start, end, mask, PROGRESS);
1392 1392 else if (header.data_type == envi_header::float64)
1393   - ((bip<double>*)file)->convolve(outfile, C, start, end, PROGRESS);
  1393 + ((bip<double>*)file)->convolve(outfile, C, start, end, mask, PROGRESS);
1394 1394 else{
1395 1395 std::cout << "ERROR: unidentified data type" << std::endl;
1396 1396 exit(1);
... ... @@ -1403,13 +1403,13 @@ public:
1403 1403 /// @param outfile is the file where the derivative approximation will be saved
1404 1404 /// @n is the derivative to be calculated
1405 1405 /// @order is the order of the error (must be even)
1406   - void deriv(std::string outfile, size_t d, size_t order, bool PROGRESS = false){
  1406 + void deriv(std::string outfile, size_t d, size_t order, unsigned char* mask = NULL, bool PROGRESS = false){
1407 1407 header.save(outfile + ".hdr");
1408 1408 if (header.interleave == envi_header::BSQ){
1409 1409 if (header.data_type == envi_header::float32)
1410   - ((bsq<float>*)file)->deriv(outfile, d, order, header.wavelength, PROGRESS);
  1410 + ((bsq<float>*)file)->deriv(outfile, d, order, header.wavelength, mask, PROGRESS);
1411 1411 else if (header.data_type == envi_header::float64)
1412   - ((bsq<double>*)file)->deriv(outfile, d, order, header.wavelength, PROGRESS);
  1412 + ((bsq<double>*)file)->deriv(outfile, d, order, header.wavelength, mask, PROGRESS);
1413 1413 else{
1414 1414 std::cout << "ERROR: unidentified data type" << std::endl;
1415 1415 exit(1);
... ... @@ -1418,9 +1418,9 @@ public:
1418 1418  
1419 1419 else if (header.interleave == envi_header::BIL){
1420 1420 if (header.data_type == envi_header::float32)
1421   - ((bil<float>*)file)->deriv(outfile, d, order, header.wavelength, PROGRESS);
  1421 + ((bil<float>*)file)->deriv(outfile, d, order, header.wavelength, mask, PROGRESS);
1422 1422 else if (header.data_type == envi_header::float64)
1423   - ((bil<double>*)file)->deriv(outfile, d, order, header.wavelength, PROGRESS);
  1423 + ((bil<double>*)file)->deriv(outfile, d, order, header.wavelength, mask, PROGRESS);
1424 1424 else{
1425 1425 std::cout << "ERROR: unidentified data type" << std::endl;
1426 1426 exit(1);
... ... @@ -1429,9 +1429,9 @@ public:
1429 1429  
1430 1430 else if (header.interleave == envi_header::BIP){
1431 1431 if (header.data_type == envi_header::float32)
1432   - ((bip<float>*)file)->deriv(outfile, d, order, header.wavelength, PROGRESS);
  1432 + ((bip<float>*)file)->deriv(outfile, d, order, header.wavelength, mask, PROGRESS);
1433 1433 else if (header.data_type == envi_header::float64)
1434   - ((bip<double>*)file)->deriv(outfile, d, order, header.wavelength, PROGRESS);
  1434 + ((bip<double>*)file)->deriv(outfile, d, order, header.wavelength, mask, PROGRESS);
1435 1435 else{
1436 1436 std::cout << "ERROR: unidentified data type" << std::endl;
1437 1437 exit(1);
... ...