IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 25616


Ignore:
Timestamp:
Sep 27, 2009, 11:15:32 AM (17 years ago)
Author:
eugene
Message:

clarify some comments; drop some deprecated code (old Unthreaded); skip DEFECTs in ApResid analysis; complete rework of trend analysis (use pmVectorSystematicError); drop old SKYBIAS, SKYSAT references; be careful to mask the desired aperture

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/20090715/psphot/src/psphotApResid.c

    r25536 r25616  
    3333    if (!measureAptrend) {
    3434        // save nan values since these were not calculated
    35         psMetadataAdd (recipe, PS_LIST_TAIL, "SKYBIAS",  PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   NAN);
    36         psMetadataAdd (recipe, PS_LIST_TAIL, "SKYSAT",   PS_DATA_F32 | PS_META_REPLACE, "aperture-determined saturation",   NAN);
    3735        psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   NAN);
    3836        psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", NAN);
     37        psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", 0);
    3938        psMetadataAdd (recipe, PS_LIST_TAIL, "APLOSS",   PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", NAN);
    40         psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", 0);
    4139        return true;
    4240    }
    … …  
    5351    maskVal |= markVal;
    5452
    55     // S/N limit to perform full non-linear fits
     53    // clipping for extreme outliers
     54    // XXX this is not currently defined in the recipe
    5655    float MAX_AP_OFFSET = psMetadataLookupF32 (&status, recipe, "MAX_AP_OFFSET");
    5756
    58     // this is the smallest radius allowed: need to at least extend growth curve down to this...
     57    // options for how the photometry is calculated
     58    // XXX are these sensible?
    5959    bool IGNORE_GROWTH = psMetadataLookupBool (&status, recipe, "IGNORE_GROWTH");
    6060    bool INTERPOLATE_AP = psMetadataLookupBool (&status, recipe, "INTERPOLATE_AP");
    … …  
    100100            PS_ARRAY_ADD_SCALAR(job->args, photMode, PS_TYPE_S32);
    101101            PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
     102            PS_ARRAY_ADD_SCALAR(job->args, markVal,  PS_TYPE_IMAGE_MASK);
    102103
    103104            PS_ARRAY_ADD_SCALAR(job->args, 0,        PS_TYPE_S32); // this is used as a return value for Nskip
    … …  
    110111            }
    111112            psFree(job);
    112 
    113 # if (0)
    114                 int nskip = 0;
    115                 int nfail = 0;
    116 
    117                 if (!psphotApResidMags_Unthreaded (&nskip, &nfail, sources, psf, photMode, maskVal)) {
    118                     psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    119                     return false;
    120                 }
    121                 Nskip += nskip;
    122                 Nfail += nfail;
    123 # endif
    124 
    125113        }
    126114
    … …  
    138126            } else {
    139127                psScalar *scalar = NULL;
    140                 scalar = job->args->data[4];
     128                scalar = job->args->data[5];
    141129                Nskip += scalar->data.S32;
    142                 scalar = job->args->data[5];
     130                scalar = job->args->data[6];
    143131                Nfail += scalar->data.S32;
    144132            }
    … …  
    150138
    151139    // gather the stats to assess the aperture residuals
    152     psVector *mask    = psVectorAllocEmpty (300, PS_TYPE_VECTOR_MASK);
    153140    psVector *mag     = psVectorAllocEmpty (300, PS_TYPE_F32);
    154141    psVector *xPos    = psVectorAllocEmpty (300, PS_TYPE_F32);
    … …  
    158145    Npsf = 0;
    159146
     147    FILE *f = fopen ("apresid.dat", "w");
     148    psAssert (f, "failed open");
     149
    160150    for (int i = 0; i < sources->n; i++) {
    161151        source = sources->data[i];
    … …  
    170160        if (source->mode &  PM_SOURCE_MODE_EXT_LIMIT) SKIPSTAR ("EXTENDED");
    171161        if (source->mode &  PM_SOURCE_MODE_CR_LIMIT) SKIPSTAR ("COSMIC RAY");
     162        if (source->mode &  PM_SOURCE_MODE_DEFECT) SKIPSTAR ("DEFECT");
    172163           
    173164        if (!isfinite(source->apMag) || !isfinite(source->psfMag)) {
    174165            continue;
    175166        }
     167
     168        // XXX make this user-configurable?
     169        if (source->errMag > 0.01) continue;
    176170
    177171        // aperture residual for this source
    … …  
    185179        }
    186180
    187         mag->data.F32[Npsf]     = source->psfMag;
    188         apResid->data.F32[Npsf] = dap;
    189         xPos->data.F32[Npsf]    = model->params->data.F32[PM_PAR_XPOS];
    190         yPos->data.F32[Npsf]    = model->params->data.F32[PM_PAR_YPOS];
    191 
    192         mask->data.PS_TYPE_VECTOR_MASK_DATA[Npsf] = 0;
    193 
    194         // XXX don't I have errMag at this point??
    195         dMag->data.F32[Npsf] = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
    196 
    197         psVectorExtend (mag,     100, 1);
    198         psVectorExtend (mask,    100, 1);
    199         psVectorExtend (xPos,    100, 1);
    200         psVectorExtend (yPos,    100, 1);
    201         psVectorExtend (dMag,    100, 1);
    202         psVectorExtend (apResid, 100, 1);
     181        fprintf (f, "%6.1f %6.1f : %6.1f %6.1f : %8.3f %8.3f %8.3f : %f : %f %f %f : %f\n",
     182                 source->peak->xf, source->peak->yf,
     183                 source->modelPSF->params->data.F32[PM_PAR_XPOS], source->modelPSF->params->data.F32[PM_PAR_YPOS],
     184                 source->psfMag, source->apMag, source->errMag,
     185                 source->modelPSF->params->data.F32[PM_PAR_I0],
     186                 source->modelPSF->params->data.F32[PM_PAR_SXX], source->modelPSF->params->data.F32[PM_PAR_SXY], source->modelPSF->params->data.F32[PM_PAR_SYY],
     187                 source->modelPSF->params->data.F32[PM_PAR_7]);
     188
     189        psVectorAppend (mag, source->psfMag);
     190        psVectorAppend (dMag,source->errMag);
     191        psVectorAppend (apResid, dap);
     192        psVectorAppend (xPos, model->params->data.F32[PM_PAR_XPOS]);
     193        psVectorAppend (yPos, model->params->data.F32[PM_PAR_YPOS]);
    203194        Npsf ++;
    204195    }
    … …  
    206197    psLogMsg ("psphot.apresid", PS_LOG_DETAIL, "measure aperture residuals for %d objects (%d skipped, %d failed, %ld invalid)\n",
    207198              Npsf, Nskip, Nfail, sources->n - Npsf - Nskip - Nfail);
     199
     200    fclose (f);
    208201
    209202    // XXX choose a better value here?
    … …  
    213206    }
    214207
    215     // XXX deprecating the old code which allowed the ApResid to be fitted as a function of flux and r^2/flux
    216     // XXX is this asymmetric clipping still needed?  this analysis should come after neighbors are subtracted...
    217     // 3hi/1lo sigma clipping on the rflux vs metric fit
    218     // systematic error information
    219     float errorScale = 0.0;
    220     float errorFloor = 0.0;
    221 
     208    // XXX set the min number of needed source more carefully
     209    if ((Npsf < 15) && (APTREND_ORDER_MAX >= 4)) APTREND_ORDER_MAX = 3;
     210    if ((Npsf < 11) && (APTREND_ORDER_MAX >= 3)) APTREND_ORDER_MAX = 2;
     211    if ((Npsf <  8) && (APTREND_ORDER_MAX >= 2)) APTREND_ORDER_MAX = 1;
     212
     213    psFree (psf->ApTrend);
     214    psf->ApTrend = NULL;
    222215    float errorFloorMin = FLT_MAX;
    223     int entryMin = -1;
    224 
    225     // Fit out the dap vs mag trend, iterate over spatial scale until error Floor increases.
    226     // Stop if Npsf / (Nx * Ny) < 3
     216
     217    // as we loop over orders, we need to refer to the initial selection, but we modify the
     218    // option values to match the current guess: save the max values here:
     219    int NX = readout->image->numCols;
     220    int NY = readout->image->numRows;
    227221    for (int i = 1; i <= APTREND_ORDER_MAX; i++) {
    228 
    229         if (!psphotApResidTrend (readout, psf, Npsf, i, &errorScale, &errorFloor, mask, xPos, yPos, apResid, dMag)) {
    230             break;
    231         }
    232 
    233         // store the resulting errorFloor values and the scales, redo the best
     222       
     223        int Nx, Ny;
     224        if (NX > NY) {
     225            Nx = i;
     226            Ny = PS_MAX (1, (int)(i * (NY / (float)(NX)) + 0.5));
     227        } else {
     228            Ny = i;
     229            Nx = PS_MAX (1, (int)(i * (NX / (float)(NY)) + 0.5));
     230        }
     231
     232        float errorFloor;
     233        pmTrend2D *apTrend = psphotApResidTrend (&errorFloor, readout, Nx, Ny, xPos, yPos, apResid, dMag);
     234        if (!apTrend) {
     235            continue;
     236        }
     237
     238        // store the minimum errorFloor and best ApTrend to keep
    234239        if (errorFloor < errorFloorMin) {
    235240            errorFloorMin = errorFloor;
    236             entryMin = i;
     241            psFree (psf->ApTrend);
     242            psf->ApTrend = psMemIncrRefCounter(apTrend);
    237243        }
    238     }
    239     if (entryMin == -1) {
     244        psFree (apTrend);
     245    }
     246    if (psf->ApTrend == NULL) {
    240247        psWarning("Failed to find a valid aperture residual value");
    241248        goto escape;
    242249    }
    243250
    244     // XXX catch error condition
    245     psphotApResidTrend (readout, psf, Npsf, entryMin, &errorScale, &errorFloor, mask, xPos, yPos, apResid, dMag);
     251    // apply ApTrend results
     252    float xc = 0.5*readout->image->numCols + readout->image->col0 + 0.5;
     253    float yc = 0.5*readout->image->numRows + readout->image->row0 + 0.5;
     254
     255    psf->ApResid  = pmTrend2DEval (psf->ApTrend, xc, yc); // ap-fit at chip center
     256    psf->dApResid = errorFloorMin;
     257    psf->nApResid = Npsf;
     258
     259    // save results for later output
     260    psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   psf->ApResid);
     261    psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->dApResid);
     262    psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", psf->nApResid);
     263    psMetadataAdd (recipe, PS_LIST_TAIL, "APLOSS",   PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", psf->growth->apLoss);
     264
     265    psLogMsg ("psphot.apresid", PS_LOG_DETAIL, "aperture residual: %f +/- %f\n", psf->ApResid, psf->dApResid);
     266    psLogMsg ("psphot.apresid", PS_LOG_INFO, "measure full-frame aperture residuals for %d sources: %f sec\n", Npsf, psTimerMark ("psphot.apresid"));
     267
     268    psFree (xPos);
     269    psFree (yPos);
     270    psFree (apResid);
     271    psFree (mag);
     272    psFree (dMag);
     273
     274    psphotVisualPlotApResid (sources, psf->ApResid, psf->dApResid);
     275
     276    return true;
     277
     278escape:
     279    // save nan values since these were not calculated
     280    psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   NAN);
     281    psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", NAN);
     282    psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", 0);
     283    psMetadataAdd (recipe, PS_LIST_TAIL, "APLOSS",   PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", NAN);
     284
     285    psFree (xPos);
     286    psFree (yPos);
     287    psFree (apResid);
     288    psFree (mag);
     289    psFree (dMag);
     290    return false;
     291}
     292
     293pmTrend2D *psphotApResidTrend (float *apResidSysErr, pmReadout *readout, int Nx, int Ny, psVector *xPos, psVector *yPos, psVector *apResid, psVector *dMag) {
     294
     295    // the mask marks the values not used to calculate the ApTrend
     296    psVector *mask = psVectorAlloc(xPos->n, PS_TYPE_VECTOR_MASK);
     297    psVectorInit (mask, 0);
     298
     299    // XXX allow user to set this optionally?
     300    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
     301
     302    // measure Trend2D for the current spatial scale
     303    pmTrend2D *apTrend = pmTrend2DAlloc (PM_TREND_MAP, readout->image, Nx, Ny, stats);
     304
     305    // XXX somewhat arbitrary: soften the errors so the few bright stars do not totally dominate:
     306    // XXX use this or not?  probably not, since this is the point of the systematic error analysis
     307    psVector *dMagSoft = psVectorAlloc (dMag->n, PS_TYPE_F32);
     308    for (int i = 0; i < dMag->n; i++) {
     309        dMagSoft->data.F32[i] = hypot(dMag->data.F32[i], 0.005);
     310    }
     311
     312    // XXX test for errors here
     313    if (!pmTrend2DFit (apTrend, mask, 0xff, xPos, yPos, apResid, dMagSoft)) {
     314        psWarning("Failed to fit trend for %d x %d map", Nx, Ny);
     315        psFree (apTrend);
     316        return NULL;
     317    }
    246318
    247319    // construct the fitted values and the residuals
    248     psVector *apResidFit = pmTrend2DEvalVector (psf->ApTrend, mask, 0xff, xPos, yPos);
     320    psVector *apResidFit = pmTrend2DEvalVector (apTrend, mask, 0xff, xPos, yPos);
    249321    psVector *apResidRes = (psVector *) psBinaryOp (NULL, (void *) apResid, "-", (void *) apResidFit);
    250     psVector *dMagSys = (psVector *) psBinaryOp (NULL, (void *) dMag, "*", (void *) psScalarAlloc(errorScale, PS_TYPE_F32));
    251 
    252     if (psTraceGetLevel("psphot") >= 2) {
    253         FILE *dumpFile = fopen ("apresid.dat", "w");
     322
     323    // measure systematic error
     324    *apResidSysErr = psVectorSystematicError (apResidRes, dMag, 0.10);
     325    if (!isfinite(*apResidSysErr)) {
     326        psWarning("Failed to find systematic error for %d x %d map", Nx, Ny);
     327        psFree (apTrend);
     328        return NULL;
     329    }
     330
     331    psLogMsg ("psphot.apresid", PS_LOG_INFO, "result of %d x %d grid\n", Nx, Ny);
     332    psLogMsg ("psphot.apresid", PS_LOG_INFO, "systematic scatter floor: %f\n", *apResidSysErr);
     333
     334    if (psTraceGetLevel("psphot") >= 4) {
     335        char filename[64];
     336        snprintf (filename, 64, "apresid.%dx%d.dat", Nx, Ny);
     337        FILE *dumpFile = fopen (filename, "w");
    254338        for (int i = 0; i < xPos->n; i++) {
    255             fprintf (dumpFile, "%f %f  %f %f %f %f %f  %x\n",
     339            fprintf (dumpFile, "%f %f  %f %f %f %f  %x\n",
    256340                     xPos->data.F32[i], yPos->data.F32[i],
    257                      mag->data.F32[i], dMag->data.F32[i], dMagSys->data.F32[i],
     341                     dMag->data.F32[i], hypot(dMag->data.F32[i], *apResidSysErr),
    258342                     apResid->data.F32[i], apResidRes->data.F32[i],
    259343                     mask->data.PS_TYPE_VECTOR_MASK_DATA[i]);
    … …  
    262346    }
    263347
    264     // apply ApTrend results
    265     float xc = 0.5*readout->image->numCols + readout->image->col0 + 0.5;
    266     float yc = 0.5*readout->image->numRows + readout->image->row0 + 0.5;
    267 
    268     psf->ApResid  = pmTrend2DEval (psf->ApTrend, xc, yc); // ap-fit at chip center
    269     psf->dApResid = errorFloor;
    270     psf->nApResid = Npsf;
    271 
    272     // save results for later output
    273     psMetadataAdd (recipe, PS_LIST_TAIL, "SKYBIAS",  PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   0.0);
    274     psMetadataAdd (recipe, PS_LIST_TAIL, "SKYSAT",   PS_DATA_F32 | PS_META_REPLACE, "aperture-determined saturation",   0.0);
    275     psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   psf->ApResid);
    276     psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->dApResid);
    277     psMetadataAdd (recipe, PS_LIST_TAIL, "APLOSS",   PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", psf->growth->apLoss);
    278     psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", psf->nApResid);
    279 
    280     psLogMsg ("psphot.apresid", PS_LOG_DETAIL, "aperture residual: %f +/- %f\n", psf->ApResid, psf->dApResid);
    281     psLogMsg ("psphot.apresid", PS_LOG_INFO, "measure full-frame aperture residuals for %d sources: %f sec\n", Npsf, psTimerMark ("psphot.apresid"));
    282 
    283     psFree (mag);
    284348    psFree (mask);
    285     psFree (xPos);
    286     psFree (yPos);
    287 
    288     psFree (apResid);
    289     psFree (apResidFit);
    290     psFree (apResidRes);
    291 
    292     psFree (dMagSys);
    293     psFree (dMag);
    294 
    295     psphotVisualPlotApResid (sources);
    296 
    297     return true;
    298 
    299 escape:
    300     // save nan values since these were not calculated
    301     psMetadataAdd (recipe, PS_LIST_TAIL, "SKYBIAS",  PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   NAN);
    302     psMetadataAdd (recipe, PS_LIST_TAIL, "SKYSAT",   PS_DATA_F32 | PS_META_REPLACE, "aperture-determined saturation",   NAN);
    303     psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   NAN);
    304     psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", NAN);
    305     psMetadataAdd (recipe, PS_LIST_TAIL, "APLOSS",   PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", NAN);
    306     psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", 0);
    307 
    308     psFree (mag);
    309     psFree (mask);
    310     psFree (xPos);
    311     psFree (yPos);
    312     psFree (apResid);
    313     psFree (dMag);
    314     return false;
    315 }
    316 
    317 /*
    318   (aprMag' - fitMag) = rflux*skyBias + ApTrend(x,y)
    319   (aprMag - rflux*skyBias) - fitMag = ApTrend(x,y)
    320   (aprMag - rflux*skyBias) = fitMag + ApTrend(x,y)
    321 */
    322 
    323 // XXX this still sucks...  need a better way to estimate the error floor...
    324 bool psphotMagErrorScale (float *errorScale, float *errorFloor, psVector *dMag, psVector *dap, psVector *mask, int nGroup) {
    325 
    326     psStats *statsS = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    327     psStats *statsM = psStatsAlloc (PS_STAT_SAMPLE_MEAN);
    328 
    329     // measure the trend in bins with 10 values each; if < 10 total, use them all
    330     int nBin = PS_MAX (dMag->n / nGroup, 1);
    331 
    332     // output vectors for ApResid trend
    333     psVector *dSo = psVectorAlloc (nBin, PS_TYPE_F32);
    334     psVector *dMo = psVectorAlloc (nBin, PS_TYPE_F32);
    335     psVector *dRo = psVectorAlloc (nBin, PS_TYPE_F32);
    336 
    337     // use dMag to group the dMag and dap vectors
    338     psVector *index = psVectorSortIndex (NULL, dMag);
    339 
    340     // subset vectors for dMag and dap values within the given range
    341     psVector *dMSubset = psVectorAllocEmpty (nGroup, PS_TYPE_F32);
    342     psVector *dASubset = psVectorAllocEmpty (nGroup, PS_TYPE_F32);
    343     psVector *mkSubset = psVectorAllocEmpty (nGroup, PS_TYPE_VECTOR_MASK);
    344 
    345     int n = 0;
    346     for (int i = 0; i < dMo->n; i++) {
    347         int j;
    348         for (j = 0; (j < nGroup) && (n < dMag->n); j++, n++) {
    349             int N = index->data.U32[n];
    350             dMSubset->data.F32[j] = dMag->data.F32[N];
    351             dASubset->data.F32[j] = dap->data.F32[N];
    352             mkSubset->data.PS_TYPE_VECTOR_MASK_DATA[j] = mask->data.PS_TYPE_VECTOR_MASK_DATA[N];
    353         }
    354         dMSubset->n = j;
    355         dASubset->n = j;
    356         mkSubset->n = j;
    357 
    358         psStatsInit (statsS);
    359         psStatsInit (statsM);
    360 
    361         if (j > 2) {
    362             if (!psVectorStats (statsS, dASubset, NULL, mkSubset, 0xff)) {
    363                 psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
    364                 return false;
    365             }
    366             if (!psVectorStats (statsM, dMSubset, NULL, mkSubset, 0xff)) {
    367                 psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
    368                 return false;
    369             }
    370             dSo->data.F32[i] = statsS->robustStdev;
    371             dMo->data.F32[i] = statsM->sampleMean;
    372             dRo->data.F32[i] = statsS->robustStdev / statsM->sampleMean;
    373         } else {
    374             dSo->data.F32[i] = NAN;
    375             dMo->data.F32[i] = NAN;
    376             dRo->data.F32[i] = NAN;
    377         }
    378     }
    379     psFree (dMSubset);
    380     psFree (dASubset);
    381     psFree (mkSubset);
    382 
    383     psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
    384     if (!psVectorStats (stats, dRo, NULL, NULL, 0)) {
    385         psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
    386         return false;
    387     }
    388 
    389     *errorScale = stats->sampleMedian;
    390     for (int i = 0; i < dSo->n; i++) {
    391         *errorFloor = dSo->data.F32[i];
    392         if (fabs(*errorFloor) <= FLT_EPSILON) continue;
    393         if (isfinite(*errorFloor)) break;
    394     }
    395 
    396     psFree (stats);
    397     psFree (index);
    398 
    399     psFree (dRo);
    400     psFree (dMo);
    401     psFree (dSo);
    402 
    403     psFree (statsS);
    404     psFree (statsM);
    405 
    406     return true;
    407 }
    408 
    409 bool psphotApResidTrend (pmReadout *readout, pmPSF *psf, int Npsf, int scale, float *errorScale, float *errorFloor, psVector *mask, psVector *xPos, psVector *yPos, psVector *apResid, psVector *dMag) {
    410 
    411     int Nx, Ny;
    412 
    413     // XXX use the psf trend scale?
    414     if (readout->image->numCols > readout->image->numRows) {
    415         Nx = scale;
    416         float AR = readout->image->numRows / (float) readout->image->numCols;
    417         Ny = (int) (Nx * AR + 0.5);
    418         Ny = PS_MAX (1, Ny);
    419     } else {
    420         Ny = scale;
    421         float AR = readout->image->numRows / (float) readout->image->numCols;
    422         Nx = (int) (Ny * AR + 0.5);
    423         Nx = PS_MAX (1, Nx);
    424     }
    425 
    426     // the mask marks the values not used to calculate the ApTrend
    427     psVectorInit (mask, 0);
    428 
    429     // XXX stats structure for use by ApTrend : make parameters user setable
    430     // XXX another stat?
    431     psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    432     stats->min = 2.0;
    433     stats->max = 3.0;
    434 
    435     // measure Trend2D for the current spatial scale
    436     psFree (psf->ApTrend);
    437     psf->ApTrend = pmTrend2DAlloc (PM_TREND_MAP, readout->image, Nx, Ny, stats);
    438 
    439     // XXX somewhat arbitrary: soften the errors so the few bright stars do not totally dominate:
    440     psVector *dMagSoft = psVectorAlloc (dMag->n, PS_TYPE_F32);
    441     for (int i = 0; i < dMag->n; i++) {
    442         dMagSoft->data.F32[i] = hypot(dMag->data.F32[i], 0.01);
    443     }
    444 
    445     // XXX test for errors here
    446     pmTrend2DFit (psf->ApTrend, mask, 0xff, xPos, yPos, apResid, dMagSoft);
    447 
    448     // construct the fitted values and the residuals
    449     psVector *apResidFit = pmTrend2DEvalVector (psf->ApTrend, mask, 0xff, xPos, yPos);
    450     psVector *apResidRes = (psVector *) psBinaryOp (NULL, (void *) apResid, "-", (void *) apResidFit);
    451 
    452     // measure systematic errorFloor & systematic / photon scale factor
    453     // XXX this is a bit arbitrary, but it forces ~3 stars from the bright bin per spatial bin
    454     int nGroup = PS_MAX (3*Nx*Ny, 10);
    455     psphotMagErrorScale (errorScale, errorFloor, dMag, apResidRes, mask, nGroup);
    456 
    457     psLogMsg ("psphot.apresid", PS_LOG_INFO, "result of %d x %d grid (%d stars per bin)\n", Nx, Ny, nGroup);
    458     psLogMsg ("psphot.apresid", PS_LOG_INFO, "systematic error / photon error: %f\n", *errorScale);
    459     psLogMsg ("psphot.apresid", PS_LOG_INFO, "systematic scatter floor: %f\n", *errorFloor);
    460 
    461349    psFree (stats);
    462350    psFree (dMagSoft);
    … …  
    464352    psFree (apResidRes);
    465353
    466     return true;
     354    return apTrend;
    467355}
    468356
    … …  
    478366    pmSourcePhotometryMode photMode = PS_SCALAR_VALUE(job->args->data[2],S32);
    479367    psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[3],PS_TYPE_IMAGE_MASK_DATA);
     368    psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
    480369
    481370    for (int i = 0; i < sources->n; i++) {
    … …  
    488377        if (source->mode &  PM_SOURCE_MODE_POOR) SKIPSTAR ("POOR STAR");
    489378
    490         if (!pmSourceMagnitudes (source, psf, photMode, maskVal)) {
    491             Nskip ++;
    492             psTrace ("psphot", 3, "skip : bad source mag");
    493             continue;
     379        // replace object in image
     380        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
     381            pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    494382        }
    495383
    496         if (!isfinite(source->apMag) || !isfinite(source->psfMag)) {
    497             Nfail ++;
    498             psTrace ("psphot", 3, "fail : nan mags : %f %f", source->apMag, source->psfMag);
    499             continue;
    500         }
    501         source->mode |= PM_SOURCE_MODE_AP_MAGS;
     384        // clear the mask bit and set the circular mask pixels
     385        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
     386        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", markVal);
     387
     388        bool status = pmSourceMagnitudes (source, psf, photMode, maskVal);
     389
     390        // clear the mask bit
     391        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
     392
     393        // re-subtract the object, leave local sky
     394        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     395
     396        if (!status) {
     397            Nskip ++;
     398            psTrace ("psphot", 3, "skip : bad source mag");
     399            continue;
     400        }
     401   
     402        if (!isfinite(source->apMag) || !isfinite(source->psfMag)) {
     403            Nfail ++;
     404            psTrace ("psphot", 3, "fail : nan mags : %f %f", source->apMag, source->psfMag);
     405            continue;
     406        }
     407        source->mode |= PM_SOURCE_MODE_AP_MAGS;
    502408    }
    503409
    504410    // change the value of a scalar on the array (wrap this and put it in psArray.h)
    505     scalar = job->args->data[4];
     411    scalar = job->args->data[5];
    506412    scalar->data.S32 = Nskip;
    507413
    508     scalar = job->args->data[5];
     414    scalar = job->args->data[6];
    509415    scalar->data.S32 = Nfail;
    510416
    511417    return true;
    512418}
    513 
    514 bool psphotApResidTrend_old (pmReadout *readout, pmPSF *psf, int Npsf, int scale, float *errorScale, float *errorFloor, psVector *mask, psVector *xPos, psVector *yPos, psVector *apResid, psVector *dMag) {
    515 
    516     int Nx, Ny;
    517 
    518     if (readout->image->numCols > readout->image->numRows) {
    519         Nx = scale;
    520         float AR = readout->image->numRows / (float) readout->image->numCols;
    521         Ny = (int) (Nx * AR + 0.5);
    522         Ny = PS_MAX (1, Ny);
    523     } else {
    524         Ny = scale;
    525         float AR = readout->image->numRows / (float) readout->image->numCols;
    526         Nx = (int) (Ny * AR + 0.5);
    527         Nx = PS_MAX (1, Nx);
    528     }
    529 
    530     // require at least 10 stars per spatial bin
    531     if (Npsf < 10*Nx*Ny) {
    532         return false;
    533     }
    534 
    535     // the mask marks the values not used to calculate the ApTrend
    536     psVectorInit (mask, 0);
    537 
    538     // XXX stats structure for use by ApTrend : make parameters user setable
    539     psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    540     stats->min = 2.0;
    541     stats->max = 3.0;
    542 
    543     // measure Trend2D for the current spatial scale
    544     psFree (psf->ApTrend);
    545     psf->ApTrend = pmTrend2DAlloc (PM_TREND_MAP, readout->image, Nx, Ny, stats);
    546 
    547     // XXX somewhat arbitrary: soften the errors so the few bright stars do not totally dominate:
    548     psVector *dMagSoft = psVectorAlloc (dMag->n, PS_TYPE_F32);
    549     for (int i = 0; i < dMag->n; i++) {
    550         dMagSoft->data.F32[i] = hypot(dMag->data.F32[i], 0.01);
    551     }
    552 
    553     // XXX test for errors here
    554     pmTrend2DFit (psf->ApTrend, mask, 0xff, xPos, yPos, apResid, dMagSoft);
    555 
    556     // construct the fitted values and the residuals
    557     psVector *apResidFit = pmTrend2DEvalVector (psf->ApTrend, mask, 0xff, xPos, yPos);
    558     psVector *apResidRes = (psVector *) psBinaryOp (NULL, (void *) apResid, "-", (void *) apResidFit);
    559 
    560     // measure systematic errorFloor & systematic / photon scale factor
    561     // XXX this is a bit arbitrary, but it forces ~3 stars from the bright bin per spatial bin
    562     int nGroup = PS_MAX (3*Nx*Ny, 10);
    563     psphotMagErrorScale (errorScale, errorFloor, dMag, apResidRes, mask, nGroup);
    564 
    565     psLogMsg ("psphot.apresid", PS_LOG_INFO, "result of %d x %d grid (%d stars per bin)\n", Nx, Ny, nGroup);
    566     psLogMsg ("psphot.apresid", PS_LOG_INFO, "systematic error / photon error: %f\n", *errorScale);
    567     psLogMsg ("psphot.apresid", PS_LOG_INFO, "systematic scatter floor: %f\n", *errorFloor);
    568 
    569     psFree (stats);
    570     psFree (dMagSoft);
    571     psFree (apResidFit);
    572     psFree (apResidRes);
    573 
    574     return true;
    575 }
    576 
Note: See TracChangeset for help on using the changeset viewer.