Changeset 28796 for branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtraction.c

Timestamp:

Jul 30, 2010, 9:38:40 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/eam_branches/ipp-20100621/psModules

Files:

• . (modified) (1 prop)
• src/imcombine/pmSubtraction.c (modified) (8 diffs)

branches/eam_branches/ipp-20100621/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);