IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jul 17, 2014, 12:30:45 PM (12 years ago)
Author:
eugene
Message:

merge changes (from past YEAR) into trunk

Location:
branches/eam_branches/ipp-ops-20130712/psphot
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-ops-20130712/psphot

    • Property svn:mergeinfo deleted

  • branches/eam_branches/ipp-ops-20130712/psphot/src

    • Property svn:ignore

      •  

        old new  
        2424psphotModelTest
        2525psphotMinimal
         26psphotFullForce
         27psmakecff
         28psphotFullForceSummary
    • Property svn:mergeinfo deleted

  • branches/eam_branches/ipp-ops-20130712/psphot/src/psphotStackReadout.c

    r34721 r37066  
    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);
     333
     334    // decide which source(s) are to be fitted with the extended source analysis code.
     335    psphotChooseAnalysisOptionsByObject (config, view, STACK_RAW, objects);
    377336
    378337    // measure elliptical apertures, petrosians (objects sorted by S/N)
    379338    // psphotExtendedSourceAnalysisByObject (config, objects, view, STACK_SRC); // pass 1 (detections->allSources)
    380     psphotExtendedSourceAnalysis (config, view, STACK_SRC); // pass 1 (detections->allSources)
     339    psphotExtendedSourceAnalysis (config, view, STACK_RAW); // pass 1 (detections->allSources)
    381340
    382341    // measure non-linear extended source models (exponential, deVaucouleur, Sersic) (sources sorted by S/N)
    383     psphotExtendedSourceFits (config, view, STACK_SRC); // pass 1 (detections->allSources)
     342    psphotExtendedSourceFits (config, view, STACK_RAW); // pass 1 (detections->allSources)
    384343
    385344    // create source children for the OUT filerule (for radial aperture photometry and output)
    386     psArray *objectsOut = psphotSourceChildrenByObject (config, view, STACK_OUT, objects);
     345    // NOTE: The new source children have image arrays pointing to the readout associated with
     346    // STACK_OUT.  in psphotStackMatchPSFsetup, we copy the current pixel values from RAW to OUT,
     347    // but keep the pointers the same so we do not break these source image references
     348    // XXX NOTE : if we use the pre-20130914 psphotStackReadout code, we need to use 'false' for the
     349    // sourcesSubtracted argument
     350    psArray *objectsOut = psphotSourceChildrenByObject (config, view, STACK_OUT, objects, true);
    387351    if (!objectsOut) {
    388352        psFree(objects);
    389353        psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
    390         return psphotReadoutCleanup (config, view, STACK_SRC);
     354        return psphotReadoutCleanup (config, view, STACK_RAW);
    391355    }
    392356
    … …  
    396360        // this forces photometry on the undetected sources from other images
    397361
    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 
     362        // set up the FWHM vector
     363        psphotStackMatchPSFsetup (config, view, STACK_OUT, STACK_RAW);
     364        psphotDumpImages (config, view, STACK_RAW, "raw.t0");
     365        psphotDumpImages (config, view, STACK_OUT, "out.t0");
    417366
    418367        int nRadialEntries = psphotStackMatchPSFsEntries(config, view, STACK_OUT);
    419         for (int entry = 1; entry < nRadialEntries; entry++) {
     368
     369        for (int entry = 0; entry < nRadialEntries; entry++) {
    420370            // NOTE: entry 0 is the unmatched image set
    421371
    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);
     372            char line[256];
    430373
    431374            // measure circular, radial apertures (objects sorted by S/N)
    432375            psphotRadialApertures (config, view, STACK_OUT, entry);
     376            snprintf (line, 256, "%s.%d", "out.t1", entry);
     377            psphotDumpImages (config, view, STACK_OUT, line);
    433378
    434379            // replace the flux in the image so it is returned to its original state
    435380            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");
     381            snprintf (line, 256, "%s.%d", "out.t2", entry);
     382            psphotDumpImages (config, view, STACK_OUT, line);
     383
     384            if (entry < nRadialEntries - 1) {
     385                // smooth to the next FWHM
     386                // this function does nothing if the targetFWHM is smaller than the currentFWHM
     387                psphotStackMatchPSFsNext(config, view, STACK_OUT, entry);
     388                snprintf (line, 256, "%s.%d", "out.t3", entry);
     389                psphotDumpImages (config, view, STACK_OUT, line);
     390                psMemDump("matched");
     391
     392                // re-measure the PSF for the smoothed image (using entries in 'allSources')
     393                if (!psphotChoosePSF (config, view, STACK_OUT, false)) {
     394                    psLogMsg ("psphot", 3, "failure to construct a psf model in radial aperture loop for entry :%d", entry);
     395                    return psphotReadoutCleanup (config, view, STACK_RAW);
     396                }
     397
     398                // this is necessary to update the models based on the new PSF
     399                psphotResetModels (config, view, STACK_OUT);
     400
     401                // this is necessary to get the right normalization for the new models
     402                // and to subtract the sources
     403                psphotFitSourcesLinear (config, view, STACK_OUT, false, false);
     404                snprintf (line, 256, "%s.%d", "out.t4", entry);
     405                psphotDumpImages (config, view, STACK_OUT, line);
     406            }
    440407        }
    441408    }
    442     psphotStackSetInputsToSkip(config, view, STACK_CNV, false);
    443     psphotStackSetInputsToSkip(config, view, STACK_OUT, false);
    444409
    445410    // measure aperture photometry corrections
    446     if (!psphotApResid (config, view, STACK_SRC)) {
     411    if (!psphotApResid (config, view, STACK_RAW)) {
    447412        psFree (objects);
    448413        psFree (objectsOut);
    449414        psLogMsg ("psphot", 3, "failed on psphotApResid");
    450         return psphotReadoutCleanup (config, view, STACK_SRC);
     415        return psphotReadoutCleanup (config, view, STACK_RAW);
    451416    }
    452417
    453418    // 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)) {
     419    psphotMagnitudes(config, view, STACK_RAW);
     420
     421    // calculate lensing parameters
     422    if (!psphotLensing(config, view, STACK_RAW)) {
     423        psErrorStackPrint(stderr, "Unable to do lensing parameters.");
     424        psErrorClear();
     425    }
     426
     427    if (!psphotEfficiency(config, view, STACK_RAW)) {
    465428        psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources");
    466429        psErrorClear();
    467430    }
    468     psphotCopyEfficiency (config, view, STACK_OUT, STACK_DET);
     431    psphotCopyEfficiency (config, view, STACK_OUT, STACK_RAW);
    469432
    470433    logMemStats("final");
    471434#if (1)
    472     psphotSourceMemory(config, view, STACK_SRC);
     435    psphotSourceMemory(config, view, STACK_RAW);
    473436    psphotSourceMemory(config, view, STACK_OUT);
    474437#endif
    475438
    476     // replace failed sources?
    477     // psphotReplaceUnfitSources (sources);
    478 
    479439    // replace background in residual image
    480     psphotSkyReplace (config, view, STACK_DET);
     440    psphotSkyReplace (config, view, STACK_RAW);
    481441
    482442    // drop the references to the image pixels held by each source
     443    psphotSourceFreePixels (config, view, STACK_RAW);
    483444    psphotSourceFreePixels (config, view, STACK_OUT);
    484     psphotSourceFreePixels (config, view, STACK_SRC);
    485445
    486446#ifdef MAKE_CHISQ_IMAGE
    487447    // 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     }
     448    psphotStackRemoveChisqFromInputs(config, STACK_RAW);
    492449#endif
    493450
    … …  
    496453
    497454    // create the exported-metadata and free local data
    498     return psphotReadoutCleanup (config, view, STACK_SRC);
     455    return psphotReadoutCleanup (config, view, STACK_RAW);
    499456}
    500457
    … …  
    552509}
    553510
    554 
    555 
    556511/* here is the process:
    557512
    558  * we have three(*) images:
    559  * RAW : unconvolved image stack
    560  * CNV : input convolved image
     513 * we have two image sets:
     514 * RAW : unconvolved image stacks
    561515
    562516 * OUT : psf-matched output image (there may be more than one of
    … …  
    649603
    650604   */
     605
     606
     607// generate a vector fwhmValues where the first has the fwhm of the raw image and the
     608// successive entries have the target values
     609
     610bool psphotStackMatchPSFsetup (pmConfig *config, const pmFPAview *view, const char *filerule, const char *filerulePSF) {
     611
     612    bool status;
     613
     614    int num = psphotFileruleCount(config, filerule);
     615
     616    // skip the chisq image (optionally?)
     617    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
     618    if (!status) chisqNum = -1;
     619
     620    // loop over the available readouts
     621    for (int i = 0; i < num; i++) {
     622        if (i == chisqNum) continue; // skip chisq image
     623
     624        if (!psphotStackMatchPSFsetupReadout (config, view, filerule, filerulePSF, i)) {
     625            psError (PSPHOT_ERR_CONFIG, false, "failed to define target PSF sizes");
     626            return false;
     627        }
     628    }
     629
     630    return true;
     631}
     632
     633float psphotPSFseeing (pmPSF *psf, pmReadout *readout, int index);
     634
     635// copy the pixels from RAW to OUT (
     636
     637bool psphotStackMatchPSFsetupReadout (pmConfig *config, const pmFPAview *view, const char *fileruleOut, const char *fileruleRaw, int index) {
     638
     639    bool status;
     640
     641    // select the appropriate recipe information
     642    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     643
     644    // find the currently selected readout
     645    pmFPAfile *fileOut = pmFPAfileSelectSingle(config->files, fileruleOut, index); // File of interest
     646    psAssert (fileOut, "missing file?");
     647
     648    pmReadout *readoutOut = pmFPAviewThisReadout(view, fileOut->fpa);
     649    psAssert (readoutOut, "missing readout?");
     650
     651    // find the currently selected readout
     652    pmFPAfile *fileRaw = pmFPAfileSelectSingle(config->files, fileruleRaw, index); // File of interest
     653    psAssert (fileRaw, "missing file?");
     654
     655    pmReadout *readoutRaw = pmFPAviewThisReadout(view, fileRaw->fpa);
     656    psAssert (readoutRaw, "missing readout?");
     657
     658    readoutOut->image = psImageCopy(readoutOut->image, readoutRaw->image, PS_TYPE_F32);
     659    if (readoutRaw->variance) {
     660        readoutOut->variance = psImageCopy(readoutOut->variance, readoutRaw->variance, PS_TYPE_F32);
     661    }
     662    if (readoutRaw->mask) {
     663        readoutOut->mask = psImageCopy(readoutOut->mask, readoutRaw->mask, PS_TYPE_IMAGE_MASK);
     664    }
     665
     666    // pmChip *chipRaw = pmFPAviewThisChip(view, fileRaw->fpa); // The chip holds the PSF
     667    // psAssert (chipRaw, "missing chip");
     668
     669    pmPSF *psf = psMetadataLookupPtr(&status, readoutRaw->analysis, "PSPHOT.PSF"); // PSF
     670    if (!psf) {
     671        // we should have a PSF by this point in psphot
     672        psError(PSPHOT_ERR_PROG, true, "Unable to find PSF.");
     673        return false;
     674    }
     675
     676    float fwhmRaw = psphotPSFseeing (psf, readoutRaw, index);
     677
     678    psVector *fwhmValues = psVectorAllocEmpty(10, PS_TYPE_F32);
     679    psVectorAppend(fwhmValues, fwhmRaw);
     680
     681    // is a single target FWHM specified, or a set of values?  set up the vector options->targetSeeing and the local 1st value
     682    float targetSeeing = psMetadataLookupF32 (&status, recipe, "PSPHOT.STACK.TARGET.PSF.FWHM");
     683    if (!status) {
     684        psVector *targetSeeing = psMetadataLookupVector(&status, recipe, "PSPHOT.STACK.TARGET.PSF.FWHM"); // Magnitude offsets
     685        psAssert (status, "missing psphot recipe value PSPHOT.STACK.TARGET.PSF.FWHM");
     686        for (int i = 0; i < targetSeeing->n; i++) {
     687            psVectorAppend(fwhmValues, targetSeeing->data.F32[i]);
     688        }           
     689    } else {
     690        psVectorAppend(fwhmValues, targetSeeing);
     691    }
     692
     693    psMetadataAddVector(readoutOut->analysis, PS_LIST_TAIL, "STACK.PSF.FWHM.VALUES", PS_META_REPLACE, "PSF sizes", fwhmValues);
     694    psFree (fwhmValues);
     695
     696    return true;
     697}
     698
     699float psphotPSFseeing (pmPSF *psf, pmReadout *readout, int index) {
     700
     701    psImage *image = readout->image;
     702
     703    int Nx = image->numCols;
     704    int Ny = image->numRows;
     705
     706    float sumFWHM = 0.0;                  // FWHM for image
     707    int numFWHM = 0;                      // Number of FWHM measurements
     708    for (float x = 0; x < Nx; x += 0.25*Nx) {
     709        for (float y = 0; y < Ny; y += 0.25*Ny) {
     710            float fwhm = pmPSFtoFWHM(psf, x, y);
     711            if (isfinite(fwhm)) {
     712                sumFWHM += fwhm;
     713                numFWHM++;
     714            }
     715        }
     716    }
     717    if (numFWHM == 0) {
     718        psLogMsg("ppStack", PS_LOG_INFO, "Unable to measure PSF FWHM for image %d --- rejected.", index);
     719        return NAN;
     720    }
     721
     722    float fwhm = sumFWHM / (float) numFWHM;
     723    psLogMsg ("psphotStack", PS_LOG_INFO, "Input Seeing for %d: %f\n", index, fwhm);
     724
     725    return fwhm;
     726}
Note: See TracChangeset for help on using the changeset viewer.