IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 9, 2013, 4:14:19 PM (13 years ago)
Author:
eugene
Message:

merge from trunk

Location:
branches/eam_branches/ipp-20130904/psphot
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20130904/psphot

  • branches/eam_branches/ipp-20130904/psphot/src

  • branches/eam_branches/ipp-20130904/psphot/src/psphotStackReadout.c

    r34721 r36198  
    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);
    121     }
     124        return psphotReadoutCleanup (config, view, STACK_RAW);
     125    }
     126
     127// XXX TEST for background:
     128    if (!psphotModelBackground (config, view, STACK_RAW)) {
     129        return psphotReadoutCleanup (config, view, STACK_RAW);
     130    }
     131    if (!psphotSubtractBackground (config, view, STACK_RAW)) {
     132        return psphotReadoutCleanup (config, view, STACK_RAW);
     133    }
     134    if (!psphotModelBackground (config, view, STACK_RAW)) {
     135        return psphotReadoutCleanup (config, view, STACK_RAW);
     136    }
     137    if (!psphotSubtractBackground (config, view, STACK_RAW)) {
     138        return psphotReadoutCleanup (config, view, STACK_RAW);
     139    }
     140// XXX TEST END
    122141
    123142#ifdef MAKE_CHISQ_IMAGE
    124143    // also make the chisq detection image
    125     if (!psphotStackChisqImage(config, view, STACK_DET, STACK_SRC)) {
     144    if (!psphotStackChisqImage(config, view, STACK_RAW, STACK_RAW)) {
    126145        psError (PSPHOT_ERR_UNKNOWN, false, "failure to generate chisq image");
    127         return psphotReadoutCleanup (config, view, STACK_SRC);
     146        return psphotReadoutCleanup (config, view, STACK_RAW);
     147    }
     148    if (!strcasecmp (breakPt, "CHISQ")) {
     149        return psphotReadoutCleanup (config, view, STACK_RAW);
    128150    }
    129151#endif
    130     if (!strcasecmp (breakPt, "CHISQ")) {
    131         return psphotReadoutCleanup (config, view, STACK_SRC);
    132     }
    133152
    134153    // find the detections (by peak and/or footprint) in the image.
    135154    // 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
     155    if (!psphotFindDetections (config, view, STACK_RAW, true)) { // pass 1
    138156        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         }
     157        return psphotReadoutCleanup (config, view, STACK_RAW);
    150158    }
    151159
    152160    // construct sources and measure basic stats (saved on detections->newSources)
    153     if (!psphotSourceStats (config, view, STACK_SRC, true)) { // pass 1
     161    if (!psphotSourceStats (config, view, STACK_RAW, true)) { // pass 1
    154162        psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
    155         return psphotReadoutCleanup (config, view, STACK_SRC);
     163        return psphotReadoutCleanup (config, view, STACK_RAW);
    156164    }
    157165    if (!strcasecmp (breakPt, "PEAKS")) {
    158         return psphotReadoutCleanup (config, view, STACK_SRC);
    159     }
    160     // psphotDumpTest (config, view, STACK_SRC);
     166        return psphotReadoutCleanup (config, view, STACK_RAW);
     167    }
     168    // psphotDumpTest (config, view, STACK_RAW);
    161169    psMemDump("sourcestats");
    162170    logMemStats("sourcestats");
    … …  
    164172    // classify sources based on moments, brightness
    165173    // only run this on detections from the input images, not chisq image
    166     if (!psphotRoughClass (config, view, STACK_SRC)) {
     174    if (!psphotRoughClass (config, view, STACK_RAW)) {
    167175        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)) {
     176        return psphotReadoutCleanup (config, view, STACK_RAW);
     177    }
     178
     179    // find and subtract radial profile models for saturated stars (XXX change name eventually)
     180    if (!psphotDeblendSatstars (config, view, STACK_RAW)) {
    179181        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);
     182        return psphotReadoutCleanup (config, view, STACK_RAW);
    187183    }
    188184
    189185    // if we were not supplied a PSF model, determine the IQ stats here (detections->newSources)
    190186    // only run this on detections from the input images, not chisq image
    191     if (!psphotImageQuality (config, view, STACK_SRC)) { // pass 1
     187    if (!psphotImageQuality (config, view, STACK_RAW)) { // pass 1
    192188        psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality");
    193         return psphotReadoutCleanup (config, view, STACK_SRC);
     189        return psphotReadoutCleanup (config, view, STACK_RAW);
    194190    }
    195191    if (!strcasecmp (breakPt, "MOMENTS")) {
    196         return psphotReadoutCleanup (config, view, STACK_SRC);
     192        return psphotReadoutCleanup (config, view, STACK_RAW);
    197193    }
    198194
    199195    // use bright stellar objects to measure PSF
    200     if (!psphotChoosePSF (config, view, STACK_SRC, true)) { // pass 1
     196    if (!psphotChoosePSF (config, view, STACK_RAW, true)) { // pass 1
    201197        psLogMsg ("psphot", 3, "failure to construct a psf model");
    202         return psphotReadoutCleanup (config, view, STACK_SRC);
     198        return psphotReadoutCleanup (config, view, STACK_RAW);
    203199    }
    204200    if (!strcasecmp (breakPt, "PSFMODEL")) {
    205         return psphotReadoutCleanup (config, view, STACK_SRC);
     201        return psphotReadoutCleanup (config, view, STACK_RAW);
    206202    }
    207203
    208204    // merge the newly selected sources into the existing list
    209205    // NOTE: merge OLD and NEW
    210     psphotMergeSources (config, view, STACK_SRC);
     206    psphotMergeSources (config, view, STACK_RAW);
    211207
    212208    // Construct an initial model for each object, set the radius to fitRadius, set circular
    213209    // fit mask.  NOTE: only applied to sources without guess models
    214     psphotGuessModels (config, view, STACK_SRC);
     210    psphotGuessModels (config, view, STACK_RAW);
    215211
    216212    // 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);
     213    psphotFitSourcesLinear (config, view, STACK_RAW, false, false);
     214    psphotStackVisualFilerule(config, view, STACK_RAW);
    219215
    220216    // measure the radial profiles to the sky
    221     psphotRadialProfileWings (config, view, STACK_SRC);
     217    psphotRadialProfileWings (config, view, STACK_RAW);
    222218
    223219    // re-measure the kron mags with models subtracted.  this pass starts with a circular
    … …  
    225221    // but iterates to an appropriately larger size
    226222    logMemStats("before.kron.1");
    227     psphotKronIterate(config, view, STACK_SRC, 1);
     223    psphotKronIterate(config, view, STACK_RAW, 1);
    228224    logMemStats("after.kron.1");
    229225       
    230226    // identify CRs and extended sources
    231     psphotSourceSize (config, view, STACK_SRC, true);
     227    psphotSourceSize (config, view, STACK_RAW, true);
    232228
    233229    // non-linear PSF and EXT fit to brighter sources
    234230    // replace model flux, adjust mask as needed, fit, subtract the models (full stamp)
    235     psphotBlendFit (config, view, STACK_SRC); // pass 1 (detections->allSources)
     231    psphotBlendFit (config, view, STACK_RAW); // pass 1 (detections->allSources)
    236232
    237233    // replace all sources (do NOT ignore subtraction state)
    238     psphotReplaceAllSources (config, view, STACK_SRC, false); // pass 1 (detections->allSources)
     234    psphotReplaceAllSources (config, view, STACK_RAW, false); // pass 1 (detections->allSources)
    239235
    240236    logMemStats("pass1");
    … …  
    245241    // linear fit to include all sources (subtract again)
    246242    // 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)
     243    // NOTE 2 : this function subtracts the models from the given filerule
     244    psphotFitSourcesLinear (config, view, STACK_RAW, true, false); // pass 2 (detections->allSources)
    249245
    250246    // NOTE: possibly re-measure background model here with objects subtracted / or masked
    … …  
    252248    // NOTE: this block performs the 2nd pass low-significance PSF detection stage
    253249    {
    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 
    266250        // add noise for subtracted objects
    267         psphotAddNoise (config, view, STACK_DET); // pass 1 (detections->allSources)
     251        psphotAddNoise (config, view, STACK_RAW); // pass 1 (detections->allSources)
    268252
    269253        // find fainter sources
    270254        // 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)
     255        psphotFindDetections (config, view, STACK_RAW, false); // pass 2 (detections->peaks, detections->footprints)
    272256
    273257        // remove noise for subtracted objects (ie, return to normal noise level)
    … …  
    276260        bool footprintsUseUnsubtracted = psMetadataLookupBool(NULL, recipe, "FOOTPRINT_USE_UNSUBTRACTED");
    277261        if (!footprintsUseUnsubtracted) {
    278             psphotSubNoise (config, view, STACK_DET); // pass 1 (detections->allSources)
     262            psphotSubNoise (config, view, STACK_RAW); // pass 1 (detections->allSources)
    279263        }
    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         }
    293264
    294265        // define new sources based on only the new peaks & measure moments
    295266        // NOTE: new sources are saved on detections->newSources
    296         psphotSourceStats (config, view, STACK_SRC, false); // pass 2 (detections->newSources)
     267        psphotSourceStats (config, view, STACK_RAW, false); // pass 2 (detections->newSources)
    297268
    298269        // set source type
    299270        // NOTE: apply only to detections->newSources
    300         if (!psphotRoughClass (config, view, STACK_SRC)) { // pass 2 (detections->newSources)
     271        if (!psphotRoughClass (config, view, STACK_RAW)) { // pass 2 (detections->newSources)
    301272            psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
    302             return psphotReadoutCleanup (config, view, STACK_SRC);
     273            return psphotReadoutCleanup (config, view, STACK_RAW);
    303274        }
    304275
    305276        // replace all sources so fit below applies to all at once
    306277        // NOTE: apply only to OLD sources (which have been subtracted)
    307         psphotReplaceAllSources (config, view, STACK_SRC, false); // pass 2
     278        psphotReplaceAllSources (config, view, STACK_RAW, false); // pass 2
    308279
    309280        // merge the newly selected sources into the existing list
    310281        // NOTE: merge OLD and NEW
    311282        // XXX check on free of sources...
    312         psphotMergeSources (config, view, STACK_SRC); // (detections->newSources + detections->allSources -> detections->allSources)
     283        psphotMergeSources (config, view, STACK_RAW); // (detections->newSources + detections->allSources -> detections->allSources)
    313284
    314285        // Construct an initial model for each object, set the radius to fitRadius, set circular
    315286        // fit mask.  NOTE: only applied to sources without guess models
    316         psphotGuessModels (config, view, STACK_SRC);
     287        psphotGuessModels (config, view, STACK_RAW);
    317288    }
    318289
    … …  
    325296    if (splitLinearFit) {
    326297        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.
     298        // Fit the detected sources separately from matched ones that we are about to create.
    328299        // 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)
     300        psphotFitSourcesLinear (config, view, STACK_RAW, true, true); // pass 3 (detections->allSources)
    330301    }
    331302
    … …  
    335306    // this just match the detections for the chisq image, and not bother measuring the source
    336307    // stats in that case...?
    337     objects = psphotMatchSources (config, view, STACK_SRC);
     308    objects = psphotMatchSources (config, view, STACK_RAW);
    338309    psMemDump("matchsources");
    339310
    … …  
    344315    // Construct an initial model for each object, set the radius to fitRadius, set circular
    345316    // fit mask.  NOTE: only applied to sources without guess models
    346     psphotGuessModels (config, view, STACK_SRC);
     317    psphotGuessModels (config, view, STACK_RAW);
    347318
    348319    psphotStackObjectsUnifyPosition (objects);
    349320
    350     psphotStackObjectsSelectForAnalysis (config, view, STACK_SRC, objects);
     321    psphotStackObjectsSelectForAnalysis (config, view, STACK_RAW, objects);
    351322
    352323    // final linear fit. NOTE: if splitLinearFit is true above, this pass will only fit
    353324    // the unsubtracted (matched) sources (the sources that we fit above are subtracted)
    354     psphotFitSourcesLinear (config, view, STACK_SRC, true, false); // pass 4 (detections->allSources)
     325    psphotFitSourcesLinear (config, view, STACK_RAW, true, false); // pass 4 (detections->allSources)
    355326
    356327    // measure the radial profiles to the sky (only measures new objects)
    357     psphotRadialProfileWings (config, view, STACK_SRC);
     328    psphotRadialProfileWings (config, view, STACK_RAW);
    358329
    359330    // re-measure the kron mags with models subtracted
    360331    // psphotKronMasked(config, view, STACK_SRC);
    361332    logMemStats("before.kron.2");
    362     psphotKronIterate(config, view, STACK_SRC, 2);
     333    psphotKronIterate(config, view, STACK_RAW, 2);
    363334    logMemStats("after.kron.2");
    364335
    365336    // measure source size for the remaining sources
    366337    // NOTE: applies only to NEW (unmeasured) sources
    367     psphotSourceSize (config, view, STACK_SRC, false); // pass 2 (detections->allSources)
     338    psphotSourceSize (config, view, STACK_RAW, false); // pass 2 (detections->allSources)
    368339
    369340    psMemDump("psfstats");
    … …  
    371342    // drop matched sources without any useful measurements and set kron radii for the ones
    372343    // we decide to keep
    373     psphotFilterMatchedSources (config, view, STACK_SRC, objects);
     344    psphotFilterMatchedSources (config, view, STACK_RAW, objects);
    374345
    375346    // measure kron fluxes for the matched sources only
    376     psphotKronIterate(config, view, STACK_SRC, 3);
     347    psphotKronIterate(config, view, STACK_RAW, 3);
    377348
    378349    // measure elliptical apertures, petrosians (objects sorted by S/N)
    379350    // psphotExtendedSourceAnalysisByObject (config, objects, view, STACK_SRC); // pass 1 (detections->allSources)
    380     psphotExtendedSourceAnalysis (config, view, STACK_SRC); // pass 1 (detections->allSources)
     351    psphotExtendedSourceAnalysis (config, view, STACK_RAW); // pass 1 (detections->allSources)
    381352
    382353    // measure non-linear extended source models (exponential, deVaucouleur, Sersic) (sources sorted by S/N)
    383     psphotExtendedSourceFits (config, view, STACK_SRC); // pass 1 (detections->allSources)
     354    psphotExtendedSourceFits (config, view, STACK_RAW); // pass 1 (detections->allSources)
    384355
    385356    // create source children for the OUT filerule (for radial aperture photometry and output)
    386     psArray *objectsOut = psphotSourceChildrenByObject (config, view, STACK_OUT, objects);
     357    // NOTE: The new source children have image arrays pointing to the readout associated with
     358    // STACK_OUT.  in psphotStackMatchPSFsetup, we copy the current pixel values from RAW to OUT,
     359    // but keep the pointers the same so we do not break these source image references
     360    // XXX NOTE : if we use the pre-20130914 psphotStackReadout code, we need to use 'false' for the
     361    // sourcesSubtracted argument
     362    psArray *objectsOut = psphotSourceChildrenByObject (config, view, STACK_OUT, objects, true);
    387363    if (!objectsOut) {
    388364        psFree(objects);
    389365        psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
    390         return psphotReadoutCleanup (config, view, STACK_SRC);
     366        return psphotReadoutCleanup (config, view, STACK_RAW);
    391367    }
    392368
    … …  
    396372        // this forces photometry on the undetected sources from other images
    397373
    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 
     374        // set up the FWHM vector
     375        psphotStackMatchPSFsetup (config, view, STACK_OUT, STACK_RAW);
     376        psphotDumpImages (config, view, STACK_RAW, "raw.t0");
     377        psphotDumpImages (config, view, STACK_OUT, "out.t0");
    417378
    418379        int nRadialEntries = psphotStackMatchPSFsEntries(config, view, STACK_OUT);
    419         for (int entry = 1; entry < nRadialEntries; entry++) {
     380
     381        for (int entry = 0; entry < nRadialEntries; entry++) {
    420382            // NOTE: entry 0 is the unmatched image set
    421383
    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);
     384            char line[256];
    430385
    431386            // measure circular, radial apertures (objects sorted by S/N)
    432387            psphotRadialApertures (config, view, STACK_OUT, entry);
     388            snprintf (line, 256, "%s.%d", "out.t1", entry);
     389            psphotDumpImages (config, view, STACK_OUT, line);
    433390
    434391            // replace the flux in the image so it is returned to its original state
    435392            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");
     393            snprintf (line, 256, "%s.%d", "out.t2", entry);
     394            psphotDumpImages (config, view, STACK_OUT, line);
     395
     396            if (entry < nRadialEntries - 1) {
     397                // smooth to the next FWHM
     398                // this function does nothing if the targetFWHM is smaller than the currentFWHM
     399                psphotStackMatchPSFsNext(config, view, STACK_OUT, entry);
     400                snprintf (line, 256, "%s.%d", "out.t3", entry);
     401                psphotDumpImages (config, view, STACK_OUT, line);
     402                psMemDump("matched");
     403
     404                // re-measure the PSF for the smoothed image (using entries in 'allSources')
     405                psphotChoosePSF (config, view, STACK_OUT, false);
     406
     407                // this is necessary to update the models based on the new PSF
     408                psphotResetModels (config, view, STACK_OUT);
     409
     410                // this is necessary to get the right normalization for the new models
     411                // and to subtract the sources
     412                psphotFitSourcesLinear (config, view, STACK_OUT, false, false);
     413                snprintf (line, 256, "%s.%d", "out.t4", entry);
     414                psphotDumpImages (config, view, STACK_OUT, line);
     415            }
    440416        }
    441417    }
    442     psphotStackSetInputsToSkip(config, view, STACK_CNV, false);
    443     psphotStackSetInputsToSkip(config, view, STACK_OUT, false);
    444418
    445419    // measure aperture photometry corrections
    446     if (!psphotApResid (config, view, STACK_SRC)) {
     420    if (!psphotApResid (config, view, STACK_RAW)) {
    447421        psFree (objects);
    448422        psFree (objectsOut);
    449423        psLogMsg ("psphot", 3, "failed on psphotApResid");
    450         return psphotReadoutCleanup (config, view, STACK_SRC);
     424        return psphotReadoutCleanup (config, view, STACK_RAW);
    451425    }
    452426
    453427    // 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)) {
     428    psphotMagnitudes(config, view, STACK_RAW);
     429
     430    if (!psphotEfficiency(config, view, STACK_RAW)) {
    465431        psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources");
    466432        psErrorClear();
    467433    }
    468     psphotCopyEfficiency (config, view, STACK_OUT, STACK_DET);
     434    psphotCopyEfficiency (config, view, STACK_OUT, STACK_RAW);
    469435
    470436    logMemStats("final");
    471437#if (1)
    472     psphotSourceMemory(config, view, STACK_SRC);
     438    psphotSourceMemory(config, view, STACK_RAW);
    473439    psphotSourceMemory(config, view, STACK_OUT);
    474440#endif
    475441
    476     // replace failed sources?
    477     // psphotReplaceUnfitSources (sources);
    478 
    479442    // replace background in residual image
    480     psphotSkyReplace (config, view, STACK_DET);
     443    psphotSkyReplace (config, view, STACK_RAW);
    481444
    482445    // drop the references to the image pixels held by each source
     446    psphotSourceFreePixels (config, view, STACK_RAW);
    483447    psphotSourceFreePixels (config, view, STACK_OUT);
    484     psphotSourceFreePixels (config, view, STACK_SRC);
    485448
    486449#ifdef MAKE_CHISQ_IMAGE
    487450    // 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     }
     451    psphotStackRemoveChisqFromInputs(config, STACK_RAW);
    492452#endif
    493453
    … …  
    496456
    497457    // create the exported-metadata and free local data
    498     return psphotReadoutCleanup (config, view, STACK_SRC);
     458    return psphotReadoutCleanup (config, view, STACK_RAW);
    499459}
    500460
    … …  
    552512}
    553513
    554 
    555 
    556514/* here is the process:
    557515
    558  * we have three(*) images:
    559  * RAW : unconvolved image stack
    560  * CNV : input convolved image
     516 * we have two image sets:
     517 * RAW : unconvolved image stacks
    561518
    562519 * OUT : psf-matched output image (there may be more than one of
    … …  
    649606
    650607   */
     608
     609
     610// generate a vector fwhmValues where the first has the fwhm of the raw image and the
     611// successive entries have the target values
     612
     613bool psphotStackMatchPSFsetup (pmConfig *config, const pmFPAview *view, const char *filerule, const char *filerulePSF) {
     614
     615    bool status;
     616
     617    int num = psphotFileruleCount(config, filerule);
     618
     619    // skip the chisq image (optionally?)
     620    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
     621    if (!status) chisqNum = -1;
     622
     623    // loop over the available readouts
     624    for (int i = 0; i < num; i++) {
     625        if (i == chisqNum) continue; // skip chisq image
     626
     627        if (!psphotStackMatchPSFsetupReadout (config, view, filerule, filerulePSF, i)) {
     628            psError (PSPHOT_ERR_CONFIG, false, "failed to define target PSF sizes");
     629            return false;
     630        }
     631    }
     632
     633    return true;
     634}
     635
     636float psphotPSFseeing (pmPSF *psf, pmReadout *readout, int index);
     637
     638// copy the pixels from RAW to OUT (
     639
     640bool psphotStackMatchPSFsetupReadout (pmConfig *config, const pmFPAview *view, const char *fileruleOut, const char *fileruleRaw, int index) {
     641
     642    bool status;
     643
     644    // select the appropriate recipe information
     645    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     646
     647    // find the currently selected readout
     648    pmFPAfile *fileOut = pmFPAfileSelectSingle(config->files, fileruleOut, index); // File of interest
     649    psAssert (fileOut, "missing file?");
     650
     651    pmReadout *readoutOut = pmFPAviewThisReadout(view, fileOut->fpa);
     652    psAssert (readoutOut, "missing readout?");
     653
     654    // find the currently selected readout
     655    pmFPAfile *fileRaw = pmFPAfileSelectSingle(config->files, fileruleRaw, index); // File of interest
     656    psAssert (fileRaw, "missing file?");
     657
     658    pmReadout *readoutRaw = pmFPAviewThisReadout(view, fileRaw->fpa);
     659    psAssert (readoutRaw, "missing readout?");
     660
     661    readoutOut->image = psImageCopy(readoutOut->image, readoutRaw->image, PS_TYPE_F32);
     662    if (readoutRaw->variance) {
     663        readoutOut->variance = psImageCopy(readoutOut->variance, readoutRaw->variance, PS_TYPE_F32);
     664    }
     665    if (readoutRaw->mask) {
     666        readoutOut->mask = psImageCopy(readoutOut->mask, readoutRaw->mask, PS_TYPE_IMAGE_MASK);
     667    }
     668
     669    // pmChip *chipRaw = pmFPAviewThisChip(view, fileRaw->fpa); // The chip holds the PSF
     670    // psAssert (chipRaw, "missing chip");
     671
     672    pmPSF *psf = psMetadataLookupPtr(&status, readoutRaw->analysis, "PSPHOT.PSF"); // PSF
     673    if (!psf) {
     674        // we should have a PSF by this point in psphot
     675        psError(PSPHOT_ERR_PROG, true, "Unable to find PSF.");
     676        return false;
     677    }
     678
     679    float fwhmRaw = psphotPSFseeing (psf, readoutRaw, index);
     680
     681    psVector *fwhmValues = psVectorAllocEmpty(10, PS_TYPE_F32);
     682    psVectorAppend(fwhmValues, fwhmRaw);
     683
     684    // is a single target FWHM specified, or a set of values?  set up the vector options->targetSeeing and the local 1st value
     685    float targetSeeing = psMetadataLookupF32 (&status, recipe, "PSPHOT.STACK.TARGET.PSF.FWHM");
     686    if (!status) {
     687        psVector *targetSeeing = psMetadataLookupVector(&status, recipe, "PSPHOT.STACK.TARGET.PSF.FWHM"); // Magnitude offsets
     688        psAssert (status, "missing psphot recipe value PSPHOT.STACK.TARGET.PSF.FWHM");
     689        for (int i = 0; i < targetSeeing->n; i++) {
     690            psVectorAppend(fwhmValues, targetSeeing->data.F32[i]);
     691        }           
     692    } else {
     693        psVectorAppend(fwhmValues, targetSeeing);
     694    }
     695
     696    psMetadataAddVector(readoutOut->analysis, PS_LIST_TAIL, "STACK.PSF.FWHM.VALUES", PS_META_REPLACE, "PSF sizes", fwhmValues);
     697    psFree (fwhmValues);
     698
     699    return true;
     700}
     701
     702float psphotPSFseeing (pmPSF *psf, pmReadout *readout, int index) {
     703
     704    psImage *image = readout->image;
     705
     706    int Nx = image->numCols;
     707    int Ny = image->numRows;
     708
     709    float sumFWHM = 0.0;                  // FWHM for image
     710    int numFWHM = 0;                      // Number of FWHM measurements
     711    for (float x = 0; x < Nx; x += 0.25*Nx) {
     712        for (float y = 0; y < Ny; y += 0.25*Ny) {
     713            float fwhm = pmPSFtoFWHM(psf, x, y);
     714            if (isfinite(fwhm)) {
     715                sumFWHM += fwhm;
     716                numFWHM++;
     717            }
     718        }
     719    }
     720    if (numFWHM == 0) {
     721        psLogMsg("ppStack", PS_LOG_INFO, "Unable to measure PSF FWHM for image %d --- rejected.", index);
     722        return NAN;
     723    }
     724
     725    float fwhm = sumFWHM / (float) numFWHM;
     726    psLogMsg ("psphotStack", PS_LOG_INFO, "Input Seeing for %d: %f\n", index, fwhm);
     727
     728    return fwhm;
     729}
Note: See TracChangeset for help on using the changeset viewer.