Changeset 36086 for trunk/psphot/src/psphotExtendedSourceFits.c

Timestamp:

Aug 31, 2013, 6:02:26 AM (13 years ago)

Author:

eugene

Message:

merge changes from EAM branch (update tap script for galaxies; count deblended satstars for log message; scale down EXT window a bit for speed (OK?); enable interactive mode for PCM fitting; allow NONCONVERGE extended fit models to succeed; init EXT model fits with EXT_AND_SKY just in case; extra verbosity and various experiments on the index iterations)

Location:

trunk/psphot

Files:

• . (modified) (1 prop)
• src (modified) (1 prop)
• src/psphotExtendedSourceFits.c (modified) (10 diffs)

trunk/psphot/src/psphotExtendedSourceFits.c

@@ -65 +65 @@
     psphotVisualShowImage(readout);
 
+    // psphotSaveImage (NULL, readout->image, "test.01.fits");
+
     pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     psAssert (detections, "missing detections?");
@@ -123 +125 @@
     }
 
+    // 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
@@ -134 +140 @@
     fitOptions->gainFactorMode   = gainFactorMode;
     fitOptions->chisqConvergence = chisqConvergence;
+    fitOptions->isInteractive    = isInteractive;
 
     // maskVal is used to test for rejected pixels, and must include markVal
@@ -296 +303 @@
             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
         }
 
@@ -378 +386 @@
     psphotSersicModelClassCleanup();
 
+    // psphotSaveImage (NULL, readout->image, "test.02.fits");
+
     psphotVisualShowResidualImage (readout, false);
 
@@ -521 +531 @@
                   psFree (source->modelFlux);
                   source->modelFlux = NULL;
+                  // if ((fabs(source->peak->xf - 2572) < 20) && (fabs(source->peak->yf - 5874) < 20)) {
+                  //   fprintf (stderr, "test object\n");
+                  // }
+
                   modelFit = psphotFitEXT (modelFit, readout, source, fitOptions, modelType, maskVal, markVal);
                   if (!modelFit) {
@@ -544 +558 @@
                       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, 
@@ -579 +591 @@
 
             if (model->flags & (PM_MODEL_STATUS_BADARGS)) continue;
-            if (model->flags & (PM_MODEL_STATUS_NONCONVERGE)) continue;
+            // if (model->flags & (PM_MODEL_STATUS_NONCONVERGE)) continue;
             if (model->flags & (PM_MODEL_STATUS_OFFIMAGE)) continue;
 
@@ -598 +610 @@
           pmSourceCacheModel (source, maskVal);
 
+          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);
           pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
 
@@ -615 +634 @@
         // cache the model flux
         if (source->modelEXT->isPCM) {
+            // 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);
         } 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);
         }