Index: trunk/psLib/src/fits/psFitsScale.c =================================================================== --- trunk/psLib/src/fits/psFitsScale.c (revision 27417) +++ trunk/psLib/src/fits/psFitsScale.c (revision 30636) @@ -100,4 +100,5 @@ } + // Determine appropriate BSCALE and BZERO for an image, mapping the standard deviation to the nominated number // of bits @@ -212,4 +213,5 @@ switch (options->scaling) { case PS_FITS_SCALE_STDEV_POSITIVE: + case PS_FITS_SCALE_LOG_STDEV_POSITIVE: // Put (mean - N sigma) at the lowest possible value: predominantly positive images imageVal = mean - options->stdevNum * stdev; @@ -217,4 +219,5 @@ break; case PS_FITS_SCALE_STDEV_NEGATIVE: + case PS_FITS_SCALE_LOG_STDEV_NEGATIVE: // Put (mean + N sigma) at the highest possible value: predominantly negative images imageVal = mean + options->stdevNum * stdev; @@ -222,4 +225,5 @@ break; case PS_FITS_SCALE_STDEV_BOTH: + case PS_FITS_SCALE_LOG_STDEV_BOTH: // Put mean right in the middle: images with an equal abundance of positive and negative values imageVal = mean; @@ -237,13 +241,230 @@ } + +static bool logscaleStdev(double *bscale, // Scaling, to return + double *bzero, // Zero point, to return + double *boffset, // Log offset, 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(image, "impossible"); + psAssert(options, "impossible"); + + psTrace("psLib.fits", 3, "Scaling image by logarithm statistics"); + int numCols = image->numCols, numRows = image->numRows; // Size of image + + psImage *copy; + + *boffset = 99e99; + + // Make a copy of the image to pass to get the scaling parameters + + 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; + } + + // Determine the minimum value on this image. + switch (image->type.type) { + case PS_TYPE_F32: + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + psF32 value = image->data.F32[y][x]; + if (isfinite(value)) { + if (value < *boffset) { + *boffset = value; + } + } + } + } + break; + case PS_TYPE_F64: + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + psF64 value = image->data.F64[y][x]; + if (isfinite(value)) { + if (value < *boffset) { + *boffset = value; + } + } + } + } + break; + default: + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Target type is not a float: %d",image->type.type); + return NULL; + break; + } + // We only need to offset images that go negative. + if (*boffset > 0.0) { + *boffset = 0.0; + } + // Write offset to header + // How? + // psMetadataAddF32(header,PS_LIST_TAIL,"LOGZERO",0,"Flux offset subtracted before taking logarithm.",offset); + // Take the logarithm of the image, applying the offset + 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] = (log10( image->data.F32[y][x] - *boffset)); + } + } + 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.F64[y][x] = (log10( image->data.F64[y][x] - *boffset)); + } + } + 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 logscaleRange(double *bscale, // Scaling, to return + double *bzero, // Zero point, to return + double *boffset, // Log offset, to return + const psImage *image, // Image to scale + const psFitsOptions *options // FITS options + ) +{ + psAssert(bscale, "impossible"); + psAssert(bzero, "impossible"); + psAssert(boffset, "impossible"); + psAssert(image, "impossible"); + psAssert(options, "impossible"); + + psTrace("psLib.fits", 3, "Scaling image by logarithm statistics"); + int numCols = image->numCols, numRows = image->numRows; // Size of image + + psImage *copy; + + *boffset = 99e99; + + // Make a copy of the image to pass to get the scaling parameters + + 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; + } + + // Determine the minimum value on this image. + switch (image->type.type) { + case PS_TYPE_F32: + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + psF32 value = image->data.F32[y][x]; + if (!isfinite(value)) { + if (value < *boffset) { + *boffset = value; + } + } + } + } + break; + case PS_TYPE_F64: + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + psF64 value = image->data.F64[y][x]; + if (!isfinite(value)) { + if (value < *boffset) { + *boffset = value; + } + } + } + } + break; + default: + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Target type is not a float: %d",image->type.type); + return NULL; + break; + } + // We only need to offset images that go negative. + if (*boffset > 0.0) { + *boffset = 0.0; + } + // Write offset to header + // How? + // psMetadataAddF32(header,PS_LIST_TAIL,"LOGZERO",0,"Flux offset subtracted before taking logarithm.",offset); + // Take the logarithm of the image, applying the offset + switch (image->type.type) { + case PS_TYPE_F32: + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + copy->data.F32[y][x] = (log10( image->data.F32[y][x] - *boffset)); + } + } + break; + case PS_TYPE_F64: + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + copy->data.F64[y][x] = (log10( image->data.F64[y][x] - *boffset)); + } + } + 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; +} + + ////////////////////////////////////////////////////////////////////////////////////////////////////////////// // Public functions ////////////////////////////////////////////////////////////////////////////////////////////////////////////// -bool psFitsScaleDetermine(double *bscale, double *bzero, long *blank, const psImage *image, +bool psFitsScaleDetermine(double *bscale, double *bzero, double *boffset, long *blank, const psImage *image, const psImage *mask, psImageMaskType maskVal, const psFits *fits) { PS_ASSERT_PTR_NON_NULL(bscale, false); PS_ASSERT_PTR_NON_NULL(bzero, false); + PS_ASSERT_PTR_NON_NULL(boffset, false); PS_ASSERT_PTR_NON_NULL(blank, false); PS_ASSERT_IMAGE_NON_NULL(image, false); @@ -256,4 +477,5 @@ *bscale = NAN; *bzero = NAN; + *boffset = 0; *blank = 0; @@ -299,7 +521,26 @@ } 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_MANUAL: *bscale = options->bscale; *bzero = options->bzero; + break; + case PS_FITS_SCALE_LOG_MANUAL: + *bscale = options->bscale; + *bzero = options->bzero; + *boffset = options->boffset; break; default: @@ -316,10 +557,11 @@ - psTrace("psLib.fits", 3, "BSCALE = %.10lf, BZERO = %.10lf, BLANK = %ld\n", *bscale, *bzero, *blank); + psTrace("psLib.fits", 3, "BSCALE = %.10lf, BZERO = %.10lf, BOFFSET = %.10lf, BLANK = %ld\n", + *bscale, *bzero, *boffset, *blank); return true; } -psImage *psFitsScaleForDisk(const psImage *image, const psFits *fits, double bscale, double bzero, +psImage *psFitsScaleForDisk(const psImage *image, const psFits *fits, double bscale, double bzero, double boffset, psRandom *rng) { @@ -377,5 +619,15 @@ for (int y = 0; y < numRows; y++) { \ for (int x = 0; x < numCols; x++) { \ - ps##INTYPE value = (IN)->data.INTYPE[y][x]; \ + ps##INTYPE value; \ + if ((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)) { \ + value = log10( (IN)->data.INTYPE[y][x] - boffset ); \ + } \ + else { \ + value = (IN)->data.INTYPE[y][x]; \ + } \ if (!isfinite(value)) { \ /* This choice of "max" for non-finite pixels is mainly cosmetic --- it has to be */ \ @@ -423,16 +675,11 @@ -#if 0 + // This function to apply BSCALE and BZERO to an image read immediately from disk should not be necessary at // 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(psFits *fits, psImage *image) +psImage *psFitsScaleFromDisk(const psImage *image, double boffset) { PS_ASSERT_IMAGE_NON_NULL(image, NULL); - - if (bscale == 0.0) { - // BSCALE = 0 means don't apply anything - return psMemIncrRefCounter(image); - } psElemType inType = image->type.type; // Type for input image @@ -471,5 +718,5 @@ for (int y = 0; y < numRows; y++) { \ for (int x = 0; x < numCols; x++) { \ - out->data.OUTTYPE[y][x] = image->data.INTYPE[y][x] * bscale + bzero; \ + out->data.OUTTYPE[y][x] = pow(10,image->data.INTYPE[y][x]) + boffset;; \ } \ } \ @@ -505,7 +752,4 @@ return out; } -#endif - - psFitsScaling psFitsScalingFromString(const char *string) @@ -517,4 +761,9 @@ if (strcasecmp(string, "STDEV_NEGATIVE") == 0) return PS_FITS_SCALE_STDEV_NEGATIVE; if (strcasecmp(string, "STDEV_BOTH") == 0) return PS_FITS_SCALE_STDEV_BOTH; + if (strcasecmp(string, "LOG_RANGE") == 0) return PS_FITS_SCALE_LOG_RANGE; + if (strcasecmp(string, "LOG_MANUAL") == 0) return PS_FITS_SCALE_LOG_MANUAL; + if (strcasecmp(string, "LOG_STDEV_POSITIVE") == 0) return PS_FITS_SCALE_LOG_STDEV_POSITIVE; + 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, "MANUAL") == 0) return PS_FITS_SCALE_MANUAL;