Changeset 18043

Timestamp:

Jun 9, 2008, 4:43:01 PM (18 years ago)

Author:

Paul Price

Message:

Gene has determined that what I had before --- normalising by the sum
of kernel squares for the variance --- is the correct thing to do, so
I'm restoring what I had before. Here's Gene's explanation:

since the errors in the output image are correlated by the smoothing
kernel, the output variance image actually should not change
normalization under the smoothing. thus, you should in fact the
normalizing to have the sum(psf^{2) = 1.0.
Note that in the output smoothed image, the per-pixel variance [sum of
(v[i][j] - <v>)}2] will decrease as you smooth, but the variance of
apertures comparable to your smoothing scale will not be decreased.
thus, the output weight image you generate is appropriate for
measuring extended parameters, but not per-pixel stats. If you are
making the assumption that the per-pixel variance should be consistent
with the mean of the variance images, you will have an error.

Location:

trunk

Files:

• psLib/src/imageops/psImageInterpolate.c (modified) (5 diffs)
• psModules/src/imcombine/pmSubtraction.c (modified) (2 diffs)

trunk/psLib/src/imageops/psImageInterpolate.c

@@ -7 +7 @@
  *  @author Paul Price, IfA
  *
- *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2008-02-07 04:03:22 $
+ *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2008-06-10 02:42:41 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -312 +312 @@
         #define KERNEL_VARIANCE_CASE(TYPE) \
             case PS_TYPE_##TYPE: { \
+                double sumKernel2 = 0.0; /* Sum of kernel squares */ \
                 for (int j = 0, yPix = yCentral - (yNum - 1) / 2; j < yNum; j++, yPix++) { \
                     for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) { \
-                        *varianceValue += PS_SQR(kernel[j][i]) * variance->data.TYPE[yPix][xPix]; \
+                        double kernel2 = PS_SQR(kernel[j][i]); /* Kernel squared */ \
+                        sumKernel2 += kernel2; \
+                        *varianceValue += kernel2 * variance->data.TYPE[yPix][xPix]; \
                     } \
                 } \
+                *varianceValue /= sumKernel2; /* Normalise so that sum of kernel squares is unity */ \
                 break; \
             }
@@ -542 +546 @@
         #define SEPARATE_VARIANCE_CASE(TYPE) \
           case PS_TYPE_##TYPE: { \
+            double ySumKernel2 = 0.0; /* Sum of kernel squared in y */ \
             for (int j = 0, yPix = yCentral - (yNum - 1) / 2; j < yNum; j++, yPix++) { \
+                double xSumKernel2 = 0.0; /* Sum of kernel squared in x */ \
                 double xInterpValue = 0.0; /* Interpolation in x */ \
                 for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) { \
-                    xInterpValue += PS_SQR(xKernel[i]) * variance->data.TYPE[yPix][xPix]; \
+                    double kernel2 = PS_SQR(xKernel[i]); /* Kernel squared */ \
+                    xSumKernel2 += kernel2; \
+                    xInterpValue += kernel2 * variance->data.TYPE[yPix][xPix]; \
                 } \
-                *varianceValue += xInterpValue * PS_SQR(yKernel[j]); /* Interpolating in y */ \
+                double kernel2 = PS_SQR(yKernel[j]); /* Kernel squared */ \
+                ySumKernel2 += xSumKernel2 * kernel2; \
+                *varianceValue += xInterpValue * kernel2; /* Interpolating in y */ \
             } \
             break; \
@@ -576 +586 @@
 
 psImageInterpolateStatus psImageInterpolate(double *imageValue, double *varianceValue, psMaskType *maskValue,
-                        float x, float y, const psImageInterpolateOptions *options)
+                                            float x, float y, const psImageInterpolateOptions *options)
 {
     PS_ASSERT_PTR_NON_NULL(options, PS_INTERPOLATE_STATUS_ERROR);
@@ -588 +598 @@
         PS_ASSERT_IMAGE_NON_NULL(variance, PS_INTERPOLATE_STATUS_ERROR);
         PS_ASSERT_IMAGE_TYPE(variance, image->type.type, PS_INTERPOLATE_STATUS_ERROR);
-        PS_ASSERT_IMAGES_SIZE_EQUAL(variance, image, PS_INTERPOLATE_STATUS_ERROR);
+        psAssert(image->numCols == variance->numCols && image->numRows == variance->numRows,
+                 "Image and variance sizes");
     }
     if (maskValue && mask) {
         PS_ASSERT_IMAGE_NON_NULL(mask, PS_INTERPOLATE_STATUS_ERROR);
         PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_MASK, PS_INTERPOLATE_STATUS_ERROR);
-        if ((image->numCols != mask->numCols) || (image->numRows != mask->numRows)) {
-          psAbort ("programming error");
-        }
-        PS_ASSERT_IMAGES_SIZE_EQUAL(mask, image, PS_INTERPOLATE_STATUS_ERROR);
-        // XXX these should probably be asserts, not PS_ASSERTS
+        psAssert(image->numCols == mask->numCols && image->numRows == mask->numRows, "Image and mask sizes");
     }
 

trunk/psModules/src/imcombine/pmSubtraction.c

@@ -4 +4 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.92 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2008-05-24 23:24:15 $
+ *  @version $Revision: 1.93 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2008-06-10 02:43:01 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -50 +50 @@
 
     // Take the square of the normal kernel
+    double sumNormal = 0.0, sumVariance = 0.0; // Sum of the normal and variance kernels
     for (int v = yMin; v <= yMax; v++) {
         for (int u = xMin; u <= xMax; u++) {
-            out->kernel[v][u] = PS_SQR(normalKernel->kernel[v][u]);
-        }
-    }
+            float value = normalKernel->kernel[v][u]; // Value of interest
+            float value2 = PS_SQR(value); // Value squared
+            sumNormal += value;
+            sumVariance += value2;
+            out->kernel[v][u] = value2;
+        }
+    }
+
+    // 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
+    psBinaryOp(out->image, out->image, "*", psScalarAlloc(PS_SQR(sumNormal) / sumVariance, PS_TYPE_F32));
 
     return out;