Changeset 25766 for branches/eam_branches/20090820/psModules/src/camera

Timestamp:

Oct 2, 2009, 5:10:19 PM (17 years ago)

Author:

eugene

Message:

merge changes from trunk

Location:

branches/eam_branches/20090820

Files:

• . (modified) (1 prop)
• psModules/src/camera/pmFPA.h (modified) (1 diff)
• psModules/src/camera/pmFPAMaskWeight.c (modified) (10 diffs)
• psModules/src/camera/pmFPAMaskWeight.h (modified) (1 diff)
• psModules/src/camera/pmFPAMosaic.c (modified) (13 diffs)
• psModules/src/camera/pmFPARead.c (modified) (1 diff)
• psModules/src/camera/pmHDUGenerate.c (modified) (3 diffs)
• psModules/src/camera/pmReadoutFake.c (modified) (7 diffs, 1 prop)
• psModules/src/camera/pmReadoutFake.h (modified) (1 diff)

branches/eam_branches/20090820/psModules/src/camera/pmFPA.h

@@ -48 +48 @@
     psPlaneTransform *toTPA;  ///< Transformation from focal plane to tangent plane, or NULL
     psProjection *toSky;         ///< Projection from tangent plane to sky, or NULL
+    bool wcsCDkeys;
     // Information
     psMetadata *concepts;               ///< FPA-level concepts

branches/eam_branches/20090820/psModules/src/camera/pmFPAMaskWeight.c

@@ -111 +111 @@
         // 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 +118 @@
         // 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 +140 @@
     // 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 +230 @@
         // 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 +237 @@
         // 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 +248 @@
     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 +262 @@
 
         // 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 +276 @@
     // 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 +362 @@
 
 
-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 +371 @@
     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 covar = sqrtf(psImageCovarianceFactor(readout->covariance)); // Covariance factor
+    float correction = stats->robustStdev / covar; // 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);
+        psFree(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;
-}
 
 

branches/eam_branches/20090820/psModules/src/camera/pmFPAMaskWeight.h

@@ -80 +80 @@
 ///
 /// The variance map is adjusted so that the mean matches the actual pixel variance in the image
-bool pmReadoutVarianceRenormPixels(
-    const pmReadout *readout,           ///< Readout to normalise
-    psImageMaskType maskVal,                 ///< Value to mask
-    psStatsOptions meanStat,            ///< Statistic to measure the mean (of the variance map)
-    psStatsOptions stdevStat,           ///< Statistic to measure the stdev (of the image)
-    psRandom *rng                       ///< Random number generator
-    );
-
-/// Renormalise the variance map to match the actual photometry variance
-///
-/// The variance map is adjusted so that the actual significance of fake sources matches the
-/// guestimated significance
-bool pmReadoutVarianceRenormPhot(
+bool pmReadoutVarianceRenormalise(
     const pmReadout *readout,           ///< Readout to normalise
     psImageMaskType maskVal,            ///< Value to mask
-    int num,                            ///< Number of instances to measure over the image
-    float width,                        ///< Photometry width
-    psStatsOptions meanStat,            ///< Statistic to measure the mean
-    psStatsOptions stdevStat,           ///< Statistic to measure the stdev
-    psRandom *rng                       ///< Random number generator
-    );
-
-/// Renormalise the variance map to match the actual variance
-///
-/// The variance in the image is measured in patches, and the variance map is adjusted so that the mean for
-/// that patch corresponds.
-bool pmReadoutVarianceRenorm(const pmReadout *readout, // Readout to normalise
-                             psImageMaskType maskVal, // Value to mask
-                             psStatsOptions meanStat, // Statistic to measure the mean (of the variance map)
-                             psStatsOptions stdevStat, // Statistic to measure the stdev (of the image)
-                             int width,   // Width of patch (pixels)
-                             psRandom *rng // Random number generator (for sub-sampling images)
+    int sample,                         ///< Sample size
+    float minValid,                     ///< Minimum valid renormalisation, or NAN
+    float maxValid                      ///< Maximum valid renormalisation, or NAN
     );
 

branches/eam_branches/20090820/psModules/src/camera/pmFPAMosaic.c

@@ -740 +740 @@
 static bool chipMosaic(psImage **mosaicImage, // The mosaic image, to be returned
                        psImage **mosaicMask, // The mosaic mask, to be returned
-                       psImage **mosaicVariances, // The mosaic variances, to be returned
+                       psImage **mosaicVariance, // The mosaic variance, to be returned
                        int *xBinChip, int *yBinChip, // The binning in x and y, to be returned
                        const pmChip *chip, // Chip to mosaic
@@ -749 +749 @@
     assert(mosaicImage);
     assert(mosaicMask);
-    assert(mosaicVariances);
+    assert(mosaicVariance);
     assert(xBinChip);
     assert(yBinChip);
@@ -826 +826 @@
     if (allGood) {
         *mosaicImage = imageMosaic(images, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0, BLANK_VALUE);
-        *mosaicVariances = imageMosaic(variances, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0, BLANK_VALUE);
+        *mosaicVariance = imageMosaic(variances, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0, BLANK_VALUE);
         *mosaicMask = imageMosaic(masks, xFlip, yFlip, xBin, yBin, *xBinChip, *yBinChip, x0, y0, blank);
     }
@@ -847 +847 @@
 static bool fpaMosaic(psImage **mosaicImage, // The mosaic image, to be returned
                       psImage **mosaicMask, // The mosaic mask, to be returned
-                      psImage **mosaicVariances, // The mosaic variances, to be returned
+                      psImage **mosaicVariance, // The mosaic variance, to be returned
                       int *xBinFPA, int *yBinFPA, // The binning in x and y, to be returned
                       const pmFPA *fpa,  // FPA to mosaic
@@ -857 +857 @@
     assert(mosaicImage);
     assert(mosaicMask);
-    assert(mosaicVariances);
+    assert(mosaicVariance);
     assert(xBinFPA);
     assert(yBinFPA);
@@ -960 +960 @@
     if (allGood) {
         *mosaicImage = imageMosaic(images, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0, BLANK_VALUE);
-        *mosaicVariances = imageMosaic(variances, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0, BLANK_VALUE);
+        *mosaicVariance = imageMosaic(variances, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0, BLANK_VALUE);
         *mosaicMask = imageMosaic(masks, xFlip, yFlip, xBin, yBin, *xBinFPA, *yBinFPA, x0, y0, blank);
     }
@@ -1025 +1025 @@
     psImage *mosaicImage   = NULL;      // The mosaic image
     psImage *mosaicMask    = NULL;      // The mosaic mask
-    psImage *mosaicVariances = NULL;      // The mosaic variances
+    psImage *mosaicVariance = NULL;      // The mosaic variances
 
     // Find the HDU
@@ -1052 +1052 @@
         }
         if (hdu->variances) {
-            mosaicVariances = psImageSubset(hdu->variances->data[0], bounds);
-            if (!mosaicVariances) {
+            mosaicVariance = psImageSubset(hdu->variances->data[0], bounds);
+            if (!mosaicVariance) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to select variance pixels.\n");
                 return false;
@@ -1061 +1061 @@
         // Case 2 --- we need to mosaic by cut and paste
         psTrace("psModules.camera", 1, "Case 2 mosaicking: cut and paste.\n");
-        if (!chipMosaic(&mosaicImage, &mosaicMask, &mosaicVariances, &xBin, &yBin, source, targetCell, blank)) {
+        if (!chipMosaic(&mosaicImage, &mosaicMask, &mosaicVariance, &xBin, &yBin, source, targetCell, blank)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to mosaic cells.\n");
             return false;
@@ -1069 +1069 @@
     }
     psTrace("psModules.camera", 1, "xBin,yBin: %d,%d\n", xBin, yBin);
+
 
     // Set the concepts for the target cell
@@ -1090 +1091 @@
     target->parent->concepts = psMetadataCopy(target->parent->concepts, source->parent->concepts); // FPA lvl
 
+    // Average the covariances
+    psList *covariances = psListAlloc(NULL); // Input covariance matrices
+    for (int i = 0; i < source->cells->n; i++) {
+        pmCell *cell = source->cells->data[i]; // Cell of interest
+        if (!cell || !cell->data_exists) {
+            continue;
+        }
+        pmReadout *ro = cell->readouts->data[0]; // Readout of interest
+        if (!ro || !ro->covariance) {
+            continue;
+        }
+        psListAdd(covariances, PS_LIST_TAIL, ro->covariance);
+    }
+    psKernel *mosaicCovariance = NULL;  // Covariance for mosaic
+    if (psListLength(covariances) > 0) {
+        psArray *covarArray = psListToArray(covariances); // Array with covariances
+        mosaicCovariance = psImageCovarianceAverage(covarArray);
+        psFree(covarArray);
+    }
+    psFree(covariances);
+
     // Now make a new readout to go in the target cell
     pmReadout *newReadout = pmReadoutAlloc(targetCell); // New readout
     newReadout->image  = mosaicImage;
     newReadout->mask   = mosaicMask;
-    newReadout->variance = mosaicVariances;
+    newReadout->variance = mosaicVariance;
+    newReadout->covariance = mosaicCovariance;
     psFree(newReadout);                 // Drop reference
 
@@ -1334 +1357 @@
     target->concepts = psMetadataCopy(target->concepts, source->concepts);
 
+    // Average the covariances
+    psList *covariances = psListAlloc(NULL); // Input covariance matrices
+    for (int i = 0; i < covariances->n; i++) {
+        pmChip *chip = chips->data[i];  // Chip of interest
+        if (!chip || !chip->data_exists) {
+            continue;
+        }
+        psArray *cells = chip->cells;   // Cells in chip
+        for (long j = 0; j < cells->n; j++) {
+            pmCell *cell = cells->data[i]; // Cell of interest
+            if (!cell || !cell->data_exists) {
+                continue;
+            }
+            pmReadout *ro = cell->readouts->data[0]; // Readout of interest
+            if (!ro || !ro->covariance) {
+                continue;
+            }
+            psListAdd(covariances, PS_LIST_TAIL, ro->covariance);
+        }
+    }
+    psKernel *mosaicCovariances = NULL; // Covariance for mosaic
+    if (psListLength(covariances) > 0) {
+        psArray *covarArray = psListToArray(covariances); // Array with covariances
+        mosaicCovariances = psImageCovarianceAverage(covarArray);
+        psFree(covarArray);
+    }
+    psFree(covariances);
+
     // Now make a new readout to go in the new cell
     pmReadout *newReadout = pmReadoutAlloc(targetCell); // New readout
@@ -1339 +1390 @@
     newReadout->mask   = mosaicMask;
     newReadout->variance = mosaicVariances;
+    newReadout->covariance = mosaicCovariances;
     psFree(newReadout);                 // Drop reference
 

branches/eam_branches/20090820/psModules/src/camera/pmFPARead.c

@@ -387 +387 @@
                 psFree(regionString);
                 psFree(readout);
+                psFree(iter);
                 return false;
             }

branches/eam_branches/20090820/psModules/src/camera/pmHDUGenerate.c

@@ -611 +611 @@
             if (cells->n == 0) {
                 // Nothing to do
+                psFree (cells);
                 return true;
             }
@@ -660 +661 @@
             if (cells->n == 0) {
                 // Nothing to do
+                psFree (cells);
                 return true;
             }
@@ -710 +712 @@
                 return true;
             }
-
             bool status = generateHDU(hdu, cells);
             psFree(cells);

branches/eam_branches/20090820/psModules/src/camera/pmReadoutFake.c

@@ -28 +28 @@
 #define MODEL_TYPE "PS_MODEL_RGAUSS"    // Type of model to use
 #define MAX_AXIS_RATIO 20.0             // Maximum axis ratio for PSF model
-#define SOURCE_MASK (PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_CR_LIMIT) // Mask to apply to input sources
 #define MODEL_MASK (PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE | \
                     PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_LIMITS) // Mask to apply to models
@@ -50 +49 @@
 }
 
-bool pmReadoutFakeFromSources(pmReadout *readout, int numCols, int numRows, const psArray *sources,
-                              const psVector *xOffset, const psVector *yOffset, const pmPSF *psf,
-                              float minFlux, int radius, bool circularise, bool normalisePeak)
+
+bool pmReadoutFakeFromVectors(pmReadout *readout, int numCols, int numRows,
+                              const psVector *x, const psVector *y, const psVector *mag,
+                              const psVector *xOffset, const psVector *yOffset,
+                              const pmPSF *psf, float minFlux, int radius,
+                              bool circularise, bool normalisePeak)
 {
     PS_ASSERT_PTR_NON_NULL(readout, false);
     PS_ASSERT_INT_LARGER_THAN(numCols, 0, false);
     PS_ASSERT_INT_LARGER_THAN(numRows, 0, false);
-    PS_ASSERT_ARRAY_NON_NULL(sources, false);
-
+    PS_ASSERT_VECTOR_NON_NULL(x, false);
+    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, false);
+    PS_ASSERT_VECTOR_NON_NULL(y, false);
+    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(y, x, false);
+    PS_ASSERT_VECTOR_NON_NULL(mag, false);
+    PS_ASSERT_VECTOR_TYPE(mag, PS_TYPE_F32, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(mag, x, false);
+    long numSources = x->n;              // Number of sources
     if (xOffset || yOffset) {
         PS_ASSERT_VECTOR_NON_NULL(xOffset, false);
@@ -64 +73 @@
         PS_ASSERT_VECTORS_SIZE_EQUAL(xOffset, yOffset, false);
         PS_ASSERT_VECTOR_TYPE(xOffset, PS_TYPE_S32, false);
-        PS_ASSERT_VECTOR_TYPE_EQUAL(xOffset, yOffset, false);
-        if (xOffset->n != sources->n) {
+        PS_ASSERT_VECTOR_TYPE(yOffset, PS_TYPE_S32, false);
+        if (xOffset->n != numSources) {
             psError(PS_ERR_BAD_PARAMETER_SIZE, true,
                     "Number of offset vectors (%ld) and sources (%ld) doesn't match",
-                    xOffset->n, sources->n);
+                    xOffset->n, numSources);
             return false;
         }
     }
     PS_ASSERT_PTR_NON_NULL(psf, false);
-    if (radius > 0 && isfinite(minFlux) && minFlux > 0.0) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Cannot define both minimum flux and fixed radius.");
-        return false;
-    }
 
     readout->image = psImageRecycle(readout->image, numCols, numRows, PS_TYPE_F32);
     psImageInit(readout->image, 0);
 
-    int numSources = sources->n;          // Number of stars
-    for (int i = 0; i < numSources; i++) {
-        pmSource *source = sources->data[i]; // Source of interest
-        if (!source) {
-            continue;
-        }
-        if (source->mode & SOURCE_MASK) {
-            continue;
-        }
-        if (!isfinite(source->psfMag)) {
-            continue;
-        }
-        float x, y;                     // Coordinates of source
-        if (source->modelPSF) {
-            x = source->modelPSF->params->data.F32[PM_PAR_XPOS];
-            y = source->modelPSF->params->data.F32[PM_PAR_YPOS];
-        } else {
-            x = source->peak->xf;
-            y = source->peak->yf;
-        }
-
-        float flux = powf(10.0, -0.4 * source->psfMag); // Flux of source
+    for (long i = 0; i < numSources; i++) {
+        float flux = powf(10.0, -0.4 * mag->data.F32[i]); // Flux of source
+        float xSrc = x->data.F32[i], ySrc = y->data.F32[i]; // Coordinates of source
 
         if (normalisePeak) {
             // Normalise flux
-            pmModel *normModel = pmModelFromPSFforXY(psf, x, y, 1.0); // Model for normalisation
+            pmModel *normModel = pmModelFromPSFforXY(psf, xSrc, ySrc, 1.0); // Model for normalisation
             if (!normModel || (normModel->flags & MODEL_MASK)) {
                 psFree(normModel);
                 continue;
             }
+            // check that all params are valid:
+            bool validParams = true;
+            for (int n = 0; validParams && (n < normModel->params->n); n++) {
+                if (n == PM_PAR_SKY) continue;
+                if (n == PM_PAR_I0) continue;
+                if (n == PM_PAR_XPOS) continue;
+                if (n == PM_PAR_YPOS) continue;
+                if (!isfinite(normModel->params->data.F32[n])) validParams = false;
+            }
+            if (!validParams) {
+                psFree(normModel);
+                continue;
+            }           
             if (circularise && !circulariseModel(normModel)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to circularise PSF model.");
@@ -121 +120 @@
         }
 
-        pmModel *fakeModel = pmModelFromPSFforXY(psf, x, y, flux);
+        pmModel *fakeModel = pmModelFromPSFforXY(psf, xSrc, ySrc, flux);
         if (!fakeModel || (fakeModel->flags & MODEL_MASK)) {
             psFree(fakeModel);
             continue;
         }
+        // check that all params are valid:
+        bool validParams = true;
+        for (int n = 0; validParams && (n < fakeModel->params->n); n++) {
+            if (n == PM_PAR_SKY) continue;
+            if (n == PM_PAR_I0) continue;
+            if (n == PM_PAR_XPOS) continue;
+            if (n == PM_PAR_YPOS) continue;
+            if (!isfinite(fakeModel->params->data.F32[n])) validParams = false;
+        }
+        if (!validParams) {
+            psFree(fakeModel);
+            continue;
+        }               
         if (circularise && !circulariseModel(fakeModel)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to circularise PSF model.");
@@ -137 +149 @@
 
         pmSource *fakeSource = pmSourceAlloc(); // Fake source to generate
-        fakeSource->peak = pmPeakAlloc(x, y, fakeModel->params->data.F32[PM_PAR_I0], PM_PEAK_LONE);
-        float fakeRadius = radius > 0 ? radius :
-            PS_MAX(1.0, fakeModel->modelRadius(fakeModel->params, minFlux)); // Radius of fake source
+        fakeSource->peak = pmPeakAlloc(xSrc, ySrc, fakeModel->params->data.F32[PM_PAR_I0], PM_PEAK_LONE);
+        float fakeRadius = 1.0;         // Radius of fake source
+        if (isfinite(minFlux)) {
+            fakeRadius = PS_MAX(fakeRadius, fakeModel->modelRadius(fakeModel->params, minFlux));
+        }
+        if (radius > 0) {
+            fakeRadius = PS_MAX(fakeRadius, radius);
+        }
 
         if (xOffset) {
-            if (!pmSourceDefinePixels(fakeSource, readout, x + xOffset->data.S32[i],
-                                      y + yOffset->data.S32[i], fakeRadius)) {
+            if (!pmSourceDefinePixels(fakeSource, readout, xSrc + xOffset->data.S32[i],
+                                      ySrc + yOffset->data.S32[i], fakeRadius)) {
                 psErrorClear();
                 continue;
@@ -153 +170 @@
             }
         } else {
-            if (!pmSourceDefinePixels(fakeSource, readout, x, y, fakeRadius)) {
+            if (!pmSourceDefinePixels(fakeSource, readout, xSrc, ySrc, fakeRadius)) {
                 psErrorClear();
                 continue;
@@ -168 +185 @@
     return true;
 }
+
+
+bool pmReadoutFakeFromSources(pmReadout *readout, int numCols, int numRows, const psArray *sources,
+                              pmSourceMode sourceMask, const psVector *xOffset, const psVector *yOffset,
+                              const pmPSF *psf, float minFlux, int radius,
+                              bool circularise, bool normalisePeak)
+{
+    PS_ASSERT_ARRAY_NON_NULL(sources, false);
+
+    int numSources = sources->n;          // Number of stars
+    psVector *x = psVectorAllocEmpty(numSources, PS_TYPE_F32);
+    psVector *y = psVectorAllocEmpty(numSources, PS_TYPE_F32);
+    psVector *mag = psVectorAllocEmpty(numSources, PS_TYPE_F32);
+
+    int numGood = 0;                    // Number of good sources
+    for (int i = 0; i < numSources; i++) {
+        pmSource *source = sources->data[i]; // Source of interest
+        if (!source) {
+            continue;
+        }
+        if (source->mode & sourceMask) {
+            continue;
+        }
+        if (!isfinite(source->psfMag)) {
+            continue;
+        }
+        float xSrc, ySrc;                     // Coordinates of source
+        if (source->modelPSF) {
+            xSrc = source->modelPSF->params->data.F32[PM_PAR_XPOS];
+            ySrc = source->modelPSF->params->data.F32[PM_PAR_YPOS];
+        } else {
+            xSrc = source->peak->xf;
+            ySrc = source->peak->yf;
+        }
+
+        x->data.F32[numGood] = xSrc;
+        y->data.F32[numGood] = ySrc;
+        mag->data.F32[numGood] = source->psfMag;
+        numGood++;
+    }
+    x->n = numGood;
+    y->n = numGood;
+    mag->n = numGood;
+
+    bool status = pmReadoutFakeFromVectors(readout, numCols, numRows, x, y, mag, xOffset, yOffset, psf,
+                                           minFlux, radius, circularise, normalisePeak);
+    psFree(x);
+    psFree(y);
+    psFree(mag);
+
+    return status;
+}

branches/eam_branches/20090820/psModules/src/camera/pmReadoutFake.h

@@ -11 +11 @@
 #include <pmTrend2D.h>
 #include <pmPSF.h>
+#include <pmSourceMasks.h>
+
+/// Generate a fake readout from vectors
+bool pmReadoutFakeFromVectors(pmReadout *readout, ///< Output readout
+                              int numCols, int numRows, ///< Dimension of image
+                              const psVector *x, const psVector *y, ///< Source coordinates
+                              const psVector *mag, ///< Source magnitudes
+                              const psVector *xOffset, ///< x offsets for sources (source -> img), or NULL
+                              const psVector *yOffset, ///< y offsets for sources (source -> img), or NULL
+                              const pmPSF *psf, ///< PSF for sources
+                              float minFlux, ///< Minimum flux to bother about; for setting source radius
+                              int radius, ///< Fixed radius for sources
+                              bool circularise, ///< Circularise PSF model?
+                              bool normalisePeak ///< Normalise the peak value?
+    );
 
 /// Generate a fake readout from an array of sources
-bool pmReadoutFakeFromSources(pmReadout *readout, ///< Output readout, or NULL
+bool pmReadoutFakeFromSources(pmReadout *readout, ///< Output readout
                               int numCols, int numRows, ///< Dimension of image
                               const psArray *sources, ///< Array of pmSource
+                              pmSourceMode sourceMask, ///< Mask for sources
                               const psVector *xOffset, ///< x offsets for sources (source -> img), or NULL
                               const psVector *yOffset, ///< y offsets for sources (source -> img), or NULL