Changeset 21476
- Timestamp:
- Feb 13, 2009, 1:52:14 PM (17 years ago)
- Location:
- trunk/psModules/src/imcombine
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/imcombine/pmStack.c
r21363 r21476 8 8 * @author GLG, MHPCC 9 9 * 10 * @version $Revision: 1.4 7$ $Name: not supported by cvs2svn $11 * @date $Date: 2009-02- 06 02:31:25$10 * @version $Revision: 1.48 $ $Name: not supported by cvs2svn $ 11 * @date $Date: 2009-02-13 23:52:14 $ 12 12 * Copyright 2004-2007 Institute for Astronomy, University of Hawaii 13 13 * … … 30 30 #define PIXEL_LIST_BUFFER 100 // Number of entries to add to pixel list at a time 31 31 #define PIXEL_MAP_BUFFER 2 // Number of entries to add to pixel map at a time 32 #define VARIANCE_FACTORS // Use variance factors when calculating the variances? 33 //#define ADD_VARIANCE // Allow additional variance (besides variance factor)? 32 #define ADD_VARIANCE // Allow additional variance (besides variance factor)? 34 33 #define NUM_DIRECT_STDEV 5 // For less than this number of values, measure stdev directly 34 35 //#define TESTING // Enable test output 36 //#define TEST_X 2318 // x coordinate to examine 37 //#define TEST_Y 2306 // y coordinate to examine 35 38 36 39 … … 190 193 } 191 194 195 #if 0 196 // Return the weighted median for the pixels 197 // This does not appear to produce as clean images as the weighted Olympic mean 198 static float combinationWeightedMedian(const psVector *values, // Values to combine 199 const psVector *weights, // Weights to combine 200 const psVector *masks, // Mask to apply 201 psVector *sortBuffer // Buffer for sorting 202 ) 203 { 204 double sumWeight = 0.0; // Sum of weights 205 for (int j = 0; j < values->n; j++) { 206 if (masks->data.PS_TYPE_VECTOR_MASK_DATA[j]) { 207 continue; 208 } 209 sumWeight += weights->data.F32[j]; 210 } 211 212 sortBuffer = psVectorSortIndex(sortBuffer, values); 213 double target = sumWeight / 2.0; // Target weight 214 215 int dominant = -1; // Index of dominant value, if any 216 double cumulativeWeight = 0.0; // Sum of weights 217 for (int j = 0; j < values->n; j++) { 218 if (masks->data.PS_TYPE_VECTOR_MASK_DATA[j]) { 219 continue; 220 } 221 int index = sortBuffer->data.S32[j]; // Index of value of interest 222 float weight = weights->data.F32[index]; // Weight for value of interest 223 if (weight >= target) { 224 // Get the weighted median of the rest 225 dominant = index; 226 sumWeight -= weight; 227 target = sumWeight / 2.0; 228 continue; 229 } 230 cumulativeWeight += weight; 231 if (cumulativeWeight >= target) { 232 float median = values->data.F32[index]; // Weighted median median 233 if (dominant != -1) { 234 // In the case that a single value contains a disproportionate weight compared to the rest, 235 // we use a weighted mean between that dominant value and the weighted median of the rest. 236 return (values->data.F32[dominant] * weights->data.F32[dominant] + median * sumWeight) / 237 (weights->data.F32[dominant] + sumWeight); 238 } 239 return median; 240 } 241 } 242 243 return NAN; 244 } 245 #endif 246 247 // Return the weighted Olympic mean for the pixels 248 static float combinationWeightedOlympic(const psVector *values, // Values to combine 249 const psVector *weights, // Weights to combine 250 const psVector *masks, // Mask to apply 251 float frac, // Fraction to discard 252 psVector *sortBuffer // Buffer for sorting 253 ) 254 { 255 int numGood = 0; // Number of good values 256 for (int i = 0; i < values->n; i++) { 257 if (masks->data.PS_TYPE_VECTOR_MASK_DATA[i]) { 258 continue; 259 } 260 numGood++; 261 } 262 263 int numBad = frac * numGood + 0.5; // Number of bad values 264 int low = numBad / 2, high = low + numGood; // Indices (modulo masked pixels) delimiting range of interest 265 266 sortBuffer = psVectorSortIndex(sortBuffer, values); 267 268 double sumValues = 0.0, sumWeight = 0.0; // Sums for weighted mean 269 for (int i = 0, j = 0; i < values->n; i++) { 270 int index = sortBuffer->data.S32[i]; // Index of interest 271 if (masks->data.PS_TYPE_VECTOR_MASK_DATA[index]) { 272 continue; 273 } 274 j++; 275 if (j > high) { 276 break; 277 } 278 if (j <= low) { 279 continue; 280 } 281 sumValues += values->data.F32[index] * weights->data.F32[index]; 282 sumWeight += weights->data.F32[index]; 283 } 284 285 return sumValues / sumWeight; 286 } 192 287 193 288 // Mark a pixel for inspection … … 213 308 const psArray *inputs, // Stack data 214 309 const psVector *weights, // Global (single value) weights for data, or NULL 215 const psVector * varFactors, // Variance factorsfor data310 const psVector *addVariance, // Additional variance for data 216 311 const psVector *reject, // Indices of pixels to reject, or NULL 217 312 int x, int y, // Coordinates of interest; frame of output image … … 221 316 float rej, // Number of standard deviations at which to reject 222 317 float sys, // Relative systematic error 318 float discard,// Fraction of values to discard (Olympic weighted mean) 223 319 bool useVariance, // Use variance for rejection when combining? 224 320 bool safe, // Combine safely? … … 232 328 assert(numIter >= 0); 233 329 assert(buffer); 234 assert( varFactors);330 assert(addVariance); 235 331 assert((useVariance && variance) || !useVariance); 236 332 … … 267 363 pixelData->data.F32[num] = image->data.F32[yIn][xIn]; 268 364 if (variance) { 269 pixelVariances->data.F32[num] = variance->data.F32[yIn][xIn] ;365 pixelVariances->data.F32[num] = variance->data.F32[yIn][xIn] * addVariance->data.F32[i]; 270 366 } 271 367 pixelWeights->data.F32[num] = data->weight; … … 280 376 pixelSources->n = num; 281 377 378 #ifdef TESTING 379 if (x == TEST_X && y == TEST_Y) { 380 for (int i = 0; i < num; i++) { 381 fprintf(stderr, "Input %d (%" PRIu16 "): %f %f %f\n", 382 i, pixelSources->data.U16[i], pixelData->data.F32[i], pixelVariances->data.F32[i], 383 pixelWeights->data.F32[i]); 384 } 385 } 386 #endif 282 387 283 388 // The sensible thing to do varies according to how many good pixels there are. … … 288 393 case 0: 289 394 // Nothing to combine: it's bad 395 #ifdef TESTING 396 if (x == TEST_X && y == TEST_Y) { 397 fprintf(stderr, "No inputs to combine, pixel is bad.\n"); 398 } 399 #endif 290 400 break; 291 401 case 1: { 292 402 // Accept the single pixel unless we have to be safe 293 403 if (!safe) { 404 #ifdef TESTING 405 if (x == TEST_X && y == TEST_Y) { 406 fprintf(stderr, "Single input to combine, safety off.\n"); 407 } 408 #endif 294 409 imageValue = pixelData->data.F32[0]; 295 410 if (variance) { … … 298 413 maskValue = 0; 299 414 } 415 #ifdef TESTING 416 else if (x == TEST_X && y == TEST_Y) { 417 fprintf(stderr, "Single input to combine, safety on, pixel is bad.\n"); 418 } 419 #endif 300 420 break; 301 421 } … … 309 429 varianceValue = var; 310 430 maskValue = 0; 431 #ifdef TESTING 432 if (x == TEST_X && y == TEST_Y) { 433 fprintf(stderr, "Two inputs to combine using variance/unsafe --> %f %f\n", 434 mean, var); 435 } 436 #endif 311 437 } 312 438 } … … 316 442 float var1 = pixelVariances->data.F32[0]; // Variance of first 317 443 float var2 = pixelVariances->data.F32[1]; // Variance of second 318 #ifdef VARIANCE_FACTORS319 // Variance factor contributes to the rejection level320 var1 *= varFactors->data.F32[pixelSources->data.U16[0]];321 var2 *= varFactors->data.F32[pixelSources->data.U16[1]];322 #endif323 444 // Systematic error contributes to the rejection level 324 445 var1 += PS_SQR(sys * pixelData->data.F32[0]); 325 446 var2 += PS_SQR(sys * pixelData->data.F32[1]); 326 447 327 float sigma2 = var1 + var2; // 448 float sigma2 = var1 + var2; // Combined variance 328 449 if (PS_SQR(diff) > PS_SQR(rej) * sigma2) { 329 450 // Not consistent: mark both for inspection 330 451 combineInspect(inputs, x, y, pixelSources->data.U16[0]); 331 452 combineInspect(inputs, x, y, pixelSources->data.U16[1]); 453 #ifdef TESTING 454 if (x == TEST_X && y == TEST_Y) { 455 fprintf(stderr, "Both pixels marked for inspection (%f > %f x %f\n)", 456 diff, rej, sqrtf(sigma2)); 457 } 458 #endif 332 459 } 333 460 } … … 344 471 varianceValue = var; 345 472 maskValue = 0; 473 #ifdef TESTING 474 if (x == TEST_X && y == TEST_Y) { 475 fprintf(stderr, "Combined inputs: %f %f\n", mean, var); 476 } 477 #endif 346 478 347 479 // Prepare for clipping iteration … … 353 485 // Using squared rejection limit because it's cheaper than sqrts 354 486 float rej2 = PS_SQR(rej); // Rejection level squared 487 double sumWeights = 0.0; 488 for (int i = 0; i < num; i++) { 489 sumWeights += pixelWeights->data.F32[i]; 490 } 355 491 for (int i = 0; i < num; i++) { 356 492 // Systematic error contributes to the rejection level 357 pixelVariances->data.F32[i] += PS_SQR(sys * pixelData->data.F32[i]); 358 #ifdef VARIANCE_FACTORS 359 // Variance factor contributes to the rejection level 360 pixelVariances->data.F32[i] *= varFactors->data.F32[pixelSources->data.U16[i]]; 361 #endif 362 pixelVariances->data.F32[i] *= rej2; 493 float sysVar = PS_SQR(sys * pixelData->data.F32[i]); 494 #ifdef TESTING 495 // Correct variance for comparison against weighted mean including itself 496 float compare = 1.0 - pixelWeights->data.F32[i] / sumWeights; 497 if (x == TEST_X && y == TEST_Y) { 498 fprintf(stderr, "Variance %d (%d): %f %f %f\n", i, pixelSources->data.U16[i], 499 pixelVariances->data.F32[i], sysVar, compare); 500 } 501 #endif 502 503 pixelVariances->data.F32[i] = rej2 * (pixelVariances->data.F32[i] + sysVar); 363 504 } 364 505 } … … 371 512 for (int i = 0; i < numIter && numClipped > 0 && num - totalClipped > 2; i++) { 372 513 numClipped = 0; 373 float median, stdev; // Median and stdev of the combination, for rejection 374 375 if (!combinationMedianStdev(&median, useVariance ? NULL : &stdev, 376 pixelData, pixelMasks, sort)) { 377 psWarning("Bad median/stdev at %d,%d", x, y); 378 break; 514 float median = NAN; // Middle of distribution 515 float limit = NAN; // Rejection limit 516 if (!useVariance) { 517 float stdev; // Median and stdev of the combination, for rejection 518 if (!combinationMedianStdev(&median, useVariance ? NULL : &stdev, 519 pixelData, pixelMasks, sort)) { 520 psWarning("Bad median/stdev at %d,%d", x, y); 521 break; 522 } 523 limit = rej * stdev; 524 #ifdef TESTING 525 if (x == TEST_X && y == TEST_Y) { 526 fprintf(stderr, "Rejection limit: %f\n", limit); 527 } 528 #endif 529 } else { 530 median = combinationWeightedOlympic(pixelData, pixelWeights, pixelMasks, discard, sort); 379 531 } 380 532 381 float limit = rej * stdev; // Rejection limit, when rejecting based on data properties 533 #ifdef TESTING 534 if (x == TEST_X && y == TEST_Y) { 535 fprintf(stderr, "Median: %f\n", median); 536 } 537 #endif 538 382 539 383 540 // Mask a pixel for inspection … … 395 552 if (useVariance) { 396 553 // Comparing squares --- cheaper than lots of sqrts 397 // pixelVariances includes the variance factor and therejection limit, from above554 // pixelVariances includes the rejection limit, from above 398 555 if (PS_SQR(diff) > pixelVariances->data.F32[j]) { 399 556 MASK_PIXEL_FOR_INSPECTION(); 557 #ifdef TESTING 558 if (x == TEST_X && y == TEST_Y) { 559 fprintf(stderr, "Rejecting input %d based on variance: %f > %f\n", 560 j, diff, sqrtf(pixelVariances->data.F32[j])); 561 } 562 #endif 400 563 } 401 564 } else if (fabsf(diff) > limit) { 402 565 MASK_PIXEL_FOR_INSPECTION(); 566 #ifdef TESTING 567 if (x == TEST_X && y == TEST_Y) { 568 fprintf(stderr, "Rejecting input %d based on distribution: %f > %f\n", 569 j, diff, limit); 570 } 571 #endif 403 572 } 404 573 } … … 564 733 /// Stack input images 565 734 bool pmStackCombine(pmReadout *combined, psArray *input, psImageMaskType maskVal, psImageMaskType bad, 566 int kernelSize, int numIter, float rej, float sys, 735 int kernelSize, int numIter, float rej, float sys, float discard, 567 736 bool entire, bool useVariance, bool safe) 568 737 { … … 596 765 } 597 766 598 psVector * varFactors = psVectorAlloc(num, PS_TYPE_F32); // Variance factorsfor each image767 psVector *addVariance = psVectorAlloc(num, PS_TYPE_F32); // Additional variance for each image 599 768 psVector *weights = psVectorAlloc(num, PS_TYPE_F32); // Relative weighting for each image 600 769 psArray *stack = psArrayAlloc(num); // Stack of readouts … … 606 775 } 607 776 weights->data.F32[i] = data->weight; 608 stack->data[i] = psMemIncrRefCounter(data->readout);609 varFactors->data.F32[i] = data->readout->covariance ?610 psImageCovarianceFactor(data->readout->covariance) : 1.0;611 #if 0777 pmReadout *ro = data->readout; // Readout of interest 778 stack->data[i] = psMemIncrRefCounter(ro); 779 addVariance->data.F32[i] = ro->covariance ? psImageCovarianceFactor(ro->covariance) : 1.0; 780 #ifdef ADD_VARIANCE 612 781 if (isfinite(data->addVariance)) { 613 varFactors->data.F32[i] *= data->addVariance;782 addVariance->data.F32[i] *= data->addVariance; 614 783 } 615 784 #endif … … 667 836 psVector *reject = pixelMapQuery(pixelMap, minInputCols, minInputRows, 668 837 x, y); // Reject these images 669 combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, varFactors, 670 reject, x, y, maskVal, bad, numIter, rej, sys, useVariance, safe, buffer); 838 combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, 839 addVariance, reject, x, y, maskVal, bad, numIter, rej, sys, discard, 840 useVariance, safe, buffer); 671 841 } 672 842 } … … 681 851 psVector *reject = pixelMapQuery(pixelMap, minInputCols, minInputRows, 682 852 x, y); // Reject these images 683 combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, varFactors, 684 reject, x, y, maskVal, bad, numIter, rej, sys, useVariance, safe, buffer); 853 combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, 854 addVariance, reject, x, y, maskVal, bad, numIter, rej, sys, discard, 855 useVariance, safe, buffer); 685 856 } 686 857 } … … 708 879 for (int y = minInputRows; y < maxInputRows; y++) { 709 880 for (int x = minInputCols; x < maxInputCols; x++) { 710 combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, varFactors, 711 NULL, x, y, maskVal, bad, numIter, rej, sys, useVariance, safe, buffer); 881 combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, 882 addVariance, NULL, x, y, maskVal, bad, numIter, rej, sys, discard, 883 useVariance, safe, buffer); 712 884 } 713 885 } -
trunk/psModules/src/imcombine/pmStack.h
r21183 r21476 8 8 * @author GLG, MHPCC 9 9 * 10 * @version $Revision: 1.1 0$ $Name: not supported by cvs2svn $11 * @date $Date: 2009-0 1-27 06:39:38$10 * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $ 11 * @date $Date: 2009-02-13 23:52:14 $ 12 12 * 13 13 * Copyright 2004-2007 Institute for Astronomy, University of Hawaii … … 49 49 float rej, ///< Rejection limit (standard deviations) 50 50 float sys, ///< Relative systematic error 51 float discard, ///< Fraction of values to discard for Olympic weighted mean 51 52 bool entire, ///< Combine entire image even if rejection lists provided? 52 53 bool useVariance, ///< Use variance values for rejection?
Note:
See TracChangeset
for help on using the changeset viewer.
