Changeset 31618 for branches/czw_branch/20110406/psModules/src/imcombine/pmStack.c

Timestamp:

Jun 8, 2011, 2:41:13 PM (15 years ago)

Author:

watersc1

Message:

compilable version of updated stack code (doesn't correctly reject) and asinh scaling (completely untested).

File:

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

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

@@ -22 +22 @@
 #include <pslib.h>
 
+#include <gsl/gsl_cdf.h>
+
 #include "pmHDU.h"
 #include "pmFPA.h"
@@ -35 +37 @@
 
 
-# if (0)
+# if (1)
 #define TESTING                         // Enable test output
-#define TEST_X 3145                       // x coordinate to examine
-#define TEST_Y 2334                       // y coordinate to examine
+#define TEST_X 5745                       // x coordinate to examine
+#define TEST_Y 5331                       // y coordinate to examine
 #define TEST_RADIUS 0.5                 // Radius to examine
 # endif
@@ -332 +334 @@
     memset (nGoodBits, 0, 16*sizeof(int));
 
-    // Extract the pixel and mask data
+
+    // Extract the pixel and mask data    
     int numGood = 0;                    // Number of good pixels
     for (int i = 0, j = 0; i < inputs->n; i++) {
@@ -364 +367 @@
             continue;
         }
-
+        
         int xIn = x - data->readout->col0, yIn = y - data->readout->row0; // Coordinates on input readout
         psImage *mask = data->readout->mask; // Mask of interest
@@ -984 +987 @@
 
 static void KMMcalculate(const psVector *values,
-                         float *Punimodal,
+                         float *Punimodal,int *iter,
                          float *pi1, float *m1, float *s1,
                          float *pi2, float *m2, float *s2) {
@@ -1012 +1015 @@
   float dL = 0;
   float oldL = -999;
-  int k = 0;
+  *iter = 0;
   logL_bimodal = logL_unimodal;
   *m1 = mU - 3 * sU;
@@ -1024 +1027 @@
   float w1,w2;
 
-  while (((dL > KMM_TOLERANCE)||(k < 3))&&(k < KMM_MAX_ITERATIONS)) {
-    k++;
+  float KMM_TOLERANCE = 1e-6;
+  int KMM_MAX_ITERATIONS = 500;
+  float KMM_SMALL_NUMBER = 1e-5;
+  while (((dL > KMM_TOLERANCE)||(*iter < 3))&&(*iter < KMM_MAX_ITERATIONS)) {
+    *iter += 1;
     dL = fabs(logL_bimodal - oldL);
     oldL = logL_bimodal;
@@ -1055 +1061 @@
       *m2 += values->data.F32[i] * P2->data.F32[i];
 
-      w1 += P1[i];
-      w2 += P2[i];
+      w1 += P1->data.F32[i];
+      w2 += P2->data.F32[i];
     }
     *m1 /= w1;
@@ -1095 +1101 @@
 }
 
-static void KMMrejectUnpopular(const psVector *values, psArray *reject) {
+static void KMMRejectUnpopular(const psArray *inputs, int x, int y) {
+  float KMM_MINIMUM_PVALUE = 0.05;
   float Punimodal,pi1,m1,s1,pi2,m2,s2;
-  KMMcalculate(values,&Punimodal,
+  int iter;
+  int j;
+
+  psVector *values = psVectorAlloc(inputs->n, PS_TYPE_F32);
+  for (j = 0; j < inputs->n; j++) {
+    pmStackData *data = inputs->data[j]; // Stack data of interest
+    psImage *image = data->readout->image; // Image of interest
+    int xIn = x - data->readout->col0, yIn = y - data->readout->row0; // Coordinates on input readout
+    values->data.F32[j] = image->data.F32[yIn][xIn];
+  }
+  
+  KMMcalculate(values,&Punimodal,&iter,
                &pi1,&m1,&s1,
                &pi2,&m2,&s2);
+  //  fprintf(stderr,
+  psTrace("psModules.imcombine",3,
+          "KMM Unpopular Test: %d,%d: Puni: %f in %d",x,y,Punimodal,iter);
+  //  CHECKPIX(x, y, "
+  
   if (Punimodal < KMM_MINIMUM_PVALUE) {
     int i;
-    float g1,g2;
+    float g1,g2,norm;
     float P1,P2;
 
@@ -1112 +1135 @@
       P2 = (pi2 * g2) / norm;
 
+      CHECKPIX(x, y, "KMM Unpopular Rejection: %d,%d: %d %f %f:(%f %f %f ) %f:(%f %f %f) rejection? %d %d\n",
+               x, y, i, values->data.F32[i],
+               P1,m1,s1,pi1,
+               P2,m2,s2,pi2,
+               (pi1 > pi2)&&(P1 < P2),
+               (pi1 < pi2)&&(P1 > P2));
       if ((pi1 > pi2)&&(P1 < P2)) { // mode 1 is more popular, but this element belongs to mode 2
-        reject_input(reject,i);
+        combineMarkReject(inputs,x,y,i);
       }
       if ((pi1 < pi2)&&(P1 > P2)) { // mode 2 is more popular, but this element belongs to mode 1
-        reject_input(reject,i);
+        combineMarkReject(inputs,x,y,i);
       }
     }
   }
+  psFree(values);
   // else do nothing.
 }
 
-static void KMMrejectBright(const psVector *values, psArray *reject) {
-  KMMcalculate(values,&Punimodal,
+static void KMMRejectBright(const psArray *inputs, int x, int y) {
+  float KMM_MINIMUM_PVALUE = 0.05;
+  float Punimodal,pi1,m1,s1,pi2,m2,s2;
+  int iter;
+  int j;
+
+  psVector *values = psVectorAlloc(inputs->n, PS_TYPE_F32);
+  for (j = 0; j < inputs->n; j++) {
+    pmStackData *data = inputs->data[j]; // Stack data of interest
+    psImage *image = data->readout->image; // Image of interest
+    int xIn = x - data->readout->col0, yIn = y - data->readout->row0; // Coordinates on input readout
+    values->data.F32[j] = image->data.F32[yIn][xIn];
+  }
+  
+  KMMcalculate(values,&Punimodal,&iter,
                &pi1,&m1,&s1,
                &pi2,&m2,&s2);
   if (Punimodal < KMM_MINIMUM_PVALUE) {
     int i;
-    float g1,g2;
+    float g1,g2,norm;
     float P1,P2;
 
@@ -1140 +1183 @@
 
       if ((m1 > m2)&&(P1 > P2)) { // m1 is larger, and this element belongs to mode 1
-        reject_input(reject,i);
+        combineMarkReject(inputs,x,y,i);
       }
       if ((m1 < m2)&&(P1 < P2)) { // m2 is larger, and this element belongs to mode 2
-        reject_input(reject,i);
+        combineMarkReject(inputs,x,y,i);
       }
     }
@@ -1310 +1353 @@
     }
 
+    // Pre-reject inputs using KMM bimodality test.
+    if (1)  {
+/*       KMMRejectUnpopular(input,x,y); */
+      rejection = true;
+    }
+   
     // Set up rejection list
     psArray *pixelMap = NULL;           // Map of pixels to source
@@ -1319 +1368 @@
     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);
+
+            // Pre-reject inputs using KMM bimodality test.
+          if (1)  {
+            KMMRejectUnpopular(input,x,y);
+/*          rejection = true; */
+          }
+          else {
+            KMMRejectBright(input,x,y);
+          }
 
 #ifdef TESTING