- Timestamp:
- Jun 8, 2011, 2:41:13 PM (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/czw_branch/20110406/psLib/src/fits/psFitsScale.c
r30636 r31618 242 242 243 243 244 245 244 246 static bool logscaleStdev(double *bscale, // Scaling, to return 245 247 double *bzero, // Zero point, to return … … 456 458 } 457 459 460 static bool asinhStdev(double *bscale, // Scaling, to return 461 double *bzero, // Zero point, to return 462 double *bsoften, // asinh softening parameter, to return 463 const psImage *image, // Image to scale 464 const psImage *mask, // Mask image 465 psImageMaskType maskVal, // value to mask 466 const psFitsOptions *options // FITS options 467 ) 468 { 469 psAssert(bscale, "impossible"); 470 psAssert(bzero, "impossible"); 471 psAssert(bsoften, "impossible"); 472 psAssert(image, "impossible"); 473 psAssert(options, "impossible"); 474 475 psTrace("psLib.fits", 3, "Scaling image by asinh method"); 476 int numCols = image->numCols, numRows = image->numRows; // Size of image 477 478 // Measure the mean and stdev 479 // psImageBackground automatically excludes pixels that are non-finite, so we don't need to bother about a 480 // mask. 481 psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); 482 psStats *stats = psStatsAlloc(MEAN_STAT | STDEV_STAT); // Statistics object 483 double mean, stdev; // Mean and standard deviation 484 if (!psImageBackground(stats, NULL, image, mask, maskVal, rng)) { 485 // It could be because the image is entirely masked, in which case we don't want to error 486 bool good = false; // Any good pixels? 487 488 489 // Find good pixels in an image, by image type 490 #define GOOD_PIXELS_CASE(TYPE) \ 491 case PS_TYPE_##TYPE: \ 492 for (int y = 0; y < image->numRows && !good; y++) { \ 493 for (int x = 0; x < image->numCols && !good; x++) { \ 494 if (mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal)) { \ 495 continue; \ 496 } \ 497 if (!isfinite(image->data.TYPE[y][x])) { \ 498 continue; \ 499 } \ 500 good = true; \ 501 } \ 502 } \ 503 break; 504 505 switch (image->type.type) { 506 GOOD_PIXELS_CASE(F32); 507 GOOD_PIXELS_CASE(F64); 508 default: 509 psAbort("Unsupported case: %x", image->type.type); 510 } 511 512 if (!good) { 513 psLogMsg("psLib.fits", PS_LOG_DETAIL, "Image has no good pixels, setting BSCALE = 1, BZERO = 0"); 514 psErrorClear(); 515 *bscale = 1.0; 516 *bzero = 0.0; 517 *bsoften = NAN; 518 psFree(rng); 519 psFree(stats); 520 return true; 521 } 522 523 // There are some good pixels in there somewhere; psImageBackground just didn't find them 524 psLogMsg("psLib.fits", PS_LOG_DETAIL, 525 "Couldn't measure background statistics for image quantisation; retrying."); 526 psErrorClear(); 527 // Retry using all the available pixels 528 stats->nSubsample = image->numCols * image->numRows + 1; 529 if (!psImageStats(stats, image, mask, maskVal)) { 530 psLogMsg("psLib.fits", PS_LOG_DETAIL, 531 "Couldn't measure background statistics for image quantisation (attempt 2); retrying."); 532 psErrorClear(); 533 // Retry with desperate statistic 534 stats->options = DESPERATE_MEAN_STAT | DESPERATE_STDEV_STAT; 535 if (!psImageStats(stats, image, mask, maskVal)) { 536 psError(PS_ERR_UNKNOWN, false, "Unable to measure background statistics for image"); 537 psFree(rng); 538 psFree(stats); 539 return false; 540 } else { 541 // Desperate retry 542 mean = psStatsGetValue(stats, DESPERATE_MEAN_STAT); 543 stdev = psStatsGetValue(stats, DESPERATE_STDEV_STAT); 544 } 545 } else { 546 // Retry with all available pixels 547 mean = psStatsGetValue(stats, MEAN_STAT); 548 stdev = psStatsGetValue(stats, STDEV_STAT); 549 } 550 } else { 551 // First attempt 552 mean = psStatsGetValue(stats, MEAN_STAT); 553 stdev = psStatsGetValue(stats, STDEV_STAT); 554 } 555 psFree(rng); 556 psFree(stats); 557 if (!isfinite(mean) || !isfinite(stdev)) { 558 psError(PS_ERR_BAD_PARAMETER_VALUE, true, 559 "Mean (%f) or stdev (%f) is non-finite.", mean, stdev); 560 return false; 561 } 562 563 psImage *copy; 564 565 switch (image->type.type) { 566 case PS_TYPE_F32: 567 copy = psImageCopy(NULL,image,PS_TYPE_F32); 568 break; 569 case PS_TYPE_F64: 570 copy = psImageCopy(NULL,image,PS_TYPE_F64); 571 break; 572 default: 573 psError(PS_ERR_UNKNOWN, true, "Target type is not a float: %d", image->type.type); 574 return NULL; 575 break; 576 } 577 // Do scaling 578 float a = 1.0857362; // 2.5 * log(e); 579 *bsoften = sqrt(a) * stdev; 580 // float m0 = 0; // Can I just arbitrarily set this? 581 582 switch (image->type.type) { 583 case PS_TYPE_F32: 584 for (int y = 0; y < numRows; y++) { 585 for (int x = 0; x < numCols; x++) { 586 /* if (x == 2331 && y == 2843) { */ 587 /* fprintf(stderr,"psFS32: %d %d %g %g %g\n",x,y,offset,image->data.F32[y][x],log10(image->data.F32[y][x] - offset)); */ 588 /* } */ 589 copy->data.F32[y][x] = -1.0 * a * asinh( image->data.F32[y][x] / (2 * *bsoften));// - 2.5 * log10(b) + m0; 590 } 591 } 592 break; 593 case PS_TYPE_F64: 594 for (int y = 0; y < numRows; y++) { 595 for (int x = 0; x < numCols; x++) { 596 // fprintf(stderr,"psFS64: %d %d %g %g %g\n",x,y,offset,image->data.F64[y][x],log10(image->data.F64[y][x] - offset)); 597 copy->data.F32[y][x] = -1.0 * a * asinh( image->data.F32[y][x] / (2 * *bsoften));// - 2.5 * log10(b) + m0; 598 } 599 } 600 break; 601 default: 602 psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Target type is not a float: %d",image->type.type); 603 return NULL; 604 break; 605 } 606 607 // Do regular scaling on the logarithm image 608 if (!scaleStdev(bscale, bzero, copy, mask, maskVal, options)) { 609 psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev"); 610 return false; 611 } 612 psFree(copy); 613 return true; 614 } 615 616 static bool asinhRange(double *bscale, // Scaling, to return 617 double *bzero, // Zero point, to return 618 double *bsoften, // asinh softening parameter, to return 619 const psImage *image, // Image to scale 620 const psImage *mask, // Mask image 621 psImageMaskType maskVal, // value to mask 622 const psFitsOptions *options // FITS options 623 ) 624 { 625 psAssert(bscale, "impossible"); 626 psAssert(bzero, "impossible"); 627 psAssert(bsoften, "impossible"); 628 psAssert(image, "impossible"); 629 psAssert(options, "impossible"); 630 631 psTrace("psLib.fits", 3, "Scaling image by asinh method"); 632 int numCols = image->numCols, numRows = image->numRows; // Size of image 633 634 // Measure the mean and stdev 635 // psImageBackground automatically excludes pixels that are non-finite, so we don't need to bother about a 636 // mask. 637 psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); 638 psStats *stats = psStatsAlloc(MEAN_STAT | STDEV_STAT); // Statistics object 639 double mean, stdev; // Mean and standard deviation 640 if (!psImageBackground(stats, NULL, image, mask, maskVal, rng)) { 641 // It could be because the image is entirely masked, in which case we don't want to error 642 bool good = false; // Any good pixels? 643 644 645 // Find good pixels in an image, by image type 646 #define GOOD_PIXELS_CASE(TYPE) \ 647 case PS_TYPE_##TYPE: \ 648 for (int y = 0; y < image->numRows && !good; y++) { \ 649 for (int x = 0; x < image->numCols && !good; x++) { \ 650 if (mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal)) { \ 651 continue; \ 652 } \ 653 if (!isfinite(image->data.TYPE[y][x])) { \ 654 continue; \ 655 } \ 656 good = true; \ 657 } \ 658 } \ 659 break; 660 661 switch (image->type.type) { 662 GOOD_PIXELS_CASE(F32); 663 GOOD_PIXELS_CASE(F64); 664 default: 665 psAbort("Unsupported case: %x", image->type.type); 666 } 667 668 if (!good) { 669 psLogMsg("psLib.fits", PS_LOG_DETAIL, "Image has no good pixels, setting BSCALE = 1, BZERO = 0"); 670 psErrorClear(); 671 *bscale = 1.0; 672 *bzero = 0.0; 673 *bsoften = NAN; 674 psFree(rng); 675 psFree(stats); 676 return true; 677 } 678 679 // There are some good pixels in there somewhere; psImageBackground just didn't find them 680 psLogMsg("psLib.fits", PS_LOG_DETAIL, 681 "Couldn't measure background statistics for image quantisation; retrying."); 682 psErrorClear(); 683 // Retry using all the available pixels 684 stats->nSubsample = image->numCols * image->numRows + 1; 685 if (!psImageStats(stats, image, mask, maskVal)) { 686 psLogMsg("psLib.fits", PS_LOG_DETAIL, 687 "Couldn't measure background statistics for image quantisation (attempt 2); retrying."); 688 psErrorClear(); 689 // Retry with desperate statistic 690 stats->options = DESPERATE_MEAN_STAT | DESPERATE_STDEV_STAT; 691 if (!psImageStats(stats, image, mask, maskVal)) { 692 psError(PS_ERR_UNKNOWN, false, "Unable to measure background statistics for image"); 693 psFree(rng); 694 psFree(stats); 695 return false; 696 } else { 697 // Desperate retry 698 mean = psStatsGetValue(stats, DESPERATE_MEAN_STAT); 699 stdev = psStatsGetValue(stats, DESPERATE_STDEV_STAT); 700 } 701 } else { 702 // Retry with all available pixels 703 mean = psStatsGetValue(stats, MEAN_STAT); 704 stdev = psStatsGetValue(stats, STDEV_STAT); 705 } 706 } else { 707 // First attempt 708 mean = psStatsGetValue(stats, MEAN_STAT); 709 stdev = psStatsGetValue(stats, STDEV_STAT); 710 } 711 psFree(rng); 712 psFree(stats); 713 if (!isfinite(mean) || !isfinite(stdev)) { 714 psError(PS_ERR_BAD_PARAMETER_VALUE, true, 715 "Mean (%f) or stdev (%f) is non-finite.", mean, stdev); 716 return false; 717 } 718 719 psImage *copy; 720 721 switch (image->type.type) { 722 case PS_TYPE_F32: 723 copy = psImageCopy(NULL,image,PS_TYPE_F32); 724 break; 725 case PS_TYPE_F64: 726 copy = psImageCopy(NULL,image,PS_TYPE_F64); 727 break; 728 default: 729 psError(PS_ERR_UNKNOWN, true, "Target type is not a float: %d", image->type.type); 730 return NULL; 731 break; 732 } 733 // Do scaling 734 float a = 1.0857362; // 2.5 * log(e); 735 *bsoften = sqrt(a) * stdev; 736 // float m0 = 0; // Can I just arbitrarily set this? 737 738 switch (image->type.type) { 739 case PS_TYPE_F32: 740 for (int y = 0; y < numRows; y++) { 741 for (int x = 0; x < numCols; x++) { 742 /* if (x == 2331 && y == 2843) { */ 743 /* fprintf(stderr,"psFS32: %d %d %g %g %g\n",x,y,offset,image->data.F32[y][x],log10(image->data.F32[y][x] - offset)); */ 744 /* } */ 745 copy->data.F32[y][x] = -1.0 * a * asinh( image->data.F32[y][x] / (2 * *bsoften));// - 2.5 * log10(b) + m0; 746 } 747 } 748 break; 749 case PS_TYPE_F64: 750 for (int y = 0; y < numRows; y++) { 751 for (int x = 0; x < numCols; x++) { 752 // fprintf(stderr,"psFS64: %d %d %g %g %g\n",x,y,offset,image->data.F64[y][x],log10(image->data.F64[y][x] - offset)); 753 copy->data.F32[y][x] = -1.0 * a * asinh( image->data.F32[y][x] / (2 * *bsoften));// - 2.5 * log10(b) + m0; 754 } 755 } 756 break; 757 default: 758 psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Target type is not a float: %d",image->type.type); 759 return NULL; 760 break; 761 } 762 763 // Do regular scaling on the logarithm image 764 if (!scaleRange(bscale, bzero, copy, options)) { 765 psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev"); 766 return false; 767 } 768 psFree(copy); 769 return true; 770 } 771 458 772 459 773 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// … … 461 775 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 462 776 463 bool psFitsScaleDetermine(double *bscale, double *bzero, double *boffset, long *blank, const psImage *image, 777 bool psFitsScaleDetermine(double *bscale, double *bzero, double *boffset, double *bsoften, 778 long *blank, const psImage *image, 464 779 const psImage *mask, psImageMaskType maskVal, const psFits *fits) 465 780 { … … 521 836 } 522 837 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; 838 case PS_FITS_SCALE_LOG_RANGE: 839 if (!logscaleRange(bscale,bzero,boffset,image,options)) { 840 psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from range"); 841 return false; 842 } 843 break; 844 case PS_FITS_SCALE_LOG_STDEV_POSITIVE: 845 case PS_FITS_SCALE_LOG_STDEV_NEGATIVE: 846 case PS_FITS_SCALE_LOG_STDEV_BOTH: 847 if (!logscaleStdev(bscale, bzero,boffset, image, mask, maskVal, options)) { 848 psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev"); 849 return false; 850 } 851 break; 852 case PS_FITS_SCALE_ASINH_RANGE: 853 if (!asinhRange(bscale,bzero,bsoften,image,mask,maskVal,options)) { 854 psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from range"); 855 return false; 856 } 857 break; 858 case PS_FITS_SCALE_ASINH_STDEV_POSITIVE: 859 case PS_FITS_SCALE_ASINH_STDEV_NEGATIVE: 860 case PS_FITS_SCALE_ASINH_STDEV_BOTH: 861 if (!asinhStdev(bscale, bzero,bsoften, image, mask, maskVal, options)) { 862 psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev"); 863 return false; 864 } 865 break; 866 867 537 868 case PS_FITS_SCALE_MANUAL: 538 869 *bscale = options->bscale; … … 563 894 564 895 565 psImage *psFitsScaleForDisk(const psImage *image, const psFits *fits, double bscale, double bzero, double boffset, 896 psImage *psFitsScaleForDisk(const psImage *image, const psFits *fits, double bscale, double bzero, double boffset, double bsoften, 566 897 psRandom *rng) 567 898 { … … 626 957 (options->scaling == PS_FITS_SCALE_LOG_STDEV_BOTH)) { \ 627 958 value = log10( (IN)->data.INTYPE[y][x] - boffset ); \ 959 } \ 960 else if ((options->scaling == PS_FITS_SCALE_ASINH_RANGE)|| \ 961 (options->scaling == PS_FITS_SCALE_ASINH_MANUAL)|| \ 962 (options->scaling == PS_FITS_SCALE_ASINH_STDEV_POSITIVE)|| \ 963 (options->scaling == PS_FITS_SCALE_ASINH_STDEV_NEGATIVE)|| \ 964 (options->scaling == PS_FITS_SCALE_ASINH_STDEV_BOTH)) { \ 965 value = -1.0 * 1.0857362 * (asinh( (IN)->data.INTYPE[y][x] - bsoften)); \ 628 966 } \ 629 967 else { \ … … 652 990 case PS_TYPE_##INTYPE: { \ 653 991 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);\992 SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, U8);; \ 993 SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S16);; \ 994 SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S32);; \ 995 SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S64);; \ 658 996 default: \ 659 997 psAbort("Should be unreachable."); \ … … 664 1002 switch (image->type.type) { 665 1003 SCALE_WRITE_IN_CASE(image, F32, out); 666 SCALE_WRITE_IN_CASE(image, F64, out); 1004 SCALE_WRITE_IN_CASE(image, F64, out); 667 1005 default: 668 1006 psAbort("Should be unreachable."); … … 679 1017 // the present time, since cfitsio should apply the scaling itself in the process of reading. However, we may 680 1018 // later desire it (e.g., if we ever make our own FITS implementation). 681 psImage *psFitsScaleFromDisk(const psImage *image, double boffset )1019 psImage *psFitsScaleFromDisk(const psImage *image, double boffset, double bsoften) 682 1020 { 683 1021 PS_ASSERT_IMAGE_NON_NULL(image, NULL); … … 718 1056 for (int y = 0; y < numRows; y++) { \ 719 1057 for (int x = 0; x < numCols; x++) { \ 720 out->data.OUTTYPE[y][x] = pow(10,image->data.INTYPE[y][x]) + boffset;; \ 1058 if (boffset) { \ 1059 out->data.OUTTYPE[y][x] = pow(10,image->data.INTYPE[y][x]) + boffset;; \ 1060 } \ 1061 else if (bsoften) { \ 1062 out->data.OUTTYPE[y][x] = sinh(image->data.INTYPE[y][x] * 2 * bsoften) / (-1.0 * 1.0857362);; \ 1063 } \ 721 1064 } \ 722 1065 } \ … … 766 1109 if (strcasecmp(string, "LOG_STDEV_NEGATIVE") == 0) return PS_FITS_SCALE_LOG_STDEV_NEGATIVE; 767 1110 if (strcasecmp(string, "LOG_STDEV_BOTH") == 0) return PS_FITS_SCALE_LOG_STDEV_BOTH; 1111 if (strcasecmp(string, "ASINH_RANGE") == 0) return PS_FITS_SCALE_ASINH_RANGE; 1112 if (strcasecmp(string, "ASINH_MANUAL") == 0) return PS_FITS_SCALE_ASINH_MANUAL; 1113 if (strcasecmp(string, "ASINH_STDEV_POSITIVE") == 0) return PS_FITS_SCALE_ASINH_STDEV_POSITIVE; 1114 if (strcasecmp(string, "ASINH_STDEV_NEGATIVE") == 0) return PS_FITS_SCALE_ASINH_STDEV_NEGATIVE; 1115 if (strcasecmp(string, "ASINH_STDEV_BOTH") == 0) return PS_FITS_SCALE_ASINH_STDEV_BOTH; 768 1116 if (strcasecmp(string, "MANUAL") == 0) return PS_FITS_SCALE_MANUAL; 769 1117
Note:
See TracChangeset
for help on using the changeset viewer.
