Changeset 28657 for branches/eam_branches/ipp-20100621/psphot/src/psphotSourceFits.c

Timestamp:

Jul 11, 2010, 3:21:36 PM (16 years ago)

Author:

eugene

Message:

add maxTol and maxChisqDOF to psMin to allow for better control over fit success and convergence; for SERSIC model, first do a grid search for the index

File:

• branches/eam_branches/ipp-20100621/psphot/src/psphotSourceFits.c (modified) (3 diffs)

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

@@ -1 +1 @@
 # include "psphotInternal.h"
+bool psphotFitSersicIndex (pmSource *source, pmModel *model, pmSourceFitOptions *fitOptions, psImageMaskType maskVal);
 
 // given a source with an existing modelPSF, attempt a full PSF fit, subtract if successful
@@ -460 +461 @@
 
     // for sersic models, get the initial guess more carefully
-    // if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
-    //  psphotGuessSersic ();
-    //  // test and return NULL on failure?
-    // }
+    if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
+        psphotFitSersicIndex (source, EXT, fitOptions, maskVal);
+    }
 
     // fit EXT (not PSF) model (set/unset the pixel mask)
@@ -473 +473 @@
 }
 
+float indexGuess[] = {0.5, 0.25, 0.125};
+
 // A sersic model is very sensitive to the index.  attempt to find the index first by grid search in just the index
-// bool psphotFitSersic (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal) {
-// 
-//     assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC"));
-// 
-// 
-//     
-//     if (!model->modelGuess(model, source)) {
-//     }
-//     
-// 
-// }
+// for a sersic model, attempt to fit just the index and normalization with a modest number of iterations
+bool psphotFitSersicIndex (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[3];
+    int iMin;
+
+    for (int i = 0; i < 3; i++) {
+        model->params->data.F32[PM_PAR_7] = indexGuess[i];
+        model->modelGuess(model, source);
+        pmSourceFitModel (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;
+}