Changeset 17340

Timestamp:

Apr 6, 2008, 10:14:53 AM (18 years ago)

Author:

eugene

Message:

test output

File:

• branches/eam_branch_20080324/psphot/src/psphotExtendedSourceFits.c (modified) (11 diffs)

branches/eam_branch_20080324/psphot/src/psphotExtendedSourceFits.c

@@ -1 @@
-# include "psphot.h"
+# include "psphotInternal.h"
 
 // non-linear model fitting for extended sources
@@ -8 +8 @@
     int Nconvolve = 0;
     int Nplain = 0;
+    bool savePics = true;
 
     // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
@@ -40 +41 @@
       }
 
-      psMetadata *model = (psMetadata *) item->data;
+      psMetadata *model = (psMetadata *) item->data.md;
 
       // check on the model type
       char *modelName = psMetadataLookupStr (&status, model, "MODEL");
-      int modelType = psModelClassGetType (modelName);
+      int modelType = pmModelClassGetType (modelName);
       if (modelType < 0) {
-        psAbort ("Unknown model class for EXTENDED_SOURCE_MODEL entry %s: %s", item->name, modelType);
+        psAbort ("Unknown model class for EXTENDED_SOURCE_MODEL entry %s: %s", item->name, modelName);
       } 
       psMetadataAddS32 (model, PS_LIST_TAIL, "MODEL_TYPE", PS_META_REPLACE, "", modelType);
@@ -61 +62 @@
       char *convolvedWord = psMetadataLookupStr (&status, model, "PSF_CONVOLVED");
       if (!status || (strcasecmp (convolvedWord, "true") && strcasecmp (convolvedWord, "false"))) {
-        psAbort ("PSF_CONVOLVED entry invalid or missing for EXTENDED_SOURCE_MODEL entry %s", item->name, modelType);
+        psAbort ("PSF_CONVOLVED entry invalid or missing for EXTENDED_SOURCE_MODEL entry %s", item->name);
       } 
       bool convolved = !strcasecmp (convolvedWord, "true");
       psMetadataAddBool (model, PS_LIST_TAIL, "PSF_CONVOLVED_VALUE", PS_META_REPLACE, "", convolved);
     }
+    psFree (iter);
 
     // option to limit analysis to a specific region
@@ -72 +74 @@
     if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
 
-    // S/N limit to perform full non-linear fits
-    float SN_LIM = psMetadataLookupF32 (&status, recipe, "EXTENDED_SOURCE_SN_LIM");
-
     // what fraction of the PSF is used? (radius in pixels : 2 -> 5x5 box)
     int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE");
-    if (!status) {
-        psfSize = 2;
-    }
+    assert (status);
 
     // source analysis is done in S/N order (brightest first)
@@ -110 +107 @@
         Next ++;
 
+        if (savePics) {
+          psphotSaveImage (NULL, readout->image, "image.xp.fits");
+        }
+
         // loop here over the models chosen for each source (exclude by S/N)
         psMetadataIterator *iter = psMetadataIteratorAlloc (models, PS_LIST_HEAD, NULL);
@@ -117 +118 @@
           // XXX this should have been forced above
           assert (item->type == PS_DATA_METADATA);
-          psMetadata *model = (psMetadata *) item->data;
+          psMetadata *model = (psMetadata *) item->data.md;
 
           // check the SNlim and skip model if source is too faint
@@ -135 +136 @@
           assert (status);
 
+          psTraceSetLevel ("psLib.math.psMinimizeLMChi2", 6);
+
           // fit the model as convolved or not
-          // XXX use the values above, rework this function
+          pmModel *modelFit = NULL;
           if (convolved) {
-            pmModel *EXT = psphotPSFConvModel (readout, source, modelType, maskVal, psfSize);
+            modelFit = psphotPSFConvModel (readout, source, modelType, maskVal, psfSize);
+            psTrace ("psphot", 4, "fit psf-conv model for %f, %f : %s chisq = %f\n", source->moments->x, source->moments->y, pmModelClassGetName (modelFit->type), modelFit->chisq);
             Nconvolve ++;
           } else {
-            pmModel *EXT = psphotFitExt (readout, source, modelType, maskVal);
+            modelFit = psphotFitEXT (readout, source, modelType, maskVal);
+            psTrace ("psphot", 4, "fit plain model for %f, %f : %s chisq = %f\n", source->moments->x, source->moments->y, pmModelClassGetName (modelFit->type), modelFit->chisq);
             Nplain ++;
           }
-          if (!EXT) {
-            psError(PSPHOT_ERR_UNKNOWN, false, "failure to fit extended model");
-            return false;
+          if (!modelFit) {
+            psTrace ("psphot", 5, "failed to generate model guess for object at %f, %f", source->moments->x, source->moments->y);
+            continue;
           } 
+
+          if (savePics) {
+            pmModel *saveModel = source->modelEXT;
+            source->modelEXT = modelFit;
+            pmSourceCacheModel (source, maskVal); // XXX put this in the source model function?
+            pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+
+            psphotSaveImage (NULL, readout->image, "image.xf.fits");
+
+            pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+            source->modelEXT = saveModel;
+
+            char key[10];
+            fprintf (stdout, "continue? ");
+            fgets (key, 8, stdin);
+            if (key[0] == 'n') {
+              savePics = false;
+            }
+          }
 
           if (source->modelFits == NULL) {
@@ -154 +178 @@
           
           // test for fit quality / result
-          psArrayAdd (source->modelFits, 5, EXT);
+          psArrayAdd (source->modelFits, 5, modelFit);
+          psFree (modelFit);
         }
+        psFree (iter);
 
         // evaluate the relative quality of the models, choose one
-        float minChisq = 0;0;
+        float minChisq = NAN;
         int minModel = -1;
         for (int i = 0; i < source->modelFits->n; i++) {
@@ -176 +202 @@
 
         if (minModel == -1) {
-            psAbort ("failed to fit extended source model");
+          // re-subtract the object, leave local sky
+          psTrace ("psphot", 5, "failed to fit extended source model to object at %f, %f", source->moments->x, source->moments->y);
+          pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+          source->mode |= PM_SOURCE_MODE_SUBTRACTED;
+
+          if (savePics) {
+            psphotSaveImage (NULL, readout->image, "image.xp.fits");
+            char key[10];
+            fprintf (stdout, "continue? ");
+            fgets (key, 8, stdin);
+            if (key[0] == 'n') {
+              savePics = false;
+            }
+          }
+
+          continue;
         }       
 
         // XXX probably need to free the current value
         psFree (source->modelEXT);
-        source->modelEXT = source->modelFits->data[minModel];
+        source->modelEXT = psMemIncrRefCounter (source->modelFits->data[minModel]);
+
+        psTrace ("psphot", 4, "best ext model for %f, %f : %s chisq = %f\n", source->moments->x, source->moments->y, pmModelClassGetName (source->modelEXT->type), source->modelEXT->chisq);
 
         // cache only the chosen model
@@ -190 +233 @@
         pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
         source->mode |= PM_SOURCE_MODE_SUBTRACTED;
+
+        if (savePics) {
+          psphotSaveImage (NULL, readout->image, "image.xm.fits");
+          char key[10];
+          fprintf (stdout, "continue? ");
+          fgets (key, 8, stdin);
+          if (key[0] == 'n') {
+            savePics = false;
+          }
+        }
     }