IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 32312


Ignore:
Timestamp:
Sep 5, 2011, 8:47:07 AM (15 years ago)
Author:
eugene
Message:

targetSeeing is now a vector; clean up and remove ppStack kruft unrelated to psphotStack (no photometric renorm, no input rejections)

File:
1 edited

Legend:

Unmodified
Added
Removed

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

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