IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Sep 20, 2009, 10:08:28 AM (17 years ago)
Author:
eugene
Message:

cleanup and attempt to improve psf model 2D variation analysis

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/20090715/psModules/src/objects/pmPSFtry.c

    r25452 r25455  
    3737#include "pmSourceVisual.h"
    3838
    39 bool printTrendMap (pmTrend2D *trend) {
    40 
    41     if (!trend->map) return false;
    42     if (!trend->map->map) return false;
    43 
    44     for (int j = 0; j < trend->map->map->numRows; j++) {
    45         for (int i = 0; i < trend->map->map->numCols; i++) {
    46             fprintf (stderr, "%5.2f  ", trend->map->map->data.F32[j][i]);
    47         }
    48         fprintf (stderr, "\t\t\t");
    49         for (int i = 0; i < trend->map->map->numCols; i++) {
    50             fprintf (stderr, "%5.2f  ", trend->map->error->data.F32[j][i]);
    51         }
    52         fprintf (stderr, "\n");
    53     }
    54     return true;
    55 }
    56 
    57 bool psImageMapCleanup (psImageMap *map) {
    58 
    59     if ((map->map->numRows == 1) && (map->map->numCols == 1)) return true;
    60 
    61     // find the weighted average of all pixels
    62     float Sum = 0.0;
    63     float Wt = 0.0;
    64     for (int j = 0; j < map->map->numRows; j++) {
    65         for (int i = 0; i < map->map->numCols; i++) {
    66             if (!isfinite(map->map->data.F32[j][i])) continue;
    67             Sum += map->map->data.F32[j][i] * map->error->data.F32[j][i];
    68             Wt += map->error->data.F32[j][i];
    69         }
    70     }
    71 
    72     float Mean = Sum / Wt;
    73 
    74     // do any of the pixels in the map need to be repaired?
    75     // XXX for now, we are just replacing bad pixels with the Mean
    76     for (int j = 0; j < map->map->numRows; j++) {
    77         for (int i = 0; i < map->map->numCols; i++) {
    78             if (isfinite(map->map->data.F32[j][i]) &&
    79                 (map->error->data.F32[j][i] > 0.0)) continue;
    80             map->map->data.F32[j][i] = Mean;
    81         }
    82     }
    83     return true;
    84 }
    85 
    8639// ********  pmPSFtry functions  **************************************************
    8740// * pmPSFtry holds a single pmPSF model test, with the input sources, the freely
    … …  
    170123
    171124    // stage 1:  fit an EXT model to all candidates PSF sources -- this is independent of the modeled 2D variations in the PSF
    172     psTimerStart ("psf.fit");
    173     for (int i = 0; i < psfTry->sources->n; i++) {
    174 
    175         pmSource *source = psfTry->sources->data[i];
    176         if (!source->moments) {
    177             psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
    178             continue;
    179         }
    180         if (!source->moments->nPixels) {
    181             psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
    182             continue;
    183         }
    184 
    185         source->modelEXT = pmSourceModelGuess (source, psfTry->psf->type);
    186         if (source->modelEXT == NULL) {
    187             psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
    188             psTrace ("psModules.objects", 4, "masking %d (%d,%d) : failed to generate model guess\n", i, source->peak->x, source->peak->y);
    189             continue;
    190         }
    191 
    192         // set object mask to define valid pixels
    193         // XXX 0.5 PIX: is the circle symmetric about the peak coordinate (given 0.5,0.5 center)?
    194         psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", markVal);
    195 
    196         // fit model as EXT, not PSF
    197         status = pmSourceFitModel (source, source->modelEXT, PM_SOURCE_FIT_EXT, maskVal);
    198 
    199         // clear object mask to define valid pixels
    200         psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal));
    201 
    202         // exclude the poor fits
    203         if (!status) {
    204             psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
    205             psTrace ("psModules.objects", 4, "masking %d (%d,%d) : status is poor\n", i, source->peak->x, source->peak->y);
    206             continue;
    207         }
    208         Next ++;
    209     }
    210     psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit ext:   %f sec for %d of %ld sources\n", psTimerMark ("psf.fit"), Next, sources->n);
    211     psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (EXT)\n", Next, sources->n);
    212 
    213     if (Next == 0) {
    214         psError(PS_ERR_UNKNOWN, false, "No sources with good extended fits from which to determine PSF.");
     125    if (!pmPSFtryFitEXT(psfTry, options, maskVal, markVal)) {
     126        psError(PS_ERR_UNKNOWN, false, "failed to fit EXT models to sources for psf model");
    215127        psFree(psfTry);
    216128        return NULL;
    217     }
    218 
    219     // stage 2: construct a psf (pmPSF) from this collection of model fits, including the 2D variation
    220     if (!pmPSFFromPSFtry (psfTry)) {
    221         psError(PS_ERR_UNKNOWN, false, "failed to construct a psf model from collection of sources");
    222         psFree(psfTry);
    223         return NULL;
    224     }
    225 
    226     // stage 3: refit with fixed shape parameters
    227     psTimerStart ("psf.fit");
    228     for (int i = 0; i < psfTry->sources->n; i++) {
    229 
    230         pmSource *source = psfTry->sources->data[i];
    231 
    232         // masked for: bad model fit, outlier in parameters
    233         if (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) {
    234             psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : source is masked\n", i, source->peak->x, source->peak->y);
    235             continue;
    236         }
    237 
    238         // set shape for this model based on PSF
    239         source->modelPSF = pmModelFromPSF (source->modelEXT, psfTry->psf);
    240         if (source->modelPSF == NULL) {
    241             psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_MODEL;
    242             abort();
    243             continue;
    244         }
    245         source->modelPSF->radiusFit = options->radius;
    246 
    247         // XXXX use a different radius for the aperture magnitude than for the PSF fit?
    248 
    249         // set object mask to define valid pixels
    250         psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", markVal);
    251 
    252         // fit the PSF model to the source
    253         status = pmSourceFitModel (source, source->modelPSF, PM_SOURCE_FIT_PSF, maskVal);
    254 
    255         // skip poor fits
    256         if (!status) {
    257             psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal));
    258             psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_PSF_FAIL;
    259             psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : failed PSF fit\n", i, source->peak->x, source->peak->y);
    260             continue;
    261         }
    262 
    263         status = pmSourceMagnitudes (source, psfTry->psf, PM_SOURCE_PHOT_INTERP, maskVal);
    264         if (!status || isnan(source->apMag)) {
    265             psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal));
    266             psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_PHOT;
    267             psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : poor photometry\n", i, source->peak->x, source->peak->y);
    268             continue;
    269         }
    270 
    271         // clear object mask to define valid pixels
    272         psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal));
    273 
    274         psfTry->fitMag->data.F32[i] = source->psfMag;
    275         psfTry->metric->data.F32[i] = source->apMag - source->psfMag;
    276         psfTry->metricErr->data.F32[i] = source->errMag;
    277 
    278         psTrace ("psModules.object", 6, "keeping source %d (%d) of %ld\n", i, Npsf, psfTry->sources->n);
    279         Npsf ++;
    280     }
    281     psfTry->psf->nPSFstars = Npsf;
    282 
    283     psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit psf:   %f sec for %d of %ld sources\n", psTimerMark ("psf.fit"), Npsf, sources->n);
    284     psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (PSF)\n", Npsf, sources->n);
    285 
    286     if (Npsf == 0) {
    287         psError(PS_ERR_UNKNOWN, false, "No sources with good PSF fits after model is built.");
    288         psFree(psfTry);
    289         return NULL;
    290     }
     129    }     
     130
     131    for (int i = 0; i < Norder; i++) {
     132        // stage 2: construct a psf (pmPSF) from this collection of model fits, including the 2D variation
     133        if (!pmPSFFromPSFtry (psfTry, Nx, Ny)) {
     134            psError(PS_ERR_UNKNOWN, false, "failed to construct a psf model from collection of sources");
     135            psFree(psfTry);
     136            return NULL;
     137        }
     138
     139        // stage 3: refit with fixed shape parameters, measure pmPSFtryMetric
     140        if (!pmPSFtryFitPSF (psfTry, Nx, Ny)) {
     141            psError(PS_ERR_UNKNOWN, false, "failed to construct a psf model from collection of sources");
     142            psFree(psfTry);
     143            return NULL;
     144        }
     145    }
     146    // XXXXX this is probably not used any more.  Are the chisq of the fits so bad? can we
     147    // fix them by softening the errors on the brightest pixels?
    291148
    292149    // measure the chi-square trend as a function of flux (PAR[PM_PAR_I0])
    … …  
    485342  (aprMag - rflux*skyBias) = fitMag + ApTrend(x,y)
    486343*/
    487 
    488 /*****************************************************************************
    489 pmPSFFromPSFtry (psfTry): build a PSF model from a collection of
    490 source->modelEXT entries.  The PSF ignores the first 4 (independent) model
    491 parameters and constructs a polynomial fit to the remaining as a function of
    492 image coordinate.
    493     input: psfTry with fitted source->modelEXT collection, pre-allocated psf
    494 Note: some of the array entries may be NULL (failed fits); ignore them.
    495  *****************************************************************************/
    496 bool pmPSFFromPSFtry (pmPSFtry *psfTry)
    497 {
    498     PS_ASSERT_PTR_NON_NULL(psfTry, false);
    499     PS_ASSERT_PTR_NON_NULL(psfTry->sources, false);
    500 
    501     pmPSF *psf = psfTry->psf;
    502 
    503     // construct the fit vectors from the collection of objects
    504     psVector *x  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
    505     psVector *y  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
    506     psVector *z  = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32);
    507 
    508     // construct the x,y terms
    509     for (int i = 0; i < psfTry->sources->n; i++) {
    510         pmSource *source = psfTry->sources->data[i];
    511         if (source->modelEXT == NULL)
    512             continue;
    513 
    514         x->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_XPOS];
    515         y->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_YPOS];
    516     }
    517 
    518     if (!pmPSFFitShapeParams (psf, psfTry->sources, x, y, psfTry->mask)) {
    519         psFree(x);
    520         psFree(y);
    521         psFree(z);
    522         return false;
    523     }
    524 
    525     // skip the unfitted parameters (X, Y, Io, Sky) and the shape parameters (SXX, SYY, SXY)
    526     // apply the values of Nx, Ny determined above for E0,E1,E2 to the remaining parameters
    527     for (int i = 0; i < psf->params->n; i++) {
    528         switch (i) {
    529           case PM_PAR_SKY:
    530           case PM_PAR_I0:
    531           case PM_PAR_XPOS:
    532           case PM_PAR_YPOS:
    533           case PM_PAR_SXX:
    534           case PM_PAR_SYY:
    535           case PM_PAR_SXY:
    536             continue;
    537           default:
    538             break;
    539         }
    540 
    541         // select the per-object fitted data for this parameter
    542         for (int j = 0; j < psfTry->sources->n; j++) {
    543             pmSource *source = psfTry->sources->data[j];
    544             if (source->modelEXT == NULL) continue;
    545             z->data.F32[j] = source->modelEXT->params->data.F32[i];
    546         }
    547 
    548         psImageBinning *binning = psImageBinningAlloc();
    549         binning->nXruff = psf->trendNx;
    550         binning->nYruff = psf->trendNy;
    551         binning->nXfine = psf->fieldNx;
    552         binning->nYfine = psf->fieldNy;
    553 
    554         if (psf->psfTrendMode == PM_TREND_MAP) {
    555             psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER);
    556             psImageBinningSetSkipByOffset (binning, psf->fieldXo, psf->fieldYo);
    557         }
    558 
    559         // free existing trend, re-alloc
    560         psFree (psf->params->data[i]);
    561         psf->params->data[i] = pmTrend2DNoImageAlloc (psf->psfTrendMode, binning, psf->psfTrendStats);
    562         psFree (binning);
    563 
    564         // fit the collection of measured parameters to the PSF 2D model
    565         // the mask is carried from previous steps and updated with this operation
    566         // the weight is either the flux error or NULL, depending on 'psf->poissonErrorParams'
    567         if (!pmTrend2DFit (psf->params->data[i], psfTry->mask, 0xff, x, y, z, NULL)) {
    568             psError(PS_ERR_UNKNOWN, false, "failed to build psf model for parameter %d", i);
    569             psFree(x);
    570             psFree(y);
    571             psFree(z);
    572             return false;
    573         }
    574     }
    575 
    576     // test dump of star parameters vs position (compare with fitted values)
    577     if (psTraceGetLevel("psModules.objects") >= 4) {
    578         FILE *f = fopen ("params.dat", "w");
    579 
    580         for (int j = 0; j < psfTry->sources->n; j++) {
    581             pmSource *source = psfTry->sources->data[j];
    582             if (source == NULL) continue;
    583             if (source->modelEXT == NULL) continue;
    584 
    585             pmModel *modelPSF = pmModelFromPSF (source->modelEXT, psf);
    586 
    587             fprintf (f, "%f %f : ", source->modelEXT->params->data.F32[PM_PAR_XPOS], source->modelEXT->params->data.F32[PM_PAR_YPOS]);
    588 
    589             for (int i = 0; i < psf->params->n; i++) {
    590                 if (psf->params->data[i] == NULL) continue;
    591                 fprintf (f, "%f %f : ", source->modelEXT->params->data.F32[i], modelPSF->params->data.F32[i]);
    592             }
    593             fprintf (f, "%f %d\n", source->modelEXT->chisq, source->modelEXT->nIter);
    594 
    595             psFree(modelPSF);
    596         }
    597         fclose (f);
    598     }
    599 
    600     psFree (x);
    601     psFree (y);
    602     psFree (z);
    603     return true;
    604 }
    605344
    606345bool pmPSFFitShapeParams (pmPSF *psf, psArray *sources, psVector *x, psVector *y, psVector *srcMask) {
    … …  
    900639        stdev = psStatsGetValue (trend->stats, stdevOption);
    901640        psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0obs_fit->n);
    902         // printTrendMap (trend);
     641        // pmTrend2DPrintMap (trend);
    903642        psImageMapCleanup (trend->map);
    904         // printTrendMap (trend);
     643        // pmTrend2DPrintMap (trend);
    905644        pmSourceVisualPSFModelResid (trend, x, y, e0obs, mask);
    906645
    … …  
    910649        stdev = psStatsGetValue (trend->stats, stdevOption);
    911650        psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1obs_fit->n);
    912         // printTrendMap (trend);
     651        // pmTrend2DPrintMap (trend);
    913652        psImageMapCleanup (trend->map);
    914         // printTrendMap (trend);
     653        // pmTrend2DPrintMap (trend);
    915654        pmSourceVisualPSFModelResid (trend, x, y, e1obs, mask);
    916655
    … …  
    920659        stdev = psStatsGetValue (trend->stats, stdevOption);
    921660        psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2obs->n);
    922         // printTrendMap (trend);
     661        // pmTrend2DPrintMap (trend);
    923662        psImageMapCleanup (trend->map);
    924         // printTrendMap (trend);
     663        // pmTrend2DPrintMap (trend);
    925664        pmSourceVisualPSFModelResid (trend, x, y, e2obs, mask);
    926665    }
Note: See TracChangeset for help on using the changeset viewer.