Changeset 30677

Timestamp:

Feb 17, 2011, 3:22:41 PM (15 years ago)

Author:

eugene

Message:

merge changes from the trunk

Location:

branches/eam_branches/ipp-20110213/psLib/src/fits

Files:

• psFits.h (modified) (2 diffs)
• psFitsImage.c (modified) (13 diffs)
• psFitsScale.c (modified) (13 diffs)
• psFitsScale.h (modified) (2 diffs)

branches/eam_branches/ipp-20110213/psLib/src/fits/psFits.h

@@ -50 +50 @@
     PS_FITS_SCALE_STDEV_NEGATIVE,       ///< Auto-scale to sample stdev, place mean at upper limit
     PS_FITS_SCALE_STDEV_BOTH,           ///< Auto-scale to sample stdev, place mean at middle
-    PS_FITS_SCALE_MANUAL                ///< Manual scaling (use specified BSCALE and BZERO)
+    PS_FITS_SCALE_MANUAL,                ///< Manual scaling (use specified BSCALE and BZERO)
+    PS_FITS_SCALE_LOG_RANGE,            ///< Take logarithm, Auto-scale to preserve dynamic range
+    PS_FITS_SCALE_LOG_STDEV_POSITIVE,   ///< Take logarithm, Auto-scale to sample stdev, place mean at lower limit
+    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)
 } psFitsScaling;
 
@@ -66 +71 @@
     bool fuzz;                          ///< Fuzz the values when quantising floating-point values?
     double bscale, bzero;               ///< Manually specified BSCALE and BZERO (for SCALE_MANUAL)
+    double boffset;                     ///< Manually specified BOFFSET (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_*)

branches/eam_branches/ipp-20110213/psLib/src/fits/psFitsImage.c

@@ -54 +54 @@
     int fitsDatatype;                   // cfitsio data type
     int psDatatype;                     // psLib data type
+    bool is_logscaled;                  // is this image log scaled using BOFFSET?
 } p_psFitsReadInfo;
 
@@ -128 +129 @@
 
     // 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 +138 @@
         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 +261 @@
 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 +274 @@
     psAssert(bscale, "impossible");
     psAssert(bzero, "impossible");
+    psAssert(boffset, "impossible");
     psAssert(floatType, "impossible");
     psAssert(fits, "impossible");
@@ -305 +322 @@
         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 +348 @@
         }
 
-        return psFitsScaleForDisk(image, fits, *bscale, *bzero, rng);
+        return psFitsScaleForDisk(image, fits, *bscale, *bzero, *boffset, rng);
     }
 
@@ -459 +476 @@
         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 +598 @@
 
     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 +637 @@
 
     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 +690 @@
     }
 
+    // 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_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)))) {
+      if (psMetadataLookup(header,"BOFFSET")) {
+        psMetadataRemoveKey(header,"BOFFSET");
+      }
+    }   
+
     // write the header, if any.
     if (header && !psFitsWriteHeaderImage(fits, header, createPHU)) {
@@ -668 +718 @@
         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_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)))) {
+          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 +827 @@
     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) {

branches/eam_branches/ipp-20110213/psLib/src/fits/psFitsScale.c

@@ -100 +100 @@
 }
 
+    
 // Determine appropriate BSCALE and BZERO for an image, mapping the standard deviation to the nominated number
 // of bits
@@ -212 +213 @@
     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 +219 @@
         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 +225 @@
         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 +241 @@
 }
 
+
+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 +477 @@
     *bscale = NAN;
     *bzero = NAN;
+    *boffset = 0;
     *blank = 0;
 
@@ -299 +521 @@
         }
         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 +557 @@
 
 
-    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 +619 @@
         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 +675 @@
 
 
-#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 +718 @@
       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 +752 @@
     return out;
 }
-#endif
-
-
 
 psFitsScaling psFitsScalingFromString(const char *string)
@@ -517 +761 @@
     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;
 

branches/eam_branches/ipp-20110213/psLib/src/fits/psFitsScale.h

@@ -10 +10 @@
 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 +28 @@
                             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