IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 30331


Ignore:
Timestamp:
Jan 20, 2011, 6:35:02 PM (15 years ago)
Author:
watersc1
Message:

need to do some checks, but it looks like this fully implements log(flux) read/writes

Location:
branches/czw_branch/20101203/psLib/src/fits
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/czw_branch/20101203/psLib/src/fits/psFitsImage.c

    r30118 r30331  
    5454    int fitsDatatype;                   // cfitsio data type
    5555    int psDatatype;                     // psLib data type
     56    bool is_logscaled;                  // is this image log scaled using BOFFSET?
    5657} p_psFitsReadInfo;
    5758
    … …  
    128129
    129130    // Check scale and zero
    130     double bscale = 0.0, bzero = 0.0;    // Scale and zero point
     131    double bscale = 0.0, bzero = 0.0, boffset = NAN;    // Scale and zero point
    131132    if (fits_read_key_dbl(fits->fd, "BSCALE", &bscale, NULL, &status) && status != KEY_NO_EXIST) {
    132133        psFitsError(status, true, "Unable to read header.");
    … …  
    137138        psFitsError(status, true, "Unable to read header.");
    138139        goto bad;
     140    }
     141    status = 0;
     142    if (fits_read_key_dbl(fits->fd, "BOFFSET", &boffset, NULL, &status) && status != KEY_NO_EXIST) {
     143        psFitsError(status, true, "Unable to read header.");
     144        goto bad;
     145    }
     146    if (status == KEY_NO_EXIST) {
     147      info->is_logscaled = false;
     148    }
     149    else if (isfinite(boffset)) {
     150      info->is_logscaled = true;
     151    }
     152    else {
     153      info->is_logscaled = false;
    139154    }
    140155    status = 0;
    … …  
    246261static psImage *imageToDiskRepresentation(double *bscale, // Scaling applied
    247262                                          double *bzero, // Zero point applied
     263                                          double *boffset, // Log offset applied
    248264                                          long *blank, // Blank value (integer data)
    249265                                          psFitsFloat *floatType, // Type of custom floating-point
    … …  
    258274    psAssert(bscale, "impossible");
    259275    psAssert(bzero, "impossible");
     276    psAssert(boffset, "impossible");
    260277    psAssert(floatType, "impossible");
    261278    psAssert(fits, "impossible");
    … …  
    305322        if (newScaleZero) {
    306323            // Choose an appropriate BSCALE and BZERO
    307             if (!psFitsScaleDetermine(bscale, bzero, blank, image, mask, maskVal, fits)) {
     324          if (!psFitsScaleDetermine(bscale, bzero, boffset, blank, image, mask, maskVal, fits)) {
    308325                // We can't have the write dying for this reason --- try to save it somehow!
    309326                psWarning("Unable to determine BSCALE and BZERO for image --- refusing to quantise.");
    … …  
    331348        }
    332349
    333         return psFitsScaleForDisk(image, fits, *bscale, *bzero, rng);
     350        return psFitsScaleForDisk(image, fits, *bscale, *bzero, *boffset, rng);
    334351    }
    335352
    … …  
    459476        return NULL;
    460477    }
    461     psFree(info);
    462478
    463479    if (floatType != PS_FITS_FLOAT_NONE) {
    464480        outImage = psFitsFloatImageFromDisk(outImage, inImage, floatType);
    465481    }
     482
     483    // Need to apply BOFFSET if info->is_logscaled is true
     484    if (info->is_logscaled) {
     485      double boffset;
     486      int status;
     487      fits_read_key_dbl(fits->fd, "BOFFSET", &boffset, NULL, &status);
     488      outImage = psFitsScaleFromDisk(outImage,boffset);
     489    }
     490    psFree(info);
     491
    466492    psFree(inImage);
    467493
    … …  
    572598
    573599    double bscale = NAN, bzero = NAN;   // Scale and zero point to put in header
     600    double boffset = NAN;               // Log offset to put into header.
    574601    long blank = 0;                     // Blank (undefined) value for image
    575602    psFitsFloat floatType;              // Custom floating-point convention type
    576     psImage *diskImage = imageToDiskRepresentation(&bscale, &bzero, &blank, &floatType, fits, image,
     603    psImage *diskImage = imageToDiskRepresentation(&bscale, &bzero, &boffset, &blank, &floatType, fits, image,
    577604                                                   mask, maskVal, NULL, true); // Image to write out
    578605    if (!diskImage) {
    … …  
    610637
    611638    psFitsOptions *options = fits->options; // FITS I/O options
     639/*     if (options) { */
     640/*       if (options->scaling == PS_FITS_SCALE_LOG_RANGE) { */
     641/*      fprintf(stderr,"it has the scaling I expect\n"); */
     642/*       } */
     643/*       else { */
     644/*      fprintf(stderr,"it does nto have the scaling I expect\n"); */
     645/*       } */
     646/*     } */
     647/*     else { */
     648/*       fprintf(stderr,"options is null, apparently? \n"); */
     649/*     } */
     650     
    612651    psAssert(!useRequestedScale || !options || bitPix == options->bitpix || options->bitpix == 0,
    613652             "Something's not consistent");
    … …  
    651690    }
    652691
     692    // Remove any BOFFSET values that exist in the header if we are not using that scaling anymore
     693    if (options&&(!((options->scaling == PS_FITS_SCALE_LOG_RANGE)||
     694                   (options->scaling == PS_FITS_SCALE_LOG_STDEV_POSITIVE)||
     695                   (options->scaling == PS_FITS_SCALE_LOG_STDEV_NEGATIVE)||
     696                   (options->scaling == PS_FITS_SCALE_LOG_STDEV_BOTH)))) {
     697      if (psMetadataLookup(header,"BOFFSET")) {
     698        psMetadataRemoveKey(header,"BOFFSET");
     699      }
     700    }   
     701
    653702    // write the header, if any.
    654703    if (header && !psFitsWriteHeaderImage(fits, header, createPHU)) {
    … …  
    668717        fits_write_key_dbl(fits->fd, "BSCALE", bscale, 12,
    669718                           "Scaling: TRUE = BZERO + BSCALE * DISK", &status);
     719        if (options&&(((options->scaling == PS_FITS_SCALE_LOG_RANGE)||
     720                       (options->scaling == PS_FITS_SCALE_LOG_STDEV_POSITIVE)||
     721                       (options->scaling == PS_FITS_SCALE_LOG_STDEV_NEGATIVE)||
     722                       (options->scaling == PS_FITS_SCALE_LOG_STDEV_BOTH)))) {
     723          fits_write_key_dbl(fits->fd, "BOFFSET", boffset, 12,
     724                             "Scaling: TRUE = BZERO + BSCALE * 10**(DISK) + BOFFSET)", &status);
     725        }       
    670726        if (psFitsError(status, true, "Could not write BSCALE/BZERO headers to file.")) {
    671727            success = false;
    … …  
    769825    bool success = true;                // Successful update?
    770826    double bscale = NAN, bzero = NAN;   // Scale and zero point to put in header
     827    double boffset = NAN;               // Log offset to put in header
    771828    long blank = 0;                     // Blank (undefined) value for image
    772829    psFitsFloat floatType;              // Custom floating-point convention type
    773     psImage *diskImage = imageToDiskRepresentation(&bscale, &bzero, &blank, &floatType, fits, input,
     830    psImage *diskImage = imageToDiskRepresentation(&bscale, &bzero, &boffset, &blank, &floatType, fits, input,
    774831                                                   mask, maskVal, NULL, false); // Image to write out
    775832    if (!diskImage) {

  • branches/czw_branch/20101203/psLib/src/fits/psFitsScale.c

    r30255 r30331  
    244244static bool logscaleStdev(double *bscale, // Scaling, to return
    245245                          double *bzero, // Zero point, to return
     246                          double *boffset, // Log offset, to return
    246247                          const psImage *image, // Image to scale
    247248                          const psImage *mask, // Mask image
    … …  
    252253    psAssert(bscale, "impossible");
    253254    psAssert(bzero, "impossible");
     255    psAssert(boffset, "impossible");
    254256    psAssert(image, "impossible");
    255257    psAssert(options, "impossible");
    … …  
    257259    psTrace("psLib.fits", 3, "Scaling image by logarithm statistics");
    258260    int numCols = image->numCols, numRows = image->numRows; // Size of image
     261
     262    psImage *copy;
    259263   
    260     double offset = 99e99;
    261 
     264    *boffset = 99e99;
     265
     266    // Make a copy of the image to pass to get the scaling parameters
     267
     268    switch (image->type.type) {
     269    case PS_TYPE_F32:
     270      copy = psImageCopy(NULL,image,PS_TYPE_F32);
     271      break;
     272    case PS_TYPE_F64:
     273      copy = psImageCopy(NULL,image,PS_TYPE_F64);
     274      break;
     275    default:
     276      psError(PS_ERR_UNKNOWN, true, "Target type is not a float: %d", image->type.type);
     277      return NULL;
     278      break;
     279    }
     280   
    262281    // Determine the minimum value on this image.
    263282    switch (image->type.type) {
    … …  
    267286          psF32 value = image->data.F32[y][x];
    268287          if (isfinite(value)) {
    269             if (value < offset) {
    270               offset = value;
     288            if (value < *boffset) {
     289              *boffset = value;
    271290            }
    272291          }
    … …  
    279298          psF64 value = image->data.F64[y][x];
    280299          if (isfinite(value)) {
    281             if (value < offset) {
    282               offset = value;
     300            if (value < *boffset) {
     301              *boffset = value;
    283302            }
    284303          }
    … …  
    292311    }
    293312    // We only need to offset images that go negative.
    294 /*     if (offset > 0.0) { */
    295 /*       offset = 0.0; */
    296 /*     } */
     313    if (*boffset > 0.0) {
     314      *boffset = 0.0;
     315    }
    297316    // Write offset to header
    298317    // How?
    … …  
    303322      for (int y = 0; y < numRows; y++) {
    304323        for (int x = 0; x < numCols; x++) {
    305           if (x == 2331 && y == 2843) {
    306             fprintf(stderr,"psFS32: %d %d %g %g %g\n",x,y,offset,image->data.F32[y][x],log10(image->data.F32[y][x] - offset));
    307           }
    308           image->data.F32[y][x] = (log10( image->data.F32[y][x] - offset));
     324/*        if (x == 2331 && y == 2843) { */
     325/*          fprintf(stderr,"psFS32: %d %d %g %g %g\n",x,y,offset,image->data.F32[y][x],log10(image->data.F32[y][x] - offset)); */
     326/*        } */
     327          copy->data.F32[y][x] = (log10( image->data.F32[y][x] - *boffset));
    309328        }
    310329      }
    … …  
    314333          for (int x = 0; x < numCols; x++) {
    315334            //      fprintf(stderr,"psFS64: %d %d %g %g %g\n",x,y,offset,image->data.F64[y][x],log10(image->data.F64[y][x] - offset));
    316             image->data.F64[y][x] = (log10( image->data.F64[y][x] - offset));
     335            copy->data.F64[y][x] = (log10( image->data.F64[y][x] - *boffset));
    317336          }
    318337        }
    … …  
    325344     
    326345    // Do regular scaling on the logarithm image
    327     if (!scaleStdev(bscale, bzero, image, mask, maskVal, options)) {
     346    if (!scaleStdev(bscale, bzero, copy, mask, maskVal, options)) {
    328347      psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev");
    329348      return false;
    330349    }
     350    psFree(copy);
    331351    return true;
    332352}
    … …  
    334354static bool logscaleRange(double *bscale, // Scaling, to return
    335355                          double *bzero, // Zero point, to return
     356                          double *boffset, // Log offset, to return
    336357                          const psImage *image, // Image to scale
    337358                          const psFitsOptions *options // FITS options
    … …  
    340361    psAssert(bscale, "impossible");
    341362    psAssert(bzero, "impossible");
     363    psAssert(boffset, "impossible");
    342364    psAssert(image, "impossible");
    343365    psAssert(options, "impossible");
    … …  
    345367    psTrace("psLib.fits", 3, "Scaling image by logarithm statistics");
    346368    int numCols = image->numCols, numRows = image->numRows; // Size of image
     369
     370    psImage *copy;
    347371   
    348     double offset = 99e99;
    349 
     372    *boffset = 99e99;
     373
     374    // Make a copy of the image to pass to get the scaling parameters
     375
     376    switch (image->type.type) {
     377    case PS_TYPE_F32:
     378      copy = psImageCopy(NULL,image,PS_TYPE_F32);
     379      break;
     380    case PS_TYPE_F64:
     381      copy = psImageCopy(NULL,image,PS_TYPE_F64);
     382      break;
     383    default:
     384      psError(PS_ERR_UNKNOWN, true, "Target type is not a float: %d", image->type.type);
     385      return NULL;
     386      break;
     387    }
     388   
    350389    // Determine the minimum value on this image.
    351390    switch (image->type.type) {
    … …  
    355394          psF32 value = image->data.F32[y][x];
    356395          if (!isfinite(value)) {
    357             if (value < offset) {
    358               offset = value;
     396            if (value < *boffset) {
     397              *boffset = value;
    359398            }
    360399          }
    … …  
    367406          psF64 value = image->data.F64[y][x];
    368407          if (!isfinite(value)) {
    369             if (value < offset) {
    370               offset = value;
     408            if (value < *boffset) {
     409              *boffset = value;
    371410            }
    372411          }
    … …  
    380419    }
    381420    // We only need to offset images that go negative.
    382     if (offset > 0.0) {
    383       offset = 0.0;
     421    if (*boffset > 0.0) {
     422      *boffset = 0.0;
    384423    }
    385424    // Write offset to header
    … …  
    391430      for (int y = 0; y < numRows; y++) {
    392431        for (int x = 0; x < numCols; x++) {
    393           image->data.F32[y][x] = (log10( image->data.F32[y][x] - offset));
     432          copy->data.F32[y][x] = (log10( image->data.F32[y][x] - *boffset));
    394433        }
    395434      }
    … …  
    398437        for (int y = 0; y < numRows; y++) {
    399438          for (int x = 0; x < numCols; x++) {
    400             image->data.F64[y][x] = (log10( image->data.F64[y][x] - offset));
     439            copy->data.F64[y][x] = (log10( image->data.F64[y][x] - *boffset));
    401440          }
    402441        }
    … …  
    409448     
    410449    // Do regular scaling on the logarithm image
    411     if (!scaleRange(bscale, bzero, image, options)) {
     450    if (!scaleRange(bscale, bzero, copy, options)) {
    412451      psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev");
    413452      return false;
    414453    }
     454    psFree(copy);
    415455    return true;
    416456}
    … …  
    421461//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    422462
    423 bool psFitsScaleDetermine(double *bscale, double *bzero, long *blank, const psImage *image,
     463bool psFitsScaleDetermine(double *bscale, double *bzero, double *boffset, long *blank, const psImage *image,
    424464                          const psImage *mask, psImageMaskType maskVal, const psFits *fits)
    425465{
    426466    PS_ASSERT_PTR_NON_NULL(bscale, false);
    427467    PS_ASSERT_PTR_NON_NULL(bzero, false);
     468    PS_ASSERT_PTR_NON_NULL(boffset, false);
    428469    PS_ASSERT_PTR_NON_NULL(blank, false);
    429470    PS_ASSERT_IMAGE_NON_NULL(image, false);
    … …  
    436477    *bscale = NAN;
    437478    *bzero = NAN;
     479    *boffset = 0;
    438480    *blank = 0;
    439481
    … …  
    480522        break;
    481523    case PS_FITS_SCALE_LOG_RANGE:
    482       if (!logscaleRange(bscale,bzero,image,options)) {
     524      if (!logscaleRange(bscale,bzero,boffset,image,options)) {
    483525        psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from range");
    484526        return false;
    … …  
    488530    case PS_FITS_SCALE_LOG_STDEV_NEGATIVE:
    489531    case PS_FITS_SCALE_LOG_STDEV_BOTH:
    490       if (!logscaleStdev(bscale, bzero, image, mask, maskVal, options)) {
     532      if (!logscaleStdev(bscale, bzero,boffset, image, mask, maskVal, options)) {
    491533            psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev");
    492534            return false;
    … …  
    510552
    511553
    512     psTrace("psLib.fits", 3, "BSCALE = %.10lf, BZERO = %.10lf, BLANK = %ld\n", *bscale, *bzero, *blank);
     554    psTrace("psLib.fits", 3, "BSCALE = %.10lf, BZERO = %.10lf, BOFFSET = %.10lf, BLANK = %ld\n",
     555            *bscale, *bzero, *boffset, *blank);
    513556    return true;
    514557}
    515558
    516559
    517 psImage *psFitsScaleForDisk(const psImage *image, const psFits *fits, double bscale, double bzero,
     560psImage *psFitsScaleForDisk(const psImage *image, const psFits *fits, double bscale, double bzero, double boffset,
    518561                            psRandom *rng)
    519562{
    … …  
    571614        for (int y = 0; y < numRows; y++) { \
    572615            for (int x = 0; x < numCols; x++) { \
    573                 ps##INTYPE value = (IN)->data.INTYPE[y][x]; \
     616              ps##INTYPE value; \
     617              if ((options->scaling == PS_FITS_SCALE_LOG_RANGE)||       \
     618                  (options->scaling == PS_FITS_SCALE_LOG_STDEV_POSITIVE)|| \
     619                  (options->scaling == PS_FITS_SCALE_LOG_STDEV_NEGATIVE)|| \
     620                  (options->scaling == PS_FITS_SCALE_LOG_STDEV_BOTH)) { \
     621                value = log10( (IN)->data.INTYPE[y][x] - boffset ); \
     622              }                                                         \
     623                else { \
     624                  value = (IN)->data.INTYPE[y][x];      \
     625                }                                           \
    574626                if (!isfinite(value)) { \
    575627                    /* This choice of "max" for non-finite pixels is mainly cosmetic --- it has to be */ \
    … …  
    617669
    618670
    619 #if 0
     671
    620672// This function to apply BSCALE and BZERO to an image read immediately from disk should not be necessary at
    621673// the present time, since cfitsio should apply the scaling itself in the process of reading.  However, we may
    622674// later desire it (e.g., if we ever make our own FITS implementation).
    623 psImage *psFitsScaleFromDisk(psFits *fits, psImage *image)
     675psImage *psFitsScaleFromDisk(const psImage *image, double boffset)
    624676{
    625677    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
    626 
    627     if (bscale == 0.0) {
    628         // BSCALE = 0 means don't apply anything
    629         return psMemIncrRefCounter(image);
    630     }
    631678
    632679    psElemType inType = image->type.type; // Type for input image
    … …  
    665712      for (int y = 0; y < numRows; y++) { \
    666713          for (int x = 0; x < numCols; x++) { \
    667               out->data.OUTTYPE[y][x] = image->data.INTYPE[y][x] * bscale + bzero; \
     714            out->data.OUTTYPE[y][x] = pow(10,image->data.INTYPE[y][x]) + boffset;; \
    668715          } \
    669716      } \
    … …  
    699746    return out;
    700747}
    701 #endif
    702 
    703 
    704748
    705749psFitsScaling psFitsScalingFromString(const char *string)

  • branches/czw_branch/20101203/psLib/src/fits/psFitsScale.h

    r21183 r30331  
    1010bool psFitsScaleDetermine(double *bscale, ///< Scaling, to return
    1111                          double *bzero, ///< Zero point, to return
     12                          double *boffset, ///< Log offset, to return
    1213                          long *blank,  ///< Blank value, to return
    1314                          const psImage *image, ///< Image to scale
    … …  
    2728                            double bscale, ///< Scaling
    2829                            double bzero, ///< Zero point
     30                            double boffset, ///< Log offset
    2931                            psRandom *rng ///< Random number generator (for the "fuzz"), or NULL
    3032    );
    31 
     33psImage *psFitsScaleFromDisk(const psImage *image, ///< Image to to unapply BOFFSET
     34                             double boffset        ///< Log offset
     35                             );
    3236/// Interpret a string as a scaling method
    3337psFitsScaling psFitsScalingFromString(const char *string ///< String to interpret
Note: See TracChangeset for help on using the changeset viewer.