Changeset 28687 for branches/eam_branches/ipp-20100621/psphot

Timestamp:

Jul 18, 2010, 2:59:19 PM (16 years ago)

Author:

eugene

Message:

update psf-convolved model fitting to be a bit more consistent with other psphot fitting, and a bit more flexible; fix the diff stats

Location:

branches/eam_branches/ipp-20100621/psphot/src

Files:

• psphotExtendedSourceFits.c (modified) (1 diff)
• psphotSourceFits.c (modified) (3 diffs)

branches/eam_branches/ipp-20100621/psphot/src/psphotExtendedSourceFits.c

@@ -210 +210 @@
           pmModel *modelFit = NULL;
           if (convolved) {
-              modelFit = psphotPSFConvModel (readout, source, modelType, maskVal, markVal, psfSize);
+              modelFit = psphotFitPCM (readout, source, fitOptions, modelType, maskVal, markVal, psfSize);
               if (!modelFit) {
                   psTrace ("psphot", 5, "failed to fit psf-conv model for object at %f, %f", source->moments->Mx, source->moments->My);

branches/eam_branches/ipp-20100621/psphot/src/psphotSourceFits.c

@@ -476 +476 @@
 }
 
+pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal) {
+
+    pmSourceFitOptions options = *fitOptions;
+
+    NfitPCM ++;
+
+    // maskVal is used to test for rejected pixels, and must include markVal
+    maskVal |= markVal;
+
+    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
+
+    // use the source moments, etc to guess basic model parameters
+    pmModel *EXT = pmSourceModelGuess (source, modelType);
+    if (!EXT) {
+        psTrace ("psphot", 5, "failed to generate a model for source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy);
+        return NULL;
+    }
+
+    if ((source->moments->Mxx < 1e-3) || (source->moments->Myy < 1e-3)) {
+        psTrace ("psphot", 5, "problem source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy);
+    }
+
+    // for sersic models, use a grid search to choose an index, then float the params there
+    if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
+        psphotFitSersicIndexPCM (source, EXT, fitOptions, maskVal);
+    }
+
+    if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
+        options.mode = PM_SOURCE_FIT_NO_INDEX;
+    } else {
+        options.mode = PM_SOURCE_FIT_EXT;
+    }
+    pmSourceFitPCM (source, EXT, &options, maskVal);
+
+    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
+    return (EXT);
+}
+
 // note that these should be 1/2n of the standard sersic index
 float indexGuess[] = {0.5, 0.33, 0.25, 0.167, 0.125, 0.083};
+define N_INDEX_GUESS 6
 
 // A sersic model is very sensitive to the index.  attempt to find the index first by grid search in just the index
@@ -491 +530 @@
     options.nIter = 3;
 
-    float xMin, chiSquare[3];
+    float xMin, chiSquare[N_INDEX_GUESS];
     int iMin;
 
-    for (int i = 0; i < 3; i++) {
+    for (int i = 0; i < N_INDEX_GUESS; i++) {
         model->params->data.F32[PM_PAR_7] = indexGuess[i];
         model->modelGuess(model, source);
@@ -516 +555 @@
     return true;
 }
+
+// A sersic model is very sensitive to the index.  attempt to find the index first by grid search in just the index
+// for a sersic model, attempt to fit just the index and normalization with a modest number of iterations
+bool psphotFitSersicIndexPCM (pmSource *source, pmModel *model, pmSourceFitOptions *fitOptions, psImageMaskType maskVal) {
+
+    assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC"));
+
+    pmSourceFitOptions options = *fitOptions;
+    
+    // fit EXT (not PSF) model (set/unset the pixel mask)
+    options.mode = PM_SOURCE_FIT_NO_INDEX;
+    options.nIter = 3;
+
+    float xMin, chiSquare[N_INDEX_GUESS];
+    int iMin;
+
+    for (int i = 0; i < N_INDEX_GUESS; i++) {
+        model->params->data.F32[PM_PAR_7] = indexGuess[i];
+        model->modelGuess(model, source);
+        pmSourceFitPCM (source, model, &options, maskVal);
+        chiSquare[i] = model->chisq;
+        if (i == 0) {
+            xMin = chiSquare[i];
+            iMin = i;
+        } else {
+            if (chiSquare[i] < xMin) {
+                xMin = chiSquare[i];
+                iMin = i;
+            }
+        }
+    }
+
+    model->flags = PM_MODEL_STATUS_NONE; // do not attempt to handle failures here, let the next iteration deal with it
+    model->params->data.F32[PM_PAR_7] = indexGuess[iMin];
+    model->modelGuess(model, source);
+
+    return true;
+}