Changeset 24244 for branches/cnb_branches/cnb_branch_20090301/psModules/src/imcombine/pmSubtractionStamps.c
- Timestamp:
- May 26, 2009, 1:59:32 PM (17 years ago)
- Location:
- branches/cnb_branches/cnb_branch_20090301
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/cnb_branches/cnb_branch_20090301
- Property svn:mergeinfo changed
-
branches/cnb_branches/cnb_branch_20090301/psModules
- Property svn:mergeinfo changed
-
branches/cnb_branches/cnb_branch_20090301/psModules/src/imcombine/pmSubtractionStamps.c
r21363 r24244 20 20 #include "pmModel.h" 21 21 #include "pmSource.h" 22 22 #include "pmErrorCodes.h" 23 23 24 24 #include "pmSubtraction.h" … … 78 78 } 79 79 80 // Is this position unmasked? 81 static bool checkStampMask(int x, int y, // Coordinates of stamp 82 const psImage *mask, // Mask 83 pmSubtractionMode mode, // Mode for subtraction 84 int footprint // Footprint to check for Bad Stuff 85 ) 86 { 87 if (!mask) { 88 return true; 89 } 90 91 bool clean = true; // Is the footprint clean? 92 int numCols = mask->numCols, numRows = mask->numRows; // Size of image 93 94 // Check the footprint bounds 95 if (x < footprint || x >= numCols - footprint || y < footprint || y >= numRows - footprint) { 96 clean = false; 97 } 98 99 // Determine mask value 100 psImageMaskType maskVal = PM_SUBTRACTION_MASK_BORDER | PM_SUBTRACTION_MASK_BAD_1 | PM_SUBTRACTION_MASK_BAD_2; 101 switch (mode) { 102 case PM_SUBTRACTION_MODE_1: 103 maskVal |= PM_SUBTRACTION_MASK_CONVOLVE_1; 104 break; 105 case PM_SUBTRACTION_MODE_2: 106 maskVal |= PM_SUBTRACTION_MASK_CONVOLVE_2; 107 break; 108 case PM_SUBTRACTION_MODE_UNSURE: 109 case PM_SUBTRACTION_MODE_DUAL: 110 maskVal |= PM_SUBTRACTION_MASK_CONVOLVE_1 | PM_SUBTRACTION_MASK_CONVOLVE_2; 111 break; 112 default: 113 psAbort("Unsupported subtraction mode: %x", mode); 114 } 115 116 // Check the immediate pixel 117 if (clean && (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & (maskVal | PM_SUBTRACTION_MASK_REJ))) { 118 clean = false; 119 } 120 121 int xMin = PS_MAX(x - footprint, 0), xMax = PS_MIN(x + footprint, numCols - 1); // Bounds in x 122 int yMin = PS_MAX(y - footprint, 0), yMax = PS_MIN(y + footprint, numRows - 1); // Bounds in y 123 124 // Check the footprint 125 if (clean) { 126 for (int j = yMin; j <= yMax; j++) { 127 for (int i = xMin; i <= xMax; i++) { 128 if (mask->data.PS_TYPE_IMAGE_MASK_DATA[j][i] & maskVal) { 129 clean = false; 130 goto CHECK_STAMP_MASK_DONE; 131 } 80 81 // Search a region for a suitable stamp 82 bool stampSearch(int *xStamp, int *yStamp, // Coordinates of stamp, to return 83 float *fluxStamp, // Flux of stamp, to return 84 const psImage *image1, const psImage *image2, // Images to search 85 float thresh1, float thresh2, // Thresholds for images 86 const psImage *subMask, // Subtraction mask 87 int xMin, int xMax, int yMin, int yMax, // Bounds of search 88 int numCols, int numRows, // Size of images 89 int border // Border around image 90 ) 91 { 92 bool found = false; // Found a suitable stamp? 93 *fluxStamp = -INFINITY; // Flux of best stamp 94 95 // Ensure we're not going to go outside the bounds of the image 96 xMin = PS_MAX(border, xMin); 97 xMax = PS_MIN(numCols - border - 1, xMax); 98 yMin = PS_MAX(border, yMin); 99 yMax = PS_MIN(numRows - border - 1, yMax); 100 101 for (int y = yMin; y <= yMax; y++) { 102 for (int x = xMin; x <= xMax; x++) { 103 if ((image1 && image1->data.F32[y][x] < thresh1) || 104 (image2 && image2->data.F32[y][x] < thresh2)) { 105 continue; 132 106 } 133 } 134 } 135 CHECK_STAMP_MASK_DONE: 136 137 if (!clean) { 138 // Mask the footprint, so we don't select something near it again 139 for (int j = yMin; j <= yMax; j++) { 140 for (int i = xMin; i <= xMax; i++) { 141 mask->data.PS_TYPE_IMAGE_MASK_DATA[j][i] |= PM_SUBTRACTION_MASK_REJ; 107 108 if (subMask && subMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 109 (PM_SUBTRACTION_MASK_BORDER | PM_SUBTRACTION_MASK_REJ)) { 110 return false; 142 111 } 143 } 144 } 145 146 return clean; 147 } 112 113 // We take the MIN to attempt to avoid transients in both images 114 float flux = (image1 && image2) ? PS_MIN(image1->data.F32[y][x], image2->data.F32[y][x]) : 115 ((image1) ? image1->data.F32[y][x] : image2->data.F32[y][x]); // Flux at pixel 116 if (flux > *fluxStamp) { 117 *xStamp = x; 118 *yStamp = y; 119 *fluxStamp = flux; 120 found = true; 121 } 122 } 123 } 124 125 return found; 126 } 127 128 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 129 // Functions for generating ds9 region files 130 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 131 132 static bool ds9regions = false; // Save ds9 region files? 133 134 void pmSubtractionRegions(bool state) 135 { 136 ds9regions = state; 137 } 138 139 FILE *pmSubtractionStampsFile(const pmSubtractionStampList *stamps, const char *filename, 140 const char *description) 141 { 142 if (!ds9regions || !stamps || !filename) { 143 return NULL; 144 } 145 146 psLogMsg("psModules.imcombine", PS_LOG_INFO, "Writing %s to ds9 region file: %s", 147 description, filename); 148 149 FILE *file = fopen(filename, "w"); 150 151 // Outline the stamps 152 for (int i = 0; i < stamps->num; i++) { 153 psRegion *region = stamps->regions->data[i]; // Region of interest 154 float xCentre = 0.5 * (region->x0 + region->x1), yCentre = 0.5 * (region->y0 + region->y1); 155 fprintf(file, "image;box(%f,%f,%f,%f,0) # color=blue\nimage;text(%f,%f,{%d}) # color=blue\n", 156 xCentre, yCentre, region->x1 - region->x0, region->y1 - region->y0, 157 xCentre, yCentre, i); 158 } 159 160 return file; 161 } 162 163 void pmSubtractionStampPrint(FILE *ds9, float x, float y, float size, const char *color) 164 { 165 if (!ds9regions || !ds9) { 166 return; 167 } 168 fprintf(ds9, "image;circle(%f,%f,%f)", x, y, size); 169 if (color && strlen(color) > 0) { 170 fprintf(ds9, " # color=%s", color); 171 } 172 fprintf(ds9, "\n"); 173 return; 174 } 175 148 176 149 177 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// … … 225 253 226 254 227 pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, const psImage *image, 228 const psImage *subMask, const psRegion *region, 229 float threshold, int footprint, float spacing, 255 pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, const psImage *image1, 256 const psImage *image2, const psImage *subMask, 257 const psRegion *region, float thresh1, float thresh2, 258 int size, int footprint, float spacing, 230 259 pmSubtractionMode mode) 231 260 { 232 PS_ASSERT_IMAGE_NON_NULL(image, NULL); 233 PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL); 261 if (!image1 && !image2) { 262 psError(PS_ERR_UNEXPECTED_NULL, true, "Must specify an image"); 263 return NULL; 264 } 265 int numCols = 0, numRows = 0; // Size of images 266 if (image1) { 267 PS_ASSERT_IMAGE_NON_NULL(image1, NULL); 268 PS_ASSERT_IMAGE_TYPE(image1, PS_TYPE_F32, NULL); 269 if (subMask) { 270 PS_ASSERT_IMAGES_SIZE_EQUAL(image1, subMask, NULL); 271 } 272 numCols = image1->numCols; 273 numRows = image1->numRows; 274 } 275 if (image2) { 276 PS_ASSERT_IMAGE_NON_NULL(image2, NULL); 277 PS_ASSERT_IMAGE_TYPE(image2, PS_TYPE_F32, NULL); 278 if (subMask) { 279 PS_ASSERT_IMAGES_SIZE_EQUAL(image2, subMask, NULL); 280 } 281 numCols = image2->numCols; 282 numRows = image2->numRows; 283 } 284 if (image1 && image2) { 285 PS_ASSERT_IMAGES_SIZE_EQUAL(image1, image2, NULL); 286 } 234 287 if (subMask) { 235 288 PS_ASSERT_IMAGE_NON_NULL(subMask, NULL); 236 PS_ASSERT_IMAGES_SIZE_EQUAL(image, subMask, NULL);237 289 PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_IMAGE_MASK, NULL); 238 } 290 PS_ASSERT_IMAGE_SIZE(subMask, numCols, numRows, NULL); 291 } 292 PS_ASSERT_INT_POSITIVE(size, NULL); 239 293 PS_ASSERT_INT_NONNEGATIVE(footprint, NULL); 240 294 PS_ASSERT_FLOAT_LARGER_THAN(spacing, 0.0, NULL); … … 246 300 return false; 247 301 } 248 if (region->x0 < 0 || region->x1 > image->numCols ||249 region->y0 < 0 || region->y1 > image->numRows) {302 if (region->x0 < 0 || region->x1 > numCols || 303 region->y0 < 0 || region->y1 > numRows) { 250 304 psString string = psRegionToString(*region); 251 305 psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Input region (%s) does not fit in image (%dx%d)", 252 string, image->numCols, image->numRows);306 string, numCols, numRows); 253 307 psFree(string); 254 308 return false; … … 256 310 } 257 311 258 int numRows = image->numRows, numCols = image->numCols; // Size of image312 int border = size + footprint; // Border size 259 313 260 314 if (!stamps) { … … 281 335 numSearch++; 282 336 283 float xStamp = 0, yStamp = 0;// Coordinates of stamp284 float fluxStamp = NAN;// Flux of stamp285 bool goodStamp = false; 337 int xStamp = 0, yStamp = 0; // Coordinates of stamp 338 float fluxStamp = -INFINITY; // Flux of stamp 339 bool goodStamp = false; // Found a good stamp? 286 340 287 341 // A couple different ways of finding stamps: … … 291 345 psVector *fluxList = stamps->flux->data[i]; // List of stamp fluxes 292 346 293 // Take stamp off the top of the (sorted) list 294 if (xList->n > 0) { 295 int index = xList->n - 1; // Index of new stamp 296 xStamp = xList->data.F32[index]; 297 yStamp = yList->data.F32[index]; 298 fluxStamp = fluxList->data.F32[index]; 347 // Take stamps off the top of the (sorted) list 348 for (int j = xList->n - 1; j >= 0 && !goodStamp; j--) { 349 int xCentre = xList->data.F32[j] + 0.5, yCentre = yList->data.F32[j] + 0.5;// Stamp centre 299 350 300 351 // Chop off the top of the list 301 xList->n = index; 302 yList->n = index; 303 fluxList->n = index; 304 305 goodStamp = (fluxStamp > threshold) ? true : false; 306 } else { 307 psTrace("psModules.imcombine", 9, "No sources in subregion %d", i); 352 xList->n = j; 353 yList->n = j; 354 fluxList->n = j; 355 356 // Fish around a bit to see if we can find a pixel that isn't masked 357 psTrace("psModules.imcombine", 7, "Searching for stamp %d around %d,%d\n", 358 i, xCentre, yCentre); 359 360 // Search bounds 361 int search = footprint - size; // Search radius 362 int xMin = PS_MAX(border, xCentre - search); 363 int xMax = PS_MIN(numCols - border -1, xCentre + search); 364 int yMin = PS_MAX(border, yCentre - search); 365 int yMax = PS_MIN(numRows - border - 1, yCentre + search); 366 367 goodStamp = stampSearch(&xStamp, &yStamp, &fluxStamp, image1, image2, thresh1, thresh2, 368 subMask, xMin, xMax, yMin, yMax, numCols, numRows, border); 308 369 } 309 370 } else { 310 371 // Use a simple method of automatically finding stamps --- take the highest pixel in the 311 372 // subregion 312 fluxStamp = threshold;313 373 psRegion *subRegion = stamps->regions->data[i]; // Sub-region of interest 314 for (int y = subRegion->y0; y <= subRegion->y1; y++) { 315 for (int x = subRegion->x0; x <= subRegion->x1; x++) { 316 if (checkStampMask(x, y, subMask, mode, footprint) && 317 image->data.F32[y][x] > fluxStamp) { 318 fluxStamp = image->data.F32[y][x]; 319 xStamp = x; 320 yStamp = y; 321 goodStamp = true; 322 } 323 } 324 } 374 375 goodStamp = stampSearch(&xStamp, &yStamp, &fluxStamp, image1, image2, thresh1, thresh2, 376 subMask, subRegion->x0, subRegion->x1, subRegion->y0, subRegion->y1, 377 numCols, numRows, border); 325 378 } 326 379 … … 355 408 if (numGood == 0 && numFound == 0) { 356 409 // No good stamps 410 psError(PM_ERR_STAMPS, true, "Unable to find suitable stamps"); 357 411 psFree(stamps); 358 412 return NULL; … … 363 417 364 418 419 365 420 pmSubtractionStampList *pmSubtractionStampsSet(const psVector *x, const psVector *y, 366 421 const psImage *image, const psImage *subMask, 367 const psRegion *region, int footprint, float spacing,368 pmSubtractionMode mode)422 const psRegion *region, int size, int footprint, 423 float spacing, pmSubtractionMode mode) 369 424 370 425 { … … 378 433 PS_ASSERT_IMAGE_NON_NULL(subMask, NULL); 379 434 PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_IMAGE_MASK, NULL); 380 if (image) { 381 PS_ASSERT_IMAGE_NON_NULL(image, NULL); 382 PS_ASSERT_IMAGES_SIZE_EQUAL(image, subMask, NULL); 383 } 384 } 435 PS_ASSERT_IMAGES_SIZE_EQUAL(image, subMask, NULL); 436 } 437 PS_ASSERT_INT_POSITIVE(size, NULL); 438 PS_ASSERT_INT_POSITIVE(footprint, NULL); 385 439 PS_ASSERT_FLOAT_LARGER_THAN(spacing, 0.0, NULL); 386 440 … … 389 443 region, footprint, spacing); // Stamp list 390 444 int numStamps = stamps->num; // Number of stamps 445 446 psString ds9name = NULL; // Filename for ds9 region file 447 static int ds9num = 0; // File number for ds9 region file 448 psStringAppend(&ds9name, "stamps_all_%d.ds9", ds9num); 449 FILE *ds9 = pmSubtractionStampsFile(stamps, ds9name, "all stamps"); // ds9 region file 450 ds9num++; 391 451 392 452 // Initialise the lists … … 408 468 psTrace("psModules.imcombine", 9, "Rejecting input stamp (%d,%d) because outside region", 409 469 xPix, yPix); 410 continue; 411 } 412 if (!checkStampMask(xPix, yPix, subMask, mode, footprint)) { 413 // Not a good stamp 414 psTrace("psModules.imcombine", 9, "Rejecting input stamp (%d,%d) because bad mask", 415 xPix, yPix); 470 pmSubtractionStampPrint(ds9, xPix, yPix, footprint, "red"); 416 471 continue; 417 472 } … … 437 492 psTrace("psModules.imcombine", 9, "Putting input stamp (%d,%d) into subregion %d", 438 493 xPix, yPix, j); 494 pmSubtractionStampPrint(ds9, xPix, yPix, footprint, "green"); 439 495 } 440 496 } … … 443 499 psTrace("psModules.imcombine", 9, "Unable to find subregion for stamp (%d,%d)", 444 500 xPix, yPix); 445 } 501 pmSubtractionStampPrint(ds9, xPix, yPix, footprint, "yellow"); 502 } 503 } 504 505 if (ds9) { 506 fclose(ds9); 446 507 } 447 508 … … 474 535 } 475 536 476 477 537 return stamps; 478 538 } … … 544 604 } 545 605 546 #if 0547 bool pmSubtractionStampsGenerate(pmSubtractionStampList *stamps, float fwhm, int kernelSize)548 {549 PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);550 PS_ASSERT_FLOAT_LARGER_THAN(fwhm, 0.0, false);551 PS_ASSERT_INT_NONNEGATIVE(kernelSize, false);552 553 int size = kernelSize + stamps->footprint; // Size of postage stamps554 int num = stamps->num; // Number of stamps555 float sigma = fwhm / (2.0 * sqrtf(2.0 * log(2.0))); // Gaussian sigma556 557 for (int i = 0; i < num; i++) {558 pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest559 if (!(stamp->status & PM_SUBTRACTION_STAMP_CALCULATE)) {560 continue;561 }562 563 float x = stamp->x, y = stamp->y; // Coordinates of stamp564 float flux = stamp->flux; // Flux of star565 if (!isfinite(flux)) {566 psWarning("Unable to generate PSF for stamp %d --- bad flux.", i);567 stamp->status = PM_SUBTRACTION_STAMP_REJECTED;568 continue;569 }570 571 float xStamp = x - (int)(x + 0.5); // x coordinate of star in stamp frame572 float yStamp = y - (int)(y + 0.5); // y coordinate of star in stamp frame573 574 psFree(stamp->image2);575 stamp->image2 = psKernelAlloc(-size, size, -size, size);576 psKernel *target = stamp->image2; // Target stamp577 578 // Put in a Waussian, just for fun!579 for (int v = -size; v <= size; v++) {580 for (int u = -size; u <= size; u++) {581 float z = (PS_SQR(u + xStamp) + PS_SQR(v + yStamp)) / (2.0 * PS_SQR(sigma));582 target->kernel[v][u] = flux / sigma * 0.5 * M_2_SQRTPI * M_SQRT1_2 / (1.0 + z + PS_SQR(z));583 }584 }585 586 }587 588 return true;589 }590 #endif591 592 606 pmSubtractionStampList *pmSubtractionStampsSetFromSources(const psArray *sources, const psImage *image, 593 607 const psImage *subMask, const psRegion *region, 594 int footprint, float spacing,608 int size, int footprint, float spacing, 595 609 pmSubtractionMode mode) 596 610 { … … 621 635 y->n = numOut; 622 636 623 pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region, 637 pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region, size, 624 638 footprint, spacing, mode); // Stamps to return 625 639 psFree(x); … … 636 650 pmSubtractionStampList *pmSubtractionStampsSetFromFile(const char *filename, const psImage *image, 637 651 const psImage *subMask, const psRegion *region, 638 int footprint, float spacing, pmSubtractionMode mode) 652 int size, int footprint, float spacing, 653 pmSubtractionMode mode) 639 654 { 640 655 PS_ASSERT_STRING_NON_EMPTY(filename, NULL); … … 652 667 psBinaryOp(y, y, "-", psScalarAlloc(1.0, PS_TYPE_F32)); 653 668 654 pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region, footprint,669 pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region, size, footprint, 655 670 spacing, mode); 656 671 psFree(data);
Note:
See TracChangeset
for help on using the changeset viewer.
