IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Sep 6, 2011, 1:32:31 PM (15 years ago)
Author:
eugene
Message:
  • clarified the output to show the stages and their timing more clearly
  • psphotStackReadout now uses the same functions as psphotReadout (deprecated FitLinearStack, ExtendedAnalysisByObject, RadialAnalysisByObject, etc)
  • psphotStack now does 2nd (5 sigma) detection pass and radial photometry of the un-matched image (as well as the matched images)
  • for speed, fitting uses a smaller subset of pixels than subtraction: the model is extended when the subtraction is performed to avoid a discontinuity
  • extended source analysis (petrosians) is threaded
  • iterative kron radius / magnitude analysis with down-weighted neighbors replaces the masked neighbor kron analysis. the result of this analysis is used by fitting guesses and is the value written to the PSF and XFIT tables.
  • the target matched PSFs are determined in advance and the 0 element of this array is always the unmatched image (with target PSF of NAN)
  • dQ, the difference between the upper and lower quartile, is measured (but not yet used to adjust the sky level. sky should be FITTED_MEAN - 2.5*dQ
  • fake sources inserted for the efficiency analysis are removed at the end of the analysis
  • radial apertures is now threaded
  • edge case petrosian mags are now handled without aborts (may still have a problem with the measured error of the petrosian radius)
  • new sersic model guess based on the psf size, the moments, and the kron flux. the psf - kron mag is used to guess the index, then the index is used to convert measured moments into model guess parameters. this uses empirical relationships between the quantities generated from fake images.
  • define the fit region and window region based on the 1st radial moments
  • use the kron S/N to determine if analysis is performed for extended sources, not the psf S/N
  • remove the old 'basic deblend' which did not really do any deblending and is superceeded by the deblending in KronIterate and the footprint analysis
Location:
trunk/psphot
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot

  • trunk/psphot/src

    • Property svn:ignore

      •  

        old new  
        2222psphotMakePSF
        2323psphotStack
         24psphotModelTest

  • trunk/psphot/src/psphotSourceFits.c

    r31452 r32348  
    211211bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *newSources, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) {
    212212
    213     float radius;
     213    float fitRadius, windowRadius;
    214214    bool okEXT, okDBL;
    215215    pmModel *ONE = NULL;
    … …  
    217217    pmModel *EXT = NULL;
    218218    psArray *DBL = NULL;
    219     pmMoments psfMoments;
     219    // pmMoments psfMoments;
    220220
    221221    // skip the source if we don't think it is extended
    … …  
    225225    if (source->type == PM_SOURCE_TYPE_SATURATED) return false;
    226226
     227# define TEST_X -420.0
     228# define TEST_Y 300.0
     229   
     230    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
     231        fprintf (stderr, "test galaxy\n");
     232    }
     233
     234# undef TEST_X
     235# undef TEST_Y
     236
    227237    // set the radius based on the footprint (also sets the mask pixels)
    228     if (!psphotSetRadiusFootprint(&radius, readout, source, markVal, 1.0)) return false;
    229 
    230     // XXX note that this changes the source moments that are published...
    231     // XXX all published moments should use the same measurement
    232     // recalculate the source moments using the larger extended-source moments radius
    233     // at this stage, skip Gaussian windowing, and do not clip pixels by S/N
    234     // this uses the footprint to judge both radius and aperture?
    235     // XXX save the psf-based moments for output
    236     psfMoments = *source->moments;
    237     if (!pmSourceMoments (source, radius, 0.0, 0.5, 0.0, maskVal)) {
    238       *source->moments = psfMoments;
    239       return false;
    240     }
     238    if (!psphotSetRadiusMoments(&fitRadius, &windowRadius, readout, source, markVal)) return false;
     239    // fprintf (stderr, "rad: %6.1f %6.1f  | %5.2f %5.2f %5.2f  ", source->peak->xf, source->peak->yf, source->moments->Mrf, fitRadius, windowRadius);
    241240
    242241    psTrace ("psphot", 5, "trying blob...\n");
    … …  
    263262        ONE = DBL->data[0];
    264263        if (ONE) {
    265             if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
     264            psAssert (isfinite(ONE->params->data.F32[PM_PAR_I0]), "nan in fit");
    266265            chiDBL = ONE->chisqNorm; // save chisq for double-star/galaxy comparison
    267             ONE->fitRadius = radius;
     266            ONE->fitRadius = fitRadius;
    268267        }
    269268
    … …  
    271270        ONE = DBL->data[1];
    272271        if (ONE) {
    273             if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
    274             ONE->fitRadius = radius;
     272            psAssert (isfinite(ONE->params->data.F32[PM_PAR_I0]), "nan in fit");
     273            ONE->fitRadius = fitRadius;
    275274        }
    276275    }
    … …  
    284283        okEXT = psphotEvalEXT (tmpSrc, EXT);
    285284        chiEXT = EXT ? EXT->chisqNorm : NAN;
     285        EXT->fitRadius = fitRadius;
    286286    }
    287287
    … …  
    294294
    295295    if (okEXT && okDBL) {
    296         psTrace ("psphot", 5, "blob chisq: %f vs %f\n", chiEXT, chiDBL);
    297296        // XXX EAM : a bogus bias: need to examine this better
    298297        if (3*chiEXT > chiDBL) goto keepDBL;
    … …  
    303302    if (!okEXT && okDBL) goto keepDBL;
    304303
     304    psTrace ("psphot", 4, "both failed: blob chisq: %f vs %f for %f,%f\n", chiEXT, chiDBL, source->peak->xf, source->peak->yf);
     305
    305306    // both models failed; reject them both
    306307    // XXX -- change type flags to psf in this case, and make sure we subtract it?
    307308    // reset the psf moments
    308     *source->moments = psfMoments;
     309    // XXX *source->moments = psfMoments;
    309310
    310311    psFree (EXT);
    … …  
    313314
    314315keepEXT:
     316    psTrace ("psphot", 4, "goto EXT : blob chisq: %f vs %f for %f,%f\n", chiEXT, chiDBL, source->peak->xf, source->peak->yf);
    315317    // sub EXT
    316318    psFree (DBL);
    … …  
    338340
    339341    // reset the psf moments
    340     *source->moments = psfMoments;
     342    // XXX *source->moments = psfMoments;
    341343    return true;
    342344
    343345keepDBL:
     346    psTrace ("psphot", 4, "goto DBL : blob chisq: %f vs %f for %f,%f\n", chiEXT, chiDBL, source->peak->xf, source->peak->yf);
    344347    // sub DLB
    345348    psFree (EXT);
    … …  
    372375            psLogMsg ("psphot", 1, "PAR %d : %f +/- %f\n", i, source->modelPSF->params->data.F32[i], source->modelPSF->dparams->data.F32[i]);
    373376        }
    374         psphotVisualShowResidualImage (readout, false);
    375377    }
    376378# endif
    377 
    378     // reset the (original) psf moments
    379     *source->moments = psfMoments;
    380     *newSrc->moments = psfMoments;
    381379
    382380    psArrayAdd (newSources, 100, newSrc);
    … …  
    386384
    387385escape:
    388     // reset the psf moments
    389     *source->moments = psfMoments;
    390386    psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
    391387    psFree (tmpSrc);
    … …  
    476472    // for sersic models, use a grid search to choose an index, then float the params there
    477473    if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
    478         // for the test fits, use a somewhat smaller radius
    479         psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5);
    480474        psphotFitSersicIndex (model, readout, source, fitOptions, maskVal, markVal);
    481     }
    482 
    483     if (!psphotSetRadiusModel (model, readout, source, markVal, true)) {
    484         psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 1.0);
    485475    }
    486476
    … …  
    494484    pmSourceFitModel (source, model, &options, maskVal);
    495485    // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
    496 
    497486    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
     487
    498488    return (model);
    499489}
     490
     491# define TIMING 0
    500492
    501493pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
    … …  
    517509        return NULL;
    518510    }
     511
     512# define TEST_X -540.0
     513# define TEST_Y 540.0
     514   
     515    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
     516        psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 5);
     517    }
     518
     519    float t1, t2, t4, t5;
     520    if (TIMING) { psTimerStart ("psphotFitPCM"); }
    519521
    520522    pmPCMdata *pcm = pmPCMinit (source, &options, model, maskVal, psfSize);
    … …  
    524526        return model;
    525527    }
    526 
    527     // use the source moments, etc to guess basic model parameters
    528     if (!pmSourceModelGuessPCM (pcm, source, maskVal, markVal)) {
    529         psFree (pcm);
    530         model->flags |= PM_MODEL_STATUS_BADARGS;
    531         return model;
    532     }
    533 
    534     // for sersic models, use a grid search to choose an index, then float the params there
     528    if (TIMING) { t1 = psTimerMark ("psphotFitPCM"); }
     529
     530    // get the guess for sersic models
    535531    if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
    536         // for the test fits, use a somewhat smaller radius
    537         psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5);
    538 
    539         if (!psphotFitSersicIndexPCM (pcm, readout, source, fitOptions, maskVal, markVal, psfSize)) {
     532        // use the source moments, etc to guess basic model parameters
     533        if (!psphotSersicModelClassGuessPCM (pcm, source)) {
    540534            psFree (pcm);
    541535            model->flags |= PM_MODEL_STATUS_BADARGS;
    542536            return model;
    543537        }
    544     }
    545 
    546     if (!psphotSetRadiusModel (model, readout, source, markVal, true)) {
    547         psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 1.0);
    548     }
     538    } else {
     539        // use the source moments, etc to guess basic model parameters
     540        if (!pmSourceModelGuessPCM (pcm, source, maskVal, markVal)) {
     541            psFree (pcm);
     542            model->flags |= PM_MODEL_STATUS_BADARGS;
     543            return model;
     544        }
     545    }
     546
     547    if (TIMING) { t2 = psTimerMark ("psphotFitPCM"); }
    549548
    550549    if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
    … …  
    555554    // update the pcm elements if we have changed the circumstance (options.mode or source->pixels)
    556555    pmPCMupdate(pcm, source, &options, model);
     556    if (TIMING) { t4 = psTimerMark ("psphotFitPCM"); }
    557557
    558558    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
    559559    pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
    560     // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
     560    if (TIMING) { t5 = psTimerMark ("psphotFitPCM"); }
     561
     562    if (TIMING) {
     563        int nPixBig = source->pixels->numCols * source->pixels->numRows;
     564        fprintf (stderr, "psphotFitPCM : nIter: %2d, radius: %6.1f, npix: %5d of %5d, t1: %6.4f, t2: %6.4f, t4: %6.4f, t5: %6.4f\n", model->nIter, model->fitRadius, model->nPix, nPixBig, t1, t2, t4, t5);
     565    }
     566
     567    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
     568        psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 0);
     569    }
    561570
    562571    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
    … …  
    566575}
    567576
     577# undef TEST_X
     578# undef TEST_Y
     579
    568580// note that these should be 1/2n of the standard sersic index
    569 float indexGuess[] = {0.5, 0.33, 0.25, 0.167, 0.125, 0.083};
    570 # define N_INDEX_GUESS 6
     581// float indexGuess[] = {0.8, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0};
     582// float indexGuess[] = {0.5, 0.33, 0.25, 0.167, 0.125, 0.083};
     583float indexGuess[] = {1.0, 2.0, 3.0, 4.0};
     584# define N_INDEX_GUESS 4
    571585
    572586// A sersic model is very sensitive to the index.  attempt to find the index first by grid search in just the index
    … …  
    586600    float chiSquare[N_INDEX_GUESS];
    587601
     602# define TEST_X -540.0
     603# define TEST_Y 540.0
     604   
     605    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
     606        psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 5);
     607    }
     608
    588609    for (int i = 0; i < N_INDEX_GUESS; i++) {
    589         model->params->data.F32[PM_PAR_7] = indexGuess[i];
     610        model->params->data.F32[PM_PAR_7] = 0.5/indexGuess[i];
    590611
    591612        if (!model->modelGuess(model, source)) {
    … …  
    594615        }
    595616
    596         // each time we change the model guess, we need to adjust the radius
    597         // XXX this did not work : we do not need such a large radius -- just uses moments-based radius
    598         if (false && !psphotSetRadiusModel (model, readout, source, markVal, false)) {
    599             psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5);
    600         }
    601        
    602617        pmSourceFitModel (source, model, &options, maskVal);
    603         // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
     618        // fprintf (stderr, "index: %f, chisq: %f, nIter: %d, radius: %f, npix: %d\n", indexGuess[i], model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
    604619
    605620        chiSquare[i] = model->chisqNorm;
    … …  
    616631    assert (iMin >= 0);
    617632
     633    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
     634        psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 0);
     635    }
     636
    618637    model->flags = PM_MODEL_STATUS_NONE; // do not attempt to handle failures here, let the next iteration deal with it
    619     model->params->data.F32[PM_PAR_7] = indexGuess[iMin];
     638    model->params->data.F32[PM_PAR_7] = 0.5/indexGuess[iMin];
    620639    model->modelGuess(model, source);
    621 
    622     // each time we change the model guess, we need to adjust the radius
    623     // if (!psphotSetRadiusModel (model, readout, source, markVal, true)) {
    624     //  psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal);
    625     // }
    626640
    627641    return true;
    … …  
    648662    float xMin = NAN;
    649663    float chiSquare[N_INDEX_GUESS];
     664
     665    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
     666        psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 5);
     667    }
    650668
    651669    for (int i = 0; i < N_INDEX_GUESS; i++) {
    … …  
    657675        }
    658676
    659         // each time we change the model guess, we need to adjust the radius
    660         // XXX this did not work : we do not need such a large radius -- just uses moments-based radius
    661         if (false && !psphotSetRadiusModel (model, readout, source, markVal, false)) {
    662             psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5);
    663         }
    664        
    665677# if (0)
     678        // this block is to test the relative speed of straight and PCM fits
    666679        pmSourceFitModel (source, model, &options, maskVal);
    667680# else
    … …  
    669682        pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
    670683# endif
     684        fprintf (stderr, "index: %f, chisq: %f, nIter: %d, radius: %f, npix: %d\n", indexGuess[i], model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
    671685        // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
    672686
    … …  
    684698    assert (iMin >= 0);
    685699   
     700    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
     701        psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 0);
     702    }
     703
    686704    model->flags = PM_MODEL_STATUS_NONE; // do not attempt to handle failures here, let the next iteration deal with it
    687705    model->params->data.F32[PM_PAR_7] = indexGuess[iMin];
Note: See TracChangeset for help on using the changeset viewer.