Changeset 25373

Timestamp:

Sep 14, 2009, 4:15:46 PM (17 years ago)

Author:

Paul Price

Message:

Using pmReadoutVarianceRenormalise

File:

• trunk/ppSub/src/ppSubVarianceRescale.c (modified) (2 diffs)

trunk/ppSub/src/ppSubVarianceRescale.c

@@ -28 +28 @@
     bool mdok; // Status of metadata lookups
 
-    psMetadata *ppSubRecipe = psMetadataLookupPtr(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSub
-    psAssert(ppSubRecipe, "Need PPSUB recipe");
+    psMetadata *recipe = psMetadataLookupPtr(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSub
+    psAssert(recipe, "Need PPSUB recipe");
 
-    if (!psMetadataLookupBool(&mdok, ppSubRecipe, "VARIANCE.RESCALE")) return true; 
+    if (!psMetadataLookupBool(&mdok, recipe, "RENORM")) return true;
+
+    int num = psMetadataLookupS32(&mdok, recipe, "RENORM.NUM");
+    if (!mdok) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "RENORM.NUM is not set in the recipe");
+        return false;
+    }
+    float minValid = psMetadataLookupF32(&mdok, recipe, "RENORM.MIN");
+    if (!mdok) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "RENORM.MIN is not set in the recipe");
+        return false;
+    }
+    float maxValid = psMetadataLookupF32(&mdok, recipe, "RENORM.MAX");
+    if (!mdok) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "RENORM.MAX is not set in the recipe");
+        return false;
+    }
 
     psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask
@@ -38 +54 @@
     pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT"); // Output image
 
-    psImage *image    = outRO->image;    // Image of interest
-    psImage *variance = outRO->variance; // Variance image
-    psImage *mask     = outRO->mask;     // Mask of interest
-
-    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
-
-    psLogMsg("psLib.psImageBackground", PS_LOG_DETAIL, "Generating the signal-to-noise image");
-
-    // generate an image representing the signal/noise:
-    psImage *SN = psImageCopy (NULL, outRO->image, PS_TYPE_F32); 
-
-    int nx = image->numCols; // Size of image
-    int ny = image->numRows; // Size of image
-
-    for (int y = 0; y < ny; y++) {
-        for (int x = 0; x < nx; x++) {
-            if (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskBad) {
-                SN->data.F32[y][x] = 0.0;
-                continue;
-            } 
-            if (!isfinite(image->data.F32[y][x])) {
-                SN->data.F32[y][x] = 0.0;
-                continue;
-            } 
-            if (!isfinite(variance->data.F32[y][x]) || (variance->data.F32[y][x] < 0.0)) {
-                SN->data.F32[y][x] = 0.0;
-                continue;
-            } 
-            SN->data.F32[y][x] = image->data.F32[y][x] / sqrt(variance->data.F32[y][x]);
-        }
-    }
-
-    // these should not be chosen by the user: only a ROBUST version makes sense
-    psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
-
-    int Nsubset = psMetadataLookupS32(&mdok, ppSubRecipe, "VARIANCE.RESCALE.NSAMPLE");
-    psAssert (mdok, "missing VARIANCE.RESCALE.NSAMPLE in ppSub recipe");
-    
-    const int Npixels = nx*ny; // Total number of pixels
-    Nsubset = PS_MIN(Nsubset, Npixels); // Number of pixels in subset
-    psLogMsg("psLib.psImageBackground", PS_LOG_DETAIL, "Searching for %d pixels in S/N image", Nsubset);
-
-    psVector *values = psVectorAllocEmpty(Nsubset, PS_TYPE_F32); // Vector containing subsample
-
-    // Minimum and maximum values
-    float min = +PS_MAX_F32;
-    float max = -PS_MAX_F32;
-
-    // select a subset of the image pixels to measure the stats
-    long n = 0;                         // Number of actual pixels in subset
-    if (Nsubset >= Npixels) {
-        // if we have an image smaller than Nsubset, just loop over the image pixels
-        for (int iy = 0; iy < ny; iy++) {
-            for (int ix = 0; ix < nx; ix++) {
-                if (!isfinite(image->data.F32[iy][ix]) || (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskBad)) {
-                    continue;
-                }
-
-                float value = SN->data.F32[iy][ix];
-                min = PS_MIN(value, min);
-                max = PS_MAX(value, max);
-                values->data.F32[n] = value;
-                n++;
-            }
-        }
-    } else {
-        for (long i = 0; i < Nsubset; i++) {
-            double frnd = psRandomUniform(rng);
-            int pixel = Npixels * frnd;
-            int ix = pixel % nx;
-            int iy = pixel / nx;
-
-            if (!isfinite(image->data.F32[iy][ix]) || (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskBad)) {
-                continue;
-            }
-
-            float value = SN->data.F32[iy][ix];
-            min = PS_MIN(value, min);
-            max = PS_MAX(value, max);
-            values->data.F32[n] = value;
-            n++;
-        }
-    }
-    if (n < 0.01*Nsubset) {
-        psLogMsg("psLib.psImageBackground", PS_LOG_DETAIL, "Unable to measure image background: too few data points (%ld)", n);
-        psFree(values);
-        return false;
-    }
-    values->n = n;
-    psFree (SN);
-
-    psLogMsg("psLib.psImageBackground", PS_LOG_DETAIL, "Using for %ld pixels in S/N image", values->n);
-
-    if (!psVectorStats (stats, values, NULL, NULL, 0)) {
-        if (psTraceGetLevel("psLib.imageops") >= 5) {
-            FILE *f = fopen ("vector.dat", "w");
-            int fd = fileno(f);
-            p_psVectorPrint (fd, values, "values");
-            fclose (f);
-        }
-        psError(PS_ERR_UNKNOWN, false, "Unable to measure statistics for image background "
-                "(%dx%d, (row0,col0) = (%d,%d)",
-                image->numRows, image->numCols, image->row0, image->col0);
-        psFree(values);
-        return false;
-    }
-    if (psTraceGetLevel("psLib.imageops") >= 6) {
-        FILE *f = fopen ("vector.dat", "w");
-        int fd = fileno(f);
-        p_psVectorPrint (fd, values, "values");
-        fclose (f);
-    }
-    psFree (values);
-    psLogMsg ("ppSub", PS_LOG_INFO, "background stdev %f\n", stats->robustStdev);
-
-    float minValid = psMetadataLookupF32(&mdok, ppSubRecipe, "VARIANCE.RESCALE.MIN.VALID");
-    psAssert (mdok, "missing VARIANCE.RESCALE.MIN.VALID in ppSub recipe");
-
-    float maxValid = psMetadataLookupF32(&mdok, ppSubRecipe, "VARIANCE.RESCALE.MAX.VALID");
-    psAssert (mdok, "missing VARIANCE.RESCALE.MAX.VALID in ppSub recipe");
-
-    psLogMsg("psLib.psImageBackground", PS_LOG_DETAIL, "Stdev of S/N image: %f (mean: %f)", stats->robustStdev, stats->robustMedian);
-
-    // valid range of correction; 0.66 to 1.5 (skip if in range 0.98 to 1.02)
-    if ((stats->robustStdev < minValid) || (stats->robustStdev > maxValid)) {
-        psWarning ("background stdev (%f) is out of allowed range; not correcting\n", stats->robustStdev);
-        // XXX set quality to poor
-        psFree (stats);
-        return true;
-    }
-
-    // valid range of correction; 0.66 to 1.5 (skip if in range 0.98 to 1.02)
-    // if ((stats->robustStdev > 0.98) && (stats->robustStdev > 1.02)) {
-    //  psLogMsg ("ppSub", PS_LOG_INFO, "background stdev (%f) is close to 1.0; do not modify variance\n", stats->robustStdev);
-    //  // XXX set quality to poor
-    //  psFree (stats);
-    //  return true;
-    // }
-
-    float varianceFactor = PS_SQR(stats->robustStdev);
-    psLogMsg ("ppSub", PS_LOG_INFO, "background stdev is not 1.0: multiplying variance by %f\n", varianceFactor);
-    for (int y = 0; y < ny; y++) {
-        for (int x = 0; x < nx; x++) {
-            if (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskBad) {
-                continue;
-            } 
-            if (!isfinite(image->data.F32[y][x])) {
-                continue;
-            } 
-            if (!isfinite(variance->data.F32[y][x]) || (variance->data.F32[y][x] < 0.0)) {
-                continue;
-            } 
-            variance->data.F32[y][x] *= varianceFactor;
-        }
-    }
-
-    pmFPA *outFPA = outRO->parent->parent->parent;
-    pmHDU *outHDU = outFPA->hdu;        // Output HDU
-
-    char history[80];
-    snprintf (history, 80, "Rescaled variance by %6.4f (stdev by %6.4f)", varianceFactor, stats->robustStdev);
-    psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, NULL, history);
-
-    psFree (stats);
-    return true;
+    return pmReadoutVarianceRenormalise(readout, maskBad, num, minValid, maxValid);
 }