Index: trunk/psModules/src/camera/pmFPAMaskWeight.c
===================================================================
--- trunk/psModules/src/camera/pmFPAMaskWeight.c	(revision 24767)
+++ trunk/psModules/src/camera/pmFPAMaskWeight.c	(revision 25372)
@@ -111,6 +111,6 @@
         // psError(PS_ERR_IO, true, "CELL.SATURATION is not set --- unable to set mask.\n");
         // return false;
-	psWarning("CELL.SATURATION is not set --- completely masking cell.\n");
-	saturation = NAN;
+        psWarning("CELL.SATURATION is not set --- completely masking cell.\n");
+        saturation = NAN;
     }
     float bad = psMetadataLookupF32(&mdok, cell->concepts, "CELL.BAD"); // Bad level
@@ -118,10 +118,10 @@
         // psError(PS_ERR_IO, true, "CELL.BAD is not set --- unable to set mask.\n");
         // return false;
-	psWarning("CELL.BAD is not set --- completely masking cell.\n");
-	bad = NAN;
+        psWarning("CELL.BAD is not set --- completely masking cell.\n");
+        bad = NAN;
     }
     psTrace("psModules.camera", 5, "Saturation: %f, bad: %f\n", saturation, bad);
 
-    // if CELL.GAIN or CELL.READNOISE are not set, then the variance will be set to NAN; 
+    // if CELL.GAIN or CELL.READNOISE are not set, then the variance will be set to NAN;
     // in this case, we have to set the mask as well
     float gain = psMetadataLookupF32(&mdok, cell->concepts, "CELL.GAIN"); // Cell gain
@@ -140,6 +140,6 @@
     // completely mask if SATURATION or BAD are invalid
     if (isnan(saturation) || isnan(bad) || isnan(gain) || isnan(readnoise)) {
-	psImageInit(mask, badMask);
-	return true;
+        psImageInit(mask, badMask);
+        return true;
     }
 
@@ -230,5 +230,5 @@
         // return false;
         psWarning("CELL.GAIN is not set --- setting variance to NAN\n");
-	gain = NAN;
+        gain = NAN;
     }
     float readnoise = psMetadataLookupF32(&mdok, cell->concepts, "CELL.READNOISE"); // Cell read noise
@@ -237,5 +237,5 @@
         // return false;
         psWarning("CELL.READNOISE is not set --- setting variance to NAN\n");
-	readnoise = NAN;
+        readnoise = NAN;
     }
     // if we have a non-NAN readnoise, then we need to ensure it has been updated (not necessary if NAN)
@@ -248,10 +248,10 @@
     if (isnan(gain) || isnan(readnoise)) {
         if (!readout->variance) {
-	    // generate the image if needed
+            // generate the image if needed
             readout->variance = psImageAlloc(readout->image->numCols, readout->image->numRows, PS_TYPE_F32);
         }
-	// XXX need to set the mask, if defined
+        // XXX need to set the mask, if defined
         psImageInit(readout->variance, NAN);
-	return true;
+        return true;
     }
 
@@ -262,5 +262,5 @@
 
         // a negative variance is non-sensical. if the image value drops below 1, the variance must be 1.
-	// XXX this calculation is wrong: limit is 1 e-, but this is in DN
+        // XXX this calculation is wrong: limit is 1 e-, but this is in DN
         readout->variance = (psImage*)psUnaryOp(readout->variance, readout->variance, "abs");
         readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "max",
@@ -276,9 +276,9 @@
     // apply a supplied readnoise map (NOTE: in DN, not electrons):
     if (noiseMap) {
-	psImage *rdVar = (psImage*)psBinaryOp(NULL, (const psPtr) noiseMap, "*", (const psPtr) noiseMap);
-	readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "+", rdVar);
-	psFree (rdVar);
+        psImage *rdVar = (psImage*)psBinaryOp(NULL, (const psPtr) noiseMap, "*", (const psPtr) noiseMap);
+        readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "+", rdVar);
+        psFree (rdVar);
     } else {
-	readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "+", psScalarAlloc(readnoise*readnoise/gain/gain, PS_TYPE_F32));
+        readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "+", psScalarAlloc(readnoise*readnoise/gain/gain, PS_TYPE_F32));
     }
 
@@ -362,6 +362,6 @@
 
 
-bool pmReadoutVarianceRenormPixels(const pmReadout *readout, psImageMaskType maskVal,
-                                 psStatsOptions meanStat, psStatsOptions stdevStat, psRandom *rng)
+bool pmReadoutVarianceRenormalise(const pmReadout *readout, psImageMaskType maskVal,
+                                  int sample, float minValid, float maxValid)
 {
     PM_ASSERT_READOUT_NON_NULL(readout, false);
@@ -371,263 +371,78 @@
     psImage *image = readout->image, *mask = readout->mask, *variance = readout->variance; // Readout parts
 
-    if (!psMemIncrRefCounter(rng)) {
-        rng = psRandomAlloc(PS_RANDOM_TAUS);
-    }
-
-    psStats *varianceStats = psStatsAlloc(meanStat);// Statistics for mean
-    if (!psImageBackground(varianceStats, NULL, variance, mask, maskVal, rng)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to measure mean variance for image");
-        psFree(varianceStats);
-        psFree(rng);
-        return false;
-    }
-    float meanVariance = varianceStats->robustMedian; // Mean variance
-    psFree(varianceStats);
-
-    psStats *imageStats = psStatsAlloc(stdevStat);// Statistics for mean
-    if (!psImageBackground(imageStats, NULL, image, mask, maskVal, rng)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to measure stdev of image");
-        psFree(imageStats);
-        psFree(rng);
-        return false;
-    }
-    float stdevImage = imageStats->robustStdev; // Standard deviation of image
-    psFree(imageStats);
-    psFree(rng);
-
-    float correction = PS_SQR(stdevImage) / meanVariance; // Correction to take variance to what it should be
-    psLogMsg("psModules.camera", PS_LOG_INFO, "Renormalising variance map by %f", correction);
-    psBinaryOp(variance, variance, "*", psScalarAlloc(correction, PS_TYPE_F32));
-
-    return true;
-}
-
-
-bool pmReadoutVarianceRenormPhot(const pmReadout *readout, psImageMaskType maskVal, int num, float width,
-                               psStatsOptions meanStat, psStatsOptions stdevStat, psRandom *rng)
-{
-    PM_ASSERT_READOUT_NON_NULL(readout, false);
-    PM_ASSERT_READOUT_IMAGE(readout, false);
-    PM_ASSERT_READOUT_VARIANCE(readout, false);
-
-    if (!psMemIncrRefCounter(rng)) {
-        rng = psRandomAlloc(PS_RANDOM_TAUS);
-    }
-
-    psImage *image = readout->image, *mask = readout->mask, *variance = readout->variance; // Readout images
-
-    // Measure background
-    psStats *bgStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);// Statistics for background
-    if (!psImageBackground(bgStats, NULL, image, mask, maskVal, rng)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to measure background for image");
-        psFree(bgStats);
-        psFree(rng);
-        return false;
-    }
-    float bgMean = bgStats->robustMedian; // Background level
-    float bgNoise = bgStats->robustStdev; // Background standard deviation
-    psFree(bgStats);
-    psTrace("psModules.camera", 5, "Background is %f +/- %f\n", bgMean, bgNoise);
-
-
-    // Construct kernels for flux measurement
-    // We use N(0,width) and N(0,width/sqrt(2)) kernels, following psphotSignificanceImage.
-    float sigFactor = 4.0 * M_PI * PS_SQR(width); // Factor for conversion from im/wt ratio to significance
-    int size = RENORM_NUM_SIGMA, fullSize = 2 * size + 1; // Half-size and full size of Gaussian
-    psVector *gauss = psVectorAlloc(fullSize, PS_TYPE_F32); // Gaussian for weighting
-    psVector *gauss2 = psVectorAlloc(fullSize, PS_TYPE_F32); // Gaussian squared
-    for (int i = 0, x = -size; i < fullSize; i++, x++) {
-        gauss->data.F32[i] = expf(-0.5 * PS_SQR(x) / PS_SQR(width));
-        gauss2->data.F32[i] = expf(-PS_SQR(x) / PS_SQR(width));
-    }
-
-    // Size of image
-    int numCols = image->numCols, numRows = image->numRows; // Size of images
-    int xSize = numCols - fullSize, ySize = numRows - fullSize; // Size of consideration
-    int xOffset = size, yOffset = size;       // Offset to region of consideration
-
-    // Measure fluxes
-    float peakFlux = RENORM_PEAK * bgNoise;     // Peak flux for fake sources
-    psVector *noise = psVectorAlloc(num, PS_TYPE_F32); // Measurements of the noise
-    psVector *source = psVectorAlloc(num, PS_TYPE_F32); // Measurements of fake sources
-    psVector *guess = psVectorAlloc(num, PS_TYPE_F32); // Guess at significance
-    psVector *photMask = psVectorAlloc(num, PS_TYPE_VECTOR_MASK); // Mask for fluxes
-    for (int i = 0; i < num; i++) {
-        // Coordinates of interest
-        int xPix = psRandomUniform(rng) * xSize + xOffset + 0.5;
-        int yPix = psRandomUniform(rng) * ySize + yOffset + 0.5;
-        psAssert(xPix - size >= 0 && xPix + size < numCols &&
-                 yPix - size >= 0 && yPix + size < numRows,
-                 "Bad pixel position: %d,%d", xPix, yPix);
-
-        // Weighted aperture photometry
-        // This has the same effect as smoothing the image by the window function
-        float sumNoise = 0.0;       // Sum for noise measurement
-        float sumSource = 0.0;      // Sum for source measurement
-        float sumVariance = 0.0;      // Sum for variance measurement
-        float sumGauss = 0.0, sumGauss2 = 0.0; // Sums of Gaussian kernels
-        for (int v = 0, y = yPix - size; v < fullSize; v++, y++) {
-            float xSumNoise = 0.0;  // Sum for noise measurement in x
-            float xSumSource = 0.0; // Sum for source measurement in x
-            float xSumVariance = 0.0; // Sum for variance measurement in x
-            float xSumGauss = 0.0, xSumGauss2 = 0.0; // Sums of Gaussian kernels in x
-            float yGauss = gauss->data.F32[v]; // Value of Gaussian in y
-            for (int u = 0, x = xPix - size; u < fullSize; u++, x++) {
-                if (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) {
+    int numCols = image->numCols, numRows = image->numRows; // Size of image
+    int numPix = numCols * numRows;                         // Number of pixels
+    int num = PS_MAX(sample, numPix);                       // Number we care about
+    psVector *signoise = psVectorAllocEmpty(num, PS_TYPE_F32);   // Signal-to-noise values
+
+    if (num >= numPix) {
+        // We have an image smaller than Nsubset, so just loop over the image pixels
+        int index = 0;                  // Index for vector
+        for (int y = 0; y < numRows; y++) {
+            for (int x = 0; x < numCols; x++) {
+                if ((mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) ||
+                    !isfinite(image->data.F32[y][x]) || !isfinite(variance->data.F32[y][x])) {
                     continue;
                 }
-                float value = image->data.F32[y][x] - bgMean; // Value of image
-                float xGauss = gauss->data.F32[u]; // Value of Gaussian in x
-                float xGauss2 = gauss2->data.F32[u]; // Value of Gaussian^2 in x
-                xSumNoise += value * xGauss;
-                xSumSource += (value + peakFlux * xGauss * yGauss) * xGauss;
-                xSumVariance += variance->data.F32[y][x] * xGauss2;
-                xSumGauss += xGauss;
-                xSumGauss2 += xGauss2;
-            }
-            float yGauss2 = gauss2->data.F32[v]; // Value of Gaussian^2 in y
-            sumNoise += xSumNoise * yGauss;
-            sumSource += xSumSource * yGauss;
-            sumVariance += xSumVariance * yGauss2;
-            sumGauss += xSumGauss * yGauss;
-            sumGauss2 += xSumGauss2 * yGauss2;
-        }
-
-        photMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = ((isfinite(sumNoise) && isfinite(sumSource) &&
-                                                isfinite(sumVariance) && sumGauss > 0 && sumGauss2 > 0) ?
-                                               0 : 0xFF);
-
-        float smoothImageNoise = sumNoise / sumGauss; // Value of smoothed image pixel for noise
-        float smoothImageSource = sumSource / sumGauss; // Value of smoothed image pixel for source
-        float smoothVariance = sumVariance / sumGauss2; // Value of smoothed variance pixel
-
-        noise->data.F32[i] = smoothImageNoise;
-        source->data.F32[i] = smoothImageSource;
-        guess->data.F32[i] = (sumSource > 0) ? sigFactor * PS_SQR(smoothImageSource) / smoothVariance : 0.0;
-        psTrace("psModules.camera", 10, "Flux %d (%d,%d): %f, %f, %f\n",
-                i, xPix, yPix, smoothImageNoise, smoothImageSource, smoothVariance);
-    }
-    psFree(gauss);
-    psFree(gauss2);
-    psFree(rng);
-
-    // Standard deviation of fluxes gives us the real significance
-    psStats *stdevStats = psStatsAlloc(stdevStat); // Statistics
-    if (!psVectorStats(stdevStats, noise, NULL, photMask, 0xFF)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to measure standard deviation of fluxes");
-        psFree(stdevStats);
-        psFree(noise);
-        psFree(source);
-        psFree(guess);
-        psFree(photMask);
+
+                signoise->data.F32[index] = image->data.F32[y][x] / sqrtf(variance->data.F32[y][x]);
+                index++;
+            }
+        }
+        signoise->n = index;
+    } else {
+        psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
+        int index = 0;                  // Index for vector
+        for (long i = 0; i < num; i++) {
+            // Pixel coordinates
+            int pixel = numPix * psRandomUniform(rng);
+            int x = pixel % numCols;
+            int y = pixel / numCols;
+
+            if ((mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) ||
+                !isfinite(image->data.F32[y][x]) || !isfinite(variance->data.F32[y][x])) {
+                continue;
+            }
+
+            signoise->data.F32[index] = image->data.F32[y][x] / sqrtf(variance->data.F32[y][x]);
+            index++;
+        }
+        signoise->n = index;
+        psFree(rng);
+    }
+
+    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_STDEV); // Statistics
+
+    if (!psVectorStats(stats, signoise, NULL, NULL, 0)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to measure statistics on S/N image");
+        psFree(signoise);
         return false;
     }
-    float stdev = psStatsGetValue(stdevStats, stdevStat); // Standard deviation of fluxes
-    psFree(stdevStats);
-    psFree(noise);
-    psTrace("psModules.camera", 5, "Standard deviation of fluxes is %f\n", stdev);
-
-    // Ratio of measured significance to guessed significance
-    psVector *ratio = psVectorAlloc(num, PS_TYPE_F32); // Ratio of measured to guess
-    for (int i = 0; i < num; i++) {
-        float measuredSig = PS_SQR(source->data.F32[i] / stdev); // Measured significance
-        ratio->data.F32[i] = measuredSig / guess->data.F32[i];
-        if (guess->data.F32[i] <= 0.0 || source->data.F32[i] <= 0.0 || !isfinite(ratio->data.F32[i])) {
-            photMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xFF;
-        }
-        psTrace("psModules.camera", 9, "Ratio %d: %f, %f, %f\n",
-                i, guess->data.F32[i], measuredSig, ratio->data.F32[i]);
-    }
-    psFree(source);
-    psFree(guess);
-
-    psStats *meanStats = psStatsAlloc(meanStat | stdevStat); // Statistics
-    if (!psVectorStats(meanStats, ratio, NULL, photMask, 0xFF)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to measure mean ratio");
-        psFree(meanStats);
-        psFree(ratio);
-        psFree(photMask);
-        return false;
-    }
-    float meanRatio = psStatsGetValue(meanStats, meanStat); // Mean ratio
-    psTrace("psModules.camera", 5, "Mean significance ratio is %f +/- %f\n",
-            meanRatio, psStatsGetValue(meanStats, stdevStat));
-    psFree(meanStats);
-    psFree(ratio);
-    psFree(photMask);
-
-    psLogMsg("psModules.camera", PS_LOG_INFO, "Renormalising variance map by %f", meanRatio);
-    psBinaryOp(variance, variance, "*", psScalarAlloc(meanRatio, PS_TYPE_F32));
+    psFree(signoise);
+
+    float correction = stats->robustStdev; // Correction factor
+    psFree(stats);
+    psLogMsg("psModules.camera", PS_LOG_DETAIL, "Variance renormalisation factor is %f", correction);
+
+    // Check valid range of correction factor
+    if ((isfinite(minValid) && correction < minValid) || (isfinite(maxValid) && correction > maxValid)) {
+        psWarning("Variance renormalisation is outside valid range: %f vs %f:%f --- no correction made",
+                  correction, minValid, maxValid);
+        return true;
+    }
+
+    psBinaryOp(variance, variance, "*", psScalarAlloc(PS_SQR(correction), PS_TYPE_F32));
+
+    pmHDU *hdu = pmHDUFromReadout(readout); // HDU for readout
+    if (hdu)  {
+        psString history = NULL;
+        psStringAppend(&history, "Rescaled variance by %6.4f (stdev by %6.4f)",
+                       PS_SQR(correction), correction);
+        psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, NULL, history);
+    }
 
     return true;
 }
 
-
-bool pmReadoutVarianceRenorm(const pmReadout *readout, psImageMaskType maskVal, psStatsOptions meanStat,
-                           psStatsOptions stdevStat, int width, psRandom *rng)
-{
-    PM_ASSERT_READOUT_NON_NULL(readout, false);
-    PM_ASSERT_READOUT_IMAGE(readout, false);
-    PM_ASSERT_READOUT_VARIANCE(readout, false);
-    PS_ASSERT_INT_POSITIVE(width, false);
-
-    if (!psMemIncrRefCounter(rng)) {
-        rng = psRandomAlloc(PS_RANDOM_TAUS);
-    }
-
-    psImage *image = readout->image, *mask = readout->mask, *variance = readout->variance; // Readout images
-    int numCols = image->numCols, numRows = image->numRows; // Size of images
-    int xNum = numCols / width + 1, yNum = numRows / width + 1; // Number of renormalisation regions
-    float xSize = numCols / (float)xNum, ySize = numRows / (float)yNum; // Size of renormalisation regions
-
-    psStats *meanStats = psStatsAlloc(meanStat), *stdevStats = psStatsAlloc(stdevStat); // Statistics
-    psVector *buffer = NULL;
-
-    for (int j = 0; j < yNum; j++) {
-        // Bounds in y
-        int yMin = j * ySize;
-        int yMax = (j + 1) * ySize;
-        for (int i = 0; i < xNum; i++) {
-            // Bounds in x
-            int xMin = i * xSize;
-            int xMax = (i + 1) * xSize;
-
-            psRegion region = psRegionSet(xMin, xMax, yMin, yMax); // Region of interest
-            psImage *subImage = psImageSubset(image, region); // Sub-image of the image pixels
-            psImage *subVariance = psImageSubset(variance, region); // Sub image of the variance pixels
-            psImage *subMask = mask ? psImageSubset(mask, region) : NULL; // Sub-image of the mask pixels
-
-            if (!psImageBackground(stdevStats, &buffer, subImage, subMask, maskVal, rng) ||
-                !psImageBackground(meanStats, &buffer, subVariance, subMask, maskVal, rng)) {
-                // Nothing we can do about it, but don't want to keel over and die, so do our best to flag it.
-                psString regionStr = psRegionToString(region); // String with region
-                psWarning("Unable to measure statistics over %s", regionStr);
-                psFree(regionStr);
-                psErrorClear();
-                psImageInit(subVariance, NAN);
-                if (subMask) {
-                    psImageInit(subMask, maskVal);
-                }
-            } else {
-                double meanVar = psStatsGetValue(meanStats, meanStat); // Mean of variance map
-                double stdev = psStatsGetValue(stdevStats, stdevStat); // Standard deviation of image
-                psTrace("psModules.camera", 3,
-                        "Region [%d:%d,%d:%d] has variance %lf, but mean of variance map is %lf\n",
-                        xMin, xMax, yMin, yMax, PS_SQR(stdev), meanVar);
-                psBinaryOp(subVariance, subVariance, "*", psScalarAlloc(PS_SQR(stdev) / meanVar, PS_TYPE_F32));
-            }
-
-            psFree(subImage);
-            psFree(subVariance);
-            psFree(subMask);
-        }
-    }
-    psFree(meanStats);
-    psFree(stdevStats);
-    psFree(rng);
-    psFree(buffer);
-
-    return true;
-}