IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 28702


Ignore:
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)

Location:
branches/eam_branches/ipp-20100621
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20100621/ippconfig/gpc1/ppMerge.config

    r26902 r28702  
    6060# END
    6161
     62# DARK.ORDINATES        METADATA
     63#       TEMP    METADATA
     64#               ORDER   S32     1
     65#               RULE    STR     CHIP.TEMP
     66#       END
     67#       TEMP2   METADATA
     68#               ORDER   S32     1
     69#               RULE    STR     CHIP.TEMP * CHIP.TEMP
     70#       END
     71# END
     72
    6273# Ordinates for fitting dark current as function of darktime and fpa temp:
    6374# Counts = C0 + (CT0 + CT1*Temp)*exptime

  • branches/eam_branches/ipp-20100621/ippconfig/recipes/psphot.config

    r28687 r28702  
    9393PSF_FIT_MAX_TOL                     F32   2.00            # Fit tolerance for PSF
    9494
    95 PSF_FIT_MIN_VALID_FLUX              F32   1e-8            # minimum allow flux for fitted source
    96 PSF_FIT_MAX_VALID_FLUX              F32   1e+8            # maximum allow flux for fitted source
     95PSF_FIT_MIN_VALID_FLUX              F32  -100000000.0     # minimum allow flux for fitted source
     96PSF_FIT_MAX_VALID_FLUX              F32  +100000000.0     # maximum allow flux for fitted source
    9797
    9898# the following is used to require a minimum quality of fit before

  • branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCM_MinimizeChisq.c

    r28692 r28702  
    291291
    292292    // convolve model and dmodel arrays with PSF
     293    // XXX speed this up by saving the FFTed psf (for each source, obviously)
     294
    293295    psImageConvolveDirect (pcm->modelConvFlux, pcm->modelFlux, pcm->psf);
    294296    for (int n = 0; n < pcm->dmodelsFlux->n; n++) {

  • branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCMdata.c

    r28692 r28702  
    141141}
    142142
    143 pmPCMdata *pmPCMinit(pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, float psfSize) {
    144 
    145     // make sure we savep a cached copy of the psf flux
     143pmPCMdata *pmPCMinit(pmSource *source, pmSourceFitOptions *fitOptions, pmModel *model, psImageMaskType maskVal, float psfSize) {
     144
     145    // make sure we save a cached copy of the psf flux
    146146    pmSourceCachePSF (source, maskVal);
    147147
    148148    // convert the cached cached psf model for this source to a psKernel
    149149    psKernel *psf = pmPCMkernelFromPSF (source, psfSize);
    150     if (!psf) return NULL;
     150    if (!psf) {
     151        // NOTE: this only happens if the source is too close to an edge
     152        model->flags |= PM_MODEL_STATUS_BADARGS;
     153        return NULL;
     154    }
    151155
    152156# if (USE_DELTA_PSF)
    … …  
    154158    psf->image->data.F32[(int)(0.5*psf->image->numRows)][(int)(0.5*psf->image->numCols)] = 1.0;
    155159# endif
    156 
    157     // allocate the model
    158     pmModel *modelConv = pmModelAlloc(modelType);
    159     if (!modelConv) {
    160         psFree (psf);
    161         return NULL;
    162     }
    163160
    164161    // count the number of unmasked pixels:
    … …  
    183180    }   
    184181
    185     psVector *params  = modelConv->params;
     182    psVector *params  = model->params;
    186183
    187184    // create the minimization constraints
    188185    psMinConstraint *constraint = psMinConstraintAlloc();
    189186    constraint->paramMask = psVectorAlloc (params->n, PS_TYPE_VECTOR_MASK);
    190     constraint->checkLimits = modelConv->modelLimits;
     187    constraint->checkLimits = model->modelLimits;
    191188
    192189    // set parameter mask based on fitting mode
    … …  
    239236    }
    240237
     238    if (nPix <  nParams + 1) {
     239        psTrace ("psModules.objects", 4, "insufficient valid pixels\n");
     240        psFree (psf);
     241        psFree (constraint);
     242        model->flags |= PM_MODEL_STATUS_BADARGS;
     243        return NULL;
     244    }
     245
    241246    // generate PCM data storage structure
    242247    pmPCMdata *pcm = pmPCMdataAlloc (params, constraint->paramMask, source);
    243248
    244     pcm->modelConv = modelConv;
    245249    pcm->psf = psf;
     250    pcm->modelConv = psMemIncrRefCounter(model);
    246251    pcm->constraint = constraint;
    247252
    248     if (nPix <  nParams + 1) {
    249         psTrace ("psModules.objects", 4, "insufficient valid pixels\n");
    250         pcm->modelConv->flags |= PM_MODEL_STATUS_BADARGS;
    251         abort ();
    252         // XXX This should not be an abort!!
    253     }
    254253    pcm->nPix = nPix;
    255254    pcm->nDOF = nPix - nParams - 1;

  • branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCMdata.h

    r28692 r28702  
    6161    pmSource *source);
    6262
    63 pmPCMdata *pmPCMinit(pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, float psfSize);
     63pmPCMdata *pmPCMinit(pmSource *source, pmSourceFitOptions *fitOptions, pmModel *model, psImageMaskType maskVal, float psfSize);
    6464
    6565psImage *pmPCMdataSaveImage (pmPCMdata *pcm);

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

    r28692 r28702  
    163163        psphotOutput.c                 \
    164164        psphotFakeSources.c            \
    165         psphotModelWithPSF.c           \
    166165        psphotExtendedSourceAnalysis.c \
    167166        psphotExtendedSourceAnalysisByObject.c \

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

    r28692 r28702  
    426426bool psphotStackObjectsUnifyPosition (psArray *objects);
    427427
    428 bool psphotFitSersicIndexPCM (pmSource *source, pmModel *model, pmSourceFitOptions *fitOptions, psImageMaskType maskVal);
     428bool psphotFitSersicIndexPCM (pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize);
    429429pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize);
    430430

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

    r28686 r28702  
    102102    }
    103103   
    104     float MIN_VALID_FLUX = psMetadataLookupBool(&status, recipe, "PSF_FIT_MIN_VALID_FLUX");
     104    float MIN_VALID_FLUX = psMetadataLookupF32(&status, recipe, "PSF_FIT_MIN_VALID_FLUX");
    105105    if (!status) {
    106106        MIN_VALID_FLUX = 1e-8;
    107107    }
    108     float MAX_VALID_FLUX = psMetadataLookupBool(&status, recipe, "PSF_FIT_MAX_VALID_FLUX");
    109     if (!status) {
    110         MAX_VALID_FLUX = 1e-8;
     108    float MAX_VALID_FLUX = psMetadataLookupF32(&status, recipe, "PSF_FIT_MAX_VALID_FLUX");
     109    if (!status) {
     110        MAX_VALID_FLUX = 1e+8;
    111111    }
    112112

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

    r28692 r28702  
    479479pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
    480480
     481    // allocate the model
     482    pmModel *model = pmModelAlloc(modelType);
     483    if (!model) {
     484        return NULL;
     485    }
     486
     487    pmSourceFitOptions options = *fitOptions;
     488
    481489    if ((source->moments->Mxx < 1e-3) || (source->moments->Myy < 1e-3)) {
    482490        psTrace ("psphot", 5, "problem source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy);
    … …  
    489497    // at this stage, skip Gaussian windowing, and do not clip pixels by S/N
    490498    // this uses the footprint to judge both radius and aperture?
    491     if (!pmSourceMoments (source, radius, 0.0, 0.0, maskVal)) return false;
    492 
    493     pmSourceFitOptions options = *fitOptions;
     499    if (!pmSourceMoments (source, radius, 0.0, 0.0, maskVal)) {
     500        // XXX set some mask bit/
     501        model->flags |= PM_MODEL_STATUS_BADARGS;
     502        return model;
     503    }
    494504
    495505    NfitPCM ++;
    … …  
    500510    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
    501511
    502     pmPCMdata *pcm = pmPCMinit (source, &options, modelType, maskVal, psfSize);
     512    pmPCMdata *pcm = pmPCMinit (source, &options, model, maskVal, psfSize);
    503513    if (!pcm) {
    504514        psTrace ("psphot", 5, "failed to generate a model for source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy);
    505         return NULL;
    506     }
    507     // XXX check for nDOF too small
     515        model->flags |= PM_MODEL_STATUS_BADARGS; // XXX this is probably already set in pmPCMinit
     516        return model;
     517    }
    508518
    509519    // use the source moments, etc to guess basic model parameters
    … …  
    512522    // for sersic models, use a grid search to choose an index, then float the params there
    513523    if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
    514         psphotFitSersicIndexPCM (pcm, source, fitOptions, maskVal, markVal);
     524        psphotFitSersicIndexPCM (pcm, source, fitOptions, maskVal, markVal, psfSize);
    515525    }
    516526
    … …  
    520530        options.mode = PM_SOURCE_FIT_EXT;
    521531    }
    522     pmSourceFitPCM (source, PCM, &options, maskVal);
     532    pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
    523533
    524534    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
    525     return (PCM);
     535    psFree (pcm);
     536
     537    return model;
    526538}
    527539
    528540// note that these should be 1/2n of the standard sersic index
    529541float indexGuess[] = {0.5, 0.33, 0.25, 0.167, 0.125, 0.083};
    530 define N_INDEX_GUESS 6
     542# define N_INDEX_GUESS 6
    531543
    532544// A sersic model is very sensitive to the index.  attempt to find the index first by grid search in just the index
    533545// for a sersic model, attempt to fit just the index and normalization with a modest number of iterations
    534 bool psphotFitSersicIndex (pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) {
    535 
    536     pmModel *model = pcm->modelConv;
     546bool psphotFitSersicIndex (pmSource *source, pmModel *model, pmSourceFitOptions *fitOptions, psImageMaskType maskVal) {
    537547
    538548    assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC"));
    … …  
    544554    options.nIter = 3;
    545555
    546     float xMin, chiSquare[N_INDEX_GUESS];
    547     int iMin;
     556    int iMin = -1;
     557    float xMin = NAN;
     558    float chiSquare[N_INDEX_GUESS];
    548559
    549560    for (int i = 0; i < N_INDEX_GUESS; i++) {
    550561        model->params->data.F32[PM_PAR_7] = indexGuess[i];
    551         pmSourceModelGuessPCM (pcm, source, maskVal, markVal);
    552 
    553         pmSourceFitPCM (pcm, source, &options, maskVal);
     562
     563        model->modelGuess(model, source);
     564        pmSourceFitModel (source, model, &options, maskVal);
     565
    554566        chiSquare[i] = model->chisq;
    555567        if (i == 0) {
    … …  
    573585// A sersic model is very sensitive to the index.  attempt to find the index first by grid search in just the index
    574586// for a sersic model, attempt to fit just the index and normalization with a modest number of iterations
    575 bool psphotFitSersicIndexPCM (pmSource *source, pmModel *model, pmSourceFitOptions *fitOptions, psImageMaskType maskVal) {
     587bool psphotFitSersicIndexPCM (pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
     588
     589    pmModel *model = pcm->modelConv;
    576590
    577591    assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC"));
    … …  
    583597    options.nIter = 3;
    584598
    585     float xMin, chiSquare[N_INDEX_GUESS];
    586     int iMin;
    587 
    588     // XXX we probably cannot be calling model->modelGuess() : this does not include the psf sigma
     599    int iMin = -1;
     600    float xMin = NAN;
     601    float chiSquare[N_INDEX_GUESS];
    589602
    590603    for (int i = 0; i < N_INDEX_GUESS; i++) {
    591604        model->params->data.F32[PM_PAR_7] = indexGuess[i];
     605       
    592606        model->modelGuess(model, source);
    593         pmSourceFitPCM (source, model, &options, maskVal);
     607        pmSourceFitModel (source, model, &options, maskVal);
     608
     609        // pmSourceModelGuessPCM(pcm, source, maskVal, markVal);
     610        // pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
     611
    594612        chiSquare[i] = model->chisq;
    595613        if (i == 0) {
    … …  
    604622    }
    605623
     624    assert (iMin >= 0);
     625   
    606626    model->flags = PM_MODEL_STATUS_NONE; // do not attempt to handle failures here, let the next iteration deal with it
    607627    model->params->data.F32[PM_PAR_7] = indexGuess[iMin];
    608     model->modelGuess(model, source);
    609 
    610     return true;
    611 }
     628
     629    pmSourceModelGuessPCM(pcm, source, maskVal, markVal);
     630
     631    return true;
     632}

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

    r28692 r28702  
    170170    psVector *ApErr = psVectorAllocEmpty (100, PS_TYPE_F32);
    171171
     172    psImageMaskType markVal = options->markVal;
    172173    psImageMaskType maskVal = options->maskVal | options->markVal;
    173174
    … …  
    295296    pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT;
    296297
     298    psImageMaskType markVal = options->markVal;
    297299    psImageMaskType maskVal = options->maskVal | options->markVal;
    298300
    … …  
    364366        float nSigmaMYY = (Myy - psfClump->Y) / hypot(psfClump->dY, psfClump->Y*psfClump->Y*source->errMag);
    365367
    366         fprintf (stderr, "%f %f : Mxx: %f, Myy: %f, dx: %f, dy: %f, psfMag: %f, apMag: %f, dMag: %f, errMag: %f, nSigmaMag: %f\n",
     368        fprintf (stderr, "%f %f : Mxx: %f, Myy: %f, dx: %f, dy: %f, psfMag: %f, apMag: %f, dMag: %f, errMag: %f, nSigmaMag: %f, nSigmaMxx: %f, nSigmaMyy: %f\n",
    367369                 source->peak->xf, source->peak->yf, Mxx, Myy, source->peak->xf - source->moments->Mx, source->peak->yf - source->moments->My,
    368                  source->psfMag, apMag, dMag, source->errMag, nSigmaMAG);
     370                 source->psfMag, apMag, dMag, source->errMag, nSigmaMAG, nSigmaMXX, nSigmaMYY);
    369371
    370372        // partially-masked sources are more likely to be mis-measured PSFs
Note: See TracChangeset for help on using the changeset viewer.