IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 14786


Ignore:
Timestamp:
Sep 7, 2007, 10:59:43 AM (19 years ago)
Author:
eugene
Message:

convert to ApTrend, simply selection options

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branch_20070830/psphot/src/psphotApResid.c

    r14732 r14786  
    11# include "psphotInternal.h"
    2 // XXXX this code fails if there are too few sources to measure the aperture residual
    3 // the larger problem is that the rules for accepting more polynomial terms are weak.
    42
    5 static pmPSFApTrendOptions DEFAULT_OPTION = PM_PSF_APTREND_SKYBIAS;
    6 
    7 // measure the aperture residual statistics
     3// measure the aperture residual statistics and 2D variations
    84bool psphotApResid (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psMaskType maskVal, psMaskType mark)
    95{
    … …  
    1612
    1713    PS_ASSERT_PTR_NON_NULL(psf, false);
    18     PS_ASSERT_PTR_NON_NULL(psf->ApTrend, false);
     14    // XXX drop this? PS_ASSERT_PTR_NON_NULL(psf->ApTrend, false);
    1915    PS_ASSERT_PTR_NON_NULL(readout, false);
    2016    PS_ASSERT_PTR_NON_NULL(sources, false);
    … …  
    3026    bool IGNORE_GROWTH = psMetadataLookupBool (&status, recipe, "IGNORE_GROWTH");
    3127    bool INTERPOLATE_AP = psMetadataLookupBool (&status, recipe, "INTERPOLATE_AP");
    32     int NSTAR_APERTURE_CORRECTION_MIN =
    33         psMetadataLookupS32(&status, recipe, "NSTAR_APERTURE_CORRECTION_MIN");
     28
     29    // XXX is this still needed?  the pmTrend2D stuff should be auto-adjusting...
     30    int NSTAR_APERTURE_CORRECTION_MIN = psMetadataLookupS32(&status, recipe, "NSTAR_APERTURE_CORRECTION_MIN");
    3431    if (!status) {
    3532        NSTAR_APERTURE_CORRECTION_MIN = 5;
    … …  
    5451
    5552    psVector *mask    = psVectorAllocEmpty (300, PS_TYPE_U8);
    56     psVector *xPos    = psVectorAllocEmpty (300, PS_TYPE_F64);
    57     psVector *yPos    = psVectorAllocEmpty (300, PS_TYPE_F64);
    58     psVector *apResid = psVectorAllocEmpty (300, PS_TYPE_F64);
    59     psVector *dMag    = psVectorAllocEmpty (300, PS_TYPE_F64);
     53    psVector *xPos    = psVectorAllocEmpty (300, PS_TYPE_F32);
     54    psVector *yPos    = psVectorAllocEmpty (300, PS_TYPE_F32);
     55    psVector *apResid = psVectorAllocEmpty (300, PS_TYPE_F32);
     56    psVector *dMag    = psVectorAllocEmpty (300, PS_TYPE_F32);
    6057    Npsf = 0;
    6158
    … …  
    9794        }
    9895
    99         apResid->data.F64[Npsf] = dap
    100         xPos->data.F64[Npsf]    = model->params->data.F32[PM_PAR_XPOS];
    101         yPos->data.F64[Npsf]    = model->params->data.F32[PM_PAR_YPOS];
     96        apResid->data.F32[Npsf] = dap;
     97        xPos->data.F32[Npsf]    = model->params->data.F32[PM_PAR_XPOS];
     98        yPos->data.F32[Npsf]    = model->params->data.F32[PM_PAR_YPOS];
    10299
    103100        mask->data.U8[Npsf] = 0;
    104101
    105         dMag->data.F64[Npsf] = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
     102        dMag->data.F32[Npsf] = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
    106103
    107104        if (psTraceGetLevel("psphot") >= 5) {
    108105            fprintf (dumpFile, "%f %f  %f %f %f  %x\n",
    109                      xPos->data.F64[Npsf], yPos->data.F64[Npsf],
    110                      source->psfMag, dMag->data.F64[Npsf], apResid->data.F64[Npsf],
     106                     xPos->data.F32[Npsf], yPos->data.F32[Npsf],
     107                     source->psfMag, dMag->data.F32[Npsf], apResid->data.F32[Npsf],
    111108                     mask->data.U8[Npsf]);
    112109        }
    … …  
    135132
    136133    // XXX deprecating the old code which allowed the ApResid to be fitted as a function of flux and r^2/flux
    137     // APTREND options : NONE CONSTANT SKYBIAS XY_LIN XY_QUAD SKY_XY_LIN FULL
    138     // APTREND options are used in the switch block below
    139     pmPSFApTrendOptions ApTrendOption = DEFAULT_OPTION;
    140     char *optionName = psMetadataLookupPtr (&status, recipe, "APTREND");
    141     if (status) ApTrendOption = pmPSFApTrendOptionFromName (optionName);
    142     if (ApTrendOption == PM_PSF_APTREND_ERROR) {
    143         psError(PSPHOT_ERR_APERTURE, true, "invalid aperture residual trend %s", optionName);
    144         return false;
    145     }
    146 
    147134    // XXX is this asymmetric clipping still needed?  this analysis should come after neighbors are subtracted...
    148135    // 3hi/1lo sigma clipping on the rflux vs metric fit
    149     psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
     136    psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN);
    150137    stats->min = 2.0;
    151138    stats->max = 3.0;
    152139
    153 #define P_APTREND_SWITCH_CLEANUP        /* Cleanup memory on error in ApTrendOption switch */ \
    154     psFree(psf->growth); psf->growth = NULL; \
    155     psFree(mask); \
    156     psFree(xPos); \
    157     psFree(yPos); \
    158     psFree(flux); \
    159     psFree(r2rflux); \
    160     psFree(apResid); \
    161     psFree(dMag); \
    162     psFree(stats)
     140    // XXX set the pmTrend2D mode, nXtrend, nYtrend based on the user inputs
     141    // XXX psf->ApTrend = pmTrend2DAlloc (mode, readout->image->numCols, readout->image->numRows, nXtrend, nYtrend, stats);
    163142
    164     // no correction
    165     switch (ApTrendOption) {
    166       case PM_PSF_APTREND_NONE:
    167         // remove ApTrend fit from pmPSFtry
    168         psf->ApTrend->coeff[0][0][0][0] = 0;
    169         break;
    170       case PM_PSF_APTREND_CONSTANT:
    171         stats->clipIter = 2;
    172         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    173         if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    174             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
    175             P_APTREND_SWITCH_CLEANUP;
    176             return false;
    177         }
    178         // apply the fit
    179         stats->clipIter = 3;
    180         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    181         if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    182             psError(PSPHOT_ERR_PHOTOM, false, "Fitting aperture correction");
    183             P_APTREND_SWITCH_CLEANUP;
    184             return false;
    185         }
    186         break;
    187       case PM_PSF_APTREND_SKYBIAS:
    188         // first clip out objects which are too far from the median
    189         stats->clipIter = 2;
    190         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    191         if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    192             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
    193             P_APTREND_SWITCH_CLEANUP;
    194             return false;
    195         }
    196         // apply the fit
    197         stats->clipIter = 3;
    198         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
    199         if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    200             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias");
    201             P_APTREND_SWITCH_CLEANUP;
    202             return false;
    203         }
    204         break;
    205       case PM_PSF_APTREND_SKYSAT:
    206         // first clip out objects which are too far from the median
    207         stats->clipIter = 2;
    208         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    209         if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    210             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
    211             P_APTREND_SWITCH_CLEANUP;
    212             return false;
    213         }
    214         // apply the fit
    215         stats->clipIter = 2;
    216         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
    217         if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    218             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias");
    219             P_APTREND_SWITCH_CLEANUP;
    220             return false;
    221         }
    222         // apply the fit
    223         stats->clipIter = 3;
    224         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYSAT);
    225         if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    226             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting skysat");
    227             P_APTREND_SWITCH_CLEANUP;
    228             return false;
    229         }
    230         break;
    231       case PM_PSF_APTREND_XY_LIN:
    232         // first clip out objects which are too far from the median
    233         stats->clipIter = 2;
    234         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    235         if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    236             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
    237             P_APTREND_SWITCH_CLEANUP;
    238             return false;
    239         }
    240         // apply the fit
    241         stats->clipIter = 3;
    242         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_XY_LIN);
    243         if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    244             psError(PSPHOT_ERR_PHOTOM, false, "fitting, XY_LIN");
    245             P_APTREND_SWITCH_CLEANUP;
    246             return false;
    247         }
    248         break;
    249       case PM_PSF_APTREND_XY_QUAD:
    250         // first clip out objects which are too far from the median
    251         stats->clipIter = 2;
    252         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    253         if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    254             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
    255             P_APTREND_SWITCH_CLEANUP;
    256             return false;
    257         }
    258         // apply the fit
    259         stats->clipIter = 3;
    260         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_XY_QUAD);
    261         if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    262             psError(PSPHOT_ERR_PHOTOM, false, "Fitting XY_QUAD");
    263             P_APTREND_SWITCH_CLEANUP;
    264             return false;
    265         }
    266         break;
    267       case PM_PSF_APTREND_SKY_XY_LIN:
    268         // first clip out objects which are too far from the median
    269         stats->clipIter = 2;
    270         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    271         if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    272             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
    273             P_APTREND_SWITCH_CLEANUP;
    274             return false;
    275         }
    276         // apply the fit
    277         stats->clipIter = 3;
    278         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKY_XY_LIN);
    279         if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    280             psError(PSPHOT_ERR_PHOTOM, false, "Fitting sky xy_lin");
    281             P_APTREND_SWITCH_CLEANUP;
    282             return false;
    283         }
    284         break;
    285       case PM_PSF_APTREND_SKYSAT_XY_LIN:
    286         // first clip out objects which are too far from the median
    287         stats->clipIter = 2;
    288         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    289         if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    290             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
    291             P_APTREND_SWITCH_CLEANUP;
    292             return false;
    293         }
    294         // apply the fit
    295         stats->clipIter = 3;
    296         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
    297         if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    298             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias");
    299             P_APTREND_SWITCH_CLEANUP;
    300             return false;
    301         }
    302         // apply the fit
    303         stats->clipIter = 3;
    304         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYSAT_XY_LIN);
    305         if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    306             psError(PSPHOT_ERR_PHOTOM, false, "Fitting skyset xy_lin");
    307             P_APTREND_SWITCH_CLEANUP;
    308             return false;
    309         }
    310         break;
    311       case PM_PSF_APTREND_ALL:
    312         // first clip out objects which are too far from the median
    313         stats->clipIter = 2;
    314         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    315         if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    316             psError(PSPHOT_ERR_PHOTOM, false, "Failed to measure apTrend");
    317             P_APTREND_SWITCH_CLEANUP;
    318             return false;
    319         }
    320         // fit just SkyBias and clip out objects which are too far from the median
    321         stats->clipIter = 2;
    322         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
    323         if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    324             psError(PSPHOT_ERR_PHOTOM, false, "fitting skyBias");
    325             P_APTREND_SWITCH_CLEANUP;
    326             return false;
    327         }
    328         // finally, fit x, y, SkyBias and clip out objects which are too far from the median
    329         stats->clipIter = 3;
    330         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_ALL);
    331         if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    332             psError(PSPHOT_ERR_PHOTOM, false, "fitting all");
    333             P_APTREND_SWITCH_CLEANUP;
    334             return false;
    335         }
    336         break;
    337       default:
    338         psError(PSPHOT_ERR_PHOTOM, true, "Unknown APTREND value: %s", optionName);
    339         return false;
    340     }
    341 #undef P_APTREND_SWITCH_CLEANUP
     143    psf->ApTrend = pmTrend2DAlloc (PM_TREND_MAP, readout->image, 5, 5, stats);
     144    pmTrend2DFit (psf->ApTrend, xPos, yPos, apResid, dMag);
     145   
     146    psphotSaveImage (NULL, psf->ApTrend->map->map, "apresid.map.fits");
    342147
    343148    // construct the fitted values and the residuals
    344     psVector *fit   = psPolynomial4DEvalVector (psf->ApTrend, xPos, yPos, r2rflux, flux);
     149    psVector *fit   = pmTrend2DEvalVector (psf->ApTrend, xPos, yPos);
    345150    psVector *resid = (psVector *) psBinaryOp (NULL, (void *) apResid, "-", (void *) fit);
    346151
    … …  
    349154    psStats *residStats = psStatsAlloc (PS_STAT_SAMPLE_STDEV);
    350155    for (int i = 0; i < dMag->n; i++) {
    351         if (dMag->data.F64[i] > 0.01) {
     156        if (dMag->data.F32[i] > 0.01) {
    352157            mask->data.U8[i] |= 0x02;
    353158        }
    … …  
    357162
    358163    // apply ApTrend results
    359     psf->skyBias  = psf->ApTrend->coeff[0][0][1][0];
    360     psf->skySat   = psf->ApTrend->coeff[0][0][0][1];
    361     psf->ApResid  = psf->ApTrend->coeff[0][0][0][0];
     164    psf->ApResid  = pmTrend2DEval (psf->ApTrend, 0.0, 0.0); // XXX use chip center?
    362165    psf->dApResid = residStats->sampleStdev;
    363     psf->ApTrend->coeff[0][0][1][0] = 0;
    364     psf->ApTrend->coeff[0][0][0][1] = 0;
    365166    psf->nApResid = Npsf;
    366167
    367168    // save results for later output
    368     psMetadataAdd (recipe, PS_LIST_TAIL, "SKYBIAS",  PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   psf->skyBias);
    369     psMetadataAdd (recipe, PS_LIST_TAIL, "SKYSAT",   PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   psf->skySat);
     169    psMetadataAdd (recipe, PS_LIST_TAIL, "SKYBIAS",  PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   0.0);
     170    psMetadataAdd (recipe, PS_LIST_TAIL, "SKYSAT",   PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   0.0);
    370171    psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   psf->ApResid);
    371172    psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->dApResid);
    … …  
    373174    psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "ap residual scatter", psf->nApResid);
    374175
    375     psLogMsg ("psphot.apresid", PS_LOG_INFO, "measure full-frame aperture residuals for %d of %d objects: %f sec\n",
    376               Nkeep, Npsf, psTimerMark ("psphot"));
    377     psLogMsg ("psphot.apresid", PS_LOG_DETAIL, "aperture residual: %f +/- %f : %f bias, %f skysat\n",
    378               psf->ApResid, psf->dApResid, psf->skyBias, psf->skySat);
    379     psLogMsg ("psphot.apresid", PS_LOG_MINUTIA, "apresid trends: %f %f %f %f %f\n",
    380               1e3*psf->ApTrend->coeff[1][0][0][0],
    381               1e6*psf->ApTrend->coeff[2][0][0][0],
    382               1e6*psf->ApTrend->coeff[1][1][0][0],
    383               1e3*psf->ApTrend->coeff[0][1][0][0],
    384               1e6*psf->ApTrend->coeff[0][2][0][0]);
     176    psLogMsg ("psphot.apresid", PS_LOG_INFO, "measure full-frame aperture residuals for %d of %d objects: %f sec\n", Nkeep, Npsf, psTimerMark ("psphot"));
     177    psLogMsg ("psphot.apresid", PS_LOG_DETAIL, "aperture residual: %f +/- %f\n", psf->ApResid, psf->dApResid);
    385178
    386179    psFree (mask);
    387180    psFree (xPos);
    388181    psFree (yPos);
    389     psFree (flux);
    390     psFree (r2rflux);
    391182    psFree (apResid);
    392183    psFree (dMag);
Note: See TracChangeset for help on using the changeset viewer.