IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Feb 1, 2009, 11:55:34 AM (17 years ago)
Author:
eugene
Message:

re-factored ppSubReadout into a sequence of logical steps; defined convolved images as their own pmFPAfiles; switched to new psphot functions (psphotReadoutFindPSF and psphotReadoutMinimal); simplified mask and psphot visualization handling

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppSub/src/ppSubReadout.c

    r21183 r21257  
    1 #ifdef HAVE_CONFIG_H
    2 #include <config.h>
    3 #endif
    4 
    5 #include <stdio.h>
    6 #include <string.h>
    7 #include <pslib.h>
    8 #include <psmodules.h>
    9 #include <psphot.h>
    10 
    111#include "ppSub.h"
    12 
    13 #define WCS_TOLERANCE 0.001             // Tolerance for WCS
    14 //#define TESTING                         // For test output
    15 
    16 #define SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_SATSTAR | PM_SOURCE_MODE_BLEND | \
    17                      PM_SOURCE_MODE_BADPSF | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | \
    18                      PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to input sources
    19 
    20 
    21 // Copy every instance of a single keyword from one metadata to another
    22 static bool metadataCopySingle(psMetadata *target, psMetadata *source, const char *name)
    23 {
    24     PS_ASSERT_METADATA_NON_NULL(target, false);
    25     PS_ASSERT_METADATA_NON_NULL(source, false);
    26     PS_ASSERT_STRING_NON_EMPTY(name, false);
    27 
    28     psString regex = NULL;      // Regular expression
    29     psStringAppend(&regex, "^%s$", name);
    30     psMetadataIterator *iter = psMetadataIteratorAlloc(source, PS_LIST_HEAD, regex); // Iterator
    31     psFree(regex);
    32     psMetadataItem *item;       // Item from iteration
    33     while ((item = psMetadataGetAndIncrement(iter))) {
    34         psMetadataAddItem(target, item, PS_LIST_TAIL, PS_META_DUPLICATE_OK);
    35     }
    36     psFree(iter);
    37 
    38     return true;
    39 }
    40 
    412
    423bool ppSubReadout(pmConfig *config, psMetadata *stats, const pmFPAview *view)
    434{
    445    psTimerStart("PPSUB_MATCH");
    45     pmReadout *inRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT"); // Input readout
    46     pmReadout *refRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF"); // Reference readout
    47     pmReadout *sourcesRO = pmFPAfileThisReadout(config->files, view, "PPSUB.SOURCES"); // Readout with sources
    48     pmReadout *inConv = pmReadoutAlloc(NULL); // Convolved version of input
    49     pmReadout *refConv = pmReadoutAlloc(NULL); // Convolved version of reference
    50     pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSUB.OUTPUT"); // Output cell
    51     pmReadout *outRO = pmReadoutAlloc(outCell); // Output readout: subtraction
    52     pmFPA *outFPA = outCell->parent->parent; // Output FPA
    53     pmHDU *outHDU = outFPA->hdu; // Output HDU
    54     if (!outHDU->header) {
    55         outHDU->header = psMetadataAlloc();
    56     }
    576
    58     psImage *input = inRO->image;       // Input image
    59     psImage *reference = refRO->image;  // Reference image
    60     PS_ASSERT_IMAGES_SIZE_EQUAL(input, reference, false);
    61 
    62     int threads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads
    63     if (threads > 0) {
    64         if (!psThreadPoolInit(threads)) {
    65             psError(PS_ERR_UNKNOWN, false, "Unable to setup %d threads", threads);
    66             return false;
    67         }
    68         pmSubtractionThreadsInit(inRO, refRO);
    69     }
    70 
    71     // Look up recipe values
    72     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSim
    73     psAssert(recipe, "We checked this earlier, so it should be here.");
    74 
    75     const char *typeStr = psMetadataLookupStr(NULL, recipe, "KERNEL.TYPE"); // Kernel type
    76     psAssert(typeStr, "We put it here in ppSubArguments.c");
    77     pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(typeStr); // Type of kernel
    78     if (type == PM_SUBTRACTION_KERNEL_NONE) {
    79         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unrecognised kernel type: %s", typeStr);
    80         psFree(inConv);
    81         psFree(refConv);
    82         psFree(outRO);
     7    if (!ppSubSetMasks (config, view)) {
     8        psError(PS_ERR_UNKNOWN, false, "Unable to subtract background.");
    839        return false;
    8410    }
    8511
    86     bool mdok;                          // Status of MD lookup
    87     int size = psMetadataLookupS32(NULL, recipe, "KERNEL.SIZE"); // Kernel half-size
    88     int order = psMetadataLookupS32(NULL, recipe, "SPATIAL.ORDER"); // Spatial polynomial order
    89     float regionSize = psMetadataLookupF32(NULL, recipe, "REGION.SIZE"); // Size of iso-kernel regs
    90     float spacing = psMetadataLookupF32(NULL, recipe, "STAMP.SPACING"); // Typical stamp spacing
    91     int footprint = psMetadataLookupS32(NULL, recipe, "STAMP.FOOTPRINT"); // Stamp half-size
    92     float threshold = psMetadataLookupF32(NULL, recipe, "STAMP.THRESHOLD"); // Threshold for stmps
    93     int stride = psMetadataLookupS32(NULL, recipe, "STRIDE"); // Size of convolution patches
    94     int iter = psMetadataLookupS32(NULL, recipe, "ITER"); // Rejection iterations
    95     float rej = psMetadataLookupF32(NULL, recipe, "REJ"); // Rejection threshold
    96     float sys = psMetadataLookupF32(NULL, recipe, "SYS"); // Relative systematic error in kernel
    97     bool reverse = psMetadataLookupBool(NULL, config->arguments, "REVERSE"); // Reverse sense of subtraction?
    98     psVector *widths = psMetadataLookupPtr(NULL, recipe, "ISIS.WIDTHS"); // ISIS Gaussian widths
    99     psVector *orders = psMetadataLookupPtr(NULL, recipe, "ISIS.ORDERS"); // ISIS Polynomial orders
    100     int inner = psMetadataLookupS32(NULL, recipe, "INNER"); // Inner radius
    101     int ringsOrder = psMetadataLookupS32(NULL, recipe, "RINGS.ORDER"); // RINGS polynomial order
    102     int binning = psMetadataLookupS32(NULL, recipe, "SPAM.BINNING"); // Binning for SPAM kernel
    103     float penalty = psMetadataLookupF32(NULL, recipe, "PENALTY"); // Penalty for wideness
    104     psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.IN"); // Name of bits to mask going in
    105     psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
    106     psString maskPoorStr = psMetadataLookupStr(NULL, recipe, "MASK.POOR"); // Name of bits to mask for poor
    107     psImageMaskType maskPoor = pmConfigMaskGet(maskPoorStr, config); // Bits to mask for poor pixels
    108     psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad
    109     psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
    110     float badFrac = psMetadataLookupF32(NULL, recipe, "BADFRAC"); // Maximum bad fraction
    111     const char *stampsName = psMetadataLookupStr(&mdok, config->arguments, "STAMPS"); // Filename for stamps
    112 
    113     bool optimum = psMetadataLookupBool(&mdok, recipe, "OPTIMUM"); // Derive optimum parameters?
    114     float optMin = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.MIN"); // Minimum width for search
    115     float optMax = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.MAX"); // Maximum width for search
    116     float optStep = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.STEP"); // Step for search
    117     float optThresh = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.TOL"); // Tolerance for search
    118     int optOrder = psMetadataLookupS32(&mdok, recipe, "OPTIMUM.ORDER"); // Order for search
    119     bool dual = psMetadataLookupBool(&mdok, recipe, "DUAL"); // Dual convolution?
    120 
    121     // Statistics for renormalisation
    122     psStatsOptions renormMean = psStatsOptionFromString(psMetadataLookupStr(NULL, recipe, "RENORM.MEAN"));
    123     psStatsOptions renormStdev = psStatsOptionFromString(psMetadataLookupStr(NULL, recipe, "RENORM.STDEV"));
    124     if (renormMean == PS_STAT_NONE || renormStdev == PS_STAT_NONE) {
    125         psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    126                 "Unable to parse renormalisation statistics from recipe.");
    127         return false;
    128     }
    129     int renormNum = psMetadataLookupS32(&mdok, recipe, "RENORM.NUM"); // Number of samples for renormalisation
    130     float renormWidth = psMetadataLookupS32(&mdok, recipe, "RENORM.WIDTH"); // Width of Gaussian phot
    131 
    132     psString interpModeStr = psMetadataLookupStr(&mdok, recipe, "INTERPOLATION"); // Interpolation mode
    133     psImageInterpolateMode interpMode = psImageInterpolateModeFromString(interpModeStr); // Interp
    134     if (interpMode == PS_INTERPOLATE_NONE) {
    135         psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unknown interpolation mode: %s", interpModeStr);
    136         return false;
    137     }
    138     float poorFrac = psMetadataLookupF32(&mdok, recipe, "POOR.FRACTION"); // Fraction for "poor"
    139 
    140     pmSubtractionMode subMode = dual ? PM_SUBTRACTION_MODE_DUAL : PM_SUBTRACTION_MODE_UNSURE; // Subtracn mode
    141 
    142     // Generate masks if they don't exist
    143     int numCols = input->numCols, numRows = input->numRows; // Image dimensions
    144     if (!inRO->mask) {
    145         if (psMetadataLookupBool(NULL, recipe, "MASK.GENERATE")) {
    146             pmReadoutSetMask(inRO, pmConfigMaskGet("SAT", config), pmConfigMaskGet("BAD", config));
    147         } else {
    148             inRO->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
    149             psImageInit(inRO->mask, 0);
    150         }
    151     }
    152     if (!refRO->mask) {
    153         if (psMetadataLookupBool(NULL, recipe, "MASK.GENERATE")) {
    154             pmReadoutSetMask(refRO, pmConfigMaskGet("SAT", config), pmConfigMaskGet("BAD", config));
    155         } else {
    156             refRO->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
    157             psImageInit(refRO->mask, 0);
    158         }
    159     }
    160 
    161     if (!pmReadoutMaskNonfinite(inRO, pmConfigMaskGet("SAT", config))) {
    162         psError(PS_ERR_UNKNOWN, false, "Unable to mask non-finite pixels in input.");
    163         return false;
    164     }
    165     if (!pmReadoutMaskNonfinite(refRO, pmConfigMaskGet("SAT", config))) {
    166         psError(PS_ERR_UNKNOWN, false, "Unable to mask non-finite pixels in reference.");
     12    if (!ppSubMatchPSFs (config, view)) {
     13        psError(PS_ERR_UNKNOWN, false, "Unable to subtract background.");
    16714        return false;
    16815    }
    16916
    170     psVector *optWidths = NULL;         // Vector with FWHMs for optimum search
    171     if (optimum) {
    172         optWidths = psVectorCreate(optWidths, optMin, optMax, optStep, PS_TYPE_F32);
    173     }
    174 
    175     psArray *sources = NULL;            // Sources in image
    176     if (sourcesRO) {
    177         sources = psMetadataLookupPtr(&mdok, sourcesRO->analysis, "PSPHOT.SOURCES");
    178     }
    179 
    180 #if 0
    181     // Interpolate over bad pixels, so the bad pixels don't explode
    182     if (!pmReadoutInterpolateBadPixels(inRO, maskVal, interpMode, poorFrac, maskPoor, maskBad)) {
    183         psError(PS_ERR_UNKNOWN, false, "Unable to interpolate bad pixels for input image.");
    184         return false;
    185     }
    186     if (!pmReadoutInterpolateBadPixels(refRO, maskVal, interpMode, poorFrac, maskPoor, maskBad)) {
    187         psError(PS_ERR_UNKNOWN, false, "Unable to interpolate bad pixels for reference image.");
    188         return false;
    189     }
    190     maskVal |= maskBad;
    191 #endif
    192 
    193     // Match the PSFs
    194     if (!pmSubtractionMatch(inConv, refConv, inRO, refRO, footprint, stride, regionSize, spacing, threshold,
    195                             sources, stampsName, type, size, order, widths, orders, inner, ringsOrder,
    196                             binning, penalty, optimum, optWidths, optOrder, optThresh, iter, rej, sys,
    197                             maskVal, maskBad, maskPoor, poorFrac, badFrac, subMode)) {
    198         psError(PS_ERR_UNKNOWN, false, "Unable to match images.");
    199         psFree(inConv);
    200         psFree(refConv);
    201         psFree(outRO);
    202         return false;
    203     }
    204     psFree(optWidths);
    205 
    206     pmSubtractionThreadsFinalize(inRO, refRO);
    207 
    208     // Add kernel descrption to header
    209     pmSubtractionKernels *kernels = psMetadataLookupPtr(&mdok, inConv->analysis,
    210                                                         PM_SUBTRACTION_ANALYSIS_KERNEL); // Subtraction kernel
    211     if (!kernels) {
    212         kernels = psMetadataLookupPtr(&mdok, refConv->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL);
    213     }
    214     if (!kernels) {
    215         psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find SUBTRACTION.KERNEL");
    216         psFree(inConv);
    217         psFree(refConv);
    218         psFree(outRO);
    219         return false;
    220     }
    221     psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "PPSUB.KERNEL", 0,
    222                      "Subtraction kernel", kernels->description);
    223 
    224     {
    225         if (psMetadataLookup(inConv->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL)) {
    226             outRO->analysis = psMetadataCopy(outRO->analysis, inConv->analysis);
    227         } else if (psMetadataLookup(refConv->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL)) {
    228             outRO->analysis = psMetadataCopy(outRO->analysis, refConv->analysis);
    229         } else {
    230             psWarning("Unable to find subtraction kernel in analysis metadata.");
    231         }
    232 
    233         // Update variance factors
    234         // It's not possible to do this perfectly, since we're combining different images:
    235         // vf_sum sigma_sum^2 = vf_1 * sigma_1^2 + vf_2 * sigma_2^2
    236         // It's not possible to write vf_sum as a function of vf_1 and vf_2 with no reference to the sigmas.
    237         // Instead, we're going to cheat.
    238         bool mdok;                      // Status of MD lookup
    239         pmSubtractionKernels *kernels = psMetadataLookupPtr(&mdok, outRO->analysis,
    240                                                             PM_SUBTRACTION_ANALYSIS_KERNEL); // Kernels
    241         if (!mdok) {
    242             psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find subtraction kernels.");
    243             psFree(inConv);
    244             psFree(refConv);
    245             psFree(outRO);
    246             return false;
    247         }
    248         float vfIn = psMetadataLookupF32(NULL, inRO->parent->concepts,
    249                                          "CELL.VARFACTOR"); // Variance factor for input
    250         if (!isfinite(vfIn)) {
    251             vfIn = 1.0;
    252         }
    253         float vfRef = psMetadataLookupF32(NULL, refRO->parent->concepts,
    254                                           "CELL.VARFACTOR"); // Variance factor for ref
    255         if (!isfinite(vfRef)) {
    256             vfRef = 1.0;
    257         }
    258         if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
    259             vfIn *= pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false);
    260         }
    261         if (kernels->mode == PM_SUBTRACTION_MODE_2) {
    262             vfRef *= pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false);
    263         } else if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
    264             vfRef *= pmSubtractionVarianceFactor(kernels, 0.0, 0.0, true);
    265         }
    266 
    267         psStats *vfStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics
    268         psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator
    269         if (!psImageBackground(vfStats, NULL, inConv->weight, inConv->mask, maskVal | maskBad, rng)) {
    270             psError(PS_ERR_UNKNOWN, false, "Unable to measure mean variance for convolved input");
    271             psFree(inConv);
    272             psFree(refConv);
    273             psFree(outRO);
    274             psFree(vfStats);
    275             psFree(rng);
    276             return false;
    277         }
    278         float inMeanVar = vfStats->robustMedian; // Mean variance of input
    279         if (!psImageBackground(vfStats, NULL, refConv->weight, refConv->mask, maskVal | maskBad, rng)) {
    280             psError(PS_ERR_UNKNOWN, false, "Unable to measure mean variance for convolved reference");
    281             psFree(inConv);
    282             psFree(refConv);
    283             psFree(outRO);
    284             psFree(vfStats);
    285             psFree(rng);
    286             return false;
    287         }
    288         float refMeanVar = vfStats->robustMedian; // Mean variance of reference
    289         psFree(rng);
    290         psFree(vfStats);
    291 
    292         float vf = (vfIn * inMeanVar + vfRef * refMeanVar) / (inMeanVar + refMeanVar); // Resulting var factor
    293         psMetadataItem *vfItem = psMetadataLookup(outRO->parent->concepts, "CELL.VARFACTOR");
    294         vfItem->data.F32 = vf;
    295 
    296         // Statistics on the matching
    297         if (stats) {
    298             metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MODE);
    299             metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_STAMPS);
    300             metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_DEV_MEAN);
    301             metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_DEV_RMS);
    302             metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_NORM);
    303             metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_BGDIFF);
    304             metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MX);
    305             metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MY);
    306             metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MXX);
    307             metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MXY);
    308             metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MYY);
    309 
    310             psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_MATCH", 0, "Time to match PSFs",
    311                              psTimerClear("PPSUB_MATCH"));
    312         }
    313     }
    314 
    315 
    316 #ifdef TESTING
    317     psImage *kernelImage = psMetadataLookupPtr(&mdok, inConv->analysis,
    318                                                "SUBTRACTION.KERNEL.IMAGE"); // Image of the kernels
    319     if (!kernelImage) {
    320         kernelImage = psMetadataLookupPtr(&mdok, refConv->analysis, "SUBTRACTION.KERNEL.IMAGE");
    321     }
    322     psFits *fits = psFitsOpen("kernel.fits", "w");
    323     psFitsWriteImage(fits, NULL, kernelImage, 0, NULL);
    324     psFitsClose(fits);
    325 #endif
    326 
    327     // Subtraction is: minuend - subtrahend
    328     pmReadout *minuend = inConv;
    329     pmReadout *subtrahend = refConv;
    330 
    331     if (reverse) {
    332         pmReadout *temp = subtrahend;
    333         subtrahend = minuend;
    334         minuend = temp;
    335     }
    336 
    337 #ifdef TESTING
    338     {
    339         pmReadoutMaskApply(minuend, maskVal);
    340         psFits *fits = psFitsOpen("minuend.fits", "w");
    341         psFitsWriteImage(fits, NULL, minuend->image, 0, NULL);
    342         psFitsClose(fits);
    343     }
    344     {
    345         pmReadoutMaskApply(subtrahend, maskVal);
    346         psFits *fits = psFitsOpen("subtrahend.fits", "w");
    347         psFitsWriteImage(fits, NULL, subtrahend->image, 0, NULL);
    348             psFitsClose(fits);
    349     }
    350 #endif
    351 
    352     // Photometry is to be performed in two stages:
    353     // 1. Measure the PSF using the PSF-matched images
    354     // 2. Find and measure sources on the subtracted image
    355     psTimerStart("PPSUB_PHOT");
    356 
    357     // Photometry stage 1: measure the PSF
    358     pmPSF *psf = NULL;                  // PSF for photometry
    359     if (psMetadataLookupBool(NULL, recipe, "PHOTOMETRY")) {
    360         outRO->image = psImageCopy(outRO->image, minuend->image, PS_TYPE_F32);
    361         if (minuend->weight) {
    362             outRO->weight = psImageCopy(outRO->weight, minuend->weight, PS_TYPE_F32);
    363         }
    364         if (minuend->mask) {
    365             outRO->mask = psImageCopy(outRO->mask, minuend->mask, PS_TYPE_IMAGE_MASK);
    366         }
    367         outRO->data_exists = outCell->data_exists = outCell->parent->data_exists = true;
    368 
    369         if (psMetadataLookupBool(&mdok, recipe, "RENORM")) {
    370             psImageMaskType maskValue = pmConfigMaskGet("BLANK", config); // Bits to mask
    371             if (!pmReadoutWeightRenormPhot(outRO, maskValue, renormNum, renormWidth,
    372                                            renormMean, renormStdev, NULL)) {
    373                 psError(PS_ERR_UNKNOWN, false, "Unable to renormalise variances.");
    374                 psFree(outRO);
    375                 return false;
    376             }
    377         }
    378 
    379         pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT");
    380         pmFPACopy(photFile->fpa, outRO->parent->parent->parent);
    381 
    382         // We have a list of sources, so we could use those to redetermine the PSF model.
    383         // Until Gene makes the necessary adaptations to psphot, we will simply redetermine the PSF model from
    384         // scratch.
    385 
    386         if (psMetadataLookupBool(&mdok, recipe, "PSPHOT.VISUAL")) {
    387             psphotSetVisual(true);
    388         }
    389 
    390         // Need to ensure aperture residual is not calculated
    391         psMetadata *psphotRecipe = psMetadataLookupMetadata(NULL, config->recipes, PSPHOT_RECIPE); // Recipe
    392         if (!psphotRecipe) {
    393             psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find %s recipe.", PSPHOT_RECIPE);
    394             return false;
    395         }
    396         const char *breakpoint = psMetadataLookupStr(&mdok, psphotRecipe, "BREAK_POINT"); // Current break point
    397         psMetadataAddStr(psphotRecipe, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE,
    398                          "ALTERED break point for psphot operations", "PSFMODEL");
    399 
    400         // set maskValue and markValue in the psphot recipe
    401         psImageMaskType maskValue = maskVal;
    402         psImageMaskType markValue = pmConfigMaskGet("MARK.VALUE", config); // Bits to use for marking
    403         psMetadataAddImageMask(psphotRecipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "Bits to mask", maskValue);
    404         psMetadataAddImageMask(psphotRecipe, PS_LIST_TAIL, "MARK.PSPHOT", PS_META_REPLACE, "Bits to use for marking",
    405                         markValue);
    406 
    407         if (!psphotReadout(config, view)) {
    408             psWarning("Unable to perform photometry on subtracted image.");
    409             psErrorStackPrint(stderr, "Error stack from photometry:");
    410             psErrorClear();
    411         }
    412 
    413         if (!pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL") ||
    414             !pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV") ||
    415             !pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND")) {
    416             psError(PS_ERR_UNKNOWN, false, "Unable to drop PSPHOT internal files.");
    417             return false;
    418         }
    419 
    420         if (breakpoint) {
    421             psMetadataAddStr(psphotRecipe, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE,
    422                              "RESTORED break point for psphot operations", breakpoint);
    423         }
    424 
    425         // Blow away the sources psphot found --- they're irrelevant for the subtraction
    426         pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout with sources
    427         psMetadataRemoveKey(photRO->analysis, "PSPHOT.SOURCES");
    428         psMetadataRemoveKey(photRO->analysis, "PSPHOT.HEADER");
    429 
    430         psf = psMemIncrRefCounter(psMetadataLookupPtr(NULL, photRO->parent->parent->analysis, "PSPHOT.PSF"));
    431         if (!psf) {
    432             psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find PSF from psphot");
    433             return false;
    434         }
    435 
    436         // Record the FWHM
    437         pmHDU *hdu = pmHDUFromCell(outRO->parent); // HDU with header
    438         psMetadataItemSupplement(hdu->header, psphotRecipe, "FWHM_MAJ");
    439         psMetadataItemSupplement(hdu->header, psphotRecipe, "FWHM_MIN");
    440 
    441 
    442         pmCell *photCell = pmFPAfileThisCell(config->files, view, "PSPHOT.INPUT");
    443         pmCellFreeReadouts(photCell);
    444     }
    445 
    446     // Do the actual subtraction
    447     outRO->image = (psImage*)psBinaryOp(outRO->image, minuend->image, "-", subtrahend->image);
    448     if (inConv->weight && refConv->weight) {
    449         outRO->weight = (psImage*)psBinaryOp(outRO->weight, inConv->weight, "+", refConv->weight);
    450     }
    451     outRO->mask = (psImage*)psBinaryOp(outRO->mask, inConv->mask, "|", refConv->mask);
    452 
    453     psFree(inConv);
    454     psFree(refConv);
    455 
    456     outRO->data_exists = outCell->data_exists = outCell->parent->data_exists = true;
    457 
    458     // Copy astrometry over
    459     // It should get into the output images and photometry
    460     pmFPA *refFPA = refRO->parent->parent->parent; // Reference FPA
    461     pmHDU *refHDU = refFPA->hdu;        // Reference HDU
    462     if (!outHDU || !refHDU) {
    463         psWarning("Unable to find HDU at FPA level to copy astrometry.");
    464     } else if (!pmAstromReadWCS(outFPA, outCell->parent, refHDU->header, 1.0)) {
    465         psWarning("Unable to read WCS astrometry from reference FPA.");
    466         psErrorClear();
    467     } else if (!pmAstromWriteWCS(outHDU->header, outFPA, outCell->parent, WCS_TOLERANCE)) {
    468         psWarning("Unable to write WCS astrometry to output FPA.");
    469         psErrorClear();
    470     }
    471 
    472 #if 0
    473     pmReadoutMaskApply(outRO, maskBad);
    474 #endif
    475 
    476 #ifdef TESTING
    477     for (int y = 0; y < outRO->image->numRows; y++) {
    478         for (int x = 0; x < outRO->image->numCols; x++) {
    479             if (isnan(outRO->image->data.F32[y][x]) && !(outRO->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal)) {
    480                 printf("Unmasked NAN at %d %d --> %d\n", x, y, outRO->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x]);
    481             }
    482         }
    483     }
    484 #endif
    485 
    486     if (!pmFPACopyConcepts(outCell->parent->parent, inRO->parent->parent->parent)) {
    487         psError(PS_ERR_UNKNOWN, false, "Unable to copy concepts from input to output.");
    488         psFree(outRO);
     17    if (!ppSubDefineOutput (config, view)) {
     18        psError(PS_ERR_UNKNOWN, false, "Unable to subtract background.");
    48919        return false;
    49020    }
    49121
    492 #if 0
    493     // XXX Under development
     22    if (!ppSubVarianceFactors (config, stats, view)) {
     23        psError(PS_ERR_UNKNOWN, false, "Unable to subtract background.");
     24        return false;
     25    }
    49426
    495     // QA: photometry of known sources with measured PSF
    496     if (sourcesRO && psMetadataLookupBool(NULL, recipe, "PHOTOMETRY")) {
    497         sources = psMetadataLookupPtr(&mdok, sourcesRO->analysis, "PSPHOT.SOURCES");
     27    if (!ppSubMakePSF(config, view)) {
     28        psError(PS_ERR_UNKNOWN, false, "Unable to subtract background.");
     29        return false;
     30    }
    49831
    499 
    500 
    501         // For each source, need to:
    502         // * generate appropriate model
    503         // * select pixels
    504         // * define fit radius
    505         // Don't need to measure sky (psphotFitSourcesLinear does that)
    506 
    507         psArray *dummy = psArrayAlloc(1); // Dummy array of sources, for psphotFitSourcesLinear with 1 source
    508         for (long i = 0; i < sources->n; i++) {
    509             pmSource *source = sources->data[i];
    510             if (!source || (source->mode & SOURCE_MASK)) {
    511                 continue;
    512             }
    513 
    514             float origMag = source->psfMag; // Original magnitude
    515 
    516             // Coordinates of source
    517             float x = source->modelPSF->params->data.F32[PM_PAR_XPOS];
    518             float y = source->modelPSF->params->data.F32[PM_PAR_YPOS];
    519 
    520             model = pmModelFromPSFforXY(psf, x, y, NAN);
    521             if (!fakeModel) {
    522                 psWarning("Can't generate model");
    523                 continue;
    524             }
    525             if (!pmSourceDefinePixels(fakeSource, readout, x, y, fakeRadius)) {
    526                 psErrorClear();
    527                 continue;
    528             }
    529 pmModel *pmModelFromPSFforXY (
    530     const pmPSF *psf,
    531     float Xo,
    532     float Yo,
    533     float Io
    534     );
    535 
    536 bool pmSourceDefinePixels(
    537     pmSource *mySource,                 ///< source to be re-defined
    538     const pmReadout *readout,  ///< base the source on this readout
    539     psF32 x,                            ///< center coords of source
    540     psF32 y,                            ///< center coords of source
    541     psF32 Radius                        ///< size of box on source
    542 );
    543 
    544 
    545 bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final) {
    546 
    547     }
    548 #endif
    549 
    550     // Photometry stage 2: find and measure sources on the subtracted image
    551     if (psMetadataLookupBool(NULL, recipe, "PHOTOMETRY")) {
    552         // The PSF should already be stored for the readout
    553         pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT");
    554         pmFPACopy(photFile->fpa, outFPA);
    555 
    556         pmReadout *psfRO = pmFPAfileThisReadout(config->files, view, "PSPHOT.PSF.LOAD");
    557         if (!psfRO) {
    558             psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find readout for file PSPHOT.PSF.LOAD");
    559             return false;
    560         }
    561         psMetadataAddPtr(psfRO->parent->parent->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN,
    562                          "PSF from matched addition", psf);
    563         psFree(psf);
    564 
    565         if (psMetadataLookupBool(&mdok, recipe, "PSPHOT.VISUAL")) {
    566             psphotSetVisual(true);
    567         }
    568 
    569         if (psMetadataLookupBool(&mdok, recipe, "RENORM")) {
    570             psImageMaskType maskValue = pmConfigMaskGet("BLANK", config); // Bits to mask
    571             if (!pmReadoutWeightRenormPhot(outRO, maskValue, renormNum, renormWidth,
    572                                            renormMean, renormStdev, NULL)) {
    573                 psError(PS_ERR_UNKNOWN, false, "Unable to renormalise variances.");
    574                 psFree(outRO);
    575                 return false;
    576             }
    577         }
    578 
    579         // Need to ensure aperture residual is not calculated
    580         psMetadata *psphotRecipe = psMetadataLookupMetadata(NULL, config->recipes, PSPHOT_RECIPE); // Recipe
    581         psMetadataItem *item = psMetadataLookup(psphotRecipe, "MEASURE.APTREND"); // Item determining aptrend
    582         if (!item) {
    583             psWarning("Unable to find MEASURE.APTREND in psphot recipe");
    584             psErrorClear();
    585         } else {
    586             item->data.B = false;
    587         }
    588 
    589         if (!psphotReadout(config, view)) {
    590             psWarning("Unable to perform photometry on subtracted image.");
    591             psErrorStackPrint(stderr, "Error stack from photometry:");
    592             psErrorClear();
    593         }
    594 
    595         if (!pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL") ||
    596             !pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV") ||
    597             !pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND")) {
    598             psError(PS_ERR_UNKNOWN, false, "Unable to drop PSPHOT internal files.");
    599             return false;
    600         }
    601 
    602         pmFPAfileActivate(config->files, false, "PSPHOT.INPUT");
    603         pmFPAfileActivate(config->files, false, "PSPHOT.LOAD.PSF");
    604 
    605         if (stats) {
    606             pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout with the sources
    607             psArray *sources = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.SOURCES"); // Sources
    608             psMetadataAddS32(stats, PS_LIST_TAIL, "NUM_SOURCES", 0, "Number of sources detected",
    609                              sources ? sources->n : 0);
    610             psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_PHOT", 0, "Time to do photometry",
    611                              psTimerClear("PPSUB_PHOT"));
    612         }
    613 
    614 #ifdef TESTING
    615         // Record data about sources: not everything gets into the output CMF files
    616         {
    617             pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout with the sources
    618             psArray *sources = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.SOURCES"); // Sources
    619             FILE *sourceFile = fopen("sources.dat", "w"); // File for sources
    620             fprintf(sourceFile,
    621                     "# x y mag mag_err psf_chisq cr_nsigma ext_nsigma psf_qf flags m_x m_y m_xx m_xy m_yy\n");
    622             for (int i = 0; i < sources->n; i++) {
    623                 pmSource *source = sources->data[i];
    624                 if (!source) {
    625                     continue;
    626                 }
    627 
    628                 float x, y;             // Position of source
    629                 float chi2;             // chi^2 for source
    630                 if (source->modelPSF) {
    631                     x = source->modelPSF->params->data.F32[PM_PAR_XPOS];
    632                     y = source->modelPSF->params->data.F32[PM_PAR_YPOS];
    633                     chi2 = source->modelPSF->chisq;
    634                 } else if (source->peak) {
    635                     x = source->peak->xf;
    636                     y = source->peak->yf;
    637                     chi2 = NAN;
    638                 } else {
    639                     psWarning("No position available for source.");
    640                     continue;
    641                 }
    642 
    643                 float xMoment = NAN, yMoment = NAN, xxMoment = NAN, xyMoment = NAN, yyMoment = NAN;
    644                 if (source->moments) {
    645                     xMoment = source->moments->Mx;
    646                     yMoment = source->moments->My;
    647                     xxMoment = source->moments->Mxx;
    648                     xyMoment = source->moments->Mxy;
    649                     yyMoment = source->moments->Myy;
    650                 }
    651 
    652                 fprintf(sourceFile, "%f %f %f %f %f %f %f %f %d %f %f %f %f %f\n",
    653                         x, y, source->psfMag, source->errMag, chi2, source->crNsigma, source->extNsigma,
    654                         source->pixWeight, source->mode, xMoment, yMoment, xxMoment, xyMoment, yyMoment);
    655             }
    656             fclose(sourceFile);
    657         }
    658 #endif
    659 
     32    if (!ppSubReadoutSubtract (config, view)) {
     33        psError(PS_ERR_UNKNOWN, false, "Unable to subtract background.");
     34        return false;
    66035    }
    66136
    66237    // Higher order background subtraction using psphot
    663     if (!ppSubBackground(outRO, config, view)) {
     38    if (!ppSubBackground(config, view)) {
    66439        psError(PS_ERR_UNKNOWN, false, "Unable to subtract background.");
    665         psFree(outRO);
     40        return false;
     41    }
     42 
     43    if (!ppSubReadoutPhotometry (config, stats, view)) {
     44        psError(PS_ERR_UNKNOWN, false, "Unable to subtract background.");
    66645        return false;
    66746    }
    66847
    669     // Renormalising for pixels, because that's what magic desires
    670     if (psMetadataLookupBool(&mdok, recipe, "RENORM")) {
    671         psImageMaskType maskValue = pmConfigMaskGet("BLANK", config); // Bits to mask
    672         if (!pmReadoutWeightRenormPixels(outRO, maskValue, renormMean, renormStdev, NULL)) {
    673             psError(PS_ERR_UNKNOWN, false, "Unable to renormalise variances.");
    674             psFree(outRO);
    675             return false;
    676         }
     48    if (!ppSubReadoutUpdate (config, view)) {
     49        psError(PS_ERR_UNKNOWN, false, "Unable to subtract background.");
     50        return false;
    67751    }
    678 
    679     // Add additional data to the header
    680     pmFPAfile *refFile = psMetadataLookupPtr(NULL, config->files, "PPSUB.REF"); // Reference file
    681     pmFPAfile *inFile = psMetadataLookupPtr(NULL, config->files, "PPSUB.INPUT"); // Input file
    682     psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "PPSUB.REFERENCE", 0,
    683                      "Subtraction reference", refFile->filename);
    684     psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "PPSUB.INPUT", 0,
    685                      "Subtraction input", inFile->filename);
    686 
    687 
    688     // Generate binned JPEGs
    689     {
    690         int bin1 = psMetadataLookupS32(NULL, recipe, "BIN1"); // First binning level
    691         int bin2 = psMetadataLookupS32(NULL, recipe, "BIN2"); // Second binning level
    692 
    693         // Target cells
    694         pmCell *cell1 = pmFPAfileThisCell(config->files, view, "PPSUB.OUTPUT.JPEG1");
    695         pmCell *cell2 = pmFPAfileThisCell(config->files, view, "PPSUB.OUTPUT.JPEG2");
    696 
    697         pmReadout *ro1 = pmReadoutAlloc(cell1), *ro2 = pmReadoutAlloc(cell2); // Binned readouts
    698         if (!pmReadoutRebin(ro1, outRO, maskBad, bin1, bin1) || !pmReadoutRebin(ro2, ro1, 0, bin2, bin2)) {
    699             psError(PS_ERR_UNKNOWN, false, "Unable to bin output.");
    700             psFree(ro1);
    701             psFree(ro2);
    702             psFree(outRO);
    703             return false;
    704         }
    705         psFree(ro1);
    706         psFree(ro2);
    707     }
    708 
    709 #ifdef TESTING
    710     // Significance image
    711     {
    712         psImage *sig = (psImage*)psBinaryOp(NULL, outRO->image, "*", outRO->image);
    713         psBinaryOp(sig, sig, "/", outRO->weight);
    714         psFits *fits = psFitsOpen("significance.fits", "w");
    715         psFitsWriteImage(fits, NULL, sig, 0, NULL);
    716         psFitsClose(fits);
    717         psFree(sig);
    718     }
    719 #endif
    720 
    721     psFree(outRO);
    72252
    72353    return true;
Note: See TracChangeset for help on using the changeset viewer.