Changeset 21363 for trunk/psModules/src/imcombine/pmSubtraction.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/imcombine/pmSubtraction.c (modified) (21 diffs)

trunk/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);
@@ -915 +915 @@
 }
 
-psImage *pmSubtractionKernelImage(const pmSubtractionKernels *kernels, float x, float y, bool wantDual)
+psKernel *pmSubtractionKernel(const pmSubtractionKernels *kernels, float x, float y, bool wantDual)
 {
     PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NULL);
@@ -922 +922 @@
     PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, NULL);
 
-    // Precalulate polynomial values
-    psImage *polyValues = p_pmSubtractionPolynomial(NULL, kernels->spatialOrder, x, y);
-
-    // The appropriate kernel
-    psKernel *kernel = solvedKernel(NULL, kernels, polyValues, wantDual);
-
+    psImage *polyValues = p_pmSubtractionPolynomial(NULL, kernels->spatialOrder, x, y); // Solved polynomial
+    psKernel *kernel = solvedKernel(NULL, kernels, polyValues, wantDual); // The appropriate kernel
     psFree(polyValues);
 
+    return kernel;
+}
+
+psImage *pmSubtractionKernelImage(const pmSubtractionKernels *kernels, float x, float y, bool wantDual)
+{
+    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NULL);
+    PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NULL);
+    PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, NULL);
+    PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, NULL);
+
+    psKernel *kernel = pmSubtractionKernel(kernels, x, y, wantDual); // Convolution kernel
     psImage *image = psMemIncrRefCounter(kernel->image); // Image of the kernel
     psFree(kernel);
@@ -989 +996 @@
 
 
-// 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 +1017 @@
 
     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 +1027 @@
 
     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 +1155 @@
             // XXX }
         }
-        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);
                 // XXX if (threaded) {
-                // XXX     psMutexInit(out1->weight);
+                // XXX     psMutexInit(out1->variance);
                 // XXX }
             }
-            psImageInit(out1->weight, 0.0);
+            psImageInit(out1->variance, 0.0);
         }
     }
@@ -1165 +1172 @@
             // XXX }
         }
-        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);
                 // XXX if (threaded) {
-                // XXX     psMutexInit(out2->weight);
+                // XXX     psMutexInit(out2->variance);
                 // XXX }
             }
-            psImageInit(out2->weight, 0.0);
+            psImageInit(out2->variance, 0.0);
         }
     }
@@ -1232 +1239 @@
     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
 
@@ -1299 +1306 @@
                 psFree(job);
             } else {
-                subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2,
-                                         sys1, sys2, subMask, maskBad, maskPoor, poorFrac, subRegion,
-                                         kernels, doBG, useFFT);
+                subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2, sys1, sys2,
+                                         subMask, maskBad, maskPoor, poorFrac, subRegion, kernels, doBG,
+                                         useFFT);
             }
             psFree(subRegion);
@@ -1334 +1341 @@
         // XXX }
     }
-
     psImageConvolveSetThreads(oldThreads);
 
     psFree(sys1);
     psFree(sys2);
+
+    // Calculate covariances
+    // This can take a while, so we only do it for a single instance
+    // XXX psImageCovarianceCalculate could be multithreaded
+    if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
+        psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel
+        out1->covariance = psImageCovarianceCalculate(kernel, ro1->covariance);
+        psFree(kernel);
+    }
+    if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
+        psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0,
+                                               kernels->mode == PM_SUBTRACTION_MODE_DUAL); // Conv. kernel
+        out2->covariance = psImageCovarianceCalculate(kernel, ro2->covariance);
+        psFree(kernel);
+    }
 
     // Copy anything that wasn't convolved
@@ -1345 +1366 @@
         if (out2) {
             out2->image = psMemIncrRefCounter(ro2->image);
-            out2->weight = psMemIncrRefCounter(ro2->weight);
+            out2->variance = psMemIncrRefCounter(ro2->variance);
             out2->mask = psMemIncrRefCounter(ro2->mask);
+            out2->covariance = psMemIncrRefCounter(ro2->covariance);
         }
         break;
@@ -1352 +1374 @@
         if (out1) {
             out1->image = psMemIncrRefCounter(ro1->image);
-            out1->weight = psMemIncrRefCounter(ro1->weight);
+            out1->variance = psMemIncrRefCounter(ro1->variance);
             out1->mask = psMemIncrRefCounter(ro1->mask);
+            out1->covariance = psMemIncrRefCounter(ro1->covariance);
         }
         break;