Changeset 13735
- Timestamp:
- Jun 8, 2007, 3:04:02 PM (19 years ago)
- Location:
- trunk/psModules/src/imcombine
- Files:
-
- 4 edited
-
pmSubtraction.c (modified) (9 diffs)
-
pmSubtraction.h (modified) (4 diffs)
-
pmSubtractionStamps.c (modified) (4 diffs)
-
pmSubtractionStamps.h (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/imcombine/pmSubtraction.c
r13511 r13735 4 4 * @author GLG, MHPCC 5 5 * 6 * @version $Revision: 1.1 2$ $Name: not supported by cvs2svn $7 * @date $Date: 2007-0 5-24 21:17:01$6 * @version $Revision: 1.13 $ $Name: not supported by cvs2svn $ 7 * @date $Date: 2007-06-09 01:04:02 $ 8 8 * 9 9 * Copyright 2004-2007 Institute for Astronomy, University of Hawaii … … 333 333 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 334 334 335 psImage *pmSubtractionMask(const psImage *inMask, const psImage *refMask, psMaskType maskVal, 336 int size, int footprint) 337 { 338 PS_ASSERT_IMAGE_NON_NULL(inMask, NULL); 339 PS_ASSERT_IMAGE_TYPE(inMask, PS_TYPE_MASK, NULL); 340 PS_ASSERT_IMAGE_NON_NULL(refMask, NULL); 341 PS_ASSERT_IMAGE_TYPE(refMask, PS_TYPE_MASK, NULL); 342 PS_ASSERT_IMAGES_SIZE_EQUAL(inMask, refMask, NULL); 343 PS_ASSERT_INT_NONNEGATIVE(size, NULL); 344 PS_ASSERT_INT_NONNEGATIVE(footprint, NULL); 345 346 // Size of the images 347 int numCols = inMask->numCols; 348 int numRows = inMask->numRows; 349 350 // Worried about the masks for bad pixels and bad stamps colliding, so make our own mask 351 psImage *mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); // The global mask 352 psImageInit(mask, 0); 353 354 // Dereference for convenience 355 psMaskType **maskData = mask->data.PS_TYPE_MASK_DATA; 356 psMaskType **inData = inMask->data.PS_TYPE_MASK_DATA; 357 psMaskType **refData = refMask->data.PS_TYPE_MASK_DATA; 358 359 // Block out a border around the edge of the image 360 361 // Bottom stripe 362 for (int y = 0; y < PS_MIN(size + footprint, numRows); y++) { 363 for (int x = 0; x < numCols; x++) { 364 maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER; 365 } 366 } 367 // Either side 368 for (int y = PS_MIN(size + footprint, numRows); y < numRows - size - footprint; y++) { 369 for (int x = 0; x < PS_MIN(size + footprint, numCols); x++) { 370 maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER; 371 } 372 for (int x = PS_MAX(numCols - size - footprint, 0); x < numCols; x++) { 373 maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER; 374 } 375 } 376 // Top stripe 377 for (int y = PS_MAX(numRows - size - footprint, 0); y < numRows; y++) { 378 for (int x = 0; x < numCols; x++) { 379 maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER; 380 } 381 } 382 383 // Mask the bad pixels, and around them 384 for (int y = footprint + size; y < numRows - footprint - size; y++) { 385 for (int x = 0; x < numCols; x++) { 386 if (inData[y][x] & maskVal) { 387 maskData[y][x] |= PM_SUBTRACTION_MASK_INPUT; 388 // Block out the entire stamp footprint around this pixel 389 for (int v = y - footprint; v <= y + footprint; v++) { 390 for (int u = x - footprint; u <= x + footprint; u++) { 391 maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT; 392 } 393 } 394 } 395 if (refData[y][x] & maskVal) { 396 maskData[y][x] |= PM_SUBTRACTION_MASK_REF; 397 398 // We want to block out with the CONVOLVE mask anything that would be bad if we convolved with 399 // this bad pixel (within 'size'). Then we want to block out with the FOOTPRINT mask 400 // everything within a footprint's distance of those (within 'footprint'). 401 402 // Bottom stripe 403 for (int v = PS_MAX(y - footprint - size, 0); v < PS_MAX(y - size, 0); v++) { 404 for (int u = PS_MAX(x - footprint - size, 0); 405 u <= PS_MIN(x + footprint + size, numCols - 1); u++) { 406 maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT; 407 } 408 } 409 // Middle 410 for (int v = PS_MAX(y - size, 0); v <= PS_MIN(y + size, numRows - 1); v++) { 411 // Left side 412 for (int u = PS_MAX(x - footprint - size, 0); u < x - size; u++) { 413 maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT; 414 } 415 // Centre 416 for (int u = PS_MAX(x - size, 0); u <= PS_MIN(x + size, numCols - 1); u++) { 417 maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT | PM_SUBTRACTION_MASK_CONVOLVE; 418 } 419 // Right side 420 for (int u = PS_MIN(x + size + 1, numCols - 1); 421 u <= PS_MIN(x + footprint + size, numCols - 1); u++) { 422 maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT; 423 } 424 } 425 // Top stripe 426 for (int v = PS_MIN(y + size + 1, numRows - 1); 427 v <= PS_MIN(y + footprint + size, numRows - 1); v++) { 428 for (int u = PS_MAX(x - footprint - size, 0); 429 u <= PS_MIN(x + footprint + size, numCols - 1); u++) { 430 maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT; 431 } 432 } 433 } 434 } 435 } 436 437 return mask; 438 } 439 335 440 bool pmSubtractionCalculateEquation(psArray *stamps, const psImage *reference, const psImage *input, 336 441 const psImage *weight, const pmSubtractionKernels *kernels, int footprint) … … 550 655 551 656 int pmSubtractionRejectStamps(psArray *stamps, const psImage *refImage, psImage *inImage, 552 psImage * mask, psMaskType badStampMaskVal, const psVector *solution,553 int footprint, float sigmaRej,const pmSubtractionKernels *kernels)657 psImage *subMask, const psVector *solution, int footprint, float sigmaRej, 658 const pmSubtractionKernels *kernels) 554 659 { 555 660 PS_ASSERT_ARRAY_NON_NULL(stamps, -1); … … 559 664 PS_ASSERT_IMAGE_TYPE(inImage, PS_TYPE_F32, -1); 560 665 PS_ASSERT_IMAGES_SIZE_EQUAL(refImage, inImage, -1); 561 PS_ASSERT_IMAGE_NON_EMPTY( mask, -1);562 PS_ASSERT_IMAGE_TYPE( mask, PS_TYPE_MASK, -1);563 PS_ASSERT_IMAGES_SIZE_EQUAL(refImage, mask, -1);666 PS_ASSERT_IMAGE_NON_EMPTY(subMask, -1); 667 PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, -1); 668 PS_ASSERT_IMAGES_SIZE_EQUAL(refImage, subMask, -1); 564 669 PS_ASSERT_VECTOR_NON_NULL(solution, -1); 565 670 PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, -1); … … 654 759 for (int y = stamp->y - footprint; y <= stamp->y + footprint; y++) { 655 760 for (int x = stamp->x - footprint; x <= stamp->x + footprint; x++) { 656 mask->data.PS_TYPE_MASK_DATA[y][x] |= badStampMaskVal;761 subMask->data.PS_TYPE_MASK_DATA[y][x] |= PM_SUBTRACTION_MASK_REJ; 657 762 } 658 763 } … … 707 812 708 813 bool pmSubtractionConvolve(psImage **outImage, psImage **outWeight, psImage **outMask, 709 const psImage *inImage, const psImage *inWeight, const psImage *inMask, 710 psMaskType maskVal, psMaskType blank, 711 const psVector *solution, const pmSubtractionKernels *kernels) 814 const psImage *inImage, const psImage *inWeight, const psImage *subMask, 815 psMaskType blank, const psVector *solution, const pmSubtractionKernels *kernels) 712 816 { 713 817 PS_ASSERT_IMAGE_NON_NULL(inImage, false); 714 818 PS_ASSERT_IMAGE_TYPE(inImage, PS_TYPE_F32, false); 715 if ( inMask) {716 PS_ASSERT_IMAGE_NON_NULL( inMask, false);717 PS_ASSERT_IMAGE_TYPE( inMask, PS_TYPE_MASK, false);718 PS_ASSERT_IMAGES_SIZE_EQUAL( inMask, inImage, false);819 if (subMask) { 820 PS_ASSERT_IMAGE_NON_NULL(subMask, false); 821 PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, false); 822 PS_ASSERT_IMAGES_SIZE_EQUAL(subMask, inImage, false); 719 823 } 720 824 if (inWeight) { … … 733 837 PS_ASSERT_IMAGES_SIZE_EQUAL(*outMask, inImage, false); 734 838 PS_ASSERT_IMAGE_TYPE(*outMask, PS_TYPE_MASK, false); 735 PS_ASSERT_IMAGE_NON_NULL( inMask, false);839 PS_ASSERT_IMAGE_NON_NULL(subMask, false); 736 840 } 737 841 if (outWeight && *outWeight) { … … 763 867 } 764 868 psImage *convMask = NULL; // Convolved mask image 765 if (outMask && inMask) {869 if (outMask && subMask) { 766 870 if (!*outMask) { 767 871 *outMask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); … … 812 916 813 917 // Check and propagate the kernel footprint, if required 814 if ( inMask) {815 for (int v = -size; v <= size; v++) {816 for (int u = -size; u <= size; u++) {817 convMask->data.PS_TYPE_MASK_DATA[y][x] |=818 inMask->data.PS_TYPE_MASK_DATA[y + v][x + u];819 }918 if (subMask && (subMask->data.PS_TYPE_MASK_DATA[y][x] & 919 (PM_SUBTRACTION_MASK_INPUT | PM_SUBTRACTION_MASK_CONVOLVE))) { 920 convMask->data.PS_TYPE_MASK_DATA[y][x] |= blank; 921 convImage->data.F32[y][x] = NAN; 922 if (inWeight) { 923 convWeight->data.F32[y][x] = NAN; 820 924 } 821 if (convMask->data.PS_TYPE_MASK_DATA[y][x] & maskVal) { 822 convImage->data.F32[y][x] = NAN; 823 if (inWeight) { 824 convWeight->data.F32[y][x] = NAN; 825 } 826 continue; 827 } 925 continue; 828 926 } 829 927 -
trunk/psModules/src/imcombine/pmSubtraction.h
r13340 r13735 6 6 * @author GLG, MHPCC 7 7 * 8 * @version $Revision: 1. 1$ $Name: not supported by cvs2svn $9 * @date $Date: 2007-0 5-11 02:21:01$8 * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2007-06-09 01:04:02 $ 10 10 * Copyright 2004-207 Institute for Astronomy, University of Hawaii 11 11 */ … … 19 19 /// @addtogroup imcombine Image Combinations 20 20 /// @{ 21 22 typedef enum { 23 PM_SUBTRACTION_MASK_CLEAR = 0x00, // No masking 24 PM_SUBTRACTION_MASK_REF = 0x01, // Reference image is bad 25 PM_SUBTRACTION_MASK_INPUT = 0x02, // Input image is bad 26 PM_SUBTRACTION_MASK_CONVOLVE = 0x04, // If convolved, would be bad 27 PM_SUBTRACTION_MASK_FOOTPRINT = 0x08, // Bad pixel within the stamp footprint 28 PM_SUBTRACTION_MASK_BORDER = 0x10, // Image border 29 PM_SUBTRACTION_MASK_REJ = 0x20, // Previously tried as a stamp, and rejected 30 } pmSubtractionMasks; 31 32 /// Generate a mask for use in the subtraction process 33 psImage *pmSubtractionMask(const psImage *inMask, ///< Mask for the input image 34 const psImage *refMask, ///< Mask for the reference image (will be convolved) 35 psMaskType maskVal, ///< Value to mask out 36 int size, ///< Half-size of the kernel (pmSubtractionKernels.size) 37 int footprint ///< Half-size of the kernel footprint 38 ); 21 39 22 40 /// Calculate the least-squares equation to match the image quality … … 38 56 const psImage *refImage, ///< Reference image 39 57 psImage *inImage, ///< Input image 40 psImage *mask, ///< Mask image 41 psMaskType badStampMaskVal, ///< Value to use in mask for bad stamp 58 psImage *subMask, ///< Subtraction mask 42 59 const psVector *solution, ///< Solution vector 43 60 int footprint, ///< Region to mask if stamp is bad … … 58 75 const psImage *inImage, ///< Input image 59 76 const psImage *inWeight, ///< Input weight map (or NULL) 60 const psImage *inMask, ///< Input mask (or NULL) 61 psMaskType maskVal, ///< Value to mask 77 const psImage *subMask, ///< Subtraction mask (or NULL) 62 78 psMaskType blank, ///< Mask value for blank regions 63 79 const psVector *solution, ///< The solution vector -
trunk/psModules/src/imcombine/pmSubtractionStamps.c
r13510 r13735 6 6 #include <pslib.h> 7 7 8 #include "pmSubtraction.h" 8 9 #include "pmSubtractionStamps.h" 9 10 … … 39 40 } 40 41 41 42 psArray *pmSubtractionFindStamps(psArray *stamps, const psImage *image, const psImage *mask, 43 psMaskType maskVal, psMaskType used, float threshold, 44 float spacing, int border) 42 psArray *pmSubtractionFindStamps(psArray *stamps, const psImage *image, const psImage *subMask, 43 float threshold, float spacing) 45 44 { 46 45 PS_ASSERT_IMAGE_NON_NULL(image, NULL); 47 46 PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL); 48 if (mask) { 49 PS_ASSERT_IMAGE_NON_NULL(mask, NULL); 50 PS_ASSERT_IMAGES_SIZE_EQUAL(image, mask, NULL); 51 PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_MASK, NULL); 52 } 53 if (used == 0) { 54 psError(PS_ERR_BAD_PARAMETER_VALUE, true, "The mask value for used stamps cannot be zero."); 55 return NULL; 47 if (subMask) { 48 PS_ASSERT_IMAGE_NON_NULL(subMask, NULL); 49 PS_ASSERT_IMAGES_SIZE_EQUAL(image, subMask, NULL); 50 PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, NULL); 56 51 } 57 52 PS_ASSERT_FLOAT_LARGER_THAN(spacing, 0.0, NULL); 58 PS_ASSERT_INT_NONNEGATIVE(border, NULL);59 PS_ASSERT_INT_LARGER_THAN(image->numCols, 2 * border, NULL);60 PS_ASSERT_INT_LARGER_THAN(image->numRows, 2 * border, NULL);61 53 62 psMaskType badCentral = maskVal | used; // For central pixel of stamp, don't get a stamp already used 63 int xNumStamps = (image->numCols - 2 * border) / spacing + 1; // Number of stamps in x 64 int yNumStamps = (image->numRows - 2 * border) / spacing + 1; // Number of stamps in y 54 // Size of image 55 int numRows = image->numRows; 56 int numCols = image->numCols; 57 58 // Number of stamps 59 int xNumStamps = numCols / spacing + 1; 60 int yNumStamps = numRows / spacing + 1; 65 61 66 62 if (stamps) { … … 78 74 } 79 75 } 80 // Footprint of image81 int numRows = image->numRows;82 int numCols = image->numCols;83 76 84 77 for (int j = 0, index = 0; j < yNumStamps; j++) { … … 96 89 97 90 // Bounds of region to search for stamp 98 int yMin = border + j * (numRows - border)/ (yNumStamps + 1);99 int yMax = border + (j + 1) * (numRows - border)/ (yNumStamps + 1) - 1;100 int xMin = border + i * (numCols - border)/ (xNumStamps + 1);101 int xMax = border + (i + 1) * (numCols - border)/ (xNumStamps + 1) - 1;102 assert(yMax < image->numRows - border && xMax < image->numCols - border&&103 yMin >= border && xMin >= border);91 int yMin = j * numRows / (yNumStamps + 1); 92 int yMax = (j + 1) * numRows / (yNumStamps + 1) - 1; 93 int xMin = i * numCols / (xNumStamps + 1); 94 int xMax = (i + 1) * numCols / (xNumStamps + 1) - 1; 95 assert(yMax < image->numRows && xMax < image->numCols && 96 yMin >= 0 && xMin >= 0); 104 97 105 98 for (int y = yMin; y <= yMax ; y++) { 106 99 for (int x = xMin; x <= xMax ; x++) { 107 if (image->data.F32[y][x] > fluxBest) { 108 bool ok = true; 109 // Check kernel footprint for bad pixels 110 if (mask && (mask->data.PS_TYPE_MASK_DATA[y][x] & badCentral)) { 111 ok = false; 112 } else { 113 for (int v = -border; v <= border && ok; v++) { 114 for (int u = -border; u <= border && ok; u++) { 115 if ((mask && 116 (mask->data.PS_TYPE_MASK_DATA[y + v][x + u] & maskVal)) || 117 !isfinite(image->data.F32[y + v][x + u])) { 118 ok = false; 119 if (mask) { 120 // Mark it so we don't have to look so hard next time 121 mask->data.PS_TYPE_MASK_DATA[y][x] |= used; 122 } 123 } 124 } 125 } 126 } 127 128 if (ok) { 129 fluxBest = image->data.F32[y][x]; 130 xBest = x; 131 yBest = y; 132 } 100 if ((!subMask || !(subMask->data.PS_TYPE_MASK_DATA[y][x] & 101 (PM_SUBTRACTION_MASK_BORDER | PM_SUBTRACTION_MASK_FOOTPRINT | 102 PM_SUBTRACTION_MASK_REJ))) && 103 image->data.F32[y][x] > fluxBest) { 104 fluxBest = image->data.F32[y][x]; 105 xBest = x; 106 yBest = y; 133 107 } 134 108 } -
trunk/psModules/src/imcombine/pmSubtractionStamps.h
r13365 r13735 3 3 4 4 #include <pslib.h> 5 6 5 7 6 /// Status of stamp … … 26 25 const psImage *image, ///< Image for which to find stamps 27 26 const psImage *mask, ///< Mask 28 psMaskType maskVal, ///< Value for mask29 psMaskType bad, ///< Mask value for bad stamps30 27 float threshold, ///< Threshold for stamps in the image 31 float spacing, ///< Rough spacing for stamps 32 int border ///< Size of border around image; no stamps in border 28 float spacing ///< Rough spacing for stamps 33 29 ); 34 30
Note:
See TracChangeset
for help on using the changeset viewer.
