IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 2, 2009, 3:12:47 PM (17 years ago)
Author:
eugene
Message:

check in changes from gene@development branch : extensive changes to moments calculation, psf model generation, aperture residuals; improvements to extended source analysis

Location:
trunk/psphot/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src

    • Property svn:ignore

      •  

        old new  
        1818psphotVersionDefinitions.h
        1919psphotMomentsStudy
         20psphotPetrosianStudy

  • trunk/psphot/src/psphotApResid.c

    r24878 r25755  
    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# ifdef DEBUG   
     148    FILE *f = fopen ("apresid.dat", "w");
     149    psAssert (f, "failed open");
     150# endif
     151
    160152    for (int i = 0; i < sources->n; i++) {
    161153        source = sources->data[i];
    … …  
    170162        if (source->mode &  PM_SOURCE_MODE_EXT_LIMIT) SKIPSTAR ("EXTENDED");
    171163        if (source->mode &  PM_SOURCE_MODE_CR_LIMIT) SKIPSTAR ("COSMIC RAY");
     164        if (source->mode &  PM_SOURCE_MODE_DEFECT) SKIPSTAR ("DEFECT");
    172165           
    173166        if (!isfinite(source->apMag) || !isfinite(source->psfMag)) {
    174167            continue;
    175168        }
     169
     170        // XXX make this user-configurable?
     171        if (source->errMag > 0.01) continue;
    176172
    177173        // aperture residual for this source
    … …  
    185181        }
    186182
    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         dMag->data.F32[Npsf] = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
    195 
    196         psVectorExtend (mag,     100, 1);
    197         psVectorExtend (mask,    100, 1);
    198         psVectorExtend (xPos,    100, 1);
    199         psVectorExtend (yPos,    100, 1);
    200         psVectorExtend (dMag,    100, 1);
    201         psVectorExtend (apResid, 100, 1);
     183# ifdef DEBUG
     184        fprintf (f, "%6.1f %6.1f : %6.1f %6.1f : %8.3f %8.3f %8.3f : %f : %f %f %f : %f\n",
     185                 source->peak->xf, source->peak->yf,
     186                 source->modelPSF->params->data.F32[PM_PAR_XPOS], source->modelPSF->params->data.F32[PM_PAR_YPOS],
     187                 source->psfMag, source->apMag, source->errMag,
     188                 source->modelPSF->params->data.F32[PM_PAR_I0],
     189                 source->modelPSF->params->data.F32[PM_PAR_SXX], source->modelPSF->params->data.F32[PM_PAR_SXY], source->modelPSF->params->data.F32[PM_PAR_SYY],
     190                 source->modelPSF->params->data.F32[PM_PAR_7]);
     191# endif
     192        if (!isfinite(source->psfMag)) psAbort ("nan in psfMag");
     193        if (!isfinite(source->errMag)) psAbort ("nan in errMag");
     194        if (!isfinite(source->apMag)) psAbort ("nan in apMag");
     195        if (!isfinite(model->params->data.F32[PM_PAR_XPOS])) psAbort ("nan in xPos");
     196        if (!isfinite(model->params->data.F32[PM_PAR_YPOS])) psAbort ("nan in yPos");
     197
     198        psVectorAppend (mag, source->psfMag);
     199        psVectorAppend (dMag,source->errMag);
     200        psVectorAppend (apResid, dap);
     201        psVectorAppend (xPos, model->params->data.F32[PM_PAR_XPOS]);
     202        psVectorAppend (yPos, model->params->data.F32[PM_PAR_YPOS]);
    202203        Npsf ++;
    203204    }
    … …  
    205206    psLogMsg ("psphot.apresid", PS_LOG_DETAIL, "measure aperture residuals for %d objects (%d skipped, %d failed, %ld invalid)\n",
    206207              Npsf, Nskip, Nfail, sources->n - Npsf - Nskip - Nfail);
     208
     209# ifdef DEBUG
     210    fclose (f);
     211# endif
    207212
    208213    // XXX choose a better value here?
    … …  
    212217    }
    213218
    214     // XXX deprecating the old code which allowed the ApResid to be fitted as a function of flux and r^2/flux
    215     // XXX is this asymmetric clipping still needed?  this analysis should come after neighbors are subtracted...
    216     // 3hi/1lo sigma clipping on the rflux vs metric fit
    217     // systematic error information
    218     float errorScale = 0.0;
    219     float errorFloor = 0.0;
    220 
     219    // XXX set the min number of needed source more carefully
     220    if ((Npsf < 15) && (APTREND_ORDER_MAX >= 4)) APTREND_ORDER_MAX = 3;
     221    if ((Npsf < 11) && (APTREND_ORDER_MAX >= 3)) APTREND_ORDER_MAX = 2;
     222    if ((Npsf <  8) && (APTREND_ORDER_MAX >= 2)) APTREND_ORDER_MAX = 1;
     223
     224    psFree (psf->ApTrend);
     225    psf->ApTrend = NULL;
    221226    float errorFloorMin = FLT_MAX;
    222     int entryMin = -1;
    223 
    224     // Fit out the dap vs mag trend, iterate over spatial scale until error Floor increases.
    225     // Stop if Npsf / (Nx * Ny) < 3
     227
     228    // as we loop over orders, we need to refer to the initial selection, but we modify the
     229    // option values to match the current guess: save the max values here:
     230    int NX = readout->image->numCols;
     231    int NY = readout->image->numRows;
    226232    for (int i = 1; i <= APTREND_ORDER_MAX; i++) {
    227 
    228         if (!psphotApResidTrend (readout, psf, Npsf, i, &errorScale, &errorFloor, mask, xPos, yPos, apResid, dMag)) {
    229             break;
    230         }
    231 
    232         // store the resulting errorFloor values and the scales, redo the best
     233       
     234        int Nx, Ny;
     235        if (NX > NY) {
     236            Nx = i;
     237            Ny = PS_MAX (1, (int)(i * (NY / (float)(NX)) + 0.5));
     238        } else {
     239            Ny = i;
     240            Nx = PS_MAX (1, (int)(i * (NX / (float)(NY)) + 0.5));
     241        }
     242
     243        float errorFloor;
     244        pmTrend2D *apTrend = psphotApResidTrend (&errorFloor, readout, Nx, Ny, xPos, yPos, apResid, dMag);
     245        if (!apTrend) {
     246            continue;
     247        }
     248
     249        // apply ApTrend results
     250        // float xc = 0.5*readout->image->numCols + readout->image->col0 + 0.5;
     251        // float yc = 0.5*readout->image->numRows + readout->image->row0 + 0.5;
     252        // float ApResid = pmTrend2DEval (psf->ApTrend, xc, yc); // ap-fit at chip center
     253        // if (!isfinite(ApResid)) psAbort("nan apresid @ center");
     254
     255        // store the minimum errorFloor and best ApTrend to keep
    233256        if (errorFloor < errorFloorMin) {
    234257            errorFloorMin = errorFloor;
    235             entryMin = i;
     258            psFree (psf->ApTrend);
     259            psf->ApTrend = psMemIncrRefCounter(apTrend);
    236260        }
    237     }
    238     if (entryMin == -1) {
     261        psFree (apTrend);
     262    }
     263    if (psf->ApTrend == NULL) {
    239264        psWarning("Failed to find a valid aperture residual value");
    240265        goto escape;
    241266    }
    242267
    243     // XXX catch error condition
    244     psphotApResidTrend (readout, psf, Npsf, entryMin, &errorScale, &errorFloor, mask, xPos, yPos, apResid, dMag);
     268    // apply ApTrend results
     269    float xc = 0.5*readout->image->numCols + readout->image->col0 + 0.5;
     270    float yc = 0.5*readout->image->numRows + readout->image->row0 + 0.5;
     271
     272    psf->ApResid  = pmTrend2DEval (psf->ApTrend, xc, yc); // ap-fit at chip center
     273    psf->dApResid = errorFloorMin;
     274    psf->nApResid = Npsf;
     275
     276    // save results for later output
     277    psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   psf->ApResid);
     278    psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->dApResid);
     279    psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", psf->nApResid);
     280    psMetadataAdd (recipe, PS_LIST_TAIL, "APLOSS",   PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", psf->growth->apLoss);
     281
     282    psLogMsg ("psphot.apresid", PS_LOG_DETAIL, "aperture residual: %f +/- %f\n", psf->ApResid, psf->dApResid);
     283    psLogMsg ("psphot.apresid", PS_LOG_INFO, "measure full-frame aperture residuals for %d sources: %f sec\n", Npsf, psTimerMark ("psphot.apresid"));
     284
     285    psFree (xPos);
     286    psFree (yPos);
     287    psFree (apResid);
     288    psFree (mag);
     289    psFree (dMag);
     290
     291    psphotVisualPlotApResid (sources, psf->ApResid, psf->dApResid);
     292
     293    return true;
     294
     295escape:
     296    // save nan values since these were not calculated
     297    psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   NAN);
     298    psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", NAN);
     299    psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", 0);
     300    psMetadataAdd (recipe, PS_LIST_TAIL, "APLOSS",   PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", NAN);
     301
     302    psFree (xPos);
     303    psFree (yPos);
     304    psFree (apResid);
     305    psFree (mag);
     306    psFree (dMag);
     307    return false;
     308}
     309
     310pmTrend2D *psphotApResidTrend (float *apResidSysErr, pmReadout *readout, int Nx, int Ny, psVector *xPos, psVector *yPos, psVector *apResid, psVector *dMag) {
     311
     312    // the mask marks the values not used to calculate the ApTrend
     313    psVector *mask = psVectorAlloc(xPos->n, PS_TYPE_VECTOR_MASK);
     314    psVectorInit (mask, 0);
     315
     316    // XXX allow user to set this optionally?
     317    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
     318
     319    // measure Trend2D for the current spatial scale
     320    pmTrend2D *apTrend = pmTrend2DAlloc (PM_TREND_MAP, readout->image, Nx, Ny, stats);
     321
     322    // XXX somewhat arbitrary: soften the errors so the few bright stars do not totally dominate:
     323    // XXX use this or not?  probably not, since this is the point of the systematic error analysis
     324    psVector *dMagSoft = psVectorAlloc (dMag->n, PS_TYPE_F32);
     325    for (int i = 0; i < dMag->n; i++) {
     326        dMagSoft->data.F32[i] = hypot(dMag->data.F32[i], 0.005);
     327    }
     328
     329    // XXX test for errors here
     330    if (!pmTrend2DFit (apTrend, mask, 0xff, xPos, yPos, apResid, dMagSoft)) {
     331        psWarning("Failed to fit trend for %d x %d map", Nx, Ny);
     332        psFree (apTrend);
     333        return NULL;
     334    }
     335    if (apTrend->mode == PM_TREND_MAP) {
     336        // p_psImagePrint (2, apTrend->map->map, "ApTrend Before"); // XXX TEST:
     337        psImageMapRepair (apTrend->map->map);
     338        // p_psImagePrint (2, apTrend->map->map, "ApTrend After"); // XXX TEST:
     339    }
    245340
    246341    // construct the fitted values and the residuals
    247     psVector *apResidFit = pmTrend2DEvalVector (psf->ApTrend, mask, 0xff, xPos, yPos);
     342    psVector *apResidFit = pmTrend2DEvalVector (apTrend, mask, 0xff, xPos, yPos);
    248343    psVector *apResidRes = (psVector *) psBinaryOp (NULL, (void *) apResid, "-", (void *) apResidFit);
    249     psVector *dMagSys = (psVector *) psBinaryOp (NULL, (void *) dMag, "*", (void *) psScalarAlloc(errorScale, PS_TYPE_F32));
    250 
    251     if (psTraceGetLevel("psphot") >= 2) {
    252         FILE *dumpFile = fopen ("apresid.dat", "w");
     344
     345    // measure systematic error
     346    *apResidSysErr = psVectorSystematicError (apResidRes, dMag, 0.10);
     347    if (!isfinite(*apResidSysErr)) {
     348        psWarning("Failed to find systematic error for %d x %d map", Nx, Ny);
     349        psFree (apTrend);
     350        return NULL;
     351    }
     352
     353    psLogMsg ("psphot.apresid", PS_LOG_INFO, "result of %d x %d grid\n", Nx, Ny);
     354    psLogMsg ("psphot.apresid", PS_LOG_INFO, "systematic scatter floor: %f\n", *apResidSysErr);
     355
     356    if (psTraceGetLevel("psphot") >= 4) {
     357        char filename[64];
     358        snprintf (filename, 64, "apresid.%dx%d.dat", Nx, Ny);
     359        FILE *dumpFile = fopen (filename, "w");
    253360        for (int i = 0; i < xPos->n; i++) {
    254             fprintf (dumpFile, "%f %f  %f %f %f %f %f  %x\n",
     361            fprintf (dumpFile, "%f %f  %f %f %f %f  %x\n",
    255362                     xPos->data.F32[i], yPos->data.F32[i],
    256                      mag->data.F32[i], dMag->data.F32[i], dMagSys->data.F32[i],
     363                     dMag->data.F32[i], hypot(dMag->data.F32[i], *apResidSysErr),
    257364                     apResid->data.F32[i], apResidRes->data.F32[i],
    258365                     mask->data.PS_TYPE_VECTOR_MASK_DATA[i]);
    … …  
    261368    }
    262369
    263     // apply ApTrend results
    264     float xc = 0.5*readout->image->numCols + readout->image->col0 + 0.5;
    265     float yc = 0.5*readout->image->numRows + readout->image->row0 + 0.5;
    266 
    267     psf->ApResid  = pmTrend2DEval (psf->ApTrend, xc, yc); // ap-fit at chip center
    268     psf->dApResid = errorFloor;
    269     psf->nApResid = Npsf;
    270 
    271     // save results for later output
    272     psMetadataAdd (recipe, PS_LIST_TAIL, "SKYBIAS",  PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   0.0);
    273     psMetadataAdd (recipe, PS_LIST_TAIL, "SKYSAT",   PS_DATA_F32 | PS_META_REPLACE, "aperture-determined saturation",   0.0);
    274     psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   psf->ApResid);
    275     psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->dApResid);
    276     psMetadataAdd (recipe, PS_LIST_TAIL, "APLOSS",   PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", psf->growth->apLoss);
    277     psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", psf->nApResid);
    278 
    279     psLogMsg ("psphot.apresid", PS_LOG_DETAIL, "aperture residual: %f +/- %f\n", psf->ApResid, psf->dApResid);
    280     psLogMsg ("psphot.apresid", PS_LOG_INFO, "measure full-frame aperture residuals for %d sources: %f sec\n", Npsf, psTimerMark ("psphot.apresid"));
    281 
    282     psFree (mag);
    283370    psFree (mask);
    284     psFree (xPos);
    285     psFree (yPos);
    286 
    287     psFree (apResid);
    288     psFree (apResidFit);
    289     psFree (apResidRes);
    290 
    291     psFree (dMagSys);
    292     psFree (dMag);
    293 
    294     psphotVisualPlotApResid (sources);
    295 
    296     return true;
    297 
    298 escape:
    299     // save nan values since these were not calculated
    300     psMetadataAdd (recipe, PS_LIST_TAIL, "SKYBIAS",  PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   NAN);
    301     psMetadataAdd (recipe, PS_LIST_TAIL, "SKYSAT",   PS_DATA_F32 | PS_META_REPLACE, "aperture-determined saturation",   NAN);
    302     psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   NAN);
    303     psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", NAN);
    304     psMetadataAdd (recipe, PS_LIST_TAIL, "APLOSS",   PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", NAN);
    305     psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", 0);
    306 
    307     psFree (mag);
    308     psFree (mask);
    309     psFree (xPos);
    310     psFree (yPos);
    311     psFree (apResid);
    312     psFree (dMag);
    313     return false;
    314 }
    315 
    316 /*
    317   (aprMag' - fitMag) = rflux*skyBias + ApTrend(x,y)
    318   (aprMag - rflux*skyBias) - fitMag = ApTrend(x,y)
    319   (aprMag - rflux*skyBias) = fitMag + ApTrend(x,y)
    320 */
    321 
    322 // XXX this still sucks...  need a better way to estimate the error floor...
    323 bool psphotMagErrorScale (float *errorScale, float *errorFloor, psVector *dMag, psVector *dap, psVector *mask, int nGroup) {
    324 
    325     psStats *statsS = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    326     psStats *statsM = psStatsAlloc (PS_STAT_SAMPLE_MEAN);
    327 
    328     // measure the trend in bins with 10 values each; if < 10 total, use them all
    329     int nBin = PS_MAX (dMag->n / nGroup, 1);
    330 
    331     // output vectors for ApResid trend
    332     psVector *dSo = psVectorAlloc (nBin, PS_TYPE_F32);
    333     psVector *dMo = psVectorAlloc (nBin, PS_TYPE_F32);
    334     psVector *dRo = psVectorAlloc (nBin, PS_TYPE_F32);
    335 
    336     // use dMag to group the dMag and dap vectors
    337     psVector *index = psVectorSortIndex (NULL, dMag);
    338 
    339     // subset vectors for dMag and dap values within the given range
    340     psVector *dMSubset = psVectorAllocEmpty (nGroup, PS_TYPE_F32);
    341     psVector *dASubset = psVectorAllocEmpty (nGroup, PS_TYPE_F32);
    342     psVector *mkSubset = psVectorAllocEmpty (nGroup, PS_TYPE_VECTOR_MASK);
    343 
    344     int n = 0;
    345     for (int i = 0; i < dMo->n; i++) {
    346         int j;
    347         for (j = 0; (j < nGroup) && (n < dMag->n); j++, n++) {
    348             int N = index->data.U32[n];
    349             dMSubset->data.F32[j] = dMag->data.F32[N];
    350             dASubset->data.F32[j] = dap->data.F32[N];
    351             mkSubset->data.PS_TYPE_VECTOR_MASK_DATA[j] = mask->data.PS_TYPE_VECTOR_MASK_DATA[N];
    352         }
    353         dMSubset->n = j;
    354         dASubset->n = j;
    355         mkSubset->n = j;
    356 
    357         psStatsInit (statsS);
    358         psStatsInit (statsM);
    359 
    360         if (j > 2) {
    361             if (!psVectorStats (statsS, dASubset, NULL, mkSubset, 0xff)) {
    362                 psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
    363                 return false;
    364             }
    365             if (!psVectorStats (statsM, dMSubset, NULL, mkSubset, 0xff)) {
    366                 psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
    367                 return false;
    368             }
    369             dSo->data.F32[i] = statsS->robustStdev;
    370             dMo->data.F32[i] = statsM->sampleMean;
    371             dRo->data.F32[i] = statsS->robustStdev / statsM->sampleMean;
    372         } else {
    373             dSo->data.F32[i] = NAN;
    374             dMo->data.F32[i] = NAN;
    375             dRo->data.F32[i] = NAN;
    376         }
    377     }
    378     psFree (dMSubset);
    379     psFree (dASubset);
    380     psFree (mkSubset);
    381 
    382     psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
    383     if (!psVectorStats (stats, dRo, NULL, NULL, 0)) {
    384         psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
    385         return false;
    386     }
    387 
    388     *errorScale = stats->sampleMedian;
    389     for (int i = 0; i < dSo->n; i++) {
    390         *errorFloor = dSo->data.F32[i];
    391         if (fabs(*errorFloor) <= FLT_EPSILON) continue;
    392         if (isfinite(*errorFloor)) break;
    393     }
    394 
    395     psFree (stats);
    396     psFree (index);
    397 
    398     psFree (dRo);
    399     psFree (dMo);
    400     psFree (dSo);
    401 
    402     psFree (statsS);
    403     psFree (statsM);
    404 
    405     return true;
    406 }
    407 
    408 bool psphotApResidTrend (pmReadout *readout, pmPSF *psf, int Npsf, int scale, float *errorScale, float *errorFloor, psVector *mask, psVector *xPos, psVector *yPos, psVector *apResid, psVector *dMag) {
    409 
    410     int Nx, Ny;
    411 
    412     if (readout->image->numCols > readout->image->numRows) {
    413         Nx = scale;
    414         float AR = readout->image->numRows / (float) readout->image->numCols;
    415         Ny = (int) (Nx * AR + 0.5);
    416         Ny = PS_MAX (1, Ny);
    417     } else {
    418         Ny = scale;
    419         float AR = readout->image->numRows / (float) readout->image->numCols;
    420         Nx = (int) (Ny * AR + 0.5);
    421         Nx = PS_MAX (1, Nx);
    422     }
    423 
    424     // require at least 10 stars per spatial bin
    425     if (Npsf < 10*Nx*Ny) {
    426         return false;
    427     }
    428 
    429     // the mask marks the values not used to calculate the ApTrend
    430     psVectorInit (mask, 0);
    431 
    432     // XXX stats structure for use by ApTrend : make parameters user setable
    433     psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    434     stats->min = 2.0;
    435     stats->max = 3.0;
    436 
    437     // measure Trend2D for the current spatial scale
    438     psFree (psf->ApTrend);
    439     psf->ApTrend = pmTrend2DAlloc (PM_TREND_MAP, readout->image, Nx, Ny, stats);
    440 
    441     // XXX somewhat arbitrary: soften the errors so the few bright stars do not totally dominate:
    442     psVector *dMagSoft = psVectorAlloc (dMag->n, PS_TYPE_F32);
    443     for (int i = 0; i < dMag->n; i++) {
    444         dMagSoft->data.F32[i] = hypot(dMag->data.F32[i], 0.01);
    445     }
    446 
    447     // XXX test for errors here
    448     pmTrend2DFit (psf->ApTrend, mask, 0xff, xPos, yPos, apResid, dMagSoft);
    449 
    450     // construct the fitted values and the residuals
    451     psVector *apResidFit = pmTrend2DEvalVector (psf->ApTrend, mask, 0xff, xPos, yPos);
    452     psVector *apResidRes = (psVector *) psBinaryOp (NULL, (void *) apResid, "-", (void *) apResidFit);
    453 
    454     // measure systematic errorFloor & systematic / photon scale factor
    455     // XXX this is a bit arbitrary, but it forces ~3 stars from the bright bin per spatial bin
    456     int nGroup = PS_MAX (3*Nx*Ny, 10);
    457     psphotMagErrorScale (errorScale, errorFloor, dMag, apResidRes, mask, nGroup);
    458 
    459     psLogMsg ("psphot.apresid", PS_LOG_INFO, "result of %d x %d grid (%d stars per bin)\n", Nx, Ny, nGroup);
    460     psLogMsg ("psphot.apresid", PS_LOG_INFO, "systematic error / photon error: %f\n", *errorScale);
    461     psLogMsg ("psphot.apresid", PS_LOG_INFO, "systematic scatter floor: %f\n", *errorFloor);
    462 
    463371    psFree (stats);
    464372    psFree (dMagSoft);
    … …  
    466374    psFree (apResidRes);
    467375
    468     return true;
     376    return apTrend;
    469377}
    470378
    … …  
    480388    pmSourcePhotometryMode photMode = PS_SCALAR_VALUE(job->args->data[2],S32);
    481389    psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[3],PS_TYPE_IMAGE_MASK_DATA);
     390    psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
    482391
    483392    for (int i = 0; i < sources->n; i++) {
    … …  
    490399        if (source->mode &  PM_SOURCE_MODE_POOR) SKIPSTAR ("POOR STAR");
    491400
    492         if (!pmSourceMagnitudes (source, psf, photMode, maskVal)) {
    493             Nskip ++;
    494             psTrace ("psphot", 3, "skip : bad source mag");
    495             continue;
     401        // replace object in image
     402        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
     403            pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    496404        }
    497405
    498         if (!isfinite(source->apMag) || !isfinite(source->psfMag)) {
    499             Nfail ++;
    500             psTrace ("psphot", 3, "fail : nan mags : %f %f", source->apMag, source->psfMag);
    501             continue;
    502         }
    503         source->mode |= PM_SOURCE_MODE_AP_MAGS;
     406        // clear the mask bit and set the circular mask pixels
     407        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
     408        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", markVal);
     409
     410        bool status = pmSourceMagnitudes (source, psf, photMode, maskVal);
     411
     412        // clear the mask bit
     413        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
     414
     415        // re-subtract the object, leave local sky
     416        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     417
     418        if (!status) {
     419            Nskip ++;
     420            psTrace ("psphot", 3, "skip : bad source mag");
     421            continue;
     422        }
     423   
     424        if (!isfinite(source->apMag) || !isfinite(source->psfMag)) {
     425            Nfail ++;
     426            psTrace ("psphot", 3, "fail : nan mags : %f %f", source->apMag, source->psfMag);
     427            continue;
     428        }
     429        source->mode |= PM_SOURCE_MODE_AP_MAGS;
    504430    }
    505431
    506432    // change the value of a scalar on the array (wrap this and put it in psArray.h)
    507     scalar = job->args->data[4];
     433    scalar = job->args->data[5];
    508434    scalar->data.S32 = Nskip;
    509435
    510     scalar = job->args->data[5];
     436    scalar = job->args->data[6];
    511437    scalar->data.S32 = Nfail;
    512438
    513439    return true;
    514440}
    515 
    516 # if (0)
    517 bool psphotApResidMags_Unthreaded (int *nskip, int *nfail, psArray *sources, pmPSF *psf, pmSourcePhotometryMode photMode, psImageMaskType maskVal) {
    518 
    519     int Nskip = 0;
    520     int Nfail = 0;
    521 
    522     for (int i = 0; i < sources->n; i++) {
    523         pmSource *source = (pmSource *) sources->data[i];
    524 
    525         if (source->type != PM_SOURCE_TYPE_STAR) SKIPSTAR ("NOT STAR");
    526         if (source->mode &  PM_SOURCE_MODE_SATSTAR) SKIPSTAR ("SATSTAR");
    527         if (source->mode &  PM_SOURCE_MODE_BLEND) SKIPSTAR ("BLEND");
    528         if (source->mode &  PM_SOURCE_MODE_FAIL) SKIPSTAR ("FAIL STAR");
    529         if (source->mode &  PM_SOURCE_MODE_POOR) SKIPSTAR ("POOR STAR");
    530 
    531         if (!pmSourceMagnitudes (source, psf, photMode, maskVal)) {
    532             Nskip ++;
    533             psTrace ("psphot", 3, "skip : bad source mag");
    534             continue;
    535         }
    536 
    537         if (!isfinite(source->apMag) || !isfinite(source->psfMag)) {
    538             Nfail ++;
    539             psTrace ("psphot", 3, "fail : nan mags : %f %f", source->apMag, source->psfMag);
    540             continue;
    541         }
    542     }
    543 
    544     // change the value of a scalar on the array (wrap this and put it in psArray.h)
    545     *nskip = Nskip;
    546     *nfail = Nfail;
    547 
    548     return true;
    549 }
    550 # endif
Note: See TracChangeset for help on using the changeset viewer.