IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 26007


Ignore:
Timestamp:
Nov 2, 2009, 5:08:31 PM (17 years ago)
Author:
Paul Price
Message:

Serious rework of pmStackCombine for clarity and some additional functionality. Split the combinePixels function into two operating modes (test and combine) to make clear what we're doing; pmStackCombine could be split into two separate functions (e.g., pmStackTest and pmStackCombine) in the future, since the testing part doesn't need to produce any output image, but this involves moving a lot of stuff around in ppStack, so not doing it yet. Verified that variance scaling is correct, and softening of variance is correct. Big improvement in minimising loss of good pixels comes from throwing out only the most variant pixel on each rejection pass (consequently changed iterations parameter to iterations per input); this involved another rework of the core combinePixels function (moving the switch on the number of valid inputs inside the iteration). Much happier with the outlier rejection now.

Location:
branches/pap
Files:
12 edited

Legend:

Unmodified
Added
Removed

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

    r23841 r26007  
    148148    psMetadataAddStr(arguments, PS_LIST_TAIL, "-stamps", 0, "Stamps file with x,y,flux per line", NULL);
    149149    psMetadataAddStr(arguments, PS_LIST_TAIL, "-stats", 0, "Statistics file", NULL);
    150     psMetadataAddS32(arguments, PS_LIST_TAIL, "-iter", 0, "Number of rejection iterations", 0);
     150    psMetadataAddF32(arguments, PS_LIST_TAIL, "-combine-iter", 0, "Number of rejection iterations per input", NAN);
    151151    psMetadataAddF32(arguments, PS_LIST_TAIL, "-combine-rej", 0,
    152152                     "Combination rejection thresold (sigma)", NAN);
    … …  
    250250    }
    251251
    252     VALUE_ARG_RECIPE_INT("-iter",              "ITER",            S32, 0);
     252    VALUE_ARG_RECIPE_FLOAT("-combine-iter",    "COMBINE.ITER",    F32);
    253253    VALUE_ARG_RECIPE_FLOAT("-combine-rej",     "COMBINE.REJ",     F32);
    254254    VALUE_ARG_RECIPE_FLOAT("-combine-sys",     "COMBINE.SYS",     F32);
    … …  
    260260    VALUE_ARG_RECIPE_FLOAT("-poor-frac",       "POOR.FRACTION",   F32);
    261261
    262     valueArgRecipeStr(arguments, recipe, "-mask-val",  "MASK.VAL",  recipe);
     262    valueArgRecipeStr(arguments, recipe, "-mask-val",  "MASK.IN",   recipe);
    263263    valueArgRecipeStr(arguments, recipe, "-mask-bad",  "MASK.BAD",  recipe);
    264264    valueArgRecipeStr(arguments, recipe, "-mask-poor", "MASK.POOR", recipe);

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

    r25964 r26007  
    1010#include "ppStackLoop.h"
    1111
    12 #define TESTING                         // Enable test output
     12//#define TESTING                         // Enable test output
    1313
    1414bool ppStackCombineFinal(pmReadout *target, ppStackThreadData *stack, psArray *covariances,

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

    r25964 r26007  
    1010#include "ppStackLoop.h"
    1111
    12 #define TESTING                         // Enable test output
     12//#define TESTING                         // Enable test output
    1313
    1414bool ppStackCombineInitial(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config)

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

    r25955 r26007  
    99#include "ppStack.h"
    1010#include "ppStackLoop.h"
     11
     12//#define TESTING
     13
    1114
    1215// Update the value of a concept
    … …  
    130133        ppStackWriteImage(options->convMasks->data[i], maskHeader, readout->mask, config);
    131134        psFree(maskHeader);
    132         psImageCovarianceTransfer(readout->variance, readout->covariance);
    133135        ppStackWriteImage(options->convVariances->data[i], hdu->header, readout->variance, config);
    134136#ifdef TESTING
    … …  
    139141            pmStackVisualPlotTestImage(readout->covariance->image, name);
    140142            psFree(name);
     143        }
     144        {
     145            int numCols = readout->image->numCols, numRows = readout->image->numRows;
     146            psImage *sn = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     147            for (int y = 0; y < numRows; y++) {
     148                for (int x = 0; x < numCols; x++) {
     149                    sn->data.F32[y][x] = readout->image->data.F32[y][x] /
     150                        sqrtf(readout->variance->data.F32[y][x]);
     151                }
     152            }
     153            psString name = NULL;
     154            psStringAppend(&name, "signoise_%d.fits", i);
     155            ppStackWriteImage(name, hdu->header, sn, config);
     156            psFree(name);
     157            psFree(sn);
    141158        }
    142159#endif
    … …  
    221238            numGood = 0;                    // Number of good images
    222239            for (int i = 0; i < num; i++) {
    223               if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_ALL) {
     240                if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_ALL) {
    224241                    continue;
    225242                }
    … …  
    245262    // Correct chi^2 for renormalisation
    246263    psBinaryOp(options->matchChi2, options->matchChi2, "/", renorms);
     264    for (int i = 0; i < num; i++) {
     265        psLogMsg("ppStack", PS_LOG_INFO, "Additional variance for image %d: %f\n",
     266                 i, options->matchChi2->data.F32[i]);
     267    }
    247268    psFree(renorms);
    248269

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

    r25964 r26007  
    130130        }
    131131        psTrace("ppStack", 2, "Stack of unconvolved images....\n");
    132         if (!ppStackCombineFinal(options->unconvRO, stack, options->origCovars, options, config, true, true)) {
     132        if (!ppStackCombineFinal(options->unconvRO, stack, options->origCovars, options, config, false, true)) {
    133133            psError(PS_ERR_UNKNOWN, false, "Unable to perform unconvolved combination.");
    134134            psFree(stack);

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

    r25959 r26007  
    1818#define COVAR_FRAC 0.01                 // Truncation fraction for covariance matrix
    1919
    20 #define TESTING                         // Enable debugging output
     20//#define TESTING                         // Enable debugging output
    2121
    2222#ifdef TESTING
    … …  
    167167    )
    168168{
     169#if 1
    169170    bool mdok; // Status of metadata lookups
    170171
    … …  
    192193    psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask
    193194
     195    psImageCovarianceTransfer(readout->variance, readout->covariance);
    194196    return pmReadoutVarianceRenormalise(readout, maskBad, num, minValid, maxValid);
     197#else
     198    return true;
     199#endif
    195200}
    196201
    … …  
    212217    int size = psMetadataLookupS32(NULL, ppsub, "KERNEL.SIZE"); // Kernel half-size
    213218
    214     psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in
     219    psString maskValStr = psMetadataLookupStr(NULL, ppsub, "MASK.VAL"); // Name of bits to mask going in
    215220    psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
    216221    psString maskPoorStr = psMetadataLookupStr(NULL, recipe, "MASK.POOR"); // Name of bits to mask for poor
    … …  
    377382            }
    378383#endif
     384
     385            fprintf(stderr, "vf = %f\n", psImageCovarianceFactor(readout->covariance));
     386
    379387
    380388            if (threads > 0) {
    … …  
    516524            psFree(iter);
    517525            options->matchChi2->data.F32[index] = sum / (psImageCovarianceFactor(readout->covariance) * num);
     526            fprintf(stderr, "chi2 = %f ; vf = %f\n", sum/num, psImageCovarianceFactor(readout->covariance));
    518527        }
    519528
    … …  
    542551    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
    543552    if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) {
    544       psWarning("Can't measure background for image.");
    545       psErrorClear();
     553        psWarning("Can't measure background for image.");
     554        psErrorClear();
    546555    } else {
    547       if (!psMetadataLookupBool(NULL, config->arguments, "PPSTACK.SKIP.BG.SUB")) {
    548         psLogMsg("ppStack", PS_LOG_INFO, "Correcting convolved image background by %lf (+/- %lf)",
    549                  psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), psStatsGetValue(bg, PS_STAT_ROBUST_STDEV));
    550         (void)psBinaryOp(readout->image, readout->image, "-",
    551                          psScalarAlloc(psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), PS_TYPE_F32));
    552       }
     556        if (!psMetadataLookupBool(NULL, config->arguments, "PPSTACK.SKIP.BG.SUB")) {
     557            psLogMsg("ppStack", PS_LOG_INFO, "Correcting convolved image background by %lf (+/- %lf)",
     558                     psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), psStatsGetValue(bg, PS_STAT_ROBUST_STDEV));
     559            (void)psBinaryOp(readout->image, readout->image, "-",
     560                             psScalarAlloc(psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), PS_TYPE_F32));
     561        }
    553562    }
    554563

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

    r25968 r26007  
    116116
    117117    bool mdok;                          // Status of MD lookup
    118     int iter = psMetadataLookupS32(NULL, recipe, "ITER"); // Rejection iterations
     118    float iter = psMetadataLookupF32(NULL, recipe, "COMBINE.ITER"); // Rejection iterations
    119119    float combineRej = psMetadataLookupF32(NULL, recipe, "COMBINE.REJ"); // Combination threshold
    120120    float combineSys = psMetadataLookupF32(NULL, recipe, "COMBINE.SYS"); // Combination systematic error
    … …  
    126126    int kernelSize = psMetadataLookupS32(NULL, ppsub, "KERNEL.SIZE"); // Kernel half-size
    127127
    128     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
    129129    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
    130132    psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad
    131133    psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
    … …  
    161163    }
    162164
    163     if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, kernelSize, iter,
     165    if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskSuspect, maskBad, kernelSize, iter,
    164166                        combineRej, combineSys, combineDiscard, useVariance, safe, false)) {
    165167        psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts with rejection.");
    … …  
    217219    bool safe = psMetadataLookupBool(&mdok, recipe, "SAFE"); // Be safe when combining small numbers of pixels
    218220
    219     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
    220222    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
    221225    psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad
    222226    psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
    … …  
    267271    }
    268272
    269     if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, 0, iter, combineRej,
     273    if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskSuspect, maskBad, 0, iter, combineRej,
    270274                        combineSys, combineDiscard, useVariance, safe, true)) {
    271275        psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts.");

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

    r25964 r26007  
    1010#include "ppStackLoop.h"
    1111
    12 #define TESTING
     12//#define TESTING
    1313
    1414bool ppStackReject(ppStackOptions *options, pmConfig *config)

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

    r25975 r26007  
    3434
    3535
    36 #define TESTING                         // Enable test output
    37 #define TEST_X 4177-1                     // x coordinate to examine
    38 #define TEST_Y 2964-1                     // y coordinate to examine
     36//#define TESTING                         // Enable test output
     37//#define TEST_X 3122-1                     // x coordinate to examine
     38//#define TEST_Y 1028-1                     // y coordinate to examine
    3939
    4040
    … …  
    4343typedef struct {
    4444    psVector *pixels;                   // Pixel values
    45     psVector *masks;                    // Pixel masks
    4645    psVector *variances;                // Pixel variances
    4746    psVector *weights;                  // Pixel weightings
    4847    psVector *sources;                  // Pixel sources (which image did they come from?)
    4948    psVector *limits;                   // Rejection limits
     49    psVector *suspects;                 // Pixel is suspect?
    5050    psVector *sort;                     // Buffer for sorting (to get a robust estimator of the standard dev)
    5151} combineBuffer;
    … …  
    5454{
    5555    psFree(buffer->pixels);
    56     psFree(buffer->masks);
    5756    psFree(buffer->variances);
    5857    psFree(buffer->weights);
    5958    psFree(buffer->sources);
    6059    psFree(buffer->limits);
     60    psFree(buffer->suspects);
    6161    psFree(buffer->sort);
    6262    return;
    … …  
    7070
    7171    buffer->pixels = psVectorAlloc(numImages, PS_TYPE_F32);
    72     buffer->masks = psVectorAlloc(numImages, PS_TYPE_VECTOR_MASK);
    7372    buffer->variances = psVectorAlloc(numImages, PS_TYPE_F32);
    7473    buffer->weights = psVectorAlloc(numImages, PS_TYPE_F32);
    7574    buffer->sources = psVectorAlloc(numImages, PS_TYPE_U16);
    7675    buffer->limits = psVectorAlloc(numImages, PS_TYPE_F32);
     76    buffer->suspects = psVectorAlloc(numImages, PS_TYPE_U8);
    7777    buffer->sort = psVectorAlloc(numImages, PS_TYPE_F32);
    7878    return buffer;
    … …  
    144144                                   float *stdev, // Standard deviation value, to return
    145145                                   const psVector *values, // Values to combine
    146                                    const psVector *masks, // Mask to apply
    147146                                   psVector *sortBuffer // Buffer for sorting
    148147                                   )
    149148{
    150149    assert(values);
    151     assert(!masks || values->n == masks->n);
    152150    assert(values->type.type == PS_TYPE_F32);
    153     assert(!masks || masks->type.type == PS_TYPE_VECTOR_MASK);
    154151    assert(sortBuffer && sortBuffer->nalloc >= values->n && sortBuffer->type.type == PS_TYPE_F32);
    155152
    156     // Need to filter out clipped values
    157     int num = 0;            // Number of valid values
    158     for (int i = 0; i < values->n; i++) {
    159         if (!masks || !masks->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
    160             sortBuffer->data.F32[num++] = values->data.F32[i];
    161         }
    162     }
    163     sortBuffer->n = num;
    164     if (!psVectorSortInPlace(sortBuffer)) {
     153    int num = values->n;                // Number of values
     154    sortBuffer = psVectorSortIndex(sortBuffer, values);
     155    if (!sortBuffer) {
     156        *median = NAN;
     157        *stdev = NAN;
    165158        return false;
    166159    }
    … …  
    168161    if (num == 3) {
    169162        // Attempt to measure standard deviation with only three values (and one of those possibly corrupted)
    170         *median = sortBuffer->data.F32[1];
     163        *median = values->data.F32[sortBuffer->data.S32[1]];
    171164        if (stdev) {
    172             float diff1 = sortBuffer->data.F32[0] - *median;
    173             float diff2 = sortBuffer->data.F32[2] - *median;
     165            float diff1 = values->data.F32[sortBuffer->data.S32[0]] - *median;
     166            float diff2 = values->data.F32[sortBuffer->data.S32[2]] - *median;
    174167            // This factor of sqrt(2) might not be exact, but it's about right
    175168            *stdev = M_SQRT2 * PS_MIN(fabsf(diff1), fabsf(diff2));
    176169        }
    177170    } else {
    178         *median = num % 2 ? sortBuffer->data.F32[num / 2] :
    179             (sortBuffer->data.F32[num / 2 - 1] + sortBuffer->data.F32[num / 2]) / 2.0;
     171        *median = num % 2 ? values->data.F32[sortBuffer->data.S32[num / 2]] :
     172            (values->data.F32[sortBuffer->data.S32[num / 2 - 1]] +
     173             values->data.F32[sortBuffer->data.S32[num / 2]]) / 2.0;
    180174        if (stdev) {
    181175            if (num <= NUM_DIRECT_STDEV) {
    … …  
    183177                double sum = 0.0;
    184178                for (int i = 0; i < num; i++) {
    185                     sum += PS_SQR(sortBuffer->data.F32[i] - *median);
     179                    sum += PS_SQR(values->data.F32[sortBuffer->data.S32[i]] - *median);
    186180                }
    187181                *stdev = sqrt(sum / (double)(num - 1));
    188182            } else {
    189183                // Standard deviation from the interquartile range
    190                 *stdev = 0.74 * (sortBuffer->data.F32[(int)(0.75 * num)] -
    191                                  sortBuffer->data.F32[(int)(0.25 * num)]);
     184                *stdev = 0.74 * (values->data.F32[sortBuffer->data.S32[(int)(0.75 * num)]] -
     185                                 values->data.F32[sortBuffer->data.S32[(int)(0.25 * num)]]);
    192186            }
    193187        }
    … …  
    196190    return true;
    197191}
    198 
    199 #if 0
    200 // Return the weighted median for the pixels
    201 // This does not appear to produce as clean images as the weighted Olympic mean
    202 static float combinationWeightedMedian(const psVector *values, // Values to combine
    203                                        const psVector *weights, // Weights to combine
    204                                        const psVector *masks, // Mask to apply
    205                                        psVector *sortBuffer // Buffer for sorting
    206                                        )
    207 {
    208     double sumWeight = 0.0;             // Sum of weights
    209     for (int j = 0; j < values->n; j++) {
    210         if (masks->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
    211             continue;
    212         }
    213         sumWeight += weights->data.F32[j];
    214     }
    215 
    216     sortBuffer = psVectorSortIndex(sortBuffer, values);
    217     double target = sumWeight / 2.0;    // Target weight
    218 
    219     int dominant = -1;                  // Index of dominant value, if any
    220     double cumulativeWeight = 0.0;      // Sum of weights
    221     for (int j = 0; j < values->n; j++) {
    222         if (masks->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
    223             continue;
    224         }
    225         int index = sortBuffer->data.S32[j];  // Index of value of interest
    226         float weight = weights->data.F32[index]; // Weight for value of interest
    227         if (weight >= target) {
    228             // Get the weighted median of the rest
    229             dominant = index;
    230             sumWeight -= weight;
    231             target = sumWeight / 2.0;
    232             continue;
    233         }
    234         cumulativeWeight += weight;
    235         if (cumulativeWeight >= target) {
    236             float median = values->data.F32[index]; // Weighted median median
    237             if (dominant != -1) {
    238                 // In the case that a single value contains a disproportionate weight compared to the rest,
    239                 // we use a weighted mean between that dominant value and the weighted median of the rest.
    240                 return (values->data.F32[dominant] * weights->data.F32[dominant] + median * sumWeight) /
    241                     (weights->data.F32[dominant] + sumWeight);
    242             }
    243             return median;
    244         }
    245     }
    246 
    247     return NAN;
    248 }
    249 #endif
    250192
    251193// Return the weighted Olympic mean for the pixels
    252194static float combinationWeightedOlympic(const psVector *values, // Values to combine
    253195                                        const psVector *weights, // Weights to combine
    254                                         const psVector *masks, // Mask to apply
    255196                                        float frac, // Fraction to discard
    256197                                        psVector *sortBuffer // Buffer for sorting
    257198                                        )
    258199{
    259     int numGood = 0;                    // Number of good values
    260     for (int i = 0; i < values->n; i++) {
    261         if (masks->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
    262             continue;
    263         }
    264         numGood++;
    265     }
     200    int numGood = values->n;            // Number of good values
    266201
    267202    int numBad = frac * numGood + 0.5;  // Number of bad values
    … …  
    273208    for (int i = 0, j = 0; i < values->n; i++) {
    274209        int index = sortBuffer->data.S32[i]; // Index of interest
    275         if (masks->data.PS_TYPE_VECTOR_MASK_DATA[index]) {
    276             continue;
    277         }
    278210        j++;
    279211        if (j > high) {
    … …  
    297229                                      )
    298230{
     231#ifdef TESTING
     232    if (x == TEST_X && y == TEST_Y) {
     233        fprintf(stderr, "Marking image %d for inspection\n", source);
     234    }
     235#endif
    299236    pmStackData *data = inputs->data[source]; // Stack data of interest
    300237    if (!data) {
    … …  
    314251                                     )
    315252{
     253#ifdef TESTING
     254    if (x == TEST_X && y == TEST_Y) {
     255        fprintf(stderr, "Marking pixel image %d for rejection\n", source);
     256    }
     257#endif
    316258    pmStackData *data = inputs->data[source]; // Stack data of interest
    317259    if (!data) {
    … …  
    321263    data->reject = psPixelsAdd(data->reject, data->reject->nalloc, x, y);
    322264    return;
    323 }
    324 
    325 // General test for multiple pixels
    326 // Returns false if we need to re-run without suspect pixels
    327 static bool combineTestGeneral(int num,      // Number of good pixels
    328                                bool suspect, // Are there suspect pixels in the list?
    329                                psArray *inputs,       // Original inputs (for flagging)
    330                                combineBuffer *buffer, // Buffer with vectors
    331                                int x, int y, // Coordinates of interest; frame of output image
    332                                int numIter, // Number of rejection iterations
    333                                float rej, // Number of standard deviations at which to reject
    334                                float sys,    // Relative systematic error
    335                                float olympic,// Fraction of values to discard (Olympic weighted mean)
    336                                bool useVariance, // Use variance for rejection when combining?
    337                                bool safe,    // Combine safely?
    338                                bool allowSuspect // Allow suspect values?
    339                                )
    340 {
    341     psVector *pixelData = buffer->pixels; // Values for the pixel of interest
    342     psVector *pixelMasks = buffer->masks; // Masks for the pixel of interest
    343     psVector *pixelVariances = buffer->variances; // Variances for the pixel of interest
    344     psVector *pixelWeights = buffer->weights; // Image weights for the pixel of interest
    345     psVector *pixelSources = buffer->sources; // Sources for the pixel of interest
    346     psVector *pixelLimits = buffer->limits; // Limits for the pixel of interest
    347     psVector *sort = buffer->sort;      // Sort buffer
    348 
    349 
    350     pixelMasks->n = num;
    351     psVectorInit(pixelMasks, 0);
    352 
    353     // Set up rejection limits
    354     if (useVariance) {
    355         // Convert to rejection limits --- saves doing it later.
    356         // Using squared rejection limit because it's cheaper than sqrts
    357         float rej2 = PS_SQR(rej); // Rejection level squared
    358         double sumWeights = 0.0;
    359         for (int i = 0; i < num; i++) {
    360             sumWeights += pixelWeights->data.F32[i];
    361         }
    362         for (int i = 0; i < num; i++) {
    363             // Systematic error contributes to the rejection level
    364             float sysVar = PS_SQR(sys * pixelData->data.F32[i]);
    365 #ifdef TESTING
    366             // Correct variance for comparison against weighted mean including itself
    367             float compare = 1.0 - pixelWeights->data.F32[i] / sumWeights;
    368             if (x == TEST_X && y == TEST_Y) {
    369                 fprintf(stderr, "Variance %d (%d): %f %f %f\n", i, pixelSources->data.U16[i],
    370                         pixelVariances->data.F32[i], sysVar, compare);
    371             }
    372 #endif
    373             pixelLimits->data.F32[i] = rej2 * (pixelVariances->data.F32[i] + sysVar);
    374         }
    375     }
    376 
    377     int numClipped = INT_MAX;     // Number of pixels clipped per iteration
    378     int totalClipped = 0;         // Total number of pixels clipped
    379     for (int i = 0; i < numIter && numClipped > 0 && num - totalClipped > 2; i++) {
    380         numClipped = 0;
    381         float median = NAN;       // Middle of distribution
    382         float limit = NAN;        // Rejection limit
    383         if (!useVariance) {
    384             float stdev;  // Median and stdev of the combination, for rejection
    385             if (!combinationMedianStdev(&median, useVariance ? NULL : &stdev,
    386                                         pixelData, pixelMasks, sort)) {
    387                 if (i == 0 && suspect) {
    388                     // Something's not right --- repeat
    389                     return false;
    390                 }
    391                 for (int j = 0; j < num; j++) {
    392                     combineMarkReject(inputs, x, y, pixelSources->data.U16[j]);
    393                 }
    394                 return true;
    395             }
    396             limit = rej * stdev;
    397 #ifdef TESTING
    398             if (x == TEST_X && y == TEST_Y) {
    399                 fprintf(stderr, "Flag without variance; limit: %f\n", limit);
    400             }
    401 #endif
    402         } else {
    403 #ifdef TESTING
    404             if (x == TEST_X && y == TEST_Y) {
    405                 fprintf(stderr, "Flag with variance...\n");
    406             }
    407 #endif
    408             median = combinationWeightedOlympic(pixelData, pixelWeights, pixelMasks, olympic, sort);
    409         }
    410 
    411 #ifdef TESTING
    412         if (x == TEST_X && y == TEST_Y) {
    413             fprintf(stderr, "Median: %f\n", median);
    414         }
    415 #endif
    416 
    417         // Mask a pixel for inspection
    418 #define MASK_PIXEL_FOR_INSPECTION()                                     \
    419         if (i == 0 && suspect) {                                        \
    420             /* Something's inconsistent, so want repeat, throwing out all suspect pixels */ \
    421             return false; \
    422         }                                                               \
    423         pixelMasks->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xFF;            \
    424         combineMarkInspect(inputs, x, y, pixelSources->data.U16[j]);    \
    425         numClipped++;                                                   \
    426         totalClipped++;
    427         // End
    428 
    429         for (int j = 0; j < num; j++) {
    430             if (pixelMasks->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
    431                 continue;
    432             }
    433             float diff = pixelData->data.F32[j] - median; // Difference from expected
    434 #ifdef TESTING
    435             if (x == TEST_X && y == TEST_Y) {
    436                 fprintf(stderr, "Testing input %d: %f\n", j, diff);
    437             }
    438 #endif
    439             if (useVariance) {
    440                 // Comparing squares --- cheaper than lots of sqrts
    441                 // pixelVariances includes the rejection limit, from above
    442                 if (PS_SQR(diff) > pixelLimits->data.F32[j]) {
    443                     MASK_PIXEL_FOR_INSPECTION();
    444 #ifdef TESTING
    445                     if (x == TEST_X && y == TEST_Y) {
    446                         fprintf(stderr, "Flagging input %d based on variance: %f > %f\n",
    447                                 j, diff, sqrtf(pixelLimits->data.F32[j]));
    448                     }
    449 #endif
    450                 }
    451             } else if (fabsf(diff) > limit) {
    452                 MASK_PIXEL_FOR_INSPECTION();
    453 #ifdef TESTING
    454                 if (x == TEST_X && y == TEST_Y) {
    455                     fprintf(stderr, "Flagging input %d based on distribution: %f > %f\n",
    456                             j, diff, limit);
    457                 }
    458 #endif
    459             }
    460         }
    461     }
    462 
    463     return true;
    464265}
    465266
    … …  
    478279                           int x, int y, // Coordinates of interest; frame of output image
    479280                           psImageMaskType maskVal, // Value to mask
    480                            psImageMaskType maskSuspect, // Value to suspect
    481                            bool allowSuspect        // Allow suspect values?
     281                           psImageMaskType maskSuspect // Value to suspect
    482282                           )
    483283{
    … …  
    493293    psVector *pixelSources = buffer->sources; // Sources for the pixel of interest
    494294    psVector *pixelLimits = buffer->limits; // Limits for the pixel of interest
     295    psVector *pixelSuspects = buffer->suspects; // Is the pixel suspect?
     296
     297    if (suspect) {
     298        *suspect = false;
     299    }
    495300
    496301    // Extract the pixel and mask data
    … …  
    514319            continue;
    515320        }
    516         if (mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & maskSuspect) {
    517             if (!allowSuspect) {
    518                 combineMarkReject(inputs, x, y, i);
    519                 continue;
    520             }
    521             if (suspect) {
    522                 *suspect = true;
    523             }
    524         }
     321
     322        pixelSuspects->data.U8[numGood] = mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & maskSuspect ?
     323            true : false;
    525324
    526325        psImage *image = data->readout->image; // Image of interest
    … …  
    541340    pixelSources->n = numGood;
    542341    pixelLimits->n = numGood;
     342    pixelSuspects->n = numGood;
    543343    *num = numGood;
    544344
    … …  
    546346    if (x == TEST_X && y == TEST_Y) {
    547347        for (int i = 0; i < numGood; i++) {
    548             fprintf(stderr, "Input %d (%" PRIu16 "): %f %f (%f) %f\n",
     348            fprintf(stderr, "Input %d (%" PRIu16 "): %f %f (%f) %f %d\n",
    549349                    i, pixelSources->data.U16[i], pixelData->data.F32[i], pixelVariances->data.F32[i],
    550                     addVariance->data.F32[i], pixelWeights->data.F32[i]);
     350                    addVariance->data.F32[i], pixelWeights->data.F32[i], pixelSuspects->data.U8[i]);
    551351        }
    552352    }
    … …  
    654454    }
    655455
     456#ifdef TESTING
     457    mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = num;
     458#endif
     459
     460
    656461    return;
    657462}
    … …  
    665470                        combineBuffer *buffer, // Buffer with vectors
    666471                        int x, int y, // Coordinates of interest; frame of output image
    667                         int numIter, // Number of rejection iterations
     472                        float iter, // Number of rejection iterations per input
    668473                        float rej, // Number of standard deviations at which to reject
    669474                        float sys,    // Relative systematic error
    670475                        float olympic,// Fraction of values to discard (Olympic weighted mean)
    671476                        bool useVariance, // Use variance for rejection when combining?
    672                         bool safe,    // Combine safely?
    673                         bool allowSuspect    // Allow suspect pixels in the combination?
     477                        bool safe    // Combine safely?
    674478                        )
    675479{
    676     if (numIter <= 0) {
     480    if (iter <= 0) {
    677481        return true;
    678482    }
    679483
     484    int numIter = PS_MAX(iter * num, 1); // Number of iterations
     485
     486#ifdef TESTING
     487    if (x == TEST_X && y == TEST_Y) {
     488        fprintf(stderr, "Testing pixel %d,%d: %d %f %f %f %d %d\n",
     489                x, y, numIter, rej, sys, olympic, useVariance, safe);
     490    }
     491#endif
     492
    680493    psVector *pixelData = buffer->pixels; // Values for the pixel of interest
     494    psVector *pixelWeights = buffer->weights; // Is the pixel suspect?
    681495    psVector *pixelVariances = buffer->variances; // Variances for the pixel of interest
    682496    psVector *pixelSources = buffer->sources; // Sources for the pixel of interest
    683 
    684     switch (num) {
    685       case 0:
    686         break;
    687       case 1:
    688           if (safe) {
    689               combineMarkReject(inputs, x, y, pixelSources->data.U16[0]);
     497    psVector *pixelSuspects = buffer->suspects; // Is the pixel suspect?
     498    psVector *pixelLimits = buffer->limits; // Is the pixel suspect?
     499
     500    // Set up rejection limits
     501    float rej2 = PS_SQR(rej); // Rejection level squared
     502    if (num > 2 && useVariance) {
     503        // Convert rejection limits --- saves doing it later multiple times
     504        // Using squared rejection limit because it's cheaper than sqrts
     505        double sumWeights = 0.0;
     506        for (int i = 0; i < num; i++) {
     507            sumWeights += pixelWeights->data.F32[i];
     508        }
     509        for (int i = 0; i < num; i++) {
     510            // Systematic error contributes to the rejection level
     511            float sysVar = PS_SQR(sys * pixelData->data.F32[i]);
     512#ifdef TESTING
     513            // Correct variance for comparison against weighted mean including itself
     514            float compare = 1.0 - pixelWeights->data.F32[i] / sumWeights;
     515            if (x == TEST_X && y == TEST_Y) {
     516                fprintf(stderr, "Variance %d (%d): %f %f %f\n", i, pixelSources->data.U16[i],
     517                        pixelVariances->data.F32[i], sysVar, compare);
     518            }
     519#endif
     520            pixelLimits->data.F32[i] = rej2 * (pixelVariances->data.F32[i] + sysVar);
     521        }
     522    }
     523
     524    int maskIndex = 0;                  // Index of pixel to mask
     525    int totalClipped = 0;               // Total number of pixels clipped
     526    for (int i = 0; i < numIter && maskIndex >= 0; i++) {
     527        maskIndex = -1;                 // Nothing to reject
     528
     529        switch (num) {
     530          case 0:
     531            break;
     532          case 1:
     533            if (i == 0 && safe) {
     534                combineMarkReject(inputs, x, y, pixelSources->data.U16[0]);
     535            }
     536            break;
     537          case 2: {
     538              if (useVariance) {
     539                  // Use variance to check that the two are consistent
     540                  float diff = 0.5 * (pixelData->data.F32[0] - pixelData->data.F32[1]); // Mean flux difference
     541                  float var1 = pixelVariances->data.F32[0]; // Variance of first
     542                  float var2 = pixelVariances->data.F32[1]; // Variance of second
     543                  // Systematic error contributes to the rejection level
     544                  var1 += PS_SQR(sys * pixelData->data.F32[0]);
     545                  var2 += PS_SQR(sys * pixelData->data.F32[1]);
     546
     547                  float sigma2 = var1 + var2; // Combined variance
     548                  if (PS_SQR(diff) > rej2 * sigma2) {
     549                      // Not consistent: don't believe either!
     550                      if (i == 0 && suspect) {
     551                          combineMarkReject(inputs, x, y, pixelSources->data.U16[0]);
     552                          combineMarkReject(inputs, x, y, pixelSources->data.U16[1]);
     553                      } else {
     554                          combineMarkInspect(inputs, x, y, pixelSources->data.U16[0]);
     555                          combineMarkInspect(inputs, x, y, pixelSources->data.U16[1]);
     556                      }
     557#ifdef TESTING
     558                      if (x == TEST_X && y == TEST_Y) {
     559                          fprintf(stderr, "Flagged both pixels (%f > %f x %f\n)",
     560                                  diff, rej, sqrtf(sigma2));
     561                      }
     562#endif
     563                  }
     564              } else if (i == 0 && safe) {
     565                  // Can't test them, and we want to be safe, so reject
     566                  combineMarkReject(inputs, x, y, pixelSources->data.U16[0]);
     567                  combineMarkReject(inputs, x, y, pixelSources->data.U16[1]);
     568              }
     569              break;
    690570          }
    691           break;
    692       case 2: {
    693           if (useVariance) {
    694               // Use variance to check that the two are consistent
    695               float diff = 0.5 * (pixelData->data.F32[0] - pixelData->data.F32[1]); // Mean flux difference
    696               float var1 = pixelVariances->data.F32[0]; // Variance of first
    697               float var2 = pixelVariances->data.F32[1]; // Variance of second
    698               // Systematic error contributes to the rejection level
    699               var1 += PS_SQR(sys * pixelData->data.F32[0]);
    700               var2 += PS_SQR(sys * pixelData->data.F32[1]);
    701 
    702               float sigma2 = var1 + var2; // Combined variance
    703               if (PS_SQR(diff) > PS_SQR(rej) * sigma2) {
    704                   // Not consistent: don't believe either!
    705                   if (allowSuspect && suspect) {
    706                       combineMarkReject(inputs, x, y, pixelSources->data.U16[0]);
    707                       combineMarkReject(inputs, x, y, pixelSources->data.U16[1]);
    708                   } else {
     571#if 0
     572          case 3: {
     573              // Want to be a bit careful on the rejection than for a larger number of inputs
     574              if (!useVariance) {
     575                  return combineTestGeneral(num, suspect, inputs, buffer, x, y, numIter, rej, sys,
     576                                            olympic, useVariance, safe, allowSuspect);
     577              }
     578
     579              // Differences between pixel values
     580              float diff01 = pixelData->data.F32[0] - pixelData->data.F32[1];
     581              float diff12 = pixelData->data.F32[1] - pixelData->data.F32[2];
     582              float diff20 = pixelData->data.F32[2] - pixelData->data.F32[0];
     583              // Variance for each pixel
     584              float var0 = pixelVariances->data.F32[0] + PS_SQR(sys * pixelData->data.F32[0]);
     585              float var1 = pixelVariances->data.F32[1] + PS_SQR(sys * pixelData->data.F32[1]);
     586              float var2 = pixelVariances->data.F32[2] + PS_SQR(sys * pixelData->data.F32[2]);
     587              // Errors in pixel differences
     588              float err01 = var0 + var1;
     589              float err12 = var1 + var2;
     590              float err20 = var2 + var0;
     591
     592#ifdef TESTING
     593              if (x == TEST_X && y == TEST_Y) {
     594                  fprintf(stderr, "Diff 0-1: %f %f\n", diff01, err01);
     595                  fprintf(stderr, "Diff 1-2: %f %f\n", diff12, err12);
     596                  fprintf(stderr, "Diff 2-0: %f %f\n", diff20, err20);
     597              }
     598#endif
     599
     600              int badPairs = 0;         // Number of bad pairs
     601              bool bad01 = false, bad12 = false, bad20 = false; // Pair is bad?
     602              if (PS_SQR(diff01) > rej2 * err01) {
     603                  bad01 = true;
     604                  badPairs++;
     605              }
     606              if (PS_SQR(diff12) > rej2 * err12) {
     607                  bad12 = true;
     608                  badPairs++;
     609              }
     610              if (PS_SQR(diff20) > rej2 * err20) {
     611                  bad20 = true;
     612                  badPairs++;
     613              }
     614
     615              if (badPairs > 0 && allowSuspect && suspect) {
     616                  return false;
     617              }
     618
     619              switch (badPairs) {
     620                case 0:
     621                  // Nothing to worry about!
     622                  break;
     623                case 1:
     624                  // Can't tell which image is bad, so be sure to get it
     625                  if (bad01) {
    709626                      combineMarkInspect(inputs, x, y, pixelSources->data.U16[0]);
    710627                      combineMarkInspect(inputs, x, y, pixelSources->data.U16[1]);
     628                      break;
    711629                  }
    712 #ifdef TESTING
    713                   if (x == TEST_X && y == TEST_Y) {
    714                       fprintf(stderr, "Flagged both pixels (%f > %f x %f\n)",
    715                               diff, rej, sqrtf(sigma2));
     630                  if (bad12) {
     631                      combineMarkInspect(inputs, x, y, pixelSources->data.U16[1]);
     632                      combineMarkInspect(inputs, x, y, pixelSources->data.U16[2]);
     633                      break;
    716634                  }
    717 #endif
    718               }
    719           } else if (safe) {
    720               // Can't test them, and we want to be safe, so reject
    721               combineMarkReject(inputs, x, y, pixelSources->data.U16[0]);
    722               combineMarkReject(inputs, x, y, pixelSources->data.U16[1]);
    723           }
    724           break;
    725       }
    726       case 3: {
    727           // Want to be a bit careful on the rejection than for a larger number of inputs
    728           if (!useVariance) {
    729               return combineTestGeneral(num, suspect, inputs, buffer, x, y, numIter, rej, sys,
    730                                         olympic, useVariance, safe, allowSuspect);
    731           }
    732 
    733           // Differences between pixel values
    734           float diff01 = pixelData->data.F32[0] - pixelData->data.F32[1];
    735           float diff12 = pixelData->data.F32[1] - pixelData->data.F32[2];
    736           float diff20 = pixelData->data.F32[2] - pixelData->data.F32[0];
    737           // Variance for each pixel
    738           float var0 = pixelVariances->data.F32[0] + PS_SQR(sys * pixelData->data.F32[0]);
    739           float var1 = pixelVariances->data.F32[1] + PS_SQR(sys * pixelData->data.F32[1]);
    740           float var2 = pixelVariances->data.F32[2] + PS_SQR(sys * pixelData->data.F32[2]);
    741           float rej2 = PS_SQR(rej); // Rejection level squared
    742           // Errors in pixel differences
    743           float err01 = var0 + var1;
    744           float err12 = var1 + var2;
    745           float err20 = var2 + var0;
    746 
    747           int badPairs = 0;         // Number of bad pairs
    748           bool bad01 = false, bad12 = false, bad20 = false; // Pair is bad?
    749           if (PS_SQR(diff01) > rej2 * err01) {
    750               bad01 = true;
    751               badPairs++;
    752           }
    753           if (PS_SQR(diff12) > rej2 * err12) {
    754               bad12 = true;
    755               badPairs++;
    756           }
    757           if (PS_SQR(diff20) > rej2 * err20) {
    758               bad20 = true;
    759               badPairs++;
    760           }
    761 
    762           if (badPairs > 0 && allowSuspect && suspect) {
    763               return false;
    764           }
    765 
    766           switch (badPairs) {
    767             case 0:
    768               // Nothing to worry about!
    769               break;
    770             case 1:
    771               // Can't tell which image is bad, so be sure to get it
    772               if (bad01) {
     635                  if (bad20) {
     636                      combineMarkInspect(inputs, x, y, pixelSources->data.U16[2]);
     637                      combineMarkInspect(inputs, x, y, pixelSources->data.U16[0]);
     638                      break;
     639                  }
     640                  psAbort("Should never get here");
     641                case 2:
     642                  if (bad01 && bad12) {
     643                      // 2 and 0 are good
     644                      combineMarkInspect(inputs, x, y, pixelSources->data.U16[1]);
     645                      break;
     646                  }
     647                  if (bad12 && bad20) {
     648                      // 0 and 1 are good
     649                      combineMarkInspect(inputs, x, y, pixelSources->data.U16[2]);
     650                      break;
     651                  }
     652                  if (bad20 && bad01) {
     653                      // 1 and 2 are good
     654                      combineMarkInspect(inputs, x, y, pixelSources->data.U16[0]);
     655                      break;
     656                  }
     657                  psAbort("Should never get here");
     658                case 3:
     659                  // Everything's bad
    773660                  combineMarkInspect(inputs, x, y, pixelSources->data.U16[0]);
    774                   combineMarkInspect(inputs, x, y, pixelSources->data.U16[1]);
    775                   break;
    776               }
    777               if (bad12) {
    778661                  combineMarkInspect(inputs, x, y, pixelSources->data.U16[1]);
    779662                  combineMarkInspect(inputs, x, y, pixelSources->data.U16[2]);
    780663                  break;
    781664              }
    782               if (bad20) {
    783                   combineMarkInspect(inputs, x, y, pixelSources->data.U16[2]);
    784                   combineMarkInspect(inputs, x, y, pixelSources->data.U16[0]);
    785                   break;
    786               }
    787               psAbort("Should never get here");
    788             case 2:
    789               if (bad01 && bad12) {
    790                   // 2 and 0 are good
    791                   combineMarkInspect(inputs, x, y, pixelSources->data.U16[1]);
    792                   break;
    793               }
    794               if (bad12 && bad20) {
    795                   // 0 and 1 are good
    796                   combineMarkInspect(inputs, x, y, pixelSources->data.U16[2]);
    797                   break;
    798               }
    799               if (bad20 && bad01) {
    800                   // 1 and 2 are good
    801                   combineMarkInspect(inputs, x, y, pixelSources->data.U16[0]);
    802                   break;
    803               }
    804               psAbort("Should never get here");
    805             case 3:
    806               // Everything's bad
    807               combineMarkInspect(inputs, x, y, pixelSources->data.U16[0]);
    808               combineMarkInspect(inputs, x, y, pixelSources->data.U16[1]);
    809               combineMarkInspect(inputs, x, y, pixelSources->data.U16[2]);
    810665              break;
    811666          }
    812           break;
    813       }
    814       default: {
    815           return combineTestGeneral(num, suspect, inputs, buffer, x, y, numIter, rej, sys,
    816                                     olympic, useVariance, safe, allowSuspect);
    817       }
     667#endif
     668          default: {
     669              if (useVariance) {
     670                  float median = combinationWeightedOlympic(pixelData, pixelWeights,
     671                                                            olympic, buffer->sort); // Median for stack
     672#ifdef TESTING
     673                  if (x == TEST_X && y == TEST_Y) {
     674                      fprintf(stderr, "Flag with variance, median = %f\n", median);
     675                  }
     676#endif
     677                  float worst = -INFINITY; // Largest deviation
     678                  for (int j = 0; j < num; j++) {
     679                      float diff = pixelData->data.F32[j] - median; // Difference from expected
     680#ifdef TESTING
     681                      if (x == TEST_X && y == TEST_Y) {
     682                          fprintf(stderr, "Testing input %d: %f\n", j, diff);
     683                      }
     684#endif
     685
     686                      // Comparing squares --- cheaper than lots of sqrts
     687                      // pixelVariances includes the rejection limit, from above
     688                      float diff2 = PS_SQR(diff); // Square difference
     689                      if (diff2 > pixelLimits->data.F32[j]) {
     690                          float dev = diff2 / pixelLimits->data.F32[j]; // Deviation
     691                          if (dev > worst) {
     692                              worst = dev;
     693                              maskIndex = j;
     694                          }
     695                      }
     696                  }
     697              } else {
     698                  float median = NAN, stdev = NAN;  // Median and stdev of the combination, for rejection
     699                  combinationMedianStdev(&median, &stdev, pixelData, buffer->sort);
     700                  float limit = rej * stdev; // Rejection limit
     701#ifdef TESTING
     702                  if (x == TEST_X && y == TEST_Y) {
     703                      fprintf(stderr, "Flag without variance; median = %f, stdev = %f, limit = %f\n",
     704                              median, stdev, limit);
     705                  }
     706#endif
     707                  float worst = -INFINITY; // Largest deviation
     708                  for (int j = 0; j < num; j++) {
     709                      float diff = fabsf(pixelData->data.F32[j] - median); // Difference from expected
     710
     711                      if (diff > limit) {
     712                          float dev = diff / limit; // Deviation
     713                          if (dev > worst) {
     714                              worst = dev;
     715                              maskIndex = j;
     716                          }
     717                      }
     718                  }
     719              }
     720          }
     721        }
     722
     723        // Do the actual rejection of the pixel
     724        if (maskIndex >= 0) {
     725            if (suspect) {
     726#ifdef TESTING
     727                if (x == TEST_X && y == TEST_Y) {
     728                    fprintf(stderr, "Throwing out all suspect pixels\n");
     729                }
     730#endif
     731                // Throw out all suspect pixels
     732                int numGood = 0;        // Number of good pixels
     733                for (int j = 0; j < num; j++) {
     734                    if (pixelSuspects->data.U8[j]) {
     735                        combineMarkReject(inputs, x, y, pixelSources->data.U16[j]);
     736                        continue;
     737                    }
     738                    if (numGood == j) {
     739                        numGood++;
     740                        continue;
     741                    }
     742                    pixelData->data.F32[numGood] = pixelData->data.F32[j];
     743                    pixelWeights->data.F32[numGood] = pixelWeights->data.F32[j];
     744                    pixelSources->data.U16[numGood] = pixelSources->data.U16[j];
     745                    pixelLimits->data.F32[numGood] = pixelLimits->data.F32[j];
     746                    pixelVariances->data.F32[numGood] = pixelVariances->data.F32[j];
     747                    numGood++;
     748                }
     749                pixelData->n = numGood;
     750                pixelWeights->n = numGood;
     751                pixelSources->n = numGood;
     752                pixelLimits->n = numGood;
     753                pixelVariances->n = numGood;
     754                totalClipped += num - numGood;
     755                num = numGood;
     756                suspect = false;
     757            } else {
     758                // Throw out masked pixel
     759#ifdef TESTING
     760                if (x == TEST_X && y == TEST_Y) {
     761                    fprintf(stderr, "Throwing out input %d\n", maskIndex);
     762                }
     763#endif
     764                combineMarkInspect(inputs, x, y, pixelSources->data.U16[maskIndex]);
     765                int numGood = 0;        // Number of good pixels
     766                for (int j = 0; j < num; j++) {
     767                    if (j == maskIndex) {
     768                        continue;
     769                    }
     770                    if (numGood == j) {
     771                        numGood++;
     772                        continue;
     773                    }
     774                    pixelData->data.F32[numGood] = pixelData->data.F32[j];
     775                    pixelWeights->data.F32[numGood] = pixelWeights->data.F32[j];
     776                    pixelSources->data.U16[numGood] = pixelSources->data.U16[j];
     777                    pixelLimits->data.F32[numGood] = pixelLimits->data.F32[j];
     778                    pixelVariances->data.F32[numGood] = pixelVariances->data.F32[j];
     779                    numGood++;
     780                }
     781                pixelData->n = numGood;
     782                pixelWeights->n = numGood;
     783                pixelSources->n = numGood;
     784                pixelLimits->n = numGood;
     785                pixelVariances->n = numGood;
     786                totalClipped++;
     787                num--;
     788            }
     789        }
    818790    }
    819791
    … …  
    821793}
    822794
    823 
    824 
    825 #if 0
    826 // Given a stack of images, combine with optional rejection.
    827 // Pixels in the stack that are rejected are marked for subsequent inspection
    828 static void combinePixels(psImage *image, // Combined image, for output
    829                           psImage *mask, // Combined mask, for output
    830                           psImage *variance, // Combined variance map, for output
    831                           const psArray *inputs, // Stack data
    832                           const psVector *weights, // Global (single value) weights for data, or NULL
    833                           const psVector *addVariance, // Additional variance for data
    834                           const psVector *reject, // Indices of pixels to reject, or NULL
    835                           int x, int y, // Coordinates of interest; frame of output image
    836                           psImageMaskType maskVal, // Value to mask
    837                           psImageMaskType suspect, // Value to suspect
    838                           psImageMaskType bad, // Value to give bad pixels
    839                           int numIter, // Number of rejection iterations
    840                           float rej, // Number of standard deviations at which to reject
    841                           float sys,    // Relative systematic error
    842                           float olympic,// Fraction of values to discard (Olympic weighted mean)
    843                           bool useVariance, // Use variance for rejection when combining?
    844                           bool safe,    // Combine safely?
    845                           bool rejection, // Reject values marked for inspection from combination?
    846                           combineBuffer *buffer, // Buffer for combination; to avoid multiple allocations
    847                           bool allowSuspect    // Allow suspect pixels in the combination?
    848                          )
    849 {
    850     // Rudimentary error checking
    851     assert(image);
    852     assert(mask);
    853     assert(inputs);
    854     assert(numIter >= 0);
    855     assert(buffer);
    856     assert(addVariance);
    857     assert((useVariance && variance) || !useVariance);
    858 
    859     psVector *pixelData = buffer->pixels; // Values for the pixel of interest
    860     psVector *pixelMasks = buffer->masks; // Masks for the pixel of interest
    861     psVector *pixelVariances = variance ? buffer->variances : NULL; // Variances for the pixel of interest
    862     psVector *pixelWeights = buffer->weights; // Image weights for the pixel of interest
    863     psVector *pixelSources = buffer->sources; // Sources for the pixel of interest
    864     psVector *pixelLimits = buffer->limits; // Limits for the pixel of interest
    865     psVector *sort = buffer->sort;      // Sort buffer
    866 
    867 
    868     // The sensible thing to do varies according to how many good pixels there are.
    869     // Default option is that the pixel is bad
    870     float imageValue = NAN, varianceValue = NAN; // Value for combined image and variance map
    871     psImageMaskType maskValue = bad;    // Value for combined mask
    872     switch (num) {
    873       case 0:
    874         // Nothing to combine: it's bad
    875 #ifdef TESTING
    876     if (x == TEST_X && y == TEST_Y) {
    877         fprintf(stderr, "No inputs to combine, pixel is bad.\n");
    878     }
    879 #endif
    880         break;
    881       case 1: {
    882           // Accept the single pixel unless we have to be safe
    883           if (!safe) {
    884 #ifdef TESTING
    885               if (x == TEST_X && y == TEST_Y) {
    886                   fprintf(stderr, "Single input to combine, safety off.\n");
    887               }
    888 #endif
    889               imageValue = pixelData->data.F32[0];
    890               if (variance) {
    891                   varianceValue = pixelVariances->data.F32[0];
    892               }
    893               maskValue = 0;
    894           } else {
    895               if (!rejection) {
    896                   combineMarkReject(inputs, x, y, pixelSources->data.U16[0]);
    897               }
    898 #ifdef TESTING
    899               numRejected = 1;
    900               if (x == TEST_X && y == TEST_Y) {
    901                   fprintf(stderr, "Single input to combine, safety on, pixel is bad.\n");
    902               }
    903 #endif
    904           }
    905           break;
    906       }
    907       case 2: {
    908           // Accept the mean of the pixels only if we're going to reject based on the variance, or we're not
    909           // playing safe
    910           if (useVariance || !safe) {
    911               float mean, var;   // Mean and variance from combination
    912               if (combinationMeanVariance(&mean, &var, pixelData, pixelVariances, pixelWeights)) {
    913                   imageValue = mean;
    914                   varianceValue = var;
    915                   maskValue = 0;
    916 #ifdef TESTING
    917                   if (x == TEST_X && y == TEST_Y) {
    918                       fprintf(stderr, "Two inputs to combine using variance/unsafe --> %f %f\n",
    919                               mean, var);
    920                   }
    921 #endif
    922               }
    923           }
    924           if (useVariance && numIter > 0) {
    925               // Use variance to check that the two are consistent
    926               float diff = 0.5 * (pixelData->data.F32[0] - pixelData->data.F32[1]); // Mean flux difference
    927               float var1 = pixelVariances->data.F32[0]; // Variance of first
    928               float var2 = pixelVariances->data.F32[1]; // Variance of second
    929               // Systematic error contributes to the rejection level
    930               var1 += PS_SQR(sys * pixelData->data.F32[0]);
    931               var2 += PS_SQR(sys * pixelData->data.F32[1]);
    932 
    933               float sigma2 = var1 + var2; // Combined variance
    934               if (PS_SQR(diff) > PS_SQR(rej) * sigma2) {
    935                   // Not consistent: don't believe either!
    936                   if (rejection) {
    937                       imageValue = NAN;
    938                       varianceValue = NAN;
    939                       maskValue = bad;
    940                   } else if (allowSuspect && numSuspect > 0) {
    941                       combineMarkReject(inputs, x, y, pixelSources->data.U16[0]);
    942                       combineMarkReject(inputs, x, y, pixelSources->data.U16[1]);
    943                   } else {
    944                       combineMarkInspect(inputs, x, y, pixelSources->data.U16[0]);
    945                       combineMarkInspect(inputs, x, y, pixelSources->data.U16[1]);
    946                   }
    947 #ifdef TESTING
    948                   if (x == TEST_X && y == TEST_Y) {
    949                       fprintf(stderr, "Both pixels marked for inspection (%f > %f x %f\n)",
    950                               diff, rej, sqrtf(sigma2));
    951                   }
    952                   numRejected = 2;
    953 #endif
    954               }
    955           }
    956           break;
    957       }
    958       case 3: {
    959           // Can combine without too much worrying, but want to be a bit careful on the rejection
    960           float mean, var;           // Mean and variance of the combination
    961           if (!combinationMeanVariance(&mean, &var, pixelData, pixelVariances, pixelWeights)) {
    962               break;
    963           }
    964           imageValue = mean;
    965           varianceValue = var;
    966           maskValue = 0;
    967 #ifdef TESTING
    968           if (x == TEST_X && y == TEST_Y) {
    969               fprintf(stderr, "Combined inputs: %f %f\n", mean, var);
    970           }
    971 #endif
    972 
    973           if (numIter > 0) {
    974               if (!useVariance) {
    975                   combineRejectionGeneral();
    976               } else {
    977                   // Differences between pixel values
    978                   float diff01 = pixelData->data.F32[0] - pixelData->data.F32[1];
    979                   float diff12 = pixelData->data.F32[1] - pixelData->data.F32[2];
    980                   float diff20 = pixelData->data.F32[2] - pixelData->data.F32[0];
    981                   // Variance for each pixel
    982                   float var0 = pixelVariances->data.F32[0] + PS_SQR(sys * pixelData->data.F32[0]);
    983                   float var1 = pixelVariances->data.F32[1] + PS_SQR(sys * pixelData->data.F32[1]);
    984                   float var2 = pixelVariances->data.F32[2] + PS_SQR(sys * pixelData->data.F32[2]);
    985                   float rej2 = PS_SQR(rej); // Rejection level squared
    986                   // Errors in pixel differences
    987                   float err01 = var0 + var1;
    988                   float err12 = var1 + var2;
    989                   float err20 = var2 + var0;
    990 
    991                   int badPairs = 0;         // Number of bad pairs
    992                   bool bad01 = false, bad12 = false, bad20 = false; // Pair is bad?
    993                   if (PS_SQR(diff01) > rej2 * err01) {
    994                       bad01 = true;
    995                       badPairs++;
    996                   }
    997                   if (PS_SQR(diff12) > rej2 * err12) {
    998                       bad12 = true;
    999                       badPairs++;
    1000                   }
    1001                   if (PS_SQR(diff20) > rej2 * err20) {
    1002                       bad20 = true;
    1003                       badPairs++;
    1004                   }
    1005 
    1006                   switch (badPairs) {
    1007                     case 0:
    1008                       // Nothing to worry about!
    1009                       break;
    1010                     case 1:
    1011                       // Can't tell which image is bad, so be sure to get it
    1012                       if (bad01) {
    1013                           combineInspect(inputs, x, y, pixelSources->data.U16[0]);
    1014                           combineInspect(inputs, x, y, pixelSources->data.U16[1]);
    1015                           break;
    1016                       }
    1017                       if (bad12) {
    1018                           combineInspect(inputs, x, y, pixelSources->data.U16[1]);
    1019                           combineInspect(inputs, x, y, pixelSources->data.U16[2]);
    1020                           break;
    1021                       }
    1022                       if (bad20) {
    1023                           combineInspect(inputs, x, y, pixelSources->data.U16[2]);
    1024                           combineInspect(inputs, x, y, pixelSources->data.U16[0]);
    1025                           break;
    1026                       }
    1027                       psAbort("Should never get here");
    1028                     case 2:
    1029                       if (bad01 && bad12) {
    1030                           // 2 and 0 are good
    1031                           combineInspect(inputs, x, y, pixelSources->data.U16[1]);
    1032                           break;
    1033                       }
    1034                       if (bad12 && bad20) {
    1035                           // 0 and 1 are good
    1036                           combineInspect(inputs, x, y, pixelSources->data.U16[2]);
    1037                           break;
    1038                       }
    1039                       if (bad20 && bad01) {
    1040                           // 1 and 2 are good
    1041                           combineInspect(inputs, x, y, pixelSources->data.U16[0]);
    1042                           break;
    1043                       }
    1044                       psAbort("Should never get here");
    1045                     case 3:
    1046                       // Everything's bad
    1047                       combineInspect(inputs, x, y, pixelSources->data.U16[0]);
    1048                       combineInspect(inputs, x, y, pixelSources->data.U16[1]);
    1049                       combineInspect(inputs, x, y, pixelSources->data.U16[2]);
    1050                       break;
    1051                   }
    1052               }
    1053           }
    1054           break;
    1055       }
    1056       default: {
    1057           // Record the value derived with no clipping, because pixels rejected using the harsh clipping
    1058           // applied in the first pass might later be accepted.
    1059           float mean, var;           // Mean and variance of the combination
    1060           if (!combinationMeanVariance(&mean, &var, pixelData, pixelVariances, pixelWeights)) {
    1061               break;
    1062           }
    1063           imageValue = mean;
    1064           varianceValue = var;
    1065           maskValue = 0;
    1066 #ifdef TESTING
    1067           if (x == TEST_X && y == TEST_Y) {
    1068               fprintf(stderr, "Combined inputs: %f %f\n", mean, var);
    1069           }
    1070 #endif
    1071 
    1072           // Prepare for clipping iteration
    1073 
    1074           break;
    1075       }
    1076     }
    1077 
    1078     return;
    1079 }
    1080 #endif
    1081795
    1082796// Ensure the input array of pmStackData is valid, and get some details out of it
    … …  
    1235949/// Stack input images
    1236950bool pmStackCombine(pmReadout *combined, psArray *input, psImageMaskType maskVal, psImageMaskType maskSuspect,
    1237                     psImageMaskType bad, int kernelSize, int numIter, float rej, float sys, float olympic,
     951                    psImageMaskType bad, int kernelSize, float iter, float rej, float sys, float olympic,
    1238952                    bool useVariance, bool safe, bool rejection)
    1239953{
    … …  
    1246960    PS_ASSERT_INT_NONNEGATIVE(kernelSize, false);
    1247961    PS_ASSERT_INT_POSITIVE(bad, false);
    1248     PS_ASSERT_INT_NONNEGATIVE(numIter, false);
    1249962    if (isnan(rej)) {
    1250         PS_ASSERT_INT_EQUAL(numIter, 0, false);
     963        PS_ASSERT_FLOAT_EQUAL(iter, 0, false);
    1251964    } else {
     965        PS_ASSERT_FLOAT_LARGER_THAN(iter, 0, false);
    1252966        PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false);
    1253967    }
    … …  
    13401054#ifdef TESTING
    13411055            if (x == TEST_X && y == TEST_Y) {
    1342                 fprintf(stderr, "Combining pixel %d,%d: %x %x %d %f %f %f %d %d %d\n",
    1343                         x, y, maskVal, bad, numIter, rej, sys, olympic, useVariance, safe, rejection);
     1056                fprintf(stderr, "Combining pixel %d,%d: %x %x %f %f %f %f %d %d %d\n",
     1057                        x, y, maskVal, bad, iter, rej, sys, olympic, useVariance, safe, rejection);
    13441058            }
    13451059#endif
    … …  
    13651079            bool suspect;               // Suspect pixels in stack?
    13661080            combineExtract(&num, &suspect, buffer, combinedImage, combinedMask, combinedVariance,
    1367                            input, weights, addVariance, reject, x, y, maskVal, maskSuspect, true);
    1368             if (numIter == 0) {
    1369                 combinePixels(combinedImage, combinedMask, combinedVariance,
    1370                               num, buffer, x, y, bad, safe);
    1371             } else {
    1372                 if (combineTest(num, suspect, input, buffer, x, y, numIter, rej, sys, olympic,
    1373                                 useVariance, safe, true)) {
    1374                     // Need to repeat without suspect pixels
    1375                     combineExtract(&num, &suspect, buffer, combinedImage, combinedMask, combinedVariance,
    1376                            input, weights, addVariance, reject, x, y, maskVal, maskSuspect, false);
    1377                     combineTest(num, suspect, input, buffer, x, y, numIter, rej, sys, olympic,
    1378                                 useVariance, safe, false);
    1379                 }
     1081                           input, weights, addVariance, reject, x, y, maskVal, maskSuspect);
     1082            combinePixels(combinedImage, combinedMask, combinedVariance, num, buffer, x, y, bad, safe);
     1083
     1084            if (iter > 0) {
     1085                combineTest(num, suspect, input, buffer, x, y, iter, rej, sys, olympic,
     1086                            useVariance, safe);
    13801087            }
    13811088        }
    … …  
    13851092    psFree(weights);
    13861093    psFree(buffer);
     1094
     1095
    13871096
    13881097#ifndef PS_NO_TRACE

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

    r25964 r26007  
    4545                    psArray *input,     ///< Input array of pmStackData
    4646                    psImageMaskType maskVal, ///< Mask value of bad pixels
     47                    psImageMaskType suspect, ///< Mask value of suspect pixels
    4748                    psImageMaskType bad,     ///< Mask value to give rejected pixels
    4849                    int kernelSize,     ///< Half-size of the convolution kernel
    49                     int numIter,        ///< Number of iterations
     50                    float iter,         ///< Number of iterations per input
    5051                    float rej,          ///< Rejection limit (standard deviations)
    5152                    float sys,          ///< Relative systematic error

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

    r25917 r26007  
    832832    psTrace("psModules.imcombine", 1, "RMS deviation: %f\n", rms);
    833833
     834    fprintf(stderr, "Mean = %f ; chi^2 = %f for %ld d.o.f.\n",
     835            mean,
     836            mean * PS_SQR(2 * stamps->footprint + 1) * numStamps,
     837            PS_SQR(2 * stamps->footprint + 1) * numStamps - kernels->solution1->n);
     838
     839
    834840    kernels->mean = mean;
    835841    kernels->rms = rms;

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

    r25899 r26007  
    228228                    double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x];
    229229#ifdef USE_VARIANCE
    230                     aa /= variance->kernel[y][x];
    231                     bb /= variance->kernel[y][x];
    232                     ab /= variance->kernel[y][x];
     230                    float var = 1.0 / variance->kernel[y][x];
     231                    aa /= var;
     232                    bb /= var;
     233                    ab /= var;
    233234#endif
    234235                    sumAA += aa;
Note: See TracChangeset for help on using the changeset viewer.