IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
May 3, 2010, 8:50:52 AM (16 years ago)
Author:
eugene
Message:

updates from trunk

Location:
branches/simtest_nebulous_branches
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/simtest_nebulous_branches

  • branches/simtest_nebulous_branches/ppStack/src

    • Property svn:ignore

      •  

        old new  
        1010stamp-h1
        1111ppStackVersionDefinitions.h
         12ppStackErrorCodes.c
         13ppStackErrorCodes.h

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

    r23577 r27840  
    2323    psVector *mask = options->inputMask; // Mask for inputs
    2424    psVector *weightings = options->weightings; // Weightings (1/noise^2) for each image
     25    psVector *exposures = options->exposures;   // Exposure times for each image
    2526    psVector *addVariance = options->matchChi2; // Additional variance when rejecting
    2627
    27     psArray *inspect = ppStackReadoutInitial(config, outRO, thread->readouts, mask,
    28                                              weightings, addVariance);
    29 
    30     job->results = inspect;
     28    job->results = ppStackReadoutInitial(config, outRO, thread->readouts, mask,
     29                                         weightings, exposures, addVariance);
    3130    thread->busy = false;
    3231
    33     return inspect ? true : false;
     32    return job->results ? true : false;
    3433}
    3534
    … …  
    4039    psArray *args = job->args;          // Arguments
    4140    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
     41    psArray *reject = args->data[1];    // Rejected pixels for each image
     42    ppStackOptions *options = args->data[2]; // Options
     43    pmConfig *config = args->data[3];   // Configuration
     44    bool safety = PS_SCALAR_VALUE(args->data[4], U8);    // Safety switch on?
     45    bool normalise = PS_SCALAR_VALUE(args->data[5], U8); // Normalise images?
     46
    4647    psVector *mask = options->inputMask; // Mask for inputs
    47     psArray *rejected = options->rejected; // Rejected pixels
    4848    psVector *weightings = options->weightings; // Weightings (1/noise^2) for each image
     49    psVector *exposures = options->exposures;   // Exposure times for each image
    4950    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
     51    psVector *norm = normalise ? options->norm : NULL; // Normalisations to apply to images
     52
     53    bool status = ppStackReadoutFinal(config, options->outRO, options->expRO, thread->readouts, mask, reject,
     54                                      weightings, exposures, addVariance, safety, norm); // Status of operation
    5355
    5456    thread->busy = false;
    5557
     58    psAssert(status, "Stacking failed.");
     59
    5660    return status;
    5761}
    … …  
    6367
    6468    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;
     69    psArray *inspects = args->data[0]; // Array of pixel arrays
     70    psArray *rejects = args->data[1];  // Array of pixel arrays
     71    int index = PS_SCALAR_VALUE(args->data[2], S32); // Index of interest
     72
     73    psArray *inInspects = inspects->data[index]; // Array of interest
     74    psArray *inRejects = rejects->data[index]; // Array of interest
     75    psAssert(inInspects->n == inRejects->n, "Size should be the same");
     76    psPixels *outInspect = NULL, *outReject = NULL; // Output pixel lists
     77    for (int i = 0; i < inInspects->n; i++) {
     78        psPixels *inInspect = inInspects->data[i]; // Input pixel list
     79        if (inInspect && inInspect->n > 0) {
     80            outInspect = psPixelsConcatenate(outInspect, inInspect);
     81        }
     82        psPixels *inReject = inRejects->data[i]; // Input pixel list
     83        if (inReject && inReject->n > 0) {
     84            outReject = psPixelsConcatenate(outReject, inReject);
     85        }
     86    }
     87
     88    // If there are no pixels to inspect, then just fake it
     89    if (!outInspect) {
     90        outInspect = psPixelsAllocEmpty(0);
     91    }
     92    if (!outReject) {
     93        outReject = psPixelsAllocEmpty(0);
     94    }
     95
     96    psFree(inspects->data[index]);
     97    inspects->data[index] = outInspect;
     98    psFree(rejects->data[index]);
     99    rejects->data[index] = outReject;
    85100
    86101    return true;
    … …  
    89104
    90105psArray *ppStackReadoutInitial(const pmConfig *config, pmReadout *outRO, const psArray *readouts,
    91                                const psVector *mask, const psVector *weightings, const psVector *addVariance)
     106                               const psVector *mask, const psVector *weightings, const psVector *exposures,
     107                               const psVector *addVariance)
    92108{
    93109    assert(config);
    … …  
    104120
    105121    bool mdok;                          // Status of MD lookup
    106     int iter = psMetadataLookupS32(NULL, recipe, "ITER"); // Rejection iterations
     122    float iter = psMetadataLookupF32(NULL, recipe, "COMBINE.ITER"); // Rejection iterations
    107123    float combineRej = psMetadataLookupF32(NULL, recipe, "COMBINE.REJ"); // Combination threshold
    108124    float combineSys = psMetadataLookupF32(NULL, recipe, "COMBINE.SYS"); // Combination systematic error
    … …  
    116132    psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in
    117133    psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
     134    psString maskSuspectStr = psMetadataLookupStr(NULL, recipe, "MASK.SUSPECT"); // Name of suspect mask bits
     135    psImageMaskType maskSuspect = pmConfigMaskGet(maskSuspectStr, config); // Suspect bits
    118136    psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad
    119137    psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
    … …  
    129147        }
    130148
    131 #if 0
    132         // This doesn't seem to work, so getting the weightings directly from a vector
    133         float weighting = psMetadataLookupF32(&mdok, ro->analysis, "PPSTACK.WEIGHTING"); // Relative weight
    134         if (!mdok || !isfinite(weighting)) {
    135             psWarning("No weighting supplied for image %d --- set to unity.", i);
    136             weighting = 1.0;
    137         } else {
    138             psLogMsg("ppStack", PS_LOG_INFO, "Weighting for image %d is %f", i, weighting);
    139         }
    140 #endif
    141 
    142149        // Ensure there is a mask, or pmStackCombine will complain
    143150        if (!ro->mask) {
    … …  
    146153        }
    147154
    148         stack->data[i] = pmStackDataAlloc(ro, weightings->data.F32[i], addVariance->data.F32[i]);
    149     }
    150 
    151     if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, kernelSize, iter,
    152                         combineRej, combineSys, combineDiscard, true, useVariance, safe, false)) {
     155        stack->data[i] = pmStackDataAlloc(ro, weightings->data.F32[i], exposures->data.F32[i],
     156                                          addVariance->data.F32[i]);
     157    }
     158
     159    if (!pmStackCombine(outRO, NULL, stack, maskVal | maskBad, maskSuspect, maskBad, kernelSize, iter,
     160                        combineRej, combineSys, combineDiscard, useVariance, safe, false)) {
    153161        psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts with rejection.");
    154162        psFree(stack);
    … …  
    156164    }
    157165
    158     // Save list of pixels to inspect
     166    // Save lists of pixels
    159167    psArray *inspect = psArrayAlloc(num); // List of pixels to inspect
     168    psArray *reject = psArrayAlloc(num);  // List of pixels rejected
    160169    for (int i = 0; i < num; i++) {
    161170        pmStackData *data = stack->data[i]; // Data for this image
    … …  
    168177        }
    169178        inspect->data[i] = psMemIncrRefCounter(data->inspect);
     179        reject->data[i] = psMemIncrRefCounter(data->reject);
    170180    }
    171181    psFree(stack);
    … …  
    173183    sectionNum++;
    174184
    175     return inspect;
    176 }
    177 
    178 
    179 
    180 bool ppStackReadoutFinal(const pmConfig *config, pmReadout *outRO, const psArray *readouts,
     185    psArray *results = psArrayAlloc(2); // Array of results
     186    results->data[0] = inspect;
     187    results->data[1] = reject;
     188
     189    return results;
     190}
     191
     192
     193
     194bool ppStackReadoutFinal(const pmConfig *config, pmReadout *outRO, pmReadout *expRO, const psArray *readouts,
    181195                         const psVector *mask, const psArray *rejected, const psVector *weightings,
    182                          const psVector *addVariance)
     196                         const psVector *exposures, const psVector *addVariance, bool safety,
     197                         const psVector *norm)
    183198{
    184199    assert(config);
    185200    assert(outRO);
     201    assert(expRO);
    186202    assert(readouts);
    187203    assert(!rejected || readouts->n == rejected->n);
    … …  
    196212
    197213    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
    202214    bool useVariance = psMetadataLookupBool(&mdok, recipe, "VARIANCE"); // Use variance for rejection?
    203215    bool safe = psMetadataLookupBool(&mdok, recipe, "SAFE"); // Be safe when combining small numbers of pixels
    … …  
    205217    psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in
    206218    psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
     219    psString maskSuspectStr = psMetadataLookupStr(NULL, recipe, "MASK.SUSPECT"); // Name of suspect mask bits
     220    psImageMaskType maskSuspect = pmConfigMaskGet(maskSuspectStr, config); // Suspect bits
    207221    psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad
    208222    psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
    … …  
    211225    psArray *stack = psArrayAlloc(num); // Array for stacking
    212226
    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     }
     227    // We have rejection from a previous combination: combine without flagging pixels to inspect
     228    safe &= safety;
     229    int iter = 0;
     230    float combineRej = NAN;
     231    float combineSys = NAN;
     232    float combineDiscard = NAN;
    222233
    223234    for (int i = 0; i < num; i++) {
    224235        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
    231             continue;
    232         }
    233 
    234 #if 0
    235         // This doesn't seem to work, so getting the weightings directly from a vector
    236         bool mdok;                      // Status of MD lookup
    237         float weighting = psMetadataLookupF32(&mdok, ro->analysis, "PPSTACK.WEIGHTING"); // Relative weight
    238         if (!mdok || !isfinite(weighting)) {
    239             psWarning("No WEIGHTING supplied for image %d --- set to unity.", i);
    240             weighting = 1.0;
    241         }
    242 #endif
     236        if (mask->data.U8[i]) {
     237            // Image completely rejected
     238            continue;
     239        }
    243240
    244241        // Ensure there is a mask, or pmStackCombine will complain
    … …  
    248245        }
    249246
    250         pmStackData *data = pmStackDataAlloc(ro, weightings->data.F32[i],
     247        pmStackData *data = pmStackDataAlloc(ro, weightings->data.F32[i], exposures->data.F32[i],
    251248                                             addVariance ? addVariance->data.F32[i] : NAN);
    252249        data->reject = rejected ? psMemIncrRefCounter(rejected->data[i]) : NULL;
    253250        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)) {
     251
     252        if (norm) {
     253            float normalise = powf(10.0, -0.4 * norm->data.F32[i]); // Normalisation
     254            psBinaryOp(ro->image, ro->image, "*", psScalarAlloc(normalise, PS_TYPE_F32));
     255            psBinaryOp(ro->variance, ro->variance, "*", psScalarAlloc(PS_SQR(normalise), PS_TYPE_F32));
     256        }
     257    }
     258
     259    if (!pmStackCombine(outRO, expRO, stack, maskVal | maskBad, maskSuspect, maskBad, 0, iter, combineRej,
     260                        combineSys, combineDiscard, useVariance, safe, rejected)) {
    259261        psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts.");
    260262        psFree(stack);
    … …  
    264266    pmCell *outCell = outRO->parent;    // Output cell
    265267    pmChip *outChip = outCell->parent;  // Output chip
    266 
    267268    outRO->data_exists = true;
    268269    outCell->data_exists = true;
    269270    outChip->data_exists = true;
    270271
     272    pmCell *expCell = expRO->parent;    // Exposure cell
     273    pmChip *expChip = expCell->parent;  // Exposure chip
     274    expRO->data_exists = true;
     275    expCell->data_exists = true;
     276    expChip->data_exists = true;
     277
    271278    psFree(stack);
    272279
Note: See TracChangeset for help on using the changeset viewer.