IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Nov 17, 2011, 5:41:38 AM (15 years ago)
Author:
eugene
Message:

psphotStack is generating duplicate detections because the detection removal before the second pass only occurs on the SRC, not DET image; working on fixing this by modifying the removal process

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20111110/psphot/src/psphotReplaceUnfit.c

    r31452 r32685  
    11# include "psphotInternal.h"
     2
     3static int replace_pass = 0;
     4static int remove_pass = 0;
    25
    36// replace the flux for sources which failed
    … …  
    912
    1013    for (int i = 0; i < sources->n; i++) {
    11       source = sources->data[i];
    12 
    13       // replace other sources?
    14       if (source->mode & PM_SOURCE_MODE_FAIL) goto replace;
    15       continue;
     14        source = sources->data[i];
     15
     16        // replace other sources?
     17        if (source->mode & PM_SOURCE_MODE_FAIL) goto replace;
     18        continue;
    1619
    1720    replace:
    … …  
    2326
    2427// for now, let's store the detections on the readout->analysis for each readout
    25 bool psphotReplaceAllSources (pmConfig *config, const pmFPAview *view, const char *filerule)
     28bool psphotReplaceAllSources (pmConfig *config, const pmFPAview *view, const char *filerule, bool ignoreState)
    2629{
    2730    bool status = true;
    … …  
    3538    // loop over the available readouts
    3639    for (int i = 0; i < num; i++) {
    37         if (!psphotReplaceAllSourcesReadout (config, view, filerule, i, recipe)) {
     40        if (!psphotReplaceAllSourcesReadout (config, view, filerule, i, recipe, ignoreState)) {
    3841            psError (PSPHOT_ERR_CONFIG, false, "failed to replace all sources for %s entry %d", filerule, i);
    3942            return false;
    4043        }
    4144    }
    42     return true;
    43 }
    44 
    45 bool psphotReplaceAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) {
     45    replace_pass ++;
     46    return true;
     47}
     48
     49bool psphotReplaceAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool ignoreState) {
    4650
    4751    bool status;
    48     pmSource *source;
    4952
    5053    psTimerStart ("psphot.replace");
    … …  
    5760    psAssert (readout, "missing readout?");
    5861
     62    char name[128];
     63    snprintf (name, 128, "testadd.neg.%02d.%02d.fits", replace_pass, index);
     64    psphotSaveImage (NULL, readout->image, name);
     65
    5966    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
    6067    psAssert (detections, "missing detections?");
    … …  
    7582
    7683    for (int i = 0; i < sources->n; i++) {
    77       source = sources->data[i];
    78 
    79       // replace other sources?
    80       if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) continue;
    81 
    82       pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    83     }
     84        pmSource *source = sources->data[i];
     85
     86        if (ignoreState) {
     87            // rely on the type of source to decide if we subtract it or not
     88
     89            // skip non-astronomical objects (very likely defects)
     90            if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
     91            if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
     92     
     93            // do not include CRs in the full ensemble fit
     94            if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue;
     95       
     96            // do not include MOMENTS_FAILURES in the fit
     97            if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
     98        } else {
     99            // if we respect the state, only remove unsubtracted sources
     100            if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) continue;
     101        }
     102
     103        pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
     104    }
     105
     106    snprintf (name, 128, "testadd.pos.%02d.%02d.fits", replace_pass, index);
     107    psphotSaveImage (NULL, readout->image, name);
    84108
    85109    psphotVisualShowImage(readout);
    … …  
    88112}
    89113
    90 bool psphotRemoveAllSources (const psArray *sources, const psMetadata *recipe) {
     114// for now, let's store the detections on the readout->analysis for each readout
     115bool psphotRemoveAllSources (pmConfig *config, const pmFPAview *view, const char *filerule, bool ignoreState)
     116{
     117    bool status = true;
     118
     119    // select the appropriate recipe information
     120    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     121    psAssert (recipe, "missing recipe?");
     122
     123    int num = psphotFileruleCount(config, filerule);
     124
     125    // loop over the available readouts
     126    for (int i = 0; i < num; i++) {
     127        if (!psphotRemoveAllSourcesReadout (config, view, filerule, i, recipe, ignoreState)) {
     128            psError (PSPHOT_ERR_CONFIG, false, "failed to replace all sources for %s entry %d", filerule, i);
     129            return false;
     130        }
     131    }
     132    remove_pass ++;
     133    return true;
     134}
     135
     136bool psphotRemoveAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool ignoreState) {
     137
     138    bool status;
     139
     140    psTimerStart ("psphot.replace");
     141
     142    // find the currently selected readout
     143    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     144    psAssert (file, "missing file?");
     145
     146    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     147    psAssert (readout, "missing readout?");
     148
     149    char name[128];
     150    snprintf (name, 128, "testsub.pos.%02d.%02d.fits", remove_pass, index);
     151    psphotSaveImage (NULL, readout->image, name);
     152
     153    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     154    psAssert (detections, "missing detections?");
     155
     156    psArray *sources = detections->allSources;
     157    psAssert (sources, "missing sources?");
     158
     159    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
     160    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
     161    psAssert (maskVal, "missing mask value?");
     162
     163    // bit-mask to mark pixels not used in analysis
     164    psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT");
     165    assert (markVal);
     166
     167    // maskVal is used to test for rejected pixels, and must include markVal
     168    maskVal |= markVal;
     169
     170    for (int i = 0; i < sources->n; i++) {
     171        pmSource *source = sources->data[i];
     172
     173        // skip non-astronomical objects (very likely defects)
     174        if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
     175        if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
     176       
     177        // do not include CRs in the full ensemble fit
     178        if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue;
     179       
     180        // do not include MOMENTS_FAILURES in the fit
     181        if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
     182       
     183        // if we respect the state, only remove unsubtracted sources
     184        if (!ignoreState && (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) continue;
     185       
     186        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     187    }
     188
     189    snprintf (name, 128, "testsub.neg.%02d.%02d.fits", remove_pass, index);
     190    psphotSaveImage (NULL, readout->image, name);
     191
     192    psphotVisualShowImage(readout);
     193    psLogMsg ("psphot.replace", PS_LOG_INFO, "replaced models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.replace"));
     194    return true;
     195}
     196
     197bool psphotRemoveAllSourcesByArray (const psArray *sources, const psMetadata *recipe) {
    91198
    92199    bool status;
    … …  
    100207
    101208    for (int i = 0; i < sources->n; i++) {
    102       source = sources->data[i];
    103 
    104       // replace other sources?
    105       if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) continue;
    106 
    107       pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     209        source = sources->data[i];
     210
     211        // replace other sources?
     212        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) continue;
     213
     214        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    108215    }
    109216    psLogMsg ("psphot.replace", PS_LOG_INFO, "replaced models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.replace"));
    … …  
    158265
    159266    for (int i = 0; i < sources->n; i++) {
    160       source = sources->data[i];
    161 
    162       // sources have not yet been subtracted in this image (but this flag may be raised)
    163       source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
    164       if (!source->modelPSF) continue;
    165 
    166       float Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS];
    167       float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS];
    168       float radius = source->modelPSF->fitRadius;
    169 
    170       // force a redefine to this image
    171       pmSourceFreePixels(source);
    172       pmSourceRedefinePixels (source, readout, Xo, Yo, radius);
     267        source = sources->data[i];
     268
     269        // sources have not yet been subtracted in this image (but this flag may be raised)
     270        source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
     271        if (!source->modelPSF) continue;
     272
     273        float Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS];
     274        float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS];
     275        float radius = source->modelPSF->fitRadius;
     276
     277        // force a redefine to this image
     278        pmSourceFreePixels(source);
     279        pmSourceRedefinePixels (source, readout, Xo, Yo, radius);
    173280    }
    174281    return true;
    … …  
    229336
    230337    for (int i = 0; i < sources->n; i++) {
    231       source = sources->data[i];
    232 
    233       // *** we need to cache the 'best' model, and we have 3 cases:
    234       // 1) model is the psf model --> generate from the new psf
    235       // 2) model is an unconvolved extended model --> just cache the copy (not perfect)
    236       // 3) model is a convolved extended model --> re-generate
    237 
    238       // use the 'best' model to cache the model (PSF or EXT : EXT may point at one of modelFits
    239       bool isPSF = false;
    240       pmModel *model = pmSourceGetModel(&isPSF, source);
    241       if (!model) continue;
    242 
    243       float radius = model->fitRadius; // save for future use below
    244 
    245       // regenerate the PSF if the model is a PSF, or if we need the PSF for a PCM
    246       if (isPSF || model->isPCM) {
    247           // the guess central intensity comes from the peak:
    248           float Io = source->peak->rawFlux;
    249           float Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS];
    250           float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS];
    251 
    252           // generate a model for this object with Io = 1.0
    253           pmModel *modelPSF = pmModelFromPSFforXY(psf, Xo, Yo, Io);
    254           if (modelPSF == NULL) {
    255               psWarning ("Failed to determine PSF model for source at (%f,%f), skipping", Xo, Yo);
    256               continue;
    257           }
    258 
    259           // set the source PSF model
    260           psFree (source->modelPSF);
    261           source->modelPSF = modelPSF;
    262           source->modelPSF->fitRadius = radius;
    263       }
    264 
    265       if (model->isPCM) {
    266           psAssert(false, "this section is not complete");
    267 
    268           pmSourceCachePSF (source, maskVal);
    269 
    270           psKernel *psfKernel = pmPCMkernelFromPSF(source, psfSize);
    271           if (!psfKernel) {
    272               psWarning ("no psf kernel");
    273           }
    274 
    275           // generate an image of the right size
    276           psImage *rawModelFlux = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
    277           psImageInit (rawModelFlux, 0.0);
     338        source = sources->data[i];
     339
     340        // *** we need to cache the 'best' model, and we have 3 cases:
     341        // 1) model is the psf model --> generate from the new psf
     342        // 2) model is an unconvolved extended model --> just cache the copy (not perfect)
     343        // 3) model is a convolved extended model --> re-generate
     344
     345        // use the 'best' model to cache the model (PSF or EXT : EXT may point at one of modelFits
     346        bool isPSF = false;
     347        pmModel *model = pmSourceGetModel(&isPSF, source);
     348        if (!model) continue;
     349
     350        float radius = model->fitRadius; // save for future use below
     351
     352        // regenerate the PSF if the model is a PSF, or if we need the PSF for a PCM
     353        if (isPSF || model->isPCM) {
     354            // the guess central intensity comes from the peak:
     355            float Io = source->peak->rawFlux;
     356            float Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS];
     357            float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS];
     358
     359            // generate a model for this object with Io = 1.0
     360            pmModel *modelPSF = pmModelFromPSFforXY(psf, Xo, Yo, Io);
     361            if (modelPSF == NULL) {
     362                psWarning ("Failed to determine PSF model for source at (%f,%f), skipping", Xo, Yo);
     363                continue;
     364            }
     365
     366            // set the source PSF model
     367            psFree (source->modelPSF);
     368            source->modelPSF = modelPSF;
     369            source->modelPSF->fitRadius = radius;
     370        }
     371
     372        if (model->isPCM) {
     373            psAssert(false, "this section is not complete");
     374
     375            pmSourceCachePSF (source, maskVal);
     376
     377            psKernel *psfKernel = pmPCMkernelFromPSF(source, psfSize);
     378            if (!psfKernel) {
     379                psWarning ("no psf kernel");
     380            }
     381
     382            // generate an image of the right size
     383            psImage *rawModelFlux = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
     384            psImageInit (rawModelFlux, 0.0);
    278385         
    279           // insert the model image normalized to 1.0
    280           pmModelAdd (rawModelFlux, NULL, model, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM, maskVal);
    281 
    282           psImageConvolveFFT (source->modelFlux, rawModelFlux, NULL, 0, psfKernel);
     386            // insert the model image normalized to 1.0
     387            pmModelAdd (rawModelFlux, NULL, model, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM, maskVal);
     388
     389            psImageConvolveFFT (source->modelFlux, rawModelFlux, NULL, 0, psfKernel);
    283390         
    284           psFree (psfKernel);
    285           psFree (rawModelFlux);
    286       }
    287 
    288       pmSourceCacheModel (source, maskVal);  // ALLOC x14 (!)
     391            psFree (psfKernel);
     392            psFree (rawModelFlux);
     393        }
     394
     395        pmSourceCacheModel (source, maskVal);  // ALLOC x14 (!)
    289396    }
    290397
Note: See TracChangeset for help on using the changeset viewer.