Changeset 37066 for branches/eam_branches/ipp-ops-20130712/psphot/src/psphotExtendedSourceFits.c

Timestamp:

Jul 17, 2014, 12:30:45 PM (12 years ago)

Author:

eugene

Message:

merge changes (from past YEAR) into trunk

Location:

branches/eam_branches/ipp-ops-20130712/psphot

Files:

• . (modified) (1 prop)
• src (modified) (2 props)
• src/psphotExtendedSourceFits.c (modified) (20 diffs)

branches/eam_branches/ipp-ops-20130712/psphot/src

@@ -24 +24 @@
 psphotModelTest
 psphotMinimal
+psphotFullForce
+psmakecff
+psphotFullForceSummary

@@ -24 +24 @@
 psphotModelTest
 psphotMinimal
+psphotFullForce
+psmakecff
+psphotFullForceSummary

branches/eam_branches/ipp-ops-20130712/psphot/src/psphotExtendedSourceFits.c

@@ -111 +111 @@
     }
 
+    float fitNsigmaConv = psMetadataLookupF32 (&status, recipe, "EXT_FIT_NSIGMA_CONV"); // number of sigma for the convolution
+    if (!status || !isfinite(fitNsigmaConv) || fitNsigmaConv <= 0) {
+        fitNsigmaConv = 5.0;
+    }
+
     bool chisqConvergence = psMetadataLookupBool (&status, recipe, "LMM_FIT_CHISQ_CONVERGENCE"); // Fit tolerance
     if (!status) {
@@ -123 +128 @@
     }
 
+    // perform full extended source non-linear fits?
+    bool isInteractive = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_FITS_INTERACTIVE");
+    if (!status) isInteractive = false;
+
     // Define source fitting parameters for extended source fits
     pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc();
-    fitOptions->mode           = PM_SOURCE_FIT_EXT;
+    fitOptions->mode           = PM_SOURCE_FIT_EXT_AND_SKY;
     fitOptions->saveCovariance = true;  // XXX make this a user option?
     fitOptions->covarFactor    = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
@@ -131 +140 @@
     fitOptions->minTol         = fitMinTol;
     fitOptions->maxTol         = fitMaxTol;
+    fitOptions->nsigma         = fitNsigmaConv;
 
     fitOptions->gainFactorMode   = gainFactorMode;
     fitOptions->chisqConvergence = chisqConvergence;
+    fitOptions->isInteractive    = isInteractive;
+
+    // use poissonian errors or local-sky errors
+    fitOptions->poissonErrors = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_FITS_POISSON");
+    if (!status) fitOptions->poissonErrors = true;
 
     // maskVal is used to test for rejected pixels, and must include markVal
@@ -278 +293 @@
             psArrayAdd(job->args, 1, cells->data[j]); // sources
             psArrayAdd(job->args, 1, models);
+
             // Allocate a metadata iterator here because psMetadataIteratorAlloc/Free are not thread safe
             psMetadataIterator *iter = psMetadataIteratorAlloc (models, PS_LIST_HEAD, NULL);
@@ -296 +312 @@
             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail
 
-// set this to 0 to run without threading
-# if (1)            
-            if (!psThreadJobAddPending(job)) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
-                psFree(AnalysisRegion);
-                psFree (fitOptions);
-                psFree (models);
-                psphotSersicModelClassCleanup();
-                return false;
-            } 
-# else
-            if (!psphotExtendedSourceFits_Threaded(job)) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
-                psFree(AnalysisRegion);
-                psFree (fitOptions);
-                psFree (models);
-                psphotSersicModelClassCleanup();
-                return false;
+// XXX TEST 
+            if (!isInteractive) {
+                if (!psThreadJobAddPending(job)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+                    psFree(AnalysisRegion);
+                    psFree (fitOptions);
+                    psFree (models);
+                    psphotSersicModelClassCleanup();
+                    return false;
+                } 
+            } else {
+                // run without threading
+                if (!psphotExtendedSourceFits_Threaded(job)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+                    psFree(AnalysisRegion);
+                    psFree (fitOptions);
+                    psFree (models);
+                    psphotSersicModelClassCleanup();
+                    return false;
+                }
+                psScalar *scalar = NULL;
+                scalar = job->args->data[9];
+                Next += scalar->data.S32;
+                scalar = job->args->data[10];
+                Nconvolve += scalar->data.S32;
+                scalar = job->args->data[11];
+                NconvolvePass += scalar->data.S32;
+                scalar = job->args->data[12];
+                Nplain += scalar->data.S32;
+                scalar = job->args->data[13];
+                NplainPass += scalar->data.S32;
+                scalar = job->args->data[14];
+                Nfaint += scalar->data.S32;
+                scalar = job->args->data[15];
+                Nfail += scalar->data.S32;
+                psFree(job->args->data[3]); // iterator allocated above
+                psFree(job);
             }
-            psScalar *scalar = NULL;
-            scalar = job->args->data[9];
-            Next += scalar->data.S32;
-            scalar = job->args->data[10];
-            Nconvolve += scalar->data.S32;
-            scalar = job->args->data[11];
-            NconvolvePass += scalar->data.S32;
-            scalar = job->args->data[12];
-            Nplain += scalar->data.S32;
-            scalar = job->args->data[13];
-            NplainPass += scalar->data.S32;
-            scalar = job->args->data[14];
-            Nfaint += scalar->data.S32;
-            scalar = job->args->data[15];
-            Nfail += scalar->data.S32;
-            psFree(job->args->data[3]); // iterator allocated above
-            psFree(job);
-# endif
         }
 
@@ -418 +435 @@
     // psTraceSetLevel ("psLib.math.psMinimizeLMChi2_Alt", 5);
 
+    pmModelStatus badModel = PM_MODEL_STATUS_NONE;
+    badModel |= PM_MODEL_STATUS_BADARGS;
+    badModel |= PM_MODEL_STATUS_OFFIMAGE;
+    badModel |= PM_MODEL_STATUS_NAN_CHISQ;
+    badModel |= PM_MODEL_SERSIC_PCM_FAIL_GUESS;
+    badModel |= PM_MODEL_SERSIC_PCM_FAIL_GRID;
+    badModel |= PM_MODEL_PCM_FAIL_GUESS;
+
     // choose the sources of interest
     for (int i = 0; i < sources->n; i++) {
@@ -423 +448 @@
         pmSource *source = sources->data[i];
 
-        // skip PSF-like and non-astronomical objects
-        if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue;
-        if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
-        if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
-
-        // skip saturated stars modeled with a radial profile 
-        if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue;
-
-        // XXX this should use peak?
+        // rules for measuring petrosian parameters for specific objects are set in
+        // psphotChooseAnalysisOptions.c
+        if (!(source->tmpFlags & PM_SOURCE_TMPF_EXT_FIT)) continue;
+
+        // limit selection by analysis region (XXX move this into psphotChooseAnalysisOption?)
         if (source->peak->x < region->x0) continue;
         if (source->peak->y < region->y0) continue;
         if (source->peak->x > region->x1) continue;
         if (source->peak->y > region->y1) continue;
-
-
-        // XXX for a test, just do the obvious trail
-        // XXX if (source->peak->xf < 1100) continue;
-        // XXX if (source->peak->xf > 1400) continue;
-        // XXX if (source->peak->yf >  245) continue;
 
         // replace object in image
@@ -452 +467 @@
         psphotSetRadiusMomentsExact(&fitRadius, &windowRadius, readout, source, markVal); // NOTE : 6 allocs
 
-        // XXX WATCH OUT HERE!!
-        // fitRadius = 30;
-
         // UPDATE : we have changed the moments calculation.  There is now an iteration within 
         // psphotKronMasked to determine moments appropriate for a larger object.  The values
@@ -464 +476 @@
             source->modelFits = psArrayAllocEmpty (models->list->n);
         }
+
+# ifdef TEST_OBJECT
+        bool testObject = false;
+        testObject |= ((fabs(source->peak->xf -  179) < 5) && (fabs(source->peak->yf - 1138) < 5));
+        if (testObject) {
+            fprintf (stderr, "test object @ %f, %f\n", source->peak->xf, source->peak->yf);
+            psTraceSetLevel ("psModules.objects.pmPCM_MinimizeChisq", 5);
+            psTraceSetLevel ("psphot.psphotExtendedSourceFits_Threaded", 5);
+        }
+# endif 
 
         // loop here over the models chosen for each source (exclude by S/N)
@@ -479 +501 @@
           assert (status);
 
-          // limit selection to some SN limit
-          bool skipSource = false;
-          if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
-              skipSource = (source->moments->KronFlux < FIT_SN_LIM * source->moments->KronFluxErr);
-          } else {
-              skipSource = (sqrt(source->peak->detValue) < FIT_SN_LIM);
-          }
-          if (skipSource) {
+          // limit selection to some SN limit for specific models (this value only applies if > EXTENDED_SOURCE_SN_LIM)
+          if (isfinite(FIT_SN_LIM)) {
+            if (source->moments->KronFlux < FIT_SN_LIM * source->moments->KronFluxErr) {
               Nfaint ++;
               continue;
+            }
           }
+
+          source->mode2 |= PM_SOURCE_MODE2_EXT_FITS_RUN;
 
           // check on the model type
@@ -507 +527 @@
                   psTrace ("psphot", 5, "failed to fit psf-conv model for object at %f, %f", source->moments->Mx, source->moments->My);
                   Nfail ++;
+                  source->mode2 |= PM_SOURCE_MODE2_EXT_FITS_FAIL;
                   continue;
               }
@@ -512 +533 @@
                        source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq, modelFit->nPix, modelFit->nIter);
               Nconvolve ++;
-              if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) {
+              if (!(modelFit->flags & badModel)) {
                   NconvolvePass ++;
                   source->mode |= PM_SOURCE_MODE_EXTENDED_FIT;
@@ -526 +547 @@
                       Nfail ++;
                       doneFits = true;
+                      source->mode2 |= PM_SOURCE_MODE2_EXT_FITS_FAIL;
                       continue;
                   }
                   psTrace ("psphot", 4, "fit plain model for %f, %f : %s chisq = %f (npix: %d, niter: %d)\n", source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq, modelFit->nPix, modelFit->nIter);
                   Nplain ++;
-                  if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) {
+                  if (!(modelFit->flags & badModel)) {
                       NplainPass ++;
                       source->mode |= PM_SOURCE_MODE_EXTENDED_FIT;
@@ -544 +566 @@
                       float Ymax = PAR[PM_PAR_YPOS] + 0.5*PAR[PM_PAR_LENGTH]*cos(PAR[PM_PAR_THETA]);
 
-                      if (false && (source->peak->xf > 1100) && 
-                          (source->peak->xf < 1400) && 
-                          (source->peak->yf < 245)) {
+                      if ((fabs(source->peak->xf - 2572) < 20) && (fabs(source->peak->yf - 5874) < 20)) {
                           fprintf (stderr, "src vs fit : %d %d - %d %d | %f %f - %f %f\n", 
                                    source->pixels->col0, source->pixels->row0, 
@@ -556 +576 @@
                           fprintf (stderr, "update window : %f %f : %f -> %f\n", source->peak->xf, source->peak->yf, fitRadius, 2*fitRadius);
                           psphotSetWindowTrail (&fitRadius, &windowRadius, readout, source, markVal, fitRadius*2.0);
+                          source->mode2 |= PM_SOURCE_MODE2_EXT_FITS_RETRY;
                       }
                   }
               }
           }
-          // XXX really need to do this in a cleaner way:
-          if (!modelFit) continue;
+          psAssert (modelFit, "modelFit not set?");
 
           // test for fit quality / result
@@ -570 +590 @@
         }
 
+        // we are allowed to fit both stars and non-stars here -- if we have fitted
+        // something which we think is a star, we should use that model to subtract the
+        // object from the image.
+        if (source->type == PM_SOURCE_TYPE_STAR) {
+          // ensure the modelPSF is cached
+          pmSourceCacheModel (source, maskVal);
+          pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+          continue;
+        }
+
         // evaluate the relative quality of the models, choose one
+        // the PSF model might be the best fit : allow it to succeed
         float minChisq = NAN;
         int minModel = -1;
@@ -576 +607 @@
             pmModel *model = source->modelFits->data[i];
 
+            // skip the really bad fits
             if (!(model->flags & PM_MODEL_STATUS_FITTED)) continue;
-
-            if (model->flags & (PM_MODEL_STATUS_BADARGS)) continue;
-            if (model->flags & (PM_MODEL_STATUS_NONCONVERGE)) continue;
-            if (model->flags & (PM_MODEL_STATUS_OFFIMAGE)) continue;
+            if (model->flags & badModel) continue;
+
+            // if (model->flags & (PM_MODEL_STATUS_NONCONVERGE)) continue;
 
             if ((minModel < 0) || (model->chisq < minChisq)) {
@@ -598 +629 @@
           pmSourceCacheModel (source, maskVal);
 
+# if (PS_TRACE_ON)
+          pmModel *model = source->modelFits->data[0];
+          int flags = 0xffffffff;
+          if (model) {
+            flags = model->flags;
+          }
+          fprintf (stderr, "failed to fit extended source model to object %d @ %f, %f (%x)\n", source->id, source->moments->Mx, source->moments->My, flags);
+#endif
           pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
 
+          source->mode2 |= PM_SOURCE_MODE2_EXT_FITS_NONE;
           continue;
         }
@@ -615 +655 @@
         // cache the model flux
         if (source->modelEXT->isPCM) {
-            pmPCMCacheModel (source, maskVal, psfSize);
+            // fprintf (stderr, "subtract PCM extended source model for object %d @ %f, %f\n", source->id, source->moments->Mx, source->moments->My);
+            pmPCMCacheModel (source, maskVal, psfSize, fitOptions->nsigma);
         } else {
+            // fprintf (stderr, "subtract non-PCM extended source model for object %d @ %f, %f\n", source->id, source->moments->Mx, source->moments->My);
             pmSourceCacheModel (source, maskVal);
         }
+        source->modelEXT->flags |= PM_MODEL_BEST_FIT;
 
         // subtract the best fit from the object, leave local sky
@@ -625 +668 @@
         psTrace ("psphot", 4, "best ext model for %f, %f : %s chisq = %f\n", source->moments->Mx, source->moments->My, pmModelClassGetName (source->modelEXT->type), source->modelEXT->chisq);
         psTrace ("psphot", 5, "extended source model for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
+
+# ifdef TEST_OBJECT
+        if (testObject) {
+            psTraceSetLevel ("psModules.objects.pmPCM_MinimizeChisq", 0);
+            psTraceSetLevel ("psphot.psphotExtendedSourceFits_Threaded", 0);
+        }
+# endif
     }
 

@@ -24 +24 @@
 psphotModelTest
 psphotMinimal
+psphotFullForce
+psmakecff
+psphotFullForceSummary