Index: /branches/czw_branch/20101203/psLib/src/fits/psFitsImage.c =================================================================== --- /branches/czw_branch/20101203/psLib/src/fits/psFitsImage.c (revision 30330) +++ /branches/czw_branch/20101203/psLib/src/fits/psFitsImage.c (revision 30331) @@ -54,4 +54,5 @@ int fitsDatatype; // cfitsio data type int psDatatype; // psLib data type + bool is_logscaled; // is this image log scaled using BOFFSET? } p_psFitsReadInfo; @@ -128,5 +129,5 @@ // Check scale and zero - double bscale = 0.0, bzero = 0.0; // Scale and zero point + double bscale = 0.0, bzero = 0.0, boffset = 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."); @@ -137,4 +138,18 @@ psFitsError(status, true, "Unable to read header."); goto bad; + } + status = 0; + if (fits_read_key_dbl(fits->fd, "BOFFSET", &boffset, NULL, &status) && status != KEY_NO_EXIST) { + psFitsError(status, true, "Unable to read header."); + goto bad; + } + if (status == KEY_NO_EXIST) { + info->is_logscaled = false; + } + else if (isfinite(boffset)) { + info->is_logscaled = true; + } + else { + info->is_logscaled = false; } status = 0; @@ -246,4 +261,5 @@ static psImage *imageToDiskRepresentation(double *bscale, // Scaling applied double *bzero, // Zero point applied + double *boffset, // Log offset applied long *blank, // Blank value (integer data) psFitsFloat *floatType, // Type of custom floating-point @@ -258,4 +274,5 @@ psAssert(bscale, "impossible"); psAssert(bzero, "impossible"); + psAssert(boffset, "impossible"); psAssert(floatType, "impossible"); psAssert(fits, "impossible"); @@ -305,5 +322,5 @@ if (newScaleZero) { // Choose an appropriate BSCALE and BZERO - if (!psFitsScaleDetermine(bscale, bzero, blank, image, mask, maskVal, fits)) { + if (!psFitsScaleDetermine(bscale, bzero, boffset, 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."); @@ -331,5 +348,5 @@ } - return psFitsScaleForDisk(image, fits, *bscale, *bzero, rng); + return psFitsScaleForDisk(image, fits, *bscale, *bzero, *boffset, rng); } @@ -459,9 +476,18 @@ return NULL; } - psFree(info); if (floatType != PS_FITS_FLOAT_NONE) { outImage = psFitsFloatImageFromDisk(outImage, inImage, floatType); } + + // Need to apply BOFFSET if info->is_logscaled is true + if (info->is_logscaled) { + double boffset; + int status; + fits_read_key_dbl(fits->fd, "BOFFSET", &boffset, NULL, &status); + outImage = psFitsScaleFromDisk(outImage,boffset); + } + psFree(info); + psFree(inImage); @@ -572,7 +598,8 @@ double bscale = NAN, bzero = NAN; // Scale and zero point to put in header + double boffset = NAN; // Log offset to put into header. long blank = 0; // Blank (undefined) value for image psFitsFloat floatType; // Custom floating-point convention type - psImage *diskImage = imageToDiskRepresentation(&bscale, &bzero, &blank, &floatType, fits, image, + psImage *diskImage = imageToDiskRepresentation(&bscale, &bzero, &boffset, &blank, &floatType, fits, image, mask, maskVal, NULL, true); // Image to write out if (!diskImage) { @@ -610,4 +637,16 @@ psFitsOptions *options = fits->options; // FITS I/O options +/* if (options) { */ +/* if (options->scaling == PS_FITS_SCALE_LOG_RANGE) { */ +/* fprintf(stderr,"it has the scaling I expect\n"); */ +/* } */ +/* else { */ +/* fprintf(stderr,"it does nto have the scaling I expect\n"); */ +/* } */ +/* } */ +/* else { */ +/* fprintf(stderr,"options is null, apparently? \n"); */ +/* } */ + psAssert(!useRequestedScale || !options || bitPix == options->bitpix || options->bitpix == 0, "Something's not consistent"); @@ -651,4 +690,14 @@ } + // Remove any BOFFSET values that exist in the header if we are not using that scaling anymore + if (options&&(!((options->scaling == PS_FITS_SCALE_LOG_RANGE)|| + (options->scaling == PS_FITS_SCALE_LOG_STDEV_POSITIVE)|| + (options->scaling == PS_FITS_SCALE_LOG_STDEV_NEGATIVE)|| + (options->scaling == PS_FITS_SCALE_LOG_STDEV_BOTH)))) { + if (psMetadataLookup(header,"BOFFSET")) { + psMetadataRemoveKey(header,"BOFFSET"); + } + } + // write the header, if any. if (header && !psFitsWriteHeaderImage(fits, header, createPHU)) { @@ -668,4 +717,11 @@ fits_write_key_dbl(fits->fd, "BSCALE", bscale, 12, "Scaling: TRUE = BZERO + BSCALE * DISK", &status); + if (options&&(((options->scaling == PS_FITS_SCALE_LOG_RANGE)|| + (options->scaling == PS_FITS_SCALE_LOG_STDEV_POSITIVE)|| + (options->scaling == PS_FITS_SCALE_LOG_STDEV_NEGATIVE)|| + (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); + } if (psFitsError(status, true, "Could not write BSCALE/BZERO headers to file.")) { success = false; @@ -769,7 +825,8 @@ bool success = true; // Successful update? double bscale = NAN, bzero = NAN; // Scale and zero point to put in header + double boffset = NAN; // Log offset to put in header long blank = 0; // Blank (undefined) value for image psFitsFloat floatType; // Custom floating-point convention type - psImage *diskImage = imageToDiskRepresentation(&bscale, &bzero, &blank, &floatType, fits, input, + psImage *diskImage = imageToDiskRepresentation(&bscale, &bzero, &boffset, &blank, &floatType, fits, input, mask, maskVal, NULL, false); // Image to write out if (!diskImage) { Index: /branches/czw_branch/20101203/psLib/src/fits/psFitsScale.c =================================================================== --- /branches/czw_branch/20101203/psLib/src/fits/psFitsScale.c (revision 30330) +++ /branches/czw_branch/20101203/psLib/src/fits/psFitsScale.c (revision 30331) @@ -244,4 +244,5 @@ 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 @@ -252,4 +253,5 @@ psAssert(bscale, "impossible"); psAssert(bzero, "impossible"); + psAssert(boffset, "impossible"); psAssert(image, "impossible"); psAssert(options, "impossible"); @@ -257,7 +259,24 @@ psTrace("psLib.fits", 3, "Scaling image by logarithm statistics"); int numCols = image->numCols, numRows = image->numRows; // Size of image + + psImage *copy; - double offset = 99e99; - + *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) { @@ -267,6 +286,6 @@ psF32 value = image->data.F32[y][x]; if (isfinite(value)) { - if (value < offset) { - offset = value; + if (value < *boffset) { + *boffset = value; } } @@ -279,6 +298,6 @@ psF64 value = image->data.F64[y][x]; if (isfinite(value)) { - if (value < offset) { - offset = value; + if (value < *boffset) { + *boffset = value; } } @@ -292,7 +311,7 @@ } // We only need to offset images that go negative. -/* if (offset > 0.0) { */ -/* offset = 0.0; */ -/* } */ + if (*boffset > 0.0) { + *boffset = 0.0; + } // Write offset to header // How? @@ -303,8 +322,8 @@ 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)); - } - image->data.F32[y][x] = (log10( image->data.F32[y][x] - offset)); +/* 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)); } } @@ -314,5 +333,5 @@ 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)); - image->data.F64[y][x] = (log10( image->data.F64[y][x] - offset)); + copy->data.F64[y][x] = (log10( image->data.F64[y][x] - *boffset)); } } @@ -325,8 +344,9 @@ // Do regular scaling on the logarithm image - if (!scaleStdev(bscale, bzero, image, mask, maskVal, options)) { + 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; } @@ -334,4 +354,5 @@ 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 @@ -340,4 +361,5 @@ psAssert(bscale, "impossible"); psAssert(bzero, "impossible"); + psAssert(boffset, "impossible"); psAssert(image, "impossible"); psAssert(options, "impossible"); @@ -345,7 +367,24 @@ psTrace("psLib.fits", 3, "Scaling image by logarithm statistics"); int numCols = image->numCols, numRows = image->numRows; // Size of image + + psImage *copy; - double offset = 99e99; - + *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) { @@ -355,6 +394,6 @@ psF32 value = image->data.F32[y][x]; if (!isfinite(value)) { - if (value < offset) { - offset = value; + if (value < *boffset) { + *boffset = value; } } @@ -367,6 +406,6 @@ psF64 value = image->data.F64[y][x]; if (!isfinite(value)) { - if (value < offset) { - offset = value; + if (value < *boffset) { + *boffset = value; } } @@ -380,6 +419,6 @@ } // We only need to offset images that go negative. - if (offset > 0.0) { - offset = 0.0; + if (*boffset > 0.0) { + *boffset = 0.0; } // Write offset to header @@ -391,5 +430,5 @@ for (int y = 0; y < numRows; y++) { for (int x = 0; x < numCols; x++) { - image->data.F32[y][x] = (log10( image->data.F32[y][x] - offset)); + copy->data.F32[y][x] = (log10( image->data.F32[y][x] - *boffset)); } } @@ -398,5 +437,5 @@ for (int y = 0; y < numRows; y++) { for (int x = 0; x < numCols; x++) { - image->data.F64[y][x] = (log10( image->data.F64[y][x] - offset)); + copy->data.F64[y][x] = (log10( image->data.F64[y][x] - *boffset)); } } @@ -409,8 +448,9 @@ // Do regular scaling on the logarithm image - if (!scaleRange(bscale, bzero, image, options)) { + 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; } @@ -421,9 +461,10 @@ ////////////////////////////////////////////////////////////////////////////////////////////////////////////// -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); @@ -436,4 +477,5 @@ *bscale = NAN; *bzero = NAN; + *boffset = 0; *blank = 0; @@ -480,5 +522,5 @@ break; case PS_FITS_SCALE_LOG_RANGE: - if (!logscaleRange(bscale,bzero,image,options)) { + if (!logscaleRange(bscale,bzero,boffset,image,options)) { psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from range"); return false; @@ -488,5 +530,5 @@ case PS_FITS_SCALE_LOG_STDEV_NEGATIVE: case PS_FITS_SCALE_LOG_STDEV_BOTH: - if (!logscaleStdev(bscale, bzero, image, mask, maskVal, options)) { + if (!logscaleStdev(bscale, bzero,boffset, image, mask, maskVal, options)) { psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev"); return false; @@ -510,10 +552,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) { @@ -571,5 +614,14 @@ 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_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 */ \ @@ -617,16 +669,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 @@ -665,5 +712,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;; \ } \ } \ @@ -699,7 +746,4 @@ return out; } -#endif - - psFitsScaling psFitsScalingFromString(const char *string) Index: /branches/czw_branch/20101203/psLib/src/fits/psFitsScale.h =================================================================== --- /branches/czw_branch/20101203/psLib/src/fits/psFitsScale.h (revision 30330) +++ /branches/czw_branch/20101203/psLib/src/fits/psFitsScale.h (revision 30331) @@ -10,4 +10,5 @@ bool psFitsScaleDetermine(double *bscale, ///< Scaling, to return double *bzero, ///< Zero point, to return + double *boffset, ///< Log offset, to return long *blank, ///< Blank value, to return const psImage *image, ///< Image to scale @@ -27,7 +28,10 @@ double bscale, ///< Scaling double bzero, ///< Zero point + double boffset, ///< Log offset psRandom *rng ///< Random number generator (for the "fuzz"), or NULL ); - +psImage *psFitsScaleFromDisk(const psImage *image, ///< Image to to unapply BOFFSET + double boffset ///< Log offset + ); /// Interpret a string as a scaling method psFitsScaling psFitsScalingFromString(const char *string ///< String to interpret