Changeset 18043 for trunk/psLib

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.

File:

• trunk/psLib/src/imageops/psImageInterpolate.c (modified) (5 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");
     }