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

Timestamp:

Jul 29, 2010, 2:36:33 PM (16 years ago)

Author:

eugene

Message:

attempting to get a better starting guess; attempting to downweight neighbor detections

File:

• branches/eam_branches/ipp-20100621/psphot/src/psphotSourceFits.c (modified) (20 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
@@ -19 +18 @@
 
     psLogMsg ("psphot.pspsf", PS_LOG_INFO, "fitted %5d psf, %5d blend, %5d ext, %5d dbl : %6.2f sec\n",
-             NfitPSF, NfitBlend, NfitEXT, NfitDBL, psTimerMark ("psphot.fits"));
+              NfitPSF, NfitBlend, NfitEXT, NfitDBL, psTimerMark ("psphot.fits"));
     return true;
 }
@@ -205 +204 @@
 }
 
+// save a local, static copy of the EXT model type so we don't have to lookup for each object
 static pmModelType modelTypeEXT;
 
-bool psphotInitLimitsEXT (psMetadata *recipe) {
+bool psphotInitLimitsEXT (psMetadata *recipe, pmReadout *readout) {
 
     bool status;
@@ -214 +214 @@
     char *modelNameEXT = psMetadataLookupStr (&status, recipe, "EXT_MODEL");
     modelTypeEXT = pmModelClassGetType (modelNameEXT);
-    psphotInitRadiusEXT (recipe, modelTypeEXT);
+
+    psphotInitRadiusEXT (recipe, readout);
 
     return true;
@@ -221 +222 @@
 bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *newSources, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) {
 
+    float radius;
     bool okEXT, okDBL;
     float chiEXT, chiDBL;
@@ -230 +232 @@
 
     // skip the source if we don't think it is extended
+    // XXX are these robust, or do we have better info from the source-size analysis??
     if (source->type == PM_SOURCE_TYPE_UNKNOWN) return false;
     if (source->type == PM_SOURCE_TYPE_DEFECT) return false;
@@ -235 +238 @@
 
     // set the radius based on the footprint (also sets the mask pixels)
-    float radius = psphotSetRadiusEXT (readout, source, markVal);
+    if (!psphotSetRadiusFootprint(&radius, readout, source, markVal, 1.0)) return false;
 
     // XXX note that this changes the source moments that are published...
+    // XXX all published moments should use the same measurement
     // recalculate the source moments using the larger extended-source moments radius
     // at this stage, skip Gaussian windowing, and do not clip pixels by S/N
     // this uses the footprint to judge both radius and aperture?
+    // XXX save the psf-based moments for output
     if (!pmSourceMoments (source, radius, 0.0, 0.0, maskVal)) return false;
 
@@ -250 +255 @@
     // this temporary source is used as a place-holder by the psphotEval functions below
     tmpSrc = pmSourceAlloc ();
-
-    // XXX need to handle failures better here
-    EXT = psphotFitEXT (readout, source, fitOptions, modelTypeEXT, maskVal, markVal);
-    if (!EXT) goto escape;
-    if (!isfinite(EXT->params->data.F32[PM_PAR_I0])) goto escape;
-
-    okEXT = psphotEvalEXT (tmpSrc, EXT);
-    chiEXT = EXT ? EXT->chisq / EXT->nDOF : NAN;
-
-    // DBL will always be defined, but DBL->data[n] might not
-    DBL = psphotFitDBL (readout, source, fitOptions, maskVal, markVal);
-    if (!DBL) goto escape;
-    if (!DBL->n) goto escape;
-
-    okDBL  = psphotEvalDBL (tmpSrc, DBL->data[0]);
-    okDBL &= psphotEvalDBL (tmpSrc, DBL->data[1]);
-    // XXX should I keep / save the flags set in the eval functions?
+    {
+        // DBL will always be defined, but DBL->data[n] might not
+        DBL = psphotFitDBL (readout, source, fitOptions, maskVal, markVal);
+        if (!DBL) goto escape;
+        if (!DBL->n) goto escape;
+
+        okDBL  = psphotEvalDBL (tmpSrc, DBL->data[0]);
+        okDBL &= psphotEvalDBL (tmpSrc, DBL->data[1]);
+        // XXX should I keep / save the flags set in the eval functions?
+
+        // correct first model chisqs for flux trend
+        chiDBL = NAN;
+        ONE = DBL->data[0];
+        if (ONE) {
+            if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
+            chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]);
+            ONE->chisqNorm = ONE->chisq / chiTrend;
+            chiDBL = ONE->chisq / ONE->nDOF; // save chisq for double-star/galaxy comparison
+            ONE->fitRadius = radius;
+        }
+
+        // correct second model chisqs for flux trend
+        ONE = DBL->data[1];
+        if (ONE) {
+            if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
+            chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]);
+            ONE->chisqNorm = ONE->chisq / chiTrend;
+            ONE->fitRadius = radius;
+        }
+    }
+
+    { 
+        // XXX need to handle failures better here
+        EXT = psphotFitEXT (readout, source, fitOptions, modelTypeEXT, maskVal, markVal);
+        if (!EXT) goto escape;
+        if (!isfinite(EXT->params->data.F32[PM_PAR_I0])) goto escape;
+
+        okEXT = psphotEvalEXT (tmpSrc, EXT);
+        chiEXT = EXT ? EXT->chisq / EXT->nDOF : NAN;
+    }
 
     // clear the circular mask
     psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); 
-
-    // correct first model chisqs for flux trend
-    chiDBL = NAN;
-    ONE = DBL->data[0];
-    if (ONE) {
-        if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
-      chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]);
-      ONE->chisqNorm = ONE->chisq / chiTrend;
-      chiDBL = ONE->chisq / ONE->nDOF; // save chisq for double-star/galaxy comparison
-    }
-
-    // correct second model chisqs for flux trend
-    ONE = DBL->data[1];
-    if (ONE) {
-        if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
-      chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]);
-      ONE->chisqNorm = ONE->chisq / chiTrend;
-    }
 
     psFree (tmpSrc);
@@ -314 +324 @@
 
     // save new model
+    // XXX save the correct radius...
     source->modelEXT = EXT;
-    source->modelEXT->fitRadius = radius;
     source->type = PM_SOURCE_TYPE_EXTENDED;
     source->mode |= PM_SOURCE_MODE_EXTMODEL;
@@ -343 +353 @@
     source->modelPSF = psMemIncrRefCounter (DBL->data[0]);
     source->mode     |= PM_SOURCE_MODE_PAIR;
-    source->modelPSF->fitRadius = radius;
 
     // copy most data from the primary source (modelEXT, blends stay NULL)
     pmSource *newSrc = pmSourceCopyData (source);
     newSrc->modelPSF = psMemIncrRefCounter (DBL->data[1]);
-    newSrc->modelPSF->fitRadius = radius;
 
     // build cached models and subtract
@@ -452 +460 @@
     maskVal |= markVal;
 
-    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
-
     // use the source moments, etc to guess basic model parameters
-    pmModel *EXT = pmSourceModelGuess (source, modelType);
-    if (!EXT) {
+    pmModel *model = pmSourceModelGuess (source, modelType);
+    if (!model) {
         psTrace ("psphot", 5, "failed to generate a model for source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy);
         return NULL;
@@ -463 +469 @@
     // for sersic models, use a grid search to choose an index, then float the params there
     if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
-        psphotFitSersicIndex (source, EXT, fitOptions, maskVal);
+        // for the test fits, use a somewhat smaller radius
+        psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5);
+        psphotFitSersicIndex (model, readout, source, fitOptions, maskVal, markVal);
+    }
+
+    if (!psphotSetRadiusModel (model, readout, source, markVal, true)) {
+        psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 1.0);
     }
 
@@ -471 +483 @@
         options.mode = PM_SOURCE_FIT_EXT;
     }
-    pmSourceFitModel (source, EXT, &options, maskVal);
+
+    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
+    pmSourceFitModel (source, model, &options, maskVal);
+    fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
 
     // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
-    return (EXT);
+    return (model);
 }
 
 pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
+
+    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);
+    }
+
+    pmSourceFitOptions options = *fitOptions;
+
+    NfitPCM ++;
+
+    // maskVal is used to test for rejected pixels, and must include markVal
+    maskVal |= markVal;
 
     // allocate the model
@@ -484 +510 @@
         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);
-    }
-
-    float radius = psphotSetRadiusEXT (readout, source, markVal);
-
-    // XXX note that this changes the source moments that are published...
-    // recalculate the source moments using the larger extended-source moments radius
-    // 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)) {
-        // XXX set some mask bit/
-        model->flags |= PM_MODEL_STATUS_BADARGS;
-        return model;
-    }
-
-    NfitPCM ++;
-
-    // maskVal is used to test for rejected pixels, and must include markVal
-    maskVal |= markVal;
-
-    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
 
     pmPCMdata *pcm = pmPCMinit (source, &options, model, maskVal, psfSize);
@@ -518 +519 @@
 
     // use the source moments, etc to guess basic model parameters
-    pmSourceModelGuessPCM (pcm, source, maskVal, markVal);
+    if (!pmSourceModelGuessPCM (pcm, source, maskVal, markVal)) {
+        psFree (pcm);
+        model->flags |= PM_MODEL_STATUS_BADARGS;
+        return model;
+    }
 
     // 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, psfSize);
+        // for the test fits, use a somewhat smaller radius
+        psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5);
+        if (!psphotFitSersicIndexPCM (pcm, readout, source, fitOptions, maskVal, markVal, psfSize)) {
+            model->flags |= PM_MODEL_STATUS_BADARGS;
+            return model;
+        }
+    }
+
+    if (!psphotSetRadiusModel (model, readout, source, markVal, true)) {
+        psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 1.0);
     }
 
@@ -530 +544 @@
         options.mode = PM_SOURCE_FIT_EXT;
     }
+
+    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
     pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
+    fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
 
     // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
@@ -544 +561 @@
 // 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 (pmSource *source, pmModel *model, pmSourceFitOptions *fitOptions, psImageMaskType maskVal) {
+bool psphotFitSersicIndex (pmModel *model, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) {
 
     assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC"));
@@ -552 +569 @@
     // fit EXT (not PSF) model (set/unset the pixel mask)
     options.mode = PM_SOURCE_FIT_NO_INDEX;
-    options.nIter = 3;
+    options.nIter = 4;
 
     int iMin = -1;
@@ -561 +578 @@
         model->params->data.F32[PM_PAR_7] = indexGuess[i];
 
-        model->modelGuess(model, source);
+        if (!model->modelGuess(model, source)) {
+            model->flags |= PM_MODEL_STATUS_BADARGS;
+            return false;
+        }
+
+        // each time we change the model guess, we need to adjust the radius
+        // XXX this did not work : we do not need such a large radius -- just uses moments-based radius
+        if (false && !psphotSetRadiusModel (model, readout, source, markVal, false)) {
+            psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5);
+        }
+        
         pmSourceFitModel (source, model, &options, maskVal);
+        fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
+
+        chiSquare[i] = model->chisqNorm;
+        if (i == 0) {
+            xMin = chiSquare[i];
+            iMin = i;
+        } else {
+            if (chiSquare[i] < xMin) {
+                xMin = chiSquare[i];
+                iMin = i;
+            }
+        }
+    }
+    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);
+
+    // each time we change the model guess, we need to adjust the radius
+    // if (!psphotSetRadiusModel (model, readout, source, markVal, true)) {
+    //  psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal);
+    // }
+
+    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 (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
+
+    pmModel *model = pcm->modelConv;
+
+    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 = 4;
+
+    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];
+        
+        if (!model->modelGuess(model, source)) {
+            model->flags |= PM_MODEL_STATUS_BADARGS;
+            return false;
+        }
+
+        // each time we change the model guess, we need to adjust the radius
+        // XXX this did not work : we do not need such a large radius -- just uses moments-based radius
+        if (false && !psphotSetRadiusModel (model, readout, source, markVal, false)) {
+            psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5);
+        }
+        
+        pmSourceFitModel (source, model, &options, maskVal);
+        fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
+
+        // pmSourceModelGuessPCM(pcm, source, maskVal, markVal);
+        // pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
 
         chiSquare[i] = model->chisq;
@@ -575 +666 @@
         }
     }
-
-    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;
-}
-
-// 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 (pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
-
-    pmModel *model = pcm->modelConv;
-
-    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;
-
-    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);
-        pmSourceFitModel (source, model, &options, maskVal);
-
-        // pmSourceModelGuessPCM(pcm, source, maskVal, markVal);
-        // pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
-
-        chiSquare[i] = model->chisq;
-        if (i == 0) {
-            xMin = chiSquare[i];
-            iMin = i;
-        } else {
-            if (chiSquare[i] < xMin) {
-                xMin = chiSquare[i];
-                iMin = i;
-            }
-        }
-    }
-
     assert (iMin >= 0);