Changeset 35394 for trunk/ppStack/src/ppStackPrepare.c

Timestamp:

Apr 11, 2013, 3:21:06 PM (13 years ago)

Author:

watersc1

Message:

Update to FWHM rejection. I've not fully tested this yet, although it does seem to work as intended in the initial tests.

File:

• trunk/ppStack/src/ppStackPrepare.c (modified) (5 diffs)

trunk/ppStack/src/ppStackPrepare.c

@@ -253 +253 @@
     float maxFWHM = psMetadataLookupF32(&mdok, recipe, "PSF.INPUT.MAX"); // max allowed input fwhm
     float clipFWHMnSig = psMetadataLookupF32(&mdok, recipe, "PSF.INPUT.CLIP.NSIGMA"); // sigma clipping of inputs
-
+    float threshFWHM = psMetadataLookupF32(&mdok, recipe, "PSF.INPUT.THRESH"); // FWHM size that we refuse to clip below
+    float asymmetryFWHM = psMetadataLookupF32(&mdok, recipe, "PSF.INPUT.ASYMMETRY"); // max bimodal asymmetry
     psString log = psStringCopy("Input seeing FWHMs:\n"); // Log message
     bool havePSFs = false;                                // Do we have any PSFs?
@@ -286 +287 @@
             options->inputSeeing->data.F32[i] = sumFWHM / (float)numFWHM;
         }
-
+#ifdef OLD_REJECTION
         // reject any input images which exceed the specified max FWHM
         if (isfinite(maxFWHM) && (options->inputSeeing->data.F32[i] > maxFWHM)) {
@@ -292 +293 @@
             psLogMsg("ppStack", PS_LOG_INFO, "PSF FWHM for image %d is too large (%f vs %f maxFWHM) --- rejected.", i, options->inputSeeing->data.F32[i], maxFWHM);
         }
-
+#endif
         psStringAppend(&log, "Input %d: %f\n", i, options->inputSeeing->data.F32[i]);
     }
@@ -300 +301 @@
     psFree(log);
 
+    // Do GMM test with two modes to decide where to put the break.
+    double Punimodal;
+    int m = 2;
+    psVector *modes = psVectorAlloc(num,PS_TYPE_F32);
+    psVector *means = psVectorAlloc(m,PS_TYPE_F32);
+    psVector *S     = psVectorAlloc(m,PS_TYPE_F32);
+    psVector *pi    = psVectorAlloc(m,PS_TYPE_F32);
+    psImage *P      = psImageAlloc(m,num,PS_TYPE_F32);
+    float limit;
+    if (!psMM1DClass(options->inputSeeing,
+                     options->inputSeeing->n,
+                     modes,means,
+                     S,pi,P,
+                     2,
+                     &Punimodal)) {
+      // Handle error here
+      psFree(modes);
+      psFree(means);
+      psFree(S);
+      psFree(pi);
+      psFree(P);
+      return(false);
+    }
+    fprintf(stderr,"means: %g %g\n",means->data.F32[0],means->data.F32[1]);
+    fprintf(stderr,"sigma: %g %g \n",S->data.F32[0],S->data.F32[1]);
+
+    fprintf(stderr,"pi:    %g %g\n",pi->data.F32[0],pi->data.F32[1]);
+    
+    //    unwrittenGMMfunction(options->inputSeeing, options->inputMask, 0xff,
+    //           &Punimodal,
+    //           &m1,&s1,&pi1,
+    //           &m2,&s2,&pi2);
+
+    // Use Gaussian mixture model analysis of the FWHM distribution to decide an optional FWHM limit
+    if ((Punimodal > 0.5)|| // This distribution is best represented by a single mode
+        (num <= 4)) {       // Or we have a small number of inputs, making this statistic poor.
+      limit = maxFWHM;
+    }
+    else {                  // This is a bimodal distribution
+      if ((fabs(pi->data.F32[0] - pi->data.F32[1]) < asymmetryFWHM)||  // However, both modes are equally populated
+          (pi->data.F32[1] > pi->data.F32[0])) {                       // Or the larger FWHM mode is more populated
+        limit = means->data.F32[1] + clipFWHMnSig * S->data.F32[1];
+      }
+      else {                                   // The smaller FWHM mode is more populated
+        limit = means->data.F32[0] + clipFWHMnSig * S->data.F32[0];
+      }
+    }
+    if (limit > maxFWHM)    { limit = maxFWHM; }    // We should not be larger than our max
+    if (limit < threshFWHM) { limit = threshFWHM; } // Nor smaller than our min
+
+    for (int i = 0; i < num; i++) {
+      if (options->inputSeeing->data.F32[i] > limit) {
+        options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_PSF;
+        psLogMsg("ppStack", PS_LOG_INFO,
+                 "PSF FWHM for image %d is too large (%f vs %f fwhmlimit) --- rejected",
+                 i, options->inputSeeing->data.F32[i],limit);
+      }
+    }
+    psFree(means);
+    psFree(modes);
+    psFree(S);
+    psFree(pi);
+    psFree(P);
+    
+#ifdef OLD_REJECTION
     // We should have the ability to filter the input list based on the seeing:
     // * reject above some max value and/or min value
@@ -323 +389 @@
         }
     }
-
+#endif
     // Generate target PSF
     if (options->convolve) {