Index: /trunk/psLib/src/imageops/psImageInterpolate.c
===================================================================
--- /trunk/psLib/src/imageops/psImageInterpolate.c	(revision 18042)
+++ /trunk/psLib/src/imageops/psImageInterpolate.c	(revision 18043)
@@ -7,6 +7,6 @@
  *  @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,9 +312,13 @@
         #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,10 +546,16 @@
         #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,5 +586,5 @@
 
 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,14 +598,11 @@
         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");
     }
 
Index: /trunk/psModules/src/imcombine/pmSubtraction.c
===================================================================
--- /trunk/psModules/src/imcombine/pmSubtraction.c	(revision 18042)
+++ /trunk/psModules/src/imcombine/pmSubtraction.c	(revision 18043)
@@ -4,6 +4,6 @@
  *  @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,9 +50,18 @@
 
     // 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;