IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 14944


Ignore:
Timestamp:
Sep 20, 2007, 2:14:47 PM (19 years ago)
Author:
eugene
Message:

substantial changes: converting to pmTrend2d; adding functions to measure the systematic error floor; iteratively choose the best scale to minimize the systematic floor

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotApResid.c

    r13900 r14944  
    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.
    4 
    5 static pmPSFApTrendOptions DEFAULT_OPTION = PM_PSF_APTREND_SKYBIAS;
    6 
    7 // measure the aperture residual statistics
     2
     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 *flux    = psVectorAllocEmpty (300, PS_TYPE_F64);
    59     psVector *r2rflux = psVectorAllocEmpty (300, PS_TYPE_F64);
    60     psVector *apResid = psVectorAllocEmpty (300, PS_TYPE_F64);
    61     psVector *dMag    = psVectorAllocEmpty (300, PS_TYPE_F64);
     53    psVector *mag     = psVectorAllocEmpty (300, PS_TYPE_F32);
     54    psVector *xPos    = psVectorAllocEmpty (300, PS_TYPE_F32);
     55    psVector *yPos    = psVectorAllocEmpty (300, PS_TYPE_F32);
     56    psVector *apResid = psVectorAllocEmpty (300, PS_TYPE_F32);
     57    psVector *dMag    = psVectorAllocEmpty (300, PS_TYPE_F32);
    6258    Npsf = 0;
    63 
    64     FILE *dumpFile = NULL;
    65     if (psTraceGetLevel("psphot") >= 5) {
    66         dumpFile = fopen ("apresid.dat", "w");
    67     }
    6859
    6960    // select all good PM_SOURCE_TYPE_STAR entries
    … …  
    9081        }
    9182
    92         // sanity clipping : if the model is so discrepant, but your expectation in the recipe
    93         apResid->data.F64[Npsf] = source->apMag - source->psfMag;
    94 
    95         if ((MAX_AP_OFFSET > 0) && (fabs(apResid->data.F64[Npsf]) > MAX_AP_OFFSET)) {
     83        // aperture residual for this source
     84        float dap = source->apMag - source->psfMag;
     85
     86        // sanity clipping : if the model is very discrepant, put your expectation in the recipe
     87        if ((MAX_AP_OFFSET > 0) && (fabs(dap) > MAX_AP_OFFSET)) {
    9688            Nfail ++;
    9789            continue;
    9890        }
    9991
    100         xPos->data.F64[Npsf]    = model->params->data.F32[PM_PAR_XPOS];
    101         yPos->data.F64[Npsf]    = model->params->data.F32[PM_PAR_YPOS];
    102 
    103         flux->data.F64[Npsf]    = pow(10.0, -0.4*source->psfMag);
    104         r2rflux->data.F64[Npsf] = PS_SQR(model->radiusFit) / flux->data.F64[Npsf];
     92        mag->data.F32[Npsf]     = source->psfMag;
     93        apResid->data.F32[Npsf] = dap;
     94        xPos->data.F32[Npsf]    = model->params->data.F32[PM_PAR_XPOS];
     95        yPos->data.F32[Npsf]    = model->params->data.F32[PM_PAR_YPOS];
    10596
    10697        mask->data.U8[Npsf] = 0;
    10798
    108         dMag->data.F64[Npsf] = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
    109 
    110         if (psTraceGetLevel("psphot") >= 5) {
    111             fprintf (dumpFile, "%f %f  %f %f %f  %x\n",
    112                      xPos->data.F64[Npsf], yPos->data.F64[Npsf],
    113                      source->psfMag, dMag->data.F64[Npsf], apResid->data.F64[Npsf],
    114                      mask->data.U8[Npsf]);
    115         }
     99        dMag->data.F32[Npsf] = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
    116100
    117101        psVectorExtend (mask,    100, 1);
    118102        psVectorExtend (xPos,    100, 1);
    119103        psVectorExtend (yPos,    100, 1);
    120         psVectorExtend (flux,    100, 1);
    121         psVectorExtend (r2rflux, 100, 1);
    122104        psVectorExtend (dMag,    100, 1);
    123105        psVectorExtend (apResid, 100, 1);
    124106        Npsf ++;
    125     }
    126 
    127     if (psTraceGetLevel("psphot") >= 5) {
    128         fclose (dumpFile);
    129107    }
    130108   
    … …  
    139117    }
    140118
    141     // APTREND options : NONE CONSTANT SKYBIAS XY_LIN XY_QUAD SKY_XY_LIN FULL
    142     // APTREND options are used in the switch block below
    143     pmPSFApTrendOptions ApTrendOption = DEFAULT_OPTION;
    144     char *optionName = psMetadataLookupPtr (&status, recipe, "APTREND");
    145     if (status) ApTrendOption = pmPSFApTrendOptionFromName (optionName);
    146     if (ApTrendOption == PM_PSF_APTREND_ERROR) {
    147         psError(PSPHOT_ERR_APERTURE, true, "invalid aperture residual trend %s", optionName);
    148         return false;
    149     }
    150 
     119    // XXX deprecating the old code which allowed the ApResid to be fitted as a function of flux and r^2/flux
     120    // XXX is this asymmetric clipping still needed?  this analysis should come after neighbors are subtracted...
    151121    // 3hi/1lo sigma clipping on the rflux vs metric fit
    152     psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    153     stats->min = 2.0;
    154     stats->max = 3.0;
    155 
    156 #define P_APTREND_SWITCH_CLEANUP        /* Cleanup memory on error in ApTrendOption switch */ \
    157     psFree(psf->growth); psf->growth = NULL; \
    158     psFree(mask); \
    159     psFree(xPos); \
    160     psFree(yPos); \
    161     psFree(flux); \
    162     psFree(r2rflux); \
    163     psFree(apResid); \
    164     psFree(dMag); \
    165     psFree(stats)
    166 
    167     // no correction
    168     switch (ApTrendOption) {
    169       case PM_PSF_APTREND_NONE:
    170         // remove ApTrend fit from pmPSFtry
    171         psf->ApTrend->coeff[0][0][0][0] = 0;
    172         break;
    173       case PM_PSF_APTREND_CONSTANT:
    174         stats->clipIter = 2;
    175         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    176         if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    177             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
    178             P_APTREND_SWITCH_CLEANUP;
    179             return false;
    180         }
    181         // apply the fit
    182         stats->clipIter = 3;
    183         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    184         if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    185             psError(PSPHOT_ERR_PHOTOM, false, "Fitting aperture correction");
    186             P_APTREND_SWITCH_CLEANUP;
    187             return false;
    188         }
    189         break;
    190       case PM_PSF_APTREND_SKYBIAS:
    191         // first clip out objects which are too far from the median
    192         stats->clipIter = 2;
    193         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    194         if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    195             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
    196             P_APTREND_SWITCH_CLEANUP;
    197             return false;
    198         }
    199         // apply the fit
    200         stats->clipIter = 3;
    201         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
    202         if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    203             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias");
    204             P_APTREND_SWITCH_CLEANUP;
    205             return false;
    206         }
    207         break;
    208       case PM_PSF_APTREND_SKYSAT:
    209         // first clip out objects which are too far from the median
    210         stats->clipIter = 2;
    211         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    212         if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    213             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
    214             P_APTREND_SWITCH_CLEANUP;
    215             return false;
    216         }
    217         // apply the fit
    218         stats->clipIter = 2;
    219         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
    220         if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    221             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias");
    222             P_APTREND_SWITCH_CLEANUP;
    223             return false;
    224         }
    225         // apply the fit
    226         stats->clipIter = 3;
    227         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYSAT);
    228         if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    229             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting skysat");
    230             P_APTREND_SWITCH_CLEANUP;
    231             return false;
    232         }
    233         break;
    234       case PM_PSF_APTREND_XY_LIN:
    235         // first clip out objects which are too far from the median
    236         stats->clipIter = 2;
    237         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    238         if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    239             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
    240             P_APTREND_SWITCH_CLEANUP;
    241             return false;
    242         }
    243         // apply the fit
    244         stats->clipIter = 3;
    245         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_XY_LIN);
    246         if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    247             psError(PSPHOT_ERR_PHOTOM, false, "fitting, XY_LIN");
    248             P_APTREND_SWITCH_CLEANUP;
    249             return false;
    250         }
    251         break;
    252       case PM_PSF_APTREND_XY_QUAD:
    253         // first clip out objects which are too far from the median
    254         stats->clipIter = 2;
    255         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    256         if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    257             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
    258             P_APTREND_SWITCH_CLEANUP;
    259             return false;
    260         }
    261         // apply the fit
    262         stats->clipIter = 3;
    263         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_XY_QUAD);
    264         if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    265             psError(PSPHOT_ERR_PHOTOM, false, "Fitting XY_QUAD");
    266             P_APTREND_SWITCH_CLEANUP;
    267             return false;
    268         }
    269         break;
    270       case PM_PSF_APTREND_SKY_XY_LIN:
    271         // first clip out objects which are too far from the median
    272         stats->clipIter = 2;
    273         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    274         if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    275             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
    276             P_APTREND_SWITCH_CLEANUP;
    277             return false;
    278         }
    279         // apply the fit
    280         stats->clipIter = 3;
    281         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKY_XY_LIN);
    282         if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    283             psError(PSPHOT_ERR_PHOTOM, false, "Fitting sky xy_lin");
    284             P_APTREND_SWITCH_CLEANUP;
    285             return false;
    286         }
    287         break;
    288       case PM_PSF_APTREND_SKYSAT_XY_LIN:
    289         // first clip out objects which are too far from the median
    290         stats->clipIter = 2;
    291         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    292         if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    293             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
    294             P_APTREND_SWITCH_CLEANUP;
    295             return false;
    296         }
    297         // apply the fit
    298         stats->clipIter = 3;
    299         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
    300         if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    301             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias");
    302             P_APTREND_SWITCH_CLEANUP;
    303             return false;
    304         }
    305         // apply the fit
    306         stats->clipIter = 3;
    307         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYSAT_XY_LIN);
    308         if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    309             psError(PSPHOT_ERR_PHOTOM, false, "Fitting skyset xy_lin");
    310             P_APTREND_SWITCH_CLEANUP;
    311             return false;
    312         }
    313         break;
    314       case PM_PSF_APTREND_ALL:
    315         // first clip out objects which are too far from the median
    316         stats->clipIter = 2;
    317         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    318         if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    319             psError(PSPHOT_ERR_PHOTOM, false, "Failed to measure apTrend");
    320             P_APTREND_SWITCH_CLEANUP;
    321             return false;
    322         }
    323         // fit just SkyBias and clip out objects which are too far from the median
    324         stats->clipIter = 2;
    325         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
    326         if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    327             psError(PSPHOT_ERR_PHOTOM, false, "fitting skyBias");
    328             P_APTREND_SWITCH_CLEANUP;
    329             return false;
    330         }
    331         // finally, fit x, y, SkyBias and clip out objects which are too far from the median
    332         stats->clipIter = 3;
    333         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_ALL);
    334         if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
    335             psError(PSPHOT_ERR_PHOTOM, false, "fitting all");
    336             P_APTREND_SWITCH_CLEANUP;
    337             return false;
    338         }
    339         break;
    340       default:
    341         psError(PSPHOT_ERR_PHOTOM, true, "Unknown APTREND value: %s", optionName);
    342         return false;
    343     }
    344 #undef P_APTREND_SWITCH_CLEANUP
     122    // systematic error information
     123    float errorScale = 0.0;
     124    float errorFloor = 0.0;
     125
     126    float errorFloorMin = FLT_MAX;
     127    int entryMin = -1;
     128
     129    // *** iterate over spatial scale until error Floor increases
     130    // *** fit out the dap vs mag trend?
     131    // *** stop if Npsf / (Nx * Ny) < 3
     132    for (int i = 1; i < 10; i++) {
     133
     134        if (!psphotApResidTrend (readout, psf, Npsf, i, &errorScale, &errorFloor, mask, xPos, yPos, apResid, dMag)) {
     135            break;
     136        }
     137
     138        // store the resulting errorFloor values and the scales, redo the best
     139        if (errorFloor < errorFloorMin) {
     140            errorFloorMin = errorFloor;
     141            entryMin = i;
     142        }
     143    }
     144    if (entryMin == -1) {
     145        psAbort ("failed on ApResid Trend");
     146    }
     147
     148    psphotApResidTrend (readout, psf, Npsf, entryMin, &errorScale, &errorFloor, mask, xPos, yPos, apResid, dMag);
    345149
    346150    // construct the fitted values and the residuals
    347     psVector *fit   = psPolynomial4DEvalVector (psf->ApTrend, xPos, yPos, r2rflux, flux);
    348     psVector *resid = (psVector *) psBinaryOp (NULL, (void *) apResid, "-", (void *) fit);
    349 
    350     // measure scatter for sources with dMag < 0.01 (S/N = 100)
    351     int Nkeep = 0;
    352     psStats *residStats = psStatsAlloc (PS_STAT_SAMPLE_STDEV);
    353     for (int i = 0; i < dMag->n; i++) {
    354         if (dMag->data.F64[i] > 0.01) {
    355             mask->data.U8[i] |= 0x02;
    356         }
    357         if (! mask->data.U8[i]) Nkeep ++;
    358     }
    359     psVectorStats (residStats, resid, NULL, mask, 0x03);
     151    psVector *apResidFit = pmTrend2DEvalVector (psf->ApTrend, xPos, yPos);
     152    psVector *apResidRes = (psVector *) psBinaryOp (NULL, (void *) apResid, "-", (void *) apResidFit);
     153    psVector *dMagSys = (psVector *) psBinaryOp (NULL, (void *) dMag, "*", (void *) psScalarAlloc(errorScale, PS_TYPE_F32));
     154
     155    psphotSaveImage (NULL, psf->ApTrend->map->map, "apresid.map.fits");
     156
     157    if (psTraceGetLevel("psphot") >= 5) {
     158        FILE *dumpFile = fopen ("apresid.dat", "w");
     159        for (int i = 0; i < xPos->n; i++) {
     160            fprintf (dumpFile, "%f %f  %f %f %f  %f %f  %x\n",
     161                     xPos->data.F32[i], yPos->data.F32[i],
     162                     mag->data.F32[i], dMag->data.F32[i], dMagSys->data.F32[i],
     163                     apResid->data.F32[i], apResidRes->data.F32[i],
     164                     mask->data.U8[i]);
     165        }
     166        fclose (dumpFile);
     167    }
    360168
    361169    // apply ApTrend results
    362     psf->skyBias  = psf->ApTrend->coeff[0][0][1][0];
    363     psf->skySat   = psf->ApTrend->coeff[0][0][0][1];
    364     psf->ApResid  = psf->ApTrend->coeff[0][0][0][0];
    365     psf->dApResid = residStats->sampleStdev;
    366     psf->ApTrend->coeff[0][0][1][0] = 0;
    367     psf->ApTrend->coeff[0][0][0][1] = 0;
     170    float xc = 0.5*readout->image->numCols + readout->image->col0 + 0.5;
     171    float yc = 0.5*readout->image->numRows + readout->image->row0 + 0.5;
     172    psf->ApResid  = pmTrend2DEval (psf->ApTrend, xc, yc); // ap-fit at chip center
     173    psf->dApResid = errorFloor;
    368174    psf->nApResid = Npsf;
    369175
    370176    // save results for later output
    371     psMetadataAdd (recipe, PS_LIST_TAIL, "SKYBIAS",  PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   psf->skyBias);
    372     psMetadataAdd (recipe, PS_LIST_TAIL, "SKYSAT",   PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   psf->skySat);
     177    psMetadataAdd (recipe, PS_LIST_TAIL, "SKYBIAS",  PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   0.0);
     178    psMetadataAdd (recipe, PS_LIST_TAIL, "SKYSAT",   PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   0.0);
    373179    psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   psf->ApResid);
    374180    psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->dApResid);
    … …  
    376182    psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "ap residual scatter", psf->nApResid);
    377183
    378     psLogMsg ("psphot.apresid", PS_LOG_INFO, "measure full-frame aperture residuals for %d of %d objects: %f sec\n",
    379               Nkeep, Npsf, psTimerMark ("psphot"));
    380     psLogMsg ("psphot.apresid", PS_LOG_DETAIL, "aperture residual: %f +/- %f : %f bias, %f skysat\n",
    381               psf->ApResid, psf->dApResid, psf->skyBias, psf->skySat);
    382     psLogMsg ("psphot.apresid", PS_LOG_MINUTIA, "apresid trends: %f %f %f %f %f\n",
    383               1e3*psf->ApTrend->coeff[1][0][0][0],
    384               1e6*psf->ApTrend->coeff[2][0][0][0],
    385               1e6*psf->ApTrend->coeff[1][1][0][0],
    386               1e3*psf->ApTrend->coeff[0][1][0][0],
    387               1e6*psf->ApTrend->coeff[0][2][0][0]);
    388 
     184    // psLogMsg ("psphot.apresid", PS_LOG_INFO, "measure full-frame aperture residuals for %d of %d objects: %f sec\n", Nkeep, Npsf, psTimerMark ("psphot"));
     185    psLogMsg ("psphot.apresid", PS_LOG_DETAIL, "aperture residual: %f +/- %f\n", psf->ApResid, psf->dApResid);
     186
     187    psFree (mag);
    389188    psFree (mask);
    390189    psFree (xPos);
    391190    psFree (yPos);
    392     psFree (flux);
    393     psFree (r2rflux);
    394191    psFree (apResid);
    395192    psFree (dMag);
    396193
    397     psFree (fit);
    398     psFree (resid);
    399     psFree (stats);
    400     psFree (residStats);
    401194    return true;
    402195}
    … …  
    408201*/
    409202
     203bool psphotMagErrorScale (float *errorScale, float *errorFloor, psVector *dMag, psVector *dap, int nGroup) {
     204
     205    psStats *statsS = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
     206    psStats *statsM = psStatsAlloc (PS_STAT_SAMPLE_MEAN);
     207
     208    // measure the trend in bins with 10 values each; if < 10 total, use them all
     209    int nBin = PS_MAX (dMag->n / nGroup, 1);
     210
     211    // output vectors for ApResid trend
     212    psVector *dSo = psVectorAlloc (nBin, PS_TYPE_F32);
     213    psVector *dMo = psVectorAlloc (nBin, PS_TYPE_F32);
     214    psVector *dRo = psVectorAlloc (nBin, PS_TYPE_F32);
     215
     216    // use dMag to group the dMag and dap vectors
     217    psVector *index = psVectorSortIndex (NULL, dMag);
     218
     219    // subset vectors for dMag and dap values within the given range
     220    psVector *dMSubset = psVectorAllocEmpty (nGroup, PS_TYPE_F32);
     221    psVector *dASubset = psVectorAllocEmpty (nGroup, PS_TYPE_F32);
     222
     223    int n = 0;
     224    for (int i = 0; i < dMo->n; i++) {
     225        int j;
     226        for (j = 0; (j < nGroup) && (n < dMag->n); j++, n++) {
     227            int N = index->data.U32[n];
     228            dMSubset->data.F32[j] = dMag->data.F32[N];
     229            dASubset->data.F32[j] = dap->data.F32[N];
     230        }
     231        dMSubset->n = j;
     232        dASubset->n = j;
     233
     234        psStatsInit (statsS);
     235        psStatsInit (statsM);
     236
     237        psVectorStats (statsS, dASubset, NULL, NULL, 0xff);
     238        psVectorStats (statsM, dMSubset, NULL, NULL, 0xff);
     239
     240        dSo->data.F32[i] = statsS->robustStdev;
     241        dMo->data.F32[i] = statsM->sampleMean;
     242
     243        dRo->data.F32[i] = statsS->robustStdev / statsM->sampleMean;
     244    }
     245    psFree (dMSubset);
     246    psFree (dASubset);
     247   
     248    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
     249    psVectorStats (stats, dRo, NULL, NULL, 0);
     250
     251    *errorScale = stats->sampleMedian;
     252    *errorFloor = dSo->data.F32[0];
     253
     254    psFree (stats);
     255    psFree (index);
     256
     257    psFree (dRo);
     258    psFree (dMo);
     259    psFree (dSo);
     260   
     261    psFree (statsS);
     262    psFree (statsM);
     263
     264    return true;
     265}
     266
     267bool psphotApResidTrend (pmReadout *readout, pmPSF *psf, int Npsf, int scale, float *errorScale, float *errorFloor, psVector *mask, psVector *xPos, psVector *yPos, psVector *apResid, psVector *dMag) {
     268
     269    int Nx, Ny;
     270       
     271    psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
     272    stats->min = 2.0;
     273    stats->max = 3.0;
     274
     275    if (readout->image->numCols > readout->image->numRows) {
     276        Nx = scale;
     277        Ny = PS_MAX (1, Nx * (readout->image->numRows / readout->image->numCols));
     278    } else {
     279        Ny = scale;
     280        Nx = PS_MAX (1, Ny * (readout->image->numCols / readout->image->numRows));
     281    }
     282
     283    if (Npsf < 3*Nx*Ny) {
     284        return false;
     285    }
     286
     287    // measure Trend2D for the current spatial scale
     288    psFree (psf->ApTrend);
     289    psf->ApTrend = pmTrend2DAlloc (PM_TREND_MAP, readout->image, Nx, Ny, stats);
     290    pmTrend2DFit (psf->ApTrend, mask, 0xff, xPos, yPos, apResid, dMag);
     291   
     292    // construct the fitted values and the residuals
     293    psVector *apResidFit = pmTrend2DEvalVector (psf->ApTrend, xPos, yPos);
     294    psVector *apResidRes = (psVector *) psBinaryOp (NULL, (void *) apResid, "-", (void *) apResidFit);
     295
     296    // measure systematic errorFloor & systematic / photon scale factor
     297    // XXX this is a bit arbitrary, but it forces ~3 stars from the bright bin per spatial bin
     298    int nGroup = PS_MAX (3*Nx*Ny, 10);
     299    psphotMagErrorScale (errorScale, errorFloor, dMag, apResidRes, nGroup);
     300
     301    psLogMsg ("psphot.apresid", PS_LOG_INFO, "result of %d x %d grid (%d stars per bin)\n", Nx, Ny, nGroup);
     302    psLogMsg ("psphot.apresid", PS_LOG_INFO, "systematic error / photon error: %f\n", *errorScale);
     303    psLogMsg ("psphot.apresid", PS_LOG_INFO, "systematic scatter floor: %f\n", *errorFloor);
     304
     305    psFree (stats);
     306    psFree (apResidFit);
     307    psFree (apResidRes);
     308
     309    return true;
     310}
     311   
Note: See TracChangeset for help on using the changeset viewer.