IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 22, 2014, 8:48:02 AM (12 years ago)
Author:
eugene
Message:

create a new function psphotPSFstatsSources to use real sources to measure the psf stats

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotChoosePSF.c

    r36863 r37767  
    11# include "psphotInternal.h"
     2bool psphotPSFstatsSources (pmReadout *readout, psArray *sources, pmPSF *psf);
    23
    34// generate a PSF model for inputs without PSF models already loaded
    … …  
    388389    psFree (models);
    389390
    390     if (!psphotPSFstats (readout, psf)) {
     391    // if (!psphotPSFstats (readout, psf)) {
     392
     393    if (!psphotPSFstatsSources (readout, sources, psf)) {
    391394        psError(PSPHOT_ERR_PSF, false, "cannot measure PSF shape terms");
    392395        psFree(options);
    … …  
    466469
    467470            if (modelPSF->params->n > 7) {
    468               psVectorAppend (psfExtra1, modelPSF->params->data.F32[7]);
     471                psVectorAppend (psfExtra1, modelPSF->params->data.F32[7]);
    469472            }
    470473            if (modelPSF->params->n > 8) {
    471               psVectorAppend (psfExtra2, modelPSF->params->data.F32[8]);
     474                psVectorAppend (psfExtra2, modelPSF->params->data.F32[8]);
    472475            }
    473476            psFree (modelPSF);
    474477
    475             // float fwhmtest = pmPSFtoFWHM(psf, xc, yc);
    476             // fprintf (stderr, "fwhm: %f, %f : %f\n", FWHM_MAJOR, FWHM_MINOR, fwhmtest);
     478            float fwhmtest = pmPSFtoFWHM(psf, xc, yc);
     479            fprintf (stderr, "fwhm: %f, %f : %f\n", FWHM_MAJOR, FWHM_MINOR, fwhmtest);
    477480        }
    478481    }
    … …  
    522525
    523526    if (psfExtra1->n) {
    524       if (!psVectorStats (stats, psfExtra1, NULL, NULL, 0)) {
    525         psError(PS_ERR_UNKNOWN, false, "failure to measure stats for PSF EXTRA 1");
     527        if (!psVectorStats (stats, psfExtra1, NULL, NULL, 0)) {
     528            psError(PS_ERR_UNKNOWN, false, "failure to measure stats for PSF EXTRA 1");
     529            goto escape;
     530        }
     531        psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PSF_EXT1", PS_META_REPLACE, "PSF extra param 1", stats->sampleMean);
     532        psLogMsg ("psphot", PS_LOG_DETAIL, "PSF extra parameter 1: %f +/- %f", stats->sampleMean, stats->sampleStdev);
     533    }
     534
     535    if (psfExtra2->n) {
     536        if (!psVectorStats (stats, psfExtra2, NULL, NULL, 0)) {
     537            psError(PS_ERR_UNKNOWN, false, "failure to measure stats for PSF EXTRA 2");
     538            goto escape;
     539        }
     540        psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PSF_EXT2", PS_META_REPLACE, "PSF extra param 2", stats->sampleMean);
     541        psLogMsg ("psphot", PS_LOG_DETAIL, "PSF extra parameter 2: %f +/- %f", stats->sampleMean, stats->sampleStdev);
     542    }
     543
     544    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "ANGLE",    PS_META_REPLACE, "PSF angle",           axes.theta);
     545    psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "NPSFSTAR", PS_META_REPLACE, "Number of stars used to make PSF", psf->nPSFstars);
     546
     547    // psLogMsg ("psphot", PS_LOG_DETAIL, "PSF angle: %f, nstars: %d", axes.theta, psf->nPSFstars);
     548
     549    char *psfModelName = pmModelClassGetName(psf->type);
     550    psMetadataAddStr(readout->analysis,  PS_LIST_TAIL, "PSFMODEL", PS_META_REPLACE, "PSF Model Name", psfModelName);
     551    psMetadataAddBool(readout->analysis, PS_LIST_TAIL, "PSF_OK",   PS_META_REPLACE, "Valid PSF Model?", true);
     552
     553    int nParams = pmModelClassParameterCount(psf->type);
     554    psMetadataAddS32(readout->analysis, PS_LIST_TAIL, "PSF_NPAR",   PS_META_REPLACE, "Number of PSF parameters", nParams);
     555
     556    psFree (fwhmMajor);
     557    psFree (fwhmMinor);
     558    psFree (psfExtra1);
     559    psFree (psfExtra2);
     560    psFree (stats);
     561    return true;
     562
     563escape:
     564    psFree (fwhmMajor);
     565    psFree (fwhmMinor);
     566    psFree (psfExtra1);
     567    psFree (psfExtra2);
     568    psFree (stats);
     569    return false;
     570}
     571
     572// measure average parameters of the PSF model
     573bool psphotPSFstatsSources (pmReadout *readout, psArray *sources, pmPSF *psf) {
     574
     575    psEllipseShape shape;
     576    psEllipseAxes axes;
     577
     578    PS_ASSERT_PTR_NON_NULL(readout, false);
     579    PS_ASSERT_PTR_NON_NULL(psf, false);
     580
     581    psImage *image = readout->image;
     582    PS_ASSERT_PTR_NON_NULL(image, false);
     583
     584    // XXX a minor hack: measure the values for a grid of points across the image, determine mean, UQ, LQ, etc:
     585    psVector *fwhmMajor = psVectorAllocEmpty (100, PS_DATA_F32);
     586    psVector *fwhmMinor = psVectorAllocEmpty (100, PS_DATA_F32);
     587    psVector *psfExtra1 = psVectorAllocEmpty (100, PS_DATA_F32);
     588    psVector *psfExtra2 = psVectorAllocEmpty (100, PS_DATA_F32);
     589
     590    for (int i = 0; i < sources->n; i++) {
     591
     592        pmSource *source = sources->data[i];
     593
     594        if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
     595
     596        float xc = source->peak->xf;
     597        float yc = source->peak->yf;
     598
     599        // create modelPSF from this model
     600        pmModel *modelPSF = pmModelFromPSFforXY (psf, xc, yc, 1.0);
     601        if (!modelPSF) {
     602            fprintf (stderr, "?");
     603            continue;
     604        }
     605
     606        // get the model full-width at half-max
     607        float FWHM_MAJOR = 2*modelPSF->class->modelRadius (modelPSF->params, 0.5);
     608
     609        // XXX make sure this is consistent with the re-definition of PM_PAR_SXX
     610        shape.sx  = modelPSF->params->data.F32[PM_PAR_SXX];
     611        shape.sy  = modelPSF->params->data.F32[PM_PAR_SYY];
     612        shape.sxy = modelPSF->params->data.F32[PM_PAR_SXY];
     613        axes = psEllipseShapeToAxes (shape, 20.0);
     614
     615        float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
     616        if (!isfinite(FWHM_MAJOR) || !isfinite(FWHM_MINOR)) {
     617            fprintf (stderr, "!");
     618            psFree (modelPSF);
     619            continue;
     620        }
     621        psVectorAppend (fwhmMajor, FWHM_MAJOR);
     622        psVectorAppend (fwhmMinor, FWHM_MINOR);
     623
     624        if (modelPSF->params->n > 7) {
     625            psVectorAppend (psfExtra1, modelPSF->params->data.F32[7]);
     626        }
     627        if (modelPSF->params->n > 8) {
     628            psVectorAppend (psfExtra2, modelPSF->params->data.F32[8]);
     629        }
     630        psFree (modelPSF);
     631
     632        float fwhmtest = pmPSFtoFWHM(psf, xc, yc);
     633        fprintf (stderr, "fwhm: %f, %f : %f\n", FWHM_MAJOR, FWHM_MINOR, fwhmtest);
     634    }
     635
     636    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_SAMPLE_QUARTILE);
     637    if (!psVectorStats (stats, fwhmMajor, NULL, NULL, 0)) {
     638        psError(PS_ERR_UNKNOWN, false, "failure to measure stats for FWHM MAJOR");
    526639        goto escape;
    527       }
    528       psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PSF_EXT1", PS_META_REPLACE, "PSF extra param 1", stats->sampleMean);
    529       psLogMsg ("psphot", PS_LOG_DETAIL, "PSF extra parameter 1: %f +/- %f", stats->sampleMean, stats->sampleStdev);
     640    }
     641
     642    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FWHM_MAJ",   PS_META_REPLACE, "PSF FWHM Major axis (mean)", stats->sampleMean);
     643    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MJ_SG",   PS_META_REPLACE, "PSF FWHM Major axis (sigma)", stats->sampleStdev);
     644    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MJ_LQ",   PS_META_REPLACE, "PSF FWHM Major axis (lower quartile)", stats->sampleLQ);
     645    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MJ_UQ",   PS_META_REPLACE, "PSF FWHM Major axis (upper quartile)", stats->sampleUQ);
     646
     647    float fwhmMaj = stats->sampleMean;  // FWHM on major axis
     648
     649    if (!psVectorStats (stats, fwhmMinor, NULL, NULL, 0)) {
     650        psError(PS_ERR_UNKNOWN, false, "failure to measure stats for FWHM MINOR");
     651        goto escape;
     652    }
     653    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FWHM_MIN",   PS_META_REPLACE, "PSF FWHM Minor axis (mean)", stats->sampleMean);
     654    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MN_SG",   PS_META_REPLACE, "PSF FWHM Minor axis (sigma)", stats->sampleStdev);
     655    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MN_LQ",   PS_META_REPLACE, "PSF FWHM Minor axis (lower quartile)", stats->sampleLQ);
     656    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MN_UQ",   PS_META_REPLACE, "PSF FWHM Minor axis (upper quartile)", stats->sampleUQ);
     657
     658    float fwhmMin = stats->sampleMean;  // FWHM on minor axis
     659
     660    if (readout->parent) {
     661
     662        // we now have 2 possible measurements of the seeing : the PSF model based version and
     663        // the Moments based value we need to define a definitive "CHIP.SEEING" value.  This
     664        // value is used by the PSF-matching programs (ppSub and ppStack).  The moments-based
     665        // value is probably more representative of the value needed by those tools (it is also
     666        // more stable?)
     667        bool status = false;
     668        float fwhmMajorMoments = psMetadataLookupF32(&status, readout->analysis, "IQ_FW1");
     669        float fwhmMinorMoments = psMetadataLookupF32(&status, readout->analysis, "IQ_FW1");
     670
     671        pmChip *chip = readout->parent->parent; // Parent chip
     672        psAssert(chip, "Cell should be attached to a chip.");
     673        psMetadataItem *item = psMetadataLookup(chip->concepts, "CHIP.SEEING"); // Item with chip
     674        item->data.F32 = 0.5 * (fwhmMajorMoments + fwhmMinorMoments);
     675
     676        psLogMsg ("psphot", PS_LOG_DETAIL, "fwhm (psf): %f,%f (moments): %f,%f", fwhmMaj, fwhmMin, fwhmMajorMoments, fwhmMinorMoments);
     677    }
     678
     679    if (psfExtra1->n) {
     680        if (!psVectorStats (stats, psfExtra1, NULL, NULL, 0)) {
     681            psError(PS_ERR_UNKNOWN, false, "failure to measure stats for PSF EXTRA 1");
     682            goto escape;
     683        }
     684        psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PSF_EXT1", PS_META_REPLACE, "PSF extra param 1", stats->sampleMean);
     685        psLogMsg ("psphot", PS_LOG_DETAIL, "PSF extra parameter 1: %f +/- %f", stats->sampleMean, stats->sampleStdev);
    530686    }
    531687
    532688    if (psfExtra2->n) {
    533       if (!psVectorStats (stats, psfExtra2, NULL, NULL, 0)) {
    534         psError(PS_ERR_UNKNOWN, false, "failure to measure stats for PSF EXTRA 2");
    535         goto escape;
    536       }
    537       psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PSF_EXT2", PS_META_REPLACE, "PSF extra param 2", stats->sampleMean);
    538       psLogMsg ("psphot", PS_LOG_DETAIL, "PSF extra parameter 2: %f +/- %f", stats->sampleMean, stats->sampleStdev);
     689        if (!psVectorStats (stats, psfExtra2, NULL, NULL, 0)) {
     690            psError(PS_ERR_UNKNOWN, false, "failure to measure stats for PSF EXTRA 2");
     691            goto escape;
     692        }
     693        psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PSF_EXT2", PS_META_REPLACE, "PSF extra param 2", stats->sampleMean);
     694        psLogMsg ("psphot", PS_LOG_DETAIL, "PSF extra parameter 2: %f +/- %f", stats->sampleMean, stats->sampleStdev);
    539695    }
    540696
Note: See TracChangeset for help on using the changeset viewer.