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

Timestamp:

Jul 22, 2010, 6:21:31 PM (16 years ago)

Author:

eugene

Message:

pcm (psf-convolved model) fitting now works (but needs to go faster)

File:

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

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

@@ -479 +479 @@
 pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
 
+    // allocate the model
+    pmModel *model = pmModelAlloc(modelType);
+    if (!model) {
+        return NULL;
+    }
+
+    pmSourceFitOptions options = *fitOptions;
+
     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);
@@ -489 +497 @@
     // at this stage, skip Gaussian windowing, and do not clip pixels by S/N
     // this uses the footprint to judge both radius and aperture?
-    if (!pmSourceMoments (source, radius, 0.0, 0.0, maskVal)) return false;
-
-    pmSourceFitOptions options = *fitOptions;
+    if (!pmSourceMoments (source, radius, 0.0, 0.0, maskVal)) {
+        // XXX set some mask bit/
+        model->flags |= PM_MODEL_STATUS_BADARGS;
+        return model;
+    }
 
     NfitPCM ++;
@@ -500 +510 @@
     // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
 
-    pmPCMdata *pcm = pmPCMinit (source, &options, modelType, maskVal, psfSize);
+    pmPCMdata *pcm = pmPCMinit (source, &options, model, maskVal, psfSize);
     if (!pcm) {
         psTrace ("psphot", 5, "failed to generate a model for source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy);
-        return NULL;
-    }
-    // XXX check for nDOF too small
+        model->flags |= PM_MODEL_STATUS_BADARGS; // XXX this is probably already set in pmPCMinit
+        return model;
+    }
 
     // use the source moments, etc to guess basic model parameters
@@ -512 +522 @@
     // for sersic models, use a grid search to choose an index, then float the params there
     if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
-        psphotFitSersicIndexPCM (pcm, source, fitOptions, maskVal, markVal);
+        psphotFitSersicIndexPCM (pcm, source, fitOptions, maskVal, markVal, psfSize);
     }
 
@@ -520 +530 @@
         options.mode = PM_SOURCE_FIT_EXT;
     }
-    pmSourceFitPCM (source, PCM, &options, maskVal);
+    pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
 
     // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
-    return (PCM);
+    psFree (pcm);
+
+    return model;
 }
 
 // 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
+# 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
 // for a sersic model, attempt to fit just the index and normalization with a modest number of iterations
-bool psphotFitSersicIndex (pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) {
-
-    pmModel *model = pcm->modelConv;
+bool psphotFitSersicIndex (pmSource *source, pmModel *model, pmSourceFitOptions *fitOptions, psImageMaskType maskVal) {
 
     assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC"));
@@ -544 +554 @@
     options.nIter = 3;
 
-    float xMin, chiSquare[N_INDEX_GUESS];
-    int iMin;
+    int iMin = -1;
+    float xMin = NAN;
+    float chiSquare[N_INDEX_GUESS];
 
     for (int i = 0; i < N_INDEX_GUESS; i++) {
         model->params->data.F32[PM_PAR_7] = indexGuess[i];
-        pmSourceModelGuessPCM (pcm, source, maskVal, markVal);
-
-        pmSourceFitPCM (pcm, source, &options, maskVal);
+
+        model->modelGuess(model, source);
+        pmSourceFitModel (source, model, &options, maskVal);
+
         chiSquare[i] = model->chisq;
         if (i == 0) {
@@ -573 +585 @@
 // 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) {
+bool psphotFitSersicIndexPCM (pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
+
+    pmModel *model = pcm->modelConv;
 
     assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC"));
@@ -583 +597 @@
     options.nIter = 3;
 
-    float xMin, chiSquare[N_INDEX_GUESS];
-    int iMin;
-
-    // XXX we probably cannot be calling model->modelGuess() : this does not include the psf sigma
+    int iMin = -1;
+    float xMin = NAN;
+    float chiSquare[N_INDEX_GUESS];
 
     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);
+        pmSourceFitModel (source, model, &options, maskVal);
+
+        // pmSourceModelGuessPCM(pcm, source, maskVal, markVal);
+        // pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
+
         chiSquare[i] = model->chisq;
         if (i == 0) {
@@ -604 +622 @@
     }
 
+    assert (iMin >= 0);
+    
     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;
-}
+
+    pmSourceModelGuessPCM(pcm, source, maskVal, markVal);
+
+    return true;
+}