Changeset 30119

Timestamp:

Dec 20, 2010, 4:25:10 PM (15 years ago)

Author:

watersc1

Message:

log(flux) compressed image work.

Location:

branches/czw_branch/20101203

Files:

• . (modified) (1 prop)
• ippconfig/recipes/fitstypes.mdc (modified) (1 diff)
• psLib/src/fits/psFits.h (modified) (1 diff)
• psLib/src/fits/psFitsScale.c (modified) (3 diffs)
• psLib/src/imageops (modified) (1 prop)
• psLib/test/math (modified) (1 prop)
• psModules/src/config/pmConfig.c (modified) (1 diff)
• psModules/src/objects (modified) (1 prop)
• psphot/src/psphotApResid.c (modified) (1 prop)
• psphot/src/psphotMakeFluxScale.c (modified) (1 prop)

branches/czw_branch/20101203/ippconfig/recipes/fitstypes.mdc

@@ -98 +98 @@
         NOISE           S32     8
 END
+# Compressed log flux image for stacks.
+COMP_STACK     METADATA
+        BITPIX          S32     16
+        SCALING         STR     LOG_STDEV_POSITIVE
+        STDEV.BITS      S32     4
+        STDEV.NUM       F32     10
+        COMPRESSION     STR     RICE
+        TILE.X          S32     0
+        TILE.Y          S32     1
+        TILE.Z          S32     1
+        NOISE           S32     8
+END
 
 # Compressed exposure image

branches/czw_branch/20101203/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_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_MANUAL                ///< Manual scaling (use specified BSCALE and BZERO)
 } psFitsScaling;

branches/czw_branch/20101203/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
@@ -237 +238 @@
 }
 
+
+static bool logscaleStdev(double *bscale, // Scaling, to return
+                          double *bzero, // Zero point, 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(image, "impossible");
+    psAssert(options, "impossible");
+
+    psTrace("psLib.fits", 3, "Scaling image by logarithm statistics");
+    int numCols = image->numCols, numRows = image->numRows; // Size of image
+    
+    double offset = 99e99;
+
+    // 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 < offset) {
+              offset = 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 < offset) {
+              offset = 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 (offset > 0.0) {
+      offset = 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++) {
+          image->data.F32[y][x] = (log10( image->data.F32[y][x] - offset));
+        }
+      }
+      break;
+    case PS_TYPE_F64:
+        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));
+          }
+        }
+      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, image, mask, maskVal, options)) {
+      psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev");
+      return false;
+    }
+    return true;
+}
+
+static bool logscaleRange(double *bscale, // Scaling, to return
+                          double *bzero, // Zero point, to return
+                          const psImage *image, // Image to scale
+                          const psFitsOptions *options // FITS options
+    )
+{
+    psAssert(bscale, "impossible");
+    psAssert(bzero, "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
+    
+    double offset = 99e99;
+
+    // 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 < offset) {
+              offset = 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 < offset) {
+              offset = 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 (offset > 0.0) {
+      offset = 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++) {
+          image->data.F32[y][x] = (log10( image->data.F32[y][x] - offset));
+        }
+      }
+      break;
+    case PS_TYPE_F64:
+        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));
+          }
+        }
+      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, image, options)) {
+      psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev");
+      return false;
+    }
+    return true;
+}
+
+
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
 // Public functions
@@ -295 +468 @@
       case PS_FITS_SCALE_STDEV_BOTH:
         if (!scaleStdev(bscale, bzero, 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,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, image, mask, maskVal, options)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev");
             return false;

branches/czw_branch/20101203/psModules/src/config/pmConfig.c

@@ -897 +897 @@
             psMetadataAddMetadata(filerules, PS_LIST_TAIL, old, PS_META_REPLACE,
                                   "Original replaced by -F option", newRule);
+        }
+        psFree(camerasIter);
+    }
+
+    // Look for command-line options for files to replace
+    while ((argNum = psArgumentGet(*argc, argv, "-R")) > 0) {
+        psArgumentRemove(argNum, argc, argv);
+        if (argNum + 2 >= *argc) {
+            psError(PM_ERR_CONFIG, true,
+                    "Filerule element switch (-R) provided without filerule element and value.");
+            psFree(config);
+            return NULL;
+        }
+
+        const char *rulename = argv[argNum]; // The filerule, to be modified
+        psArgumentRemove(argNum, argc, argv);
+        const char *element  = argv[argNum]; // The element, to be modified
+        psArgumentRemove(argNum, argc, argv);
+        const char *value    = argv[argNum]; // The value, to be set
+        psArgumentRemove(argNum, argc, argv);
+
+        psMetadata *cameras = psMetadataLookupMetadata(NULL, config->system, "CAMERAS"); // List of cameras
+        if (!cameras) {
+            psError(PM_ERR_CONFIG, false, "Unable to find CAMERAS in the site configuration.\n");
+            return false;
+        }
+
+        psMetadataIterator *camerasIter = psMetadataIteratorAlloc(cameras, PS_LIST_HEAD, NULL); // Iterator
+        psMetadataItem *cameraItem;     // Item from iteration
+        while ((cameraItem = psMetadataGetAndIncrement(camerasIter))) {
+            // Silently ignore problems --- they will be caught later, because if the user wants the nominated
+            // file and it's not available for that camera, then they will know.
+
+            if (cameraItem->type != PS_DATA_METADATA) {
+                psTrace("psModules.config", 2,
+                        "Entry %s in CAMERAS is not of type METADATA --- ignored.", cameraItem->name);
+                continue;
+            }
+            psMetadata *camera = cameraItem->data.md; // Camera configuration
+
+            psMetadata *newRule = pmConfigFileRule(config, camera, rulename); // The rule of interest
+            if (!newRule) {
+                psTrace("psModules.config", 2,
+                        "Unable to find filerule %s in camera %s --- ignored.", rulename, cameraItem->name);
+                continue;
+            }
+
+            // By calling pmConfigFileRule, we've assured that the FILERULES is now a metadata
+            psMetadata *filerules = psMetadataLookupMetadata(NULL, camera, "FILERULES"); // File rules
+            if (!filerules) {
+                psTrace("psModules.config", 2,
+                        "Can't find FILERULES of type METADATA in camera %s --- ignored.", cameraItem->name);
+                continue;
+            }
+
+            // Convert newRule to have the element value requested.
+            if (!psMetadataLookupStr(NULL,newRule,element)) {
+              psTrace("psModules.config", 2,
+                      "Unable to find filerule element %s in filerule %s in camera %s --- ignored.",
+                      element,rulename,cameraItem->name);
+              continue;
+            }
+            psMetadataAddStr(newRule, PS_LIST_TAIL, element, PS_META_REPLACE,
+                             "Original replaced by -R option", value);
+            
+            psMetadataAddMetadata(filerules, PS_LIST_TAIL, rulename, PS_META_REPLACE,
+                                  "Original replaced by -R option", newRule);
         }
         psFree(camerasIter);