Index: /branches/czw_branch/20110406/ippconfig/recipes/fitstypes.mdc =================================================================== --- /branches/czw_branch/20110406/ippconfig/recipes/fitstypes.mdc (revision 31617) +++ /branches/czw_branch/20110406/ippconfig/recipes/fitstypes.mdc (revision 31618) @@ -101,5 +101,5 @@ COMP_STACK METADATA BITPIX S32 16 - SCALING STR LOG_STDEV_POSITIVE + SCALING STR ASINH_STDEV_POSITIVE STDEV.BITS S32 4 STDEV.NUM F32 10 Index: /branches/czw_branch/20110406/psLib/src/fits/psFits.h =================================================================== --- /branches/czw_branch/20110406/psLib/src/fits/psFits.h (revision 31617) +++ /branches/czw_branch/20110406/psLib/src/fits/psFits.h (revision 31618) @@ -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: /branches/czw_branch/20110406/psLib/src/fits/psFitsImage.c =================================================================== --- /branches/czw_branch/20110406/psLib/src/fits/psFitsImage.c (revision 31617) +++ /branches/czw_branch/20110406/psLib/src/fits/psFitsImage.c (revision 31618) @@ -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,19 @@ 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; + } + else { + info->is_asinh = false; + } + status = 0; @@ -262,4 +278,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 +292,5 @@ psAssert(bzero, "impossible"); psAssert(boffset, "impossible"); + psAssert(bsoften, "impossible"); psAssert(floatType, "impossible"); psAssert(fits, "impossible"); @@ -322,5 +340,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 +366,5 @@ } - return psFitsScaleForDisk(image, fits, *bscale, *bzero, *boffset, rng); + return psFitsScaleForDisk(image, fits, *bscale, *bzero, *boffset, *bsoften, rng); } @@ -487,5 +505,14 @@ 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; + int status; + status = 0; + fits_read_key_dbl(fits->fd, "BSOFTEN", &bsoften, NULL, &status); + psImage *newImage = psFitsScaleFromDisk(outImage,0.0,bsoften); psFree (outImage); outImage = newImage; @@ -602,7 +629,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) { @@ -831,7 +859,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: /branches/czw_branch/20110406/psLib/src/fits/psFitsScale.c =================================================================== --- /branches/czw_branch/20110406/psLib/src/fits/psFitsScale.c (revision 31617) +++ /branches/czw_branch/20110406/psLib/src/fits/psFitsScale.c (revision 31618) @@ -242,4 +242,6 @@ + + static bool logscaleStdev(double *bscale, // Scaling, to return double *bzero, // Zero point, to return @@ -456,4 +458,316 @@ } +static bool asinhStdev(double *bscale, // Scaling, to return + double *bzero, // Zero point, 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; + // 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)); */ +/* } */ + copy->data.F32[y][x] = -1.0 * a * asinh( image->data.F32[y][x] / (2 * *bsoften));// - 2.5 * log10(b) + m0; + } + } + 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)); + copy->data.F32[y][x] = -1.0 * a * asinh( image->data.F32[y][x] / (2 * *bsoften));// - 2.5 * log10(b) + m0; + } + } + 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 *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; + // 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)); */ +/* } */ + copy->data.F32[y][x] = -1.0 * a * asinh( image->data.F32[y][x] / (2 * *bsoften));// - 2.5 * log10(b) + m0; + } + } + 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)); + copy->data.F32[y][x] = -1.0 * a * asinh( image->data.F32[y][x] / (2 * *bsoften));// - 2.5 * log10(b) + m0; + } + } + 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 (!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 +775,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 +836,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,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,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; @@ -563,5 +894,5 @@ -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) { @@ -626,4 +957,11 @@ (options->scaling == PS_FITS_SCALE_LOG_STDEV_BOTH)) { \ value = log10( (IN)->data.INTYPE[y][x] - boffset ); \ + } \ + 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)) { \ + value = -1.0 * 1.0857362 * (asinh( (IN)->data.INTYPE[y][x] - bsoften)); \ } \ else { \ @@ -652,8 +990,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 +1002,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 +1017,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 +1056,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 (boffset) { \ + out->data.OUTTYPE[y][x] = pow(10,image->data.INTYPE[y][x]) + boffset;; \ + } \ + else if (bsoften) { \ + out->data.OUTTYPE[y][x] = sinh(image->data.INTYPE[y][x] * 2 * bsoften) / (-1.0 * 1.0857362);; \ + } \ } \ } \ @@ -766,4 +1109,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: /branches/czw_branch/20110406/psLib/src/fits/psFitsScale.h =================================================================== --- /branches/czw_branch/20110406/psLib/src/fits/psFitsScale.h (revision 31617) +++ /branches/czw_branch/20110406/psLib/src/fits/psFitsScale.h (revision 31618) @@ -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: /branches/czw_branch/20110406/psModules/src/camera/pmFPAfileDefine.c =================================================================== --- /branches/czw_branch/20110406/psModules/src/camera/pmFPAfileDefine.c (revision 31617) +++ /branches/czw_branch/20110406/psModules/src/camera/pmFPAfileDefine.c (revision 31618) @@ -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,14 @@ 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->bsoften = parseOptionDouble(scheme, "BSOFTEN", source); // Softening parameter + break; + default: psAbort("Should never get here."); } Index: /branches/czw_branch/20110406/psModules/src/imcombine/pmStack.c =================================================================== --- /branches/czw_branch/20110406/psModules/src/imcombine/pmStack.c (revision 31617) +++ /branches/czw_branch/20110406/psModules/src/imcombine/pmStack.c (revision 31618) @@ -22,4 +22,6 @@ #include +#include + #include "pmHDU.h" #include "pmFPA.h" @@ -35,8 +37,8 @@ -# if (0) +# if (1) #define TESTING // Enable test output -#define TEST_X 3145 // x coordinate to examine -#define TEST_Y 2334 // y coordinate to examine +#define TEST_X 5745 // x coordinate to examine +#define TEST_Y 5331 // y coordinate to examine #define TEST_RADIUS 0.5 // Radius to examine # endif @@ -332,5 +334,6 @@ memset (nGoodBits, 0, 16*sizeof(int)); - // Extract the pixel and mask data + + // Extract the pixel and mask data int numGood = 0; // Number of good pixels for (int i = 0, j = 0; i < inputs->n; i++) { @@ -364,5 +367,5 @@ continue; } - + int xIn = x - data->readout->col0, yIn = y - data->readout->row0; // Coordinates on input readout psImage *mask = data->readout->mask; // Mask of interest @@ -984,5 +987,5 @@ static void KMMcalculate(const psVector *values, - float *Punimodal, + float *Punimodal,int *iter, float *pi1, float *m1, float *s1, float *pi2, float *m2, float *s2) { @@ -1012,5 +1015,5 @@ float dL = 0; float oldL = -999; - int k = 0; + *iter = 0; logL_bimodal = logL_unimodal; *m1 = mU - 3 * sU; @@ -1024,6 +1027,9 @@ float w1,w2; - while (((dL > KMM_TOLERANCE)||(k < 3))&&(k < KMM_MAX_ITERATIONS)) { - k++; + float KMM_TOLERANCE = 1e-6; + int KMM_MAX_ITERATIONS = 500; + float KMM_SMALL_NUMBER = 1e-5; + while (((dL > KMM_TOLERANCE)||(*iter < 3))&&(*iter < KMM_MAX_ITERATIONS)) { + *iter += 1; dL = fabs(logL_bimodal - oldL); oldL = logL_bimodal; @@ -1055,6 +1061,6 @@ *m2 += values->data.F32[i] * P2->data.F32[i]; - w1 += P1[i]; - w2 += P2[i]; + w1 += P1->data.F32[i]; + w2 += P2->data.F32[i]; } *m1 /= w1; @@ -1095,12 +1101,29 @@ } -static void KMMrejectUnpopular(const psVector *values, psArray *reject) { +static void KMMRejectUnpopular(const psArray *inputs, int x, int y) { + float KMM_MINIMUM_PVALUE = 0.05; float Punimodal,pi1,m1,s1,pi2,m2,s2; - KMMcalculate(values,&Punimodal, + int iter; + int j; + + psVector *values = psVectorAlloc(inputs->n, PS_TYPE_F32); + for (j = 0; j < inputs->n; j++) { + pmStackData *data = inputs->data[j]; // Stack data of interest + psImage *image = data->readout->image; // Image of interest + int xIn = x - data->readout->col0, yIn = y - data->readout->row0; // Coordinates on input readout + values->data.F32[j] = image->data.F32[yIn][xIn]; + } + + KMMcalculate(values,&Punimodal,&iter, &pi1,&m1,&s1, &pi2,&m2,&s2); + // fprintf(stderr, + psTrace("psModules.imcombine",3, + "KMM Unpopular Test: %d,%d: Puni: %f in %d",x,y,Punimodal,iter); + // CHECKPIX(x, y, " + if (Punimodal < KMM_MINIMUM_PVALUE) { int i; - float g1,g2; + float g1,g2,norm; float P1,P2; @@ -1112,22 +1135,42 @@ P2 = (pi2 * g2) / norm; + CHECKPIX(x, y, "KMM Unpopular Rejection: %d,%d: %d %f %f:(%f %f %f ) %f:(%f %f %f) rejection? %d %d\n", + x, y, i, values->data.F32[i], + P1,m1,s1,pi1, + P2,m2,s2,pi2, + (pi1 > pi2)&&(P1 < P2), + (pi1 < pi2)&&(P1 > P2)); if ((pi1 > pi2)&&(P1 < P2)) { // mode 1 is more popular, but this element belongs to mode 2 - reject_input(reject,i); + combineMarkReject(inputs,x,y,i); } if ((pi1 < pi2)&&(P1 > P2)) { // mode 2 is more popular, but this element belongs to mode 1 - reject_input(reject,i); + combineMarkReject(inputs,x,y,i); } } } + psFree(values); // else do nothing. } -static void KMMrejectBright(const psVector *values, psArray *reject) { - KMMcalculate(values,&Punimodal, +static void KMMRejectBright(const psArray *inputs, int x, int y) { + float KMM_MINIMUM_PVALUE = 0.05; + float Punimodal,pi1,m1,s1,pi2,m2,s2; + int iter; + int j; + + psVector *values = psVectorAlloc(inputs->n, PS_TYPE_F32); + for (j = 0; j < inputs->n; j++) { + pmStackData *data = inputs->data[j]; // Stack data of interest + psImage *image = data->readout->image; // Image of interest + int xIn = x - data->readout->col0, yIn = y - data->readout->row0; // Coordinates on input readout + values->data.F32[j] = image->data.F32[yIn][xIn]; + } + + KMMcalculate(values,&Punimodal,&iter, &pi1,&m1,&s1, &pi2,&m2,&s2); if (Punimodal < KMM_MINIMUM_PVALUE) { int i; - float g1,g2; + float g1,g2,norm; float P1,P2; @@ -1140,8 +1183,8 @@ if ((m1 > m2)&&(P1 > P2)) { // m1 is larger, and this element belongs to mode 1 - reject_input(reject,i); + combineMarkReject(inputs,x,y,i); } if ((m1 < m2)&&(P1 < P2)) { // m2 is larger, and this element belongs to mode 2 - reject_input(reject,i); + combineMarkReject(inputs,x,y,i); } } @@ -1310,4 +1353,10 @@ } + // Pre-reject inputs using KMM bimodality test. + if (1) { +/* KMMRejectUnpopular(input,x,y); */ + rejection = true; + } + // Set up rejection list psArray *pixelMap = NULL; // Map of pixels to source @@ -1319,5 +1368,15 @@ for (int y = minInputRows; y < maxInputRows; y++) { for (int x = minInputCols; x < maxInputCols; x++) { + CHECKPIX(x, y, "Combining pixel %d,%d: %x %x %f %f %f %f %d %d %d\n", x, y, badMaskBits, blankMaskBits, iter, rej, sys, olympic, useVariance, safe, rejection); + + // Pre-reject inputs using KMM bimodality test. + if (1) { + KMMRejectUnpopular(input,x,y); +/* rejection = true; */ + } + else { + KMMRejectBright(input,x,y); + } #ifdef TESTING