- Timestamp:
- Jul 17, 2014, 12:36:19 PM (12 years ago)
- Location:
- branches/eam_branches/ipp-ops-20130712
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
psLib/src/fits/psFitsScale.c (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/eam_branches/ipp-ops-20130712
- Property svn:mergeinfo changed
-
branches/eam_branches/ipp-ops-20130712/psLib/src/fits/psFitsScale.c
r35674 r37068 929 929 930 930 931 psImage *psFitsScaleForDisk(const psImage *image, const psFits *fits, double bscale, double bzero, double boffset, double bsoften, 932 psRandom *rng) 931 psImage *psFitsScaleForDisk(const psImage *image, const psFits *fits, double bscale, double bzero, double boffset, double bsoften, psRandom *rng) 933 932 { 934 933 PS_ASSERT_IMAGE_NON_NULL(image, NULL); … … 971 970 psImage *out = psImageAlloc(numCols, numRows, outType); // Output image 972 971 973 if (!psMemIncrRefCounter(rng) && options->fuzz) { 974 // Don't blab about which seed we're going to get --- it's not necessary for this purpose 975 rng = psRandomAlloc(PS_RANDOM_TAUS); 976 } 977 978 #define SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, OUTTYPE) \ 979 case PS_TYPE_##OUTTYPE: { \ 980 double scale = 1.0 / bscale; \ 981 double zero = bzero; \ 972 // fprintf (stderr, "corruption check psFitsScaleForDisk 1\n"); 973 // psMemCheckCorruption (stderr, true); 974 975 if (options->fuzz && !psMemIncrRefCounter (rng)) { 976 // Don't blab about which seed we're going to get --- it's not necessary for this purpose 977 rng = psRandomAlloc(PS_RANDOM_TAUS); 978 } 979 980 #define SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, OUTTYPE) \ 981 case PS_TYPE_##OUTTYPE: { \ 982 double scale = 1.0 / bscale; \ 983 double zero = bzero; \ 982 984 /* Note: BITPIX=8 treated differently, since it uses unsigned values; the rest use signed */ \ 983 double min = bitpix == 8 ? 0 : -pow(2.0, options->bitpix - 1); \985 double min = bitpix == 8 ? 0 : -pow(2.0, options->bitpix - 1); \ 984 986 double max = bitpix == 8 ? 255 : (pow(2.0, options->bitpix - 1) - 1.0); \ 985 for (int y = 0; y < numRows; y++) { \ 986 for (int x = 0; x < numCols; x++) { \ 987 ps##INTYPE value; \ 988 if ((options->scaling == PS_FITS_SCALE_LOG_RANGE)|| \ 989 (options->scaling == PS_FITS_SCALE_LOG_MANUAL)|| \ 990 (options->scaling == PS_FITS_SCALE_LOG_STDEV_POSITIVE)|| \ 991 (options->scaling == PS_FITS_SCALE_LOG_STDEV_NEGATIVE)|| \ 992 (options->scaling == PS_FITS_SCALE_LOG_STDEV_BOTH)) { \ 993 if (isfinite((IN)->data.INTYPE[y][x])) { \ 994 value = log10( (IN)->data.INTYPE[y][x] - boffset ); \ 987 for (int y = 0; y < numRows; y++) { \ 988 for (int x = 0; x < numCols; x++) { \ 989 ps##INTYPE value; \ 990 if ((options->scaling == PS_FITS_SCALE_LOG_RANGE)|| \ 991 (options->scaling == PS_FITS_SCALE_LOG_MANUAL)|| \ 992 (options->scaling == PS_FITS_SCALE_LOG_STDEV_POSITIVE)|| \ 993 (options->scaling == PS_FITS_SCALE_LOG_STDEV_NEGATIVE)|| \ 994 (options->scaling == PS_FITS_SCALE_LOG_STDEV_BOTH)) { \ 995 if (isfinite((IN)->data.INTYPE[y][x])) { \ 996 value = log10( (IN)->data.INTYPE[y][x] - boffset ); \ 997 } \ 998 else { \ 999 value = (IN)->data.INTYPE[y][x]; \ 1000 } \ 1001 } \ 1002 else if ((options->scaling == PS_FITS_SCALE_ASINH_RANGE)|| \ 1003 (options->scaling == PS_FITS_SCALE_ASINH_MANUAL)|| \ 1004 (options->scaling == PS_FITS_SCALE_ASINH_STDEV_POSITIVE)|| \ 1005 (options->scaling == PS_FITS_SCALE_ASINH_STDEV_NEGATIVE)|| \ 1006 (options->scaling == PS_FITS_SCALE_ASINH_STDEV_BOTH)) { \ 1007 if (isfinite((IN)->data.INTYPE[y][x])) { \ 1008 value = 1.0857362 * (asinh( ((IN)->data.INTYPE[y][x] - boffset) / (2.0 * bsoften))); \ 1009 } \ 1010 else { \ 1011 value = (IN)->data.INTYPE[y][x]; \ 1012 } \ 995 1013 } \ 996 1014 else { \ 997 value = (IN)->data.INTYPE[y][x]; \1015 value = (IN)->data.INTYPE[y][x]; \ 998 1016 } \ 999 } \ 1000 else if ((options->scaling == PS_FITS_SCALE_ASINH_RANGE)|| \ 1001 (options->scaling == PS_FITS_SCALE_ASINH_MANUAL)|| \ 1002 (options->scaling == PS_FITS_SCALE_ASINH_STDEV_POSITIVE)|| \ 1003 (options->scaling == PS_FITS_SCALE_ASINH_STDEV_NEGATIVE)|| \ 1004 (options->scaling == PS_FITS_SCALE_ASINH_STDEV_BOTH)) { \ 1005 if (isfinite((IN)->data.INTYPE[y][x])) { \ 1006 value = 1.0857362 * (asinh( ((IN)->data.INTYPE[y][x] - boffset) / (2.0 * bsoften))); \ 1017 if (!isfinite(value)) { \ 1018 /* This choice of "max" for non-finite pixels is mainly cosmetic --- it has to be */ \ 1019 /* something, and "min" would produce holes in the cores of bright stars. */ \ 1020 (OUT)->data.OUTTYPE[y][x] = max; \ 1021 } else { \ 1022 value = (value - zero) * scale; \ 1023 if (options->fuzz && (value - (int)value != 0.0)) { \ 1024 /* Add random factor [0.0,1.0): adds a variance of 1/12, */ \ 1025 /* but preserves the expectation value given the floor() */ \ 1026 value += psRandomUniform(rng) ; \ 1027 } \ 1028 /* Check for underflow and overflow; set either to max */ \ 1029 (OUT)->data.OUTTYPE[y][x] = (value < min || value > max ? max : floor(value)); \ 1007 1030 } \ 1008 else { \ 1009 value = (IN)->data.INTYPE[y][x]; \ 1010 } \ 1011 if ((x == 1000)&&(y == 1000)) { \ 1012 fprintf(stderr,"ASINH: %f %f %f", \ 1013 (IN)->data.INTYPE[y][x],value,bsoften); \ 1014 } \ 1015 } \ 1016 else { \ 1017 value = (IN)->data.INTYPE[y][x]; \ 1018 } \ 1019 if (!isfinite(value)) { \ 1020 /* This choice of "max" for non-finite pixels is mainly cosmetic --- it has to be */ \ 1021 /* something, and "min" would produce holes in the cores of bright stars. */ \ 1022 (OUT)->data.OUTTYPE[y][x] = max; \ 1023 } else { \ 1024 value = (value - zero) * scale; \ 1025 if (options->fuzz && (value - (int)value != 0.0)) { \ 1026 /* Add random factor [-0.5,0.5): adds a variance of 1/12, */ \ 1027 /* but preserves the expectation value */ \ 1028 value += psRandomUniform(rng) - 0.5; \ 1029 } \ 1030 /* Check for underflow and overflow; set either to max */ \ 1031 (OUT)->data.OUTTYPE[y][x] = (value < min || value > max ? max : value); \ 1032 } \ 1033 } \ 1034 } \ 1035 break; \ 1036 } 1037 1038 #define SCALE_WRITE_IN_CASE(IN, INTYPE, OUT) \ 1039 case PS_TYPE_##INTYPE: { \ 1040 switch (outType) { \ 1041 SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, U8);; \ 1042 SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S16);; \ 1043 SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S32);; \ 1044 SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S64);; \ 1045 default: \ 1046 psAbort("Should be unreachable."); \ 1047 } \ 1048 break; \ 1049 } 1050 1031 } \ 1032 } \ 1033 break; \ 1034 } 1035 1036 #define SCALE_WRITE_IN_CASE(IN, INTYPE, OUT) \ 1037 case PS_TYPE_##INTYPE: { \ 1038 switch (outType) { \ 1039 SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, U8);; \ 1040 SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S16);; \ 1041 SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S32);; \ 1042 SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S64);; \ 1043 default: \ 1044 psAbort("Should be unreachable."); \ 1045 } \ 1046 break; \ 1047 } 1048 1051 1049 switch (image->type.type) { 1052 SCALE_WRITE_IN_CASE(image, F32, out);1053 SCALE_WRITE_IN_CASE(image, F64, out);1050 SCALE_WRITE_IN_CASE(image, F32, out); 1051 SCALE_WRITE_IN_CASE(image, F64, out); 1054 1052 default: 1055 1053 psAbort("Should be unreachable."); … … 1057 1055 1058 1056 psFree(rng); 1059 1060 1057 return out; 1061 1058 }
Note:
See TracChangeset
for help on using the changeset viewer.
