IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jul 20, 2010, 7:21:53 PM (16 years ago)
Author:
eugene
Message:

update PCM fitting (not quite ready)

Location:
branches/eam_branches/ipp-20100621/psphot/src
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20100621/psphot/src/Makefile.am

    r28683 r28692  
    168168        psphotExtendedSourceFits.c     \
    169169        psphotKernelFromPSF.c          \
    170         psphotPSFConvModel.c           \
    171170        psphotFitSet.c                 \
    172171        psphotSourceFreePixels.c       \

  • branches/eam_branches/ipp-20100621/psphot/src/psphot.h

    r28643 r28692  
    1212
    1313#define PSPHOT_RECIPE_PSF_FAKE_ALLOW "PSF.FAKE.ALLOW" // Name for recipe component permitting fake PSFs
    14 
    15 // pmPCMData : PSF Convolved Model data storage structure
    16 typedef struct {
    17     psImage *model;
    18     psArray *dmodels;
    19     psImage *modelConv;
    20     psArray *dmodelsConv;
    21 } pmPCMData;
    2214
    2315// top-level psphot functions
    … …  
    293285bool psphotRadialBins (psMetadata *recipe, pmSource *source, float radiusMax, float skynoise);
    294286
    295 // structures & functions to support psf-convolved model fitting
    296 
    297 // psf-convolved model fitting
    298 bool psphotModelWithPSF_LMM (
    299     psMinimization *min,
    300     psImage *covar,
    301     psVector *params,
    302     psMinConstraint *constraint,
    303     pmSource *source,
    304     const psKernel *psf,
    305     psMinimizeLMChi2Func func);
    306 
    307 psF32 psphotModelWithPSF_SetABX(
    308     psImage  *alpha,
    309     psVector *beta,
    310     const psVector *params,
    311     const psVector *paramMask,
    312     pmPCMData *pcm,
    313     const pmSource *source,
    314     const psKernel *psf,
    315     psMinimizeLMChi2Func func);
    316 
    317 pmPCMData *pmPCMDataAlloc (
    318     const psVector *params,
    319     const psVector *paramMask,
    320     pmSource *source);
    321 
    322 psImage *pmPCMDataSaveImage (pmPCMData *pcm);
    323 
    324287int psphotKapaOpen (void);
    325288bool psphotKapaClose (void);
    … …  
    463426bool psphotStackObjectsUnifyPosition (psArray *objects);
    464427
     428bool psphotFitSersicIndexPCM (pmSource *source, pmModel *model, pmSourceFitOptions *fitOptions, psImageMaskType maskVal);
     429pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize);
     430
    465431#endif

  • branches/eam_branches/ipp-20100621/psphot/src/psphotApResid.c

    r28405 r28692  
    459459        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", markVal);
    460460
    461         bool status = pmSourceMagnitudes (source, psf, photMode, maskVal);
     461        bool status = pmSourceMagnitudes (source, psf, photMode, maskVal, markVal);
    462462
    463463        // clear the mask bit

  • branches/eam_branches/ipp-20100621/psphot/src/psphotMagnitudes.c

    r28405 r28692  
    186186        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", markVal);
    187187
    188         status = pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal
     188        status = pmSourceMagnitudes (source, psf, photMode, maskVal, markVal);
    189189        if (status && isfinite(source->apMag)) Nap ++;
    190190

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

    r28687 r28692  
    88static int NfitDBL = 0;
    99static int NfitEXT = 0;
     10static int NfitPCM = 0;
    1011
    1112bool psphotFitInit (int nThreads) {
    … …  
    440441pmModel *psphotFitEXT (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal) {
    441442
     443    if ((source->moments->Mxx < 1e-3) || (source->moments->Myy < 1e-3)) {
     444        psTrace ("psphot", 5, "problem source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy);
     445    }
     446
    442447    pmSourceFitOptions options = *fitOptions;
    443448
    … …  
    456461    }
    457462
    458     if ((source->moments->Mxx < 1e-3) || (source->moments->Myy < 1e-3)) {
    459         psTrace ("psphot", 5, "problem source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy);
    460     }
    461 
    462463    // for sersic models, use a grid search to choose an index, then float the params there
    463464    if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
    … …  
    476477}
    477478
    478 pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal) {
     479pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
     480
     481    if ((source->moments->Mxx < 1e-3) || (source->moments->Myy < 1e-3)) {
     482        psTrace ("psphot", 5, "problem source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy);
     483    }
     484
     485    float radius = psphotSetRadiusEXT (readout, source, markVal);
     486
     487    // XXX note that this changes the source moments that are published...
     488    // recalculate the source moments using the larger extended-source moments radius
     489    // at this stage, skip Gaussian windowing, and do not clip pixels by S/N
     490    // this uses the footprint to judge both radius and aperture?
     491    if (!pmSourceMoments (source, radius, 0.0, 0.0, maskVal)) return false;
    479492
    480493    pmSourceFitOptions options = *fitOptions;
    … …  
    487500    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
    488501
    489     // use the source moments, etc to guess basic model parameters
    490     pmModel *EXT = pmSourceModelGuess (source, modelType);
    491     if (!EXT) {
     502    pmPCMdata *pcm = pmPCMinit (source, &options, modelType, maskVal, psfSize);
     503    if (!pcm) {
    492504        psTrace ("psphot", 5, "failed to generate a model for source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy);
    493505        return NULL;
    494506    }
    495 
    496     if ((source->moments->Mxx < 1e-3) || (source->moments->Myy < 1e-3)) {
    497         psTrace ("psphot", 5, "problem source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy);
    498     }
     507    // XXX check for nDOF too small
     508
     509    // use the source moments, etc to guess basic model parameters
     510    pmSourceModelGuessPCM (pcm, source, maskVal, markVal);
    499511
    500512    // for sersic models, use a grid search to choose an index, then float the params there
    501513    if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
    502         psphotFitSersicIndexPCM (source, EXT, fitOptions, maskVal);
     514        psphotFitSersicIndexPCM (pcm, source, fitOptions, maskVal, markVal);
    503515    }
    504516
    … …  
    508520        options.mode = PM_SOURCE_FIT_EXT;
    509521    }
    510     pmSourceFitPCM (source, EXT, &options, maskVal);
     522    pmSourceFitPCM (source, PCM, &options, maskVal);
    511523
    512524    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
    513     return (EXT);
     525    return (PCM);
    514526}
    515527
    … …  
    520532// A sersic model is very sensitive to the index.  attempt to find the index first by grid search in just the index
    521533// for a sersic model, attempt to fit just the index and normalization with a modest number of iterations
    522 bool psphotFitSersicIndex (pmSource *source, pmModel *model, pmSourceFitOptions *fitOptions, psImageMaskType maskVal) {
     534bool psphotFitSersicIndex (pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) {
     535
     536    pmModel *model = pcm->modelConv;
    523537
    524538    assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC"));
    … …  
    535549    for (int i = 0; i < N_INDEX_GUESS; i++) {
    536550        model->params->data.F32[PM_PAR_7] = indexGuess[i];
    537         model->modelGuess(model, source);
    538         pmSourceFitModel (source, model, &options, maskVal);
     551        pmSourceModelGuessPCM (pcm, source, maskVal, markVal);
     552
     553        pmSourceFitPCM (pcm, source, &options, maskVal);
    539554        chiSquare[i] = model->chisq;
    540555        if (i == 0) {
    … …  
    570585    float xMin, chiSquare[N_INDEX_GUESS];
    571586    int iMin;
     587
     588    // XXX we probably cannot be calling model->modelGuess() : this does not include the psf sigma
    572589
    573590    for (int i = 0; i < N_INDEX_GUESS; i++) {

  • branches/eam_branches/ipp-20100621/psphot/src/psphotSourceSize.c

    r28677 r28692  
    191191
    192192        // XXX can we test if psfMag is set and calculate only if needed?
    193         pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal
     193        pmSourceMagnitudes (source, psf, photMode, maskVal, markVal);
    194194
    195195        // clear the mask bit
    … …  
    342342
    343343        // XXX can we test if psfMag is set and calculate only if needed?
    344         pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal
     344        pmSourceMagnitudes (source, psf, photMode, maskVal, markVal);
    345345
    346346        // clear the mask bit
Note: See TracChangeset for help on using the changeset viewer.