Changeset 16351

Timestamp:

Feb 6, 2008, 6:03:22 PM (18 years ago)

Author:

Paul Price

Message:

I no longer believe that I should be normalising the variance kernel.
This isn't done by Andy Becker in hotpants, and I can't see it in the
math:

If I(x,y) is the image, V(x,y) is the variance map and k(u,v) is the
kernel, then the convolved image is: sum_u,v I(x-u,y-v).k(u,v). Then
the variance of the convolved image is: sum_u,v V(x-u,y-v).k(u,v)^{2.
If you want to normalise by the sum of the kernel, then the convolved
image is: sum_u,v I(x-u,y-v).k(u,v)/sum_u,v k(u,v) and the variance
is: sum_u,v V(x-u,y-v).k(u,v)}2/(sum_u,v k(u,v))^{2. But, of course in
practise you absorb the normalisation into the kernel itself, so you
get the first form. This has no division by the sum of the square.}

See also psModules/src/imcombine/pmSubtraction.c

File:

• trunk/psLib/src/imageops/psImageInterpolate.c (modified) (4 diffs)

trunk/psLib/src/imageops/psImageInterpolate.c

@@ -7 +7 @@
  *  @author Paul Price, IfA
  *
- *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-11-17 00:39:29 $
+ *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2008-02-07 04:03:22 $
  *
  *  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++) { \
-                        double kernel2 = PS_SQR(kernel[j][i]); /* Kernel squared */ \
-                        sumKernel2 += kernel2; \
-                        *varianceValue += kernel2 * variance->data.TYPE[yPix][xPix]; \
+                        *varianceValue += PS_SQR(kernel[j][i]) * variance->data.TYPE[yPix][xPix]; \
                     } \
                 } \
-                *varianceValue /= sumKernel2; /* Normalise so that sum of kernel squares is unity */ \
                 break; \
             }
@@ -546 +542 @@
         #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 xInterpValue = 0.0; /* Interpolation in x */ \
-                double xSumKernel2 = 0.0; /* Sum of kernel squared in x */ \
                 for (int i = 0, xPix = xCentral - (xNum - 1) / 2; i < xNum; i++, xPix++) { \
-                    double kernel2 = PS_SQR(xKernel[i]); /* Kernel squared */ \
-                    xSumKernel2 += kernel2; \
-                    xInterpValue += kernel2 * variance->data.TYPE[yPix][xPix]; \
+                    xInterpValue += PS_SQR(xKernel[i]) * variance->data.TYPE[yPix][xPix]; \
                 } \
-                double kernel2 = PS_SQR(yKernel[j]); /* Kernel squared */ \
-                ySumKernel2 += xSumKernel2 * kernel2; \
-                *varianceValue += xInterpValue * kernel2; /* Interpolating in y */ \
+                *varianceValue += xInterpValue * PS_SQR(yKernel[j]); /* Interpolating in y */ \
             } \
-            *varianceValue /= ySumKernel2; /* Normalise so that sum of kernel squares is unity */ \
             break; \
           }
@@ -604 +593 @@
         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");
-        }
+        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
+        // XXX these should probably be asserts, not PS_ASSERTS
     }