IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 25964


Ignore:
Timestamp:
Oct 28, 2009, 5:35:05 PM (17 years ago)
Author:
Paul Price
Message:

Reworking stack combination because there are *three* modes for pixels going into the final stack (after rejection), not just two: tested and good, tested and rejected, and not tested. The code did not recognise the third, which is a distinct state because we don't want these pixels grown, as we do for rejected pixels. This cannot be fixed merely by using the 'safe' combination because that would discard the 'tested and good' pixels that have only a single unrejected input but are good because they have survived the testing process. Needed to add a new state into the combination process. Now I add these pixels straight into the 'reject' list. This requires a little bit more fiddling around in ppStack. Not sure it's working yet.

Location:
branches/pap
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • branches/pap/ppStack/src/ppStack.h

    r25924 r25964  
    5959// Perform stacking on a readout
    6060//
    61 // Returns an array of pixels to inspect for each input image
     61// Returns two arrays: pixels to inspect for each input image, and pixels to reject for each input image.
    6262psArray *ppStackReadoutInitial(const pmConfig *config,   // Configuration
    6363                               pmReadout *outRO,   // Output readout
    … …  
    8484                         const psVector *weightings, // Weighting factors for each image
    8585                         const psVector *addVariance, // Additional variance for rejection
    86                          bool full,                   // Combine full image?
    8786                         bool safety,                 // Enable safety switch?
    8887                         const psVector *norm         // Normalisations to apply

  • branches/pap/ppStack/src/ppStackCombineFinal.c

    r25950 r25964  
    1313
    1414bool ppStackCombineFinal(pmReadout *target, ppStackThreadData *stack, psArray *covariances,
    15                          ppStackOptions *options, pmConfig *config, bool full, bool safe, bool normalise)
     15                         ppStackOptions *options, pmConfig *config, bool safe, bool normalise)
    1616{
    1717    psAssert(stack, "Require stack");
    … …  
    4343        psArrayAdd(job->args, 1, options);
    4444        psArrayAdd(job->args, 1, config);
    45         PS_ARRAY_ADD_SCALAR(job->args, full, PS_TYPE_U8);
    4645        PS_ARRAY_ADD_SCALAR(job->args, safe, PS_TYPE_U8);
    4746        PS_ARRAY_ADD_SCALAR(job->args, normalise, PS_TYPE_U8);

  • branches/pap/ppStack/src/ppStackCombineInitial.c

    r25950 r25964  
    6262    // Harvest the jobs, gathering the inspection lists
    6363    options->inspect = psArrayAlloc(options->num);
     64    options->rejected = psArrayAlloc(options->num);
    6465    for (int i = 0; i < options->num; i++) {
    6566        if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
    … …  
    6768        }
    6869        options->inspect->data[i] = psArrayAllocEmpty(numChunk);
     70        options->rejected->data[i] = psArrayAllocEmpty(numChunk);
    6971    }
    7072    psThreadJob *job;               // Completed job
    … …  
    7375                 "Job has incorrect type: %s", job->type);
    7476        psArray *results = job->results; // Results of job
     77        psAssert(results->n == 2, "Results array has wrong size!");
     78        psArray *inspect = results->data[0]; // Pixels to inspect
     79        psArray *reject = results->data[1];  // Pixels to reject
    7580        for (int i = 0; i < options->num; i++) {
    7681            if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
    7782                continue;
    7883            }
    79             options->inspect->data[i] = psArrayAdd(options->inspect->data[i], 1, results->data[i]);
     84            options->inspect->data[i] = psArrayAdd(options->inspect->data[i], 1, inspect->data[i]);
     85            options->rejected->data[i] = psArrayAdd(options->rejected->data[i], 1, reject->data[i]);
    8086        }
    8187        psFree(job);

  • branches/pap/ppStack/src/ppStackLoop.c

    r25950 r25964  
    9797    // Final combination
    9898    psTrace("ppStack", 2, "Final stack of convolved images....\n");
    99     if (!ppStackCombineFinal(options->outRO, stack, options->convCovars, options, config,
    100                              true, false, false)) {
     99    if (!ppStackCombineFinal(options->outRO, stack, options->convCovars, options, config, false, false)) {
    101100        psError(PS_ERR_UNKNOWN, false, "Unable to perform final combination.");
    102101        psFree(stack);
    … …  
    106105    psLogMsg("ppStack", PS_LOG_INFO, "Stage 5: Final Stack: %f sec", psTimerClear("PPSTACK_STEPS"));
    107106    ppStackMemDump("final");
    108 
    109107
    110108    // Clean up
    … …  
    121119    psFree(stack);
    122120
    123 
    124121#if 1
    125122    // Unconvolved stack --- it's cheap to calculate, compared to everything else!
    … …  
    133130        }
    134131        psTrace("ppStack", 2, "Stack of unconvolved images....\n");
    135         if (!ppStackCombineFinal(options->unconvRO, stack, options->origCovars, options, config,
    136                                  true, true, true)) {
     132        if (!ppStackCombineFinal(options->unconvRO, stack, options->origCovars, options, config, true, true)) {
    137133            psError(PS_ERR_UNKNOWN, false, "Unable to perform unconvolved combination.");
    138134            psFree(stack);
    … …  
    158154    ppStackMemDump("photometry");
    159155
    160 
    161156    // Finish up
    162157    psTrace("ppStack", 1, "Finishing up....\n");
    … …  
    170165
    171166    psFree(options);
     167
    172168    return true;
    173169}

  • branches/pap/ppStack/src/ppStackLoop.h

    r25924 r25964  
    6161    ppStackOptions *options,            // Options
    6262    pmConfig *config,                   // Configuration
    63     bool full,                          // Combine full image?
    6463    bool safe,                          // Allow safe combination?
    6564    bool norm                           // Normalise images?

  • branches/pap/ppStack/src/ppStackReadout.c

    r25924 r25964  
    2525    psVector *addVariance = options->matchChi2; // Additional variance when rejecting
    2626
    27     psArray *inspect = ppStackReadoutInitial(config, outRO, thread->readouts, mask,
    28                                              weightings, addVariance);
    29 
    30     job->results = inspect;
     27    job->results = ppStackReadoutInitial(config, outRO, thread->readouts, mask,
     28                                         weightings, addVariance);
    3129    thread->busy = false;
    3230
    33     return inspect ? true : false;
     31    return job->results ? true : false;
    3432}
    3533
    … …  
    4341    ppStackOptions *options = args->data[2]; // Options
    4442    pmConfig *config = args->data[3];   // Configuration
    45     bool full = PS_SCALAR_VALUE(args->data[4], U8); // Combine full image?
    46     bool safety = PS_SCALAR_VALUE(args->data[5], U8);    // Safety switch on?
    47     bool normalise = PS_SCALAR_VALUE(args->data[6], U8); // Normalise images?
     43    bool safety = PS_SCALAR_VALUE(args->data[4], U8);    // Safety switch on?
     44    bool normalise = PS_SCALAR_VALUE(args->data[5], U8); // Normalise images?
    4845
    4946    psVector *mask = options->inputMask; // Mask for inputs
    … …  
    5451
    5552    bool status = ppStackReadoutFinal(config, target, thread->readouts, mask, rejected,
    56                                       weightings, addVariance, full, safety, norm); // Status of operation
     53                                      weightings, addVariance, safety, norm); // Status of operation
    5754
    5855    thread->busy = false;
    … …  
    6764
    6865    psArray *args = job->args;  // Input arguments
    69     psArray *inspect = args->data[0]; // Array of pixel arrays
    70     int index = PS_SCALAR_VALUE(args->data[1], S32); // Index of interest
    71 
    72     psArray *inputs = inspect->data[index]; // Array of interest
    73     psPixels *output = NULL;    // Output pixel list
    74     for (int i = 0; i < inputs->n; i++) {
    75         psPixels *input = inputs->data[i]; // Input pixel list
    76         if (!input || input->n == 0) {
    77             continue;
    78         }
    79         output = psPixelsConcatenate(output, input);
    80     }
    81 
    82     if (!output) {
    83         // If there are no pixels to inspect, then just fake it
    84         output = psPixelsAllocEmpty(0);
    85     }
    86 
    87     psFree(inputs);
    88     inspect->data[index] = output;
     66    psArray *inspects = args->data[0]; // Array of pixel arrays
     67    psArray *rejects = args->data[1];  // Array of pixel arrays
     68    int index = PS_SCALAR_VALUE(args->data[2], S32); // Index of interest
     69
     70    psArray *inInspects = inspects->data[index]; // Array of interest
     71    psArray *inRejects = rejects->data[index]; // Array of interest
     72    psAssert(inInspects->n == inRejects->n, "Size should be the same");
     73    psPixels *outInspect = NULL, *outReject = NULL; // Output pixel lists
     74    for (int i = 0; i < inInspects->n; i++) {
     75        psPixels *inInspect = inInspects->data[i]; // Input pixel list
     76        if (inInspect && inInspect->n > 0) {
     77            outInspect = psPixelsConcatenate(outInspect, inInspect);
     78        }
     79        psPixels *inReject = inRejects->data[i]; // Input pixel list
     80        if (inReject && inReject->n > 0) {
     81            outReject = psPixelsConcatenate(outReject, inReject);
     82        }
     83    }
     84
     85    // If there are no pixels to inspect, then just fake it
     86    if (!outInspect) {
     87        outInspect = psPixelsAllocEmpty(0);
     88    }
     89    if (!outReject) {
     90        outReject = psPixelsAllocEmpty(0);
     91    }
     92
     93    psFree(inspects->data[index]);
     94    inspects->data[index] = outInspect;
     95    psFree(rejects->data[index]);
     96    rejects->data[index] = outReject;
    8997
    9098    return true;
    … …  
    154162
    155163    if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, kernelSize, iter,
    156                         combineRej, combineSys, combineDiscard, true, useVariance, safe, false)) {
     164                        combineRej, combineSys, combineDiscard, useVariance, safe, false)) {
    157165        psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts with rejection.");
    158166        psFree(stack);
    … …  
    160168    }
    161169
    162     // Save list of pixels to inspect
     170    // Save lists of pixels
    163171    psArray *inspect = psArrayAlloc(num); // List of pixels to inspect
     172    psArray *reject = psArrayAlloc(num);  // List of pixels rejected
    164173    for (int i = 0; i < num; i++) {
    165174        pmStackData *data = stack->data[i]; // Data for this image
    … …  
    172181        }
    173182        inspect->data[i] = psMemIncrRefCounter(data->inspect);
     183        reject->data[i] = psMemIncrRefCounter(data->reject);
    174184    }
    175185    psFree(stack);
    … …  
    177187    sectionNum++;
    178188
    179     return inspect;
     189    psArray *results = psArrayAlloc(2); // Array of results
     190    results->data[0] = inspect;
     191    results->data[1] = reject;
     192
     193    return results;
    180194}
    181195
    … …  
    184198bool ppStackReadoutFinal(const pmConfig *config, pmReadout *outRO, const psArray *readouts,
    185199                         const psVector *mask, const psArray *rejected, const psVector *weightings,
    186                          const psVector *addVariance, bool full, bool safety, const psVector *norm)
     200                         const psVector *addVariance, bool safety, const psVector *norm)
    187201{
    188202    assert(config);
    … …  
    225239    for (int i = 0; i < num; i++) {
    226240        pmReadout *ro = readouts->data[i];
    227         if (mask->data.U8[i] & (PPSTACK_MASK_REJECT | PPSTACK_MASK_BAD)) {
    228             // Image completely rejected since previous combination
    229             full = true;
    230             continue;
    231         } else if (mask->data.U8[i]) {
    232             // Image completely rejected before original combination
     241        if (mask->data.U8[i]) {
     242            // Image completely rejected
    233243            continue;
    234244        }
    … …  
    262272    }
    263273
    264     if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, 0,
    265                         iter, combineRej, combineSys, combineDiscard,
    266                         full, useVariance, safe, !rejected)) {
     274    if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, 0, iter, combineRej,
     275                        combineSys, combineDiscard, useVariance, safe, !rejected)) {
    267276        psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts.");
    268277        psFree(stack);

  • branches/pap/ppStack/src/ppStackReject.c

    r25950 r25964  
    2323
    2424    int num = options->num;             // Number of inputs
    25     options->rejected = psArrayAlloc(num);
    2625
    2726    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
    … …  
    5352        psThreadJob *job = psThreadJobAlloc("PPSTACK_INSPECT"); // Job to start
    5453        psArrayAdd(job->args, 1, options->inspect);
     54        psArrayAdd(job->args, 1, options->rejected);
    5555        PS_ARRAY_ADD_SCALAR(job->args, i, PS_TYPE_S32);
    5656        if (!psThreadJobAddPending(job)) {
    … …  
    9292#endif
    9393
    94         psPixels *reject = pmStackReject(options->inspect->data[i], options->numCols, options->numRows,
    95                                          threshold, poorFrac, stride, options->regions->data[i],
    96                                          options->kernels->data[i]); // Rejected pixels
    97 
    9894#ifdef TESTING
    9995        {
    100             psImage *mask = psPixelsToMask(NULL, reject,
     96            psImage *mask = psPixelsToMask(NULL, options->rejected->data[i],
    10197                                           psRegionSet(0, options->numCols - 1, 0, options->numRows - 1),
    10298                                           0xff); // Mask image
    10399            psString name = NULL;           // Name of image
    104             psStringAppend(&name, "reject_%03d.fits", i);
     100            psStringAppend(&name, "pre_reject_%03d.fits", i);
    105101            pmStackVisualPlotTestImage(mask, name);
    106102            psFits *fits = psFitsOpen(name, "w");
    … …  
    111107        }
    112108#endif
     109
     110        psPixels *reject = pmStackReject(options->inspect->data[i], options->numCols, options->numRows,
     111                                         threshold, poorFrac, stride, options->regions->data[i],
     112                                         options->kernels->data[i]); // Rejected pixels
    113113
    114114        psFree(options->inspect->data[i]);
    … …  
    127127                          "exceeds limit (%.3f)", i, frac, imageRej);
    128128                psFree(reject);
    129                 // reject == NULL means reject image completely
    130129                reject = NULL;
    131130                options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PPSTACK_MASK_BAD;
    … …  
    134133        }
    135134
    136         // Images without a list of rejected pixels (the list may be empty) are rejected completely
    137         options->rejected->data[i] = reject;
     135        if (reject) {
     136            // Add to list of pixels already rejected
     137            reject = psPixelsConcatenate(reject, options->rejected->data[i]);
     138            options->rejected->data[i] = psPixelsDuplicates(options->rejected->data[i], reject);
     139        }
     140
     141#ifdef TESTING
     142        {
     143            psImage *mask = psPixelsToMask(NULL, options->rejected->data[i],
     144                                           psRegionSet(0, options->numCols - 1, 0, options->numRows - 1),
     145                                           0xff); // Mask image
     146            psString name = NULL;           // Name of image
     147            psStringAppend(&name, "reject_%03d.fits", i);
     148            pmStackVisualPlotTestImage(mask, name);
     149            psFits *fits = psFitsOpen(name, "w");
     150            psFree(name);
     151            psFitsWriteImage(fits, NULL, mask, 0, NULL);
     152            psFree(mask);
     153            psFitsClose(fits);
     154        }
     155#endif
    138156
    139157        if (options->stats) {
    … …  
    143161                             "Number of pixels rejected", reject ? reject->n : 0);
    144162        }
     163
     164        psFree(reject);
    145165        psLogMsg("ppStack", PS_LOG_INFO, "Time to perform rejection on image %d: %f sec", i,
    146166                 psTimerClear("PPSTACK_REJECT"));

  • branches/pap/ppStack/src/ppStackThread.c

    r25924 r25964  
    275275
    276276    {
    277         psThreadTask *task = psThreadTaskAlloc("PPSTACK_INSPECT", 2);
     277        psThreadTask *task = psThreadTaskAlloc("PPSTACK_INSPECT", 3);
    278278        task->function = &ppStackInspect;
    279279        psThreadTaskAdd(task);
    … …  
    282282
    283283    {
    284         psThreadTask *task = psThreadTaskAlloc("PPSTACK_FINAL_COMBINE", 7);
     284        psThreadTask *task = psThreadTaskAlloc("PPSTACK_FINAL_COMBINE", 6);
    285285        task->function = &ppStackReadoutFinalThread;
    286286        psThreadTaskAdd(task);

  • branches/pap/psModules/src/imcombine/pmStack.c

    r25960 r25964  
    290290
    291291// Mark a pixel for inspection
     292// Value in pixel doesn't seem to agree with the stack, so need to look closer
    292293static inline void combineInspect(const psArray *inputs, // Stack data
    293294                                  int x, int y, // Pixel
    … …  
    301302    }
    302303    data->inspect = psPixelsAdd(data->inspect, data->inspect->nalloc, x, y);
     304    return;
     305}
     306
     307// Mark a pixel for rejection
     308// Cannot possibly inspect this pixel and confirm that it's good.
     309// e.g., Only a single input
     310static inline void combineReject(const psArray *inputs, // Stack data
     311                                 int x, int y, // Pixel
     312                                 int source // Source image index
     313                                 )
     314{
     315    pmStackData *data = inputs->data[source]; // Stack data of interest
     316    if (!data) {
     317        psWarning("Can't find input data for source %d", source);
     318        return;
     319    }
     320    data->reject = psPixelsAdd(data->reject, data->reject->nalloc, x, y);
    303321    return;
    304322}
    … …  
    319337                          float rej, // Number of standard deviations at which to reject
    320338                          float sys,    // Relative systematic error
    321                           float discard,// Fraction of values to discard (Olympic weighted mean)
     339                          float olympic,// Fraction of values to discard (Olympic weighted mean)
    322340                          bool useVariance, // Use variance for rejection when combining?
    323341                          bool safe,    // Combine safely?
    324                           bool rejectInspect, // Reject values marked for inspection from combination?
     342                          bool rejection, // Reject values marked for inspection from combination?
    325343                          combineBuffer *buffer // Buffer for combination; to avoid multiple allocations
    326344                         )
    … …  
    419437              }
    420438              maskValue = 0;
    421           }
    422 #ifdef TESTING
    423           else {
     439          } else {
     440              if (!rejection) {
     441                  combineReject(inputs, x, y, pixelSources->data.U16[0]);
     442              }
     443#ifdef TESTING
    424444              numRejected = 1;
    425445              if (x == TEST_X && y == TEST_Y) {
    426446                  fprintf(stderr, "Single input to combine, safety on, pixel is bad.\n");
    427447              }
     448#endif
    428449          }
    429 #endif
    430450          break;
    431451      }
    … …  
    459479              if (PS_SQR(diff) > PS_SQR(rej) * sigma2) {
    460480                  // Not consistent: mark both for inspection
    461                   if (rejectInspect) {
     481                  if (rejection) {
    462482                      imageValue = NAN;
    463483                      varianceValue = NAN;
    … …  
    550570                  }
    551571#endif
    552                   median = combinationWeightedOlympic(pixelData, pixelWeights, pixelMasks, discard, sort);
     572                  median = combinationWeightedOlympic(pixelData, pixelWeights, pixelMasks, olympic, sort);
    553573              }
    554574
    … …  
    563583#define MASK_PIXEL_FOR_INSPECTION() \
    564584    pixelMasks->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xff; \
    565     if (!rejectInspect) { \
     585    if (!rejection) { \
    566586        combineInspect(inputs, x, y, pixelSources->data.U16[j]); \
    567587    } \
    … …  
    605625          }
    606626
    607           if (rejectInspect && totalClipped > 0) {
     627          if (rejection && totalClipped > 0) {
    608628              // Get rid of the masked values
    609629              // The alternative to this is to make combinationMeanVariance() accept a mask
    … …  
    658678// Ensure the input array of pmStackData is valid, and get some details out of it
    659679static bool validateInputData(bool *haveVariances, // Do we have variance maps in the sky images?
    660                               bool *haveRejects, // Do we have lists of rejected pixels?
    661680                              int *num,    // Number of inputs
    662681                              int *numCols, int *numRows, // Size of (sky) images
    663                               psArray *input // Input array of pmStackData to validate
     682                              const psArray *input, // Input array of pmStackData to validate
     683                              const pmReadout *output // Output readout
    664684    )
    665685{
    … …  
    687707        PS_ASSERT_IMAGE_TYPE(data->readout->variance, PS_TYPE_F32, false);
    688708    }
    689     *haveRejects = (data->reject != NULL);
     709    bool haveRejects = (data->reject != NULL); // Do we have rejected pixels?
    690710
    691711    // Make sure the rest correspond with the first
    … …  
    700720            return false;
    701721        }
    702         if ((*haveRejects && !data->reject) || (data->reject && !*haveRejects)) {
     722        if ((haveRejects && !data->reject) || (data->reject && !haveRejects)) {
    703723            psError(PS_ERR_UNEXPECTED_NULL, true,
    704724                    "The rejected pixels are specified in some but not all inputs.");
    … …  
    716736            PS_ASSERT_IMAGE_TYPE(data->readout->variance, PS_TYPE_F32, false);
    717737        }
     738    }
     739
     740    PM_ASSERT_READOUT_NON_NULL(output, false);
     741    if (output->image) {
     742        PS_ASSERT_IMAGE_NON_NULL(output->image, false);
     743        PS_ASSERT_IMAGE_TYPE(output->image, PS_TYPE_F32, false);
     744        PS_ASSERT_IMAGE_NON_NULL(output->mask, false);
     745        PS_ASSERT_IMAGE_TYPE(output->mask, PS_TYPE_IMAGE_MASK, false);
     746        PS_ASSERT_IMAGES_SIZE_EQUAL(output->image, output->mask, false);
    718747    }
    719748
    … …  
    802831/// Stack input images
    803832bool pmStackCombine(pmReadout *combined, psArray *input, psImageMaskType maskVal, psImageMaskType bad,
    804                     int kernelSize, int numIter, float rej, float sys, float discard,
    805                     bool entire, bool useVariance, bool safe, bool rejectInspect)
    806 {
    807     PS_ASSERT_PTR_NON_NULL(combined, false);
     833                    int kernelSize, int numIter, float rej, float sys, float olympic,
     834                    bool useVariance, bool safe, bool rejection)
     835{
    808836    bool haveVariances;                 // Do we have the variance maps?
    809     bool haveRejects;                   // Do we have lists of rejected pixels?
    810837    int num;                            // Number of inputs
    811838    int numCols, numRows;               // Size of (sky) images
    812     if (!validateInputData(&haveVariances, &haveRejects, &num, &numCols, &numRows, input)) {
     839    if (!validateInputData(&haveVariances, &num, &numCols, &numRows, input, combined)) {
    813840        return false;
    814841    }
    … …  
    821848        PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false);
    822849    }
    823     if (haveRejects) {
    824         // This is a subsequent combination, so expect that the image and mask already exist
    825         PS_ASSERT_IMAGE_NON_NULL(combined->image, false);
    826         PS_ASSERT_IMAGE_TYPE(combined->image, PS_TYPE_F32, false);
    827         PS_ASSERT_IMAGE_NON_NULL(combined->mask, false);
    828         PS_ASSERT_IMAGE_TYPE(combined->mask, PS_TYPE_IMAGE_MASK, false);
    829         PS_ASSERT_IMAGES_SIZE_EQUAL(combined->image, combined->mask, false);
    830     }
    831850    if (useVariance && !haveVariances) {
    832851        psWarning("Unable to use variance in rejection if no variance maps supplied --- option turned off");
    … …  
    852871        }
    853872#endif
    854         if (!haveRejects && !data->inspect) {
    855             data->inspect = psPixelsAllocEmpty(PIXEL_LIST_BUFFER);
     873        if (!rejection) {
     874            // Ensure pixels can be put on the appropriate list
     875            if (!data->inspect) {
     876                data->inspect = psPixelsAllocEmpty(PIXEL_LIST_BUFFER);
     877            }
     878            if (!data->reject) {
     879                data->reject = psPixelsAllocEmpty(PIXEL_LIST_BUFFER);
     880            }
    856881        }
    857882    }
    … …  
    882907    combineBuffer *buffer = combineBufferAlloc(num);
    883908
    884     if (haveRejects) {
    885         psImage *combinedImage = combined->image; // Combined image
    886         psImage *combinedMask = combined->mask; // Combined mask
    887         psImage *combinedVariance = combined->variance; // Combined variance map
    888 
    889         psArray *pixelMap = pixelMapGenerate(input, minInputCols, maxInputCols,
    890                                              minInputRows, maxInputRows); // Map of pixels to source
    891         psPixels *pixels = NULL;            // Total list of pixels, with no duplicates
     909    // Pull the products out, allocate if necessary
     910    psImage *combinedImage = combined->image; // Combined image
     911    if (!combinedImage) {
     912        combined->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     913        combinedImage = combined->image;
     914    }
     915    psImage *combinedMask = combined->mask; // Combined mask
     916    if (!combinedMask) {
     917        combined->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
     918        combinedMask = combined->mask;
     919    }
     920
     921    psImage *combinedVariance = combined->variance; // Combined variance map
     922    if (haveVariances && !combinedVariance) {
     923        combined->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     924        combinedVariance = combined->variance;
     925    }
     926
     927    // Set up rejection list
     928    psArray *pixelMap = NULL;           // Map of pixels to source
     929    if (rejection) {
     930        pixelMap = pixelMapGenerate(input, minInputCols, maxInputCols, minInputRows, maxInputRows);
     931    }
     932
     933    // Combine each pixel
     934    for (int y = minInputRows; y < maxInputRows; y++) {
     935        for (int x = minInputCols; x < maxInputCols; x++) {
     936            psVector *reject = NULL; // Images to reject for this pixel
     937            if (rejection) {
     938                reject = pixelMapQuery(pixelMap, minInputCols, minInputRows, x, y);
     939            }
     940            combinePixels(combinedImage, combinedMask, combinedVariance, input, weights,
     941                          addVariance, reject, x, y, maskVal, bad, numIter, rej, sys, olympic,
     942                          useVariance, safe, rejection, buffer);
     943        }
     944    }
     945
     946    psFree(pixelMap);
     947    psFree(weights);
     948    psFree(buffer);
     949
     950#ifndef PS_NO_TRACE
     951    if (!rejection && psTraceGetLevel("psModules.imcombine") >= 5) {
    892952        for (int i = 0; i < num; i++) {
    893             pmStackData *data = input->data[i]; // Stacking data; contains the list of pixels
    894             if (!data) {
     953            pmStackData *data = input->data[i]; // Stack data for this input
     954            if (!data || !data->inspect) {
    895955                continue;
    896956            }
    897             pixels = psPixelsConcatenate(pixels, data->reject);
    898         }
    899         pixels = psPixelsDuplicates(pixels, pixels);
    900 
    901         if (entire) {
    902             // Combine entire image
    903             for (int y = minInputRows; y < maxInputRows; y++) {
    904                 for (int x = minInputCols; x < maxInputCols; x++) {
    905                     psVector *reject = pixelMapQuery(pixelMap, minInputCols, minInputRows,
    906                                                      x, y); // Reject these images
    907                     combinePixels(combinedImage, combinedMask, combinedVariance, input, weights,
    908                                   addVariance, reject, x, y, maskVal, bad, numIter, rej, sys, discard,
    909                                   useVariance, safe, rejectInspect, buffer);
    910                 }
    911             }
    912         } else {
    913             // Only combine previously rejected pixels
    914             for (int i = 0; i < pixels->n; i++) {
    915                 // Pixel coordinates are in the frame of the original image
    916                 int x = pixels->data[i].x, y = pixels->data[i].y; // Coordinates of interest
    917                 if (x < minInputCols || x >= maxInputCols || y < minInputRows || y >= maxInputRows) {
    918                     continue;
    919                 }
    920                 psVector *reject = pixelMapQuery(pixelMap, minInputCols, minInputRows,
    921                                                  x, y); // Reject these images
    922                 combinePixels(combinedImage, combinedMask, combinedVariance, input, weights,
    923                               addVariance, reject, x, y, maskVal, bad, numIter, rej, sys, discard,
    924                               useVariance, safe, rejectInspect, buffer);
    925             }
    926         }
    927         psFree(pixels);
    928         psFree(pixelMap);
    929     } else {
    930         // Pull the products out, allocate if necessary
    931         psImage *combinedImage = combined->image; // Combined image
    932         if (!combinedImage) {
    933             combined->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    934             combinedImage = combined->image;
    935         }
    936         psImage *combinedMask = combined->mask; // Combined mask
    937         if (!combinedMask) {
    938             combined->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
    939             combinedMask = combined->mask;
    940         }
    941 
    942         psImage *combinedVariance = combined->variance; // Combined variance map
    943         if (haveVariances && !combinedVariance) {
    944             combined->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    945             combinedVariance = combined->variance;
    946         }
    947 
    948         for (int y = minInputRows; y < maxInputRows; y++) {
    949             for (int x = minInputCols; x < maxInputCols; x++) {
    950                 combinePixels(combinedImage, combinedMask, combinedVariance, input, weights,
    951                               addVariance, NULL, x, y, maskVal, bad, numIter, rej, sys, discard,
    952                               useVariance, safe, rejectInspect, buffer);
    953             }
    954         }
    955 
    956 #ifndef PS_NO_TRACE
    957         if (psTraceGetLevel("psModules.imcombine") >= 5) {
    958             for (int i = 0; i < num; i++) {
    959                 pmStackData *data = input->data[i]; // Stack data for this input
    960                 if (!data || !data->inspect) {
    961                     continue;
    962                 }
    963                 psTrace("psModules.imcombine", 5, "Image %d: %ld pixels to inspect.\n", i, data->inspect->n);
    964             }
    965         }
    966 #endif
    967     }
    968 
    969     psFree(weights);
    970     psFree(buffer);
     957            psTrace("psModules.imcombine", 5, "Image %d: %ld pixels to inspect.\n", i, data->inspect->n);
     958        }
     959    }
     960#endif
    971961
    972962    return true;

  • branches/pap/psModules/src/imcombine/pmStack.h

    r23577 r25964  
    3030    psPixels *reject;                   ///< Pixels to reject
    3131    psPixels *inspect;                  ///< Pixels to inspect
     32    psPixels *discard;                  ///< Pixels to discard
    3233    float weight;                       ///< Relative weighting for image
    3334    float addVariance;                  ///< Additional variance when rejecting
    … …  
    4950                    float rej,          ///< Rejection limit (standard deviations)
    5051                    float sys,          ///< Relative systematic error
    51                     float discard,      ///< Fraction of values to discard for Olympic weighted mean
    52                     bool entire,        ///< Combine entire image even if rejection lists provided?
     52                    float olympic,      ///< Fraction of values to discard for Olympic weighted mean
    5353                    bool useVariance,   ///< Use variance values for rejection?
    5454                    bool safe,          ///< Play safe with small numbers of input pixels (mask if N <= 2)?
    55                     bool rejectInspect  ///< Reject pixels instead of marking them for inspection?
     55                    bool rejection      ///< Reject pixels instead of marking them for inspection?
    5656    );
    5757
Note: See TracChangeset for help on using the changeset viewer.