IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 36119


Ignore:
Timestamp:
Sep 14, 2013, 6:11:12 AM (13 years ago)
Author:
eugene
Message:

major rework of psphotStack sequence : ignore the CNV images, use the RAW image for detection and all analysis, measure radial apertures on RAW (copied to OUT), then smooth the RAW to the target fwhm, measure radial apertures, etc

Location:
trunk/psphot/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotStackImageLoop.c

    r36086 r36119  
    2222
    2323    pmFPAview *view = pmFPAviewAlloc (0);
    24     pmFPAfile *inputRaw = psMetadataLookupPtr (&status, config->files, "PSPHOT.STACK.INPUT.RAW");
    25     pmFPAfile *inputCnv = psMetadataLookupPtr (&status, config->files, "PSPHOT.STACK.INPUT.CNV");
    26     pmFPAfile *input = inputRaw ? inputRaw : inputCnv;
     24    pmFPAfile *input = psMetadataLookupPtr (&status, config->files, "PSPHOT.STACK.INPUT.RAW");
    2725
    2826    if (!input) {
    … …  
    3230
    3331    psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE);
    34     bool useRaw = psMetadataLookupBool (NULL, recipe, "PSPHOT.STACK.USE.RAW");
    35     if (useRaw && inputRaw == NULL) {
    36         psLogMsg ("psphot", 1, "PSPHOT.STACK.USE.RAW set but no raw input.");
    37         useRaw = false;
    38     }
    39 
    40     bool radial_apertures = psMetadataLookupBool(NULL, recipe, "RADIAL_APERTURES");
    41     bool match_psfs = radial_apertures;
    42     bool needConvolved = radial_apertures || !useRaw;
    43     if (!needConvolved) {
    44         pmFPAfileActivate (config->files, false, "PSPHOT.STACK.INPUT.CNV");
    45         pmFPAfileActivate (config->files, false, "PSPHOT.STACK.MASK.CNV");
    46         pmFPAfileActivate (config->files, false, "PSPHOT.STACK.VARIANCE.CNV");
    47         pmFPAfileActivate (config->files, false, "PSPHOT.STACK.PSF.CNV");
    48     }
    4932
    5033    // just load the full set of images up front except for EXPNUM which we defer
    5134    pmFPAfileActivate (config->files, false, "PSPHOT.STACK.EXPNUM.RAW");
    52     pmFPAfileActivate (config->files, false, "PSPHOT.STACK.EXPNUM.CNV");
    5335    if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for fpa in psphot.");
    5436
    … …  
    7153                psMemDump("load");
    7254
    73                 if (match_psfs) {
    74                     // Generate the 1st PSF-matched image set (larger target PSFs are generated by smoothing this image)
    75                     if (!psphotStackMatchPSFs (config, view)) {
    76                         psError(psErrorCodeLast(), false, "failure in psphotStackMatchPSFs for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout);
    77                         psFree (view);
    78                         return false;
    79                     }
    80                 } else {
    81                     if (!psphotStackAllocateOutput (config, view, recipe)) {
    82                         psError(psErrorCodeLast(), false, "failure in psphotStackAllocateOutput for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout);
    83                         psFree (view);
    84                         return false;
    85                     }
     55                if (!psphotStackAllocateOutput (config, view, recipe)) {
     56                    psError(psErrorCodeLast(), false, "failure in psphotStackAllocateOutput for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout);
     57                    psFree (view);
     58                    return false;
    8659                }
    8760                psMemDump("stackmatch");
    … …  
    12699
    127100    // Load the appropriate EXPNUM image
    128     pmFPAfileActivate (config->files, true, useRaw ? "PSPHOT.STACK.EXPNUM.RAW" : "PSPHOT.STACK.EXPNUM.CNV");
     101    pmFPAfileActivate (config->files, true, "PSPHOT.STACK.EXPNUM.RAW");
    129102
    130103    if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for fpa EXPNUM in psphot.");
    … …  
    146119                if (! readout->data_exists) { continue; }
    147120
    148                 if (!psphotSetNFrames (config, view, useRaw ? inputRaw->name : inputCnv->name)) ESCAPE ("failed to setNFrames.");
     121                if (!psphotSetNFrames (config, view, input->name)) ESCAPE ("failed to setNFrames.");
    149122            }
    150123        }
    … …  
    167140bool UpdateHeadersForReadout (pmConfig *config, pmFPAview *view) {
    168141
    169     psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE);
    170     bool useRaw = psMetadataLookupBool (NULL, recipe, "PSPHOT.STACK.USE.RAW");
    171 
    172142    int num = psphotFileruleCount(config, "PSPHOT.INPUT");
    173143
    … …  
    179149        psAssert (output, "missing file?");
    180150
    181         pmFPAfile *inputRaw = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
    182         pmFPAfile *inputCnv = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.CNV", i); // File of interest
    183         pmFPAfile *input = useRaw ? inputRaw : inputCnv;
     151        pmFPAfile *input = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
    184152        psAssert (input, "missing input file");
    185153
    … …  
    205173        psAssert (output, "missing file?");
    206174
    207         pmFPAfile *inputRaw = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
    208         pmFPAfile *inputCnv = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.CNV", i); // File of interest
    209         pmFPAfile *input = inputRaw ? inputRaw : inputCnv;
     175        pmFPAfile *input = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
    210176        psAssert (input, "missing input file");
    211177
    … …  
    242208        psAssert (output, "missing file?");
    243209
    244         pmFPAfile *inputRaw = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
    245         pmFPAfile *inputCnv = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.CNV", i); // File of interest
    246         pmFPAfile *input = inputRaw ? inputRaw : inputCnv;
     210        pmFPAfile *input = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
    247211        psAssert (input, "missing input file");
    248212

  • trunk/psphot/src/psphotStackReadout.c

    r34721 r36119  
    11# include "psphotInternal.h"
    22
     3static bool psphotStackMatchPSFsetup (pmConfig *config, const pmFPAview *view, const char *filerule, const char *fPSF);
     4static bool psphotStackMatchPSFsetupReadout (pmConfig *config, const pmFPAview *view, const char *filerule, const char *fPSF, int index);
    35static bool psphotStackLoadWCS(pmConfig *config, const pmFPAview *view, const char *filerule);
    46static void logMemStats(const char *heading);
    57
    6 // we have 3 possible real filesets:
     8// relevant filesets:
    79# define STACK_RAW "PSPHOT.STACK.INPUT.RAW"
    8 # define STACK_CNV "PSPHOT.STACK.INPUT.CNV"
    9 # define STACK_OUT "PSPHOT.STACK.OUTPUT.IMAGE"  /* the psf-matched image */
    10 
    11 // we have 3 files on which we operate:
    12 // DET (detection image)       : nominally RAW (optionally CNV?)
    13 // SRC (source analysis image) : nominally CNV (optionally RAW)
    14 // OUT (psf-matched images)    : always OUT
     10# define STACK_OUT "PSPHOT.STACK.OUTPUT.IMAGE"
     11
     12// XXX STACK_OUT currently is a copy of STACK_RAW, but should be a pointer to is as in psphot (single)
     13
     14// TEST CODE, can be removed
     15bool psphotDumpImages (pmConfig *config, const pmFPAview *view, const char *filerule, char *base) {
     16
     17    // XXX do nothing
     18    return true;
     19
     20    int num = psphotFileruleCount(config, "PSPHOT.INPUT");
     21
     22    for (int i = 0; i < num; i++) {
     23        // find the currently selected readout
     24        pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest
     25        psAssert (file, "missing file?");
     26
     27        pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     28        psAssert (readout, "missing readout?");
     29
     30        char line[256];
     31        snprintf (line, 256, "%s.%d.im.fits", base, i);
     32        psphotSaveImage (NULL, readout->image, line);
     33
     34        snprintf (line, 256, "%s.%d.wt.fits", base, i);
     35        psphotSaveImage (NULL, readout->variance, line);
     36
     37        snprintf (line, 256, "%s.%d.mk.fits", base, i);
     38        psphotSaveImage (NULL, readout->mask, line);
     39    }
     40    return true;
     41}
    1542
    1643bool psphotStackVisualFilerule(pmConfig *config, const pmFPAview *view, const char *filerule) {
    … …  
    6794    psAssert (breakPt, "configuration error: set BREAK_POINT");
    6895
    69     // we have 3 relevant files: RAW (unconvolved), CNV (convolved stack), OUT (psf-matched stack)
    70     // select which image (RAW or CNV) is used for analysis (RAW always used for detection)
    71     bool useRaw = psMetadataLookupBool (NULL, recipe, "PSPHOT.STACK.USE.RAW");
    72     char *STACK_SRC = useRaw ? STACK_RAW : STACK_CNV;
    73     char *STACK_DET = STACK_RAW;
    74 
    7596    // load WCS
    76     if (!psphotStackLoadWCS(config, view, STACK_SRC)) {
    77         psError (PSPHOT_ERR_CONFIG, false, "trouble loading WCS for %s", STACK_SRC);
     97    if (!psphotStackLoadWCS(config, view, STACK_RAW)) {
     98        psError (PSPHOT_ERR_CONFIG, false, "trouble loading WCS for %s", STACK_RAW);
    7899        return false;
    79100    }
    80101
    81102    // set the photcode for each image
    82     if (!psphotAddPhotcode (config, view, STACK_SRC)) {
     103    if (!psphotAddPhotcode (config, view, STACK_RAW)) {
    83104        psError (PSPHOT_ERR_CONFIG, false, "trouble defining the photcode");
    84105        return false;
    … …  
    86107
    87108    // Generate the mask and weight images (if not supplied) and set mask bits.
    88     // This also insures that all invalid pixels are masked (this is done for STACK_CNV in psphotStackMatchPSFs)
    89     if (!psphotSetMaskAndVariance (config, view, STACK_DET)) {
    90         return psphotReadoutCleanup (config, view, STACK_SRC);
    91     }
    92     if (!psphotSetMaskAndVariance (config, view, STACK_OUT)) {
    93         return psphotReadoutCleanup (config, view, STACK_SRC);
     109    if (!psphotSetMaskAndVariance (config, view, STACK_RAW)) {
     110        return psphotReadoutCleanup (config, view, STACK_RAW);
    94111    }
    95112    if (!strcasecmp (breakPt, "NOTHING")) {
    96         return psphotReadoutCleanup (config, view, STACK_SRC);
     113        return psphotReadoutCleanup (config, view, STACK_RAW);
    97114    }
    98115
    99116    // generate a background model (median, smoothed image)
    100     if (!psphotModelBackground (config, view, STACK_DET)) {
    101         return psphotReadoutCleanup (config, view, STACK_SRC);
    102     }
    103     if (!psphotSubtractBackground (config, view, STACK_DET)) {
    104         return psphotReadoutCleanup (config, view, STACK_SRC);
    105     }
    106     if (strcmp(STACK_SRC, STACK_DET)) {
    107 #define MODEL_BACKGROUND_SRC 1
    108 #ifdef MODEL_BACKGROUND_SRC
    109         // work around the fact that the background levels on the convolved
    110         // and unconvolved stacks can be different
    111         if (!psphotModelBackground (config, view, STACK_SRC)) {
    112             return psphotReadoutCleanup (config, view, STACK_SRC);
    113         }
    114 #endif
    115         if (!psphotSubtractBackground (config, view, STACK_SRC)) {
    116             return psphotReadoutCleanup (config, view, STACK_SRC);
    117         }
     117    if (!psphotModelBackground (config, view, STACK_RAW)) {
     118        return psphotReadoutCleanup (config, view, STACK_RAW);
     119    }
     120    if (!psphotSubtractBackground (config, view, STACK_RAW)) {
     121        return psphotReadoutCleanup (config, view, STACK_RAW);
    118122    }
    119123    if (!strcasecmp (breakPt, "BACKMDL")) {
    120         return psphotReadoutCleanup (config, view, STACK_SRC);
     124        return psphotReadoutCleanup (config, view, STACK_RAW);
    121125    }
    122126
    123127#ifdef MAKE_CHISQ_IMAGE
    124128    // also make the chisq detection image
    125     if (!psphotStackChisqImage(config, view, STACK_DET, STACK_SRC)) {
     129    if (!psphotStackChisqImage(config, view, STACK_RAW, STACK_RAW)) {
    126130        psError (PSPHOT_ERR_UNKNOWN, false, "failure to generate chisq image");
    127         return psphotReadoutCleanup (config, view, STACK_SRC);
     131        return psphotReadoutCleanup (config, view, STACK_RAW);
     132    }
     133    if (!strcasecmp (breakPt, "CHISQ")) {
     134        return psphotReadoutCleanup (config, view, STACK_RAW);
    128135    }
    129136#endif
    130     if (!strcasecmp (breakPt, "CHISQ")) {
    131         return psphotReadoutCleanup (config, view, STACK_SRC);
    132     }
    133137
    134138    // find the detections (by peak and/or footprint) in the image.
    135139    // This finds the detections on Chisq image as well as the individuals
    136     if (!psphotFindDetections (config, view, STACK_DET, true)) { // pass 1
    137         // this only happens if we had an error in psphotFindDetections
     140    if (!psphotFindDetections (config, view, STACK_RAW, true)) { // pass 1
    138141        psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
    139         return psphotReadoutCleanup (config, view, STACK_SRC);
    140     }
    141 
    142     // If DET and SRC are different images, copy the detections from DET to SRC.  This 'copy'
    143     // is just a copy of the container pointer; the sources on both DET and SRC are the same
    144     // memory objects
    145     if (strcmp(STACK_SRC, STACK_DET)) {
    146         if (!psphotCopySources (config, view, STACK_SRC, STACK_DET)) {
    147             psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
    148             return psphotReadoutCleanup (config, view, STACK_SRC);
    149         }
     142        return psphotReadoutCleanup (config, view, STACK_RAW);
    150143    }
    151144
    152145    // construct sources and measure basic stats (saved on detections->newSources)
    153     if (!psphotSourceStats (config, view, STACK_SRC, true)) { // pass 1
     146    if (!psphotSourceStats (config, view, STACK_RAW, true)) { // pass 1
    154147        psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
    155         return psphotReadoutCleanup (config, view, STACK_SRC);
     148        return psphotReadoutCleanup (config, view, STACK_RAW);
    156149    }
    157150    if (!strcasecmp (breakPt, "PEAKS")) {
    158         return psphotReadoutCleanup (config, view, STACK_SRC);
    159     }
    160     // psphotDumpTest (config, view, STACK_SRC);
     151        return psphotReadoutCleanup (config, view, STACK_RAW);
     152    }
     153    // psphotDumpTest (config, view, STACK_RAW);
    161154    psMemDump("sourcestats");
    162155    logMemStats("sourcestats");
    … …  
    164157    // classify sources based on moments, brightness
    165158    // only run this on detections from the input images, not chisq image
    166     if (!psphotRoughClass (config, view, STACK_SRC)) {
     159    if (!psphotRoughClass (config, view, STACK_RAW)) {
    167160        psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough classifications");
    168         return psphotReadoutCleanup (config, view, STACK_SRC);
    169     }
    170 
    171     // If DET and SRC are different images, subtract radial profiles for the convolved
    172     // image first.  The profiles found on the convolved image will be replaced by those
    173     // found for the unconvolved image.  Downstream, in psphotFindDetections, we will
    174     // replace the the profiles (and re-subtract them) for the detection image, so we want
    175     // to keep those versions of the profiles on the sources.
    176     if (strcmp(STACK_SRC, STACK_DET)) {
    177       // find and subtract radial profile models for saturated stars (XXX change name eventually)
    178       if (!psphotDeblendSatstars (config, view, STACK_SRC)) {
     161        return psphotReadoutCleanup (config, view, STACK_RAW);
     162    }
     163
     164    // find and subtract radial profile models for saturated stars (XXX change name eventually)
     165    if (!psphotDeblendSatstars (config, view, STACK_RAW)) {
    179166        psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis");
    180         return psphotReadoutCleanup (config, view, STACK_SRC);
    181       }
    182     }
    183     // find and subtract radial profile models for saturated stars (XXX change name eventually)
    184     if (!psphotDeblendSatstars (config, view, STACK_DET)) {
    185         psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis");
    186         return psphotReadoutCleanup (config, view, STACK_SRC);
     167        return psphotReadoutCleanup (config, view, STACK_RAW);
    187168    }
    188169
    189170    // if we were not supplied a PSF model, determine the IQ stats here (detections->newSources)
    190171    // only run this on detections from the input images, not chisq image
    191     if (!psphotImageQuality (config, view, STACK_SRC)) { // pass 1
     172    if (!psphotImageQuality (config, view, STACK_RAW)) { // pass 1
    192173        psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality");
    193         return psphotReadoutCleanup (config, view, STACK_SRC);
     174        return psphotReadoutCleanup (config, view, STACK_RAW);
    194175    }
    195176    if (!strcasecmp (breakPt, "MOMENTS")) {
    196         return psphotReadoutCleanup (config, view, STACK_SRC);
     177        return psphotReadoutCleanup (config, view, STACK_RAW);
    197178    }
    198179
    199180    // use bright stellar objects to measure PSF
    200     if (!psphotChoosePSF (config, view, STACK_SRC, true)) { // pass 1
     181    if (!psphotChoosePSF (config, view, STACK_RAW, true)) { // pass 1
    201182        psLogMsg ("psphot", 3, "failure to construct a psf model");
    202         return psphotReadoutCleanup (config, view, STACK_SRC);
     183        return psphotReadoutCleanup (config, view, STACK_RAW);
    203184    }
    204185    if (!strcasecmp (breakPt, "PSFMODEL")) {
    205         return psphotReadoutCleanup (config, view, STACK_SRC);
     186        return psphotReadoutCleanup (config, view, STACK_RAW);
    206187    }
    207188
    208189    // merge the newly selected sources into the existing list
    209190    // NOTE: merge OLD and NEW
    210     psphotMergeSources (config, view, STACK_SRC);
     191    psphotMergeSources (config, view, STACK_RAW);
    211192
    212193    // Construct an initial model for each object, set the radius to fitRadius, set circular
    213194    // fit mask.  NOTE: only applied to sources without guess models
    214     psphotGuessModels (config, view, STACK_SRC);
     195    psphotGuessModels (config, view, STACK_RAW);
    215196
    216197    // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
    217     psphotFitSourcesLinear (config, view, STACK_SRC, false, false);
    218     psphotStackVisualFilerule(config, view, STACK_SRC);
     198    psphotFitSourcesLinear (config, view, STACK_RAW, false, false);
     199    psphotStackVisualFilerule(config, view, STACK_RAW);
    219200
    220201    // measure the radial profiles to the sky
    221     psphotRadialProfileWings (config, view, STACK_SRC);
     202    psphotRadialProfileWings (config, view, STACK_RAW);
    222203
    223204    // re-measure the kron mags with models subtracted.  this pass starts with a circular
    … …  
    225206    // but iterates to an appropriately larger size
    226207    logMemStats("before.kron.1");
    227     psphotKronIterate(config, view, STACK_SRC, 1);
     208    psphotKronIterate(config, view, STACK_RAW, 1);
    228209    logMemStats("after.kron.1");
    229210       
    230211    // identify CRs and extended sources
    231     psphotSourceSize (config, view, STACK_SRC, true);
     212    psphotSourceSize (config, view, STACK_RAW, true);
    232213
    233214    // non-linear PSF and EXT fit to brighter sources
    234215    // replace model flux, adjust mask as needed, fit, subtract the models (full stamp)
    235     psphotBlendFit (config, view, STACK_SRC); // pass 1 (detections->allSources)
     216    psphotBlendFit (config, view, STACK_RAW); // pass 1 (detections->allSources)
    236217
    237218    // replace all sources (do NOT ignore subtraction state)
    238     psphotReplaceAllSources (config, view, STACK_SRC, false); // pass 1 (detections->allSources)
     219    psphotReplaceAllSources (config, view, STACK_RAW, false); // pass 1 (detections->allSources)
    239220
    240221    logMemStats("pass1");
    … …  
    245226    // linear fit to include all sources (subtract again)
    246227    // NOTE : apply to ALL sources (extended + psf)
    247     // NOTE 2 : this function subtracts the models from the given filerule (SRC), not DET
    248     psphotFitSourcesLinear (config, view, STACK_SRC, true, false); // pass 2 (detections->allSources)
     228    // NOTE 2 : this function subtracts the models from the given filerule
     229    psphotFitSourcesLinear (config, view, STACK_RAW, true, false); // pass 2 (detections->allSources)
    249230
    250231    // NOTE: possibly re-measure background model here with objects subtracted / or masked
    … …  
    252233    // NOTE: this block performs the 2nd pass low-significance PSF detection stage
    253234    {
    254         // if DET and SRC are different images, generate children sources for all sources in
    255         // the SRC image.  This operation replaces the existing DETECTION container on DET
    256         // which is currently a view to the one on SRC).  children sources go to
    257         // det->allSources
    258         if (strcmp(STACK_SRC, STACK_DET)) {
    259             psphotSourceChildren (config, view, STACK_DET, STACK_SRC);
    260 
    261             //  subtract all sources from DET (this will subtract using the psf model for SRC, which
    262             //  will somewhat oversubtract the sources -- this is OK
    263             psphotRemoveAllSources (config, view, STACK_DET, false); // do not ignore subtraction state for sources
    264         }
    265 
    266235        // add noise for subtracted objects
    267         psphotAddNoise (config, view, STACK_DET); // pass 1 (detections->allSources)
     236        psphotAddNoise (config, view, STACK_RAW); // pass 1 (detections->allSources)
    268237
    269238        // find fainter sources
    270239        // NOTE: finds new peaks and new footprints, OLD and FULL set are saved on detections
    271         psphotFindDetections (config, view, STACK_DET, false); // pass 2 (detections->peaks, detections->footprints)
     240        psphotFindDetections (config, view, STACK_RAW, false); // pass 2 (detections->peaks, detections->footprints)
    272241
    273242        // remove noise for subtracted objects (ie, return to normal noise level)
    … …  
    276245        bool footprintsUseUnsubtracted = psMetadataLookupBool(NULL, recipe, "FOOTPRINT_USE_UNSUBTRACTED");
    277246        if (!footprintsUseUnsubtracted) {
    278             psphotSubNoise (config, view, STACK_DET); // pass 1 (detections->allSources)
     247            psphotSubNoise (config, view, STACK_RAW); // pass 1 (detections->allSources)
    279248        }
    280 
    281         // if DET and SRC are different images, copy the detections from DET to SRC
    282         // (this operation just ensures the metadata container has a view on SRC as well
    283         if (strcmp(STACK_SRC, STACK_DET)) {
    284             // replace all sources in DET
    285             psphotReplaceAllSources (config, view, STACK_DET, false); // ignore subtraction state for sources
    286 
    287             // copy the newly detected peaks from DET to SRC so SourceStats below can operate on them
    288             if (!psphotCopyPeaks (config, view, STACK_SRC, STACK_DET)) {
    289                 psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
    290                 return psphotReadoutCleanup (config, view, STACK_SRC);
    291             }
    292         }
    293249
    294250        // define new sources based on only the new peaks & measure moments
    295251        // NOTE: new sources are saved on detections->newSources
    296         psphotSourceStats (config, view, STACK_SRC, false); // pass 2 (detections->newSources)
     252        psphotSourceStats (config, view, STACK_RAW, false); // pass 2 (detections->newSources)
    297253
    298254        // set source type
    299255        // NOTE: apply only to detections->newSources
    300         if (!psphotRoughClass (config, view, STACK_SRC)) { // pass 2 (detections->newSources)
     256        if (!psphotRoughClass (config, view, STACK_RAW)) { // pass 2 (detections->newSources)
    301257            psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
    302             return psphotReadoutCleanup (config, view, STACK_SRC);
     258            return psphotReadoutCleanup (config, view, STACK_RAW);
    303259        }
    304260
    305261        // replace all sources so fit below applies to all at once
    306262        // NOTE: apply only to OLD sources (which have been subtracted)
    307         psphotReplaceAllSources (config, view, STACK_SRC, false); // pass 2
     263        psphotReplaceAllSources (config, view, STACK_RAW, false); // pass 2
    308264
    309265        // merge the newly selected sources into the existing list
    310266        // NOTE: merge OLD and NEW
    311267        // XXX check on free of sources...
    312         psphotMergeSources (config, view, STACK_SRC); // (detections->newSources + detections->allSources -> detections->allSources)
     268        psphotMergeSources (config, view, STACK_RAW); // (detections->newSources + detections->allSources -> detections->allSources)
    313269
    314270        // Construct an initial model for each object, set the radius to fitRadius, set circular
    315271        // fit mask.  NOTE: only applied to sources without guess models
    316         psphotGuessModels (config, view, STACK_SRC);
     272        psphotGuessModels (config, view, STACK_RAW);
    317273    }
    318274
    … …  
    325281    if (splitLinearFit) {
    326282        psLogMsg ("psphot", 3, "splitting fit of detected and matched soures\n");
    327         // Fit the detected sources separately from matched that wea are about to create.
     283        // Fit the detected sources separately from matched ones that we are about to create.
    328284        // NOTE: apply to ALL sources but only include sources with postitive flux in the fit
    329         psphotFitSourcesLinear (config, view, STACK_SRC, true, true); // pass 3 (detections->allSources)
     285        psphotFitSourcesLinear (config, view, STACK_RAW, true, true); // pass 3 (detections->allSources)
    330286    }
    331287
    … …  
    335291    // this just match the detections for the chisq image, and not bother measuring the source
    336292    // stats in that case...?
    337     objects = psphotMatchSources (config, view, STACK_SRC);
     293    objects = psphotMatchSources (config, view, STACK_RAW);
    338294    psMemDump("matchsources");
    339295
    … …  
    344300    // Construct an initial model for each object, set the radius to fitRadius, set circular
    345301    // fit mask.  NOTE: only applied to sources without guess models
    346     psphotGuessModels (config, view, STACK_SRC);
     302    psphotGuessModels (config, view, STACK_RAW);
    347303
    348304    psphotStackObjectsUnifyPosition (objects);
    349305
    350     psphotStackObjectsSelectForAnalysis (config, view, STACK_SRC, objects);
     306    psphotStackObjectsSelectForAnalysis (config, view, STACK_RAW, objects);
    351307
    352308    // final linear fit. NOTE: if splitLinearFit is true above, this pass will only fit
    353309    // the unsubtracted (matched) sources (the sources that we fit above are subtracted)
    354     psphotFitSourcesLinear (config, view, STACK_SRC, true, false); // pass 4 (detections->allSources)
     310    psphotFitSourcesLinear (config, view, STACK_RAW, true, false); // pass 4 (detections->allSources)
    355311
    356312    // measure the radial profiles to the sky (only measures new objects)
    357     psphotRadialProfileWings (config, view, STACK_SRC);
     313    psphotRadialProfileWings (config, view, STACK_RAW);
    358314
    359315    // re-measure the kron mags with models subtracted
    360316    // psphotKronMasked(config, view, STACK_SRC);
    361317    logMemStats("before.kron.2");
    362     psphotKronIterate(config, view, STACK_SRC, 2);
     318    psphotKronIterate(config, view, STACK_RAW, 2);
    363319    logMemStats("after.kron.2");
    364320
    365321    // measure source size for the remaining sources
    366322    // NOTE: applies only to NEW (unmeasured) sources
    367     psphotSourceSize (config, view, STACK_SRC, false); // pass 2 (detections->allSources)
     323    psphotSourceSize (config, view, STACK_RAW, false); // pass 2 (detections->allSources)
    368324
    369325    psMemDump("psfstats");
    … …  
    371327    // drop matched sources without any useful measurements and set kron radii for the ones
    372328    // we decide to keep
    373     psphotFilterMatchedSources (config, view, STACK_SRC, objects);
     329    psphotFilterMatchedSources (config, view, STACK_RAW, objects);
    374330
    375331    // measure kron fluxes for the matched sources only
    376     psphotKronIterate(config, view, STACK_SRC, 3);
     332    psphotKronIterate(config, view, STACK_RAW, 3);
    377333
    378334    // measure elliptical apertures, petrosians (objects sorted by S/N)
    379335    // psphotExtendedSourceAnalysisByObject (config, objects, view, STACK_SRC); // pass 1 (detections->allSources)
    380     psphotExtendedSourceAnalysis (config, view, STACK_SRC); // pass 1 (detections->allSources)
     336    psphotExtendedSourceAnalysis (config, view, STACK_RAW); // pass 1 (detections->allSources)
    381337
    382338    // measure non-linear extended source models (exponential, deVaucouleur, Sersic) (sources sorted by S/N)
    383     psphotExtendedSourceFits (config, view, STACK_SRC); // pass 1 (detections->allSources)
     339    psphotExtendedSourceFits (config, view, STACK_RAW); // pass 1 (detections->allSources)
    384340
    385341    // create source children for the OUT filerule (for radial aperture photometry and output)
    386     psArray *objectsOut = psphotSourceChildrenByObject (config, view, STACK_OUT, objects);
     342    // NOTE: The new source children have image arrays pointing to the readout associated with
     343    // STACK_OUT.  in psphotStackMatchPSFsetup, we copy the current pixel values from RAW to OUT,
     344    // but keep the pointers the same so we do not break these source image references
     345    // XXX NOTE : if we use the pre-20130914 psphotStackReadout code, we need to use 'false' for the
     346    // sourcesSubtracted argument
     347    psArray *objectsOut = psphotSourceChildrenByObject (config, view, STACK_OUT, objects, true);
    387348    if (!objectsOut) {
    388349        psFree(objects);
    389350        psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
    390         return psphotReadoutCleanup (config, view, STACK_SRC);
     351        return psphotReadoutCleanup (config, view, STACK_RAW);
    391352    }
    392353
    … …  
    396357        // this forces photometry on the undetected sources from other images
    397358
    398         // NOTE: we always do the radial apertures analysis on the convolved image since
    399         // those are the ones that are psf matched and are the source of STACK_OUT's pixels
    400         // XXX: Actually if PSPHOT.STACK.MATCH.PSF.SOURCE were set to RAW this wouldn't be true.
    401         // but in that case we don't get past the psf matching step because there is no
    402         // target psf for the RAW inputs
    403 
    404         // If useRaw copy the sources to the convolved readout
    405         if (strcmp(STACK_SRC, STACK_CNV)) {
    406             if (!psphotCopySources (config, view, STACK_CNV, STACK_SRC)) {
    407                 psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
    408                 return psphotReadoutCleanup (config, view, STACK_SRC);
    409             }
    410         }
    411         // mark any inputs that we want to skip the matched apertures for
    412         psphotStackSetInputsToSkip(config, view, STACK_CNV, true);
    413         psphotStackSetInputsToSkip(config, view, STACK_OUT, true);
    414         psphotRadialApertures (config, view, STACK_CNV, 0); // entry 0 == unmatched
    415         psMemDump("extmeas");
    416 
     359        // set up the FWHM vector
     360        psphotStackMatchPSFsetup (config, view, STACK_OUT, STACK_RAW);
     361        psphotDumpImages (config, view, STACK_RAW, "raw.t0");
     362        psphotDumpImages (config, view, STACK_OUT, "out.t0");
    417363
    418364        int nRadialEntries = psphotStackMatchPSFsEntries(config, view, STACK_OUT);
    419         for (int entry = 1; entry < nRadialEntries; entry++) {
     365
     366        for (int entry = 0; entry < nRadialEntries; entry++) {
    420367            // NOTE: entry 0 is the unmatched image set
    421368
    422             // re-measure the PSF for the smoothed image (using entries in 'allSources')
    423             psphotChoosePSF (config, view, STACK_OUT, false);
    424 
    425             // this is necessary to update the models based on the new PSF
    426             psphotResetModels (config, view, STACK_OUT);
    427 
    428             // this is necessary to get the right normalization for the new models
    429             psphotFitSourcesLinear (config, view, STACK_OUT, false, false);
     369            char line[256];
    430370
    431371            // measure circular, radial apertures (objects sorted by S/N)
    432372            psphotRadialApertures (config, view, STACK_OUT, entry);
     373            snprintf (line, 256, "%s.%d", "out.t1", entry);
     374            psphotDumpImages (config, view, STACK_OUT, line);
    433375
    434376            // replace the flux in the image so it is returned to its original state
    435377            psphotReplaceAllSources (config, view, STACK_OUT, false);
    436 
    437             // smooth to the next FWHM, or set 'smoothAgain' to false if no more
    438             psphotStackMatchPSFsNext(config, view, STACK_OUT, entry);
    439             psMemDump("matched");
     378            snprintf (line, 256, "%s.%d", "out.t2", entry);
     379            psphotDumpImages (config, view, STACK_OUT, line);
     380
     381            if (entry < nRadialEntries - 1) {
     382                // smooth to the next FWHM
     383                // this function does nothing if the targetFWHM is smaller than the currentFWHM
     384                psphotStackMatchPSFsNext(config, view, STACK_OUT, entry);
     385                snprintf (line, 256, "%s.%d", "out.t3", entry);
     386                psphotDumpImages (config, view, STACK_OUT, line);
     387                psMemDump("matched");
     388
     389                // re-measure the PSF for the smoothed image (using entries in 'allSources')
     390                psphotChoosePSF (config, view, STACK_OUT, false);
     391
     392                // this is necessary to update the models based on the new PSF
     393                psphotResetModels (config, view, STACK_OUT);
     394
     395                // this is necessary to get the right normalization for the new models
     396                // and to subtract the sources
     397                psphotFitSourcesLinear (config, view, STACK_OUT, false, false);
     398                snprintf (line, 256, "%s.%d", "out.t4", entry);
     399                psphotDumpImages (config, view, STACK_OUT, line);
     400            }
    440401        }
    441402    }
    442     psphotStackSetInputsToSkip(config, view, STACK_CNV, false);
    443     psphotStackSetInputsToSkip(config, view, STACK_OUT, false);
    444403
    445404    // measure aperture photometry corrections
    446     if (!psphotApResid (config, view, STACK_SRC)) {
     405    if (!psphotApResid (config, view, STACK_RAW)) {
    447406        psFree (objects);
    448407        psFree (objectsOut);
    449408        psLogMsg ("psphot", 3, "failed on psphotApResid");
    450         return psphotReadoutCleanup (config, view, STACK_SRC);
     409        return psphotReadoutCleanup (config, view, STACK_RAW);
    451410    }
    452411
    453412    // calculate source magnitudes
    454     psphotMagnitudes(config, view, STACK_SRC);
    455 
    456     if (!useRaw) {
    457         // psphotEfficiency wants to have the PSF of the image, but since we are measuring on
    458         // the convolved images we need to generate PSFs for the DET images
    459         if (!psphotChoosePSF (config, view, STACK_DET, false)) {
    460             psLogMsg ("psphot", 3, "failure to construct a psf model for raw input");
    461             return psphotReadoutCleanup (config, view, STACK_DET);
    462         }
    463     }
    464     if (!psphotEfficiency(config, view, STACK_DET)) {
     413    psphotMagnitudes(config, view, STACK_RAW);
     414
     415    if (!psphotEfficiency(config, view, STACK_RAW)) {
    465416        psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources");
    466417        psErrorClear();
    467418    }
    468     psphotCopyEfficiency (config, view, STACK_OUT, STACK_DET);
     419    psphotCopyEfficiency (config, view, STACK_OUT, STACK_RAW);
    469420
    470421    logMemStats("final");
    471422#if (1)
    472     psphotSourceMemory(config, view, STACK_SRC);
     423    psphotSourceMemory(config, view, STACK_RAW);
    473424    psphotSourceMemory(config, view, STACK_OUT);
    474425#endif
    475426
    476     // replace failed sources?
    477     // psphotReplaceUnfitSources (sources);
    478 
    479427    // replace background in residual image
    480     psphotSkyReplace (config, view, STACK_DET);
     428    psphotSkyReplace (config, view, STACK_RAW);
    481429
    482430    // drop the references to the image pixels held by each source
     431    psphotSourceFreePixels (config, view, STACK_RAW);
    483432    psphotSourceFreePixels (config, view, STACK_OUT);
    484     psphotSourceFreePixels (config, view, STACK_SRC);
    485433
    486434#ifdef MAKE_CHISQ_IMAGE
    487435    // remove chisq image from config->file:PSPHOT.INPUT
    488     psphotStackRemoveChisqFromInputs(config, STACK_DET);
    489     if (strcmp(STACK_SRC, STACK_DET)) {
    490         psphotStackRemoveChisqFromInputs(config, STACK_SRC);
    491     }
     436    psphotStackRemoveChisqFromInputs(config, STACK_RAW);
    492437#endif
    493438
    … …  
    496441
    497442    // create the exported-metadata and free local data
    498     return psphotReadoutCleanup (config, view, STACK_SRC);
     443    return psphotReadoutCleanup (config, view, STACK_RAW);
    499444}
    500445
    … …  
    552497}
    553498
    554 
    555 
    556499/* here is the process:
    557500
    558  * we have three(*) images:
    559  * RAW : unconvolved image stack
    560  * CNV : input convolved image
     501 * we have two image sets:
     502 * RAW : unconvolved image stacks
    561503
    562504 * OUT : psf-matched output image (there may be more than one of
    … …  
    649591
    650592   */
     593
     594
     595// generate a vector fwhmValues where the first has the fwhm of the raw image and the
     596// successive entries have the target values
     597
     598bool psphotStackMatchPSFsetup (pmConfig *config, const pmFPAview *view, const char *filerule, const char *filerulePSF) {
     599
     600    bool status;
     601
     602    int num = psphotFileruleCount(config, filerule);
     603
     604    // skip the chisq image (optionally?)
     605    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
     606    if (!status) chisqNum = -1;
     607
     608    // loop over the available readouts
     609    for (int i = 0; i < num; i++) {
     610        if (i == chisqNum) continue; // skip chisq image
     611
     612        if (!psphotStackMatchPSFsetupReadout (config, view, filerule, filerulePSF, i)) {
     613            psError (PSPHOT_ERR_CONFIG, false, "failed to define target PSF sizes");
     614            return false;
     615        }
     616    }
     617
     618    return true;
     619}
     620
     621float psphotPSFseeing (pmPSF *psf, pmReadout *readout, int index);
     622
     623// copy the pixels from RAW to OUT (
     624
     625bool psphotStackMatchPSFsetupReadout (pmConfig *config, const pmFPAview *view, const char *fileruleOut, const char *fileruleRaw, int index) {
     626
     627    bool status;
     628
     629    // select the appropriate recipe information
     630    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     631
     632    // find the currently selected readout
     633    pmFPAfile *fileOut = pmFPAfileSelectSingle(config->files, fileruleOut, index); // File of interest
     634    psAssert (fileOut, "missing file?");
     635
     636    pmReadout *readoutOut = pmFPAviewThisReadout(view, fileOut->fpa);
     637    psAssert (readoutOut, "missing readout?");
     638
     639    // find the currently selected readout
     640    pmFPAfile *fileRaw = pmFPAfileSelectSingle(config->files, fileruleRaw, index); // File of interest
     641    psAssert (fileRaw, "missing file?");
     642
     643    pmReadout *readoutRaw = pmFPAviewThisReadout(view, fileRaw->fpa);
     644    psAssert (readoutRaw, "missing readout?");
     645
     646    readoutOut->image = psImageCopy(readoutOut->image, readoutRaw->image, PS_TYPE_F32);
     647    if (readoutRaw->variance) {
     648        readoutOut->variance = psImageCopy(readoutOut->variance, readoutRaw->variance, PS_TYPE_F32);
     649    }
     650    if (readoutRaw->mask) {
     651        readoutOut->mask = psImageCopy(readoutOut->mask, readoutRaw->mask, PS_TYPE_IMAGE_MASK);
     652    }
     653
     654    // pmChip *chipRaw = pmFPAviewThisChip(view, fileRaw->fpa); // The chip holds the PSF
     655    // psAssert (chipRaw, "missing chip");
     656
     657    pmPSF *psf = psMetadataLookupPtr(&status, readoutRaw->analysis, "PSPHOT.PSF"); // PSF
     658    if (!psf) {
     659        // we should have a PSF by this point in psphot
     660        psError(PSPHOT_ERR_PROG, true, "Unable to find PSF.");
     661        return false;
     662    }
     663
     664    float fwhmRaw = psphotPSFseeing (psf, readoutRaw, index);
     665
     666    psVector *fwhmValues = psVectorAllocEmpty(10, PS_TYPE_F32);
     667    psVectorAppend(fwhmValues, fwhmRaw);
     668
     669    // is a single target FWHM specified, or a set of values?  set up the vector options->targetSeeing and the local 1st value
     670    float targetSeeing = psMetadataLookupF32 (&status, recipe, "PSPHOT.STACK.TARGET.PSF.FWHM");
     671    if (!status) {
     672        psVector *targetSeeing = psMetadataLookupVector(&status, recipe, "PSPHOT.STACK.TARGET.PSF.FWHM"); // Magnitude offsets
     673        psAssert (status, "missing psphot recipe value PSPHOT.STACK.TARGET.PSF.FWHM");
     674        for (int i = 0; i < targetSeeing->n; i++) {
     675            psVectorAppend(fwhmValues, targetSeeing->data.F32[i]);
     676        }           
     677    } else {
     678        psVectorAppend(fwhmValues, targetSeeing);
     679    }
     680
     681    psMetadataAddVector(readoutOut->analysis, PS_LIST_TAIL, "STACK.PSF.FWHM.VALUES", PS_META_REPLACE, "PSF sizes", fwhmValues);
     682    psFree (fwhmValues);
     683
     684    return true;
     685}
     686
     687float psphotPSFseeing (pmPSF *psf, pmReadout *readout, int index) {
     688
     689    psImage *image = readout->image;
     690
     691    int Nx = image->numCols;
     692    int Ny = image->numRows;
     693
     694    float sumFWHM = 0.0;                  // FWHM for image
     695    int numFWHM = 0;                      // Number of FWHM measurements
     696    for (float x = 0; x < Nx; x += 0.25*Nx) {
     697        for (float y = 0; y < Ny; y += 0.25*Ny) {
     698            float fwhm = pmPSFtoFWHM(psf, x, y);
     699            if (isfinite(fwhm)) {
     700                sumFWHM += fwhm;
     701                numFWHM++;
     702            }
     703        }
     704    }
     705    if (numFWHM == 0) {
     706        psLogMsg("ppStack", PS_LOG_INFO, "Unable to measure PSF FWHM for image %d --- rejected.", index);
     707        return NAN;
     708    }
     709
     710    float fwhm = sumFWHM / (float) numFWHM;
     711    psLogMsg ("psphotStack", PS_LOG_INFO, "Input Seeing for %d: %f\n", index, fwhm);
     712
     713    return fwhm;
     714}
Note: See TracChangeset for help on using the changeset viewer.