IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Mar 22, 2010, 8:34:28 PM (16 years ago)
Author:
Paul Price
Message:

Adding exposure map (both exposure time and number of inputs) to the outputs of ppStack. These are calculated both for the regular (convolved) and unconvolved stacks. These should compress pretty well, but I haven't enabled the compression yet. I also plan to add a weighted exposure time output.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r27310 r27400  
    3535
    3636//#define TESTING                         // Enable test output
    37 //#define TEST_X 2148-1                     // x coordinate to examine
    38 //#define TEST_Y 248-1                     // y coordinate to examine
     37//#define TEST_X 843-1                     // x coordinate to examine
     38//#define TEST_Y 813-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
    4849    psVector *sources;                  // Pixel sources (which image did they come from?)
    4950    psVector *limits;                   // Rejection limits
    … …  
    5758    psFree(buffer->variances);
    5859    psFree(buffer->weights);
     60    psFree(buffer->exps);
    5961    psFree(buffer->sources);
    6062    psFree(buffer->limits);
    … …  
    7375    buffer->variances = psVectorAlloc(numImages, PS_TYPE_F32);
    7476    buffer->weights = psVectorAlloc(numImages, PS_TYPE_F32);
     77    buffer->exps = psVectorAlloc(numImages, PS_TYPE_F32);
    7578    buffer->sources = psVectorAlloc(numImages, PS_TYPE_U16);
    7679    buffer->limits = psVectorAlloc(numImages, PS_TYPE_F32);
    … …  
    9598static bool combinationMeanVariance(float *mean, // Mean value, to return
    9699                                    float *var, // Variance value, to return
     100                                    float *exp, // Exposure time, to return
     101                                    float *expWeight,          // Weighted exposure time, to return
    97102                                    const psVector *values, // Values to combine
    98103                                    const psVector *variances, // Pixel variances to combine
     104                                    const psVector *exps,      // Exposure times to combine
    99105                                    const psVector *weights // Weights to apply
    100106                                    )
    … …  
    121127    float sumVarianceWeight = 0.0;     // Sum of the pixel variances multiplied by the global weights
    122128    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
    123131    for (int i = 0; i < values->n; i++) {
    124132        sumValueWeight += values->data.F32[i] * weights->data.F32[i];
    … …  
    127135            sumVarianceWeight += variances->data.F32[i] * PS_SQR(weights->data.F32[i]);
    128136        }
     137        if (exps) {
     138            sumExp += exps->data.F32[i];
     139            sumExpWeight += exps->data.F32[i] * weights->data.F32[i];
     140        }
    129141    }
    130142
    … …  
    136148    if (var) {
    137149        *var = sumVarianceWeight / PS_SQR(sumWeight);
     150    }
     151    if (exp) {
     152        *exp = sumExp;
     153    }
     154    if (expWeight) {
     155        *expWeight = sumExpWeight / sumWeight;
    138156    }
    139157    return true;
    … …  
    276294                           const psArray *inputs, // Stack data
    277295                           const psVector *weights, // Global (single value) weights for data, or NULL
     296                           const psVector *exps,    // Exposures for data, or NULL
    278297                           const psVector *addVariance, // Additional variance for data
    279298                           const psVector *reject, // Indices of pixels to reject, or NULL
    … …  
    292311    psVector *pixelVariances = variance ? buffer->variances : NULL; // Variances for the pixel of interest
    293312    psVector *pixelWeights = buffer->weights; // Image weights for the pixel of interest
     313    psVector *pixelExps = buffer->exps;       // Exposure times
    294314    psVector *pixelSources = buffer->sources; // Sources for the pixel of interest
    295315    psVector *pixelLimits = buffer->limits; // Limits for the pixel of interest
    … …  
    331351        }
    332352        pixelWeights->data.F32[numGood] = data->weight;
     353        pixelExps->data.F32[numGood] = data->exp;
    333354        pixelSources->data.U16[numGood] = i;
    334355        numGood++;
    … …  
    347368    if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
    348369        for (int i = 0; i < numGood; i++) {
    349             fprintf(stderr, "Input %d, pixel %d,%d (%" PRIu16 "): %f %f (%f) %f %d\n",
     370            fprintf(stderr, "Input %d, pixel %d,%d (%" PRIu16 "): %f %f (%f) %f %f %d\n",
    350371                    i, x, y, pixelSources->data.U16[i], pixelData->data.F32[i], pixelVariances->data.F32[i],
    351                     addVariance->data.F32[i], pixelWeights->data.F32[i], pixelSuspects->data.U8[i]);
     372                    addVariance->data.F32[i], pixelWeights->data.F32[i], pixelExps->data.F32[i],
     373                    pixelSuspects->data.U8[i]);
    352374        }
    353375    }
    … …  
    362384                          psImage *mask, // Combined mask, for output
    363385                          psImage *variance, // Combined variance map, for output
     386                          psImage *exp,   // Exposure map (time), for output
     387                          psImage *expnum,       // Exposure map (number) for output
     388                          psImage *expweight,    // Exposure map (weighted time) for output
    364389                          int num,      // Number of good pixels
    365390                          combineBuffer *buffer, // Buffer with vectors
    … …  
    372397    psVector *pixelVariances = variance ? buffer->variances : NULL; // Variances for the pixel of interest
    373398    psVector *pixelWeights = buffer->weights; // Image weights for the pixel of interest
     399    psVector *pixelExps = buffer->exps;       // Exposure times
    374400
    375401    // Default option is that the pixel is bad
    376402    float imageValue = NAN, varianceValue = NAN; // Value for combined image and variance map
    377403    psImageMaskType maskValue = bad;    // Value for combined mask
     404    float expValue = 0.0, expWeightValue = NAN; // Exposure value (straight, and weighted)
     405
    378406    switch (num) {
    379407      case 0: {
    … …  
    393421                  varianceValue = pixelVariances->data.F32[0];
    394422              }
     423              if (exp) {
     424                  expValue = pixelExps->data.F32[0];
     425              }
    395426              maskValue = 0;
    396427#ifdef TESTING
    … …  
    411442          // Automatically accept the mean of the pixels only if we're not playing safe
    412443          if (!safe) {
    413               float mean, var;   // Mean and variance from combination
    414               if (combinationMeanVariance(&mean, &var, pixelData, pixelVariances, pixelWeights)) {
    415                   imageValue = mean;
    416                   varianceValue = var;
     444              if (combinationMeanVariance(&imageValue, &varianceValue, &expValue, &expWeightValue,
     445                                          pixelData, pixelVariances, pixelExps, pixelWeights)) {
    417446                  maskValue = 0;
    418447#ifdef TESTING
    419448                  if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
    420449                      fprintf(stderr, "Two inputs to combine using unsafe, pixel %d,%d --> %f %f\n",
    421                               x, y, mean, var);
     450                              x, y, imageValue, varianceValue);
    422451                  }
    423452#endif
    … …  
    435464      default: {
    436465          // Can combine without too much worrying
    437           float mean, var;           // Mean and variance of the combination
    438           if (!combinationMeanVariance(&mean, &var, pixelData, pixelVariances, pixelWeights)) {
     466          if (!combinationMeanVariance(&imageValue, &varianceValue, &expValue, &expWeightValue,
     467                                       pixelData, pixelVariances, pixelExps, pixelWeights)) {
    439468              break;
    440469          }
    441           imageValue = mean;
    442           varianceValue = var;
    443470          maskValue = 0;
    444471#ifdef TESTING
    445472          if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
    446               fprintf(stderr, "Combined inputs, pixel %d,%d --> %f %f\n", x, y, mean, var);
     473              fprintf(stderr, "Combined inputs, pixel %d,%d --> %f %f\n", x, y, imageValue, varianceValue);
    447474          }
    448475#endif
    … …  
    456483        variance->data.F32[y][x] = varianceValue;
    457484    }
    458 
    459 #ifdef TESTING
    460     mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = num;
    461 #endif
    462 
     485    if (exp) {
     486        exp->data.F32[y][x] = expValue;
     487    }
     488    if (expnum) {
     489        expnum->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = num;
     490    }
     491    if (expweight) {
     492        expweight->data.F32[y][x] = expWeightValue;
     493    }
    463494
    464495    return;
    … …  
    803834                              int *numCols, int *numRows, // Size of (sky) images
    804835                              const psArray *input, // Input array of pmStackData to validate
    805                               const pmReadout *output // Output readout
     836                              const pmReadout *output, // Output readout
     837                              const pmReadout *exp    // Exposure map
    806838    )
    807839{
    … …  
    869901    }
    870902
     903    if (exp) {
     904        PM_ASSERT_READOUT_NON_NULL(exp, false);
     905        if (exp->image) {
     906            PS_ASSERT_IMAGES_SIZE_EQUAL(exp->image, output->image, false);
     907        }
     908        if (exp->mask) {
     909            PS_ASSERT_IMAGES_SIZE_EQUAL(exp->mask, output->image, false);
     910        }
     911    }
     912
    871913    return true;
    872914}
    … …  
    937979
    938980/// Constructor
    939 pmStackData *pmStackDataAlloc(pmReadout *readout, float weight, float addVariance)
     981pmStackData *pmStackDataAlloc(pmReadout *readout, float weight, float exp, float addVariance)
    940982{
    941983    pmStackData *data = psAlloc(sizeof(pmStackData)); // Stack data, to return
    … …  
    946988    data->inspect = NULL;
    947989    data->weight = weight;
     990    data->exp = exp;
    948991    data->addVariance = addVariance;
    949992
    … …  
    952995
    953996/// Stack input images
    954 bool pmStackCombine(pmReadout *combined, psArray *input, psImageMaskType maskVal, psImageMaskType maskSuspect,
    955                     psImageMaskType bad, int kernelSize, float iter, float rej, float sys, float olympic,
     997bool pmStackCombine(pmReadout *combined, pmReadout *expmaps, psArray *input,
     998                    psImageMaskType maskVal, psImageMaskType maskSuspect,
     999                    psImageMaskType bad, int kernelSize,
     1000                    float iter, float rej, float sys, float olympic,
    9561001                    bool useVariance, bool safe, bool rejection)
    9571002{
    … …  
    9591004    int num;                            // Number of inputs
    9601005    int numCols, numRows;               // Size of (sky) images
    961     if (!validateInputData(&haveVariances, &num, &numCols, &numRows, input, combined)) {
     1006    if (!validateInputData(&haveVariances, &num, &numCols, &numRows, input, combined, expmaps)) {
    9621007        return false;
    9631008    }
    … …  
    9771022    psVector *addVariance = psVectorAlloc(num, PS_TYPE_F32); // Additional variance for each image
    9781023    psVector *weights = psVectorAlloc(num, PS_TYPE_F32); // Relative weighting for each image
     1024    psVector *exps = psVectorAlloc(num, PS_TYPE_F32);    // Exposure times for each image
    9791025    psArray *stack = psArrayAlloc(num); // Stack of readouts
    9801026    for (int i = 0; i < num; i++) {
    … …  
    9821028        if (!data) {
    9831029            weights->data.F32[i] = 0.0;
     1030            exps->data.F32[i] = NAN;
    9841031            continue;
    9851032        }
    9861033        weights->data.F32[i] = data->weight;
     1034        exps->data.F32[i] = data->exp;
    9871035        pmReadout *ro = data->readout;  // Readout of interest
    9881036        stack->data[i] = psMemIncrRefCounter(ro);
    … …  
    10451093        combined->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    10461094        combinedVariance = combined->variance;
     1095    }
     1096
     1097    psImage *exp = NULL, *expnum = NULL; // Exposure map and exposure number
     1098    if (expmaps) {
     1099        if (!expmaps->image) {
     1100            expmaps->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     1101        }
     1102        exp = expmaps->image;
     1103
     1104        if (!expmaps->mask) {
     1105            expmaps->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
     1106        }
     1107        expnum = expmaps->mask;
    10471108    }
    10481109
    … …  
    10831144            bool suspect;               // Suspect pixels in stack?
    10841145            combineExtract(&num, &suspect, buffer, combinedImage, combinedMask, combinedVariance,
    1085                            input, weights, addVariance, reject, x, y, maskVal, maskSuspect);
    1086             combinePixels(combinedImage, combinedMask, combinedVariance, num, buffer, x, y, bad, safe);
     1146                           input, weights, exps, addVariance, reject, x, y, maskVal, maskSuspect);
     1147            combinePixels(combinedImage, combinedMask, combinedVariance, exp, expnum, NULL,
     1148                          num, buffer, x, y, bad, safe);
    10871149
    10881150            if (iter > 0) {
Note: See TracChangeset for help on using the changeset viewer.