IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Feb 13, 2009, 1:53:56 PM (17 years ago)
Author:
Paul Price
Message:

Using covariance factor to set the correct variance level for the combination. Adding additional recipe parameter to provide the fraction of values to discard when performing a 'olympic weighted mean'.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppStack/src/ppStackReadout.c

    r21183 r21477  
    1010#include "ppStack.h"
    1111
    12 //#define TESTING                  // Write debugging output?
    13 
    1412bool ppStackReadoutInitialThread(psThreadJob *job)
    1513{
    … …  
    2018    pmConfig *config = args->data[1];   // Configuration
    2119    pmReadout *outRO = args->data[2];   // Output readout
    22     psArray *subRegions = args->data[3]; // Regions for PSF-matching
    23     psArray *subKernels = args->data[4]; // Kernels for PSF-matching
    24     psVector *weightings = args->data[5]; // Weightings (1/noise^2) for each image
    25     psVector *addVariance = args->data[6]; // Additional variance when rejecting
    26 
    27     psArray *inspect = ppStackReadoutInitial(config, outRO, thread->readouts, subRegions, subKernels,
     20    psVector *mask = args->data[3];     // Mask for inputs
     21    psVector *weightings = args->data[4]; // Weightings (1/noise^2) for each image
     22    psVector *addVariance = args->data[5]; // Additional variance when rejecting
     23
     24    psArray *inspect = ppStackReadoutInitial(config, outRO, thread->readouts, mask,
    2825                                             weightings, addVariance);
    2926
    … …  
    4239    pmConfig *config = args->data[1];   // Configuration
    4340    pmReadout *outRO = args->data[2];   // Output readout
    44     psArray *rejected = args->data[3];  // Rejected pixels
    45     psVector *weightings = args->data[4];  // Weighting (1/noise^2) for each image
    46 
    47     bool status = ppStackReadoutFinal(config, outRO, thread->readouts, rejected,
    48                                       weightings); // Status of operation
     41    psVector *mask = args->data[3];     // Mask for inputs
     42    psArray *rejected = args->data[4];  // Rejected pixels
     43    psVector *weightings = args->data[5];  // Weighting (1/noise^2) for each image
     44    psVector *addVariance = args->data[6];  // Weighting (1/noise^2) for each image
     45
     46    bool status = ppStackReadoutFinal(config, outRO, thread->readouts, mask, rejected,
     47                                      weightings, addVariance); // Status of operation
    4948
    5049    thread->busy = false;
    … …  
    8584
    8685psArray *ppStackReadoutInitial(const pmConfig *config, pmReadout *outRO, const psArray *readouts,
    87                                const psArray *regions, const psArray *kernels, const psVector *weightings,
    88                                const psVector *addVariance)
     86                               const psVector *mask, const psVector *weightings, const psVector *addVariance)
    8987{
    9088    assert(config);
    9189    assert(outRO);
    9290    assert(readouts);
    93     assert(regions);
    94     assert(kernels);
    95     assert(readouts->n == regions->n);
    96     assert(regions->n == kernels->n);
     91    assert(mask && mask->n == readouts->n && mask->type.type == PS_TYPE_VECTOR_MASK);
    9792    assert(weightings && weightings->n == readouts->n && weightings->type.type == PS_TYPE_F32);
    9893    assert(addVariance && addVariance->n == readouts->n && addVariance->type.type == PS_TYPE_F32);
    … …  
    107102    float combineRej = psMetadataLookupF32(NULL, recipe, "COMBINE.REJ"); // Combination threshold
    108103    float combineSys = psMetadataLookupF32(NULL, recipe, "COMBINE.SYS"); // Combination systematic error
     104    float combineDiscard = psMetadataLookupF32(NULL, recipe, "COMBINE.DISCARD"); // Olympic discard fraction
    109105    bool useVariance = psMetadataLookupBool(&mdok, recipe, "VARIANCE"); // Use variance for rejection?
    110106    bool safe = psMetadataLookupBool(&mdok, recipe, "SAFE"); // Be safe when combining small numbers of pixels
    … …  
    123119    for (int i = 0; i < num; i++) {
    124120        pmReadout *ro = readouts->data[i];
    125         if (!ro) {
     121        if (!ro || mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
    126122            // Bad image
    127123            continue;
    … …  
    148144    }
    149145
    150     if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, kernelSize, iter, combineRej, combineSys,
    151                         true, useVariance, safe)) {
     146    if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, kernelSize, iter,
     147                        combineRej, combineSys, combineDiscard, true, useVariance, safe)) {
    152148        psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts with rejection.");
    153149        psFree(stack);
    154150        return false;
    155151    }
    156 
    157 #ifdef TESTING
    158     {
    159         psString name = NULL;           // Name of image
    160         psStringAppend(&name, "combined_initial_%03d.fits", sectionNum);
    161         psFits *fits = psFitsOpen(name, "w");
    162         psFree(name);
    163         psFitsWriteImage(fits, NULL, outRO->image, 0, NULL);
    164         psFitsClose(fits);
    165     }
    166 #endif
    167152
    168153    // Save list of pixels to inspect
    … …  
    189174
    190175bool ppStackReadoutFinal(const pmConfig *config, pmReadout *outRO, const psArray *readouts,
    191                          const psArray *rejected, const psVector *weightings)
     176                         const psVector *mask, const psArray *rejected, const psVector *weightings,
     177                         const psVector *addVariance)
    192178{
    193179    assert(config);
    … …  
    196182    assert(rejected);
    197183    assert(readouts->n == rejected->n);
     184    assert(mask && mask->n == readouts->n && mask->type.type == PS_TYPE_VECTOR_MASK);
    198185    assert(weightings && weightings->n == readouts->n && weightings->type.type == PS_TYPE_F32);
    199186
    … …  
    217204    int numGood = num;                  // Number of good inputs: images that haven't been completely rejected
    218205    for (int i = 0; i < num; i++) {
    219         if (!rejected->data[i]) {
     206        pmReadout *ro = readouts->data[i];
     207        if (!ro || !rejected->data[i] || mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
    220208            // Image completely rejected
    221209            numGood--;
    222210            continue;
    223211        }
    224 
    225         pmReadout *ro = readouts->data[i];
    226         assert(ro);
    227212
    228213#if 0
    … …  
    242227        }
    243228
    244         pmStackData *data = pmStackDataAlloc(ro, weightings->data.F32[i], NAN);
     229        pmStackData *data = pmStackDataAlloc(ro, weightings->data.F32[i], addVariance->data.F32[i]);
    245230        data->reject = psMemIncrRefCounter(rejected->data[i]);
    246231        stack->data[i] = data;
    247232    }
    248233
    249     if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, 0, 0, NAN, NAN,
     234    if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, 0, 0, NAN, NAN, NAN,
    250235                        numGood != num, useVariance, false)) {
    251236        psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts.");
    … …  
    254239    }
    255240
    256 #ifdef TESTING
    257     {
    258         psString name = NULL;           // Name of image
    259         psStringAppend(&name, "combined_final_%03d.fits", sectionNum);
    260         psFits *fits = psFitsOpen(name, "w");
    261         psFree(name);
    262         psFitsWriteImage(fits, NULL, outRO->image, 0, NULL);
    263         psFitsClose(fits);
    264     }
    265 #endif
    266 
    267241    pmCell *outCell = outRO->parent;    // Output cell
    268242    pmChip *outChip = outCell->parent;  // Output chip
Note: See TracChangeset for help on using the changeset viewer.