Changeset 31631

Timestamp:

Jun 15, 2011, 4:53:30 PM (15 years ago)

Author:

watersc1

Message:

disabled debug mode, and added test case for the number of inputs so that we can be sure that we're only running KMM when we have enough inputs to believe the results (N > 6). Ready for merging back into the trunk.

File:

• branches/czw_branch/20110406/psModules/src/imcombine/pmStack.c (modified) (19 diffs)

branches/czw_branch/20110406/psModules/src/imcombine/pmStack.c

@@ -37 +37 @@
 
 
-# if (1)
+# if (0)
 #define TESTING                         // Enable test output
 /* #define TEST_X 5745                       // x coordinate to examine */
@@ -106 +106 @@
 
 // KMM functions to do bimodality rejection of pixels
-
 float gaussian(float x, float m, float s) {
   return(pow(s * sqrt(2 * M_PI),-1) * exp(-0.5 * pow( (x - m) / s, 2)));
@@ -116 +115 @@
                          float *pi1, float *m1, float *s1,
                          float *pi2, float *m2, float *s2) {
+  assert(values);
+  assert(values->type.type == PS_TYPE_F32);
+  
   double logL_bimodal = 0, logL_unimodal;
   psVector *P1 = psVectorAlloc(values->n,PS_TYPE_F32);
@@ -121 +123 @@
   int i;
 /*   int debug = 0; */
-/*   if (fabs(values->data.F32[0] - -145.624954) < 1e-4) { */
-/*     debug = 1; */
-/*   } */
   
   // Calculate unimodal properties
@@ -156 +155 @@
   float w1,w2;
 
+  // These should be options.
   float KMM_TOLERANCE = 1e-6;
   int KMM_MAX_ITERATIONS = 500;
@@ -191 +191 @@
     *s1 = 0;
     *s2 = 0;
-    for (i = 0; i < values->n; i++) { // Calculate new means
+    for (i = 0; i < values->n; i++) { // Calculate new means and weights
       *m1 += values->data.F32[i] * P1->data.F32[i];
       *m2 += values->data.F32[i] * P2->data.F32[i];
@@ -226 +226 @@
   // Calculate Punimodal
   double lambda = -2.0 * (logL_unimodal - logL_bimodal);
-  int    df     = 2 + 2 * 1;
+  int    df     = 2 + 2 * 1; // I can't find my reference on this. 
   if (lambda > 0) {
     *Punimodal = gsl_cdf_chisq_Q(lambda,df);
   }
-  else {
+  else { // If lambda <= 0, then logL_unimodal > logL_bimodal, so Punimodal must be by definition 1.0
     *Punimodal = 1.0;
   }
@@ -238 +238 @@
 
 static void KMMFindPopular(const psVector *values, float *Punimodal, float *mean, float *sigma, float *pi) {
-  float KMM_MINIMUM_PVALUE = 0.05;
+  float KMM_MINIMUM_PVALUE = 0.05; // Should be an option.
   float mU,sU;
   float pi1,m1,s1,pi2,m2,s2;
   int iter;
 
+  assert(values);
+  assert(values->type.type == PS_TYPE_F32);
+  
   KMMcalculate(values,Punimodal,&iter,
                &mU,&sU,
@@ -255 +258 @@
   }
   else {
+    // Is bimodal. Select most popular mode.
     if (pi1 >= pi2) {
       *mean = m1;
@@ -455 +459 @@
 }
 #if (0)
+// Currently unused function to reject inputs for a given pixel(x,y) based on the selecting the
+// most popular mode after running the KMM test, and rejecting all inputs that belong to the
+// least popular mode.
 static void KMMRejectUnpopular(const psArray *inputs, int x, int y) {
   float KMM_MINIMUM_PVALUE = 0.05;
@@ -479 +486 @@
                &pi1,&m1,&s1,
                &pi2,&m2,&s2);
-  //  fprintf(stderr,
-/*   psTrace("psModules.imcombine",3, */
 
    CHECKPIX(x, y, 
@@ -516 +521 @@
 }
 
+// Currently unused function to reject inputs for a given pixel(x,y) based on the selecting the
+// faintest mode as determined by the KMM test, and rejecting all inputs that belong to the brighest.
 static void KMMRejectBright(const psArray *inputs, int x, int y) {
   float KMM_MINIMUM_PVALUE = 0.05;
@@ -558 +565 @@
   // else do nothing.
 }
-#endif
+#endif // End if(0) to prevent KMMReject{Unpopular|Bright} from being defined.
 
 // Extract vectors for simple combination/rejection operations
@@ -701 +708 @@
     if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
         for (int i = 0; i < numGood; i++) {
-          fprintf(stderr,"Input %d, pixel %d,%d (%" PRIu16 "): %f %f (%f) %f %f %d %x %x -> %x %x\n",
+            fprintf(stderr,"Input %d, pixel %d,%d (%" PRIu16 "): %f %f (%f) %f %f %d %x %x -> %x %x\n",
                     i, x, y, pixelSources->data.U16[i], pixelData->data.F32[i], pixelVariances->data.F32[i],
                     addVariance->data.F32[i], pixelWeights->data.F32[i], pixelExps->data.F32[i],
@@ -843 +850 @@
 
     // KMM values;
-    float Punimodal,KMMmean,KMMsigma,KMMpi;
+    float Punimodal = 1.0,KMMmean = NAN,KMMsigma = NAN,KMMpi = NAN;
+    int KMM_MINIMUM_INPUTS = 6;
+    bool useKMM = false;
+    if (num >= KMM_MINIMUM_INPUTS) {
+      useKMM = true;
+    }
     
     // Set up rejection limits
@@ -853 +865 @@
 
         // Determine the systematic error from the most popular population in the sample
-        KMMFindPopular(pixelData,&Punimodal,&KMMmean,&KMMsigma,&KMMpi);
-        CHECKPIX(x,y,"KMM Popularity Contest: (%d,%d) Puni: %g Mean: %f Sigma %f Pi: %f\n",
-                 x,y,Punimodal,KMMmean,KMMsigma,KMMpi);
-        
+        // This should probably be an option
+        if (useKMM) {
+          KMMFindPopular(pixelData,&Punimodal,&KMMmean,&KMMsigma,&KMMpi);
+          CHECKPIX(x,y,"KMM Popularity Contest: (%d,%d) Puni: %g Mean: %f Sigma %f Pi: %f\n",
+                   x,y,Punimodal,KMMmean,KMMsigma,KMMpi);
+        }
         for (int i = 0; i < num; i++) {
             sumWeights += pixelWeights->data.F32[i];
         }
         for (int i = 0; i < num; i++) {
-            // Systematic error contributes to the rejection level
-/*             float sysVar = PS_SQR(sys * pixelData->data.F32[i]); */
-          float sysVar = PS_SQR(KMMsigma);
+          // Systematic error contributes to the rejection level
+          float sysVar;
+          if (useKMM) { // If we can trust KMM results, set the systematic variance
+            sysVar = PS_SQR(KMMsigma);
+          }
+          else { // Otherwise, use the 10% systematic variance we've done in the past.
+            sysVar = PS_SQR(sys * pixelData->data.F32[i]);
+          }
+
             CHECKPIX(x, y, "Variance %d (%d), pixel %d,%d: %f %f %f\n", 
                      i, pixelSources->data.U16[i], x, y, 
@@ -1013 +1033 @@
               if (useVariance) {
                 float median;
-                if (Punimodal > 0.05) {
+                if ((useKMM)&&(Punimodal < 0.05)) {
+                  median = KMMmean;
+                }
+                else {
                   median = combinationWeightedOlympic(pixelData, pixelWeights,
                                                       olympic, buffer->sort); // Median for stack
                 }
-                else {
-                  median = KMMmean;
-                }
+                
                   CHECKPIX(x, y, "Flag with variance pixel %d,%d: median = %f\n", x, y, median);
                   float worst = -INFINITY; // Largest deviation
                   for (int j = 0; j < num; j++) {
-                    
-                      float diff = pixelData->data.F32[j] - median; // Difference from expected
+                      float diff = pixelData->data.F32[j] - median; // Difference from expected
                       CHECKPIX(x, y, "Testing input %d for pixel %d,%d: %f\n", j, x, y, diff);
 
@@ -1273 +1293 @@
 }
 
-  
-  
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -1445 +1463 @@
     for (int y = minInputRows; y < maxInputRows; y++) {
         for (int x = minInputCols; x < maxInputCols; x++) {
-
-/*          CHECKPIX(x, y, "Combining pixel %d,%d: %x %x %f %f %f %f %d %d %d\n", x, y, badMaskBits, blankMaskBits, iter, rej, sys, olympic, useVariance, safe, rejection); */
+            CHECKPIX(x, y, "Combining pixel %d,%d: %x %x %f %f %f %f %d %d %d\n", x, y, badMaskBits, blankMaskBits, iter, rej, sys, olympic, useVariance, safe, rejection);
 
 /*          // Pre-reject inputs using KMM bimodality test. */