Changeset 14195
- Timestamp:
- Jul 13, 2007, 10:28:28 AM (19 years ago)
- Location:
- trunk/psModules/src/imcombine
- Files:
-
- 4 edited
-
pmSubtraction.c (modified) (10 diffs)
-
pmSubtraction.h (modified) (4 diffs)
-
pmSubtractionStamps.c (modified) (4 diffs)
-
pmSubtractionStamps.h (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/imcombine/pmSubtraction.c
r14125 r14195 4 4 * @author GLG, MHPCC 5 5 * 6 * @version $Revision: 1.1 6$ $Name: not supported by cvs2svn $7 * @date $Date: 2007-07-1 1 01:28:39$6 * @version $Revision: 1.17 $ $Name: not supported by cvs2svn $ 7 * @date $Date: 2007-07-13 20:28:28 $ 8 8 * 9 9 * Copyright 2004-2007 Institute for Astronomy, University of Hawaii … … 410 410 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 411 411 412 psImage *pmSubtractionMask(const psImage *inMask, const psImage *refMask, psMaskType maskVal, 413 int size, int footprint) 414 { 415 PS_ASSERT_IMAGE_NON_NULL(inMask, NULL); 416 PS_ASSERT_IMAGE_TYPE(inMask, PS_TYPE_MASK, NULL); 417 PS_ASSERT_IMAGE_NON_NULL(refMask, NULL); 418 PS_ASSERT_IMAGE_TYPE(refMask, PS_TYPE_MASK, NULL); 419 PS_ASSERT_IMAGES_SIZE_EQUAL(inMask, refMask, NULL); 420 PS_ASSERT_INT_NONNEGATIVE(size, NULL); 421 PS_ASSERT_INT_NONNEGATIVE(footprint, NULL); 422 423 // Size of the images 424 int numCols = inMask->numCols; 425 int numRows = inMask->numRows; 426 427 // Worried about the masks for bad pixels and bad stamps colliding, so make our own mask 428 psImage *mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); // The global mask 429 psImageInit(mask, 0); 430 431 // Dereference for convenience 432 psMaskType **maskData = mask->data.PS_TYPE_MASK_DATA; 433 psMaskType **inData = inMask->data.PS_TYPE_MASK_DATA; 434 psMaskType **refData = refMask->data.PS_TYPE_MASK_DATA; 435 436 // Block out a border around the edge of the image 437 438 // Bottom stripe 439 for (int y = 0; y < PS_MIN(size + footprint, numRows); y++) { 440 for (int x = 0; x < numCols; x++) { 441 maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER; 442 } 443 } 444 // Either side 445 for (int y = PS_MIN(size + footprint, numRows); y < numRows - size - footprint; y++) { 446 for (int x = 0; x < PS_MIN(size + footprint, numCols); x++) { 447 maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER; 448 } 449 for (int x = PS_MAX(numCols - size - footprint, 0); x < numCols; x++) { 450 maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER; 451 } 452 } 453 // Top stripe 454 for (int y = PS_MAX(numRows - size - footprint, 0); y < numRows; y++) { 455 for (int x = 0; x < numCols; x++) { 456 maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER; 457 } 458 } 459 460 // Mask the bad pixels, and around them 461 for (int y = footprint + size; y < numRows - footprint - size; y++) { 462 for (int x = 0; x < numCols; x++) { 463 if (inData[y][x] & maskVal) { 464 maskData[y][x] |= PM_SUBTRACTION_MASK_INPUT; 465 // Block out the entire stamp footprint around this pixel 466 for (int v = y - footprint; v <= y + footprint; v++) { 467 for (int u = x - footprint; u <= x + footprint; u++) { 468 maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT; 469 } 470 } 471 } 472 if (refData[y][x] & maskVal) { 473 maskData[y][x] |= PM_SUBTRACTION_MASK_REF; 474 475 // We want to block out with the CONVOLVE mask anything that would be bad if we convolved with 476 // this bad pixel (within 'size'). Then we want to block out with the FOOTPRINT mask 477 // everything within a footprint's distance of those (within 'footprint'). 478 479 // Bottom stripe 480 for (int v = PS_MAX(y - footprint - size, 0); v < y - size; v++) { 481 for (int u = PS_MAX(x - footprint - size, 0); 482 u <= PS_MIN(x + footprint + size, numCols - 1); u++) { 483 maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT; 484 } 485 } 486 // Middle 487 for (int v = PS_MAX(y - size, 0); v <= PS_MIN(y + size, numRows - 1); v++) { 488 // Left side 489 for (int u = PS_MAX(x - footprint - size, 0); u < x - size; u++) { 490 maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT; 491 } 492 // Centre 493 for (int u = PS_MAX(x - size, 0); u <= PS_MIN(x + size, numCols - 1); u++) { 494 maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT | PM_SUBTRACTION_MASK_CONVOLVE; 495 } 496 // Right side 497 for (int u = x + size + 1; u <= PS_MIN(x + footprint + size, numCols - 1); u++) { 498 maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT; 499 } 500 } 501 // Top stripe 502 for (int v = y + size + 1; v <= PS_MIN(y + footprint + size, numRows - 1); v++) { 503 for (int u = PS_MAX(x - footprint - size, 0); 504 u <= PS_MIN(x + footprint + size, numCols - 1); u++) { 505 maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT; 506 } 507 } 508 } 509 } 510 } 511 512 return mask; 513 } 514 515 412 516 bool pmSubtractionCalculateEquation(psArray *stamps, const psImage *reference, const psImage *input, 413 517 const psImage *weight, const pmSubtractionKernels *kernels, int footprint) … … 507 611 psFree(polyValues); 508 612 509 // Fill in lower diagonal of symmetric matrix 613 // Fill in lower diagonal of symmetric matrix, while checking for bad values 614 bool bad = false; // Are there bad values? 510 615 for (int i = 0; i < numKernels; i++) { 511 616 for (int j = 0; j < i; j++) { 512 617 stampMatrix->data.F64[i][j] = stampMatrix->data.F64[j][i]; 618 if (!isfinite(stampMatrix->data.F64[j][i])) { 619 bad = true; 620 } 513 621 } 514 622 stampMatrix->data.F64[bgIndex][i] = stampMatrix->data.F64[i][bgIndex]; 515 } 516 517 stamp->status = PM_SUBTRACTION_STAMP_USED; 623 if (!isfinite(stampMatrix->data.F64[i][bgIndex]) || 624 !isfinite(stampMatrix->data.F64[i][i]) || 625 !isfinite(stampVector->data.F64[i])) { 626 bad = true; 627 } 628 } 629 if (!isfinite(stampVector->data.F64[bgIndex])) { 630 bad = true; 631 } 632 633 if (bad) { 634 stamp->status = PM_SUBTRACTION_STAMP_REJECTED; 635 psTrace("psModules.imcombine", 3, "Rejecting stamp %d because of bad equation\n", i); 636 } else { 637 stamp->status = PM_SUBTRACTION_STAMP_USED; 638 } 518 639 519 640 if (psTraceGetLevel("psModules.imcombine.equation") >= 10) { … … 610 731 611 732 int pmSubtractionRejectStamps(psArray *stamps, const psImage *refImage, psImage *inImage, 612 psImage * mask, psMaskType badStampMaskVal, const psVector *solution,613 int footprint, float sigmaRej,const pmSubtractionKernels *kernels)733 psImage *subMask, const psVector *solution, int footprint, float sigmaRej, 734 const pmSubtractionKernels *kernels) 614 735 { 615 736 PS_ASSERT_ARRAY_NON_NULL(stamps, -1); … … 619 740 PS_ASSERT_IMAGE_TYPE(inImage, PS_TYPE_F32, -1); 620 741 PS_ASSERT_IMAGES_SIZE_EQUAL(refImage, inImage, -1); 621 PS_ASSERT_IMAGE_NON_EMPTY( mask, -1);622 PS_ASSERT_IMAGE_TYPE( mask, PS_TYPE_MASK, -1);623 PS_ASSERT_IMAGES_SIZE_EQUAL(refImage, mask, -1);742 PS_ASSERT_IMAGE_NON_EMPTY(subMask, -1); 743 PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, -1); 744 PS_ASSERT_IMAGES_SIZE_EQUAL(refImage, subMask, -1); 624 745 PS_ASSERT_VECTOR_NON_NULL(solution, -1); 625 746 PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, -1); … … 714 835 for (int y = stamp->y - footprint; y <= stamp->y + footprint; y++) { 715 836 for (int x = stamp->x - footprint; x <= stamp->x + footprint; x++) { 716 mask->data.PS_TYPE_MASK_DATA[y][x] |= badStampMaskVal;837 subMask->data.PS_TYPE_MASK_DATA[y][x] |= PM_SUBTRACTION_MASK_REJ; 717 838 } 718 839 } … … 767 888 768 889 bool pmSubtractionConvolve(psImage **outImage, psImage **outWeight, psImage **outMask, 769 const psImage *inImage, const psImage *inWeight, const psImage *inMask, 770 psMaskType maskVal, psMaskType blank, 771 const psVector *solution, const pmSubtractionKernels *kernels) 890 const psImage *inImage, const psImage *inWeight, const psImage *subMask, 891 psMaskType blank, const psVector *solution, const pmSubtractionKernels *kernels) 772 892 { 773 893 PS_ASSERT_IMAGE_NON_NULL(inImage, false); 774 894 PS_ASSERT_IMAGE_TYPE(inImage, PS_TYPE_F32, false); 775 if ( inMask) {776 PS_ASSERT_IMAGE_NON_NULL( inMask, false);777 PS_ASSERT_IMAGE_TYPE( inMask, PS_TYPE_MASK, false);778 PS_ASSERT_IMAGES_SIZE_EQUAL( inMask, inImage, false);895 if (subMask) { 896 PS_ASSERT_IMAGE_NON_NULL(subMask, false); 897 PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, false); 898 PS_ASSERT_IMAGES_SIZE_EQUAL(subMask, inImage, false); 779 899 } 780 900 if (inWeight) { … … 793 913 PS_ASSERT_IMAGES_SIZE_EQUAL(*outMask, inImage, false); 794 914 PS_ASSERT_IMAGE_TYPE(*outMask, PS_TYPE_MASK, false); 795 PS_ASSERT_IMAGE_NON_NULL( inMask, false);915 PS_ASSERT_IMAGE_NON_NULL(subMask, false); 796 916 } 797 917 if (outWeight && *outWeight) { … … 823 943 } 824 944 psImage *convMask = NULL; // Convolved mask image 825 if (outMask && inMask) {945 if (outMask && subMask) { 826 946 if (!*outMask) { 827 947 *outMask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); … … 872 992 873 993 // Check and propagate the kernel footprint, if required 874 if (inMask) { 875 for (int v = -size; v <= size; v++) { 876 for (int u = -size; u <= size; u++) { 877 convMask->data.PS_TYPE_MASK_DATA[y][x] |= 878 inMask->data.PS_TYPE_MASK_DATA[y + v][x + u]; 879 } 880 } 881 if (convMask->data.PS_TYPE_MASK_DATA[y][x] & maskVal) { 882 convImage->data.F32[y][x] = NAN; 883 if (inWeight) { 884 convWeight->data.F32[y][x] = NAN; 885 } 886 continue; 994 if (subMask && (subMask->data.PS_TYPE_MASK_DATA[y][x] & 995 (PM_SUBTRACTION_MASK_INPUT | PM_SUBTRACTION_MASK_CONVOLVE))) { 996 convMask->data.PS_TYPE_MASK_DATA[y][x] |= blank; 997 convImage->data.F32[y][x] = NAN; 998 if (inWeight) { 999 convWeight->data.F32[y][x] = NAN; 887 1000 } 888 1001 } -
trunk/psModules/src/imcombine/pmSubtraction.h
r14106 r14195 6 6 * @author GLG, MHPCC 7 7 * 8 * @version $Revision: 1. 3$ $Name: not supported by cvs2svn $9 * @date $Date: 2007-07-1 0 23:54:26$8 * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2007-07-13 20:28:28 $ 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
r14106 r14195 6 6 #include <pslib.h> 7 7 8 #include "pmSubtraction.h" 8 9 #include "pmSubtractionStamps.h" 9 10 … … 40 41 41 42 42 psArray *pmSubtractionFindStamps(psArray *stamps, const psImage *image, const psImage *mask, 43 psMaskType maskVal, psMaskType used, float threshold, 44 float spacing, int border) 43 psArray *pmSubtractionFindStamps(psArray *stamps, const psImage *image, const psImage *subMask, 44 float threshold, float spacing) 45 45 { 46 46 PS_ASSERT_IMAGE_NON_NULL(image, NULL); 47 47 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; 48 if (subMask) { 49 PS_ASSERT_IMAGE_NON_NULL(subMask, NULL); 50 PS_ASSERT_IMAGES_SIZE_EQUAL(image, subMask, NULL); 51 PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, NULL); 56 52 } 57 53 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 54 62 maskVal |= used; // Make sure we don't get stamps we've 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 55 // Size of image 56 int numRows = image->numRows; 57 int numCols = image->numCols; 58 59 // Number of stamps 60 int xNumStamps = numCols / spacing + 1; 61 int yNumStamps = numRows / spacing + 1; 65 62 66 63 if (stamps) { … … 78 75 } 79 76 } 80 // Footprint of image81 int numRows = image->numRows;82 int numCols = image->numCols;83 77 84 78 for (int j = 0, index = 0; j < yNumStamps; j++) { … … 96 90 97 91 // 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); 92 int yMin = j * numRows / (yNumStamps + 1); 93 int yMax = (j + 1) * numRows / (yNumStamps + 1) - 1; 94 int xMin = i * numCols / (xNumStamps + 1); 95 int xMax = (i + 1) * numCols / (xNumStamps + 1) - 1; 96 assert(yMax < image->numRows && xMax < image->numCols && 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 if (mask) { 110 // Check kernel footprint for bad pixels 111 if (mask->data.PS_TYPE_MASK_DATA[y][x] & maskVal) { 112 ok = false; 113 } else { 114 for (int v = -border; v <= border && ok; v++) { 115 for (int u = -border; u <= border && ok; u++) { 116 if (mask->data.PS_TYPE_MASK_DATA[y + v][x + u] & maskVal) { 117 ok = false; 118 #if 0 119 // Mark it so we don't have to look so hard next time 120 mask->data.PS_TYPE_MASK_DATA[y][x] |= maskVal; 121 #endif 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
r14106 r14195 26 26 const psImage *image, ///< Image for which to find stamps 27 27 const psImage *mask, ///< Mask 28 psMaskType maskVal, ///< Value for mask29 psMaskType bad, ///< Mask value for bad stamps30 28 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 29 float spacing ///< Rough spacing for stamps 33 30 ); 34 31
Note:
See TracChangeset
for help on using the changeset viewer.
