Index: trunk/ppSub/src/ppSubVarianceRescale.c
===================================================================
--- trunk/ppSub/src/ppSubVarianceRescale.c	(revision 24672)
+++ trunk/ppSub/src/ppSubVarianceRescale.c	(revision 25373)
@@ -28,8 +28,24 @@
     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,168 +54,4 @@
     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);
 }