IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 16626


Ignore:
Timestamp:
Feb 22, 2008, 5:00:30 PM (18 years ago)
Author:
Paul Price
Message:

Trying to get better about "variance" vs "weight". The "variance"
comes from the variance map (the unfortunately-named
"readout->weight"), while the "weight" is the relative weighting for
each image.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r16624 r16626  
    88 *  @author GLG, MHPCC
    99 *
    10  *  @version $Revision: 1.21 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2008-02-23 02:33:59 $
     10 *  @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2008-02-23 03:00:30 $
    1212 *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
    1313 *
    … …  
    4343    psVector *pixels;                   // Pixel values
    4444    psVector *masks;                    // Pixel masks
    45     psVector *weights;                  // Pixel weights
     45    psVector *variances;                // Pixel variances
     46    psVector *weights;                  // Pixel weightings
    4647    psVector *sources;                  // Pixel sources (which image did they come from?)
    4748    psVector *sort;                     // Buffer for sorting (to get a robust estimator of the standard dev)
    … …  
    5253    psFree(buffer->pixels);
    5354    psFree(buffer->masks);
     55    psFree(buffer->variances);
    5456    psFree(buffer->weights);
    5557    psFree(buffer->sources);
    … …  
    6769    buffer->pixels = psVectorAlloc(numImages, PS_TYPE_F32);
    6870    buffer->masks = psVectorAlloc(numImages, PS_TYPE_MASK);
     71    buffer->variances = psVectorAlloc(numImages, PS_TYPE_F32);
    6972    buffer->weights = psVectorAlloc(numImages, PS_TYPE_F32);
    7073    buffer->sources = psVectorAlloc(numImages, PS_TYPE_U16);
    … …  
    8487
    8588// Generate a mean value for the combination
    86 // Not using psVectorStats because it assumes that the weights are errors, and weights by 1/error^2
     89// Not using psVectorStats because it assumes that the "weights" are errors, and weights by 1/error^2
    8790static bool combinationMean(float *mean, // Mean value, to return
    8891                            const psVector *values, // Values to combine
    … …  
    216219}
    217220
    218 // Common path for a bad combined pixel
    219 static inline bool combineBad(psImage *image, // Combined image
    220                               psImage *mask, // Combined mask
    221                               psImage *weight, // Combined variance map
    222                               int x, int y, // Coordinates of interest
    223                               psMaskType bad // Value for bad pixels
    224                               )
    225 {
    226     image->data.F32[y][x] = NAN;
    227     if (weight) {
    228         weight->data.F32[y][x] = NAN;
    229     }
    230     mask->data.PS_TYPE_MASK_DATA[y][x] = bad;
    231     return false;
    232 }
    233221
    234222// Mark a pixel for inspection
    … …  
    250238static void combinePixels(psImage *image, // Combined image, for output
    251239                          psImage *mask, // Combined mask, for output
    252                           psImage *weight, // Combined variance map, for output
     240                          psImage *variance, // Combined variance map, for output
    253241                          const psArray *inputs, // Stack data
    254242                          const psVector *weights, // Global (single value) weights for data, or NULL
    … …  
    271259    assert(rej > 0);
    272260    assert(buffer);
    273     assert((useVariance && weight) || !useVariance);
     261    assert((useVariance && variance) || !useVariance);
    274262
    275263    psVector *pixelData = buffer->pixels; // Values for the pixel of interest
    276264    psVector *pixelMasks = buffer->masks; // Masks for the pixel of interest
    277     psVector *pixelWeights = buffer->weights; // Weights for the pixel of interest
     265    psVector *pixelVariances = buffer->variances; // Variances for the pixel of interest
     266    psVector *pixelWeights = buffer->weights; // Image weights for the pixel of interest
    278267    psVector *pixelSources = buffer->sources; // Sources for the pixel of interest
    279268    psVector *sort = buffer->sort;      // Sort buffer
    … …  
    301290
    302291        psImage *image = data->readout->image; // Image of interest
    303         psImage *weight = data->readout->weight; // Weight map of interest
     292        psImage *variance = data->readout->weight; // Variance ("weight") map of interest
    304293        pixelData->data.F32[num] = image->data.F32[yIn][xIn];
    305         if (weight) {
    306             pixelWeights->data.F32[num] = weight->data.F32[yIn][xIn];
    307         }
     294        if (variance) {
     295            pixelVariances->data.F32[num] = variance->data.F32[yIn][xIn];
     296        }
     297        pixelWeights->data.F32[num] = data->weight;
    308298        pixelSources->data.U16[num] = i;
    309299        num++;
    310300    }
    311301    pixelData->n = num;
    312     if (weight) {
    313         pixelWeights->n = num;
    314     }
     302    if (variance) {
     303        pixelVariances->n = num;
     304    }
     305    pixelWeights->n = num;
    315306    pixelSources->n = num;
     307
    316308
    317309    // The sensible thing to do varies according to how many good pixels there are.
    … …  
    327319          if (!safe) {
    328320              imageValue = pixelData->data.F32[0];
    329               if (weight) {
    330                   varianceValue = pixelWeights->data.F32[0];
     321              if (variance) {
     322                  varianceValue = pixelVariances->data.F32[0];
    331323              }
    332324              maskValue = 0;
    … …  
    338330          // playing safe
    339331          if (useVariance || !safe) {
    340               float mean, variance;   // Mean and variance from combination
    341               if (combinationMean(&mean, pixelData, weights, pixelMasks, maskVal) &&
    342                   (!weight || combinationVariance(&variance, pixelWeights, weights, pixelMasks, maskVal))) {
     332              float mean, var;   // Mean and variance from combination
     333              if (combinationMean(&mean, pixelData, pixelWeights, NULL, maskVal) &&
     334                  (!variance || combinationVariance(&var, pixelVariances, pixelWeights, NULL, maskVal))) {
    343335                  imageValue = mean;
    344                   varianceValue = variance;
     336                  varianceValue = var;
    345337                  maskValue = 0;
    346338              }
    … …  
    349341              // Use variance to check that the two are consistent
    350342              float diff = pixelData->data.F32[0] - pixelData->data.F32[1];
    351               float sigma = sqrtf(pixelWeights->data.F32[0] + pixelWeights->data.F32[1]);
     343              float sigma = sqrtf(pixelVariances->data.F32[0] + pixelVariances->data.F32[1]);
    352344              if (fabs(diff) > rej * sigma) {
    353345                  // Not consistent: mark both for inspection
    … …  
    361353          // Record the value derived with no clipping, because pixels rejected using the harsh clipping
    362354          // applied in the first pass might later be accepted.
    363           float mean, variance;           // Mean and variance of the combination
    364           if (!combinationMean(&mean, pixelData, weights, pixelMasks, maskVal) ||
    365               (weight && !combinationVariance(&variance, pixelWeights, weights, pixelMasks, maskVal))) {
     355          float mean, var;           // Mean and variance of the combination
     356          if (!combinationMean(&mean, pixelData, pixelWeights, pixelMasks, maskVal) ||
     357              (variance && !combinationVariance(&var, pixelVariances, pixelWeights, pixelMasks, maskVal))) {
    366358              break;
    367359          }
    368360          imageValue = mean;
    369           varianceValue = variance;
     361          varianceValue = var;
    370362          maskValue = 0;
    371363
    … …  
    377369                  // Convert to rejection limits
    378370                  for (int i = 0; i < num; i++) {
    379                       pixelWeights->data.F32[i] = rej * sqrtf(pixelWeights->data.F32[i]);
     371                      pixelVariances->data.F32[i] = rej * sqrtf(pixelVariances->data.F32[i]);
    380372                  }
    381373              }
    … …  
    406398                  }
    407399                  float diff = fabsf(pixelData->data.F32[j] - median); // Difference from expected
    408                   if (diff > (useVariance ? pixelWeights->data.F32[i] : limit)) {
     400                  if (diff > (useVariance ? pixelVariances->data.F32[i] : limit)) {
    409401                      pixelMasks->data.PS_TYPE_MASK_DATA[j] = 0xff;
    410402                      combineInspect(inputs, x, y, pixelSources->data.U16[j]);
    … …  
    423415
    424416// Ensure the input array of pmStackData is valid, and get some details out of it
    425 static bool validateInputData(bool *haveWeights, // Do we have weights in the sky images?
     417static bool validateInputData(bool *haveVariances, // Do we have variance maps in the sky images?
    426418                              bool *havePixels, // Do we have lists of pixels?
    427419                              int *num,    // Number of inputs
    … …  
    439431    PS_ASSERT_PTR_NON_NULL(data, false);
    440432    assert(psMemGetDeallocator(data) == (psFreeFunc)stackDataFree); // Ensure it's the right type
    441     *haveWeights = false;
     433    *haveVariances = false;
    442434    PS_ASSERT_IMAGE_NON_NULL(data->readout->image, false);
    443435    PS_ASSERT_IMAGE_TYPE(data->readout->image, PS_TYPE_F32, false);
    … …  
    448440    *numRows = data->readout->image->numRows;
    449441    if (data->readout->weight) {
    450         *haveWeights = true;
     442        *haveVariances = true;
    451443        PS_ASSERT_IMAGE_NON_NULL(data->readout->weight, false);
    452444        PS_ASSERT_IMAGES_SIZE_EQUAL(data->readout->image, data->readout->weight, false);
    … …  
    476468        PS_ASSERT_IMAGE_SIZE(data->readout->image, *numCols, *numRows, false);
    477469        PS_ASSERT_IMAGES_SIZE_EQUAL(data->readout->image, data->readout->mask, false);
    478         if (*haveWeights) {
     470        if (*haveVariances) {
    479471            PS_ASSERT_IMAGE_NON_NULL(data->readout->weight, false);
    480472            PS_ASSERT_IMAGES_SIZE_EQUAL(data->readout->image, data->readout->weight, false);
    … …  
    651643{
    652644    PS_ASSERT_PTR_NON_NULL(combined, false);
    653     bool haveWeights;                   // Do we have the weight maps?
     645    bool haveVariances;                 // Do we have the variance maps?
    654646    bool havePixels;                    // Do we have lists of pixels?
    655647    int num;                            // Number of inputs
    656648    int numCols, numRows;               // Size of (sky) images
    657     if (!validateInputData(&haveWeights, &havePixels, &num, &numCols, &numRows, input)) {
     649    if (!validateInputData(&haveVariances, &havePixels, &num, &numCols, &numRows, input)) {
    658650        return false;
    659651    }
    … …  
    669661        PS_ASSERT_IMAGES_SIZE_EQUAL(combined->image, combined->mask, false);
    670662    }
    671     if (useVariance && !haveWeights) {
     663    if (useVariance && !haveVariances) {
    672664        psWarning("Unable to use variance in rejection if no variance maps supplied --- option turned off");
    673665        useVariance = false;
    674666    }
    675667
    676     // Pull out the weightings
     668    // Pull out the image weightings
    677669    psVector *weights = psVectorAlloc(num, PS_TYPE_F32);
    678670    for (int i = 0; i < num; i++) {
    … …  
    712704        psImage *combinedImage = combined->image; // Combined image
    713705        psImage *combinedMask = combined->mask; // Combined mask
    714         psImage *combinedWeight = combined->weight; // Combined mask
     706        psImage *combinedVariance = combined->weight; // Combined variance map
    715707
    716708        psArray *pixelMap = pixelMapGenerate(input, maxInputCols, maxInputRows); // Map of pixels to source
    … …  
    731723            }
    732724            psVector *reject = pixelMapQuery(pixelMap, x, y); // Inspect these images closely
    733             combinePixels(combinedImage, combinedMask, combinedWeight, input, weights, reject, x, y,
     725            combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, reject, x, y,
    734726                          maskVal, bad, numIter, rej, useVariance, safe, buffer);
    735727        }
    … …  
    750742        }
    751743
    752         psImage *combinedWeights = combined->weight; // Combined weight map
    753         if (haveWeights && !combinedWeights) {
     744        psImage *combinedVariance = combined->weight; // Combined variance map
     745        if (haveVariances && !combinedVariance) {
    754746            combined->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    755             combinedWeights = combined->weight;
     747            combinedVariance = combined->weight;
    756748        }
    757749
    … …  
    769761        for (int y = minInputRows; y < maxInputRows; y++) {
    770762            for (int x = minInputCols; x < maxInputCols; x++) {
    771                 combinePixels(combinedImage, combinedMask, combinedWeights, input, weights, NULL, x, y,
     763                combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, NULL, x, y,
    772764                              maskVal, bad, numIter, rej, useVariance, safe, buffer);
    773765            }
Note: See TracChangeset for help on using the changeset viewer.