IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Sep 5, 2011, 9:02:13 AM (15 years ago)
Author:
eugene
Message:

major re-work of psphotStackReadout to be more consistent with psphotReadout; psphotRadialApertures now assumes the 0 element is the unconvolved version; for the convolved versions, the vector of target psfs is saved as STACK.PSF.FWHM.VALUES

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r31154 r32322  
    4444bool psphotStackReadout (pmConfig *config, const pmFPAview *view) {
    4545
     46    // measure the total elapsed time in psphotReadout.  dtime is the elapsed time used jointly
     47    // by the multiple threads, not the total time used by all threads.
    4648    psTimerStart ("psphotReadout");
    4749
    … …  
    5658    // optional break-point for processing
    5759    char *breakPt = psMetadataLookupStr (NULL, recipe, "BREAK_POINT");
    58     PS_ASSERT_PTR_NON_NULL (breakPt, false);
    59 
    60     // XXX or do I set OUT to be a pmFPAfile pointing at the input of interest?
     60    psAssert (breakPt, "configuration error: set BREAK_POINT");
     61
     62    // we have 3 relevant files: RAW (unconvolved), CNV (convolved stack), OUT (psf-matched stack)
     63    // select which image (RAW or CNV) is used for analysis (RAW always used for detection)
    6164    bool useRaw = psMetadataLookupBool (NULL, recipe, "PSPHOT.STACK.USE.RAW");
    6265    char *STACK_SRC = useRaw ? STACK_RAW : STACK_CNV;
    63     char *STACK_DET = STACK_RAW; // XXX optionally allow this to be CNV?
    64 
    65     // we have 3 relevant files: RAW, CNV, OUT
     66    char *STACK_DET = STACK_RAW;
    6667
    6768    // set the photcode for each image
    … …  
    7172    }
    7273
    73     // Generate the mask and weight images
    74     // XXX this should be done before we perform the convolutions
     74    // Generate the mask and weight images (if not supplied) and set mask bits
    7575    if (!psphotSetMaskAndVariance (config, view, STACK_DET)) {
    7676        return psphotReadoutCleanup (config, view, STACK_SRC);
    … …  
    8080    }
    8181
     82    // XXX I think this is not defined correctly for an array of images.
     83    // XXX I probably need to subtract the model (same model?) for both RAW and OUT.
     84    // XXX But, probably not a problem in practice since the stacks are constructed with 0.0 mean level.
     85
    8286    // generate a background model (median, smoothed image)
    83     // XXX I think this is not defined correctly for an array of images.
    84     // XXX probably need to subtract the model (same model?) for both RAW and OUT
    8587    if (!psphotModelBackground (config, view, STACK_DET)) {
    8688        return psphotReadoutCleanup (config, view, STACK_SRC);
    … …  
    9395    }
    9496
     97    // also make the chisq detection image
    9598    if (!psphotStackChisqImage(config, view, STACK_DET, STACK_SRC)) {
    9699        psError (PSPHOT_ERR_UNKNOWN, false, "failure to generate chisq image");
    … …  
    109112    }
    110113
    111     // copy the detections from DET to SRC
     114    // if DET and SRC are different images, copy the detections from DET to SRC
    112115    if (strcmp(STACK_SRC, STACK_DET)) {
    113116        if (!psphotCopySources (config, view, STACK_SRC, STACK_DET)) {
    … …  
    122125        return psphotReadoutCleanup (config, view, STACK_SRC);
    123126    }
    124 
    125     if (!strcasecmp (breakPt, "TEST1")) {
    126         return psphotReadoutCleanup (config, view, STACK_SRC);
    127     }
    128 
     127    if (!strcasecmp (breakPt, "PEAKS")) {
     128        return psphotReadoutCleanup (config, view, STACK_SRC);
     129    }
    129130    psMemDump("sourcestats");
    130131
    131     // generate the objects (object unify the sources from the different images)
     132    // classify sources based on moments, brightness
     133    // only run this on detections from the input images, not chisq image
     134    if (!psphotRoughClass (config, view, STACK_SRC)) {
     135        psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough classifications");
     136        return psphotReadoutCleanup (config, view, STACK_SRC);
     137    }
     138
     139    // if we were not supplied a PSF model, determine the IQ stats here (detections->newSources)
     140    // only run this on detections from the input images, not chisq image
     141    if (!psphotImageQuality (config, view, STACK_SRC)) { // pass 1
     142        psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality");
     143        return psphotReadoutCleanup (config, view, STACK_SRC);
     144    }
     145    if (!strcasecmp (breakPt, "MOMENTS")) {
     146        return psphotReadoutCleanup (config, view, STACK_SRC);
     147    }
     148
     149    // use bright stellar objects to measure PSF
     150    if (!psphotChoosePSF (config, view, STACK_SRC, true)) { // pass 1
     151        psLogMsg ("psphot", 3, "failure to construct a psf model");
     152        return psphotReadoutCleanup (config, view, STACK_SRC);
     153    }
     154    if (!strcasecmp (breakPt, "PSFMODEL")) {
     155        return psphotReadoutCleanup (config, view, STACK_SRC);
     156    }
     157
     158    // construct an initial model for each object, set the radius to fitRadius, set circular fit mask
     159    psphotGuessModels (config, view, STACK_SRC);
     160
     161    // merge the newly selected sources into the existing list
     162    // NOTE: merge OLD and NEW
     163    psphotMergeSources (config, view, STACK_SRC);
     164
     165    // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
     166    // XXX why do this as a stack operation?
     167    // psphotFitSourcesLinearStack (config, objects, false);
     168    psphotFitSourcesLinear (config, view, STACK_SRC, false);
     169    psphotStackVisualFilerule(config, view, STACK_SRC);
     170
     171    // re-measure the kron mags with models subtracted.  this pass starts with a circular
     172    // window of size PSF_MOMENTS_RADIUS (same window used to measure the psf-scale moments)
     173    // but iterates to an appropriately larger size
     174    psphotKronIterate(config, view, STACK_SRC);
     175
     176    // identify CRs and extended sources
     177    psphotSourceSize (config, view, STACK_SRC, true);
     178
     179    // non-linear PSF and EXT fit to brighter sources
     180    // replace model flux, adjust mask as needed, fit, subtract the models (full stamp)
     181    psphotBlendFit (config, view, STACK_SRC); // pass 1 (detections->allSources)
     182
     183    // replace all sources
     184    psphotReplaceAllSources (config, view, STACK_SRC); // pass 1 (detections->allSources)
     185
     186    // linear fit to include all sources (subtract again)
     187    // NOTE : apply to ALL sources (extended + psf)
     188    psphotFitSourcesLinear (config, view, STACK_SRC, true); // pass 2 (detections->allSources)
     189
     190    // if we only do one pass, skip to extended source analysis
     191    if (!strcasecmp (breakPt, "PASS1")) goto pass1finish;
     192
     193    // NOTE: possibly re-measure background model here with objects subtracted / or masked
     194
     195    // NOTE: this block performs the 2nd pass low-significance PSF detection stage
     196    {
     197        // add noise for subtracted objects
     198        psphotAddNoise (config, view, STACK_DET); // pass 1 (detections->allSources)
     199
     200        // find fainter sources
     201        // NOTE: finds new peaks and new footprints, OLD and FULL set are saved on detections
     202        psphotFindDetections (config, view, STACK_DET, false); // pass 2 (detections->peaks, detections->footprints)
     203
     204        // remove noise for subtracted objects (ie, return to normal noise level)
     205        // NOTE: this needs to operate only on the OLD sources
     206        psphotSubNoise (config, view, STACK_DET); // pass 1 (detections->allSources)
     207
     208        // if DET and SRC are different images, copy the detections from DET to SRC
     209        if (strcmp(STACK_SRC, STACK_DET)) {
     210            // XXX how does this handle 1st vs 2nd pass sources?
     211            if (!psphotCopySources (config, view, STACK_SRC, STACK_DET)) {
     212                psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
     213                return psphotReadoutCleanup (config, view, STACK_SRC);
     214            }
     215        }
     216
     217        // define new sources based on only the new peaks & measure moments
     218        // NOTE: new sources are saved on detections->newSources
     219        psphotSourceStats (config, view, STACK_SRC, false); // pass 2 (detections->newSources)
     220
     221        // set source type
     222        // NOTE: apply only to detections->newSources
     223        if (!psphotRoughClass (config, view, STACK_SRC)) { // pass 2 (detections->newSources)
     224            psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
     225            return psphotReadoutCleanup (config, view, STACK_SRC);
     226        }
     227
     228        // create full input models, set the radius to fitRadius, set circular fit mask
     229        // NOTE: apply only to detections->newSources
     230        psphotGuessModels (config, view, STACK_SRC); // pass 2 (detections->newSources)
     231
     232        // replace all sources so fit below applies to all at once
     233        // NOTE: apply only to OLD sources (which have been subtracted)
     234        psphotReplaceAllSources (config, view, STACK_SRC); // pass 2
     235
     236        // merge the newly selected sources into the existing list
     237        // NOTE: merge OLD and NEW
     238        // XXX check on free of sources...
     239        psphotMergeSources (config, view, STACK_SRC); // (detections->newSources + detections->allSources -> detections->allSources)
     240
     241        // NOTE: apply to ALL sources
     242        psphotFitSourcesLinear (config, view, STACK_SRC, true); // pass 3 (detections->allSources)
     243    }
     244
     245pass1finish:
     246
     247    // re-measure the kron mags with models subtracted
     248    // psphotKronMasked(config, view, STACK_SRC);
     249    psphotKronIterate(config, view, STACK_SRC);
     250
     251    // measure source size for the remaining sources
     252    // NOTE: applies only to NEW (unmeasured) sources
     253    psphotSourceSize (config, view, STACK_SRC, false); // pass 2 (detections->allSources)
     254
     255    psMemDump("psfstats");
     256
     257    // generate the objects (objects unify the sources from the different images)
    132258    // XXX this could just match the detections for the chisq image, and not bother measuring the
    133259    // source stats in that case...
    134260    psArray *objects = psphotMatchSources (config, view, STACK_SRC);
    135 
    136261    psMemDump("matchsources");
    137262
    138     if (!strcasecmp (breakPt, "TEST2")) {
    139         psFree(objects);
    140         return psphotReadoutCleanup (config, view, STACK_SRC);
    141     }
    142 
    143     // construct sources for the newly-generated sources (from other images)
    144     if (!psphotSourceStats (config, view, STACK_SRC, false)) { // pass 1
    145         psFree(objects);
    146         psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
    147         return psphotReadoutCleanup (config, view, STACK_SRC);
    148     }
    149 
    150     psMemDump("sourcestats");
    151 
    152     // find blended neighbors of very saturated stars (detections->newSources)
    153     // if (!psphotDeblendSatstars (config, view)) {
    154     //     psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis");
    155     //     return psphotReadoutCleanup (config, view, STACK_SRC);
    156     // }
    157 
    158     // mark blended peaks PS_SOURCE_BLEND (detections->newSources)
    159     // if (!psphotBasicDeblend (config, view)) {
    160     //     psError (PSPHOT_ERR_UNKNOWN, false, "failed on deblend analysis");
    161     //     return psphotReadoutCleanup (config, view, STACK_SRC);
    162     // }
    163 
    164     // classify sources based on moments, brightness
    165     // only run this on detections from the input images, not chisq image
    166     if (!psphotRoughClass (config, view, STACK_SRC)) {
    167         psFree(objects);
    168         psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough classifications");
    169         return psphotReadoutCleanup (config, view, STACK_SRC);
    170     }
    171     // if we were not supplied a PSF model, determine the IQ stats here (detections->newSources)
    172     // only run this on detections from the input images, not chisq image
    173     if (!psphotImageQuality (config, view, STACK_SRC)) { // pass 1
    174         psFree(objects);
    175         psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality");
    176         return psphotReadoutCleanup (config, view, STACK_SRC);
    177     }
    178     if (!strcasecmp (breakPt, "MOMENTS")) {
    179         psFree(objects);
    180         return psphotReadoutCleanup (config, view, STACK_SRC);
    181     }
    182 
    183     // use bright stellar objects to measure PSF if we were supplied a PSF for any input file,
    184     // this step is skipped
    185     if (!psphotChoosePSF (config, view, STACK_SRC, true)) { // pass 1
    186         psFree(objects);
    187         psLogMsg ("psphot", 3, "failure to construct a psf model");
    188         return psphotReadoutCleanup (config, view, STACK_SRC);
    189     }
    190     if (!strcasecmp (breakPt, "PSFMODEL")) {
    191         psFree(objects);
    192         return psphotReadoutCleanup (config, view, STACK_SRC);
    193     }
    194 
    195     // construct an initial model for each object, set the radius to fitRadius, set circular fit mask
    196     psphotGuessModels (config, view, STACK_SRC);
    197 
    198     // merge the newly selected sources into the existing list
    199     // NOTE: merge OLD and NEW
    200     psphotMergeSources (config, view, STACK_SRC);
    201 
    202     // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
    203     psphotFitSourcesLinearStack (config, objects, FALSE);
    204     psphotStackVisualFilerule(config, view, STACK_SRC);
    205 
    206     // identify CRs and extended sources
    207     psphotSourceSize (config, view, STACK_SRC, TRUE);
    208 
    209     // XXX do we want to do a preliminary (unconvolved) model fit here, and then
    210     // do a second detection pass? (like standard psphot)
    211 
    212     // measure aperture photometry corrections
    213     if (!psphotApResid (config, view, STACK_SRC)) {
    214         psFree (objects);
    215         psLogMsg ("psphot", 3, "failed on psphotApResid");
    216         return psphotReadoutCleanup (config, view, STACK_SRC);
    217     }
    218 
    219     psMemDump("psfstats");
    220 
    221263    psphotStackObjectsUnifyPosition (objects);
    222264
     265    psphotStackObjectsSelectForAnalysis (config, view, STACK_SRC, objects);
     266
    223267    // measure elliptical apertures, petrosians (objects sorted by S/N)
    224     psphotExtendedSourceAnalysisByObject (config, objects, view, STACK_SRC); // pass 1 (detections->allSources)
     268    // psphotExtendedSourceAnalysisByObject (config, objects, view, STACK_SRC); // pass 1 (detections->allSources)
     269    psphotExtendedSourceAnalysis (config, view, STACK_SRC); // pass 1 (detections->allSources)
    225270
    226271    // measure non-linear extended source models (exponential, deVaucouleur, Sersic) (sources sorted by S/N)
    227272    psphotExtendedSourceFits (config, view, STACK_SRC); // pass 1 (detections->allSources)
    228 
    229     // calculate source magnitudes
    230     psphotMagnitudes(config, view, STACK_SRC);
    231273
    232274    // create source children for the OUT filerule (for radial aperture photometry)
    … …  
    238280    }
    239281
     282    // measure circular, radial apertures (objects sorted by S/N)
     283    // XXX can we just use psphotRadialApertures
     284    // XXX make sure the headers are consistent with this (which PSF convolutions, ie mark 'none')
     285    // psphotRadialAperturesByObject (config, objectsRadial, view, STACK_SRC, nMatchedPSF);
     286    psphotRadialApertures (config, view, STACK_SRC, 0); // XXX entry 0 == unmatched?
    240287    psMemDump("extmeas");
    241288
    242     bool smoothAgain = true;
    243     for (int nMatchedPSF = 0; smoothAgain; nMatchedPSF++) {
     289    int nRadialEntries = psphotStackMatchPSFsEntries(config, view, STACK_OUT);
     290
     291    for (int entry = 1; entry < nRadialEntries; entry++) {
     292        // NOTE: entry 0 is the unmatched image set
    244293
    245294        // re-measure the PSF for the smoothed image (using entries in 'allSources')
    … …  
    253302
    254303        // measure circular, radial apertures (objects sorted by S/N)
    255         psphotRadialAperturesByObject (config, objectsRadial, view, STACK_OUT, nMatchedPSF);
     304        // entry 0 == unmatched? pass entry + 1?
     305        psphotRadialApertures (config, view, STACK_OUT, entry);
    256306
    257307        // replace the flux in the image so it is returned to its original state
    … …  
    259309
    260310        // smooth to the next FWHM, or set 'smoothAgain' to false if no more
    261         psphotStackMatchPSFsNext(&smoothAgain, config, view, STACK_OUT, nMatchedPSF);
     311        psphotStackMatchPSFsNext(config, view, STACK_OUT, entry);
    262312        psMemDump("matched");
    263313    }
    264314
    265     if (0 && !psphotEfficiency(config, view, STACK_OUT)) {
     315    // measure aperture photometry corrections
     316    if (!psphotApResid (config, view, STACK_SRC)) {
     317        psFree (objects);
     318        psLogMsg ("psphot", 3, "failed on psphotApResid");
     319        return psphotReadoutCleanup (config, view, STACK_SRC);
     320    }
     321
     322    // calculate source magnitudes
     323    psphotMagnitudes(config, view, STACK_SRC);
     324
     325    if (0 && !psphotEfficiency(config, view, STACK_DET)) {
    266326        psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources");
    267327        psErrorClear();
Note: See TracChangeset for help on using the changeset viewer.