IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 25964 for branches/pap/ppStack


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/ppStack/src
Files:
8 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);
Note: See TracChangeset for help on using the changeset viewer.