Changeset 17314

Timestamp:

Apr 2, 2008, 8:04:05 PM (18 years ago)

Author:

eugene

Message:

fit each model, validate model selections in recipe

File:

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

branches/eam_branch_20080324/psphot/src/psphotExtendedSourceFits.c

@@ -6 +6 @@
     bool status;
     int Next = 0;
+    int Nconvolve = 0;
+    int Nplain = 0;
 
     // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
@@ -46 +48 @@
         psAbort ("Unknown model class for EXTENDED_SOURCE_MODEL entry %s: %s", item->name, modelType);
       } 
+      psMetadataAddS32 (model, PS_LIST_TAIL, "MODEL_TYPE", PS_META_REPLACE, "", modelType);
 
       // check on the SNLIM, set a float value
@@ -72 +75 @@
     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;
+    }
+
     // source analysis is done in S/N order (brightest first)
     sources = psArraySort (sources, pmSourceSortBySN);
@@ -84 +93 @@
         if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
         if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
-
-        // limit selection to some SN limit
-        // XXX move this inside model loop
-        assert (source->peak); // how can a source not have a peak?
-        if (source->peak->SN < SN_LIM) continue;
 
         // XXX this should use peak?
@@ -119 +123 @@
           assert (status);
 
-          if (X < SNlim) {
-            continue;
-          }
+          // limit selection to some SN limit
+          assert (source->peak); // how can a source not have a peak?
+          if (source->peak->SN < SNlim) continue;
 
           // check on the model type
-          char *modelName = psMetadataLookupStr (&status, model, "MODEL");
+          pmModelType modelType = psMetadataLookupS32 (&status, model, "MODEL_TYPE");
           assert (status);
 
@@ -134 +138 @@
           // XXX use the values above, rework this function
           if (convolved) {
-            pmModel *EXT = psphotPSFConvModel (source, recipe, maskVal);
+            pmModel *EXT = psphotPSFConvModel (readout, source, modelType, maskVal, psfSize);
+            Nconvolve ++;
           } else {
-            pmModel *EXT = psphotFitExt (readout, source, maskVal);
+            pmModel *EXT = psphotFitExt (readout, source, modelType, maskVal);
+            Nplain ++;
           }
           if (!EXT) {
@@ -143 +149 @@
           } 
 
-          // XXX do something with the resulting fit
+          if (source->modelFits == NULL) {
+              source->modelFits = psArrayAllocEmpty (5);
+          }
           
-          
-
+          // test for fit quality / result
+          psArrayAdd (source->modelFits, 5, EXT);
         }
 
         // evaluate the relative quality of the models, choose one
+        float minChisq = 0;0;
+        int minModel = -1;
+        for (int i = 0; i < source->modelFits->n; i++) {
+            pmModel *model = source->modelFits->data[i];
+            
+            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 ((minModel < 0) || (model->chisq < minChisq)) {
+                minChisq = model->chisq;
+                minModel = i;
+            }
+        }           
+
+        if (minModel == -1) {
+            psAbort ("failed to fit extended source model");
+        }       
+
+        // XXX probably need to free the current value
+        psFree (source->modelEXT);
+        source->modelEXT = source->modelFits->data[minModel];
 
         // cache only the chosen model
         pmSourceCacheModel (source, maskVal); // XXX put this in the source model function?
         psTrace ("psphot", 5, "psf-convolved model for source at %7.1f, %7.1f", source->moments->x, source->moments->y);
-        Nconvolve ++;
 
         // re-subtract the object, leave local sky
@@ -163 +194 @@
     psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot"), Next);
     psLogMsg ("psphot", PS_LOG_INFO, "  %d convolved models\n", Nconvolve);
-    psLogMsg ("psphot", PS_LOG_INFO, "  %d petrosian\n", Npetro);
-    psLogMsg ("psphot", PS_LOG_INFO, "  %d isophotal\n", Nisophot);
-    psLogMsg ("psphot", PS_LOG_INFO, "  %d annuli\n", Nannuli);
-    psLogMsg ("psphot", PS_LOG_INFO, "  %d kron\n", Nkron);
+    psLogMsg ("psphot", PS_LOG_INFO, "  %d plain models\n", Nplain);
     return true;
 }