Index: /trunk/ippconfig/recipes/fitstypes.mdc =================================================================== --- /trunk/ippconfig/recipes/fitstypes.mdc (revision 31632) +++ /trunk/ippconfig/recipes/fitstypes.mdc (revision 31633) @@ -101,7 +101,8 @@ COMP_STACK METADATA BITPIX S32 16 - SCALING STR LOG_STDEV_POSITIVE - STDEV.BITS S32 4 - STDEV.NUM F32 10 + SCALING STR ASINH_RANGE +# SCALING STR LOG_RANGE +# STDEV.BITS S32 4 +# STDEV.NUM F32 100 COMPRESSION STR RICE TILE.X S32 0 Index: /trunk/ippconfig/recipes/ppStack.config =================================================================== --- /trunk/ippconfig/recipes/ppStack.config (revision 31632) +++ /trunk/ippconfig/recipes/ppStack.config (revision 31633) @@ -132,5 +132,9 @@ END -# +# Stack mode for testing stacking procedure +TESTSTACK METADATA + VARIANCE BOOL FALSE +END + STACK_DEEP METADATA Index: /trunk/ippconfig/recipes/reductionClasses.mdc =================================================================== --- /trunk/ippconfig/recipes/reductionClasses.mdc (revision 31632) +++ /trunk/ippconfig/recipes/reductionClasses.mdc (revision 31633) @@ -284,4 +284,10 @@ END +TESTSTACK METADATA + STACK_PPSTACK STR TESTSTACK + STACK_PPSUB STR STACK + STACK_PSPHOT STR STACK +END + # basic science analysis PHOTFEST.SINGLE1 METADATA Index: /trunk/psLib/src/fits/psFits.h =================================================================== --- /trunk/psLib/src/fits/psFits.h (revision 31632) +++ /trunk/psLib/src/fits/psFits.h (revision 31633) @@ -55,5 +55,10 @@ PS_FITS_SCALE_LOG_STDEV_NEGATIVE, ///< Take logarithm, Auto-scale to sample stdev, place mean at upper limit PS_FITS_SCALE_LOG_STDEV_BOTH, ///< Take logarithm, Auto-scale to sample stdev, place mean at middle - PS_FITS_SCALE_LOG_MANUAL ///< Manual scaling (use specified BSCALE, BZERO, and BOFFSET) + PS_FITS_SCALE_LOG_MANUAL, ///< Manual scaling (use specified BSCALE, BZERO, and BOFFSET) + PS_FITS_SCALE_ASINH_RANGE, ///< Do asinh scaling, auto-scale to preserve dynamic range + PS_FITS_SCALE_ASINH_STDEV_POSITIVE, ///< Do asinh scaling, auto-scale to sample stdev, place mean at lower limit + PS_FITS_SCALE_ASINH_STDEV_NEGATIVE, ///< Do asinh scaling, auto-scale to sample stdev, place mean at upper limit + PS_FITS_SCALE_ASINH_STDEV_BOTH, ///< Do asinh scaling, auto-scale to sample stdev, place mean at middle + PS_FITS_SCALE_ASINH_MANUAL ///< Manual scaling after doing asinh scaling. } psFitsScaling; @@ -72,4 +77,5 @@ double bscale, bzero; ///< Manually specified BSCALE and BZERO (for SCALE_MANUAL) double boffset; ///< Manually specified BOFFSET (for SCALE_MANUAL) + double bsoften; ///< Manuallly specified BSOFTEN (for SCALE_MANUAL) double mean, stdev; ///< Mean and standard deviation of image int stdevBits; ///< Number of bits to sample a standard deviation (for SCALE_STDEV_*) Index: /trunk/psLib/src/fits/psFitsImage.c =================================================================== --- /trunk/psLib/src/fits/psFitsImage.c (revision 31632) +++ /trunk/psLib/src/fits/psFitsImage.c (revision 31633) @@ -55,4 +55,5 @@ int psDatatype; // psLib data type bool is_logscaled; // is this image log scaled using BOFFSET? + bool is_asinh; // is this image asinh scaled using BSOFTEN? } p_psFitsReadInfo; @@ -129,5 +130,5 @@ // Check scale and zero - double bscale = 0.0, bzero = 0.0, boffset = NAN; // Scale and zero point + double bscale = 0.0, bzero = 0.0, boffset = NAN, bsoften = NAN; // Scale and zero point if (fits_read_key_dbl(fits->fd, "BSCALE", &bscale, NULL, &status) && status != KEY_NO_EXIST) { psFitsError(status, true, "Unable to read header."); @@ -153,4 +154,20 @@ info->is_logscaled = false; } + status = 0; + if (fits_read_key_dbl(fits->fd, "BSOFTEN", &bsoften, NULL, &status) && status != KEY_NO_EXIST) { + psFitsError(status, true, "Unable to read header."); + goto bad; + } + if (status == KEY_NO_EXIST) { + info->is_asinh = false; + } + else if (isfinite(bsoften)) { + info->is_asinh = true; + info->is_logscaled = false; + } + else { + info->is_asinh = false; + } + status = 0; @@ -262,4 +279,5 @@ double *bzero, // Zero point applied double *boffset, // Log offset applied + double *bsoften, // asinh offset applied long *blank, // Blank value (integer data) psFitsFloat *floatType, // Type of custom floating-point @@ -275,4 +293,5 @@ psAssert(bzero, "impossible"); psAssert(boffset, "impossible"); + psAssert(bsoften, "impossible"); psAssert(floatType, "impossible"); psAssert(fits, "impossible"); @@ -322,5 +341,5 @@ if (newScaleZero) { // Choose an appropriate BSCALE and BZERO - if (!psFitsScaleDetermine(bscale, bzero, boffset, blank, image, mask, maskVal, fits)) { + if (!psFitsScaleDetermine(bscale, bzero, boffset, bsoften, blank, image, mask, maskVal, fits)) { // We can't have the write dying for this reason --- try to save it somehow! psWarning("Unable to determine BSCALE and BZERO for image --- refusing to quantise."); @@ -348,5 +367,5 @@ } - return psFitsScaleForDisk(image, fits, *bscale, *bzero, *boffset, rng); + return psFitsScaleForDisk(image, fits, *bscale, *bzero, *boffset, *bsoften, rng); } @@ -480,6 +499,6 @@ outImage = psFitsFloatImageFromDisk(outImage, inImage, floatType); } - // Need to apply BOFFSET if info->is_logscaled is true + if (info->is_logscaled) { double boffset; @@ -487,5 +506,15 @@ status = 0; fits_read_key_dbl(fits->fd, "BOFFSET", &boffset, NULL, &status); - psImage *newImage = psFitsScaleFromDisk(outImage,boffset); + psImage *newImage = psFitsScaleFromDisk(outImage,boffset,0.0); + psFree (outImage); + outImage = newImage; + } + else if (info->is_asinh) { + double bsoften,boffset; + int status; + status = 0; + fits_read_key_dbl(fits->fd, "BSOFTEN", &bsoften, NULL, &status); + fits_read_key_dbl(fits->fd, "BOFFSET", &boffset, NULL, &status); + psImage *newImage = psFitsScaleFromDisk(outImage,boffset,bsoften); psFree (outImage); outImage = newImage; @@ -602,7 +631,8 @@ double bscale = NAN, bzero = NAN; // Scale and zero point to put in header double boffset = NAN; // Log offset to put into header. + double bsoften = NAN; // Asinh softening parameter to put into header long blank = 0; // Blank (undefined) value for image psFitsFloat floatType; // Custom floating-point convention type - psImage *diskImage = imageToDiskRepresentation(&bscale, &bzero, &boffset, &blank, &floatType, fits, image, + psImage *diskImage = imageToDiskRepresentation(&bscale, &bzero, &boffset, &bsoften, &blank, &floatType, fits, image, mask, maskVal, NULL, true); // Image to write out if (!diskImage) { @@ -696,9 +726,24 @@ if (options&&(!((options->scaling == PS_FITS_SCALE_LOG_RANGE)|| (options->scaling == PS_FITS_SCALE_LOG_MANUAL)|| - (options->scaling == PS_FITS_SCALE_LOG_STDEV_POSITIVE)|| - (options->scaling == PS_FITS_SCALE_LOG_STDEV_NEGATIVE)|| - (options->scaling == PS_FITS_SCALE_LOG_STDEV_BOTH)))) { + (options->scaling == PS_FITS_SCALE_LOG_STDEV_POSITIVE)|| + (options->scaling == PS_FITS_SCALE_LOG_STDEV_NEGATIVE)|| + (options->scaling == PS_FITS_SCALE_LOG_STDEV_BOTH)|| + (options->scaling == PS_FITS_SCALE_ASINH_RANGE)|| + (options->scaling == PS_FITS_SCALE_ASINH_MANUAL)|| + (options->scaling == PS_FITS_SCALE_ASINH_STDEV_POSITIVE)|| + (options->scaling == PS_FITS_SCALE_ASINH_STDEV_NEGATIVE)|| + (options->scaling == PS_FITS_SCALE_ASINH_STDEV_BOTH)))) { if (psMetadataLookup(header,"BOFFSET")) { psMetadataRemoveKey(header,"BOFFSET"); + } + } + // Remove any BSOFTENvalues that exist in the header if we are not using that scaling anymore + if (options&&(!((options->scaling == PS_FITS_SCALE_ASINH_RANGE)|| + (options->scaling == PS_FITS_SCALE_ASINH_MANUAL)|| + (options->scaling == PS_FITS_SCALE_ASINH_STDEV_POSITIVE)|| + (options->scaling == PS_FITS_SCALE_ASINH_STDEV_NEGATIVE)|| + (options->scaling == PS_FITS_SCALE_ASINH_STDEV_BOTH)))) { + if (psMetadataLookup(header,"BSOFTEN")) { + psMetadataRemoveKey(header,"BSOFTEN"); } } @@ -727,6 +772,19 @@ (options->scaling == PS_FITS_SCALE_LOG_STDEV_BOTH)))) { fits_write_key_dbl(fits->fd, "BOFFSET", boffset, 12, - "Scaling: TRUE = BZERO + BSCALE * 10**(DISK) + BOFFSET)", &status); + "Scaling: ORIGINAL = 10**(TRUE) + BOFFSET)", &status); } + if (options&&(((options->scaling == PS_FITS_SCALE_ASINH_RANGE)|| + (options->scaling == PS_FITS_SCALE_ASINH_MANUAL)|| + (options->scaling == PS_FITS_SCALE_ASINH_STDEV_POSITIVE)|| + (options->scaling == PS_FITS_SCALE_ASINH_STDEV_NEGATIVE)|| + (options->scaling == PS_FITS_SCALE_ASINH_STDEV_BOTH)))) { + fits_write_key_dbl(fits->fd, "BOFFSET", boffset, 12, + "Scaling: NORM = DISK + BOFFSET", &status); + fits_write_key_dbl(fits->fd, "BSOFTEN", bsoften, 12, + "Scaling: TRUE = 2 * BSOFTEN * sinh(NORM/a)",&status); + + + } + if (psFitsError(status, true, "Could not write BSCALE/BZERO headers to file.")) { success = false; @@ -831,7 +889,8 @@ double bscale = NAN, bzero = NAN; // Scale and zero point to put in header double boffset = NAN; // Log offset to put in header + double bsoften = NAN; // Asinh softening parameter to put in header long blank = 0; // Blank (undefined) value for image psFitsFloat floatType; // Custom floating-point convention type - psImage *diskImage = imageToDiskRepresentation(&bscale, &bzero, &boffset, &blank, &floatType, fits, input, + psImage *diskImage = imageToDiskRepresentation(&bscale, &bzero, &boffset, &bsoften, &blank, &floatType, fits, input, mask, maskVal, NULL, false); // Image to write out if (!diskImage) { Index: /trunk/psLib/src/fits/psFitsScale.c =================================================================== --- /trunk/psLib/src/fits/psFitsScale.c (revision 31632) +++ /trunk/psLib/src/fits/psFitsScale.c (revision 31633) @@ -214,4 +214,5 @@ case PS_FITS_SCALE_STDEV_POSITIVE: case PS_FITS_SCALE_LOG_STDEV_POSITIVE: + case PS_FITS_SCALE_ASINH_STDEV_POSITIVE: // Put (mean - N sigma) at the lowest possible value: predominantly positive images imageVal = mean - options->stdevNum * stdev; @@ -220,4 +221,5 @@ case PS_FITS_SCALE_STDEV_NEGATIVE: case PS_FITS_SCALE_LOG_STDEV_NEGATIVE: + case PS_FITS_SCALE_ASINH_STDEV_NEGATIVE: // Put (mean + N sigma) at the highest possible value: predominantly negative images imageVal = mean + options->stdevNum * stdev; @@ -226,4 +228,5 @@ case PS_FITS_SCALE_STDEV_BOTH: case PS_FITS_SCALE_LOG_STDEV_BOTH: + case PS_FITS_SCALE_ASINH_STDEV_BOTH: // Put mean right in the middle: images with an equal abundance of positive and negative values imageVal = mean; @@ -241,4 +244,6 @@ } + + static bool logscaleStdev(double *bscale, // Scaling, to return @@ -456,4 +461,347 @@ } +static bool asinhStdev(double *bscale, // Scaling, to return + double *bzero, // Zero point, to return + double *boffset, // asinh flux zeropoint, to return + double *bsoften, // asinh softening parameter, to return + const psImage *image, // Image to scale + const psImage *mask, // Mask image + psImageMaskType maskVal, // value to mask + const psFitsOptions *options // FITS options + ) +{ + psAssert(bscale, "impossible"); + psAssert(bzero, "impossible"); + psAssert(boffset, "impossible"); + psAssert(bsoften, "impossible"); + psAssert(image, "impossible"); + psAssert(options, "impossible"); + + psTrace("psLib.fits", 3, "Scaling image by asinh method"); + int numCols = image->numCols, numRows = image->numRows; // Size of image + + // Measure the mean and stdev + // psImageBackground automatically excludes pixels that are non-finite, so we don't need to bother about a + // mask. + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); + psStats *stats = psStatsAlloc(MEAN_STAT | STDEV_STAT); // Statistics object + double mean, stdev; // Mean and standard deviation + if (!psImageBackground(stats, NULL, image, mask, maskVal, rng)) { + // It could be because the image is entirely masked, in which case we don't want to error + bool good = false; // Any good pixels? + + +// Find good pixels in an image, by image type +#define GOOD_PIXELS_CASE(TYPE) \ + case PS_TYPE_##TYPE: \ + for (int y = 0; y < image->numRows && !good; y++) { \ + for (int x = 0; x < image->numCols && !good; x++) { \ + if (mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal)) { \ + continue; \ + } \ + if (!isfinite(image->data.TYPE[y][x])) { \ + continue; \ + } \ + good = true; \ + } \ + } \ + break; + + switch (image->type.type) { + GOOD_PIXELS_CASE(F32); + GOOD_PIXELS_CASE(F64); + default: + psAbort("Unsupported case: %x", image->type.type); + } + + if (!good) { + psLogMsg("psLib.fits", PS_LOG_DETAIL, "Image has no good pixels, setting BSCALE = 1, BZERO = 0"); + psErrorClear(); + *bscale = 1.0; + *bzero = 0.0; + *bsoften = NAN; + psFree(rng); + psFree(stats); + return true; + } + + // There are some good pixels in there somewhere; psImageBackground just didn't find them + psLogMsg("psLib.fits", PS_LOG_DETAIL, + "Couldn't measure background statistics for image quantisation; retrying."); + psErrorClear(); + // Retry using all the available pixels + stats->nSubsample = image->numCols * image->numRows + 1; + if (!psImageStats(stats, image, mask, maskVal)) { + psLogMsg("psLib.fits", PS_LOG_DETAIL, + "Couldn't measure background statistics for image quantisation (attempt 2); retrying."); + psErrorClear(); + // Retry with desperate statistic + stats->options = DESPERATE_MEAN_STAT | DESPERATE_STDEV_STAT; + if (!psImageStats(stats, image, mask, maskVal)) { + psError(PS_ERR_UNKNOWN, false, "Unable to measure background statistics for image"); + psFree(rng); + psFree(stats); + return false; + } else { + // Desperate retry + mean = psStatsGetValue(stats, DESPERATE_MEAN_STAT); + stdev = psStatsGetValue(stats, DESPERATE_STDEV_STAT); + } + } else { + // Retry with all available pixels + mean = psStatsGetValue(stats, MEAN_STAT); + stdev = psStatsGetValue(stats, STDEV_STAT); + } + } else { + // First attempt + mean = psStatsGetValue(stats, MEAN_STAT); + stdev = psStatsGetValue(stats, STDEV_STAT); + } + psFree(rng); + psFree(stats); + if (!isfinite(mean) || !isfinite(stdev)) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "Mean (%f) or stdev (%f) is non-finite.", mean, stdev); + return false; + } + + psImage *copy; + + switch (image->type.type) { + case PS_TYPE_F32: + copy = psImageCopy(NULL,image,PS_TYPE_F32); + break; + case PS_TYPE_F64: + copy = psImageCopy(NULL,image,PS_TYPE_F64); + break; + default: + psError(PS_ERR_UNKNOWN, true, "Target type is not a float: %d", image->type.type); + return NULL; + break; + } + // Do scaling + float a = 1.0857362; // 2.5 * log(e); + *bsoften = sqrt(a) * stdev; + *boffset = mean; + // float m0 = 0; // Can I just arbitrarily set this? + + switch (image->type.type) { + case PS_TYPE_F32: + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { +/* if (x == 266 && y == 4584) { */ +/* 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), */ +/* a * asinh((image->data.F32[y][x] - *boffset) / (2 * *bsoften))); */ +/* } */ + if (isfinite(image->data.F32[y][x])) { + copy->data.F32[y][x] = a * asinh( (image->data.F32[y][x] - *boffset) / (2 * *bsoften));// - 2.5 * log10(b) + m0; + } + else { + copy->data.F32[y][x] = image->data.F32[y][x]; + } + } + } + break; + case PS_TYPE_F64: + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + // fprintf(stderr,"psFS64: %d %d %g %g %g\n",x,y,offset,image->data.F64[y][x],log10(image->data.F64[y][x] - offset)); + if (isfinite(image->data.F64[y][x])) { + copy->data.F64[y][x] = a * asinh( (image->data.F64[y][x] - *boffset) / (2 * *bsoften));// - 2.5 * log10(b) + m0; + } + else { + copy->data.F64[y][x] = image->data.F64[y][x]; + } + } + } + break; + default: + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Target type is not a float: %d",image->type.type); + return NULL; + break; + } + + // Do regular scaling on the logarithm image + if (!scaleStdev(bscale, bzero, copy, mask, maskVal, options)) { + psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev"); + return false; + } + psFree(copy); + return true; +} + +static bool asinhRange(double *bscale, // Scaling, to return + double *bzero, // Zero point, to return + double *boffset, // asinh flux zeropoint, to return + double *bsoften, // asinh softening parameter, to return + const psImage *image, // Image to scale + const psImage *mask, // Mask image + psImageMaskType maskVal, // value to mask + const psFitsOptions *options // FITS options + ) +{ + psAssert(bscale, "impossible"); + psAssert(bzero, "impossible"); + psAssert(bsoften, "impossible"); + psAssert(image, "impossible"); + psAssert(options, "impossible"); + + psTrace("psLib.fits", 3, "Scaling image by asinh method"); + int numCols = image->numCols, numRows = image->numRows; // Size of image + + // Measure the mean and stdev + // psImageBackground automatically excludes pixels that are non-finite, so we don't need to bother about a + // mask. + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); + psStats *stats = psStatsAlloc(MEAN_STAT | STDEV_STAT); // Statistics object + double mean, stdev; // Mean and standard deviation + if (!psImageBackground(stats, NULL, image, mask, maskVal, rng)) { + // It could be because the image is entirely masked, in which case we don't want to error + bool good = false; // Any good pixels? + + +// Find good pixels in an image, by image type +#define GOOD_PIXELS_CASE(TYPE) \ + case PS_TYPE_##TYPE: \ + for (int y = 0; y < image->numRows && !good; y++) { \ + for (int x = 0; x < image->numCols && !good; x++) { \ + if (mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal)) { \ + continue; \ + } \ + if (!isfinite(image->data.TYPE[y][x])) { \ + continue; \ + } \ + good = true; \ + } \ + } \ + break; + + switch (image->type.type) { + GOOD_PIXELS_CASE(F32); + GOOD_PIXELS_CASE(F64); + default: + psAbort("Unsupported case: %x", image->type.type); + } + + if (!good) { + psLogMsg("psLib.fits", PS_LOG_DETAIL, "Image has no good pixels, setting BSCALE = 1, BZERO = 0"); + psErrorClear(); + *bscale = 1.0; + *bzero = 0.0; + *bsoften = NAN; + psFree(rng); + psFree(stats); + return true; + } + + // There are some good pixels in there somewhere; psImageBackground just didn't find them + psLogMsg("psLib.fits", PS_LOG_DETAIL, + "Couldn't measure background statistics for image quantisation; retrying."); + psErrorClear(); + // Retry using all the available pixels + stats->nSubsample = image->numCols * image->numRows + 1; + if (!psImageStats(stats, image, mask, maskVal)) { + psLogMsg("psLib.fits", PS_LOG_DETAIL, + "Couldn't measure background statistics for image quantisation (attempt 2); retrying."); + psErrorClear(); + // Retry with desperate statistic + stats->options = DESPERATE_MEAN_STAT | DESPERATE_STDEV_STAT; + if (!psImageStats(stats, image, mask, maskVal)) { + psError(PS_ERR_UNKNOWN, false, "Unable to measure background statistics for image"); + psFree(rng); + psFree(stats); + return false; + } else { + // Desperate retry + mean = psStatsGetValue(stats, DESPERATE_MEAN_STAT); + stdev = psStatsGetValue(stats, DESPERATE_STDEV_STAT); + } + } else { + // Retry with all available pixels + mean = psStatsGetValue(stats, MEAN_STAT); + stdev = psStatsGetValue(stats, STDEV_STAT); + } + } else { + // First attempt + mean = psStatsGetValue(stats, MEAN_STAT); + stdev = psStatsGetValue(stats, STDEV_STAT); + } + psFree(rng); + psFree(stats); + if (!isfinite(mean) || !isfinite(stdev)) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "Mean (%f) or stdev (%f) is non-finite.", mean, stdev); + return false; + } + + psImage *copy; + + switch (image->type.type) { + case PS_TYPE_F32: + copy = psImageCopy(NULL,image,PS_TYPE_F32); + break; + case PS_TYPE_F64: + copy = psImageCopy(NULL,image,PS_TYPE_F64); + break; + default: + psError(PS_ERR_UNKNOWN, true, "Target type is not a float: %d", image->type.type); + return NULL; + break; + } + // Do scaling + float a = 1.0857362; // 2.5 * log(e); + *bsoften = sqrt(a) * stdev; + *boffset = mean; + // float m0 = 0; // Can I just arbitrarily set this? + + switch (image->type.type) { + case PS_TYPE_F32: + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { +/* if (x == 2331 && y == 2843) { */ +/* fprintf(stderr,"psFS32: %d %d %g %g %g\n",x,y,offset,image->data.F32[y][x],log10(image->data.F32[y][x] - offset)); */ +/* } */ +/* if (x == 266 && y == 4584) { */ +/* 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), */ +/* a * asinh((image->data.F32[y][x] - *boffset) / (2 * *bsoften))); */ +/* } */ + + if (isfinite(image->data.F32[y][x])) { + copy->data.F32[y][x] = a * asinh( (image->data.F32[y][x] - *boffset) / (2 * *bsoften));// - 2.5 * log10(b) + m0; + } + else { + copy->data.F32[y][x] = image->data.F32[y][x]; + } + } + } + break; + case PS_TYPE_F64: + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + // fprintf(stderr,"psFS64: %d %d %g %g %g\n",x,y,offset,image->data.F64[y][x],log10(image->data.F64[y][x] - offset)); + if (isfinite(image->data.F64[y][x])) { + copy->data.F64[y][x] = a * asinh( (image->data.F64[y][x] - *boffset)/ (2 * *bsoften));// - 2.5 * log10(b) + m0; + } + else { + copy->data.F64[y][x] = image->data.F32[y][x]; + } + } + } + break; + default: + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Target type is not a float: %d",image->type.type); + return NULL; + break; + } + + // Do regular scaling on the asinh image + if (!scaleRange(bscale, bzero, copy, options)) { + psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev"); + return false; + } + psFree(copy); + return true; +} + ////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -461,5 +809,6 @@ ////////////////////////////////////////////////////////////////////////////////////////////////////////////// -bool psFitsScaleDetermine(double *bscale, double *bzero, double *boffset, long *blank, const psImage *image, +bool psFitsScaleDetermine(double *bscale, double *bzero, double *boffset, double *bsoften, + long *blank, const psImage *image, const psImage *mask, psImageMaskType maskVal, const psFits *fits) { @@ -521,18 +870,34 @@ } break; - case PS_FITS_SCALE_LOG_RANGE: - if (!logscaleRange(bscale,bzero,boffset,image,options)) { - psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from range"); - return false; - } - break; - case PS_FITS_SCALE_LOG_STDEV_POSITIVE: - case PS_FITS_SCALE_LOG_STDEV_NEGATIVE: - case PS_FITS_SCALE_LOG_STDEV_BOTH: - if (!logscaleStdev(bscale, bzero,boffset, image, mask, maskVal, options)) { - psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev"); - return false; - } - break; + case PS_FITS_SCALE_LOG_RANGE: + if (!logscaleRange(bscale,bzero,boffset,image,options)) { + psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from range"); + return false; + } + break; + case PS_FITS_SCALE_LOG_STDEV_POSITIVE: + case PS_FITS_SCALE_LOG_STDEV_NEGATIVE: + case PS_FITS_SCALE_LOG_STDEV_BOTH: + if (!logscaleStdev(bscale, bzero,boffset, image, mask, maskVal, options)) { + psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev"); + return false; + } + break; + case PS_FITS_SCALE_ASINH_RANGE: + if (!asinhRange(bscale,bzero,boffset,bsoften,image,mask,maskVal,options)) { + psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from range"); + return false; + } + break; + case PS_FITS_SCALE_ASINH_STDEV_POSITIVE: + case PS_FITS_SCALE_ASINH_STDEV_NEGATIVE: + case PS_FITS_SCALE_ASINH_STDEV_BOTH: + if (!asinhStdev(bscale, bzero,boffset,bsoften, image, mask, maskVal, options)) { + psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev"); + return false; + } + break; + + case PS_FITS_SCALE_MANUAL: *bscale = options->bscale; @@ -557,11 +922,11 @@ - psTrace("psLib.fits", 3, "BSCALE = %.10lf, BZERO = %.10lf, BOFFSET = %.10lf, BLANK = %ld\n", - *bscale, *bzero, *boffset, *blank); + psTrace("psLib.fits", 3, "BSCALE = %.10lf, BZERO = %.10lf, BOFFSET = %.10lf, BSOFTEN = %.10lf, BLANK = %ld\n", + *bscale, *bzero, *boffset, *bsoften, *blank); return true; } -psImage *psFitsScaleForDisk(const psImage *image, const psFits *fits, double bscale, double bzero, double boffset, +psImage *psFitsScaleForDisk(const psImage *image, const psFits *fits, double bscale, double bzero, double boffset, double bsoften, psRandom *rng) { @@ -625,5 +990,26 @@ (options->scaling == PS_FITS_SCALE_LOG_STDEV_NEGATIVE)|| \ (options->scaling == PS_FITS_SCALE_LOG_STDEV_BOTH)) { \ - value = log10( (IN)->data.INTYPE[y][x] - boffset ); \ + if (isfinite((IN)->data.INTYPE[y][x])) { \ + value = log10( (IN)->data.INTYPE[y][x] - boffset ); \ + } \ + else { \ + value = (IN)->data.INTYPE[y][x]; \ + } \ + } \ + else if ((options->scaling == PS_FITS_SCALE_ASINH_RANGE)|| \ + (options->scaling == PS_FITS_SCALE_ASINH_MANUAL)|| \ + (options->scaling == PS_FITS_SCALE_ASINH_STDEV_POSITIVE)|| \ + (options->scaling == PS_FITS_SCALE_ASINH_STDEV_NEGATIVE)|| \ + (options->scaling == PS_FITS_SCALE_ASINH_STDEV_BOTH)) { \ + if (isfinite((IN)->data.INTYPE[y][x])) { \ + value = 1.0857362 * (asinh( ((IN)->data.INTYPE[y][x] - boffset) / (2.0 * bsoften))); \ + } \ + else { \ + value = (IN)->data.INTYPE[y][x]; \ + } \ + if ((x == 1000)&&(y == 1000)) { \ + fprintf(stderr,"ASINH: %f %f %f", \ + (IN)->data.INTYPE[y][x],value,bsoften); \ + } \ } \ else { \ @@ -652,8 +1038,8 @@ case PS_TYPE_##INTYPE: { \ switch (outType) { \ - SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, U8); \ - SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S16); \ - SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S32); \ - SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S64); \ + SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, U8);; \ + SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S16);; \ + SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S32);; \ + SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S64);; \ default: \ psAbort("Should be unreachable."); \ @@ -664,5 +1050,5 @@ switch (image->type.type) { SCALE_WRITE_IN_CASE(image, F32, out); - SCALE_WRITE_IN_CASE(image, F64, out); + SCALE_WRITE_IN_CASE(image, F64, out); default: psAbort("Should be unreachable."); @@ -679,5 +1065,5 @@ // the present time, since cfitsio should apply the scaling itself in the process of reading. However, we may // later desire it (e.g., if we ever make our own FITS implementation). -psImage *psFitsScaleFromDisk(const psImage *image, double boffset) +psImage *psFitsScaleFromDisk(const psImage *image, double boffset, double bsoften) { PS_ASSERT_IMAGE_NON_NULL(image, NULL); @@ -718,5 +1104,10 @@ for (int y = 0; y < numRows; y++) { \ for (int x = 0; x < numCols; x++) { \ - out->data.OUTTYPE[y][x] = pow(10,image->data.INTYPE[y][x]) + boffset;; \ + if (bsoften) { \ + out->data.OUTTYPE[y][x] = 2 * bsoften * sinh(image->data.INTYPE[y][x] / (1.0857362)) + boffset; \ + } \ + else if (boffset) { \ + out->data.OUTTYPE[y][x] = pow(10,image->data.INTYPE[y][x]) + boffset;; \ + } \ } \ } \ @@ -766,4 +1157,9 @@ if (strcasecmp(string, "LOG_STDEV_NEGATIVE") == 0) return PS_FITS_SCALE_LOG_STDEV_NEGATIVE; if (strcasecmp(string, "LOG_STDEV_BOTH") == 0) return PS_FITS_SCALE_LOG_STDEV_BOTH; + if (strcasecmp(string, "ASINH_RANGE") == 0) return PS_FITS_SCALE_ASINH_RANGE; + if (strcasecmp(string, "ASINH_MANUAL") == 0) return PS_FITS_SCALE_ASINH_MANUAL; + if (strcasecmp(string, "ASINH_STDEV_POSITIVE") == 0) return PS_FITS_SCALE_ASINH_STDEV_POSITIVE; + if (strcasecmp(string, "ASINH_STDEV_NEGATIVE") == 0) return PS_FITS_SCALE_ASINH_STDEV_NEGATIVE; + if (strcasecmp(string, "ASINH_STDEV_BOTH") == 0) return PS_FITS_SCALE_ASINH_STDEV_BOTH; if (strcasecmp(string, "MANUAL") == 0) return PS_FITS_SCALE_MANUAL; Index: /trunk/psLib/src/fits/psFitsScale.h =================================================================== --- /trunk/psLib/src/fits/psFitsScale.h (revision 31632) +++ /trunk/psLib/src/fits/psFitsScale.h (revision 31633) @@ -11,4 +11,5 @@ double *bzero, ///< Zero point, to return double *boffset, ///< Log offset, to return + double *bsoften, ///< asinh softening parameter, to return long *blank, ///< Blank value, to return const psImage *image, ///< Image to scale @@ -29,8 +30,10 @@ double bzero, ///< Zero point double boffset, ///< Log offset + double bsoften, ///< asinh softening parameter psRandom *rng ///< Random number generator (for the "fuzz"), or NULL ); psImage *psFitsScaleFromDisk(const psImage *image, ///< Image to to unapply BOFFSET - double boffset ///< Log offset + double boffset, ///< Log offset + double bsoften ///< asinh softening parameter ); /// Interpret a string as a scaling method Index: /trunk/psModules/src/camera/pmFPAfileDefine.c =================================================================== --- /trunk/psModules/src/camera/pmFPAfileDefine.c (revision 31632) +++ /trunk/psModules/src/camera/pmFPAfileDefine.c (revision 31633) @@ -283,4 +283,5 @@ case PS_FITS_SCALE_RANGE: case PS_FITS_SCALE_LOG_RANGE: + case PS_FITS_SCALE_ASINH_RANGE: // No options required break; @@ -289,4 +290,6 @@ case PS_FITS_SCALE_LOG_STDEV_POSITIVE: case PS_FITS_SCALE_LOG_STDEV_NEGATIVE: + case PS_FITS_SCALE_ASINH_STDEV_POSITIVE: + case PS_FITS_SCALE_ASINH_STDEV_NEGATIVE: options->stdevNum = parseOptionFloat(scheme, "STDEV.NUM", source); // Padding to edge if (!isfinite(options->stdevNum)) { @@ -299,4 +302,5 @@ case PS_FITS_SCALE_STDEV_BOTH: case PS_FITS_SCALE_LOG_STDEV_BOTH: + case PS_FITS_SCALE_ASINH_STDEV_BOTH: options->stdevBits = parseOptionInt(scheme, "STDEV.BITS", source, 0); // Bits for stdev if (options->stdevBits <= 0) { @@ -312,10 +316,15 @@ options->bzero = parseOptionDouble(scheme, "BZERO", source); // Zero point break; - case PS_FITS_SCALE_LOG_MANUAL: - options->bscale = parseOptionDouble(scheme, "BSCALE", source); // Scaling - options->bzero = parseOptionDouble(scheme, "BZERO", source); // Zero point - options->boffset = parseOptionDouble(scheme, "BOFFSET", source); // Log offset - break; - default: + case PS_FITS_SCALE_LOG_MANUAL: + options->bscale = parseOptionDouble(scheme, "BSCALE", source); // Scaling + options->bzero = parseOptionDouble(scheme, "BZERO", source); // Zero point + options->boffset = parseOptionDouble(scheme, "BOFFSET", source); // Log offset + case PS_FITS_SCALE_ASINH_MANUAL: + options->bscale = parseOptionDouble(scheme, "BSCALE", source); // Scaling + options->bzero = parseOptionDouble(scheme, "BZERO", source); // Zero point + options->boffset = parseOptionDouble(scheme, "BOFFSET", source); // Log offset + options->bsoften = parseOptionDouble(scheme, "BSOFTEN", source); // Softening parameter + break; + default: psAbort("Should never get here."); }