Changeset 27402 for trunk/psModules/src/imcombine/pmStack.c
- Timestamp:
- Mar 23, 2010, 9:12:53 AM (16 years ago)
- File:
-
- 1 edited
-
trunk/psModules/src/imcombine/pmStack.c (modified) (28 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/imcombine/pmStack.c
r27400 r27402 35 35 36 36 //#define TESTING // Enable test output 37 //#define TEST_X 843-1 // x coordinate to examine38 //#define TEST_Y 813-1 // y coordinate to examine37 //#define TEST_X 2148-1 // x coordinate to examine 38 //#define TEST_Y 248-1 // y coordinate to examine 39 39 //#define TEST_RADIUS 0 // Radius to examine 40 40 … … 46 46 psVector *variances; // Pixel variances 47 47 psVector *weights; // Pixel weightings 48 psVector *exps; // Pixel exposures49 48 psVector *sources; // Pixel sources (which image did they come from?) 50 49 psVector *limits; // Rejection limits … … 58 57 psFree(buffer->variances); 59 58 psFree(buffer->weights); 60 psFree(buffer->exps);61 59 psFree(buffer->sources); 62 60 psFree(buffer->limits); … … 75 73 buffer->variances = psVectorAlloc(numImages, PS_TYPE_F32); 76 74 buffer->weights = psVectorAlloc(numImages, PS_TYPE_F32); 77 buffer->exps = psVectorAlloc(numImages, PS_TYPE_F32);78 75 buffer->sources = psVectorAlloc(numImages, PS_TYPE_U16); 79 76 buffer->limits = psVectorAlloc(numImages, PS_TYPE_F32); … … 98 95 static bool combinationMeanVariance(float *mean, // Mean value, to return 99 96 float *var, // Variance value, to return 100 float *exp, // Exposure time, to return101 float *expWeight, // Weighted exposure time, to return102 97 const psVector *values, // Values to combine 103 98 const psVector *variances, // Pixel variances to combine 104 const psVector *exps, // Exposure times to combine105 99 const psVector *weights // Weights to apply 106 100 ) … … 127 121 float sumVarianceWeight = 0.0; // Sum of the pixel variances multiplied by the global weights 128 122 float sumWeight = 0.0; // Sum of the image weights 129 float sumExp = 0.0; // Sum of the exposure time130 float sumExpWeight = 0.0; // Sum of the exposure time multiplied by the global weights131 123 for (int i = 0; i < values->n; i++) { 132 124 sumValueWeight += values->data.F32[i] * weights->data.F32[i]; … … 135 127 sumVarianceWeight += variances->data.F32[i] * PS_SQR(weights->data.F32[i]); 136 128 } 137 if (exps) {138 sumExp += exps->data.F32[i];139 sumExpWeight += exps->data.F32[i] * weights->data.F32[i];140 }141 129 } 142 130 … … 148 136 if (var) { 149 137 *var = sumVarianceWeight / PS_SQR(sumWeight); 150 }151 if (exp) {152 *exp = sumExp;153 }154 if (expWeight) {155 *expWeight = sumExpWeight / sumWeight;156 138 } 157 139 return true; … … 294 276 const psArray *inputs, // Stack data 295 277 const psVector *weights, // Global (single value) weights for data, or NULL 296 const psVector *exps, // Exposures for data, or NULL297 278 const psVector *addVariance, // Additional variance for data 298 279 const psVector *reject, // Indices of pixels to reject, or NULL … … 311 292 psVector *pixelVariances = variance ? buffer->variances : NULL; // Variances for the pixel of interest 312 293 psVector *pixelWeights = buffer->weights; // Image weights for the pixel of interest 313 psVector *pixelExps = buffer->exps; // Exposure times314 294 psVector *pixelSources = buffer->sources; // Sources for the pixel of interest 315 295 psVector *pixelLimits = buffer->limits; // Limits for the pixel of interest … … 351 331 } 352 332 pixelWeights->data.F32[numGood] = data->weight; 353 pixelExps->data.F32[numGood] = data->exp;354 333 pixelSources->data.U16[numGood] = i; 355 334 numGood++; … … 368 347 if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) { 369 348 for (int i = 0; i < numGood; i++) { 370 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", 371 350 i, x, y, pixelSources->data.U16[i], pixelData->data.F32[i], pixelVariances->data.F32[i], 372 addVariance->data.F32[i], pixelWeights->data.F32[i], pixelExps->data.F32[i], 373 pixelSuspects->data.U8[i]); 351 addVariance->data.F32[i], pixelWeights->data.F32[i], pixelSuspects->data.U8[i]); 374 352 } 375 353 } … … 384 362 psImage *mask, // Combined mask, for output 385 363 psImage *variance, // Combined variance map, for output 386 psImage *exp, // Exposure map (time), for output387 psImage *expnum, // Exposure map (number) for output388 psImage *expweight, // Exposure map (weighted time) for output389 364 int num, // Number of good pixels 390 365 combineBuffer *buffer, // Buffer with vectors … … 397 372 psVector *pixelVariances = variance ? buffer->variances : NULL; // Variances for the pixel of interest 398 373 psVector *pixelWeights = buffer->weights; // Image weights for the pixel of interest 399 psVector *pixelExps = buffer->exps; // Exposure times400 374 401 375 // Default option is that the pixel is bad 402 376 float imageValue = NAN, varianceValue = NAN; // Value for combined image and variance map 403 377 psImageMaskType maskValue = bad; // Value for combined mask 404 float expValue = 0.0, expWeightValue = NAN; // Exposure value (straight, and weighted)405 406 378 switch (num) { 407 379 case 0: { … … 421 393 varianceValue = pixelVariances->data.F32[0]; 422 394 } 423 if (exp) {424 expValue = pixelExps->data.F32[0];425 }426 395 maskValue = 0; 427 396 #ifdef TESTING … … 442 411 // Automatically accept the mean of the pixels only if we're not playing safe 443 412 if (!safe) { 444 if (combinationMeanVariance(&imageValue, &varianceValue, &expValue, &expWeightValue, 445 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; 446 417 maskValue = 0; 447 418 #ifdef TESTING 448 419 if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) { 449 420 fprintf(stderr, "Two inputs to combine using unsafe, pixel %d,%d --> %f %f\n", 450 x, y, imageValue, varianceValue);421 x, y, mean, var); 451 422 } 452 423 #endif … … 464 435 default: { 465 436 // Can combine without too much worrying 466 if (!combinationMeanVariance(&imageValue, &varianceValue, &expValue, &expWeightValue,467 pixelData, pixelVariances, pixelExps, pixelWeights)) {437 float mean, var; // Mean and variance of the combination 438 if (!combinationMeanVariance(&mean, &var, pixelData, pixelVariances, pixelWeights)) { 468 439 break; 469 440 } 441 imageValue = mean; 442 varianceValue = var; 470 443 maskValue = 0; 471 444 #ifdef TESTING 472 445 if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) { 473 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); 474 447 } 475 448 #endif … … 483 456 variance->data.F32[y][x] = varianceValue; 484 457 } 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 } 458 459 #ifdef TESTING 460 mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = num; 461 #endif 462 494 463 495 464 return; … … 834 803 int *numCols, int *numRows, // Size of (sky) images 835 804 const psArray *input, // Input array of pmStackData to validate 836 const pmReadout *output, // Output readout 837 const pmReadout *exp // Exposure map 805 const pmReadout *output // Output readout 838 806 ) 839 807 { … … 901 869 } 902 870 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 913 871 return true; 914 872 } … … 979 937 980 938 /// Constructor 981 pmStackData *pmStackDataAlloc(pmReadout *readout, float weight, float exp, floataddVariance)939 pmStackData *pmStackDataAlloc(pmReadout *readout, float weight, float addVariance) 982 940 { 983 941 pmStackData *data = psAlloc(sizeof(pmStackData)); // Stack data, to return … … 988 946 data->inspect = NULL; 989 947 data->weight = weight; 990 data->exp = exp;991 948 data->addVariance = addVariance; 992 949 … … 995 952 996 953 /// Stack input images 997 bool 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, 954 bool pmStackCombine(pmReadout *combined, psArray *input, psImageMaskType maskVal, psImageMaskType maskSuspect, 955 psImageMaskType bad, int kernelSize, float iter, float rej, float sys, float olympic, 1001 956 bool useVariance, bool safe, bool rejection) 1002 957 { … … 1004 959 int num; // Number of inputs 1005 960 int numCols, numRows; // Size of (sky) images 1006 if (!validateInputData(&haveVariances, &num, &numCols, &numRows, input, combined , expmaps)) {961 if (!validateInputData(&haveVariances, &num, &numCols, &numRows, input, combined)) { 1007 962 return false; 1008 963 } … … 1022 977 psVector *addVariance = psVectorAlloc(num, PS_TYPE_F32); // Additional variance for each image 1023 978 psVector *weights = psVectorAlloc(num, PS_TYPE_F32); // Relative weighting for each image 1024 psVector *exps = psVectorAlloc(num, PS_TYPE_F32); // Exposure times for each image1025 979 psArray *stack = psArrayAlloc(num); // Stack of readouts 1026 980 for (int i = 0; i < num; i++) { … … 1028 982 if (!data) { 1029 983 weights->data.F32[i] = 0.0; 1030 exps->data.F32[i] = NAN;1031 984 continue; 1032 985 } 1033 986 weights->data.F32[i] = data->weight; 1034 exps->data.F32[i] = data->exp;1035 987 pmReadout *ro = data->readout; // Readout of interest 1036 988 stack->data[i] = psMemIncrRefCounter(ro); … … 1093 1045 combined->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); 1094 1046 combinedVariance = combined->variance; 1095 }1096 1097 psImage *exp = NULL, *expnum = NULL; // Exposure map and exposure number1098 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;1108 1047 } 1109 1048 … … 1144 1083 bool suspect; // Suspect pixels in stack? 1145 1084 combineExtract(&num, &suspect, buffer, combinedImage, combinedMask, combinedVariance, 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); 1085 input, weights, addVariance, reject, x, y, maskVal, maskSuspect); 1086 combinePixels(combinedImage, combinedMask, combinedVariance, num, buffer, x, y, bad, safe); 1149 1087 1150 1088 if (iter > 0) {
Note:
See TracChangeset
for help on using the changeset viewer.
