Changeset 28667 for trunk/psModules/src/imcombine/pmSubtraction.c

Timestamp:

Jul 13, 2010, 2:34:04 PM (16 years ago)

Author:

Paul Price

Message:

It's bothered me for a while that the noise model that we have been using is only kind of spatially variable. The variance map is intended to track the change in the noise properties as a function of position, and while it does keep track of different noise levels due to, e.g., different detectors, vignetting, etc., it has not tracked the different noise that comes from convolution kernels changing over an image. This had been an intentional oversight --- the prescription we were using for calculating the variance and covariance removed the spatial variation (by prescribing that the variance level under a convolution remain unchanged) and we were sweeping it under the carpet so we could calculate the covariance pseudo-matrix and use that for setting the absolute level of the noise.

I've now figured out how to have our cake and eat it too. By changing where the scalings are applied, the variance map can track the changing noise introduced by a varying convolution kernel, while the covariance pseudo-matrix tracks an absolute change. This fix is purely algorithmic --- it does not affect any of the variances/covariances stored on disk, nor is it affected by them. I hope that this means that any Magic recertification required is either minimal or completely unnecessary.

Again, all that has changed is where the scalings are applied.

File:

• trunk/psModules/src/imcombine/pmSubtraction.c (modified) (8 diffs)

trunk/psModules/src/imcombine/pmSubtraction.c

@@ -52 +52 @@
 
     // Take the square of the normal kernel
-    double sumVariance = 0.0; // Sum of the variance kernels
     for (int v = yMin; v <= yMax; v++) {
         for (int u = xMin; u <= xMax; u++) {
-            sumVariance += out->kernel[v][u] = PS_SQR(normalKernel->kernel[v][u]);
-        }
-    }
-
-    // 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 variance map
-    psBinaryOp(out->image, out->image, "*", psScalarAlloc(1.0 / sumVariance, PS_TYPE_F32));
-
+            out->kernel[v][u] = PS_SQR(normalKernel->kernel[v][u]);
+        }
+    }
     return out;
 }
@@ -271 +265 @@
 // Convolve an image using FFT
 static void convolveVarianceFFT(psImage *target,// Place the result in here
-                              psImage *variance, // Variance map to convolve
-                              psImage *kernelErr, // Kernel error image
-                              psImage *mask, // Mask image
-                              psImageMaskType maskVal, // Value to mask
-                              const psKernel *kernel, // Kernel by which to convolve
-                              psRegion region,// Region of interest
-                              int size        // Size of (square) kernel
-                              )
+                                psImage *variance, // Variance map to convolve
+                                psImage *kernelErr, // Kernel error image
+                                psImage *mask, // Mask image
+                                psImageMaskType maskVal, // Value to mask
+                                const psKernel *kernel, // Kernel by which to convolve
+                                psRegion region,// Region of interest
+                                int size        // Size of (square) kernel
+                                )
 {
     psRegion border = psRegionSet(region.x0 - size, region.x1 + size,
@@ -348 +342 @@
                                   psImage *image, // Image to convolve
                                   psImage *variance, // Variance map to convolve, or NULL
+                                  const psKernel *covar,               // Covariance, or NULL
                                   psImage *kernelErr, // Kernel error image, or NULL
                                   psImage *subMask, // Subtraction mask
@@ -393 +388 @@
         if (variance) {
             convolveDirect(convVariance, variance, *kernelVariance, region, 0.0, kernels->size);
+        }
+    }
+
+    if (variance && covar) {
+        // Apply covariance factor to variance map, to allow for spatial variation
+        float factor = psImageCovarianceCalculateFactor(*kernelImage, covar); // Factor to apply
+        for (int y = region.y0; y < region.y1; y++) {
+            for (int x = region.x0; x < region.x1; x++) {
+                convVariance->data.F32[y][x] *= factor;
+            }
         }
     }
@@ -1085 +1090 @@
     if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
         convolveRegion(out1->image, out1->variance, out1->mask, &kernelImage, &kernelVariance,
-                       ro1->image, ro1->variance, kernelErr1, subMask, kernels, polyValues, background,
-                       *region, maskBad, maskPoor, poorFrac, useFFT, false);
+                       ro1->image, ro1->variance, ro1->covariance, kernelErr1, 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->variance, out2->mask, &kernelImage, &kernelVariance,
-                       ro2->image, ro2->variance, kernelErr2, subMask, kernels, polyValues, background,
-                       *region, maskBad, maskPoor, poorFrac, useFFT,
+                       ro2->image, ro2->variance, ro2->covariance, kernelErr2, subMask, kernels,
+                       polyValues, background, *region, maskBad, maskPoor, poorFrac, useFFT,
                        kernels->mode == PM_SUBTRACTION_MODE_DUAL);
     }
@@ -1325 +1330 @@
 
     // Calculate covariances
-    // This can be fairly involved, so we only do it for a single instance
-    // Enable threads for covariance calculation, since we're not threading on top of it.
+    // This can be fairly involved, so we only do it for a small number of instances
     float position[NUM_COVAR_POS] = { -1.0, -0.5, 0.0, +0.5, +1.0 }; // Positions for covariance calculations
+    // Enable threads for covariance calculation, since we're not threading on top of it
     oldThreads = psImageCovarianceSetThreads(true);
     if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
@@ -1342 +1347 @@
         out1->covariance = psImageCovarianceAverage(covars);
         psFree(covars);
+        // Remove covariance factor from covariance, since we've put it in the variance map already
+        float factor = psImageCovarianceFactor(out1->covariance);
+        psBinaryOp(out1->covariance->image, out1->covariance->image, "/", psScalarAlloc(factor, PS_TYPE_F32));
     }
     if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
@@ -1356 +1364 @@
         out2->covariance = psImageCovarianceAverage(covars);
         psFree(covars);
+        // Remove covariance factor from covariance, since we've put it in the variance map already
+        float factor = psImageCovarianceFactor(out2->covariance);
+        psBinaryOp(out2->covariance->image, out2->covariance->image, "/", psScalarAlloc(factor, PS_TYPE_F32));
     }
     psImageCovarianceSetThreads(oldThreads);