implemented update from Daniel Scheffler

Author: davidfrantz

src/higher-level/spec-adjust-hl.c

@@ -41,9 +41,12 @@ int spectral_predict(ard_t ard, small *mask_, int nc, int sid);
 #define _SPECHOMO_N_DST_ 10 // number of destination bands
 #define _SPECHOMO_N_COF_  7 // number of coefficients
 #define _SPECHOMO_N_SIM_ 10 // max. number of close clusters to be used for kNN
-#define _SPECHOMO_GOOD_SAM_ 0.0698132 // accept clusters that are closer than this angle
-#define _SPECHOMO_POOR_SAM_ 0.2617995 // used to compute a weight for each cluster
+#define _SPECHOMO_GOOD_SAM_ 0.0698132 // good clusters are closer than this angle (4°)
+#define _SPECHOMO_MEDI_SAM_ 0.2094395 // ok'ish clusters are closer than this angle (12°)
+#define _SPECHOMO_POOR_SAM_ 0.2617995 // clusters further away from this angle (15°) should not be uses
+                                      // also used to compute a weight for each cluster
 #define _SPECHOMO_MIN_WEIGHT_ 1.0 - _SPECHOMO_GOOD_SAM_/_SPECHOMO_POOR_SAM_
+#define _SPECHOMO_MED_WEIGHT_ 1.0 - _SPECHOMO_MEDI_SAM_/_SPECHOMO_POOR_SAM_
 
 const char _SPECHOMO_SENSOR_[_SPECHOMO_N_SEN_][NPOW_04] = {
   "LND04", "LND05", "LND07", "LND08", "MOD01", "MOD02" };
@@ -736,7 +739,7 @@ double pred, wpred[_SPECHOMO_N_DST_], wsum;
         // weight for each cluster center
         // not exactly as in the paper, pragmatic suggestion by D. Scheffler
         if (sam > _SPECHOMO_POOR_SAM_){
-          weight[s] = 0.01;
+          weight[s] = 0.00001;
         } else {
           weight[s] = 1.0 - sam / _SPECHOMO_POOR_SAM_;
         }
@@ -750,27 +753,26 @@ double pred, wpred[_SPECHOMO_N_DST_], wsum;
         printf("lib: ");
         for (b=0; b<_SPECHOMO_N_SRC_; b++) printf("%05d ", _SPECHOMO_CENTER_[sid][b][s]);
         printf("\n");
-        printf("xx: %.2f, yy: %.2f, xy: %.2f, sam: %.2f, weight: %.2f\n", xx, yy, xy, sam, weight[s]);
+        printf("xx: %.2f, yy: %.2f, xy: %.2f, sam: %.5f, weight: %.5f\n", xx, yy, xy, sam, weight[s]);
         #endif
 
         // maximum weight of close clusters (-1 if no close cluster)
-        if (weight[s] >= _SPECHOMO_MIN_WEIGHT_ && 
-            weight[s] >   max_weight && 
+        if (weight[s] > max_weight && 
             s != _SPECHOMO_N_SIM_-1){
 
-            max_weight  = weight[s];
-            max_cluster = s;
+          max_weight  = weight[s];
+          max_cluster = s;
 
         }
 
         #ifdef FORCE_DEBUG
-        printf("max_weight: %.2f, max_cluster: %02d\n", max_weight, max_cluster);
+        printf("max_weight: %.5f, max_cluster: %02d\n", max_weight, max_cluster);
         #endif
 
       }
 
       #ifdef FORCE_DEBUG
-      print_dvector(weight,  "all weights", _SPECHOMO_N_CLS_, 2, 4);
+      print_dvector(weight,  "all weights", _SPECHOMO_N_CLS_, 2, 5);
       #endif
 
 
@@ -779,29 +781,30 @@ double pred, wpred[_SPECHOMO_N_DST_], wsum;
       // find closest clusters, fill with global
       for (c=0; c<_SPECHOMO_N_SIM_; c++){
 
-        // if no more close cluster can be found, add global cluster, then stop
-        if (max_weight < 0){
-          cluster_select[c] = _SPECHOMO_N_CLS_-1;
-          weight_select[c]  = weight[cluster_select[c]];
-          break;
-        }
+        // init with global cluster
+        cluster_select[c] = _SPECHOMO_N_CLS_-1;
+        weight_select[c]  = weight[_SPECHOMO_N_CLS_-1];
 
-        // copy closest cluster
-        cluster_select[c] = max_cluster;
-        weight_select[c]  = weight[cluster_select[c]];
+        if (max_weight > _SPECHOMO_MIN_WEIGHT_ ||
+           (max_weight > _SPECHOMO_MED_WEIGHT_ &&
+            max_weight > weight[_SPECHOMO_N_CLS_-1])){
+        
+          // copy closest cluster
+          cluster_select[c] = max_cluster;
+          weight_select[c]  = max_weight;
 
-        // remove weight from closest cluster
-        weight[cluster_select[c]] = 0.0;
+          // remove weight from closest cluster
+          weight[max_cluster] = 0.0;
+
+        } else break;
 
         // find the next closest cluster
         max_weight = -1.0;
         for (s=0; s<(_SPECHOMO_N_CLS_-1); s++){
-          if (weight[s] >= _SPECHOMO_MIN_WEIGHT_ && 
-              weight[s] >  max_weight){
+          if (weight[s] >  max_weight){
             max_weight  = weight[s];
             max_cluster = s;
           }
-
         }
 
       }
@@ -811,8 +814,8 @@ double pred, wpred[_SPECHOMO_N_DST_], wsum;
 
       #ifdef FORCE_DEBUG
       printf("found %d clusters\n", n_cluster);
-      print_ivector(cluster_select, "cluster", _SPECHOMO_N_SIM_, 7);
-      print_dvector(weight_select,  "weights", _SPECHOMO_N_SIM_, 2, 4);
+      print_ivector(cluster_select, "cluster", _SPECHOMO_N_SIM_, 8);
+      print_dvector(weight_select,  "weights", _SPECHOMO_N_SIM_, 2, 5);
       #endif
 
       // predict the target reflectance by using a weighted average of