Changeset 31633
- Timestamp:
- Jun 15, 2011, 6:09:25 PM (15 years ago)
- Location:
- trunk
- Files:
-
- 8 edited
-
ippconfig/recipes/fitstypes.mdc (modified) (1 diff)
-
ippconfig/recipes/ppStack.config (modified) (1 diff)
-
ippconfig/recipes/reductionClasses.mdc (modified) (1 diff, 1 prop)
-
psLib/src/fits/psFits.h (modified) (2 diffs)
-
psLib/src/fits/psFitsImage.c (modified) (13 diffs)
-
psLib/src/fits/psFitsScale.c (modified) (14 diffs)
-
psLib/src/fits/psFitsScale.h (modified) (2 diffs)
-
psModules/src/camera/pmFPAfileDefine.c (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/ippconfig/recipes/fitstypes.mdc
r30636 r31633 101 101 COMP_STACK METADATA 102 102 BITPIX S32 16 103 SCALING STR LOG_STDEV_POSITIVE 104 STDEV.BITS S32 4 105 STDEV.NUM F32 10 103 SCALING STR ASINH_RANGE 104 # SCALING STR LOG_RANGE 105 # STDEV.BITS S32 4 106 # STDEV.NUM F32 100 106 107 COMPRESSION STR RICE 107 108 TILE.X S32 0 -
trunk/ippconfig/recipes/ppStack.config
r31422 r31633 132 132 END 133 133 134 # 134 # Stack mode for testing stacking procedure 135 TESTSTACK METADATA 136 VARIANCE BOOL FALSE 137 END 138 135 139 136 140 STACK_DEEP METADATA -
trunk/ippconfig/recipes/reductionClasses.mdc
- Property svn:mergeinfo changed
/branches/czw_branch/20110406/ippconfig/recipes/reductionClasses.mdc merged: 31606,31618,31632
r31456 r31633 284 284 END 285 285 286 TESTSTACK METADATA 287 STACK_PPSTACK STR TESTSTACK 288 STACK_PPSUB STR STACK 289 STACK_PSPHOT STR STACK 290 END 291 286 292 # basic science analysis 287 293 PHOTFEST.SINGLE1 METADATA - Property svn:mergeinfo changed
-
trunk/psLib/src/fits/psFits.h
r30636 r31633 55 55 PS_FITS_SCALE_LOG_STDEV_NEGATIVE, ///< Take logarithm, Auto-scale to sample stdev, place mean at upper limit 56 56 PS_FITS_SCALE_LOG_STDEV_BOTH, ///< Take logarithm, Auto-scale to sample stdev, place mean at middle 57 PS_FITS_SCALE_LOG_MANUAL ///< Manual scaling (use specified BSCALE, BZERO, and BOFFSET) 57 PS_FITS_SCALE_LOG_MANUAL, ///< Manual scaling (use specified BSCALE, BZERO, and BOFFSET) 58 PS_FITS_SCALE_ASINH_RANGE, ///< Do asinh scaling, auto-scale to preserve dynamic range 59 PS_FITS_SCALE_ASINH_STDEV_POSITIVE, ///< Do asinh scaling, auto-scale to sample stdev, place mean at lower limit 60 PS_FITS_SCALE_ASINH_STDEV_NEGATIVE, ///< Do asinh scaling, auto-scale to sample stdev, place mean at upper limit 61 PS_FITS_SCALE_ASINH_STDEV_BOTH, ///< Do asinh scaling, auto-scale to sample stdev, place mean at middle 62 PS_FITS_SCALE_ASINH_MANUAL ///< Manual scaling after doing asinh scaling. 58 63 } psFitsScaling; 59 64 … … 72 77 double bscale, bzero; ///< Manually specified BSCALE and BZERO (for SCALE_MANUAL) 73 78 double boffset; ///< Manually specified BOFFSET (for SCALE_MANUAL) 79 double bsoften; ///< Manuallly specified BSOFTEN (for SCALE_MANUAL) 74 80 double mean, stdev; ///< Mean and standard deviation of image 75 81 int stdevBits; ///< Number of bits to sample a standard deviation (for SCALE_STDEV_*) -
trunk/psLib/src/fits/psFitsImage.c
r31597 r31633 55 55 int psDatatype; // psLib data type 56 56 bool is_logscaled; // is this image log scaled using BOFFSET? 57 bool is_asinh; // is this image asinh scaled using BSOFTEN? 57 58 } p_psFitsReadInfo; 58 59 … … 129 130 130 131 // Check scale and zero 131 double bscale = 0.0, bzero = 0.0, boffset = NAN ; // Scale and zero point132 double bscale = 0.0, bzero = 0.0, boffset = NAN, bsoften = NAN; // Scale and zero point 132 133 if (fits_read_key_dbl(fits->fd, "BSCALE", &bscale, NULL, &status) && status != KEY_NO_EXIST) { 133 134 psFitsError(status, true, "Unable to read header."); … … 153 154 info->is_logscaled = false; 154 155 } 156 status = 0; 157 if (fits_read_key_dbl(fits->fd, "BSOFTEN", &bsoften, NULL, &status) && status != KEY_NO_EXIST) { 158 psFitsError(status, true, "Unable to read header."); 159 goto bad; 160 } 161 if (status == KEY_NO_EXIST) { 162 info->is_asinh = false; 163 } 164 else if (isfinite(bsoften)) { 165 info->is_asinh = true; 166 info->is_logscaled = false; 167 } 168 else { 169 info->is_asinh = false; 170 } 171 155 172 status = 0; 156 173 … … 262 279 double *bzero, // Zero point applied 263 280 double *boffset, // Log offset applied 281 double *bsoften, // asinh offset applied 264 282 long *blank, // Blank value (integer data) 265 283 psFitsFloat *floatType, // Type of custom floating-point … … 275 293 psAssert(bzero, "impossible"); 276 294 psAssert(boffset, "impossible"); 295 psAssert(bsoften, "impossible"); 277 296 psAssert(floatType, "impossible"); 278 297 psAssert(fits, "impossible"); … … 322 341 if (newScaleZero) { 323 342 // Choose an appropriate BSCALE and BZERO 324 if (!psFitsScaleDetermine(bscale, bzero, boffset, b lank, image, mask, maskVal, fits)) {343 if (!psFitsScaleDetermine(bscale, bzero, boffset, bsoften, blank, image, mask, maskVal, fits)) { 325 344 // We can't have the write dying for this reason --- try to save it somehow! 326 345 psWarning("Unable to determine BSCALE and BZERO for image --- refusing to quantise."); … … 348 367 } 349 368 350 return psFitsScaleForDisk(image, fits, *bscale, *bzero, *boffset, rng);369 return psFitsScaleForDisk(image, fits, *bscale, *bzero, *boffset, *bsoften, rng); 351 370 } 352 371 … … 480 499 outImage = psFitsFloatImageFromDisk(outImage, inImage, floatType); 481 500 } 482 483 501 // Need to apply BOFFSET if info->is_logscaled is true 502 484 503 if (info->is_logscaled) { 485 504 double boffset; … … 487 506 status = 0; 488 507 fits_read_key_dbl(fits->fd, "BOFFSET", &boffset, NULL, &status); 489 psImage *newImage = psFitsScaleFromDisk(outImage,boffset); 508 psImage *newImage = psFitsScaleFromDisk(outImage,boffset,0.0); 509 psFree (outImage); 510 outImage = newImage; 511 } 512 else if (info->is_asinh) { 513 double bsoften,boffset; 514 int status; 515 status = 0; 516 fits_read_key_dbl(fits->fd, "BSOFTEN", &bsoften, NULL, &status); 517 fits_read_key_dbl(fits->fd, "BOFFSET", &boffset, NULL, &status); 518 psImage *newImage = psFitsScaleFromDisk(outImage,boffset,bsoften); 490 519 psFree (outImage); 491 520 outImage = newImage; … … 602 631 double bscale = NAN, bzero = NAN; // Scale and zero point to put in header 603 632 double boffset = NAN; // Log offset to put into header. 633 double bsoften = NAN; // Asinh softening parameter to put into header 604 634 long blank = 0; // Blank (undefined) value for image 605 635 psFitsFloat floatType; // Custom floating-point convention type 606 psImage *diskImage = imageToDiskRepresentation(&bscale, &bzero, &boffset, &b lank, &floatType, fits, image,636 psImage *diskImage = imageToDiskRepresentation(&bscale, &bzero, &boffset, &bsoften, &blank, &floatType, fits, image, 607 637 mask, maskVal, NULL, true); // Image to write out 608 638 if (!diskImage) { … … 696 726 if (options&&(!((options->scaling == PS_FITS_SCALE_LOG_RANGE)|| 697 727 (options->scaling == PS_FITS_SCALE_LOG_MANUAL)|| 698 (options->scaling == PS_FITS_SCALE_LOG_STDEV_POSITIVE)|| 699 (options->scaling == PS_FITS_SCALE_LOG_STDEV_NEGATIVE)|| 700 (options->scaling == PS_FITS_SCALE_LOG_STDEV_BOTH)))) { 728 (options->scaling == PS_FITS_SCALE_LOG_STDEV_POSITIVE)|| 729 (options->scaling == PS_FITS_SCALE_LOG_STDEV_NEGATIVE)|| 730 (options->scaling == PS_FITS_SCALE_LOG_STDEV_BOTH)|| 731 (options->scaling == PS_FITS_SCALE_ASINH_RANGE)|| 732 (options->scaling == PS_FITS_SCALE_ASINH_MANUAL)|| 733 (options->scaling == PS_FITS_SCALE_ASINH_STDEV_POSITIVE)|| 734 (options->scaling == PS_FITS_SCALE_ASINH_STDEV_NEGATIVE)|| 735 (options->scaling == PS_FITS_SCALE_ASINH_STDEV_BOTH)))) { 701 736 if (psMetadataLookup(header,"BOFFSET")) { 702 737 psMetadataRemoveKey(header,"BOFFSET"); 738 } 739 } 740 // Remove any BSOFTENvalues that exist in the header if we are not using that scaling anymore 741 if (options&&(!((options->scaling == PS_FITS_SCALE_ASINH_RANGE)|| 742 (options->scaling == PS_FITS_SCALE_ASINH_MANUAL)|| 743 (options->scaling == PS_FITS_SCALE_ASINH_STDEV_POSITIVE)|| 744 (options->scaling == PS_FITS_SCALE_ASINH_STDEV_NEGATIVE)|| 745 (options->scaling == PS_FITS_SCALE_ASINH_STDEV_BOTH)))) { 746 if (psMetadataLookup(header,"BSOFTEN")) { 747 psMetadataRemoveKey(header,"BSOFTEN"); 703 748 } 704 749 } … … 727 772 (options->scaling == PS_FITS_SCALE_LOG_STDEV_BOTH)))) { 728 773 fits_write_key_dbl(fits->fd, "BOFFSET", boffset, 12, 729 "Scaling: TRUE = BZERO + BSCALE * 10**(DISK) + BOFFSET)", &status);774 "Scaling: ORIGINAL = 10**(TRUE) + BOFFSET)", &status); 730 775 } 776 if (options&&(((options->scaling == PS_FITS_SCALE_ASINH_RANGE)|| 777 (options->scaling == PS_FITS_SCALE_ASINH_MANUAL)|| 778 (options->scaling == PS_FITS_SCALE_ASINH_STDEV_POSITIVE)|| 779 (options->scaling == PS_FITS_SCALE_ASINH_STDEV_NEGATIVE)|| 780 (options->scaling == PS_FITS_SCALE_ASINH_STDEV_BOTH)))) { 781 fits_write_key_dbl(fits->fd, "BOFFSET", boffset, 12, 782 "Scaling: NORM = DISK + BOFFSET", &status); 783 fits_write_key_dbl(fits->fd, "BSOFTEN", bsoften, 12, 784 "Scaling: TRUE = 2 * BSOFTEN * sinh(NORM/a)",&status); 785 786 787 } 788 731 789 if (psFitsError(status, true, "Could not write BSCALE/BZERO headers to file.")) { 732 790 success = false; … … 831 889 double bscale = NAN, bzero = NAN; // Scale and zero point to put in header 832 890 double boffset = NAN; // Log offset to put in header 891 double bsoften = NAN; // Asinh softening parameter to put in header 833 892 long blank = 0; // Blank (undefined) value for image 834 893 psFitsFloat floatType; // Custom floating-point convention type 835 psImage *diskImage = imageToDiskRepresentation(&bscale, &bzero, &boffset, &b lank, &floatType, fits, input,894 psImage *diskImage = imageToDiskRepresentation(&bscale, &bzero, &boffset, &bsoften, &blank, &floatType, fits, input, 836 895 mask, maskVal, NULL, false); // Image to write out 837 896 if (!diskImage) { -
trunk/psLib/src/fits/psFitsScale.c
r30636 r31633 214 214 case PS_FITS_SCALE_STDEV_POSITIVE: 215 215 case PS_FITS_SCALE_LOG_STDEV_POSITIVE: 216 case PS_FITS_SCALE_ASINH_STDEV_POSITIVE: 216 217 // Put (mean - N sigma) at the lowest possible value: predominantly positive images 217 218 imageVal = mean - options->stdevNum * stdev; … … 220 221 case PS_FITS_SCALE_STDEV_NEGATIVE: 221 222 case PS_FITS_SCALE_LOG_STDEV_NEGATIVE: 223 case PS_FITS_SCALE_ASINH_STDEV_NEGATIVE: 222 224 // Put (mean + N sigma) at the highest possible value: predominantly negative images 223 225 imageVal = mean + options->stdevNum * stdev; … … 226 228 case PS_FITS_SCALE_STDEV_BOTH: 227 229 case PS_FITS_SCALE_LOG_STDEV_BOTH: 230 case PS_FITS_SCALE_ASINH_STDEV_BOTH: 228 231 // Put mean right in the middle: images with an equal abundance of positive and negative values 229 232 imageVal = mean; … … 241 244 } 242 245 246 247 243 248 244 249 static bool logscaleStdev(double *bscale, // Scaling, to return … … 456 461 } 457 462 463 static bool asinhStdev(double *bscale, // Scaling, to return 464 double *bzero, // Zero point, to return 465 double *boffset, // asinh flux zeropoint, to return 466 double *bsoften, // asinh softening parameter, to return 467 const psImage *image, // Image to scale 468 const psImage *mask, // Mask image 469 psImageMaskType maskVal, // value to mask 470 const psFitsOptions *options // FITS options 471 ) 472 { 473 psAssert(bscale, "impossible"); 474 psAssert(bzero, "impossible"); 475 psAssert(boffset, "impossible"); 476 psAssert(bsoften, "impossible"); 477 psAssert(image, "impossible"); 478 psAssert(options, "impossible"); 479 480 psTrace("psLib.fits", 3, "Scaling image by asinh method"); 481 int numCols = image->numCols, numRows = image->numRows; // Size of image 482 483 // Measure the mean and stdev 484 // psImageBackground automatically excludes pixels that are non-finite, so we don't need to bother about a 485 // mask. 486 psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); 487 psStats *stats = psStatsAlloc(MEAN_STAT | STDEV_STAT); // Statistics object 488 double mean, stdev; // Mean and standard deviation 489 if (!psImageBackground(stats, NULL, image, mask, maskVal, rng)) { 490 // It could be because the image is entirely masked, in which case we don't want to error 491 bool good = false; // Any good pixels? 492 493 494 // Find good pixels in an image, by image type 495 #define GOOD_PIXELS_CASE(TYPE) \ 496 case PS_TYPE_##TYPE: \ 497 for (int y = 0; y < image->numRows && !good; y++) { \ 498 for (int x = 0; x < image->numCols && !good; x++) { \ 499 if (mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal)) { \ 500 continue; \ 501 } \ 502 if (!isfinite(image->data.TYPE[y][x])) { \ 503 continue; \ 504 } \ 505 good = true; \ 506 } \ 507 } \ 508 break; 509 510 switch (image->type.type) { 511 GOOD_PIXELS_CASE(F32); 512 GOOD_PIXELS_CASE(F64); 513 default: 514 psAbort("Unsupported case: %x", image->type.type); 515 } 516 517 if (!good) { 518 psLogMsg("psLib.fits", PS_LOG_DETAIL, "Image has no good pixels, setting BSCALE = 1, BZERO = 0"); 519 psErrorClear(); 520 *bscale = 1.0; 521 *bzero = 0.0; 522 *bsoften = NAN; 523 psFree(rng); 524 psFree(stats); 525 return true; 526 } 527 528 // There are some good pixels in there somewhere; psImageBackground just didn't find them 529 psLogMsg("psLib.fits", PS_LOG_DETAIL, 530 "Couldn't measure background statistics for image quantisation; retrying."); 531 psErrorClear(); 532 // Retry using all the available pixels 533 stats->nSubsample = image->numCols * image->numRows + 1; 534 if (!psImageStats(stats, image, mask, maskVal)) { 535 psLogMsg("psLib.fits", PS_LOG_DETAIL, 536 "Couldn't measure background statistics for image quantisation (attempt 2); retrying."); 537 psErrorClear(); 538 // Retry with desperate statistic 539 stats->options = DESPERATE_MEAN_STAT | DESPERATE_STDEV_STAT; 540 if (!psImageStats(stats, image, mask, maskVal)) { 541 psError(PS_ERR_UNKNOWN, false, "Unable to measure background statistics for image"); 542 psFree(rng); 543 psFree(stats); 544 return false; 545 } else { 546 // Desperate retry 547 mean = psStatsGetValue(stats, DESPERATE_MEAN_STAT); 548 stdev = psStatsGetValue(stats, DESPERATE_STDEV_STAT); 549 } 550 } else { 551 // Retry with all available pixels 552 mean = psStatsGetValue(stats, MEAN_STAT); 553 stdev = psStatsGetValue(stats, STDEV_STAT); 554 } 555 } else { 556 // First attempt 557 mean = psStatsGetValue(stats, MEAN_STAT); 558 stdev = psStatsGetValue(stats, STDEV_STAT); 559 } 560 psFree(rng); 561 psFree(stats); 562 if (!isfinite(mean) || !isfinite(stdev)) { 563 psError(PS_ERR_BAD_PARAMETER_VALUE, true, 564 "Mean (%f) or stdev (%f) is non-finite.", mean, stdev); 565 return false; 566 } 567 568 psImage *copy; 569 570 switch (image->type.type) { 571 case PS_TYPE_F32: 572 copy = psImageCopy(NULL,image,PS_TYPE_F32); 573 break; 574 case PS_TYPE_F64: 575 copy = psImageCopy(NULL,image,PS_TYPE_F64); 576 break; 577 default: 578 psError(PS_ERR_UNKNOWN, true, "Target type is not a float: %d", image->type.type); 579 return NULL; 580 break; 581 } 582 // Do scaling 583 float a = 1.0857362; // 2.5 * log(e); 584 *bsoften = sqrt(a) * stdev; 585 *boffset = mean; 586 // float m0 = 0; // Can I just arbitrarily set this? 587 588 switch (image->type.type) { 589 case PS_TYPE_F32: 590 for (int y = 0; y < numRows; y++) { 591 for (int x = 0; x < numCols; x++) { 592 /* if (x == 266 && y == 4584) { */ 593 /* fprintf(stderr,"psFS32: %d %d %g %g %g %g %g\n",x,y,*boffset,*bsoften,image->data.F32[y][x],log10(image->data.F32[y][x] - *boffset), */ 594 /* a * asinh((image->data.F32[y][x] - *boffset) / (2 * *bsoften))); */ 595 /* } */ 596 if (isfinite(image->data.F32[y][x])) { 597 copy->data.F32[y][x] = a * asinh( (image->data.F32[y][x] - *boffset) / (2 * *bsoften));// - 2.5 * log10(b) + m0; 598 } 599 else { 600 copy->data.F32[y][x] = image->data.F32[y][x]; 601 } 602 } 603 } 604 break; 605 case PS_TYPE_F64: 606 for (int y = 0; y < numRows; y++) { 607 for (int x = 0; x < numCols; x++) { 608 // fprintf(stderr,"psFS64: %d %d %g %g %g\n",x,y,offset,image->data.F64[y][x],log10(image->data.F64[y][x] - offset)); 609 if (isfinite(image->data.F64[y][x])) { 610 copy->data.F64[y][x] = a * asinh( (image->data.F64[y][x] - *boffset) / (2 * *bsoften));// - 2.5 * log10(b) + m0; 611 } 612 else { 613 copy->data.F64[y][x] = image->data.F64[y][x]; 614 } 615 } 616 } 617 break; 618 default: 619 psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Target type is not a float: %d",image->type.type); 620 return NULL; 621 break; 622 } 623 624 // Do regular scaling on the logarithm image 625 if (!scaleStdev(bscale, bzero, copy, mask, maskVal, options)) { 626 psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev"); 627 return false; 628 } 629 psFree(copy); 630 return true; 631 } 632 633 static bool asinhRange(double *bscale, // Scaling, to return 634 double *bzero, // Zero point, to return 635 double *boffset, // asinh flux zeropoint, to return 636 double *bsoften, // asinh softening parameter, to return 637 const psImage *image, // Image to scale 638 const psImage *mask, // Mask image 639 psImageMaskType maskVal, // value to mask 640 const psFitsOptions *options // FITS options 641 ) 642 { 643 psAssert(bscale, "impossible"); 644 psAssert(bzero, "impossible"); 645 psAssert(bsoften, "impossible"); 646 psAssert(image, "impossible"); 647 psAssert(options, "impossible"); 648 649 psTrace("psLib.fits", 3, "Scaling image by asinh method"); 650 int numCols = image->numCols, numRows = image->numRows; // Size of image 651 652 // Measure the mean and stdev 653 // psImageBackground automatically excludes pixels that are non-finite, so we don't need to bother about a 654 // mask. 655 psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); 656 psStats *stats = psStatsAlloc(MEAN_STAT | STDEV_STAT); // Statistics object 657 double mean, stdev; // Mean and standard deviation 658 if (!psImageBackground(stats, NULL, image, mask, maskVal, rng)) { 659 // It could be because the image is entirely masked, in which case we don't want to error 660 bool good = false; // Any good pixels? 661 662 663 // Find good pixels in an image, by image type 664 #define GOOD_PIXELS_CASE(TYPE) \ 665 case PS_TYPE_##TYPE: \ 666 for (int y = 0; y < image->numRows && !good; y++) { \ 667 for (int x = 0; x < image->numCols && !good; x++) { \ 668 if (mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal)) { \ 669 continue; \ 670 } \ 671 if (!isfinite(image->data.TYPE[y][x])) { \ 672 continue; \ 673 } \ 674 good = true; \ 675 } \ 676 } \ 677 break; 678 679 switch (image->type.type) { 680 GOOD_PIXELS_CASE(F32); 681 GOOD_PIXELS_CASE(F64); 682 default: 683 psAbort("Unsupported case: %x", image->type.type); 684 } 685 686 if (!good) { 687 psLogMsg("psLib.fits", PS_LOG_DETAIL, "Image has no good pixels, setting BSCALE = 1, BZERO = 0"); 688 psErrorClear(); 689 *bscale = 1.0; 690 *bzero = 0.0; 691 *bsoften = NAN; 692 psFree(rng); 693 psFree(stats); 694 return true; 695 } 696 697 // There are some good pixels in there somewhere; psImageBackground just didn't find them 698 psLogMsg("psLib.fits", PS_LOG_DETAIL, 699 "Couldn't measure background statistics for image quantisation; retrying."); 700 psErrorClear(); 701 // Retry using all the available pixels 702 stats->nSubsample = image->numCols * image->numRows + 1; 703 if (!psImageStats(stats, image, mask, maskVal)) { 704 psLogMsg("psLib.fits", PS_LOG_DETAIL, 705 "Couldn't measure background statistics for image quantisation (attempt 2); retrying."); 706 psErrorClear(); 707 // Retry with desperate statistic 708 stats->options = DESPERATE_MEAN_STAT | DESPERATE_STDEV_STAT; 709 if (!psImageStats(stats, image, mask, maskVal)) { 710 psError(PS_ERR_UNKNOWN, false, "Unable to measure background statistics for image"); 711 psFree(rng); 712 psFree(stats); 713 return false; 714 } else { 715 // Desperate retry 716 mean = psStatsGetValue(stats, DESPERATE_MEAN_STAT); 717 stdev = psStatsGetValue(stats, DESPERATE_STDEV_STAT); 718 } 719 } else { 720 // Retry with all available pixels 721 mean = psStatsGetValue(stats, MEAN_STAT); 722 stdev = psStatsGetValue(stats, STDEV_STAT); 723 } 724 } else { 725 // First attempt 726 mean = psStatsGetValue(stats, MEAN_STAT); 727 stdev = psStatsGetValue(stats, STDEV_STAT); 728 } 729 psFree(rng); 730 psFree(stats); 731 if (!isfinite(mean) || !isfinite(stdev)) { 732 psError(PS_ERR_BAD_PARAMETER_VALUE, true, 733 "Mean (%f) or stdev (%f) is non-finite.", mean, stdev); 734 return false; 735 } 736 737 psImage *copy; 738 739 switch (image->type.type) { 740 case PS_TYPE_F32: 741 copy = psImageCopy(NULL,image,PS_TYPE_F32); 742 break; 743 case PS_TYPE_F64: 744 copy = psImageCopy(NULL,image,PS_TYPE_F64); 745 break; 746 default: 747 psError(PS_ERR_UNKNOWN, true, "Target type is not a float: %d", image->type.type); 748 return NULL; 749 break; 750 } 751 // Do scaling 752 float a = 1.0857362; // 2.5 * log(e); 753 *bsoften = sqrt(a) * stdev; 754 *boffset = mean; 755 // float m0 = 0; // Can I just arbitrarily set this? 756 757 switch (image->type.type) { 758 case PS_TYPE_F32: 759 for (int y = 0; y < numRows; y++) { 760 for (int x = 0; x < numCols; x++) { 761 /* if (x == 2331 && y == 2843) { */ 762 /* fprintf(stderr,"psFS32: %d %d %g %g %g\n",x,y,offset,image->data.F32[y][x],log10(image->data.F32[y][x] - offset)); */ 763 /* } */ 764 /* if (x == 266 && y == 4584) { */ 765 /* fprintf(stderr,"psFS32: %d %d %g %g %g %g %g\n",x,y,*boffset,*bsoften,image->data.F32[y][x],log10(image->data.F32[y][x] - *boffset), */ 766 /* a * asinh((image->data.F32[y][x] - *boffset) / (2 * *bsoften))); */ 767 /* } */ 768 769 if (isfinite(image->data.F32[y][x])) { 770 copy->data.F32[y][x] = a * asinh( (image->data.F32[y][x] - *boffset) / (2 * *bsoften));// - 2.5 * log10(b) + m0; 771 } 772 else { 773 copy->data.F32[y][x] = image->data.F32[y][x]; 774 } 775 } 776 } 777 break; 778 case PS_TYPE_F64: 779 for (int y = 0; y < numRows; y++) { 780 for (int x = 0; x < numCols; x++) { 781 // fprintf(stderr,"psFS64: %d %d %g %g %g\n",x,y,offset,image->data.F64[y][x],log10(image->data.F64[y][x] - offset)); 782 if (isfinite(image->data.F64[y][x])) { 783 copy->data.F64[y][x] = a * asinh( (image->data.F64[y][x] - *boffset)/ (2 * *bsoften));// - 2.5 * log10(b) + m0; 784 } 785 else { 786 copy->data.F64[y][x] = image->data.F32[y][x]; 787 } 788 } 789 } 790 break; 791 default: 792 psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Target type is not a float: %d",image->type.type); 793 return NULL; 794 break; 795 } 796 797 // Do regular scaling on the asinh image 798 if (!scaleRange(bscale, bzero, copy, options)) { 799 psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev"); 800 return false; 801 } 802 psFree(copy); 803 return true; 804 } 805 458 806 459 807 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// … … 461 809 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 462 810 463 bool psFitsScaleDetermine(double *bscale, double *bzero, double *boffset, long *blank, const psImage *image, 811 bool psFitsScaleDetermine(double *bscale, double *bzero, double *boffset, double *bsoften, 812 long *blank, const psImage *image, 464 813 const psImage *mask, psImageMaskType maskVal, const psFits *fits) 465 814 { … … 521 870 } 522 871 break; 523 case PS_FITS_SCALE_LOG_RANGE: 524 if (!logscaleRange(bscale,bzero,boffset,image,options)) { 525 psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from range"); 526 return false; 527 } 528 break; 529 case PS_FITS_SCALE_LOG_STDEV_POSITIVE: 530 case PS_FITS_SCALE_LOG_STDEV_NEGATIVE: 531 case PS_FITS_SCALE_LOG_STDEV_BOTH: 532 if (!logscaleStdev(bscale, bzero,boffset, image, mask, maskVal, options)) { 533 psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev"); 534 return false; 535 } 536 break; 872 case PS_FITS_SCALE_LOG_RANGE: 873 if (!logscaleRange(bscale,bzero,boffset,image,options)) { 874 psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from range"); 875 return false; 876 } 877 break; 878 case PS_FITS_SCALE_LOG_STDEV_POSITIVE: 879 case PS_FITS_SCALE_LOG_STDEV_NEGATIVE: 880 case PS_FITS_SCALE_LOG_STDEV_BOTH: 881 if (!logscaleStdev(bscale, bzero,boffset, image, mask, maskVal, options)) { 882 psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev"); 883 return false; 884 } 885 break; 886 case PS_FITS_SCALE_ASINH_RANGE: 887 if (!asinhRange(bscale,bzero,boffset,bsoften,image,mask,maskVal,options)) { 888 psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from range"); 889 return false; 890 } 891 break; 892 case PS_FITS_SCALE_ASINH_STDEV_POSITIVE: 893 case PS_FITS_SCALE_ASINH_STDEV_NEGATIVE: 894 case PS_FITS_SCALE_ASINH_STDEV_BOTH: 895 if (!asinhStdev(bscale, bzero,boffset,bsoften, image, mask, maskVal, options)) { 896 psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev"); 897 return false; 898 } 899 break; 900 901 537 902 case PS_FITS_SCALE_MANUAL: 538 903 *bscale = options->bscale; … … 557 922 558 923 559 psTrace("psLib.fits", 3, "BSCALE = %.10lf, BZERO = %.10lf, BOFFSET = %.10lf, B LANK = %ld\n",560 *bscale, *bzero, *boffset, *b lank);924 psTrace("psLib.fits", 3, "BSCALE = %.10lf, BZERO = %.10lf, BOFFSET = %.10lf, BSOFTEN = %.10lf, BLANK = %ld\n", 925 *bscale, *bzero, *boffset, *bsoften, *blank); 561 926 return true; 562 927 } 563 928 564 929 565 psImage *psFitsScaleForDisk(const psImage *image, const psFits *fits, double bscale, double bzero, double boffset, 930 psImage *psFitsScaleForDisk(const psImage *image, const psFits *fits, double bscale, double bzero, double boffset, double bsoften, 566 931 psRandom *rng) 567 932 { … … 625 990 (options->scaling == PS_FITS_SCALE_LOG_STDEV_NEGATIVE)|| \ 626 991 (options->scaling == PS_FITS_SCALE_LOG_STDEV_BOTH)) { \ 627 value = log10( (IN)->data.INTYPE[y][x] - boffset ); \ 992 if (isfinite((IN)->data.INTYPE[y][x])) { \ 993 value = log10( (IN)->data.INTYPE[y][x] - boffset ); \ 994 } \ 995 else { \ 996 value = (IN)->data.INTYPE[y][x]; \ 997 } \ 998 } \ 999 else if ((options->scaling == PS_FITS_SCALE_ASINH_RANGE)|| \ 1000 (options->scaling == PS_FITS_SCALE_ASINH_MANUAL)|| \ 1001 (options->scaling == PS_FITS_SCALE_ASINH_STDEV_POSITIVE)|| \ 1002 (options->scaling == PS_FITS_SCALE_ASINH_STDEV_NEGATIVE)|| \ 1003 (options->scaling == PS_FITS_SCALE_ASINH_STDEV_BOTH)) { \ 1004 if (isfinite((IN)->data.INTYPE[y][x])) { \ 1005 value = 1.0857362 * (asinh( ((IN)->data.INTYPE[y][x] - boffset) / (2.0 * bsoften))); \ 1006 } \ 1007 else { \ 1008 value = (IN)->data.INTYPE[y][x]; \ 1009 } \ 1010 if ((x == 1000)&&(y == 1000)) { \ 1011 fprintf(stderr,"ASINH: %f %f %f", \ 1012 (IN)->data.INTYPE[y][x],value,bsoften); \ 1013 } \ 628 1014 } \ 629 1015 else { \ … … 652 1038 case PS_TYPE_##INTYPE: { \ 653 1039 switch (outType) { \ 654 SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, U8);\655 SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S16);\656 SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S32);\657 SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S64);\1040 SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, U8);; \ 1041 SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S16);; \ 1042 SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S32);; \ 1043 SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S64);; \ 658 1044 default: \ 659 1045 psAbort("Should be unreachable."); \ … … 664 1050 switch (image->type.type) { 665 1051 SCALE_WRITE_IN_CASE(image, F32, out); 666 SCALE_WRITE_IN_CASE(image, F64, out); 1052 SCALE_WRITE_IN_CASE(image, F64, out); 667 1053 default: 668 1054 psAbort("Should be unreachable."); … … 679 1065 // the present time, since cfitsio should apply the scaling itself in the process of reading. However, we may 680 1066 // later desire it (e.g., if we ever make our own FITS implementation). 681 psImage *psFitsScaleFromDisk(const psImage *image, double boffset )1067 psImage *psFitsScaleFromDisk(const psImage *image, double boffset, double bsoften) 682 1068 { 683 1069 PS_ASSERT_IMAGE_NON_NULL(image, NULL); … … 718 1104 for (int y = 0; y < numRows; y++) { \ 719 1105 for (int x = 0; x < numCols; x++) { \ 720 out->data.OUTTYPE[y][x] = pow(10,image->data.INTYPE[y][x]) + boffset;; \ 1106 if (bsoften) { \ 1107 out->data.OUTTYPE[y][x] = 2 * bsoften * sinh(image->data.INTYPE[y][x] / (1.0857362)) + boffset; \ 1108 } \ 1109 else if (boffset) { \ 1110 out->data.OUTTYPE[y][x] = pow(10,image->data.INTYPE[y][x]) + boffset;; \ 1111 } \ 721 1112 } \ 722 1113 } \ … … 766 1157 if (strcasecmp(string, "LOG_STDEV_NEGATIVE") == 0) return PS_FITS_SCALE_LOG_STDEV_NEGATIVE; 767 1158 if (strcasecmp(string, "LOG_STDEV_BOTH") == 0) return PS_FITS_SCALE_LOG_STDEV_BOTH; 1159 if (strcasecmp(string, "ASINH_RANGE") == 0) return PS_FITS_SCALE_ASINH_RANGE; 1160 if (strcasecmp(string, "ASINH_MANUAL") == 0) return PS_FITS_SCALE_ASINH_MANUAL; 1161 if (strcasecmp(string, "ASINH_STDEV_POSITIVE") == 0) return PS_FITS_SCALE_ASINH_STDEV_POSITIVE; 1162 if (strcasecmp(string, "ASINH_STDEV_NEGATIVE") == 0) return PS_FITS_SCALE_ASINH_STDEV_NEGATIVE; 1163 if (strcasecmp(string, "ASINH_STDEV_BOTH") == 0) return PS_FITS_SCALE_ASINH_STDEV_BOTH; 768 1164 if (strcasecmp(string, "MANUAL") == 0) return PS_FITS_SCALE_MANUAL; 769 1165 -
trunk/psLib/src/fits/psFitsScale.h
r30636 r31633 11 11 double *bzero, ///< Zero point, to return 12 12 double *boffset, ///< Log offset, to return 13 double *bsoften, ///< asinh softening parameter, to return 13 14 long *blank, ///< Blank value, to return 14 15 const psImage *image, ///< Image to scale … … 29 30 double bzero, ///< Zero point 30 31 double boffset, ///< Log offset 32 double bsoften, ///< asinh softening parameter 31 33 psRandom *rng ///< Random number generator (for the "fuzz"), or NULL 32 34 ); 33 35 psImage *psFitsScaleFromDisk(const psImage *image, ///< Image to to unapply BOFFSET 34 double boffset ///< Log offset 36 double boffset, ///< Log offset 37 double bsoften ///< asinh softening parameter 35 38 ); 36 39 /// Interpret a string as a scaling method -
trunk/psModules/src/camera/pmFPAfileDefine.c
r30636 r31633 283 283 case PS_FITS_SCALE_RANGE: 284 284 case PS_FITS_SCALE_LOG_RANGE: 285 case PS_FITS_SCALE_ASINH_RANGE: 285 286 // No options required 286 287 break; … … 289 290 case PS_FITS_SCALE_LOG_STDEV_POSITIVE: 290 291 case PS_FITS_SCALE_LOG_STDEV_NEGATIVE: 292 case PS_FITS_SCALE_ASINH_STDEV_POSITIVE: 293 case PS_FITS_SCALE_ASINH_STDEV_NEGATIVE: 291 294 options->stdevNum = parseOptionFloat(scheme, "STDEV.NUM", source); // Padding to edge 292 295 if (!isfinite(options->stdevNum)) { … … 299 302 case PS_FITS_SCALE_STDEV_BOTH: 300 303 case PS_FITS_SCALE_LOG_STDEV_BOTH: 304 case PS_FITS_SCALE_ASINH_STDEV_BOTH: 301 305 options->stdevBits = parseOptionInt(scheme, "STDEV.BITS", source, 0); // Bits for stdev 302 306 if (options->stdevBits <= 0) { … … 312 316 options->bzero = parseOptionDouble(scheme, "BZERO", source); // Zero point 313 317 break; 314 case PS_FITS_SCALE_LOG_MANUAL: 315 options->bscale = parseOptionDouble(scheme, "BSCALE", source); // Scaling 316 options->bzero = parseOptionDouble(scheme, "BZERO", source); // Zero point 317 options->boffset = parseOptionDouble(scheme, "BOFFSET", source); // Log offset 318 break; 319 default: 318 case PS_FITS_SCALE_LOG_MANUAL: 319 options->bscale = parseOptionDouble(scheme, "BSCALE", source); // Scaling 320 options->bzero = parseOptionDouble(scheme, "BZERO", source); // Zero point 321 options->boffset = parseOptionDouble(scheme, "BOFFSET", source); // Log offset 322 case PS_FITS_SCALE_ASINH_MANUAL: 323 options->bscale = parseOptionDouble(scheme, "BSCALE", source); // Scaling 324 options->bzero = parseOptionDouble(scheme, "BZERO", source); // Zero point 325 options->boffset = parseOptionDouble(scheme, "BOFFSET", source); // Log offset 326 options->bsoften = parseOptionDouble(scheme, "BSOFTEN", source); // Softening parameter 327 break; 328 default: 320 329 psAbort("Should never get here."); 321 330 }
Note:
See TracChangeset
for help on using the changeset viewer.
