IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Mar 23, 2010, 3:23:07 PM (16 years ago)
Author:
Paul Price
Message:

Merging development branch branches/pap_stack/ (exposure maps and their compression).

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk

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

    r27402 r27417  
    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
     131    int numGood = 0;                    // Number of good exposures
    123132    for (int i = 0; i < values->n; i++) {
    124133        sumValueWeight += values->data.F32[i] * weights->data.F32[i];
    … …  
    127136            sumVarianceWeight += variances->data.F32[i] * PS_SQR(weights->data.F32[i]);
    128137        }
     138        if (exps) {
     139            sumExp += exps->data.F32[i];
     140            sumExpWeight += exps->data.F32[i] * weights->data.F32[i];
     141            numGood++;
     142        }
    129143    }
    130144
    … …  
    136150    if (var) {
    137151        *var = sumVarianceWeight / PS_SQR(sumWeight);
     152    }
     153    if (exp) {
     154        *exp = sumExp;
     155    }
     156    if (expWeight) {
     157        *expWeight = sumExpWeight;
    138158    }
    139159    return true;
    … …  
    276296                           const psArray *inputs, // Stack data
    277297                           const psVector *weights, // Global (single value) weights for data, or NULL
     298                           const psVector *exps,    // Exposures for data, or NULL
    278299                           const psVector *addVariance, // Additional variance for data
    279300                           const psVector *reject, // Indices of pixels to reject, or NULL
    … …  
    292313    psVector *pixelVariances = variance ? buffer->variances : NULL; // Variances for the pixel of interest
    293314    psVector *pixelWeights = buffer->weights; // Image weights for the pixel of interest
     315    psVector *pixelExps = buffer->exps;       // Exposure times
    294316    psVector *pixelSources = buffer->sources; // Sources for the pixel of interest
    295317    psVector *pixelLimits = buffer->limits; // Limits for the pixel of interest
    … …  
    331353        }
    332354        pixelWeights->data.F32[numGood] = data->weight;
     355        pixelExps->data.F32[numGood] = data->exp;
    333356        pixelSources->data.U16[numGood] = i;
    334357        numGood++;
    … …  
    347370    if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
    348371        for (int i = 0; i < numGood; i++) {
    349             fprintf(stderr, "Input %d, pixel %d,%d (%" PRIu16 "): %f %f (%f) %f %d\n",
     372            fprintf(stderr, "Input %d, pixel %d,%d (%" PRIu16 "): %f %f (%f) %f %f %d\n",
    350373                    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]);
     374                    addVariance->data.F32[i], pixelWeights->data.F32[i], pixelExps->data.F32[i],
     375                    pixelSuspects->data.U8[i]);
    352376        }
    353377    }
    … …  
    362386                          psImage *mask, // Combined mask, for output
    363387                          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
    364391                          int num,      // Number of good pixels
    365392                          combineBuffer *buffer, // Buffer with vectors
    366393                          int x, int y, // Coordinates of interest; frame of output image
    367394                          psImageMaskType bad, // Value for bad pixels
    368                           bool safe             // Safe combination?
     395                          bool safe,           // Safe combination?
     396                          float invTotalWeight    // Inverse of total weight for all inputs
    369397                          )
    370398{
    … …  
    372400    psVector *pixelVariances = variance ? buffer->variances : NULL; // Variances for the pixel of interest
    373401    psVector *pixelWeights = buffer->weights; // Image weights for the pixel of interest
     402    psVector *pixelExps = buffer->exps;       // Exposure times
    374403
    375404    // Default option is that the pixel is bad
    376405    float imageValue = NAN, varianceValue = NAN; // Value for combined image and variance map
    377406    psImageMaskType maskValue = bad;    // Value for combined mask
     407    float expValue = 0.0, expWeightValue = NAN; // Exposure value (straight, and weighted)
     408
    378409    switch (num) {
    379410      case 0: {
    … …  
    393424                  varianceValue = pixelVariances->data.F32[0];
    394425              }
     426              if (exp) {
     427                  expValue = pixelExps->data.F32[0];
     428                  expWeightValue = pixelExps->data.F32[0];
     429              }
    395430              maskValue = 0;
    396431#ifdef TESTING
    … …  
    411446          // Automatically accept the mean of the pixels only if we're not playing safe
    412447          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;
     448              if (combinationMeanVariance(&imageValue, &varianceValue, &expValue, &expWeightValue,
     449                                          pixelData, pixelVariances, pixelExps, pixelWeights)) {
    417450                  maskValue = 0;
    418451#ifdef TESTING
    419452                  if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
    420453                      fprintf(stderr, "Two inputs to combine using unsafe, pixel %d,%d --> %f %f\n",
    421                               x, y, mean, var);
     454                              x, y, imageValue, varianceValue);
    422455                  }
    423456#endif
    … …  
    435468      default: {
    436469          // Can combine without too much worrying
    437           float mean, var;           // Mean and variance of the combination
    438           if (!combinationMeanVariance(&mean, &var, pixelData, pixelVariances, pixelWeights)) {
     470          if (!combinationMeanVariance(&imageValue, &varianceValue, &expValue, &expWeightValue,
     471                                       pixelData, pixelVariances, pixelExps, pixelWeights)) {
    439472              break;
    440473          }
    441           imageValue = mean;
    442           varianceValue = var;
    443474          maskValue = 0;
    444475#ifdef TESTING
    445476          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);
     477              fprintf(stderr, "Combined inputs, pixel %d,%d --> %f %f\n", x, y, imageValue, varianceValue);
    447478          }
    448479#endif
    … …  
    456487        variance->data.F32[y][x] = varianceValue;
    457488    }
    458 
    459 #ifdef TESTING
    460     mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = num;
    461 #endif
    462 
     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    }
    463498
    464499    return;
    … …  
    803838                              int *numCols, int *numRows, // Size of (sky) images
    804839                              const psArray *input, // Input array of pmStackData to validate
    805                               const pmReadout *output // Output readout
     840                              const pmReadout *output, // Output readout
     841                              const pmReadout *exp    // Exposure map
    806842    )
    807843{
    … …  
    869905    }
    870906
     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
    871917    return true;
    872918}
    … …  
    937983
    938984/// Constructor
    939 pmStackData *pmStackDataAlloc(pmReadout *readout, float weight, float addVariance)
     985pmStackData *pmStackDataAlloc(pmReadout *readout, float weight, float exp, float addVariance)
    940986{
    941987    pmStackData *data = psAlloc(sizeof(pmStackData)); // Stack data, to return
    … …  
    946992    data->inspect = NULL;
    947993    data->weight = weight;
     994    data->exp = exp;
    948995    data->addVariance = addVariance;
    949996
    … …  
    952999
    9531000/// 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,
     1001bool 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,
    9561005                    bool useVariance, bool safe, bool rejection)
    9571006{
    … …  
    9591008    int num;                            // Number of inputs
    9601009    int numCols, numRows;               // Size of (sky) images
    961     if (!validateInputData(&haveVariances, &num, &numCols, &numRows, input, combined)) {
     1010    if (!validateInputData(&haveVariances, &num, &numCols, &numRows, input, combined, expmaps)) {
    9621011        return false;
    9631012    }
    … …  
    9771026    psVector *addVariance = psVectorAlloc(num, PS_TYPE_F32); // Additional variance for each image
    9781027    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
    9791029    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
    9801032    for (int i = 0; i < num; i++) {
    9811033        pmStackData *data = input->data[i]; // Stack data for this input
    9821034        if (!data) {
    9831035            weights->data.F32[i] = 0.0;
     1036            exps->data.F32[i] = NAN;
    9841037            continue;
    9851038        }
    9861039        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];
    9871043        pmReadout *ro = data->readout;  // Readout of interest
    9881044        stack->data[i] = psMemIncrRefCounter(ro);
    … …  
    10031059        }
    10041060    }
     1061    totalExpWeight = totalExp / totalExpWeight;    // Convert to inverse
    10051062
    10061063    int minInputCols, maxInputCols, minInputRows, maxInputRows; // Smallest and largest values to combine
    … …  
    10451102        combined->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    10461103        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;
    10471122    }
    10481123
    … …  
    10831158            bool suspect;               // Suspect pixels in stack?
    10841159            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);
     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);
    10871163
    10881164            if (iter > 0) {
Note: See TracChangeset for help on using the changeset viewer.