IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Sep 6, 2011, 1:32:31 PM (15 years ago)
Author:
eugene
Message:
  • clarified the output to show the stages and their timing more clearly
  • psphotStackReadout now uses the same functions as psphotReadout (deprecated FitLinearStack, ExtendedAnalysisByObject, RadialAnalysisByObject, etc)
  • psphotStack now does 2nd (5 sigma) detection pass and radial photometry of the un-matched image (as well as the matched images)
  • for speed, fitting uses a smaller subset of pixels than subtraction: the model is extended when the subtraction is performed to avoid a discontinuity
  • extended source analysis (petrosians) is threaded
  • iterative kron radius / magnitude analysis with down-weighted neighbors replaces the masked neighbor kron analysis. the result of this analysis is used by fitting guesses and is the value written to the PSF and XFIT tables.
  • the target matched PSFs are determined in advance and the 0 element of this array is always the unmatched image (with target PSF of NAN)
  • dQ, the difference between the upper and lower quartile, is measured (but not yet used to adjust the sky level. sky should be FITTED_MEAN - 2.5*dQ
  • fake sources inserted for the efficiency analysis are removed at the end of the analysis
  • radial apertures is now threaded
  • edge case petrosian mags are now handled without aborts (may still have a problem with the measured error of the petrosian radius)
  • new sersic model guess based on the psf size, the moments, and the kron flux. the psf - kron mag is used to guess the index, then the index is used to convert measured moments into model guess parameters. this uses empirical relationships between the quantities generated from fake images.
  • define the fit region and window region based on the 1st radial moments
  • use the kron S/N to determine if analysis is performed for extended sources, not the psf S/N
  • remove the old 'basic deblend' which did not really do any deblending and is superceeded by the deblending in KronIterate and the footprint analysis
Location:
trunk/psphot
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot

  • trunk/psphot/src

    • Property svn:ignore

      •  

        old new  
        2222psphotMakePSF
        2323psphotStack
         24psphotModelTest

  • trunk/psphot/src/psphotStackReadout.c

    r31154 r32348  
    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.