Changeset 35529

Timestamp:

May 7, 2013, 2:30:23 PM (13 years ago)

Author:

watersc1

Message:

Added more rejection options to help with refstacks.

Location:

trunk

Files:

• ippconfig/recipes/ppStack.config (modified) (1 diff)
• ppStack/src/ppStackPrepare.c (modified) (3 diffs)

trunk/ippconfig/recipes/ppStack.config

@@ -90 +90 @@
 
 
-PSF.INPUT.MAX           F32  NAN
-PSF.INPUT.CLIP.NSIGMA   F32  NAN
-PSF.INPUT.THRESH        F32  NAN
-PSF.INPUT.ASYMMETRY     F32  NAN
+PSF.INPUT.CLIP.SIMPLE   BOOL F          # Should we do a simple sigma clip (T, requires CLIP.NSIGMA) or the mixture model (F, requires all parameters)
+PSF.INPUT.MAX           F32  NAN        # Set maximum limit on input FWHM (defaults to inf)
+PSF.INPUT.CLIP.NSIGMA   F32  NAN        # Set n-sigma cut limit for input rejection
+PSF.INPUT.THRESH        F32  NAN        # Set minimum limit below which we do not exclude an input (defaults to 0.0)
+PSF.INPUT.ASYMMETRY     F32  NAN        # Set difference in mixture model populations to consider equal.
 
 TEMP.IMAGE      STR     conv.im.fits    # Suffix for temporary convolved images

trunk/ppStack/src/ppStackPrepare.c

@@ -257 +257 @@
 
     bool mdok = false;
+    bool  simpleClip = psMetadataLookupF32(&mdok, recipe, "PSF.INPUT.CLIP.SIMPLE");
     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?
@@ -309 +311 @@
 
 #ifdef NEW_REJECTION
-    // 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
+    // Begin new rejection code.
+    float limit = INFINITY; // Default to a big value.
+    if (!isfinite(threshFWHM)) { // If this isn't set, initialize it to zero.  No inputs may be clipped if smaller than this
+      threshFWHM = 0.0;
+    }
+    if (!isfinite(maxFWHM)) { // No maxFWHM was set, accept all inputs
+      maxFWHM = INFINITY;
+    }
+    limit = maxFWHM; // Initialize the limit to be the maxFWHM
+    if (limit < threshFWHM) {
+      limit = threshFWHM;
+    }
+    if (!isfinite(clipFWHMnSig)) { // We cannot do a sigma clip with nSig, so use the maxFWHM to set the limit.
+      psLogMsg("ppStack",PS_LOG_INFO,
+               "PSF FWHM distribution: Unable to determine limit based on recipe. Setting to %f. (clipFWHMnSig: %f maxFWHM: %f threshFWHM: %f asymmetryFWHM: %f)",
+               limit,
+               clipFWHMnSig,maxFWHM,threshFWHM,asymmetryFWHM);
+    }
+    if (simpleClip) { // Do a sigma clip like the old rejection code.
+      psStats *fwhmStats = psStatsAlloc(PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
+      psVectorStats (fwhmStats, options->inputSeeing, NULL, options->inputMask, 0xff);
+      psLogMsg("ppStack", PS_LOG_INFO, "Input FWHMs : %f +/- %f", fwhmStats->clippedMean, fwhmStats->clippedStdev);
+      
+      options->clippedMean =  fwhmStats->clippedMean;
+      options->clippedStdev = fwhmStats->clippedStdev;
+      limit = options->clippedMean + clipFWHMnSig * options->clippedStdev;
+      if (limit < threshFWHM) {
+        limit = threshFWHM;
+      }
+      if (limit > maxFWHM) {
+        limit = maxFWHM;
+      }
+      psLogMsg("ppStack",PS_LOG_INFO,
+               "PSF FWHM distribution: Used simple clip method.  Limit set to %f.",
+               limit);
+      psFree(fwhmStats);
+    }
+    else { // Do the mixture model rejection.
+      if(!isfinite(asymmetryFWHM)) { // Or not, because the parameters aren't set.
+        psLogMsg("ppStack",PS_LOG_INFO,
+                 "PSF FWHM distribution: Unable to determine limit based on recipe. Setting to %f. (clipFWHMnSig: %f maxFWHM: %f threshFWHM: %f asymmetryFWHM: %f)",
+                 limit,
+                 clipFWHMnSig,maxFWHM,threshFWHM,asymmetryFWHM);
+      }
+      // 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);
+      
+      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]);
+      
+      // 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
+      psLogMsg("ppStack",PS_LOG_INFO,
+               "PSF FWHM distribution: limit: %f (%f %f %f) (%f %f %f) %f",
+               limit,means->data.F32[0],S->data.F32[0],pi->data.F32[0],
+               means->data.F32[1],S->data.F32[1],pi->data.F32[1],
+               Punimodal);
+      psFree(means);
       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]);
-
-    // 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
-    psLogMsg("ppStack",PS_LOG_INFO,
-             "PSF FWHM distribution: limit: %f (%f %f %f) (%f %f %f) %f",
-             limit,means->data.F32[0],S->data.F32[0],pi->data.F32[0],
-             means->data.F32[1],S->data.F32[1],pi->data.F32[1],
-             Punimodal);
+    } // End mixture model case
+
+    // Perform rejection using the limit set
     for (int i = 0; i < num; i++) {
       if (options->inputSeeing->data.F32[i] > limit) {
@@ -365 +419 @@
       }
     }
-    psFree(means);
-    psFree(modes);
-    psFree(S);
-    psFree(pi);
-    psFree(P);
     // End new rejection code
 #endif