IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 27420


Ignore:
Timestamp:
Mar 23, 2010, 4:38:36 PM (16 years ago)
Author:
bills
Message:

Grabbed this version of pmStack.c from Paul's branch. This builds

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r27417 r27420  
    3535
    3636//#define TESTING                         // Enable test output
    37 //#define TEST_X 843-1                     // x coordinate to examine
    38 //#define TEST_Y 813-1                     // y coordinate to examine
     37//#define TEST_X 2148-1                     // x coordinate to examine
     38//#define TEST_Y 248-1                     // y coordinate to examine
    3939//#define TEST_RADIUS 0                    // Radius to examine
    4040
    … …  
    4646    psVector *variances;                // Pixel variances
    4747    psVector *weights;                  // Pixel weightings
    48     psVector *exps;                     // Pixel exposures
    4948    psVector *sources;                  // Pixel sources (which image did they come from?)
    5049    psVector *limits;                   // Rejection limits
    … …  
    5857    psFree(buffer->variances);
    5958    psFree(buffer->weights);
    60     psFree(buffer->exps);
    6159    psFree(buffer->sources);
    6260    psFree(buffer->limits);
    … …  
    7573    buffer->variances = psVectorAlloc(numImages, PS_TYPE_F32);
    7674    buffer->weights = psVectorAlloc(numImages, PS_TYPE_F32);
    77     buffer->exps = psVectorAlloc(numImages, PS_TYPE_F32);
    7875    buffer->sources = psVectorAlloc(numImages, PS_TYPE_U16);
    7976    buffer->limits = psVectorAlloc(numImages, PS_TYPE_F32);
    … …  
    9895static bool combinationMeanVariance(float *mean, // Mean value, to return
    9996                                    float *var, // Variance value, to return
    100                                     float *exp, // Exposure time, to return
    101                                     float *expWeight,          // Weighted exposure time, to return
    10297                                    const psVector *values, // Values to combine
    10398                                    const psVector *variances, // Pixel variances to combine
    104                                     const psVector *exps,      // Exposure times to combine
    10599                                    const psVector *weights // Weights to apply
    106100                                    )
    … …  
    127121    float sumVarianceWeight = 0.0;     // Sum of the pixel variances multiplied by the global weights
    128122    float sumWeight = 0.0;              // Sum of the image weights
    129     float sumExp = 0.0;                 // Sum of the exposure time
    130     float sumExpWeight = 0.0;           // Sum of the exposure time multiplied by the global weights
    131     int numGood = 0;                    // Number of good exposures
    132123    for (int i = 0; i < values->n; i++) {
    133124        sumValueWeight += values->data.F32[i] * weights->data.F32[i];
    … …  
    136127            sumVarianceWeight += variances->data.F32[i] * PS_SQR(weights->data.F32[i]);
    137128        }
    138         if (exps) {
    139             sumExp += exps->data.F32[i];
    140             sumExpWeight += exps->data.F32[i] * weights->data.F32[i];
    141             numGood++;
    142         }
    143129    }
    144130
    … …  
    150136    if (var) {
    151137        *var = sumVarianceWeight / PS_SQR(sumWeight);
    152     }
    153     if (exp) {
    154         *exp = sumExp;
    155     }
    156     if (expWeight) {
    157         *expWeight = sumExpWeight;
    158138    }
    159139    return true;
    … …  
    296276                           const psArray *inputs, // Stack data
    297277                           const psVector *weights, // Global (single value) weights for data, or NULL
    298                            const psVector *exps,    // Exposures for data, or NULL
    299278                           const psVector *addVariance, // Additional variance for data
    300279                           const psVector *reject, // Indices of pixels to reject, or NULL
    … …  
    313292    psVector *pixelVariances = variance ? buffer->variances : NULL; // Variances for the pixel of interest
    314293    psVector *pixelWeights = buffer->weights; // Image weights for the pixel of interest
    315     psVector *pixelExps = buffer->exps;       // Exposure times
    316294    psVector *pixelSources = buffer->sources; // Sources for the pixel of interest
    317295    psVector *pixelLimits = buffer->limits; // Limits for the pixel of interest
    … …  
    353331        }
    354332        pixelWeights->data.F32[numGood] = data->weight;
    355         pixelExps->data.F32[numGood] = data->exp;
    356333        pixelSources->data.U16[numGood] = i;
    357334        numGood++;
    … …  
    370347    if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
    371348        for (int i = 0; i < numGood; i++) {
    372             fprintf(stderr, "Input %d, pixel %d,%d (%" PRIu16 "): %f %f (%f) %f %f %d\n",
     349            fprintf(stderr, "Input %d, pixel %d,%d (%" PRIu16 "): %f %f (%f) %f %d\n",
    373350                    i, x, y, pixelSources->data.U16[i], pixelData->data.F32[i], pixelVariances->data.F32[i],
    374                     addVariance->data.F32[i], pixelWeights->data.F32[i], pixelExps->data.F32[i],
    375                     pixelSuspects->data.U8[i]);
     351                    addVariance->data.F32[i], pixelWeights->data.F32[i], pixelSuspects->data.U8[i]);
    376352        }
    377353    }
    … …  
    386362                          psImage *mask, // Combined mask, for output
    387363                          psImage *variance, // Combined variance map, for output
    388                           psImage *exp,   // Exposure map (time), for output
    389                           psImage *expnum,       // Exposure map (number) for output
    390                           psImage *expweight,    // Exposure map (weighted time) for output
    391364                          int num,      // Number of good pixels
    392365                          combineBuffer *buffer, // Buffer with vectors
    393366                          int x, int y, // Coordinates of interest; frame of output image
    394367                          psImageMaskType bad, // Value for bad pixels
    395                           bool safe,           // Safe combination?
    396                           float invTotalWeight    // Inverse of total weight for all inputs
     368                          bool safe             // Safe combination?
    397369                          )
    398370{
    … …  
    400372    psVector *pixelVariances = variance ? buffer->variances : NULL; // Variances for the pixel of interest
    401373    psVector *pixelWeights = buffer->weights; // Image weights for the pixel of interest
    402     psVector *pixelExps = buffer->exps;       // Exposure times
    403374
    404375    // Default option is that the pixel is bad
    405376    float imageValue = NAN, varianceValue = NAN; // Value for combined image and variance map
    406377    psImageMaskType maskValue = bad;    // Value for combined mask
    407     float expValue = 0.0, expWeightValue = NAN; // Exposure value (straight, and weighted)
    408 
    409378    switch (num) {
    410379      case 0: {
    … …  
    424393                  varianceValue = pixelVariances->data.F32[0];
    425394              }
    426               if (exp) {
    427                   expValue = pixelExps->data.F32[0];
    428                   expWeightValue = pixelExps->data.F32[0];
    429               }
    430395              maskValue = 0;
    431396#ifdef TESTING
    … …  
    446411          // Automatically accept the mean of the pixels only if we're not playing safe
    447412          if (!safe) {
    448               if (combinationMeanVariance(&imageValue, &varianceValue, &expValue, &expWeightValue,
    449                                           pixelData, pixelVariances, pixelExps, pixelWeights)) {
     413              float mean, var;   // Mean and variance from combination
     414              if (combinationMeanVariance(&mean, &var, pixelData, pixelVariances, pixelWeights)) {
     415                  imageValue = mean;
     416                  varianceValue = var;
    450417                  maskValue = 0;
    451418#ifdef TESTING
    452419                  if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
    453420                      fprintf(stderr, "Two inputs to combine using unsafe, pixel %d,%d --> %f %f\n",
    454                               x, y, imageValue, varianceValue);
     421                              x, y, mean, var);
    455422                  }
    456423#endif
    … …  
    468435      default: {
    469436          // Can combine without too much worrying
    470           if (!combinationMeanVariance(&imageValue, &varianceValue, &expValue, &expWeightValue,
    471                                        pixelData, pixelVariances, pixelExps, pixelWeights)) {
     437          float mean, var;           // Mean and variance of the combination
     438          if (!combinationMeanVariance(&mean, &var, pixelData, pixelVariances, pixelWeights)) {
    472439              break;
    473440          }
     441          imageValue = mean;
     442          varianceValue = var;
    474443          maskValue = 0;
    475444#ifdef TESTING
    476445          if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
    477               fprintf(stderr, "Combined inputs, pixel %d,%d --> %f %f\n", x, y, imageValue, varianceValue);
     446              fprintf(stderr, "Combined inputs, pixel %d,%d --> %f %f\n", x, y, mean, var);
    478447          }
    479448#endif
    … …  
    487456        variance->data.F32[y][x] = varianceValue;
    488457    }
    489     if (exp) {
    490         exp->data.F32[y][x] = expValue;
    491     }
    492     if (expnum) {
    493         expnum->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = num;
    494     }
    495     if (expweight) {
    496         expweight->data.F32[y][x] = expWeightValue * invTotalWeight;
    497     }
     458
     459#ifdef TESTING
     460    mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = num;
     461#endif
     462
    498463
    499464    return;
    … …  
    838803                              int *numCols, int *numRows, // Size of (sky) images
    839804                              const psArray *input, // Input array of pmStackData to validate
    840                               const pmReadout *output, // Output readout
    841                               const pmReadout *exp    // Exposure map
     805                              const pmReadout *output // Output readout
    842806    )
    843807{
    … …  
    905869    }
    906870
    907     if (exp) {
    908         PM_ASSERT_READOUT_NON_NULL(exp, false);
    909         if (exp->image) {
    910             PS_ASSERT_IMAGES_SIZE_EQUAL(exp->image, output->image, false);
    911         }
    912         if (exp->mask) {
    913             PS_ASSERT_IMAGES_SIZE_EQUAL(exp->mask, output->image, false);
    914         }
    915     }
    916 
    917871    return true;
    918872}
    … …  
    983937
    984938/// Constructor
    985 pmStackData *pmStackDataAlloc(pmReadout *readout, float weight, float exp, float addVariance)
     939pmStackData *pmStackDataAlloc(pmReadout *readout, float weight, float addVariance)
    986940{
    987941    pmStackData *data = psAlloc(sizeof(pmStackData)); // Stack data, to return
    … …  
    992946    data->inspect = NULL;
    993947    data->weight = weight;
    994     data->exp = exp;
    995948    data->addVariance = addVariance;
    996949
    … …  
    999952
    1000953/// Stack input images
    1001 bool pmStackCombine(pmReadout *combined, pmReadout *expmaps, psArray *input,
    1002                     psImageMaskType maskVal, psImageMaskType maskSuspect,
    1003                     psImageMaskType bad, int kernelSize,
    1004                     float iter, float rej, float sys, float olympic,
     954bool pmStackCombine(pmReadout *combined, psArray *input, psImageMaskType maskVal, psImageMaskType maskSuspect,
     955                    psImageMaskType bad, int kernelSize, float iter, float rej, float sys, float olympic,
    1005956                    bool useVariance, bool safe, bool rejection)
    1006957{
    … …  
    1008959    int num;                            // Number of inputs
    1009960    int numCols, numRows;               // Size of (sky) images
    1010     if (!validateInputData(&haveVariances, &num, &numCols, &numRows, input, combined, expmaps)) {
     961    if (!validateInputData(&haveVariances, &num, &numCols, &numRows, input, combined)) {
    1011962        return false;
    1012963    }
    … …  
    1026977    psVector *addVariance = psVectorAlloc(num, PS_TYPE_F32); // Additional variance for each image
    1027978    psVector *weights = psVectorAlloc(num, PS_TYPE_F32); // Relative weighting for each image
    1028     psVector *exps = psVectorAlloc(num, PS_TYPE_F32);    // Exposure times for each image
    1029979    psArray *stack = psArrayAlloc(num); // Stack of readouts
    1030     float totalExpWeight = 0.0;           // Total value of all weighted exposure times
    1031     float totalExp = 0.0;                 // Total exposure time
    1032980    for (int i = 0; i < num; i++) {
    1033981        pmStackData *data = input->data[i]; // Stack data for this input
    1034982        if (!data) {
    1035983            weights->data.F32[i] = 0.0;
    1036             exps->data.F32[i] = NAN;
    1037984            continue;
    1038985        }
    1039986        weights->data.F32[i] = data->weight;
    1040         exps->data.F32[i] = data->exp;
    1041         totalExp += exps->data.F32[i];
    1042         totalExpWeight += exps->data.F32[i] * weights->data.F32[i];
    1043987        pmReadout *ro = data->readout;  // Readout of interest
    1044988        stack->data[i] = psMemIncrRefCounter(ro);
    … …  
    10591003        }
    10601004    }
    1061     totalExpWeight = totalExp / totalExpWeight;    // Convert to inverse
    10621005
    10631006    int minInputCols, maxInputCols, minInputRows, maxInputRows; // Smallest and largest values to combine
    … …  
    10651008    if (!pmReadoutStackValidate(&minInputCols, &maxInputCols, &minInputRows, &maxInputRows, &xSize, &ySize,
    10661009                                stack)) {
    1067         psError(psErrorCodeLast(), false, "Input stack is not valid.");
     1010        psError(PS_ERR_UNKNOWN, false, "Input stack is not valid.");
    10681011        psFree(stack);
    10691012        return false;
    … …  
    11021045        combined->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    11031046        combinedVariance = combined->variance;
    1104     }
    1105 
    1106     psImage *exp = NULL, *expnum = NULL, *expweight = NULL; // Exposure map and exposure number
    1107     if (expmaps) {
    1108         if (!expmaps->image) {
    1109             expmaps->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    1110         }
    1111         exp = expmaps->image;
    1112 
    1113         if (!expmaps->mask) {
    1114             expmaps->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
    1115         }
    1116         expnum = expmaps->mask;
    1117 
    1118         if (!expmaps->variance) {
    1119             expmaps->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    1120         }
    1121         expweight = expmaps->variance;
    11221047    }
    11231048
    … …  
    11581083            bool suspect;               // Suspect pixels in stack?
    11591084            combineExtract(&num, &suspect, buffer, combinedImage, combinedMask, combinedVariance,
    1160                            input, weights, exps, addVariance, reject, x, y, maskVal, maskSuspect);
    1161             combinePixels(combinedImage, combinedMask, combinedVariance, exp, expnum, expweight,
    1162                           num, buffer, x, y, bad, safe, totalExpWeight);
     1085                           input, weights, addVariance, reject, x, y, maskVal, maskSuspect);
     1086            combinePixels(combinedImage, combinedMask, combinedVariance, num, buffer, x, y, bad, safe);
    11631087
    11641088            if (iter > 0) {
    … …  
    11721096    psFree(weights);
    11731097    psFree(buffer);
    1174     psFree(addVariance);
     1098
    11751099
    11761100
Note: See TracChangeset for help on using the changeset viewer.