IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 26681


Ignore:
Timestamp:
Jan 25, 2010, 7:52:41 PM (16 years ago)
Author:
eugene
Message:

various API updates to work with multiple INPUT files (not yet finished)

Location:
branches/eam_branches/psphot.stack.20100120/src
Files:
20 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/psphot.stack.20100120/src/psphot.h

    r26635 r26681  
    7575bool            psphotReplaceUnfitSources (psArray *sources, psImageMaskType maskVal);
    7676
     77bool psphotSetMomentsWindow (psMetadata *recipe, psMetadata *analysis, psArray *sources);
     78
    7779bool            psphotReplaceAllSources (psArray *sources, psMetadata *recipe);
    7880bool            psphotRemoveAllSources (const psArray *sources, const psMetadata *recipe);
    … …  
    126128bool            psphotSetMaskAndVarianceReadout (pmConfig *config, pmReadout *readout, psMetadata *recipe);
    127129void            psphotSourceFreePixels (psArray *sources);
     130
     131bool psphotRoughClassRegion (int nRegion, psRegion *region, psArray *sources, psMetadata *target, psMetadata *recipe, const bool havePSF);
    128132
    129133// functions to set the correct source pixels

  • branches/eam_branches/psphot.stack.20100120/src/psphotAddNoise.c

    r25755 r26681  
    11# include "psphotInternal.h"
    22
    3 bool psphotAddNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe) {
    4   return psphotAddOrSubNoise (readout, sources, recipe, true);
     3bool psphotAddNoise (pmConfig *config, const pmFPAview *view) {
     4    return psphotAddOrSubNoise (config, view, true);
    55}
    66
    7 bool psphotSubNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe) {
    8   return psphotAddOrSubNoise (readout, sources, recipe, false);
     7bool psphotSubNoise (pmConfig *config, const pmFPAview *view) {
     8    return psphotAddOrSubNoise (config, view, false);
    99}
    1010
    11 bool psphotAddOrSubNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe, bool add) {
     11// for now, let's store the detections on the readout->analysis for each readout
     12bool psphotAddOrSubNoise (pmConfig *config, const pmFPAview *view, bool add)
     13{
     14    bool status = true;
     15
     16    // select the appropriate recipe information
     17    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     18    psAssert (recipe, "missing recipe?");
     19
     20    int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     21    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     22
     23    // loop over the available readouts
     24    for (int i = 0; i < num; i++) {
     25        if (!psphotAddOrSubNoiseReadout (config, view, "PSPHOT.INPUT", i, recipe, add)) {
     26            psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for PSPHOT.INPUT entry %d", i);
     27            return false;
     28        }
     29    }
     30    return true;
     31}
     32
     33bool psphotAddOrSubNoise (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool add) {
    1234
    1335    bool status = false;
    … …  
    1638    psEllipseAxes axes;
    1739
    18     PS_ASSERT (readout, false);
    19     PS_ASSERT (readout->parent, false);
    20     PS_ASSERT (readout->parent->concepts, false);
     40    // find the currently selected readout
     41    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     42    psAssert (readout, "missing file?");
     43
     44    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     45    psAssert (readout, "missing readout?");
     46    psAssert (readout->parent, "missing cell?");
     47    psAssert (readout->parent->concepts, "missing concepts?");
     48
     49    psArray *sources = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES");
     50    psAssert (sources, "missing sources?");
    2151
    2252    psTimerStart ("psphot.noise");
    … …  
    2454    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    2555    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
    26     assert (maskVal);
     56    psAssert (maskVal, "missing mask value?");
    2757
    2858    // increase variance by factor*(object noise):

  • branches/eam_branches/psphot.stack.20100120/src/psphotBasicDeblend.c

    r26643 r26681  
    11# include "psphotInternal.h"
     2
     3// for now, let's store the detections on the readout->analysis for each readout
     4bool psphotBasicDeblend (pmConfig *config, const pmFPAview *view)
     5{
     6    bool status = true;
     7
     8    int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     9    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     10
     11    // loop over the available readouts
     12    for (int i = 0; i < num; i++) {
     13        if (!psphotBasicDeblendReadout (config, view, "PSPHOT.INPUT", i)) {
     14            psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for PSPHOT.INPUT entry %d", i);
     15            return false;
     16        }
     17    }
     18    return true;
     19}
    220
    321bool psphotBasicDeblend (pmConfig *config, const pmFPAview *view, const char *filename, int index) {
    … …  
    2139    psArray *sources = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES");
    2240    psAssert (sources, "missing sources?");
     41
     42    if (!sources->n) {
     43        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping basic deblend");
     44        return true;
     45    }
    2346
    2447    // select the appropriate recipe information
    … …  
    142165    return true;
    143166}
    144 
    145 // for now, let's store the detections on the readout->analysis for each readout
    146 bool psphotBasicDeblend (pmConfig *config, const pmFPAview *view)
    147 {
    148     bool status = true;
    149 
    150     int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
    151     psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
    152 
    153     // loop over the available readouts
    154     for (int i = 0; i < num; i++) {
    155         if (!psphotBasicDeblendReadout (config, view, "PSPHOT.INPUT", i)) {
    156             psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for PSPHOT.INPUT entry %d", i);
    157             return false;
    158         }
    159     }
    160     return true;
    161 }

  • branches/eam_branches/psphot.stack.20100120/src/psphotBlendFit.c

    r25755 r26681  
    11# include "psphotInternal.h"
    22
     3// for now, let's store the detections on the readout->analysis for each readout
     4bool psphotBlendFit (pmConfig *config, const pmFPAview *view)
     5{
     6    bool status = true;
     7
     8    // select the appropriate recipe information
     9    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     10    psAssert (recipe, "missing recipe?");
     11
     12    int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     13    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     14
     15    // loop over the available readouts
     16    for (int i = 0; i < num; i++) {
     17        if (!psphotBlendFitReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
     18            psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for PSPHOT.INPUT entry %d", i);
     19            return false;
     20        }
     21    }
     22    return true;
     23}
     24
    325// XXX I don't like this name
    4 bool psphotBlendFit (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf) {
     26bool psphotBlendFitReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) {
    527
    628    int Nfit = 0;
    … …  
    1234    psTimerStart ("psphot.fit.nonlinear");
    1335
    14     // select the appropriate recipe information
    15     psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
    16     assert (recipe);
     36    // find the currently selected readout
     37    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     38    psAssert (readout, "missing file?");
     39
     40    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     41    psAssert (readout, "missing readout?");
     42
     43    psArray *sources = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES");
     44    psAssert (sources, "missing sources?");
     45
     46    if (!sources->n) {
     47        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend fit");
     48        return true;
     49    }
     50
     51    pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
     52    psAssert (psf, "missing psf?");
    1753
    1854    // determine the number of allowed threads
    … …  
    4783    // source analysis is done in S/N order (brightest first)
    4884    sources = psArraySort (sources, pmSourceSortBySN);
     85    if (!sources->n) {
     86        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend");
     87        return true;
     88    }
    4989
    5090    // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)

  • branches/eam_branches/psphot.stack.20100120/src/psphotChoosePSF.c

    r26643 r26681  
    11# include "psphotInternal.h"
    22
     3// generate a PSF model for inputs without PSF models already loaded
     4bool psphotChoosePSF (pmConfig *config, const pmFPAview *view)
     5{
     6    bool status = true;
     7
     8    // select the appropriate recipe information
     9    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     10    psAssert (recipe, "missing recipe?");
     11
     12    int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     13    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     14
     15    // loop over the available readouts
     16    for (int i = 0; i < num; i++) {
     17        if (!psphotChoosePSFReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
     18            psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for PSPHOT.INPUT entry %d", i);
     19            return false;
     20        }
     21    }
     22    return true;
     23}
     24
    325// try PSF models and select best option
    4 bool psphotChoosePSFReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index) {
     26bool psphotChoosePSFReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) {
    527
    628    bool status;
     29
     30    // do not generate a PSF if we already were supplied one
     31    if (psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF")) {
     32        psLogMsg ("psphot", PS_LOG_DETAIL, "psf model supplied for input file %d", i);
     33        return true;
     34    }
    735
    836    psTimerStart ("psphot.choose.psf");
    … …  
    1846    psAssert (sources, "missing sources?");
    1947
    20     // select the appropriate recipe information
    21     psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
    22     psAssert (recipe, "missing recipe?");
     48    if (!sources->n) {
     49        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping PSF model");
     50        return true;
     51    }
    2352
    2453    // bit-masks to test for good/bad pixels
    … …  
    322351        return false;
    323352    }
     353    psFree (psf); // XXX double-check
     354
     355    // move into psphotChoosePSF
     356    psphotVisualShowPSFModel (readout, psf);
    324357
    325358    return true;
    … …  
    466499    return true;
    467500}
    468 
    469 // for now, let's store the detections on the readout->analysis for each readout
    470 bool psphotChoosePSF (pmConfig *config, const pmFPAview *view)
    471 {
    472     bool status = true;
    473 
    474     int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
    475     psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
    476 
    477     // loop over the available readouts
    478     for (int i = 0; i < num; i++) {
    479         if (!psphotChoosePSFReadout (config, view, "PSPHOT.INPUT", i, havePSF)) {
    480             psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for PSPHOT.INPUT entry %d", i);
    481             return false;
    482         }
    483     }
    484     return true;
    485 }

  • branches/eam_branches/psphot.stack.20100120/src/psphotDeblendSatstars.c

    r26643 r26681  
    11# include "psphotInternal.h"
     2
     3// for now, let's store the detections on the readout->analysis for each readout
     4bool psphotDeblendSatstars (pmConfig *config, const pmFPAview *view)
     5{
     6    bool status = true;
     7
     8    int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     9    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     10
     11    // loop over the available readouts
     12    for (int i = 0; i < num; i++) {
     13        if (!psphotDeblendSatstarsReadout (config, view, "PSPHOT.INPUT", i)) {
     14            psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for PSPHOT.INPUT entry %d", i);
     15            return false;
     16        }
     17    }
     18    return true;
     19}
    220
    321bool psphotDeblendSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index) {
    … …  
    2038    psArray *sources = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES");
    2139    psAssert (sources, "missing sources?");
     40
     41    if (!sources->n) {
     42        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping satstar blend");
     43        return true;
     44    }
    2245
    2346    // select the appropriate recipe information
    … …  
    183206    return true;
    184207}
    185 
    186 // for now, let's store the detections on the readout->analysis for each readout
    187 bool psphotDeblendSatstars (pmConfig *config, const pmFPAview *view)
    188 {
    189     bool status = true;
    190 
    191     int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
    192     psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
    193 
    194     // loop over the available readouts
    195     for (int i = 0; i < num; i++) {
    196         if (!psphotDeblendSatstarsReadout (config, view, "PSPHOT.INPUT", i)) {
    197             psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for PSPHOT.INPUT entry %d", i);
    198             return false;
    199         }
    200     }
    201     return true;
    202 }

  • branches/eam_branches/psphot.stack.20100120/src/psphotFitSourcesLinear.c

    r25990 r26681  
    1212static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER, psImageMaskType markVal);
    1313
    14 bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, const psMetadata *recipe, const pmPSF *psf, bool final) {
     14// for now, let's store the detections on the readout->analysis for each readout
     15bool psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, bool final)
     16{
     17    bool status = true;
     18
     19    int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     20    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     21
     22    // loop over the available readouts
     23    for (int i = 0; i < num; i++) {
     24        if (!psphotFitSourcesLinearReadout (config, view, "PSPHOT.INPUT", i, final)) {
     25            psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for PSPHOT.INPUT entry %d", i);
     26            return false;
     27        }
     28    }
     29    return true;
     30}
     31
     32bool psphotFitSourcesLinearReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index) {
    1533
    1634    bool status;
    … …  
    1937    float f;
    2038    // float r;
     39
     40    // select the appropriate recipe information
     41    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     42    assert (recipe);
     43
     44    // find the currently selected readout
     45    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     46    psAssert (readout, "missing file?");
     47
     48    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     49    psAssert (readout, "missing readout?");
     50
     51    psArray *sources = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES");
     52    psAssert (sources, "missing sources?");
     53
     54    if (!sources->n) {
     55        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping linear fit");
     56        return true;
     57    }
     58
     59    pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
     60    psAssert (sources, "missing psf?");
    2161
    2262    psTimerStart ("psphot.linear");
    … …  
    248288    psphotVisualPlotChisq (sources);
    249289    // psphotVisualShowFlags (sources);
     290
     291    // We have to place this visualization here because the models are not realized until
     292    // psphotGuessModels or fitted until psphotFitSourcesLinear.
     293    psphotVisualShowPSFStars (recipe, psf, sources);
    250294
    251295    return true;

  • branches/eam_branches/psphot.stack.20100120/src/psphotGuessModels.c

    r25755 r26681  
    77// 3) define the threaded function to work with sources for a given cell
    88
     9// for now, let's store the detections on the readout->analysis for each readout
     10bool psphotGuessModels (pmConfig *config, const pmFPAview *view)
     11{
     12    bool status = true;
     13
     14    int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     15    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     16
     17    // loop over the available readouts
     18    for (int i = 0; i < num; i++) {
     19        if (!psphotGuessModelsReadout (config, view, "PSPHOT.INPUT", i)) {
     20            psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for PSPHOT.INPUT entry %d", i);
     21            return false;
     22        }
     23    }
     24    return true;
     25}
     26
    927// construct an initial PSF model for each object
    10 bool psphotGuessModels (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf) {
     28bool psphotGuessModelsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index) {
    1129
    1230    bool status;
    1331
    1432    psTimerStart ("psphot.models");
     33
     34    // find the currently selected readout
     35    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     36    psAssert (readout, "missing file?");
     37
     38    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     39    psAssert (readout, "missing readout?");
     40
     41    psArray *sources = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES");
     42    psAssert (sources, "missing sources?");
     43
     44    if (!sources->n) {
     45        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping model guess");
     46        return true;
     47    }
     48
     49    pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
     50    psAssert (sources, "missing psf?");
    1551
    1652    // select the appropriate recipe information

  • branches/eam_branches/psphot.stack.20100120/src/psphotImageQuality.c

    r25755 r26681  
    44// XXX Lots of code duplication here
    55
     6// for now, let's store the detections on the readout->analysis for each readout
     7bool psphotImageQuality (pmConfig *config, const pmFPAview *view)
     8{
     9    bool status = true;
     10
     11    // select the appropriate recipe information
     12    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     13    psAssert (recipe, "missing recipe?");
     14
     15    int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     16    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     17
     18    // loop over the available readouts
     19    for (int i = 0; i < num; i++) {
     20        if (!psphotImageQualityReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
     21            psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for PSPHOT.INPUT entry %d", i);
     22            return false;
     23        }
     24    }
     25    return true;
     26}
     27
    628// selecting the 'good' stars (likely to be psf stars), measure the M_cn, M_sn terms for n = 2,3,4
    7 bool psphotImageQuality(psMetadata *recipe, psArray *sources)
     29bool psphotImageQualityReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe)
    830{
     31    // find the currently selected readout
     32    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     33    psAssert (readout, "missing file?");
     34
     35    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     36    psAssert (readout, "missing readout?");
     37
     38    // if we have a PSF, skip this analysis
     39    if (psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF")) {
     40        return true;
     41    }
     42
     43    psArray *sources = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES");
     44    psAssert (sources, "missing sources?");
     45
     46    if (!sources->n) {
     47        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping image quality");
     48        return true;
     49    }
     50
    951    psVector *FWHM_MAJOR = psVectorAllocEmpty(100, PS_TYPE_F32);
    1052    psVector *FWHM_MINOR = psVectorAllocEmpty(100, PS_TYPE_F32);

  • branches/eam_branches/psphot.stack.20100120/src/psphotLoadPSF.c

    r18002 r26681  
    22
    33// load an externally supplied psf model
    4 pmPSF *psphotLoadPSF (pmConfig *config, const pmFPAview *view, psMetadata *recipe) {
     4pmPSF *psphotLoadPSFReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) {
     5
     6    // find the currently selected readout
     7    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     8    psAssert (readout, "missing file?");
    59
    610    // find the currently selected chip
    7     pmChip *chip = pmFPAfileThisChip (config->files, view, "PSPHOT.PSF.LOAD");
    8     if (!chip) return NULL;
     11    pmChip *chip = pmFPAviewThisChip (view, file->fpa);
     12    if (!chip) return false;
    913
    1014    // find the currently selected readout
    11     pmReadout *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.PSF.LOAD");
    12     if (!readout) return NULL;
     15    pmReadout *readout = pmFPAfileThisReadout (view, file->fpa);
     16    if (!readout) return false;
    1317
    1418    // check if a PSF model is supplied by the user
    … …  
    1620    if (psf == NULL) {
    1721        psLogMsg ("psphot", 3, "no psf supplied for this chip");
    18         return NULL;
     22        return true;
    1923    }
    2024
    21     if (!psphotPSFstats (readout, recipe, psf)) psAbort("cannot measure PSF shape terms");
     25    if (!psphotPSFstats (readout, recipe, psf)) {
     26        psAbort("cannot measure PSF shape terms");
     27    }
    2228
    2329    psLogMsg ("psphot", 3, "using externally supplied PSF model for this readout");
    2430
    25     // we return a psf which can be free (and which may be removed from the metadata)
    26     psMemIncrRefCounter (psf);
    27     return psf;
     31    // save PSF on readout->analysis
     32    if (!psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_PTR, "psphot psf model", psf)) {
     33        psError (PSPHOT_ERR_UNKNOWN, false, "problem saving sources on readout");
     34        return false;
     35    }
     36
     37    return true;
    2838}
     39
     40bool psphotLoadPSF (pmConfig *config, const pmFPAview *view) {
     41
     42    bool status = false;
     43
     44    // select the appropriate recipe information
     45    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     46    psAssert (recipe, "missing recipe?");
     47
     48    // XXX PSPHOT.PSF.LOAD vs PSPHOT.INPUT??
     49    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     50    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     51
     52    // loop over the available readouts
     53    for (int i = 0; i < num; i++) {
     54
     55        // Generate the mask and weight images, including the user-defined analysis region of interest
     56        if (!psphotLoadPSFReadout (config, view, "PSPHOT.PSF.LOAD", i, recipe)) {
     57            psError (PSPHOT_ERR_CONFIG, false, "failed to load PSF model for PSPHOT.PSF.LOAD entry %d", i);
     58            return false;
     59        }
     60    }
     61    return true;
     62}

  • branches/eam_branches/psphot.stack.20100120/src/psphotMaskReadout.c

    r26542 r26681  
    11# include "psphotInternal.h"
    22
     3bool psphotSetMaskAndVariance (pmConfig *config, const pmFPAview *view) {
     4
     5    bool status = false;
     6
     7    // select the appropriate recipe information
     8    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     9    psAssert (recipe, "missing recipe?");
     10
     11    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     12    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     13
     14    // loop over the available readouts
     15    for (int i = 0; i < num; i++) {
     16
     17        // Generate the mask and weight images, including the user-defined analysis region of interest
     18        if (!psphotSetMaskAndVarianceReadout (config, readout, recipe)) {
     19            psError (PSPHOT_ERR_CONFIG, false, "failed to generate mask for PSPHOT.INPUT entry %d", i);
     20            return false;
     21        }
     22
     23        // display the image, weight, mask (ch 1,2,3)
     24        psphotVisualShowImage (readout);
     25    }
     26    return true;
     27}
     28
    329// generate mask and variance if not defined, additional mask for restricted subregion
    4 bool psphotSetMaskAndVarianceReadout (pmConfig *config, pmReadout *readout, psMetadata *recipe) {
     30bool psphotSetMaskAndVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) {
    531
    632    bool status;
     33
     34    pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", index); // File of interest
     35    psAssert (file, "missing file?");
     36
     37    // find the currently selected readout
     38    pmReadout  *readout = pmFPAviewThisReadout (view, file->fpa);
     39    psAssert (readout, "missing readout?");
    740
    841    // ** Interpret the mask values:
    … …  
    73106    }
    74107
     108    // display the image, weight, mask (ch 1,2,3)
     109    psphotVisualShowImage (readout);
     110
    75111    return true;
    76112}
    77 
    78 bool psphotSetMaskAndVariance (pmConfig *config, const pmFPAview *view, psMetadata *recipe) {
    79 
    80     bool status = false;
    81 
    82     int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
    83     psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
    84 
    85     // loop over the available readouts
    86     for (int i = 0; i < num; i++) {
    87 
    88         pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", i); // File of interest
    89         PS_ASSERT_PTR_NON_NULL (file, false);
    90 
    91         // find the currently selected readout
    92         pmReadout  *readout = pmFPAviewThisReadout (view, file->fpa);
    93         PS_ASSERT_PTR_NON_NULL (readout, false);
    94 
    95         // Generate the mask and weight images, including the user-defined analysis region of interest
    96         if (!psphotSetMaskAndVarianceReadout (config, readout, recipe)) {
    97             psError (PSPHOT_ERR_CONFIG, false, "failed to generate mask for PSPHOT.INPUT entry %d", i);
    98             return false;
    99         }
    100 
    101         // display the image, weight, mask (ch 1,2,3)
    102         psphotVisualShowImage (readout);
    103     }
    104     return true;
    105 }

  • branches/eam_branches/psphot.stack.20100120/src/psphotModelBackground.c

    r26542 r26681  
    356356    // find the currently selected readout
    357357    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
    358     PS_ASSERT_PTR_NON_NULL (file, false);
     358    psAssert (file, "missing file?");
    359359
    360360    pmFPA *inFPA = file->fpa;
    361361    pmReadout *readout = pmFPAviewThisReadout(view, inFPA);
     362    psAssert (readout, "missing readout?");
    362363
    363364    psImageBinning *binning = psphotBackgroundBinning(readout->image, config); // Image binning parameters

  • branches/eam_branches/psphot.stack.20100120/src/psphotOutput.c

    r26542 r26681  
    161161}
    162162
    163 bool psphotSetHeaderNstars (psMetadata *recipe, psArray *sources) {
     163bool psphotSetHeaderNstars (psMetadata *header, psArray *sources) {
    164164
    165165    int nSrc = 0;
    … …  
    190190    }
    191191
    192     psMetadataAddS32 (recipe, PS_LIST_TAIL, "NSTARS",   PS_META_REPLACE, "Number of sources", nSrc);
    193     psMetadataAddS32 (recipe, PS_LIST_TAIL, "NDET_EXT", PS_META_REPLACE, "Number of extended sources", nEXT);
    194     psMetadataAddS32 (recipe, PS_LIST_TAIL, "NDET_CR",  PS_META_REPLACE, "Number of cosmic rays", nCR);
    195     return true;
    196 }
    197 
    198 // these values are saved in an output header stub - they are added to either the
    199 // PHU header (CMP) or the MEF table header (CMF)
    200 // XXX these are currently carried by the recipe -- this will not work in a Stack context
    201 // XXX they should be place in the readout->analysis of the given image
     192    psMetadataAddS32 (header, PS_LIST_TAIL, "NSTARS",   PS_META_REPLACE, "Number of sources", nSrc);
     193    psMetadataAddS32 (header, PS_LIST_TAIL, "NDET_EXT", PS_META_REPLACE, "Number of extended sources", nEXT);
     194    psMetadataAddS32 (header, PS_LIST_TAIL, "NDET_CR",  PS_META_REPLACE, "Number of cosmic rays", nCR);
     195    return true;
     196}
     197
     198// these values are saved in an output header stub - they are added to either the PHU header
     199// (CMP) or the MEF table header (CMF) XXX these are currently carried by the recipe -- this
     200// will not work in a Stack context XXX they should be place in the readout->analysis of the
     201// given image or directly on the output header
    202202psMetadata *psphotDefineHeader (psMetadata *recipe) {
    203203

  • branches/eam_branches/psphot.stack.20100120/src/psphotReadout.c

    r26643 r26681  
    33// this should be called by every program that links against libpsphot
    44bool psphotInit (void) {
    5 
    65    psphotErrorRegister();              // register our error codes/messages
    76    psphotModelClassInit ();            // load implementation-specific models
    … …  
    1211bool psphotReadout(pmConfig *config, const pmFPAview *view) {
    1312
    14     // measure the total elapsed time in psphotReadout.  XXX the current threading plan
    15     // for psphot envisions threading within psphotReadout, not multiple threads calling
    16     // the same psphotReadout.  In the current plan, this dtime is the elapsed time used
    17     // jointly by the multiple threads, not the total time used by all threads.
     13    // measure the total elapsed time in psphotReadout.  dtime is the elapsed time used jointly
     14    // by the multiple threads, not the total time used by all threads.
    1815    psTimerStart ("psphotReadout");
    1916
    … …  
    2623        return false;
    2724    }
     25    // optional break-point for processing
     26    char *breakPt = psMetadataLookupStr (NULL, recipe, "BREAK_POINT");
     27    psAssert (breakPt, "configuration error: set BREAK_POINT");
    2828
    2929    // set the photcode for this image
    … …  
    3333    }
    3434
    35     // find the currently selected readout
    36     pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT");
    37     PS_ASSERT_PTR_NON_NULL (readout, false);
    38 
    39     // optional break-point for processing
    40     char *breakPt = psMetadataLookupStr (NULL, recipe, "BREAK_POINT");
    41     PS_ASSERT_PTR_NON_NULL (breakPt, false);
    42 
    4335    // Generate the mask and weight images, including the user-defined analysis region of interest
    44     psphotSetMaskAndVariance (config, view, recipe);
     36    psphotSetMaskAndVariance (config, view);
    4537    if (!strcasecmp (breakPt, "NOTHING")) {
    46         return psphotReadoutCleanup(config, readout, recipe, NULL, NULL, NULL);
    47     }
    48 
    49     // display the image, weight, mask (ch 1,2,3)
    50     // XXX move this into the loop above
    51     psphotVisualShowImage (readout);
     38        return psphotReadoutCleanup(config, view);
     39    }
    5240
    5341    // generate a background model (median, smoothed image)
    5442    if (!psphotModelBackground (config, view)) {
    55         return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL);
     43        return psphotReadoutCleanup (config, view);
    5644    }
    5745    if (!psphotSubtractBackground (config, view)) {
    58         return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL);
     46        return psphotReadoutCleanup (config, view);
    5947    }
    6048    if (!strcasecmp (breakPt, "BACKMDL")) {
    61         return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL);
    62     }
    63 
    64     // display the backsub and backgnd images
    65     // move this inthe the subtract background loop
    66     psphotVisualShowBackground (config, view, readout);
     49        return psphotReadoutCleanup (config, view);
     50    }
    6751
    6852    // run a single-model test if desired (exits from here if test is run)
    69     psphotModelTest (config, view, recipe);
     53    // XXX drop this and only keep the stand-alone program?
     54    // psphotModelTest (config, view, recipe);
    7055
    7156    // load the psf model, if suppled.  FWHM_X,FWHM_Y,etc are saved in the recipe
    72     pmPSF *psf = psphotLoadPSF (config, view, recipe);
    73 
    74     // several functions below behave differently if we have a PSF model already
    75     bool havePSF = (psf != NULL);
    76 
     57    // XXX this needs to save the PSF on the chip->analysis or readout->analysis
     58    if (!psphotLoadPSF (config, view, recipe)) {
     59        // this only happens if we had a programming error in psphotLoadPSF
     60        psError (PSPHOT_ERR_UNKNOWN, false, "error loading psf model");
     61        return psphotReadoutCleanup (config, view);
     62    }
     63       
    7764    // find the detections (by peak and/or footprint) in the image.
    7865    if (!psphotFindDetections (config, view)) {
    7966        // this only happens if we had an error in psphotFindDetections
    8067        psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
    81         return psphotReadoutCleanup (config, readout, recipe, detections, psf, NULL);
    82     }
     68        return psphotReadoutCleanup (config, view);
     69    }
     70
     71XXX: // need to handle or ignore this !!
    8372    if (!detections->peaks->n) {
    8473        psLogMsg ("psphot", 3, "unable to find detections in this image");
    85         return psphotReadoutCleanup (config, readout, recipe, detections, psf, NULL);
     74        return psphotReadoutCleanup (config, view);
    8675    }
    8776
    8877    // construct sources and measure basic stats
    8978    if (!psphotSourceStats (config, view, true)) {
    90         // XXX do I need to raise an error here?
    91         return false;
     79        psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
     80        return psphotReadoutCleanup (config, view);
    9281    }
    9382    if (!strcasecmp (breakPt, "PEAKS")) {
    94         return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources);
     83        return psphotReadoutCleanup(config, view);
    9584    }
    9685
    9786    // find blended neighbors of very saturated stars
    98     // XXX merge this with Basic Deblend?
    9987    if (!psphotDeblendSatstars (config, view)) {
    100         psLogMsg ("psphot", 3, "failed on satstar deblend analysis");
    101         return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
     88        psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis");
     89        return psphotReadoutCleanup (config, view);
    10290    }
    10391
    10492    // mark blended peaks PS_SOURCE_BLEND
    10593    if (!psphotBasicDeblend (config, view)) {
    106         psLogMsg ("psphot", 3, "failed on deblend analysis");
    107         return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
    108     }
    109 
    110     // classify sources based on moments, brightness
    111     if (!psphotRoughClass (config, view, havePSF)) {
    112         psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
    113         return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
     94        psError (PSPHOT_ERR_UNKNOWN, false, "failed on deblend analysis");
     95        return psphotReadoutCleanup (config, view);
     96    }
     97
     98    // classify sources based on moments, brightness.  if a PSF model has been loaded, the PSF
     99    // clump defined for it is used not measured
     100    if (!psphotRoughClass (config, view)) {
     101        psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough classifications");
     102        return psphotReadoutCleanup (config, view);
    114103    }
    115104    if (!strcasecmp (breakPt, "MOMENTS")) {
    116         return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources);
    117     }
    118 
    119     if (!havePSF && !psphotImageQuality (recipe, sources)) {
     105        return psphotReadoutCleanup(config, view);
     106    }
     107
     108    // if we were not supplied a PSF model, determine the IQ stats here
     109    if (!psphotImageQuality (config, view)) {
    120110        psLogMsg("psphot", 3, "failed to measure image quality");
    121         return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources);
    122     }
    123 
    124     // if we were not supplied a PSF, choose one here
    125     if (psf == NULL) {
    126         // use bright stellar objects to measure PSF
    127         // XXX if we do not have enough stars to generate the PSF, build one
    128         // from the SEEING guess and model class
    129         if (!psphotChoosePSF (config, view)) {
    130             psLogMsg ("psphot", 3, "failure to construct a psf model");
    131             return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
    132         }
    133         havePSF = true;
     111        return psphotReadoutCleanup(config, view);
     112    }
     113
     114    // use bright stellar objects to measure PSF if we were supplied a PSF for any input file,
     115    // this step is skipped
     116    if (!psphotChoosePSF (config, view)) {
     117        psLogMsg ("psphot", 3, "failure to construct a psf model");
     118        return psphotReadoutCleanup (config, view);
    134119    }
    135120    if (!strcasecmp (breakPt, "PSFMODEL")) {
    136         return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
    137     }
    138     // move into psphotChoosePSF
    139     psphotVisualShowPSFModel (readout, psf);
     121        return psphotReadoutCleanup (config, view);
     122    }
    140123
    141124    // include externally-supplied sources
     125    // XXX hmm...
    142126    psphotLoadExtSources (config, view, sources);
    143127
    144128    // construct an initial model for each object, set the radius to fitRadius, set circular fit mask
    145     psphotGuessModels (config, readout, sources, psf);
     129    psphotGuessModels (config, view);
    146130
    147131    // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
    148     psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE);
    149 
    150     // We have to place these visualizations here because the models are not realized until
    151     // psphotGuessModels or fitted until psphotFitSourcesLinear.
    152     psphotVisualShowPSFStars (recipe, psf, sources);
     132    psphotFitSourcesLinear (config, view, FALSE);
    153133
    154134    // identify CRs and extended sources
    … …  
    157137        goto finish;
    158138    }
    159     psphotVisualShowSatStars (recipe, psf, sources);
    160139
    161140    // non-linear PSF and EXT fit to brighter sources
    162141    // replace model flux, adjust mask as needed, fit, subtract the models (full stamp)
    163     psphotBlendFit (config, readout, sources, psf);
     142    psphotBlendFit (config, view);
    164143
    165144    // replace all sources
    166     psphotReplaceAllSources (sources, recipe);
     145    psphotReplaceAllSources (config, view);
    167146
    168147    // linear fit to include all sources (subtract again)
    169     psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE);
     148    psphotFitSourcesLinear (config, view, TRUE);
    170149
    171150    // if we only do one pass, skip to extended source analysis
    … …  
    176155
    177156    // add noise for subtracted objects
    178     psphotAddNoise (readout, sources, recipe);
     157    psphotAddNoise (config, view);
    179158
    180159    // find fainter sources (pass 2)
    181     detections = psphotFindDetections (detections, readout, recipe);
     160    // XXX need to distinguish old from new sources
     161    psphotFindDetections (config, view);
    182162
    183163    // remove noise for subtracted objects (ie, return to normal noise level)
    184     psphotSubNoise (readout, sources, recipe);
     164    // XXX this needs to operate only on the OLD sources
     165    psphotSubNoise (config, view);
    185166
    186167    // define new sources based on only the new peaks
    187     psArray *newSources = psphotSourceStats (config, readout, detections, false);
     168    // XXX need to distinguish old from new sources
     169    psphotSourceStats (config, view, false);
    188170
    189171    // set source type
    190     if (!psphotRoughClass (config, view, havePSF)) {
     172    if (!psphotRoughClass (config, view)) {
    191173        psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
    192         return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
     174        return psphotReadoutCleanup (config, view);
    193175    }
    194176
    195177    // create full input models, set the radius to fitRadius, set circular fit mask
    196     psphotGuessModels (config, readout, newSources, psf);
     178    // XXX need to distinguish old from new sources
     179    psphotGuessModels (config, view);
    197180
    198181    // replace all sources so fit below applies to all at once
    199     psphotReplaceAllSources (sources, recipe);
     182    psphotReplaceAllSources (config, view);
    200183
    201184    // merge the newly selected sources into the existing list
    202     psphotMergeSources (sources, newSources);
     185    // XXX old vs new????
     186FIX:    psphotMergeSources (sources, newSources);
    203187    psFree (newSources);
    204188
    205189    // linear fit to all sources
    206     psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE);
     190    psphotFitSourcesLinear (config, view, TRUE);
    207191
    208192pass1finish:
    … …  
    211195    psphotSourceSize (config, view);
    212196
    213     psphotExtendedSourceAnalysis (readout, sources, recipe);
    214 
    215     psphotExtendedSourceFits (readout, sources, recipe);
     197FIX:    psphotExtendedSourceAnalysis (readout, sources, recipe);
     198FIX:    psphotExtendedSourceFits (readout, sources, recipe);
    216199
    217200finish:
    … …  
    221204
    222205    // measure aperture photometry corrections
    223     if (!psphotApResid (config, readout, sources, psf)) {
     206FIX:    if (!psphotApResid (config, readout, sources, psf)) {
    224207        psLogMsg ("psphot", 3, "failed on psphotApResid");
    225         return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
     208        return psphotReadoutCleanup (config, view);
    226209    }
    227210
    228211    // calculate source magnitudes
    229     psphotMagnitudes(config, readout, view, sources, psf);
    230 
    231     if (!psphotEfficiency(config, readout, view, psf, recipe, sources)) {
     212FIX:    psphotMagnitudes(config, readout, view, sources, psf);
     213
     214FIX:    if (!psphotEfficiency(config, readout, view, psf, recipe, sources)) {
    232215        psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources");
    233216        psErrorClear();
    … …  
    238221
    239222    // replace background in residual image
    240     psphotSkyReplace (config, view);
     223FIX:    psphotSkyReplace (config, view);
    241224
    242225    // drop the references to the image pixels held by each source
    243     psphotSourceFreePixels (sources);
     226FIX:    psphotSourceFreePixels (sources);
    244227
    245228    // create the exported-metadata and free local data
    246     return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources);
     229    return psphotReadoutCleanup(config, view);
    247230}
    248 

  • branches/eam_branches/psphot.stack.20100120/src/psphotReadoutCleanup.c

    r26643 r26681  
    11# include "psphotInternal.h"
    22
    3 // psphotReadoutCleanup is called on exit from psphotReadout.  If the last raised error is
    4 // not a DATA error, then there was a serious problem.  Only in this case, or if the fail
    5 // on the stats measurement, do we return false
    6 bool psphotReadoutCleanup (pmConfig *config, pmReadout *readout, psMetadata *recipe, pmDetections *detections, pmPSF *psf, psArray *sources) {
     3// for now, let's store the detections on the readout->analysis for each readout
     4bool psphotReadoutCleanup (pmConfig *config, const pmFPAview *view)
     5{
     6    bool status = true;
    77
    88    // remove internal pmFPAfiles, if created
    … …  
    1212    }
    1313    if (psErrorCodeLast() != PS_ERR_NONE) {
    14         psFree (psf);
    15         psFree (sources);
    16         psFree (detections);
    1714        return false;
    1815    }
     16
     17    int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     18    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     19
     20    // loop over the available readouts
     21    for (int i = 0; i < num; i++) {
     22        if (!psphotReadoutCleanupReadout (config, view, "PSPHOT.INPUT", i)) {
     23            psError (PSPHOT_ERR_CONFIG, false, "failed on psphotReadoutCleanup for PSPHOT.INPUT entry %d", i);
     24            return false;
     25        }
     26    }
     27
     28    // XXX move this to top of loop?
     29    pmKapaClose ();
     30
     31    return true;
     32}
     33
     34// psphotReadoutCleanup is called on exit from psphotReadout.  If the last raised error is
     35// not a DATA error, then there was a serious problem.  Only in this case, or if the fail
     36// on the stats measurement, do we return false
     37bool psphotReadoutCleanupReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index) {
     38
     39    pmReadout *readout, psMetadata *recipe, pmDetections *detections, pmPSF *psf, psArray *sources;
     40
     41    // select the appropriate recipe information
     42    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     43    psAssert (recipe, "missing recipe?");
     44
     45    // find the currently selected readout
     46    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     47    psAssert (readout, "missing file?");
     48
     49    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     50    psAssert (readout, "missing readout?");
     51
     52    // when psphotReadoutCleanup is called, these are not necessarily defined
     53    psArray *sources         = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES");
     54    pmPSF *psf               = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
     55    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     56
     57    // XXX where do we free these, in here (psMetadataRemove?)
    1958
    2059    // use the psf-model to measure FWHM stats
    … …  
    2261        if (!psphotPSFstats (readout, recipe, psf)) {
    2362            psError(PSPHOT_ERR_PROG, false, "Failed to measure PSF shape parameters");
    24             psFree (psf);
    25             psFree (sources);
    26             psFree (detections);
    2763            return false;
    2864        }
    … …  
    3268        if (!psphotMomentsStats (readout, recipe, sources)) {
    3369            psError(PSPHOT_ERR_PROG, false, "Failed to measure Moment shape parameters");
    34             psFree (psf);
    35             psFree (sources);
    36             psFree (detections);
    3770            return false;
    3871        }
    … …  
    4073
    4174    // Check to see if the image quality was measured
     75    // XXX not sure we want / need this test
    4276    if (!psf) {
    4377        bool mdok;                      // Status of MD lookup
    … …  
    4579        if (!mdok || nIQ <= 0) {
    4680            psError(PSPHOT_ERR_DATA, false, "Unable to measure image quality");
    47             psFree (psf);
    48             psFree (sources);
    49             psFree (detections);
    5081            return false;
    5182        }
    5283    }
    5384
     85    // create an output header with stats results XXX this needs to be reworked :
     86    // psphotImageQuality and others need to write these values to a header.  that needs to be
     87    // stored on the readout->analysis
     88    psMetadata *header = psphotDefineHeader (recipe);
    5489
    5590    // write NSTARS to the image header
    56     psphotSetHeaderNstars (recipe, sources);
    57 
    58     // create an output header with stats results
    59     psMetadata *header = psphotDefineHeader (recipe);
     91    psphotSetHeaderNstars (header, sources);
    6092
    6193    // save the results of the analysis
     94    // this should happen way up stream (when needed?)
    6295    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.HEADER",  PS_DATA_METADATA, "header stats", header);
    6396
    6497    if (psf) {
     98        // XXX this seems a little silly : we saved the psf on readout->analysis above, but now
     99        // we are moving it to chip->analysis.
    65100        // save the psf for possible output.  if there was already an entry, it was loaded from external sources
    66101        // the new one may have been updated or modified, so replace the existing entry.  We
    … …  
    77112    }
    78113
    79     // XXX move this to top of loop?
    80     pmKapaClose ();
    81 
    82     psFree (detections);
    83     psFree (psf);
    84114    psFree (header);
    85     psFree (sources);
    86 
    87115    return true;
    88116}

  • branches/eam_branches/psphot.stack.20100120/src/psphotReplaceUnfit.c

    r25755 r26681  
    2222}
    2323
    24 bool psphotReplaceAllSources (psArray *sources, psMetadata *recipe) {
     24// for now, let's store the detections on the readout->analysis for each readout
     25bool psphotReplaceAllSources (pmConfig *config, const pmFPAview *view)
     26{
     27    bool status = true;
     28
     29    // select the appropriate recipe information
     30    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     31    psAssert (recipe, "missing recipe?");
     32
     33    int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     34    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     35
     36    // loop over the available readouts
     37    for (int i = 0; i < num; i++) {
     38        if (!psphotReplaceAllSourcesReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
     39            psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for PSPHOT.INPUT entry %d", i);
     40            return false;
     41        }
     42    }
     43    return true;
     44}
     45
     46bool psphotReplaceAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) {
    2547
    2648    bool status;
    … …  
    2951    psTimerStart ("psphot.replace");
    3052
     53    // find the currently selected readout
     54    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     55    psAssert (readout, "missing file?");
     56
     57    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     58    psAssert (readout, "missing readout?");
     59
     60    psArray *sources = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES");
     61    psAssert (sources, "missing sources?");
     62
    3163    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    3264    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
    33     assert (maskVal);
     65    psAssert (maskVal, "missing mask value?");
    3466
    3567    for (int i = 0; i < sources->n; i++) {

  • branches/eam_branches/psphot.stack.20100120/src/psphotRoughClass.c

    r26643 r26681  
    77        } }
    88
    9 bool psphotRoughClassRegion (int nRegion, psRegion *region, psArray *sources, psMetadata *target, psMetadata *recipe, const bool havePSF);
     9// for now, let's store the detections on the readout->analysis for each readout
     10bool psphotRoughClass (pmConfig *config, const pmFPAview *view)
     11{
     12    bool status = true;
    1013
    11 // 2006.02.02 : no leaks
    12 bool psphotRoughClassReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, const bool havePSF) {
     14    // select the appropriate recipe information
     15    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     16    psAssert (recipe, "missing recipe?");
     17
     18    int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     19    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     20
     21    // loop over the available readouts
     22    for (int i = 0; i < num; i++) {
     23        if (!psphotRoughClassReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
     24            psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for PSPHOT.INPUT entry %d", i);
     25            return false;
     26        }
     27    }
     28    return true;
     29}
     30
     31bool psphotRoughClassReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) {
    1332
    1433    bool status;
    … …  
    2342    psAssert (readout, "missing readout?");
    2443
     44    // if we have a PSF, use the existing PSF clump region below
     45    bool havePSF = false;
     46    if (psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF")) {
     47        havePSF = true;
     48    }
     49
    2550    psArray *sources = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES");
    2651    psAssert (sources, "missing sources?");
    2752
    28     // select the appropriate recipe information
    29     psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
    30     psAssert (recipe, "missing recipe?");
     53    if (!sources->n) {
     54        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping rough classification");
     55        return true;
     56    }
    3157
    3258    // we make this measurement on a NxM grid of regions across the readout
    … …  
    125151    return true;
    126152}
    127 
    128 // for now, let's store the detections on the readout->analysis for each readout
    129 bool psphotRoughClass (pmConfig *config, const pmFPAview *view, const bool havePSF)
    130 {
    131     bool status = true;
    132 
    133     int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
    134     psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
    135 
    136     // loop over the available readouts
    137     for (int i = 0; i < num; i++) {
    138         if (!psphotRoughClassReadout (config, view, "PSPHOT.INPUT", i, havePSF)) {
    139             psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for PSPHOT.INPUT entry %d", i);
    140             return false;
    141         }
    142     }
    143     return true;
    144 }

  • branches/eam_branches/psphot.stack.20100120/src/psphotSourceSize.c

    r26643 r26681  
    1515} psphotSourceSizeOptions;
    1616
     17// local functions:
     18bool psphotSourceSizeReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index);
    1719bool psphotSourceSizePSF (psphotSourceSizeOptions *options, psArray *sources, pmPSF *psf);
    1820bool psphotSourceClass (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options);
    1921bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options);
    2022bool psphotSourceSizeCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options);
    21 bool psphotMaskCosmicRay (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask);
    22 bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask);
     23bool psphotMaskCosmicRay (pmReadout *readout, pmSource *source, psImageMaskType maskVal);
     24bool psphotMaskCosmicRayFootprintCheck (psArray *sources);
    2325
    2426// we need to call this function after sources have been fitted to the PSF model and
    … …  
    3335    bool status = true;
    3436
     37    // select the appropriate recipe information
     38    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     39    psAssert (recipe, "missing recipe?");
     40
    3541    int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
    3642    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
    … …  
    3844    // loop over the available readouts
    3945    for (int i = 0; i < num; i++) {
    40         if (!psphotSourceSizeReadout (config, view, "PSPHOT.INPUT", i)) {
     46        if (!psphotSourceSizeReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
    4147            psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for PSPHOT.INPUT entry %d", i);
    4248            return false;
    … …  
    4753
    4854// XXX use an internal flag to mark sources which have already been measured
    49 // how to track the sources already measured?
    50 bool psphotSourceSize(pmConfig *config, const pmFPAview *view, const char *filename, int index)
     55bool psphotSourceSizeReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe)
    5156{
    5257    bool status;
    … …  
    6570    psAssert (sources, "missing sources?");
    6671
     72    if (!sources->n) {
     73        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source size");
     74        return true;
     75    }
     76
    6777    pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
    6878    psAssert (psf, "missing psf?");
    69 
    70     // select the appropriate recipe information
    71     psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
    72     psAssert (recipe, "missing recipe?");
    7379
    7480    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    … …  
    106112
    107113    // We are using the value psfMag - 2.5*log10(moment.Sum) as a measure of the extendedness
    108     // of and object.  We need to model this distribution for the PSF stars before we can test
     114    // of an object.  We need to model this distribution for the PSF stars before we can test
    109115    // the significance for a specific object
    110116    // XXX move this to the code that generates the PSF?
    … …  
    122128    psphotVisualShowSourceSize (readout, sources);
    123129    psphotVisualPlotApResid (sources, options.ApResid, options.ApSysErr);
    124 
    125     return true;
    126 }
    127 
    128 bool psphotMaskCosmicRay (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) {
     130    psphotVisualShowSatStars (recipe, psf, sources);
     131
     132    return true;
     133}
     134
     135// model the apmifit distribution for the psf stars:
     136bool psphotSourceSizePSF (psphotSourceSizeOptions *options, psArray *sources, pmPSF *psf) {
     137
     138    // select stats from the psf stars
     139    psVector *Ap = psVectorAllocEmpty (100, PS_TYPE_F32);
     140    psVector *ApErr = psVectorAllocEmpty (100, PS_TYPE_F32);
     141
     142    psImageMaskType maskVal = options->maskVal | options->markVal;
     143
     144    // XXX  why PHOT_WEIGHT??
     145    pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT;
     146
     147    for (int i = 0; i < sources->n; i++) {
     148        pmSource *source = sources->data[i];
     149        if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
     150
     151        // replace object in image
     152        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
     153            pmSourceAdd (source, PM_MODEL_OP_FULL, options->maskVal);
     154        }
     155
     156        // clear the mask bit and set the circular mask pixels
     157        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal));
     158        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", options->markVal);
     159
     160        // XXX can we test if psfMag is set and calculate only if needed?
     161        pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal
     162
     163        // clear the mask bit
     164        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal));
     165
     166        // re-subtract the object, leave local sky
     167        pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal);
     168
     169        float apMag = -2.5*log10(source->moments->Sum);
     170        float dMag = source->psfMag - apMag;
     171
     172        psVectorAppend (Ap, 100, dMag);
     173        psVectorAppend (ApErr, 100, source->errMag);
     174    }
     175
     176    // model the distribution as a mean or median value and a systematic error from that value:
     177    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
     178    psVectorStats (stats, Ap, NULL, NULL, 0);
     179
     180    psVector *dAp = psVectorAlloc (Ap->n, PS_TYPE_F32);
     181    for (int i = 0; i < Ap->n; i++) {
     182        dAp->data.F32[i] = Ap->data.F32[i] - stats->robustMedian;
     183    }
     184
     185    options->ApResid = stats->robustMedian;
     186    options->ApSysErr = psVectorSystematicError(dAp, ApErr, 0.05);
     187    psLogMsg ("psphot", PS_LOG_DETAIL, "psf - Sum: %f +/- %f\n", options->ApResid, options->ApSysErr);
     188
     189    psFree (Ap);
     190    psFree (ApErr);
     191    psFree (stats);
     192    psFree (dAp);
     193
     194    return true;
     195}
     196
     197// classify sources based on the combination of psf-mag, Mxx, Myy
     198bool psphotSourceClass (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options) {
     199
     200    bool status;
     201    pmPSFClump psfClump;
     202    char regionName[64];
     203
     204    psLogMsg("psModules.objects", PS_LOG_INFO, "Source Size classifications: %4s %4s %4s %4s %4s %4s", "Npsf", "Next", "Nsat", "Ncr", "Nmiss", "Nskip");
     205
     206    int nRegions = psMetadataLookupS32 (&status, readout->analysis, "PSF.CLUMP.NREGIONS");
     207    for (int i = 0; i < nRegions; i ++) {
     208        snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i);
     209        psMetadata *regionMD = psMetadataLookupPtr (&status, readout->analysis, regionName);
     210        psAssert (regionMD, "regions must be defined by earlier call to psphotRoughClassRegion");
     211
     212        psRegion *region = psMetadataLookupPtr (&status, regionMD, "REGION");
     213        psAssert (region, "regions must be defined by earlier call to psphotRoughClassRegion");
     214
     215        // pull FWHM_X,Y from the recipe, use to define psfClump.X,Y
     216        psfClump.X  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");   psAssert (status, "missing PSF.CLUMP.X");
     217        psfClump.Y  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");   psAssert (status, "missing PSF.CLUMP.Y");
     218        psfClump.dX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");  psAssert (status, "missing PSF.CLUMP.DX");
     219        psfClump.dY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");  psAssert (status, "missing PSF.CLUMP.DY");
     220
     221        if ((psfClump.X < 0) || (psfClump.Y < 0) || !psfClump.X || !psfClump.Y || isnan(psfClump.X) || isnan(psfClump.Y)) {
     222            psLogMsg ("psphot", 4, "Failed to find a valid PSF clump for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1);
     223            continue;
     224        }
     225
     226        if (!psphotSourceClassRegion (region, &psfClump, sources, recipe, psf, options)) {
     227            psLogMsg ("psphot", 4, "Failed to determine source classification for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1);
     228            continue;
     229        }
     230    }
     231
     232    return true;
     233}
     234
     235bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options) {
     236
     237    PS_ASSERT_PTR_NON_NULL(sources, false);
     238    PS_ASSERT_PTR_NON_NULL(recipe, false);
     239
     240    int Nsat  = 0;
     241    int Next  = 0;
     242    int Npsf  = 0;
     243    int Ncr   = 0;
     244    int Nmiss = 0;
     245    int Nskip = 0;
     246
     247    pmSourceMode noMoments = PM_SOURCE_MODE_MOMENTS_FAILURE | PM_SOURCE_MODE_SKYVAR_FAILURE | PM_SOURCE_MODE_SKY_FAILURE | PM_SOURCE_MODE_BELOW_MOMENTS_SN;
     248    pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT;
     249
     250    psImageMaskType maskVal = options->maskVal | options->markVal;
     251
     252    for (psS32 i = 0 ; i < sources->n ; i++) {
     253
     254        pmSource *source = (pmSource *) sources->data[i];
     255
     256        // psfClumps are found for image subregions:
     257        // skip sources not in this region
     258        if (source->peak->x <  region->x0) continue;
     259        if (source->peak->x >= region->x1) continue;
     260        if (source->peak->y <  region->y0) continue;
     261        if (source->peak->y >= region->y1) continue;
     262
     263        // skip source if it was already measured
     264        if (source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED) {
     265            psTrace("psphot", 7, "Not calculating source size since it has already been measured\n");
     266            continue;
     267        }
     268
     269        // source must have been subtracted
     270        if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) {
     271            source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED;
     272            psTrace("psphot", 7, "Not calculating source size since source is not subtracted\n");
     273            continue;
     274        }
     275
     276        // we are basically classifying by moments; use the default if not found
     277        psAssert (source->moments, "why is this source missing moments?");
     278        if (source->mode & noMoments) {
     279            Nskip ++;
     280            continue;
     281        }
     282
     283        psF32 Mxx = source->moments->Mxx;
     284        psF32 Myy = source->moments->Myy;
     285
     286        // replace object in image
     287        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
     288            pmSourceAdd (source, PM_MODEL_OP_FULL, options->maskVal);
     289        }
     290
     291        // clear the mask bit and set the circular mask pixels
     292        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal));
     293        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", options->markVal);
     294
     295        // XXX can we test if psfMag is set and calculate only if needed?
     296        pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal
     297
     298        // clear the mask bit
     299        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal));
     300
     301        // re-subtract the object, leave local sky
     302        pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal);
     303
     304        float apMag = -2.5*log10(source->moments->Sum);
     305        float dMag = source->psfMag - apMag;
     306        float nSigma = (dMag - options->ApResid) / hypot(source->errMag, options->ApSysErr);
     307
     308        source->extNsigma = nSigma;
     309        source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
     310
     311        // Anything within this region is a probably PSF-like object. Saturated stars may land
     312        // in this region, but are detected elsewhere on the basis of their peak value.
     313        bool isPSF = ((fabs(nSigma) < options->nSigmaApResid) &&
     314                      (fabs(Mxx - psfClump->X) < options->nSigmaMoments*psfClump->dX) &&
     315                      (fabs(Myy - psfClump->Y) < options->nSigmaMoments*psfClump->dY));
     316        if (isPSF) {
     317          psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g PSF\t%g %g\n",
     318                  source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma,
     319                  options->nSigmaApResid,options->nSigmaMoments);
     320          Npsf ++;
     321          continue;
     322        }
     323
     324        // Defects may not always match CRs from peak curvature analysis
     325        // Defects may also be marked as SATSTAR -- XXX deactivate this flag?
     326        // XXX this rule is not great
     327        if ((Mxx < psfClump->X) || (Myy < psfClump->Y)) {
     328
     329          psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g CR\t%g %g\n",
     330                  source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma,
     331                  options->nSigmaApResid,options->nSigmaMoments);
     332            source->mode |= PM_SOURCE_MODE_DEFECT;
     333            Ncr ++;
     334            continue;
     335        }
     336
     337        // saturated star (determined in PSF fit).  These may also be saturated galaxies, or
     338        // just large saturated regions.
     339        if (source->mode & PM_SOURCE_MODE_SATSTAR) {
     340          psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g SAT\t%g %g\n",
     341                  source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma,
     342                  options->nSigmaApResid,options->nSigmaMoments);
     343         
     344            Nsat ++;
     345            continue;
     346        }
     347
     348        // XXX allow the Mxx, Myy to be less than psfClump->X,Y (by some nSigma)?
     349        bool isEXT = (nSigma > options->nSigmaApResid) || ((Mxx > psfClump->X) && (Myy > psfClump->Y));
     350        if (isEXT) {
     351          psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g Ext\t%g %g\n",
     352                  source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma,
     353                  options->nSigmaApResid,options->nSigmaMoments);
     354         
     355            source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
     356            Next ++;
     357            continue;
     358        }
     359        psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g Unk\t%g %g\n",
     360                source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma,
     361                options->nSigmaApResid,options->nSigmaMoments);
     362       
     363        psWarning ("sourse size was missed for %f,%f : %f %f -- %f\n", source->peak->xf, source->peak->yf, Mxx, Myy, nSigma);
     364        Nmiss ++;
     365    }
     366
     367    psLogMsg("psModules.objects", PS_LOG_INFO, "Source Size classifications: %4d %4d %4d %4d %4d %4d", Npsf, Next, Nsat, Ncr, Nmiss, Nskip);
     368
     369    return true;
     370}
     371
     372// given an object suspected to be a defect, generate a pixel mask using the Lapacian transform
     373// if enough of the object is detected as 'sharp', consider the object a cosmic ray
     374bool psphotSourceSizeCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options) {
     375
     376    for (int i = 0; i < sources->n; i++) {
     377        pmSource *source = sources->data[i];
     378
     379        // Integer position of peak
     380        int xPeak = source->peak->xf - source->pixels->col0 + 0.5;
     381        int yPeak = source->peak->yf - source->pixels->row0 + 0.5;
     382
     383        // Skip sources which are too close to a boundary.  These are mostly caught as DEFECT
     384        if (xPeak < 1 || xPeak > source->pixels->numCols - 2 ||
     385            yPeak < 1 || yPeak > source->pixels->numRows - 2) {
     386            psTrace("psphot", 7, "Not calculating crNsigma due to edge\n");
     387            //      psTrace("psphot.czw", 2, "Not calculating crNsigma due to edge\n");
     388            continue;
     389        }
     390       
     391        // this source is thought to be a cosmic ray.  flag the detection and mask the pixels
     392        if (source->mode & PM_SOURCE_MODE_DEFECT) {
     393            // XXX this is running slowly and is too agressive, but it more-or-less works
     394            // psphotMaskCosmicRayCZW(readout, source, options->crMask);
     395           
     396            // XXX these are the old versions which we are not using
     397            // psphotMaskCosmicRay (readout->mask, source, maskVal, crMask);
     398            // psphotMaskCosmicRay_Old (source, maskVal, crMask);
     399        }
     400    }
     401   
     402    // now that we have masked pixels associated with CRs, we can grow the mask
     403    if (options->grow > 0) {
     404        bool oldThreads = psImageConvolveSetThreads(true); // Old value of threading for psImageConvolveMask
     405        psImage *newMask = psImageConvolveMask(NULL, readout->mask, options->crMask, options->crMask, -options->grow, options->grow, -options->grow, options->grow);
     406        psImageConvolveSetThreads(oldThreads);
     407        if (!newMask) {
     408            psError(PS_ERR_UNKNOWN, false, "Unable to grow CR mask");
     409            return false;
     410        }
     411        psFree(readout->mask);
     412        readout->mask = newMask;
     413    }
     414
     415    // XXX test : save the mask image
     416    if (0) {
     417        psphotSaveImage (NULL, readout->mask,   "mask.fits");
     418    }
     419    return true;
     420}
     421
     422// Comments by CZW 20091209 : Mechanics of how to identify CR pixels taken from "Cosmic-Ray
     423// Rejection by Laplacian Edge Detection" by Pieter van Dokkum, arXiv:astro-ph/0108003.  This
     424// does no repair or recovery of the CR pixels, it only masks them out.  My test code can be
     425// found at /data/ipp031.0/watersc1/psphot.20091209/algo_check.c
     426bool psphotMaskCosmicRay (pmReadout *readout, pmSource *source, psImageMaskType maskVal) {
     427
     428    // Get the actual images and information about the peak.
     429    psImage *mask = readout->mask;
     430    pmPeak *peak = source->peak;
     431    pmFootprint *footprint = peak->footprint;
     432
     433    int xm = footprint->bbox.x0;
     434    int xM = footprint->bbox.x1;
     435    int ym = footprint->bbox.y0;
     436    int yM = footprint->bbox.y1;
     437
     438    if (xm < 0) { xm = 0; }
     439    if (ym < 0) { ym = 0; }
     440    if (xM > mask->numCols) { xM = mask->numCols; }
     441    if (yM > mask->numRows) { yM = mask->numRows; }
     442    int dx = xM - xm;
     443    int dy = yM - ym;
     444
     445    // Bounding boxes are inclusive of final pixel
     446    // XXX ensure dx,dy do not become too large here
     447    dx++;
     448    dy++;
     449
     450    psImage *image= readout->image;
     451    psImage *variance = readout->variance;
     452     
     453    int binning = 1;
     454    float sigma_thresh = 2.0;
     455    int iteration = 0;
     456    int max_iter = 5;
     457    float noise_factor = 5.0 / 4.0;  // Intrinsic to the Laplacian making noise spikes spikier.
     458
     459    // Temporary images.
     460    psImage *mypix  = psImageAlloc(dx,dy,image->type.type);
     461    psImage *myvar  = psImageAlloc(dx,dy,image->type.type);
     462    psImage *binned = psImageAlloc(dx * binning,dy * binning,image->type.type);
     463    psImage *conved = psImageAlloc(dx * binning,dy * binning,image->type.type);
     464    psImage *edges  = psImageAlloc(dx,dy,image->type.type);
     465    psImage *mymask = psImageAlloc(dx,dy,PS_TYPE_IMAGE_MASK);
     466     
     467    int x,y;
     468    // Load my copy of things.
     469    for (x = 0; x < dx; x++) {
     470        for (y = 0; y < dy; y++) {
     471            psImageSet(mypix,x,y,psImageGet(image,x+xm,y+ym));
     472            psImageSet(myvar,x,y,psImageGet(variance,x+xm,y+ym));
     473            mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0x00;
     474        }
     475    }
     476    // Mask so I can see on the output image where the footprint is.
     477    for (int i = 0; i < footprint->spans->n; i++) {
     478        pmSpan *sp = footprint->spans->data[i];
     479        for (int j = sp->x0; j <= sp->x1; j++) {
     480            y = sp->y - ym;
     481            x = j - xm;
     482            mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= 0x01;
     483        }
     484    }
     485
     486    int CRpix_count = 0;     
     487    do {
     488        CRpix_count = 0;
     489        // Zero out my temp images.
     490        for (x = 0; x < binning * dx; x++) {
     491            for (y = 0; y < binning * dy; y++) {
     492                psImageSet(binned,x,y,0.0);
     493                psImageSet(conved,x,y,0.0);
     494                psImageSet(edges,x/binning,y/ binning,0.0);
     495            }
     496        }     
     497        // Make subsampled image. Maybe this should be called "unbinned" or something
     498        for (x = 0; x < binning * dx; x++) {
     499            for (y = 0; y < binning * dy; y++) {
     500                psImageSet(binned,x,y,psImageGet(mypix,x / binning,y / binning));
     501            }
     502        }
     503        // Apply Laplace transform (kernel = [[0 -0.25 0][-0.25 1 -0.25][0 -0.25 0]]), clipping at zero
     504        for (x = 1; x < dx - 1; x++) {
     505            for (y = 1; y < dy - 1; y++) {
     506                psImageSet(conved,x,y,psImageGet(binned,x,y) - 0.25 *
     507                           (psImageGet(binned,x-1,y) + psImageGet(binned,x+1,y) +
     508                            psImageGet(binned,x,y-1) + psImageGet(binned,x,y+1)));
     509                if (psImageGet(conved,x,y) < 0.0) {
     510                    psImageSet(conved,x,y,0.0);
     511                }
     512            }
     513        }
     514        // Create an edge map by rebinning
     515        for (x = 0; x < binning * dx; x++) {
     516            for (y = 0; y < binning * dy; y++) {
     517                psImageSet(edges,x / binning, y / binning,
     518                           psImageGet(edges, x / binning, y / binning) +
     519                           psImageGet(conved,x,y));
     520            }
     521        }
     522        // Modify my mask if we're above the significance threshold
     523        for (x = 0; x < dx; x++) {
     524            for (y = 0; y < dy; y++) {
     525                if ( psImageGet(edges,x,y) / (binning * sqrt(noise_factor * psImageGet(myvar,x,y))) > sigma_thresh ) {
     526                    if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x01) {
     527                        mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= 0x40;
     528                        CRpix_count++;
     529                    }
     530                }
     531            }
     532        }
     533
     534        // "Repair" Masked pixels for the next round.
     535        for (x = 1; x < dx - 1; x++) {
     536            for (y = 1; y < dy - 1; y++) {
     537                if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40) {
     538                    psImageSet(mypix,x,y,0.25 *
     539                               (psImageGet(mypix,x-1,y) + psImageGet(mypix,x+1,y) +
     540                                psImageGet(mypix,x,y-1) + psImageGet(mypix,x,y+1)));
     541                }
     542            }
     543        }
     544       
     545# if 0
     546        if ((psTraceGetLevel("psphot.czw") >= 2)&&(abs(xm - 2770) < 5)&&(abs(ym - 2581) < 5)&&(iteration == 0)) {
     547          psTrace("psphot.czw",2,"TRACEMOTRON %d %d %d %d %d\n",xm,ym,dx,dy,iteration);
     548          psphotSaveImage (NULL, mypix,   "czw.pix.fits");
     549          psphotSaveImage (NULL, myvar,   "czw.var.fits");
     550          psphotSaveImage (NULL, binned,  "czw.binn.fits");
     551          psphotSaveImage (NULL, conved,  "czw.conv.fits");
     552          psphotSaveImage (NULL, edges,   "czw.edge.fits");
     553          psphotSaveImage (NULL, mymask,  "czw.mask.fits");
     554        }
     555# endif
     556
     557        psTrace("psphot.czw",2,"Iter: %d Count: %d",iteration,CRpix_count);
     558        iteration++;
     559    } while ((iteration < max_iter)&&(CRpix_count > 0));
     560
     561    // A solitary masked pixel is likely a lie. Remove those.
     562    for (x = 0; x < dx; x++) {
     563        for (y = 0; y < dy; y++) {
     564            if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40) {
     565                if ((x-1 >= 0)&&(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x-1] & 0x40)) {
     566                    continue;
     567                }
     568                if ((y-1 >= 0)&&(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y-1][x] & 0x40)) {
     569                    continue;
     570                }
     571                if ((x+1 < dx)&&(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x+1] & 0x40)) {
     572                    continue;
     573                }
     574                if ((y+1 < dy)&&(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y+1][x] & 0x40)) {
     575                    continue;
     576                }
     577                mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] ^= 0x40;
     578            }
     579        }
     580    }
     581
     582    // transfer temporary mask to real mask
     583    for (x = 0; x < dx; x++) {
     584        for (y = 0; y < dy; y++) {
     585            if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40) {
     586                mask->data.PS_TYPE_IMAGE_MASK_DATA[y+ym+mask->row0][x+xm+mask->col0] |= maskVal;
     587            }
     588        }
     589    }
     590     
     591    // XXX if we decide this REALLY is a cosmic ray, set the CR_LIMIT bit
     592    source->mode |= PM_SOURCE_MODE_CR_LIMIT;
     593
     594    psFree(mypix);
     595    psFree(myvar);
     596    psFree(binned);
     597    psFree(conved);
     598    psFree(edges);
     599    psFree(mymask);
     600     
     601    return true;
     602}
     603
     604bool psphotMaskCosmicRayFootprintCheck (psArray *sources) {
     605
     606    for (int i = 0; i < sources->n; i++) {
     607        pmSource *source = sources->data[i];
     608        pmPeak *peak = source->peak;
     609        pmFootprint *footprint = peak->footprint;
     610        if (!footprint) continue;
     611        for (int j = 0; j < footprint->spans->n; j++) {
     612            pmSpan *sp = footprint->spans->data[j];
     613            psAssert (sp, "missing span");
     614        }
     615    }
     616    return true;
     617}
     618
     619/**** ------ old versions of cosmic ray masking ----- ****/
     620
     621// This attempt to mask the cosmic rays used the isophotal boundary
     622bool psphotMaskCosmicRay_V1 (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) {
    129623
    130624    // replace the source flux
    … …  
    139633    pmFootprint *footprint = peak->footprint;
    140634    if (!footprint) {
     635      psTrace("psphot.czw",2,"Using isophot CR mask code.");
     636     
    141637        // if we have not footprint, use the old code to mask by isophot
    142638        psphotMaskCosmicRayIsophot (source, maskVal, crMask);
    … …  
    145641
    146642    if (!footprint->spans) {
     643      psTrace("psphot.czw",2,"Using isophot CR mask code.");
     644     
    147645        // if we have no footprint, use the old code to mask by isophot
    148646        psphotMaskCosmicRayIsophot (source, maskVal, crMask);
    149647        return true;
    150648    }
    151 
     649    psphotMaskCosmicRayIsophot (source, maskVal, crMask);
    152650    // mask all of the pixels covered by the spans of the footprint
    153     for (int j = 1; j < footprint->spans->n; j++) {
    154         pmSpan *span1 = footprint->spans->data[j];
    155 
    156         int iy = span1->y;
    157         int xs = span1->x0;
    158         int xe = span1->x1;
    159 
    160         for (int ix = xs; ix < xe; ix++) {
    161             mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
    162         }
    163     }
     651    for (int j = 1; j < footprint->spans->n; j++) { 
     652        pmSpan *span1 = footprint->spans->data[j];
     653
     654        int iy = span1->y;
     655        int xs = span1->x0;
     656        int xe = span1->x1;
     657
     658        for (int ix = xs; ix < xe; ix++) {
     659            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
     660        }
     661    } 
    164662    return true;
    165663}
    … …  
    228726        }
    229727    }
    230     return true;
    231 }
    232 
    233 // model the apmifit distribution for the psf stars:
    234 bool psphotSourceSizePSF (psphotSourceSizeOptions *options, psArray *sources, pmPSF *psf) {
    235 
    236     // select stats from the psf stars
    237     psVector *Ap = psVectorAllocEmpty (100, PS_TYPE_F32);
    238     psVector *ApErr = psVectorAllocEmpty (100, PS_TYPE_F32);
    239 
    240     psImageMaskType maskVal = options->maskVal | options->markVal;
    241 
    242     // XXX  why PHOT_WEIGHT??
    243     pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT;
    244 
    245     for (int i = 0; i < sources->n; i++) {
    246         pmSource *source = sources->data[i];
    247         if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
    248 
    249         // replace object in image
    250         if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
    251             pmSourceAdd (source, PM_MODEL_OP_FULL, options->maskVal);
    252         }
    253 
    254         // clear the mask bit and set the circular mask pixels
    255         psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal));
    256         psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", options->markVal);
    257 
    258         // XXX can we test if psfMag is set and calculate only if needed?
    259         pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal
    260 
    261         // clear the mask bit
    262         psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal));
    263 
    264         // re-subtract the object, leave local sky
    265         pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal);
    266 
    267         float apMag = -2.5*log10(source->moments->Sum);
    268         float dMag = source->psfMag - apMag;
    269 
    270         psVectorAppend (Ap, dMag);
    271         psVectorAppend (ApErr, source->errMag);
    272     }
    273 
    274     // model the distribution as a mean or median value and a systematic error from that value:
    275     psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
    276     psVectorStats (stats, Ap, NULL, NULL, 0);
    277 
    278     psVector *dAp = psVectorAlloc (Ap->n, PS_TYPE_F32);
    279     for (int i = 0; i < Ap->n; i++) {
    280         dAp->data.F32[i] = Ap->data.F32[i] - stats->robustMedian;
    281     }
    282 
    283     options->ApResid = stats->robustMedian;
    284     options->ApSysErr = psVectorSystematicError(dAp, ApErr, 0.05);
    285     // XXX this is quite arbitrary...
    286     if (!isfinite(options->ApSysErr)) options->ApSysErr = 0.01;
    287     psLogMsg ("psphot", PS_LOG_DETAIL, "psf - Sum: %f +/- %f\n", options->ApResid, options->ApSysErr);
    288 
    289     psFree (Ap);
    290     psFree (ApErr);
    291     psFree (stats);
    292     psFree (dAp);
    293 
    294     return true;
    295 }
    296 
    297 // classify sources based on the combination of psf-mag, Mxx, Myy
    298 bool psphotSourceClass (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options) {
    299 
    300     bool status;
    301     pmPSFClump psfClump;
    302     char regionName[64];
    303 
    304     psLogMsg("psModules.objects", PS_LOG_INFO, "Source Size classifications: %4s %4s %4s %4s %4s %4s", "Npsf", "Next", "Nsat", "Ncr", "Nmiss", "Nskip");
    305 
    306     int nRegions = psMetadataLookupS32 (&status, readout->analysis, "PSF.CLUMP.NREGIONS");
    307     for (int i = 0; i < nRegions; i ++) {
    308         snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i);
    309         psMetadata *regionMD = psMetadataLookupPtr (&status, readout->analysis, regionName);
    310         psAssert (regionMD, "regions must be defined by earlier call to psphotRoughClassRegion");
    311 
    312         psRegion *region = psMetadataLookupPtr (&status, regionMD, "REGION");
    313         psAssert (region, "regions must be defined by earlier call to psphotRoughClassRegion");
    314 
    315         // pull FWHM_X,Y from the recipe, use to define psfClump.X,Y
    316         psfClump.X  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");   psAssert (status, "missing PSF.CLUMP.X");
    317         psfClump.Y  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");   psAssert (status, "missing PSF.CLUMP.Y");
    318         psfClump.dX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");  psAssert (status, "missing PSF.CLUMP.DX");
    319         psfClump.dY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");  psAssert (status, "missing PSF.CLUMP.DY");
    320 
    321         if ((psfClump.X < 0) || (psfClump.Y < 0) || !psfClump.X || !psfClump.Y || isnan(psfClump.X) || isnan(psfClump.Y)) {
    322             psLogMsg ("psphot", 4, "Failed to find a valid PSF clump for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1);
    323             continue;
    324         }
    325 
    326         if (!psphotSourceClassRegion (region, &psfClump, sources, recipe, psf, options)) {
    327             psLogMsg ("psphot", 4, "Failed to determine source classification for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1);
    328             continue;
    329         }
    330     }
    331 
    332     return true;
    333 }
    334 
    335 bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options) {
    336 
    337     PS_ASSERT_PTR_NON_NULL(sources, false);
    338     PS_ASSERT_PTR_NON_NULL(recipe, false);
    339 
    340     int Nsat  = 0;
    341     int Next  = 0;
    342     int Npsf  = 0;
    343     int Ncr   = 0;
    344     int Nmiss = 0;
    345     int Nskip = 0;
    346 
    347     pmSourceMode noMoments = PM_SOURCE_MODE_MOMENTS_FAILURE | PM_SOURCE_MODE_SKYVAR_FAILURE | PM_SOURCE_MODE_SKY_FAILURE | PM_SOURCE_MODE_BELOW_MOMENTS_SN;
    348     pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT;
    349 
    350     psImageMaskType maskVal = options->maskVal | options->markVal;
    351 
    352     for (psS32 i = 0 ; i < sources->n ; i++) {
    353 
    354         pmSource *source = (pmSource *) sources->data[i];
    355 
    356         // psfClumps are found for image subregions:
    357         // skip sources not in this region
    358         if (source->peak->x <  region->x0) continue;
    359         if (source->peak->x >= region->x1) continue;
    360         if (source->peak->y <  region->y0) continue;
    361         if (source->peak->y >= region->y1) continue;
    362 
    363         // skip source if it was already measured
    364         if (source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED) {
    365             psTrace("psphot", 7, "Not calculating source size since it has already been measured\n");
    366             continue;
    367         }
    368 
    369         // source must have been subtracted
    370         if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) {
    371             source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED;
    372             psTrace("psphot", 7, "Not calculating source size since source is not subtracted\n");
    373             continue;
    374         }
    375 
    376         // we are basically classifying by moments; use the default if not found
    377         psAssert (source->moments, "why is this source missing moments?");
    378         if (source->mode & noMoments) {
    379             Nskip ++;
    380             continue;
    381         }
    382 
    383         psF32 Mxx = source->moments->Mxx;
    384         psF32 Myy = source->moments->Myy;
    385 
    386         // replace object in image
    387         if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
    388             pmSourceAdd (source, PM_MODEL_OP_FULL, options->maskVal);
    389         }
    390 
    391         // clear the mask bit and set the circular mask pixels
    392         psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal));
    393         psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", options->markVal);
    394 
    395         // XXX can we test if psfMag is set and calculate only if needed?
    396         pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal
    397 
    398         // clear the mask bit
    399         psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal));
    400 
    401         // re-subtract the object, leave local sky
    402         pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal);
    403 
    404         float apMag = -2.5*log10(source->moments->Sum);
    405         float dMag = source->psfMag - apMag;
    406         float nSigma = (dMag - options->ApResid) / hypot(source->errMag, options->ApSysErr);
    407 
    408         source->extNsigma = nSigma;
    409         source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
    410 
    411         // Anything within this region is a probably PSF-like object. Saturated stars may land
    412         // in this region, but are detected elsewhere on the basis of their peak value.
    413         bool isPSF = (fabs(nSigma) < options->nSigmaApResid) && (fabs(Mxx - psfClump->X) < options->nSigmaMoments*psfClump->dX) && (fabs(Myy - psfClump->Y) < options->nSigmaMoments*psfClump->dY);
    414         if (isPSF) {
    415             Npsf ++;
    416             continue;
    417         }
    418 
    419         // Defects may not always match CRs from peak curvature analysis
    420         // Defects may also be marked as SATSTAR -- XXX deactivate this flag?
    421         // XXX this rule is not great
    422         if ((Mxx < psfClump->X) || (Myy < psfClump->Y)) {
    423             source->mode |= PM_SOURCE_MODE_DEFECT;
    424             Ncr ++;
    425             continue;
    426         }
    427 
    428         // saturated star (determined in PSF fit).  These may also be saturated galaxies, or
    429         // just large saturated regions.
    430         if (source->mode & PM_SOURCE_MODE_SATSTAR) {
    431             Nsat ++;
    432             continue;
    433         }
    434 
    435         // XXX allow the Mxx, Myy to be less than psfClump->X,Y (by some nSigma)?
    436         bool isEXT = (nSigma > options->nSigmaApResid) || ((Mxx > psfClump->X) && (Myy > psfClump->Y));
    437         if (isEXT) {
    438             source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
    439             Next ++;
    440             continue;
    441         }
    442 
    443         psWarning ("sourse size was missed for %f,%f : %f %f -- %f\n", source->peak->xf, source->peak->yf, Mxx, Myy, nSigma);
    444         Nmiss ++;
    445     }
    446 
    447     psLogMsg("psModules.objects", PS_LOG_INFO, "Source Size classifications: %4d %4d %4d %4d %4d %4d", Npsf, Next, Nsat, Ncr, Nmiss, Nskip);
    448 
    449728    return true;
    450729}
    … …  
    524803}
    525804
    526 bool psphotSourceSizeCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options) {
     805// this was an old attempt to identify cosmic rays based on the peak curvature
     806bool psphotSourcePeakCurvature (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options) {
    527807
    528808    // classify the sources based on the CR test (place this in a function?)
    … …  
    656936    return true;
    657937}
     938

  • branches/eam_branches/psphot.stack.20100120/src/psphotSourceStats.c

    r26643 r26681  
    11# include "psphotInternal.h"
    22
    3 bool psphotSetMomentsWindow (psMetadata *recipe, psMetadata *analysis, psArray *sources);
    4 
    5 bool psphotSourceStatsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, bool setWindow) {
    6 
    7     bool status = false;
    8     psArray *sources = NULL;
    9 
    10     psTimerStart ("psphot.stats");
    11 
    12     // find the currently selected readout
    13     pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
    14     psAssert (readout, "missing file?");
    15 
    16     pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
    17     psAssert (readout, "missing readout?");
    18 
    19     pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
    20     psAssert (detections, "missing detections?");
     3// convert detections to sources and measure their basic properties (moments, local sky, sky variance)
     4bool psphotSourceStats (pmConfig *config, const pmFPAview *view, bool setWindow)
     5{
     6    bool status = true;
    217
    228    // select the appropriate recipe information
    239    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
    2410    psAssert (recipe, "missing recipe?");
     11
     12    int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     13    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     14
     15    // loop over the available readouts
     16    for (int i = 0; i < num; i++) {
     17        if (!psphotSourceStatsReadout (config, view, "PSPHOT.INPUT", i, recipe, setWindow)) {
     18            psError (PSPHOT_ERR_CONFIG, false, "failed to find initial detections for PSPHOT.INPUT entry %d", i);
     19            return false;
     20        }
     21    }
     22    return true;
     23}
     24
     25bool psphotSourceStatsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool setWindow) {
     26
     27    bool status = false;
     28    psArray *sources = NULL;
     29
     30    psTimerStart ("psphot.stats");
     31
     32    // find the currently selected readout
     33    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     34    psAssert (file, "missing file?");
     35
     36    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     37    psAssert (readout, "missing readout?");
     38
     39    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     40    psAssert (detections, "missing detections?");
    2541
    2642    // determine the number of allowed threads
    … …  
    3854
    3955    char *breakPt  = psMetadataLookupStr (&status, recipe, "BREAK_POINT");
    40     if (!status) return false;
     56    psAssert (status, "missing BREAK_POINT?");
    4157
    4258    psArray *peaks = detections->peaks;
    … …  
    4864    // generate the array of sources, define the associated pixel
    4965    sources = psArrayAllocEmpty (peaks->n);
     66
     67    // if there are no peaks, we save the empty source array and return
     68    if (!peaks->n) {
     69        // save sources on the readout->analysis
     70        if (!psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY, "psphot sources", sources)) {
     71            psError (PSPHOT_ERR_UNKNOWN, false, "problem saving sources on readout");
     72            psFree(sources);
     73            return false;
     74        }
     75        psFree(sources);
     76        return true;
     77    }
    5078
    5179    for (int i = 0; i < peaks->n; i++) {
    … …  
    485513// otherwise it only contains the new peaks.
    486514
    487 // for now, let's store the detections on the readout->analysis for each readout
    488 bool psphotSourceStats (pmConfig *config, const pmFPAview *view, bool setWindow)
    489 {
    490     bool status = true;
    491 
    492     int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
    493     psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
    494 
    495     // loop over the available readouts
    496     for (int i = 0; i < num; i++) {
    497         if (!psphotFindDetectionsReadout (config, view, "PSPHOT.INPUT", i, setWindow)) {
    498             psError (PSPHOT_ERR_CONFIG, false, "failed to find initial detections for PSPHOT.INPUT entry %d", i);
    499             return false;
    500         }
    501     }
    502     return true;
    503 }

  • branches/eam_branches/psphot.stack.20100120/src/psphotSubtractBackground.c

    r26542 r26681  
    44// generate the median in NxN boxes, clipping heavily
    55// linear interpolation to generate full-scale model
    6 bool psphotSubtractBackgroundReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index)
     6bool psphotSubtractBackgroundReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe)
    77{
    88    bool status = true;
    … …  
    2525    pmReadout *model = pmFPAviewThisReadout (view, modelFile->fpa);
    2626    assert (model);
    27 
    28     // select the appropriate recipe information
    29     psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
    30     assert (recipe);
    3127
    3228    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    … …  
    114110    // the pmReadout selected in this function are all view on entries in config->files
    115111
     112    // display the backsub and backgnd images
     113    // move this inthe the subtract background loop
     114    psphotVisualShowBackground (config, view, readout);
     115
    116116    npass ++;
    117117    return true;
    … …  
    122122    bool status = false;
    123123
     124    // select the appropriate recipe information
     125    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     126    psAssert (recipe, "missing recipe?");
     127
    124128    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
    125129    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
    … …  
    127131    // loop over the available readouts
    128132    for (int i = 0; i < num; i++) {
    129         if (!psphotSubtractBackgroundReadout (config, view, "PSPHOT.INPUT", i)) {
     133        if (!psphotSubtractBackgroundReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
    130134            psError (PSPHOT_ERR_CONFIG, false, "failed to subtract background for PSPHOT.INPUT entry %d", i);
    131135            return false;
Note: See TracChangeset for help on using the changeset viewer.