IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Feb 13, 2009, 1:52:14 PM (17 years ago)
Author:
Paul Price
Message:

Using covariance factor to set the correct variance level. Median does not do a good job of selecting the middle of the distribution when there's a variety of input weights (esp. when one image has dominant weight), so now using an 'olympic weighted mean' to provide the mean value for rejection. Tried a 'weighted median', but it doesn't do as good a job.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/imcombine/pmStack.c

    r21363 r21476  
    88 *  @author GLG, MHPCC
    99 *
    10  *  @version $Revision: 1.47 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2009-02-06 02:31:25 $
     10 *  @version $Revision: 1.48 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2009-02-13 23:52:14 $
    1212 *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
    1313 *
    … …  
    3030#define PIXEL_LIST_BUFFER 100           // Number of entries to add to pixel list at a time
    3131#define PIXEL_MAP_BUFFER 2              // Number of entries to add to pixel map at a time
    32 #define VARIANCE_FACTORS                // Use variance factors when calculating the variances?
    33 //#define ADD_VARIANCE                  // Allow additional variance (besides variance factor)?
     32#define ADD_VARIANCE                    // Allow additional variance (besides variance factor)?
    3433#define NUM_DIRECT_STDEV 5              // For less than this number of values, measure stdev directly
     34
     35//#define TESTING                         // Enable test output
     36//#define TEST_X 2318                     // x coordinate to examine
     37//#define TEST_Y 2306                     // y coordinate to examine
    3538
    3639
    … …  
    190193}
    191194
     195#if 0
     196// Return the weighted median for the pixels
     197// This does not appear to produce as clean images as the weighted Olympic mean
     198static float combinationWeightedMedian(const psVector *values, // Values to combine
     199                                       const psVector *weights, // Weights to combine
     200                                       const psVector *masks, // Mask to apply
     201                                       psVector *sortBuffer // Buffer for sorting
     202                                       )
     203{
     204    double sumWeight = 0.0;             // Sum of weights
     205    for (int j = 0; j < values->n; j++) {
     206        if (masks->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
     207            continue;
     208        }
     209        sumWeight += weights->data.F32[j];
     210    }
     211
     212    sortBuffer = psVectorSortIndex(sortBuffer, values);
     213    double target = sumWeight / 2.0;    // Target weight
     214
     215    int dominant = -1;                  // Index of dominant value, if any
     216    double cumulativeWeight = 0.0;      // Sum of weights
     217    for (int j = 0; j < values->n; j++) {
     218        if (masks->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
     219            continue;
     220        }
     221        int index = sortBuffer->data.S32[j];  // Index of value of interest
     222        float weight = weights->data.F32[index]; // Weight for value of interest
     223        if (weight >= target) {
     224            // Get the weighted median of the rest
     225            dominant = index;
     226            sumWeight -= weight;
     227            target = sumWeight / 2.0;
     228            continue;
     229        }
     230        cumulativeWeight += weight;
     231        if (cumulativeWeight >= target) {
     232            float median = values->data.F32[index]; // Weighted median median
     233            if (dominant != -1) {
     234                // In the case that a single value contains a disproportionate weight compared to the rest,
     235                // we use a weighted mean between that dominant value and the weighted median of the rest.
     236                return (values->data.F32[dominant] * weights->data.F32[dominant] + median * sumWeight) /
     237                    (weights->data.F32[dominant] + sumWeight);
     238            }
     239            return median;
     240        }
     241    }
     242
     243    return NAN;
     244}
     245#endif
     246
     247// Return the weighted Olympic mean for the pixels
     248static float combinationWeightedOlympic(const psVector *values, // Values to combine
     249                                        const psVector *weights, // Weights to combine
     250                                        const psVector *masks, // Mask to apply
     251                                        float frac, // Fraction to discard
     252                                        psVector *sortBuffer // Buffer for sorting
     253                                        )
     254{
     255    int numGood = 0;                    // Number of good values
     256    for (int i = 0; i < values->n; i++) {
     257        if (masks->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
     258            continue;
     259        }
     260        numGood++;
     261    }
     262
     263    int numBad = frac * numGood + 0.5;  // Number of bad values
     264    int low = numBad / 2, high = low + numGood; // Indices (modulo masked pixels) delimiting range of interest
     265
     266    sortBuffer = psVectorSortIndex(sortBuffer, values);
     267
     268    double sumValues = 0.0, sumWeight = 0.0; // Sums for weighted mean
     269    for (int i = 0, j = 0; i < values->n; i++) {
     270        int index = sortBuffer->data.S32[i]; // Index of interest
     271        if (masks->data.PS_TYPE_VECTOR_MASK_DATA[index]) {
     272            continue;
     273        }
     274        j++;
     275        if (j > high) {
     276            break;
     277        }
     278        if (j <= low) {
     279            continue;
     280        }
     281        sumValues += values->data.F32[index] * weights->data.F32[index];
     282        sumWeight += weights->data.F32[index];
     283    }
     284
     285    return sumValues / sumWeight;
     286}
    192287
    193288// Mark a pixel for inspection
    … …  
    213308                          const psArray *inputs, // Stack data
    214309                          const psVector *weights, // Global (single value) weights for data, or NULL
    215                           const psVector *varFactors, // Variance factors for data
     310                          const psVector *addVariance, // Additional variance for data
    216311                          const psVector *reject, // Indices of pixels to reject, or NULL
    217312                          int x, int y, // Coordinates of interest; frame of output image
    … …  
    221316                          float rej, // Number of standard deviations at which to reject
    222317                          float sys,    // Relative systematic error
     318                          float discard,// Fraction of values to discard (Olympic weighted mean)
    223319                          bool useVariance, // Use variance for rejection when combining?
    224320                          bool safe,    // Combine safely?
    … …  
    232328    assert(numIter >= 0);
    233329    assert(buffer);
    234     assert(varFactors);
     330    assert(addVariance);
    235331    assert((useVariance && variance) || !useVariance);
    236332
    … …  
    267363        pixelData->data.F32[num] = image->data.F32[yIn][xIn];
    268364        if (variance) {
    269             pixelVariances->data.F32[num] = variance->data.F32[yIn][xIn];
     365            pixelVariances->data.F32[num] = variance->data.F32[yIn][xIn] * addVariance->data.F32[i];
    270366        }
    271367        pixelWeights->data.F32[num] = data->weight;
    … …  
    280376    pixelSources->n = num;
    281377
     378#ifdef TESTING
     379    if (x == TEST_X && y == TEST_Y) {
     380        for (int i = 0; i < num; i++) {
     381            fprintf(stderr, "Input %d (%" PRIu16 "): %f %f %f\n",
     382                    i, pixelSources->data.U16[i], pixelData->data.F32[i], pixelVariances->data.F32[i],
     383                    pixelWeights->data.F32[i]);
     384        }
     385    }
     386#endif
    282387
    283388    // The sensible thing to do varies according to how many good pixels there are.
    … …  
    288393      case 0:
    289394        // Nothing to combine: it's bad
     395#ifdef TESTING
     396    if (x == TEST_X && y == TEST_Y) {
     397        fprintf(stderr, "No inputs to combine, pixel is bad.\n");
     398    }
     399#endif
    290400        break;
    291401      case 1: {
    292402          // Accept the single pixel unless we have to be safe
    293403          if (!safe) {
     404#ifdef TESTING
     405              if (x == TEST_X && y == TEST_Y) {
     406                  fprintf(stderr, "Single input to combine, safety off.\n");
     407              }
     408#endif
    294409              imageValue = pixelData->data.F32[0];
    295410              if (variance) {
    … …  
    298413              maskValue = 0;
    299414          }
     415#ifdef TESTING
     416          else if (x == TEST_X && y == TEST_Y) {
     417              fprintf(stderr, "Single input to combine, safety on, pixel is bad.\n");
     418          }
     419#endif
    300420          break;
    301421      }
    … …  
    309429                  varianceValue = var;
    310430                  maskValue = 0;
     431#ifdef TESTING
     432                  if (x == TEST_X && y == TEST_Y) {
     433                      fprintf(stderr, "Two inputs to combine using variance/unsafe --> %f %f\n",
     434                              mean, var);
     435                  }
     436#endif
    311437              }
    312438          }
    … …  
    316442              float var1 = pixelVariances->data.F32[0]; // Variance of first
    317443              float var2 = pixelVariances->data.F32[1]; // Variance of second
    318 #ifdef VARIANCE_FACTORS
    319               // Variance factor contributes to the rejection level
    320               var1 *= varFactors->data.F32[pixelSources->data.U16[0]];
    321               var2 *= varFactors->data.F32[pixelSources->data.U16[1]];
    322 #endif
    323444              // Systematic error contributes to the rejection level
    324445              var1 += PS_SQR(sys * pixelData->data.F32[0]);
    325446              var2 += PS_SQR(sys * pixelData->data.F32[1]);
    326447
    327               float sigma2 = var1 + var2; //
     448              float sigma2 = var1 + var2; // Combined variance
    328449              if (PS_SQR(diff) > PS_SQR(rej) * sigma2) {
    329450                  // Not consistent: mark both for inspection
    330451                  combineInspect(inputs, x, y, pixelSources->data.U16[0]);
    331452                  combineInspect(inputs, x, y, pixelSources->data.U16[1]);
     453#ifdef TESTING
     454                  if (x == TEST_X && y == TEST_Y) {
     455                      fprintf(stderr, "Both pixels marked for inspection (%f > %f x %f\n)",
     456                              diff, rej, sqrtf(sigma2));
     457                  }
     458#endif
    332459              }
    333460          }
    … …  
    344471          varianceValue = var;
    345472          maskValue = 0;
     473#ifdef TESTING
     474          if (x == TEST_X && y == TEST_Y) {
     475              fprintf(stderr, "Combined inputs: %f %f\n", mean, var);
     476          }
     477#endif
    346478
    347479          // Prepare for clipping iteration
    … …  
    353485                  // Using squared rejection limit because it's cheaper than sqrts
    354486                  float rej2 = PS_SQR(rej); // Rejection level squared
     487                  double sumWeights = 0.0;
     488                  for (int i = 0; i < num; i++) {
     489                      sumWeights += pixelWeights->data.F32[i];
     490                  }
    355491                  for (int i = 0; i < num; i++) {
    356492                      // Systematic error contributes to the rejection level
    357                       pixelVariances->data.F32[i] += PS_SQR(sys * pixelData->data.F32[i]);
    358 #ifdef VARIANCE_FACTORS
    359                       // Variance factor contributes to the rejection level
    360                       pixelVariances->data.F32[i] *= varFactors->data.F32[pixelSources->data.U16[i]];
    361 #endif
    362                       pixelVariances->data.F32[i] *= rej2;
     493                      float sysVar = PS_SQR(sys * pixelData->data.F32[i]);
     494#ifdef TESTING
     495                      // Correct variance for comparison against weighted mean including itself
     496                      float compare = 1.0 - pixelWeights->data.F32[i] / sumWeights;
     497                      if (x == TEST_X && y == TEST_Y) {
     498                          fprintf(stderr, "Variance %d (%d): %f %f %f\n", i, pixelSources->data.U16[i],
     499                                  pixelVariances->data.F32[i], sysVar, compare);
     500                      }
     501#endif
     502
     503                      pixelVariances->data.F32[i] = rej2 * (pixelVariances->data.F32[i] + sysVar);
    363504                  }
    364505              }
    … …  
    371512          for (int i = 0; i < numIter && numClipped > 0 && num - totalClipped > 2; i++) {
    372513              numClipped = 0;
    373               float median, stdev;    // Median and stdev of the combination, for rejection
    374 
    375               if (!combinationMedianStdev(&median, useVariance ? NULL : &stdev,
    376                                           pixelData, pixelMasks, sort)) {
    377                   psWarning("Bad median/stdev at %d,%d", x, y);
    378                   break;
     514              float median = NAN;       // Middle of distribution
     515              float limit = NAN;        // Rejection limit
     516              if (!useVariance) {
     517                  float stdev;  // Median and stdev of the combination, for rejection
     518                  if (!combinationMedianStdev(&median, useVariance ? NULL : &stdev,
     519                                              pixelData, pixelMasks, sort)) {
     520                      psWarning("Bad median/stdev at %d,%d", x, y);
     521                      break;
     522                  }
     523                  limit = rej * stdev;
     524#ifdef TESTING
     525                  if (x == TEST_X && y == TEST_Y) {
     526                      fprintf(stderr, "Rejection limit: %f\n", limit);
     527                  }
     528#endif
     529              } else {
     530                  median = combinationWeightedOlympic(pixelData, pixelWeights, pixelMasks, discard, sort);
    379531              }
    380532
    381               float limit = rej * stdev; // Rejection limit, when rejecting based on data properties
     533#ifdef TESTING
     534              if (x == TEST_X && y == TEST_Y) {
     535                  fprintf(stderr, "Median: %f\n", median);
     536              }
     537#endif
     538
    382539
    383540// Mask a pixel for inspection
    … …  
    395552                  if (useVariance) {
    396553                      // Comparing squares --- cheaper than lots of sqrts
    397                       // pixelVariances includes the variance factor and the rejection limit, from above
     554                      // pixelVariances includes the rejection limit, from above
    398555                      if (PS_SQR(diff) > pixelVariances->data.F32[j]) {
    399556                          MASK_PIXEL_FOR_INSPECTION();
     557#ifdef TESTING
     558                          if (x == TEST_X && y == TEST_Y) {
     559                              fprintf(stderr, "Rejecting input %d based on variance: %f > %f\n",
     560                                      j, diff, sqrtf(pixelVariances->data.F32[j]));
     561                          }
     562#endif
    400563                      }
    401564                  } else if (fabsf(diff) > limit) {
    402565                      MASK_PIXEL_FOR_INSPECTION();
     566#ifdef TESTING
     567                      if (x == TEST_X && y == TEST_Y) {
     568                          fprintf(stderr, "Rejecting input %d based on distribution: %f > %f\n",
     569                                  j, diff, limit);
     570                      }
     571#endif
    403572                  }
    404573              }
    … …  
    564733/// Stack input images
    565734bool pmStackCombine(pmReadout *combined, psArray *input, psImageMaskType maskVal, psImageMaskType bad,
    566                     int kernelSize, int numIter, float rej, float sys,
     735                    int kernelSize, int numIter, float rej, float sys, float discard,
    567736                    bool entire, bool useVariance, bool safe)
    568737{
    … …  
    596765    }
    597766
    598     psVector *varFactors = psVectorAlloc(num, PS_TYPE_F32); // Variance factors for each image
     767    psVector *addVariance = psVectorAlloc(num, PS_TYPE_F32); // Additional variance for each image
    599768    psVector *weights = psVectorAlloc(num, PS_TYPE_F32); // Relative weighting for each image
    600769    psArray *stack = psArrayAlloc(num); // Stack of readouts
    … …  
    606775        }
    607776        weights->data.F32[i] = data->weight;
    608         stack->data[i] = psMemIncrRefCounter(data->readout);
    609         varFactors->data.F32[i] = data->readout->covariance ?
    610             psImageCovarianceFactor(data->readout->covariance) : 1.0;
    611 #if 0
     777        pmReadout *ro = data->readout;  // Readout of interest
     778        stack->data[i] = psMemIncrRefCounter(ro);
     779        addVariance->data.F32[i] = ro->covariance ? psImageCovarianceFactor(ro->covariance) : 1.0;
     780#ifdef ADD_VARIANCE
    612781        if (isfinite(data->addVariance)) {
    613             varFactors->data.F32[i] *= data->addVariance;
     782            addVariance->data.F32[i] *= data->addVariance;
    614783        }
    615784#endif
    … …  
    667836                    psVector *reject = pixelMapQuery(pixelMap, minInputCols, minInputRows,
    668837                                                     x, y); // Reject these images
    669                     combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, varFactors,
    670                                   reject, x, y, maskVal, bad, numIter, rej, sys, useVariance, safe, buffer);
     838                    combinePixels(combinedImage, combinedMask, combinedVariance, input, weights,
     839                                  addVariance, reject, x, y, maskVal, bad, numIter, rej, sys, discard,
     840                                  useVariance, safe, buffer);
    671841                }
    672842            }
    … …  
    681851                psVector *reject = pixelMapQuery(pixelMap, minInputCols, minInputRows,
    682852                                                 x, y); // Reject these images
    683                 combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, varFactors,
    684                               reject, x, y, maskVal, bad, numIter, rej, sys, useVariance, safe, buffer);
     853                combinePixels(combinedImage, combinedMask, combinedVariance, input, weights,
     854                              addVariance, reject, x, y, maskVal, bad, numIter, rej, sys, discard,
     855                              useVariance, safe, buffer);
    685856            }
    686857        }
    … …  
    708879        for (int y = minInputRows; y < maxInputRows; y++) {
    709880            for (int x = minInputCols; x < maxInputCols; x++) {
    710                 combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, varFactors,
    711                               NULL, x, y, maskVal, bad, numIter, rej, sys, useVariance, safe, buffer);
     881                combinePixels(combinedImage, combinedMask, combinedVariance, input, weights,
     882                              addVariance, NULL, x, y, maskVal, bad, numIter, rej, sys, discard,
     883                              useVariance, safe, buffer);
    712884            }
    713885        }
Note: See TracChangeset for help on using the changeset viewer.