Changeset 21363 for trunk/psModules/src/camera/pmFPAMaskWeight.c

Timestamp:

Feb 5, 2009, 4:31:25 PM (17 years ago)

Author:

Paul Price

Message:

Merging pap_branch_20090128. Resolved a small number of conflicts. Compiles, but not tested in detail.

File:

• trunk/psModules/src/camera/pmFPAMaskWeight.c (modified) (37 diffs)

trunk/psModules/src/camera/pmFPAMaskWeight.c

@@ -49 +49 @@
 }
 
-// Create the parent weight images that reside in the HDU
-static void createParentWeights(pmHDU *hdu // The HDU for which to create
+// Create the parent variance images that reside in the HDU
+static void createParentVariances(pmHDU *hdu // The HDU for which to create
                                )
 {
@@ -58 +58 @@
     // Generate the parent mask images
     psArray *images = hdu->images;      // Array of images
-    psArray *weights = hdu->weights;    // Array of weight images
-    if (!weights) {
-        weights = psArrayAlloc(images->n);
-        hdu->weights = weights;
+    psArray *variances = hdu->variances;    // Array of variance images
+    if (!variances) {
+        variances = psArrayAlloc(images->n);
+        hdu->variances = variances;
     }
 
     for (long i = 0; i < images->n; i++) {
         psImage *image = images->data[i]; // The image for this readout
-        if (!image || weights->data[i]) {
+        if (!image || variances->data[i]) {
             continue;
         }
-        weights->data[i] = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32);
+        variances->data[i] = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32);
     }
 
@@ -199 +199 @@
 }
 
-bool pmReadoutSetWeight(pmReadout *readout, bool poisson)
+bool pmReadoutSetVariance(pmReadout *readout, bool poisson)
 {
     PS_ASSERT_PTR_NON_NULL(readout, false);
@@ -209 +209 @@
     float gain = psMetadataLookupF32(&mdok, cell->concepts, "CELL.GAIN"); // Cell gain
     if (!mdok || isnan(gain)) {
-        psError(PS_ERR_IO, true, "CELL.GAIN is not set --- unable to set weight.\n");
+        psError(PS_ERR_IO, true, "CELL.GAIN is not set --- unable to set variance.\n");
         return false;
     }
     float readnoise = psMetadataLookupF32(&mdok, cell->concepts, "CELL.READNOISE"); // Cell read noise
     if (!mdok || isnan(readnoise)) {
-        psError(PS_ERR_IO, true, "CELL.READNOISE is not set --- unable to set weight.\n");
+        psError(PS_ERR_IO, true, "CELL.READNOISE is not set --- unable to set variance.\n");
         return false;
     }
@@ -223 +223 @@
 
     if (poisson) {
-        // Set weight image to the variance in ADU = f/g + rn^2
+        // Set variance image to the variance in ADU = f/g + rn^2
         psImage *image = readout->image;    // The image pixels
-        readout->weight = (psImage*)psBinaryOp(readout->weight, image, "/", psScalarAlloc(gain, PS_TYPE_F32));
-
-        // a negative weight is non-sensical. if the image value drops below 1, the weight must be 1.
-        readout->weight = (psImage*)psUnaryOp(readout->weight, readout->weight, "abs");
-        readout->weight = (psImage*)psBinaryOp(readout->weight, readout->weight, "max",
+        readout->variance = (psImage*)psBinaryOp(readout->variance, image, "/", psScalarAlloc(gain, PS_TYPE_F32));
+
+        // a negative variance is non-sensical. if the image value drops below 1, the variance must be 1.
+        readout->variance = (psImage*)psUnaryOp(readout->variance, readout->variance, "abs");
+        readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "max",
                                                psScalarAlloc(1, PS_TYPE_F32));
     } else {
         // Just use the read noise
-        if (!readout->weight) {
-            readout->weight = psImageAlloc(readout->image->numCols, readout->image->numRows, PS_TYPE_F32);
-        }
-        psImageInit(readout->weight, 0.0);
-    }
-
-    readout->weight = (psImage*)psBinaryOp(readout->weight, readout->weight, "+",
+        if (!readout->variance) {
+            readout->variance = psImageAlloc(readout->image->numCols, readout->image->numRows, PS_TYPE_F32);
+        }
+        psImageInit(readout->variance, 0.0);
+    }
+
+    readout->variance = (psImage*)psBinaryOp(readout->variance, readout->variance, "+",
                                            psScalarAlloc(readnoise*readnoise/gain/gain, PS_TYPE_F32));
 
@@ -245 +245 @@
 }
 
-// this function creates the weight pixels, or uses the existing weight pixels.  it will set
-// the noise pixel values only if the weight image is not supplied
-bool pmReadoutGenerateWeight(pmReadout *readout, bool poisson)
+// this function creates the variance pixels, or uses the existing variance pixels.  it will set
+// the noise pixel values only if the variance image is not supplied
+bool pmReadoutGenerateVariance(pmReadout *readout, bool poisson)
 {
     PS_ASSERT_PTR_NON_NULL(readout, false);
@@ -254 +254 @@
     bool mdok = true;                   // Status of MD lookup
 
-    // Create the weight image if required
-    if (readout->weight)
+    // Create the variance image if required
+    if (readout->variance)
         return true;
 
     psRegion *trimsec = psMetadataLookupPtr(&mdok, cell->concepts, "CELL.TRIMSEC"); // Trim section
     if (!mdok || psRegionIsNaN(*trimsec)) {
-        psError(PS_ERR_IO, true, "CELL.TRIMSEC is not set --- unable to set weight.\n");
+        psError(PS_ERR_IO, true, "CELL.TRIMSEC is not set --- unable to set variance.\n");
         return false;
     }
@@ -271 +271 @@
     }
 
-    createParentWeights(hdu);
+    createParentVariances(hdu);
 
     // Need to identify which readout we're working with....
@@ -280 +280 @@
     }
 
-    psImage *weight = psImageSubset(hdu->weights->data[index], *trimsec); // The weight pixels
-    if (!weight) {
+    psImage *variance = psImageSubset(hdu->variances->data[index], *trimsec); // The variance pixels
+    if (!variance) {
         psString trimsecString = psRegionToString(*trimsec);
-        psError(PS_ERR_UNKNOWN, false, "Unable to set weight from HDU with trimsec: %s.\n",
+        psError(PS_ERR_UNKNOWN, false, "Unable to set variance from HDU with trimsec: %s.\n",
                 trimsecString);
         psFree(trimsecString);
         return false;
     }
-    psImageInit(weight, 0);
-    readout->weight = weight;
-
-    return pmReadoutSetWeight(readout, poisson);
-}
-
-bool pmReadoutGenerateMaskWeight(pmReadout *readout, psImageMaskType satMask, psImageMaskType badMask, bool poisson)
+    psImageInit(variance, 0);
+    readout->variance = variance;
+
+    return pmReadoutSetVariance(readout, poisson);
+}
+
+bool pmReadoutGenerateMaskVariance(pmReadout *readout, psImageMaskType satMask, psImageMaskType badMask, bool poisson)
 {
     PS_ASSERT_PTR_NON_NULL(readout, false);
@@ -301 +301 @@
 
     success &= pmReadoutGenerateMask(readout, satMask, badMask);
-    success &= pmReadoutGenerateWeight(readout, poisson);
+    success &= pmReadoutGenerateVariance(readout, poisson);
 
     return success;
 }
 
-bool pmCellGenerateMaskWeight(pmCell *cell, psImageMaskType satMask, psImageMaskType badMask, bool poisson)
+bool pmCellGenerateMaskVariance(pmCell *cell, psImageMaskType satMask, psImageMaskType badMask, bool poisson)
 {
     PS_ASSERT_PTR_NON_NULL(cell, false);
@@ -314 +314 @@
     for (int i = 0; i < readouts->n; i++) {
         pmReadout *readout = readouts->data[i]; // The readout
-        success &= pmReadoutGenerateMaskWeight(readout, poisson, satMask, badMask);
+        success &= pmReadoutGenerateMaskVariance(readout, poisson, satMask, badMask);
     }
 
@@ -321 +321 @@
 
 
-bool pmReadoutWeightRenormPixels(const pmReadout *readout, psImageMaskType maskVal,
+bool pmReadoutVarianceRenormPixels(const pmReadout *readout, psImageMaskType maskVal,
                                  psStatsOptions meanStat, psStatsOptions stdevStat, psRandom *rng)
 {
     PM_ASSERT_READOUT_NON_NULL(readout, false);
     PM_ASSERT_READOUT_IMAGE(readout, false);
-    PM_ASSERT_READOUT_WEIGHT(readout, false);
-
-    psImage *image = readout->image, *mask = readout->mask, *weight = readout->weight; // Readout parts
+    PM_ASSERT_READOUT_VARIANCE(readout, false);
+
+    psImage *image = readout->image, *mask = readout->mask, *variance = readout->variance; // Readout parts
 
     if (!psMemIncrRefCounter(rng)) {
@@ -334 +334 @@
     }
 
-    psStats *weightStats = psStatsAlloc(meanStat);// Statistics for mean
-    if (!psImageBackground(weightStats, NULL, weight, mask, maskVal, rng)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to measure mean weight for image");
-        psFree(weightStats);
+    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 = weightStats->robustMedian; // Mean variance
-    psFree(weightStats);
+    float meanVariance = varianceStats->robustMedian; // Mean variance
+    psFree(varianceStats);
 
     psStats *imageStats = psStatsAlloc(stdevStat);// Statistics for mean
@@ -356 +356 @@
 
     float correction = PS_SQR(stdevImage) / meanVariance; // Correction to take variance to what it should be
-    psLogMsg("psModules.camera", PS_LOG_INFO, "Renormalising weight map by %f", correction);
-    psBinaryOp(weight, weight, "*", psScalarAlloc(correction, PS_TYPE_F32));
+    psLogMsg("psModules.camera", PS_LOG_INFO, "Renormalising variance map by %f", correction);
+    psBinaryOp(variance, variance, "*", psScalarAlloc(correction, PS_TYPE_F32));
 
     return true;
@@ -363 +363 @@
 
 
-bool pmReadoutWeightRenormPhot(const pmReadout *readout, psImageMaskType maskVal, int num, float width,
+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_WEIGHT(readout, false);
+    PM_ASSERT_READOUT_VARIANCE(readout, false);
 
     if (!psMemIncrRefCounter(rng)) {
@@ -374 +374 @@
     }
 
-    psImage *image = readout->image, *mask = readout->mask, *weight = readout->weight; // Readout images
+    psImage *image = readout->image, *mask = readout->mask, *variance = readout->variance; // Readout images
 
     // Measure background
@@ -424 +424 @@
         float sumNoise = 0.0;       // Sum for noise measurement
         float sumSource = 0.0;      // Sum for source measurement
-        float sumWeight = 0.0;      // Sum for weight 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 xSumWeight = 0.0; // Sum for weight 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
@@ -441 +441 @@
                 xSumNoise += value * xGauss;
                 xSumSource += (value + peakFlux * xGauss * yGauss) * xGauss;
-                xSumWeight += weight->data.F32[y][x] * xGauss2;
+                xSumVariance += variance->data.F32[y][x] * xGauss2;
                 xSumGauss += xGauss;
                 xSumGauss2 += xGauss2;
@@ -448 +448 @@
             sumNoise += xSumNoise * yGauss;
             sumSource += xSumSource * yGauss;
-            sumWeight += xSumWeight * yGauss2;
+            sumVariance += xSumVariance * yGauss2;
             sumGauss += xSumGauss * yGauss;
             sumGauss2 += xSumGauss2 * yGauss2;
@@ -454 +454 @@
 
         photMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = ((isfinite(sumNoise) && isfinite(sumSource) &&
-                                                isfinite(sumWeight) && sumGauss > 0 && sumGauss2 > 0) ?
+                                                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 smoothWeight = sumWeight / sumGauss2; // Value of smoothed weight pixel
+        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) / smoothWeight : 0.0;
+        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, smoothWeight);
+                i, xPix, yPix, smoothImageNoise, smoothImageSource, smoothVariance);
     }
     psFree(gauss);
@@ -516 +516 @@
     psFree(photMask);
 
-    psLogMsg("psModules.camera", PS_LOG_INFO, "Renormalising weight map by %f", meanRatio);
-    psBinaryOp(weight, weight, "*", psScalarAlloc(meanRatio, PS_TYPE_F32));
+    psLogMsg("psModules.camera", PS_LOG_INFO, "Renormalising variance map by %f", meanRatio);
+    psBinaryOp(variance, variance, "*", psScalarAlloc(meanRatio, PS_TYPE_F32));
 
     return true;
@@ -523 +523 @@
 
 
-bool pmReadoutWeightRenorm(const pmReadout *readout, psImageMaskType maskVal, psStatsOptions meanStat,
+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_WEIGHT(readout, false);
+    PM_ASSERT_READOUT_VARIANCE(readout, false);
     PS_ASSERT_INT_POSITIVE(width, false);
 
@@ -535 +535 @@
     }
 
-    psImage *image = readout->image, *mask = readout->mask, *weight = readout->weight; // Readout images
+    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
@@ -554 +554 @@
             psRegion region = psRegionSet(xMin, xMax, yMin, yMax); // Region of interest
             psImage *subImage = psImageSubset(image, region); // Sub-image of the image pixels
-            psImage *subWeight = psImageSubset(weight, region); // Sub image of the weight 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, subWeight, 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
@@ -564 +564 @@
                 psFree(regionStr);
                 psErrorClear();
-                psImageInit(subWeight, NAN);
+                psImageInit(subVariance, NAN);
                 if (subMask) {
                     psImageInit(subMask, maskVal);
@@ -574 +574 @@
                         "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(subWeight, subWeight, "*", psScalarAlloc(PS_SQR(stdev) / meanVar, PS_TYPE_F32));
+                psBinaryOp(subVariance, subVariance, "*", psScalarAlloc(PS_SQR(stdev) / meanVar, PS_TYPE_F32));
             }
 
             psFree(subImage);
-            psFree(subWeight);
+            psFree(subVariance);
             psFree(subMask);
         }
@@ -597 +597 @@
 
     psImage *image = readout->image;    // Readout's image
-    psImage *weight = readout->weight;  // Readout's weight
+    psImage *variance = readout->variance;  // Readout's variance
     int numCols = image->numCols, numRows = image->numRows; // Size of image
 
@@ -607 +607 @@
     for (int y = 0; y < numRows; y++) {
         for (int x = 0; x < numCols; x++) {
-            if (!isfinite(image->data.F32[y][x]) || (weight && !isfinite(weight->data.F32[y][x]))) {
+            if (!isfinite(image->data.F32[y][x]) || (variance && !isfinite(variance->data.F32[y][x]))) {
                 mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= maskVal;
             }
@@ -627 +627 @@
     psImageMaskType **maskData = readout->mask->data.PS_TYPE_IMAGE_MASK_DATA; // Dereference mask
     psF32 **imageData = readout->image->data.F32;// Dereference image
-    psF32 **weightData = readout->weight ? readout->weight->data.F32 : NULL; // Dereference weight map
+    psF32 **varianceData = readout->variance ? readout->variance->data.F32 : NULL; // Dereference variance map
 
     for (int y = 0; y < numRows; y++) {
@@ -633 +633 @@
             if (maskData[y][x] & maskVal) {
                 imageData[y][x] = NAN;
-                if (weightData) {
-                    weightData[y][x] = NAN;
+                if (varianceData) {
+                    varianceData[y][x] = NAN;
                 }
             }
@@ -656 +656 @@
     psImage *image = readout->image;    // Image
     psImage *mask = readout->mask;      // Mask
-    psImage *weight = readout->weight;  // Weight map
-
-    psImageInterpolation *interp = psImageInterpolationAlloc(mode, image, weight, mask, maskVal,
+    psImage *variance = readout->variance;  // Variance map
+
+    psImageInterpolation *interp = psImageInterpolationAlloc(mode, image, variance, mask, maskVal,
                                                              NAN, NAN, maskBad, maskPoor, poorFrac, 0);
     interp->shifting = false;           // Turn off "exact shifts" so we get proper interpolation
@@ -666 +666 @@
     psPixels *pixels = psPixelsAllocEmpty(PIXELS_BUFFER); // Pixels that have been interpolated
     psVector *imagePix = psVectorAllocEmpty(PIXELS_BUFFER, PS_TYPE_F32); // Corresponding values for image
-    psVector *weightPix = psVectorAllocEmpty(PIXELS_BUFFER, PS_TYPE_F32); // Corresponding values for weight
+    psVector *variancePix = psVectorAllocEmpty(PIXELS_BUFFER, PS_TYPE_F32); // Corresponding values for variance
     psVector *maskPix = psVectorAllocEmpty(PIXELS_BUFFER, PS_TYPE_IMAGE_MASK); // Corresponding values for mask
     // NOTE: maskPix carries the actual image mask values -- do NOT use
@@ -675 +675 @@
         for (int x = 0; x < numCols; x++) {
             if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) {
-                double imageValue, weightValue; // Image and weight value from interpolation
+                double imageValue, varianceValue; // Image and variance value from interpolation
                 psImageMaskType maskValue = 0; // Mask value from interpolation
-                psImageInterpolateStatus status = psImageInterpolate(&imageValue, &weightValue, &maskValue, x, y, interp);
+                psImageInterpolateStatus status = psImageInterpolate(&imageValue, &varianceValue, &maskValue, x, y, interp);
                 if (status == PS_INTERPOLATE_STATUS_ERROR || status == PS_INTERPOLATE_STATUS_OFF) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to interpolate readout at %d,%d", x, y);
@@ -683 +683 @@
                     psFree(pixels);
                     psFree(imagePix);
-                    psFree(weightPix);
+                    psFree(variancePix);
                     psFree(maskPix);
                     return false;
@@ -694 +694 @@
                 pixels = psPixelsAdd(pixels, PIXELS_BUFFER, x, y);
                 imagePix = psVectorExtend(imagePix, PIXELS_BUFFER, 1);
-                weightPix = psVectorExtend(weightPix, PIXELS_BUFFER, 1);
+                variancePix = psVectorExtend(variancePix, PIXELS_BUFFER, 1);
                 maskPix = psVectorExtend(maskPix, PIXELS_BUFFER, 1);
                 imagePix->data.F32[numBad] = imageValue;
-                weightPix->data.F32[numBad] = weightValue;
+                variancePix->data.F32[numBad] = varianceValue;
                 maskPix->data.PS_TYPE_IMAGE_MASK_DATA[numBad] = maskValue;
                 numBad++;
@@ -709 +709 @@
         int x = pixels->data[i].x, y = pixels->data[i].y; // Coordinates of pixel
         image->data.F32[y][x] = imagePix->data.F32[i];
-        weight->data.F32[y][x] = weightPix->data.F32[i];
+        variance->data.F32[y][x] = variancePix->data.F32[i];
         mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = maskPix->data.PS_TYPE_IMAGE_MASK_DATA[i];
     }
@@ -715 +715 @@
     psFree(pixels);
     psFree(imagePix);
-    psFree(weightPix);
+    psFree(variancePix);
     psFree(maskPix);