Changeset 29060 for branches/sc_branches/trunkTest/psphot/src/psphotSourceFits.c

Timestamp:

Aug 26, 2010, 9:18:39 AM (16 years ago)

Author:

Serge CHASTEL

Message:

Merging trunk in branch

Location:

branches/sc_branches/trunkTest

Files:

• . (modified) (1 prop)
• psphot (modified) (1 prop)
• psphot/src/psphotSourceFits.c (modified) (23 diffs)

branches/sc_branches/trunkTest/psphot/src/psphotSourceFits.c

@@ -2 +2 @@
 
 // given a source with an existing modelPSF, attempt a full PSF fit, subtract if successful
-// XXX this function does not call pmSourceFitModelInit : fix this?
 
 static int NfitPSF = 0;
@@ -8 +7 @@
 static int NfitDBL = 0;
 static int NfitEXT = 0;
+static int NfitPCM = 0;
 
 bool psphotFitInit (int nThreads) {
@@ -18 +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"));
-    return true;
-}
-
-bool psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf, psImageMaskType maskVal, psImageMaskType markVal) {
+              NfitPSF, NfitBlend, NfitEXT, NfitDBL, psTimerMark ("psphot.fits"));
+    return true;
+}
+
+bool psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) {
 
     float x, y, dR;
+
+    pmSourceFitOptions options = *fitOptions;
 
     // maskVal is used to test for rejected pixels, and must include markVal
@@ -31 +33 @@
     // if this source is not a possible blend, just fit as PSF
     if ((source->blends == NULL) || (source->mode & PM_SOURCE_MODE_SATSTAR)) {
-        bool status = psphotFitPSF (readout, source, psf, maskVal, markVal);
+        bool status = psphotFitPSF (readout, source, psf, fitOptions, maskVal, markVal);
         return status;
     }
@@ -91 +93 @@
 
     // fit PSF model
-    pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF, maskVal);
+    options.mode = PM_SOURCE_FIT_PSF;
+    pmSourceFitSet (source, modelSet, &options, maskVal);
 
     // clear the circular mask
@@ -154 +157 @@
 }
 
-bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf, psImageMaskType maskVal, psImageMaskType markVal) {
+bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) {
 
     double chiTrend;
+    pmSourceFitOptions options = *fitOptions;
 
     // maskVal is used to test for rejected pixels, and must include markVal
@@ -171 +175 @@
 
     // fit PSF model (set/unset the pixel mask)
-    pmSourceFitModel (source, PSF, PM_SOURCE_FIT_PSF, maskVal);
+    options.mode = PM_SOURCE_FIT_PSF;
+    pmSourceFitModel (source, PSF, &options, maskVal);
 
     if (!isfinite(PSF->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
@@ -199 +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;
@@ -208 +214 @@
     char *modelNameEXT = psMetadataLookupStr (&status, recipe, "EXT_MODEL");
     modelTypeEXT = pmModelClassGetType (modelNameEXT);
-    psphotInitRadiusEXT (recipe, modelTypeEXT);
-
-    return true;
-}
-
-bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *newSources, pmPSF *psf, psImageMaskType maskVal, psImageMaskType markVal) {
-
+
+    psphotInitRadiusEXT (recipe, readout);
+
+    return true;
+}
+
+bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *newSources, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) {
+
+    float radius;
     bool okEXT, okDBL;
     float chiEXT, chiDBL;
@@ -222 +230 @@
     pmModel *EXT = NULL;
     psArray *DBL = NULL;
+    pmMoments psfMoments;
 
     // 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;
@@ -229 +239 @@
 
     // 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?
-    if (!pmSourceMoments (source, radius, 0.0, 0.0, maskVal)) return false;
+    // XXX save the psf-based moments for output
+    psfMoments = *source->moments;
+    if (!pmSourceMoments (source, radius, 0.0, 0.5, maskVal)) {
+      *source->moments = psfMoments;
+      return false;
+    }
 
     psTrace ("psphot", 5, "trying blob...\n");
@@ -244 +260 @@
     // 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, 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, 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);
@@ -298 +319 @@
 
     // both models failed; reject them both
-    // XXX -- change type flags to psf in this case and keep original moments?
+    // XXX -- change type flags to psf in this case, and make sure we subtract it?
+    // reset the psf moments
+    *source->moments = psfMoments;
+
     psFree (EXT);
     psFree (DBL);
@@ -308 +332 @@
 
     // 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;
@@ -327 +351 @@
 # endif
 
+    // reset the psf moments
+    *source->moments = psfMoments;
     return true;
 
@@ -337 +363 @@
     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 = pmSourceCopy (source);
+    pmSource *newSrc = pmSourceCopyData (source);
     newSrc->modelPSF = psMemIncrRefCounter (DBL->data[1]);
-    newSrc->modelPSF->fitRadius = radius;
 
     // build cached models and subtract
@@ -365 +389 @@
 # endif
 
+    // reset the (original) psf moments
+    *source->moments = psfMoments;
+    *newSrc->moments = psfMoments;
+
     psArrayAdd (newSources, 100, newSrc);
     psFree (newSrc);
@@ -371 +399 @@
 
 escape:
+    // reset the psf moments
+    *source->moments = psfMoments;
     psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); 
     psFree (tmpSrc);
@@ -379 +409 @@
 
 // fit a double PSF source to an extended blob
-psArray *psphotFitDBL (pmReadout *readout, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal) {
+psArray *psphotFitDBL (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) {
 
     float dx, dy;
@@ -387 +417 @@
     psEllipseMoments moments;
     psArray *modelSet;
+
+    pmSourceFitOptions options = *fitOptions;
 
     NfitDBL ++;
@@ -403 +435 @@
     axes = psEllipseMomentsToAxes (moments, 20.0);
 
+    if (isnan(axes.major)) return NULL;
+    if (isnan(axes.minor)) return NULL;
+    if (isnan(axes.theta)) return NULL;
+
     // XXX this is really arbitrary: 4 pixel separation?
     dx = 2 * cos (axes.theta);
@@ -426 +462 @@
 
     // fit PSF model (set/unset the pixel mask)
-    pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF, maskVal);
+    options.mode = PM_SOURCE_FIT_PSF;
+    pmSourceFitSet (source, modelSet, &options, maskVal);
     return (modelSet);
 }
 
-pmModel *psphotFitEXT (pmReadout *readout, pmSource *source, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal) {
+pmModel *psphotFitEXT (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal) {
+
+    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;
 
     NfitEXT ++;
@@ -438 +481 @@
 
     // 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;
     }
 
+    // for sersic models, use a grid search to choose an index, then float the params there
+    if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
+        // 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);
+    }
+
+    if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
+        options.mode = PM_SOURCE_FIT_NO_INDEX;
+    } else {
+        options.mode = PM_SOURCE_FIT_EXT;
+    }
+
+    // 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 (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
+    pmModel *model = pmModelAlloc(modelType);
+    if (!model) {
+        return NULL;
+    }
+
+    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);
+        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
+    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")) {
+        // 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);
+    }
+
+    if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
+        options.mode = PM_SOURCE_FIT_NO_INDEX;
+    } else {
+        options.mode = PM_SOURCE_FIT_EXT;
+    }
+    // update the pcm elements if we have changed the circumstance (options.mode or source->pixels)
+    pmPCMupdate(pcm, source, &options, model);
+
+    // 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);
+    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
+
+// 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 (pmModel *model, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) {
+
+    assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC"));
+
+    pmSourceFitOptions options = *fitOptions;
+    
     // fit EXT (not PSF) model (set/unset the pixel mask)
-    pmSourceFitModel (source, EXT, PM_SOURCE_FIT_EXT, maskVal);
-    return (EXT);
-}
+    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);
+
+        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;
+
+    // update the pcm elements if we have changed the circumstance (here, options.mode)
+    pmPCMupdate(pcm, source, &options, model);
+
+    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);
+        }
+        
+# if (0)
+        pmSourceFitModel (source, model, &options, maskVal);
+# else
+        pmSourceModelGuessPCM(pcm, source, maskVal, markVal);
+        pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
+# endif
+        // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
+
+        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);
+    
+    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];
+
+    pmSourceModelGuessPCM(pcm, source, maskVal, markVal);
+
+    return true;
+}