Changeset 35394

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.

Location:

trunk

Files:

• ippconfig/gpc1/ppStack.config (modified) (2 diffs)
• ippconfig/recipes/ppStack.config (modified) (1 diff)
• ppStack/src/ppStackPrepare.c (modified) (5 diffs)
• psLib/src/math/Makefile.am (modified) (2 diffs)
• psLib/src/math/psMixtureModels.c (added)
• psLib/src/math/psMixtureModels.h (added)
• psLib/src/pslib_strict.h (modified) (1 diff)

trunk/ippconfig/gpc1/ppStack.config

@@ -1 +1 @@
 PSF.MODEL       STR     PS_MODEL_PS1_V1 # Model for PSF generation
+BACKGROUND.MODEL BOOL   F
 
 ZP.AIRMASS      METADATA                # Airmass terms by filter
@@ -30 +31 @@
     OUTPUT.LOGFLUX          BOOL  FALSE
     OUTPUT.REPLICATE        BOOL  TRUE
-    PSF.INPUT.MAX           F32   7.5
+    PSF.INPUT.MAX           F32   10.0
     PSF.INPUT.CLIP.NSIGMA   F32   1.5
+    PSF.INPUT.THRESH        F32   5.0
+    PSF.INPUT.ASYMMETRY     F32   0.2
 END
 

trunk/ippconfig/recipes/ppStack.config

@@ -92 +92 @@
 PSF.INPUT.MAX           F32  NAN
 PSF.INPUT.CLIP.NSIGMA   F32  NAN
+PSF.INPUT.THRESH        F32  NAN
+PSF.INPUT.ASYMMETRY     F32  NAN
 
 TEMP.IMAGE      STR     conv.im.fits    # Suffix for temporary convolved images

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) {

trunk/psLib/src/math/Makefile.am

@@ -16 +16 @@
         psMinimizePowell.c \
         psMinimizePolyFit.c \
+        psMixtureModels.c \
         psPolynomial.c \
         psPolynomialMetadata.c \
@@ -45 +46 @@
         psMinimizePowell.h \
         psMinimizePolyFit.h \
+        psMixtureModels.h \
         psPolynomial.h \
         psPolynomialMetadata.h \

trunk/psLib/src/pslib_strict.h

@@ -76 +76 @@
 #include "psMinimizePowell.h"
 #include "psMinimizePolyFit.h"
+#include "psMixtureModels.h"
 #include "psMutex.h"
 #include "psRandom.h"