Changeset 27838 for branches/tap_branches/psLib/src/fits/psFitsScale.c

Timestamp:

May 3, 2010, 8:41:49 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/tap_branches

Files:

• . (modified) (1 prop)
• psLib (modified) (1 prop)
• psLib/src/fits/psFitsScale.c (modified) (4 diffs)

branches/tap_branches/psLib/src/fits/psFitsScale.c

@@ -15 +15 @@
 #include "psImage.h"
 #include "psFits.h"
+#include "psStats.h"
+#include "psImageStats.h"
 #include "psImageBackground.h"
-#include "psStats.h"
 #include "psImageStructManip.h"
 
@@ -29 +30 @@
 #define MEAN_STAT PS_STAT_ROBUST_MEDIAN // Statistic to use for mean
 #define STDEV_STAT PS_STAT_ROBUST_STDEV // Statistic to use for stdev
+
+#define DESPERATE_MEAN_STAT PS_STAT_SAMPLE_MEDIAN // Statistic to use for mean when deperate
+#define DESPERATE_STDEV_STAT PS_STAT_SAMPLE_QUARTILE // Statistic to use for stdev when desperate
 
 
@@ -118 +122 @@
     psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
     psStats *stats = psStatsAlloc(MEAN_STAT | STDEV_STAT); // Statistics object
+    double mean, stdev;                                    // Mean and standard deviation
     if (!psImageBackground(stats, NULL, image, mask, maskVal, rng)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to perform statistics on image");
-        psFree(rng);
-        psFree(stats);
-        return false;
+        // It could be because the image is entirely masked, in which case we don't want to error
+        bool good = false;              // Any good pixels?
+
+
+// Find good pixels in an image, by image type
+#define GOOD_PIXELS_CASE(TYPE) \
+      case PS_TYPE_##TYPE: \
+        for (int y = 0; y < image->numRows && !good; y++) { \
+            for (int x = 0; x < image->numCols && !good; x++) { \
+                if (mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal)) { \
+                    continue; \
+                } \
+                if (!isfinite(image->data.TYPE[y][x])) { \
+                    continue; \
+                } \
+                good = true; \
+            } \
+        } \
+        break;
+
+        switch (image->type.type) {
+            GOOD_PIXELS_CASE(F32);
+            GOOD_PIXELS_CASE(F64);
+          default:
+            psAbort("Unsupported case: %x", image->type.type);
+        }
+
+        if (!good) {
+            psLogMsg("psLib.fits", PS_LOG_DETAIL, "Image has no good pixels, setting BSCALE = 1, BZERO = 0");
+            psErrorClear();
+            *bscale = 1.0;
+            *bzero = 0.0;
+            psFree(rng);
+            psFree(stats);
+            return true;
+        }
+
+        // There are some good pixels in there somewhere; psImageBackground just didn't find them
+        psLogMsg("psLib.fits", PS_LOG_DETAIL,
+                 "Couldn't measure background statistics for image quantisation; retrying.");
+        psErrorClear();
+        // Retry using all the available pixels
+        stats->nSubsample = image->numCols * image->numRows + 1;
+        if (!psImageStats(stats, image, mask, maskVal)) {
+            psLogMsg("psLib.fits", PS_LOG_DETAIL,
+                     "Couldn't measure background statistics for image quantisation (attempt 2); retrying.");
+            psErrorClear();
+            // Retry with desperate statistic
+            stats->options = DESPERATE_MEAN_STAT | DESPERATE_STDEV_STAT;
+            if (!psImageStats(stats, image, mask, maskVal)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to measure background statistics for image");
+                psFree(rng);
+                psFree(stats);
+                return false;
+            } else {
+                // Desperate retry
+                mean = psStatsGetValue(stats, DESPERATE_MEAN_STAT);
+                stdev = psStatsGetValue(stats, DESPERATE_STDEV_STAT);
+            }
+        } else {
+            // Retry with all available pixels
+            mean = psStatsGetValue(stats, MEAN_STAT);
+            stdev = psStatsGetValue(stats, STDEV_STAT);
+        }
+    } else {
+        // First attempt
+        mean = psStatsGetValue(stats, MEAN_STAT);
+        stdev = psStatsGetValue(stats, STDEV_STAT);
     }
     psFree(rng);
-
-    double mean = psStatsGetValue(stats, MEAN_STAT); // Mean
-    double stdev = psStatsGetValue(stats, STDEV_STAT); // Standard deviation
     psFree(stats);
     if (!isfinite(mean) || !isfinite(stdev)) {
@@ -318 +384 @@
                 } else { \
                     value = (value - zero) * scale; \
-                    if (options->fuzz) { \
+                    if (options->fuzz && (value - (int)value != 0.0)) { \
                        /* Add random factor [-0.5,0.5): adds a variance of 1/12, */ \
                        /* but preserves the expectation value */ \