IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Nov 9, 2009, 2:53:12 PM (17 years ago)
Author:
Paul Price
Message:

Merging branches/pap (unconvolved stacks, reworked combinePixels function into clearer roles, only throw out most variant pixel on each iteration, throw suspect pixels out as first rejection step) stack development. I think I've got everything, but not entirely sure, since I've already merged this branch once before (for dual convolution).

Location:
trunk/ppStack
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppStack

  • trunk/ppStack/src/ppStackReadout.c

    r23577 r26076  
    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
    … …  
    3937
    4038    psArray *args = job->args;          // Arguments
    41     ppStackThread *thread = args->data[0]; // Thread
    42     ppStackOptions *options = args->data[1]; // Options
    43     pmConfig *config = args->data[2];   // Configuration
    44 
    45     pmReadout *outRO = options->outRO;  // Output readout
     39    pmReadout *target = args->data[0];  // Output readout
     40    ppStackThread *thread = args->data[1]; // Thread
     41    ppStackOptions *options = args->data[2]; // Options
     42    pmConfig *config = args->data[3];   // Configuration
     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?
     45
    4646    psVector *mask = options->inputMask; // Mask for inputs
    4747    psArray *rejected = options->rejected; // Rejected pixels
    4848    psVector *weightings = options->weightings; // Weightings (1/noise^2) for each image
    4949    psVector *addVariance = options->matchChi2; // Additional variance when rejecting
    50 
    51     bool status = ppStackReadoutFinal(config, outRO, thread->readouts, mask, rejected,
    52                                       weightings, addVariance); // Status of operation
     50    psVector *norm = normalise ? options->norm : NULL; // Normalisations to apply to images
     51
     52    bool status = ppStackReadoutFinal(config, target, thread->readouts, mask, rejected,
     53                                      weightings, addVariance, safety, norm); // Status of operation
    5354
    5455    thread->busy = false;
    … …  
    6364
    6465    psArray *args = job->args;  // Input arguments
    65     psArray *inspect = args->data[0]; // Array of pixel arrays
    66     int index = PS_SCALAR_VALUE(args->data[1], S32); // Index of interest
    67 
    68     psArray *inputs = inspect->data[index]; // Array of interest
    69     psPixels *output = NULL;    // Output pixel list
    70     for (int i = 0; i < inputs->n; i++) {
    71         psPixels *input = inputs->data[i]; // Input pixel list
    72         if (!input || input->n == 0) {
    73             continue;
    74         }
    75         output = psPixelsConcatenate(output, input);
    76     }
    77 
    78     if (!output) {
    79         // If there are no pixels to inspect, then just fake it
    80         output = psPixelsAllocEmpty(0);
    81     }
    82 
    83     psFree(inputs);
    84     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;
    8597
    8698    return true;
    … …  
    104116
    105117    bool mdok;                          // Status of MD lookup
    106     int iter = psMetadataLookupS32(NULL, recipe, "ITER"); // Rejection iterations
     118    float iter = psMetadataLookupF32(NULL, recipe, "COMBINE.ITER"); // Rejection iterations
    107119    float combineRej = psMetadataLookupF32(NULL, recipe, "COMBINE.REJ"); // Combination threshold
    108120    float combineSys = psMetadataLookupF32(NULL, recipe, "COMBINE.SYS"); // Combination systematic error
    … …  
    114126    int kernelSize = psMetadataLookupS32(NULL, ppsub, "KERNEL.SIZE"); // Kernel half-size
    115127
    116     psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in
     128    psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.IN"); // Name of bits to mask going in
    117129    psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
     130    psString maskSuspectStr = psMetadataLookupStr(NULL, recipe, "MASK.SUSPECT"); // Name of suspect mask bits
     131    psImageMaskType maskSuspect = pmConfigMaskGet(maskSuspectStr, config); // Suspect bits
    118132    psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad
    119133    psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
    … …  
    149163    }
    150164
    151     if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, kernelSize, iter,
    152                         combineRej, combineSys, combineDiscard, true, useVariance, safe, false)) {
     165    if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskSuspect, maskBad, kernelSize, iter,
     166                        combineRej, combineSys, combineDiscard, useVariance, safe, false)) {
    153167        psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts with rejection.");
    154168        psFree(stack);
    … …  
    156170    }
    157171
    158     // Save list of pixels to inspect
     172    // Save lists of pixels
    159173    psArray *inspect = psArrayAlloc(num); // List of pixels to inspect
     174    psArray *reject = psArrayAlloc(num);  // List of pixels rejected
    160175    for (int i = 0; i < num; i++) {
    161176        pmStackData *data = stack->data[i]; // Data for this image
    … …  
    168183        }
    169184        inspect->data[i] = psMemIncrRefCounter(data->inspect);
     185        reject->data[i] = psMemIncrRefCounter(data->reject);
    170186    }
    171187    psFree(stack);
    … …  
    173189    sectionNum++;
    174190
    175     return inspect;
     191    psArray *results = psArrayAlloc(2); // Array of results
     192    results->data[0] = inspect;
     193    results->data[1] = reject;
     194
     195    return results;
    176196}
    177197
    … …  
    180200bool ppStackReadoutFinal(const pmConfig *config, pmReadout *outRO, const psArray *readouts,
    181201                         const psVector *mask, const psArray *rejected, const psVector *weightings,
    182                          const psVector *addVariance)
     202                         const psVector *addVariance, bool safety, const psVector *norm)
    183203{
    184204    assert(config);
    … …  
    196216
    197217    bool mdok;                          // Status of MD lookup
    198     int iter = psMetadataLookupS32(NULL, recipe, "ITER"); // Rejection iterations
    199     float combineRej = psMetadataLookupF32(NULL, recipe, "COMBINE.REJ"); // Combination threshold
    200     float combineSys = psMetadataLookupF32(NULL, recipe, "COMBINE.SYS"); // Combination systematic error
    201     float combineDiscard = psMetadataLookupF32(NULL, recipe, "COMBINE.DISCARD"); // Olympic discard fraction
    202218    bool useVariance = psMetadataLookupBool(&mdok, recipe, "VARIANCE"); // Use variance for rejection?
    203219    bool safe = psMetadataLookupBool(&mdok, recipe, "SAFE"); // Be safe when combining small numbers of pixels
    204220
    205     psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in
     221    psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.IN"); // Name of bits to mask going in
    206222    psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
     223    psString maskSuspectStr = psMetadataLookupStr(NULL, recipe, "MASK.SUSPECT"); // Name of suspect mask bits
     224    psImageMaskType maskSuspect = pmConfigMaskGet(maskSuspectStr, config); // Suspect bits
    207225    psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad
    208226    psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
    … …  
    211229    psArray *stack = psArrayAlloc(num); // Array for stacking
    212230
    213     bool entire = true;                 // Combine entire image?
    214     if (rejected) {
    215         // We have rejection from a previous combination: combine without flagging pixels to inspect
    216         entire = false;
    217         safe = false;
    218         iter = 0;
    219         combineRej = NAN;
    220         combineSys = NAN;
    221     }
     231    // We have rejection from a previous combination: combine without flagging pixels to inspect
     232    safe &= safety;
     233    int iter = 0;
     234    float combineRej = NAN;
     235    float combineSys = NAN;
     236    float combineDiscard = NAN;
    222237
    223238    for (int i = 0; i < num; i++) {
    224239        pmReadout *ro = readouts->data[i];
    225         if (mask->data.U8[i] & (PPSTACK_MASK_REJECT | PPSTACK_MASK_BAD)) {
    226             // Image completely rejected since previous combination
    227             entire = true;
    228             continue;
    229         } else if (mask->data.U8[i]) {
    230             // Image completely rejected before original combination
     240        if (mask->data.U8[i]) {
     241            // Image completely rejected
    231242            continue;
    232243        }
    … …  
    252263        data->reject = rejected ? psMemIncrRefCounter(rejected->data[i]) : NULL;
    253264        stack->data[i] = data;
    254     }
    255 
    256     if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, 0,
    257                         iter, combineRej, combineSys, combineDiscard,
    258                         entire, useVariance, safe, !rejected)) {
     265
     266        if (norm) {
     267            float normalise = powf(10.0, -0.4 * norm->data.F32[i]); // Normalisation
     268            psBinaryOp(ro->image, ro->image, "*", psScalarAlloc(normalise, PS_TYPE_F32));
     269            psBinaryOp(ro->variance, ro->variance, "*", psScalarAlloc(PS_SQR(normalise), PS_TYPE_F32));
     270        }
     271    }
     272
     273    if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskSuspect, maskBad, 0, iter, combineRej,
     274                        combineSys, combineDiscard, useVariance, safe, true)) {
    259275        psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts.");
    260276        psFree(stack);
Note: See TracChangeset for help on using the changeset viewer.