Changeset 21211 for branches/pap_branch_20090128/psModules/src/imcombine/pmSubtraction.c

Timestamp:

Jan 28, 2009, 2:33:51 PM (17 years ago)

Author:

Paul Price

Message:

Changing pmReadout.weight to variance. Adding pmReadout.covariance which will carry around a covariance pseudo-matrix, allowing us to calculate the pixel-to-pixel variance. pmReadout.covariance now exists and is set, but no mechanism yet to read/write, or use it.

File:

• branches/pap_branch_20090128/psModules/src/imcombine/pmSubtraction.c (modified) (17 diffs)

branches/pap_branch_20090128/psModules/src/imcombine/pmSubtraction.c

@@ -60 +60 @@
 
     // Normalise so that the sum of the variance kernel is the square of the sum of the normal kernel
-    // This is required to keep the relative scaling between the image and the weight map
+    // This is required to keep the relative scaling between the image and the variance map
     psBinaryOp(out->image, out->image, "*", psScalarAlloc(PS_SQR(sumNormal) / sumVariance, PS_TYPE_F32));
 
@@ -287 +287 @@
 
 // Convolve an image using FFT
-static void convolveWeightFFT(psImage *target,// Place the result in here
-                              psImage *weight, // Weight map to convolve
+static void convolveVarianceFFT(psImage *target,// Place the result in here
+                              psImage *variance, // Variance map to convolve
                               psImage *sys, // Systematic error image
                               psImage *mask, // Mask image
@@ -302 +302 @@
     bool threaded = pmSubtractionThreaded(); // Are we running threaded?
 
-    psImage *subWeight = convolveSubsetAlloc(weight, border, threaded); // Weight map
+    psImage *subVariance = convolveSubsetAlloc(variance, border, threaded); // Variance map
     psImage *subSys = convolveSubsetAlloc(sys, border, threaded); // Systematic error image
     psImage *subMask = convolveSubsetAlloc(mask, border, threaded); // Mask
 
     // XXX Can trim this a little by combining the convolution: only have to take the FFT of the kernel once
-    psImage *convWeight = psImageConvolveFFT(NULL, subWeight, subMask, maskVal, kernel); // Convolved weight
+    psImage *convVariance = psImageConvolveFFT(NULL, subVariance, subMask, maskVal, kernel); // Convolved variance
     psImage *convSys = subSys ? psImageConvolveFFT(NULL, subSys, subMask, maskVal, kernel) : NULL; // Conv sys
 
-    convolveSubsetFree(weight, subWeight, threaded);
+    convolveSubsetFree(variance, subVariance, threaded);
     convolveSubsetFree(sys, subSys, threaded);
     convolveSubsetFree(mask, subMask, threaded);
@@ -319 +319 @@
         for (int yTarget = yMin, ySource = size; yTarget < yMax; yTarget++, ySource++) {
             for (int xTarget = xMin, xSource = size; xTarget < xMax; xTarget++, xSource++) {
-                target->data.F32[yTarget][xTarget] = convWeight->data.F32[ySource][xSource] +
+                target->data.F32[yTarget][xTarget] = convVariance->data.F32[ySource][xSource] +
                     convSys->data.F32[ySource][xSource];
             }
@@ -326 +326 @@
         int numBytes = (xMax - xMin) * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes to copy
         for (int yTarget = yMin, ySource = size; yTarget < yMax; yTarget++, ySource++) {
-            memcpy(&target->data.F32[yTarget][xMin], &convWeight->data.F32[ySource][size], numBytes);
-        }
-    }
-
-    psFree(convWeight);
+            memcpy(&target->data.F32[yTarget][xMin], &convVariance->data.F32[ySource][size], numBytes);
+        }
+    }
+
+    psFree(convVariance);
     psFree(convSys);
 
@@ -361 +361 @@
 // Convolve a region of an image
 static inline void convolveRegion(psImage *convImage, // Convolved image (output)
-                                  psImage *convWeight, // Convolved weight map (output), or NULL
+                                  psImage *convVariance, // Convolved variance map (output), or NULL
                                   psImage *convMask, // Convolve mask (output), or NULL
                                   psKernel **kernelImage, // Convolution kernel for the image
-                                  psKernel **kernelWeight, // Convolution kernel for the weight map, or NULL
+                                  psKernel **kernelVariance, // Convolution kernel for the variance map, or NULL
                                   psImage *image, // Image to convolve
-                                  psImage *weight, // Weight map to convolve, or NULL
+                                  psImage *variance, // Variance map to convolve, or NULL
                                   psImage *sys, // Systematic error image, or NULL
                                   psImage *subMask, // Subtraction mask
@@ -381 +381 @@
 {
     *kernelImage = solvedKernel(*kernelImage, kernels, polyValues, wantDual);
-    if (weight || subMask) {
-        *kernelWeight = varianceKernel(*kernelWeight, *kernelImage);
+    if (variance || subMask) {
+        *kernelVariance = varianceKernel(*kernelVariance, *kernelImage);
     }
 
@@ -398 +398 @@
     }
 
-    // Convolve the image and weight
+    // Convolve the image and variance
     if (useFFT) {
         // Use Fast Fourier Transform to do the convolution
         // This provides a big speed-up for large kernels
         convolveFFT(convImage, image, subMask, subBad, *kernelImage, region, background, kernels->size);
-        if (weight) {
-            convolveWeightFFT(convWeight, weight, sys, subMask, subBad, *kernelWeight, region, kernels->size);
+        if (variance) {
+            convolveVarianceFFT(convVariance, variance, sys, subMask, subBad, *kernelVariance, region, kernels->size);
         }
     } else {
         // XXX Direct convolution doesn't account for bad pixels yet
         convolveDirect(convImage, image, *kernelImage, region, background, kernels->size);
-        if (weight) {
-            convolveDirect(convWeight, weight, *kernelWeight, region, 0.0, kernels->size);
+        if (variance) {
+            convolveDirect(convVariance, variance, *kernelVariance, region, 0.0, kernels->size);
         }
     }
@@ -885 +885 @@
                 psFree(stamp->image1);
                 psFree(stamp->image2);
-                psFree(stamp->weight);
-                stamp->image1 = stamp->image2 = stamp->weight = NULL;
+                psFree(stamp->variance);
+                stamp->image1 = stamp->image2 = stamp->variance = NULL;
                 psFree(stamp->matrix1);
                 psFree(stamp->matrix2);
@@ -989 +989 @@
 
 
-// XXX Put kernelImage, kernelWeight and polyValues on thread-dependent data
+// XXX Put kernelImage, kernelVariance and polyValues on thread-dependent data
 static bool subtractionConvolvePatch(int numCols, int numRows, // Size of image
                                      int x0, int y0, // Offsets for image
@@ -1010 +1010 @@
 
     psKernel *kernelImage = NULL;       // Kernel for the images
-    psKernel *kernelWeight = NULL;      // Kernel for the weight maps
+    psKernel *kernelVariance = NULL;      // Kernel for the variance maps
 
     // Only generate polynomial values every kernel footprint, since we have already assumed
@@ -1020 +1020 @@
 
     if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-        convolveRegion(out1->image, out1->weight, convMask, &kernelImage, &kernelWeight,
-                       ro1->image, ro1->weight, sys1, subMask, kernels, polyValues, background, *region,
+        convolveRegion(out1->image, out1->variance, convMask, &kernelImage, &kernelVariance,
+                       ro1->image, ro1->variance, sys1, subMask, kernels, polyValues, background, *region,
                        maskBad, maskPoor, poorFrac, useFFT, false);
     }
     if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-        convolveRegion(out2->image, out2->weight, convMask, &kernelImage, &kernelWeight,
-                       ro2->image, ro2->weight, sys2, subMask, kernels, polyValues, background, *region,
+        convolveRegion(out2->image, out2->variance, convMask, &kernelImage, &kernelVariance,
+                       ro2->image, ro2->variance, sys2, subMask, kernels, polyValues, background, *region,
                        maskBad, maskPoor, poorFrac, useFFT, kernels->mode == PM_SUBTRACTION_MODE_DUAL);
     }
 
     psFree(kernelImage);
-    psFree(kernelWeight);
+    psFree(kernelVariance);
     psFree(polyValues);
 
@@ -1148 +1148 @@
             }
         }
-        if (ro1->weight) {
-            if (!out1->weight) {
-                out1->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+        if (ro1->variance) {
+            if (!out1->variance) {
+                out1->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
                 if (threaded) {
-                    psMutexInit(out1->weight);
+                    psMutexInit(out1->variance);
                 }
             }
-            psImageInit(out1->weight, 0.0);
+            psImageInit(out1->variance, 0.0);
         }
     }
@@ -1165 +1165 @@
             }
         }
-        if (ro2->weight) {
-            if (!out2->weight) {
-                out2->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+        if (ro2->variance) {
+            if (!out2->variance) {
+                out2->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
                 if (threaded) {
-                    psMutexInit(out2->weight);
+                    psMutexInit(out2->variance);
                 }
             }
-            psImageInit(out2->weight, 0.0);
+            psImageInit(out2->variance, 0.0);
         }
     }
@@ -1232 +1232 @@
     psImage *polyValues = NULL;         // Pre-calculated polynomial values
     psKernel *kernelImage = NULL;       // Kernel for the images
-    psKernel *kernelWeight = NULL;      // Kernel for the weight maps
+    psKernel *kernelVariance = NULL;      // Kernel for the variance maps
 #endif
 
@@ -1345 +1345 @@
         if (out2) {
             out2->image = psMemIncrRefCounter(ro2->image);
-            out2->weight = psMemIncrRefCounter(ro2->weight);
+            out2->variance = psMemIncrRefCounter(ro2->variance);
             out2->mask = psMemIncrRefCounter(ro2->mask);
         }
@@ -1352 +1352 @@
         if (out1) {
             out1->image = psMemIncrRefCounter(ro1->image);
-            out1->weight = psMemIncrRefCounter(ro1->weight);
+            out1->variance = psMemIncrRefCounter(ro1->variance);
             out1->mask = psMemIncrRefCounter(ro1->mask);
         }