Changeset 14106
- Timestamp:
- Jul 10, 2007, 1:54:26 PM (19 years ago)
- Location:
- trunk/psModules/src/imcombine
- Files:
-
- 6 edited
-
pmSubtraction.c (modified) (28 diffs)
-
pmSubtraction.h (modified) (4 diffs)
-
pmSubtractionKernels.c (modified) (4 diffs)
-
pmSubtractionKernels.h (modified) (3 diffs)
-
pmSubtractionStamps.c (modified) (4 diffs)
-
pmSubtractionStamps.h (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/imcombine/pmSubtraction.c
r13736 r14106 4 4 * @author GLG, MHPCC 5 5 * 6 * @version $Revision: 1.1 4$ $Name: not supported by cvs2svn $7 * @date $Date: 2007-0 6-09 01:06:48$6 * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $ 7 * @date $Date: 2007-07-10 23:54:26 $ 8 8 * 9 9 * Copyright 2004-2007 Institute for Astronomy, University of Hawaii … … 24 24 25 25 #include "pmSubtraction.h" 26 27 //#define CAN_PASS_STATIC_INLINE // Can I pass the address of a static inline function? 26 28 27 29 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// … … 53 55 } 54 56 57 #ifndef CAN_PASS_STATIC_INLINE 55 58 // Generate an image of the solved kernel 56 59 // Meant to replace the following function declaration: … … 61 64 // double (*weightFunc)(double value) // Function for weighting 62 65 // ); 63 // because even if the weightFunc can't be "static inline", they still appearin the assembly.66 // because even though the weightFunc is "static inline", it still appears in the assembly. 64 67 // TARGET: Kernel, to return (psKernel*) 65 68 // SOLUTION: Solution to the least-squares problem (psVector*) … … 75 78 \ 76 79 /* Ensure the subIndex for POIS kernels is what is expected */ \ 77 assert((KERNELS)->type != PM_SUBTRACTION_KERNEL_POIS || \ 80 assert(((KERNELS)->type != PM_SUBTRACTION_KERNEL_POIS && \ 81 (KERNELS)->type != PM_SUBTRACTION_KERNEL_SPAM && \ 82 (KERNELS)->type != PM_SUBTRACTION_KERNEL_FRIES) || \ 78 83 ((KERNELS)->u->data.S32[(KERNELS)->subIndex] == 0 && \ 79 84 (KERNELS)->v->data.S32[(KERNELS)->subIndex] == 0 && \ … … 103 108 break; \ 104 109 } \ 110 /* SPAM and FRIES use the same method */ \ 111 case PM_SUBTRACTION_KERNEL_SPAM: \ 112 case PM_SUBTRACTION_KERNEL_FRIES: { \ 113 int uStart = (KERNELS)->u->data.S32[i]; \ 114 int uStop = (KERNELS)->uStop->data.S32[i]; \ 115 int vStart = (KERNELS)->v->data.S32[i]; \ 116 int vStop = (KERNELS)->vStop->data.S32[i]; \ 117 for (int v = vStart; v <= vStop; v++) { \ 118 for (int u = uStart; u <= uStop; u++) { \ 119 (TARGET)->kernel[v][u] += value; \ 120 } \ 121 } \ 122 if ((KERNELS)->spatialOrder > 0 && i != (KERNELS)->subIndex) { \ 123 /* The (0,0) element is subtracted from most kernels to preserve photometric scaling */ \ 124 (TARGET)->kernel[0][0] += subValue; \ 125 } \ 126 break; \ 127 } \ 105 128 case PM_SUBTRACTION_KERNEL_ISIS: { \ 106 129 psKernel *preCalc = (KERNELS)->preCalc->data[i]; /* Precalculated values */ \ … … 122 145 } \ 123 146 } 124 125 #if 0 147 #else 126 148 // Generate an image of the solved kernel 127 149 static inline psKernel *solvedKernel(psKernel *kernel, // Kernel, to return … … 139 161 140 162 // Ensure the subIndex for POIS kernels is what is expected 141 assert(kernels->type != PM_SUBTRACTION_KERNEL_POIS || 163 assert((kernels->type != PM_SUBTRACTION_KERNEL_POIS && 164 kernels->type != PM_SUBTRACTION_KERNEL_SPAM && 165 kernels->type != PM_SUBTRACTION_KERNEL_FRIES) || 142 166 (kernels->u->data.S32[kernels->subIndex] == 0 && kernels->v->data.S32[kernels->subIndex] == 0 && 143 167 kernels->xOrder->data.S32[kernels->subIndex] == 0 && … … 167 191 break; 168 192 } 193 /* SPAM and FRIES use the same method */ 194 case PM_SUBTRACTION_KERNEL_SPAM: 195 case PM_SUBTRACTION_KERNEL_FRIES: 196 int uStart = kernels->u->data.S32[i]; 197 int uStop = kernels->uStop->data.S32[i]; 198 int vStart = kernels->v->data.S32[i]; 199 int vStop = kernels->vStop->data.S32[i]; 200 for (int v = vStart; v <= vStop; v++) { 201 for (int u = uStart; u <= uStop; u++) { 202 kernel->kernel[v][u] += value; 203 } 204 } 205 if (kernels->spatialOrder > 0 && i != kernels->subIndex) { 206 /* The (0,0) element is subtracted from most kernels to preserve photometric scaling */ 207 kernel->kernel[0][0] += subValue; 208 } 209 break; 169 210 case PM_SUBTRACTION_KERNEL_ISIS: { 170 211 psKernel *preCalc = kernels->preCalc->data[i];// Precalculated values … … 190 231 #endif 191 232 233 #ifndef CAN_PASS_STATIC_INLINE 192 234 // Generate the convolved pixel value 193 235 // Meant to replace the following function declaration: … … 199 241 // double (*weightFunc)(double value) // Function for weighting 200 242 // ); 201 // because even if the weightFunc can't be "static inline", they still appearin the assembly.243 // because even though weightFunc is "static inline", it still appears in the assembly. 202 244 // 203 245 // TARGET: The value to 'return' (double) … … 225 267 break; \ 226 268 } \ 227 case PM_SUBTRACTION_KERNEL_ISIS: { \ 269 /* Method for SPAM and FRIES is the same */ \ 270 case PM_SUBTRACTION_KERNEL_SPAM: \ 271 case PM_SUBTRACTION_KERNEL_FRIES: { \ 272 int uStart = (KERNELS)->u->data.S32[(INDEX)]; \ 273 int uStop = (KERNELS)->uStop->data.S32[(INDEX)]; \ 274 int vStart = (KERNELS)->v->data.S32[(INDEX)]; \ 275 int vStop = (KERNELS)->vStop->data.S32[(INDEX)]; \ 276 for (int v = vStart; v <= vStop; v++) { \ 277 for (int u = uStart; u <= uStop; u++) { \ 278 (TARGET) += FUNC(polyValue) * (IMAGE)->data.F32[(Y) + v][(X) + u]; \ 279 } \ 280 } \ 281 if ((KERNELS)->spatialOrder > 0 && (INDEX) != (KERNELS)->subIndex) { \ 282 /* The (0,0) element is subtracted from most kernels to preserve photometric scaling */ \ 283 (TARGET) += FUNC(-1.0) * (IMAGE)->data.F32[(Y)][(X)]; \ 284 } \ 285 break; \ 286 } \ 287 case PM_SUBTRACTION_KERNEL_ISIS: { \ 228 288 psKernel *kernel = (KERNELS)->preCalc->data[(INDEX)]; /* The convolution kernel */ \ 229 289 psKernel *subKernel = (KERNELS)->preCalc->data[(KERNELS)->subIndex]; /* Kernel to subtract */ \ … … 231 291 double sum = 0.0; /* Accumulated sum from convolution */ \ 232 292 double sub = 0.0; /* Accumulated sum to subtract */ \ 233 for (int v = -size ; v <= size; v++) { \234 for (int u = -size ; u <= size; u++) { \235 sum += FUNC(kernel->kernel[v][u]) * (IMAGE)->data.F32[ (Y) + v][(X) + u]; \293 for (int v = -size, y = (Y) - size; v <= size; v++, y++) { \ 294 for (int u = -size, x = (X) - size; u <= size; u++, x++) { \ 295 sum += FUNC(kernel->kernel[v][u]) * (IMAGE)->data.F32[y][x]; \ 236 296 /* The (0,0) kernel is subtracted from other kernels to preserve photometric scaling */ \ 237 297 if ((KERNELS)->spatialOrder > 0 && (INDEX) != (KERNELS)->subIndex) { \ 238 sub += FUNC(subKernel->kernel[v][u]) * (IMAGE)->data.F32[ (Y) + v][(X) + u]; \298 sub += FUNC(subKernel->kernel[v][u]) * (IMAGE)->data.F32[y][x]; \ 239 299 } \ 240 300 } \ … … 247 307 } \ 248 308 } 249 250 #if 0 309 #else 251 310 // Generate the convolved pixel value 252 311 static inline double convolvePixel(const pmSubtractionKernels *kernels, // Kernel basis functions … … 273 332 } 274 333 return value; 334 } 335 /* Method for SPAM and FRIES is the same */ 336 case PM_SUBTRACTION_KERNEL_SPAM: 337 case PM_SUBTRACTION_KERNEL_FRIES: { 338 int uStart = kernels->u->data.S32[index]; 339 int uStop = kernels->uStop->data.S32[index]; 340 int vStart = kernels->v->data.S32[index]; 341 int vStop = kernels->vStop->data.S32[index]; 342 for (int v = vStart; v <= vStop; v++) { 343 for (int u = uStart; u <= uStop; u++) { 344 value += weightFunc(polyValue) * image->data.F32[y + v][x + u]; 345 } 346 } 347 if (kernels->spatialOrder > 0 && index != kernels->subIndex) { 348 /* The (0,0) element is subtracted from most kernels to preserve photometric scaling */ 349 value += weightFunc(-1.0) * image->data.F32[y][x]; 350 } 351 break; 275 352 } 276 353 case PM_SUBTRACTION_KERNEL_ISIS: { … … 333 410 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 334 411 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 images347 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 mask351 psImage *mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); // The global mask352 psImageInit(mask, 0);353 354 // Dereference for convenience355 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 image360 361 // Bottom stripe362 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 side368 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 stripe377 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 them384 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 pixel389 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 with399 // this bad pixel (within 'size'). Then we want to block out with the FOOTPRINT mask400 // everything within a footprint's distance of those (within 'footprint').401 402 // Bottom stripe403 for (int v = PS_MAX(y - footprint - size, 0); v < y - size; 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 // Middle410 for (int v = PS_MAX(y - size, 0); v <= PS_MIN(y + size, numRows - 1); v++) {411 // Left side412 for (int u = PS_MAX(x - footprint - size, 0); u < x - size; u++) {413 maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT;414 }415 // Centre416 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 side420 for (int u = x + size + 1; u <= PS_MIN(x + footprint + size, numCols - 1); u++) {421 maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT;422 }423 }424 // Top stripe425 for (int v = y + size + 1; v <= PS_MIN(y + footprint + size, numRows - 1); v++) {426 for (int u = PS_MAX(x - footprint - size, 0);427 u <= PS_MIN(x + footprint + size, numCols - 1); u++) {428 maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT;429 }430 }431 }432 }433 }434 435 return mask;436 }437 438 412 bool pmSubtractionCalculateEquation(psArray *stamps, const psImage *reference, const psImage *input, 439 413 const psImage *weight, const pmSubtractionKernels *kernels, int footprint) … … 504 478 // Generate the convolutions 505 479 for (int i = 0; i < numKernels; i++) { 506 #if 0480 #ifdef CAN_PASS_STATIC_INLINE 507 481 convolutions->data.F64[i] = convolvePixel(kernels, i, x, y, reference, polyValues, 508 482 imageWeighting); … … 533 507 psFree(polyValues); 534 508 535 // Fill in lower diagonal of symmetric matrix, while checking for bad values 536 bool bad = false; // Are there bad values? 509 // Fill in lower diagonal of symmetric matrix 537 510 for (int i = 0; i < numKernels; i++) { 538 511 for (int j = 0; j < i; j++) { 539 if (!isfinite(stampMatrix->data.F64[j][i])) {540 bad = true;541 }542 512 stampMatrix->data.F64[i][j] = stampMatrix->data.F64[j][i]; 543 513 } 544 if (!isfinite(stampMatrix->data.F64[i][bgIndex]) ||545 !isfinite(stampMatrix->data.F64[i][i]) ||546 !isfinite(stampVector->data.F64[i])) {547 bad = true;548 }549 514 stampMatrix->data.F64[bgIndex][i] = stampMatrix->data.F64[i][bgIndex]; 550 515 } 551 if (!isfinite(stampVector->data.F64[bgIndex])) { 552 bad = true; 553 } 554 555 if (bad) { 556 stamp->status = PM_SUBTRACTION_STAMP_REJECTED; 557 psTrace("psModules.imcombine", 3, "Rejecting stamp %d because of bad equation\n", i); 558 } else { 559 stamp->status = PM_SUBTRACTION_STAMP_USED; 560 } 516 517 stamp->status = PM_SUBTRACTION_STAMP_USED; 561 518 562 519 if (psTraceGetLevel("psModules.imcombine.equation") >= 10) { … … 653 610 654 611 int pmSubtractionRejectStamps(psArray *stamps, const psImage *refImage, psImage *inImage, 655 psImage * subMask, const psVector *solution, int footprint, float sigmaRej,656 const pmSubtractionKernels *kernels)612 psImage *mask, psMaskType badStampMaskVal, const psVector *solution, 613 int footprint, float sigmaRej, const pmSubtractionKernels *kernels) 657 614 { 658 615 PS_ASSERT_ARRAY_NON_NULL(stamps, -1); … … 662 619 PS_ASSERT_IMAGE_TYPE(inImage, PS_TYPE_F32, -1); 663 620 PS_ASSERT_IMAGES_SIZE_EQUAL(refImage, inImage, -1); 664 PS_ASSERT_IMAGE_NON_EMPTY( subMask, -1);665 PS_ASSERT_IMAGE_TYPE( subMask, PS_TYPE_MASK, -1);666 PS_ASSERT_IMAGES_SIZE_EQUAL(refImage, subMask, -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); 667 624 PS_ASSERT_VECTOR_NON_NULL(solution, -1); 668 625 PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, -1); … … 707 664 for (int x = footprint; x <= footprint; x++) { 708 665 for (int k = 0; k < numKernels; k++) { 709 #if 0666 #ifdef CAN_PASS_STATIC_INLINE 710 667 convolvedStamp->data.F32[y + footprint][x + footprint] += 711 668 convolvePixel(kernels, k, xStamp + x, yStamp + y, refImage, polyValues, … … 753 710 if (stamp->status == PM_SUBTRACTION_STAMP_USED && deviations->data.F32[i] > limit) { 754 711 // Mask out the stamp in the image so you it's not found again 755 psTrace("psModules.imcombine", 3, " Clipping stamp %d\n", i);712 psTrace("psModules.imcombine", 3, "Rejecting stamp %d\n", i); 756 713 numRejected++; 757 714 for (int y = stamp->y - footprint; y <= stamp->y + footprint; y++) { 758 715 for (int x = stamp->x - footprint; x <= stamp->x + footprint; x++) { 759 subMask->data.PS_TYPE_MASK_DATA[y][x] |= PM_SUBTRACTION_MASK_REJ;716 mask->data.PS_TYPE_MASK_DATA[y][x] |= badStampMaskVal; 760 717 } 761 718 } … … 794 751 795 752 // The appropriate kernel 796 #if 0753 #ifdef CAN_PASS_STATIC_INLINE 797 754 psKernel *kernel = solvedKernel(NULL, solution, kernels, polyValues, imageWeighting); 798 755 #else … … 810 767 811 768 bool pmSubtractionConvolve(psImage **outImage, psImage **outWeight, psImage **outMask, 812 const psImage *inImage, const psImage *inWeight, const psImage *subMask, 813 psMaskType blank, const psVector *solution, const pmSubtractionKernels *kernels) 769 const psImage *inImage, const psImage *inWeight, const psImage *inMask, 770 psMaskType maskVal, psMaskType blank, 771 const psVector *solution, const pmSubtractionKernels *kernels) 814 772 { 815 773 PS_ASSERT_IMAGE_NON_NULL(inImage, false); 816 774 PS_ASSERT_IMAGE_TYPE(inImage, PS_TYPE_F32, false); 817 if ( subMask) {818 PS_ASSERT_IMAGE_NON_NULL( subMask, false);819 PS_ASSERT_IMAGE_TYPE( subMask, PS_TYPE_MASK, false);820 PS_ASSERT_IMAGES_SIZE_EQUAL( subMask, inImage, 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); 821 779 } 822 780 if (inWeight) { … … 835 793 PS_ASSERT_IMAGES_SIZE_EQUAL(*outMask, inImage, false); 836 794 PS_ASSERT_IMAGE_TYPE(*outMask, PS_TYPE_MASK, false); 837 PS_ASSERT_IMAGE_NON_NULL( subMask, false);795 PS_ASSERT_IMAGE_NON_NULL(inMask, false); 838 796 } 839 797 if (outWeight && *outWeight) { … … 865 823 } 866 824 psImage *convMask = NULL; // Convolved mask image 867 if (outMask && subMask) {825 if (outMask && inMask) { 868 826 if (!*outMask) { 869 827 *outMask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); … … 896 854 2.0 * (float)(i + size - numCols/2.0) / (float)numCols, 897 855 2.0 * (float)(j + size - numRows/2.0) / (float)numRows); 898 #if 0856 #ifdef CAN_PASS_STATIC_INLINE 899 857 kernelImage = solvedKernel(kernelImage, solution, kernels, polyValues, imageWeighting); 900 858 if (inWeight) { … … 914 872 915 873 // Check and propagate the kernel footprint, if required 916 if ( subMask && (subMask->data.PS_TYPE_MASK_DATA[y][x] &917 (PM_SUBTRACTION_MASK_INPUT | PM_SUBTRACTION_MASK_CONVOLVE))) {918 convMask->data.PS_TYPE_MASK_DATA[y][x] |= blank;919 convImage->data.F32[y][x] = NAN;920 if (inWeight) {921 convWeight->data.F32[y][x] = NAN;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 } 922 880 } 923 continue; 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; 887 } 924 888 } 925 889 -
trunk/psModules/src/imcombine/pmSubtraction.h
r13735 r14106 6 6 * @author GLG, MHPCC 7 7 * 8 * @version $Revision: 1. 2$ $Name: not supported by cvs2svn $9 * @date $Date: 2007-0 6-09 01:04:02$8 * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2007-07-10 23:54:26 $ 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 masking24 PM_SUBTRACTION_MASK_REF = 0x01, // Reference image is bad25 PM_SUBTRACTION_MASK_INPUT = 0x02, // Input image is bad26 PM_SUBTRACTION_MASK_CONVOLVE = 0x04, // If convolved, would be bad27 PM_SUBTRACTION_MASK_FOOTPRINT = 0x08, // Bad pixel within the stamp footprint28 PM_SUBTRACTION_MASK_BORDER = 0x10, // Image border29 PM_SUBTRACTION_MASK_REJ = 0x20, // Previously tried as a stamp, and rejected30 } pmSubtractionMasks;31 32 /// Generate a mask for use in the subtraction process33 psImage *pmSubtractionMask(const psImage *inMask, ///< Mask for the input image34 const psImage *refMask, ///< Mask for the reference image (will be convolved)35 psMaskType maskVal, ///< Value to mask out36 int size, ///< Half-size of the kernel (pmSubtractionKernels.size)37 int footprint ///< Half-size of the kernel footprint38 );39 21 40 22 /// Calculate the least-squares equation to match the image quality … … 56 38 const psImage *refImage, ///< Reference image 57 39 psImage *inImage, ///< Input image 58 psImage *subMask, ///< Subtraction mask 40 psImage *mask, ///< Mask image 41 psMaskType badStampMaskVal, ///< Value to use in mask for bad stamp 59 42 const psVector *solution, ///< Solution vector 60 43 int footprint, ///< Region to mask if stamp is bad … … 75 58 const psImage *inImage, ///< Input image 76 59 const psImage *inWeight, ///< Input weight map (or NULL) 77 const psImage *subMask, ///< Subtraction mask (or NULL) 60 const psImage *inMask, ///< Input mask (or NULL) 61 psMaskType maskVal, ///< Value to mask 78 62 psMaskType blank, ///< Mask value for blank regions 79 63 const psVector *solution, ///< The solution vector -
trunk/psModules/src/imcombine/pmSubtractionKernels.c
r13390 r14106 15 15 psFree(kernels->v); 16 16 psFree(kernels->sigma); 17 psFree(kernels->uStop); 18 psFree(kernels->vStop); 17 19 psFree(kernels->xOrder); 18 20 psFree(kernels->yOrder); … … 48 50 kernels->u = psVectorAlloc(numBasisFunctions, PS_TYPE_S32); 49 51 kernels->v = psVectorAlloc(numBasisFunctions, PS_TYPE_S32); 52 kernels->sigma = NULL; 53 kernels->uStop = NULL; 54 kernels->vStop = NULL; 50 55 kernels->xOrder = psVectorAlloc(numBasisFunctions, PS_TYPE_S32); 51 56 kernels->yOrder = psVectorAlloc(numBasisFunctions, PS_TYPE_S32); 52 kernels->sigma = NULL;53 57 kernels->subIndex = 0; 54 58 kernels->preCalc = NULL; … … 187 191 } 188 192 193 /// Generate SPAM kernels 194 pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, ///< Half-size of the kernel 195 int spatialOrder, ///< Order of spatial variations 196 int inner, ///< Inner radius to preserve unbinned 197 int binning ///< Kernel binning factor 198 ) 199 { 200 PS_ASSERT_INT_POSITIVE(size, NULL); 201 PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL); 202 PS_ASSERT_INT_NONNEGATIVE(inner, NULL); 203 PS_ASSERT_INT_LARGER_THAN(size, inner, NULL); 204 PS_ASSERT_INT_POSITIVE(binning, NULL); 205 206 // The outer region should be divisible by the "binning"; otherwise allocate remainder to the inner region 207 int numOuter = (size - inner) / binning; // Number of summed pixels in the outer region 208 int numInner = inner + (size - inner) % binning; // Number of pixels in the inner region 209 assert(numOuter * binning + numInner == size); 210 int numTotal = numOuter + numInner; // Total number of summed pixels 211 212 psTrace("psModules.imcombine", 3, "Inner: %d Outer: %d\n", numInner, numOuter); 213 214 int num = PS_SQR(2 * numTotal + 1) * 215 (spatialOrder + 1) * (spatialOrder + 2) / 2; // Number of basis functions 216 217 psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num); 218 219 pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM, 220 size, spatialOrder); // The kernels 221 222 kernels->uStop = psVectorAlloc(num, PS_TYPE_F32); 223 kernels->vStop = psVectorAlloc(num, PS_TYPE_F32); 224 225 psVector *locations = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); // Locations for each kernel element 226 psVector *widths = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); // Widths for each kernel element 227 locations->data.S32[numTotal] = 0; 228 widths->data.S32[numTotal] = 0; 229 for (int i = 1; i <= numInner; i++) { 230 locations->data.S32[numTotal + i] = i; 231 widths->data.S32[numTotal + i] = 0; 232 locations->data.S32[numTotal - i] = - i; 233 widths->data.S32[numTotal - i] = 0; 234 } 235 for (int i = numInner + 1; i <= numTotal; i++) { 236 locations->data.S32[numTotal + i] = locations->data.S32[numTotal + i - 1] + 237 widths->data.S32[numTotal + i - 1] + 1; 238 widths->data.S32[numTotal + i] = binning - 1; 239 locations->data.S32[numTotal - i] = locations->data.S32[numTotal - i + 1] - 240 widths->data.S32[numTotal - i + 1] - binning; 241 widths->data.S32[numTotal - i] = binning - 1; 242 } 243 244 if (psTraceGetLevel("psModules.imcombine") >= 10) { 245 for (int i = 0; i < 2 * numTotal + 1; i++) { 246 psTrace("psModules.imcombine", 10, "%d: %d -> %d\n", i, locations->data.S32[i], 247 locations->data.S32[i] + widths->data.S32[i]); 248 } 249 } 250 251 // Set the kernel parameters 252 for (int i = - numTotal, index = 0; i <= numTotal; i++) { 253 int u = locations->data.S32[numTotal + i]; // Location of pixel 254 int uStop = u + widths->data.S32[numTotal + i]; // Width of pixel 255 256 for (int j = - numTotal; j <= numTotal; j++) { 257 int v = locations->data.S32[numTotal + j]; // Location of pixel 258 int vStop = v + widths->data.S32[numTotal + j]; // Width of pixel 259 260 // Iterate over spatial order. This loop creates the terms for 261 // x^xOrder * y^yOrder such that (xOrder+yOrder) <= spatialOrder. 262 for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) { 263 for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) { 264 kernels->u->data.S32[index] = u; 265 kernels->v->data.S32[index] = v; 266 kernels->uStop->data.S32[index] = uStop; 267 kernels->vStop->data.S32[index] = vStop; 268 kernels->xOrder->data.S32[index] = xOrder; 269 kernels->yOrder->data.S32[index] = yOrder; 270 271 psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d %d %d\n", index, 272 u, uStop, v, vStop, xOrder, yOrder); 273 } 274 } 275 } 276 } 277 278 kernels->subIndex = (num - (spatialOrder + 1) * (spatialOrder + 2) / 2) / 2; 279 assert(kernels->u->data.S32[kernels->subIndex] == 0 && 280 kernels->v->data.S32[kernels->subIndex] == 0 && 281 kernels->uStop->data.S32[kernels->subIndex] == 0 && 282 kernels->vStop->data.S32[kernels->subIndex] == 0 && 283 kernels->xOrder->data.S32[kernels->subIndex] == 0 && 284 kernels->yOrder->data.S32[kernels->subIndex] == 0); 285 286 psFree(locations); 287 psFree(widths); 288 289 return kernels; 290 } 291 292 293 /// Generate FRIES kernels 294 pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, ///< Half-size of the kernel 295 int spatialOrder, ///< Order of spatial variations 296 int inner ///< Inner radius to preserve unbinned 297 ) 298 { 299 PS_ASSERT_INT_POSITIVE(size, NULL); 300 PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL); 301 PS_ASSERT_INT_NONNEGATIVE(inner, NULL); 302 PS_ASSERT_INT_LARGER_THAN(size, inner, NULL); 303 304 int fibNum = 0; // Number of Fibonacci values 305 int fibLast = 1, fibTotal = 2; // Fibonacci sequence 306 while (fibTotal < size - inner) { 307 int temp = fibTotal; 308 fibTotal += fibLast; 309 fibLast = temp; 310 fibNum++; 311 } 312 313 int numInner = inner; // Number of pixels in the inner region 314 int numOuter = fibNum; // Number of summed pixels in the outer region 315 int numTotal = numOuter + numInner; // Total number of summed pixels 316 317 psTrace("psModules.imcombine", 3, "Inner: %d Outer: %d\n", numInner, numOuter); 318 319 int num = PS_SQR(2 * numTotal + 1) * 320 (spatialOrder + 1) * (spatialOrder + 2) / 2; // Number of basis functions 321 322 psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num); 323 324 pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM, 325 size, spatialOrder); // The kernels 326 kernels->uStop = psVectorAlloc(num, PS_TYPE_F32); 327 kernels->vStop = psVectorAlloc(num, PS_TYPE_F32); 328 329 psVector *start = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); 330 psVector *stop = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); 331 start->data.S32[numTotal] = 0; 332 stop->data.S32[numTotal] = 0; 333 for (int i = 1; i <= numInner; i++) { 334 start->data.S32[numTotal + i] = i; 335 stop->data.S32[numTotal + i] = i; 336 start->data.S32[numTotal - i] = -i; 337 stop->data.S32[numTotal - i] = -i; 338 } 339 for (int i = numInner + 1, fibLast = 1, fib = 2, temp; i <= numTotal; 340 i++, fib = (temp = fib) + fibLast, fibLast = temp) { 341 start->data.S32[numTotal + i] = stop->data.S32[numTotal + i - 1] + 1; 342 stop->data.S32[numTotal + i] = PS_MIN(start->data.S32[numTotal + i] + fib - 1, size); 343 start->data.S32[numTotal - i] = - stop->data.S32[numTotal + i]; 344 stop->data.S32[numTotal - i] = - start->data.S32[numTotal + i]; 345 } 346 347 if (psTraceGetLevel("psModules.imcombine") >= 10) { 348 for (int i = 0; i < 2 * numTotal + 1; i++) { 349 psTrace("psModules.imcombine", 10, "%d: %d -> %d\n", i, start->data.S32[i], stop->data.S32[i]); 350 } 351 } 352 353 // Set the kernel parameters 354 for (int i = - numTotal, index = 0; i <= numTotal; i++) { 355 int u = start->data.S32[numTotal + i]; // Location of pixel 356 int uStop = stop->data.S32[numTotal + i]; // Width of pixel 357 for (int j = - numTotal; j <= numTotal; j++) { 358 int v = start->data.S32[numTotal + j]; // Location of pixel 359 int vStop = stop->data.S32[numTotal + j]; // Width of pixel 360 361 // Iterate over spatial order. This loop creates the terms for 362 // x^xOrder * y^yOrder such that (xOrder+yOrder) <= spatialOrder. 363 for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) { 364 for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) { 365 kernels->u->data.S32[index] = u; 366 kernels->v->data.S32[index] = v; 367 kernels->uStop->data.S32[index] = uStop; 368 kernels->vStop->data.S32[index] = vStop; 369 kernels->xOrder->data.S32[index] = xOrder; 370 kernels->yOrder->data.S32[index] = yOrder; 371 372 psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d %d %d\n", index, 373 u, uStop, v, vStop, xOrder, yOrder); 374 } 375 } 376 } 377 } 378 379 kernels->subIndex = (num - (spatialOrder + 1) * (spatialOrder + 2) / 2) / 2; 380 assert(kernels->u->data.S32[kernels->subIndex] == 0 && 381 kernels->v->data.S32[kernels->subIndex] == 0 && 382 kernels->uStop->data.S32[kernels->subIndex] == 0 && 383 kernels->vStop->data.S32[kernels->subIndex] == 0 && 384 kernels->xOrder->data.S32[kernels->subIndex] == 0 && 385 kernels->yOrder->data.S32[kernels->subIndex] == 0); 386 387 psFree(start); 388 psFree(stop); 389 390 return kernels; 391 } 392 393 189 394 pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, int size, int spatialOrder, 190 const psVector *sigmas, const psVector *orders) 395 const psVector *sigmas, const psVector *orders, int inner, 396 int binning) 191 397 { 192 398 switch (type) { … … 195 401 case PM_SUBTRACTION_KERNEL_ISIS: 196 402 return pmSubtractionKernelsISIS(size, spatialOrder, sigmas, orders); 403 case PM_SUBTRACTION_KERNEL_SPAM: 404 return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning); 405 case PM_SUBTRACTION_KERNEL_FRIES: 406 return pmSubtractionKernelsFRIES(size, spatialOrder, inner); 197 407 default: 198 408 psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unknown kernel type: %x", type); -
trunk/psModules/src/imcombine/pmSubtractionKernels.h
r13340 r14106 6 6 /// Type of subtraction kernel 7 7 typedef enum { 8 PM_SUBTRACTION_KERNEL_POIS, ///< POIS kernel --- delta functions 9 PM_SUBTRACTION_KERNEL_ISIS ///< ISIS kernel --- gaussians modified by polynomials 8 PM_SUBTRACTION_KERNEL_POIS, ///< Pan-STARRS Optimal Image Subtraction --- delta functions 9 PM_SUBTRACTION_KERNEL_ISIS, ///< Traditional kernel --- gaussians modified by polynomials 10 PM_SUBTRACTION_KERNEL_SPAM, ///< Summed Pixels for Advanced Matching --- summed delta functions 11 PM_SUBTRACTION_KERNEL_FRIES, ///< Fibonacci Radius Increases Excellence of Subtraction 10 12 } pmSubtractionKernelsType; 11 13 12 14 /// Kernels specification 13 15 typedef struct { 14 pmSubtractionKernelsType type; 16 pmSubtractionKernelsType type; ///< Type of kernels --- allowing the use of multiple kernels 15 17 psVector *u, *v; ///< Offset (for POIS) or polynomial order (for ISIS) 16 psVector *sigma; ///< Width of Gaussian (for ISIS) 18 psVector *sigma; ///< Gaussian widths (ISIS only) 19 psVector *uStop, *vStop; ///< Width of kernel element (SPAM,FRIES only) 17 20 psVector *xOrder, *yOrder; ///< Spatial Polynomial order (for all) 18 21 int subIndex; ///< Index of kernel to be subtracted (to maintain flux conservation) … … 44 47 ); 45 48 49 /// Generate SPAM kernels 50 pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, ///< Half-size of the kernel 51 int spatialOrder, ///< Order of spatial variations 52 int inner, ///< Inner radius to preserve unbinned 53 int binning ///< Kernel binning factor 54 ); 55 56 /// Generate FRIES kernels 57 pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, ///< Half-size of the kernel 58 int spatialOrder, ///< Order of spatial variations 59 int inner ///< Inner radius to preserve unbinned 60 ); 61 46 62 /// Generate a kernel of a specified type 47 63 pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, ///< Kernel type … … 49 65 int spatialOrder, ///< Order of spatial variations 50 66 const psVector *sigmas, ///< Gaussian widths 51 const psVector *orders ///< Polynomial order of gaussians 52 ); 67 const psVector *orders, ///< Polynomial order of gaussians 68 int inner, ///< Inner radius to preserve unbinned 69 int binning ///< Kernel binning factor 70 ); 53 71 54 72 #endif -
trunk/psModules/src/imcombine/pmSubtractionStamps.c
r13735 r14106 6 6 #include <pslib.h> 7 7 8 #include "pmSubtraction.h"9 8 #include "pmSubtractionStamps.h" 10 9 … … 40 39 } 41 40 42 psArray *pmSubtractionFindStamps(psArray *stamps, const psImage *image, const psImage *subMask, 43 float threshold, float spacing) 41 42 psArray *pmSubtractionFindStamps(psArray *stamps, const psImage *image, const psImage *mask, 43 psMaskType maskVal, psMaskType used, float threshold, 44 float spacing, int border) 44 45 { 45 46 PS_ASSERT_IMAGE_NON_NULL(image, NULL); 46 47 PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, 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); 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; 51 56 } 52 57 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); 53 61 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; 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 61 65 62 66 if (stamps) { … … 74 78 } 75 79 } 80 // Footprint of image 81 int numRows = image->numRows; 82 int numCols = image->numCols; 76 83 77 84 for (int j = 0, index = 0; j < yNumStamps; j++) { … … 89 96 90 97 // Bounds of region to search for stamp 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);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); 97 104 98 105 for (int y = yMin; y <= yMax ; y++) { 99 106 for (int x = xMin; x <= xMax ; x++) { 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; 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 } 107 133 } 108 134 } -
trunk/psModules/src/imcombine/pmSubtractionStamps.h
r13735 r14106 3 3 4 4 #include <pslib.h> 5 5 6 6 7 /// Status of stamp … … 25 26 const psImage *image, ///< Image for which to find stamps 26 27 const psImage *mask, ///< Mask 28 psMaskType maskVal, ///< Value for mask 29 psMaskType bad, ///< Mask value for bad stamps 27 30 float threshold, ///< Threshold for stamps in the image 28 float spacing ///< Rough spacing for stamps 31 float spacing, ///< Rough spacing for stamps 32 int border ///< Size of border around image; no stamps in border 29 33 ); 30 34
Note:
See TracChangeset
for help on using the changeset viewer.
