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/psphotStackMatchPSFsUtils.c

    r31154 r32348  
    22# define ARRAY_BUFFER 16                 // Number to add to array at a time
    33
    4 // XXX better name
    5 bool readImage(psImage **target, // Target for image
    6                const char *name, // Name of FITS file
    7                const pmConfig *config // Configuration
    8     )
    9 {
    10     psString resolved = pmConfigConvertFilename(name, config, false, false); // Resolved filename
    11     psFits *fits = psFitsOpen(resolved, "r");
    12     psFree(resolved);
    13     if (!fits) {
    14         psError(PSPHOT_ERR_IO, false, "Unable to open previously produced image: %s", name);
    15         return false;
    16     }
    17     psImage *image = psFitsReadImage(fits, psRegionSet(0,0,0,0), 0); // Image of interest
    18     if (!image) {
    19         psError(PSPHOT_ERR_IO, false, "Unable to read previously produced image: %s", name);
    20         psFitsClose(fits);
    21         return false;
    22     }
    23     psFitsClose(fits);
    24 
    25     psFree(*target);
    26     *target = image;
    27 
    28     return true;
    29 }
     4psVector *SetOptWidths (bool *optimum, psMetadata *recipe);
     5pmReadout *makeFakeReadout(pmConfig *config, pmReadout *raw, psArray *sources, pmPSF *psf, psImageMaskType maskVal, int fullSize);
     6bool saveMatchData (pmReadout *readout, psphotStackOptions *options, int index);
     7bool matchKernel(pmConfig *config, pmReadout *cnv, pmReadout *raw, psphotStackOptions *options, int index);
     8bool dumpImageDiff(pmReadout *readoutConv, pmReadout *readoutFake, pmReadout *readoutRef, int index, char *rootname);
     9bool dumpImage(pmReadout *readoutOut, pmReadout *readoutRef, int index, char *rootname);
    3010
    3111// Get coordinates from a source
    … …  
    148128
    149129// Renormalise a readout's variance map
    150 bool stackRenormaliseReadout(const pmConfig *config, // Configuration
    151                              pmReadout *readout      // Readout to renormalise
     130bool psphotStackRenormaliseVariance(const pmConfig *config, // Configuration
     131                              pmReadout *readout      // Readout to renormalise
    152132    )
    153133{
    … …  
    180160    return pmReadoutVarianceRenormalise(readout, maskBad, num, minValid, maxValid);
    181161}
    182 
    183 // This is a hack to use the temporary convolved images and kernel generated previously.
    184 // This makes the 'matching' operation much faster, allowing debugging of the stack process easier.
    185 // It implicitly assumes the output root name is the same between invocations.
    186 
    187 # if (0)
    188 bool loadKernel (pmConfig *config, pmReadout *readoutCnv, psphotStackOptions *options, int index) {
    189 
    190     // Read the convolution kernel from the saved file
    191     pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.CONV.KERNEL", index);
    192     psAssert(file, "Require file");
    193 
    194     pmFPAview *view = pmFPAviewAlloc(0); // View to readout of interest
    195     view->chip = view->cell = view->readout = 0;
    196     psString filename = pmFPAfileNameFromRule(file->filerule, file, view); // Filename of interest
    197 
    198     // Read convolution kernel data
    199     psString resolved = pmConfigConvertFilename(filename, config, false, false); // Resolved filename
    200     psFree(filename);
    201     psFits *fits = psFitsOpen(resolved, "r"); // FITS file for subtraction kernel
    202     psFree(resolved);
    203     if (!fits || !pmReadoutReadSubtractionKernels(readoutCnv, fits)) {
    204         psError(PSPHOT_ERR_IO, false, "Unable to read previously produced kernel");
    205         psFitsClose(fits);
    206         return false;
    207     }
    208     psFitsClose(fits);
    209 
    210     // read the convolved pixels (image, mask, variance) -- names are pre-defined
    211     if (!readImage(&readoutCnv->image,    options->convImages->data[index],    config) ||
    212         !readImage(&readoutCnv->mask,     options->convMasks->data[index],     config) ||
    213         !readImage(&readoutCnv->variance, options->convVariances->data[index], config)) {
    214         psError(PSPHOT_ERR_IO, false, "Unable to read previously produced image.");
    215         return false;
    216     }
    217 
    218     // XXX ??? not sure what is happening here -- consult Paul
    219     psRegion *region = psMetadataLookupPtr(NULL, readoutCnv->analysis, PM_SUBTRACTION_ANALYSIS_REGION); // Convolution region
    220     pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, readoutCnv->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL);
    221 
    222     pmSubtractionAnalysis(readoutCnv->analysis, NULL, kernels, region, readoutCnv->image->numCols, readoutCnv->image->numRows);
    223 
    224     psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel
    225 
    226     // update the covariance matrix
    227     // XXX why is this needed if we have correctly read the saved data?
    228     bool oldThreads = psImageCovarianceSetThreads(true);              // Old thread setting
    229     psKernel *covar = psImageCovarianceCalculate(kernel, readoutCnv->covariance); // Covariance matrix
    230     psImageCovarianceSetThreads(oldThreads);
    231     psFree(readoutCnv->covariance);
    232     readoutCnv->covariance = covar;
    233     psFree(kernel);
    234     return true;
    235 }
    236 # endif
    237162
    238163bool dumpImage(pmReadout *readoutOut, pmReadout *readoutRef, int index, char *rootname) {
    … …  
    362287        widthsCopy = psVectorCopy(NULL, widths, PS_TYPE_F32); // Copy of kernel widths
    363288
    364         // we need to register the FWHM values for use downstream
    365         pmSubtractionSetFWHMs(options->targetSeeing, options->inputSeeing->data.F32[index]);
     289        // we need to register the FWHM values for use by pmSubtraction code
     290        pmSubtractionSetFWHMs(options->targetSeeing->data.F32[0], options->inputSeeing->data.F32[index]);
    366291
    367292        pmSubtractionParamScaleOptions(scale, scaleRef, scaleMin, scaleMax);
    … …  
    469394}
    470395
    471 // Kernel normalisation for convolved readout
    472 bool renormKernel(pmReadout *readout, psphotStackOptions *options, int index) {
    473 
    474     double sum = 0.0;           // Sum of chi^2
    475     int num = 0;                // Number of measurements of chi^2
    476     psString regex = NULL;      // Regular expression
    477     psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_NORM);
    478     psMetadataIterator *iter = psMetadataIteratorAlloc(readout->analysis, PS_LIST_HEAD, regex);
    479     psFree(regex);
    480     psMetadataItem *item = NULL;// Item from iteration
    481     while ((item = psMetadataGetAndIncrement(iter))) {
    482         assert(item->type == PS_TYPE_F32);
    483         float norm = item->data.F32; // Normalisation
    484         sum += norm;
    485         num++;
    486     }
    487     psFree(iter);
    488     float conv = sum/num;       // Mean normalisation from convolution
    489     float stars = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation from stars
    490     float renorm =  stars / conv; // Renormalisation to apply
    491     psLogMsg("psphotStack", PS_LOG_INFO, "Renormalising image %d by %f (kernel: %f, stars: %f)\n", index, renorm, conv, stars);
    492 
    493     psBinaryOp(readout->image, readout->image, "*", psScalarAlloc(renorm, PS_TYPE_F32));
    494     psBinaryOp(readout->variance, readout->variance, "*", psScalarAlloc(PS_SQR(renorm), PS_TYPE_F32));
    495     return true;
    496 }
    497 
    498 // adjust scaling for readout (remove background, ..., determine weighting)
    499 bool rescaleData(pmReadout *readout, pmConfig *config, psphotStackOptions *options, int index) {
    500 
    501     psMetadata *stackRecipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSTACK"); // ppStack recipe
    502     psAssert(stackRecipe, "We've thrown an error on this before.");
    503 
    504     // Look up appropriate values from the ppSub recipe
    505     psMetadata *subRecipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe
    506     psAssert(subRecipe, "recipe missing");
    507 
    508     psString maskValStr = psMetadataLookupStr(NULL, subRecipe, "MASK.VAL"); // Name of bits to mask going in
    509     psString maskBadStr = psMetadataLookupStr(NULL, stackRecipe, "MASK.BAD"); // Name of bits to mask for bad
    510 
    511     psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
    512     psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
    513 
    514     // Ensure the background value is zero
    515     psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for background
    516     psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
    517 
    518     // XXX why is this in config->arguments and not recipe?
    519     if (!psMetadataLookupBool(NULL, config->arguments, "PPSTACK.SKIP.BG.SUB")) {
    520         if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) {
    521             psAbort("Can't measure background for image.");
    522             // XXX we used to clear error: why is this acceptable? psErrorClear();
    523         }
    524 
    525         float value = psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN);
    526         float stdev = psStatsGetValue(bg, PS_STAT_ROBUST_STDEV);
    527 
    528         psLogMsg("psphotStack", PS_LOG_INFO, "Correcting convolved image background by %lf (+/- %lf)", value, stdev);
    529         psBinaryOp(readout->image, readout->image, "-", psScalarAlloc(value, PS_TYPE_F32));
    530     }
    531 
    532     if (!stackRenormaliseReadout(config, readout)) {
    533         psFree(rng);
    534         psFree(bg);
    535         return false;
    536     }
    537 
    538     // Measure the variance level for the weighting
    539     if (psMetadataLookupBool(NULL, stackRecipe, "WEIGHTS")) {
    540         if (!psImageBackground(bg, NULL, readout->variance, readout->mask, maskVal | maskBad, rng)) {
    541             psError(PSPHOT_ERR_DATA, false, "Can't measure mean variance for image.");
    542             psFree(rng);
    543             psFree(bg);
    544             return false;
    545         }
    546         options->weightings->data.F32[index] = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) * psImageCovarianceFactor(readout->covariance));
    547     } else {
    548         options->weightings->data.F32[index] = 1.0;
    549     }
    550     psLogMsg("psphotStack", PS_LOG_INFO, "Weighting for image %d is %f\n", index, options->weightings->data.F32[index]);
    551 
    552     psFree(rng);
    553     psFree(bg);
    554     return true;
    555 }
    556 
    557396# define NOISE_FRACTION 0.01             // Set minimum flux to this fraction of noise
    558397# define SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to input sources
Note: See TracChangeset for help on using the changeset viewer.